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