]>
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 | ||
702ab87e | 18 | /* History of cvs commits: |
19 | * | |
20 | * $Log$ | |
af885e0f | 21 | * Revision 1.111 2007/05/04 14:49:29 policheh |
22 | * AliPHOSRecPoint inheritance from AliCluster | |
23 | * | |
9ee9f78d | 24 | * Revision 1.110 2007/04/24 10:08:03 kharlov |
25 | * Vertex extraction from GenHeader | |
26 | * | |
2c06dc7a | 27 | * Revision 1.109 2007/04/18 09:34:05 kharlov |
28 | * Geometry bug fixes | |
29 | * | |
e5b7b511 | 30 | * Revision 1.108 2007/04/16 09:03:37 kharlov |
31 | * Incedent angle correction fixed | |
32 | * | |
753b19cd | 33 | * Revision 1.107 2007/04/02 15:00:16 cvetan |
34 | * No more calls to gAlice in the reconstruction | |
35 | * | |
999f9a8f | 36 | * Revision 1.106 2007/04/01 15:40:15 kharlov |
37 | * Correction for actual vertex position implemented | |
38 | * | |
407d15b3 | 39 | * Revision 1.105 2007/03/06 06:57:46 kharlov |
40 | * DP:calculation of distance to CPV done in TSM | |
41 | * | |
26aa7e4a | 42 | * Revision 1.104 2006/12/15 10:46:26 hristov |
43 | * Using TMath::Abs instead of fabs | |
44 | * | |
e4df4b30 | 45 | * Revision 1.103 2006/09/07 18:31:08 kharlov |
46 | * Effective c++ corrections (T.Pocheptsov) | |
47 | * | |
3f7dbdb7 | 48 | * Revision 1.102 2006/01/23 17:51:48 hristov |
49 | * Using the recommended way of forward declarations for TVector and TMatrix (see v5-08-00 release notes). Additional clean-up | |
50 | * | |
e8d02863 | 51 | * Revision 1.101 2005/05/28 14:19:04 schutz |
52 | * Compilation warnings fixed by T.P. | |
53 | * | |
702ab87e | 54 | */ |
55 | ||
6ad0bfa0 | 56 | //_________________________________________________________________________ |
b2a60966 | 57 | // Implementation version v1 of the PHOS particle identifier |
7acf6008 | 58 | // Particle identification based on the |
148b2bba | 59 | // - RCPV: distance from CPV recpoint to EMCA recpoint. |
60 | // - TOF | |
61 | // - PCA: Principal Components Analysis.. | |
62 | // The identified particle has an identification number corresponding | |
63 | // to a 9 bits number: | |
bc0c084c | 64 | // -Bit 0 to 2: bit set if RCPV > CpvEmcDistance (each bit corresponds |
148b2bba | 65 | // to a different efficiency-purity point of the photon identification) |
bc0c084c | 66 | // -Bit 3 to 5: bit set if TOF < TimeGate (each bit corresponds |
148b2bba | 67 | // to a different efficiency-purity point of the photon identification) |
68 | // -Bit 6 to 9: bit set if Principal Components are | |
50739f15 | 69 | // inside an ellipse defined by the parameters a, b, c, x0 and y0. |
148b2bba | 70 | // (each bit corresponds to a different efficiency-purity point of the |
50739f15 | 71 | // photon identification) |
72 | // The PCA (Principal components analysis) needs a file that contains | |
73 | // a previous analysis of the correlations between the particles. This | |
bc0c084c | 74 | // file is $ALICE_ROOT/PHOS/PCA8pa15_0.5-100.root. Analysis done for |
50739f15 | 75 | // energies between 0.5 and 100 GeV. |
9fa5f1d0 | 76 | // A calibrated energy is calculated. The energy of the reconstructed |
50739f15 | 77 | // cluster is corrected with the formula A + B * E + C * E^2, whose |
bc0c084c | 78 | // parameters where obtained through the study of the reconstructed |
50739f15 | 79 | // energy distribution of monoenergetic photons. |
a4e98857 | 80 | // |
bc0c084c | 81 | // All the parameters (RCPV(2 rows-3 columns),TOF(1r-3c),PCA(5r-4c) |
50739f15 | 82 | // and calibration(1r-3c))are stored in a file called |
83 | // $ALICE_ROOT/PHOS/Parameters.dat. Each time that AliPHOSPIDv1 is | |
bc0c084c | 84 | // initialized, this parameters are copied to a Matrix (9,4), a |
50739f15 | 85 | // TMatrixD object. |
7acf6008 | 86 | // |
a4e98857 | 87 | // use case: |
50739f15 | 88 | // root [0] AliPHOSPIDv1 * p = new AliPHOSPIDv1("galice1.root") |
a4e98857 | 89 | // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated |
50739f15 | 90 | // // reading headers from file galice1.root and create RecParticles |
91 | // TrackSegments and RecPoints are used | |
92 | // // set file name for the branch RecParticles | |
f0a4c9e9 | 93 | // root [1] p->ExecuteTask("deb all time") |
50739f15 | 94 | // // available options |
95 | // // "deb" - prints # of reconstructed particles | |
96 | // // "deb all" - prints # and list of RecParticles | |
97 | // // "time" - prints benchmarking results | |
7acf6008 | 98 | // |
50739f15 | 99 | // root [2] AliPHOSPIDv1 * p2 = new AliPHOSPIDv1("galice1.root","v1",kTRUE) |
148b2bba | 100 | // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated |
50739f15 | 101 | // //Split mode. |
f0a4c9e9 | 102 | // root [3] p2->ExecuteTask() |
103 | // | |
50739f15 | 104 | |
f0a4c9e9 | 105 | |
7acf6008 | 106 | //*-- Author: Yves Schutz (SUBATECH) & Gines Martinez (SUBATECH) & |
148b2bba | 107 | // Gustavo Conesa April 2002 |
50739f15 | 108 | // PCA redesigned by Gustavo Conesa October 2002: |
109 | // The way of using the PCA has changed. Instead of 2 | |
110 | // files with the PCA, each one with different energy ranges | |
111 | // of application, we use the wide one (0.5-100 GeV), and instead | |
bc0c084c | 112 | // of fixing 3 ellipses for different ranges of energy, it has been |
50739f15 | 113 | // studied the dependency of the ellipses parameters with the |
114 | // energy, and they are implemented in the code as a funtion | |
115 | // of the energy. | |
116 | // | |
117 | // | |
118 | // | |
6ad0bfa0 | 119 | // --- ROOT system --- |
c947e71a | 120 | |
121 | ||
122 | // --- Standard library --- | |
e8d02863 | 123 | #include <TMatrixF.h> |
acb5beb7 | 124 | #include "TFormula.h" |
7acf6008 | 125 | #include "TBenchmark.h" |
148b2bba | 126 | #include "TPrincipal.h" |
c947e71a | 127 | #include "TFile.h" |
e3817e5f | 128 | #include "TSystem.h" |
407d15b3 | 129 | #include "TVector3.h" |
148b2bba | 130 | |
6ad0bfa0 | 131 | // --- AliRoot header files --- |
c947e71a | 132 | //#include "AliLog.h" |
e3817e5f | 133 | #include "AliPHOS.h" |
26d4b141 | 134 | #include "AliPHOSPIDv1.h" |
7b7c1533 | 135 | #include "AliPHOSGetter.h" |
af885e0f | 136 | #include "AliESDEvent.h" |
407d15b3 | 137 | #include "AliESDVertex.h" |
2c06dc7a | 138 | #include "AliHeader.h" |
139 | #include "AliGenEventHeader.h" | |
6ad0bfa0 | 140 | |
26d4b141 | 141 | ClassImp( AliPHOSPIDv1) |
6ad0bfa0 | 142 | |
1cb7c1ee | 143 | //____________________________________________________________________________ |
3f7dbdb7 | 144 | AliPHOSPIDv1::AliPHOSPIDv1() : |
145 | fBayesian(kFALSE), | |
146 | fDefaultInit(kFALSE), | |
147 | fWrite(kFALSE), | |
148 | fNEvent(0), | |
149 | fFileNamePrincipalPhoton(), | |
150 | fFileNamePrincipalPi0(), | |
151 | fFileNameParameters(), | |
152 | fPrincipalPhoton(0), | |
153 | fPrincipalPi0(0), | |
154 | fX(0), | |
155 | fPPhoton(0), | |
156 | fPPi0(0), | |
157 | fRecParticlesInRun(0), | |
158 | fParameters(0), | |
159 | fTFphoton(0), | |
160 | fTFpiong(0), | |
161 | fTFkaong(0), | |
162 | fTFkaonl(0), | |
163 | fTFhhadrong(0), | |
164 | fTFhhadronl(0), | |
165 | fDFmuon(0), | |
166 | fERecWeight(0), | |
167 | fChargedNeutralThreshold(0.), | |
168 | fTOFEnThreshold(0), | |
169 | fDispEnThreshold(0), | |
170 | fDispMultThreshold(0) | |
1cb7c1ee | 171 | { |
a4e98857 | 172 | // default ctor |
148b2bba | 173 | |
8d0f3f77 | 174 | InitParameters() ; |
92f521a9 | 175 | fDefaultInit = kTRUE ; |
7acf6008 | 176 | } |
177 | ||
581354c5 | 178 | //____________________________________________________________________________ |
3f7dbdb7 | 179 | AliPHOSPIDv1::AliPHOSPIDv1(const AliPHOSPIDv1 & pid ) : |
180 | AliPHOSPID(pid), | |
181 | fBayesian(kFALSE), | |
182 | fDefaultInit(kFALSE), | |
183 | fWrite(kFALSE), | |
184 | fNEvent(0), | |
185 | fFileNamePrincipalPhoton(), | |
186 | fFileNamePrincipalPi0(), | |
187 | fFileNameParameters(), | |
188 | fPrincipalPhoton(0), | |
189 | fPrincipalPi0(0), | |
190 | fX(0), | |
191 | fPPhoton(0), | |
192 | fPPi0(0), | |
193 | fRecParticlesInRun(0), | |
194 | fParameters(0), | |
195 | fTFphoton(0), | |
196 | fTFpiong(0), | |
197 | fTFkaong(0), | |
198 | fTFkaonl(0), | |
199 | fTFhhadrong(0), | |
200 | fTFhhadronl(0), | |
201 | fDFmuon(0), | |
202 | fERecWeight(0), | |
203 | fChargedNeutralThreshold(0.), | |
204 | fTOFEnThreshold(0), | |
205 | fDispEnThreshold(0), | |
206 | fDispMultThreshold(0) | |
207 | ||
581354c5 | 208 | { |
386aef34 | 209 | // ctor |
581354c5 | 210 | InitParameters() ; |
581354c5 | 211 | Init() ; |
581354c5 | 212 | |
213 | } | |
214 | ||
7acf6008 | 215 | //____________________________________________________________________________ |
3f7dbdb7 | 216 | AliPHOSPIDv1::AliPHOSPIDv1(const TString alirunFileName, const TString eventFolderName) : |
217 | AliPHOSPID(alirunFileName, eventFolderName), | |
218 | fBayesian(kFALSE), | |
219 | fDefaultInit(kFALSE), | |
220 | fWrite(kFALSE), | |
221 | fNEvent(0), | |
222 | fFileNamePrincipalPhoton(), | |
223 | fFileNamePrincipalPi0(), | |
224 | fFileNameParameters(), | |
225 | fPrincipalPhoton(0), | |
226 | fPrincipalPi0(0), | |
227 | fX(0), | |
228 | fPPhoton(0), | |
229 | fPPi0(0), | |
230 | fRecParticlesInRun(0), | |
231 | fParameters(0), | |
232 | fTFphoton(0), | |
233 | fTFpiong(0), | |
234 | fTFkaong(0), | |
235 | fTFkaonl(0), | |
236 | fTFhhadrong(0), | |
237 | fTFhhadronl(0), | |
238 | fDFmuon(0), | |
239 | fERecWeight(0), | |
240 | fChargedNeutralThreshold(0.), | |
241 | fTOFEnThreshold(0), | |
242 | fDispEnThreshold(0), | |
243 | fDispMultThreshold(0) | |
244 | ||
7acf6008 | 245 | { |
a4e98857 | 246 | //ctor with the indication on where to look for the track segments |
7b7c1533 | 247 | |
8d0f3f77 | 248 | InitParameters() ; |
2bd5457f | 249 | Init() ; |
92f521a9 | 250 | fDefaultInit = kFALSE ; |
7acf6008 | 251 | } |
7b7c1533 | 252 | |
7acf6008 | 253 | //____________________________________________________________________________ |
254 | AliPHOSPIDv1::~AliPHOSPIDv1() | |
255 | { | |
79bb1b62 | 256 | // dtor |
acb5beb7 | 257 | fPrincipalPhoton = 0; |
258 | fPrincipalPi0 = 0; | |
9fa5f1d0 | 259 | |
e3817e5f | 260 | delete [] fX ; // Principal input |
261 | delete [] fPPhoton ; // Photon Principal components | |
262 | delete [] fPPi0 ; // Pi0 Principal components | |
acb5beb7 | 263 | |
264 | delete fParameters; | |
265 | delete fTFphoton; | |
266 | delete fTFpiong; | |
267 | delete fTFkaong; | |
268 | delete fTFkaonl; | |
269 | delete fTFhhadrong; | |
270 | delete fTFhhadronl; | |
271 | delete fDFmuon; | |
7acf6008 | 272 | } |
a496c46c | 273 | //____________________________________________________________________________ |
274 | const TString AliPHOSPIDv1::BranchName() const | |
275 | { | |
88cb7938 | 276 | |
277 | return GetName() ; | |
a496c46c | 278 | } |
279 | ||
148b2bba | 280 | //____________________________________________________________________________ |
281 | void AliPHOSPIDv1::Init() | |
282 | { | |
283 | // Make all memory allocations that are not possible in default constructor | |
284 | // Add the PID task to the list of PHOS tasks | |
a496c46c | 285 | |
adcca1e6 | 286 | AliPHOSGetter * gime = AliPHOSGetter::Instance() ; |
287 | if(!gime) | |
288 | gime = AliPHOSGetter::Instance(GetTitle(), fEventFolderName.Data()) ; | |
88cb7938 | 289 | |
290 | if ( !gime->PID() ) | |
291 | gime->PostPID(this) ; | |
148b2bba | 292 | } |
8d0f3f77 | 293 | |
294 | //____________________________________________________________________________ | |
295 | void AliPHOSPIDv1::InitParameters() | |
296 | { | |
e3817e5f | 297 | // Initialize PID parameters |
adcca1e6 | 298 | fWrite = kTRUE ; |
8d0f3f77 | 299 | fRecParticlesInRun = 0 ; |
8d0f3f77 | 300 | fNEvent = 0 ; |
8d0f3f77 | 301 | fRecParticlesInRun = 0 ; |
35adb638 | 302 | fBayesian = kTRUE ; |
9fa5f1d0 | 303 | SetParameters() ; // fill the parameters matrix from parameters file |
eabde521 | 304 | SetEventRange(0,-1) ; |
35adb638 | 305 | |
cc1fe362 | 306 | // initialisation of response function parameters |
307 | // Tof | |
2924941c | 308 | |
309 | // // Photons | |
310 | // fTphoton[0] = 0.218 ; | |
311 | // fTphoton[1] = 1.55E-8 ; | |
312 | // fTphoton[2] = 5.05E-10 ; | |
313 | // fTFphoton = new TFormula("ToF response to photons" , "gaus") ; | |
314 | // fTFphoton->SetParameters( fTphoton[0], fTphoton[1], fTphoton[2]) ; | |
315 | ||
316 | // // Pions | |
317 | // //Gaus (0 to max probability) | |
318 | // fTpiong[0] = 0.0971 ; | |
319 | // fTpiong[1] = 1.58E-8 ; | |
320 | // fTpiong[2] = 5.69E-10 ; | |
321 | // fTFpiong = new TFormula("ToF response to pions" , "gaus") ; | |
322 | // fTFpiong->SetParameters( fTpiong[0], fTpiong[1], fTpiong[2]) ; | |
323 | ||
324 | // // Kaons | |
325 | // //Gaus (0 to max probability) | |
326 | // fTkaong[0] = 0.0542 ; | |
327 | // fTkaong[1] = 1.64E-8 ; | |
328 | // fTkaong[2] = 6.07E-10 ; | |
329 | // fTFkaong = new TFormula("ToF response to kaon" , "gaus") ; | |
330 | // fTFkaong->SetParameters( fTkaong[0], fTkaong[1], fTkaong[2]) ; | |
331 | // //Landau (max probability to inf) | |
332 | // fTkaonl[0] = 0.264 ; | |
333 | // fTkaonl[1] = 1.68E-8 ; | |
334 | // fTkaonl[2] = 4.10E-10 ; | |
335 | // fTFkaonl = new TFormula("ToF response to kaon" , "landau") ; | |
336 | // fTFkaonl->SetParameters( fTkaonl[0], fTkaonl[1], fTkaonl[2]) ; | |
337 | ||
338 | // //Heavy Hadrons | |
339 | // //Gaus (0 to max probability) | |
340 | // fThhadrong[0] = 0.0302 ; | |
341 | // fThhadrong[1] = 1.73E-8 ; | |
342 | // fThhadrong[2] = 9.52E-10 ; | |
343 | // fTFhhadrong = new TFormula("ToF response to heavy hadrons" , "gaus") ; | |
344 | // fTFhhadrong->SetParameters( fThhadrong[0], fThhadrong[1], fThhadrong[2]) ; | |
345 | // //Landau (max probability to inf) | |
346 | // fThhadronl[0] = 0.139 ; | |
347 | // fThhadronl[1] = 1.745E-8 ; | |
348 | // fThhadronl[2] = 1.00E-9 ; | |
349 | // fTFhhadronl = new TFormula("ToF response to heavy hadrons" , "landau") ; | |
350 | // fTFhhadronl->SetParameters( fThhadronl[0], fThhadronl[1], fThhadronl[2]) ; | |
351 | ||
cc1fe362 | 352 | // Photons |
2924941c | 353 | fTphoton[0] = 7.83E8 ; |
35adb638 | 354 | fTphoton[1] = 1.55E-8 ; |
2924941c | 355 | fTphoton[2] = 5.09E-10 ; |
35adb638 | 356 | fTFphoton = new TFormula("ToF response to photons" , "gaus") ; |
cc1fe362 | 357 | fTFphoton->SetParameters( fTphoton[0], fTphoton[1], fTphoton[2]) ; |
35adb638 | 358 | |
359 | // Pions | |
360 | //Gaus (0 to max probability) | |
2924941c | 361 | fTpiong[0] = 6.73E8 ; |
35adb638 | 362 | fTpiong[1] = 1.58E-8 ; |
2924941c | 363 | fTpiong[2] = 5.87E-10 ; |
35adb638 | 364 | fTFpiong = new TFormula("ToF response to pions" , "gaus") ; |
365 | fTFpiong->SetParameters( fTpiong[0], fTpiong[1], fTpiong[2]) ; | |
35adb638 | 366 | |
367 | // Kaons | |
368 | //Gaus (0 to max probability) | |
2924941c | 369 | fTkaong[0] = 3.93E8 ; |
35adb638 | 370 | fTkaong[1] = 1.64E-8 ; |
2924941c | 371 | fTkaong[2] = 6.07E-10 ; |
35adb638 | 372 | fTFkaong = new TFormula("ToF response to kaon" , "gaus") ; |
373 | fTFkaong->SetParameters( fTkaong[0], fTkaong[1], fTkaong[2]) ; | |
374 | //Landau (max probability to inf) | |
2924941c | 375 | fTkaonl[0] = 2.0E9 ; |
35adb638 | 376 | fTkaonl[1] = 1.68E-8 ; |
377 | fTkaonl[2] = 4.10E-10 ; | |
378 | fTFkaonl = new TFormula("ToF response to kaon" , "landau") ; | |
379 | fTFkaonl->SetParameters( fTkaonl[0], fTkaonl[1], fTkaonl[2]) ; | |
380 | ||
381 | //Heavy Hadrons | |
382 | //Gaus (0 to max probability) | |
2924941c | 383 | fThhadrong[0] = 2.02E8 ; |
35adb638 | 384 | fThhadrong[1] = 1.73E-8 ; |
385 | fThhadrong[2] = 9.52E-10 ; | |
386 | fTFhhadrong = new TFormula("ToF response to heavy hadrons" , "gaus") ; | |
387 | fTFhhadrong->SetParameters( fThhadrong[0], fThhadrong[1], fThhadrong[2]) ; | |
388 | //Landau (max probability to inf) | |
2924941c | 389 | fThhadronl[0] = 1.10E9 ; |
390 | fThhadronl[1] = 1.74E-8 ; | |
391 | fThhadronl[2] = 1.00E-9 ; | |
35adb638 | 392 | fTFhhadronl = new TFormula("ToF response to heavy hadrons" , "landau") ; |
393 | fTFhhadronl->SetParameters( fThhadronl[0], fThhadronl[1], fThhadronl[2]) ; | |
394 | ||
35adb638 | 395 | |
396 | ||
397 | // Shower shape: dispersion gaussian parameters | |
398 | // Photons | |
2924941c | 399 | |
400 | // fDphoton[0] = 4.62e-2; fDphoton[1] = 1.39e-2 ; fDphoton[2] = -3.80e-2;//constant | |
401 | // fDphoton[3] = 1.53 ; fDphoton[4] =-6.62e-2 ; fDphoton[5] = 0.339 ;//mean | |
402 | // fDphoton[6] = 6.89e-2; fDphoton[7] =-6.59e-2 ; fDphoton[8] = 0.194 ;//sigma | |
403 | ||
404 | // fDpi0[0] = 0.0586 ; fDpi0[1] = 1.06E-3 ; fDpi0[2] = 0. ;//constant | |
405 | // fDpi0[3] = 2.67 ; fDpi0[4] =-2.00E-2 ; fDpi0[5] = 9.37E-5 ;//mean | |
406 | // fDpi0[6] = 0.153 ; fDpi0[7] = 9.34E-4 ; fDpi0[8] =-1.49E-5 ;//sigma | |
407 | ||
408 | // fDhadron[0] = 1.61E-2 ; fDhadron[1] = 3.03E-3 ; fDhadron[2] = 1.01E-2 ;//constant | |
409 | // fDhadron[3] = 3.81 ; fDhadron[4] = 0.232 ; fDhadron[5] =-1.25 ;//mean | |
410 | // fDhadron[6] = 0.897 ; fDhadron[7] = 0.0987 ; fDhadron[8] =-0.534 ;//sigma | |
411 | ||
412 | fDphoton[0] = 1.5 ; fDphoton[1] = 0.49 ; fDphoton[2] =-1.7E-2 ;//constant | |
413 | fDphoton[3] = 1.5 ; fDphoton[4] = 4.0E-2 ; fDphoton[5] = 0.21 ;//mean | |
414 | fDphoton[6] = 4.8E-2 ; fDphoton[7] =-0.12 ; fDphoton[8] = 0.27 ;//sigma | |
415 | fDphoton[9] = 16.; //for E> fDphoton[9] parameters calculated at fDphoton[9] | |
416 | ||
417 | fDpi0[0] = 0.25 ; fDpi0[1] = 3.3E-2 ; fDpi0[2] =-1.0e-5 ;//constant | |
418 | fDpi0[3] = 1.50 ; fDpi0[4] = 398. ; fDpi0[5] = 12. ;//mean | |
419 | fDpi0[6] =-7.0E-2 ; fDpi0[7] =-524. ; fDpi0[8] = 22. ;//sigma | |
420 | fDpi0[9] = 110.; //for E> fDpi0[9] parameters calculated at fDpi0[9] | |
421 | ||
422 | fDhadron[0] = 6.5 ; fDhadron[1] =-5.3 ; fDhadron[2] = 1.5 ;//constant | |
423 | fDhadron[3] = 3.8 ; fDhadron[4] = 0.23 ; fDhadron[5] =-1.2 ;//mean | |
424 | fDhadron[6] = 0.88 ; fDhadron[7] = 9.3E-2 ; fDhadron[8] =-0.51 ;//sigma | |
425 | fDhadron[9] = 2.; //for E> fDhadron[9] parameters calculated at fDhadron[9] | |
426 | ||
427 | fDmuon[0] = 0.0631 ; | |
428 | fDmuon[1] = 1.4 ; | |
35adb638 | 429 | fDmuon[2] = 0.0557 ; |
430 | fDFmuon = new TFormula("Shower shape response to muons" , "landau") ; | |
431 | fDFmuon->SetParameters( fDmuon[0], fDmuon[1], fDmuon[2]) ; | |
432 | ||
35adb638 | 433 | |
c947e71a | 434 | // x(CPV-EMC) distance gaussian parameters |
435 | ||
2924941c | 436 | // fXelectron[0] = 8.06e-2 ; fXelectron[1] = 1.00e-2; fXelectron[2] =-5.14e-2;//constant |
437 | // fXelectron[3] = 0.202 ; fXelectron[4] = 8.15e-3; fXelectron[5] = 4.55 ;//mean | |
438 | // fXelectron[6] = 0.334 ; fXelectron[7] = 0.186 ; fXelectron[8] = 4.32e-2;//sigma | |
c947e71a | 439 | |
2924941c | 440 | // //charged hadrons gaus |
441 | // fXcharged[0] = 6.43e-3 ; fXcharged[1] =-4.19e-5; fXcharged[2] = 1.42e-3;//constant | |
442 | // fXcharged[3] = 2.75 ; fXcharged[4] =-0.40 ; fXcharged[5] = 1.68 ;//mean | |
443 | // fXcharged[6] = 3.135 ; fXcharged[7] =-9.41e-2; fXcharged[8] = 1.31e-2;//sigma | |
c947e71a | 444 | |
2924941c | 445 | // // z(CPV-EMC) distance gaussian parameters |
446 | ||
447 | // fZelectron[0] = 8.22e-2 ; fZelectron[1] = 5.11e-3; fZelectron[2] =-3.05e-2;//constant | |
448 | // fZelectron[3] = 3.09e-2 ; fZelectron[4] = 5.87e-2; fZelectron[5] =-9.49e-2;//mean | |
449 | // fZelectron[6] = 0.263 ; fZelectron[7] =-9.02e-3; fZelectron[8] = 0.151 ;//sigma | |
c947e71a | 450 | |
2924941c | 451 | // //charged hadrons gaus |
c947e71a | 452 | |
2924941c | 453 | // fZcharged[0] = 1.00e-2 ; fZcharged[1] = 2.82E-4 ; fZcharged[2] = 2.87E-3 ;//constant |
454 | // fZcharged[3] =-4.68e-2 ; fZcharged[4] =-9.21e-3 ; fZcharged[5] = 4.91e-2 ;//mean | |
455 | // fZcharged[6] = 1.425 ; fZcharged[7] =-5.90e-2 ; fZcharged[8] = 5.07e-2 ;//sigma | |
456 | ||
457 | ||
458 | fXelectron[0] =-1.6E-2 ; fXelectron[1] = 0.77 ; fXelectron[2] =-0.15 ;//constant | |
459 | fXelectron[3] = 0.35 ; fXelectron[4] = 0.25 ; fXelectron[5] = 4.12 ;//mean | |
460 | fXelectron[6] = 0.30 ; fXelectron[7] = 0.11 ; fXelectron[8] = 0.16 ;//sigma | |
461 | fXelectron[9] = 3.; //for E> fXelectron[9] parameters calculated at fXelectron[9] | |
462 | ||
c947e71a | 463 | //charged hadrons gaus |
2924941c | 464 | fXcharged[0] = 0.14 ; fXcharged[1] =-3.0E-2 ; fXcharged[2] = 0 ;//constant |
465 | fXcharged[3] = 1.4 ; fXcharged[4] =-9.3E-2 ; fXcharged[5] = 1.4 ;//mean | |
466 | fXcharged[6] = 5.7 ; fXcharged[7] = 0.27 ; fXcharged[8] =-1.8 ;//sigma | |
467 | fXcharged[9] = 1.2; //for E> fXcharged[9] parameters calculated at fXcharged[9] | |
468 | ||
469 | // z(CPV-EMC) distance gaussian parameters | |
c947e71a | 470 | |
2924941c | 471 | fZelectron[0] = 0.49 ; fZelectron[1] = 0.53 ; fZelectron[2] =-9.8E-2 ;//constant |
472 | fZelectron[3] = 2.8E-2 ; fZelectron[4] = 5.0E-2 ; fZelectron[5] =-8.2E-2 ;//mean | |
473 | fZelectron[6] = 0.25 ; fZelectron[7] =-1.7E-2 ; fZelectron[8] = 0.17 ;//sigma | |
474 | fZelectron[9] = 3.; //for E> fZelectron[9] parameters calculated at fZelectron[9] | |
475 | ||
476 | //charged hadrons gaus | |
c947e71a | 477 | |
2924941c | 478 | fZcharged[0] = 0.46 ; fZcharged[1] =-0.65 ; fZcharged[2] = 0.52 ;//constant |
479 | fZcharged[3] = 1.1E-2 ; fZcharged[4] = 0. ; fZcharged[5] = 0. ;//mean | |
480 | fZcharged[6] = 0.60 ; fZcharged[7] =-8.2E-2 ; fZcharged[8] = 0.45 ;//sigma | |
481 | fZcharged[9] = 1.2; //for E> fXcharged[9] parameters calculated at fXcharged[9] | |
482 | ||
fb7b51ad | 483 | //Threshold to differentiate between charged and neutral |
484 | fChargedNeutralThreshold = 1e-5; | |
2924941c | 485 | fTOFEnThreshold = 2; //Maximum energy to use TOF |
486 | fDispEnThreshold = 0.5; //Minimum energy to use shower shape | |
487 | fDispMultThreshold = 3; //Minimum multiplicity to use shower shape | |
fb7b51ad | 488 | |
489 | //Weight to hadrons recontructed energy | |
490 | ||
491 | fERecWeightPar[0] = 0.32 ; | |
492 | fERecWeightPar[1] = 3.8 ; | |
493 | fERecWeightPar[2] = 5.4E-3 ; | |
494 | fERecWeightPar[3] = 5.6E-2 ; | |
495 | fERecWeight = new TFormula("Weight for hadrons" , "[0]*exp(-x*[1])+[2]*exp(-x*[3])") ; | |
496 | fERecWeight ->SetParameters(fERecWeightPar[0],fERecWeightPar[1] ,fERecWeightPar[2] ,fERecWeightPar[3]) ; | |
497 | ||
498 | ||
304864ab | 499 | for (Int_t i =0; i< AliPID::kSPECIESN ; i++) |
35adb638 | 500 | fInitPID[i] = 1.; |
fb7b51ad | 501 | |
8d0f3f77 | 502 | } |
503 | ||
88cb7938 | 504 | //________________________________________________________________________ |
eabde521 | 505 | void AliPHOSPIDv1::Exec(Option_t *option) |
88cb7938 | 506 | { |
eabde521 | 507 | // Steering method to perform particle reconstruction and identification |
508 | // for the event range from fFirstEvent to fLastEvent. | |
509 | // This range is optionally set by SetEventRange(). | |
510 | // if fLastEvent=-1 (by default), then process events until the end. | |
61d3d6aa | 511 | |
88cb7938 | 512 | if(strstr(option,"tim")) |
513 | gBenchmark->Start("PHOSPID"); | |
514 | ||
515 | if(strstr(option,"print")) { | |
516 | Print() ; | |
517 | return ; | |
518 | } | |
519 | ||
520 | ||
adcca1e6 | 521 | AliPHOSGetter * gime = AliPHOSGetter::Instance() ; |
88cb7938 | 522 | |
71cee46d | 523 | if (fLastEvent == -1) |
524 | fLastEvent = gime->MaxEvent() - 1 ; | |
525 | else | |
526 | fLastEvent = TMath::Min(fLastEvent,gime->MaxEvent()); | |
eabde521 | 527 | Int_t nEvents = fLastEvent - fFirstEvent + 1; |
88cb7938 | 528 | |
71cee46d | 529 | Int_t ievent ; |
530 | for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) { | |
88cb7938 | 531 | gime->Event(ievent,"TR") ; |
532 | if(gime->TrackSegments() && //Skip events, where no track segments made | |
533 | gime->TrackSegments()->GetEntriesFast()) { | |
7fb46731 | 534 | |
407d15b3 | 535 | GetVertex() ; |
88cb7938 | 536 | MakeRecParticles() ; |
35adb638 | 537 | |
538 | if(fBayesian) | |
539 | MakePID() ; | |
540 | ||
90cceaf6 | 541 | WriteRecParticles(); |
88cb7938 | 542 | if(strstr(option,"deb")) |
543 | PrintRecParticles(option) ; | |
544 | //increment the total number of rec particles per run | |
545 | fRecParticlesInRun += gime->RecParticles()->GetEntriesFast() ; | |
546 | } | |
547 | } | |
ff417097 | 548 | if(strstr(option,"deb")) |
549 | PrintRecParticles(option); | |
88cb7938 | 550 | if(strstr(option,"tim")){ |
551 | gBenchmark->Stop("PHOSPID"); | |
351dd634 | 552 | AliInfo(Form("took %f seconds for PID %f seconds per event", |
88cb7938 | 553 | gBenchmark->GetCpuTime("PHOSPID"), |
351dd634 | 554 | gBenchmark->GetCpuTime("PHOSPID")/nEvents)) ; |
88cb7938 | 555 | } |
adcca1e6 | 556 | if(fWrite) |
557 | Unload(); | |
88cb7938 | 558 | } |
559 | ||
35adb638 | 560 | //________________________________________________________________________ |
17323043 | 561 | Double_t AliPHOSPIDv1::GausF(Double_t x, Double_t y, Double_t * par) |
35adb638 | 562 | { |
c947e71a | 563 | //Given the energy x and the parameter y (tof, shower dispersion or cpv-emc distance), |
564 | //this method returns a density probability of this parameter, given by a gaussian | |
565 | //function whose parameters depend with the energy with a function: a/(x*x)+b/x+b | |
2924941c | 566 | //Float_t xorg = x; |
567 | if (x > par[9]) x = par[9]; | |
568 | ||
569 | //Double_t cnt = par[1] / (x*x) + par[2] / x + par[0] ; | |
570 | Double_t cnt = par[0] + par[1] * x + par[2] * x * x ; | |
35adb638 | 571 | Double_t mean = par[4] / (x*x) + par[5] / x + par[3] ; |
572 | Double_t sigma = par[7] / (x*x) + par[8] / x + par[6] ; | |
c947e71a | 573 | |
2924941c | 574 | // if(xorg > 30) |
575 | // cout<<"En_in = "<<xorg<<"; En_out = "<<x<<"; cnt = "<<cnt | |
576 | // <<"; mean = "<<mean<<"; sigma = "<<sigma<<endl; | |
577 | ||
35adb638 | 578 | // Double_t arg = - (y-mean) * (y-mean) / (2*sigma*sigma) ; |
579 | // return cnt * TMath::Exp(arg) ; | |
c947e71a | 580 | if(TMath::Abs(sigma) > 1.e-10){ |
acb5beb7 | 581 | return cnt*TMath::Gaus(y,mean,sigma); |
35adb638 | 582 | } |
583 | else | |
584 | return 0.; | |
c947e71a | 585 | |
35adb638 | 586 | } |
587 | //________________________________________________________________________ | |
17323043 | 588 | Double_t AliPHOSPIDv1::GausPol2(Double_t x, Double_t y, Double_t * par) |
35adb638 | 589 | { |
c947e71a | 590 | //Given the energy x and the parameter y (tof, shower dispersion or cpv-emc distance), |
591 | //this method returns a density probability of this parameter, given by a gaussian | |
592 | //function whose parameters depend with the energy like second order polinomial | |
593 | ||
35adb638 | 594 | Double_t cnt = par[0] + par[1] * x + par[2] * x * x ; |
595 | Double_t mean = par[3] + par[4] * x + par[5] * x * x ; | |
596 | Double_t sigma = par[6] + par[7] * x + par[8] * x * x ; | |
597 | ||
c947e71a | 598 | if(TMath::Abs(sigma) > 1.e-10){ |
acb5beb7 | 599 | return cnt*TMath::Gaus(y,mean,sigma); |
35adb638 | 600 | } |
601 | else | |
602 | return 0.; | |
c947e71a | 603 | |
604 | ||
605 | ||
35adb638 | 606 | } |
607 | ||
69183710 | 608 | //____________________________________________________________________________ |
e3817e5f | 609 | const TString AliPHOSPIDv1::GetFileNamePrincipal(TString particle) const |
148b2bba | 610 | { |
e3817e5f | 611 | //Get file name that contains the PCA for a particle ("photon or pi0") |
612 | particle.ToLower(); | |
613 | TString name; | |
351dd634 | 614 | if (particle=="photon") |
615 | name = fFileNamePrincipalPhoton ; | |
616 | else if (particle=="pi0" ) | |
617 | name = fFileNamePrincipalPi0 ; | |
618 | else | |
619 | AliError(Form("Wrong particle name: %s (choose from pi0/photon)\n", | |
620 | particle.Data())); | |
e3817e5f | 621 | return name; |
622 | } | |
bc0c084c | 623 | |
e3817e5f | 624 | //____________________________________________________________________________ |
fc7e2f43 | 625 | Float_t AliPHOSPIDv1::GetParameterCalibration(Int_t i) const |
e3817e5f | 626 | { |
627 | // Get the i-th parameter "Calibration" | |
628 | Float_t param = 0.; | |
351dd634 | 629 | if (i>2 || i<0) { |
630 | AliError(Form("Invalid parameter number: %d",i)); | |
631 | } else | |
e3817e5f | 632 | param = (*fParameters)(0,i); |
633 | return param; | |
634 | } | |
bc0c084c | 635 | |
88cb7938 | 636 | //____________________________________________________________________________ |
fc7e2f43 | 637 | Float_t AliPHOSPIDv1::GetCalibratedEnergy(Float_t e) const |
88cb7938 | 638 | { |
639 | // It calibrates Energy depending on the recpoint energy. | |
640 | // The energy of the reconstructed cluster is corrected with | |
641 | // the formula A + B* E + C* E^2, whose parameters where obtained | |
642 | // through the study of the reconstructed energy distribution of | |
643 | // monoenergetic photons. | |
644 | ||
645 | Float_t p[]={0.,0.,0.}; | |
646 | for (Int_t i=0; i<3; i++) p[i] = GetParameterCalibration(i); | |
647 | Float_t enerec = p[0] + p[1]*e + p[2]*e*e; | |
648 | return enerec ; | |
649 | ||
650 | } | |
651 | ||
e3817e5f | 652 | //____________________________________________________________________________ |
fc7e2f43 | 653 | Float_t AliPHOSPIDv1::GetParameterCpv2Emc(Int_t i, TString axis) const |
e3817e5f | 654 | { |
655 | // Get the i-th parameter "CPV-EMC distance" for the specified axis | |
656 | Float_t param = 0.; | |
351dd634 | 657 | if(i>2 || i<0) { |
658 | AliError(Form("Invalid parameter number: %d",i)); | |
659 | } else { | |
e3817e5f | 660 | axis.ToLower(); |
351dd634 | 661 | if (axis == "x") |
662 | param = (*fParameters)(1,i); | |
663 | else if (axis == "z") | |
664 | param = (*fParameters)(2,i); | |
665 | else { | |
666 | AliError(Form("Invalid axis name: %s",axis.Data())); | |
667 | } | |
e3817e5f | 668 | } |
669 | return param; | |
670 | } | |
671 | ||
672 | //____________________________________________________________________________ | |
fc7e2f43 | 673 | Float_t AliPHOSPIDv1::GetCpv2EmcDistanceCut(TString axis, Float_t e) const |
e3817e5f | 674 | { |
88cb7938 | 675 | // Get CpvtoEmcDistance Cut depending on the cluster energy, axis and |
676 | // Purity-Efficiency point | |
677 | ||
678 | axis.ToLower(); | |
679 | Float_t p[]={0.,0.,0.}; | |
680 | for (Int_t i=0; i<3; i++) p[i] = GetParameterCpv2Emc(i,axis); | |
681 | Float_t sig = p[0] + TMath::Exp(p[1] - p[2]*e); | |
682 | return sig; | |
e3817e5f | 683 | } |
684 | ||
88cb7938 | 685 | //____________________________________________________________________________ |
fc7e2f43 | 686 | Float_t AliPHOSPIDv1::GetEllipseParameter(TString particle, TString param, Float_t e) const |
88cb7938 | 687 | { |
688 | // Calculates the parameter param of the ellipse | |
e3817e5f | 689 | |
690 | particle.ToLower(); | |
691 | param. ToLower(); | |
88cb7938 | 692 | Float_t p[4]={0.,0.,0.,0.}; |
693 | Float_t value = 0.0; | |
694 | for (Int_t i=0; i<4; i++) p[i] = GetParameterToCalculateEllipse(particle,param,i); | |
695 | if (particle == "photon") { | |
696 | if (param.Contains("a")) e = TMath::Min((Double_t)e,70.); | |
697 | else if (param.Contains("b")) e = TMath::Min((Double_t)e,70.); | |
698 | else if (param.Contains("x0")) e = TMath::Max((Double_t)e,1.1); | |
699 | } | |
e3817e5f | 700 | |
443caba9 | 701 | if (particle == "photon") |
702 | value = p[0]/TMath::Sqrt(e) + p[1]*e + p[2]*e*e + p[3]; | |
703 | else if (particle == "pi0") | |
704 | value = p[0] + p[1]*e + p[2]*e*e; | |
705 | ||
88cb7938 | 706 | return value; |
e3817e5f | 707 | } |
708 | ||
709 | //_____________________________________________________________________________ | |
fc7e2f43 | 710 | Float_t AliPHOSPIDv1::GetParameterPhotonBoundary (Int_t i) const |
e3817e5f | 711 | { |
712 | // Get the parameter "i" to calculate the boundary on the moment M2x | |
713 | // for photons at high p_T | |
714 | Float_t param = 0; | |
351dd634 | 715 | if (i>3 || i<0) { |
716 | AliError(Form("Wrong parameter number: %d\n",i)); | |
717 | } else | |
e3817e5f | 718 | param = (*fParameters)(14,i) ; |
719 | return param; | |
148b2bba | 720 | } |
e3817e5f | 721 | |
148b2bba | 722 | //____________________________________________________________________________ |
fc7e2f43 | 723 | Float_t AliPHOSPIDv1::GetParameterPi0Boundary (Int_t i) const |
e3817e5f | 724 | { |
725 | // Get the parameter "i" to calculate the boundary on the moment M2x | |
726 | // for pi0 at high p_T | |
727 | Float_t param = 0; | |
351dd634 | 728 | if (i>2 || i<0) { |
729 | AliError(Form("Wrong parameter number: %d\n",i)); | |
730 | } else | |
e3817e5f | 731 | param = (*fParameters)(15,i) ; |
732 | return param; | |
733 | } | |
148b2bba | 734 | |
e3817e5f | 735 | //____________________________________________________________________________ |
fc7e2f43 | 736 | Float_t AliPHOSPIDv1::GetParameterTimeGate(Int_t i) const |
e3817e5f | 737 | { |
88cb7938 | 738 | // Get TimeGate parameter depending on Purity-Efficiency i: |
739 | // i=0 - Low purity, i=1 - Medium purity, i=2 - High purity | |
740 | Float_t param = 0.; | |
351dd634 | 741 | if(i>2 || i<0) { |
742 | AliError(Form("Invalid Efficiency-Purity choice %d",i)); | |
743 | } else | |
88cb7938 | 744 | param = (*fParameters)(3,i) ; |
745 | return param; | |
e3817e5f | 746 | } |
747 | ||
148b2bba | 748 | //_____________________________________________________________________________ |
fc7e2f43 | 749 | Float_t AliPHOSPIDv1::GetParameterToCalculateEllipse(TString particle, TString param, Int_t i) const |
88cb7938 | 750 | { |
751 | // Get the parameter "i" that is needed to calculate the ellipse | |
752 | // parameter "param" for the particle "particle" ("photon" or "pi0") | |
753 | ||
e3817e5f | 754 | particle.ToLower(); |
755 | param. ToLower(); | |
88cb7938 | 756 | Int_t offset = -1; |
351dd634 | 757 | if (particle == "photon") |
758 | offset=0; | |
759 | else if (particle == "pi0") | |
760 | offset=5; | |
e3817e5f | 761 | else |
351dd634 | 762 | AliError(Form("Wrong particle name: %s (choose from pi0/photon)\n", |
763 | particle.Data())); | |
88cb7938 | 764 | |
765 | Int_t p= -1; | |
766 | Float_t par = 0; | |
e3817e5f | 767 | |
768 | if (param.Contains("a")) p=4+offset; | |
769 | else if(param.Contains("b")) p=5+offset; | |
770 | else if(param.Contains("c")) p=6+offset; | |
771 | else if(param.Contains("x0"))p=7+offset; | |
772 | else if(param.Contains("y0"))p=8+offset; | |
12022e83 | 773 | |
351dd634 | 774 | if (i>4 || i<0) { |
775 | AliError(Form("No parameter with index %d", i)) ; | |
776 | } else if (p==-1) { | |
777 | AliError(Form("No parameter with name %s", param.Data() )) ; | |
778 | } else | |
88cb7938 | 779 | par = (*fParameters)(p,i) ; |
780 | ||
781 | return par; | |
12022e83 | 782 | } |
783 | ||
12022e83 | 784 | |
26aa7e4a | 785 | //DP____________________________________________________________________________ |
786 | //Float_t AliPHOSPIDv1::GetDistance(AliPHOSEmcRecPoint * emc,AliPHOSCpvRecPoint * cpv, Option_t * axis)const | |
787 | //{ | |
788 | // // Calculates the distance between the EMC RecPoint and the PPSD RecPoint | |
789 | // | |
790 | // const AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ; | |
791 | // TVector3 vecEmc ; | |
792 | // TVector3 vecCpv ; | |
793 | // if(cpv){ | |
794 | // emc->GetLocalPosition(vecEmc) ; | |
795 | // cpv->GetLocalPosition(vecCpv) ; | |
796 | // | |
797 | // if(emc->GetPHOSMod() == cpv->GetPHOSMod()){ | |
798 | // // Correct to difference in CPV and EMC position due to different distance to center. | |
799 | // // we assume, that particle moves from center | |
800 | // Float_t dCPV = geom->GetIPtoOuterCoverDistance(); | |
801 | // Float_t dEMC = geom->GetIPtoCrystalSurface() ; | |
802 | // dEMC = dEMC / dCPV ; | |
803 | // vecCpv = dEMC * vecCpv - vecEmc ; | |
804 | // if (axis == "X") return vecCpv.X(); | |
805 | // if (axis == "Y") return vecCpv.Y(); | |
806 | // if (axis == "Z") return vecCpv.Z(); | |
807 | // if (axis == "R") return vecCpv.Mag(); | |
808 | // } | |
809 | // return 100000000 ; | |
810 | // } | |
811 | // return 100000000 ; | |
812 | //} | |
12022e83 | 813 | //____________________________________________________________________________ |
26aa7e4a | 814 | Int_t AliPHOSPIDv1::GetCPVBit(AliPHOSTrackSegment * ts, Int_t effPur, Float_t e) const |
bc0c084c | 815 | { |
c947e71a | 816 | //Calculates the pid bit for the CPV selection per each purity. |
e3817e5f | 817 | if(effPur>2 || effPur<0) |
351dd634 | 818 | AliError(Form("Invalid Efficiency-Purity choice %d",effPur)); |
26aa7e4a | 819 | |
407d15b3 | 820 | //DP if(ts->GetCpvIndex()<0) |
821 | //DP return 1 ; //no CPV cluster | |
bc0c084c | 822 | |
e3817e5f | 823 | Float_t sigX = GetCpv2EmcDistanceCut("X",e); |
824 | Float_t sigZ = GetCpv2EmcDistanceCut("Z",e); | |
bc0c084c | 825 | |
26aa7e4a | 826 | Float_t deltaX = TMath::Abs(ts->GetCpvDistance("X")); |
827 | Float_t deltaZ = TMath::Abs(ts->GetCpvDistance("Z")); | |
407d15b3 | 828 | // Info("GetCPVBit"," xdist %f, sigx %f, zdist %f, sigz %f",deltaX, sigX, deltaZ,sigZ) ; |
7fb46731 | 829 | |
830 | //if(deltaX>sigX*(effPur+1)) | |
831 | //if((deltaX>sigX*(effPur+1)) || (deltaZ>sigZ*(effPur+1))) | |
832 | if((deltaX>sigX*(effPur+1)) && (deltaZ>sigZ*(effPur+1))) | |
bc0c084c | 833 | return 1;//Neutral |
834 | else | |
835 | return 0;//Charged | |
bc0c084c | 836 | } |
69183710 | 837 | |
6ad0bfa0 | 838 | //____________________________________________________________________________ |
fc7e2f43 | 839 | Int_t AliPHOSPIDv1::GetPrincipalBit(TString particle, const Double_t* p, Int_t effPur, Float_t e)const |
148b2bba | 840 | { |
50739f15 | 841 | //Is the particle inside de PCA ellipse? |
581354c5 | 842 | |
e3817e5f | 843 | particle.ToLower(); |
844 | Int_t prinbit = 0 ; | |
7fb46731 | 845 | Float_t a = GetEllipseParameter(particle,"a" , e); |
846 | Float_t b = GetEllipseParameter(particle,"b" , e); | |
847 | Float_t c = GetEllipseParameter(particle,"c" , e); | |
e3817e5f | 848 | Float_t x0 = GetEllipseParameter(particle,"x0", e); |
849 | Float_t y0 = GetEllipseParameter(particle,"y0", e); | |
850 | ||
851 | Float_t r = TMath::Power((p[0] - x0)/a,2) + | |
852 | TMath::Power((p[1] - y0)/b,2) + | |
853 | c*(p[0] - x0)*(p[1] - y0)/(a*b) ; | |
50739f15 | 854 | //3 different ellipses defined |
e3817e5f | 855 | if((effPur==2) && (r<1./2.)) prinbit= 1; |
856 | if((effPur==1) && (r<2. )) prinbit= 1; | |
857 | if((effPur==0) && (r<9./2.)) prinbit= 1; | |
50739f15 | 858 | |
581354c5 | 859 | if(r<0) |
351dd634 | 860 | AliError("Negative square?") ; |
1f0e7ccd | 861 | |
862 | return prinbit; | |
148b2bba | 863 | |
148b2bba | 864 | } |
1f0e7ccd | 865 | //____________________________________________________________________________ |
fc7e2f43 | 866 | Int_t AliPHOSPIDv1::GetHardPhotonBit(AliPHOSEmcRecPoint * emc) const |
1f0e7ccd | 867 | { |
e3817e5f | 868 | // Set bit for identified hard photons (E > 30 GeV) |
869 | // if the second moment M2x is below the boundary | |
870 | ||
871 | Float_t e = emc->GetEnergy(); | |
872 | if (e < 30.0) return 0; | |
873 | Float_t m2x = emc->GetM2x(); | |
874 | Float_t m2xBoundary = GetParameterPhotonBoundary(0) * | |
875 | TMath::Exp(-TMath::Power(e-GetParameterPhotonBoundary(1),2)/2.0/ | |
876 | TMath::Power(GetParameterPhotonBoundary(2),2)) + | |
877 | GetParameterPhotonBoundary(3); | |
351dd634 | 878 | AliDebug(1, Form("GetHardPhotonBit","E=%f, m2x=%f, boundary=%f", |
879 | e,m2x,m2xBoundary)); | |
e3817e5f | 880 | if (m2x < m2xBoundary) |
881 | return 1;// A hard photon | |
882 | else | |
883 | return 0;// Not a hard photon | |
1f0e7ccd | 884 | } |
92f521a9 | 885 | |
e3817e5f | 886 | //____________________________________________________________________________ |
fc7e2f43 | 887 | Int_t AliPHOSPIDv1::GetHardPi0Bit(AliPHOSEmcRecPoint * emc) const |
e3817e5f | 888 | { |
889 | // Set bit for identified hard pi0 (E > 30 GeV) | |
890 | // if the second moment M2x is above the boundary | |
891 | ||
892 | Float_t e = emc->GetEnergy(); | |
893 | if (e < 30.0) return 0; | |
894 | Float_t m2x = emc->GetM2x(); | |
895 | Float_t m2xBoundary = GetParameterPi0Boundary(0) + | |
896 | e * GetParameterPi0Boundary(1); | |
351dd634 | 897 | AliDebug(1,Form("E=%f, m2x=%f, boundary=%f",e,m2x,m2xBoundary)); |
e3817e5f | 898 | if (m2x > m2xBoundary) |
899 | return 1;// A hard pi0 | |
bc0c084c | 900 | else |
e3817e5f | 901 | return 0;// Not a hard pi0 |
f0a4c9e9 | 902 | } |
e3817e5f | 903 | |
904 | //____________________________________________________________________________ | |
8d4608b5 | 905 | TVector3 AliPHOSPIDv1::GetMomentumDirection(AliPHOSEmcRecPoint * emc, AliPHOSCpvRecPoint * )const |
88cb7938 | 906 | { |
907 | // Calculates the momentum direction: | |
908 | // 1. if only a EMC RecPoint, direction is given by IP and this RecPoint | |
909 | // 2. if a EMC RecPoint and CPV RecPoint, direction is given by the line through the 2 recpoints | |
910 | // However because of the poor position resolution of PPSD the direction is always taken as if we were | |
911 | // in case 1. | |
f0a4c9e9 | 912 | |
407d15b3 | 913 | TVector3 local ; |
914 | emc->GetLocalPosition(local) ; | |
915 | ||
916 | AliPHOSGeometry * phosgeom = AliPHOSGeometry::GetInstance() ; | |
917 | //Correct for the non-perpendicular incidence | |
918 | // Correction for the depth of the shower starting point (TDR p 127) | |
919 | Float_t para = 0.925 ; | |
920 | Float_t parb = 6.52 ; | |
921 | ||
922 | //Remove Old correction (vertex at 0,0,0) | |
923 | TVector3 vtxOld(0.,0.,0.) ; | |
924 | TVector3 vInc ; | |
925 | Float_t x=local.X() ; | |
926 | Float_t z=local.Z() ; | |
927 | phosgeom->GetIncidentVector(vtxOld,emc->GetPHOSMod(),x,z,vInc) ; | |
928 | Float_t depthxOld = 0.; | |
929 | Float_t depthzOld = 0.; | |
930 | Float_t energy = emc->GetEnergy() ; | |
931 | if (energy > 0 && vInc.Y()!=0.) { | |
753b19cd | 932 | depthxOld = ( para * TMath::Log(energy) + parb ) * vInc.X()/TMath::Abs(vInc.Y()) ; |
933 | depthzOld = ( para * TMath::Log(energy) + parb ) * vInc.Z()/TMath::Abs(vInc.Y()) ; | |
7fb46731 | 934 | } |
adcca1e6 | 935 | else{ |
407d15b3 | 936 | AliError("Cluster with zero energy \n"); |
2c06dc7a | 937 | } |
407d15b3 | 938 | //Apply Real vertex |
939 | phosgeom->GetIncidentVector(fVtx,emc->GetPHOSMod(),x,z,vInc) ; | |
940 | Float_t depthx = 0.; | |
941 | Float_t depthz = 0.; | |
942 | if (energy > 0 && vInc.Y()!=0.) { | |
e5b7b511 | 943 | depthx = ( para * TMath::Log(energy) + parb ) * vInc.X()/TMath::Abs(vInc.Y()) ; |
944 | depthz = ( para * TMath::Log(energy) + parb ) * vInc.Z()/TMath::Abs(vInc.Y()) ; | |
adcca1e6 | 945 | } |
407d15b3 | 946 | |
e5b7b511 | 947 | //Correct for the vertex position and shower depth |
948 | Double_t xd=x+(depthxOld-depthx) ; | |
949 | Double_t zd=z+(depthzOld-depthz) ; | |
950 | TVector3 dir(0,0,0) ; | |
951 | phosgeom->Local2Global(emc->GetPHOSMod(),xd,zd,dir) ; | |
407d15b3 | 952 | |
e5b7b511 | 953 | dir-=fVtx ; |
7fb46731 | 954 | dir.SetMag(1.) ; |
e3817e5f | 955 | |
88cb7938 | 956 | return dir ; |
7b7c1533 | 957 | } |
958 | ||
35adb638 | 959 | //________________________________________________________________________ |
17323043 | 960 | Double_t AliPHOSPIDv1::LandauF(Double_t x, Double_t y, Double_t * par) |
35adb638 | 961 | { |
c947e71a | 962 | //Given the energy x and the parameter y (tof, shower dispersion or cpv-emc distance), |
963 | //this method returns a density probability of this parameter, given by a landau | |
964 | //function whose parameters depend with the energy with a function: a/(x*x)+b/x+b | |
965 | ||
2924941c | 966 | if (x > par[9]) x = par[9]; |
967 | ||
968 | //Double_t cnt = par[1] / (x*x) + par[2] / x + par[0] ; | |
969 | Double_t cnt = par[0] + par[1] * x + par[2] * x * x ; | |
35adb638 | 970 | Double_t mean = par[4] / (x*x) + par[5] / x + par[3] ; |
971 | Double_t sigma = par[7] / (x*x) + par[8] / x + par[6] ; | |
c947e71a | 972 | |
973 | if(TMath::Abs(sigma) > 1.e-10){ | |
acb5beb7 | 974 | return cnt*TMath::Landau(y,mean,sigma); |
35adb638 | 975 | } |
976 | else | |
977 | return 0.; | |
978 | ||
979 | } | |
980 | //________________________________________________________________________ | |
17323043 | 981 | Double_t AliPHOSPIDv1::LandauPol2(Double_t x, Double_t y, Double_t * par) |
35adb638 | 982 | { |
c947e71a | 983 | |
984 | //Given the energy x and the parameter y (tof, shower dispersion or cpv-emc distance), | |
985 | //this method returns a density probability of this parameter, given by a landau | |
986 | //function whose parameters depend with the energy like second order polinomial | |
987 | ||
35adb638 | 988 | Double_t cnt = par[2] * (x*x) + par[1] * x + par[0] ; |
c947e71a | 989 | Double_t mean = par[5] * (x*x) + par[4] * x + par[3] ; |
990 | Double_t sigma = par[8] * (x*x) + par[7] * x + par[6] ; | |
35adb638 | 991 | |
c947e71a | 992 | if(TMath::Abs(sigma) > 1.e-10){ |
acb5beb7 | 993 | return cnt*TMath::Landau(y,mean,sigma); |
35adb638 | 994 | } |
995 | else | |
996 | return 0.; | |
c947e71a | 997 | |
998 | ||
35adb638 | 999 | } |
1000 | // //________________________________________________________________________ | |
1001 | // Double_t AliPHOSPIDv1::ChargedHadronDistProb(Double_t x, Double_t y, Double_t * parg, Double_t * parl) | |
1002 | // { | |
1003 | // Double_t cnt = 0.0 ; | |
1004 | // Double_t mean = 0.0 ; | |
1005 | // Double_t sigma = 0.0 ; | |
1006 | // Double_t arg = 0.0 ; | |
1007 | // if (y < parl[4] / (x*x) + parl[5] / x + parl[3]){ | |
1008 | // cnt = parg[1] / (x*x) + parg[2] / x + parg[0] ; | |
1009 | // mean = parg[4] / (x*x) + parg[5] / x + parg[3] ; | |
1010 | // sigma = parg[7] / (x*x) + parg[8] / x + parg[6] ; | |
1011 | // TF1 * f = new TF1("gaus","gaus",0.,100.); | |
1012 | // f->SetParameters(cnt,mean,sigma); | |
1013 | // arg = f->Eval(y) ; | |
1014 | // } | |
1015 | // else{ | |
1016 | // cnt = parl[1] / (x*x) + parl[2] / x + parl[0] ; | |
1017 | // mean = parl[4] / (x*x) + parl[5] / x + parl[3] ; | |
1018 | // sigma = parl[7] / (x*x) + parl[8] / x + parl[6] ; | |
1019 | // TF1 * f = new TF1("landau","landau",0.,100.); | |
1020 | // f->SetParameters(cnt,mean,sigma); | |
1021 | // arg = f->Eval(y) ; | |
1022 | // } | |
1023 | // // Double_t mean = par[3] + par[4] * x + par[5] * x * x ; | |
1024 | // // Double_t sigma = par[6] + par[7] * x + par[8] * x * x ; | |
1025 | ||
1026 | // //Double_t arg = -(y-mean)*(y-mean)/(2*sigma*sigma) ; | |
1027 | // //return cnt * TMath::Exp(arg) ; | |
1028 | ||
1029 | // return arg; | |
1030 | ||
1031 | // } | |
2cc71c1e | 1032 | //____________________________________________________________________________ |
1033 | void AliPHOSPIDv1::MakePID() | |
1034 | { | |
1035 | // construct the PID weight from a Bayesian Method | |
c947e71a | 1036 | |
304864ab | 1037 | const Int_t kSPECIES = AliPID::kSPECIESN ; |
7fb46731 | 1038 | |
1039 | AliPHOSGetter * gime = AliPHOSGetter::Instance() ; | |
1040 | ||
1041 | Int_t nparticles = gime->RecParticles()->GetEntriesFast() ; | |
1042 | ||
7fb46731 | 1043 | TObjArray * emcRecPoints = gime->EmcRecPoints() ; |
1044 | TObjArray * cpvRecPoints = gime->CpvRecPoints() ; | |
1045 | TClonesArray * trackSegments = gime->TrackSegments() ; | |
1046 | if ( !emcRecPoints || !cpvRecPoints || !trackSegments ) { | |
1047 | AliFatal("RecPoints or TrackSegments not found !") ; | |
1048 | } | |
1049 | TIter next(trackSegments) ; | |
1050 | AliPHOSTrackSegment * ts ; | |
1051 | Int_t index = 0 ; | |
1052 | ||
35adb638 | 1053 | Double_t * stof[kSPECIES] ; |
1054 | Double_t * sdp [kSPECIES] ; | |
1055 | Double_t * scpv[kSPECIES] ; | |
fb7b51ad | 1056 | Double_t * sw [kSPECIES] ; |
35adb638 | 1057 | //Info("MakePID","Begin MakePID"); |
1058 | ||
1059 | for (Int_t i =0; i< kSPECIES; i++){ | |
1060 | stof[i] = new Double_t[nparticles] ; | |
1061 | sdp [i] = new Double_t[nparticles] ; | |
1062 | scpv[i] = new Double_t[nparticles] ; | |
fb7b51ad | 1063 | sw [i] = new Double_t[nparticles] ; |
35adb638 | 1064 | } |
1065 | ||
7fb46731 | 1066 | |
1067 | while ( (ts = (AliPHOSTrackSegment *)next()) ) { | |
1068 | ||
1069 | //cout<<">>>>>> Bayesian Index "<<index<<endl; | |
1070 | ||
1071 | AliPHOSEmcRecPoint * emc = 0 ; | |
1072 | if(ts->GetEmcIndex()>=0) | |
1073 | emc = (AliPHOSEmcRecPoint *) emcRecPoints->At(ts->GetEmcIndex()) ; | |
1074 | ||
407d15b3 | 1075 | // AliPHOSCpvRecPoint * cpv = 0 ; |
1076 | // if(ts->GetCpvIndex()>=0) | |
1077 | // cpv = (AliPHOSCpvRecPoint *) cpvRecPoints->At(ts->GetCpvIndex()) ; | |
1078 | // | |
1079 | //// Int_t track = 0 ; | |
1080 | //// track = ts->GetTrackIndex() ; //TPC tracks ? | |
cc1fe362 | 1081 | |
7fb46731 | 1082 | if (!emc) { |
1083 | AliFatal(Form("-> emc(%d) = %d", ts->GetEmcIndex(), emc )) ; | |
1084 | } | |
c947e71a | 1085 | |
2924941c | 1086 | |
7fb46731 | 1087 | // ############Tof############################# |
c947e71a | 1088 | |
7fb46731 | 1089 | // Info("MakePID", "TOF"); |
fb7b51ad | 1090 | Float_t en = emc->GetEnergy(); |
7fb46731 | 1091 | Double_t time = emc->GetTime() ; |
1092 | // cout<<">>>>>>>Energy "<<en<<"Time "<<time<<endl; | |
fb7b51ad | 1093 | |
35adb638 | 1094 | // now get the signals probability |
1095 | // s(pid) in the Bayesian formulation | |
cc1fe362 | 1096 | |
304864ab | 1097 | stof[AliPID::kPhoton][index] = 1.; |
1098 | stof[AliPID::kElectron][index] = 1.; | |
304864ab | 1099 | stof[AliPID::kEleCon][index] = 1.; |
7fb46731 | 1100 | //We assing the same prob to charged hadrons, sum is 1 |
1101 | stof[AliPID::kPion][index] = 1./3.; | |
1102 | stof[AliPID::kKaon][index] = 1./3.; | |
1103 | stof[AliPID::kProton][index] = 1./3.; | |
1104 | //We assing the same prob to neutral hadrons, sum is 1 | |
1105 | stof[AliPID::kNeutron][index] = 1./2.; | |
1106 | stof[AliPID::kKaon0][index] = 1./2.; | |
304864ab | 1107 | stof[AliPID::kMuon][index] = 1.; |
2924941c | 1108 | |
1109 | if(en < fTOFEnThreshold) { | |
7fb46731 | 1110 | |
1111 | Double_t pTofPion = fTFpiong ->Eval(time) ; //gaus distribution | |
1112 | Double_t pTofKaon = 0; | |
1113 | ||
35adb638 | 1114 | if(time < fTkaonl[1]) |
fb7b51ad | 1115 | pTofKaon = fTFkaong ->Eval(time) ; //gaus distribution |
35adb638 | 1116 | else |
fb7b51ad | 1117 | pTofKaon = fTFkaonl ->Eval(time) ; //landau distribution |
7fb46731 | 1118 | |
1119 | Double_t pTofNucleon = 0; | |
1120 | ||
35adb638 | 1121 | if(time < fThhadronl[1]) |
7fb46731 | 1122 | pTofNucleon = fTFhhadrong ->Eval(time) ; //gaus distribution |
35adb638 | 1123 | else |
7fb46731 | 1124 | pTofNucleon = fTFhhadronl ->Eval(time) ; //landau distribution |
7fb46731 | 1125 | //We assing the same prob to neutral hadrons, sum is the average prob |
fb7b51ad | 1126 | Double_t pTofNeHadron = (pTofKaon + pTofNucleon)/2. ; |
1127 | //We assing the same prob to charged hadrons, sum is the average prob | |
1128 | Double_t pTofChHadron = (pTofPion + pTofKaon + pTofNucleon)/3. ; | |
7fb46731 | 1129 | |
fb7b51ad | 1130 | stof[AliPID::kPhoton][index] = fTFphoton ->Eval(time) ; |
1131 | //gaus distribution | |
1132 | stof[AliPID::kEleCon][index] = stof[AliPID::kPhoton][index] ; | |
1133 | //a conversion electron has the photon ToF | |
304864ab | 1134 | stof[AliPID::kMuon][index] = stof[AliPID::kPhoton][index] ; |
7fb46731 | 1135 | |
fb7b51ad | 1136 | stof[AliPID::kElectron][index] = pTofPion ; |
7fb46731 | 1137 | |
1138 | stof[AliPID::kPion][index] = pTofChHadron ; | |
1139 | stof[AliPID::kKaon][index] = pTofChHadron ; | |
1140 | stof[AliPID::kProton][index] = pTofChHadron ; | |
1141 | ||
1142 | stof[AliPID::kKaon0][index] = pTofNeHadron ; | |
1143 | stof[AliPID::kNeutron][index] = pTofNeHadron ; | |
cc1fe362 | 1144 | } |
c947e71a | 1145 | |
1146 | // Info("MakePID", "Dispersion"); | |
cc1fe362 | 1147 | |
7fb46731 | 1148 | // ###########Shower shape: Dispersion#################### |
35adb638 | 1149 | Float_t dispersion = emc->GetDispersion(); |
407d15b3 | 1150 | //DP: Correct for non-perpendicular incidence |
1151 | //DP: still to be done | |
1152 | ||
35adb638 | 1153 | //dispersion is not well defined if the cluster is only in few crystals |
cc1fe362 | 1154 | |
304864ab | 1155 | sdp[AliPID::kPhoton][index] = 1. ; |
1156 | sdp[AliPID::kElectron][index] = 1. ; | |
1157 | sdp[AliPID::kPion][index] = 1. ; | |
1158 | sdp[AliPID::kKaon][index] = 1. ; | |
1159 | sdp[AliPID::kProton][index] = 1. ; | |
1160 | sdp[AliPID::kNeutron][index] = 1. ; | |
1161 | sdp[AliPID::kEleCon][index] = 1. ; | |
1162 | sdp[AliPID::kKaon0][index] = 1. ; | |
1163 | sdp[AliPID::kMuon][index] = 1. ; | |
c947e71a | 1164 | |
2924941c | 1165 | if(en > fDispEnThreshold && emc->GetMultiplicity() > fDispMultThreshold){ |
7fb46731 | 1166 | sdp[AliPID::kPhoton][index] = GausF(en , dispersion, fDphoton) ; |
304864ab | 1167 | sdp[AliPID::kElectron][index] = sdp[AliPID::kPhoton][index] ; |
1168 | sdp[AliPID::kPion][index] = LandauF(en , dispersion, fDhadron ) ; | |
1169 | sdp[AliPID::kKaon][index] = sdp[AliPID::kPion][index] ; | |
1170 | sdp[AliPID::kProton][index] = sdp[AliPID::kPion][index] ; | |
1171 | sdp[AliPID::kNeutron][index] = sdp[AliPID::kPion][index] ; | |
1172 | sdp[AliPID::kEleCon][index] = sdp[AliPID::kPhoton][index]; | |
1173 | sdp[AliPID::kKaon0][index] = sdp[AliPID::kPion][index] ; | |
fb7b51ad | 1174 | sdp[AliPID::kMuon][index] = fDFmuon ->Eval(dispersion) ; |
1175 | //landau distribution | |
35adb638 | 1176 | } |
cc1fe362 | 1177 | |
7fb46731 | 1178 | // Info("MakePID","multiplicity %d, dispersion %f", emc->GetMultiplicity(), dispersion); |
1179 | // Info("MakePID","ss: photon %f, hadron %f ", sdp[AliPID::kPhoton][index], sdp[AliPID::kPion][index]); | |
c947e71a | 1180 | // cout<<">>>>>multiplicity "<<emc->GetMultiplicity()<<", dispersion "<< dispersion<<endl ; |
304864ab | 1181 | // cout<<"<<<<<ss: photon "<<sdp[AliPID::kPhoton][index]<<", hadron "<<sdp[AliPID::kPion][index]<<endl; |
c947e71a | 1182 | |
7fb46731 | 1183 | //########## CPV-EMC Distance####################### |
1184 | // Info("MakePID", "Distance"); | |
fb7b51ad | 1185 | |
26aa7e4a | 1186 | Float_t x = TMath::Abs(ts->GetCpvDistance("X")) ; |
1187 | Float_t z = ts->GetCpvDistance("Z") ; | |
fb7b51ad | 1188 | |
7fb46731 | 1189 | Double_t pcpv = 0 ; |
c947e71a | 1190 | Double_t pcpvneutral = 0. ; |
7fb46731 | 1191 | |
1192 | Double_t elprobx = GausF(en , x, fXelectron) ; | |
1193 | Double_t elprobz = GausF(en , z, fZelectron) ; | |
1194 | Double_t chprobx = GausF(en , x, fXcharged) ; | |
1195 | Double_t chprobz = GausF(en , z, fZcharged) ; | |
c947e71a | 1196 | Double_t pcpvelectron = elprobx * elprobz; |
1197 | Double_t pcpvcharged = chprobx * chprobz; | |
7fb46731 | 1198 | |
1199 | // cout<<">>>>energy "<<en<<endl; | |
c947e71a | 1200 | // cout<<">>>>electron : x "<<x<<" xprob "<<elprobx<<" z "<<z<<" zprob "<<elprobz<<endl; |
1201 | // cout<<">>>>hadron : x "<<x<<" xprob "<<chprobx<<" z "<<z<<" zprob "<<chprobz<<endl; | |
1202 | // cout<<">>>>electron : px*pz "<<pcpvelectron <<" hadron: px*pz "<<pcpvcharged<<endl; | |
1203 | ||
2924941c | 1204 | // Is neutral or charged? |
35adb638 | 1205 | if(pcpvelectron >= pcpvcharged) |
1206 | pcpv = pcpvelectron ; | |
1207 | else | |
1208 | pcpv = pcpvcharged ; | |
1209 | ||
fb7b51ad | 1210 | if(pcpv < fChargedNeutralThreshold) |
35adb638 | 1211 | { |
1212 | pcpvneutral = 1. ; | |
1213 | pcpvcharged = 0. ; | |
1214 | pcpvelectron = 0. ; | |
1215 | } | |
c947e71a | 1216 | // else |
1217 | // cout<<">>>>>>>>>>>CHARGED>>>>>>>>>>>"<<endl; | |
35adb638 | 1218 | |
304864ab | 1219 | scpv[AliPID::kPion][index] = pcpvcharged ; |
1220 | scpv[AliPID::kKaon][index] = pcpvcharged ; | |
1221 | scpv[AliPID::kProton][index] = pcpvcharged ; | |
7fb46731 | 1222 | |
1223 | scpv[AliPID::kMuon][index] = pcpvelectron ; | |
304864ab | 1224 | scpv[AliPID::kElectron][index] = pcpvelectron ; |
304864ab | 1225 | scpv[AliPID::kEleCon][index] = pcpvelectron ; |
7fb46731 | 1226 | |
1227 | scpv[AliPID::kPhoton][index] = pcpvneutral ; | |
1228 | scpv[AliPID::kNeutron][index] = pcpvneutral ; | |
304864ab | 1229 | scpv[AliPID::kKaon0][index] = pcpvneutral ; |
7fb46731 | 1230 | |
35adb638 | 1231 | |
1232 | // Info("MakePID", "CPV passed"); | |
c947e71a | 1233 | |
7fb46731 | 1234 | //############## Pi0 ############################# |
304864ab | 1235 | stof[AliPID::kPi0][index] = 0. ; |
1236 | scpv[AliPID::kPi0][index] = 0. ; | |
1237 | sdp [AliPID::kPi0][index] = 0. ; | |
c947e71a | 1238 | |
35adb638 | 1239 | if(en > 30.){ |
1240 | // pi0 are detected via decay photon | |
2924941c | 1241 | stof[AliPID::kPi0][index] = stof[AliPID::kPhoton][index]; |
304864ab | 1242 | scpv[AliPID::kPi0][index] = pcpvneutral ; |
2924941c | 1243 | if(emc->GetMultiplicity() > fDispMultThreshold) |
1244 | sdp [AliPID::kPi0][index] = GausF(en , dispersion, fDpi0) ; | |
1245 | //sdp [AliPID::kPi0][index] = GausPol2(en , dispersion, fDpi0) ; | |
1246 | // cout<<"E = "<<en<<" GeV; disp = "<<dispersion<<"; mult = " | |
1247 | // <<emc->GetMultiplicity()<<endl; | |
1248 | // cout<<"PDF: photon = "<<sdp [AliPID::kPhoton][index]<<"; pi0 = " | |
1249 | // <<sdp [AliPID::kPi0][index]<<endl; | |
35adb638 | 1250 | } |
1251 | ||
2924941c | 1252 | |
7fb46731 | 1253 | |
2924941c | 1254 | |
7fb46731 | 1255 | //############## muon ############################# |
1256 | ||
35adb638 | 1257 | if(en > 0.5){ |
1258 | //Muons deposit few energy | |
304864ab | 1259 | scpv[AliPID::kMuon][index] = 0 ; |
1260 | stof[AliPID::kMuon][index] = 0 ; | |
1261 | sdp [AliPID::kMuon][index] = 0 ; | |
c947e71a | 1262 | } |
1263 | ||
fb7b51ad | 1264 | //Weight to apply to hadrons due to energy reconstruction |
1265 | ||
1266 | Float_t weight = fERecWeight ->Eval(en) ; | |
1267 | ||
1268 | sw[AliPID::kPhoton][index] = 1. ; | |
1269 | sw[AliPID::kElectron][index] = 1. ; | |
1270 | sw[AliPID::kPion][index] = weight ; | |
1271 | sw[AliPID::kKaon][index] = weight ; | |
1272 | sw[AliPID::kProton][index] = weight ; | |
1273 | sw[AliPID::kNeutron][index] = weight ; | |
1274 | sw[AliPID::kEleCon][index] = 1. ; | |
1275 | sw[AliPID::kKaon0][index] = weight ; | |
1276 | sw[AliPID::kMuon][index] = weight ; | |
1277 | sw[AliPID::kPi0][index] = 1. ; | |
1278 | ||
7fb46731 | 1279 | // if(en > 0.5){ |
1280 | // cout<<"######################################################"<<endl; | |
1281 | // //cout<<"MakePID: energy "<<en<<", tof "<<time<<", distance "<<distance<<", dispersion "<<dispersion<<endl ; | |
1282 | // cout<<"MakePID: energy "<<en<<", tof "<<time<<", dispersion "<<dispersion<<", x "<<x<<", z "<<z<<endl ; | |
1283 | // cout<<">>>>>multiplicity "<<emc->GetMultiplicity()<<endl; | |
1284 | // cout<<">>>>electron : xprob "<<elprobx<<" zprob "<<elprobz<<endl; | |
1285 | // cout<<">>>>hadron : xprob "<<chprobx<<" zprob "<<chprobz<<endl; | |
1286 | // cout<<">>>>electron : px*pz "<<pcpvelectron <<" hadron: px*pz "<<pcpvcharged<<endl; | |
c947e71a | 1287 | |
2924941c | 1288 | // cout<<"Photon , pid "<< fInitPID[AliPID::kPhoton]<<" tof "<<stof[AliPID::kPhoton][index] |
1289 | // <<", cpv "<<scpv[AliPID::kPhoton][index]<<", ss "<<sdp[AliPID::kPhoton][index]<<endl; | |
7fb46731 | 1290 | // cout<<"EleCon , pid "<< fInitPID[AliPID::kEleCon]<<", tof "<<stof[AliPID::kEleCon][index] |
1291 | // <<", cpv "<<scpv[AliPID::kEleCon][index]<<" ss "<<sdp[AliPID::kEleCon][index]<<endl; | |
1292 | // cout<<"Electron , pid "<< fInitPID[AliPID::kElectron]<<", tof "<<stof[AliPID::kElectron][index] | |
1293 | // <<", cpv "<<scpv[AliPID::kElectron][index]<<" ss "<<sdp[AliPID::kElectron][index]<<endl; | |
1294 | // cout<<"Muon , pid "<< fInitPID[AliPID::kMuon]<<", tof "<<stof[AliPID::kMuon][index] | |
1295 | // <<", cpv "<<scpv[AliPID::kMuon][index]<<" ss "<<sdp[AliPID::kMuon][index]<<endl; | |
2924941c | 1296 | // cout<<"Pi0 , pid "<< fInitPID[AliPID::kPi0]<<", tof "<<stof[AliPID::kPi0][index] |
1297 | // <<", cpv "<<scpv[AliPID::kPi0][index]<<" ss "<<sdp[AliPID::kPi0][index]<<endl; | |
7fb46731 | 1298 | // cout<<"Pion , pid "<< fInitPID[AliPID::kPion]<<", tof "<<stof[AliPID::kPion][index] |
1299 | // <<", cpv "<<scpv[AliPID::kPion][index]<<" ss "<<sdp[AliPID::kPion][index]<<endl; | |
1300 | // cout<<"Kaon0 , pid "<< fInitPID[AliPID::kKaon0]<<", tof "<<stof[AliPID::kKaon0][index] | |
1301 | // <<", cpv "<<scpv[AliPID::kKaon0][index]<<" ss "<<sdp[AliPID::kKaon0][index]<<endl; | |
1302 | // cout<<"Kaon , pid "<< fInitPID[AliPID::kKaon]<<", tof "<<stof[AliPID::kKaon][index] | |
1303 | // <<", cpv "<<scpv[AliPID::kKaon][index]<<" ss "<<sdp[AliPID::kKaon][index]<<endl; | |
1304 | // cout<<"Neutron , pid "<< fInitPID[AliPID::kNeutron]<<", tof "<<stof[AliPID::kNeutron][index] | |
1305 | // <<", cpv "<<scpv[AliPID::kNeutron][index]<<" ss "<<sdp[AliPID::kNeutron][index]<<endl; | |
1306 | // cout<<"Proton , pid "<< fInitPID[AliPID::kProton]<<", tof "<<stof[AliPID::kProton][index] | |
1307 | // <<", cpv "<<scpv[AliPID::kProton][index]<<" ss "<<sdp[AliPID::kProton][index]<<endl; | |
1308 | // cout<<"######################################################"<<endl; | |
1309 | // } | |
1310 | index++; | |
cc1fe362 | 1311 | } |
35adb638 | 1312 | |
1313 | //for (index = 0 ; index < kSPECIES ; index++) | |
1314 | // pid[index] /= nparticles ; | |
1315 | ||
7fb46731 | 1316 | |
35adb638 | 1317 | // Info("MakePID", "Total Probability calculation"); |
1318 | ||
cc1fe362 | 1319 | for(index = 0 ; index < nparticles ; index ++) { |
2924941c | 1320 | |
1321 | AliPHOSRecParticle * recpar = gime->RecParticle(index) ; | |
1322 | ||
1323 | //Conversion electron? | |
1324 | ||
1325 | if(recpar->IsEleCon()){ | |
1326 | fInitPID[AliPID::kEleCon] = 1. ; | |
1327 | fInitPID[AliPID::kPhoton] = 0. ; | |
1328 | fInitPID[AliPID::kElectron] = 0. ; | |
1329 | } | |
1330 | else{ | |
1331 | fInitPID[AliPID::kEleCon] = 0. ; | |
1332 | fInitPID[AliPID::kPhoton] = 1. ; | |
1333 | fInitPID[AliPID::kElectron] = 1. ; | |
1334 | } | |
1335 | // fInitPID[AliPID::kEleCon] = 0. ; | |
1336 | ||
1337 | ||
cc1fe362 | 1338 | // calculates the Bayesian weight |
7fb46731 | 1339 | |
cc1fe362 | 1340 | Int_t jndex ; |
1341 | Double_t wn = 0.0 ; | |
1342 | for (jndex = 0 ; jndex < kSPECIES ; jndex++) | |
2924941c | 1343 | wn += stof[jndex][index] * sdp[jndex][index] * scpv[jndex][index] * |
1344 | sw[jndex][index] * fInitPID[jndex] ; | |
1345 | ||
7fb46731 | 1346 | // cout<<"*************wn "<<wn<<endl; |
e74ea0e9 | 1347 | if (TMath::Abs(wn)>0) |
1348 | for (jndex = 0 ; jndex < kSPECIES ; jndex++) { | |
35adb638 | 1349 | //cout<<"jndex "<<jndex<<" wn "<<wn<<" SetPID * wn" |
1350 | //<<stof[jndex][index] * sdp[jndex][index] * pid[jndex] << endl; | |
1351 | //cout<<" tof "<<stof[jndex][index] << " disp " <<sdp[jndex][index] << " pid "<< fInitPID[jndex] << endl; | |
2924941c | 1352 | // if(jndex == AliPID::kPi0 || jndex == AliPID::kPhoton){ |
1353 | // cout<<"Particle "<<jndex<<" final prob * wn " | |
1354 | // <<stof[jndex][index] * sdp[jndex][index] * scpv[jndex][index] * | |
1355 | // fInitPID[jndex] <<" wn "<< wn<<endl; | |
1356 | // cout<<"pid "<< fInitPID[jndex]<<", tof "<<stof[jndex][index] | |
1357 | // <<", cpv "<<scpv[jndex][index]<<" ss "<<sdp[jndex][index]<<endl; | |
1358 | // } | |
1359 | recpar->SetPID(jndex, stof[jndex][index] * sdp[jndex][index] * | |
1360 | sw[jndex][index] * scpv[jndex][index] * | |
1361 | fInitPID[jndex] / wn) ; | |
e74ea0e9 | 1362 | } |
2cc71c1e | 1363 | } |
35adb638 | 1364 | // Info("MakePID", "Delete"); |
1365 | ||
2924941c | 1366 | for (Int_t i =0; i< kSPECIES; i++){ |
1367 | delete [] stof[i]; | |
1368 | delete [] sdp [i]; | |
1369 | delete [] scpv[i]; | |
1370 | delete [] sw [i]; | |
1371 | } | |
35adb638 | 1372 | // Info("MakePID","End MakePID"); |
2cc71c1e | 1373 | } |
1374 | ||
7acf6008 | 1375 | //____________________________________________________________________________ |
e3817e5f | 1376 | void AliPHOSPIDv1::MakeRecParticles() |
1377 | { | |
b2a60966 | 1378 | // Makes a RecParticle out of a TrackSegment |
148b2bba | 1379 | |
88cb7938 | 1380 | AliPHOSGetter * gime = AliPHOSGetter::Instance() ; |
fbf811ec | 1381 | TObjArray * emcRecPoints = gime->EmcRecPoints() ; |
1382 | TObjArray * cpvRecPoints = gime->CpvRecPoints() ; | |
1383 | TClonesArray * trackSegments = gime->TrackSegments() ; | |
148b2bba | 1384 | if ( !emcRecPoints || !cpvRecPoints || !trackSegments ) { |
351dd634 | 1385 | AliFatal("RecPoints or TrackSegments not found !") ; |
148b2bba | 1386 | } |
fbf811ec | 1387 | TClonesArray * recParticles = gime->RecParticles() ; |
01a599c9 | 1388 | recParticles->Clear(); |
148b2bba | 1389 | |
7b7c1533 | 1390 | TIter next(trackSegments) ; |
7acf6008 | 1391 | AliPHOSTrackSegment * ts ; |
6ad0bfa0 | 1392 | Int_t index = 0 ; |
09fc14a0 | 1393 | AliPHOSRecParticle * rp ; |
7acf6008 | 1394 | while ( (ts = (AliPHOSTrackSegment *)next()) ) { |
7fb46731 | 1395 | // cout<<">>>>>>>>>>>>>>>PCA Index "<<index<<endl; |
7b7c1533 | 1396 | new( (*recParticles)[index] ) AliPHOSRecParticle() ; |
1397 | rp = (AliPHOSRecParticle *)recParticles->At(index) ; | |
f0a4c9e9 | 1398 | rp->SetTrackSegment(index) ; |
9688c1dd | 1399 | rp->SetIndexInList(index) ; |
148b2bba | 1400 | |
7acf6008 | 1401 | AliPHOSEmcRecPoint * emc = 0 ; |
1402 | if(ts->GetEmcIndex()>=0) | |
7b7c1533 | 1403 | emc = (AliPHOSEmcRecPoint *) emcRecPoints->At(ts->GetEmcIndex()) ; |
fad3e5b9 | 1404 | |
8d4608b5 | 1405 | AliPHOSCpvRecPoint * cpv = 0 ; |
7acf6008 | 1406 | if(ts->GetCpvIndex()>=0) |
8d4608b5 | 1407 | cpv = (AliPHOSCpvRecPoint *) cpvRecPoints->At(ts->GetCpvIndex()) ; |
fad3e5b9 | 1408 | |
bd76890a | 1409 | Int_t track = 0 ; |
1410 | track = ts->GetTrackIndex() ; | |
1411 | ||
148b2bba | 1412 | // Now set type (reconstructed) of the particle |
1413 | ||
1414 | // Choose the cluster energy range | |
9fa5f1d0 | 1415 | |
fbf811ec | 1416 | if (!emc) { |
351dd634 | 1417 | AliFatal(Form("-> emc(%d) = %d", ts->GetEmcIndex(), emc )) ; |
fbf811ec | 1418 | } |
50739f15 | 1419 | |
e3817e5f | 1420 | Float_t e = emc->GetEnergy() ; |
bc0c084c | 1421 | |
6f969528 | 1422 | Float_t lambda[2] ; |
1423 | emc->GetElipsAxis(lambda) ; | |
407d15b3 | 1424 | |
50739f15 | 1425 | if((lambda[0]>0.01) && (lambda[1]>0.01)){ |
1426 | // Looking PCA. Define and calculate the data (X), | |
bc0c084c | 1427 | // introduce in the function X2P that gives the components (P). |
1428 | ||
c947e71a | 1429 | Float_t spher = 0. ; |
1430 | Float_t emaxdtotal = 0. ; | |
50739f15 | 1431 | |
bc0c084c | 1432 | if((lambda[0]+lambda[1])!=0) |
e4df4b30 | 1433 | spher=TMath::Abs(lambda[0]-lambda[1])/(lambda[0]+lambda[1]); |
50739f15 | 1434 | |
c947e71a | 1435 | emaxdtotal=emc->GetMaximalEnergy()/emc->GetEnergy(); |
50739f15 | 1436 | |
1437 | fX[0] = lambda[0] ; | |
1438 | fX[1] = lambda[1] ; | |
1439 | fX[2] = emc->GetDispersion() ; | |
c947e71a | 1440 | fX[3] = spher ; |
50739f15 | 1441 | fX[4] = emc->GetMultiplicity() ; |
c947e71a | 1442 | fX[5] = emaxdtotal ; |
50739f15 | 1443 | fX[6] = emc->GetCoreEnergy() ; |
1444 | ||
e3817e5f | 1445 | fPrincipalPhoton->X2P(fX,fPPhoton); |
1446 | fPrincipalPi0 ->X2P(fX,fPPi0); | |
1f0e7ccd | 1447 | |
50739f15 | 1448 | } |
1449 | else{ | |
e3817e5f | 1450 | fPPhoton[0]=-100.0; //We do not accept clusters with |
1451 | fPPhoton[1]=-100.0; //one cell as a photon-like | |
1452 | fPPi0[0] =-100.0; | |
1453 | fPPi0[1] =-100.0; | |
50739f15 | 1454 | } |
1455 | ||
2cc71c1e | 1456 | Float_t time = emc->GetTime() ; |
1457 | rp->SetTof(time) ; | |
9fa5f1d0 | 1458 | |
bc0c084c | 1459 | // Loop of Efficiency-Purity (the 3 points of purity or efficiency |
1460 | // are taken into account to set the particle identification) | |
e3817e5f | 1461 | for(Int_t effPur = 0; effPur < 3 ; effPur++){ |
50739f15 | 1462 | |
bc0c084c | 1463 | // Looking at the CPV detector. If RCPV greater than CpvEmcDistance, |
1464 | // 1st,2nd or 3rd bit (depending on the efficiency-purity point ) | |
1465 | // is set to 1 | |
26aa7e4a | 1466 | if(GetCPVBit(ts, effPur,e) == 1 ){ |
e3817e5f | 1467 | rp->SetPIDBit(effPur) ; |
35adb638 | 1468 | //cout<<"CPV bit "<<effPur<<endl; |
1469 | } | |
50739f15 | 1470 | // Looking the TOF. If TOF smaller than gate, 4th, 5th or 6th |
1471 | // bit (depending on the efficiency-purity point )is set to 1 | |
2cc71c1e | 1472 | if(time< (*fParameters)(3,effPur)) |
e3817e5f | 1473 | rp->SetPIDBit(effPur+3) ; |
2cc71c1e | 1474 | |
e3817e5f | 1475 | //Photon PCA |
50739f15 | 1476 | //If we are inside the ellipse, 7th, 8th or 9th |
1477 | // bit (depending on the efficiency-purity point )is set to 1 | |
e3817e5f | 1478 | if(GetPrincipalBit("photon",fPPhoton,effPur,e) == 1) |
1479 | rp->SetPIDBit(effPur+6) ; | |
1f0e7ccd | 1480 | |
e3817e5f | 1481 | //Pi0 PCA |
1f0e7ccd | 1482 | //If we are inside the ellipse, 10th, 11th or 12th |
1483 | // bit (depending on the efficiency-purity point )is set to 1 | |
e3817e5f | 1484 | if(GetPrincipalBit("pi0" ,fPPi0 ,effPur,e) == 1) |
1485 | rp->SetPIDBit(effPur+9) ; | |
f0a4c9e9 | 1486 | } |
e3817e5f | 1487 | if(GetHardPhotonBit(emc)) |
1488 | rp->SetPIDBit(12) ; | |
1489 | if(GetHardPi0Bit (emc)) | |
1490 | rp->SetPIDBit(13) ; | |
1f0e7ccd | 1491 | |
bd76890a | 1492 | if(track >= 0) |
1493 | rp->SetPIDBit(14) ; | |
1494 | ||
9fa5f1d0 | 1495 | //Set momentum, energy and other parameters |
50739f15 | 1496 | Float_t encal = GetCalibratedEnergy(e); |
9fa5f1d0 | 1497 | TVector3 dir = GetMomentumDirection(emc,cpv) ; |
1498 | dir.SetMag(encal) ; | |
1499 | rp->SetMomentum(dir.X(),dir.Y(),dir.Z(),encal) ; | |
1500 | rp->SetCalcMass(0); | |
e0ed2e49 | 1501 | rp->Name(); //If photon sets the particle pdg name to gamma |
407d15b3 | 1502 | rp->SetProductionVertex(fVtx.X(),fVtx.Y(),fVtx.Z(),0); |
e747b8da | 1503 | rp->SetFirstMother(-1); |
1504 | rp->SetLastMother(-1); | |
1505 | rp->SetFirstDaughter(-1); | |
1506 | rp->SetLastDaughter(-1); | |
1507 | rp->SetPolarisation(0,0,0); | |
d956e9b7 | 1508 | //Set the position in global coordinate system from the RecPoint |
1509 | AliPHOSGeometry * geom = gime->PHOSGeometry() ; | |
1510 | AliPHOSTrackSegment * ts = gime->TrackSegment(rp->GetPHOSTSIndex()) ; | |
1511 | AliPHOSEmcRecPoint * erp = gime->EmcRecPoint(ts->GetEmcIndex()) ; | |
1512 | TVector3 pos ; | |
9ee9f78d | 1513 | geom->GetGlobalPHOS(erp, pos) ; |
d956e9b7 | 1514 | rp->SetPos(pos); |
6ad0bfa0 | 1515 | index++ ; |
1516 | } | |
6ad0bfa0 | 1517 | } |
e3817e5f | 1518 | |
09fc14a0 | 1519 | //____________________________________________________________________________ |
702ab87e | 1520 | void AliPHOSPIDv1::Print(const Option_t *) const |
09fc14a0 | 1521 | { |
b2a60966 | 1522 | // Print the parameters used for the particle type identification |
bc0c084c | 1523 | |
351dd634 | 1524 | AliInfo("=============== AliPHOSPIDv1 ================") ; |
88cb7938 | 1525 | printf("Making PID\n") ; |
1526 | printf(" Pricipal analysis file from 0.5 to 100 %s\n", fFileNamePrincipalPhoton.Data() ) ; | |
1527 | printf(" Name of parameters file %s\n", fFileNameParameters.Data() ) ; | |
1528 | printf(" Matrix of Parameters: 14x4\n") ; | |
1529 | printf(" Energy Calibration 1x3 [3 parametres to calibrate energy: A + B* E + C * E^2]\n") ; | |
1530 | printf(" RCPV 2x3 rows x and z, columns function cut parameters\n") ; | |
1531 | printf(" TOF 1x3 [High Eff-Low Pur,Medium Eff-Pur, Low Eff-High Pur]\n") ; | |
1532 | printf(" PCA 5x4 [5 ellipse parametres and 4 parametres to calculate them: A/Sqrt(E) + B* E + C * E^2 + D]\n") ; | |
407d15b3 | 1533 | printf(" Pi0 PCA 5x3 [5 ellipse parametres and 3 parametres to calculate them: A + B* E + C * E^2]\n") ; |
50739f15 | 1534 | fParameters->Print() ; |
09fc14a0 | 1535 | } |
1536 | ||
a496c46c | 1537 | |
69183710 | 1538 | |
7acf6008 | 1539 | //____________________________________________________________________________ |
a4e98857 | 1540 | void AliPHOSPIDv1::PrintRecParticles(Option_t * option) |
1541 | { | |
dd5c4038 | 1542 | // Print table of reconstructed particles |
1543 | ||
88cb7938 | 1544 | AliPHOSGetter *gime = AliPHOSGetter::Instance() ; |
bf8f1fbd | 1545 | |
88cb7938 | 1546 | TClonesArray * recParticles = gime->RecParticles() ; |
21cd0c07 | 1547 | |
1548 | TString message ; | |
3bf72d32 | 1549 | message = "\nevent " ; |
999f9a8f | 1550 | message += gime->EventNumber(); |
3bf72d32 | 1551 | message += " found " ; |
1552 | message += recParticles->GetEntriesFast(); | |
1553 | message += " RecParticles\n" ; | |
1554 | ||
7acf6008 | 1555 | if(strstr(option,"all")) { // printing found TS |
3bf72d32 | 1556 | message += "\n PARTICLE Index \n" ; |
7acf6008 | 1557 | |
1558 | Int_t index ; | |
7b7c1533 | 1559 | for (index = 0 ; index < recParticles->GetEntries() ; index++) { |
21cd0c07 | 1560 | AliPHOSRecParticle * rp = (AliPHOSRecParticle * ) recParticles->At(index) ; |
3bf72d32 | 1561 | message += "\n" ; |
1562 | message += rp->Name().Data() ; | |
1563 | message += " " ; | |
1564 | message += rp->GetIndexInList() ; | |
1565 | message += " " ; | |
1566 | message += rp->GetType() ; | |
7acf6008 | 1567 | } |
3bf72d32 | 1568 | } |
351dd634 | 1569 | AliInfo(message.Data() ) ; |
69183710 | 1570 | } |
88cb7938 | 1571 | |
1572 | //____________________________________________________________________________ | |
1573 | void AliPHOSPIDv1::SetParameters() | |
1574 | { | |
1575 | // PCA : To do the Principal Components Analysis it is necessary | |
1576 | // the Principal file, which is opened here | |
1577 | fX = new double[7]; // Data for the PCA | |
1578 | fPPhoton = new double[7]; // Eigenvalues of the PCA | |
1579 | fPPi0 = new double[7]; // Eigenvalues of the Pi0 PCA | |
1580 | ||
1581 | // Read photon principals from the photon file | |
1582 | ||
1583 | fFileNamePrincipalPhoton = "$ALICE_ROOT/PHOS/PCA8pa15_0.5-100.root" ; | |
1584 | TFile f( fFileNamePrincipalPhoton.Data(), "read" ) ; | |
1585 | fPrincipalPhoton = dynamic_cast<TPrincipal*> (f.Get("principal")) ; | |
1586 | f.Close() ; | |
1587 | ||
1588 | // Read pi0 principals from the pi0 file | |
1589 | ||
1590 | fFileNamePrincipalPi0 = "$ALICE_ROOT/PHOS/PCA_pi0_40-120.root" ; | |
1591 | TFile fPi0( fFileNamePrincipalPi0.Data(), "read" ) ; | |
1592 | fPrincipalPi0 = dynamic_cast<TPrincipal*> (fPi0.Get("principal")) ; | |
1593 | fPi0.Close() ; | |
1594 | ||
1595 | // Open parameters file and initialization of the Parameters matrix. | |
1596 | // In the File Parameters.dat are all the parameters. These are introduced | |
1597 | // in a matrix of 16x4 | |
1598 | // | |
1599 | // All the parameters defined in this file are, in order of row: | |
1600 | // line 0 : calibration | |
1601 | // lines 1,2 : CPV rectangular cat for X and Z | |
1602 | // line 3 : TOF cut | |
1603 | // lines 4-8 : parameters to calculate photon PCA ellipse | |
1604 | // lines 9-13: parameters to calculate pi0 PCA ellipse | |
1605 | // lines 14-15: parameters to calculate border for high-pt photons and pi0 | |
1606 | ||
1607 | fFileNameParameters = gSystem->ExpandPathName("$ALICE_ROOT/PHOS/Parameters.dat"); | |
e8d02863 | 1608 | fParameters = new TMatrixF(16,4) ; |
c947e71a | 1609 | const Int_t kMaxLeng=255; |
1610 | char string[kMaxLeng]; | |
88cb7938 | 1611 | |
1612 | // Open a text file with PID parameters | |
1613 | FILE *fd = fopen(fFileNameParameters.Data(),"r"); | |
1614 | if (!fd) | |
351dd634 | 1615 | AliFatal(Form("File %s with a PID parameters cannot be opened\n", |
1616 | fFileNameParameters.Data())); | |
88cb7938 | 1617 | |
1618 | Int_t i=0; | |
1619 | // Read parameter file line-by-line and skip empty line and comments | |
c947e71a | 1620 | while (fgets(string,kMaxLeng,fd) != NULL) { |
88cb7938 | 1621 | if (string[0] == '\n' ) continue; |
1622 | if (string[0] == '!' ) continue; | |
1623 | sscanf(string, "%f %f %f %f", | |
1624 | &(*fParameters)(i,0), &(*fParameters)(i,1), | |
1625 | &(*fParameters)(i,2), &(*fParameters)(i,3)); | |
1626 | i++; | |
351dd634 | 1627 | AliDebug(1, Form("SetParameters", "line %d: %s",i,string)); |
88cb7938 | 1628 | } |
1629 | fclose(fd); | |
1630 | } | |
1631 | ||
1632 | //____________________________________________________________________________ | |
1633 | void AliPHOSPIDv1::SetParameterCalibration(Int_t i,Float_t param) | |
1634 | { | |
1635 | // Set parameter "Calibration" i to a value param | |
351dd634 | 1636 | if(i>2 || i<0) { |
1637 | AliError(Form("Invalid parameter number: %d",i)); | |
1638 | } else | |
88cb7938 | 1639 | (*fParameters)(0,i) = param ; |
1640 | } | |
1641 | ||
1642 | //____________________________________________________________________________ | |
1643 | void AliPHOSPIDv1::SetParameterCpv2Emc(Int_t i, TString axis, Float_t cut) | |
1644 | { | |
1645 | // Set the parameters to calculate Cpv-to-Emc Distance Cut depending on | |
1646 | // Purity-Efficiency point i | |
1647 | ||
351dd634 | 1648 | if(i>2 || i<0) { |
1649 | AliError(Form("Invalid parameter number: %d",i)); | |
1650 | } else { | |
88cb7938 | 1651 | axis.ToLower(); |
1652 | if (axis == "x") (*fParameters)(1,i) = cut; | |
1653 | else if (axis == "z") (*fParameters)(2,i) = cut; | |
351dd634 | 1654 | else { |
1655 | AliError(Form("Invalid axis name: %s",axis.Data())); | |
1656 | } | |
88cb7938 | 1657 | } |
1658 | } | |
1659 | ||
1660 | //____________________________________________________________________________ | |
1661 | void AliPHOSPIDv1::SetParameterPhotonBoundary(Int_t i,Float_t param) | |
1662 | { | |
1663 | // Set parameter "Hard photon boundary" i to a value param | |
351dd634 | 1664 | if(i>4 || i<0) { |
1665 | AliError(Form("Invalid parameter number: %d",i)); | |
1666 | } else | |
88cb7938 | 1667 | (*fParameters)(14,i) = param ; |
1668 | } | |
1669 | ||
1670 | //____________________________________________________________________________ | |
1671 | void AliPHOSPIDv1::SetParameterPi0Boundary(Int_t i,Float_t param) | |
1672 | { | |
1673 | // Set parameter "Hard pi0 boundary" i to a value param | |
351dd634 | 1674 | if(i>1 || i<0) { |
1675 | AliError(Form("Invalid parameter number: %d",i)); | |
1676 | } else | |
88cb7938 | 1677 | (*fParameters)(15,i) = param ; |
1678 | } | |
1679 | ||
1680 | //_____________________________________________________________________________ | |
1681 | void AliPHOSPIDv1::SetParameterTimeGate(Int_t i, Float_t gate) | |
1682 | { | |
1683 | // Set the parameter TimeGate depending on Purity-Efficiency point i | |
351dd634 | 1684 | if (i>2 || i<0) { |
1685 | AliError(Form("Invalid Efficiency-Purity choice %d",i)); | |
1686 | } else | |
88cb7938 | 1687 | (*fParameters)(3,i)= gate ; |
1688 | } | |
1689 | ||
1690 | //_____________________________________________________________________________ | |
1691 | void AliPHOSPIDv1::SetParameterToCalculateEllipse(TString particle, TString param, Int_t i, Float_t par) | |
1692 | { | |
1693 | // Set the parameter "i" that is needed to calculate the ellipse | |
1694 | // parameter "param" for a particle "particle" | |
1695 | ||
1696 | particle.ToLower(); | |
1697 | param. ToLower(); | |
1698 | Int_t p= -1; | |
1699 | Int_t offset=0; | |
1700 | ||
1701 | if (particle == "photon") offset=0; | |
1702 | else if (particle == "pi0") offset=5; | |
1703 | else | |
351dd634 | 1704 | AliError(Form("Wrong particle name: %s (choose from pi0/photon)\n", |
1705 | particle.Data())); | |
88cb7938 | 1706 | |
1707 | if (param.Contains("a")) p=4+offset; | |
1708 | else if(param.Contains("b")) p=5+offset; | |
1709 | else if(param.Contains("c")) p=6+offset; | |
1710 | else if(param.Contains("x0"))p=7+offset; | |
1711 | else if(param.Contains("y0"))p=8+offset; | |
351dd634 | 1712 | if((i>4)||(i<0)) { |
1713 | AliError(Form("No parameter with index %d", i)) ; | |
1714 | } else if(p==-1) { | |
1715 | AliError(Form("No parameter with name %s", param.Data() )) ; | |
1716 | } else | |
88cb7938 | 1717 | (*fParameters)(p,i) = par ; |
1718 | } | |
1719 | ||
1720 | //____________________________________________________________________________ | |
1721 | void AliPHOSPIDv1::Unload() | |
1722 | { | |
c947e71a | 1723 | //Unloads RecPoints, Tracks and RecParticles |
88cb7938 | 1724 | AliPHOSGetter * gime = AliPHOSGetter::Instance() ; |
1725 | gime->PhosLoader()->UnloadRecPoints() ; | |
1726 | gime->PhosLoader()->UnloadTracks() ; | |
1727 | gime->PhosLoader()->UnloadRecParticles() ; | |
1728 | } | |
1729 | ||
1730 | //____________________________________________________________________________ | |
90cceaf6 | 1731 | void AliPHOSPIDv1::WriteRecParticles() |
88cb7938 | 1732 | { |
c947e71a | 1733 | //It writes reconstructed particles and pid to file |
1734 | ||
88cb7938 | 1735 | AliPHOSGetter *gime = AliPHOSGetter::Instance() ; |
1736 | ||
1737 | TClonesArray * recParticles = gime->RecParticles() ; | |
1738 | recParticles->Expand(recParticles->GetEntriesFast() ) ; | |
adcca1e6 | 1739 | if(fWrite){ |
1740 | TTree * treeP = gime->TreeP(); | |
1741 | ||
1742 | //First rp | |
1743 | Int_t bufferSize = 32000 ; | |
1744 | TBranch * rpBranch = treeP->Branch("PHOSRP",&recParticles,bufferSize); | |
1745 | rpBranch->SetTitle(BranchName()); | |
1746 | ||
1747 | rpBranch->Fill() ; | |
1748 | ||
1749 | gime->WriteRecParticles("OVERWRITE"); | |
1750 | gime->WritePID("OVERWRITE"); | |
1751 | } | |
88cb7938 | 1752 | } |
407d15b3 | 1753 | //____________________________________________________________________________ |
1754 | void AliPHOSPIDv1::GetVertex(void) | |
1755 | { //extract vertex either using ESD or generator | |
1756 | ||
1757 | //Try to extract vertex from data | |
1758 | if(fESD){ | |
1759 | const AliESDVertex *esdVtx = fESD->GetVertex() ; | |
e5b7b511 | 1760 | if(esdVtx && esdVtx->GetChi2()!=0.){ |
407d15b3 | 1761 | fVtx.SetXYZ(esdVtx->GetXv(),esdVtx->GetYv(),esdVtx->GetZv()) ; |
1762 | return ; | |
1763 | } | |
1764 | } | |
2c06dc7a | 1765 | if(gAlice && gAlice->GetHeader() && gAlice->GetHeader()->GenEventHeader()){ |
1766 | AliGenEventHeader *eh = gAlice->GetHeader()->GenEventHeader() ; | |
1767 | TArrayF ftx ; | |
1768 | eh->PrimaryVertex(ftx); | |
1769 | fVtx.SetXYZ(ftx[0],ftx[1],ftx[2]) ; | |
e5b7b511 | 1770 | return ; |
1771 | } | |
407d15b3 | 1772 | |
1773 | AliWarning("Can not read vertex from data, use fixed \n") ; | |
1774 | fVtx.SetXYZ(0.,0.,0.) ; | |
1775 | ||
1776 | } | |
35adb638 | 1777 | //_______________________________________________________________________ |
1778 | void AliPHOSPIDv1::SetInitPID(const Double_t *p) { | |
1779 | // Sets values for the initial population of each particle type | |
304864ab | 1780 | for (Int_t i=0; i<AliPID::kSPECIESN; i++) fInitPID[i] = p[i]; |
35adb638 | 1781 | } |
1782 | //_______________________________________________________________________ | |
1783 | void AliPHOSPIDv1::GetInitPID(Double_t *p) const { | |
1784 | // Gets values for the initial population of each particle type | |
304864ab | 1785 | for (Int_t i=0; i<AliPID::kSPECIESN; i++) p[i] = fInitPID[i]; |
35adb638 | 1786 | } |