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.h>
  84 #include <fstream.h>
  85 #include <iomanip.h>
  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 }
 109
 110 //____________________________________________________________________________
 111 AliPHOSPIDv1::AliPHOSPIDv1(const char * headerFile,const char * name, const char * from) : AliPHOSPID(headerFile, name)
 112
 113
 114 {
 115   //ctor with the indication on where to look for the track segments
 116
 117   InitParameters() ;
 118
 119   if ( from == 0 )
 120     fFrom = name ;
 121   else
 122     fFrom = from ;
 123
 124   Init() ;
 125
 126 }
 127
 128 //____________________________________________________________________________
 129 AliPHOSPIDv1::~AliPHOSPIDv1()
 130 {
 131   // dtor
 132   // gime=0 if PID created by default ctor (to get just the parameters)
 133
 134
 135   delete [] fX ; // Principal input
 136   delete [] fP ; // Principal components
 137   delete fParameters ; // Matrix of Parameters
 138   delete fParameters5 ; // Matrix of Parameters
 139   delete fParameters100 ; // Matrix of Parameters
 140
 141
 142   AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
 143
 144   if(gime){
 145   // remove the task from the folder list
 146   gime->RemoveTask("P",GetName()) ;
 147   TString name(GetName()) ;
 148   name.ReplaceAll("pid", "clu") ;
 149   gime->RemoveTask("C",name) ;
 150
 151   // remove the data from the folder list
 152   name = GetName() ;
 153   name.Remove(name.Index(":")) ;
 154   gime->RemoveObjects("RE", name) ; // EMCARecPoints
 155   gime->RemoveObjects("RC", name) ; // CPVRecPoints
 156   gime->RemoveObjects("T", name) ;  // TrackSegments
 157   gime->RemoveObjects("P", name) ;  // RecParticles
 158
 159   // Delete gAlice
 160   gime->CloseFile() ;
 161   fSplitFile = 0 ;
 162   }
 163 }
 164
 165 //____________________________________________________________________________
 166 const TString AliPHOSPIDv1::BranchName() const
 167 {
 168   TString branchName(GetName() ) ;
 169   branchName.Remove(branchName.Index(Version())-1) ;
 170   return branchName ;
 171 }
 172
 173 //____________________________________________________________________________
 174 void AliPHOSPIDv1::Init()
 175 {
 176   // Make all memory allocations that are not possible in default constructor
 177   // Add the PID task to the list of PHOS tasks
 178
 179   if ( strcmp(GetTitle(), "") == 0 )
 180     SetTitle("galice.root") ;
 181
 182   AliPHOSGetter * gime = AliPHOSGetter::GetInstance(GetTitle(), fFrom.Data()) ;
 183
 184   gime->SetRecParticlesTitle(BranchName()) ;
 185   if ( gime == 0 ) {
 186     cerr << "ERROR: AliPHOSPIDv1::Init -> Could not obtain the Getter object !" << endl ;
 187     return ;
 188   }
 189
 190   gime->PostPID(this) ;
 191   // create a folder on the white board //YSAlice/WhiteBoard/RecParticles/PHOS/recparticlesName
 192   gime->PostRecParticles(BranchName()) ;
 193
 194 }
 195
 196 //____________________________________________________________________________
 197 void AliPHOSPIDv1::InitParameters()
 198 {
 199   fFrom               = "" ;
 200   fHeaderFileName     = GetTitle() ;
 201   TString name(GetName()) ;
 202   if (name.IsNull())
 203     name = "Default" ;
 204   fTrackSegmentsTitle = name ;
 205   fRecPointsTitle     = name ;
 206   fRecParticlesTitle  = name ;
 207   name.Append(":") ;
 208   name.Append(Version()) ;
 209   SetName(name) ;
 210   fRecParticlesInRun = 0 ;
 211   fNEvent            = 0 ;
 212   fClusterizer       = 0 ;
 213   fTSMaker           = 0 ;
 214   fRecParticlesInRun = 0 ;
 215   SetParameters() ; // fill the parameters matrix from parameters file
 216 }
 217
 218 //____________________________________________________________________________
 219 Double_t  AliPHOSPIDv1::GetCpvtoEmcDistanceCut(const Float_t Cluster_En, const TString Eff_Pur)
 220 {
 221   // Get CpvtoEmcDistanceCut parameter depending on the cluster energy and
 222   // Purity-Efficiency point (possible options "HIGH EFFICIENCY"
 223   // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing
 224   // EFFICIENCY by PURITY)
 225
 226   Int_t eff_pur = GetEffPurOption(Eff_Pur);
 227
 228   GetAnalysisParameters(Cluster_En) ;
 229   if((fClusterrcpv!= -1)&&(eff_pur != -1))
 230     return (*fParameters)(fClusterrcpv,eff_pur) ;
 231   else
 232     return 0.0;
 233 }
 234 //____________________________________________________________________________
 235
 236 Double_t  AliPHOSPIDv1::GetTimeGate(const Float_t Cluster_En, const TString Eff_Pur)
 237 {
 238   // Get TimeGate parameter depending on the cluster energy and
 239   // Purity-Efficiency point (possible options "HIGH EFFICIENCY"
 240   // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing
 241   // EFFICIENCY by PURITY)
 242
 243   Int_t eff_pur = GetEffPurOption(Eff_Pur);
 244   GetAnalysisParameters(Cluster_En) ;
 245
 246   if((fCluster!= -1)&&(eff_pur != -1))
 247     return (*fParameters)(fCluster+3+fMatrixExtraRow,eff_pur) ;
 248   else
 249     return 0.0;
 250
 251 }
 252 //_____________________________________________________________________________
 253 Float_t  AliPHOSPIDv1::GetDistance(AliPHOSEmcRecPoint * emc,AliPHOSRecPoint * cpv, Option_t *  Axis)const
 254 {
 255   // Calculates the distance between the EMC RecPoint and the PPSD RecPoint
 256
 257   const AliPHOSGeometry * geom = AliPHOSGetter::GetInstance()->PHOSGeometry() ;
 258   TVector3 vecEmc ;
 259   TVector3 vecCpv ;
 260   if(cpv){
 261     emc->GetLocalPosition(vecEmc) ;
 262     cpv->GetLocalPosition(vecCpv) ;
 263     if(emc->GetPHOSMod() == cpv->GetPHOSMod()){
 264       // Correct to difference in CPV and EMC position due to different distance to center.
 265       // we assume, that particle moves from center
 266       Float_t dCPV = geom->GetIPtoOuterCoverDistance();
 267       Float_t dEMC = geom->GetIPtoCrystalSurface() ;
 268       dEMC         = dEMC / dCPV ;
 269       vecCpv = dEMC * vecCpv  - vecEmc ;
 270       if (Axis == "X") return vecCpv.X();
 271       if (Axis == "Y") return vecCpv.Y();
 272       if (Axis == "Z") return vecCpv.Z();
 273       if (Axis == "R") return vecCpv.Mag();
 274   }
 275
 276     return 100000000 ;
 277   }
 278   return 100000000 ;
 279 }
 280
 281 //____________________________________________________________________________
 282 Double_t  AliPHOSPIDv1::CalibratedEnergy(Float_t e){
 283   //It calibrates Energy depending on the recpoint energy.
 284 //      The energy of the reconstructed
 285 //      cluster is corrected with the formula A + B* E  + C* E^2, whose parameters
 286 //      where obtained through the study of the reconstructed energy
 287 //      distribution of monoenergetic photons.
 288   Double_t enerec;
 289     enerec = fACalParameter + fBCalParameter * e+ fCCalParameter * e * e;
 290   return enerec ;
 291
 292 }
 293 //____________________________________________________________________________
 294 Int_t  AliPHOSPIDv1::GetPrincipalSign(Double_t* P, Int_t cluster, Int_t eff_pur)const
 295 {
 296   //This method gives if the PCA of the particle are inside a defined ellipse
 297   // Get the parameters that define the ellipse stored in the
 298   // fParameters matrix.
 299   Double_t X_center = (*fParameters)(cluster+6,eff_pur) ;
 300   Double_t Y_center = (*fParameters)(cluster+9,eff_pur) ;
 301   Double_t A        = (*fParameters)(cluster+12,eff_pur) ;
 302   Double_t B        = (*fParameters)(cluster+15,eff_pur) ;
 303   Double_t Angle    = (*fParameters)(cluster+18,eff_pur) ;
 304
 305   Int_t      prinsign;
 306   Double_t   Dx        = 0. ;
 307   Double_t   Delta     = 0. ;
 308   Double_t   Y         = 0. ;
 309   Double_t   Y_1       = 0. ;
 310   Double_t   Y_2       = 0. ;
 311   Double_t   Pi        = TMath::Pi() ;
 312   Double_t   Cos_Theta = TMath::Cos(Pi*Angle/180.) ;
 313   Double_t   Sin_Theta = TMath::Sin(Pi*Angle/180.) ;
 314
 315   Dx = P[0] - X_center ;
 316   Delta = 4.*A*A*B*B* (A*A*Cos_Theta*Cos_Theta
 317                         + B*B*Sin_Theta*Sin_Theta - Dx*Dx) ;
 318   if (Delta < 0.)
 319     {prinsign=0;}
 320
 321   else if (Delta == 0.)
 322     {
 323       Y = Cos_Theta*Sin_Theta*(A*A - B*B)*Dx /
 324         (A*A*Cos_Theta*Cos_Theta + B*B*Sin_Theta*Sin_Theta) ;
 325       Y += Y_center ;
 326       if(P[1]==Y )
 327         {prinsign=1;}
 328       else
 329         {prinsign=0;}
 330     }
 331   else
 332     {
 333       Y_1 = (Cos_Theta*Sin_Theta*(A*A - B*B) *Dx +
 334              TMath::Sqrt(Delta)/2.)/(A*A*Cos_Theta*Cos_Theta +
 335                                      B*B*Sin_Theta*Sin_Theta) ;
 336       Y_2 = (Cos_Theta*Sin_Theta*(A*A - B*B) *Dx -
 337              TMath::Sqrt(Delta)/2.)/(A*A*Cos_Theta*Cos_Theta
 338                                       + B*B*Sin_Theta*Sin_Theta) ;
 339       Y_1 += Y_center ;
 340       Y_2 += Y_center ;
 341       if ((P[1]<=Y_1) && (P[1]>=Y_2))
 342         {prinsign=1;}
 343       else
 344         {prinsign=0;}
 345     }
 346   return prinsign;
 347 }
 348
 349 //____________________________________________________________________________
 350 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)
 351 {
 352
 353   // Set all ellipse parameters depending on the cluster energy and
 354   // Purity-Efficiency point (possible options "HIGH EFFICIENCY"
 355   // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing
 356   // EFFICIENCY by PURITY)
 357
 358   Int_t eff_pur = GetEffPurOption(Eff_Pur);
 359   GetAnalysisParameters(Cluster_En) ;
 360   if((fCluster!= -1)&&(eff_pur != -1)){
 361     (*fParameters)(fCluster+6 +fMatrixExtraRow,eff_pur) = x ;
 362     (*fParameters)(fCluster+9 +fMatrixExtraRow,eff_pur) = y ;
 363     (*fParameters)(fCluster+12+fMatrixExtraRow,eff_pur) = a ;
 364     (*fParameters)(fCluster+15+fMatrixExtraRow,eff_pur) = b ;
 365     (*fParameters)(fCluster+18+fMatrixExtraRow,eff_pur) = angle ;
 366   }
 367
 368 }
 369 //__________________________________________________________________________
 370 void  AliPHOSPIDv1::SetEllipseXCenter(Float_t Cluster_En, TString Eff_Pur, Float_t x)
 371 {
 372   // Set the ellipse parameter x_center depending on the custer energy and
 373   // Purity-Efficiency point (possible options "HIGH EFFICIENCY"
 374   // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing
 375   // EFFICIENCY by PURITY)
 376   Int_t eff_pur = GetEffPurOption(Eff_Pur);
 377   GetAnalysisParameters(Cluster_En) ;
 378   if((fCluster!= -1)&&(eff_pur != -1))
 379     (*fParameters)(fCluster+6+fMatrixExtraRow,eff_pur) = x ;
 380 }
 381 //_________________________________________________________________________
 382 void  AliPHOSPIDv1::SetEllipseYCenter(Float_t Cluster_En, TString Eff_Pur, Float_t y)
 383 {
 384
 385   // Set the ellipse parameter y_center depending on the cluster energy and
 386   // Purity-Efficiency point (possible options "HIGH EFFICIENCY"
 387   // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing
 388   // EFFICIENCY by PURITY)
 389
 390   Int_t eff_pur = GetEffPurOption(Eff_Pur);
 391   GetAnalysisParameters(Cluster_En) ;
 392   if((fCluster!= -1)&&(eff_pur != -1))
 393     (*fParameters)(fCluster+9+fMatrixExtraRow,eff_pur) = y ;
 394 }
 395 //_________________________________________________________________________
 396 void  AliPHOSPIDv1::SetEllipseAParameter(Float_t Cluster_En, TString Eff_Pur, Float_t a)
 397 {
 398   // Set the ellipse parameter a depending on the cluster energy and
 399   // Purity-Efficiency point (possible options "HIGH EFFICIENCY"
 400   // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing
 401   // EFFICIENCY by PURITY)
 402
 403   Int_t eff_pur = GetEffPurOption(Eff_Pur);
 404   GetAnalysisParameters(Cluster_En) ;
 405   if((fCluster!= -1)&&(eff_pur != -1))
 406     (*fParameters)(fCluster+12+fMatrixExtraRow,eff_pur) = a ;
 407 }
 408 //________________________________________________________________________
 409 void  AliPHOSPIDv1::SetEllipseBParameter(Float_t Cluster_En, TString Eff_Pur, Float_t b)
 410 {
 411   // Set the ellipse parameter b depending on the cluster energy and
 412   // Purity-Efficiency point (possible options "HIGH EFFICIENCY"
 413   // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing
 414   // EFFICIENCY by PURITY)
 415
 416   Int_t eff_pur = GetEffPurOption(Eff_Pur);
 417   GetAnalysisParameters(Cluster_En) ;
 418   if((fCluster!= -1)&&(eff_pur != -1))
 419     (*fParameters)(fCluster+15+fMatrixExtraRow,eff_pur) = b ;
 420 }
 421 //________________________________________________________________________
 422 void  AliPHOSPIDv1::SetEllipseAngle(Float_t Cluster_En, TString Eff_Pur, Float_t angle)
 423 {
 424
 425   // Set the ellipse parameter angle depending on the cluster energy and
 426   // Purity-Efficiency point (possible options "HIGH EFFICIENCY"
 427   // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing
 428   // EFFICIENCY by PURITY)
 429
 430   Int_t eff_pur = GetEffPurOption(Eff_Pur);
 431   GetAnalysisParameters(Cluster_En) ;
 432   if((fCluster!= -1)&&(eff_pur != -1))
 433     (*fParameters)(fCluster+18+fMatrixExtraRow,eff_pur) = angle ;
 434 }
 435 //_____________________________________________________________________________
 436 void  AliPHOSPIDv1::SetCpvtoEmcDistanceCut(Float_t Cluster_En, TString Eff_Pur, Float_t cut)
 437 {
 438
 439   // Set the parameter Cpvto EmcDistanceCut depending on the cluster energy and
 440   // Purity-Efficiency point (possible options "HIGH EFFICIENCY"
 441   // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing
 442   // EFFICIENCY by PURITY)
 443
 444
 445   Int_t eff_pur = GetEffPurOption(Eff_Pur);
 446   GetAnalysisParameters(Cluster_En) ;
 447   if((fClusterrcpv!= -1)&&(eff_pur != -1))
 448     (*fParameters)(fClusterrcpv,eff_pur) = cut ;
 449 }
 450 //_____________________________________________________________________________
 451 void  AliPHOSPIDv1::SetTimeGate(Float_t Cluster_En, TString Eff_Pur, Float_t gate)
 452 {
 453
 454   // Set the parameter TimeGate depending on the cluster energy and
 455   // Purity-Efficiency point (possible options "HIGH EFFICIENCY"
 456   // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing
 457   // EFFICIENCY by PURITY)
 458
 459   Int_t eff_pur = GetEffPurOption(Eff_Pur);
 460   GetAnalysisParameters(Cluster_En) ;
 461   if((fCluster!= -1)&&(eff_pur != -1))
 462     (*fParameters)(fCluster+3+fMatrixExtraRow,eff_pur) = gate ;
 463 }
 464 //_____________________________________________________________________________
 465 void  AliPHOSPIDv1::SetParameters()
 466                                   //TString OptFileName)
 467 {
 468   // PCA : To do the Principal Components Analysis it is necessary
 469   // the Principal file, which is opened here
 470   fX         = new double[7]; // Data for the PCA
 471   fP         = new double[7]; // Eigenvalues of the PCA
 472
 473
 474   // Set the principal and parameters files to be used
 475   fFileName5  = "$ALICE_ROOT/PHOS/PCA8pa15_0.5-5.root" ;
 476   fFileNamePar5 = gSystem->ExpandPathName("$ALICE_ROOT/PHOS/Parameters_0.5_5.dat");
 477   fFileName100  = "$ALICE_ROOT/PHOS/PCA8pa15_0.5-100.root" ;
 478   fFileNamePar100 = gSystem->ExpandPathName("$ALICE_ROOT/PHOS/Parameters_0.5_100.dat");
 479
 480   //SetPrincipalFileOptions();
 481   //fOptFileName);
 482   TFile f5( fFileName5.Data(), "read" ) ;
 483   fPrincipal5 = dynamic_cast<TPrincipal*> (f5.Get("principal")) ;
 484   f5.Close() ;
 485   TFile f100( fFileName100.Data(), "read" ) ;
 486   fPrincipal100 = dynamic_cast<TPrincipal*> (f100.Get("principal")) ;
 487   f100.Close() ;
 488   TFile f( fFileName100.Data(), "read" ) ;
 489   fPrincipal = dynamic_cast<TPrincipal*> (f.Get("principal")) ;
 490   f.Close() ;
 491   // Initialization of the Parameters matrix. In the File ParametersXX.dat
 492   // are all the parameters. These are introduced in a matrix of 21x3 or 22x3
 493   // elements (depending on the principal file 21 rows for 0.5-5 GeV and 22
 494   // rows for 5-100).
 495   // All the parameters defined in this file are, in order of row (there are
 496   // 3 rows per parameter): CpvtoEmcDistanceCut(if the principal file is 5-100
 497   // GeV then 4 rows), TimeGate and the ellipse parameters, X_center, Y_center,
 498   // a, b, angle. Each row of a given parameter depends on the cluster energy range
 499   // (wich depends on the chosen principal file)
 500   // Each column designs the parameters for a point in the Efficiency-Purity
 501   // of the photon identification P1(96%,63%), P2(87%,0.88%) and P3(68%,94%)
 502   // for the principal file from 0.5-5 GeV and for the other one P1(95%,79%),
 503   // P2(89%,90%) and P3(72%,96%)
 504
 505
 506   fParameters5 = new TMatrixD(21,3) ;
 507   fParameters100 = new TMatrixD(22,3) ;
 508   fParameters = new TMatrixD(22,3) ;
 509
 510   ifstream paramFile5(fFileNamePar5, ios::in) ;
 511
 512   Int_t i,j ;
 513
 514   for(i = 0; i< 21; i++){
 515     for(j = 0; j< 3; j++){
 516       paramFile5 >> (*fParameters5)(i,j) ;
 517     }
 518   }
 519   paramFile5.close();
 520
 521   ifstream paramFile100(fFileNamePar100, ios::in) ;
 522
 523   Int_t l,k ;
 524
 525   for(l = 0; l< 22; l++){
 526     for(k = 0; k< 3; k++){
 527       paramFile100 >> (*fParameters100)(l,k) ;
 528     }
 529   }
 530   paramFile100.close();
 531
 532   ifstream paramFile(fFileNamePar100, ios::in) ;
 533   Int_t h,n;
 534   for(h = 0; h< 22; h++){
 535     for(n = 0; n< 3; n++){
 536       paramFile >> (*fParameters)(h,n) ;
 537     }
 538   }
 539   paramFile.close();
 540
 541   fCluster = -1;
 542   fClusterrcpv = -1;
 543   fMatrixExtraRow = 0;
 544
 545   //Calibration parameters Encal = C * E^2 + B * E + A  (E is the energy from cluster)
 546   fACalParameter = 0.0241  ;
 547   fBCalParameter = 1.0504  ;
 548   fCCalParameter = 0.000249 ;
 549
 550   // fParameters->Print();
 551 }
 552 //_____________________________________________________________________________
 553 void  AliPHOSPIDv1::GetAnalysisParameters(Float_t Cluster_En)
 554 {
 555   if(Cluster_En <= 5.){
 556     fPrincipal  = fPrincipal5;
 557     fParameters = fParameters5;
 558     fMatrixExtraRow = 0;
 559     GetClusterOption(Cluster_En,kFALSE) ;
 560   }
 561   else{
 562     fPrincipal  = fPrincipal100;
 563     fParameters = fParameters100;
 564     fMatrixExtraRow = 1;
 565     GetClusterOption(Cluster_En,kTRUE) ;
 566   }
 567 }
 568
 569 //_____________________________________________________________________________
 570 void  AliPHOSPIDv1::GetClusterOption(const Float_t Cluster_En, const Bool_t range)
 571 {
 572
 573   // Gives the cluster energy range.
 574   // range = kFALSE Default analisys range from 0.5 to 5 GeV
 575   // range = kTRUE  analysis range from 0.5 to 100 GeV
 576
 577
 578   //Int_t cluster = -1 ;
 579
 580   if((range == kFALSE)){
 581     if((Cluster_En > 0.3)&&(Cluster_En <= 1.0)){
 582       fCluster = 0 ;
 583       fClusterrcpv = 0 ;
 584     }
 585     if((Cluster_En > 1.0)&&(Cluster_En <= 2.0)){
 586       fCluster = 1 ;
 587       fClusterrcpv = 1 ;
 588     }
 589     if( Cluster_En > 2.0){
 590       fCluster = 2 ;
 591       fClusterrcpv = 2 ;
 592     }
 593   }
 594   else if(range == kTRUE){
 595     if((Cluster_En > 0.5 )&&(Cluster_En <= 20.0)) fCluster = 0 ;
 596     if((Cluster_En > 20.0)&&(Cluster_En <= 50.0)) fCluster = 1 ;
 597     if( Cluster_En > 50.0)                        fCluster = 2 ;
 598     if((Cluster_En > 5.0 )&&(Cluster_En <= 10.0)) fClusterrcpv = 0 ;
 599     if((Cluster_En > 10.0)&&(Cluster_En <= 20.0)) fClusterrcpv = 1 ;
 600     if((Cluster_En > 20.0)&&(Cluster_En <= 30.0)) fClusterrcpv = 2 ;
 601     if( Cluster_En > 30.0)                        fClusterrcpv = 3 ;
 602   }
 603   else {
 604     fCluster = -1 ;
 605     fClusterrcpv = -1;
 606     cout<<"Invalid Energy option"<<endl;
 607   }
 608
 609   //return cluster;
 610 }
 611 //____________________________________________________________________________
 612 Int_t  AliPHOSPIDv1::GetEffPurOption(const TString Eff_Pur) const
 613 {
 614
 615   // Looks for the Purity-Efficiency point (possible options "HIGH EFFICIENCY"
 616   // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing
 617   // EFFICIENCY by PURITY)
 618
 619   Int_t eff_pur = -1 ;
 620
 621   if(Eff_Pur.Contains("HIGH EFFICIENCY") ||Eff_Pur.Contains("LOW PURITY") )
 622     eff_pur = 0 ;
 623   else if(Eff_Pur.Contains("MEDIUM EFFICIENCY") ||Eff_Pur.Contains("MEDIUM PURITY") )
 624     eff_pur = 1 ;
 625   else if(Eff_Pur.Contains("LOW EFFICIENCY")||Eff_Pur.Contains("HIGH PURITY") )
 626     eff_pur = 2 ;
 627   else{
 628     eff_pur = -1;
 629     cout<<"Invalid Efficiency-Purity option"<<endl;
 630     cout<<"Possible options: HIGH EFFICIENCY =    LOW PURITY"<<endl;
 631     cout<<"                MEDIUM EFFICIENCY = MEDIUM PURITY"<<endl;
 632     cout<<"                   LOW EFFICIENCY =   HIGH PURITY"<<endl;
 633   }
 634
 635   return eff_pur;
 636 }
 637 //____________________________________________________________________________
 638
 639 void  AliPHOSPIDv1::Exec(Option_t * option)
 640 {
 641   //Steering method
 642
 643   if( strcmp(GetName(), "")== 0 )
 644     Init() ;
 645
 646   if(strstr(option,"tim"))
 647     gBenchmark->Start("PHOSPID");
 648
 649   if(strstr(option,"print")) {
 650     Print("") ;
 651     return ;
 652   }
 653
 654   //cout << gDirectory->GetName() << endl ;
 655
 656   gAlice->GetEvent(0) ;
 657
 658   //check, if the branch with name of this" already exits?
 659   if (gAlice->TreeR()) {
 660     TObjArray * lob = (TObjArray*)gAlice->TreeR()->GetListOfBranches() ;
 661     TIter next(lob) ;
 662     TBranch * branch = 0 ;
 663     Bool_t phospidfound = kFALSE, pidfound = kFALSE ;
 664
 665     TString taskName(GetName()) ;
 666     taskName.Remove(taskName.Index(Version())-1) ;
 667
 668     while ( (branch = (TBranch*)next()) && (!phospidfound || !pidfound) ) {
 669       if ( (strcmp(branch->GetName(), "PHOSPID")==0) && (strcmp(branch->GetTitle(), taskName.Data())==0) )
 670         phospidfound = kTRUE ;
 671
 672       else if ( (strcmp(branch->GetName(), "AliPHOSPID")==0) && (strcmp(branch->GetTitle(), taskName.Data())==0) )
 673         pidfound = kTRUE ;
 674     }
 675
 676     if ( phospidfound || pidfound ) {
 677       cerr << "WARNING: AliPHOSPIDv1::Exec -> RecParticles and/or PIDtMaker branch with name "
 678            << taskName.Data() << " already exits" << endl ;
 679       return ;
 680     }
 681   }
 682
 683   Int_t nevents = (Int_t) gAlice->TreeE()->GetEntries() ;
 684   Int_t ievent ;
 685   AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
 686   for(ievent = 0; ievent < nevents; ievent++){
 687     gime->Event(ievent,"R") ;
 688
 689     MakeRecParticles() ;
 690
 691     WriteRecParticles(ievent);
 692
 693     if(strstr(option,"deb"))
 694       PrintRecParticles(option) ;
 695
 696     //increment the total number of rec particles per run
 697     fRecParticlesInRun += gime->RecParticles(BranchName())->GetEntriesFast() ;
 698
 699   }
 700
 701   if(strstr(option,"tim")){
 702     gBenchmark->Stop("PHOSPID");
 703     cout << "AliPHOSPID:" << endl ;
 704     cout << "  took " << gBenchmark->GetCpuTime("PHOSPID") << " seconds for PID "
 705          <<  gBenchmark->GetCpuTime("PHOSPID")/nevents << " seconds per event " << endl ;
 706     cout << endl ;
 707   }
 708
 709 }
 710
 711 //____________________________________________________________________________
 712 void  AliPHOSPIDv1::MakeRecParticles(){
 713
 714   // Makes a RecParticle out of a TrackSegment
 715
 716   AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
 717   TObjArray * emcRecPoints = gime->EmcRecPoints(fFrom) ;
 718   TObjArray * cpvRecPoints = gime->CpvRecPoints(fFrom) ;
 719   TClonesArray * trackSegments = gime->TrackSegments(fFrom) ;
 720   if ( !emcRecPoints || !cpvRecPoints || !trackSegments ) {
 721     cerr << "ERROR:  AliPHOSPIDv1::MakeRecParticles -> RecPoints or TrackSegments with name "
 722          << fFrom << " not found ! " << endl ;
 723     abort() ;
 724   }
 725   TClonesArray * recParticles  = gime->RecParticles(BranchName()) ;
 726   recParticles->Clear();
 727
 728
 729   TIter next(trackSegments) ;
 730   AliPHOSTrackSegment * ts ;
 731   Int_t index = 0 ;
 732   AliPHOSRecParticle * rp ;
 733
 734   while ( (ts = (AliPHOSTrackSegment *)next()) ) {
 735
 736     new( (*recParticles)[index] ) AliPHOSRecParticle() ;
 737     rp = (AliPHOSRecParticle *)recParticles->At(index) ;
 738     rp->SetTrackSegment(index) ;
 739     rp->SetIndexInList(index) ;
 740
 741     AliPHOSEmcRecPoint * emc = 0 ;
 742     if(ts->GetEmcIndex()>=0)
 743       emc = (AliPHOSEmcRecPoint *) emcRecPoints->At(ts->GetEmcIndex()) ;
 744
 745     AliPHOSRecPoint    * cpv = 0 ;
 746     if(ts->GetCpvIndex()>=0)
 747       cpv = (AliPHOSRecPoint *)   cpvRecPoints->At(ts->GetCpvIndex()) ;
 748
 749     // Now set type (reconstructed) of the particle
 750
 751     // Choose the cluster energy range
 752
 753     Float_t    e = emc->GetEnergy() ;
 754
 755
 756     GetAnalysisParameters(e);
 757
 758     if((fCluster== -1)||(fClusterrcpv == -1)) continue ;
 759
 760     // Ellipse and rcpv cut in function of the cluster energy
 761
 762     // Loop of Efficiency-Purity (the 3 points of purity or efficiency are taken
 763     // into account to set the particle identification)
 764     for(Int_t eff_pur = 0; eff_pur < 3 ; eff_pur++){
 765
 766       Float_t  lambda[2] ;
 767       emc->GetElipsAxis(lambda) ;
 768       Float_t time =emc->GetTime() ;
 769
 770       if((lambda[0]>0.01) && (lambda[1]>0.01) && time > 0.){
 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