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