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