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: |
26 | // -Bit 0 to 2: bit set if RCPV > fCpvEmcDistance (each bit corresponds |
27 | // to a different efficiency-purity point of the photon identification) |
28 | // -Bit 3 to 5: bit set if TOF < fTimeGate (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 fX_center, fY_center, fA, fB, fAngle |
32 | // (each bit corresponds to a different efficiency-purity point of the |
33 | // photon identification) |
a4e98857 |
34 | // |
7acf6008 |
35 | // |
a4e98857 |
36 | // use case: |
37 | // root [0] AliPHOSPIDv1 * p1 = new AliPHOSPIDv1("galice.root") |
38 | // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated |
39 | // root [1] p1->SetIdentificationMethod("disp ellipse") |
40 | // root [2] p1->ExecuteTask() |
148b2bba |
41 | // root [3] AliPHOSPIDv1 * p2 = new AliPHOSPIDv1("galice1.root","v1") |
a4e98857 |
42 | // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated |
148b2bba |
43 | // // reading headers from file galice1.root and create RecParticles with title v1 |
44 | // TrackSegments and RecPoints with title "v1" are used |
7acf6008 |
45 | // // set file name for the branch RecParticles |
148b2bba |
46 | // root [4] p2->ExecuteTask("deb all time") |
7acf6008 |
47 | // // available options |
48 | // // "deb" - prints # of reconstructed particles |
49 | // // "deb all" - prints # and list of RecParticles |
50 | // // "time" - prints benchmarking results |
51 | // |
148b2bba |
52 | // root [5] AliPHOSPIDv1 * p3 = new AliPHOSPIDv1("galice1.root","v1","v0") |
53 | // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated |
54 | // // reading headers from file galice1.root and create RecParticles with title v1 |
55 | // RecPoints and TrackSegments with title "v0" are used |
56 | // root [6] p3->ExecuteTask() |
7acf6008 |
57 | //*-- Author: Yves Schutz (SUBATECH) & Gines Martinez (SUBATECH) & |
148b2bba |
58 | // Gustavo Conesa April 2002 |
6ad0bfa0 |
59 | |
60 | // --- ROOT system --- |
7acf6008 |
61 | #include "TROOT.h" |
62 | #include "TTree.h" |
63 | #include "TFile.h" |
64 | #include "TF2.h" |
65 | #include "TFormula.h" |
66 | #include "TCanvas.h" |
67 | #include "TFolder.h" |
68 | #include "TSystem.h" |
69 | #include "TBenchmark.h" |
148b2bba |
70 | #include "TMatrixD.h" |
71 | #include "TPrincipal.h" |
72 | #include "TSystem.h" |
73 | |
6ad0bfa0 |
74 | // --- Standard library --- |
75 | |
de9ec31b |
76 | #include <iostream.h> |
148b2bba |
77 | #include <fstream.h> |
7acf6008 |
78 | #include <iomanip.h> |
6ad0bfa0 |
79 | |
80 | // --- AliRoot header files --- |
81 | |
7acf6008 |
82 | #include "AliRun.h" |
83 | #include "AliGenerator.h" |
84 | #include "AliPHOS.h" |
26d4b141 |
85 | #include "AliPHOSPIDv1.h" |
7acf6008 |
86 | #include "AliPHOSClusterizerv1.h" |
6ad0bfa0 |
87 | #include "AliPHOSTrackSegment.h" |
7acf6008 |
88 | #include "AliPHOSTrackSegmentMakerv1.h" |
6ad0bfa0 |
89 | #include "AliPHOSRecParticle.h" |
7b7c1533 |
90 | #include "AliPHOSGeometry.h" |
91 | #include "AliPHOSGetter.h" |
6ad0bfa0 |
92 | |
26d4b141 |
93 | ClassImp( AliPHOSPIDv1) |
6ad0bfa0 |
94 | |
1cb7c1ee |
95 | //____________________________________________________________________________ |
96 | AliPHOSPIDv1::AliPHOSPIDv1():AliPHOSPID() |
97 | { |
a4e98857 |
98 | // default ctor |
148b2bba |
99 | |
a496c46c |
100 | fFileName = "" ; |
101 | fFileNamePar = "" ; |
102 | fFrom = "" ; |
103 | fHeaderFileName = "" ; |
e0ed2e49 |
104 | fOptFileName = "Default" ; |
a496c46c |
105 | fTrackSegmentsTitle = "" ; |
106 | fRecPointsTitle = "" ; |
107 | fRecParticlesTitle = "" ; |
108 | |
109 | fNEvent = 0 ; |
110 | fClusterizer = 0 ; |
111 | fTSMaker = 0 ; |
2685bf00 |
112 | fRecParticlesInRun = 0 ; |
a496c46c |
113 | fX = 0 ; |
114 | fP = 0 ; |
115 | fParameters = 0 ; |
116 | |
7acf6008 |
117 | } |
118 | |
119 | //____________________________________________________________________________ |
148b2bba |
120 | AliPHOSPIDv1::AliPHOSPIDv1(const char * headerFile,const char * name, const char * from) : AliPHOSPID(headerFile, name) |
e0ed2e49 |
121 | |
122 | |
7acf6008 |
123 | { |
a4e98857 |
124 | //ctor with the indication on where to look for the track segments |
7b7c1533 |
125 | |
126 | fHeaderFileName = GetTitle() ; |
127 | fTrackSegmentsTitle = GetName() ; |
128 | fRecPointsTitle = GetName() ; |
148b2bba |
129 | fRecParticlesTitle = GetName() ; |
7b7c1533 |
130 | TString tempo(GetName()) ; |
bf8f1fbd |
131 | tempo.Append(":") ; |
7b7c1533 |
132 | tempo.Append(Version()) ; |
bf8f1fbd |
133 | SetName(tempo) ; |
2b60655b |
134 | fRecParticlesInRun = 0 ; |
148b2bba |
135 | if ( from == 0 ) |
136 | fFrom = name ; |
137 | else |
138 | fFrom = from ; |
e0ed2e49 |
139 | fOptFileName = "Default" ; |
2bd5457f |
140 | Init() ; |
7acf6008 |
141 | |
142 | } |
7b7c1533 |
143 | |
7acf6008 |
144 | //____________________________________________________________________________ |
145 | AliPHOSPIDv1::~AliPHOSPIDv1() |
146 | { |
a496c46c |
147 | |
148 | delete [] fX ; // Principal input |
148b2bba |
149 | delete [] fP ; // Principal components |
150 | delete fParameters ; // Matrix of Parameters |
7acf6008 |
151 | } |
2bd5457f |
152 | |
a496c46c |
153 | //____________________________________________________________________________ |
154 | const TString AliPHOSPIDv1::BranchName() const |
155 | { |
156 | TString branchName(GetName() ) ; |
157 | branchName.Remove(branchName.Index(Version())-1) ; |
158 | return branchName ; |
159 | } |
160 | |
148b2bba |
161 | //____________________________________________________________________________ |
162 | void AliPHOSPIDv1::Init() |
163 | { |
164 | // Make all memory allocations that are not possible in default constructor |
165 | // Add the PID task to the list of PHOS tasks |
a496c46c |
166 | |
148b2bba |
167 | if ( strcmp(GetTitle(), "") == 0 ) |
168 | SetTitle("galice.root") ; |
169 | |
a496c46c |
170 | SetParameters() ; // fill the parameters matrix from parameters file |
e0ed2e49 |
171 | |
148b2bba |
172 | AliPHOSGetter * gime = AliPHOSGetter::GetInstance(GetTitle(), fFrom.Data()) ; |
173 | |
a496c46c |
174 | gime->SetRecParticlesTitle(BranchName()) ; |
148b2bba |
175 | if ( gime == 0 ) { |
176 | cerr << "ERROR: AliPHOSPIDv1::Init -> Could not obtain the Getter object !" << endl ; |
177 | return ; |
178 | } |
a496c46c |
179 | |
148b2bba |
180 | gime->PostPID(this) ; |
181 | // create a folder on the white board //YSAlice/WhiteBoard/RecParticles/PHOS/recparticlesName |
a496c46c |
182 | gime->PostRecParticles(BranchName()) ; |
148b2bba |
183 | |
184 | } |
69183710 |
185 | //____________________________________________________________________________ |
148b2bba |
186 | Double_t AliPHOSPIDv1::GetCpvtoEmcDistanceCut(const Float_t Cluster_En, const TString Eff_Pur)const |
187 | { |
188 | // Get CpvtoEmcDistanceCut parameter depending on the cluster energy and |
189 | // Purity-Efficiency point (possible options "HIGH EFFICIENCY" |
190 | // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing |
191 | // EFFICIENCY by PURITY) |
192 | |
193 | Int_t cluster = -1 ; |
194 | Int_t eff_pur = -1 ; |
195 | |
196 | // Check the cluster energy range |
197 | if((Cluster_En > 0.3)&&(Cluster_En <= 1.0)) cluster = 0 ; |
198 | if((Cluster_En > 1.0)&&(Cluster_En <= 2.0)) cluster = 1 ; |
199 | if( Cluster_En > 2.0) cluster = 2 ; |
200 | |
e0ed2e49 |
201 | if(Eff_Pur.Contains("HIGH EFFICIENCY") ||Eff_Pur.Contains("LOW PURITY") ) eff_pur = 0 ; |
202 | if(Eff_Pur.Contains("MEDIUM EFFICIENCY") ||Eff_Pur.Contains("MEDIUM PURITY") ) eff_pur = 1 ; |
203 | if(Eff_Pur.Contains("LOW EFFICIENCY")||Eff_Pur.Contains("HIGH PURITY") ) eff_pur = 2 ; |
148b2bba |
204 | |
205 | if(cluster ==-1){ |
206 | cout<<"Invalid Cluster Energy option"<<endl; |
207 | } |
208 | else if(eff_pur ==-1){ |
209 | cout<<"Invalid Efficiency-Purity option"<<endl; |
210 | } |
211 | else{ |
212 | return (*fParameters)(cluster,eff_pur) ; |
213 | } |
214 | } |
215 | //____________________________________________________________________________ |
216 | |
217 | Double_t AliPHOSPIDv1::GetTimeGate(const Float_t Cluster_En, const TString Eff_Pur) const |
218 | { |
219 | // Get TimeGate parameter depending on the cluster energy and |
220 | // Purity-Efficiency point (possible options "HIGH EFFICIENCY" |
221 | // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing |
222 | // EFFICIENCY by PURITY) |
223 | Int_t cluster = -1 ; |
224 | Int_t eff_pur = -1 ; |
225 | |
226 | if((Cluster_En > 0.3)&&(Cluster_En <= 1.0)) cluster = 0 ; |
227 | if((Cluster_En > 1.0)&&(Cluster_En <= 2.0)) cluster = 1 ; |
228 | if( Cluster_En > 2.0) cluster = 2 ; |
229 | |
e0ed2e49 |
230 | if(Eff_Pur.Contains("HIGH EFFICIENCY") ||Eff_Pur.Contains("LOW PURITY") ) eff_pur = 0 ; |
231 | if(Eff_Pur.Contains("MEDIUM EFFICIENCY") ||Eff_Pur.Contains("MEDIUM PURITY") ) eff_pur = 1 ; |
232 | if(Eff_Pur.Contains("LOW EFFICIENCY")||Eff_Pur.Contains("HIGH PURITY") ) eff_pur = 2 ; |
233 | |
148b2bba |
234 | if(cluster ==-1){ |
235 | cout<<"Invalid Cluster Energy option"<<endl; |
236 | } |
237 | else if(eff_pur ==-1){ |
238 | cout<<"Invalid Efficiency-Purity option"<<endl; |
239 | } |
240 | else{ |
241 | return (*fParameters)(cluster+3,eff_pur) ; |
242 | } |
243 | } |
244 | //_____________________________________________________________________________ |
7acf6008 |
245 | Float_t AliPHOSPIDv1::GetDistance(AliPHOSEmcRecPoint * emc,AliPHOSRecPoint * cpv, Option_t * Axis)const |
69183710 |
246 | { |
247 | // Calculates the distance between the EMC RecPoint and the PPSD RecPoint |
148b2bba |
248 | |
7b7c1533 |
249 | const AliPHOSGeometry * geom = AliPHOSGetter::GetInstance()->PHOSGeometry() ; |
69183710 |
250 | TVector3 vecEmc ; |
7acf6008 |
251 | TVector3 vecCpv ; |
148b2bba |
252 | if(cpv){ |
253 | emc->GetLocalPosition(vecEmc) ; |
254 | cpv->GetLocalPosition(vecCpv) ; |
255 | if(emc->GetPHOSMod() == cpv->GetPHOSMod()){ |
256 | // Correct to difference in CPV and EMC position due to different distance to center. |
257 | // we assume, that particle moves from center |
258 | Float_t dCPV = geom->GetIPtoOuterCoverDistance(); |
259 | Float_t dEMC = geom->GetIPtoCrystalSurface() ; |
260 | dEMC = dEMC / dCPV ; |
261 | vecCpv = dEMC * vecCpv - vecEmc ; |
262 | if (Axis == "X") return vecCpv.X(); |
263 | if (Axis == "Y") return vecCpv.Y(); |
264 | if (Axis == "Z") return vecCpv.Z(); |
265 | if (Axis == "R") return vecCpv.Mag(); |
7acf6008 |
266 | } |
148b2bba |
267 | |
268 | return 100000000 ; |
269 | } |
7acf6008 |
270 | return 100000000 ; |
69183710 |
271 | } |
272 | |
6ad0bfa0 |
273 | //____________________________________________________________________________ |
148b2bba |
274 | Int_t AliPHOSPIDv1::GetPrincipalSign(Double_t* P, Int_t cluster, Int_t eff_pur )const |
275 | { |
276 | //This method gives if the PCA of the particle are inside a defined ellipse |
277 | |
278 | // Get the parameters that define the ellipse stored in the |
279 | // fParameters matrix. |
280 | Double_t X_center = (*fParameters)(cluster+6,eff_pur) ; |
281 | Double_t Y_center = (*fParameters)(cluster+9,eff_pur) ; |
282 | Double_t A = (*fParameters)(cluster+12,eff_pur) ; |
283 | Double_t B = (*fParameters)(cluster+15,eff_pur) ; |
284 | Double_t Angle = (*fParameters)(cluster+18,eff_pur) ; |
285 | |
286 | Int_t prinsign; |
287 | Double_t Dx = 0. ; |
288 | Double_t Delta = 0. ; |
289 | Double_t Y = 0. ; |
290 | Double_t Y_1 = 0. ; |
291 | Double_t Y_2 = 0. ; |
292 | Double_t Pi = TMath::Pi() ; |
293 | Double_t Cos_Theta = TMath::Cos(Pi*Angle/180.) ; |
294 | Double_t Sin_Theta = TMath::Sin(Pi*Angle/180.) ; |
295 | |
296 | Dx = P[0] - X_center ; |
297 | Delta = 4.*A*A*A*B* (A*A*Cos_Theta*Cos_Theta |
298 | + B*B*Sin_Theta*Sin_Theta - Dx*Dx) ; |
299 | if (Delta < 0.) |
300 | {prinsign=0;} |
301 | |
302 | else if (Delta == 0.) |
303 | { |
304 | Y = Cos_Theta*Sin_Theta*(A*A - B*B)*Dx / |
305 | (A*A*Cos_Theta*Cos_Theta + B*B*Sin_Theta*Sin_Theta) ; |
306 | Y += Y_center ; |
307 | if(P[1]==Y ) |
308 | {prinsign=1;} |
309 | else |
310 | {prinsign=0;} |
311 | } |
312 | else |
313 | { |
314 | Y_1 = (Cos_Theta*Sin_Theta*(A*A - B*B) *Dx + |
315 | TMath::Sqrt(Delta)/2.)/(A*A*Cos_Theta*Cos_Theta + |
316 | B*B*Sin_Theta*Sin_Theta) ; |
317 | Y_2 = (Cos_Theta*Sin_Theta*(A*A - B*B) *Dx - |
318 | TMath::Sqrt(Delta)/2.)/(A*A*Cos_Theta*Cos_Theta |
319 | + B*B*Sin_Theta*Sin_Theta) ; |
320 | Y_1 += Y_center ; |
321 | Y_2 += Y_center ; |
322 | if ((P[1]<=Y_1) && (P[1]>=Y_2)) |
323 | {prinsign=1;} |
324 | else |
325 | {prinsign=0;} |
326 | } |
327 | return prinsign; |
328 | } |
329 | |
330 | //____________________________________________________________________________ |
e0ed2e49 |
331 | void AliPHOSPIDv1::SetPrincipalFileOptions(TString OptFileName) { |
332 | |
333 | if(OptFileName.Contains("Small energy range")||OptFileName.Contains("Default")){ |
334 | fFileName = "$ALICE_ROOT/PHOS/PCA8pa15_0.5-5.root" ; |
335 | fFileNamePar = gSystem->ExpandPathName("$ALICE_ROOT/PHOS/Parameters_0.5_5.dat"); |
336 | } |
337 | |
338 | if(OptFileName.Contains("Wide energy range")){ |
339 | fFileName = "$ALICE_ROOT/PHOS/PCA8pa15_0.5-100.root" ; |
340 | fFileNamePar = gSystem->ExpandPathName("$ALICE_ROOT/PHOS/Parameters_0.5_100.dat"); |
341 | } |
342 | } |
343 | |
344 | //____________________________________________________________________________ |
345 | void AliPHOSPIDv1::SetEllipseParameters(Float_t Cluster_En, TString Eff_Pur, Float_t x, Float_t y,Float_t a, Float_t b,Float_t angle) |
148b2bba |
346 | { |
347 | |
348 | // Set all ellipse parameters depending on the cluster energy and |
349 | // Purity-Efficiency point (possible options "HIGH EFFICIENCY" |
350 | // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing |
351 | // EFFICIENCY by PURITY) |
352 | |
353 | Int_t cluster = -1 ; |
354 | Int_t eff_pur = -1 ; |
e0ed2e49 |
355 | |
356 | if((Cluster_En > 0.3)&&(Cluster_En <= 1.0)) cluster = 0 ; |
357 | if((Cluster_En > 1.0)&&(Cluster_En <= 2.0)) cluster = 1 ; |
358 | if( Cluster_En > 2.0) cluster= 2 ; |
359 | |
360 | if(Eff_Pur.Contains("HIGH EFFICIENCY") ||Eff_Pur.Contains("LOW PURITY") ) eff_pur = 0 ; |
361 | if(Eff_Pur.Contains("MEDIUM EFFICIENCY") ||Eff_Pur.Contains("MEDIUM PURITY") ) eff_pur = 1 ; |
362 | if(Eff_Pur.Contains("LOW EFFICIENCY")||Eff_Pur.Contains("HIGH PURITY") ) eff_pur = 2 ; |
363 | |
364 | if(cluster ==-1){ |
148b2bba |
365 | cout<<"Invalid Cluster Energy option"<<endl; |
366 | } |
367 | else if(eff_pur ==-1){ |
368 | cout<<"Invalid Efficiency-Purity option"<<endl; |
369 | } |
370 | else{ |
371 | (*fParameters)(cluster+6,eff_pur) = x ; |
372 | (*fParameters)(cluster+9,eff_pur) = y ; |
373 | (*fParameters)(cluster+12,eff_pur) = a ; |
374 | (*fParameters)(cluster+15,eff_pur) = b ; |
375 | (*fParameters)(cluster+18,eff_pur) = angle ; |
376 | } |
377 | |
378 | } |
379 | //__________________________________________________________________________ |
380 | void AliPHOSPIDv1::SetEllipseXCenter(Float_t Cluster_En, TString Eff_Pur, Float_t x) |
381 | { |
382 | // Set the ellipse parameter x_center depending on the custer energy and |
383 | // Purity-Efficiency point (possible options "HIGH EFFICIENCY" |
384 | // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing |
385 | // EFFICIENCY by PURITY) |
386 | |
387 | Int_t cluster = -1 ; |
388 | Int_t eff_pur = -1 ; |
389 | |
390 | if((Cluster_En > 0.3)&&(Cluster_En <= 1.0)) cluster = 0 ; |
391 | if((Cluster_En > 1.0)&&(Cluster_En <= 2.0)) cluster = 1 ; |
392 | if( Cluster_En > 2.0) cluster = 2 ; |
393 | |
e0ed2e49 |
394 | if(Eff_Pur.Contains("HIGH EFFICIENCY") ||Eff_Pur.Contains("LOW PURITY") ) eff_pur = 0 ; |
395 | if(Eff_Pur.Contains("MEDIUM EFFICIENCY") ||Eff_Pur.Contains("MEDIUM PURITY") ) eff_pur = 1 ; |
396 | if(Eff_Pur.Contains("LOW EFFICIENCY")||Eff_Pur.Contains("HIGH PURITY") ) eff_pur = 2 ; |
148b2bba |
397 | |
398 | if(cluster ==-1){ |
399 | cout<<"Invalid Cluster Energy option"<<endl; |
400 | } |
401 | else if(eff_pur ==-1){ |
402 | cout<<"Invalid Efficiency-Purity option"<<endl; |
403 | } |
404 | else{ |
405 | (*fParameters)(cluster+6,eff_pur) = x ; |
406 | } |
407 | } |
408 | //_________________________________________________________________________ |
409 | void AliPHOSPIDv1::SetEllipseYCenter(Float_t Cluster_En, TString Eff_Pur, Float_t y) |
410 | { |
411 | |
412 | // Set the ellipse parameter y_center depending on the cluster energy and |
413 | // Purity-Efficiency point (possible options "HIGH EFFICIENCY" |
414 | // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing |
415 | // EFFICIENCY by PURITY) |
416 | |
417 | Int_t cluster = -1 ; |
418 | Int_t eff_pur = -1 ; |
419 | |
420 | if((Cluster_En > 0.3)&&(Cluster_En <= 1.0)) cluster = 0 ; |
421 | if((Cluster_En > 1.0)&&(Cluster_En <= 2.0)) cluster = 1 ; |
422 | if( Cluster_En > 2.0) cluster = 2 ; |
423 | |
e0ed2e49 |
424 | if(Eff_Pur.Contains("HIGH EFFICIENCY") ||Eff_Pur.Contains("LOW PURITY") ) eff_pur = 0 ; |
425 | if(Eff_Pur.Contains("MEDIUM EFFICIENCY") ||Eff_Pur.Contains("MEDIUM PURITY") ) eff_pur = 1 ; |
426 | if(Eff_Pur.Contains("LOW EFFICIENCY")||Eff_Pur.Contains("HIGH PURITY") ) eff_pur = 2 ; |
148b2bba |
427 | |
428 | if(cluster ==-1){ |
429 | cout<<"Invalid Cluster Energy option"<<endl; |
430 | } |
431 | else if(eff_pur ==-1){ |
432 | cout<<"Invalid Efficiency-Purity option"<<endl ; |
433 | } |
434 | else{ |
435 | (*fParameters)(cluster+9,eff_pur) = y ; |
436 | } |
437 | } |
438 | //_________________________________________________________________________ |
439 | void AliPHOSPIDv1::SetEllipseAParameter(Float_t Cluster_En, TString Eff_Pur, Float_t a) |
440 | { |
441 | // Set the ellipse parameter a depending on the cluster energy and |
442 | // Purity-Efficiency point (possible options "HIGH EFFICIENCY" |
443 | // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing |
444 | // EFFICIENCY by PURITY) |
445 | |
446 | Int_t cluster = -1 ; |
447 | Int_t eff_pur = -1 ; |
448 | |
449 | if((Cluster_En > 0.3)&&(Cluster_En <= 1.0)) cluster = 0 ; |
450 | if((Cluster_En > 1.0)&&(Cluster_En <= 2.0)) cluster = 1 ; |
451 | if( Cluster_En > 2.0) cluster = 2 ; |
452 | |
e0ed2e49 |
453 | if(Eff_Pur.Contains("HIGH EFFICIENCY") ||Eff_Pur.Contains("LOW PURITY") ) eff_pur = 0 ; |
454 | if(Eff_Pur.Contains("MEDIUM EFFICIENCY") ||Eff_Pur.Contains("MEDIUM PURITY") ) eff_pur = 1 ; |
455 | if(Eff_Pur.Contains("LOW EFFICIENCY")||Eff_Pur.Contains("HIGH PURITY") ) eff_pur = 2 ; |
148b2bba |
456 | |
457 | if(cluster ==-1){ |
458 | cout<<"Invalid Cluster Energy option"<<endl; |
459 | } |
460 | else if(eff_pur ==-1){ |
461 | cout<<"Invalid Efficiency-Purity option"<<endl; |
462 | } |
463 | else{ |
464 | (*fParameters)(cluster+12,eff_pur) = a ; |
465 | } |
466 | } |
467 | //________________________________________________________________________ |
468 | void AliPHOSPIDv1::SetEllipseBParameter(Float_t Cluster_En, TString Eff_Pur, Float_t b) |
469 | { |
470 | // Set the ellipse parameter b depending on the cluster energy and |
471 | // Purity-Efficiency point (possible options "HIGH EFFICIENCY" |
472 | // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing |
473 | // EFFICIENCY by PURITY) |
474 | |
475 | Int_t cluster = -1 ; |
476 | Int_t eff_pur = -1 ; |
e0ed2e49 |
477 | |
148b2bba |
478 | if((Cluster_En > 0.3)&&(Cluster_En <= 1.0)) cluster = 0 ; |
479 | if((Cluster_En > 1.0)&&(Cluster_En <= 2.0)) cluster = 1 ; |
480 | if( Cluster_En > 2.0) cluster = 2 ; |
481 | |
e0ed2e49 |
482 | if(Eff_Pur.Contains("HIGH EFFICIENCY") ||Eff_Pur.Contains("LOW PURITY") ) eff_pur = 0 ; |
483 | if(Eff_Pur.Contains("MEDIUM EFFICIENCY") ||Eff_Pur.Contains("MEDIUM PURITY") ) eff_pur = 1 ; |
484 | if(Eff_Pur.Contains("LOW EFFICIENCY")||Eff_Pur.Contains("HIGH PURITY") ) eff_pur = 2 ; |
485 | |
486 | if(cluster ==-1){ |
487 | cout<<"Invalid Cluster Energy option"<<endl; |
488 | } |
489 | else if(eff_pur ==-1){ |
148b2bba |
490 | cout<<"Invalid Efficiency-Purity option"<<endl; |
e0ed2e49 |
491 | } |
492 | else{ |
493 | (*fParameters)(cluster+15,eff_pur) = b ; |
494 | } |
148b2bba |
495 | } |
496 | //________________________________________________________________________ |
497 | void AliPHOSPIDv1::SetEllipseAngle(Float_t Cluster_En, TString Eff_Pur, Float_t angle) |
498 | { |
499 | |
500 | // Set the ellipse parameter angle depending on the cluster energy and |
501 | // Purity-Efficiency point (possible options "HIGH EFFICIENCY" |
502 | // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing |
503 | // EFFICIENCY by PURITY) |
504 | |
505 | Int_t cluster = -1 ; |
506 | Int_t eff_pur = -1 ; |
507 | |
508 | if((Cluster_En > 0.3)&&(Cluster_En <= 1.0)) cluster = 0 ; |
509 | if((Cluster_En > 1.0)&&(Cluster_En <= 2.0)) cluster = 1 ; |
510 | if( Cluster_En > 2.0) cluster = 2 ; |
511 | |
e0ed2e49 |
512 | if(Eff_Pur.Contains("HIGH EFFICIENCY") ||Eff_Pur.Contains("LOW PURITY") ) eff_pur = 0 ; |
513 | if(Eff_Pur.Contains("MEDIUM EFFICIENCY") ||Eff_Pur.Contains("MEDIUM PURITY") ) eff_pur = 1 ; |
514 | if(Eff_Pur.Contains("LOW EFFICIENCY")||Eff_Pur.Contains("HIGH PURITY") ) eff_pur = 2 ; |
515 | |
148b2bba |
516 | if(cluster ==-1){ |
517 | cout<<"Invalid Cluster Energy option"<<endl; |
518 | } |
519 | else if(eff_pur ==-1){ |
520 | cout<<"Invalid Efficiency-Purity option"<<endl; |
521 | } |
522 | else{ |
523 | (*fParameters)(cluster+18,eff_pur) = angle ; |
524 | } |
525 | } |
526 | //_____________________________________________________________________________ |
527 | void AliPHOSPIDv1::SetCpvtoEmcDistanceCut(Float_t Cluster_En, TString Eff_Pur, Float_t cut) |
528 | { |
529 | |
530 | // Set the parameter Cpvto EmcDistanceCut depending on the cluster energy and |
531 | // Purity-Efficiency point (possible options "HIGH EFFICIENCY" |
532 | // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing |
533 | // EFFICIENCY by PURITY) |
534 | |
535 | Int_t cluster = -1 ; |
536 | Int_t eff_pur = -1 ; |
537 | |
538 | if((Cluster_En > 0.3)&&(Cluster_En <= 1.0)) cluster = 0 ; |
539 | if((Cluster_En > 1.0)&&(Cluster_En <= 2.0)) cluster = 1 ; |
540 | if( Cluster_En > 2.0) cluster = 2 ; |
541 | |
e0ed2e49 |
542 | if(Eff_Pur.Contains("HIGH EFFICIENCY") ||Eff_Pur.Contains("LOW PURITY") ) eff_pur = 0 ; |
543 | if(Eff_Pur.Contains("MEDIUM EFFICIENCY") ||Eff_Pur.Contains("MEDIUM PURITY") ) eff_pur = 1 ; |
544 | if(Eff_Pur.Contains("LOW EFFICIENCY")||Eff_Pur.Contains("HIGH PURITY") ) eff_pur = 2 ; |
545 | |
148b2bba |
546 | if(cluster ==-1){ |
547 | cout<<"Invalid Cluster Energy option"<<endl; |
548 | } |
549 | else if(eff_pur ==-1){ |
550 | cout<<"Invalid Efficiency-Purity option"<<endl; |
551 | } |
552 | else{ |
553 | (*fParameters)(cluster,eff_pur) = cut ; |
554 | } |
555 | } |
a496c46c |
556 | |
557 | //_____________________________________________________________________________ |
558 | void AliPHOSPIDv1::SetParameters() |
559 | { |
560 | // PCA : To do the Principal Components Analysis it is necessary |
561 | // the Principal file, which is opened here |
562 | fX = new double[7]; // Data for the PCA |
563 | fP = new double[7]; // Eigenvalues of the PCA |
e0ed2e49 |
564 | |
565 | SetPrincipalFileOptions(fOptFileName); |
a496c46c |
566 | TFile f( fFileName.Data(), "read" ) ; |
567 | fPrincipal = dynamic_cast<TPrincipal*> (f.Get("principal")) ; |
568 | f.Close() ; |
569 | |
570 | // Initialization of the Parameters matrix. In the File Parameters.dat |
571 | // are all the parameters. These are introduced in a matrix of 21x3 elements. |
572 | // All the parameters defined in this file are, in order of row (there are |
573 | // 3 rows per parameter): CpvtoEmcDistanceCut, TimeGate (and the ellipse |
574 | // parameters), X_center, Y_center, a, b, angle. Each row of a given parameter |
575 | // depends on the cluster energy range (0.3-1,1-2, >2 ) |
576 | // Each column designes the parameters for a point in the Efficiency-Purity |
577 | // of the photon identification P1(0.959,0.625), P2(0.919,0.835), P3(0.833,0.901). |
578 | |
579 | fParameters = new TMatrixD(21,3) ; |
580 | |
a496c46c |
581 | ifstream paramFile(fFileNamePar, ios::in) ; |
582 | |
583 | Int_t i,j ; |
584 | |
585 | for(i = 0; i< 21; i++){ |
586 | for(j = 0; j< 3; j++){ |
587 | paramFile >> (*fParameters)(i,j) ; |
588 | } |
589 | } |
590 | paramFile.close(); |
591 | } |
592 | |
148b2bba |
593 | //_____________________________________________________________________________ |
594 | void AliPHOSPIDv1::SetTimeGate(Float_t Cluster_En, TString Eff_Pur, Float_t gate) |
595 | { |
596 | |
597 | // Set the parameter TimeGate depending on the cluster energy and |
598 | // Purity-Efficiency point (possible options "HIGH EFFICIENCY" |
599 | // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing |
600 | // EFFICIENCY by PURITY) |
601 | |
602 | Int_t cluster = -1 ; |
603 | Int_t eff_pur = -1 ; |
604 | |
605 | if((Cluster_En > 0.3)&&(Cluster_En <= 1.0)) cluster = 0 ; |
606 | if((Cluster_En > 1.0)&&(Cluster_En <= 2.0)) cluster = 1 ; |
607 | if( Cluster_En > 2.0) cluster = 2 ; |
608 | |
e0ed2e49 |
609 | if(Eff_Pur.Contains("HIGH EFFICIENCY") ||Eff_Pur.Contains("LOW PURITY") ) eff_pur = 0 ; |
610 | if(Eff_Pur.Contains("MEDIUM EFFICIENCY") ||Eff_Pur.Contains("MEDIUM PURITY") ) eff_pur = 1 ; |
611 | if(Eff_Pur.Contains("LOW EFFICIENCY")||Eff_Pur.Contains("HIGH PURITY") ) eff_pur = 2 ; |
148b2bba |
612 | |
613 | if(cluster ==-1){ |
614 | cout<<"Invalid Cluster Energy option"<<endl; |
615 | } |
616 | else if(eff_pur ==-1){ |
617 | cout<<"Invalid Efficiency-Purity option"<<endl; |
618 | } |
619 | else{ |
620 | (*fParameters)(cluster+3,eff_pur) = gate ; |
621 | } |
622 | } |
623 | //____________________________________________________________________________ |
624 | |
7acf6008 |
625 | void AliPHOSPIDv1::Exec(Option_t * option) |
6ad0bfa0 |
626 | { |
f035f6ce |
627 | //Steering method |
9688c1dd |
628 | |
bf8f1fbd |
629 | if( strcmp(GetName(), "")== 0 ) |
7acf6008 |
630 | Init() ; |
bf8f1fbd |
631 | |
632 | if(strstr(option,"tim")) |
7acf6008 |
633 | gBenchmark->Start("PHOSPID"); |
bf8f1fbd |
634 | |
635 | if(strstr(option,"print")) { |
7b7c1533 |
636 | Print("") ; |
637 | return ; |
bf8f1fbd |
638 | } |
9688c1dd |
639 | |
148b2bba |
640 | //cout << gDirectory->GetName() << endl ; |
641 | |
bf8f1fbd |
642 | gAlice->GetEvent(0) ; |
148b2bba |
643 | |
7b7c1533 |
644 | //check, if the branch with name of this" already exits? |
645 | TObjArray * lob = (TObjArray*)gAlice->TreeR()->GetListOfBranches() ; |
646 | TIter next(lob) ; |
647 | TBranch * branch = 0 ; |
648 | Bool_t phospidfound = kFALSE, pidfound = kFALSE ; |
649 | |
650 | TString taskName(GetName()) ; |
bf8f1fbd |
651 | taskName.Remove(taskName.Index(Version())-1) ; |
9688c1dd |
652 | |
7b7c1533 |
653 | while ( (branch = (TBranch*)next()) && (!phospidfound || !pidfound) ) { |
654 | if ( (strcmp(branch->GetName(), "PHOSPID")==0) && (strcmp(branch->GetTitle(), taskName.Data())==0) ) |
655 | phospidfound = kTRUE ; |
656 | |
657 | else if ( (strcmp(branch->GetName(), "AliPHOSPID")==0) && (strcmp(branch->GetTitle(), taskName.Data())==0) ) |
658 | pidfound = kTRUE ; |
659 | } |
9688c1dd |
660 | |
7b7c1533 |
661 | if ( phospidfound || pidfound ) { |
662 | cerr << "WARNING: AliPHOSPIDv1::Exec -> RecParticles and/or PIDtMaker branch with name " |
663 | << taskName.Data() << " already exits" << endl ; |
664 | return ; |
665 | } |
bf8f1fbd |
666 | |
7b7c1533 |
667 | Int_t nevents = (Int_t) gAlice->TreeE()->GetEntries() ; |
668 | Int_t ievent ; |
148b2bba |
669 | AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ; |
7b7c1533 |
670 | for(ievent = 0; ievent < nevents; ievent++){ |
a496c46c |
671 | gime->Event(ievent,"R") ; |
148b2bba |
672 | |
7acf6008 |
673 | MakeRecParticles() ; |
9688c1dd |
674 | |
7b7c1533 |
675 | WriteRecParticles(ievent); |
9688c1dd |
676 | |
7acf6008 |
677 | if(strstr(option,"deb")) |
678 | PrintRecParticles(option) ; |
94de8339 |
679 | |
680 | //increment the total number of rec particles per run |
a496c46c |
681 | fRecParticlesInRun += gime->RecParticles(BranchName())->GetEntriesFast() ; |
94de8339 |
682 | |
7acf6008 |
683 | } |
9688c1dd |
684 | |
7acf6008 |
685 | if(strstr(option,"tim")){ |
686 | gBenchmark->Stop("PHOSPID"); |
687 | cout << "AliPHOSPID:" << endl ; |
688 | cout << " took " << gBenchmark->GetCpuTime("PHOSPID") << " seconds for PID " |
7b7c1533 |
689 | << gBenchmark->GetCpuTime("PHOSPID")/nevents << " seconds per event " << endl ; |
7acf6008 |
690 | cout << endl ; |
691 | } |
9688c1dd |
692 | |
7b7c1533 |
693 | } |
694 | |
7acf6008 |
695 | //____________________________________________________________________________ |
696 | void AliPHOSPIDv1::MakeRecParticles(){ |
697 | |
b2a60966 |
698 | // Makes a RecParticle out of a TrackSegment |
148b2bba |
699 | |
7b7c1533 |
700 | AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ; |
148b2bba |
701 | TObjArray * emcRecPoints = gime->EmcRecPoints(fFrom) ; |
702 | TObjArray * cpvRecPoints = gime->CpvRecPoints(fFrom) ; |
703 | TClonesArray * trackSegments = gime->TrackSegments(fFrom) ; |
704 | if ( !emcRecPoints || !cpvRecPoints || !trackSegments ) { |
705 | cerr << "ERROR: AliPHOSPIDv1::MakeRecParticles -> RecPoints or TrackSegments with name " |
706 | << fFrom << " not found ! " << endl ; |
707 | abort() ; |
708 | } |
a496c46c |
709 | TClonesArray * recParticles = gime->RecParticles(BranchName()) ; |
01a599c9 |
710 | recParticles->Clear(); |
148b2bba |
711 | |
712 | |
7b7c1533 |
713 | TIter next(trackSegments) ; |
7acf6008 |
714 | AliPHOSTrackSegment * ts ; |
6ad0bfa0 |
715 | Int_t index = 0 ; |
09fc14a0 |
716 | AliPHOSRecParticle * rp ; |
7acf6008 |
717 | |
7acf6008 |
718 | while ( (ts = (AliPHOSTrackSegment *)next()) ) { |
fad3e5b9 |
719 | |
7b7c1533 |
720 | new( (*recParticles)[index] ) AliPHOSRecParticle() ; |
721 | rp = (AliPHOSRecParticle *)recParticles->At(index) ; |
7acf6008 |
722 | rp->SetTraskSegment(index) ; |
9688c1dd |
723 | rp->SetIndexInList(index) ; |
148b2bba |
724 | |
7acf6008 |
725 | AliPHOSEmcRecPoint * emc = 0 ; |
726 | if(ts->GetEmcIndex()>=0) |
7b7c1533 |
727 | emc = (AliPHOSEmcRecPoint *) emcRecPoints->At(ts->GetEmcIndex()) ; |
fad3e5b9 |
728 | |
7acf6008 |
729 | AliPHOSRecPoint * cpv = 0 ; |
730 | if(ts->GetCpvIndex()>=0) |
7b7c1533 |
731 | cpv = (AliPHOSRecPoint *) cpvRecPoints->At(ts->GetCpvIndex()) ; |
fad3e5b9 |
732 | |
7acf6008 |
733 | //set momentum and energy first |
734 | Float_t e = emc->GetEnergy() ; |
9688c1dd |
735 | TVector3 dir = GetMomentumDirection(emc,cpv) ; |
7acf6008 |
736 | dir.SetMag(e) ; |
737 | |
738 | rp->SetMomentum(dir.X(),dir.Y(),dir.Z(),e) ; |
739 | rp->SetCalcMass(0); |
fad3e5b9 |
740 | |
148b2bba |
741 | // Now set type (reconstructed) of the particle |
742 | |
743 | // Choose the cluster energy range |
744 | Int_t cluster = 0 ; // Ellipse and rcpv cut in function of the cluster energy |
745 | if((e > 0.3)&&(e <= 1.0)) cluster = 0 ; |
746 | if((e > 1.0)&&(e <= 2.0)) cluster = 1 ; |
747 | if( e > 2.0) cluster = 2 ; |
748 | |
749 | // Loop of Efficiency-Purity (the 3 points of purity or efficiency are taken |
750 | // into account to set the particle identification) |
751 | for(Int_t eff_pur = 0; eff_pur < 3 ; eff_pur++){ |
752 | |
753 | // Looking at the CPV detector. If RCPV greater than CpvEmcDistance, 1st, |
754 | // 2nd or 3rd bit (depending on the efficiency-purity point )is set to 1 . |
755 | if(GetDistance(emc, cpv, "R") > (*fParameters)(cluster,eff_pur) ) |
756 | rp->SetPIDBit(eff_pur) ; |
757 | |
758 | // Looking the TOF. If TOF smaller than gate, 4th, 5th or 6th |
759 | // bit (depending on the efficiency-purity point )is set to 1 |
760 | if(emc->GetTime()< (*fParameters)(cluster+3,eff_pur)) |
761 | rp->SetPIDBit(eff_pur+3) ; |
762 | |
763 | // Looking PCA. Define and calculate the data (X), introduce in the function |
764 | // X2P that gives the components (P). |
765 | Float_t fSpher = 0. ; |
766 | Float_t fEmaxdtotal = 0. ; |
767 | Float_t lambda[2] ; |
768 | |
769 | emc->GetElipsAxis(lambda) ; |
770 | if((lambda[0]+lambda[1])!=0) fSpher=fabs(lambda[0]-lambda[1])/(lambda[0]+lambda[1]); |
fad3e5b9 |
771 | |
148b2bba |
772 | fEmaxdtotal=emc->GetMaximalEnergy()/emc->GetEnergy(); |
7acf6008 |
773 | |
148b2bba |
774 | fX[0] = lambda[0] ; |
775 | fX[1] = lambda[1] ; |
776 | fX[2] = emc->GetDispersion() ; |
777 | fX[3] = fSpher ; |
778 | fX[4] = emc->GetMultiplicity() ; |
779 | fX[5] = fEmaxdtotal ; |
780 | fX[6] = emc->GetCoreEnergy() ; |
7acf6008 |
781 | |
148b2bba |
782 | fPrincipal->X2P(fX,fP); |
783 | |
784 | //If we are inside the ellipse, 7th, 8th or 9th |
785 | // bit (depending on the efficiency-purity point )is set to 1 |
786 | if(GetPrincipalSign(fP,cluster,eff_pur) == 1) |
787 | rp->SetPIDBit(eff_pur+6) ; |
788 | |
789 | } |
e0ed2e49 |
790 | |
791 | rp->Name(); //If photon sets the particle pdg name to gamma |
e747b8da |
792 | rp->SetProductionVertex(0,0,0,0); |
793 | rp->SetFirstMother(-1); |
794 | rp->SetLastMother(-1); |
795 | rp->SetFirstDaughter(-1); |
796 | rp->SetLastDaughter(-1); |
797 | rp->SetPolarisation(0,0,0); |
6ad0bfa0 |
798 | index++ ; |
799 | } |
7acf6008 |
800 | |
6ad0bfa0 |
801 | } |
802 | |
09fc14a0 |
803 | //____________________________________________________________________________ |
a496c46c |
804 | void AliPHOSPIDv1:: Print() |
09fc14a0 |
805 | { |
b2a60966 |
806 | // Print the parameters used for the particle type identification |
f035f6ce |
807 | cout << "=============== AliPHOSPID1 ================" << endl ; |
808 | cout << "Making PID "<< endl ; |
809 | cout << " Headers file: " << fHeaderFileName.Data() << endl ; |
810 | cout << " RecPoints branch title: " << fRecPointsTitle.Data() << endl ; |
7b7c1533 |
811 | cout << " TrackSegments Branch title: " << fTrackSegmentsTitle.Data() << endl ; |
812 | cout << " RecParticles Branch title " << fRecParticlesTitle.Data() << endl; |
e0ed2e49 |
813 | cout << " Matrix of Parameters: "<<endl; |
814 | cout << " 3 Columns [High Eff-Low Pur,Medium Eff-Pur, Low Eff-High Pur]"<<endl; |
815 | cout << " 21 Rows, each 3 [ RCPV, TOF, X_Center, Y_Center, A, B, Angle ]"<<endl; |
816 | SetParameters() ; |
a496c46c |
817 | fParameters->Print() ; |
f035f6ce |
818 | cout << "============================================" << endl ; |
09fc14a0 |
819 | } |
820 | |
7acf6008 |
821 | //____________________________________________________________________________ |
7b7c1533 |
822 | void AliPHOSPIDv1::WriteRecParticles(Int_t event) |
09fc14a0 |
823 | { |
7b7c1533 |
824 | |
825 | AliPHOSGetter *gime = AliPHOSGetter::GetInstance() ; |
a496c46c |
826 | |
827 | TClonesArray * recParticles = gime->RecParticles(BranchName()) ; |
bf8f1fbd |
828 | recParticles->Expand(recParticles->GetEntriesFast() ) ; |
9688c1dd |
829 | |
7b7c1533 |
830 | //Make branch in TreeR for RecParticles |
7acf6008 |
831 | char * filename = 0; |
832 | if(gSystem->Getenv("CONFIG_SPLIT_FILE")!=0){ //generating file name |
833 | filename = new char[strlen(gAlice->GetBaseFile())+20] ; |
834 | sprintf(filename,"%s/PHOS.Reco.root",gAlice->GetBaseFile()) ; |
835 | } |
bf8f1fbd |
836 | |
7acf6008 |
837 | TDirectory *cwd = gDirectory; |
838 | |
839 | //First rp |
840 | Int_t bufferSize = 32000 ; |
7b7c1533 |
841 | TBranch * rpBranch = gAlice->TreeR()->Branch("PHOSRP",&recParticles,bufferSize); |
bf8f1fbd |
842 | rpBranch->SetTitle(fRecParticlesTitle); |
7acf6008 |
843 | if (filename) { |
844 | rpBranch->SetFile(filename); |
845 | TIter next( rpBranch->GetListOfBranches()); |
761e34c0 |
846 | TBranch * sb ; |
847 | while ((sb=(TBranch*)next())) { |
848 | sb->SetFile(filename); |
7acf6008 |
849 | } |
850 | cwd->cd(); |
851 | } |
9688c1dd |
852 | |
7acf6008 |
853 | //second, pid |
854 | Int_t splitlevel = 0 ; |
855 | AliPHOSPIDv1 * pid = this ; |
7b7c1533 |
856 | TBranch * pidBranch = gAlice->TreeR()->Branch("AliPHOSPID","AliPHOSPIDv1",&pid,bufferSize,splitlevel); |
9688c1dd |
857 | pidBranch->SetTitle(fRecParticlesTitle.Data()); |
7acf6008 |
858 | if (filename) { |
859 | pidBranch->SetFile(filename); |
860 | TIter next( pidBranch->GetListOfBranches()); |
761e34c0 |
861 | TBranch * sb ; |
862 | while ((sb=(TBranch*)next())) { |
863 | sb->SetFile(filename); |
7acf6008 |
864 | } |
865 | cwd->cd(); |
866 | } |
867 | |
761e34c0 |
868 | rpBranch->Fill() ; |
a496c46c |
869 | pidBranch->Fill() ; |
9688c1dd |
870 | |
7acf6008 |
871 | gAlice->TreeR()->Write(0,kOverwrite) ; |
7b7c1533 |
872 | |
148b2bba |
873 | delete [] filename ; |
7acf6008 |
874 | } |
69183710 |
875 | |
7acf6008 |
876 | //____________________________________________________________________________ |
9688c1dd |
877 | TVector3 AliPHOSPIDv1::GetMomentumDirection(AliPHOSEmcRecPoint * emc, AliPHOSRecPoint * cpv)const |
7acf6008 |
878 | { |
879 | // Calculates the momentum direction: |
880 | // 1. if only a EMC RecPoint, direction is given by IP and this RecPoint |
9688c1dd |
881 | // 2. if a EMC RecPoint and CPV RecPoint, direction is given by the line through the 2 recpoints |
882 | // However because of the poor position resolution of PPSD the direction is always taken as if we were |
7acf6008 |
883 | // in case 1. |
884 | |
885 | TVector3 dir(0,0,0) ; |
886 | |
887 | TVector3 emcglobalpos ; |
888 | TMatrix dummy ; |
889 | |
890 | emc->GetGlobalPosition(emcglobalpos, dummy) ; |
891 | |
148b2bba |
892 | |
7acf6008 |
893 | dir = emcglobalpos ; |
894 | dir.SetZ( -dir.Z() ) ; // why ? |
895 | dir.SetMag(1.) ; |
896 | |
897 | //account correction to the position of IP |
898 | Float_t xo,yo,zo ; //Coordinates of the origin |
899 | gAlice->Generator()->GetOrigin(xo,yo,zo) ; |
900 | TVector3 origin(xo,yo,zo); |
901 | dir = dir - origin ; |
902 | |
903 | return dir ; |
904 | } |
905 | //____________________________________________________________________________ |
a4e98857 |
906 | void AliPHOSPIDv1::PrintRecParticles(Option_t * option) |
907 | { |
dd5c4038 |
908 | // Print table of reconstructed particles |
909 | |
7b7c1533 |
910 | AliPHOSGetter *gime = AliPHOSGetter::GetInstance() ; |
bf8f1fbd |
911 | |
a496c46c |
912 | TClonesArray * recParticles = gime->RecParticles(BranchName()) ; |
9688c1dd |
913 | |
01a599c9 |
914 | cout << "AliPHOSPIDv1: event "<<gAlice->GetEvNumber() << endl ; |
7b7c1533 |
915 | cout << " found " << recParticles->GetEntriesFast() << " RecParticles " << endl ; |
9688c1dd |
916 | |
7acf6008 |
917 | if(strstr(option,"all")) { // printing found TS |
918 | |
e0ed2e49 |
919 | cout << " PARTICLE " |
7acf6008 |
920 | << " Index " << endl ; |
7acf6008 |
921 | |
922 | Int_t index ; |
7b7c1533 |
923 | for (index = 0 ; index < recParticles->GetEntries() ; index++) { |
148b2bba |
924 | AliPHOSRecParticle * rp = (AliPHOSRecParticle * ) recParticles->At(index) ; |
e0ed2e49 |
925 | |
926 | cout << setw(10) << rp->Name() << " " |
927 | << setw(5) << rp->GetIndexInList() << " " <<endl; |
928 | cout << "Type "<< rp->GetType() << endl; |
7acf6008 |
929 | } |
930 | cout << "-------------------------------------------" << endl ; |
931 | } |
932 | |
69183710 |
933 | } |
934 | |
7acf6008 |
935 | |
936 | |