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