PHOS/AliPHOSPIDv1.cxx

   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
  16 /* $Id$ */
  17
  18 //_________________________________________________________________________
  19 // Implementation version v1 of the PHOS particle identifier
  20 // Particle identification based on the
  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
  31 //      inside an ellipse defined by fX_center, fY_center, fA, fB, fAngle
  32 //      (each bit corresponds to a different efficiency-purity point of the
  33 //      photon identification)
  34 //      A calibrated energy is calculated. The energy of the reconstructed
  35 //      cluster is corrected with the formula A + B * E  + C * E^2, whose parameters
  36 //      where obtained thourgh the study of the reconstructed energy
  37 //      distribution of monoenergetic photons.
  38 //
  39 //
  40 //
  41 // use case:
  42 //  root [0] AliPHOSPIDv1 * p = new AliPHOSPIDv1("galice1.root","v1")
  43 //  Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
  44 //                // reading headers from file galice1.root and create  RecParticles with title v1
  45                   // TrackSegments and RecPoints with title "v1" are used
  46 //                // set file name for the branch RecParticles
  47 //  root [1] p->ExecuteTask("deb all time")
  48 //                // available options
  49 //                // "deb" - prints # of reconstructed particles
  50 //                // "deb all" -  prints # and list of RecParticles
  51 //                // "time" - prints benchmarking results
  52 //
  53 //  root [2] AliPHOSPIDv1 * p2 = new AliPHOSPIDv1("galice1.root","v1","v0")
  54 //  Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
  55 //                // reading headers from file galice1.root and create  RecParticles with title v1
  56                   // RecPoints and TrackSegments with title "v0" are used
  57 //  root [3] p2->ExecuteTask()
  58 //
  59 //  There are two possible principal files available to do the analysis.
  60 //  One for energy ranges from 0.5 to 5 GeV, and another
  61 //  one from 5 to 100 GeV. This files are automatically called in function
  62 //  of the cluster energy.
  63
  64 //*-- Author: Yves Schutz (SUBATECH)  & Gines Martinez (SUBATECH) &
  65 //            Gustavo Conesa April 2002
  66
  67 // --- ROOT system ---
  68 #include "TROOT.h"
  69 #include "TTree.h"
  70 #include "TFile.h"
  71 #include "TF2.h"
  72 #include "TFormula.h"
  73 #include "TCanvas.h"
  74 #include "TFolder.h"
  75 #include "TSystem.h"
  76 #include "TBenchmark.h"
  77 #include "TMatrixD.h"
  78 #include "TPrincipal.h"
  79 #include "TSystem.h"
  80
  81 // --- Standard library ---
  82
  83 #include <iostream>
  84 #include <fstream>
  85 #include <iomanip>
  86
  87 // --- AliRoot header files ---
  88
  89 #include "AliRun.h"
  90 #include "AliGenerator.h"
  91 #include "AliPHOS.h"
  92 #include "AliPHOSPIDv1.h"
  93 #include "AliPHOSClusterizerv1.h"
  94 #include "AliPHOSTrackSegment.h"
  95 #include "AliPHOSTrackSegmentMakerv1.h"
  96 #include "AliPHOSRecParticle.h"
  97 #include "AliPHOSGeometry.h"
  98 #include "AliPHOSGetter.h"
  99
 100 ClassImp( AliPHOSPIDv1)
 101
 102 //____________________________________________________________________________
 103 AliPHOSPIDv1::AliPHOSPIDv1():AliPHOSPID()
 104 {
 105   // default ctor
 106
 107   InitParameters() ;
 108   fDefaultInit = kTRUE ;
 109
 110 }
 111
 112 //____________________________________________________________________________
 113 AliPHOSPIDv1::AliPHOSPIDv1(const char * headerFile,const char * name, const Bool_t toSplit)
 114 :AliPHOSPID(headerFile, name,toSplit)
 115 {
 116   //ctor with the indication on where to look for the track segments
 117
 118   InitParameters() ;
 119
 120   Init() ;
 121   fDefaultInit = kFALSE ;
 122
 123 }
 124
 125 //____________________________________________________________________________
 126 AliPHOSPIDv1::~AliPHOSPIDv1()
 127 {
 128   // dtor
 129   // fDefaultInit = kTRUE if PID created by default ctor (to get just the parameters)
 130
 131   delete [] fX ; // Principal input
 132   delete [] fP ; // Principal components
 133 //  delete fParameters ; // Matrix of Parameters
 134 //  delete fParameters5 ; // Matrix of Parameters
 135 //  delete fParameters100 ; // Matrix of Parameters
 136
 137
 138   if (!fDefaultInit) {
 139 //    AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
 140     // remove the task from the folder list
 141 //     gime->RemoveTask("P",GetName()) ;
 142 //     TString name(GetName()) ;
 143 //     name.ReplaceAll("pid", "clu") ;
 144 //     gime->RemoveTask("C",name) ;
 145
 146 //     // remove the data from the folder list
 147 //     name = GetName() ;
 148 //     name.Remove(name.Index(":")) ;
 149 //     gime->RemoveObjects("RE", name) ; // EMCARecPoints
 150 //     gime->RemoveObjects("RC", name) ; // CPVRecPoints
 151 //     gime->RemoveObjects("T", name) ;  // TrackSegments
 152 //     gime->RemoveObjects("P", name) ;  // RecParticles
 153
 154 //     // Delete gAlice
 155 //     gime->CloseFile() ;
 156
 157     fSplitFile = 0 ;
 158   }
 159 }
 160
 161 //____________________________________________________________________________
 162 const TString AliPHOSPIDv1::BranchName() const
 163 {
 164   TString branchName(GetName() ) ;
 165   branchName.Remove(branchName.Index(Version())-1) ;
 166   return branchName ;
 167 }
 168
 169 //____________________________________________________________________________
 170 void AliPHOSPIDv1::Init()
 171 {
 172   // Make all memory allocations that are not possible in default constructor
 173   // Add the PID task to the list of PHOS tasks
 174
 175   if ( strcmp(GetTitle(), "") == 0 )
 176     SetTitle("galice.root") ;
 177
 178   TString branchname(GetName()) ;
 179   branchname.Remove(branchname.Index(Version())-1) ;
 180   AliPHOSGetter * gime = AliPHOSGetter::GetInstance(GetTitle(),branchname.Data(),fToSplit ) ;
 181
 182   //  gime->SetRecParticlesTitle(BranchName()) ;
 183   if ( gime == 0 ) {
 184     cerr << "ERROR: AliPHOSPIDv1::Init -> Could not obtain the Getter object !" << endl ;
 185     return ;
 186   }
 187
 188   fSplitFile = 0 ;
 189   if(fToSplit){
 190     //First - extract full path if necessary
 191     TString fileName(GetTitle()) ;
 192     Ssiz_t islash = fileName.Last('/') ;
 193     if(islash<fileName.Length())
 194       fileName.Remove(islash+1,fileName.Length()) ;
 195     else
 196       fileName="" ;
 197     fileName+="PHOS.RecData." ;
 198     if((strcmp(branchname.Data(),"Default")!=0)&&(strcmp(branchname.Data(),"")!=0)){
 199       fileName+=branchname.Data() ;
 200       fileName+="." ;
 201     }
 202     fileName+="root" ;
 203     fSplitFile = static_cast<TFile*>(gROOT->GetFile(fileName.Data()));
 204     if(!fSplitFile)
 205       fSplitFile =  TFile::Open(fileName.Data(),"update") ;
 206   }
 207
 208   gime->PostPID(this) ;
 209   // create a folder on the white board //YSAlice/WhiteBoard/RecParticles/PHOS/recparticlesName
 210   gime->PostRecParticles(branchname) ;
 211
 212 }
 213
 214 //____________________________________________________________________________
 215 void AliPHOSPIDv1::InitParameters()
 216 {
 217 //   fFrom               = "" ;
 218 //   fHeaderFileName     = GetTitle() ;
 219 //   TString name(GetName()) ;
 220 //   if (name.IsNull())
 221 //     name = "Default" ;
 222 //   fTrackSegmentsTitle = name ;
 223 //   fRecPointsTitle     = name ;
 224 //   fRecParticlesTitle  = name ;
 225 //   name.Append(":") ;
 226 //   name.Append(Version()) ;
 227 //   SetName(name) ;
 228   fRecParticlesInRun = 0 ;
 229   fNEvent            = 0 ;
 230   //  fClusterizer       = 0 ;
 231   //  fTSMaker           = 0 ;
 232   fRecParticlesInRun = 0 ;
 233   TString pidName( GetName()) ;
 234   if (pidName.IsNull() )
 235     pidName = "Default" ;
 236   pidName.Append(":") ;
 237   pidName.Append(Version()) ;
 238   SetName(pidName) ;
 239   SetParameters() ; // fill the parameters matrix from parameters file
 240 }
 241
 242 //____________________________________________________________________________
 243 Double_t  AliPHOSPIDv1::GetCpvtoEmcDistanceCut(const Float_t Cluster_En, const TString Eff_Pur)
 244 {
 245   // Get CpvtoEmcDistanceCut parameter depending on the cluster energy and
 246   // Purity-Efficiency point (possible options "HIGH EFFICIENCY"
 247   // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing
 248   // EFFICIENCY by PURITY)
 249
 250   Int_t eff_pur = GetEffPurOption(Eff_Pur);
 251
 252   GetAnalysisParameters(Cluster_En) ;
 253   if((fClusterrcpv!= -1)&&(eff_pur != -1))
 254     return (*fParameters)(fClusterrcpv,eff_pur) ;
 255   else
 256     return 0.0;
 257 }
 258 //____________________________________________________________________________
 259
 260 Double_t  AliPHOSPIDv1::GetTimeGate(const Float_t Cluster_En, const TString Eff_Pur)
 261 {
 262   // Get TimeGate parameter depending on the cluster energy and
 263   // Purity-Efficiency point (possible options "HIGH EFFICIENCY"
 264   // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing
 265   // EFFICIENCY by PURITY)
 266
 267   Int_t eff_pur = GetEffPurOption(Eff_Pur);
 268   GetAnalysisParameters(Cluster_En) ;
 269
 270   if((fCluster!= -1)&&(eff_pur != -1))
 271     return (*fParameters)(fCluster+3+fMatrixExtraRow,eff_pur) ;
 272   else
 273     return 0.0;
 274
 275 }
 276 //_____________________________________________________________________________
 277 Float_t  AliPHOSPIDv1::GetDistance(AliPHOSEmcRecPoint * emc,AliPHOSRecPoint * cpv, Option_t *  Axis)const
 278 {
 279   // Calculates the distance between the EMC RecPoint and the PPSD RecPoint
 280
 281   const AliPHOSGeometry * geom = AliPHOSGetter::GetInstance()->PHOSGeometry() ;
 282   TVector3 vecEmc ;
 283   TVector3 vecCpv ;
 284   if(cpv){
 285     emc->GetLocalPosition(vecEmc) ;
 286     cpv->GetLocalPosition(vecCpv) ;
 287     if(emc->GetPHOSMod() == cpv->GetPHOSMod()){
 288       // Correct to difference in CPV and EMC position due to different distance to center.
 289       // we assume, that particle moves from center
 290       Float_t dCPV = geom->GetIPtoOuterCoverDistance();
 291       Float_t dEMC = geom->GetIPtoCrystalSurface() ;
 292       dEMC         = dEMC / dCPV ;
 293       vecCpv = dEMC * vecCpv  - vecEmc ;
 294       if (Axis == "X") return vecCpv.X();
 295       if (Axis == "Y") return vecCpv.Y();
 296       if (Axis == "Z") return vecCpv.Z();
 297       if (Axis == "R") return vecCpv.Mag();
 298   }
 299
 300     return 100000000 ;
 301   }
 302   return 100000000 ;
 303 }
 304
 305 //____________________________________________________________________________
 306 Double_t  AliPHOSPIDv1::CalibratedEnergy(Float_t e){
 307   //It calibrates Energy depending on the recpoint energy.
 308 //      The energy of the reconstructed
 309 //      cluster is corrected with the formula A + B* E  + C* E^2, whose parameters
 310 //      where obtained through the study of the reconstructed energy
 311 //      distribution of monoenergetic photons.
 312   Double_t enerec;
 313     enerec = fACalParameter + fBCalParameter * e+ fCCalParameter * e * e;
 314   return enerec ;
 315
 316 }
 317 //____________________________________________________________________________
 318 Int_t  AliPHOSPIDv1::GetPrincipalSign(Double_t* P, Int_t cluster, Int_t eff_pur)const
 319 {
 320   //This method gives if the PCA of the particle are inside a defined ellipse
 321   // Get the parameters that define the ellipse stored in the
 322   // fParameters matrix.
 323   Double_t X_center = (*fParameters)(cluster+6,eff_pur) ;
 324   Double_t Y_center = (*fParameters)(cluster+9,eff_pur) ;
 325   Double_t A        = (*fParameters)(cluster+12,eff_pur) ;
 326   Double_t B        = (*fParameters)(cluster+15,eff_pur) ;
 327   Double_t Angle    = (*fParameters)(cluster+18,eff_pur) ;
 328
 329   Int_t      prinsign;
 330   Double_t   Dx        = 0. ;
 331   Double_t   Delta     = 0. ;
 332   Double_t   Y         = 0. ;
 333   Double_t   Y_1       = 0. ;
 334   Double_t   Y_2       = 0. ;
 335   Double_t   Pi        = TMath::Pi() ;
 336   Double_t   Cos_Theta = TMath::Cos(Pi*Angle/180.) ;
 337   Double_t   Sin_Theta = TMath::Sin(Pi*Angle/180.) ;
 338
 339   Dx = P[0] - X_center ;
 340   Delta = 4.*A*A*B*B* (A*A*Cos_Theta*Cos_Theta
 341                         + B*B*Sin_Theta*Sin_Theta - Dx*Dx) ;
 342   if (Delta < 0.)
 343     {prinsign=0;}
 344
 345   else if (Delta == 0.)
 346     {
 347       Y = Cos_Theta*Sin_Theta*(A*A - B*B)*Dx /
 348         (A*A*Cos_Theta*Cos_Theta + B*B*Sin_Theta*Sin_Theta) ;
 349       Y += Y_center ;
 350       if(P[1]==Y )
 351         {prinsign=1;}
 352       else
 353         {prinsign=0;}
 354     }
 355   else
 356     {
 357       Y_1 = (Cos_Theta*Sin_Theta*(A*A - B*B) *Dx +
 358              TMath::Sqrt(Delta)/2.)/(A*A*Cos_Theta*Cos_Theta +
 359                                      B*B*Sin_Theta*Sin_Theta) ;
 360       Y_2 = (Cos_Theta*Sin_Theta*(A*A - B*B) *Dx -
 361              TMath::Sqrt(Delta)/2.)/(A*A*Cos_Theta*Cos_Theta
 362                                       + B*B*Sin_Theta*Sin_Theta) ;
 363       Y_1 += Y_center ;
 364       Y_2 += Y_center ;
 365       if ((P[1]<=Y_1) && (P[1]>=Y_2))
 366         {prinsign=1;}
 367       else
 368         {prinsign=0;}
 369     }
 370   return prinsign;
 371 }
 372
 373 //____________________________________________________________________________
 374 void  AliPHOSPIDv1::SetEllipseParameters(Float_t Cluster_En, TString Eff_Pur, Float_t x, Float_t y,Float_t a, Float_t b,Float_t angle)
 375 {
 376
 377   // Set all ellipse parameters depending on the cluster energy and
 378   // Purity-Efficiency point (possible options "HIGH EFFICIENCY"
 379   // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing
 380   // EFFICIENCY by PURITY)
 381
 382   Int_t eff_pur = GetEffPurOption(Eff_Pur);
 383   GetAnalysisParameters(Cluster_En) ;
 384   if((fCluster!= -1)&&(eff_pur != -1)){
 385     (*fParameters)(fCluster+6 +fMatrixExtraRow,eff_pur) = x ;
 386     (*fParameters)(fCluster+9 +fMatrixExtraRow,eff_pur) = y ;
 387     (*fParameters)(fCluster+12+fMatrixExtraRow,eff_pur) = a ;
 388     (*fParameters)(fCluster+15+fMatrixExtraRow,eff_pur) = b ;
 389     (*fParameters)(fCluster+18+fMatrixExtraRow,eff_pur) = angle ;
 390   }
 391
 392 }
 393 //__________________________________________________________________________
 394 void  AliPHOSPIDv1::SetEllipseXCenter(Float_t Cluster_En, TString Eff_Pur, Float_t x)
 395 {
 396   // Set the ellipse parameter x_center depending on the custer energy and
 397   // Purity-Efficiency point (possible options "HIGH EFFICIENCY"
 398   // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing
 399   // EFFICIENCY by PURITY)
 400   Int_t eff_pur = GetEffPurOption(Eff_Pur);
 401   GetAnalysisParameters(Cluster_En) ;
 402   if((fCluster!= -1)&&(eff_pur != -1))
 403     (*fParameters)(fCluster+6+fMatrixExtraRow,eff_pur) = x ;
 404 }
 405 //_________________________________________________________________________
 406 void  AliPHOSPIDv1::SetEllipseYCenter(Float_t Cluster_En, TString Eff_Pur, Float_t y)
 407 {
 408
 409   // Set the ellipse parameter y_center depending on the cluster energy and
 410   // Purity-Efficiency point (possible options "HIGH EFFICIENCY"
 411   // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing
 412   // EFFICIENCY by PURITY)
 413
 414   Int_t eff_pur = GetEffPurOption(Eff_Pur);
 415   GetAnalysisParameters(Cluster_En) ;
 416   if((fCluster!= -1)&&(eff_pur != -1))
 417     (*fParameters)(fCluster+9+fMatrixExtraRow,eff_pur) = y ;
 418 }
 419 //_________________________________________________________________________
 420 void  AliPHOSPIDv1::SetEllipseAParameter(Float_t Cluster_En, TString Eff_Pur, Float_t a)
 421 {
 422   // Set the ellipse parameter a depending on the cluster energy and
 423   // Purity-Efficiency point (possible options "HIGH EFFICIENCY"
 424   // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing
 425   // EFFICIENCY by PURITY)
 426
 427   Int_t eff_pur = GetEffPurOption(Eff_Pur);
 428   GetAnalysisParameters(Cluster_En) ;
 429   if((fCluster!= -1)&&(eff_pur != -1))
 430     (*fParameters)(fCluster+12+fMatrixExtraRow,eff_pur) = a ;
 431 }
 432 //________________________________________________________________________
 433 void  AliPHOSPIDv1::SetEllipseBParameter(Float_t Cluster_En, TString Eff_Pur, Float_t b)
 434 {
 435   // Set the ellipse parameter b depending on the cluster energy and
 436   // Purity-Efficiency point (possible options "HIGH EFFICIENCY"
 437   // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing
 438   // EFFICIENCY by PURITY)
 439
 440   Int_t eff_pur = GetEffPurOption(Eff_Pur);
 441   GetAnalysisParameters(Cluster_En) ;
 442   if((fCluster!= -1)&&(eff_pur != -1))
 443     (*fParameters)(fCluster+15+fMatrixExtraRow,eff_pur) = b ;
 444 }
 445 //________________________________________________________________________
 446 void  AliPHOSPIDv1::SetEllipseAngle(Float_t Cluster_En, TString Eff_Pur, Float_t angle)
 447 {
 448
 449   // Set the ellipse parameter angle depending on the cluster energy and
 450   // Purity-Efficiency point (possible options "HIGH EFFICIENCY"
 451   // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing
 452   // EFFICIENCY by PURITY)
 453
 454   Int_t eff_pur = GetEffPurOption(Eff_Pur);
 455   GetAnalysisParameters(Cluster_En) ;
 456   if((fCluster!= -1)&&(eff_pur != -1))
 457     (*fParameters)(fCluster+18+fMatrixExtraRow,eff_pur) = angle ;
 458 }
 459 //_____________________________________________________________________________
 460 void  AliPHOSPIDv1::SetCpvtoEmcDistanceCut(Float_t Cluster_En, TString Eff_Pur, Float_t cut)
 461 {
 462
 463   // Set the parameter Cpvto EmcDistanceCut depending on the cluster energy and
 464   // Purity-Efficiency point (possible options "HIGH EFFICIENCY"
 465   // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing
 466   // EFFICIENCY by PURITY)
 467
 468
 469   Int_t eff_pur = GetEffPurOption(Eff_Pur);
 470   GetAnalysisParameters(Cluster_En) ;
 471   if((fClusterrcpv!= -1)&&(eff_pur != -1))
 472     (*fParameters)(fClusterrcpv,eff_pur) = cut ;
 473 }
 474 //_____________________________________________________________________________
 475 void  AliPHOSPIDv1::SetTimeGate(Float_t Cluster_En, TString Eff_Pur, Float_t gate)
 476 {
 477
 478   // Set the parameter TimeGate depending on the cluster energy and
 479   // Purity-Efficiency point (possible options "HIGH EFFICIENCY"
 480   // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing
 481   // EFFICIENCY by PURITY)
 482
 483   Int_t eff_pur = GetEffPurOption(Eff_Pur);
 484   GetAnalysisParameters(Cluster_En) ;
 485   if((fCluster!= -1)&&(eff_pur != -1))
 486     (*fParameters)(fCluster+3+fMatrixExtraRow,eff_pur) = gate ;
 487 }
 488 //_____________________________________________________________________________
 489 void  AliPHOSPIDv1::SetParameters()
 490                                   //TString OptFileName)
 491 {
 492   // PCA : To do the Principal Components Analysis it is necessary
 493   // the Principal file, which is opened here
 494   fX         = new double[7]; // Data for the PCA
 495   fP         = new double[7]; // Eigenvalues of the PCA
 496
 497
 498   // Set the principal and parameters files to be used
 499   fFileName5  = "$ALICE_ROOT/PHOS/PCA8pa15_0.5-5.root" ;
 500   fFileNamePar5 = gSystem->ExpandPathName("$ALICE_ROOT/PHOS/Parameters_0.5_5.dat");
 501   fFileName100  = "$ALICE_ROOT/PHOS/PCA8pa15_0.5-100.root" ;
 502   fFileNamePar100 = gSystem->ExpandPathName("$ALICE_ROOT/PHOS/Parameters_0.5_100.dat");
 503
 504   //SetPrincipalFileOptions();
 505   //fOptFileName);
 506   TFile f5( fFileName5.Data(), "read" ) ;
 507   fPrincipal5 = dynamic_cast<TPrincipal*> (f5.Get("principal")) ;
 508   f5.Close() ;
 509   TFile f100( fFileName100.Data(), "read" ) ;
 510   fPrincipal100 = dynamic_cast<TPrincipal*> (f100.Get("principal")) ;
 511   f100.Close() ;
 512   TFile f( fFileName100.Data(), "read" ) ;
 513   fPrincipal = dynamic_cast<TPrincipal*> (f.Get("principal")) ;
 514   f.Close() ;
 515   // Initialization of the Parameters matrix. In the File ParametersXX.dat
 516   // are all the parameters. These are introduced in a matrix of 21x3 or 22x3
 517   // elements (depending on the principal file 21 rows for 0.5-5 GeV and 22
 518   // rows for 5-100).
 519   // All the parameters defined in this file are, in order of row (there are
 520   // 3 rows per parameter): CpvtoEmcDistanceCut(if the principal file is 5-100
 521   // GeV then 4 rows), TimeGate and the ellipse parameters, X_center, Y_center,
 522   // a, b, angle. Each row of a given parameter depends on the cluster energy range
 523   // (wich depends on the chosen principal file)
 524   // Each column designs the parameters for a point in the Efficiency-Purity
 525   // of the photon identification P1(96%,63%), P2(87%,0.88%) and P3(68%,94%)
 526   // for the principal file from 0.5-5 GeV and for the other one P1(95%,79%),
 527   // P2(89%,90%) and P3(72%,96%)
 528
 529   fEnergyAnalysisCut = 5.; // Energy cut to change PCA
 530
 531   fParameters5 = new TMatrixD(21,3) ;
 532   fParameters100 = new TMatrixD(22,3) ;
 533   fParameters = new TMatrixD(22,3) ;
 534
 535   ifstream paramFile5(fFileNamePar5) ;
 536
 537   Int_t i,j ;
 538
 539   for(i = 0; i< 21; i++){
 540     for(j = 0; j< 3; j++){
 541       paramFile5 >> (*fParameters5)(i,j) ;
 542     }
 543   }
 544   paramFile5.close();
 545
 546   ifstream paramFile100(fFileNamePar100) ;
 547
 548   Int_t l,k ;
 549
 550   for(l = 0; l< 22; l++){
 551     for(k = 0; k< 3; k++){
 552       paramFile100 >> (*fParameters100)(l,k) ;
 553     }
 554   }
 555   paramFile100.close();
 556
 557   ifstream paramFile(fFileNamePar100) ;
 558   Int_t h,n;
 559   for(h = 0; h< 22; h++){
 560     for(n = 0; n< 3; n++){
 561       paramFile >> (*fParameters)(h,n) ;
 562     }
 563   }
 564   paramFile.close();
 565
 566   fCluster = -1;
 567   fClusterrcpv = -1;
 568   fMatrixExtraRow = 0;
 569
 570   //Calibration parameters Encal = C * E^2 + B * E + A  (E is the energy from cluster)
 571   fACalParameter = 0.0241  ;
 572   fBCalParameter = 1.0504  ;
 573   fCCalParameter = 0.000249 ;
 574
 575   // fParameters->Print();
 576 }
 577 //_____________________________________________________________________________
 578 void  AliPHOSPIDv1::GetAnalysisParameters(Float_t Cluster_En)
 579 {
 580   if(Cluster_En <=  fEnergyAnalysisCut){
 581     fPrincipal  = fPrincipal5;
 582     fParameters = fParameters5;
 583     fMatrixExtraRow = 0;
 584     GetClusterOption(Cluster_En,kFALSE) ;
 585   }
 586   else{
 587     fPrincipal  = fPrincipal100;
 588     fParameters = fParameters100;
 589     fMatrixExtraRow = 1;
 590     GetClusterOption(Cluster_En,kTRUE) ;
 591   }
 592 }
 593
 594 //_____________________________________________________________________________
 595 void  AliPHOSPIDv1::GetClusterOption(const Float_t Cluster_En, const Bool_t range)
 596 {
 597
 598   // Gives the cluster energy range.
 599   // range = kFALSE Default analysis range from 0.5 to 5 GeV
 600   // range = kTRUE  analysis range from 0.5 to 100 GeV
 601
 602
 603   //Int_t cluster = -1 ;
 604
 605   if((range == kFALSE)){
 606     if((Cluster_En > 0.3)&&(Cluster_En <= 1.0)){
 607       fCluster = 0 ;
 608       fClusterrcpv = 0 ;
 609     }
 610     if((Cluster_En > 1.0)&&(Cluster_En <= 2.0)){
 611       fCluster = 1 ;
 612       fClusterrcpv = 1 ;
 613     }
 614     if( Cluster_En > 2.0){
 615       fCluster = 2 ;
 616       fClusterrcpv = 2 ;
 617     }
 618   }
 619   else if(range == kTRUE){
 620     if((Cluster_En > 0.5 )&&(Cluster_En <= 20.0)) fCluster = 0 ;
 621     if((Cluster_En > 20.0)&&(Cluster_En <= 50.0)) fCluster = 1 ;
 622     if( Cluster_En > 50.0)                        fCluster = 2 ;
 623     if((Cluster_En > 5.0 )&&(Cluster_En <= 10.0)) fClusterrcpv = 0 ;
 624     if((Cluster_En > 10.0)&&(Cluster_En <= 20.0)) fClusterrcpv = 1 ;
 625     if((Cluster_En > 20.0)&&(Cluster_En <= 30.0)) fClusterrcpv = 2 ;
 626     if( Cluster_En > 30.0)                        fClusterrcpv = 3 ;
 627   }
 628   else {
 629     fCluster = -1 ;
 630     fClusterrcpv = -1;
 631     cout<<"Invalid Energy option"<<endl;
 632   }
 633
 634   //return cluster;
 635 }
 636 //____________________________________________________________________________
 637 Int_t  AliPHOSPIDv1::GetEffPurOption(const TString Eff_Pur) const
 638 {
 639
 640   // Looks for the Purity-Efficiency point (possible options "HIGH EFFICIENCY"
 641   // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing
 642   // EFFICIENCY by PURITY)
 643
 644   Int_t eff_pur = -1 ;
 645
 646   if(Eff_Pur.Contains("HIGH EFFICIENCY") ||Eff_Pur.Contains("LOW PURITY") )
 647     eff_pur = 0 ;
 648   else if(Eff_Pur.Contains("MEDIUM EFFICIENCY") ||Eff_Pur.Contains("MEDIUM PURITY") )
 649     eff_pur = 1 ;
 650   else if(Eff_Pur.Contains("LOW EFFICIENCY")||Eff_Pur.Contains("HIGH PURITY") )
 651     eff_pur = 2 ;
 652   else{
 653     eff_pur = -1;
 654     cout<<"Invalid Efficiency-Purity option"<<endl;
 655     cout<<"Possible options: HIGH EFFICIENCY =    LOW PURITY"<<endl;
 656     cout<<"                MEDIUM EFFICIENCY = MEDIUM PURITY"<<endl;
 657     cout<<"                   LOW EFFICIENCY =   HIGH PURITY"<<endl;
 658   }
 659
 660   return eff_pur;
 661 }
 662 //____________________________________________________________________________
 663
 664 void  AliPHOSPIDv1::Exec(Option_t * option)
 665 {
 666   //Steering method
 667
 668   if( strcmp(GetName(), "")== 0 )
 669     Init() ;
 670
 671   if(strstr(option,"tim"))
 672     gBenchmark->Start("PHOSPID");
 673
 674   if(strstr(option,"print")) {
 675     Print("") ;
 676     return ;
 677   }
 678
 679
 680 //   gAlice->GetEvent(0) ;
 681
 682 //   //check, if the branch with name of this" already exits?
 683 //   if (gAlice->TreeR()) {
 684 //     TObjArray * lob = (TObjArray*)gAlice->TreeR()->GetListOfBranches() ;
 685 //     TIter next(lob) ;
 686 //     TBranch * branch = 0 ;
 687 //     Bool_t phospidfound = kFALSE, pidfound = kFALSE ;
 688
 689 //     TString taskName(GetName()) ;
 690 //     taskName.Remove(taskName.Index(Version())-1) ;
 691
 692 //     while ( (branch = (TBranch*)next()) && (!phospidfound || !pidfound) ) {
 693 //       if ( (strcmp(branch->GetName(), "PHOSPID")==0) && (strcmp(branch->GetTitle(), taskName.Data())==0) )
 694 //      phospidfound = kTRUE ;
 695
 696 //       else if ( (strcmp(branch->GetName(), "AliPHOSPID")==0) && (strcmp(branch->GetTitle(), taskName.Data())==0) )
 697 //      pidfound = kTRUE ;
 698 //     }
 699
 700 //     if ( phospidfound || pidfound ) {
 701 //       cerr << "WARNING: AliPHOSPIDv1::Exec -> RecParticles and/or PIDtMaker branch with name "
 702 //         << taskName.Data() << " already exits" << endl ;
 703 //       return ;
 704 //     }
 705 //   }
 706
 707 //   Int_t nevents = (Int_t) gAlice->TreeE()->GetEntries() ;
 708 //   Int_t ievent ;
 709 //   AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
 710
 711   AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
 712   if(gime->BranchExists("RecParticles") )
 713     return ;
 714   Int_t nevents = gime->MaxEvent() ;       //(Int_t) gAlice->TreeE()->GetEntries() ;
 715   Int_t ievent ;
 716
 717
 718   for(ievent = 0; ievent < nevents; ievent++){
 719     gime->Event(ievent,"R") ;
 720
 721     MakeRecParticles() ;
 722
 723     WriteRecParticles(ievent);
 724
 725     if(strstr(option,"deb"))
 726       PrintRecParticles(option) ;
 727
 728     //increment the total number of rec particles per run
 729     fRecParticlesInRun += gime->RecParticles(BranchName())->GetEntriesFast() ;
 730
 731   }
 732
 733   if(strstr(option,"tim")){
 734     gBenchmark->Stop("PHOSPID");
 735     cout << "AliPHOSPID:" << endl ;
 736     cout << "  took " << gBenchmark->GetCpuTime("PHOSPID") << " seconds for PID "
 737          <<  gBenchmark->GetCpuTime("PHOSPID")/nevents << " seconds per event " << endl ;
 738     cout << endl ;
 739   }
 740
 741 }
 742
 743 //____________________________________________________________________________
 744 void  AliPHOSPIDv1::MakeRecParticles(){
 745
 746   // Makes a RecParticle out of a TrackSegment
 747
 748   AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
 749   TObjArray * emcRecPoints = gime->EmcRecPoints() ;
 750   TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
 751   TClonesArray * trackSegments = gime->TrackSegments() ;
 752   if ( !emcRecPoints || !cpvRecPoints || !trackSegments ) {
 753     cerr << "ERROR:  AliPHOSPIDv1::MakeRecParticles -> RecPoints or TrackSegments not found ! " << endl ;
 754     abort() ;
 755   }
 756   TClonesArray * recParticles  = gime->RecParticles() ;
 757   recParticles->Clear();
 758
 759
 760   TIter next(trackSegments) ;
 761   AliPHOSTrackSegment * ts ;
 762   Int_t index = 0 ;
 763   AliPHOSRecParticle * rp ;
 764   while ( (ts = (AliPHOSTrackSegment *)next()) ) {
 765
 766     new( (*recParticles)[index] ) AliPHOSRecParticle() ;
 767     rp = (AliPHOSRecParticle *)recParticles->At(index) ;
 768     rp->SetTrackSegment(index) ;
 769     rp->SetIndexInList(index) ;
 770
 771     AliPHOSEmcRecPoint * emc = 0 ;
 772     if(ts->GetEmcIndex()>=0)
 773       emc = (AliPHOSEmcRecPoint *) emcRecPoints->At(ts->GetEmcIndex()) ;
 774
 775     AliPHOSRecPoint    * cpv = 0 ;
 776     if(ts->GetCpvIndex()>=0)
 777       cpv = (AliPHOSRecPoint *)   cpvRecPoints->At(ts->GetCpvIndex()) ;
 778
 779     // Now set type (reconstructed) of the particle
 780
 781     // Choose the cluster energy range
 782
 783     // YK: check if (emc != 0) !!!
 784     if (!emc) {
 785       cerr << "ERROR:  AliPHOSPIDv1::MakeRecParticles -> emc("
 786            <<ts->GetEmcIndex()<<") = "         <<emc<< endl;
 787       abort();
 788     }
 789     Float_t    e = emc->GetEnergy() ;
 790
 791     GetAnalysisParameters(e);// Gives value to fCluster, fClusterrcpv, fMatrixExtraRow, and to fPrincipal and fParameters depending on the energy.
 792
 793     if((fCluster== -1)||(fClusterrcpv == -1)) continue ;
 794
 795     Float_t  lambda[2] ;
 796     emc->GetElipsAxis(lambda) ;
 797     Float_t time =emc->GetTime() ;
 798
 799     if((lambda[0]>0.01) && (lambda[1]>0.01) && time > 0.){
 800
 801       // Loop of Efficiency-Purity (the 3 points of purity or efficiency are taken
 802       // into account to set the particle identification)
 803       for(Int_t eff_pur = 0; eff_pur < 3 ; eff_pur++){
 804
 805         // Looking at the CPV detector. If RCPV greater than CpvEmcDistance, 1st,
 806         // 2nd or 3rd bit (depending on the efficiency-purity point )is set to 1 .
 807
 808         if(GetDistance(emc, cpv,  "R") > (*fParameters)(fClusterrcpv,eff_pur) )
 809           rp->SetPIDBit(eff_pur) ;
 810
 811         // Looking the TOF. If TOF smaller than gate,  4th, 5th or 6th
 812         // bit (depending on the efficiency-purity point )is set to 1
 813         if(time< (*fParameters)(fCluster+3+fMatrixExtraRow,eff_pur))
 814           rp->SetPIDBit(eff_pur+3) ;
 815
 816         // Looking PCA. Define and calculate the data (X), introduce in the function
 817         // X2P that gives the components (P).
 818         Float_t  Spher = 0. ;
 819         Float_t  Emaxdtotal = 0. ;
 820
 821         if((lambda[0]+lambda[1])!=0) Spher=fabs(lambda[0]-lambda[1])/(lambda[0]+lambda[1]);
 822
 823         Emaxdtotal=emc->GetMaximalEnergy()/emc->GetEnergy();
 824
 825         fX[0] = lambda[0] ;
 826         fX[1] = lambda[1] ;
 827         fX[2] = emc->GetDispersion() ;
 828         fX[3] = Spher ;
 829         fX[4] = emc->GetMultiplicity() ;
 830         fX[5] = Emaxdtotal ;
 831         fX[6] = emc->GetCoreEnergy() ;
 832
 833         fPrincipal->X2P(fX,fP);
 834
 835         //If we are inside the ellipse, 7th, 8th or 9th
 836         // bit (depending on the efficiency-purity point )is set to 1
 837         if(GetPrincipalSign(fP,fCluster+fMatrixExtraRow,eff_pur) == 1)
 838           rp->SetPIDBit(eff_pur+6) ;
 839
 840       }
 841     }
 842
 843     //Set momentum, energy and other parameters
 844     Float_t  encal = CalibratedEnergy(e);
 845     TVector3 dir   = GetMomentumDirection(emc,cpv) ;
 846     dir.SetMag(encal) ;
 847     rp->SetMomentum(dir.X(),dir.Y(),dir.Z(),encal) ;
 848     rp->SetCalcMass(0);
 849     rp->Name(); //If photon sets the particle pdg name to gamma
 850     rp->SetProductionVertex(0,0,0,0);
 851     rp->SetFirstMother(-1);
 852     rp->SetLastMother(-1);
 853     rp->SetFirstDaughter(-1);
 854     rp->SetLastDaughter(-1);
 855     rp->SetPolarisation(0,0,0);
 856     index++ ;
 857   }
 858
 859 }
 860
 861 //____________________________________________________________________________
 862 void  AliPHOSPIDv1:: Print()
 863 {
 864   // Print the parameters used for the particle type identification
 865     cout <<  "=============== AliPHOSPID1 ================" << endl ;
 866     cout <<  "Making PID "<< endl ;
 867 //     cout <<  "    Headers file:               " << fHeaderFileName.Data() << endl ;
 868 //     cout <<  "    RecPoints branch title:     " << fRecPointsTitle.Data() << endl ;
 869 //     cout <<  "    TrackSegments Branch title: " << fTrackSegmentsTitle.Data() << endl ;
 870 //     cout <<  "    RecParticles Branch title   " << fRecParticlesTitle.Data() << endl;
 871
 872     cout <<  "    Pricipal analysis file from 0.5 to 5 " << fFileName5.Data() << endl;
 873     cout <<  "    Name of parameters file     "<<fFileNamePar5.Data() << endl ;
 874     cout <<  "    Matrix of Parameters: "<<endl;
 875     cout <<  "           3  Columns [High Eff-Low Pur,Medium Eff-Pur, Low Eff-High Pur]"<<endl;
 876     cout <<  "           21 Rows, each 3 [ RCPV, TOF, X_Center, Y_Center, A, B, Angle ]"<<endl;
 877     fParameters5->Print() ;
 878
 879     cout <<  "    Pricipal analysis file from 5 to 100 " << fFileName100.Data() << endl;
 880     cout <<  "    Name of parameters file     "<<fFileNamePar100.Data() << endl ;
 881     cout <<  "    Matrix of Parameters: "<<endl;
 882     cout <<  "           3  Columns [High Eff-Low Pur,Medium Eff-Pur, Low Eff-High Pur]"<<endl;
 883     cout <<  "           22 Rows, [ 4 RCPV, 3 TOF, 3 X_Center, 3 Y_Center, 3 A, 3 B, 3 Angle ]"<<endl;
 884     fParameters100->Print() ;
 885
 886     cout <<  "    Energy Calibration Parameters  A + B* E + C * E^2"<<endl;
 887     cout <<  "    E is the energy from the cluster "<<endl;
 888     cout <<  "           A = "<< fACalParameter  << endl;
 889     cout <<  "           B = "<< fBCalParameter  << endl;
 890     cout <<  "           C = "<< fCCalParameter  << endl;
 891     cout <<  "============================================" << endl ;
 892 }
 893
 894 //____________________________________________________________________________
 895 void  AliPHOSPIDv1::WriteRecParticles(Int_t event)
 896 {
 897
 898   AliPHOSGetter *gime = AliPHOSGetter::GetInstance() ;
 899
 900   TClonesArray * recParticles = gime->RecParticles() ;
 901   recParticles->Expand(recParticles->GetEntriesFast() ) ;
 902   TTree * treeR ;
 903
 904   if(fToSplit){
 905     if(!fSplitFile)
 906       return ;
 907     fSplitFile->cd() ;
 908     char name[10] ;
 909     sprintf(name,"%s%d", "TreeR",event) ;
 910     treeR = dynamic_cast<TTree*>(fSplitFile->Get(name));
 911   }
 912   else{
 913     treeR = gAlice->TreeR();
 914   }
 915
 916   if(!treeR){
 917     gAlice->MakeTree("R", fSplitFile);
 918     treeR = gAlice->TreeR() ;
 919   }
 920
 921   //First rp
 922   Int_t bufferSize = 32000 ;
 923   TBranch * rpBranch = treeR->Branch("PHOSRP",&recParticles,bufferSize);
 924   rpBranch->SetTitle(BranchName());
 925
 926
 927   //second, pid
 928   Int_t splitlevel = 0 ;
 929   AliPHOSPIDv1 * pid = this ;
 930   TBranch * pidBranch = treeR->Branch("AliPHOSPID","AliPHOSPIDv1",&pid,bufferSize,splitlevel);
 931   pidBranch->SetTitle(BranchName());
 932
 933   rpBranch->Fill() ;
 934   pidBranch->Fill() ;
 935
 936   treeR->AutoSave() ; //Write(0,kOverwrite) ;
 937   if(gAlice->TreeR()!=treeR){
 938     treeR->Delete();
 939   }
 940 }
 941
 942 //____________________________________________________________________________
 943 TVector3 AliPHOSPIDv1::GetMomentumDirection(AliPHOSEmcRecPoint * emc, AliPHOSRecPoint * cpv)const
 944 {
 945   // Calculates the momentum direction:
 946   //   1. if only a EMC RecPoint, direction is given by IP and this RecPoint
 947   //   2. if a EMC RecPoint and CPV RecPoint, direction is given by the line through the 2 recpoints
 948   //  However because of the poor position resolution of PPSD the direction is always taken as if we were
 949   //  in case 1.
 950
 951   TVector3 dir(0,0,0) ;
 952
 953   TVector3 emcglobalpos ;
 954   TMatrix  dummy ;
 955
 956   emc->GetGlobalPosition(emcglobalpos, dummy) ;
 957
 958
 959   dir = emcglobalpos ;
 960   dir.SetZ( -dir.Z() ) ;   // why ?
 961   dir.SetMag(1.) ;
 962
 963   //account correction to the position of IP
 964   Float_t xo,yo,zo ; //Coordinates of the origin
 965   gAlice->Generator()->GetOrigin(xo,yo,zo) ;
 966   TVector3 origin(xo,yo,zo);
 967   dir = dir - origin ;
 968
 969   return dir ;
 970 }
 971 //____________________________________________________________________________
 972 void AliPHOSPIDv1::PrintRecParticles(Option_t * option)
 973 {
 974   // Print table of reconstructed particles
 975
 976   AliPHOSGetter *gime = AliPHOSGetter::GetInstance() ;
 977
 978   TClonesArray * recParticles = gime->RecParticles(BranchName()) ;
 979
 980   cout << "AliPHOSPIDv1: event "<<gAlice->GetEvNumber()  << endl ;
 981   cout << "       found " << recParticles->GetEntriesFast() << " RecParticles " << endl ;
 982
 983   if(strstr(option,"all")) {  // printing found TS
 984
 985     cout << "  PARTICLE       "
 986          << "  Index    "  << endl ;
 987
 988     Int_t index ;
 989     for (index = 0 ; index < recParticles->GetEntries() ; index++) {
 990        AliPHOSRecParticle * rp = (AliPHOSRecParticle * ) recParticles->At(index) ;
 991
 992        cout << setw(10) << rp->Name() << "  "
 993             << setw(5) <<  rp->GetIndexInList() << " " <<endl;
 994        cout << "Type "<<  rp->GetType() << endl;
 995     }
 996     cout << "-------------------------------------------" << endl ;
 997   }
 998
 999 }
1000
1001
1002