Giter Club home page Giter Club logo

cern-physics's Introduction

$\rho'$ analysis notes

GIRD Selection criteria

These criteria were applied on the selection stage(GRID):

Events:

  • >= 4 tracks

Tracks:

  • Has Point On inner OR outer ITS Layer
  • Not ITS SA
  • |dca1| < 3 && |dca0| < 3;

Data info

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 110004 entries, 0 to 110003
Data columns (total 49 columns):
 #   Column            Non-Null Count   Dtype  
---  ------            --------------   -----  
 0   RunNum            110004 non-null  int32  
 1   PeriodNumber      110004 non-null  int64  
 2   OrbitNumber       110004 non-null  int64  
 3   BunchCrossNumber  110004 non-null  int64  
 4   Mass              110004 non-null  float32
 5   Pt                110004 non-null  float32
 6   Q                 110004 non-null  int32  
 7   Rapidity          110004 non-null  float32
 8   Phi               110004 non-null  float32
 9   ZNAenergy         110004 non-null  float32
 10  ZNCenergy         110004 non-null  float32
 11  ZPAenergy         110004 non-null  float32
 12  ZPCenergy         110004 non-null  float32
 13  VtxX              110004 non-null  float32
 14  VtxY              110004 non-null  float32
 15  VtxZ              110004 non-null  float32
 16  VtxContrib        110004 non-null  int32  
 17  VtxChi2           110004 non-null  float32
 18  VtxNDF            110004 non-null  float32
 19  SpdVtxX           110004 non-null  float32
 20  SpdVtxY           110004 non-null  float32
 21  SpdVtxZ           110004 non-null  float32
 22  SpdVtxContrib     110004 non-null  int32  
 23  V0Adecision       110004 non-null  int32  
 24  V0Cdecision       110004 non-null  int32  
 25  ADAdecision       110004 non-null  int32  
 26  ADCdecision       110004 non-null  int32  
 27  V0Afired          110004 non-null  bool   
 28  V0Cfired          110004 non-null  bool   
 29  ADAfired          110004 non-null  bool   
 30  ADCfired          110004 non-null  bool   
 31  STPfired          110004 non-null  bool   
 32  SMBfired          110004 non-null  bool   
 33  SM2fired          110004 non-null  bool   
 34  SH1fired          110004 non-null  bool   
 35  OM2fired          110004 non-null  bool   
 36  OMUfired          110004 non-null  bool   
 37  IsTriggered       110004 non-null  bool   
 38  nTracklets        110004 non-null  int32  
 39  nTracks           110004 non-null  int32  
 40  ZDCAtime_0        110004 non-null  float32
 41  ZDCAtime_1        110004 non-null  float32
 42  ZDCAtime_2        110004 non-null  float32
 43  ZDCAtime_3        110004 non-null  float32
 44  ZDCCtime_0        110004 non-null  float32
 45  ZDCCtime_1        110004 non-null  float32
 46  ZDCCtime_2        110004 non-null  float32
 47  ZDCCtime_3        110004 non-null  float32
 48  FORChip           110004 non-null  object 
dtypes: bool(11), float32(24), int32(10), int64(3), object(1)
memory usage: 18.8+ MB
<class 'pandas.core.frame.DataFrame'>
MultiIndex: 26978093 entries, (0, 0) to (110003, 3)
Data columns (total 13 columns):
 #   Column                  Dtype  
---  ------                  -----  
 0   T_NumberOfSigmaTPCPion  float32
 1   T_Eta                   float32
 2   T_Phi                   float32
 3   T_Px                    float32
 4   T_Py                    float32
 5   T_Pz                    float32
 6   T_Q                     int32  
 7   T_HasPointOnITSLayer0   bool   
 8   T_HasPointOnITSLayer1   bool   
 9   T_ITSModuleInner        int32  
 10  T_ITSModuleOuter        int32  
 11  T_TPCNCls               int32  
 12  T_TPCRefit              bool   
dtypes: bool(3), float32(6), int32(4)
memory usage: 1.2 GB

Transverse momenta of initial data set

<style scoped> .dataframe tbody tr th:only-of-type { vertical-align: middle; }
.dataframe tbody tr th {
    vertical-align: top;
}

.dataframe thead th {
    text-align: right;
}
</style>
count mean std min 25% 50% 75% max
zero net charge event $p_T$ 19456.0 0.837667 12.995872 0.000275 0.215188 0.483686 0.846086 1726.771484
non-zero net charge event $p_T$ 13762.0 1.163582 12.508871 0.010555 0.386721 0.648763 1.049536 1328.583740

png

C:\Users\bdrum\AppData\Roaming\Python\Python39\site-packages\mplhep\plot.py:238: RuntimeWarning: invalid value encountered in sqrt
  _yerr = np.sqrt(h)

png

Analysis cuts

Global tracks

It's known that global tracks consist from ITS and TPC identification, in our case we can implicitly add checks for TPC identification track and see what happens with the data.

Let's aply further conditions for the tracks:

  • |NumberOfSigmaTPCPion| < 3
  • Number of TPC Clusters > 50
  • TPCRefit

Low energy tracks and TPC

There is an idea about that tracks with small energies (low pt) not able to reach TPC. Idea is that addding such condition will decrease our signal and background level.

png

As we can see above adding of gobal tracks will decrease statistics level, but it is correct for both signal and background. Let's estimate what number of global tracks in 4 tracks is enough:

png

Here we see such construction each row contains three plots:

  • starting point or what we have (first cell correspond our initial plot at the very begining of this notebook)
  • what we will throw
  • the difference between 1 and 2

As we can see in case of transition from zero global tracks to one we will lose only background. I guess it's easy to see that the best case when we throw almost only background event is more than 2 global tracks in event.

In further analysis I will be use this case. 'ITS & (>= 2TPC)'

False triggering

There are some situations when CCUP9 trigger could be fired false: It may occured when some fake or random track fires FOR and trigger will provide.

We can check list of FORs of event and what chipkey has each of four tracks. In case it has matches and produce back to back topology this means correct trigger state.

img

See debugging details in one of the issue

Let's do the same thing as we do few cell above and let's try to understand what we will throw after splitting event by fake triggered or not:

png

Let's see how looks masses for correct and fake trigger events:

png

As we can see above the shape of masses for fake trigger and correct trigger the same. Here we use pdf view for a normalization purposes. Despite of this fact we saw on $p_t$ distribution that fake triggers contains signal part, but most of the fake trigger data is bacgrkound events.

In further analysis we will use only correct trigger events. Also we hope that it will make our data closer to MC quality.

ZDC cuts

ZDC again allows us to make signal more clear. Neutrons in ZDC could be a markers about peripheral events.

No handles with labels found to put in legend.

png

No handles with labels found to put in legend.

png

No handles with labels found to put in legend.

png

Beside energy distirbution we also have to make corrections for ZDC timing:

No handles with labels found to put in legend.

png

Event called passed ZDC cuts if satisfied such criteria:

  • For only one of the two sides
  • neutrons energy less or equal than given parameter ZDC_En value
  • or
  • one of ZDC times not in an interval of +- ZDC_Time_Delta parameter value
<style scoped> .dataframe tbody tr th:only-of-type { vertical-align: middle; }
.dataframe tbody tr th {
    vertical-align: top;
}

.dataframe thead th {
    text-align: right;
}
</style>
count mean std min 25% 50% 75% max
ZDC < 5TeV | abs(ZDC_times) >100 6605.0 0.408039 2.210536 0.000275 0.071012 0.281112 0.548197 176.975204
ZDC > 5TeV & abs(ZDC_times) <=100 3032.0 0.596154 0.887427 0.002845 0.259760 0.495125 0.785955 42.316540

png

Transversal momentum distribution after all cuts

<style scoped> .dataframe tbody tr th:only-of-type { vertical-align: middle; }
.dataframe tbody tr th {
    vertical-align: top;
}

.dataframe thead th {
    text-align: right;
}
</style>
count mean std min 25% 50% 75% max
orig $p_t$ 19456.0 0.837667 12.995872 0.000275 0.215188 0.483686 0.846086 1726.771484
$p_t$ all cuts 6605.0 0.408039 2.210536 0.000275 0.071012 0.281112 0.548197 176.975204

png

<style scoped> .dataframe tbody tr th:only-of-type { vertical-align: middle; }
.dataframe tbody tr th {
    vertical-align: top;
}

.dataframe thead th {
    text-align: right;
}
</style>
count mean std min 25% 50% 75% max
zero net charge event $p_t$ 6605.0 0.408039 2.210536 0.000275 0.071012 0.281112 0.548197 176.975204
non-zero net charge event $p_t$ 13762.0 1.163582 12.508871 0.010555 0.386721 0.648763 1.049536 1328.583740

png

Now let's try to see what tracks we lost from signal area and what contribution they have:

TPC and ITS has different coverage for polar angle:

img1

Perhaps we have tracks that not only can't reach TPC, but also has $\theta$ values that TPC doesn't cover.

Below we can see polar angle distribution for tracks that covers three cases:

  1. All tracks from events were reconstructed by ITS and TPC
  2. Only ITS tracks from events with only part TPC tracks. Here tracks that not reconstructed by TPC
  3. All tracks from events were reconstructed by ITS or TPC

We can see small gaps with for the second case, that allow to speak about correctness of the suggestion, but anyway low energy of tracks is the main reason why TPC can't reconstructed tracks.

Mass

Let's see on the mass distribution of the events

<style scoped> .dataframe tbody tr th:only-of-type { vertical-align: middle; }
.dataframe tbody tr th {
    vertical-align: top;
}

.dataframe thead th {
    text-align: right;
}
</style>
count mean std min 25% 50% 75% max
ft_zq_all_cuts 2418.0 1.613113 0.390436 0.740946 1.383315 1.549550 1.747386 4.907002
ft_nzq_pt_cut 663.0 1.413462 0.563853 0.637530 0.978157 1.344274 1.708745 4.250375

png

Pions subsystems

In our process 4 pions were producted. The most probably intermediate state including two pions and $\rho$ i.e. $$\rho' \rightarrow \rho \ \pi^+ \pi^- \rightarrow \pi^+ \pi^- \pi^+ \pi^-$$

We can see this on distribution of mass that can be obtained as all combinations of pairs from intial four tracks, i.e. only four pairs:

img

Here we can plot two distirbutions:

  1. Make all possible(4) combinations of pairs. Then take lightest and pair that belong to one combination with that. Plot masses of these two pairs.
  2. Plot masses of masses from possible combinations.
<style scoped> .dataframe tbody tr th:only-of-type { vertical-align: middle; }
.dataframe tbody tr th {
    vertical-align: top;
}

.dataframe thead th {
    text-align: right;
}
</style>
count mean std min 25% 50% 75% max
heavy 4836.0 0.812236 0.259066 0.291725 0.674686 0.778129 0.902377 3.544038
light 4836.0 0.525252 0.164924 0.279545 0.399250 0.498259 0.627641 2.129988

png

As we can see above, for second case (all possible pairs) we've got stronger signal in comparison with light-recoil pair as it made in STAR work. Let's build 2d distirbuition and marginals component separately:

png

Data modeling

First of all let's try to fit $\rho(0)$ that we can see as heavy pair of final 4 pions state:

png

C:\Users\bdrum\Anaconda3\envs\hep\lib\site-packages\lmfit\minimizer.py:857: RuntimeWarning: invalid value encountered in sqrt
  (par.stderr * np.sqrt(self.result.covar[jvar, jvar])))
C:\Users\bdrum\Anaconda3\envs\hep\lib\site-packages\lmfit\minimizer.py:850: RuntimeWarning: invalid value encountered in sqrt
  par.stderr = np.sqrt(self.result.covar[ivar, ivar])


[[Model]]
    (Model(bw, prefix='bw_') + Model(polynomial, prefix='bckg_'))
[[Fit Statistics]]
    # fitting method   = Nelder-Mead
    # function evals   = 766
    # data points      = 50
    # variables        = 9
    chi-square         = 40.7026433
    reduced chi-square = 0.99274740
    Akaike info crit   = 7.71350158
    Bayesian info crit = 24.9217086
##  Warning: uncertainties could not be estimated:
[[Variables]]
    bw_M:     0.76257570 +/- 0.00278473 (0.37%) (init = 0.77)
    bw_G:     0.17054290 +/- 0.00407858 (2.39%) (init = 0.17)
    bw_amp:   220.155833 +/- 0.07016377 (0.03%) (init = 230)
    bckg_c0:  375.237149 +/- 0.33459047 (0.09%) (init = 374.961)
    bckg_c1: -2511.03059 +/- 0.03314420 (0.00%) (init = -2493.255)
    bckg_c2:  7312.39507 +/- 0.03578421 (0.00%) (init = 7288.046)
    bckg_c3: -9574.79936 +/-        nan (nan%) (init = -9558.787)
    bckg_c4:  5677.15511 +/-        nan (nan%) (init = 5664.544)
    bckg_c5: -1250.88619 +/- 0.00275394 (0.00%) (init = -1244.663)
[[Correlations]] (unreported correlations are < 0.100)
    C(bw_G, bckg_c0) = -0.334

png

C:\Users\bdrum\Anaconda3\envs\hep\lib\site-packages\lmfit\minimizer.py:850: RuntimeWarning: invalid value encountered in sqrt
  par.stderr = np.sqrt(self.result.covar[ivar, ivar])
C:\Users\bdrum\Anaconda3\envs\hep\lib\site-packages\lmfit\minimizer.py:857: RuntimeWarning: invalid value encountered in sqrt
  (par.stderr * np.sqrt(self.result.covar[jvar, jvar])))


[[Model]]
    ((Model(bw, prefix='bw1_') + Model(bw, prefix='bw2_')) + Model(polynomial, prefix='bckg_'))
[[Fit Statistics]]
    # fitting method   = Nelder-Mead
    # function evals   = 1485
    # data points      = 30
    # variables        = 11
    chi-square         = 25.8808301
    reduced chi-square = 1.36214895
    Akaike info crit   = 17.5691549
    Bayesian info crit = 32.9823261
##  Warning: uncertainties could not be estimated:
[[Variables]]
    bw1_M:    1.38115837 +/-        nan (nan%) (init = 1.4)
    bw1_G:    0.26730196 +/-        nan (nan%) (init = 0.2)
    bw1_amp:  115.902720 +/-        nan (nan%) (init = 150)
    bw2_M:    1.58211568 +/-        nan (nan%) (init = 1.55)
    bw2_G:    0.28454224 +/-        nan (nan%) (init = 0.3)
    bw2_amp:  95.5333644 +/-        nan (nan%) (init = 135)
    bckg_c0:  135.876946 +/-        nan (nan%) (init = 159.0529)
    bckg_c1: -531.778166 +/-        nan (nan%) (init = -502.5231)
    bckg_c2:  693.204812 +/- 0.03743163 (0.01%) (init = 626.8528)
    bckg_c3: -341.989204 +/-        nan (nan%) (init = -321.2902)
    bckg_c4:  56.3015457 +/- 0.01892836 (0.03%) (init = 57.279)

png

Decay angle for $\rho_0$

Here we would like to study angle decay for $\rho_0 \rightarrow \pi^+\pi^-$ between the momentum of one of the pions in the rest frame of $\rho_0$ and original(lab frame) momentum of $\rho_0$.

For the sake of simplicity first we will create new lab system with OZ axis directed along momentum of $\rho_0$ in lab frame.

The transition from original lab frame to the new one could be obtain by two rotation:

  • First, around x axis:

Rotation matrix for this case: $$ R_{x,\alpha}= \Bigg( \begin{matrix} 1 & 0 & 0\ 0 & \cos{\alpha} & -\sin{\alpha} \ 0 & \sin{\alpha} & \cos{\alpha} \end{matrix} \Bigg) $$

Rotation angle for x axis:

$$ \begin{matrix} \sin{\alpha} = \frac{p_y}{\sqrt{p_y^2+p_z^2}} \\ \cos{\alpha} = \frac{p_z}{\sqrt{p_y^2+p_z^2}} \\ \end{matrix} $$

Now, momentum vector of $\rho_0$ in this system will looks like

$$ \boldsymbol{p'}=R_{x,\alpha}\boldsymbol{p}= \Bigg( \begin{matrix} p_x\\ p_y\cos{\alpha} - p_z\sin{\alpha}\\ p_y\sin{\alpha} + p_z\cos{\alpha}) \end{matrix} \Bigg) $$

  • Second, around y axis:

Rotation matrix for this case will looks like $$ R = R_{y,-\beta}= \Bigg( \begin{matrix} \cos{\beta} & 0 & -\sin{\beta}\ 0 & 1 & 0 \ \sin{\beta} & 0 & \cos{\beta} \end{matrix} \Bigg) $$

We have to apply it for already rotated vector around x axis:

$$ \boldsymbol{p''} =R_{y,-\beta}\boldsymbol{p'} = \Bigg( \begin{matrix} \cos{\beta} & 0 & -\sin{\beta}\ 0 & 1 & 0 \ \sin{\beta} & 0 & \cos{\beta} \end{matrix} \Bigg) \Bigg( \begin{matrix} {p'}_x\ {p'}_y\ {p'}_z \end{matrix} \Bigg)

\Bigg( \begin{matrix} {p'}_x\cos{\beta}-{p'}_z\sin{\beta}\ {p'}_y\ {p'}_x\sin{\beta}+{p'}_z\cos{\beta} \end{matrix} \Bigg) $$

Where rotation angle based on the new rotated vector $\boldsymbol{p'}$ coordinates:

$$ \begin{matrix} \sin{\beta} = \frac{{p'}_x}{\sqrt{{p'_x}^2+{p'}_z^2}} \\ \cos{\beta} = \frac{{p'}_z}{\sqrt{{p'_x}^2+{p'}_z^2}} \\ \end{matrix} $$

As a result final transition looks like:


$$ \boldsymbol{p''} = \Bigg( \begin{matrix} {p'}_x\cos{\beta}-{p'}_z\sin{\beta}\ {p'}_y\ {p'}_x\sin{\beta}+{p'}_z\cos{\beta} \end{matrix} \Bigg)= \Bigg( \begin{matrix} p_x\cos{\beta}-p_y\sin{\alpha}\sin{\beta}-p_z\cos{\alpha}\sin{\beta} \ p_y\cos{\alpha}-p_z\sin{\alpha} \ p_x\sin{\beta}+p_y\sin{\alpha}\cos{\beta}+p_z\cos{\alpha}\cos{\beta} \end{matrix} \Bigg) $$ where rotation angles are:

$$ \begin{matrix} \sin{\alpha} = \frac{p_y}{\sqrt{p_y^2+p_z^2}} \\ \cos{\alpha} = \frac{p_z}{\sqrt{p_y^2+p_z^2}} \\ \sin{\beta} = \frac{p_x}{\sqrt{p_x^2+(p_y\sin{\alpha}+p_z\cos{\alpha})^2}} \\ \cos{\beta} = \frac{p_y\sin{\alpha}+p_z\cos{\alpha}}{\sqrt{p_x^2+(p_y\sin{\alpha}+p_z\cos{\alpha})^2}} \\ \end{matrix} $$

Now let's consider moving coordinate system with $\rho_0$ so that OZ axis direct along $\boldsymbol{p}_{\rho_0}$

We know components of original momentum of $\pi^+$ in the such system and now let's boost their via Lorentz Transormation:

$$ \begin{matrix} \ {E'} = \gamma E - \Gamma p_z \\ \ {p'}_x= p_x \\ \ {p'}_y= p_y \\ \ {p'}_z= \gamma p_z - \Gamma E \\ \end{matrix} $$

where $$ \boldsymbol{\beta} = \frac{\boldsymbol{p}}{E} $$ $$ \gamma = \frac{E}{m}$$ $$ \Gamma = \gamma \beta = \frac{p}{m}$$

$$ \boldsymbol{\beta} = \frac{\boldsymbol{p}}{E} $$ $$ \gamma = \frac{1}{\sqrt{1-\frac{p^2}{E^2}}}$$

Now the searched angle can be obtain from scalar multiplication of $\pi^+$ momentum in the rest frame of $\rho_0$ and momentum of $\rho_0$ in the lab frame:

$$\cos{\theta}=\frac{\boldsymbol{{p'}{\pi^+}}\boldsymbol{p{\rho_0}}} {{p'}{\pi^+}p{\rho_0}}$$

Cross section

Let's see to cross section of my events.

For this we should take luminosity of runs.

Unfortunately file that I have to use for getting luminosity have a reference to special class AliTriggerInfo and moreover it packed into TObjArray, so I can't read it via uproot4. This is the reason why I used pure root again. Here is the script that I used.

Then let's see how much events do we have in each run.

$$L = \frac{1}{\sigma} \frac{\delta N}{\delta t}$$, this means that

$$\sigma \approx \frac{N}{L}$$

Cross section of phenomena should be flat and independent from runs.

Let's check it:

count    117.000000
mean       3.555185
std        1.101315
min        0.716081
25%        2.823498
50%        3.477399
75%        4.276723
max        6.551968
Name: sigma, dtype: float64

png

Other decays

In PDG I've seen also other interesting modes for $\rho'$:

    1. $\rho' \rightarrow \eta_0 \rho_0$ | ?
    • 1.1. $\rho_0 \rightarrow 4 \pi$ | $2*10^{-5}%$
    • 1.2. $\rho_0 \rightarrow \pi^+ \pi^-$ | $10^{-2}%$
    • 1.3. $\eta_0' \rightarrow \pi^+ \pi^- \gamma$ | $4%$
    • 1.4. $\eta_0' \rightarrow \pi^+ \pi^- \pi^0$ | $23%$
    1. $\rho' \rightarrow 4 \pi$ | ?

What if $\rho' \rightarrow \rho_0 \rho_0$ possible?

cern-physics's People

Contributors

bdrum avatar

Stargazers

 avatar

Watchers

 avatar

Forkers

abylinkin

cern-physics's Issues

Reorganizing project dir

Now in my project dir all things are together.

I prefer split notebooks, cpp project, and readme should be related with main notebook

Produce Monte Carlo data

Ask Michal Broz to run simulation for the:

Run anchored: Get Run number with maximum of 'good' events
Period: LHC15o
Process: rho(1720) → rho pi+ pi- → 4pi
Number of generated events: 100 000

GeneratorCode

AliGenerator* GeneratorCustom() 
{
  AliGenCocktail* genCocktail = GeneratorCocktail("StarlightEvtGen");
  
  AliGenStarLight* genStarLight = (AliGenStarLight*)GeneratorStarlight();
  genStarLight->SetParameter("PROD_PID     =   999    #Channel of interest (not relevant for photonuclear processes)");
  genStarLight->SetParameter("rho0PrimeMass     =   1.720");
  genStarLight->SetParameter("rho0PrimeWidth     =   0.249");
  genStarLight->SetParameter("W_MAX     =   5.0");
  genCocktail->AddGenerator(genStarLight, "Starlight", 1.);

  AliGenSLEvtGen *genEvtGen = new AliGenSLEvtGen();
  TString decayTable="RHOPRIME.RHOPIPI.DEC"; // default
  genEvtGen->SetUserDecayTable(gSystem->ExpandPathName(Form("$ALICE_ROOT/STARLIGHT/AliStarLight/DecayTables/%s",decayTable.Data())));
  genEvtGen->SetEvtGenPartNumber(30113);

  genEvtGen->SetPolarization(1);           // Set polarization: transversal(1),longitudinal(-1),unpolarized(0)

  genCocktail->AddGenerator(genEvtGen,"EvtGen",1.);
  return genCocktail;
}

Here is the question about decay channel. Now it specified as 999 that correspond to 4 pions final state, but in startlight one more exists - 913 that produces a final state consisting of two pions, also it adds the amplitudes (with interference) for photoproduction of the rho, and for photoproduction of direct pi^+pi^-., but again the final state is only two pions.

Comparison runs quality

I've processed 125 runs, but only ~80 were recommended by DPG.

I can compare the quality of such runs and the rest via visible c-s.

Unbinned fitting (Likelihood)

For fitting in #44 I have standard $\chi^2$ test for deterring quality of the fit, but for area with small statistics (Puasson) it's wrong.
E.g. I've calculated error as $\sqrt(N)$ for each bin, but in case of small count of events in bin it provide us the situation when 2 events in bin and 1.4 sigma.

Here is the fragment from lection of Cambridge prof. Mark Thomson about the question:

image

Move code from cell to scripts

Keeping a code in jupyter notebook cells make notebook oversized.
I think the good solution is keeping scripts in the 'scripts' folder and use their as modules.

Generalize computations

Now, when I want to plot mass with total charge zero and non-zero, I can't do it because I don't have a functions or objects that could change the state. So I should repeat a code. This is bad solution.

I have to parametrize calculations and states.

Comparison with Valeri data

Now me and Valeri have discrepancies in the data.

He has less number of processed runs, but more statistics than me.

What is the reason of such discrepancies?

He'll provide to me file with event info and we can find out the difference...

Discrepancies in the data

In my analysis two the same criteria, but exposed on the different manner gives difference result:

selectOnlyStandard    = data['T_TPCRefit']  * (data['T_TPCNCls'] > 50) * \
            (data['T_ITSNCls'] > 3) * \
            (~newT_ITSsa) * (np.abs(data['T_NumberOfSigmaTPCPion']) < 3)
GoodEventsStandard = np.argwhere(selectOnlyStandard.sum()==4).flatten() # get events with 4 good tracks
GoodEventsStandard = GoodEventsStandard[np.argwhere(data['T_Q'][selectOnlyStandard][GoodEventsStandard].sum()==0).flatten()].flatten()  

pxstd = data['T_Px'][selectOnlyStandard][GoodEventsStandard]
pystd = data['T_Py'][selectOnlyStandard][GoodEventsStandard]
    
ptstd = np.sqrt(pxstd.sum()**2  + pystd.sum()**2)

fig = plt.figure(figsize=(15, 7))
ax = fig.add_axes([0,0,1,1])
fig.suptitle(f'4pr pt standard criteria', fontsize=32)
plt.style.use(hep.style.ROOT)
counts, bins = np.histogram(ptstd, bins=100, range=(0,2))
_ = ax.hist(ptstd, bins=bins, color='black', histtype='step', label=f'Q = 0;Entries {np.sum(counts)}')
plt.xlabel('Pt, GeV')
plt.ylabel('# events')
ax.legend()

GoodEvents = GetGoodEvents(WithGoodNTpcTracks=4)
counts, bins = np.histogram(GetPt(Draw=False), bins=100, range=(0,2))
_ = ax.hist(GetPt(Draw=False), bins=bins, color='red', histtype='step', label=f'4 TPC tracks;Entries {np.sum(counts)}', linewidth=2)
plt.xlabel('Pt, GeV')
plt.ylabel('# events')
ax.legend()

gives:

First criteria direct combination of second criteria.

Obviously for direct criteria condition to total charge doesn't work:

We can see the sum of zero charge events and non zero:

Picture for the second case, but the meaning is understood

What is the best way to access UPC events?

In the slides of the one presentation I've seen:

We have two types of UPC events Triggered by MUON Triggered in the Central Barrel In both cases the events are almost empty The average size of AOD UPC events in the central barrel for 2015 PbPb is around 0.5 kB We have some 30 M triggered events, which means some 15 GB of data. The purity of this sample is small, meaning that a simple pre-selection may reduce this number by a large factor (between 2 and 3).

In the past we have used a LEGO train to access the events and produced trees to perform local analyses. The code in the LEGO train is quite stable Latchezar commented that this approach may not be optimal for such small data set distributed over all the PbPb data set and recommended that we had a look at nanoAODs.

The code to select only those branches we use and create with them ‘a nanoAOD’ is ready. As it is so general, we do not expect the code to be changed frequently. Latchezar mentioned that there exist AOD trains in the central framework.

Could it be a good option to use them to produce UPC AODs?

In Summary, what is the best way to run over the UPC triggers in the PbPb sample so that we do not have to run over all the PbPb events to find our few gigas of data?

ZDC energy distribution

Plot histogram of zdc energy distribution.

I've to see smth like abs(sinc)....

Also demonstrate that basic part (let say first half of the plot) contains events that form Mass spectra meanwhile second half mass spectra contains background events.

Move data preparation from pandas to numpy

Till this moment I've used uproot4 for reading root files and output result was pandas dataframe.
Now my file has grown from 500mb to 1.5gb, so current version of script doesn't work because it takes forever and takes all memory.

I've tested parsing to numpy arrays and it works fine.

Now I have to move the main script form pandas usage to numpy.

Aliases for most common visualizing scenarios

I have many common scenarios for data visualizing.

E.g.

  • vertical comparison or horizontal comparison on the different plots
  • adding data to the same plots with labels
  • histogramming with errors bars
  • custom styles for histograms

How can I simplify usage of such behavior?

First of all I know that in matplotlib I can specify default parameters for some functions. I can use this e.g. for changing default type of histograms to 'step', change defalut colors, applying root style as default.

Then I can prepare functions that can be wrap described common scenarios with additional parameters like arrays of colors and labels, and extend standard arguments of matplotlib via args and kwargs.

Influence of ITSNCls >=3 or >3

I've seen that expanding condition from >3 to >=3 for ITSNCls parameter gives me about 600 events to statistics.

What are these events, background or signal?

Overlapping problem

The same track can be triggered with different names (CCUP2 and CCUP9 and so on).

I have to estimate this influence.

Mass from pairs

I've got mass from 4 tracks, but what if I will calculate it from pion subsystems?

Events losses after splitting by trigger

I faced a problem that after splitting event to triggered and untriggered I have a difference between origin total event count and
triggered event plus untriggered:

ind nTPC total_cnt trig_cnt untrig_cnt trig % untrig % trig_cnt_plus_untrig_cnt diff
0 0 19194 10349 8631 53.9 45 18980 214
1 1 17641 10138 7334 57.5 41.6 17472 169
2 2 16506 9572 6780 58 41.1 16352 154
3 3 12952 7583 5249 58.5 40.5 12832 120
4 4 5916 3506 2359 59.3 39.9 5865 51

Definitely this is wrong

New structure for analysis scripts

In the solution process #22 I realized that a problem of structure is not only in oversized. Also reinitializing of variables in case of restarting the kernel could be a headache. Moreover jumping from cell to cell without tags is really slowly.

I've a plan to create package for my analysis based on python modules and in notebooks just import modules and call functions.

Histogramming problems

ROOT has excellent histogram tools with different visualizations modes (error bars, markers), operations for histogram such as ratio or difference.

What I need is

  • Stats table, such as entries, mean, std, overflow, integral
  • Operations with histograms such as division and difference
  • Histogram with Error bars
  • Visualization with comparison purposes it could be n hist at one plot, subplots, so on

I've seen many of histogramming tools including ROOT, e.g.
scikit-hep/hist
scikit-hep/aghast
scikit-hep/boost-histogram
numpy.histogram

Perhaps many of these tick could be easy solved via let say standard packages, e.g. stats could be provided by pd.describe (except overflow and integral)

Migrate drawing to pandas

Now I use such constructions for drawing:

bw, = ax.plot(np.linspace(x.min(), x.max(),1000), fit.bw(x=np.linspace(x.min(), x.max(),1000), M = result.best_values['bw_M'], G=result.best_values['bw_G'], amp=result.best_values['bw_amp'] ), 'g--',label=r'Breit-Wigner')
bkg, = plt.plot(np.linspace(x.min(), x.max(),1000), fit.polyn(np.linspace(x.min(), x.max(),1000), result.best_values['bckg_c0'], result.best_values['bckg_c1'], result.best_values['bckg_c2'], result.best_values['bckg_c3'], result.best_values['bckg_c4'],0,0,0, result.best_values['bckg_amp']), 'b--', label='Bckgr fit')
# plt.plot(bcg_x, bcg_y, 'r+', label='bkg data')
result.plot_fit(datafmt='ko',fitfmt='r-',numpoints=100,data_kws={'xerr':0.022})
ax.set_title('')
_ = ax.set_ylabel(r'Counts per 20 Mev/$c^2$')
_ = ax.set_xlabel(r'M$(\pi^+\pi^-\pi^+\pi^-)$ (GeV/$c^2$)')
_ = plt.legend([dat,ft, bw, bkg],['Data','Fit','Breit-Wigner','Bkg'])

Instead of this, I can combine data to DataFrames and use pandas plot for drawing purposes

Influence of ITS only tracks to 4 tracks events

Let's see how does further parameters for tracks are distributed for TPC+ITS and ITS only

  • full momentum
  • pt
  • pseudorapidity

Mass based on ITS only tracks should be biased to the soft side.

image

Cross section

Now I have number of events with $\rho'$ for each runs and I can get luminosity of these runs.

This means that I can get cross section of the $\rho'$

In fact this is not true, but why?

$$\sigma=\frac{N}{L}$$

Obviously this should be flat. But for such definition I can't confirm that:

img

Actually in case we will take std as $\frac{1}{\sqrt{nev}}$ such calculation looks strange because of scale doesn't match.

Perhaps, we have to normalize the data. In this case we can see flat curve(last picture):

img

What eta distribution for pion subsytems?

The idea is that we can check $\eta$ distribution for the pairs and check, perhaps lightest mass condition for selecting pair is to hard and we can increase accuracy if make condition for $\eta$

Fit mass of resonance

Here is ρ' resonance.

img

According with PDG description
The best fit could be obtained with the hypothesis of ρ-like resonances at 1420 and 1770 MeV, with widths of about 250 MeV.

I can create two BW model and fit resonance by them.

Trigger influence investigation

Estimate triggers influence: selected triggers are CCUP9, CCUP2, CCUP4, C1ZED, CCUP11. what difference between them according for my events?

Investigate other decay modes for $\rho'$

Until now I was studying only one decay mode:

$\rho' \rightarrow 4 \pi$

In PDG I've seen also other interesting modes for $\rho'$:

    1. $\rho' \rightarrow \eta_0 \rho_0$ | ?
    • 1.1. $\rho_0 \rightarrow 4 \pi$ | $2*10^{-5}%$
    • 1.2. $\rho_0 \rightarrow \pi^+ \pi^-$ | $10^{-2}%$
    • 1.3. $\eta_0' \rightarrow \pi^+ \pi^- \gamma$ | $4%$
    • 1.4. $\eta_0' \rightarrow \pi^+ \pi^- \pi^0$ | $23%$
    1. $\rho' \rightarrow 4 \pi$ | ?

What about $\rho' \rightarrow \rho_0 \rho_0$ is it possible?

$\rho_0$ decay angle for the rest frame of $\rho_0$

4tracks_analysis

During 4$\pi$ photoproduction in reaction of the $\rho'$ decay:

$$\rho'\rightarrow \pi^+\pi^-\pi^+\pi^-$$

I have a intermediate two-body decay state when $\rho_0\rightarrow \pi^+\pi^-$.

This allows us to check $\rho_0$ via determining of the decay angle for rest frame of $\rho_0$.

Track duplicates with only ITS criteria

I see in my events plenty of duplicated tracks:

There are only in case of ITS selection

((data['T_HasPointOnITSLayer0']) + (data['T_HasPointOnITSLayer1'])) * \
data['T_ITSRefit']  * (np.abs(data['T_NumberOfSigmaITSPion']) < 3)

img

In case of adding TPC criteria looks like (cause of I can't get it) that such tracks will be thrown...

4tr mass distr

Plot the Mass for 4 prongs events with zero and non-zero total charge for different pt intervals:

  • (0,150]
  • (150,800]
  • (800,2000)

Influence of limitation for tracks number

In the selection script PVN limits the number of THE GOOD(passed by some criteria) TRACKStracks and fill predefined arrays.

I have chosen another strategy and fill vectors with dynamic size.

That's interesting to compare data losses in PVN case.

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. 📊📈🎉

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google ❤️ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.