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