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" |
85 | #include "TF2.h" |
86 | #include "TFormula.h" |
87 | #include "TCanvas.h" |
88 | #include "TFolder.h" |
89 | #include "TSystem.h" |
90 | #include "TBenchmark.h" |
148b2bba |
91 | #include "TMatrixD.h" |
92 | #include "TPrincipal.h" |
93 | #include "TSystem.h" |
94 | |
6ad0bfa0 |
95 | // --- Standard library --- |
96 | |
75a6835b |
97 | #include <Riostream.h> |
98 | |
6ad0bfa0 |
99 | // --- AliRoot header files --- |
100 | |
7acf6008 |
101 | #include "AliRun.h" |
102 | #include "AliGenerator.h" |
103 | #include "AliPHOS.h" |
26d4b141 |
104 | #include "AliPHOSPIDv1.h" |
7acf6008 |
105 | #include "AliPHOSClusterizerv1.h" |
6ad0bfa0 |
106 | #include "AliPHOSTrackSegment.h" |
7acf6008 |
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 | |
124 | //____________________________________________________________________________ |
fbf811ec |
125 | AliPHOSPIDv1::AliPHOSPIDv1(const char * headerFile,const char * name, const Bool_t toSplit) |
126 | :AliPHOSPID(headerFile, name,toSplit) |
7acf6008 |
127 | { |
a4e98857 |
128 | //ctor with the indication on where to look for the track segments |
7b7c1533 |
129 | |
8d0f3f77 |
130 | InitParameters() ; |
131 | |
2bd5457f |
132 | Init() ; |
92f521a9 |
133 | fDefaultInit = kFALSE ; |
7acf6008 |
134 | |
135 | } |
7b7c1533 |
136 | |
7acf6008 |
137 | //____________________________________________________________________________ |
138 | AliPHOSPIDv1::~AliPHOSPIDv1() |
139 | { |
79bb1b62 |
140 | // dtor |
92f521a9 |
141 | // fDefaultInit = kTRUE if PID created by default ctor (to get just the parameters) |
9fa5f1d0 |
142 | |
a496c46c |
143 | delete [] fX ; // Principal input |
148b2bba |
144 | delete [] fP ; // Principal components |
8d0f3f77 |
145 | |
92f521a9 |
146 | if (!fDefaultInit) { |
fbf811ec |
147 | // AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ; |
92f521a9 |
148 | // remove the task from the folder list |
fbf811ec |
149 | // gime->RemoveTask("P",GetName()) ; |
150 | // TString name(GetName()) ; |
151 | // name.ReplaceAll("pid", "clu") ; |
152 | // gime->RemoveTask("C",name) ; |
92f521a9 |
153 | |
fbf811ec |
154 | // // remove the data from the folder list |
155 | // name = GetName() ; |
156 | // name.Remove(name.Index(":")) ; |
157 | // gime->RemoveObjects("RE", name) ; // EMCARecPoints |
158 | // gime->RemoveObjects("RC", name) ; // CPVRecPoints |
159 | // gime->RemoveObjects("T", name) ; // TrackSegments |
160 | // gime->RemoveObjects("P", name) ; // RecParticles |
92f521a9 |
161 | |
fbf811ec |
162 | // // Delete gAlice |
163 | // gime->CloseFile() ; |
92f521a9 |
164 | |
165 | fSplitFile = 0 ; |
79bb1b62 |
166 | } |
7acf6008 |
167 | } |
2bd5457f |
168 | |
a496c46c |
169 | //____________________________________________________________________________ |
170 | const TString AliPHOSPIDv1::BranchName() const |
171 | { |
172 | TString branchName(GetName() ) ; |
173 | branchName.Remove(branchName.Index(Version())-1) ; |
174 | return branchName ; |
175 | } |
176 | |
148b2bba |
177 | //____________________________________________________________________________ |
178 | void AliPHOSPIDv1::Init() |
179 | { |
180 | // Make all memory allocations that are not possible in default constructor |
181 | // Add the PID task to the list of PHOS tasks |
a496c46c |
182 | |
148b2bba |
183 | if ( strcmp(GetTitle(), "") == 0 ) |
184 | SetTitle("galice.root") ; |
148b2bba |
185 | |
fbf811ec |
186 | TString branchname(GetName()) ; |
187 | branchname.Remove(branchname.Index(Version())-1) ; |
188 | AliPHOSGetter * gime = AliPHOSGetter::GetInstance(GetTitle(),branchname.Data(),fToSplit ) ; |
189 | |
190 | // gime->SetRecParticlesTitle(BranchName()) ; |
148b2bba |
191 | if ( gime == 0 ) { |
21cd0c07 |
192 | Error("Init", "Could not obtain the Getter object !" ) ; |
148b2bba |
193 | return ; |
194 | } |
fbf811ec |
195 | |
196 | fSplitFile = 0 ; |
197 | if(fToSplit){ |
198 | //First - extract full path if necessary |
199 | TString fileName(GetTitle()) ; |
200 | Ssiz_t islash = fileName.Last('/') ; |
201 | if(islash<fileName.Length()) |
202 | fileName.Remove(islash+1,fileName.Length()) ; |
203 | else |
204 | fileName="" ; |
205 | fileName+="PHOS.RecData." ; |
206 | if((strcmp(branchname.Data(),"Default")!=0)&&(strcmp(branchname.Data(),"")!=0)){ |
207 | fileName+=branchname.Data() ; |
208 | fileName+="." ; |
209 | } |
210 | fileName+="root" ; |
211 | fSplitFile = static_cast<TFile*>(gROOT->GetFile(fileName.Data())); |
212 | if(!fSplitFile) |
213 | fSplitFile = TFile::Open(fileName.Data(),"update") ; |
214 | } |
a496c46c |
215 | |
148b2bba |
216 | gime->PostPID(this) ; |
217 | // create a folder on the white board //YSAlice/WhiteBoard/RecParticles/PHOS/recparticlesName |
fbf811ec |
218 | gime->PostRecParticles(branchname) ; |
148b2bba |
219 | |
220 | } |
8d0f3f77 |
221 | |
222 | //____________________________________________________________________________ |
223 | void AliPHOSPIDv1::InitParameters() |
224 | { |
fbf811ec |
225 | // fFrom = "" ; |
226 | // fHeaderFileName = GetTitle() ; |
227 | // TString name(GetName()) ; |
228 | // if (name.IsNull()) |
229 | // name = "Default" ; |
230 | // fTrackSegmentsTitle = name ; |
231 | // fRecPointsTitle = name ; |
232 | // fRecParticlesTitle = name ; |
233 | // name.Append(":") ; |
234 | // name.Append(Version()) ; |
235 | // SetName(name) ; |
8d0f3f77 |
236 | fRecParticlesInRun = 0 ; |
8d0f3f77 |
237 | fNEvent = 0 ; |
fbf811ec |
238 | // fClusterizer = 0 ; |
239 | // fTSMaker = 0 ; |
8d0f3f77 |
240 | fRecParticlesInRun = 0 ; |
fbf811ec |
241 | TString pidName( GetName()) ; |
242 | if (pidName.IsNull() ) |
243 | pidName = "Default" ; |
244 | pidName.Append(":") ; |
245 | pidName.Append(Version()) ; |
246 | SetName(pidName) ; |
bc0c084c |
247 | fPi0Analysis = kFALSE ; |
9fa5f1d0 |
248 | SetParameters() ; // fill the parameters matrix from parameters file |
8d0f3f77 |
249 | } |
250 | |
69183710 |
251 | //____________________________________________________________________________ |
bc0c084c |
252 | const Float_t AliPHOSPIDv1::GetCpvtoEmcDistanceCut(const Float_t e, const TString Axis) const |
148b2bba |
253 | { |
bc0c084c |
254 | // Get CpvtoEmcDistance Cut depending on the cluster energy, axis and |
255 | // Purity-Efficiency point |
256 | |
257 | Int_t i = -1; |
258 | if(Axis.Contains("X"))i = 0; |
259 | else if(Axis.Contains("Z"))i = 1; |
9fa5f1d0 |
260 | else |
bc0c084c |
261 | Error("GetCpvtoEmcDistanceCut"," Invalid axis option "); |
262 | |
263 | Float_t a = (*fParameters)(i,0) ; |
264 | Float_t b = (*fParameters)(i,1) ; |
265 | Float_t c = (*fParameters)(i,2) ; |
266 | |
267 | Float_t sig = a + TMath::Exp(b-c*e); |
268 | return sig; |
148b2bba |
269 | } |
270 | //____________________________________________________________________________ |
271 | |
bc0c084c |
272 | const Double_t AliPHOSPIDv1::GetTimeGate(const Int_t Eff_Pur) const |
148b2bba |
273 | { |
bc0c084c |
274 | // Get TimeGate parameter depending on Purity-Efficiency point |
9fa5f1d0 |
275 | |
bc0c084c |
276 | if(Eff_Pur>2 || Eff_Pur<0) |
277 | Error("GetTimeGate","Invalid Efficiency-Purity choice %d",Eff_Pur); |
278 | return (*fParameters)(2,Eff_Pur) ; |
f0a4c9e9 |
279 | |
148b2bba |
280 | } |
281 | //_____________________________________________________________________________ |
50739f15 |
282 | const Float_t AliPHOSPIDv1::GetDistance(AliPHOSEmcRecPoint * emc,AliPHOSRecPoint * cpv, Option_t * Axis)const |
69183710 |
283 | { |
284 | // Calculates the distance between the EMC RecPoint and the PPSD RecPoint |
148b2bba |
285 | |
7b7c1533 |
286 | const AliPHOSGeometry * geom = AliPHOSGetter::GetInstance()->PHOSGeometry() ; |
69183710 |
287 | TVector3 vecEmc ; |
7acf6008 |
288 | TVector3 vecCpv ; |
148b2bba |
289 | if(cpv){ |
290 | emc->GetLocalPosition(vecEmc) ; |
291 | cpv->GetLocalPosition(vecCpv) ; |
292 | if(emc->GetPHOSMod() == cpv->GetPHOSMod()){ |
293 | // Correct to difference in CPV and EMC position due to different distance to center. |
294 | // we assume, that particle moves from center |
295 | Float_t dCPV = geom->GetIPtoOuterCoverDistance(); |
296 | Float_t dEMC = geom->GetIPtoCrystalSurface() ; |
297 | dEMC = dEMC / dCPV ; |
298 | vecCpv = dEMC * vecCpv - vecEmc ; |
299 | if (Axis == "X") return vecCpv.X(); |
300 | if (Axis == "Y") return vecCpv.Y(); |
301 | if (Axis == "Z") return vecCpv.Z(); |
302 | if (Axis == "R") return vecCpv.Mag(); |
7acf6008 |
303 | } |
148b2bba |
304 | |
305 | return 100000000 ; |
306 | } |
7acf6008 |
307 | return 100000000 ; |
69183710 |
308 | } |
bc0c084c |
309 | //____________________________________________________________________________ |
310 | const Int_t AliPHOSPIDv1::GetCPVBit(AliPHOSEmcRecPoint * emc,AliPHOSRecPoint * cpv,const Int_t EffPur, const Float_t e) const |
311 | { |
312 | if(EffPur>2 || EffPur<0) |
313 | Error("GetCPVBit","Invalid Efficiency-Purity choice %d",EffPur); |
314 | |
315 | Float_t sigX = GetCpvtoEmcDistanceCut(e,"X"); |
316 | Float_t sigZ = GetCpvtoEmcDistanceCut(e,"Z"); |
317 | |
318 | Float_t deltaX = TMath::Abs(GetDistance(emc, cpv, "X")); |
319 | Float_t deltaZ = TMath::Abs(GetDistance(emc, cpv, "Z")); |
320 | |
321 | if((deltaX>sigX*(EffPur+1))|(deltaZ>sigZ*(EffPur+1))) |
322 | return 1;//Neutral |
323 | else |
324 | return 0;//Charged |
325 | |
326 | } |
69183710 |
327 | |
6ad0bfa0 |
328 | //____________________________________________________________________________ |
50739f15 |
329 | const Double_t AliPHOSPIDv1::GetCalibratedEnergy(const Float_t e) const |
330 | { |
331 | // It calibrates Energy depending on the recpoint energy. |
332 | // The energy of the reconstructed cluster is corrected with |
333 | // the formula A + B* E + C* E^2, whose parameters where obtained |
334 | // through the study of the reconstructed energy distribution of |
335 | // monoenergetic photons. |
336 | |
337 | Double_t p[]={0.,0.,0.}; |
338 | Int_t i; |
bc0c084c |
339 | for(i=0;i<3;i++) p[i]= (*fParameters)(8,i); |
50739f15 |
340 | Double_t enerec = p[0] + p[1]* e+ p[2] * e * e; |
9fa5f1d0 |
341 | return enerec ; |
342 | |
343 | } |
344 | //____________________________________________________________________________ |
bc0c084c |
345 | const Int_t AliPHOSPIDv1::GetPrincipalBit(const Double_t* P,const Int_t eff_pur, const Float_t E)const |
148b2bba |
346 | { |
50739f15 |
347 | //Is the particle inside de PCA ellipse? |
348 | |
349 | Int_t prinsign= 0 ; |
350 | Double_t A = GetEllipseParameter("a", E); |
351 | Double_t B = GetEllipseParameter("b", E); |
352 | Double_t C = GetEllipseParameter("c", E); |
353 | Double_t X_center = GetEllipseParameter("x0", E); |
354 | Double_t Y_center = GetEllipseParameter("y0", E); |
148b2bba |
355 | |
50739f15 |
356 | Double_t R = TMath::Power((P[0] - X_center)/A,2) + |
357 | TMath::Power((P[1] - Y_center)/B,2) + |
358 | C*(P[0] - X_center)*(P[1] - Y_center)/(A*B) ; |
359 | //3 different ellipses defined |
360 | if((eff_pur==2)&&(R <1./2.)) prinsign= 1; |
361 | if((eff_pur==1)&&(R <2. )) prinsign= 1; |
362 | if((eff_pur==0)&&(R <9./2.)) prinsign= 1; |
363 | |
21cd0c07 |
364 | if(R<0) |
365 | Error("GetPrincipalSign", "Negative square?") ; |
148b2bba |
366 | return prinsign; |
148b2bba |
367 | |
148b2bba |
368 | } |
148b2bba |
369 | |
148b2bba |
370 | //_____________________________________________________________________________ |
bc0c084c |
371 | void AliPHOSPIDv1::SetCpvtoEmcDistanceCutParameters(Float_t e, Int_t Eff_Pur, TString Axis,Float_t cut) |
148b2bba |
372 | { |
bc0c084c |
373 | // Set the parameters to calculate Cpvto EmcDistanceCut depending on the cluster energy and |
374 | // Purity-Efficiency point |
148b2bba |
375 | |
bc0c084c |
376 | if(Eff_Pur>2 || Eff_Pur<0) |
377 | Error("SetCpvtoEmcDistanceCutParameters","Invalid Efficiency-Purity choice %d",Eff_Pur); |
92f521a9 |
378 | |
bc0c084c |
379 | Int_t i = -1; |
380 | if(Axis.Contains("X"))i = 0; |
381 | else if(Axis.Contains("Z"))i = 1; |
382 | else |
383 | Error("SetCpvtoEmcDistanceCutParameters"," Invalid axis option"); |
384 | |
385 | (*fParameters)(i,Eff_Pur) = cut ; |
f0a4c9e9 |
386 | } |
387 | //_____________________________________________________________________________ |
bc0c084c |
388 | void AliPHOSPIDv1::SetTimeGate(Int_t Eff_Pur, Float_t gate) |
f0a4c9e9 |
389 | { |
f0a4c9e9 |
390 | // Set the parameter TimeGate depending on the cluster energy and |
bc0c084c |
391 | // Purity-Efficiency point |
392 | if(Eff_Pur>2 || Eff_Pur<0) |
393 | Error("SetTimeGate","Invalid Efficiency-Purity choice %d",Eff_Pur); |
394 | |
395 | (*fParameters)(2,Eff_Pur)= gate ; |
f0a4c9e9 |
396 | } |
a496c46c |
397 | //_____________________________________________________________________________ |
50739f15 |
398 | void AliPHOSPIDv1::SetParameters() |
a496c46c |
399 | { |
400 | // PCA : To do the Principal Components Analysis it is necessary |
401 | // the Principal file, which is opened here |
402 | fX = new double[7]; // Data for the PCA |
f0a4c9e9 |
403 | fP = new double[7]; // Eigenvalues of the PCA |
e0ed2e49 |
404 | |
bc0c084c |
405 | // Open principal file |
50739f15 |
406 | |
bc0c084c |
407 | fFileName = "$ALICE_ROOT/PHOS/PCA8pa15_0.5-100.root" ; |
50739f15 |
408 | TFile f( fFileName.Data(), "read" ) ; |
a496c46c |
409 | fPrincipal = dynamic_cast<TPrincipal*> (f.Get("principal")) ; |
410 | f.Close() ; |
a496c46c |
411 | |
bc0c084c |
412 | // Open parameters file and initialization of the Parameters matrix. |
413 | // In the File Parameters.dat are all the parameters. These are introduced |
414 | // in a matrix of 9x4 |
50739f15 |
415 | // |
416 | // All the parameters defined in this file are, in order of row: |
bc0c084c |
417 | // -CpvtoEmcDistanceCut (2 row (x and z) and 3 columns, each one depending |
418 | // on the parameter of the funtion that sets the cut in x or z. |
419 | // -TimeGate, 1 row and 3 columns (3 efficiency-purty cuts) |
420 | // -PCA, parameters of the functions that |
50739f15 |
421 | // calculate the ellipse parameters, x0,y0,a, b, c. These 5 parameters |
bc0c084c |
422 | // (5 rows) depend on 4 parameters (columns). |
423 | // -Finally there is a row with the energy calibration parameters, |
424 | // 3 parameters. |
425 | |
426 | fFileNamePar = gSystem->ExpandPathName("$ALICE_ROOT/PHOS/Parameters.dat"); |
427 | fParameters = new TMatrixD(9,4) ; |
9fa5f1d0 |
428 | |
bc0c084c |
429 | FILE *par = fopen(fFileNamePar,"r"); |
430 | for(int i = 0;i<9;i++){ |
431 | fscanf(par, "%lf %lf %lf %lf", &(*fParameters)(i,0), |
432 | &(*fParameters)(i,1), |
433 | &(*fParameters)(i,2), &(*fParameters)(i,3)); |
a496c46c |
434 | } |
bc0c084c |
435 | fclose(par); |
436 | //fParameters->Print(); |
f0a4c9e9 |
437 | } |
f0a4c9e9 |
438 | |
f0a4c9e9 |
439 | |
50739f15 |
440 | //________________________________________________________________________ |
441 | void AliPHOSPIDv1::SetEllipseParameter(TString Param, Int_t i, Double_t par) |
442 | { |
443 | // Set the parameter "i" that is needed to calculate the ellipse |
444 | // parameter "Param". |
bc0c084c |
445 | |
50739f15 |
446 | Int_t p= -1; |
447 | |
bc0c084c |
448 | if(Param.Contains("a")) p=3; |
449 | if(Param.Contains("b")) p=4; |
450 | if(Param.Contains("c")) p=5; |
451 | if(Param.Contains("x0"))p=6; |
452 | if(Param.Contains("y0"))p=7; |
50739f15 |
453 | if((i>4)||(i<0)) |
21cd0c07 |
454 | Error("SetEllipseParameter", "No parameter with index %d", i) ; |
50739f15 |
455 | else if(p==-1) |
21cd0c07 |
456 | Error("SetEllipseParameter", "No parameter with name %s", Param.Data() ) ; |
50739f15 |
457 | else |
458 | (*fParameters)(p,i) = par ; |
459 | } |
460 | //________________________________________________________________________ |
461 | const Double_t AliPHOSPIDv1::GetParameterToCalculateEllipse(const TString Param, const Int_t i) const |
462 | { |
463 | // Get the parameter "i" that is needed to calculate the ellipse |
464 | // parameter "Param". |
465 | |
466 | Int_t p= -1; |
467 | Double_t par = -1; |
468 | |
bc0c084c |
469 | if(Param.Contains("a")) p=3; |
470 | if(Param.Contains("b")) p=4; |
471 | if(Param.Contains("c")) p=5; |
472 | if(Param.Contains("x0"))p=6; |
473 | if(Param.Contains("y0"))p=7; |
50739f15 |
474 | |
475 | if((i>4)||(i<0)) |
21cd0c07 |
476 | Error("GetParameterToCalculateEllipse", "No parameter with index", i) ; |
50739f15 |
477 | else if(p==-1) |
21cd0c07 |
478 | Error("GetParameterToCalculateEllipse", "No parameter with name %s", Param.Data() ) ; |
50739f15 |
479 | else |
480 | par = (*fParameters)(p,i) ; |
481 | |
482 | return par; |
483 | |
484 | } |
485 | //____________________________________________________________________________ |
486 | void AliPHOSPIDv1::SetCalibrationParameter(Int_t i,Double_t param) |
487 | { |
488 | (*fParameters)(17,i) = param ; |
489 | } |
490 | //____________________________________________________________________________ |
491 | const Double_t AliPHOSPIDv1::GetCalibrationParameter(const Int_t i) const |
492 | { |
bc0c084c |
493 | Float_t param = (*fParameters)(8,i); |
50739f15 |
494 | return param; |
495 | } |
496 | //____________________________________________________________________________ |
497 | const Double_t AliPHOSPIDv1::GetEllipseParameter(const TString Param,Float_t E) const |
498 | { |
bc0c084c |
499 | // Calculates the parameter Param of the ellipse |
500 | |
50739f15 |
501 | Double_t p[4]={0.,0.,0.,0.}; |
502 | Double_t value = 0.0; |
503 | Int_t i; |
504 | |
505 | if(Param.Contains("a")){ |
bc0c084c |
506 | for(i=0;i<4;i++)p[i]=(*fParameters)(3,i); |
50739f15 |
507 | if(E>70.)E=70.; |
508 | } |
509 | |
510 | if(Param.Contains("b")){ |
bc0c084c |
511 | for(i=0;i<4;i++)p[i]=(*fParameters)(4,i); |
50739f15 |
512 | if(E>70.)E=70.; |
513 | } |
514 | |
515 | if(Param.Contains("c")) |
bc0c084c |
516 | for(i=0;i<4;i++)p[i]=(*fParameters)(5,i); |
50739f15 |
517 | |
518 | if(Param.Contains("x0")){ |
bc0c084c |
519 | for(i=0;i<4;i++)p[i]=(*fParameters)(6,i); |
50739f15 |
520 | if(E<1.)E=1.1; |
521 | } |
522 | if(Param.Contains("y0")) |
bc0c084c |
523 | for(i=0;i<4;i++)p[i]=(*fParameters)(7,i); |
50739f15 |
524 | |
525 | value = p[0]/TMath::Sqrt(E)+p[1]*E+p[2]*E*E+p[3]; |
526 | return value; |
527 | } |
148b2bba |
528 | //____________________________________________________________________________ |
529 | |
7acf6008 |
530 | void AliPHOSPIDv1::Exec(Option_t * option) |
6ad0bfa0 |
531 | { |
f035f6ce |
532 | //Steering method |
9688c1dd |
533 | |
bf8f1fbd |
534 | if( strcmp(GetName(), "")== 0 ) |
7acf6008 |
535 | Init() ; |
bf8f1fbd |
536 | |
537 | if(strstr(option,"tim")) |
7acf6008 |
538 | gBenchmark->Start("PHOSPID"); |
bf8f1fbd |
539 | |
540 | if(strstr(option,"print")) { |
7b7c1533 |
541 | Print("") ; |
542 | return ; |
bf8f1fbd |
543 | } |
9688c1dd |
544 | |
148b2bba |
545 | |
fbf811ec |
546 | // gAlice->GetEvent(0) ; |
148b2bba |
547 | |
fbf811ec |
548 | // //check, if the branch with name of this" already exits? |
549 | // if (gAlice->TreeR()) { |
550 | // TObjArray * lob = (TObjArray*)gAlice->TreeR()->GetListOfBranches() ; |
551 | // TIter next(lob) ; |
552 | // TBranch * branch = 0 ; |
553 | // Bool_t phospidfound = kFALSE, pidfound = kFALSE ; |
8d0f3f77 |
554 | |
fbf811ec |
555 | // TString taskName(GetName()) ; |
556 | // taskName.Remove(taskName.Index(Version())-1) ; |
8d0f3f77 |
557 | |
fbf811ec |
558 | // while ( (branch = (TBranch*)next()) && (!phospidfound || !pidfound) ) { |
559 | // if ( (strcmp(branch->GetName(), "PHOSPID")==0) && (strcmp(branch->GetTitle(), taskName.Data())==0) ) |
560 | // phospidfound = kTRUE ; |
8d0f3f77 |
561 | |
fbf811ec |
562 | // else if ( (strcmp(branch->GetName(), "AliPHOSPID")==0) && (strcmp(branch->GetTitle(), taskName.Data())==0) ) |
563 | // pidfound = kTRUE ; |
564 | // } |
7b7c1533 |
565 | |
fbf811ec |
566 | // if ( phospidfound || pidfound ) { |
21cd0c07 |
567 | // Error("Exec", "RecParticles and/or PIDtMaker branch with name %s already exists", taskName.Data() ) ; |
fbf811ec |
568 | // return ; |
569 | // } |
570 | // } |
571 | |
572 | // Int_t nevents = (Int_t) gAlice->TreeE()->GetEntries() ; |
573 | // Int_t ievent ; |
574 | // AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ; |
9688c1dd |
575 | |
fbf811ec |
576 | AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ; |
577 | if(gime->BranchExists("RecParticles") ) |
578 | return ; |
579 | Int_t nevents = gime->MaxEvent() ; //(Int_t) gAlice->TreeE()->GetEntries() ; |
7b7c1533 |
580 | Int_t ievent ; |
fbf811ec |
581 | |
582 | |
7b7c1533 |
583 | for(ievent = 0; ievent < nevents; ievent++){ |
a496c46c |
584 | gime->Event(ievent,"R") ; |
148b2bba |
585 | |
7acf6008 |
586 | MakeRecParticles() ; |
9688c1dd |
587 | |
7b7c1533 |
588 | WriteRecParticles(ievent); |
9688c1dd |
589 | |
7acf6008 |
590 | if(strstr(option,"deb")) |
591 | PrintRecParticles(option) ; |
94de8339 |
592 | |
593 | //increment the total number of rec particles per run |
a496c46c |
594 | fRecParticlesInRun += gime->RecParticles(BranchName())->GetEntriesFast() ; |
94de8339 |
595 | |
7acf6008 |
596 | } |
9688c1dd |
597 | |
7acf6008 |
598 | if(strstr(option,"tim")){ |
599 | gBenchmark->Stop("PHOSPID"); |
21cd0c07 |
600 | Info("Exec", "took %f seconds for PID %f seconds per event", |
601 | gBenchmark->GetCpuTime("PHOSPID"), |
602 | gBenchmark->GetCpuTime("PHOSPID")/nevents) ; |
603 | } |
7b7c1533 |
604 | } |
605 | |
7acf6008 |
606 | //____________________________________________________________________________ |
607 | void AliPHOSPIDv1::MakeRecParticles(){ |
608 | |
b2a60966 |
609 | // Makes a RecParticle out of a TrackSegment |
148b2bba |
610 | |
7b7c1533 |
611 | AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ; |
fbf811ec |
612 | TObjArray * emcRecPoints = gime->EmcRecPoints() ; |
613 | TObjArray * cpvRecPoints = gime->CpvRecPoints() ; |
614 | TClonesArray * trackSegments = gime->TrackSegments() ; |
148b2bba |
615 | if ( !emcRecPoints || !cpvRecPoints || !trackSegments ) { |
f1611b7c |
616 | Fatal("MakeRecParticles", "RecPoints or TrackSegments not found !") ; |
148b2bba |
617 | } |
fbf811ec |
618 | TClonesArray * recParticles = gime->RecParticles() ; |
01a599c9 |
619 | recParticles->Clear(); |
148b2bba |
620 | |
7b7c1533 |
621 | TIter next(trackSegments) ; |
7acf6008 |
622 | AliPHOSTrackSegment * ts ; |
6ad0bfa0 |
623 | Int_t index = 0 ; |
09fc14a0 |
624 | AliPHOSRecParticle * rp ; |
7acf6008 |
625 | while ( (ts = (AliPHOSTrackSegment *)next()) ) { |
fad3e5b9 |
626 | |
7b7c1533 |
627 | new( (*recParticles)[index] ) AliPHOSRecParticle() ; |
628 | rp = (AliPHOSRecParticle *)recParticles->At(index) ; |
f0a4c9e9 |
629 | rp->SetTrackSegment(index) ; |
9688c1dd |
630 | rp->SetIndexInList(index) ; |
148b2bba |
631 | |
7acf6008 |
632 | AliPHOSEmcRecPoint * emc = 0 ; |
633 | if(ts->GetEmcIndex()>=0) |
7b7c1533 |
634 | emc = (AliPHOSEmcRecPoint *) emcRecPoints->At(ts->GetEmcIndex()) ; |
fad3e5b9 |
635 | |
7acf6008 |
636 | AliPHOSRecPoint * cpv = 0 ; |
637 | if(ts->GetCpvIndex()>=0) |
7b7c1533 |
638 | cpv = (AliPHOSRecPoint *) cpvRecPoints->At(ts->GetCpvIndex()) ; |
fad3e5b9 |
639 | |
148b2bba |
640 | // Now set type (reconstructed) of the particle |
641 | |
642 | // Choose the cluster energy range |
9fa5f1d0 |
643 | |
fbf811ec |
644 | // YK: check if (emc != 0) !!! |
645 | if (!emc) { |
f1611b7c |
646 | Fatal("MakeRecParticles", "-> emc(%d) = %d", ts->GetEmcIndex(), emc ) ; |
fbf811ec |
647 | } |
50739f15 |
648 | |
bc0c084c |
649 | Float_t e = emc->GetEnergy() ; |
650 | |
6f969528 |
651 | Float_t lambda[2] ; |
652 | emc->GetElipsAxis(lambda) ; |
50739f15 |
653 | |
654 | if((lambda[0]>0.01) && (lambda[1]>0.01)){ |
655 | // Looking PCA. Define and calculate the data (X), |
bc0c084c |
656 | // introduce in the function X2P that gives the components (P). |
657 | |
50739f15 |
658 | Float_t Spher = 0. ; |
659 | Float_t Emaxdtotal = 0. ; |
660 | |
bc0c084c |
661 | if((lambda[0]+lambda[1])!=0) |
662 | Spher=fabs(lambda[0]-lambda[1])/(lambda[0]+lambda[1]); |
50739f15 |
663 | |
664 | Emaxdtotal=emc->GetMaximalEnergy()/emc->GetEnergy(); |
665 | |
666 | fX[0] = lambda[0] ; |
667 | fX[1] = lambda[1] ; |
668 | fX[2] = emc->GetDispersion() ; |
669 | fX[3] = Spher ; |
670 | fX[4] = emc->GetMultiplicity() ; |
671 | fX[5] = Emaxdtotal ; |
672 | fX[6] = emc->GetCoreEnergy() ; |
673 | |
674 | fPrincipal->X2P(fX,fP); |
675 | } |
676 | else{ |
677 | fP[0]=-100.0; //We do not accept clusters with |
678 | fP[1]=-100.0; //one cell as a photon-like |
679 | } |
680 | |
6f969528 |
681 | Float_t time =emc->GetTime() ; |
9fa5f1d0 |
682 | |
bc0c084c |
683 | // Loop of Efficiency-Purity (the 3 points of purity or efficiency |
684 | // are taken into account to set the particle identification) |
50739f15 |
685 | for(Int_t eff_pur = 0; eff_pur < 3 ; eff_pur++){ |
686 | |
bc0c084c |
687 | // Looking at the CPV detector. If RCPV greater than CpvEmcDistance, |
688 | // 1st,2nd or 3rd bit (depending on the efficiency-purity point ) |
689 | // is set to 1 |
690 | if(GetCPVBit(emc, cpv, eff_pur,e) == 1 ) |
50739f15 |
691 | rp->SetPIDBit(eff_pur) ; |
f0a4c9e9 |
692 | |
50739f15 |
693 | // Looking the TOF. If TOF smaller than gate, 4th, 5th or 6th |
694 | // bit (depending on the efficiency-purity point )is set to 1 |
bc0c084c |
695 | if(time< (*fParameters)(2,eff_pur)) { |
50739f15 |
696 | rp->SetPIDBit(eff_pur+3) ; |
9fa5f1d0 |
697 | } |
50739f15 |
698 | |
699 | //If we are inside the ellipse, 7th, 8th or 9th |
700 | // bit (depending on the efficiency-purity point )is set to 1 |
bc0c084c |
701 | if(GetPrincipalBit(fP,eff_pur,e) == 1) |
50739f15 |
702 | rp->SetPIDBit(eff_pur+6) ; |
f0a4c9e9 |
703 | } |
50739f15 |
704 | |
9fa5f1d0 |
705 | //Set momentum, energy and other parameters |
50739f15 |
706 | Float_t encal = GetCalibratedEnergy(e); |
9fa5f1d0 |
707 | TVector3 dir = GetMomentumDirection(emc,cpv) ; |
708 | dir.SetMag(encal) ; |
709 | rp->SetMomentum(dir.X(),dir.Y(),dir.Z(),encal) ; |
710 | rp->SetCalcMass(0); |
e0ed2e49 |
711 | rp->Name(); //If photon sets the particle pdg name to gamma |
e747b8da |
712 | rp->SetProductionVertex(0,0,0,0); |
713 | rp->SetFirstMother(-1); |
714 | rp->SetLastMother(-1); |
715 | rp->SetFirstDaughter(-1); |
716 | rp->SetLastDaughter(-1); |
717 | rp->SetPolarisation(0,0,0); |
6ad0bfa0 |
718 | index++ ; |
719 | } |
7acf6008 |
720 | |
6ad0bfa0 |
721 | } |
722 | |
09fc14a0 |
723 | //____________________________________________________________________________ |
a496c46c |
724 | void AliPHOSPIDv1:: Print() |
09fc14a0 |
725 | { |
b2a60966 |
726 | // Print the parameters used for the particle type identification |
bc0c084c |
727 | |
21cd0c07 |
728 | TString message ; |
729 | message = "\n=============== AliPHOSPID1 ================\n" ; |
730 | message += "Making PID\n"; |
731 | message += " Pricipal analysis file from 0.5 to 100 %s\n" ; |
732 | message += " Name of parameters file %s\n" ; |
bc0c084c |
733 | message += " Matrix of Parameters: 9x4\n" ; |
734 | message += " RCPV 2x3 rows x and z, columns function cut parameters\n" ; |
735 | message += " TOF 1x3 [High Eff-Low Pur,Medium Eff-Pur, Low Eff-High Pur]\n" ; |
21cd0c07 |
736 | message += " PCA 5x4 [5 ellipse parametres and 4 parametres to calculate them: A/Sqrt(E) + B* E + C * E^2 + D]\n" ; |
737 | message += " Energy Calibration 1x3 [3 parametres to calibrate energy: A + B* E + C * E^2]\n" ; |
738 | Info("Print", message.Data(), fFileName.Data(), fFileNamePar.Data() ) ; |
50739f15 |
739 | fParameters->Print() ; |
09fc14a0 |
740 | } |
741 | |
7acf6008 |
742 | //____________________________________________________________________________ |
7b7c1533 |
743 | void AliPHOSPIDv1::WriteRecParticles(Int_t event) |
09fc14a0 |
744 | { |
7b7c1533 |
745 | |
746 | AliPHOSGetter *gime = AliPHOSGetter::GetInstance() ; |
a496c46c |
747 | |
fbf811ec |
748 | TClonesArray * recParticles = gime->RecParticles() ; |
bf8f1fbd |
749 | recParticles->Expand(recParticles->GetEntriesFast() ) ; |
fbf811ec |
750 | TTree * treeR ; |
751 | |
752 | if(fToSplit){ |
753 | if(!fSplitFile) |
754 | return ; |
755 | fSplitFile->cd() ; |
756 | char name[10] ; |
757 | sprintf(name,"%s%d", "TreeR",event) ; |
758 | treeR = dynamic_cast<TTree*>(fSplitFile->Get(name)); |
759 | } |
760 | else{ |
761 | treeR = gAlice->TreeR(); |
762 | } |
763 | |
764 | if(!treeR){ |
765 | gAlice->MakeTree("R", fSplitFile); |
766 | treeR = gAlice->TreeR() ; |
767 | } |
7acf6008 |
768 | |
769 | //First rp |
770 | Int_t bufferSize = 32000 ; |
8d0f3f77 |
771 | TBranch * rpBranch = treeR->Branch("PHOSRP",&recParticles,bufferSize); |
fbf811ec |
772 | rpBranch->SetTitle(BranchName()); |
8d0f3f77 |
773 | |
9688c1dd |
774 | |
7acf6008 |
775 | //second, pid |
776 | Int_t splitlevel = 0 ; |
777 | AliPHOSPIDv1 * pid = this ; |
8d0f3f77 |
778 | TBranch * pidBranch = treeR->Branch("AliPHOSPID","AliPHOSPIDv1",&pid,bufferSize,splitlevel); |
fbf811ec |
779 | pidBranch->SetTitle(BranchName()); |
7acf6008 |
780 | |
761e34c0 |
781 | rpBranch->Fill() ; |
a496c46c |
782 | pidBranch->Fill() ; |
9688c1dd |
783 | |
fbf811ec |
784 | treeR->AutoSave() ; //Write(0,kOverwrite) ; |
785 | if(gAlice->TreeR()!=treeR){ |
786 | treeR->Delete(); |
787 | } |
7acf6008 |
788 | } |
69183710 |
789 | |
7acf6008 |
790 | //____________________________________________________________________________ |
9688c1dd |
791 | TVector3 AliPHOSPIDv1::GetMomentumDirection(AliPHOSEmcRecPoint * emc, AliPHOSRecPoint * cpv)const |
7acf6008 |
792 | { |
793 | // Calculates the momentum direction: |
794 | // 1. if only a EMC RecPoint, direction is given by IP and this RecPoint |
9688c1dd |
795 | // 2. if a EMC RecPoint and CPV RecPoint, direction is given by the line through the 2 recpoints |
796 | // However because of the poor position resolution of PPSD the direction is always taken as if we were |
7acf6008 |
797 | // in case 1. |
798 | |
799 | TVector3 dir(0,0,0) ; |
800 | |
801 | TVector3 emcglobalpos ; |
802 | TMatrix dummy ; |
803 | |
804 | emc->GetGlobalPosition(emcglobalpos, dummy) ; |
805 | |
148b2bba |
806 | |
7acf6008 |
807 | dir = emcglobalpos ; |
808 | dir.SetZ( -dir.Z() ) ; // why ? |
809 | dir.SetMag(1.) ; |
810 | |
811 | //account correction to the position of IP |
812 | Float_t xo,yo,zo ; //Coordinates of the origin |
813 | gAlice->Generator()->GetOrigin(xo,yo,zo) ; |
814 | TVector3 origin(xo,yo,zo); |
815 | dir = dir - origin ; |
816 | |
817 | return dir ; |
818 | } |
819 | //____________________________________________________________________________ |
a4e98857 |
820 | void AliPHOSPIDv1::PrintRecParticles(Option_t * option) |
821 | { |
dd5c4038 |
822 | // Print table of reconstructed particles |
823 | |
7b7c1533 |
824 | AliPHOSGetter *gime = AliPHOSGetter::GetInstance() ; |
bf8f1fbd |
825 | |
a496c46c |
826 | TClonesArray * recParticles = gime->RecParticles(BranchName()) ; |
21cd0c07 |
827 | |
828 | TString message ; |
3bf72d32 |
829 | message = "\nevent " ; |
830 | message += gAlice->GetEvNumber() ; |
831 | message += " found " ; |
832 | message += recParticles->GetEntriesFast(); |
833 | message += " RecParticles\n" ; |
834 | |
7acf6008 |
835 | if(strstr(option,"all")) { // printing found TS |
3bf72d32 |
836 | message += "\n PARTICLE Index \n" ; |
7acf6008 |
837 | |
838 | Int_t index ; |
7b7c1533 |
839 | for (index = 0 ; index < recParticles->GetEntries() ; index++) { |
21cd0c07 |
840 | AliPHOSRecParticle * rp = (AliPHOSRecParticle * ) recParticles->At(index) ; |
3bf72d32 |
841 | message += "\n" ; |
842 | message += rp->Name().Data() ; |
843 | message += " " ; |
844 | message += rp->GetIndexInList() ; |
845 | message += " " ; |
846 | message += rp->GetType() ; |
7acf6008 |
847 | } |
3bf72d32 |
848 | } |
849 | Info("Print", message.Data() ) ; |
69183710 |
850 | } |
851 | |
7acf6008 |
852 | |
853 | |