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