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