Updated PID with photon/pi0 recognition by M2x at high-pt
[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
6ad0bfa0 114//____________________________________________________________________________
1cb7c1ee 115AliPHOSPIDv1::AliPHOSPIDv1():AliPHOSPID()
116{
a4e98857 117 // default ctor
148b2bba 118
8d0f3f77 119 InitParameters() ;
92f521a9 120 fDefaultInit = kTRUE ;
121
7acf6008 122}
123
124//____________________________________________________________________________
581354c5 125AliPHOSPIDv1::AliPHOSPIDv1(AliPHOSPIDv1 & pid ):AliPHOSPID(pid)
126{
386aef34 127 // ctor
581354c5 128 InitParameters() ;
129
130 Init() ;
131 fDefaultInit = kFALSE ;
132
133}
134
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
148b2bba 162//____________________________________________________________________________
a496c46c 163const TString AliPHOSPIDv1::BranchName() const
164{
165 TString branchName(GetName() ) ;
166 branchName.Remove(branchName.Index(Version())-1) ;
167 return branchName ;
168}
169
170//____________________________________________________________________________
148b2bba 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
1cb7c1ee 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//____________________________________________________________________________
243const Float_t AliPHOSPIDv1::GetParameterCalibration(const Int_t i) const
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//____________________________________________________________________________
271const Float_t AliPHOSPIDv1::GetParameterTimeGate(const Int_t i) const
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}
452//________________________________________________________________________
453const Float_t AliPHOSPIDv1::GetDistance(AliPHOSEmcRecPoint * emc,AliPHOSRecPoint * cpv, Option_t * axis)const
69183710 454{
455 // Calculates the distance between the EMC RecPoint and the PPSD RecPoint
148b2bba 456
7b7c1533 457 const AliPHOSGeometry * geom = AliPHOSGetter::GetInstance()->PHOSGeometry() ;
69183710 458 TVector3 vecEmc ;
7acf6008 459 TVector3 vecCpv ;
148b2bba 460 if(cpv){
461 emc->GetLocalPosition(vecEmc) ;
462 cpv->GetLocalPosition(vecCpv) ;
463 if(emc->GetPHOSMod() == cpv->GetPHOSMod()){
464 // Correct to difference in CPV and EMC position due to different distance to center.
465 // we assume, that particle moves from center
466 Float_t dCPV = geom->GetIPtoOuterCoverDistance();
467 Float_t dEMC = geom->GetIPtoCrystalSurface() ;
468 dEMC = dEMC / dCPV ;
469 vecCpv = dEMC * vecCpv - vecEmc ;
e3817e5f 470 if (axis == "X") return vecCpv.X();
471 if (axis == "Y") return vecCpv.Y();
472 if (axis == "Z") return vecCpv.Z();
473 if (axis == "R") return vecCpv.Mag();
474 }
148b2bba 475 return 100000000 ;
476 }
7acf6008 477 return 100000000 ;
69183710 478}
bc0c084c 479//____________________________________________________________________________
e3817e5f 480const Int_t AliPHOSPIDv1::GetCPVBit(AliPHOSEmcRecPoint * emc,AliPHOSRecPoint * cpv,const Int_t effPur, const Float_t e) const
bc0c084c 481{
e3817e5f 482 if(effPur>2 || effPur<0)
483 Error("GetCPVBit","Invalid Efficiency-Purity choice %d",effPur);
bc0c084c 484
e3817e5f 485 Float_t sigX = GetCpv2EmcDistanceCut("X",e);
486 Float_t sigZ = GetCpv2EmcDistanceCut("Z",e);
bc0c084c 487
488 Float_t deltaX = TMath::Abs(GetDistance(emc, cpv, "X"));
489 Float_t deltaZ = TMath::Abs(GetDistance(emc, cpv, "Z"));
490
e3817e5f 491 if((deltaX>sigX*(effPur+1))|(deltaZ>sigZ*(effPur+1)))
bc0c084c 492 return 1;//Neutral
493 else
494 return 0;//Charged
bc0c084c 495}
69183710 496
69183710 497//____________________________________________________________________________
e3817e5f 498const Int_t AliPHOSPIDv1::GetPrincipalBit(TString particle, const Double_t* p,const Int_t effPur, const Float_t e)const
148b2bba 499{
50739f15 500 //Is the particle inside de PCA ellipse?
581354c5 501
e3817e5f 502 particle.ToLower();
503 Int_t prinbit = 0 ;
504 Float_t a = GetEllipseParameter(particle,"a" , e);
505 Float_t b = GetEllipseParameter(particle,"b" , e);
506 Float_t c = GetEllipseParameter(particle,"c" , e);
507 Float_t x0 = GetEllipseParameter(particle,"x0", e);
508 Float_t y0 = GetEllipseParameter(particle,"y0", e);
509
510 Float_t r = TMath::Power((p[0] - x0)/a,2) +
511 TMath::Power((p[1] - y0)/b,2) +
512 c*(p[0] - x0)*(p[1] - y0)/(a*b) ;
50739f15 513 //3 different ellipses defined
e3817e5f 514 if((effPur==2) && (r<1./2.)) prinbit= 1;
515 if((effPur==1) && (r<2. )) prinbit= 1;
516 if((effPur==0) && (r<9./2.)) prinbit= 1;
50739f15 517
581354c5 518 if(r<0)
e3817e5f 519 Error("GetPrincipalBit", "Negative square?") ;
1f0e7ccd 520
521 return prinbit;
148b2bba 522
148b2bba 523}
1f0e7ccd 524//____________________________________________________________________________
e3817e5f 525const Int_t AliPHOSPIDv1::GetHardPhotonBit(AliPHOSEmcRecPoint * emc) const
1f0e7ccd 526{
e3817e5f 527 // Set bit for identified hard photons (E > 30 GeV)
528 // if the second moment M2x is below the boundary
529
530 Float_t e = emc->GetEnergy();
531 if (e < 30.0) return 0;
532 Float_t m2x = emc->GetM2x();
533 Float_t m2xBoundary = GetParameterPhotonBoundary(0) *
534 TMath::Exp(-TMath::Power(e-GetParameterPhotonBoundary(1),2)/2.0/
535 TMath::Power(GetParameterPhotonBoundary(2),2)) +
536 GetParameterPhotonBoundary(3);
537 Info("GetHardPhotonBit","E=%f, m2x=%f, boundary=%f",e,m2x,m2xBoundary);
538 if (m2x < m2xBoundary)
539 return 1;// A hard photon
540 else
541 return 0;// Not a hard photon
1f0e7ccd 542}
92f521a9 543
e3817e5f 544//____________________________________________________________________________
545const Int_t AliPHOSPIDv1::GetHardPi0Bit(AliPHOSEmcRecPoint * emc) const
546{
547 // Set bit for identified hard pi0 (E > 30 GeV)
548 // if the second moment M2x is above the boundary
549
550 Float_t e = emc->GetEnergy();
551 if (e < 30.0) return 0;
552 Float_t m2x = emc->GetM2x();
553 Float_t m2xBoundary = GetParameterPi0Boundary(0) +
554 e * GetParameterPi0Boundary(1);
555 Info("GetHardPi0Bit","E=%f, m2x=%f, boundary=%f",e,m2x,m2xBoundary);
556 if (m2x > m2xBoundary)
557 return 1;// A hard pi0
bc0c084c 558 else
e3817e5f 559 return 0;// Not a hard pi0
f0a4c9e9 560}
e3817e5f 561
562//____________________________________________________________________________
50739f15 563void AliPHOSPIDv1::SetParameters()
a496c46c 564{
565 // PCA : To do the Principal Components Analysis it is necessary
566 // the Principal file, which is opened here
e3817e5f 567 fX = new double[7]; // Data for the PCA
568 fPPhoton = new double[7]; // Eigenvalues of the PCA
569 fPPi0 = new double[7]; // Eigenvalues of the Pi0 PCA
1f0e7ccd 570
571 // Read photon principals from the photon file
e0ed2e49 572
e3817e5f 573 fFileNamePrincipalPhoton = "$ALICE_ROOT/PHOS/PCA8pa15_0.5-100.root" ;
574 TFile f( fFileNamePrincipalPhoton.Data(), "read" ) ;
575 fPrincipalPhoton = dynamic_cast<TPrincipal*> (f.Get("principal")) ;
a496c46c 576 f.Close() ;
1f0e7ccd 577
578 // Read pi0 principals from the pi0 file
579
e3817e5f 580 fFileNamePrincipalPi0 = "$ALICE_ROOT/PHOS/PCA_pi0_40-120.root" ;
581 TFile fPi0( fFileNamePrincipalPi0.Data(), "read" ) ;
582 fPrincipalPi0 = dynamic_cast<TPrincipal*> (fPi0.Get("principal")) ;
1f0e7ccd 583 fPi0.Close() ;
584
bc0c084c 585 // Open parameters file and initialization of the Parameters matrix.
586 // In the File Parameters.dat are all the parameters. These are introduced
e3817e5f 587 // in a matrix of 16x4
50739f15 588 //
589 // All the parameters defined in this file are, in order of row:
e3817e5f 590 // line 0 : calibration
591 // lines 1,2 : CPV rectangular cat for X and Z
592 // line 3 : TOF cut
593 // lines 4-8 : parameters to calculate photon PCA ellipse
594 // lines 9-13: parameters to calculate pi0 PCA ellipse
595 // lines 14-15: parameters to calculate border for high-pt photons and pi0
596
597 fFileNameParameters = gSystem->ExpandPathName("$ALICE_ROOT/PHOS/Parameters.dat");
598 fParameters = new TMatrix(16,4) ;
599 const Int_t maxLeng=255;
600 char string[maxLeng];
1f0e7ccd 601
602 // Open a text file with PID parameters
e3817e5f 603 FILE *fd = fopen(fFileNameParameters.Data(),"r");
1f0e7ccd 604 if (!fd)
605 Fatal("SetParameter","File %s with a PID parameters cannot be opened\n",
e3817e5f 606 fFileNameParameters.Data());
1f0e7ccd 607
608 Int_t i=0;
609 // Read parameter file line-by-line and skip empty line and comments
e3817e5f 610 while (fgets(string,maxLeng,fd) != NULL) {
1f0e7ccd 611 if (string[0] == '\n' ) continue;
612 if (string[0] == '!' ) continue;
e3817e5f 613 sscanf(string, "%f %f %f %f",
1f0e7ccd 614 &(*fParameters)(i,0), &(*fParameters)(i,1),
bc0c084c 615 &(*fParameters)(i,2), &(*fParameters)(i,3));
1f0e7ccd 616 i++;
e3817e5f 617 printf("line %d: %s",i,string);
a496c46c 618 }
1f0e7ccd 619 fclose(fd);
f0a4c9e9 620}
f0a4c9e9 621
50739f15 622//________________________________________________________________________
7acf6008 623void AliPHOSPIDv1::Exec(Option_t * option)
6ad0bfa0 624{
f035f6ce 625 //Steering method
9688c1dd 626
bf8f1fbd 627 if( strcmp(GetName(), "")== 0 )
7acf6008 628 Init() ;
bf8f1fbd 629
630 if(strstr(option,"tim"))
7acf6008 631 gBenchmark->Start("PHOSPID");
bf8f1fbd 632
633 if(strstr(option,"print")) {
7b7c1533 634 Print("") ;
635 return ;
bf8f1fbd 636 }
9688c1dd 637
148b2bba 638
e3817e5f 639// gAlice->GetEvent(0) ;
640
641// //check, if the branch with name of this" already exits?
642// if (gAlice->TreeR()) {
643// TObjArray * lob = (TObjArray*)gAlice->TreeR()->GetListOfBranches() ;
644// TIter next(lob) ;
645// TBranch * branch = 0 ;
646// Bool_t phospidfound = kFALSE, pidfound = kFALSE ;
647
648// TString taskName(GetName()) ;
649// taskName.Remove(taskName.Index(Version())-1) ;
650
651// while ( (branch = (TBranch*)next()) && (!phospidfound || !pidfound) ) {
652// if ( (strcmp(branch->GetName(), "PHOSPID")==0) && (strcmp(branch->GetTitle(), taskName.Data())==0) )
653// phospidfound = kTRUE ;
654
655// else if ( (strcmp(branch->GetName(), "AliPHOSPID")==0) && (strcmp(branch->GetTitle(), taskName.Data())==0) )
656// pidfound = kTRUE ;
657// }
658
659// if ( phospidfound || pidfound ) {
660// Error("Exec", "RecParticles and/or PIDtMaker branch with name %s already exists", taskName.Data() ) ;
661// return ;
662// }
663// }
664
665// Int_t nevents = (Int_t) gAlice->TreeE()->GetEntries() ;
666// Int_t ievent ;
667// AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
668
fbf811ec 669 AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
670 if(gime->BranchExists("RecParticles") )
671 return ;
e3817e5f 672 Int_t nevents = gime->MaxEvent() ; //(Int_t) gAlice->TreeE()->GetEntries() ;
7b7c1533 673 Int_t ievent ;
fbf811ec 674
e3817e5f 675
7b7c1533 676 for(ievent = 0; ievent < nevents; ievent++){
a496c46c 677 gime->Event(ievent,"R") ;
148b2bba 678
e3817e5f 679 MakeRecParticles() ;
680
681 WriteRecParticles(ievent);
682
683 if(strstr(option,"deb"))
684 PrintRecParticles(option) ;
685
686 //increment the total number of rec particles per run
687 fRecParticlesInRun += gime->RecParticles(BranchName())->GetEntriesFast() ;
688
7acf6008 689 }
9688c1dd 690
7acf6008 691 if(strstr(option,"tim")){
692 gBenchmark->Stop("PHOSPID");
21cd0c07 693 Info("Exec", "took %f seconds for PID %f seconds per event",
694 gBenchmark->GetCpuTime("PHOSPID"),
695 gBenchmark->GetCpuTime("PHOSPID")/nevents) ;
696 }
7acf6008 697}
7b7c1533 698
699//____________________________________________________________________________
e3817e5f 700void AliPHOSPIDv1::MakeRecParticles()
701{
b2a60966 702 // Makes a RecParticle out of a TrackSegment
148b2bba 703
7b7c1533 704 AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
fbf811ec 705 TObjArray * emcRecPoints = gime->EmcRecPoints() ;
706 TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
707 TClonesArray * trackSegments = gime->TrackSegments() ;
148b2bba 708 if ( !emcRecPoints || !cpvRecPoints || !trackSegments ) {
f1611b7c 709 Fatal("MakeRecParticles", "RecPoints or TrackSegments not found !") ;
148b2bba 710 }
fbf811ec 711 TClonesArray * recParticles = gime->RecParticles() ;
01a599c9 712 recParticles->Clear();
148b2bba 713
7b7c1533 714 TIter next(trackSegments) ;
7acf6008 715 AliPHOSTrackSegment * ts ;
6ad0bfa0 716 Int_t index = 0 ;
09fc14a0 717 AliPHOSRecParticle * rp ;
7acf6008 718 while ( (ts = (AliPHOSTrackSegment *)next()) ) {
fad3e5b9 719
7b7c1533 720 new( (*recParticles)[index] ) AliPHOSRecParticle() ;
721 rp = (AliPHOSRecParticle *)recParticles->At(index) ;
f0a4c9e9 722 rp->SetTrackSegment(index) ;
9688c1dd 723 rp->SetIndexInList(index) ;
148b2bba 724
7acf6008 725 AliPHOSEmcRecPoint * emc = 0 ;
726 if(ts->GetEmcIndex()>=0)
7b7c1533 727 emc = (AliPHOSEmcRecPoint *) emcRecPoints->At(ts->GetEmcIndex()) ;
fad3e5b9 728
7acf6008 729 AliPHOSRecPoint * cpv = 0 ;
730 if(ts->GetCpvIndex()>=0)
7b7c1533 731 cpv = (AliPHOSRecPoint *) cpvRecPoints->At(ts->GetCpvIndex()) ;
fad3e5b9 732
148b2bba 733 // Now set type (reconstructed) of the particle
734
735 // Choose the cluster energy range
9fa5f1d0 736
fbf811ec 737 if (!emc) {
f1611b7c 738 Fatal("MakeRecParticles", "-> emc(%d) = %d", ts->GetEmcIndex(), emc ) ;
fbf811ec 739 }
50739f15 740
e3817e5f 741 Float_t e = emc->GetEnergy() ;
bc0c084c 742
6f969528 743 Float_t lambda[2] ;
744 emc->GetElipsAxis(lambda) ;
50739f15 745
746 if((lambda[0]>0.01) && (lambda[1]>0.01)){
747 // Looking PCA. Define and calculate the data (X),
bc0c084c 748 // introduce in the function X2P that gives the components (P).
749
e3817e5f 750 Float_t Spher = 0. ;
751 Float_t Emaxdtotal = 0. ;
50739f15 752
bc0c084c 753 if((lambda[0]+lambda[1])!=0)
e3817e5f 754 Spher=fabs(lambda[0]-lambda[1])/(lambda[0]+lambda[1]);
50739f15 755
e3817e5f 756 Emaxdtotal=emc->GetMaximalEnergy()/emc->GetEnergy();
50739f15 757
758 fX[0] = lambda[0] ;
759 fX[1] = lambda[1] ;
760 fX[2] = emc->GetDispersion() ;
e3817e5f 761 fX[3] = Spher ;
50739f15 762 fX[4] = emc->GetMultiplicity() ;
e3817e5f 763 fX[5] = Emaxdtotal ;
50739f15 764 fX[6] = emc->GetCoreEnergy() ;
765
e3817e5f 766 fPrincipalPhoton->X2P(fX,fPPhoton);
767 fPrincipalPi0 ->X2P(fX,fPPi0);
1f0e7ccd 768
50739f15 769 }
770 else{
e3817e5f 771 fPPhoton[0]=-100.0; //We do not accept clusters with
772 fPPhoton[1]=-100.0; //one cell as a photon-like
773 fPPi0[0] =-100.0;
774 fPPi0[1] =-100.0;
50739f15 775 }
776
6f969528 777 Float_t time =emc->GetTime() ;
9fa5f1d0 778
bc0c084c 779 // Loop of Efficiency-Purity (the 3 points of purity or efficiency
780 // are taken into account to set the particle identification)
e3817e5f 781 for(Int_t effPur = 0; effPur < 3 ; effPur++){
50739f15 782
bc0c084c 783 // Looking at the CPV detector. If RCPV greater than CpvEmcDistance,
784 // 1st,2nd or 3rd bit (depending on the efficiency-purity point )
785 // is set to 1
e3817e5f 786 if(GetCPVBit(emc, cpv, effPur,e) == 1 )
787 rp->SetPIDBit(effPur) ;
f0a4c9e9 788
50739f15 789 // Looking the TOF. If TOF smaller than gate, 4th, 5th or 6th
790 // bit (depending on the efficiency-purity point )is set to 1
e3817e5f 791 if(time< (*fParameters)(2,effPur))
792 rp->SetPIDBit(effPur+3) ;
50739f15 793
e3817e5f 794 //Photon PCA
50739f15 795 //If we are inside the ellipse, 7th, 8th or 9th
796 // bit (depending on the efficiency-purity point )is set to 1
e3817e5f 797 if(GetPrincipalBit("photon",fPPhoton,effPur,e) == 1)
798 rp->SetPIDBit(effPur+6) ;
1f0e7ccd 799
e3817e5f 800 //Pi0 PCA
1f0e7ccd 801 //If we are inside the ellipse, 10th, 11th or 12th
802 // bit (depending on the efficiency-purity point )is set to 1
e3817e5f 803 if(GetPrincipalBit("pi0" ,fPPi0 ,effPur,e) == 1)
804 rp->SetPIDBit(effPur+9) ;
f0a4c9e9 805 }
e3817e5f 806 if(GetHardPhotonBit(emc))
807 rp->SetPIDBit(12) ;
808 if(GetHardPi0Bit (emc))
809 rp->SetPIDBit(13) ;
1f0e7ccd 810
9fa5f1d0 811 //Set momentum, energy and other parameters
50739f15 812 Float_t encal = GetCalibratedEnergy(e);
9fa5f1d0 813 TVector3 dir = GetMomentumDirection(emc,cpv) ;
814 dir.SetMag(encal) ;
815 rp->SetMomentum(dir.X(),dir.Y(),dir.Z(),encal) ;
816 rp->SetCalcMass(0);
e0ed2e49 817 rp->Name(); //If photon sets the particle pdg name to gamma
e747b8da 818 rp->SetProductionVertex(0,0,0,0);
819 rp->SetFirstMother(-1);
820 rp->SetLastMother(-1);
821 rp->SetFirstDaughter(-1);
822 rp->SetLastDaughter(-1);
823 rp->SetPolarisation(0,0,0);
6ad0bfa0 824 index++ ;
825 }
6ad0bfa0 826}
e3817e5f 827
09fc14a0 828//____________________________________________________________________________
1f0e7ccd 829void AliPHOSPIDv1::Print()
09fc14a0 830{
b2a60966 831 // Print the parameters used for the particle type identification
bc0c084c 832
21cd0c07 833 TString message ;
834 message = "\n=============== AliPHOSPID1 ================\n" ;
835 message += "Making PID\n";
836 message += " Pricipal analysis file from 0.5 to 100 %s\n" ;
837 message += " Name of parameters file %s\n" ;
1f0e7ccd 838 message += " Matrix of Parameters: 14x4\n" ;
839 message += " Energy Calibration 1x3 [3 parametres to calibrate energy: A + B* E + C * E^2]\n" ;
bc0c084c 840 message += " RCPV 2x3 rows x and z, columns function cut parameters\n" ;
841 message += " TOF 1x3 [High Eff-Low Pur,Medium Eff-Pur, Low Eff-High Pur]\n" ;
21cd0c07 842 message += " PCA 5x4 [5 ellipse parametres and 4 parametres to calculate them: A/Sqrt(E) + B* E + C * E^2 + D]\n" ;
1f0e7ccd 843 message += " Pi0 PCA 5x3 [5 ellipse parametres and 3 parametres to calculate them: A + B* E + C * E^2]\n" ;
e3817e5f 844 Info("Print", message.Data(), fFileNamePrincipalPhoton.Data(), fFileNameParameters.Data() ) ;
50739f15 845 fParameters->Print() ;
09fc14a0 846}
847
848//____________________________________________________________________________
7b7c1533 849void AliPHOSPIDv1::WriteRecParticles(Int_t event)
09fc14a0 850{
e3817e5f 851
7b7c1533 852 AliPHOSGetter *gime = AliPHOSGetter::GetInstance() ;
a496c46c 853
fbf811ec 854 TClonesArray * recParticles = gime->RecParticles() ;
bf8f1fbd 855 recParticles->Expand(recParticles->GetEntriesFast() ) ;
fbf811ec 856 TTree * treeR ;
857
858 if(fToSplit){
859 if(!fSplitFile)
860 return ;
861 fSplitFile->cd() ;
862 char name[10] ;
863 sprintf(name,"%s%d", "TreeR",event) ;
864 treeR = dynamic_cast<TTree*>(fSplitFile->Get(name));
865 }
866 else{
867 treeR = gAlice->TreeR();
868 }
869
870 if(!treeR){
871 gAlice->MakeTree("R", fSplitFile);
872 treeR = gAlice->TreeR() ;
873 }
7acf6008 874
875 //First rp
876 Int_t bufferSize = 32000 ;
8d0f3f77 877 TBranch * rpBranch = treeR->Branch("PHOSRP",&recParticles,bufferSize);
fbf811ec 878 rpBranch->SetTitle(BranchName());
8d0f3f77 879
9688c1dd 880
7acf6008 881 //second, pid
882 Int_t splitlevel = 0 ;
883 AliPHOSPIDv1 * pid = this ;
8d0f3f77 884 TBranch * pidBranch = treeR->Branch("AliPHOSPID","AliPHOSPIDv1",&pid,bufferSize,splitlevel);
fbf811ec 885 pidBranch->SetTitle(BranchName());
7acf6008 886
761e34c0 887 rpBranch->Fill() ;
a496c46c 888 pidBranch->Fill() ;
9688c1dd 889
fbf811ec 890 treeR->AutoSave() ; //Write(0,kOverwrite) ;
891 if(gAlice->TreeR()!=treeR){
892 treeR->Delete();
893 }
7acf6008 894}
69183710 895
7acf6008 896//____________________________________________________________________________
9688c1dd 897TVector3 AliPHOSPIDv1::GetMomentumDirection(AliPHOSEmcRecPoint * emc, AliPHOSRecPoint * cpv)const
7acf6008 898{
899 // Calculates the momentum direction:
900 // 1. if only a EMC RecPoint, direction is given by IP and this RecPoint
9688c1dd 901 // 2. if a EMC RecPoint and CPV RecPoint, direction is given by the line through the 2 recpoints
902 // However because of the poor position resolution of PPSD the direction is always taken as if we were
7acf6008 903 // in case 1.
904
905 TVector3 dir(0,0,0) ;
906
907 TVector3 emcglobalpos ;
908 TMatrix dummy ;
909
910 emc->GetGlobalPosition(emcglobalpos, dummy) ;
e3817e5f 911
148b2bba 912
7acf6008 913 dir = emcglobalpos ;
914 dir.SetZ( -dir.Z() ) ; // why ?
e3817e5f 915 dir.SetMag(1.) ;
7acf6008 916
917 //account correction to the position of IP
918 Float_t xo,yo,zo ; //Coordinates of the origin
919 gAlice->Generator()->GetOrigin(xo,yo,zo) ;
920 TVector3 origin(xo,yo,zo);
921 dir = dir - origin ;
e3817e5f 922
7acf6008 923 return dir ;
924}
925//____________________________________________________________________________
a4e98857 926void AliPHOSPIDv1::PrintRecParticles(Option_t * option)
927{
dd5c4038 928 // Print table of reconstructed particles
929
7b7c1533 930 AliPHOSGetter *gime = AliPHOSGetter::GetInstance() ;
bf8f1fbd 931
a496c46c 932 TClonesArray * recParticles = gime->RecParticles(BranchName()) ;
21cd0c07 933
934 TString message ;
3bf72d32 935 message = "\nevent " ;
936 message += gAlice->GetEvNumber() ;
937 message += " found " ;
938 message += recParticles->GetEntriesFast();
939 message += " RecParticles\n" ;
940
7acf6008 941 if(strstr(option,"all")) { // printing found TS
3bf72d32 942 message += "\n PARTICLE Index \n" ;
7acf6008 943
944 Int_t index ;
7b7c1533 945 for (index = 0 ; index < recParticles->GetEntries() ; index++) {
21cd0c07 946 AliPHOSRecParticle * rp = (AliPHOSRecParticle * ) recParticles->At(index) ;
3bf72d32 947 message += "\n" ;
948 message += rp->Name().Data() ;
949 message += " " ;
950 message += rp->GetIndexInList() ;
951 message += " " ;
952 message += rp->GetType() ;
7acf6008 953 }
3bf72d32 954 }
955 Info("Print", message.Data() ) ;
69183710 956}