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