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