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