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