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 /* History of cvs commits:
  19  *
  20  * $Log$
  21  * Revision 1.106  2007/04/01 15:40:15  kharlov
  22  * Correction for actual vertex position implemented
  23  *
  24  * Revision 1.105  2007/03/06 06:57:46  kharlov
  25  * DP:calculation of distance to CPV done in TSM
  26  *
  27  * Revision 1.104  2006/12/15 10:46:26  hristov
  28  * Using TMath::Abs instead of fabs
  29  *
  30  * Revision 1.103  2006/09/07 18:31:08  kharlov
  31  * Effective c++ corrections (T.Pocheptsov)
  32  *
  33  * Revision 1.102  2006/01/23 17:51:48  hristov
  34  * Using the recommended way of forward declarations for TVector and TMatrix (see v5-08-00 release notes). Additional clean-up
  35  *
  36  * Revision 1.101  2005/05/28 14:19:04  schutz
  37  * Compilation warnings fixed by T.P.
  38  *
  39  */
  40
  41 //_________________________________________________________________________
  42 // Implementation version v1 of the PHOS particle identifier
  43 // Particle identification based on the
  44 //     - RCPV: distance from CPV recpoint to EMCA recpoint.
  45 //     - TOF
  46 //     - PCA: Principal Components Analysis..
  47 // The identified particle has an identification number corresponding
  48 // to a 9 bits number:
  49 //     -Bit 0 to 2: bit set if RCPV > CpvEmcDistance (each bit corresponds
  50 //      to a different efficiency-purity point of the photon identification)
  51 //     -Bit 3 to 5: bit set if TOF  < TimeGate (each bit corresponds
  52 //      to a different efficiency-purity point of the photon identification)
  53 //     -Bit 6 to 9: bit set if Principal Components are
  54 //      inside an ellipse defined by the parameters a, b, c, x0 and y0.
  55 //      (each bit corresponds to a different efficiency-purity point of the
  56 //      photon identification)
  57 //      The PCA (Principal components analysis) needs a file that contains
  58 //      a previous analysis of the correlations between the particles. This
  59 //      file is $ALICE_ROOT/PHOS/PCA8pa15_0.5-100.root. Analysis done for
  60 //      energies between 0.5 and 100 GeV.
  61 //      A calibrated energy is calculated. The energy of the reconstructed
  62 //      cluster is corrected with the formula A + B * E  + C * E^2, whose
  63 //      parameters where obtained through the study of the reconstructed
  64 //      energy distribution of monoenergetic photons.
  65 //
  66 //      All the parameters (RCPV(2 rows-3 columns),TOF(1r-3c),PCA(5r-4c)
  67 //      and calibration(1r-3c))are stored in a file called
  68 //      $ALICE_ROOT/PHOS/Parameters.dat. Each time that AliPHOSPIDv1 is
  69 //      initialized, this parameters are copied to a Matrix (9,4), a
  70 //      TMatrixD object.
  71 //
  72 // use case:
  73 //  root [0] AliPHOSPIDv1 * p = new AliPHOSPIDv1("galice1.root")
  74 //  Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
  75 //          // reading headers from file galice1.root and create  RecParticles
  76             // TrackSegments and RecPoints are used
  77 //          // set file name for the branch RecParticles
  78 //  root [1] p->ExecuteTask("deb all time")
  79 //          // available options
  80 //          // "deb" - prints # of reconstructed particles
  81 //          // "deb all" -  prints # and list of RecParticles
  82 //          // "time" - prints benchmarking results
  83 //
  84 //  root [2] AliPHOSPIDv1 * p2 = new AliPHOSPIDv1("galice1.root","v1",kTRUE)
  85 //  Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
  86 //                //Split mode.
  87 //  root [3] p2->ExecuteTask()
  88 //
  89
  90
  91 //*-- Author: Yves Schutz (SUBATECH)  & Gines Martinez (SUBATECH) &
  92 //            Gustavo Conesa April 2002
  93 //            PCA redesigned by Gustavo Conesa October 2002:
  94 //            The way of using the PCA has changed. Instead of 2
  95 //            files with the PCA, each one with different energy ranges
  96 //            of application, we use the wide one (0.5-100 GeV), and instead
  97 //            of fixing 3 ellipses for different ranges of energy, it has been
  98 //            studied the dependency of the ellipses parameters with the
  99 //            energy, and they are implemented in the code as a funtion
 100 //            of the energy.
 101 //
 102 //
 103 //
 104 // --- ROOT system ---
 105
 106
 107 // --- Standard library ---
 108 #include <TMatrixF.h>
 109 #include "TFormula.h"
 110 #include "TBenchmark.h"
 111 #include "TPrincipal.h"
 112 #include "TFile.h"
 113 #include "TSystem.h"
 114 #include "TVector3.h"
 115
 116 // --- AliRoot header files ---
 117               //#include "AliLog.h"
 118 #include "AliPHOS.h"
 119 #include "AliPHOSPIDv1.h"
 120 #include "AliPHOSGetter.h"
 121 #include "AliESD.h"
 122 #include "AliESDVertex.h"
 123
 124 ClassImp( AliPHOSPIDv1)
 125
 126 //____________________________________________________________________________
 127 AliPHOSPIDv1::AliPHOSPIDv1() :
 128   fBayesian(kFALSE),
 129   fDefaultInit(kFALSE),
 130   fWrite(kFALSE),
 131   fNEvent(0),
 132   fFileNamePrincipalPhoton(),
 133   fFileNamePrincipalPi0(),
 134   fFileNameParameters(),
 135   fPrincipalPhoton(0),
 136   fPrincipalPi0(0),
 137   fX(0),
 138   fPPhoton(0),
 139   fPPi0(0),
 140   fRecParticlesInRun(0),
 141   fParameters(0),
 142   fTFphoton(0),
 143   fTFpiong(0),
 144   fTFkaong(0),
 145   fTFkaonl(0),
 146   fTFhhadrong(0),
 147   fTFhhadronl(0),
 148   fDFmuon(0),
 149   fERecWeight(0),
 150   fChargedNeutralThreshold(0.),
 151   fTOFEnThreshold(0),
 152   fDispEnThreshold(0),
 153   fDispMultThreshold(0)
 154 {
 155   // default ctor
 156
 157   InitParameters() ;
 158   fDefaultInit = kTRUE ;
 159 }
 160
 161 //____________________________________________________________________________
 162 AliPHOSPIDv1::AliPHOSPIDv1(const AliPHOSPIDv1 & pid ) :
 163   AliPHOSPID(pid),
 164   fBayesian(kFALSE),
 165   fDefaultInit(kFALSE),
 166   fWrite(kFALSE),
 167   fNEvent(0),
 168   fFileNamePrincipalPhoton(),
 169   fFileNamePrincipalPi0(),
 170   fFileNameParameters(),
 171   fPrincipalPhoton(0),
 172   fPrincipalPi0(0),
 173   fX(0),
 174   fPPhoton(0),
 175   fPPi0(0),
 176   fRecParticlesInRun(0),
 177   fParameters(0),
 178   fTFphoton(0),
 179   fTFpiong(0),
 180   fTFkaong(0),
 181   fTFkaonl(0),
 182   fTFhhadrong(0),
 183   fTFhhadronl(0),
 184   fDFmuon(0),
 185   fERecWeight(0),
 186   fChargedNeutralThreshold(0.),
 187   fTOFEnThreshold(0),
 188   fDispEnThreshold(0),
 189   fDispMultThreshold(0)
 190
 191 {
 192   // ctor
 193   InitParameters() ;
 194   Init() ;
 195
 196 }
 197
 198 //____________________________________________________________________________
 199 AliPHOSPIDv1::AliPHOSPIDv1(const TString alirunFileName, const TString eventFolderName) :
 200   AliPHOSPID(alirunFileName, eventFolderName),
 201   fBayesian(kFALSE),
 202   fDefaultInit(kFALSE),
 203   fWrite(kFALSE),
 204   fNEvent(0),
 205   fFileNamePrincipalPhoton(),
 206   fFileNamePrincipalPi0(),
 207   fFileNameParameters(),
 208   fPrincipalPhoton(0),
 209   fPrincipalPi0(0),
 210   fX(0),
 211   fPPhoton(0),
 212   fPPi0(0),
 213   fRecParticlesInRun(0),
 214   fParameters(0),
 215   fTFphoton(0),
 216   fTFpiong(0),
 217   fTFkaong(0),
 218   fTFkaonl(0),
 219   fTFhhadrong(0),
 220   fTFhhadronl(0),
 221   fDFmuon(0),
 222   fERecWeight(0),
 223   fChargedNeutralThreshold(0.),
 224   fTOFEnThreshold(0),
 225   fDispEnThreshold(0),
 226   fDispMultThreshold(0)
 227
 228 {
 229   //ctor with the indication on where to look for the track segments
 230
 231   InitParameters() ;
 232   Init() ;
 233   fDefaultInit = kFALSE ;
 234 }
 235
 236 //____________________________________________________________________________
 237 AliPHOSPIDv1::~AliPHOSPIDv1()
 238 {
 239   // dtor
 240   fPrincipalPhoton = 0;
 241   fPrincipalPi0 = 0;
 242
 243   delete [] fX ;       // Principal input
 244   delete [] fPPhoton ; // Photon Principal components
 245   delete [] fPPi0 ;    // Pi0 Principal components
 246
 247   delete fParameters;
 248   delete fTFphoton;
 249   delete fTFpiong;
 250   delete fTFkaong;
 251   delete fTFkaonl;
 252   delete fTFhhadrong;
 253   delete fTFhhadronl;
 254   delete fDFmuon;
 255 }
 256 //____________________________________________________________________________
 257 const TString AliPHOSPIDv1::BranchName() const
 258 {
 259
 260   return GetName() ;
 261 }
 262
 263 //____________________________________________________________________________
 264 void AliPHOSPIDv1::Init()
 265 {
 266   // Make all memory allocations that are not possible in default constructor
 267   // Add the PID task to the list of PHOS tasks
 268
 269   AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
 270   if(!gime)
 271     gime = AliPHOSGetter::Instance(GetTitle(), fEventFolderName.Data()) ;
 272
 273   if ( !gime->PID() )
 274     gime->PostPID(this) ;
 275 }
 276
 277 //____________________________________________________________________________
 278 void AliPHOSPIDv1::InitParameters()
 279 {
 280   // Initialize PID parameters
 281   fWrite                   = kTRUE ;
 282   fRecParticlesInRun = 0 ;
 283   fNEvent            = 0 ;
 284   fRecParticlesInRun = 0 ;
 285   fBayesian          = kTRUE ;
 286   SetParameters() ; // fill the parameters matrix from parameters file
 287   SetEventRange(0,-1) ;
 288
 289   // initialisation of response function parameters
 290   // Tof
 291
 292 //   // Photons
 293 //   fTphoton[0] = 0.218    ;
 294 //   fTphoton[1] = 1.55E-8  ;
 295 //   fTphoton[2] = 5.05E-10 ;
 296 //   fTFphoton = new TFormula("ToF response to photons" , "gaus") ;
 297 //   fTFphoton->SetParameters( fTphoton[0], fTphoton[1], fTphoton[2]) ;
 298
 299 //   // Pions
 300 //   //Gaus (0 to max probability)
 301 //   fTpiong[0] = 0.0971    ;
 302 //   fTpiong[1] = 1.58E-8  ;
 303 //   fTpiong[2] = 5.69E-10 ;
 304 //   fTFpiong = new TFormula("ToF response to pions" , "gaus") ;
 305 //   fTFpiong->SetParameters( fTpiong[0], fTpiong[1], fTpiong[2]) ;
 306
 307 //   // Kaons
 308 //   //Gaus (0 to max probability)
 309 //   fTkaong[0] = 0.0542  ;
 310 //   fTkaong[1] = 1.64E-8 ;
 311 //   fTkaong[2] = 6.07E-10 ;
 312 //   fTFkaong = new TFormula("ToF response to kaon" , "gaus") ;
 313 //   fTFkaong->SetParameters( fTkaong[0], fTkaong[1], fTkaong[2]) ;
 314 //   //Landau (max probability to inf)
 315 //   fTkaonl[0] = 0.264   ;
 316 //   fTkaonl[1] = 1.68E-8  ;
 317 //   fTkaonl[2] = 4.10E-10 ;
 318 //   fTFkaonl = new TFormula("ToF response to kaon" , "landau") ;
 319 //   fTFkaonl->SetParameters( fTkaonl[0], fTkaonl[1], fTkaonl[2]) ;
 320
 321 //   //Heavy Hadrons
 322 //   //Gaus (0 to max probability)
 323 //   fThhadrong[0] = 0.0302   ;
 324 //   fThhadrong[1] = 1.73E-8  ;
 325 //   fThhadrong[2] = 9.52E-10 ;
 326 //   fTFhhadrong = new TFormula("ToF response to heavy hadrons" , "gaus") ;
 327 //   fTFhhadrong->SetParameters( fThhadrong[0], fThhadrong[1], fThhadrong[2]) ;
 328 //   //Landau (max probability to inf)
 329 //   fThhadronl[0] = 0.139    ;
 330 //   fThhadronl[1] = 1.745E-8  ;
 331 //   fThhadronl[2] = 1.00E-9  ;
 332 //   fTFhhadronl = new TFormula("ToF response to heavy hadrons" , "landau") ;
 333 //   fTFhhadronl->SetParameters( fThhadronl[0], fThhadronl[1], fThhadronl[2]) ;
 334
 335   // Photons
 336   fTphoton[0] = 7.83E8   ;
 337   fTphoton[1] = 1.55E-8  ;
 338   fTphoton[2] = 5.09E-10 ;
 339   fTFphoton = new TFormula("ToF response to photons" , "gaus") ;
 340   fTFphoton->SetParameters( fTphoton[0], fTphoton[1], fTphoton[2]) ;
 341
 342   // Pions
 343   //Gaus (0 to max probability)
 344   fTpiong[0] = 6.73E8    ;
 345   fTpiong[1] = 1.58E-8  ;
 346   fTpiong[2] = 5.87E-10 ;
 347   fTFpiong = new TFormula("ToF response to pions" , "gaus") ;
 348   fTFpiong->SetParameters( fTpiong[0], fTpiong[1], fTpiong[2]) ;
 349
 350   // Kaons
 351   //Gaus (0 to max probability)
 352   fTkaong[0] = 3.93E8  ;
 353   fTkaong[1] = 1.64E-8 ;
 354   fTkaong[2] = 6.07E-10 ;
 355   fTFkaong = new TFormula("ToF response to kaon" , "gaus") ;
 356   fTFkaong->SetParameters( fTkaong[0], fTkaong[1], fTkaong[2]) ;
 357   //Landau (max probability to inf)
 358   fTkaonl[0] = 2.0E9    ;
 359   fTkaonl[1] = 1.68E-8  ;
 360   fTkaonl[2] = 4.10E-10 ;
 361   fTFkaonl = new TFormula("ToF response to kaon" , "landau") ;
 362   fTFkaonl->SetParameters( fTkaonl[0], fTkaonl[1], fTkaonl[2]) ;
 363
 364   //Heavy Hadrons
 365   //Gaus (0 to max probability)
 366   fThhadrong[0] = 2.02E8   ;
 367   fThhadrong[1] = 1.73E-8  ;
 368   fThhadrong[2] = 9.52E-10 ;
 369   fTFhhadrong = new TFormula("ToF response to heavy hadrons" , "gaus") ;
 370   fTFhhadrong->SetParameters( fThhadrong[0], fThhadrong[1], fThhadrong[2]) ;
 371   //Landau (max probability to inf)
 372   fThhadronl[0] = 1.10E9    ;
 373   fThhadronl[1] = 1.74E-8   ;
 374   fThhadronl[2] = 1.00E-9   ;
 375   fTFhhadronl = new TFormula("ToF response to heavy hadrons" , "landau") ;
 376   fTFhhadronl->SetParameters( fThhadronl[0], fThhadronl[1], fThhadronl[2]) ;
 377
 378
 379
 380   // Shower shape: dispersion gaussian parameters
 381   // Photons
 382
 383 //   fDphoton[0] = 4.62e-2;  fDphoton[1] = 1.39e-2 ; fDphoton[2] = -3.80e-2;//constant
 384 //   fDphoton[3] = 1.53   ;  fDphoton[4] =-6.62e-2 ; fDphoton[5] = 0.339   ;//mean
 385 //   fDphoton[6] = 6.89e-2;  fDphoton[7] =-6.59e-2 ; fDphoton[8] = 0.194   ;//sigma
 386
 387 //   fDpi0[0] = 0.0586  ;  fDpi0[1] = 1.06E-3 ; fDpi0[2] = 0.      ;//constant
 388 //   fDpi0[3] = 2.67    ;  fDpi0[4] =-2.00E-2 ; fDpi0[5] = 9.37E-5 ;//mean
 389 //   fDpi0[6] = 0.153   ;  fDpi0[7] = 9.34E-4 ; fDpi0[8] =-1.49E-5 ;//sigma
 390
 391 //   fDhadron[0] = 1.61E-2 ;  fDhadron[1] = 3.03E-3 ; fDhadron[2] = 1.01E-2 ;//constant
 392 //   fDhadron[3] = 3.81    ;  fDhadron[4] = 0.232   ; fDhadron[5] =-1.25    ;//mean
 393 //   fDhadron[6] = 0.897   ;  fDhadron[7] = 0.0987  ; fDhadron[8] =-0.534   ;//sigma
 394
 395   fDphoton[0] = 1.5    ;  fDphoton[1] = 0.49    ; fDphoton[2] =-1.7E-2 ;//constant
 396   fDphoton[3] = 1.5    ;  fDphoton[4] = 4.0E-2  ; fDphoton[5] = 0.21   ;//mean
 397   fDphoton[6] = 4.8E-2 ;  fDphoton[7] =-0.12    ; fDphoton[8] = 0.27   ;//sigma
 398   fDphoton[9] = 16.; //for E>  fDphoton[9] parameters calculated at  fDphoton[9]
 399
 400   fDpi0[0] = 0.25      ;  fDpi0[1] = 3.3E-2     ; fDpi0[2] =-1.0e-5    ;//constant
 401   fDpi0[3] = 1.50      ;  fDpi0[4] = 398.       ; fDpi0[5] = 12.       ;//mean
 402   fDpi0[6] =-7.0E-2    ;  fDpi0[7] =-524.       ; fDpi0[8] = 22.       ;//sigma
 403   fDpi0[9] = 110.; //for E>  fDpi0[9] parameters calculated at  fDpi0[9]
 404
 405   fDhadron[0] = 6.5    ;  fDhadron[1] =-5.3     ; fDhadron[2] = 1.5    ;//constant
 406   fDhadron[3] = 3.8    ;  fDhadron[4] = 0.23    ; fDhadron[5] =-1.2    ;//mean
 407   fDhadron[6] = 0.88   ;  fDhadron[7] = 9.3E-2  ; fDhadron[8] =-0.51   ;//sigma
 408   fDhadron[9] = 2.; //for E>  fDhadron[9] parameters calculated at  fDhadron[9]
 409
 410   fDmuon[0] = 0.0631 ;
 411   fDmuon[1] = 1.4    ;
 412   fDmuon[2] = 0.0557 ;
 413   fDFmuon = new TFormula("Shower shape response to muons" , "landau") ;
 414   fDFmuon->SetParameters( fDmuon[0], fDmuon[1], fDmuon[2]) ;
 415
 416
 417   // x(CPV-EMC) distance gaussian parameters
 418
 419 //   fXelectron[0] = 8.06e-2 ;  fXelectron[1] = 1.00e-2; fXelectron[2] =-5.14e-2;//constant
 420 //   fXelectron[3] = 0.202   ;  fXelectron[4] = 8.15e-3; fXelectron[5] = 4.55   ;//mean
 421 //   fXelectron[6] = 0.334   ;  fXelectron[7] = 0.186  ; fXelectron[8] = 4.32e-2;//sigma
 422
 423 //   //charged hadrons gaus
 424 //   fXcharged[0] = 6.43e-3 ;  fXcharged[1] =-4.19e-5; fXcharged[2] = 1.42e-3;//constant
 425 //   fXcharged[3] = 2.75    ;  fXcharged[4] =-0.40   ; fXcharged[5] = 1.68   ;//mean
 426 //   fXcharged[6] = 3.135   ;  fXcharged[7] =-9.41e-2; fXcharged[8] = 1.31e-2;//sigma
 427
 428 //   // z(CPV-EMC) distance gaussian parameters
 429
 430 //   fZelectron[0] = 8.22e-2 ;  fZelectron[1] = 5.11e-3; fZelectron[2] =-3.05e-2;//constant
 431 //   fZelectron[3] = 3.09e-2 ;  fZelectron[4] = 5.87e-2; fZelectron[5] =-9.49e-2;//mean
 432 //   fZelectron[6] = 0.263   ;  fZelectron[7] =-9.02e-3; fZelectron[8] = 0.151 ;//sigma
 433
 434 //   //charged hadrons gaus
 435
 436 //   fZcharged[0] = 1.00e-2 ;  fZcharged[1] = 2.82E-4 ; fZcharged[2] = 2.87E-3 ;//constant
 437 //   fZcharged[3] =-4.68e-2 ;  fZcharged[4] =-9.21e-3 ; fZcharged[5] = 4.91e-2 ;//mean
 438 //   fZcharged[6] = 1.425   ;  fZcharged[7] =-5.90e-2 ; fZcharged[8] = 5.07e-2 ;//sigma
 439
 440
 441   fXelectron[0] =-1.6E-2 ;  fXelectron[1] = 0.77  ; fXelectron[2] =-0.15 ;//constant
 442   fXelectron[3] = 0.35   ;  fXelectron[4] = 0.25  ; fXelectron[5] = 4.12 ;//mean
 443   fXelectron[6] = 0.30   ;  fXelectron[7] = 0.11  ; fXelectron[8] = 0.16 ;//sigma
 444   fXelectron[9] = 3.; //for E>  fXelectron[9] parameters calculated at  fXelectron[9]
 445
 446   //charged hadrons gaus
 447   fXcharged[0] = 0.14    ;  fXcharged[1] =-3.0E-2 ; fXcharged[2] = 0     ;//constant
 448   fXcharged[3] = 1.4     ;  fXcharged[4] =-9.3E-2 ; fXcharged[5] = 1.4   ;//mean
 449   fXcharged[6] = 5.7     ;  fXcharged[7] = 0.27   ; fXcharged[8] =-1.8   ;//sigma
 450   fXcharged[9] = 1.2; //for E>  fXcharged[9] parameters calculated at  fXcharged[9]
 451
 452   // z(CPV-EMC) distance gaussian parameters
 453
 454   fZelectron[0] = 0.49   ;  fZelectron[1] = 0.53   ; fZelectron[2] =-9.8E-2 ;//constant
 455   fZelectron[3] = 2.8E-2 ;  fZelectron[4] = 5.0E-2 ; fZelectron[5] =-8.2E-2 ;//mean
 456   fZelectron[6] = 0.25   ;  fZelectron[7] =-1.7E-2 ; fZelectron[8] = 0.17   ;//sigma
 457   fZelectron[9] = 3.; //for E>  fZelectron[9] parameters calculated at  fZelectron[9]
 458
 459   //charged hadrons gaus
 460
 461   fZcharged[0] = 0.46    ;  fZcharged[1] =-0.65    ; fZcharged[2] = 0.52    ;//constant
 462   fZcharged[3] = 1.1E-2  ;  fZcharged[4] = 0.      ; fZcharged[5] = 0.      ;//mean
 463   fZcharged[6] = 0.60    ;  fZcharged[7] =-8.2E-2  ; fZcharged[8] = 0.45    ;//sigma
 464   fZcharged[9] = 1.2; //for E>  fXcharged[9] parameters calculated at  fXcharged[9]
 465
 466   //Threshold to differentiate between charged and neutral
 467   fChargedNeutralThreshold = 1e-5;
 468   fTOFEnThreshold          = 2;          //Maximum energy to use TOF
 469   fDispEnThreshold         = 0.5;       //Minimum energy to use shower shape
 470   fDispMultThreshold       = 3;       //Minimum multiplicity to use shower shape
 471
 472   //Weight to hadrons recontructed energy
 473
 474   fERecWeightPar[0] = 0.32 ;
 475   fERecWeightPar[1] = 3.8  ;
 476   fERecWeightPar[2] = 5.4E-3 ;
 477   fERecWeightPar[3] = 5.6E-2 ;
 478   fERecWeight = new TFormula("Weight for hadrons" , "[0]*exp(-x*[1])+[2]*exp(-x*[3])") ;
 479   fERecWeight ->SetParameters(fERecWeightPar[0],fERecWeightPar[1] ,fERecWeightPar[2] ,fERecWeightPar[3]) ;
 480
 481
 482   for (Int_t i =0; i<  AliPID::kSPECIESN ; i++)
 483     fInitPID[i] = 1.;
 484
 485 }
 486
 487 //________________________________________________________________________
 488 void  AliPHOSPIDv1::Exec(Option_t *option)
 489 {
 490   // Steering method to perform particle reconstruction and identification
 491   // for the event range from fFirstEvent to fLastEvent.
 492   // This range is optionally set by SetEventRange().
 493   // if fLastEvent=-1 (by default), then process events until the end.
 494
 495   if(strstr(option,"tim"))
 496     gBenchmark->Start("PHOSPID");
 497
 498   if(strstr(option,"print")) {
 499     Print() ;
 500     return ;
 501   }
 502
 503
 504   AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
 505
 506   if (fLastEvent == -1)
 507     fLastEvent = gime->MaxEvent() - 1 ;
 508   else
 509     fLastEvent = TMath::Min(fLastEvent,gime->MaxEvent());
 510   Int_t nEvents   = fLastEvent - fFirstEvent + 1;
 511
 512   Int_t ievent ;
 513   for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
 514     gime->Event(ievent,"TR") ;
 515     if(gime->TrackSegments() && //Skip events, where no track segments made
 516        gime->TrackSegments()->GetEntriesFast()) {
 517
 518       GetVertex() ;
 519       MakeRecParticles() ;
 520
 521       if(fBayesian)
 522         MakePID() ;
 523
 524       WriteRecParticles();
 525       if(strstr(option,"deb"))
 526         PrintRecParticles(option) ;
 527       //increment the total number of rec particles per run
 528       fRecParticlesInRun += gime->RecParticles()->GetEntriesFast() ;
 529     }
 530   }
 531   if(strstr(option,"deb"))
 532       PrintRecParticles(option);
 533   if(strstr(option,"tim")){
 534     gBenchmark->Stop("PHOSPID");
 535     AliInfo(Form("took %f seconds for PID %f seconds per event",
 536          gBenchmark->GetCpuTime("PHOSPID"),
 537          gBenchmark->GetCpuTime("PHOSPID")/nEvents)) ;
 538   }
 539   if(fWrite)
 540     Unload();
 541 }
 542
 543 //________________________________________________________________________
 544 Double_t  AliPHOSPIDv1::GausF(Double_t  x, Double_t  y, Double_t * par)
 545 {
 546   //Given the energy x and the parameter y (tof, shower dispersion or cpv-emc distance),
 547   //this method returns a density probability of this parameter, given by a gaussian
 548   //function whose parameters depend with the energy  with a function: a/(x*x)+b/x+b
 549   //Float_t xorg = x;
 550   if (x > par[9]) x = par[9];
 551
 552   //Double_t cnt    = par[1] / (x*x) + par[2] / x + par[0] ;
 553   Double_t cnt    = par[0] + par[1] * x + par[2] * x * x ;
 554   Double_t mean   = par[4] / (x*x) + par[5] / x + par[3] ;
 555   Double_t sigma  = par[7] / (x*x) + par[8] / x + par[6] ;
 556
 557 //   if(xorg > 30)
 558 //     cout<<"En_in = "<<xorg<<"; En_out = "<<x<<"; cnt = "<<cnt
 559 //      <<"; mean = "<<mean<<"; sigma = "<<sigma<<endl;
 560
 561   //  Double_t arg    = - (y-mean) * (y-mean) / (2*sigma*sigma) ;
 562   //  return cnt * TMath::Exp(arg) ;
 563   if(TMath::Abs(sigma) > 1.e-10){
 564     return cnt*TMath::Gaus(y,mean,sigma);
 565   }
 566   else
 567     return 0.;
 568
 569 }
 570 //________________________________________________________________________
 571 Double_t  AliPHOSPIDv1::GausPol2(Double_t  x, Double_t y, Double_t * par)
 572 {
 573   //Given the energy x and the parameter y (tof, shower dispersion or cpv-emc distance),
 574   //this method returns a density probability of this parameter, given by a gaussian
 575   //function whose parameters depend with the energy like second order polinomial
 576
 577   Double_t cnt    = par[0] + par[1] * x + par[2] * x * x ;
 578   Double_t mean   = par[3] + par[4] * x + par[5] * x * x ;
 579   Double_t sigma  = par[6] + par[7] * x + par[8] * x * x ;
 580
 581   if(TMath::Abs(sigma) > 1.e-10){
 582     return cnt*TMath::Gaus(y,mean,sigma);
 583   }
 584   else
 585     return 0.;
 586
 587
 588
 589 }
 590
 591 //____________________________________________________________________________
 592 const TString AliPHOSPIDv1::GetFileNamePrincipal(TString particle) const
 593 {
 594   //Get file name that contains the PCA for a particle ("photon or pi0")
 595   particle.ToLower();
 596   TString name;
 597   if      (particle=="photon")
 598     name = fFileNamePrincipalPhoton ;
 599   else if (particle=="pi0"   )
 600     name = fFileNamePrincipalPi0    ;
 601   else
 602     AliError(Form("Wrong particle name: %s (choose from pi0/photon)\n",
 603                   particle.Data()));
 604   return name;
 605 }
 606
 607 //____________________________________________________________________________
 608 Float_t  AliPHOSPIDv1::GetParameterCalibration(Int_t i) const
 609 {
 610   // Get the i-th parameter "Calibration"
 611   Float_t param = 0.;
 612   if (i>2 || i<0) {
 613     AliError(Form("Invalid parameter number: %d",i));
 614   } else
 615     param = (*fParameters)(0,i);
 616   return param;
 617 }
 618
 619 //____________________________________________________________________________
 620 Float_t  AliPHOSPIDv1::GetCalibratedEnergy(Float_t e) const
 621 {
 622 //      It calibrates Energy depending on the recpoint energy.
 623 //      The energy of the reconstructed cluster is corrected with
 624 //      the formula A + B* E  + C* E^2, whose parameters where obtained
 625 //      through the study of the reconstructed energy distribution of
 626 //      monoenergetic photons.
 627
 628   Float_t p[]={0.,0.,0.};
 629   for (Int_t i=0; i<3; i++) p[i] = GetParameterCalibration(i);
 630   Float_t enerec = p[0] +  p[1]*e + p[2]*e*e;
 631   return enerec ;
 632
 633 }
 634
 635 //____________________________________________________________________________
 636 Float_t  AliPHOSPIDv1::GetParameterCpv2Emc(Int_t i, TString axis) const
 637 {
 638   // Get the i-th parameter "CPV-EMC distance" for the specified axis
 639   Float_t param = 0.;
 640   if(i>2 || i<0) {
 641     AliError(Form("Invalid parameter number: %d",i));
 642   } else {
 643     axis.ToLower();
 644     if      (axis == "x")
 645       param = (*fParameters)(1,i);
 646     else if (axis == "z")
 647       param = (*fParameters)(2,i);
 648     else {
 649       AliError(Form("Invalid axis name: %s",axis.Data()));
 650     }
 651   }
 652   return  param;
 653 }
 654
 655 //____________________________________________________________________________
 656 Float_t  AliPHOSPIDv1::GetCpv2EmcDistanceCut(TString axis, Float_t e) const
 657 {
 658   // Get CpvtoEmcDistance Cut depending on the cluster energy, axis and
 659   // Purity-Efficiency point
 660
 661   axis.ToLower();
 662   Float_t p[]={0.,0.,0.};
 663   for (Int_t i=0; i<3; i++) p[i] = GetParameterCpv2Emc(i,axis);
 664   Float_t sig = p[0] + TMath::Exp(p[1] - p[2]*e);
 665   return sig;
 666 }
 667
 668 //____________________________________________________________________________
 669 Float_t  AliPHOSPIDv1::GetEllipseParameter(TString particle, TString param, Float_t e) const
 670 {
 671   // Calculates the parameter param of the ellipse
 672
 673   particle.ToLower();
 674   param.   ToLower();
 675   Float_t p[4]={0.,0.,0.,0.};
 676   Float_t value = 0.0;
 677   for (Int_t i=0; i<4; i++) p[i] = GetParameterToCalculateEllipse(particle,param,i);
 678   if (particle == "photon") {
 679     if      (param.Contains("a"))  e = TMath::Min((Double_t)e,70.);
 680     else if (param.Contains("b"))  e = TMath::Min((Double_t)e,70.);
 681     else if (param.Contains("x0")) e = TMath::Max((Double_t)e,1.1);
 682   }
 683
 684  if (particle == "photon")
 685     value = p[0]/TMath::Sqrt(e) + p[1]*e + p[2]*e*e + p[3];
 686   else if (particle == "pi0")
 687     value = p[0] + p[1]*e + p[2]*e*e;
 688
 689   return value;
 690 }
 691
 692 //_____________________________________________________________________________
 693 Float_t  AliPHOSPIDv1::GetParameterPhotonBoundary (Int_t i) const
 694 {
 695   // Get the parameter "i" to calculate the boundary on the moment M2x
 696   // for photons at high p_T
 697   Float_t param = 0;
 698   if (i>3 || i<0) {
 699     AliError(Form("Wrong parameter number: %d\n",i));
 700   } else
 701     param = (*fParameters)(14,i) ;
 702   return param;
 703 }
 704
 705 //____________________________________________________________________________
 706 Float_t  AliPHOSPIDv1::GetParameterPi0Boundary (Int_t i) const
 707 {
 708   // Get the parameter "i" to calculate the boundary on the moment M2x
 709   // for pi0 at high p_T
 710   Float_t param = 0;
 711   if (i>2 || i<0) {
 712     AliError(Form("Wrong parameter number: %d\n",i));
 713   } else
 714     param = (*fParameters)(15,i) ;
 715   return param;
 716 }
 717
 718 //____________________________________________________________________________
 719 Float_t  AliPHOSPIDv1::GetParameterTimeGate(Int_t i) const
 720 {
 721   // Get TimeGate parameter depending on Purity-Efficiency i:
 722   // i=0 - Low purity, i=1 - Medium purity, i=2 - High purity
 723   Float_t param = 0.;
 724   if(i>2 || i<0) {
 725     AliError(Form("Invalid Efficiency-Purity choice %d",i));
 726   } else
 727     param = (*fParameters)(3,i) ;
 728   return param;
 729 }
 730
 731 //_____________________________________________________________________________
 732 Float_t  AliPHOSPIDv1::GetParameterToCalculateEllipse(TString particle, TString param, Int_t i) const
 733 {
 734   // Get the parameter "i" that is needed to calculate the ellipse
 735   // parameter "param" for the particle "particle" ("photon" or "pi0")
 736
 737   particle.ToLower();
 738   param.   ToLower();
 739   Int_t offset = -1;
 740   if      (particle == "photon")
 741     offset=0;
 742   else if (particle == "pi0")
 743     offset=5;
 744   else
 745     AliError(Form("Wrong particle name: %s (choose from pi0/photon)\n",
 746                   particle.Data()));
 747
 748   Int_t p= -1;
 749   Float_t par = 0;
 750
 751   if     (param.Contains("a")) p=4+offset;
 752   else if(param.Contains("b")) p=5+offset;
 753   else if(param.Contains("c")) p=6+offset;
 754   else if(param.Contains("x0"))p=7+offset;
 755   else if(param.Contains("y0"))p=8+offset;
 756
 757   if      (i>4 || i<0) {
 758     AliError(Form("No parameter with index %d", i)) ;
 759   } else if (p==-1) {
 760     AliError(Form("No parameter with name %s", param.Data() )) ;
 761   } else
 762     par = (*fParameters)(p,i) ;
 763
 764   return par;
 765 }
 766
 767
 768 //DP____________________________________________________________________________
 769 //Float_t  AliPHOSPIDv1::GetDistance(AliPHOSEmcRecPoint * emc,AliPHOSCpvRecPoint * cpv, Option_t *  axis)const
 770 //{
 771 //  // Calculates the distance between the EMC RecPoint and the PPSD RecPoint
 772 //
 773 //  const AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ;
 774 //  TVector3 vecEmc ;
 775 //  TVector3 vecCpv ;
 776 //  if(cpv){
 777 //    emc->GetLocalPosition(vecEmc) ;
 778 //    cpv->GetLocalPosition(vecCpv) ;
 779 //
 780 //    if(emc->GetPHOSMod() == cpv->GetPHOSMod()){
 781 //      // Correct to difference in CPV and EMC position due to different distance to center.
 782 //      // we assume, that particle moves from center
 783 //      Float_t dCPV = geom->GetIPtoOuterCoverDistance();
 784 //      Float_t dEMC = geom->GetIPtoCrystalSurface() ;
 785 //      dEMC         = dEMC / dCPV ;
 786 //      vecCpv = dEMC * vecCpv  - vecEmc ;
 787 //      if (axis == "X") return vecCpv.X();
 788 //      if (axis == "Y") return vecCpv.Y();
 789 //      if (axis == "Z") return vecCpv.Z();
 790 //      if (axis == "R") return vecCpv.Mag();
 791 //    }
 792 //    return 100000000 ;
 793 //  }
 794 //  return 100000000 ;
 795 //}
 796 //____________________________________________________________________________
 797 Int_t  AliPHOSPIDv1::GetCPVBit(AliPHOSTrackSegment * ts, Int_t effPur, Float_t e) const
 798 {
 799   //Calculates the pid bit for the CPV selection per each purity.
 800   if(effPur>2 || effPur<0)
 801     AliError(Form("Invalid Efficiency-Purity choice %d",effPur));
 802
 803 //DP  if(ts->GetCpvIndex()<0)
 804 //DP    return 1 ; //no CPV cluster
 805
 806   Float_t sigX = GetCpv2EmcDistanceCut("X",e);
 807   Float_t sigZ = GetCpv2EmcDistanceCut("Z",e);
 808
 809   Float_t deltaX = TMath::Abs(ts->GetCpvDistance("X"));
 810   Float_t deltaZ = TMath::Abs(ts->GetCpvDistance("Z"));
 811 //  Info("GetCPVBit"," xdist %f, sigx %f, zdist %f, sigz %f",deltaX, sigX, deltaZ,sigZ) ;
 812
 813   //if(deltaX>sigX*(effPur+1))
 814   //if((deltaX>sigX*(effPur+1)) || (deltaZ>sigZ*(effPur+1)))
 815   if((deltaX>sigX*(effPur+1)) && (deltaZ>sigZ*(effPur+1)))
 816     return 1;//Neutral
 817   else
 818     return 0;//Charged
 819 }
 820
 821 //____________________________________________________________________________
 822 Int_t  AliPHOSPIDv1::GetPrincipalBit(TString particle, const Double_t* p, Int_t effPur, Float_t e)const
 823 {
 824   //Is the particle inside de PCA ellipse?
 825
 826   particle.ToLower();
 827   Int_t    prinbit  = 0 ;
 828   Float_t a  = GetEllipseParameter(particle,"a" , e);
 829   Float_t b  = GetEllipseParameter(particle,"b" , e);
 830   Float_t c  = GetEllipseParameter(particle,"c" , e);
 831   Float_t x0 = GetEllipseParameter(particle,"x0", e);
 832   Float_t y0 = GetEllipseParameter(particle,"y0", e);
 833
 834   Float_t r = TMath::Power((p[0] - x0)/a,2) +
 835               TMath::Power((p[1] - y0)/b,2) +
 836             c*(p[0] -  x0)*(p[1] - y0)/(a*b) ;
 837   //3 different ellipses defined
 838   if((effPur==2) && (r<1./2.)) prinbit= 1;
 839   if((effPur==1) && (r<2.   )) prinbit= 1;
 840   if((effPur==0) && (r<9./2.)) prinbit= 1;
 841
 842   if(r<0)
 843     AliError("Negative square?") ;
 844
 845   return prinbit;
 846
 847 }
 848 //____________________________________________________________________________
 849 Int_t  AliPHOSPIDv1::GetHardPhotonBit(AliPHOSEmcRecPoint * emc) const
 850 {
 851   // Set bit for identified hard photons (E > 30 GeV)
 852   // if the second moment M2x is below the boundary
 853
 854   Float_t e   = emc->GetEnergy();
 855   if (e < 30.0) return 0;
 856   Float_t m2x = emc->GetM2x();
 857   Float_t m2xBoundary = GetParameterPhotonBoundary(0) *
 858     TMath::Exp(-TMath::Power(e-GetParameterPhotonBoundary(1),2)/2.0/
 859                 TMath::Power(GetParameterPhotonBoundary(2),2)) +
 860     GetParameterPhotonBoundary(3);
 861   AliDebug(1, Form("GetHardPhotonBit","E=%f, m2x=%f, boundary=%f",
 862                        e,m2x,m2xBoundary));
 863   if (m2x < m2xBoundary)
 864     return 1;// A hard photon
 865   else
 866     return 0;// Not a hard photon
 867 }
 868
 869 //____________________________________________________________________________
 870 Int_t  AliPHOSPIDv1::GetHardPi0Bit(AliPHOSEmcRecPoint * emc) const
 871 {
 872   // Set bit for identified hard pi0  (E > 30 GeV)
 873   // if the second moment M2x is above the boundary
 874
 875   Float_t e   = emc->GetEnergy();
 876   if (e < 30.0) return 0;
 877   Float_t m2x = emc->GetM2x();
 878   Float_t m2xBoundary = GetParameterPi0Boundary(0) +
 879                     e * GetParameterPi0Boundary(1);
 880   AliDebug(1,Form("E=%f, m2x=%f, boundary=%f",e,m2x,m2xBoundary));
 881   if (m2x > m2xBoundary)
 882     return 1;// A hard pi0
 883   else
 884     return 0;// Not a hard pi0
 885 }
 886
 887 //____________________________________________________________________________
 888 TVector3 AliPHOSPIDv1::GetMomentumDirection(AliPHOSEmcRecPoint * emc, AliPHOSCpvRecPoint * )const
 889 {
 890   // Calculates the momentum direction:
 891   //   1. if only a EMC RecPoint, direction is given by IP and this RecPoint
 892   //   2. if a EMC RecPoint and CPV RecPoint, direction is given by the line through the 2 recpoints
 893   //  However because of the poor position resolution of PPSD the direction is always taken as if we were
 894   //  in case 1.
 895
 896   TVector3 dir(0,0,0) ;
 897   TMatrixF  dummy ;
 898
 899   emc->GetGlobalPosition(dir, dummy) ;
 900
 901   TVector3 local ;
 902   emc->GetLocalPosition(local) ;
 903
 904   AliPHOSGeometry * phosgeom = AliPHOSGeometry::GetInstance() ;
 905   //Correct for the non-perpendicular incidence
 906   // Correction for the depth of the shower starting point (TDR p 127)
 907   Float_t para = 0.925 ;
 908   Float_t parb = 6.52 ;
 909
 910   //Remove Old correction (vertex at 0,0,0)
 911   TVector3 vtxOld(0.,0.,0.) ;
 912   TVector3 vInc ;
 913   Float_t x=local.X() ;
 914   Float_t z=local.Z() ;
 915   phosgeom->GetIncidentVector(vtxOld,emc->GetPHOSMod(),x,z,vInc) ;
 916   Float_t depthxOld = 0.;
 917   Float_t depthzOld = 0.;
 918   Float_t energy = emc->GetEnergy() ;
 919   if (energy > 0 && vInc.Y()!=0.) {
 920     depthxOld = ( para * TMath::Log(energy) + parb ) * vInc.X()/vInc.Y() ;
 921     depthzOld = ( para * TMath::Log(energy) + parb ) * vInc.Z()/vInc.Y() ;
 922   }
 923   else{
 924     AliError("Cluster with zero energy \n");
 925   }
 926   //Apply Real vertex
 927   phosgeom->GetIncidentVector(fVtx,emc->GetPHOSMod(),x,z,vInc) ;
 928   Float_t depthx = 0.;
 929   Float_t depthz = 0.;
 930   if (energy > 0 && vInc.Y()!=0.) {
 931     depthx = ( para * TMath::Log(energy) + parb ) * vInc.X()/vInc.Y() ;
 932     depthz = ( para * TMath::Log(energy) + parb ) * vInc.Z()/vInc.Y() ;
 933   }
 934
 935   dir.SetXYZ(dir.X()-(depthxOld-depthx)*TMath::Sin(dir.Phi()),
 936              dir.Y()-(depthxOld-depthx)*TMath::Cos(dir.Phi()),
 937              dir.Z()+depthzOld-depthz) ;
 938
 939   //Correct for the vertex position
 940   dir = dir - fVtx ;
 941   dir.SetMag(1.) ;
 942
 943   return dir ;
 944 }
 945
 946 //________________________________________________________________________
 947 Double_t  AliPHOSPIDv1::LandauF(Double_t  x, Double_t y, Double_t * par)
 948 {
 949   //Given the energy x and the parameter y (tof, shower dispersion or cpv-emc distance),
 950   //this method returns a density probability of this parameter, given by a landau
 951   //function whose parameters depend with the energy  with a function: a/(x*x)+b/x+b
 952
 953   if (x > par[9]) x = par[9];
 954
 955   //Double_t cnt    = par[1] / (x*x) + par[2] / x + par[0] ;
 956   Double_t cnt    = par[0] + par[1] * x + par[2] * x * x ;
 957   Double_t mean   = par[4] / (x*x) + par[5] / x + par[3] ;
 958   Double_t sigma  = par[7] / (x*x) + par[8] / x + par[6] ;
 959
 960   if(TMath::Abs(sigma) > 1.e-10){
 961     return cnt*TMath::Landau(y,mean,sigma);
 962   }
 963   else
 964     return 0.;
 965
 966 }
 967 //________________________________________________________________________
 968 Double_t  AliPHOSPIDv1::LandauPol2(Double_t  x, Double_t y, Double_t * par)
 969 {
 970
 971   //Given the energy x and the parameter y (tof, shower dispersion or cpv-emc distance),
 972   //this method returns a density probability of this parameter, given by a landau
 973   //function whose parameters depend with the energy like second order polinomial
 974
 975   Double_t cnt    = par[2] * (x*x) + par[1] * x + par[0] ;
 976   Double_t mean   = par[5] * (x*x) + par[4] * x + par[3] ;
 977   Double_t sigma  = par[8] * (x*x) + par[7] * x + par[6] ;
 978
 979    if(TMath::Abs(sigma) > 1.e-10){
 980     return cnt*TMath::Landau(y,mean,sigma);
 981   }
 982   else
 983     return 0.;
 984
 985
 986 }
 987 // //________________________________________________________________________
 988 // Double_t  AliPHOSPIDv1::ChargedHadronDistProb(Double_t  x, Double_t y, Double_t * parg, Double_t * parl)
 989 // {
 990 //   Double_t cnt   = 0.0 ;
 991 //   Double_t mean  = 0.0 ;
 992 //   Double_t sigma = 0.0 ;
 993 //   Double_t arg   = 0.0 ;
 994 //   if (y < parl[4] / (x*x) + parl[5] / x + parl[3]){
 995 //     cnt    = parg[1] / (x*x) + parg[2] / x + parg[0] ;
 996 //     mean   = parg[4] / (x*x) + parg[5] / x + parg[3] ;
 997 //     sigma  = parg[7] / (x*x) + parg[8] / x + parg[6] ;
 998 //     TF1 * f = new TF1("gaus","gaus",0.,100.);
 999 //     f->SetParameters(cnt,mean,sigma);
1000 //     arg  = f->Eval(y) ;
1001 //   }
1002 //   else{
1003 //     cnt    = parl[1] / (x*x) + parl[2] / x + parl[0] ;
1004 //     mean   = parl[4] / (x*x) + parl[5] / x + parl[3] ;
1005 //     sigma  = parl[7] / (x*x) + parl[8] / x + parl[6] ;
1006 //     TF1 * f = new TF1("landau","landau",0.,100.);
1007 //     f->SetParameters(cnt,mean,sigma);
1008 //     arg  = f->Eval(y) ;
1009 //   }
1010 //   //  Double_t mean   = par[3] + par[4] * x + par[5] * x * x ;
1011 //   //   Double_t sigma  = par[6] + par[7] * x + par[8] * x * x ;
1012
1013 //   //Double_t arg    = -(y-mean)*(y-mean)/(2*sigma*sigma) ;
1014 //   //return cnt * TMath::Exp(arg) ;
1015
1016 //   return arg;
1017
1018 // }
1019 //____________________________________________________________________________
1020 void  AliPHOSPIDv1::MakePID()
1021 {
1022   // construct the PID weight from a Bayesian Method
1023
1024   const Int_t kSPECIES = AliPID::kSPECIESN ;
1025
1026   AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
1027
1028   Int_t nparticles = gime->RecParticles()->GetEntriesFast() ;
1029
1030   TObjArray * emcRecPoints = gime->EmcRecPoints() ;
1031   TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
1032   TClonesArray * trackSegments = gime->TrackSegments() ;
1033   if ( !emcRecPoints || !cpvRecPoints || !trackSegments ) {
1034     AliFatal("RecPoints or TrackSegments not found !") ;
1035   }
1036   TIter next(trackSegments) ;
1037   AliPHOSTrackSegment * ts ;
1038   Int_t index = 0 ;
1039
1040   Double_t * stof[kSPECIES] ;
1041   Double_t * sdp [kSPECIES]  ;
1042   Double_t * scpv[kSPECIES] ;
1043   Double_t * sw  [kSPECIES] ;
1044   //Info("MakePID","Begin MakePID");
1045
1046   for (Int_t i =0; i< kSPECIES; i++){
1047     stof[i] = new Double_t[nparticles] ;
1048     sdp [i] = new Double_t[nparticles] ;
1049     scpv[i] = new Double_t[nparticles] ;
1050     sw  [i] = new Double_t[nparticles] ;
1051   }
1052
1053
1054   while ( (ts = (AliPHOSTrackSegment *)next()) ) {
1055
1056     //cout<<">>>>>> Bayesian Index "<<index<<endl;
1057
1058     AliPHOSEmcRecPoint * emc = 0 ;
1059     if(ts->GetEmcIndex()>=0)
1060       emc = (AliPHOSEmcRecPoint *) emcRecPoints->At(ts->GetEmcIndex()) ;
1061
1062 //    AliPHOSCpvRecPoint * cpv = 0 ;
1063 //    if(ts->GetCpvIndex()>=0)
1064 //      cpv = (AliPHOSCpvRecPoint *) cpvRecPoints->At(ts->GetCpvIndex()) ;
1065 //
1066 ////     Int_t track = 0 ;
1067 ////     track = ts->GetTrackIndex() ; //TPC tracks ?
1068
1069     if (!emc) {
1070       AliFatal(Form("-> emc(%d) = %d", ts->GetEmcIndex(), emc )) ;
1071     }
1072
1073
1074     // ############Tof#############################
1075
1076     //    Info("MakePID", "TOF");
1077     Float_t  en   = emc->GetEnergy();
1078     Double_t time = emc->GetTime() ;
1079     //    cout<<">>>>>>>Energy "<<en<<"Time "<<time<<endl;
1080
1081     // now get the signals probability
1082     // s(pid) in the Bayesian formulation
1083
1084     stof[AliPID::kPhoton][index]   = 1.;
1085     stof[AliPID::kElectron][index] = 1.;
1086     stof[AliPID::kEleCon][index]   = 1.;
1087     //We assing the same prob to charged hadrons, sum is 1
1088     stof[AliPID::kPion][index]     = 1./3.;
1089     stof[AliPID::kKaon][index]     = 1./3.;
1090     stof[AliPID::kProton][index]   = 1./3.;
1091     //We assing the same prob to neutral hadrons, sum is 1
1092     stof[AliPID::kNeutron][index]  = 1./2.;
1093     stof[AliPID::kKaon0][index]    = 1./2.;
1094     stof[AliPID::kMuon][index]     = 1.;
1095
1096     if(en <  fTOFEnThreshold) {
1097
1098       Double_t pTofPion = fTFpiong ->Eval(time) ; //gaus distribution
1099       Double_t pTofKaon = 0;
1100
1101       if(time < fTkaonl[1])
1102         pTofKaon = fTFkaong  ->Eval(time) ; //gaus distribution
1103       else
1104         pTofKaon = fTFkaonl  ->Eval(time) ; //landau distribution
1105
1106       Double_t pTofNucleon = 0;
1107
1108       if(time < fThhadronl[1])
1109         pTofNucleon = fTFhhadrong   ->Eval(time) ; //gaus distribution
1110       else
1111         pTofNucleon = fTFhhadronl   ->Eval(time) ; //landau distribution
1112       //We assing the same prob to neutral hadrons, sum is the average prob
1113       Double_t pTofNeHadron =  (pTofKaon + pTofNucleon)/2. ;
1114       //We assing the same prob to charged hadrons, sum is the average prob
1115       Double_t pTofChHadron =  (pTofPion + pTofKaon + pTofNucleon)/3. ;
1116
1117       stof[AliPID::kPhoton][index]   = fTFphoton     ->Eval(time) ;
1118       //gaus distribution
1119       stof[AliPID::kEleCon][index]   = stof[AliPID::kPhoton][index] ;
1120       //a conversion electron has the photon ToF
1121       stof[AliPID::kMuon][index]     = stof[AliPID::kPhoton][index] ;
1122
1123       stof[AliPID::kElectron][index] = pTofPion  ;
1124
1125       stof[AliPID::kPion][index]     =  pTofChHadron ;
1126       stof[AliPID::kKaon][index]     =  pTofChHadron ;
1127       stof[AliPID::kProton][index]   =  pTofChHadron ;
1128
1129       stof[AliPID::kKaon0][index]    =  pTofNeHadron ;
1130       stof[AliPID::kNeutron][index]  =  pTofNeHadron ;
1131     }
1132
1133     //    Info("MakePID", "Dispersion");
1134
1135     // ###########Shower shape: Dispersion####################
1136     Float_t dispersion = emc->GetDispersion();
1137     //DP: Correct for non-perpendicular incidence
1138     //DP: still to be done
1139
1140     //dispersion is not well defined if the cluster is only in few crystals
1141
1142     sdp[AliPID::kPhoton][index]   = 1. ;
1143     sdp[AliPID::kElectron][index] = 1. ;
1144     sdp[AliPID::kPion][index]     = 1. ;
1145     sdp[AliPID::kKaon][index]     = 1. ;
1146     sdp[AliPID::kProton][index]   = 1. ;
1147     sdp[AliPID::kNeutron][index]  = 1. ;
1148     sdp[AliPID::kEleCon][index]   = 1. ;
1149     sdp[AliPID::kKaon0][index]    = 1. ;
1150     sdp[AliPID::kMuon][index]     = 1. ;
1151
1152     if(en > fDispEnThreshold && emc->GetMultiplicity() >  fDispMultThreshold){
1153       sdp[AliPID::kPhoton][index]   = GausF(en , dispersion, fDphoton) ;
1154       sdp[AliPID::kElectron][index] = sdp[AliPID::kPhoton][index] ;
1155       sdp[AliPID::kPion][index]     = LandauF(en , dispersion, fDhadron ) ;
1156       sdp[AliPID::kKaon][index]     = sdp[AliPID::kPion][index]  ;
1157       sdp[AliPID::kProton][index]   = sdp[AliPID::kPion][index]  ;
1158       sdp[AliPID::kNeutron][index]  = sdp[AliPID::kPion][index]  ;
1159       sdp[AliPID::kEleCon][index]   = sdp[AliPID::kPhoton][index];
1160       sdp[AliPID::kKaon0][index]    = sdp[AliPID::kPion][index]  ;
1161       sdp[AliPID::kMuon][index]     = fDFmuon ->Eval(dispersion) ;
1162       //landau distribution
1163     }
1164
1165 //      Info("MakePID","multiplicity %d, dispersion %f", emc->GetMultiplicity(), dispersion);
1166 //      Info("MakePID","ss: photon %f, hadron %f ",  sdp[AliPID::kPhoton][index],  sdp[AliPID::kPion][index]);
1167 //       cout<<">>>>>multiplicity "<<emc->GetMultiplicity()<<", dispersion "<< dispersion<<endl ;
1168 //       cout<<"<<<<<ss: photon   "<<sdp[AliPID::kPhoton][index]<<", hadron    "<<sdp[AliPID::kPion][index]<<endl;
1169
1170     //########## CPV-EMC  Distance#######################
1171     //     Info("MakePID", "Distance");
1172
1173     Float_t x             = TMath::Abs(ts->GetCpvDistance("X")) ;
1174     Float_t z             = ts->GetCpvDistance("Z") ;
1175
1176     Double_t pcpv         = 0 ;
1177     Double_t pcpvneutral  = 0. ;
1178
1179     Double_t elprobx      = GausF(en , x, fXelectron) ;
1180     Double_t elprobz      = GausF(en , z, fZelectron) ;
1181     Double_t chprobx      = GausF(en , x, fXcharged)  ;
1182     Double_t chprobz      = GausF(en , z, fZcharged)  ;
1183     Double_t pcpvelectron = elprobx * elprobz;
1184     Double_t pcpvcharged  = chprobx * chprobz;
1185
1186 //     cout<<">>>>energy "<<en<<endl;
1187 //     cout<<">>>>electron : x "<<x<<" xprob "<<elprobx<<" z "<<z<<" zprob "<<elprobz<<endl;
1188 //     cout<<">>>>hadron   : x "<<x<<" xprob "<<chprobx<<" z "<<z<<" zprob "<<chprobz<<endl;
1189 //     cout<<">>>>electron : px*pz "<<pcpvelectron <<" hadron: px*pz "<<pcpvcharged<<endl;
1190
1191     // Is neutral or charged?
1192     if(pcpvelectron >= pcpvcharged)
1193       pcpv = pcpvelectron ;
1194     else
1195       pcpv = pcpvcharged ;
1196
1197     if(pcpv < fChargedNeutralThreshold)
1198       {
1199         pcpvneutral  = 1. ;
1200         pcpvcharged  = 0. ;
1201         pcpvelectron = 0. ;
1202       }
1203     //    else
1204     //      cout<<">>>>>>>>>>>CHARGED>>>>>>>>>>>"<<endl;
1205
1206     scpv[AliPID::kPion][index]     =  pcpvcharged  ;
1207     scpv[AliPID::kKaon][index]     =  pcpvcharged  ;
1208     scpv[AliPID::kProton][index]   =  pcpvcharged  ;
1209
1210     scpv[AliPID::kMuon][index]     =  pcpvelectron ;
1211     scpv[AliPID::kElectron][index] =  pcpvelectron ;
1212     scpv[AliPID::kEleCon][index]   =  pcpvelectron ;
1213
1214     scpv[AliPID::kPhoton][index]   =  pcpvneutral  ;
1215     scpv[AliPID::kNeutron][index]  =  pcpvneutral  ;
1216     scpv[AliPID::kKaon0][index]    =  pcpvneutral  ;
1217
1218
1219     //   Info("MakePID", "CPV passed");
1220
1221     //############## Pi0 #############################
1222     stof[AliPID::kPi0][index]      = 0. ;
1223     scpv[AliPID::kPi0][index]      = 0. ;
1224     sdp [AliPID::kPi0][index]      = 0. ;
1225
1226     if(en > 30.){
1227       // pi0 are detected via decay photon
1228       stof[AliPID::kPi0][index]  =   stof[AliPID::kPhoton][index];
1229       scpv[AliPID::kPi0][index]  = pcpvneutral  ;
1230       if(emc->GetMultiplicity() >  fDispMultThreshold)
1231         sdp [AliPID::kPi0][index]  = GausF(en , dispersion, fDpi0) ;
1232         //sdp [AliPID::kPi0][index]  = GausPol2(en , dispersion, fDpi0) ;
1233 //       cout<<"E = "<<en<<" GeV; disp = "<<dispersion<<"; mult = "
1234 //        <<emc->GetMultiplicity()<<endl;
1235 //       cout<<"PDF: photon = "<<sdp [AliPID::kPhoton][index]<<"; pi0 = "
1236 //        <<sdp [AliPID::kPi0][index]<<endl;
1237     }
1238
1239
1240
1241
1242     //############## muon #############################
1243
1244     if(en > 0.5){
1245       //Muons deposit few energy
1246       scpv[AliPID::kMuon][index]     =  0 ;
1247       stof[AliPID::kMuon][index]     =  0 ;
1248       sdp [AliPID::kMuon][index]     =  0 ;
1249     }
1250
1251     //Weight to apply to hadrons due to energy reconstruction
1252
1253     Float_t weight = fERecWeight ->Eval(en) ;
1254
1255     sw[AliPID::kPhoton][index]   = 1. ;
1256     sw[AliPID::kElectron][index] = 1. ;
1257     sw[AliPID::kPion][index]     = weight ;
1258     sw[AliPID::kKaon][index]     = weight ;
1259     sw[AliPID::kProton][index]   = weight ;
1260     sw[AliPID::kNeutron][index]  = weight ;
1261     sw[AliPID::kEleCon][index]   = 1. ;
1262     sw[AliPID::kKaon0][index]    = weight ;
1263     sw[AliPID::kMuon][index]     = weight ;
1264     sw[AliPID::kPi0][index]      = 1. ;
1265
1266 //     if(en > 0.5){
1267 //       cout<<"######################################################"<<endl;
1268 //       //cout<<"MakePID: energy "<<en<<", tof "<<time<<", distance "<<distance<<", dispersion "<<dispersion<<endl ;
1269 //       cout<<"MakePID: energy "<<en<<", tof "<<time<<", dispersion "<<dispersion<<", x "<<x<<", z "<<z<<endl ;
1270 //       cout<<">>>>>multiplicity "<<emc->GetMultiplicity()<<endl;
1271 //       cout<<">>>>electron : xprob "<<elprobx<<" zprob "<<elprobz<<endl;
1272 //       cout<<">>>>hadron   : xprob "<<chprobx<<" zprob "<<chprobz<<endl;
1273 //       cout<<">>>>electron : px*pz "<<pcpvelectron <<" hadron: px*pz "<<pcpvcharged<<endl;
1274
1275 //        cout<<"Photon   , pid "<< fInitPID[AliPID::kPhoton]<<" tof "<<stof[AliPID::kPhoton][index]
1276 //        <<", cpv "<<scpv[AliPID::kPhoton][index]<<", ss "<<sdp[AliPID::kPhoton][index]<<endl;
1277 //       cout<<"EleCon   , pid "<< fInitPID[AliPID::kEleCon]<<", tof "<<stof[AliPID::kEleCon][index]
1278 //        <<", cpv "<<scpv[AliPID::kEleCon][index]<<" ss "<<sdp[AliPID::kEleCon][index]<<endl;
1279 //       cout<<"Electron , pid "<< fInitPID[AliPID::kElectron]<<", tof "<<stof[AliPID::kElectron][index]
1280 //        <<", cpv "<<scpv[AliPID::kElectron][index]<<" ss "<<sdp[AliPID::kElectron][index]<<endl;
1281 //       cout<<"Muon     , pid "<< fInitPID[AliPID::kMuon]<<", tof "<<stof[AliPID::kMuon][index]
1282 //        <<", cpv "<<scpv[AliPID::kMuon][index]<<" ss "<<sdp[AliPID::kMuon][index]<<endl;
1283 //        cout<<"Pi0      , pid "<< fInitPID[AliPID::kPi0]<<", tof "<<stof[AliPID::kPi0][index]
1284 //        <<", cpv "<<scpv[AliPID::kPi0][index]<<" ss "<<sdp[AliPID::kPi0][index]<<endl;
1285 //       cout<<"Pion     , pid "<< fInitPID[AliPID::kPion]<<", tof "<<stof[AliPID::kPion][index]
1286 //        <<", cpv "<<scpv[AliPID::kPion][index]<<" ss "<<sdp[AliPID::kPion][index]<<endl;
1287 //       cout<<"Kaon0    , pid "<< fInitPID[AliPID::kKaon0]<<", tof "<<stof[AliPID::kKaon0][index]
1288 //        <<", cpv "<<scpv[AliPID::kKaon0][index]<<" ss "<<sdp[AliPID::kKaon0][index]<<endl;
1289 //       cout<<"Kaon     , pid "<< fInitPID[AliPID::kKaon]<<", tof "<<stof[AliPID::kKaon][index]
1290 //        <<", cpv "<<scpv[AliPID::kKaon][index]<<" ss "<<sdp[AliPID::kKaon][index]<<endl;
1291 //       cout<<"Neutron  , pid "<< fInitPID[AliPID::kNeutron]<<", tof "<<stof[AliPID::kNeutron][index]
1292 //        <<", cpv "<<scpv[AliPID::kNeutron][index]<<" ss "<<sdp[AliPID::kNeutron][index]<<endl;
1293 //       cout<<"Proton   , pid "<< fInitPID[AliPID::kProton]<<", tof "<<stof[AliPID::kProton][index]
1294 //        <<", cpv "<<scpv[AliPID::kProton][index]<<" ss "<<sdp[AliPID::kProton][index]<<endl;
1295 //       cout<<"######################################################"<<endl;
1296 //     }
1297       index++;
1298   }
1299
1300   //for (index = 0 ; index < kSPECIES ; index++)
1301   // pid[index] /= nparticles ;
1302
1303
1304   //  Info("MakePID", "Total Probability calculation");
1305
1306   for(index = 0 ; index < nparticles ; index ++) {
1307
1308     AliPHOSRecParticle * recpar = gime->RecParticle(index) ;
1309
1310     //Conversion electron?
1311
1312     if(recpar->IsEleCon()){
1313       fInitPID[AliPID::kEleCon]   = 1. ;
1314       fInitPID[AliPID::kPhoton]   = 0. ;
1315       fInitPID[AliPID::kElectron] = 0. ;
1316     }
1317     else{
1318       fInitPID[AliPID::kEleCon]   = 0. ;
1319       fInitPID[AliPID::kPhoton]   = 1. ;
1320       fInitPID[AliPID::kElectron] = 1. ;
1321     }
1322     //  fInitPID[AliPID::kEleCon]   = 0. ;
1323
1324
1325     // calculates the Bayesian weight
1326
1327     Int_t jndex ;
1328     Double_t wn = 0.0 ;
1329     for (jndex = 0 ; jndex < kSPECIES ; jndex++)
1330       wn += stof[jndex][index] * sdp[jndex][index]  * scpv[jndex][index] *
1331         sw[jndex][index] * fInitPID[jndex] ;
1332
1333     //    cout<<"*************wn "<<wn<<endl;
1334     if (TMath::Abs(wn)>0)
1335       for (jndex = 0 ; jndex < kSPECIES ; jndex++) {
1336         //cout<<"jndex "<<jndex<<" wn "<<wn<<" SetPID * wn"
1337         //<<stof[jndex][index] * sdp[jndex][index] * pid[jndex]  << endl;
1338         //cout<<" tof "<<stof[jndex][index] << " disp " <<sdp[jndex][index] << " pid "<< fInitPID[jndex] << endl;
1339         //      if(jndex ==  AliPID::kPi0 || jndex ==  AliPID::kPhoton){
1340         //        cout<<"Particle "<<jndex<<"  final prob * wn   "
1341         //            <<stof[jndex][index] * sdp[jndex][index] * scpv[jndex][index] *
1342         //          fInitPID[jndex] <<"  wn  "<< wn<<endl;
1343         //        cout<<"pid "<< fInitPID[jndex]<<", tof "<<stof[jndex][index]
1344         //            <<", cpv "<<scpv[jndex][index]<<" ss "<<sdp[jndex][index]<<endl;
1345         //      }
1346         recpar->SetPID(jndex, stof[jndex][index] * sdp[jndex][index] *
1347                        sw[jndex][index] * scpv[jndex][index] *
1348                        fInitPID[jndex] / wn) ;
1349       }
1350   }
1351   //  Info("MakePID", "Delete");
1352
1353   for (Int_t i =0; i< kSPECIES; i++){
1354     delete [] stof[i];
1355     delete [] sdp [i];
1356     delete [] scpv[i];
1357     delete [] sw  [i];
1358   }
1359   //  Info("MakePID","End MakePID");
1360 }
1361
1362 //____________________________________________________________________________
1363 void  AliPHOSPIDv1::MakeRecParticles()
1364 {
1365   // Makes a RecParticle out of a TrackSegment
1366
1367   AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
1368   TObjArray * emcRecPoints = gime->EmcRecPoints() ;
1369   TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
1370   TClonesArray * trackSegments = gime->TrackSegments() ;
1371   if ( !emcRecPoints || !cpvRecPoints || !trackSegments ) {
1372     AliFatal("RecPoints or TrackSegments not found !") ;
1373   }
1374   TClonesArray * recParticles  = gime->RecParticles() ;
1375   recParticles->Clear();
1376
1377   TIter next(trackSegments) ;
1378   AliPHOSTrackSegment * ts ;
1379   Int_t index = 0 ;
1380   AliPHOSRecParticle * rp ;
1381   while ( (ts = (AliPHOSTrackSegment *)next()) ) {
1382     //  cout<<">>>>>>>>>>>>>>>PCA Index "<<index<<endl;
1383     new( (*recParticles)[index] ) AliPHOSRecParticle() ;
1384     rp = (AliPHOSRecParticle *)recParticles->At(index) ;
1385     rp->SetTrackSegment(index) ;
1386     rp->SetIndexInList(index) ;
1387
1388     AliPHOSEmcRecPoint * emc = 0 ;
1389     if(ts->GetEmcIndex()>=0)
1390       emc = (AliPHOSEmcRecPoint *) emcRecPoints->At(ts->GetEmcIndex()) ;
1391
1392     AliPHOSCpvRecPoint * cpv = 0 ;
1393     if(ts->GetCpvIndex()>=0)
1394       cpv = (AliPHOSCpvRecPoint *) cpvRecPoints->At(ts->GetCpvIndex()) ;
1395
1396     Int_t track = 0 ;
1397     track = ts->GetTrackIndex() ;
1398
1399     // Now set type (reconstructed) of the particle
1400
1401     // Choose the cluster energy range
1402
1403     if (!emc) {
1404       AliFatal(Form("-> emc(%d) = %d", ts->GetEmcIndex(), emc )) ;
1405     }
1406
1407     Float_t e = emc->GetEnergy() ;
1408
1409     Float_t  lambda[2] ;
1410     emc->GetElipsAxis(lambda) ;
1411
1412     if((lambda[0]>0.01) && (lambda[1]>0.01)){
1413       // Looking PCA. Define and calculate the data (X),
1414       // introduce in the function X2P that gives the components (P).
1415
1416       Float_t  spher = 0. ;
1417       Float_t  emaxdtotal = 0. ;
1418
1419       if((lambda[0]+lambda[1])!=0)
1420         spher=TMath::Abs(lambda[0]-lambda[1])/(lambda[0]+lambda[1]);
1421
1422       emaxdtotal=emc->GetMaximalEnergy()/emc->GetEnergy();
1423
1424       fX[0] = lambda[0] ;
1425       fX[1] = lambda[1] ;
1426       fX[2] = emc->GetDispersion() ;
1427       fX[3] = spher ;
1428       fX[4] = emc->GetMultiplicity() ;
1429       fX[5] = emaxdtotal ;
1430       fX[6] = emc->GetCoreEnergy() ;
1431
1432       fPrincipalPhoton->X2P(fX,fPPhoton);
1433       fPrincipalPi0   ->X2P(fX,fPPi0);
1434
1435     }
1436     else{
1437       fPPhoton[0]=-100.0;  //We do not accept clusters with
1438       fPPhoton[1]=-100.0;  //one cell as a photon-like
1439       fPPi0[0]   =-100.0;
1440       fPPi0[1]   =-100.0;
1441     }
1442
1443     Float_t time = emc->GetTime() ;
1444     rp->SetTof(time) ;
1445
1446     // Loop of Efficiency-Purity (the 3 points of purity or efficiency
1447     // are taken into account to set the particle identification)
1448     for(Int_t effPur = 0; effPur < 3 ; effPur++){
1449
1450       // Looking at the CPV detector. If RCPV greater than CpvEmcDistance,
1451       // 1st,2nd or 3rd bit (depending on the efficiency-purity point )
1452       // is set to 1
1453       if(GetCPVBit(ts, effPur,e) == 1 ){
1454         rp->SetPIDBit(effPur) ;
1455         //cout<<"CPV bit "<<effPur<<endl;
1456       }
1457       // Looking the TOF. If TOF smaller than gate,  4th, 5th or 6th
1458       // bit (depending on the efficiency-purity point )is set to 1
1459       if(time< (*fParameters)(3,effPur))
1460         rp->SetPIDBit(effPur+3) ;
1461
1462       //Photon PCA
1463       //If we are inside the ellipse, 7th, 8th or 9th
1464       // bit (depending on the efficiency-purity point )is set to 1
1465       if(GetPrincipalBit("photon",fPPhoton,effPur,e) == 1)
1466         rp->SetPIDBit(effPur+6) ;
1467
1468       //Pi0 PCA
1469       //If we are inside the ellipse, 10th, 11th or 12th
1470       // bit (depending on the efficiency-purity point )is set to 1
1471       if(GetPrincipalBit("pi0"   ,fPPi0   ,effPur,e) == 1)
1472         rp->SetPIDBit(effPur+9) ;
1473     }
1474     if(GetHardPhotonBit(emc))
1475       rp->SetPIDBit(12) ;
1476     if(GetHardPi0Bit   (emc))
1477       rp->SetPIDBit(13) ;
1478
1479     if(track >= 0)
1480       rp->SetPIDBit(14) ;
1481
1482     //Set momentum, energy and other parameters
1483     Float_t  encal = GetCalibratedEnergy(e);
1484     TVector3 dir   = GetMomentumDirection(emc,cpv) ;
1485     dir.SetMag(encal) ;
1486     rp->SetMomentum(dir.X(),dir.Y(),dir.Z(),encal) ;
1487     rp->SetCalcMass(0);
1488     rp->Name(); //If photon sets the particle pdg name to gamma
1489     rp->SetProductionVertex(fVtx.X(),fVtx.Y(),fVtx.Z(),0);
1490     rp->SetFirstMother(-1);
1491     rp->SetLastMother(-1);
1492     rp->SetFirstDaughter(-1);
1493     rp->SetLastDaughter(-1);
1494     rp->SetPolarisation(0,0,0);
1495     //Set the position in global coordinate system from the RecPoint
1496     AliPHOSGeometry * geom = gime->PHOSGeometry() ;
1497     AliPHOSTrackSegment * ts  = gime->TrackSegment(rp->GetPHOSTSIndex()) ;
1498     AliPHOSEmcRecPoint  * erp = gime->EmcRecPoint(ts->GetEmcIndex()) ;
1499     TVector3 pos ;
1500     geom->GetGlobal(erp, pos) ;
1501     rp->SetPos(pos);
1502     index++ ;
1503   }
1504 }
1505
1506 //____________________________________________________________________________
1507 void  AliPHOSPIDv1::Print(const Option_t *) const
1508 {
1509   // Print the parameters used for the particle type identification
1510
1511     AliInfo("=============== AliPHOSPIDv1 ================") ;
1512     printf("Making PID\n") ;
1513     printf("    Pricipal analysis file from 0.5 to 100 %s\n", fFileNamePrincipalPhoton.Data() )   ;
1514     printf("    Name of parameters file     %s\n", fFileNameParameters.Data() )  ;
1515     printf("    Matrix of Parameters: 14x4\n") ;
1516     printf("        Energy Calibration  1x3 [3 parametres to calibrate energy: A + B* E + C * E^2]\n") ;
1517     printf("        RCPV 2x3 rows x and z, columns function cut parameters\n") ;
1518     printf("        TOF  1x3 [High Eff-Low Pur,Medium Eff-Pur, Low Eff-High Pur]\n") ;
1519     printf("        PCA  5x4 [5 ellipse parametres and 4 parametres to calculate them: A/Sqrt(E) + B* E + C * E^2 + D]\n") ;
1520     printf("    Pi0 PCA  5x3 [5 ellipse parametres and 3 parametres to calculate them: A + B* E + C * E^2]\n") ;
1521     fParameters->Print() ;
1522 }
1523
1524
1525
1526 //____________________________________________________________________________
1527 void AliPHOSPIDv1::PrintRecParticles(Option_t * option)
1528 {
1529   // Print table of reconstructed particles
1530
1531   AliPHOSGetter *gime = AliPHOSGetter::Instance() ;
1532
1533   TClonesArray * recParticles = gime->RecParticles() ;
1534
1535   TString message ;
1536   message  = "\nevent " ;
1537   message += gime->EventNumber();
1538   message += "       found " ;
1539   message += recParticles->GetEntriesFast();
1540   message += " RecParticles\n" ;
1541
1542   if(strstr(option,"all")) {  // printing found TS
1543     message += "\n  PARTICLE         Index    \n" ;
1544
1545     Int_t index ;
1546     for (index = 0 ; index < recParticles->GetEntries() ; index++) {
1547       AliPHOSRecParticle * rp = (AliPHOSRecParticle * ) recParticles->At(index) ;
1548       message += "\n" ;
1549       message += rp->Name().Data() ;
1550       message += " " ;
1551       message += rp->GetIndexInList() ;
1552       message += " " ;
1553       message += rp->GetType()  ;
1554     }
1555   }
1556   AliInfo(message.Data() ) ;
1557 }
1558
1559 //____________________________________________________________________________
1560 void  AliPHOSPIDv1::SetParameters()
1561 {
1562   // PCA : To do the Principal Components Analysis it is necessary
1563   // the Principal file, which is opened here
1564   fX       = new double[7]; // Data for the PCA
1565   fPPhoton = new double[7]; // Eigenvalues of the PCA
1566   fPPi0    = new double[7]; // Eigenvalues of the Pi0 PCA
1567
1568   // Read photon principals from the photon file
1569
1570   fFileNamePrincipalPhoton = "$ALICE_ROOT/PHOS/PCA8pa15_0.5-100.root" ;
1571   TFile f( fFileNamePrincipalPhoton.Data(), "read" ) ;
1572   fPrincipalPhoton = dynamic_cast<TPrincipal*> (f.Get("principal")) ;
1573   f.Close() ;
1574
1575   // Read pi0 principals from the pi0 file
1576
1577   fFileNamePrincipalPi0    = "$ALICE_ROOT/PHOS/PCA_pi0_40-120.root" ;
1578   TFile fPi0( fFileNamePrincipalPi0.Data(), "read" ) ;
1579   fPrincipalPi0    = dynamic_cast<TPrincipal*> (fPi0.Get("principal")) ;
1580   fPi0.Close() ;
1581
1582   // Open parameters file and initialization of the Parameters matrix.
1583   // In the File Parameters.dat are all the parameters. These are introduced
1584   // in a matrix of 16x4
1585   //
1586   // All the parameters defined in this file are, in order of row:
1587   // line   0   : calibration
1588   // lines  1,2 : CPV rectangular cat for X and Z
1589   // line   3   : TOF cut
1590   // lines  4-8 : parameters to calculate photon PCA ellipse
1591   // lines  9-13: parameters to calculate pi0 PCA ellipse
1592   // lines 14-15: parameters to calculate border for high-pt photons and pi0
1593
1594   fFileNameParameters = gSystem->ExpandPathName("$ALICE_ROOT/PHOS/Parameters.dat");
1595   fParameters = new TMatrixF(16,4) ;
1596   const Int_t kMaxLeng=255;
1597   char string[kMaxLeng];
1598
1599   // Open a text file with PID parameters
1600   FILE *fd = fopen(fFileNameParameters.Data(),"r");
1601   if (!fd)
1602     AliFatal(Form("File %s with a PID parameters cannot be opened\n",
1603           fFileNameParameters.Data()));
1604
1605   Int_t i=0;
1606   // Read parameter file line-by-line and skip empty line and comments
1607   while (fgets(string,kMaxLeng,fd) != NULL) {
1608     if (string[0] == '\n' ) continue;
1609     if (string[0] == '!'  ) continue;
1610     sscanf(string, "%f %f %f %f",
1611            &(*fParameters)(i,0), &(*fParameters)(i,1),
1612            &(*fParameters)(i,2), &(*fParameters)(i,3));
1613     i++;
1614     AliDebug(1, Form("SetParameters", "line %d: %s",i,string));
1615   }
1616   fclose(fd);
1617 }
1618
1619 //____________________________________________________________________________
1620 void  AliPHOSPIDv1::SetParameterCalibration(Int_t i,Float_t param)
1621 {
1622   // Set parameter "Calibration" i to a value param
1623   if(i>2 || i<0) {
1624     AliError(Form("Invalid parameter number: %d",i));
1625   } else
1626     (*fParameters)(0,i) = param ;
1627 }
1628
1629 //____________________________________________________________________________
1630 void  AliPHOSPIDv1::SetParameterCpv2Emc(Int_t i, TString axis, Float_t cut)
1631 {
1632   // Set the parameters to calculate Cpv-to-Emc Distance Cut depending on
1633   // Purity-Efficiency point i
1634
1635   if(i>2 || i<0) {
1636     AliError(Form("Invalid parameter number: %d",i));
1637   } else {
1638     axis.ToLower();
1639     if      (axis == "x") (*fParameters)(1,i) = cut;
1640     else if (axis == "z") (*fParameters)(2,i) = cut;
1641     else {
1642       AliError(Form("Invalid axis name: %s",axis.Data()));
1643     }
1644   }
1645 }
1646
1647 //____________________________________________________________________________
1648 void  AliPHOSPIDv1::SetParameterPhotonBoundary(Int_t i,Float_t param)
1649 {
1650   // Set parameter "Hard photon boundary" i to a value param
1651   if(i>4 || i<0) {
1652     AliError(Form("Invalid parameter number: %d",i));
1653   } else
1654     (*fParameters)(14,i) = param ;
1655 }
1656
1657 //____________________________________________________________________________
1658 void  AliPHOSPIDv1::SetParameterPi0Boundary(Int_t i,Float_t param)
1659 {
1660   // Set parameter "Hard pi0 boundary" i to a value param
1661   if(i>1 || i<0) {
1662     AliError(Form("Invalid parameter number: %d",i));
1663   } else
1664     (*fParameters)(15,i) = param ;
1665 }
1666
1667 //_____________________________________________________________________________
1668 void  AliPHOSPIDv1::SetParameterTimeGate(Int_t i, Float_t gate)
1669 {
1670   // Set the parameter TimeGate depending on Purity-Efficiency point i
1671   if (i>2 || i<0) {
1672     AliError(Form("Invalid Efficiency-Purity choice %d",i));
1673   } else
1674     (*fParameters)(3,i)= gate ;
1675 }
1676
1677 //_____________________________________________________________________________
1678 void  AliPHOSPIDv1::SetParameterToCalculateEllipse(TString particle, TString param, Int_t i, Float_t par)
1679 {
1680   // Set the parameter "i" that is needed to calculate the ellipse
1681   // parameter "param" for a particle "particle"
1682
1683   particle.ToLower();
1684   param.   ToLower();
1685   Int_t p= -1;
1686   Int_t offset=0;
1687
1688   if      (particle == "photon") offset=0;
1689   else if (particle == "pi0")    offset=5;
1690   else
1691     AliError(Form("Wrong particle name: %s (choose from pi0/photon)\n",
1692                   particle.Data()));
1693
1694   if     (param.Contains("a")) p=4+offset;
1695   else if(param.Contains("b")) p=5+offset;
1696   else if(param.Contains("c")) p=6+offset;
1697   else if(param.Contains("x0"))p=7+offset;
1698   else if(param.Contains("y0"))p=8+offset;
1699   if((i>4)||(i<0)) {
1700     AliError(Form("No parameter with index %d", i)) ;
1701   } else if(p==-1) {
1702     AliError(Form("No parameter with name %s", param.Data() )) ;
1703   } else
1704     (*fParameters)(p,i) = par ;
1705 }
1706
1707 //____________________________________________________________________________
1708 void AliPHOSPIDv1::Unload()
1709 {
1710   //Unloads RecPoints, Tracks and RecParticles
1711   AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
1712   gime->PhosLoader()->UnloadRecPoints() ;
1713   gime->PhosLoader()->UnloadTracks() ;
1714   gime->PhosLoader()->UnloadRecParticles() ;
1715 }
1716
1717 //____________________________________________________________________________
1718 void  AliPHOSPIDv1::WriteRecParticles()
1719 {
1720   //It writes reconstructed particles and pid to file
1721
1722   AliPHOSGetter *gime = AliPHOSGetter::Instance() ;
1723
1724   TClonesArray * recParticles = gime->RecParticles() ;
1725   recParticles->Expand(recParticles->GetEntriesFast() ) ;
1726   if(fWrite){
1727     TTree * treeP =  gime->TreeP();
1728
1729     //First rp
1730     Int_t bufferSize = 32000 ;
1731     TBranch * rpBranch = treeP->Branch("PHOSRP",&recParticles,bufferSize);
1732     rpBranch->SetTitle(BranchName());
1733
1734     rpBranch->Fill() ;
1735
1736     gime->WriteRecParticles("OVERWRITE");
1737     gime->WritePID("OVERWRITE");
1738   }
1739 }
1740 //____________________________________________________________________________
1741 void AliPHOSPIDv1::GetVertex(void)
1742 { //extract vertex either using ESD or generator
1743
1744   //Try to extract vertex from data
1745   if(fESD){
1746     const AliESDVertex *esdVtx = fESD->GetVertex() ;
1747     if(esdVtx){
1748       fVtx.SetXYZ(esdVtx->GetXv(),esdVtx->GetYv(),esdVtx->GetZv()) ;
1749       return ;
1750     }
1751   }
1752
1753   AliWarning("Can not read vertex from data, use fixed \n") ;
1754   fVtx.SetXYZ(0.,0.,0.) ;
1755
1756 }
1757 //_______________________________________________________________________
1758 void AliPHOSPIDv1::SetInitPID(const Double_t *p) {
1759   // Sets values for the initial population of each particle type
1760   for (Int_t i=0; i<AliPID::kSPECIESN; i++) fInitPID[i] = p[i];
1761 }
1762 //_______________________________________________________________________
1763 void AliPHOSPIDv1::GetInitPID(Double_t *p) const {
1764   // Gets values for the initial population of each particle type
1765   for (Int_t i=0; i<AliPID::kSPECIESN; i++) p[i] = fInitPID[i];
1766 }