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