]>
Commit | Line | Data |
---|---|---|
6ad0bfa0 | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | ||
b2a60966 | 16 | /* $Id$ */ |
17 | ||
6ad0bfa0 | 18 | //_________________________________________________________________________ |
b2a60966 | 19 | // Implementation version v1 of the PHOS particle identifier |
7acf6008 | 20 | // Particle identification based on the |
148b2bba | 21 | // - RCPV: distance from CPV recpoint to EMCA recpoint. |
22 | // - TOF | |
23 | // - PCA: Principal Components Analysis.. | |
24 | // The identified particle has an identification number corresponding | |
25 | // to a 9 bits number: | |
bc0c084c | 26 | // -Bit 0 to 2: bit set if RCPV > CpvEmcDistance (each bit corresponds |
148b2bba | 27 | // to a different efficiency-purity point of the photon identification) |
bc0c084c | 28 | // -Bit 3 to 5: bit set if TOF < TimeGate (each bit corresponds |
148b2bba | 29 | // to a different efficiency-purity point of the photon identification) |
30 | // -Bit 6 to 9: bit set if Principal Components are | |
50739f15 | 31 | // inside an ellipse defined by the parameters a, b, c, x0 and y0. |
148b2bba | 32 | // (each bit corresponds to a different efficiency-purity point of the |
50739f15 | 33 | // photon identification) |
34 | // The PCA (Principal components analysis) needs a file that contains | |
35 | // a previous analysis of the correlations between the particles. This | |
bc0c084c | 36 | // file is $ALICE_ROOT/PHOS/PCA8pa15_0.5-100.root. Analysis done for |
50739f15 | 37 | // energies between 0.5 and 100 GeV. |
9fa5f1d0 | 38 | // A calibrated energy is calculated. The energy of the reconstructed |
50739f15 | 39 | // cluster is corrected with the formula A + B * E + C * E^2, whose |
bc0c084c | 40 | // parameters where obtained through the study of the reconstructed |
50739f15 | 41 | // energy distribution of monoenergetic photons. |
a4e98857 | 42 | // |
bc0c084c | 43 | // All the parameters (RCPV(2 rows-3 columns),TOF(1r-3c),PCA(5r-4c) |
50739f15 | 44 | // and calibration(1r-3c))are stored in a file called |
45 | // $ALICE_ROOT/PHOS/Parameters.dat. Each time that AliPHOSPIDv1 is | |
bc0c084c | 46 | // initialized, this parameters are copied to a Matrix (9,4), a |
50739f15 | 47 | // TMatrixD object. |
7acf6008 | 48 | // |
a4e98857 | 49 | // use case: |
50739f15 | 50 | // root [0] AliPHOSPIDv1 * p = new AliPHOSPIDv1("galice1.root") |
a4e98857 | 51 | // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated |
50739f15 | 52 | // // reading headers from file galice1.root and create RecParticles |
53 | // TrackSegments and RecPoints are used | |
54 | // // set file name for the branch RecParticles | |
f0a4c9e9 | 55 | // root [1] p->ExecuteTask("deb all time") |
50739f15 | 56 | // // available options |
57 | // // "deb" - prints # of reconstructed particles | |
58 | // // "deb all" - prints # and list of RecParticles | |
59 | // // "time" - prints benchmarking results | |
7acf6008 | 60 | // |
50739f15 | 61 | // root [2] AliPHOSPIDv1 * p2 = new AliPHOSPIDv1("galice1.root","v1",kTRUE) |
148b2bba | 62 | // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated |
50739f15 | 63 | // //Split mode. |
f0a4c9e9 | 64 | // root [3] p2->ExecuteTask() |
65 | // | |
50739f15 | 66 | |
f0a4c9e9 | 67 | |
7acf6008 | 68 | //*-- Author: Yves Schutz (SUBATECH) & Gines Martinez (SUBATECH) & |
148b2bba | 69 | // Gustavo Conesa April 2002 |
50739f15 | 70 | // PCA redesigned by Gustavo Conesa October 2002: |
71 | // The way of using the PCA has changed. Instead of 2 | |
72 | // files with the PCA, each one with different energy ranges | |
73 | // of application, we use the wide one (0.5-100 GeV), and instead | |
bc0c084c | 74 | // of fixing 3 ellipses for different ranges of energy, it has been |
50739f15 | 75 | // studied the dependency of the ellipses parameters with the |
76 | // energy, and they are implemented in the code as a funtion | |
77 | // of the energy. | |
78 | // | |
79 | // | |
80 | // | |
6ad0bfa0 | 81 | // --- ROOT system --- |
c947e71a | 82 | |
83 | ||
84 | // --- Standard library --- | |
85 | #include "TF1.h" | |
7acf6008 | 86 | #include "TBenchmark.h" |
148b2bba | 87 | #include "TPrincipal.h" |
c947e71a | 88 | #include "TFile.h" |
e3817e5f | 89 | #include "TSystem.h" |
148b2bba | 90 | |
6ad0bfa0 | 91 | // --- AliRoot header files --- |
c947e71a | 92 | //#include "AliLog.h" |
7acf6008 | 93 | #include "AliGenerator.h" |
e3817e5f | 94 | #include "AliPHOS.h" |
26d4b141 | 95 | #include "AliPHOSPIDv1.h" |
7b7c1533 | 96 | #include "AliPHOSGetter.h" |
6ad0bfa0 | 97 | |
26d4b141 | 98 | ClassImp( AliPHOSPIDv1) |
6ad0bfa0 | 99 | |
1cb7c1ee | 100 | //____________________________________________________________________________ |
101 | AliPHOSPIDv1::AliPHOSPIDv1():AliPHOSPID() | |
102 | { | |
a4e98857 | 103 | // default ctor |
148b2bba | 104 | |
8d0f3f77 | 105 | InitParameters() ; |
92f521a9 | 106 | fDefaultInit = kTRUE ; |
7acf6008 | 107 | } |
108 | ||
581354c5 | 109 | //____________________________________________________________________________ |
88cb7938 | 110 | AliPHOSPIDv1::AliPHOSPIDv1(const AliPHOSPIDv1 & pid ):AliPHOSPID(pid) |
581354c5 | 111 | { |
386aef34 | 112 | // ctor |
581354c5 | 113 | InitParameters() ; |
581354c5 | 114 | Init() ; |
581354c5 | 115 | |
116 | } | |
117 | ||
7acf6008 | 118 | //____________________________________________________________________________ |
88cb7938 | 119 | AliPHOSPIDv1::AliPHOSPIDv1(const TString alirunFileName, const TString eventFolderName):AliPHOSPID(alirunFileName, eventFolderName) |
7acf6008 | 120 | { |
a4e98857 | 121 | //ctor with the indication on where to look for the track segments |
7b7c1533 | 122 | |
8d0f3f77 | 123 | InitParameters() ; |
2bd5457f | 124 | Init() ; |
92f521a9 | 125 | fDefaultInit = kFALSE ; |
7acf6008 | 126 | } |
7b7c1533 | 127 | |
7acf6008 | 128 | //____________________________________________________________________________ |
129 | AliPHOSPIDv1::~AliPHOSPIDv1() | |
130 | { | |
79bb1b62 | 131 | // dtor |
9fa5f1d0 | 132 | |
e3817e5f | 133 | delete [] fX ; // Principal input |
134 | delete [] fPPhoton ; // Photon Principal components | |
135 | delete [] fPPi0 ; // Pi0 Principal components | |
7acf6008 | 136 | } |
a496c46c | 137 | //____________________________________________________________________________ |
138 | const TString AliPHOSPIDv1::BranchName() const | |
139 | { | |
88cb7938 | 140 | |
141 | return GetName() ; | |
a496c46c | 142 | } |
143 | ||
148b2bba | 144 | //____________________________________________________________________________ |
145 | void AliPHOSPIDv1::Init() | |
146 | { | |
147 | // Make all memory allocations that are not possible in default constructor | |
148 | // Add the PID task to the list of PHOS tasks | |
a496c46c | 149 | |
adcca1e6 | 150 | AliPHOSGetter * gime = AliPHOSGetter::Instance() ; |
151 | if(!gime) | |
152 | gime = AliPHOSGetter::Instance(GetTitle(), fEventFolderName.Data()) ; | |
88cb7938 | 153 | |
154 | if ( !gime->PID() ) | |
155 | gime->PostPID(this) ; | |
148b2bba | 156 | } |
8d0f3f77 | 157 | |
158 | //____________________________________________________________________________ | |
159 | void AliPHOSPIDv1::InitParameters() | |
160 | { | |
e3817e5f | 161 | // Initialize PID parameters |
adcca1e6 | 162 | fWrite = kTRUE ; |
8d0f3f77 | 163 | fRecParticlesInRun = 0 ; |
8d0f3f77 | 164 | fNEvent = 0 ; |
8d0f3f77 | 165 | fRecParticlesInRun = 0 ; |
35adb638 | 166 | fBayesian = kTRUE ; |
9fa5f1d0 | 167 | SetParameters() ; // fill the parameters matrix from parameters file |
eabde521 | 168 | SetEventRange(0,-1) ; |
35adb638 | 169 | |
cc1fe362 | 170 | // initialisation of response function parameters |
171 | // Tof | |
172 | // Photons | |
35adb638 | 173 | fTphoton[0] = 0.218 ; |
174 | //fTphoton[0] = 1. ; | |
175 | fTphoton[1] = 1.55E-8 ; | |
176 | fTphoton[2] = 5.05E-10 ; | |
177 | fTFphoton = new TFormula("ToF response to photons" , "gaus") ; | |
cc1fe362 | 178 | fTFphoton->SetParameters( fTphoton[0], fTphoton[1], fTphoton[2]) ; |
35adb638 | 179 | // // Electrons |
180 | // fTelectron[0] = 0.2 ; | |
181 | // fTelectron[1] = 1.55E-8 ; | |
182 | // fTelectron[2] = 5.35E-10 ; | |
183 | // fTFelectron = new TFormula("ToF response to electrons" , "gaus") ; | |
184 | // fTFelectron->SetParameters( fTelectron[0], fTelectron[1], fTelectron[2]) ; | |
185 | // // Muons | |
186 | // fTmuon[0] = 0.2 ; | |
187 | // fTmuon[1] = 1.55E-8 ; | |
188 | // fTmuon[2] = 5.1E-10 ; | |
189 | // fTFmuon = new TFormula("ToF response to muons" , "gaus") ; | |
190 | // fTFmuon->SetParameters( fTmuon[0], fTmuon[1], fTmuon[2]) ; | |
191 | ||
192 | // Pions | |
193 | //Gaus (0 to max probability) | |
194 | fTpiong[0] = 0.0971 ; | |
195 | //fTpiong[0] = 1. ; | |
196 | fTpiong[1] = 1.58E-8 ; | |
197 | fTpiong[2] = 5.69E-10 ; | |
198 | fTFpiong = new TFormula("ToF response to pions" , "gaus") ; | |
199 | fTFpiong->SetParameters( fTpiong[0], fTpiong[1], fTpiong[2]) ; | |
200 | // Landau (max probability to inf) | |
201 | // fTpionl[0] = 0.05 ; | |
202 | // //fTpionl[0] = 5.53 ; | |
203 | // fTpionl[1] = 1.68E-8 ; | |
204 | // fTpionl[2] = 5.38E-10 ; | |
205 | // fTFpionl = new TFormula("ToF response to pions" , "landau") ; | |
206 | // fTFpionl->SetParameters( fTpionl[0], fTpionl[1], fTpionl[2]) ; | |
207 | ||
208 | ||
209 | // Kaons | |
210 | //Gaus (0 to max probability) | |
211 | fTkaong[0] = 0.0542 ; | |
212 | //fTkaong[0] = 1. ; | |
213 | fTkaong[1] = 1.64E-8 ; | |
214 | fTkaong[2] = 6.07-10 ; | |
215 | fTFkaong = new TFormula("ToF response to kaon" , "gaus") ; | |
216 | fTFkaong->SetParameters( fTkaong[0], fTkaong[1], fTkaong[2]) ; | |
217 | //Landau (max probability to inf) | |
218 | fTkaonl[0] = 0.264 ; | |
219 | //fTkaonl[0] = 5.53 ; | |
220 | fTkaonl[1] = 1.68E-8 ; | |
221 | fTkaonl[2] = 4.10E-10 ; | |
222 | fTFkaonl = new TFormula("ToF response to kaon" , "landau") ; | |
223 | fTFkaonl->SetParameters( fTkaonl[0], fTkaonl[1], fTkaonl[2]) ; | |
224 | ||
225 | //Heavy Hadrons | |
226 | //Gaus (0 to max probability) | |
227 | fThhadrong[0] = 0.0302 ; | |
228 | //fThhadrong[0] = 1. ; | |
229 | fThhadrong[1] = 1.73E-8 ; | |
230 | fThhadrong[2] = 9.52E-10 ; | |
231 | fTFhhadrong = new TFormula("ToF response to heavy hadrons" , "gaus") ; | |
232 | fTFhhadrong->SetParameters( fThhadrong[0], fThhadrong[1], fThhadrong[2]) ; | |
233 | //Landau (max probability to inf) | |
234 | fThhadronl[0] = 0.139 ; | |
235 | //fThhadronl[0] = 5.53 ; | |
236 | fThhadronl[1] = 1.745E-8 ; | |
237 | fThhadronl[2] = 1.00E-9 ; | |
238 | fTFhhadronl = new TFormula("ToF response to heavy hadrons" , "landau") ; | |
239 | fTFhhadronl->SetParameters( fThhadronl[0], fThhadronl[1], fThhadronl[2]) ; | |
240 | ||
241 | /// /gaussian parametrization for pions | |
242 | // fTpion[0] = 3.93E-2 ; fTpion[1] = 0.130 ; fTpion[2] =-6.37E-2 ;//constant | |
243 | // fTpion[3] = 1.65E-8 ; fTpion[4] =-1.40E-9 ; fTpion[5] = 5.96E-10;//mean | |
244 | // fTpion[6] = 8.09E-10; fTpion[7] =-4.65E-10; fTpion[8] = 1.50E-10;//sigma | |
245 | ||
246 | // //landau parametrization for kaons | |
247 | // fTkaon[0] = 0.107 ; fTkaon[1] = 0.166 ; fTkaon[2] = 0.243 ;//constant | |
248 | // fTkaon[3] = 1.80E-8 ; fTkaon[4] =-2.96E-9 ; fTkaon[5] = 9.60E-10;//mean | |
249 | // fTkaon[6] = 1.37E-9 ; fTkaon[7] =-1.80E-9 ; fTkaon[8] = 6.74E-10;//sigma | |
250 | ||
251 | // //landau parametrization for nucleons | |
252 | // fThhadron[0] = 6.33E-2 ; fThhadron[1] = 2.52E-2 ; fThhadron[2] = 2.16E-2 ;//constant | |
253 | // fThhadron[3] = 1.94E-8 ; fThhadron[4] =-7.06E-10; fThhadron[5] =-4.69E-10;//mean | |
254 | // fThhadron[6] = 2.55E-9 ; fThhadron[7] =-1.90E-9 ; fThhadron[8] = 5.41E-10;//sigma | |
255 | ||
256 | ||
257 | // Shower shape: dispersion gaussian parameters | |
258 | // Photons | |
259 | ||
c947e71a | 260 | // fDphoton[0] = 3.84e-2; fDphoton[1] = 4.46e-3 ; fDphoton[2] = -2.36e-2;//constant |
261 | // //fDphoton[0] = 1.0 ; fDphoton[1] = 0. ; fDphoton[2] = 0. ;//constant | |
262 | // fDphoton[3] = 1.55 ; fDphoton[4] =-0.0863 ; fDphoton[5] = 0.287 ;//mean | |
263 | // fDphoton[6] = 0.0451 ; fDphoton[7] =-0.0803 ; fDphoton[8] = 0.314 ;//sigma | |
264 | ||
265 | fDphoton[0] = 4.62e-2; fDphoton[1] = 1.39e-2 ; fDphoton[2] = -3.80e-2;//constant | |
266 | //fDphoton[0] = 1.0 ; fDphoton[1] = 0. ; fDphoton[2] = 0. ;//constant | |
267 | fDphoton[3] = 1.53 ; fDphoton[4] =-6.62e-2 ; fDphoton[5] = 0.339 ;//mean | |
268 | fDphoton[6] = 6.89e-2; fDphoton[7] =-6.59e-2 ; fDphoton[8] = 0.194 ;//sigma | |
35adb638 | 269 | |
270 | fDpi0[0] = 0.0586 ; fDpi0[1] = 1.06E-3 ; fDpi0[2] = 0. ;//constant | |
271 | //fDpi0[0] = 1.0 ; fDpi0[1] = 0.0 ; fDpi0[2] = 0. ;//constant | |
272 | fDpi0[3] = 2.67 ; fDpi0[4] =-2.00E-2 ; fDpi0[5] = 9.37E-5 ;//mean | |
273 | fDpi0[6] = 0.153 ; fDpi0[7] = 9.34E-4 ; fDpi0[8] =-1.49E-5 ;//sigma | |
274 | //landau | |
275 | // fDhadron[0] = 0.007 ; fDhadron[1] = 0. ; fDhadron[2] = 0. ;//constant | |
276 | // //fDhadron[0] = 5.53 ; fDhadron[1] = 0. ; fDhadron[2] = 0. ;//constant | |
277 | // fDhadron[3] = 3.38 ; fDhadron[4] = 0.0833 ; fDhadron[5] =-0.845 ;//mean | |
278 | // fDhadron[6] = 0.627 ; fDhadron[7] = 0.012 ; fDhadron[8] =-0.170 ;//sigma | |
279 | ||
c947e71a | 280 | fDhadron[0] = 1.61E-2 ; fDhadron[1] = 3.03E-3 ; fDhadron[2] = 1.01E-2 ;//constant |
281 | fDhadron[3] = 3.81 ; fDhadron[4] = 0.232 ; fDhadron[5] =-1.25 ;//mean | |
282 | fDhadron[6] = 0.897 ; fDhadron[7] = 0.0987 ; fDhadron[8] =-0.534 ;//sigma | |
35adb638 | 283 | // Muons |
284 | fDmuon[0] = 0.0631 ; | |
285 | fDmuon[1] = 1.4 ; | |
286 | fDmuon[2] = 0.0557 ; | |
287 | fDFmuon = new TFormula("Shower shape response to muons" , "landau") ; | |
288 | fDFmuon->SetParameters( fDmuon[0], fDmuon[1], fDmuon[2]) ; | |
289 | ||
35adb638 | 290 | |
c947e71a | 291 | // x(CPV-EMC) distance gaussian parameters |
292 | ||
293 | fXelectron[0] = 8.06e-2 ; fXelectron[1] = 1.00e-2; fXelectron[2] =-5.14e-2;//constant | |
294 | //fXelectron[0] = 1.0 ; fXelectron[1] = 0. ; fXelectron[2] = 0. ;//constant | |
295 | fXelectron[3] = 0.202 ; fXelectron[4] = 8.15e-3; fXelectron[5] = 4.55 ;//mean | |
296 | fXelectron[6] = 0.334 ; fXelectron[7] = 0.186 ; fXelectron[8] = 4.32e-2;//sigma | |
297 | ||
298 | //charged hadrons gaus | |
299 | fXcharged[0] = 6.43e-3 ; fXcharged[1] =-4.19e-5; fXcharged[2] = 1.42e-3;//constant | |
300 | fXcharged[3] = 2.75 ; fXcharged[4] =-0.40 ; fXcharged[5] = 1.68 ;//mean | |
301 | fXcharged[6] = 3.135 ; fXcharged[7] =-9.41e-2; fXcharged[8] = 1.31e-2;//sigma | |
302 | ||
303 | // z(CPV-EMC) distance gaussian parameters | |
304 | ||
305 | fZelectron[0] = 8.22e-2 ; fZelectron[1] = 5.11e-3; fZelectron[2] =-3.05e-2;//constant | |
306 | //fZelectron[0] = 1.0 ; fZelectron[1] = 0. ; fZelectron[2] = 0. ;//constant | |
307 | fZelectron[3] = 3.09e-2 ; fZelectron[4] = 5.87e-2; fZelectron[5] =-9.49e-2;//mean | |
308 | fZelectron[6] = 0.263 ; fZelectron[7] =-9.02e-3; fZelectron[8] = 0.151 ;//sigma | |
309 | ||
310 | //charged hadrons gaus | |
311 | ||
312 | fZcharged[0] = 1.00e-2 ; fZcharged[1] = 2.82E-4 ; fZcharged[2] = 2.87E-3 ;//constant | |
313 | fZcharged[3] =-4.68e-2 ; fZcharged[4] =-9.21e-3 ; fZcharged[5] = 4.91e-2 ;//mean | |
314 | fZcharged[6] = 1.425 ; fZcharged[7] =-5.90e-2 ; fZcharged[8] = 5.07e-2 ;//sigma | |
315 | ||
35adb638 | 316 | for (Int_t i =0; i< AliESDtrack::kSPECIESN ; i++) |
317 | fInitPID[i] = 1.; | |
318 | ||
8d0f3f77 | 319 | } |
320 | ||
88cb7938 | 321 | //________________________________________________________________________ |
eabde521 | 322 | void AliPHOSPIDv1::Exec(Option_t *option) |
88cb7938 | 323 | { |
eabde521 | 324 | // Steering method to perform particle reconstruction and identification |
325 | // for the event range from fFirstEvent to fLastEvent. | |
326 | // This range is optionally set by SetEventRange(). | |
327 | // if fLastEvent=-1 (by default), then process events until the end. | |
61d3d6aa | 328 | |
88cb7938 | 329 | if(strstr(option,"tim")) |
330 | gBenchmark->Start("PHOSPID"); | |
331 | ||
332 | if(strstr(option,"print")) { | |
333 | Print() ; | |
334 | return ; | |
335 | } | |
336 | ||
337 | ||
adcca1e6 | 338 | AliPHOSGetter * gime = AliPHOSGetter::Instance() ; |
88cb7938 | 339 | |
71cee46d | 340 | if (fLastEvent == -1) |
341 | fLastEvent = gime->MaxEvent() - 1 ; | |
342 | else | |
343 | fLastEvent = TMath::Min(fLastEvent,gime->MaxEvent()); | |
eabde521 | 344 | Int_t nEvents = fLastEvent - fFirstEvent + 1; |
88cb7938 | 345 | |
71cee46d | 346 | Int_t ievent ; |
347 | for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) { | |
88cb7938 | 348 | gime->Event(ievent,"TR") ; |
349 | if(gime->TrackSegments() && //Skip events, where no track segments made | |
350 | gime->TrackSegments()->GetEntriesFast()) { | |
351 | MakeRecParticles() ; | |
35adb638 | 352 | |
353 | if(fBayesian) | |
354 | MakePID() ; | |
355 | ||
90cceaf6 | 356 | WriteRecParticles(); |
88cb7938 | 357 | if(strstr(option,"deb")) |
358 | PrintRecParticles(option) ; | |
359 | //increment the total number of rec particles per run | |
360 | fRecParticlesInRun += gime->RecParticles()->GetEntriesFast() ; | |
361 | } | |
362 | } | |
ff417097 | 363 | if(strstr(option,"deb")) |
364 | PrintRecParticles(option); | |
88cb7938 | 365 | if(strstr(option,"tim")){ |
366 | gBenchmark->Stop("PHOSPID"); | |
351dd634 | 367 | AliInfo(Form("took %f seconds for PID %f seconds per event", |
88cb7938 | 368 | gBenchmark->GetCpuTime("PHOSPID"), |
351dd634 | 369 | gBenchmark->GetCpuTime("PHOSPID")/nEvents)) ; |
88cb7938 | 370 | } |
adcca1e6 | 371 | if(fWrite) |
372 | Unload(); | |
88cb7938 | 373 | } |
374 | ||
35adb638 | 375 | //________________________________________________________________________ |
c947e71a | 376 | const Double_t AliPHOSPIDv1::GausF(Double_t x, Double_t y, Double_t * par) |
35adb638 | 377 | { |
c947e71a | 378 | //Given the energy x and the parameter y (tof, shower dispersion or cpv-emc distance), |
379 | //this method returns a density probability of this parameter, given by a gaussian | |
380 | //function whose parameters depend with the energy with a function: a/(x*x)+b/x+b | |
381 | Double_t cnt = par[1] / (x*x) + par[2] / x + par[0] ; | |
35adb638 | 382 | Double_t mean = par[4] / (x*x) + par[5] / x + par[3] ; |
383 | Double_t sigma = par[7] / (x*x) + par[8] / x + par[6] ; | |
c947e71a | 384 | |
35adb638 | 385 | // Double_t arg = - (y-mean) * (y-mean) / (2*sigma*sigma) ; |
386 | // return cnt * TMath::Exp(arg) ; | |
c947e71a | 387 | if(TMath::Abs(sigma) > 1.e-10){ |
35adb638 | 388 | TF1 * f = new TF1("gaus","gaus",0,100); |
389 | f->SetParameters(cnt,mean,sigma); | |
35adb638 | 390 | Double_t arg = f->Eval(y) ; |
c947e71a | 391 | //cout<<"GausF : en "<<x<<" y "<<y<<" c "<< cnt << " mean "<<mean<<" sigma "<<sigma<<endl; |
35adb638 | 392 | return arg; |
393 | } | |
394 | else | |
395 | return 0.; | |
c947e71a | 396 | |
35adb638 | 397 | } |
398 | //________________________________________________________________________ | |
c947e71a | 399 | const Double_t AliPHOSPIDv1::GausPol2(Double_t x, Double_t y, Double_t * par) |
35adb638 | 400 | { |
c947e71a | 401 | //Given the energy x and the parameter y (tof, shower dispersion or cpv-emc distance), |
402 | //this method returns a density probability of this parameter, given by a gaussian | |
403 | //function whose parameters depend with the energy like second order polinomial | |
404 | ||
35adb638 | 405 | Double_t cnt = par[0] + par[1] * x + par[2] * x * x ; |
406 | Double_t mean = par[3] + par[4] * x + par[5] * x * x ; | |
407 | Double_t sigma = par[6] + par[7] * x + par[8] * x * x ; | |
408 | ||
c947e71a | 409 | if(TMath::Abs(sigma) > 1.e-10){ |
35adb638 | 410 | TF1 * f = new TF1("gaus","gaus",0,100); |
411 | f->SetParameters(cnt,mean,sigma); | |
412 | Double_t arg = f->Eval(y) ; | |
c947e71a | 413 | //cout<<"GausPol : en "<<x<<" y "<<y<<" c "<< cnt << " mean "<<mean<<" sigma "<<sigma<<endl; |
414 | ||
35adb638 | 415 | return arg; |
416 | } | |
417 | else | |
418 | return 0.; | |
c947e71a | 419 | |
420 | ||
421 | ||
35adb638 | 422 | } |
423 | ||
69183710 | 424 | //____________________________________________________________________________ |
e3817e5f | 425 | const TString AliPHOSPIDv1::GetFileNamePrincipal(TString particle) const |
148b2bba | 426 | { |
e3817e5f | 427 | //Get file name that contains the PCA for a particle ("photon or pi0") |
428 | particle.ToLower(); | |
429 | TString name; | |
351dd634 | 430 | if (particle=="photon") |
431 | name = fFileNamePrincipalPhoton ; | |
432 | else if (particle=="pi0" ) | |
433 | name = fFileNamePrincipalPi0 ; | |
434 | else | |
435 | AliError(Form("Wrong particle name: %s (choose from pi0/photon)\n", | |
436 | particle.Data())); | |
e3817e5f | 437 | return name; |
438 | } | |
bc0c084c | 439 | |
e3817e5f | 440 | //____________________________________________________________________________ |
fc7e2f43 | 441 | Float_t AliPHOSPIDv1::GetParameterCalibration(Int_t i) const |
e3817e5f | 442 | { |
443 | // Get the i-th parameter "Calibration" | |
444 | Float_t param = 0.; | |
351dd634 | 445 | if (i>2 || i<0) { |
446 | AliError(Form("Invalid parameter number: %d",i)); | |
447 | } else | |
e3817e5f | 448 | param = (*fParameters)(0,i); |
449 | return param; | |
450 | } | |
bc0c084c | 451 | |
88cb7938 | 452 | //____________________________________________________________________________ |
fc7e2f43 | 453 | Float_t AliPHOSPIDv1::GetCalibratedEnergy(Float_t e) const |
88cb7938 | 454 | { |
455 | // It calibrates Energy depending on the recpoint energy. | |
456 | // The energy of the reconstructed cluster is corrected with | |
457 | // the formula A + B* E + C* E^2, whose parameters where obtained | |
458 | // through the study of the reconstructed energy distribution of | |
459 | // monoenergetic photons. | |
460 | ||
461 | Float_t p[]={0.,0.,0.}; | |
462 | for (Int_t i=0; i<3; i++) p[i] = GetParameterCalibration(i); | |
463 | Float_t enerec = p[0] + p[1]*e + p[2]*e*e; | |
464 | return enerec ; | |
465 | ||
466 | } | |
467 | ||
e3817e5f | 468 | //____________________________________________________________________________ |
fc7e2f43 | 469 | Float_t AliPHOSPIDv1::GetParameterCpv2Emc(Int_t i, TString axis) const |
e3817e5f | 470 | { |
471 | // Get the i-th parameter "CPV-EMC distance" for the specified axis | |
472 | Float_t param = 0.; | |
351dd634 | 473 | if(i>2 || i<0) { |
474 | AliError(Form("Invalid parameter number: %d",i)); | |
475 | } else { | |
e3817e5f | 476 | axis.ToLower(); |
351dd634 | 477 | if (axis == "x") |
478 | param = (*fParameters)(1,i); | |
479 | else if (axis == "z") | |
480 | param = (*fParameters)(2,i); | |
481 | else { | |
482 | AliError(Form("Invalid axis name: %s",axis.Data())); | |
483 | } | |
e3817e5f | 484 | } |
485 | return param; | |
486 | } | |
487 | ||
488 | //____________________________________________________________________________ | |
fc7e2f43 | 489 | Float_t AliPHOSPIDv1::GetCpv2EmcDistanceCut(TString axis, Float_t e) const |
e3817e5f | 490 | { |
88cb7938 | 491 | // Get CpvtoEmcDistance Cut depending on the cluster energy, axis and |
492 | // Purity-Efficiency point | |
493 | ||
494 | axis.ToLower(); | |
495 | Float_t p[]={0.,0.,0.}; | |
496 | for (Int_t i=0; i<3; i++) p[i] = GetParameterCpv2Emc(i,axis); | |
497 | Float_t sig = p[0] + TMath::Exp(p[1] - p[2]*e); | |
498 | return sig; | |
e3817e5f | 499 | } |
500 | ||
88cb7938 | 501 | //____________________________________________________________________________ |
fc7e2f43 | 502 | Float_t AliPHOSPIDv1::GetEllipseParameter(TString particle, TString param, Float_t e) const |
88cb7938 | 503 | { |
504 | // Calculates the parameter param of the ellipse | |
e3817e5f | 505 | |
506 | particle.ToLower(); | |
507 | param. ToLower(); | |
88cb7938 | 508 | Float_t p[4]={0.,0.,0.,0.}; |
509 | Float_t value = 0.0; | |
510 | for (Int_t i=0; i<4; i++) p[i] = GetParameterToCalculateEllipse(particle,param,i); | |
511 | if (particle == "photon") { | |
512 | if (param.Contains("a")) e = TMath::Min((Double_t)e,70.); | |
513 | else if (param.Contains("b")) e = TMath::Min((Double_t)e,70.); | |
514 | else if (param.Contains("x0")) e = TMath::Max((Double_t)e,1.1); | |
515 | } | |
e3817e5f | 516 | |
443caba9 | 517 | if (particle == "photon") |
518 | value = p[0]/TMath::Sqrt(e) + p[1]*e + p[2]*e*e + p[3]; | |
519 | else if (particle == "pi0") | |
520 | value = p[0] + p[1]*e + p[2]*e*e; | |
521 | ||
88cb7938 | 522 | return value; |
e3817e5f | 523 | } |
524 | ||
525 | //_____________________________________________________________________________ | |
fc7e2f43 | 526 | Float_t AliPHOSPIDv1::GetParameterPhotonBoundary (Int_t i) const |
e3817e5f | 527 | { |
528 | // Get the parameter "i" to calculate the boundary on the moment M2x | |
529 | // for photons at high p_T | |
530 | Float_t param = 0; | |
351dd634 | 531 | if (i>3 || i<0) { |
532 | AliError(Form("Wrong parameter number: %d\n",i)); | |
533 | } else | |
e3817e5f | 534 | param = (*fParameters)(14,i) ; |
535 | return param; | |
148b2bba | 536 | } |
e3817e5f | 537 | |
148b2bba | 538 | //____________________________________________________________________________ |
fc7e2f43 | 539 | Float_t AliPHOSPIDv1::GetParameterPi0Boundary (Int_t i) const |
e3817e5f | 540 | { |
541 | // Get the parameter "i" to calculate the boundary on the moment M2x | |
542 | // for pi0 at high p_T | |
543 | Float_t param = 0; | |
351dd634 | 544 | if (i>2 || i<0) { |
545 | AliError(Form("Wrong parameter number: %d\n",i)); | |
546 | } else | |
e3817e5f | 547 | param = (*fParameters)(15,i) ; |
548 | return param; | |
549 | } | |
148b2bba | 550 | |
e3817e5f | 551 | //____________________________________________________________________________ |
fc7e2f43 | 552 | Float_t AliPHOSPIDv1::GetParameterTimeGate(Int_t i) const |
e3817e5f | 553 | { |
88cb7938 | 554 | // Get TimeGate parameter depending on Purity-Efficiency i: |
555 | // i=0 - Low purity, i=1 - Medium purity, i=2 - High purity | |
556 | Float_t param = 0.; | |
351dd634 | 557 | if(i>2 || i<0) { |
558 | AliError(Form("Invalid Efficiency-Purity choice %d",i)); | |
559 | } else | |
88cb7938 | 560 | param = (*fParameters)(3,i) ; |
561 | return param; | |
e3817e5f | 562 | } |
563 | ||
148b2bba | 564 | //_____________________________________________________________________________ |
fc7e2f43 | 565 | Float_t AliPHOSPIDv1::GetParameterToCalculateEllipse(TString particle, TString param, Int_t i) const |
88cb7938 | 566 | { |
567 | // Get the parameter "i" that is needed to calculate the ellipse | |
568 | // parameter "param" for the particle "particle" ("photon" or "pi0") | |
569 | ||
e3817e5f | 570 | particle.ToLower(); |
571 | param. ToLower(); | |
88cb7938 | 572 | Int_t offset = -1; |
351dd634 | 573 | if (particle == "photon") |
574 | offset=0; | |
575 | else if (particle == "pi0") | |
576 | offset=5; | |
e3817e5f | 577 | else |
351dd634 | 578 | AliError(Form("Wrong particle name: %s (choose from pi0/photon)\n", |
579 | particle.Data())); | |
88cb7938 | 580 | |
581 | Int_t p= -1; | |
582 | Float_t par = 0; | |
e3817e5f | 583 | |
584 | if (param.Contains("a")) p=4+offset; | |
585 | else if(param.Contains("b")) p=5+offset; | |
586 | else if(param.Contains("c")) p=6+offset; | |
587 | else if(param.Contains("x0"))p=7+offset; | |
588 | else if(param.Contains("y0"))p=8+offset; | |
12022e83 | 589 | |
351dd634 | 590 | if (i>4 || i<0) { |
591 | AliError(Form("No parameter with index %d", i)) ; | |
592 | } else if (p==-1) { | |
593 | AliError(Form("No parameter with name %s", param.Data() )) ; | |
594 | } else | |
88cb7938 | 595 | par = (*fParameters)(p,i) ; |
596 | ||
597 | return par; | |
12022e83 | 598 | } |
599 | ||
12022e83 | 600 | |
601 | //____________________________________________________________________________ | |
8d4608b5 | 602 | Float_t AliPHOSPIDv1::GetDistance(AliPHOSEmcRecPoint * emc,AliPHOSCpvRecPoint * cpv, Option_t * axis)const |
69183710 | 603 | { |
604 | // Calculates the distance between the EMC RecPoint and the PPSD RecPoint | |
148b2bba | 605 | |
88cb7938 | 606 | const AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ; |
69183710 | 607 | TVector3 vecEmc ; |
7acf6008 | 608 | TVector3 vecCpv ; |
148b2bba | 609 | if(cpv){ |
610 | emc->GetLocalPosition(vecEmc) ; | |
611 | cpv->GetLocalPosition(vecCpv) ; | |
2bb500e5 | 612 | |
148b2bba | 613 | if(emc->GetPHOSMod() == cpv->GetPHOSMod()){ |
614 | // Correct to difference in CPV and EMC position due to different distance to center. | |
615 | // we assume, that particle moves from center | |
616 | Float_t dCPV = geom->GetIPtoOuterCoverDistance(); | |
617 | Float_t dEMC = geom->GetIPtoCrystalSurface() ; | |
618 | dEMC = dEMC / dCPV ; | |
619 | vecCpv = dEMC * vecCpv - vecEmc ; | |
e3817e5f | 620 | if (axis == "X") return vecCpv.X(); |
621 | if (axis == "Y") return vecCpv.Y(); | |
622 | if (axis == "Z") return vecCpv.Z(); | |
623 | if (axis == "R") return vecCpv.Mag(); | |
624 | } | |
148b2bba | 625 | return 100000000 ; |
626 | } | |
7acf6008 | 627 | return 100000000 ; |
69183710 | 628 | } |
bc0c084c | 629 | //____________________________________________________________________________ |
8d4608b5 | 630 | Int_t AliPHOSPIDv1::GetCPVBit(AliPHOSEmcRecPoint * emc,AliPHOSCpvRecPoint * cpv, Int_t effPur, Float_t e) const |
bc0c084c | 631 | { |
c947e71a | 632 | //Calculates the pid bit for the CPV selection per each purity. |
e3817e5f | 633 | if(effPur>2 || effPur<0) |
351dd634 | 634 | AliError(Form("Invalid Efficiency-Purity choice %d",effPur)); |
bc0c084c | 635 | |
e3817e5f | 636 | Float_t sigX = GetCpv2EmcDistanceCut("X",e); |
637 | Float_t sigZ = GetCpv2EmcDistanceCut("Z",e); | |
bc0c084c | 638 | |
639 | Float_t deltaX = TMath::Abs(GetDistance(emc, cpv, "X")); | |
640 | Float_t deltaZ = TMath::Abs(GetDistance(emc, cpv, "Z")); | |
35adb638 | 641 | //Info("GetCPVBit"," xdist %f, sigx %f, zdist %f, sigz %f",deltaX, sigX, deltaZ,sigZ ) ; |
642 | if((deltaX>sigX*(effPur+1))&&(deltaZ>sigZ*(effPur+1))) | |
bc0c084c | 643 | return 1;//Neutral |
644 | else | |
645 | return 0;//Charged | |
bc0c084c | 646 | } |
69183710 | 647 | |
6ad0bfa0 | 648 | //____________________________________________________________________________ |
fc7e2f43 | 649 | Int_t AliPHOSPIDv1::GetPrincipalBit(TString particle, const Double_t* p, Int_t effPur, Float_t e)const |
148b2bba | 650 | { |
50739f15 | 651 | //Is the particle inside de PCA ellipse? |
581354c5 | 652 | |
e3817e5f | 653 | particle.ToLower(); |
654 | Int_t prinbit = 0 ; | |
655 | Float_t a = GetEllipseParameter(particle,"a" , e); | |
656 | Float_t b = GetEllipseParameter(particle,"b" , e); | |
657 | Float_t c = GetEllipseParameter(particle,"c" , e); | |
658 | Float_t x0 = GetEllipseParameter(particle,"x0", e); | |
659 | Float_t y0 = GetEllipseParameter(particle,"y0", e); | |
660 | ||
661 | Float_t r = TMath::Power((p[0] - x0)/a,2) + | |
662 | TMath::Power((p[1] - y0)/b,2) + | |
663 | c*(p[0] - x0)*(p[1] - y0)/(a*b) ; | |
50739f15 | 664 | //3 different ellipses defined |
e3817e5f | 665 | if((effPur==2) && (r<1./2.)) prinbit= 1; |
666 | if((effPur==1) && (r<2. )) prinbit= 1; | |
667 | if((effPur==0) && (r<9./2.)) prinbit= 1; | |
50739f15 | 668 | |
581354c5 | 669 | if(r<0) |
351dd634 | 670 | AliError("Negative square?") ; |
1f0e7ccd | 671 | |
672 | return prinbit; | |
148b2bba | 673 | |
148b2bba | 674 | } |
1f0e7ccd | 675 | //____________________________________________________________________________ |
fc7e2f43 | 676 | Int_t AliPHOSPIDv1::GetHardPhotonBit(AliPHOSEmcRecPoint * emc) const |
1f0e7ccd | 677 | { |
e3817e5f | 678 | // Set bit for identified hard photons (E > 30 GeV) |
679 | // if the second moment M2x is below the boundary | |
680 | ||
681 | Float_t e = emc->GetEnergy(); | |
682 | if (e < 30.0) return 0; | |
683 | Float_t m2x = emc->GetM2x(); | |
684 | Float_t m2xBoundary = GetParameterPhotonBoundary(0) * | |
685 | TMath::Exp(-TMath::Power(e-GetParameterPhotonBoundary(1),2)/2.0/ | |
686 | TMath::Power(GetParameterPhotonBoundary(2),2)) + | |
687 | GetParameterPhotonBoundary(3); | |
351dd634 | 688 | AliDebug(1, Form("GetHardPhotonBit","E=%f, m2x=%f, boundary=%f", |
689 | e,m2x,m2xBoundary)); | |
e3817e5f | 690 | if (m2x < m2xBoundary) |
691 | return 1;// A hard photon | |
692 | else | |
693 | return 0;// Not a hard photon | |
1f0e7ccd | 694 | } |
92f521a9 | 695 | |
e3817e5f | 696 | //____________________________________________________________________________ |
fc7e2f43 | 697 | Int_t AliPHOSPIDv1::GetHardPi0Bit(AliPHOSEmcRecPoint * emc) const |
e3817e5f | 698 | { |
699 | // Set bit for identified hard pi0 (E > 30 GeV) | |
700 | // if the second moment M2x is above the boundary | |
701 | ||
702 | Float_t e = emc->GetEnergy(); | |
703 | if (e < 30.0) return 0; | |
704 | Float_t m2x = emc->GetM2x(); | |
705 | Float_t m2xBoundary = GetParameterPi0Boundary(0) + | |
706 | e * GetParameterPi0Boundary(1); | |
351dd634 | 707 | AliDebug(1,Form("E=%f, m2x=%f, boundary=%f",e,m2x,m2xBoundary)); |
e3817e5f | 708 | if (m2x > m2xBoundary) |
709 | return 1;// A hard pi0 | |
bc0c084c | 710 | else |
e3817e5f | 711 | return 0;// Not a hard pi0 |
f0a4c9e9 | 712 | } |
e3817e5f | 713 | |
714 | //____________________________________________________________________________ | |
8d4608b5 | 715 | TVector3 AliPHOSPIDv1::GetMomentumDirection(AliPHOSEmcRecPoint * emc, AliPHOSCpvRecPoint * )const |
88cb7938 | 716 | { |
717 | // Calculates the momentum direction: | |
718 | // 1. if only a EMC RecPoint, direction is given by IP and this RecPoint | |
719 | // 2. if a EMC RecPoint and CPV RecPoint, direction is given by the line through the 2 recpoints | |
720 | // However because of the poor position resolution of PPSD the direction is always taken as if we were | |
721 | // in case 1. | |
f0a4c9e9 | 722 | |
88cb7938 | 723 | TVector3 dir(0,0,0) ; |
9688c1dd | 724 | |
88cb7938 | 725 | TVector3 emcglobalpos ; |
726 | TMatrix dummy ; | |
bf8f1fbd | 727 | |
88cb7938 | 728 | emc->GetGlobalPosition(emcglobalpos, dummy) ; |
bf8f1fbd | 729 | |
fbf811ec | 730 | |
88cb7938 | 731 | dir = emcglobalpos ; |
88cb7938 | 732 | dir.SetMag(1.) ; |
e3817e5f | 733 | |
88cb7938 | 734 | //account correction to the position of IP |
735 | Float_t xo,yo,zo ; //Coordinates of the origin | |
adcca1e6 | 736 | if(gAlice && gAlice->GetMCApp() && gAlice->Generator()) |
737 | gAlice->Generator()->GetOrigin(xo,yo,zo) ; | |
738 | else{ | |
739 | xo=yo=zo=0.; | |
740 | } | |
88cb7938 | 741 | TVector3 origin(xo,yo,zo); |
742 | dir = dir - origin ; | |
e3817e5f | 743 | |
88cb7938 | 744 | return dir ; |
7b7c1533 | 745 | } |
746 | ||
35adb638 | 747 | //________________________________________________________________________ |
c947e71a | 748 | const Double_t AliPHOSPIDv1::LandauF(Double_t x, Double_t y, Double_t * par) |
35adb638 | 749 | { |
c947e71a | 750 | //Given the energy x and the parameter y (tof, shower dispersion or cpv-emc distance), |
751 | //this method returns a density probability of this parameter, given by a landau | |
752 | //function whose parameters depend with the energy with a function: a/(x*x)+b/x+b | |
753 | ||
754 | Double_t cnt = par[1] / (x*x) + par[2] / x + par[0] ; | |
35adb638 | 755 | Double_t mean = par[4] / (x*x) + par[5] / x + par[3] ; |
756 | Double_t sigma = par[7] / (x*x) + par[8] / x + par[6] ; | |
c947e71a | 757 | |
758 | if(TMath::Abs(sigma) > 1.e-10){ | |
35adb638 | 759 | TF1 * f = new TF1("landau","landau",0.,100.); |
760 | f->SetParameters(cnt,mean,sigma); | |
761 | Double_t arg = f->Eval(y) ; | |
c947e71a | 762 | // cout<<"LandauF : en "<<x<<" y "<<y<<" c "<< cnt << " mean "<<mean<<" sigma "<<sigma<<endl; |
763 | ||
35adb638 | 764 | return arg; |
765 | } | |
766 | else | |
767 | return 0.; | |
768 | ||
769 | } | |
770 | //________________________________________________________________________ | |
c947e71a | 771 | const Double_t AliPHOSPIDv1::LandauPol2(Double_t x, Double_t y, Double_t * par) |
35adb638 | 772 | { |
c947e71a | 773 | |
774 | //Given the energy x and the parameter y (tof, shower dispersion or cpv-emc distance), | |
775 | //this method returns a density probability of this parameter, given by a landau | |
776 | //function whose parameters depend with the energy like second order polinomial | |
777 | ||
35adb638 | 778 | Double_t cnt = par[2] * (x*x) + par[1] * x + par[0] ; |
c947e71a | 779 | Double_t mean = par[5] * (x*x) + par[4] * x + par[3] ; |
780 | Double_t sigma = par[8] * (x*x) + par[7] * x + par[6] ; | |
35adb638 | 781 | |
c947e71a | 782 | if(TMath::Abs(sigma) > 1.e-10){ |
35adb638 | 783 | TF1 * f = new TF1("landau","landau",0.,100.); |
784 | f->SetParameters(cnt,mean,sigma); | |
785 | Double_t arg = f->Eval(y) ; | |
c947e71a | 786 | //cout<<"LandauPol : en "<<x<<" y "<<y<<" c "<< cnt << " mean "<<mean<<" sigma "<<sigma<<endl; |
35adb638 | 787 | return arg; |
788 | } | |
789 | else | |
790 | return 0.; | |
c947e71a | 791 | |
792 | ||
35adb638 | 793 | } |
794 | // //________________________________________________________________________ | |
795 | // Double_t AliPHOSPIDv1::ChargedHadronDistProb(Double_t x, Double_t y, Double_t * parg, Double_t * parl) | |
796 | // { | |
797 | // Double_t cnt = 0.0 ; | |
798 | // Double_t mean = 0.0 ; | |
799 | // Double_t sigma = 0.0 ; | |
800 | // Double_t arg = 0.0 ; | |
801 | // if (y < parl[4] / (x*x) + parl[5] / x + parl[3]){ | |
802 | // cnt = parg[1] / (x*x) + parg[2] / x + parg[0] ; | |
803 | // mean = parg[4] / (x*x) + parg[5] / x + parg[3] ; | |
804 | // sigma = parg[7] / (x*x) + parg[8] / x + parg[6] ; | |
805 | // TF1 * f = new TF1("gaus","gaus",0.,100.); | |
806 | // f->SetParameters(cnt,mean,sigma); | |
807 | // arg = f->Eval(y) ; | |
808 | // } | |
809 | // else{ | |
810 | // cnt = parl[1] / (x*x) + parl[2] / x + parl[0] ; | |
811 | // mean = parl[4] / (x*x) + parl[5] / x + parl[3] ; | |
812 | // sigma = parl[7] / (x*x) + parl[8] / x + parl[6] ; | |
813 | // TF1 * f = new TF1("landau","landau",0.,100.); | |
814 | // f->SetParameters(cnt,mean,sigma); | |
815 | // arg = f->Eval(y) ; | |
816 | // } | |
817 | // // Double_t mean = par[3] + par[4] * x + par[5] * x * x ; | |
818 | // // Double_t sigma = par[6] + par[7] * x + par[8] * x * x ; | |
819 | ||
820 | // //Double_t arg = -(y-mean)*(y-mean)/(2*sigma*sigma) ; | |
821 | // //return cnt * TMath::Exp(arg) ; | |
822 | ||
823 | // return arg; | |
824 | ||
825 | // } | |
2cc71c1e | 826 | //____________________________________________________________________________ |
827 | void AliPHOSPIDv1::MakePID() | |
828 | { | |
829 | // construct the PID weight from a Bayesian Method | |
c947e71a | 830 | |
2cc71c1e | 831 | Int_t index ; |
cc1fe362 | 832 | const Int_t kSPECIES = AliESDtrack::kSPECIESN ; |
2cc71c1e | 833 | Int_t nparticles = AliPHOSGetter::Instance()->RecParticles()->GetEntriesFast() ; |
c947e71a | 834 | |
835 | // const Int_t kMAXPARTICLES = 2000 ; | |
836 | // if (nparticles >= kMAXPARTICLES) | |
837 | // Error("MakePID", "Change size of MAXPARTICLES") ; | |
838 | // Double_t stof[kSPECIES][kMAXPARTICLES] ; | |
839 | ||
53ab54c8 | 840 | // const Int_t kMAXPARTICLES = 2000 ; |
841 | // if (nparticles >= kMAXPARTICLES) | |
351dd634 | 842 | // AliError("Change size of MAXPARTICLES") ; |
53ab54c8 | 843 | // Double_t stof[kSPECIES][kMAXPARTICLES] ; |
53ab54c8 | 844 | |
35adb638 | 845 | Double_t * stof[kSPECIES] ; |
846 | Double_t * sdp [kSPECIES] ; | |
847 | Double_t * scpv[kSPECIES] ; | |
848 | ||
849 | //Info("MakePID","Begin MakePID"); | |
850 | ||
851 | for (Int_t i =0; i< kSPECIES; i++){ | |
852 | stof[i] = new Double_t[nparticles] ; | |
853 | sdp [i] = new Double_t[nparticles] ; | |
854 | scpv[i] = new Double_t[nparticles] ; | |
855 | } | |
856 | ||
cc1fe362 | 857 | // make the normalized distribution of pid for this event |
858 | // w(pid) in the Bayesian formulation | |
2cc71c1e | 859 | for(index = 0 ; index < nparticles ; index ++) { |
cc1fe362 | 860 | |
35adb638 | 861 | AliPHOSRecParticle * recpar = AliPHOSGetter::Instance()->RecParticle(index) ; |
862 | AliPHOSEmcRecPoint * emc = AliPHOSGetter::Instance()->EmcRecPoint(index) ; | |
863 | AliPHOSCpvRecPoint * cpv = AliPHOSGetter::Instance()->CpvRecPoint(index) ; | |
cc1fe362 | 864 | |
35adb638 | 865 | Float_t en = emc->GetEnergy(); |
c947e71a | 866 | |
35adb638 | 867 | // Tof |
c947e71a | 868 | |
869 | // Info("MakePID", "TOF"); | |
35adb638 | 870 | Double_t time = recpar->ToF() ; |
871 | //cout<<">>>>>>>Energy "<<en<<"Time "<<time<<endl; | |
872 | //Electrons initial population to be removed | |
873 | fInitPID[AliESDtrack::kEleCon] = 0. ; | |
cc1fe362 | 874 | |
35adb638 | 875 | // now get the signals probability |
876 | // s(pid) in the Bayesian formulation | |
cc1fe362 | 877 | |
c947e71a | 878 | stof[AliESDtrack::kPhoton][index] = 1.; |
879 | stof[AliESDtrack::kElectron][index] = 1.; | |
880 | stof[AliESDtrack::kPion][index] = 1.; | |
881 | stof[AliESDtrack::kKaon][index] = 1.; | |
882 | stof[AliESDtrack::kProton][index] = 1.; | |
883 | stof[AliESDtrack::kNeutron][index] = 1.; | |
884 | stof[AliESDtrack::kEleCon][index] = 1.; | |
885 | stof[AliESDtrack::kKaon0][index] = 1.; | |
886 | stof[AliESDtrack::kMuon][index] = 1.; | |
887 | ||
35adb638 | 888 | if(en < 2.) { |
c947e71a | 889 | stof[AliESDtrack::kPhoton][index] = fTFphoton ->Eval(time) ; //gaus distribution |
890 | stof[AliESDtrack::kPion][index] = fTFpiong ->Eval(time) ; //gaus distribution | |
891 | stof[AliESDtrack::kElectron][index] = stof[AliESDtrack::kPion][index] ; | |
892 | ||
35adb638 | 893 | if(time < fTkaonl[1]) |
894 | stof[AliESDtrack::kKaon][index] = fTFkaong ->Eval(time) ; //gaus distribution | |
895 | else | |
896 | stof[AliESDtrack::kKaon][index] = fTFkaonl ->Eval(time) ; //landau distribution | |
897 | if(time < fThhadronl[1]) | |
898 | stof[AliESDtrack::kProton][index] = fTFhhadrong ->Eval(time) ; //gaus distribution | |
899 | else | |
900 | stof[AliESDtrack::kProton][index] = fTFhhadronl ->Eval(time) ; //landau distribution | |
c947e71a | 901 | |
35adb638 | 902 | stof[AliESDtrack::kNeutron][index] = stof[AliESDtrack::kProton][index] ; |
903 | stof[AliESDtrack::kEleCon][index] = stof[AliESDtrack::kPhoton][index] ; | |
904 | // a conversion electron has the photon ToF | |
905 | stof[AliESDtrack::kKaon0][index] = stof[AliESDtrack::kKaon][index] ; | |
906 | stof[AliESDtrack::kMuon][index] = stof[AliESDtrack::kPhoton][index] ; | |
cc1fe362 | 907 | } |
c947e71a | 908 | |
909 | // Info("MakePID", "Dispersion"); | |
cc1fe362 | 910 | |
35adb638 | 911 | // Shower shape: Dispersion |
912 | Float_t dispersion = emc->GetDispersion(); | |
913 | //dispersion is not well defined if the cluster is only in few crystals | |
cc1fe362 | 914 | |
c947e71a | 915 | sdp[AliESDtrack::kPhoton][index] = 1. ; |
916 | sdp[AliESDtrack::kElectron][index] = 1. ; | |
917 | sdp[AliESDtrack::kPion][index] = 1. ; | |
918 | sdp[AliESDtrack::kKaon][index] = 1. ; | |
919 | sdp[AliESDtrack::kProton][index] = 1. ; | |
920 | sdp[AliESDtrack::kNeutron][index] = 1. ; | |
921 | sdp[AliESDtrack::kEleCon][index] = 1. ; | |
922 | sdp[AliESDtrack::kKaon0][index] = 1. ; | |
923 | sdp[AliESDtrack::kMuon][index] = 1. ; | |
924 | ||
925 | if(en > 0.5 && emc->GetMultiplicity() > 3){ | |
35adb638 | 926 | sdp[AliESDtrack::kPhoton][index] = GausF (en , dispersion, fDphoton) ; |
927 | sdp[AliESDtrack::kElectron][index] = sdp[AliESDtrack::kPhoton][index] ; | |
928 | sdp[AliESDtrack::kPion][index] = LandauF(en , dispersion, fDhadron ) ; | |
929 | sdp[AliESDtrack::kKaon][index] = sdp[AliESDtrack::kPion][index] ; | |
930 | sdp[AliESDtrack::kProton][index] = sdp[AliESDtrack::kPion][index] ; | |
931 | sdp[AliESDtrack::kNeutron][index] = sdp[AliESDtrack::kPion][index] ; | |
932 | sdp[AliESDtrack::kEleCon][index] = sdp[AliESDtrack::kPhoton][index]; | |
933 | sdp[AliESDtrack::kKaon0][index] = sdp[AliESDtrack::kPion][index] ; | |
934 | sdp[AliESDtrack::kMuon][index] = fDFmuon ->Eval(dispersion) ; //landau distribution | |
935 | } | |
cc1fe362 | 936 | |
c947e71a | 937 | // Info("MakePID","multiplicity %d, dispersion %f", emc->GetMultiplicity(), dispersion); |
938 | // Info("MakePID","ss: photon %f, hadron %f ", sdp[AliESDtrack::kPhoton][index], sdp[AliESDtrack::kPion][index]); | |
adcca1e6 | 939 | |
c947e71a | 940 | // cout<<">>>>>multiplicity "<<emc->GetMultiplicity()<<", dispersion "<< dispersion<<endl ; |
941 | // cout<<"<<<<<ss: photon "<<sdp[AliESDtrack::kPhoton][index]<<", hadron "<<sdp[AliESDtrack::kPion][index]<<endl; | |
942 | ||
35adb638 | 943 | // CPV-EMC Distance |
c947e71a | 944 | // Info("MakePID", "Distance"); |
945 | // Float_t distance = GetDistance(emc, cpv, "R") ; | |
946 | Float_t x = TMath::Abs(GetDistance(emc, cpv, "X")) ; | |
947 | Float_t z = GetDistance(emc, cpv, "Z") ; | |
35adb638 | 948 | // Info("MakePID", "Distance %f", distance); |
c947e71a | 949 | Double_t pcpv = 0 ; |
950 | Double_t pcpvneutral = 0. ; | |
951 | Double_t elprobx = GausF(en , x, fXelectron) ; | |
952 | Double_t elprobz = GausF(en , z, fZelectron) ; | |
953 | Double_t chprobx = GausF(en , x, fXcharged) ; | |
954 | Double_t chprobz = GausF(en , z, fZcharged) ; | |
955 | Double_t pcpvelectron = elprobx * elprobz; | |
956 | Double_t pcpvcharged = chprobx * chprobz; | |
957 | ||
958 | // cout<<">>>>electron : x "<<x<<" xprob "<<elprobx<<" z "<<z<<" zprob "<<elprobz<<endl; | |
959 | // cout<<">>>>hadron : x "<<x<<" xprob "<<chprobx<<" z "<<z<<" zprob "<<chprobz<<endl; | |
960 | // cout<<">>>>electron : px*pz "<<pcpvelectron <<" hadron: px*pz "<<pcpvcharged<<endl; | |
961 | ||
35adb638 | 962 | if(pcpvelectron >= pcpvcharged) |
963 | pcpv = pcpvelectron ; | |
964 | else | |
965 | pcpv = pcpvcharged ; | |
966 | ||
c947e71a | 967 | if(pcpv < 1e-7) |
35adb638 | 968 | { |
969 | pcpvneutral = 1. ; | |
970 | pcpvcharged = 0. ; | |
971 | pcpvelectron = 0. ; | |
972 | } | |
c947e71a | 973 | // else |
974 | // cout<<">>>>>>>>>>>CHARGED>>>>>>>>>>>"<<endl; | |
35adb638 | 975 | |
976 | scpv[AliESDtrack::kPion][index] = pcpvcharged ; | |
977 | scpv[AliESDtrack::kKaon][index] = pcpvcharged ; | |
978 | scpv[AliESDtrack::kProton][index] = pcpvcharged ; | |
979 | scpv[AliESDtrack::kPhoton][index] = pcpvneutral ; | |
980 | scpv[AliESDtrack::kElectron][index] = pcpvelectron ; | |
981 | scpv[AliESDtrack::kNeutron][index] = pcpvneutral ; | |
982 | scpv[AliESDtrack::kEleCon][index] = pcpvelectron ; | |
983 | scpv[AliESDtrack::kKaon0][index] = pcpvneutral ; | |
984 | scpv[AliESDtrack::kMuon][index] = pcpvelectron ; | |
985 | ||
986 | // Info("MakePID", "CPV passed"); | |
c947e71a | 987 | |
988 | //Pi0 | |
989 | stof[AliESDtrack::kPi0][index] = 0. ; | |
990 | scpv[AliESDtrack::kPi0][index] = 0. ; | |
991 | sdp [AliESDtrack::kPi0][index] = 0. ; | |
992 | fInitPID[AliESDtrack::kPi0] = 0. ; | |
993 | ||
35adb638 | 994 | if(en > 30.){ |
995 | // pi0 are detected via decay photon | |
996 | stof[AliESDtrack::kPi0][index] = fTFphoton ->Eval(time) ; | |
997 | scpv[AliESDtrack::kPi0][index] = pcpvneutral ; | |
c947e71a | 998 | sdp [AliESDtrack::kPi0][index] = 1. ; |
999 | if(emc->GetMultiplicity() > 3) | |
35adb638 | 1000 | sdp [AliESDtrack::kPi0][index] = GausPol2(en , dispersion, fDpi0) ; |
35adb638 | 1001 | } |
1002 | ||
1003 | if(en > 0.5){ | |
1004 | //Muons deposit few energy | |
c947e71a | 1005 | scpv[AliESDtrack::kMuon][index] = 0 ; |
1006 | stof[AliESDtrack::kMuon][index] = 0 ; | |
1007 | sdp [AliESDtrack::kMuon][index] = 0 ; | |
1008 | } | |
1009 | ||
1010 | if(en > 0.5){ | |
1011 | cout<<"######################################################"<<endl; | |
1012 | //cout<<"MakePID: energy "<<en<<", tof "<<time<<", distance "<<distance<<", dispersion "<<dispersion<<endl ; | |
1013 | cout<<"MakePID: energy "<<en<<", tof "<<time<<", dispersion "<<dispersion<<", x "<<x<<", z "<<z<<endl ; | |
1014 | cout<<">>>>>multiplicity "<<emc->GetMultiplicity()<<endl; | |
1015 | cout<<">>>>electron : xprob "<<elprobx<<" zprob "<<elprobz<<endl; | |
1016 | cout<<">>>>hadron : xprob "<<chprobx<<" zprob "<<chprobz<<endl; | |
1017 | cout<<">>>>electron : px*pz "<<pcpvelectron <<" hadron: px*pz "<<pcpvcharged<<endl; | |
1018 | ||
1019 | ||
1020 | cout<<"Photon , pid "<< fInitPID[AliESDtrack::kPhoton]<<" tof "<<stof[AliESDtrack::kPhoton][index] | |
1021 | <<", cpv "<<scpv[AliESDtrack::kPhoton][index]<<", ss "<<sdp[AliESDtrack::kPhoton][index]<<endl; | |
1022 | cout<<"EleCon , pid "<< fInitPID[AliESDtrack::kEleCon]<<", tof "<<stof[AliESDtrack::kEleCon][index] | |
1023 | <<", cpv "<<scpv[AliESDtrack::kEleCon][index]<<" ss "<<sdp[AliESDtrack::kEleCon][index]<<endl; | |
1024 | cout<<"Electron , pid "<< fInitPID[AliESDtrack::kElectron]<<", tof "<<stof[AliESDtrack::kElectron][index] | |
1025 | <<", cpv "<<scpv[AliESDtrack::kElectron][index]<<" ss "<<sdp[AliESDtrack::kElectron][index]<<endl; | |
1026 | cout<<"Muon , pid "<< fInitPID[AliESDtrack::kMuon]<<", tof "<<stof[AliESDtrack::kMuon][index] | |
1027 | <<", cpv "<<scpv[AliESDtrack::kMuon][index]<<" ss "<<sdp[AliESDtrack::kMuon][index]<<endl; | |
1028 | cout<<"Pi0 , pid "<< fInitPID[AliESDtrack::kPi0]<<", tof "<<stof[AliESDtrack::kPi0][index] | |
1029 | <<", cpv "<<scpv[AliESDtrack::kPi0][index]<<" ss "<<sdp[AliESDtrack::kPi0][index]<<endl; | |
1030 | cout<<"Pion , pid "<< fInitPID[AliESDtrack::kPion]<<", tof "<<stof[AliESDtrack::kPion][index] | |
1031 | <<", cpv "<<scpv[AliESDtrack::kPion][index]<<" ss "<<sdp[AliESDtrack::kPion][index]<<endl; | |
1032 | cout<<"Kaon0 , pid "<< fInitPID[AliESDtrack::kKaon0]<<", tof "<<stof[AliESDtrack::kKaon0][index] | |
1033 | <<", cpv "<<scpv[AliESDtrack::kKaon0][index]<<" ss "<<sdp[AliESDtrack::kKaon0][index]<<endl; | |
1034 | cout<<"Kaon , pid "<< fInitPID[AliESDtrack::kKaon]<<", tof "<<stof[AliESDtrack::kKaon][index] | |
1035 | <<", cpv "<<scpv[AliESDtrack::kKaon][index]<<" ss "<<sdp[AliESDtrack::kKaon][index]<<endl; | |
1036 | cout<<"Neutron , pid "<< fInitPID[AliESDtrack::kNeutron]<<", tof "<<stof[AliESDtrack::kNeutron][index] | |
1037 | <<", cpv "<<scpv[AliESDtrack::kNeutron][index]<<" ss "<<sdp[AliESDtrack::kNeutron][index]<<endl; | |
1038 | cout<<"Proton , pid "<< fInitPID[AliESDtrack::kProton]<<", tof "<<stof[AliESDtrack::kProton][index] | |
1039 | <<", cpv "<<scpv[AliESDtrack::kProton][index]<<" ss "<<sdp[AliESDtrack::kProton][index]<<endl; | |
1040 | cout<<"######################################################"<<endl; | |
35adb638 | 1041 | } |
cc1fe362 | 1042 | } |
35adb638 | 1043 | |
1044 | //for (index = 0 ; index < kSPECIES ; index++) | |
1045 | // pid[index] /= nparticles ; | |
1046 | ||
1047 | // Info("MakePID", "Total Probability calculation"); | |
1048 | ||
cc1fe362 | 1049 | for(index = 0 ; index < nparticles ; index ++) { |
1050 | // calculates the Bayesian weight | |
1051 | Int_t jndex ; | |
1052 | Double_t wn = 0.0 ; | |
1053 | for (jndex = 0 ; jndex < kSPECIES ; jndex++) | |
35adb638 | 1054 | //wn += stof[jndex][index] * pid[jndex] ; |
1055 | wn += stof[jndex][index] * sdp[jndex][index] * scpv[jndex][index] * fInitPID[jndex] ; | |
1056 | //cout<<"*************wn "<<wn<<endl; | |
cc1fe362 | 1057 | AliPHOSRecParticle * recpar = AliPHOSGetter::Instance()->RecParticle(index) ; |
e74ea0e9 | 1058 | if (TMath::Abs(wn)>0) |
1059 | for (jndex = 0 ; jndex < kSPECIES ; jndex++) { | |
35adb638 | 1060 | //cout<<"jndex "<<jndex<<" wn "<<wn<<" SetPID * wn" |
1061 | //<<stof[jndex][index] * sdp[jndex][index] * pid[jndex] << endl; | |
1062 | //cout<<" tof "<<stof[jndex][index] << " disp " <<sdp[jndex][index] << " pid "<< fInitPID[jndex] << endl; | |
1063 | // cout<<"Particle "<<jndex<<" final prob * wn " | |
1064 | // <<stof[jndex][index] * sdp[jndex][index] * scpv[jndex][index] * fInitPID[jndex] <<" wn "<< wn<<endl; | |
1065 | recpar->SetPID(jndex, stof[jndex][index] * sdp[jndex][index] * | |
1066 | scpv[jndex][index] * fInitPID[jndex] / wn) ; | |
1067 | // cout<<"final prob "<<stof[jndex][index] * sdp[jndex][index] * scpv[jndex][index] * fInitPID[jndex] / wn<<endl; | |
1068 | //recpar->SetPID(jndex, stof[jndex][index] * fInitPID[jndex] / wn) ; | |
1069 | //cout<<"After SetPID"<<endl; | |
1070 | //recpar->Print(); | |
e74ea0e9 | 1071 | } |
2cc71c1e | 1072 | } |
35adb638 | 1073 | // Info("MakePID", "Delete"); |
1074 | ||
1075 | // for (Int_t i =0; i< kSPECIES; i++){ | |
1076 | // delete [] stof[i]; | |
1077 | // cout<<i<<endl; | |
1078 | // delete [] sdp[i]; | |
1079 | // cout<<i<<endl; | |
1080 | // delete [] scpv[i]; | |
1081 | // cout<<i<<endl; | |
1082 | // } | |
1083 | ||
1084 | // Info("MakePID","End MakePID"); | |
2cc71c1e | 1085 | } |
1086 | ||
7acf6008 | 1087 | //____________________________________________________________________________ |
e3817e5f | 1088 | void AliPHOSPIDv1::MakeRecParticles() |
1089 | { | |
b2a60966 | 1090 | // Makes a RecParticle out of a TrackSegment |
148b2bba | 1091 | |
88cb7938 | 1092 | AliPHOSGetter * gime = AliPHOSGetter::Instance() ; |
fbf811ec | 1093 | TObjArray * emcRecPoints = gime->EmcRecPoints() ; |
1094 | TObjArray * cpvRecPoints = gime->CpvRecPoints() ; | |
1095 | TClonesArray * trackSegments = gime->TrackSegments() ; | |
148b2bba | 1096 | if ( !emcRecPoints || !cpvRecPoints || !trackSegments ) { |
351dd634 | 1097 | AliFatal("RecPoints or TrackSegments not found !") ; |
148b2bba | 1098 | } |
fbf811ec | 1099 | TClonesArray * recParticles = gime->RecParticles() ; |
01a599c9 | 1100 | recParticles->Clear(); |
148b2bba | 1101 | |
7b7c1533 | 1102 | TIter next(trackSegments) ; |
7acf6008 | 1103 | AliPHOSTrackSegment * ts ; |
6ad0bfa0 | 1104 | Int_t index = 0 ; |
09fc14a0 | 1105 | AliPHOSRecParticle * rp ; |
7acf6008 | 1106 | while ( (ts = (AliPHOSTrackSegment *)next()) ) { |
fad3e5b9 | 1107 | |
7b7c1533 | 1108 | new( (*recParticles)[index] ) AliPHOSRecParticle() ; |
1109 | rp = (AliPHOSRecParticle *)recParticles->At(index) ; | |
f0a4c9e9 | 1110 | rp->SetTrackSegment(index) ; |
9688c1dd | 1111 | rp->SetIndexInList(index) ; |
148b2bba | 1112 | |
7acf6008 | 1113 | AliPHOSEmcRecPoint * emc = 0 ; |
1114 | if(ts->GetEmcIndex()>=0) | |
7b7c1533 | 1115 | emc = (AliPHOSEmcRecPoint *) emcRecPoints->At(ts->GetEmcIndex()) ; |
fad3e5b9 | 1116 | |
8d4608b5 | 1117 | AliPHOSCpvRecPoint * cpv = 0 ; |
7acf6008 | 1118 | if(ts->GetCpvIndex()>=0) |
8d4608b5 | 1119 | cpv = (AliPHOSCpvRecPoint *) cpvRecPoints->At(ts->GetCpvIndex()) ; |
fad3e5b9 | 1120 | |
bd76890a | 1121 | Int_t track = 0 ; |
1122 | track = ts->GetTrackIndex() ; | |
1123 | ||
148b2bba | 1124 | // Now set type (reconstructed) of the particle |
1125 | ||
1126 | // Choose the cluster energy range | |
9fa5f1d0 | 1127 | |
fbf811ec | 1128 | if (!emc) { |
351dd634 | 1129 | AliFatal(Form("-> emc(%d) = %d", ts->GetEmcIndex(), emc )) ; |
fbf811ec | 1130 | } |
50739f15 | 1131 | |
e3817e5f | 1132 | Float_t e = emc->GetEnergy() ; |
bc0c084c | 1133 | |
6f969528 | 1134 | Float_t lambda[2] ; |
1135 | emc->GetElipsAxis(lambda) ; | |
50739f15 | 1136 | |
1137 | if((lambda[0]>0.01) && (lambda[1]>0.01)){ | |
1138 | // Looking PCA. Define and calculate the data (X), | |
bc0c084c | 1139 | // introduce in the function X2P that gives the components (P). |
1140 | ||
c947e71a | 1141 | Float_t spher = 0. ; |
1142 | Float_t emaxdtotal = 0. ; | |
50739f15 | 1143 | |
bc0c084c | 1144 | if((lambda[0]+lambda[1])!=0) |
c947e71a | 1145 | spher=fabs(lambda[0]-lambda[1])/(lambda[0]+lambda[1]); |
50739f15 | 1146 | |
c947e71a | 1147 | emaxdtotal=emc->GetMaximalEnergy()/emc->GetEnergy(); |
50739f15 | 1148 | |
1149 | fX[0] = lambda[0] ; | |
1150 | fX[1] = lambda[1] ; | |
1151 | fX[2] = emc->GetDispersion() ; | |
c947e71a | 1152 | fX[3] = spher ; |
50739f15 | 1153 | fX[4] = emc->GetMultiplicity() ; |
c947e71a | 1154 | fX[5] = emaxdtotal ; |
50739f15 | 1155 | fX[6] = emc->GetCoreEnergy() ; |
1156 | ||
e3817e5f | 1157 | fPrincipalPhoton->X2P(fX,fPPhoton); |
1158 | fPrincipalPi0 ->X2P(fX,fPPi0); | |
1f0e7ccd | 1159 | |
50739f15 | 1160 | } |
1161 | else{ | |
e3817e5f | 1162 | fPPhoton[0]=-100.0; //We do not accept clusters with |
1163 | fPPhoton[1]=-100.0; //one cell as a photon-like | |
1164 | fPPi0[0] =-100.0; | |
1165 | fPPi0[1] =-100.0; | |
50739f15 | 1166 | } |
1167 | ||
2cc71c1e | 1168 | Float_t time = emc->GetTime() ; |
1169 | rp->SetTof(time) ; | |
9fa5f1d0 | 1170 | |
bc0c084c | 1171 | // Loop of Efficiency-Purity (the 3 points of purity or efficiency |
1172 | // are taken into account to set the particle identification) | |
e3817e5f | 1173 | for(Int_t effPur = 0; effPur < 3 ; effPur++){ |
50739f15 | 1174 | |
bc0c084c | 1175 | // Looking at the CPV detector. If RCPV greater than CpvEmcDistance, |
1176 | // 1st,2nd or 3rd bit (depending on the efficiency-purity point ) | |
1177 | // is set to 1 | |
35adb638 | 1178 | if(GetCPVBit(emc, cpv, effPur,e) == 1 ){ |
e3817e5f | 1179 | rp->SetPIDBit(effPur) ; |
35adb638 | 1180 | //cout<<"CPV bit "<<effPur<<endl; |
1181 | } | |
50739f15 | 1182 | // Looking the TOF. If TOF smaller than gate, 4th, 5th or 6th |
1183 | // bit (depending on the efficiency-purity point )is set to 1 | |
2cc71c1e | 1184 | if(time< (*fParameters)(3,effPur)) |
e3817e5f | 1185 | rp->SetPIDBit(effPur+3) ; |
2cc71c1e | 1186 | |
e3817e5f | 1187 | //Photon PCA |
50739f15 | 1188 | //If we are inside the ellipse, 7th, 8th or 9th |
1189 | // bit (depending on the efficiency-purity point )is set to 1 | |
e3817e5f | 1190 | if(GetPrincipalBit("photon",fPPhoton,effPur,e) == 1) |
1191 | rp->SetPIDBit(effPur+6) ; | |
1f0e7ccd | 1192 | |
e3817e5f | 1193 | //Pi0 PCA |
1f0e7ccd | 1194 | //If we are inside the ellipse, 10th, 11th or 12th |
1195 | // bit (depending on the efficiency-purity point )is set to 1 | |
e3817e5f | 1196 | if(GetPrincipalBit("pi0" ,fPPi0 ,effPur,e) == 1) |
1197 | rp->SetPIDBit(effPur+9) ; | |
f0a4c9e9 | 1198 | } |
e3817e5f | 1199 | if(GetHardPhotonBit(emc)) |
1200 | rp->SetPIDBit(12) ; | |
1201 | if(GetHardPi0Bit (emc)) | |
1202 | rp->SetPIDBit(13) ; | |
1f0e7ccd | 1203 | |
bd76890a | 1204 | if(track >= 0) |
1205 | rp->SetPIDBit(14) ; | |
1206 | ||
9fa5f1d0 | 1207 | //Set momentum, energy and other parameters |
50739f15 | 1208 | Float_t encal = GetCalibratedEnergy(e); |
9fa5f1d0 | 1209 | TVector3 dir = GetMomentumDirection(emc,cpv) ; |
1210 | dir.SetMag(encal) ; | |
1211 | rp->SetMomentum(dir.X(),dir.Y(),dir.Z(),encal) ; | |
1212 | rp->SetCalcMass(0); | |
e0ed2e49 | 1213 | rp->Name(); //If photon sets the particle pdg name to gamma |
e747b8da | 1214 | rp->SetProductionVertex(0,0,0,0); |
1215 | rp->SetFirstMother(-1); | |
1216 | rp->SetLastMother(-1); | |
1217 | rp->SetFirstDaughter(-1); | |
1218 | rp->SetLastDaughter(-1); | |
1219 | rp->SetPolarisation(0,0,0); | |
d956e9b7 | 1220 | //Set the position in global coordinate system from the RecPoint |
1221 | AliPHOSGeometry * geom = gime->PHOSGeometry() ; | |
1222 | AliPHOSTrackSegment * ts = gime->TrackSegment(rp->GetPHOSTSIndex()) ; | |
1223 | AliPHOSEmcRecPoint * erp = gime->EmcRecPoint(ts->GetEmcIndex()) ; | |
1224 | TVector3 pos ; | |
1225 | geom->GetGlobal(erp, pos) ; | |
1226 | rp->SetPos(pos); | |
6ad0bfa0 | 1227 | index++ ; |
1228 | } | |
6ad0bfa0 | 1229 | } |
e3817e5f | 1230 | |
09fc14a0 | 1231 | //____________________________________________________________________________ |
88cb7938 | 1232 | void AliPHOSPIDv1::Print() const |
09fc14a0 | 1233 | { |
b2a60966 | 1234 | // Print the parameters used for the particle type identification |
bc0c084c | 1235 | |
351dd634 | 1236 | AliInfo("=============== AliPHOSPIDv1 ================") ; |
88cb7938 | 1237 | printf("Making PID\n") ; |
1238 | printf(" Pricipal analysis file from 0.5 to 100 %s\n", fFileNamePrincipalPhoton.Data() ) ; | |
1239 | printf(" Name of parameters file %s\n", fFileNameParameters.Data() ) ; | |
1240 | printf(" Matrix of Parameters: 14x4\n") ; | |
1241 | printf(" Energy Calibration 1x3 [3 parametres to calibrate energy: A + B* E + C * E^2]\n") ; | |
1242 | printf(" RCPV 2x3 rows x and z, columns function cut parameters\n") ; | |
1243 | printf(" TOF 1x3 [High Eff-Low Pur,Medium Eff-Pur, Low Eff-High Pur]\n") ; | |
1244 | printf(" PCA 5x4 [5 ellipse parametres and 4 parametres to calculate them: A/Sqrt(E) + B* E + C * E^2 + D]\n") ; | |
1245 | Printf(" Pi0 PCA 5x3 [5 ellipse parametres and 3 parametres to calculate them: A + B* E + C * E^2]\n") ; | |
50739f15 | 1246 | fParameters->Print() ; |
09fc14a0 | 1247 | } |
1248 | ||
a496c46c | 1249 | |
69183710 | 1250 | |
7acf6008 | 1251 | //____________________________________________________________________________ |
a4e98857 | 1252 | void AliPHOSPIDv1::PrintRecParticles(Option_t * option) |
1253 | { | |
dd5c4038 | 1254 | // Print table of reconstructed particles |
1255 | ||
88cb7938 | 1256 | AliPHOSGetter *gime = AliPHOSGetter::Instance() ; |
bf8f1fbd | 1257 | |
88cb7938 | 1258 | TClonesArray * recParticles = gime->RecParticles() ; |
21cd0c07 | 1259 | |
1260 | TString message ; | |
3bf72d32 | 1261 | message = "\nevent " ; |
1262 | message += gAlice->GetEvNumber() ; | |
1263 | message += " found " ; | |
1264 | message += recParticles->GetEntriesFast(); | |
1265 | message += " RecParticles\n" ; | |
1266 | ||
7acf6008 | 1267 | if(strstr(option,"all")) { // printing found TS |
3bf72d32 | 1268 | message += "\n PARTICLE Index \n" ; |
7acf6008 | 1269 | |
1270 | Int_t index ; | |
7b7c1533 | 1271 | for (index = 0 ; index < recParticles->GetEntries() ; index++) { |
21cd0c07 | 1272 | AliPHOSRecParticle * rp = (AliPHOSRecParticle * ) recParticles->At(index) ; |
3bf72d32 | 1273 | message += "\n" ; |
1274 | message += rp->Name().Data() ; | |
1275 | message += " " ; | |
1276 | message += rp->GetIndexInList() ; | |
1277 | message += " " ; | |
1278 | message += rp->GetType() ; | |
7acf6008 | 1279 | } |
3bf72d32 | 1280 | } |
351dd634 | 1281 | AliInfo(message.Data() ) ; |
69183710 | 1282 | } |
88cb7938 | 1283 | |
1284 | //____________________________________________________________________________ | |
1285 | void AliPHOSPIDv1::SetParameters() | |
1286 | { | |
1287 | // PCA : To do the Principal Components Analysis it is necessary | |
1288 | // the Principal file, which is opened here | |
1289 | fX = new double[7]; // Data for the PCA | |
1290 | fPPhoton = new double[7]; // Eigenvalues of the PCA | |
1291 | fPPi0 = new double[7]; // Eigenvalues of the Pi0 PCA | |
1292 | ||
1293 | // Read photon principals from the photon file | |
1294 | ||
1295 | fFileNamePrincipalPhoton = "$ALICE_ROOT/PHOS/PCA8pa15_0.5-100.root" ; | |
1296 | TFile f( fFileNamePrincipalPhoton.Data(), "read" ) ; | |
1297 | fPrincipalPhoton = dynamic_cast<TPrincipal*> (f.Get("principal")) ; | |
1298 | f.Close() ; | |
1299 | ||
1300 | // Read pi0 principals from the pi0 file | |
1301 | ||
1302 | fFileNamePrincipalPi0 = "$ALICE_ROOT/PHOS/PCA_pi0_40-120.root" ; | |
1303 | TFile fPi0( fFileNamePrincipalPi0.Data(), "read" ) ; | |
1304 | fPrincipalPi0 = dynamic_cast<TPrincipal*> (fPi0.Get("principal")) ; | |
1305 | fPi0.Close() ; | |
1306 | ||
1307 | // Open parameters file and initialization of the Parameters matrix. | |
1308 | // In the File Parameters.dat are all the parameters. These are introduced | |
1309 | // in a matrix of 16x4 | |
1310 | // | |
1311 | // All the parameters defined in this file are, in order of row: | |
1312 | // line 0 : calibration | |
1313 | // lines 1,2 : CPV rectangular cat for X and Z | |
1314 | // line 3 : TOF cut | |
1315 | // lines 4-8 : parameters to calculate photon PCA ellipse | |
1316 | // lines 9-13: parameters to calculate pi0 PCA ellipse | |
1317 | // lines 14-15: parameters to calculate border for high-pt photons and pi0 | |
1318 | ||
1319 | fFileNameParameters = gSystem->ExpandPathName("$ALICE_ROOT/PHOS/Parameters.dat"); | |
1320 | fParameters = new TMatrix(16,4) ; | |
c947e71a | 1321 | const Int_t kMaxLeng=255; |
1322 | char string[kMaxLeng]; | |
88cb7938 | 1323 | |
1324 | // Open a text file with PID parameters | |
1325 | FILE *fd = fopen(fFileNameParameters.Data(),"r"); | |
1326 | if (!fd) | |
351dd634 | 1327 | AliFatal(Form("File %s with a PID parameters cannot be opened\n", |
1328 | fFileNameParameters.Data())); | |
88cb7938 | 1329 | |
1330 | Int_t i=0; | |
1331 | // Read parameter file line-by-line and skip empty line and comments | |
c947e71a | 1332 | while (fgets(string,kMaxLeng,fd) != NULL) { |
88cb7938 | 1333 | if (string[0] == '\n' ) continue; |
1334 | if (string[0] == '!' ) continue; | |
1335 | sscanf(string, "%f %f %f %f", | |
1336 | &(*fParameters)(i,0), &(*fParameters)(i,1), | |
1337 | &(*fParameters)(i,2), &(*fParameters)(i,3)); | |
1338 | i++; | |
351dd634 | 1339 | AliDebug(1, Form("SetParameters", "line %d: %s",i,string)); |
88cb7938 | 1340 | } |
1341 | fclose(fd); | |
1342 | } | |
1343 | ||
1344 | //____________________________________________________________________________ | |
1345 | void AliPHOSPIDv1::SetParameterCalibration(Int_t i,Float_t param) | |
1346 | { | |
1347 | // Set parameter "Calibration" i to a value param | |
351dd634 | 1348 | if(i>2 || i<0) { |
1349 | AliError(Form("Invalid parameter number: %d",i)); | |
1350 | } else | |
88cb7938 | 1351 | (*fParameters)(0,i) = param ; |
1352 | } | |
1353 | ||
1354 | //____________________________________________________________________________ | |
1355 | void AliPHOSPIDv1::SetParameterCpv2Emc(Int_t i, TString axis, Float_t cut) | |
1356 | { | |
1357 | // Set the parameters to calculate Cpv-to-Emc Distance Cut depending on | |
1358 | // Purity-Efficiency point i | |
1359 | ||
351dd634 | 1360 | if(i>2 || i<0) { |
1361 | AliError(Form("Invalid parameter number: %d",i)); | |
1362 | } else { | |
88cb7938 | 1363 | axis.ToLower(); |
1364 | if (axis == "x") (*fParameters)(1,i) = cut; | |
1365 | else if (axis == "z") (*fParameters)(2,i) = cut; | |
351dd634 | 1366 | else { |
1367 | AliError(Form("Invalid axis name: %s",axis.Data())); | |
1368 | } | |
88cb7938 | 1369 | } |
1370 | } | |
1371 | ||
1372 | //____________________________________________________________________________ | |
1373 | void AliPHOSPIDv1::SetParameterPhotonBoundary(Int_t i,Float_t param) | |
1374 | { | |
1375 | // Set parameter "Hard photon boundary" i to a value param | |
351dd634 | 1376 | if(i>4 || i<0) { |
1377 | AliError(Form("Invalid parameter number: %d",i)); | |
1378 | } else | |
88cb7938 | 1379 | (*fParameters)(14,i) = param ; |
1380 | } | |
1381 | ||
1382 | //____________________________________________________________________________ | |
1383 | void AliPHOSPIDv1::SetParameterPi0Boundary(Int_t i,Float_t param) | |
1384 | { | |
1385 | // Set parameter "Hard pi0 boundary" i to a value param | |
351dd634 | 1386 | if(i>1 || i<0) { |
1387 | AliError(Form("Invalid parameter number: %d",i)); | |
1388 | } else | |
88cb7938 | 1389 | (*fParameters)(15,i) = param ; |
1390 | } | |
1391 | ||
1392 | //_____________________________________________________________________________ | |
1393 | void AliPHOSPIDv1::SetParameterTimeGate(Int_t i, Float_t gate) | |
1394 | { | |
1395 | // Set the parameter TimeGate depending on Purity-Efficiency point i | |
351dd634 | 1396 | if (i>2 || i<0) { |
1397 | AliError(Form("Invalid Efficiency-Purity choice %d",i)); | |
1398 | } else | |
88cb7938 | 1399 | (*fParameters)(3,i)= gate ; |
1400 | } | |
1401 | ||
1402 | //_____________________________________________________________________________ | |
1403 | void AliPHOSPIDv1::SetParameterToCalculateEllipse(TString particle, TString param, Int_t i, Float_t par) | |
1404 | { | |
1405 | // Set the parameter "i" that is needed to calculate the ellipse | |
1406 | // parameter "param" for a particle "particle" | |
1407 | ||
1408 | particle.ToLower(); | |
1409 | param. ToLower(); | |
1410 | Int_t p= -1; | |
1411 | Int_t offset=0; | |
1412 | ||
1413 | if (particle == "photon") offset=0; | |
1414 | else if (particle == "pi0") offset=5; | |
1415 | else | |
351dd634 | 1416 | AliError(Form("Wrong particle name: %s (choose from pi0/photon)\n", |
1417 | particle.Data())); | |
88cb7938 | 1418 | |
1419 | if (param.Contains("a")) p=4+offset; | |
1420 | else if(param.Contains("b")) p=5+offset; | |
1421 | else if(param.Contains("c")) p=6+offset; | |
1422 | else if(param.Contains("x0"))p=7+offset; | |
1423 | else if(param.Contains("y0"))p=8+offset; | |
351dd634 | 1424 | if((i>4)||(i<0)) { |
1425 | AliError(Form("No parameter with index %d", i)) ; | |
1426 | } else if(p==-1) { | |
1427 | AliError(Form("No parameter with name %s", param.Data() )) ; | |
1428 | } else | |
88cb7938 | 1429 | (*fParameters)(p,i) = par ; |
1430 | } | |
1431 | ||
1432 | //____________________________________________________________________________ | |
1433 | void AliPHOSPIDv1::Unload() | |
1434 | { | |
c947e71a | 1435 | //Unloads RecPoints, Tracks and RecParticles |
88cb7938 | 1436 | AliPHOSGetter * gime = AliPHOSGetter::Instance() ; |
1437 | gime->PhosLoader()->UnloadRecPoints() ; | |
1438 | gime->PhosLoader()->UnloadTracks() ; | |
1439 | gime->PhosLoader()->UnloadRecParticles() ; | |
1440 | } | |
1441 | ||
1442 | //____________________________________________________________________________ | |
90cceaf6 | 1443 | void AliPHOSPIDv1::WriteRecParticles() |
88cb7938 | 1444 | { |
c947e71a | 1445 | //It writes reconstructed particles and pid to file |
1446 | ||
88cb7938 | 1447 | AliPHOSGetter *gime = AliPHOSGetter::Instance() ; |
1448 | ||
1449 | TClonesArray * recParticles = gime->RecParticles() ; | |
1450 | recParticles->Expand(recParticles->GetEntriesFast() ) ; | |
adcca1e6 | 1451 | if(fWrite){ |
1452 | TTree * treeP = gime->TreeP(); | |
1453 | ||
1454 | //First rp | |
1455 | Int_t bufferSize = 32000 ; | |
1456 | TBranch * rpBranch = treeP->Branch("PHOSRP",&recParticles,bufferSize); | |
1457 | rpBranch->SetTitle(BranchName()); | |
1458 | ||
1459 | rpBranch->Fill() ; | |
1460 | ||
1461 | gime->WriteRecParticles("OVERWRITE"); | |
1462 | gime->WritePID("OVERWRITE"); | |
1463 | } | |
88cb7938 | 1464 | } |
1465 | ||
35adb638 | 1466 | |
1467 | //_______________________________________________________________________ | |
1468 | void AliPHOSPIDv1::SetInitPID(const Double_t *p) { | |
1469 | // Sets values for the initial population of each particle type | |
1470 | for (Int_t i=0; i<AliESDtrack::kSPECIESN; i++) fInitPID[i] = p[i]; | |
1471 | } | |
1472 | //_______________________________________________________________________ | |
1473 | void AliPHOSPIDv1::GetInitPID(Double_t *p) const { | |
1474 | // Gets values for the initial population of each particle type | |
1475 | for (Int_t i=0; i<AliESDtrack::kSPECIESN; i++) p[i] = fInitPID[i]; | |
1476 | } |