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