]>
Commit | Line | Data |
---|---|---|
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 | ||
16 | /* $Id$ */ | |
17 | ||
18 | //_________________________________________________________________________ | |
19 | // Implementation version v1 of the PHOS particle identifier | |
20 | // Particle identification based on the | |
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: | |
26 | // -Bit 0 to 2: bit set if RCPV > CpvEmcDistance (each bit corresponds | |
27 | // to a different efficiency-purity point of the photon identification) | |
28 | // -Bit 3 to 5: bit set if TOF < TimeGate (each bit corresponds | |
29 | // to a different efficiency-purity point of the photon identification) | |
30 | // -Bit 6 to 9: bit set if Principal Components are | |
31 | // inside an ellipse defined by the parameters a, b, c, x0 and y0. | |
32 | // (each bit corresponds to a different efficiency-purity point of the | |
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 | |
36 | // file is $ALICE_ROOT/PHOS/PCA8pa15_0.5-100.root. Analysis done for | |
37 | // energies between 0.5 and 100 GeV. | |
38 | // A calibrated energy is calculated. The energy of the reconstructed | |
39 | // cluster is corrected with the formula A + B * E + C * E^2, whose | |
40 | // parameters where obtained through the study of the reconstructed | |
41 | // energy distribution of monoenergetic photons. | |
42 | // | |
43 | // All the parameters (RCPV(2 rows-3 columns),TOF(1r-3c),PCA(5r-4c) | |
44 | // and calibration(1r-3c))are stored in a file called | |
45 | // $ALICE_ROOT/PHOS/Parameters.dat. Each time that AliPHOSPIDv1 is | |
46 | // initialized, this parameters are copied to a Matrix (9,4), a | |
47 | // TMatrixD object. | |
48 | // | |
49 | // use case: | |
50 | // root [0] AliPHOSPIDv1 * p = new AliPHOSPIDv1("galice1.root") | |
51 | // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated | |
52 | // // reading headers from file galice1.root and create RecParticles | |
53 | // TrackSegments and RecPoints are used | |
54 | // // set file name for the branch RecParticles | |
55 | // root [1] p->ExecuteTask("deb all time") | |
56 | // // available options | |
57 | // // "deb" - prints # of reconstructed particles | |
58 | // // "deb all" - prints # and list of RecParticles | |
59 | // // "time" - prints benchmarking results | |
60 | // | |
61 | // root [2] AliPHOSPIDv1 * p2 = new AliPHOSPIDv1("galice1.root","v1",kTRUE) | |
62 | // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated | |
63 | // //Split mode. | |
64 | // root [3] p2->ExecuteTask() | |
65 | // | |
66 | ||
67 | ||
68 | //*-- Author: Yves Schutz (SUBATECH) & Gines Martinez (SUBATECH) & | |
69 | // Gustavo Conesa April 2002 | |
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 | |
74 | // of fixing 3 ellipses for different ranges of energy, it has been | |
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 | // | |
81 | // --- ROOT system --- | |
82 | #include "TROOT.h" | |
83 | #include "TTree.h" | |
84 | #include "TFile.h" | |
85 | #include "TF2.h" | |
86 | #include "TFormula.h" | |
87 | #include "TCanvas.h" | |
88 | #include "TFolder.h" | |
89 | #include "TSystem.h" | |
90 | #include "TBenchmark.h" | |
91 | #include "TMatrixD.h" | |
92 | #include "TPrincipal.h" | |
93 | #include "TSystem.h" | |
94 | ||
95 | // --- Standard library --- | |
96 | ||
97 | ||
98 | // --- AliRoot header files --- | |
99 | ||
100 | #include "AliGenerator.h" | |
101 | #include "AliPHOS.h" | |
102 | #include "AliPHOSPIDv1.h" | |
103 | #include "AliPHOSClusterizerv1.h" | |
104 | #include "AliPHOSEmcRecPoint.h" | |
105 | #include "AliPHOSTrackSegment.h" | |
106 | #include "AliPHOSTrackSegmentMakerv1.h" | |
107 | #include "AliPHOSRecParticle.h" | |
108 | #include "AliPHOSGeometry.h" | |
109 | #include "AliPHOSGetter.h" | |
110 | ||
111 | ClassImp( AliPHOSPIDv1) | |
112 | ||
113 | //____________________________________________________________________________ | |
114 | AliPHOSPIDv1::AliPHOSPIDv1():AliPHOSPID() | |
115 | { | |
116 | // default ctor | |
117 | ||
118 | InitParameters() ; | |
119 | fDefaultInit = kTRUE ; | |
120 | } | |
121 | ||
122 | //____________________________________________________________________________ | |
123 | AliPHOSPIDv1::AliPHOSPIDv1(const AliPHOSPIDv1 & pid ):AliPHOSPID(pid) | |
124 | { | |
125 | // ctor | |
126 | InitParameters() ; | |
127 | Init() ; | |
128 | ||
129 | } | |
130 | ||
131 | //____________________________________________________________________________ | |
132 | AliPHOSPIDv1::AliPHOSPIDv1(const TString alirunFileName, const TString eventFolderName):AliPHOSPID(alirunFileName, eventFolderName) | |
133 | { | |
134 | //ctor with the indication on where to look for the track segments | |
135 | ||
136 | InitParameters() ; | |
137 | Init() ; | |
138 | fDefaultInit = kFALSE ; | |
139 | } | |
140 | ||
141 | //____________________________________________________________________________ | |
142 | AliPHOSPIDv1::~AliPHOSPIDv1() | |
143 | { | |
144 | // dtor | |
145 | ||
146 | delete [] fX ; // Principal input | |
147 | delete [] fPPhoton ; // Photon Principal components | |
148 | delete [] fPPi0 ; // Pi0 Principal components | |
149 | } | |
150 | //____________________________________________________________________________ | |
151 | const TString AliPHOSPIDv1::BranchName() const | |
152 | { | |
153 | ||
154 | return GetName() ; | |
155 | } | |
156 | ||
157 | //____________________________________________________________________________ | |
158 | void AliPHOSPIDv1::Init() | |
159 | { | |
160 | // Make all memory allocations that are not possible in default constructor | |
161 | // Add the PID task to the list of PHOS tasks | |
162 | ||
163 | AliPHOSGetter * gime = AliPHOSGetter::Instance(GetTitle(), fEventFolderName.Data()) ; | |
164 | ||
165 | if ( !gime->PID() ) | |
166 | gime->PostPID(this) ; | |
167 | } | |
168 | ||
169 | //____________________________________________________________________________ | |
170 | void AliPHOSPIDv1::InitParameters() | |
171 | { | |
172 | // Initialize PID parameters | |
173 | fRecParticlesInRun = 0 ; | |
174 | fNEvent = 0 ; | |
175 | fRecParticlesInRun = 0 ; | |
176 | SetParameters() ; // fill the parameters matrix from parameters file | |
177 | SetEventRange(0,-1) ; | |
178 | } | |
179 | ||
180 | //________________________________________________________________________ | |
181 | void AliPHOSPIDv1::Exec(Option_t *option) | |
182 | { | |
183 | // Steering method to perform particle reconstruction and identification | |
184 | // for the event range from fFirstEvent to fLastEvent. | |
185 | // This range is optionally set by SetEventRange(). | |
186 | // if fLastEvent=-1 (by default), then process events until the end. | |
187 | ||
188 | if(strstr(option,"tim")) | |
189 | gBenchmark->Start("PHOSPID"); | |
190 | ||
191 | if(strstr(option,"print")) { | |
192 | Print() ; | |
193 | return ; | |
194 | } | |
195 | ||
196 | ||
197 | AliPHOSGetter * gime = AliPHOSGetter::Instance(GetTitle()) ; | |
198 | ||
199 | if (fLastEvent == -1) | |
200 | fLastEvent = gime->MaxEvent() - 1 ; | |
201 | else | |
202 | fLastEvent = TMath::Min(fLastEvent,gime->MaxEvent()); | |
203 | Int_t nEvents = fLastEvent - fFirstEvent + 1; | |
204 | ||
205 | Int_t ievent ; | |
206 | for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) { | |
207 | gime->Event(ievent,"TR") ; | |
208 | if(gime->TrackSegments() && //Skip events, where no track segments made | |
209 | gime->TrackSegments()->GetEntriesFast()) { | |
210 | MakeRecParticles() ; | |
211 | MakePID() ; | |
212 | WriteRecParticles(); | |
213 | if(strstr(option,"deb")) | |
214 | PrintRecParticles(option) ; | |
215 | //increment the total number of rec particles per run | |
216 | fRecParticlesInRun += gime->RecParticles()->GetEntriesFast() ; | |
217 | } | |
218 | } | |
219 | if(strstr(option,"deb")) | |
220 | PrintRecParticles(option); | |
221 | if(strstr(option,"tim")){ | |
222 | gBenchmark->Stop("PHOSPID"); | |
223 | Info("Exec", "took %f seconds for PID %f seconds per event", | |
224 | gBenchmark->GetCpuTime("PHOSPID"), | |
225 | gBenchmark->GetCpuTime("PHOSPID")/nEvents) ; | |
226 | } | |
227 | Unload(); | |
228 | } | |
229 | ||
230 | //____________________________________________________________________________ | |
231 | const TString AliPHOSPIDv1::GetFileNamePrincipal(TString particle) const | |
232 | { | |
233 | //Get file name that contains the PCA for a particle ("photon or pi0") | |
234 | particle.ToLower(); | |
235 | TString name; | |
236 | if (particle=="photon") name = fFileNamePrincipalPhoton ; | |
237 | else if (particle=="pi0" ) name = fFileNamePrincipalPi0 ; | |
238 | else Error("GetFileNamePrincipal","Wrong particle name: %s (choose from pi0/photon)\n",particle.Data()); | |
239 | return name; | |
240 | } | |
241 | ||
242 | //____________________________________________________________________________ | |
243 | Float_t AliPHOSPIDv1::GetParameterCalibration(Int_t i) const | |
244 | { | |
245 | // Get the i-th parameter "Calibration" | |
246 | Float_t param = 0.; | |
247 | if (i>2 || i<0) | |
248 | Error("GetParameterCalibration","Invalid parameter number: %d",i); | |
249 | else | |
250 | param = (*fParameters)(0,i); | |
251 | return param; | |
252 | } | |
253 | ||
254 | //____________________________________________________________________________ | |
255 | Float_t AliPHOSPIDv1::GetCalibratedEnergy(Float_t e) const | |
256 | { | |
257 | // It calibrates Energy depending on the recpoint energy. | |
258 | // The energy of the reconstructed cluster is corrected with | |
259 | // the formula A + B* E + C* E^2, whose parameters where obtained | |
260 | // through the study of the reconstructed energy distribution of | |
261 | // monoenergetic photons. | |
262 | ||
263 | Float_t p[]={0.,0.,0.}; | |
264 | for (Int_t i=0; i<3; i++) p[i] = GetParameterCalibration(i); | |
265 | Float_t enerec = p[0] + p[1]*e + p[2]*e*e; | |
266 | return enerec ; | |
267 | ||
268 | } | |
269 | ||
270 | //____________________________________________________________________________ | |
271 | Float_t AliPHOSPIDv1::GetParameterCpv2Emc(Int_t i, TString axis) const | |
272 | { | |
273 | // Get the i-th parameter "CPV-EMC distance" for the specified axis | |
274 | Float_t param = 0.; | |
275 | if(i>2 || i<0) | |
276 | Error("GetParameterCpv2Emc","Invalid parameter number: %d",i); | |
277 | else { | |
278 | axis.ToLower(); | |
279 | if (axis == "x") param = (*fParameters)(1,i); | |
280 | else if (axis == "z") param = (*fParameters)(2,i); | |
281 | else Error("GetParameterCpv2Emc","Invalid axis name: %s",axis.Data()); | |
282 | } | |
283 | return param; | |
284 | } | |
285 | ||
286 | //____________________________________________________________________________ | |
287 | Float_t AliPHOSPIDv1::GetCpv2EmcDistanceCut(TString axis, Float_t e) const | |
288 | { | |
289 | // Get CpvtoEmcDistance Cut depending on the cluster energy, axis and | |
290 | // Purity-Efficiency point | |
291 | ||
292 | axis.ToLower(); | |
293 | Float_t p[]={0.,0.,0.}; | |
294 | for (Int_t i=0; i<3; i++) p[i] = GetParameterCpv2Emc(i,axis); | |
295 | Float_t sig = p[0] + TMath::Exp(p[1] - p[2]*e); | |
296 | return sig; | |
297 | } | |
298 | ||
299 | //____________________________________________________________________________ | |
300 | Float_t AliPHOSPIDv1::GetEllipseParameter(TString particle, TString param, Float_t e) const | |
301 | { | |
302 | // Calculates the parameter param of the ellipse | |
303 | ||
304 | particle.ToLower(); | |
305 | param. ToLower(); | |
306 | Float_t p[4]={0.,0.,0.,0.}; | |
307 | Float_t value = 0.0; | |
308 | for (Int_t i=0; i<4; i++) p[i] = GetParameterToCalculateEllipse(particle,param,i); | |
309 | if (particle == "photon") { | |
310 | if (param.Contains("a")) e = TMath::Min((Double_t)e,70.); | |
311 | else if (param.Contains("b")) e = TMath::Min((Double_t)e,70.); | |
312 | else if (param.Contains("x0")) e = TMath::Max((Double_t)e,1.1); | |
313 | } | |
314 | ||
315 | value = p[0]/TMath::Sqrt(e) + p[1]*e + p[2]*e*e + p[3]; | |
316 | return value; | |
317 | } | |
318 | ||
319 | //_____________________________________________________________________________ | |
320 | Float_t AliPHOSPIDv1::GetParameterPhotonBoundary (Int_t i) const | |
321 | { | |
322 | // Get the parameter "i" to calculate the boundary on the moment M2x | |
323 | // for photons at high p_T | |
324 | Float_t param = 0; | |
325 | if (i>3 || i<0) | |
326 | Error("GetParameterPhotonBoundary","Wrong parameter number: %d\n",i); | |
327 | else | |
328 | param = (*fParameters)(14,i) ; | |
329 | return param; | |
330 | } | |
331 | ||
332 | //____________________________________________________________________________ | |
333 | Float_t AliPHOSPIDv1::GetParameterPi0Boundary (Int_t i) const | |
334 | { | |
335 | // Get the parameter "i" to calculate the boundary on the moment M2x | |
336 | // for pi0 at high p_T | |
337 | Float_t param = 0; | |
338 | if (i>2 || i<0) | |
339 | Error("GetParameterPi0Boundary","Wrong parameter number: %d\n",i); | |
340 | else | |
341 | param = (*fParameters)(15,i) ; | |
342 | return param; | |
343 | } | |
344 | ||
345 | //____________________________________________________________________________ | |
346 | Float_t AliPHOSPIDv1::GetParameterTimeGate(Int_t i) const | |
347 | { | |
348 | // Get TimeGate parameter depending on Purity-Efficiency i: | |
349 | // i=0 - Low purity, i=1 - Medium purity, i=2 - High purity | |
350 | Float_t param = 0.; | |
351 | if(i>2 || i<0) | |
352 | Error("GetParameterTimeGate","Invalid Efficiency-Purity choice %d",i); | |
353 | else | |
354 | param = (*fParameters)(3,i) ; | |
355 | return param; | |
356 | } | |
357 | ||
358 | //_____________________________________________________________________________ | |
359 | Float_t AliPHOSPIDv1::GetParameterToCalculateEllipse(TString particle, TString param, Int_t i) const | |
360 | { | |
361 | // Get the parameter "i" that is needed to calculate the ellipse | |
362 | // parameter "param" for the particle "particle" ("photon" or "pi0") | |
363 | ||
364 | particle.ToLower(); | |
365 | param. ToLower(); | |
366 | Int_t offset = -1; | |
367 | if (particle == "photon") offset=0; | |
368 | else if (particle == "pi0") offset=5; | |
369 | else | |
370 | Error("GetParameterToCalculateEllipse","Wrong particle name: %s (choose from pi0/photon)\n",particle.Data()); | |
371 | ||
372 | Int_t p= -1; | |
373 | Float_t par = 0; | |
374 | ||
375 | if (param.Contains("a")) p=4+offset; | |
376 | else if(param.Contains("b")) p=5+offset; | |
377 | else if(param.Contains("c")) p=6+offset; | |
378 | else if(param.Contains("x0"))p=7+offset; | |
379 | else if(param.Contains("y0"))p=8+offset; | |
380 | ||
381 | if (i>4 || i<0) | |
382 | Error("GetParameterToCalculateEllipse", "No parameter with index", i) ; | |
383 | else if (p==-1) | |
384 | Error("GetParameterToCalculateEllipse", "No parameter with name %s", param.Data() ) ; | |
385 | else | |
386 | par = (*fParameters)(p,i) ; | |
387 | ||
388 | return par; | |
389 | } | |
390 | ||
391 | ||
392 | //____________________________________________________________________________ | |
393 | Float_t AliPHOSPIDv1::GetDistance(AliPHOSEmcRecPoint * emc,AliPHOSCpvRecPoint * cpv, Option_t * axis)const | |
394 | { | |
395 | // Calculates the distance between the EMC RecPoint and the PPSD RecPoint | |
396 | ||
397 | const AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ; | |
398 | TVector3 vecEmc ; | |
399 | TVector3 vecCpv ; | |
400 | if(cpv){ | |
401 | emc->GetLocalPosition(vecEmc) ; | |
402 | cpv->GetLocalPosition(vecCpv) ; | |
403 | ||
404 | if(emc->GetPHOSMod() == cpv->GetPHOSMod()){ | |
405 | // Correct to difference in CPV and EMC position due to different distance to center. | |
406 | // we assume, that particle moves from center | |
407 | Float_t dCPV = geom->GetIPtoOuterCoverDistance(); | |
408 | Float_t dEMC = geom->GetIPtoCrystalSurface() ; | |
409 | dEMC = dEMC / dCPV ; | |
410 | vecCpv = dEMC * vecCpv - vecEmc ; | |
411 | if (axis == "X") return vecCpv.X(); | |
412 | if (axis == "Y") return vecCpv.Y(); | |
413 | if (axis == "Z") return vecCpv.Z(); | |
414 | if (axis == "R") return vecCpv.Mag(); | |
415 | } | |
416 | return 100000000 ; | |
417 | } | |
418 | return 100000000 ; | |
419 | } | |
420 | //____________________________________________________________________________ | |
421 | Int_t AliPHOSPIDv1::GetCPVBit(AliPHOSEmcRecPoint * emc,AliPHOSCpvRecPoint * cpv, Int_t effPur, Float_t e) const | |
422 | { | |
423 | if(effPur>2 || effPur<0) | |
424 | Error("GetCPVBit","Invalid Efficiency-Purity choice %d",effPur); | |
425 | ||
426 | Float_t sigX = GetCpv2EmcDistanceCut("X",e); | |
427 | Float_t sigZ = GetCpv2EmcDistanceCut("Z",e); | |
428 | ||
429 | Float_t deltaX = TMath::Abs(GetDistance(emc, cpv, "X")); | |
430 | Float_t deltaZ = TMath::Abs(GetDistance(emc, cpv, "Z")); | |
431 | ||
432 | if((deltaX>sigX*(effPur+1))|(deltaZ>sigZ*(effPur+1))) | |
433 | return 1;//Neutral | |
434 | else | |
435 | return 0;//Charged | |
436 | } | |
437 | ||
438 | //____________________________________________________________________________ | |
439 | Int_t AliPHOSPIDv1::GetPrincipalBit(TString particle, const Double_t* p, Int_t effPur, Float_t e)const | |
440 | { | |
441 | //Is the particle inside de PCA ellipse? | |
442 | ||
443 | particle.ToLower(); | |
444 | Int_t prinbit = 0 ; | |
445 | Float_t a = GetEllipseParameter(particle,"a" , e); | |
446 | Float_t b = GetEllipseParameter(particle,"b" , e); | |
447 | Float_t c = GetEllipseParameter(particle,"c" , e); | |
448 | Float_t x0 = GetEllipseParameter(particle,"x0", e); | |
449 | Float_t y0 = GetEllipseParameter(particle,"y0", e); | |
450 | ||
451 | Float_t r = TMath::Power((p[0] - x0)/a,2) + | |
452 | TMath::Power((p[1] - y0)/b,2) + | |
453 | c*(p[0] - x0)*(p[1] - y0)/(a*b) ; | |
454 | //3 different ellipses defined | |
455 | if((effPur==2) && (r<1./2.)) prinbit= 1; | |
456 | if((effPur==1) && (r<2. )) prinbit= 1; | |
457 | if((effPur==0) && (r<9./2.)) prinbit= 1; | |
458 | ||
459 | if(r<0) | |
460 | Error("GetPrincipalBit", "Negative square?") ; | |
461 | ||
462 | return prinbit; | |
463 | ||
464 | } | |
465 | //____________________________________________________________________________ | |
466 | Int_t AliPHOSPIDv1::GetHardPhotonBit(AliPHOSEmcRecPoint * emc) const | |
467 | { | |
468 | // Set bit for identified hard photons (E > 30 GeV) | |
469 | // if the second moment M2x is below the boundary | |
470 | ||
471 | Float_t e = emc->GetEnergy(); | |
472 | if (e < 30.0) return 0; | |
473 | Float_t m2x = emc->GetM2x(); | |
474 | Float_t m2xBoundary = GetParameterPhotonBoundary(0) * | |
475 | TMath::Exp(-TMath::Power(e-GetParameterPhotonBoundary(1),2)/2.0/ | |
476 | TMath::Power(GetParameterPhotonBoundary(2),2)) + | |
477 | GetParameterPhotonBoundary(3); | |
478 | //Info("GetHardPhotonBit","E=%f, m2x=%f, boundary=%f",e,m2x,m2xBoundary); | |
479 | if (m2x < m2xBoundary) | |
480 | return 1;// A hard photon | |
481 | else | |
482 | return 0;// Not a hard photon | |
483 | } | |
484 | ||
485 | //____________________________________________________________________________ | |
486 | Int_t AliPHOSPIDv1::GetHardPi0Bit(AliPHOSEmcRecPoint * emc) const | |
487 | { | |
488 | // Set bit for identified hard pi0 (E > 30 GeV) | |
489 | // if the second moment M2x is above the boundary | |
490 | ||
491 | Float_t e = emc->GetEnergy(); | |
492 | if (e < 30.0) return 0; | |
493 | Float_t m2x = emc->GetM2x(); | |
494 | Float_t m2xBoundary = GetParameterPi0Boundary(0) + | |
495 | e * GetParameterPi0Boundary(1); | |
496 | //Info("GetHardPi0Bit","E=%f, m2x=%f, boundary=%f",e,m2x,m2xBoundary); | |
497 | if (m2x > m2xBoundary) | |
498 | return 1;// A hard pi0 | |
499 | else | |
500 | return 0;// Not a hard pi0 | |
501 | } | |
502 | ||
503 | //____________________________________________________________________________ | |
504 | TVector3 AliPHOSPIDv1::GetMomentumDirection(AliPHOSEmcRecPoint * emc, AliPHOSCpvRecPoint * )const | |
505 | { | |
506 | // Calculates the momentum direction: | |
507 | // 1. if only a EMC RecPoint, direction is given by IP and this RecPoint | |
508 | // 2. if a EMC RecPoint and CPV RecPoint, direction is given by the line through the 2 recpoints | |
509 | // However because of the poor position resolution of PPSD the direction is always taken as if we were | |
510 | // in case 1. | |
511 | ||
512 | TVector3 dir(0,0,0) ; | |
513 | ||
514 | TVector3 emcglobalpos ; | |
515 | TMatrix dummy ; | |
516 | ||
517 | emc->GetGlobalPosition(emcglobalpos, dummy) ; | |
518 | ||
519 | ||
520 | dir = emcglobalpos ; | |
521 | dir.SetMag(1.) ; | |
522 | ||
523 | //account correction to the position of IP | |
524 | Float_t xo,yo,zo ; //Coordinates of the origin | |
525 | gAlice->Generator()->GetOrigin(xo,yo,zo) ; | |
526 | TVector3 origin(xo,yo,zo); | |
527 | dir = dir - origin ; | |
528 | ||
529 | return dir ; | |
530 | } | |
531 | ||
532 | //____________________________________________________________________________ | |
533 | void AliPHOSPIDv1::MakePID() | |
534 | { | |
535 | // construct the PID weight from a Bayesian Method | |
536 | ||
537 | Int_t index ; | |
538 | Float_t pid1, pid2, pid3, pid4, pid5, pid6 ; | |
539 | pid1 = pid2 = pid3 = pid4 = pid5 = pid6 = 0 ; | |
540 | Int_t nparticles = AliPHOSGetter::Instance()->RecParticles()->GetEntriesFast() ; | |
541 | for(index = 0 ; index < nparticles ; index ++) { | |
542 | AliPHOSRecParticle * recpar = AliPHOSGetter::Instance()->RecParticle(index) ; | |
543 | if (recpar->IsPhoton() || recpar->IsHardPhoton()) pid1++ ; | |
544 | else if (recpar->IsPi0() || recpar->IsHardPi0()) pid2++ ; | |
545 | else if (recpar->IsElectron()) pid3++ ; | |
546 | else if (recpar->IsChargedHadron()) pid4++ ; | |
547 | else if (recpar->IsNeutralHadron()) pid5++ ; | |
548 | else if (recpar->IsEleCon()) pid6++ ; | |
549 | } | |
550 | } | |
551 | ||
552 | //____________________________________________________________________________ | |
553 | void AliPHOSPIDv1::MakeRecParticles() | |
554 | { | |
555 | // Makes a RecParticle out of a TrackSegment | |
556 | ||
557 | AliPHOSGetter * gime = AliPHOSGetter::Instance() ; | |
558 | TObjArray * emcRecPoints = gime->EmcRecPoints() ; | |
559 | TObjArray * cpvRecPoints = gime->CpvRecPoints() ; | |
560 | TClonesArray * trackSegments = gime->TrackSegments() ; | |
561 | if ( !emcRecPoints || !cpvRecPoints || !trackSegments ) { | |
562 | Fatal("MakeRecParticles", "RecPoints or TrackSegments not found !") ; | |
563 | } | |
564 | TClonesArray * recParticles = gime->RecParticles() ; | |
565 | recParticles->Clear(); | |
566 | ||
567 | TIter next(trackSegments) ; | |
568 | AliPHOSTrackSegment * ts ; | |
569 | Int_t index = 0 ; | |
570 | AliPHOSRecParticle * rp ; | |
571 | while ( (ts = (AliPHOSTrackSegment *)next()) ) { | |
572 | ||
573 | new( (*recParticles)[index] ) AliPHOSRecParticle() ; | |
574 | rp = (AliPHOSRecParticle *)recParticles->At(index) ; | |
575 | rp->SetTrackSegment(index) ; | |
576 | rp->SetIndexInList(index) ; | |
577 | ||
578 | AliPHOSEmcRecPoint * emc = 0 ; | |
579 | if(ts->GetEmcIndex()>=0) | |
580 | emc = (AliPHOSEmcRecPoint *) emcRecPoints->At(ts->GetEmcIndex()) ; | |
581 | ||
582 | AliPHOSCpvRecPoint * cpv = 0 ; | |
583 | if(ts->GetCpvIndex()>=0) | |
584 | cpv = (AliPHOSCpvRecPoint *) cpvRecPoints->At(ts->GetCpvIndex()) ; | |
585 | ||
586 | Int_t track = 0 ; | |
587 | track = ts->GetTrackIndex() ; | |
588 | ||
589 | // Now set type (reconstructed) of the particle | |
590 | ||
591 | // Choose the cluster energy range | |
592 | ||
593 | if (!emc) { | |
594 | Fatal("MakeRecParticles", "-> emc(%d) = %d", ts->GetEmcIndex(), emc ) ; | |
595 | } | |
596 | ||
597 | Float_t e = emc->GetEnergy() ; | |
598 | ||
599 | Float_t lambda[2] ; | |
600 | emc->GetElipsAxis(lambda) ; | |
601 | ||
602 | if((lambda[0]>0.01) && (lambda[1]>0.01)){ | |
603 | // Looking PCA. Define and calculate the data (X), | |
604 | // introduce in the function X2P that gives the components (P). | |
605 | ||
606 | Float_t Spher = 0. ; | |
607 | Float_t Emaxdtotal = 0. ; | |
608 | ||
609 | if((lambda[0]+lambda[1])!=0) | |
610 | Spher=fabs(lambda[0]-lambda[1])/(lambda[0]+lambda[1]); | |
611 | ||
612 | Emaxdtotal=emc->GetMaximalEnergy()/emc->GetEnergy(); | |
613 | ||
614 | fX[0] = lambda[0] ; | |
615 | fX[1] = lambda[1] ; | |
616 | fX[2] = emc->GetDispersion() ; | |
617 | fX[3] = Spher ; | |
618 | fX[4] = emc->GetMultiplicity() ; | |
619 | fX[5] = Emaxdtotal ; | |
620 | fX[6] = emc->GetCoreEnergy() ; | |
621 | ||
622 | fPrincipalPhoton->X2P(fX,fPPhoton); | |
623 | fPrincipalPi0 ->X2P(fX,fPPi0); | |
624 | ||
625 | } | |
626 | else{ | |
627 | fPPhoton[0]=-100.0; //We do not accept clusters with | |
628 | fPPhoton[1]=-100.0; //one cell as a photon-like | |
629 | fPPi0[0] =-100.0; | |
630 | fPPi0[1] =-100.0; | |
631 | } | |
632 | ||
633 | Float_t time = emc->GetTime() ; | |
634 | rp->SetTof(time) ; | |
635 | ||
636 | // Loop of Efficiency-Purity (the 3 points of purity or efficiency | |
637 | // are taken into account to set the particle identification) | |
638 | for(Int_t effPur = 0; effPur < 3 ; effPur++){ | |
639 | ||
640 | // Looking at the CPV detector. If RCPV greater than CpvEmcDistance, | |
641 | // 1st,2nd or 3rd bit (depending on the efficiency-purity point ) | |
642 | // is set to 1 | |
643 | if(GetCPVBit(emc, cpv, effPur,e) == 1 ) | |
644 | rp->SetPIDBit(effPur) ; | |
645 | ||
646 | // Looking the TOF. If TOF smaller than gate, 4th, 5th or 6th | |
647 | // bit (depending on the efficiency-purity point )is set to 1 | |
648 | if(time< (*fParameters)(3,effPur)) | |
649 | rp->SetPIDBit(effPur+3) ; | |
650 | ||
651 | //Photon PCA | |
652 | //If we are inside the ellipse, 7th, 8th or 9th | |
653 | // bit (depending on the efficiency-purity point )is set to 1 | |
654 | if(GetPrincipalBit("photon",fPPhoton,effPur,e) == 1) | |
655 | rp->SetPIDBit(effPur+6) ; | |
656 | ||
657 | //Pi0 PCA | |
658 | //If we are inside the ellipse, 10th, 11th or 12th | |
659 | // bit (depending on the efficiency-purity point )is set to 1 | |
660 | if(GetPrincipalBit("pi0" ,fPPi0 ,effPur,e) == 1) | |
661 | rp->SetPIDBit(effPur+9) ; | |
662 | } | |
663 | if(GetHardPhotonBit(emc)) | |
664 | rp->SetPIDBit(12) ; | |
665 | if(GetHardPi0Bit (emc)) | |
666 | rp->SetPIDBit(13) ; | |
667 | ||
668 | if(track >= 0) | |
669 | rp->SetPIDBit(14) ; | |
670 | ||
671 | //Set momentum, energy and other parameters | |
672 | Float_t encal = GetCalibratedEnergy(e); | |
673 | TVector3 dir = GetMomentumDirection(emc,cpv) ; | |
674 | dir.SetMag(encal) ; | |
675 | rp->SetMomentum(dir.X(),dir.Y(),dir.Z(),encal) ; | |
676 | rp->SetCalcMass(0); | |
677 | rp->Name(); //If photon sets the particle pdg name to gamma | |
678 | rp->SetProductionVertex(0,0,0,0); | |
679 | rp->SetFirstMother(-1); | |
680 | rp->SetLastMother(-1); | |
681 | rp->SetFirstDaughter(-1); | |
682 | rp->SetLastDaughter(-1); | |
683 | rp->SetPolarisation(0,0,0); | |
684 | //Set the position in global coordinate system from the RecPoint | |
685 | AliPHOSGeometry * geom = gime->PHOSGeometry() ; | |
686 | AliPHOSTrackSegment * ts = gime->TrackSegment(rp->GetPHOSTSIndex()) ; | |
687 | AliPHOSEmcRecPoint * erp = gime->EmcRecPoint(ts->GetEmcIndex()) ; | |
688 | TVector3 pos ; | |
689 | geom->GetGlobal(erp, pos) ; | |
690 | rp->SetPos(pos); | |
691 | index++ ; | |
692 | } | |
693 | } | |
694 | ||
695 | //____________________________________________________________________________ | |
696 | void AliPHOSPIDv1::Print() const | |
697 | { | |
698 | // Print the parameters used for the particle type identification | |
699 | ||
700 | Info("Print", "=============== AliPHOSPIDv1 ================") ; | |
701 | printf("Making PID\n") ; | |
702 | printf(" Pricipal analysis file from 0.5 to 100 %s\n", fFileNamePrincipalPhoton.Data() ) ; | |
703 | printf(" Name of parameters file %s\n", fFileNameParameters.Data() ) ; | |
704 | printf(" Matrix of Parameters: 14x4\n") ; | |
705 | printf(" Energy Calibration 1x3 [3 parametres to calibrate energy: A + B* E + C * E^2]\n") ; | |
706 | printf(" RCPV 2x3 rows x and z, columns function cut parameters\n") ; | |
707 | printf(" TOF 1x3 [High Eff-Low Pur,Medium Eff-Pur, Low Eff-High Pur]\n") ; | |
708 | printf(" PCA 5x4 [5 ellipse parametres and 4 parametres to calculate them: A/Sqrt(E) + B* E + C * E^2 + D]\n") ; | |
709 | Printf(" Pi0 PCA 5x3 [5 ellipse parametres and 3 parametres to calculate them: A + B* E + C * E^2]\n") ; | |
710 | fParameters->Print() ; | |
711 | } | |
712 | ||
713 | ||
714 | ||
715 | //____________________________________________________________________________ | |
716 | void AliPHOSPIDv1::PrintRecParticles(Option_t * option) | |
717 | { | |
718 | // Print table of reconstructed particles | |
719 | ||
720 | AliPHOSGetter *gime = AliPHOSGetter::Instance() ; | |
721 | ||
722 | TClonesArray * recParticles = gime->RecParticles() ; | |
723 | ||
724 | TString message ; | |
725 | message = "\nevent " ; | |
726 | message += gAlice->GetEvNumber() ; | |
727 | message += " found " ; | |
728 | message += recParticles->GetEntriesFast(); | |
729 | message += " RecParticles\n" ; | |
730 | ||
731 | if(strstr(option,"all")) { // printing found TS | |
732 | message += "\n PARTICLE Index \n" ; | |
733 | ||
734 | Int_t index ; | |
735 | for (index = 0 ; index < recParticles->GetEntries() ; index++) { | |
736 | AliPHOSRecParticle * rp = (AliPHOSRecParticle * ) recParticles->At(index) ; | |
737 | message += "\n" ; | |
738 | message += rp->Name().Data() ; | |
739 | message += " " ; | |
740 | message += rp->GetIndexInList() ; | |
741 | message += " " ; | |
742 | message += rp->GetType() ; | |
743 | } | |
744 | } | |
745 | Info("Print", message.Data() ) ; | |
746 | } | |
747 | ||
748 | //____________________________________________________________________________ | |
749 | void AliPHOSPIDv1::SetParameters() | |
750 | { | |
751 | // PCA : To do the Principal Components Analysis it is necessary | |
752 | // the Principal file, which is opened here | |
753 | fX = new double[7]; // Data for the PCA | |
754 | fPPhoton = new double[7]; // Eigenvalues of the PCA | |
755 | fPPi0 = new double[7]; // Eigenvalues of the Pi0 PCA | |
756 | ||
757 | // Read photon principals from the photon file | |
758 | ||
759 | fFileNamePrincipalPhoton = "$ALICE_ROOT/PHOS/PCA8pa15_0.5-100.root" ; | |
760 | TFile f( fFileNamePrincipalPhoton.Data(), "read" ) ; | |
761 | fPrincipalPhoton = dynamic_cast<TPrincipal*> (f.Get("principal")) ; | |
762 | f.Close() ; | |
763 | ||
764 | // Read pi0 principals from the pi0 file | |
765 | ||
766 | fFileNamePrincipalPi0 = "$ALICE_ROOT/PHOS/PCA_pi0_40-120.root" ; | |
767 | TFile fPi0( fFileNamePrincipalPi0.Data(), "read" ) ; | |
768 | fPrincipalPi0 = dynamic_cast<TPrincipal*> (fPi0.Get("principal")) ; | |
769 | fPi0.Close() ; | |
770 | ||
771 | // Open parameters file and initialization of the Parameters matrix. | |
772 | // In the File Parameters.dat are all the parameters. These are introduced | |
773 | // in a matrix of 16x4 | |
774 | // | |
775 | // All the parameters defined in this file are, in order of row: | |
776 | // line 0 : calibration | |
777 | // lines 1,2 : CPV rectangular cat for X and Z | |
778 | // line 3 : TOF cut | |
779 | // lines 4-8 : parameters to calculate photon PCA ellipse | |
780 | // lines 9-13: parameters to calculate pi0 PCA ellipse | |
781 | // lines 14-15: parameters to calculate border for high-pt photons and pi0 | |
782 | ||
783 | fFileNameParameters = gSystem->ExpandPathName("$ALICE_ROOT/PHOS/Parameters.dat"); | |
784 | fParameters = new TMatrix(16,4) ; | |
785 | const Int_t maxLeng=255; | |
786 | char string[maxLeng]; | |
787 | ||
788 | // Open a text file with PID parameters | |
789 | FILE *fd = fopen(fFileNameParameters.Data(),"r"); | |
790 | if (!fd) | |
791 | Fatal("SetParameter","File %s with a PID parameters cannot be opened\n", | |
792 | fFileNameParameters.Data()); | |
793 | ||
794 | Int_t i=0; | |
795 | // Read parameter file line-by-line and skip empty line and comments | |
796 | while (fgets(string,maxLeng,fd) != NULL) { | |
797 | if (string[0] == '\n' ) continue; | |
798 | if (string[0] == '!' ) continue; | |
799 | sscanf(string, "%f %f %f %f", | |
800 | &(*fParameters)(i,0), &(*fParameters)(i,1), | |
801 | &(*fParameters)(i,2), &(*fParameters)(i,3)); | |
802 | i++; | |
803 | //Info("SetParameters", "line %d: %s",i,string); | |
804 | } | |
805 | fclose(fd); | |
806 | } | |
807 | ||
808 | //____________________________________________________________________________ | |
809 | void AliPHOSPIDv1::SetParameterCalibration(Int_t i,Float_t param) | |
810 | { | |
811 | // Set parameter "Calibration" i to a value param | |
812 | if(i>2 || i<0) | |
813 | Error("SetParameterCalibration","Invalid parameter number: %d",i); | |
814 | else | |
815 | (*fParameters)(0,i) = param ; | |
816 | } | |
817 | ||
818 | //____________________________________________________________________________ | |
819 | void AliPHOSPIDv1::SetParameterCpv2Emc(Int_t i, TString axis, Float_t cut) | |
820 | { | |
821 | // Set the parameters to calculate Cpv-to-Emc Distance Cut depending on | |
822 | // Purity-Efficiency point i | |
823 | ||
824 | if(i>2 || i<0) | |
825 | Error("SetParameterCpv2Emc","Invalid parameter number: %d",i); | |
826 | else { | |
827 | axis.ToLower(); | |
828 | if (axis == "x") (*fParameters)(1,i) = cut; | |
829 | else if (axis == "z") (*fParameters)(2,i) = cut; | |
830 | else Error("SetParameterCpv2Emc","Invalid axis name: %s",axis.Data()); | |
831 | } | |
832 | } | |
833 | ||
834 | //____________________________________________________________________________ | |
835 | void AliPHOSPIDv1::SetParameterPhotonBoundary(Int_t i,Float_t param) | |
836 | { | |
837 | // Set parameter "Hard photon boundary" i to a value param | |
838 | if(i>4 || i<0) | |
839 | Error("SetParameterPhotonBoundary","Invalid parameter number: %d",i); | |
840 | else | |
841 | (*fParameters)(14,i) = param ; | |
842 | } | |
843 | ||
844 | //____________________________________________________________________________ | |
845 | void AliPHOSPIDv1::SetParameterPi0Boundary(Int_t i,Float_t param) | |
846 | { | |
847 | // Set parameter "Hard pi0 boundary" i to a value param | |
848 | if(i>1 || i<0) | |
849 | Error("SetParameterPi0Boundary","Invalid parameter number: %d",i); | |
850 | else | |
851 | (*fParameters)(15,i) = param ; | |
852 | } | |
853 | ||
854 | //_____________________________________________________________________________ | |
855 | void AliPHOSPIDv1::SetParameterTimeGate(Int_t i, Float_t gate) | |
856 | { | |
857 | // Set the parameter TimeGate depending on Purity-Efficiency point i | |
858 | if (i>2 || i<0) | |
859 | Error("SetParameterTimeGate","Invalid Efficiency-Purity choice %d",i); | |
860 | else | |
861 | (*fParameters)(3,i)= gate ; | |
862 | } | |
863 | ||
864 | //_____________________________________________________________________________ | |
865 | void AliPHOSPIDv1::SetParameterToCalculateEllipse(TString particle, TString param, Int_t i, Float_t par) | |
866 | { | |
867 | // Set the parameter "i" that is needed to calculate the ellipse | |
868 | // parameter "param" for a particle "particle" | |
869 | ||
870 | particle.ToLower(); | |
871 | param. ToLower(); | |
872 | Int_t p= -1; | |
873 | Int_t offset=0; | |
874 | ||
875 | if (particle == "photon") offset=0; | |
876 | else if (particle == "pi0") offset=5; | |
877 | else | |
878 | Error("SetParameterToCalculateEllipse","Wrong particle name: %s (choose from pi0/photon)\n",particle.Data()); | |
879 | ||
880 | if (param.Contains("a")) p=4+offset; | |
881 | else if(param.Contains("b")) p=5+offset; | |
882 | else if(param.Contains("c")) p=6+offset; | |
883 | else if(param.Contains("x0"))p=7+offset; | |
884 | else if(param.Contains("y0"))p=8+offset; | |
885 | if((i>4)||(i<0)) | |
886 | Error("SetEllipseParameter", "No parameter with index %d", i) ; | |
887 | else if(p==-1) | |
888 | Error("SetEllipseParameter", "No parameter with name %s", param.Data() ) ; | |
889 | else | |
890 | (*fParameters)(p,i) = par ; | |
891 | } | |
892 | ||
893 | //____________________________________________________________________________ | |
894 | void AliPHOSPIDv1::Unload() | |
895 | { | |
896 | AliPHOSGetter * gime = AliPHOSGetter::Instance() ; | |
897 | gime->PhosLoader()->UnloadRecPoints() ; | |
898 | gime->PhosLoader()->UnloadTracks() ; | |
899 | gime->PhosLoader()->UnloadRecParticles() ; | |
900 | } | |
901 | ||
902 | //____________________________________________________________________________ | |
903 | void AliPHOSPIDv1::WriteRecParticles() | |
904 | { | |
905 | ||
906 | AliPHOSGetter *gime = AliPHOSGetter::Instance() ; | |
907 | ||
908 | TClonesArray * recParticles = gime->RecParticles() ; | |
909 | recParticles->Expand(recParticles->GetEntriesFast() ) ; | |
910 | TTree * treeP = gime->TreeP(); | |
911 | ||
912 | //First rp | |
913 | Int_t bufferSize = 32000 ; | |
914 | TBranch * rpBranch = treeP->Branch("PHOSRP",&recParticles,bufferSize); | |
915 | rpBranch->SetTitle(BranchName()); | |
916 | ||
917 | rpBranch->Fill() ; | |
918 | ||
919 | gime->WriteRecParticles("OVERWRITE"); | |
920 | gime->WritePID("OVERWRITE"); | |
921 | } | |
922 |