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