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