]>
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" | |
e3817e5f | 85 | #include "TF2.h" |
86 | #include "TFormula.h" | |
87 | #include "TCanvas.h" | |
88 | #include "TFolder.h" | |
7acf6008 | 89 | #include "TSystem.h" |
90 | #include "TBenchmark.h" | |
148b2bba | 91 | #include "TMatrixD.h" |
92 | #include "TPrincipal.h" | |
e3817e5f | 93 | #include "TSystem.h" |
148b2bba | 94 | |
6ad0bfa0 | 95 | // --- Standard library --- |
96 | ||
75a6835b | 97 | |
6ad0bfa0 | 98 | // --- AliRoot header files --- |
99 | ||
7acf6008 | 100 | #include "AliGenerator.h" |
e3817e5f | 101 | #include "AliPHOS.h" |
26d4b141 | 102 | #include "AliPHOSPIDv1.h" |
e3817e5f | 103 | #include "AliPHOSClusterizerv1.h" |
d956e9b7 | 104 | #include "AliPHOSEmcRecPoint.h" |
6ad0bfa0 | 105 | #include "AliPHOSTrackSegment.h" |
e3817e5f | 106 | #include "AliPHOSTrackSegmentMakerv1.h" |
6ad0bfa0 | 107 | #include "AliPHOSRecParticle.h" |
7b7c1533 | 108 | #include "AliPHOSGeometry.h" |
109 | #include "AliPHOSGetter.h" | |
6ad0bfa0 | 110 | |
26d4b141 | 111 | ClassImp( AliPHOSPIDv1) |
6ad0bfa0 | 112 | |
1cb7c1ee | 113 | //____________________________________________________________________________ |
114 | AliPHOSPIDv1::AliPHOSPIDv1():AliPHOSPID() | |
115 | { | |
a4e98857 | 116 | // default ctor |
148b2bba | 117 | |
8d0f3f77 | 118 | InitParameters() ; |
92f521a9 | 119 | fDefaultInit = kTRUE ; |
7acf6008 | 120 | } |
121 | ||
581354c5 | 122 | //____________________________________________________________________________ |
88cb7938 | 123 | AliPHOSPIDv1::AliPHOSPIDv1(const AliPHOSPIDv1 & pid ):AliPHOSPID(pid) |
581354c5 | 124 | { |
386aef34 | 125 | // ctor |
581354c5 | 126 | InitParameters() ; |
581354c5 | 127 | Init() ; |
581354c5 | 128 | |
129 | } | |
130 | ||
7acf6008 | 131 | //____________________________________________________________________________ |
88cb7938 | 132 | AliPHOSPIDv1::AliPHOSPIDv1(const TString alirunFileName, const TString eventFolderName):AliPHOSPID(alirunFileName, eventFolderName) |
7acf6008 | 133 | { |
a4e98857 | 134 | //ctor with the indication on where to look for the track segments |
7b7c1533 | 135 | |
8d0f3f77 | 136 | InitParameters() ; |
2bd5457f | 137 | Init() ; |
92f521a9 | 138 | fDefaultInit = kFALSE ; |
7acf6008 | 139 | } |
7b7c1533 | 140 | |
7acf6008 | 141 | //____________________________________________________________________________ |
142 | AliPHOSPIDv1::~AliPHOSPIDv1() | |
143 | { | |
79bb1b62 | 144 | // dtor |
9fa5f1d0 | 145 | |
e3817e5f | 146 | delete [] fX ; // Principal input |
147 | delete [] fPPhoton ; // Photon Principal components | |
148 | delete [] fPPi0 ; // Pi0 Principal components | |
7acf6008 | 149 | } |
a496c46c | 150 | //____________________________________________________________________________ |
151 | const TString AliPHOSPIDv1::BranchName() const | |
152 | { | |
88cb7938 | 153 | |
154 | return GetName() ; | |
a496c46c | 155 | } |
156 | ||
148b2bba | 157 | //____________________________________________________________________________ |
158 | void AliPHOSPIDv1::Init() | |
159 | { | |
160 | // Make all memory allocations that are not possible in default constructor | |
161 | // Add the PID task to the list of PHOS tasks | |
a496c46c | 162 | |
88cb7938 | 163 | AliPHOSGetter * gime = AliPHOSGetter::Instance(GetTitle(), fEventFolderName.Data()) ; |
164 | ||
165 | if ( !gime->PID() ) | |
166 | gime->PostPID(this) ; | |
148b2bba | 167 | } |
8d0f3f77 | 168 | |
169 | //____________________________________________________________________________ | |
170 | void AliPHOSPIDv1::InitParameters() | |
171 | { | |
e3817e5f | 172 | // Initialize PID parameters |
8d0f3f77 | 173 | fRecParticlesInRun = 0 ; |
8d0f3f77 | 174 | fNEvent = 0 ; |
8d0f3f77 | 175 | fRecParticlesInRun = 0 ; |
9fa5f1d0 | 176 | SetParameters() ; // fill the parameters matrix from parameters file |
eabde521 | 177 | SetEventRange(0,-1) ; |
cc1fe362 | 178 | // initialisation of response function parameters |
179 | // Tof | |
180 | // Photons | |
181 | fTphoton[0] = 2.97E-1 ; | |
182 | fTphoton[1] = 1.55E-8 ; | |
183 | fTphoton[2] = 5.40E-10 ; | |
184 | fTFphoton = new TFormula("ToF response to photon" , "gaus") ; | |
185 | fTFphoton->SetParameters( fTphoton[0], fTphoton[1], fTphoton[2]) ; | |
186 | // Electrons | |
187 | fTelectron[0] = 2.73E-1 ; | |
188 | fTelectron[1] = 1.56E-8 ; | |
189 | fTelectron[2] = 5.70E-10 ; | |
190 | fTFelectron = new TFormula("ToF response to electron" , "gaus") ; | |
191 | fTFelectron->SetParameters( fTelectron[0], fTelectron[1], fTelectron[2]) ; | |
192 | //Charged Hadrons | |
193 | fTchargedhadron[0] = 1.58E-1 ; | |
194 | fTchargedhadron[1] = 1.64E-8 ; | |
195 | fTchargedhadron[2] = 3.59E-10 ; | |
196 | fTFchargedhadron = new TFormula("ToF response to charged hadron" , "landau") ; | |
197 | fTFchargedhadron->SetParameters( fTchargedhadron[0], fTchargedhadron[1], fTchargedhadron[2]) ; | |
198 | //Neutral Hadrons | |
199 | fTneutralhadron[0] = 9.62E-1 ; | |
200 | fTneutralhadron[1] = 1.65E-8 ; | |
201 | fTneutralhadron[2] = 6.46E-10 ; | |
202 | fTFneutralhadron = new TFormula("ToF response to neutral hadron" , "landau") ; | |
203 | fTFneutralhadron->SetParameters( fTneutralhadron[0], fTneutralhadron[1], fTneutralhadron[2]) ; | |
8d0f3f77 | 204 | } |
205 | ||
88cb7938 | 206 | //________________________________________________________________________ |
eabde521 | 207 | void AliPHOSPIDv1::Exec(Option_t *option) |
88cb7938 | 208 | { |
eabde521 | 209 | // Steering method to perform particle reconstruction and identification |
210 | // for the event range from fFirstEvent to fLastEvent. | |
211 | // This range is optionally set by SetEventRange(). | |
212 | // if fLastEvent=-1 (by default), then process events until the end. | |
61d3d6aa | 213 | |
88cb7938 | 214 | if(strstr(option,"tim")) |
215 | gBenchmark->Start("PHOSPID"); | |
216 | ||
217 | if(strstr(option,"print")) { | |
218 | Print() ; | |
219 | return ; | |
220 | } | |
221 | ||
222 | ||
fb43ada4 | 223 | AliPHOSGetter * gime = AliPHOSGetter::Instance(GetTitle()) ; |
88cb7938 | 224 | |
71cee46d | 225 | if (fLastEvent == -1) |
226 | fLastEvent = gime->MaxEvent() - 1 ; | |
227 | else | |
228 | fLastEvent = TMath::Min(fLastEvent,gime->MaxEvent()); | |
eabde521 | 229 | Int_t nEvents = fLastEvent - fFirstEvent + 1; |
88cb7938 | 230 | |
71cee46d | 231 | Int_t ievent ; |
232 | for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) { | |
88cb7938 | 233 | gime->Event(ievent,"TR") ; |
234 | if(gime->TrackSegments() && //Skip events, where no track segments made | |
235 | gime->TrackSegments()->GetEntriesFast()) { | |
236 | MakeRecParticles() ; | |
2cc71c1e | 237 | MakePID() ; |
90cceaf6 | 238 | WriteRecParticles(); |
88cb7938 | 239 | if(strstr(option,"deb")) |
240 | PrintRecParticles(option) ; | |
241 | //increment the total number of rec particles per run | |
242 | fRecParticlesInRun += gime->RecParticles()->GetEntriesFast() ; | |
243 | } | |
244 | } | |
ff417097 | 245 | if(strstr(option,"deb")) |
246 | PrintRecParticles(option); | |
88cb7938 | 247 | if(strstr(option,"tim")){ |
248 | gBenchmark->Stop("PHOSPID"); | |
249 | Info("Exec", "took %f seconds for PID %f seconds per event", | |
250 | gBenchmark->GetCpuTime("PHOSPID"), | |
eabde521 | 251 | gBenchmark->GetCpuTime("PHOSPID")/nEvents) ; |
88cb7938 | 252 | } |
88cb7938 | 253 | Unload(); |
254 | } | |
255 | ||
69183710 | 256 | //____________________________________________________________________________ |
e3817e5f | 257 | const TString AliPHOSPIDv1::GetFileNamePrincipal(TString particle) const |
148b2bba | 258 | { |
e3817e5f | 259 | //Get file name that contains the PCA for a particle ("photon or pi0") |
260 | particle.ToLower(); | |
261 | TString name; | |
262 | if (particle=="photon") name = fFileNamePrincipalPhoton ; | |
263 | else if (particle=="pi0" ) name = fFileNamePrincipalPi0 ; | |
264 | else Error("GetFileNamePrincipal","Wrong particle name: %s (choose from pi0/photon)\n",particle.Data()); | |
265 | return name; | |
266 | } | |
bc0c084c | 267 | |
e3817e5f | 268 | //____________________________________________________________________________ |
fc7e2f43 | 269 | Float_t AliPHOSPIDv1::GetParameterCalibration(Int_t i) const |
e3817e5f | 270 | { |
271 | // Get the i-th parameter "Calibration" | |
272 | Float_t param = 0.; | |
273 | if (i>2 || i<0) | |
274 | Error("GetParameterCalibration","Invalid parameter number: %d",i); | |
9fa5f1d0 | 275 | else |
e3817e5f | 276 | param = (*fParameters)(0,i); |
277 | return param; | |
278 | } | |
bc0c084c | 279 | |
88cb7938 | 280 | //____________________________________________________________________________ |
fc7e2f43 | 281 | Float_t AliPHOSPIDv1::GetCalibratedEnergy(Float_t e) const |
88cb7938 | 282 | { |
283 | // It calibrates Energy depending on the recpoint energy. | |
284 | // The energy of the reconstructed cluster is corrected with | |
285 | // the formula A + B* E + C* E^2, whose parameters where obtained | |
286 | // through the study of the reconstructed energy distribution of | |
287 | // monoenergetic photons. | |
288 | ||
289 | Float_t p[]={0.,0.,0.}; | |
290 | for (Int_t i=0; i<3; i++) p[i] = GetParameterCalibration(i); | |
291 | Float_t enerec = p[0] + p[1]*e + p[2]*e*e; | |
292 | return enerec ; | |
293 | ||
294 | } | |
295 | ||
e3817e5f | 296 | //____________________________________________________________________________ |
fc7e2f43 | 297 | Float_t AliPHOSPIDv1::GetParameterCpv2Emc(Int_t i, TString axis) const |
e3817e5f | 298 | { |
299 | // Get the i-th parameter "CPV-EMC distance" for the specified axis | |
300 | Float_t param = 0.; | |
301 | if(i>2 || i<0) | |
302 | Error("GetParameterCpv2Emc","Invalid parameter number: %d",i); | |
303 | else { | |
304 | axis.ToLower(); | |
305 | if (axis == "x") param = (*fParameters)(1,i); | |
306 | else if (axis == "z") param = (*fParameters)(2,i); | |
307 | else Error("GetParameterCpv2Emc","Invalid axis name: %s",axis.Data()); | |
308 | } | |
309 | return param; | |
310 | } | |
311 | ||
312 | //____________________________________________________________________________ | |
fc7e2f43 | 313 | Float_t AliPHOSPIDv1::GetCpv2EmcDistanceCut(TString axis, Float_t e) const |
e3817e5f | 314 | { |
88cb7938 | 315 | // Get CpvtoEmcDistance Cut depending on the cluster energy, axis and |
316 | // Purity-Efficiency point | |
317 | ||
318 | axis.ToLower(); | |
319 | Float_t p[]={0.,0.,0.}; | |
320 | for (Int_t i=0; i<3; i++) p[i] = GetParameterCpv2Emc(i,axis); | |
321 | Float_t sig = p[0] + TMath::Exp(p[1] - p[2]*e); | |
322 | return sig; | |
e3817e5f | 323 | } |
324 | ||
88cb7938 | 325 | //____________________________________________________________________________ |
fc7e2f43 | 326 | Float_t AliPHOSPIDv1::GetEllipseParameter(TString particle, TString param, Float_t e) const |
88cb7938 | 327 | { |
328 | // Calculates the parameter param of the ellipse | |
e3817e5f | 329 | |
330 | particle.ToLower(); | |
331 | param. ToLower(); | |
88cb7938 | 332 | Float_t p[4]={0.,0.,0.,0.}; |
333 | Float_t value = 0.0; | |
334 | for (Int_t i=0; i<4; i++) p[i] = GetParameterToCalculateEllipse(particle,param,i); | |
335 | if (particle == "photon") { | |
336 | if (param.Contains("a")) e = TMath::Min((Double_t)e,70.); | |
337 | else if (param.Contains("b")) e = TMath::Min((Double_t)e,70.); | |
338 | else if (param.Contains("x0")) e = TMath::Max((Double_t)e,1.1); | |
339 | } | |
e3817e5f | 340 | |
443caba9 | 341 | if (particle == "photon") |
342 | value = p[0]/TMath::Sqrt(e) + p[1]*e + p[2]*e*e + p[3]; | |
343 | else if (particle == "pi0") | |
344 | value = p[0] + p[1]*e + p[2]*e*e; | |
345 | ||
88cb7938 | 346 | return value; |
e3817e5f | 347 | } |
348 | ||
349 | //_____________________________________________________________________________ | |
fc7e2f43 | 350 | Float_t AliPHOSPIDv1::GetParameterPhotonBoundary (Int_t i) const |
e3817e5f | 351 | { |
352 | // Get the parameter "i" to calculate the boundary on the moment M2x | |
353 | // for photons at high p_T | |
354 | Float_t param = 0; | |
355 | if (i>3 || i<0) | |
356 | Error("GetParameterPhotonBoundary","Wrong parameter number: %d\n",i); | |
357 | else | |
358 | param = (*fParameters)(14,i) ; | |
359 | return param; | |
148b2bba | 360 | } |
e3817e5f | 361 | |
148b2bba | 362 | //____________________________________________________________________________ |
fc7e2f43 | 363 | Float_t AliPHOSPIDv1::GetParameterPi0Boundary (Int_t i) const |
e3817e5f | 364 | { |
365 | // Get the parameter "i" to calculate the boundary on the moment M2x | |
366 | // for pi0 at high p_T | |
367 | Float_t param = 0; | |
368 | if (i>2 || i<0) | |
369 | Error("GetParameterPi0Boundary","Wrong parameter number: %d\n",i); | |
370 | else | |
371 | param = (*fParameters)(15,i) ; | |
372 | return param; | |
373 | } | |
148b2bba | 374 | |
e3817e5f | 375 | //____________________________________________________________________________ |
fc7e2f43 | 376 | Float_t AliPHOSPIDv1::GetParameterTimeGate(Int_t i) const |
e3817e5f | 377 | { |
88cb7938 | 378 | // Get TimeGate parameter depending on Purity-Efficiency i: |
379 | // i=0 - Low purity, i=1 - Medium purity, i=2 - High purity | |
380 | Float_t param = 0.; | |
e3817e5f | 381 | if(i>2 || i<0) |
88cb7938 | 382 | Error("GetParameterTimeGate","Invalid Efficiency-Purity choice %d",i); |
e3817e5f | 383 | else |
88cb7938 | 384 | param = (*fParameters)(3,i) ; |
385 | return param; | |
e3817e5f | 386 | } |
387 | ||
148b2bba | 388 | //_____________________________________________________________________________ |
fc7e2f43 | 389 | Float_t AliPHOSPIDv1::GetParameterToCalculateEllipse(TString particle, TString param, Int_t i) const |
88cb7938 | 390 | { |
391 | // Get the parameter "i" that is needed to calculate the ellipse | |
392 | // parameter "param" for the particle "particle" ("photon" or "pi0") | |
393 | ||
e3817e5f | 394 | particle.ToLower(); |
395 | param. ToLower(); | |
88cb7938 | 396 | Int_t offset = -1; |
e3817e5f | 397 | if (particle == "photon") offset=0; |
398 | else if (particle == "pi0") offset=5; | |
399 | else | |
88cb7938 | 400 | Error("GetParameterToCalculateEllipse","Wrong particle name: %s (choose from pi0/photon)\n",particle.Data()); |
401 | ||
402 | Int_t p= -1; | |
403 | Float_t par = 0; | |
e3817e5f | 404 | |
405 | if (param.Contains("a")) p=4+offset; | |
406 | else if(param.Contains("b")) p=5+offset; | |
407 | else if(param.Contains("c")) p=6+offset; | |
408 | else if(param.Contains("x0"))p=7+offset; | |
409 | else if(param.Contains("y0"))p=8+offset; | |
12022e83 | 410 | |
88cb7938 | 411 | if (i>4 || i<0) |
412 | Error("GetParameterToCalculateEllipse", "No parameter with index", i) ; | |
413 | else if (p==-1) | |
414 | Error("GetParameterToCalculateEllipse", "No parameter with name %s", param.Data() ) ; | |
12022e83 | 415 | else |
88cb7938 | 416 | par = (*fParameters)(p,i) ; |
417 | ||
418 | return par; | |
12022e83 | 419 | } |
420 | ||
12022e83 | 421 | |
422 | //____________________________________________________________________________ | |
8d4608b5 | 423 | Float_t AliPHOSPIDv1::GetDistance(AliPHOSEmcRecPoint * emc,AliPHOSCpvRecPoint * cpv, Option_t * axis)const |
69183710 | 424 | { |
425 | // Calculates the distance between the EMC RecPoint and the PPSD RecPoint | |
148b2bba | 426 | |
88cb7938 | 427 | const AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ; |
69183710 | 428 | TVector3 vecEmc ; |
7acf6008 | 429 | TVector3 vecCpv ; |
148b2bba | 430 | if(cpv){ |
431 | emc->GetLocalPosition(vecEmc) ; | |
432 | cpv->GetLocalPosition(vecCpv) ; | |
2bb500e5 | 433 | |
148b2bba | 434 | if(emc->GetPHOSMod() == cpv->GetPHOSMod()){ |
435 | // Correct to difference in CPV and EMC position due to different distance to center. | |
436 | // we assume, that particle moves from center | |
437 | Float_t dCPV = geom->GetIPtoOuterCoverDistance(); | |
438 | Float_t dEMC = geom->GetIPtoCrystalSurface() ; | |
439 | dEMC = dEMC / dCPV ; | |
440 | vecCpv = dEMC * vecCpv - vecEmc ; | |
e3817e5f | 441 | if (axis == "X") return vecCpv.X(); |
442 | if (axis == "Y") return vecCpv.Y(); | |
443 | if (axis == "Z") return vecCpv.Z(); | |
444 | if (axis == "R") return vecCpv.Mag(); | |
445 | } | |
148b2bba | 446 | return 100000000 ; |
447 | } | |
7acf6008 | 448 | return 100000000 ; |
69183710 | 449 | } |
bc0c084c | 450 | //____________________________________________________________________________ |
8d4608b5 | 451 | Int_t AliPHOSPIDv1::GetCPVBit(AliPHOSEmcRecPoint * emc,AliPHOSCpvRecPoint * cpv, Int_t effPur, Float_t e) const |
bc0c084c | 452 | { |
e3817e5f | 453 | if(effPur>2 || effPur<0) |
454 | Error("GetCPVBit","Invalid Efficiency-Purity choice %d",effPur); | |
bc0c084c | 455 | |
e3817e5f | 456 | Float_t sigX = GetCpv2EmcDistanceCut("X",e); |
457 | Float_t sigZ = GetCpv2EmcDistanceCut("Z",e); | |
bc0c084c | 458 | |
459 | Float_t deltaX = TMath::Abs(GetDistance(emc, cpv, "X")); | |
460 | Float_t deltaZ = TMath::Abs(GetDistance(emc, cpv, "Z")); | |
8d4608b5 | 461 | |
e3817e5f | 462 | if((deltaX>sigX*(effPur+1))|(deltaZ>sigZ*(effPur+1))) |
bc0c084c | 463 | return 1;//Neutral |
464 | else | |
465 | return 0;//Charged | |
bc0c084c | 466 | } |
69183710 | 467 | |
6ad0bfa0 | 468 | //____________________________________________________________________________ |
fc7e2f43 | 469 | Int_t AliPHOSPIDv1::GetPrincipalBit(TString particle, const Double_t* p, Int_t effPur, Float_t e)const |
148b2bba | 470 | { |
50739f15 | 471 | //Is the particle inside de PCA ellipse? |
581354c5 | 472 | |
e3817e5f | 473 | particle.ToLower(); |
474 | Int_t prinbit = 0 ; | |
475 | Float_t a = GetEllipseParameter(particle,"a" , e); | |
476 | Float_t b = GetEllipseParameter(particle,"b" , e); | |
477 | Float_t c = GetEllipseParameter(particle,"c" , e); | |
478 | Float_t x0 = GetEllipseParameter(particle,"x0", e); | |
479 | Float_t y0 = GetEllipseParameter(particle,"y0", e); | |
480 | ||
481 | Float_t r = TMath::Power((p[0] - x0)/a,2) + | |
482 | TMath::Power((p[1] - y0)/b,2) + | |
483 | c*(p[0] - x0)*(p[1] - y0)/(a*b) ; | |
50739f15 | 484 | //3 different ellipses defined |
e3817e5f | 485 | if((effPur==2) && (r<1./2.)) prinbit= 1; |
486 | if((effPur==1) && (r<2. )) prinbit= 1; | |
487 | if((effPur==0) && (r<9./2.)) prinbit= 1; | |
50739f15 | 488 | |
581354c5 | 489 | if(r<0) |
e3817e5f | 490 | Error("GetPrincipalBit", "Negative square?") ; |
1f0e7ccd | 491 | |
492 | return prinbit; | |
148b2bba | 493 | |
148b2bba | 494 | } |
1f0e7ccd | 495 | //____________________________________________________________________________ |
fc7e2f43 | 496 | Int_t AliPHOSPIDv1::GetHardPhotonBit(AliPHOSEmcRecPoint * emc) const |
1f0e7ccd | 497 | { |
e3817e5f | 498 | // Set bit for identified hard photons (E > 30 GeV) |
499 | // if the second moment M2x is below the boundary | |
500 | ||
501 | Float_t e = emc->GetEnergy(); | |
502 | if (e < 30.0) return 0; | |
503 | Float_t m2x = emc->GetM2x(); | |
504 | Float_t m2xBoundary = GetParameterPhotonBoundary(0) * | |
505 | TMath::Exp(-TMath::Power(e-GetParameterPhotonBoundary(1),2)/2.0/ | |
506 | TMath::Power(GetParameterPhotonBoundary(2),2)) + | |
507 | GetParameterPhotonBoundary(3); | |
61d3d6aa | 508 | //Info("GetHardPhotonBit","E=%f, m2x=%f, boundary=%f",e,m2x,m2xBoundary); |
e3817e5f | 509 | if (m2x < m2xBoundary) |
510 | return 1;// A hard photon | |
511 | else | |
512 | return 0;// Not a hard photon | |
1f0e7ccd | 513 | } |
92f521a9 | 514 | |
e3817e5f | 515 | //____________________________________________________________________________ |
fc7e2f43 | 516 | Int_t AliPHOSPIDv1::GetHardPi0Bit(AliPHOSEmcRecPoint * emc) const |
e3817e5f | 517 | { |
518 | // Set bit for identified hard pi0 (E > 30 GeV) | |
519 | // if the second moment M2x is above the boundary | |
520 | ||
521 | Float_t e = emc->GetEnergy(); | |
522 | if (e < 30.0) return 0; | |
523 | Float_t m2x = emc->GetM2x(); | |
524 | Float_t m2xBoundary = GetParameterPi0Boundary(0) + | |
525 | e * GetParameterPi0Boundary(1); | |
61d3d6aa | 526 | //Info("GetHardPi0Bit","E=%f, m2x=%f, boundary=%f",e,m2x,m2xBoundary); |
e3817e5f | 527 | if (m2x > m2xBoundary) |
528 | return 1;// A hard pi0 | |
bc0c084c | 529 | else |
e3817e5f | 530 | return 0;// Not a hard pi0 |
f0a4c9e9 | 531 | } |
e3817e5f | 532 | |
533 | //____________________________________________________________________________ | |
8d4608b5 | 534 | TVector3 AliPHOSPIDv1::GetMomentumDirection(AliPHOSEmcRecPoint * emc, AliPHOSCpvRecPoint * )const |
88cb7938 | 535 | { |
536 | // Calculates the momentum direction: | |
537 | // 1. if only a EMC RecPoint, direction is given by IP and this RecPoint | |
538 | // 2. if a EMC RecPoint and CPV RecPoint, direction is given by the line through the 2 recpoints | |
539 | // However because of the poor position resolution of PPSD the direction is always taken as if we were | |
540 | // in case 1. | |
f0a4c9e9 | 541 | |
88cb7938 | 542 | TVector3 dir(0,0,0) ; |
9688c1dd | 543 | |
88cb7938 | 544 | TVector3 emcglobalpos ; |
545 | TMatrix dummy ; | |
bf8f1fbd | 546 | |
88cb7938 | 547 | emc->GetGlobalPosition(emcglobalpos, dummy) ; |
bf8f1fbd | 548 | |
fbf811ec | 549 | |
88cb7938 | 550 | dir = emcglobalpos ; |
88cb7938 | 551 | dir.SetMag(1.) ; |
e3817e5f | 552 | |
88cb7938 | 553 | //account correction to the position of IP |
554 | Float_t xo,yo,zo ; //Coordinates of the origin | |
555 | gAlice->Generator()->GetOrigin(xo,yo,zo) ; | |
556 | TVector3 origin(xo,yo,zo); | |
557 | dir = dir - origin ; | |
e3817e5f | 558 | |
88cb7938 | 559 | return dir ; |
7b7c1533 | 560 | } |
561 | ||
2cc71c1e | 562 | //____________________________________________________________________________ |
563 | void AliPHOSPIDv1::MakePID() | |
564 | { | |
565 | // construct the PID weight from a Bayesian Method | |
566 | ||
567 | Int_t index ; | |
cc1fe362 | 568 | const Int_t kSPECIES = AliESDtrack::kSPECIESN ; |
569 | Double_t pid[kSPECIES] = {0., 0., 0., 0., 0., 0.} ; | |
2cc71c1e | 570 | Int_t nparticles = AliPHOSGetter::Instance()->RecParticles()->GetEntriesFast() ; |
53ab54c8 | 571 | // const Int_t kMAXPARTICLES = 2000 ; |
572 | // if (nparticles >= kMAXPARTICLES) | |
573 | // Error("MakePID", "Change size of MAXPARTICLES") ; | |
574 | // Double_t stof[kSPECIES][kMAXPARTICLES] ; | |
575 | Double_t * stof[kSPECIES]; | |
576 | for (Int_t i =0; i< kSPECIES; i++) | |
577 | stof[i] = new Double_t[nparticles]; | |
578 | ||
cc1fe362 | 579 | // make the normalized distribution of pid for this event |
580 | // w(pid) in the Bayesian formulation | |
2cc71c1e | 581 | for(index = 0 ; index < nparticles ; index ++) { |
582 | AliPHOSRecParticle * recpar = AliPHOSGetter::Instance()->RecParticle(index) ; | |
cc1fe362 | 583 | |
584 | pid[AliESDtrack::kKaon0] = 0 ; | |
585 | ||
586 | if (recpar->IsPhoton() || recpar->IsHardPhoton()) | |
587 | pid[AliESDtrack::kPhoton]++ ; | |
588 | ||
589 | else if (recpar->IsPi0() || recpar->IsHardPi0()) | |
590 | pid[AliESDtrack::kPi0]++ ; | |
591 | ||
592 | else if (recpar->IsElectron()) { | |
593 | pid[AliESDtrack::kElectron]++ ; | |
594 | pid[AliESDtrack::kMuon]++ ; | |
595 | } | |
596 | ||
597 | else if (recpar->IsChargedHadron()){ | |
598 | pid[AliESDtrack::kPion]++ ; | |
599 | pid[AliESDtrack::kKaon]++ ; | |
600 | pid[AliESDtrack::kProton]++ ; | |
601 | } | |
602 | ||
603 | else if (recpar->IsNeutralHadron()) | |
604 | pid[AliESDtrack::kNeutron]++ ; | |
605 | ||
606 | else if (recpar->IsEleCon()) | |
607 | pid[AliESDtrack::kEleCon]++ ; | |
608 | ||
609 | // now get the signals probability | |
610 | // s(pid) in the Bayesian formulation | |
611 | // Tof | |
612 | stof[AliESDtrack::kPhoton][index] = fTFphoton->Eval(recpar->ToF()) ; | |
613 | stof[AliESDtrack::kPi0][index] = fTFphoton->Eval(recpar->ToF()) ; // pi0 are detected via decay photon | |
614 | stof[AliESDtrack::kElectron][index] = fTFelectron->Eval(recpar->ToF()) ; | |
615 | stof[AliESDtrack::kPion][index] = fTFchargedhadron->Eval(recpar->ToF()) ; | |
616 | stof[AliESDtrack::kKaon][index] = fTFchargedhadron->Eval(recpar->ToF()) ; | |
617 | stof[AliESDtrack::kProton][index] = fTFchargedhadron->Eval(recpar->ToF()) ; | |
618 | stof[AliESDtrack::kNeutron][index] = fTFneutralhadron->Eval(recpar->ToF()) ; | |
619 | stof[AliESDtrack::kEleCon][index] = fTFphoton->Eval(recpar->ToF()) ; // a conversion electron has the photon ToF | |
620 | stof[AliESDtrack::kKaon0][index] = 0 ; // do not know yet what to to with K0 | |
621 | ||
622 | } | |
623 | for (index = 0 ; index < kSPECIES ; index++) | |
624 | pid[index] /= nparticles ; | |
625 | ||
626 | ||
627 | for(index = 0 ; index < nparticles ; index ++) { | |
628 | // calculates the Bayesian weight | |
629 | Int_t jndex ; | |
630 | Double_t wn = 0.0 ; | |
631 | for (jndex = 0 ; jndex < kSPECIES ; jndex++) | |
632 | wn += stof[jndex][index] * pid[jndex] ; | |
633 | AliPHOSRecParticle * recpar = AliPHOSGetter::Instance()->RecParticle(index) ; | |
e74ea0e9 | 634 | if (TMath::Abs(wn)>0) |
635 | for (jndex = 0 ; jndex < kSPECIES ; jndex++) { | |
636 | recpar->SetPID(jndex, stof[jndex][index] * pid[jndex] / wn) ; | |
637 | } | |
2cc71c1e | 638 | } |
53ab54c8 | 639 | for (Int_t i =0; i< kSPECIES; i++) |
640 | delete [] stof[i]; | |
2cc71c1e | 641 | } |
642 | ||
7acf6008 | 643 | //____________________________________________________________________________ |
e3817e5f | 644 | void AliPHOSPIDv1::MakeRecParticles() |
645 | { | |
b2a60966 | 646 | // Makes a RecParticle out of a TrackSegment |
148b2bba | 647 | |
88cb7938 | 648 | AliPHOSGetter * gime = AliPHOSGetter::Instance() ; |
fbf811ec | 649 | TObjArray * emcRecPoints = gime->EmcRecPoints() ; |
650 | TObjArray * cpvRecPoints = gime->CpvRecPoints() ; | |
651 | TClonesArray * trackSegments = gime->TrackSegments() ; | |
148b2bba | 652 | if ( !emcRecPoints || !cpvRecPoints || !trackSegments ) { |
f1611b7c | 653 | Fatal("MakeRecParticles", "RecPoints or TrackSegments not found !") ; |
148b2bba | 654 | } |
fbf811ec | 655 | TClonesArray * recParticles = gime->RecParticles() ; |
01a599c9 | 656 | recParticles->Clear(); |
148b2bba | 657 | |
7b7c1533 | 658 | TIter next(trackSegments) ; |
7acf6008 | 659 | AliPHOSTrackSegment * ts ; |
6ad0bfa0 | 660 | Int_t index = 0 ; |
09fc14a0 | 661 | AliPHOSRecParticle * rp ; |
7acf6008 | 662 | while ( (ts = (AliPHOSTrackSegment *)next()) ) { |
fad3e5b9 | 663 | |
7b7c1533 | 664 | new( (*recParticles)[index] ) AliPHOSRecParticle() ; |
665 | rp = (AliPHOSRecParticle *)recParticles->At(index) ; | |
f0a4c9e9 | 666 | rp->SetTrackSegment(index) ; |
9688c1dd | 667 | rp->SetIndexInList(index) ; |
148b2bba | 668 | |
7acf6008 | 669 | AliPHOSEmcRecPoint * emc = 0 ; |
670 | if(ts->GetEmcIndex()>=0) | |
7b7c1533 | 671 | emc = (AliPHOSEmcRecPoint *) emcRecPoints->At(ts->GetEmcIndex()) ; |
fad3e5b9 | 672 | |
8d4608b5 | 673 | AliPHOSCpvRecPoint * cpv = 0 ; |
7acf6008 | 674 | if(ts->GetCpvIndex()>=0) |
8d4608b5 | 675 | cpv = (AliPHOSCpvRecPoint *) cpvRecPoints->At(ts->GetCpvIndex()) ; |
fad3e5b9 | 676 | |
bd76890a | 677 | Int_t track = 0 ; |
678 | track = ts->GetTrackIndex() ; | |
679 | ||
148b2bba | 680 | // Now set type (reconstructed) of the particle |
681 | ||
682 | // Choose the cluster energy range | |
9fa5f1d0 | 683 | |
fbf811ec | 684 | if (!emc) { |
f1611b7c | 685 | Fatal("MakeRecParticles", "-> emc(%d) = %d", ts->GetEmcIndex(), emc ) ; |
fbf811ec | 686 | } |
50739f15 | 687 | |
e3817e5f | 688 | Float_t e = emc->GetEnergy() ; |
bc0c084c | 689 | |
6f969528 | 690 | Float_t lambda[2] ; |
691 | emc->GetElipsAxis(lambda) ; | |
50739f15 | 692 | |
693 | if((lambda[0]>0.01) && (lambda[1]>0.01)){ | |
694 | // Looking PCA. Define and calculate the data (X), | |
bc0c084c | 695 | // introduce in the function X2P that gives the components (P). |
696 | ||
e3817e5f | 697 | Float_t Spher = 0. ; |
698 | Float_t Emaxdtotal = 0. ; | |
50739f15 | 699 | |
bc0c084c | 700 | if((lambda[0]+lambda[1])!=0) |
e3817e5f | 701 | Spher=fabs(lambda[0]-lambda[1])/(lambda[0]+lambda[1]); |
50739f15 | 702 | |
e3817e5f | 703 | Emaxdtotal=emc->GetMaximalEnergy()/emc->GetEnergy(); |
50739f15 | 704 | |
705 | fX[0] = lambda[0] ; | |
706 | fX[1] = lambda[1] ; | |
707 | fX[2] = emc->GetDispersion() ; | |
e3817e5f | 708 | fX[3] = Spher ; |
50739f15 | 709 | fX[4] = emc->GetMultiplicity() ; |
e3817e5f | 710 | fX[5] = Emaxdtotal ; |
50739f15 | 711 | fX[6] = emc->GetCoreEnergy() ; |
712 | ||
e3817e5f | 713 | fPrincipalPhoton->X2P(fX,fPPhoton); |
714 | fPrincipalPi0 ->X2P(fX,fPPi0); | |
1f0e7ccd | 715 | |
50739f15 | 716 | } |
717 | else{ | |
e3817e5f | 718 | fPPhoton[0]=-100.0; //We do not accept clusters with |
719 | fPPhoton[1]=-100.0; //one cell as a photon-like | |
720 | fPPi0[0] =-100.0; | |
721 | fPPi0[1] =-100.0; | |
50739f15 | 722 | } |
723 | ||
2cc71c1e | 724 | Float_t time = emc->GetTime() ; |
725 | rp->SetTof(time) ; | |
9fa5f1d0 | 726 | |
bc0c084c | 727 | // Loop of Efficiency-Purity (the 3 points of purity or efficiency |
728 | // are taken into account to set the particle identification) | |
e3817e5f | 729 | for(Int_t effPur = 0; effPur < 3 ; effPur++){ |
50739f15 | 730 | |
bc0c084c | 731 | // Looking at the CPV detector. If RCPV greater than CpvEmcDistance, |
732 | // 1st,2nd or 3rd bit (depending on the efficiency-purity point ) | |
733 | // is set to 1 | |
e3817e5f | 734 | if(GetCPVBit(emc, cpv, effPur,e) == 1 ) |
735 | rp->SetPIDBit(effPur) ; | |
f0a4c9e9 | 736 | |
50739f15 | 737 | // Looking the TOF. If TOF smaller than gate, 4th, 5th or 6th |
738 | // bit (depending on the efficiency-purity point )is set to 1 | |
2cc71c1e | 739 | if(time< (*fParameters)(3,effPur)) |
e3817e5f | 740 | rp->SetPIDBit(effPur+3) ; |
2cc71c1e | 741 | |
e3817e5f | 742 | //Photon PCA |
50739f15 | 743 | //If we are inside the ellipse, 7th, 8th or 9th |
744 | // bit (depending on the efficiency-purity point )is set to 1 | |
e3817e5f | 745 | if(GetPrincipalBit("photon",fPPhoton,effPur,e) == 1) |
746 | rp->SetPIDBit(effPur+6) ; | |
1f0e7ccd | 747 | |
e3817e5f | 748 | //Pi0 PCA |
1f0e7ccd | 749 | //If we are inside the ellipse, 10th, 11th or 12th |
750 | // bit (depending on the efficiency-purity point )is set to 1 | |
e3817e5f | 751 | if(GetPrincipalBit("pi0" ,fPPi0 ,effPur,e) == 1) |
752 | rp->SetPIDBit(effPur+9) ; | |
f0a4c9e9 | 753 | } |
e3817e5f | 754 | if(GetHardPhotonBit(emc)) |
755 | rp->SetPIDBit(12) ; | |
756 | if(GetHardPi0Bit (emc)) | |
757 | rp->SetPIDBit(13) ; | |
1f0e7ccd | 758 | |
bd76890a | 759 | if(track >= 0) |
760 | rp->SetPIDBit(14) ; | |
761 | ||
9fa5f1d0 | 762 | //Set momentum, energy and other parameters |
50739f15 | 763 | Float_t encal = GetCalibratedEnergy(e); |
9fa5f1d0 | 764 | TVector3 dir = GetMomentumDirection(emc,cpv) ; |
765 | dir.SetMag(encal) ; | |
766 | rp->SetMomentum(dir.X(),dir.Y(),dir.Z(),encal) ; | |
767 | rp->SetCalcMass(0); | |
e0ed2e49 | 768 | rp->Name(); //If photon sets the particle pdg name to gamma |
e747b8da | 769 | rp->SetProductionVertex(0,0,0,0); |
770 | rp->SetFirstMother(-1); | |
771 | rp->SetLastMother(-1); | |
772 | rp->SetFirstDaughter(-1); | |
773 | rp->SetLastDaughter(-1); | |
774 | rp->SetPolarisation(0,0,0); | |
d956e9b7 | 775 | //Set the position in global coordinate system from the RecPoint |
776 | AliPHOSGeometry * geom = gime->PHOSGeometry() ; | |
777 | AliPHOSTrackSegment * ts = gime->TrackSegment(rp->GetPHOSTSIndex()) ; | |
778 | AliPHOSEmcRecPoint * erp = gime->EmcRecPoint(ts->GetEmcIndex()) ; | |
779 | TVector3 pos ; | |
780 | geom->GetGlobal(erp, pos) ; | |
781 | rp->SetPos(pos); | |
6ad0bfa0 | 782 | index++ ; |
783 | } | |
6ad0bfa0 | 784 | } |
e3817e5f | 785 | |
09fc14a0 | 786 | //____________________________________________________________________________ |
88cb7938 | 787 | void AliPHOSPIDv1::Print() const |
09fc14a0 | 788 | { |
b2a60966 | 789 | // Print the parameters used for the particle type identification |
bc0c084c | 790 | |
88cb7938 | 791 | Info("Print", "=============== AliPHOSPIDv1 ================") ; |
792 | printf("Making PID\n") ; | |
793 | printf(" Pricipal analysis file from 0.5 to 100 %s\n", fFileNamePrincipalPhoton.Data() ) ; | |
794 | printf(" Name of parameters file %s\n", fFileNameParameters.Data() ) ; | |
795 | printf(" Matrix of Parameters: 14x4\n") ; | |
796 | printf(" Energy Calibration 1x3 [3 parametres to calibrate energy: A + B* E + C * E^2]\n") ; | |
797 | printf(" RCPV 2x3 rows x and z, columns function cut parameters\n") ; | |
798 | printf(" TOF 1x3 [High Eff-Low Pur,Medium Eff-Pur, Low Eff-High Pur]\n") ; | |
799 | printf(" PCA 5x4 [5 ellipse parametres and 4 parametres to calculate them: A/Sqrt(E) + B* E + C * E^2 + D]\n") ; | |
800 | Printf(" Pi0 PCA 5x3 [5 ellipse parametres and 3 parametres to calculate them: A + B* E + C * E^2]\n") ; | |
50739f15 | 801 | fParameters->Print() ; |
09fc14a0 | 802 | } |
803 | ||
a496c46c | 804 | |
69183710 | 805 | |
7acf6008 | 806 | //____________________________________________________________________________ |
a4e98857 | 807 | void AliPHOSPIDv1::PrintRecParticles(Option_t * option) |
808 | { | |
dd5c4038 | 809 | // Print table of reconstructed particles |
810 | ||
88cb7938 | 811 | AliPHOSGetter *gime = AliPHOSGetter::Instance() ; |
bf8f1fbd | 812 | |
88cb7938 | 813 | TClonesArray * recParticles = gime->RecParticles() ; |
21cd0c07 | 814 | |
815 | TString message ; | |
3bf72d32 | 816 | message = "\nevent " ; |
817 | message += gAlice->GetEvNumber() ; | |
818 | message += " found " ; | |
819 | message += recParticles->GetEntriesFast(); | |
820 | message += " RecParticles\n" ; | |
821 | ||
7acf6008 | 822 | if(strstr(option,"all")) { // printing found TS |
3bf72d32 | 823 | message += "\n PARTICLE Index \n" ; |
7acf6008 | 824 | |
825 | Int_t index ; | |
7b7c1533 | 826 | for (index = 0 ; index < recParticles->GetEntries() ; index++) { |
21cd0c07 | 827 | AliPHOSRecParticle * rp = (AliPHOSRecParticle * ) recParticles->At(index) ; |
3bf72d32 | 828 | message += "\n" ; |
829 | message += rp->Name().Data() ; | |
830 | message += " " ; | |
831 | message += rp->GetIndexInList() ; | |
832 | message += " " ; | |
833 | message += rp->GetType() ; | |
7acf6008 | 834 | } |
3bf72d32 | 835 | } |
836 | Info("Print", message.Data() ) ; | |
69183710 | 837 | } |
88cb7938 | 838 | |
839 | //____________________________________________________________________________ | |
840 | void AliPHOSPIDv1::SetParameters() | |
841 | { | |
842 | // PCA : To do the Principal Components Analysis it is necessary | |
843 | // the Principal file, which is opened here | |
844 | fX = new double[7]; // Data for the PCA | |
845 | fPPhoton = new double[7]; // Eigenvalues of the PCA | |
846 | fPPi0 = new double[7]; // Eigenvalues of the Pi0 PCA | |
847 | ||
848 | // Read photon principals from the photon file | |
849 | ||
850 | fFileNamePrincipalPhoton = "$ALICE_ROOT/PHOS/PCA8pa15_0.5-100.root" ; | |
851 | TFile f( fFileNamePrincipalPhoton.Data(), "read" ) ; | |
852 | fPrincipalPhoton = dynamic_cast<TPrincipal*> (f.Get("principal")) ; | |
853 | f.Close() ; | |
854 | ||
855 | // Read pi0 principals from the pi0 file | |
856 | ||
857 | fFileNamePrincipalPi0 = "$ALICE_ROOT/PHOS/PCA_pi0_40-120.root" ; | |
858 | TFile fPi0( fFileNamePrincipalPi0.Data(), "read" ) ; | |
859 | fPrincipalPi0 = dynamic_cast<TPrincipal*> (fPi0.Get("principal")) ; | |
860 | fPi0.Close() ; | |
861 | ||
862 | // Open parameters file and initialization of the Parameters matrix. | |
863 | // In the File Parameters.dat are all the parameters. These are introduced | |
864 | // in a matrix of 16x4 | |
865 | // | |
866 | // All the parameters defined in this file are, in order of row: | |
867 | // line 0 : calibration | |
868 | // lines 1,2 : CPV rectangular cat for X and Z | |
869 | // line 3 : TOF cut | |
870 | // lines 4-8 : parameters to calculate photon PCA ellipse | |
871 | // lines 9-13: parameters to calculate pi0 PCA ellipse | |
872 | // lines 14-15: parameters to calculate border for high-pt photons and pi0 | |
873 | ||
874 | fFileNameParameters = gSystem->ExpandPathName("$ALICE_ROOT/PHOS/Parameters.dat"); | |
875 | fParameters = new TMatrix(16,4) ; | |
876 | const Int_t maxLeng=255; | |
877 | char string[maxLeng]; | |
878 | ||
879 | // Open a text file with PID parameters | |
880 | FILE *fd = fopen(fFileNameParameters.Data(),"r"); | |
881 | if (!fd) | |
882 | Fatal("SetParameter","File %s with a PID parameters cannot be opened\n", | |
883 | fFileNameParameters.Data()); | |
884 | ||
885 | Int_t i=0; | |
886 | // Read parameter file line-by-line and skip empty line and comments | |
887 | while (fgets(string,maxLeng,fd) != NULL) { | |
888 | if (string[0] == '\n' ) continue; | |
889 | if (string[0] == '!' ) continue; | |
890 | sscanf(string, "%f %f %f %f", | |
891 | &(*fParameters)(i,0), &(*fParameters)(i,1), | |
892 | &(*fParameters)(i,2), &(*fParameters)(i,3)); | |
893 | i++; | |
2cc71c1e | 894 | //Info("SetParameters", "line %d: %s",i,string); |
88cb7938 | 895 | } |
896 | fclose(fd); | |
897 | } | |
898 | ||
899 | //____________________________________________________________________________ | |
900 | void AliPHOSPIDv1::SetParameterCalibration(Int_t i,Float_t param) | |
901 | { | |
902 | // Set parameter "Calibration" i to a value param | |
903 | if(i>2 || i<0) | |
904 | Error("SetParameterCalibration","Invalid parameter number: %d",i); | |
905 | else | |
906 | (*fParameters)(0,i) = param ; | |
907 | } | |
908 | ||
909 | //____________________________________________________________________________ | |
910 | void AliPHOSPIDv1::SetParameterCpv2Emc(Int_t i, TString axis, Float_t cut) | |
911 | { | |
912 | // Set the parameters to calculate Cpv-to-Emc Distance Cut depending on | |
913 | // Purity-Efficiency point i | |
914 | ||
915 | if(i>2 || i<0) | |
916 | Error("SetParameterCpv2Emc","Invalid parameter number: %d",i); | |
917 | else { | |
918 | axis.ToLower(); | |
919 | if (axis == "x") (*fParameters)(1,i) = cut; | |
920 | else if (axis == "z") (*fParameters)(2,i) = cut; | |
921 | else Error("SetParameterCpv2Emc","Invalid axis name: %s",axis.Data()); | |
922 | } | |
923 | } | |
924 | ||
925 | //____________________________________________________________________________ | |
926 | void AliPHOSPIDv1::SetParameterPhotonBoundary(Int_t i,Float_t param) | |
927 | { | |
928 | // Set parameter "Hard photon boundary" i to a value param | |
929 | if(i>4 || i<0) | |
930 | Error("SetParameterPhotonBoundary","Invalid parameter number: %d",i); | |
931 | else | |
932 | (*fParameters)(14,i) = param ; | |
933 | } | |
934 | ||
935 | //____________________________________________________________________________ | |
936 | void AliPHOSPIDv1::SetParameterPi0Boundary(Int_t i,Float_t param) | |
937 | { | |
938 | // Set parameter "Hard pi0 boundary" i to a value param | |
939 | if(i>1 || i<0) | |
940 | Error("SetParameterPi0Boundary","Invalid parameter number: %d",i); | |
941 | else | |
942 | (*fParameters)(15,i) = param ; | |
943 | } | |
944 | ||
945 | //_____________________________________________________________________________ | |
946 | void AliPHOSPIDv1::SetParameterTimeGate(Int_t i, Float_t gate) | |
947 | { | |
948 | // Set the parameter TimeGate depending on Purity-Efficiency point i | |
949 | if (i>2 || i<0) | |
950 | Error("SetParameterTimeGate","Invalid Efficiency-Purity choice %d",i); | |
951 | else | |
952 | (*fParameters)(3,i)= gate ; | |
953 | } | |
954 | ||
955 | //_____________________________________________________________________________ | |
956 | void AliPHOSPIDv1::SetParameterToCalculateEllipse(TString particle, TString param, Int_t i, Float_t par) | |
957 | { | |
958 | // Set the parameter "i" that is needed to calculate the ellipse | |
959 | // parameter "param" for a particle "particle" | |
960 | ||
961 | particle.ToLower(); | |
962 | param. ToLower(); | |
963 | Int_t p= -1; | |
964 | Int_t offset=0; | |
965 | ||
966 | if (particle == "photon") offset=0; | |
967 | else if (particle == "pi0") offset=5; | |
968 | else | |
969 | Error("SetParameterToCalculateEllipse","Wrong particle name: %s (choose from pi0/photon)\n",particle.Data()); | |
970 | ||
971 | if (param.Contains("a")) p=4+offset; | |
972 | else if(param.Contains("b")) p=5+offset; | |
973 | else if(param.Contains("c")) p=6+offset; | |
974 | else if(param.Contains("x0"))p=7+offset; | |
975 | else if(param.Contains("y0"))p=8+offset; | |
976 | if((i>4)||(i<0)) | |
977 | Error("SetEllipseParameter", "No parameter with index %d", i) ; | |
978 | else if(p==-1) | |
979 | Error("SetEllipseParameter", "No parameter with name %s", param.Data() ) ; | |
980 | else | |
981 | (*fParameters)(p,i) = par ; | |
982 | } | |
983 | ||
984 | //____________________________________________________________________________ | |
985 | void AliPHOSPIDv1::Unload() | |
986 | { | |
987 | AliPHOSGetter * gime = AliPHOSGetter::Instance() ; | |
988 | gime->PhosLoader()->UnloadRecPoints() ; | |
989 | gime->PhosLoader()->UnloadTracks() ; | |
990 | gime->PhosLoader()->UnloadRecParticles() ; | |
991 | } | |
992 | ||
993 | //____________________________________________________________________________ | |
90cceaf6 | 994 | void AliPHOSPIDv1::WriteRecParticles() |
88cb7938 | 995 | { |
996 | ||
997 | AliPHOSGetter *gime = AliPHOSGetter::Instance() ; | |
998 | ||
999 | TClonesArray * recParticles = gime->RecParticles() ; | |
1000 | recParticles->Expand(recParticles->GetEntriesFast() ) ; | |
1001 | TTree * treeP = gime->TreeP(); | |
1002 | ||
1003 | //First rp | |
1004 | Int_t bufferSize = 32000 ; | |
1005 | TBranch * rpBranch = treeP->Branch("PHOSRP",&recParticles,bufferSize); | |
1006 | rpBranch->SetTitle(BranchName()); | |
1007 | ||
1008 | rpBranch->Fill() ; | |
1009 | ||
1010 | gime->WriteRecParticles("OVERWRITE"); | |
1011 | gime->WritePID("OVERWRITE"); | |
1012 | } | |
1013 |