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