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