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::GetCalibratedEnergy(Float_t e) const
 595 {
 596 //      It calibrates Energy depending on the recpoint energy.
 597 //      The energy of the reconstructed cluster is corrected with
 598 //      the formula A + B* E  + C* E^2, whose parameters where obtained
 599 //      through the study of the reconstructed energy distribution of
 600 //      monoenergetic photons.
 601
 602   if(!fEnergyCorrectionOn) return e;
 603
 604   Float_t p[]={0.,0.,0.};
 605   for (Int_t i=0; i<3; i++) p[i] = GetParameterCalibration(i);
 606   Float_t enerec = p[0] +  p[1]*e + p[2]*e*e;
 607   return enerec ;
 608
 609 }
 610
 611 //____________________________________________________________________________
 612 Float_t  AliPHOSPIDv1::GetParameterCpv2Emc(Int_t i, TString axis) const
 613 {
 614   // Get the i-th parameter "CPV-EMC distance" for the specified axis
 615   Float_t param = 0.;
 616   if(i>2 || i<0) {
 617     AliError(Form("Invalid parameter number: %d",i));
 618   } else {
 619     axis.ToLower();
 620     if      (axis == "x")
 621       param = (*fParameters)(1,i);
 622     else if (axis == "z")
 623       param = (*fParameters)(2,i);
 624     else {
 625       AliError(Form("Invalid axis name: %s",axis.Data()));
 626     }
 627   }
 628   return  param;
 629 }
 630
 631 //____________________________________________________________________________
 632 Float_t  AliPHOSPIDv1::GetCpv2EmcDistanceCut(TString axis, Float_t e) const
 633 {
 634   // Get CpvtoEmcDistance Cut depending on the cluster energy, axis and
 635   // Purity-Efficiency point
 636
 637   axis.ToLower();
 638   Float_t p[]={0.,0.,0.};
 639   for (Int_t i=0; i<3; i++) p[i] = GetParameterCpv2Emc(i,axis);
 640   Float_t sig = p[0] + TMath::Exp(p[1] - p[2]*e);
 641   return sig;
 642 }
 643
 644 //____________________________________________________________________________
 645 Float_t  AliPHOSPIDv1::GetEllipseParameter(TString particle, TString param, Float_t e) const
 646 {
 647   // Calculates the parameter param of the ellipse
 648
 649   particle.ToLower();
 650   param.   ToLower();
 651   Float_t p[4]={0.,0.,0.,0.};
 652   Float_t value = 0.0;
 653   for (Int_t i=0; i<4; i++) p[i] = GetParameterToCalculateEllipse(particle,param,i);
 654   if (particle == "photon") {
 655     if      (param.Contains("a"))  e = TMath::Min((Double_t)e,70.);
 656     else if (param.Contains("b"))  e = TMath::Min((Double_t)e,70.);
 657     else if (param.Contains("x0")) e = TMath::Max((Double_t)e,1.1);
 658   }
 659
 660  if (particle == "photon")
 661     value = p[0]/TMath::Sqrt(e) + p[1]*e + p[2]*e*e + p[3];
 662   else if (particle == "pi0")
 663     value = p[0] + p[1]*e + p[2]*e*e;
 664
 665   return value;
 666 }
 667
 668 //_____________________________________________________________________________
 669 Float_t  AliPHOSPIDv1::GetParameterPhotonBoundary (Int_t i) const
 670 {
 671   // Get the parameter "i" to calculate the boundary on the moment M2x
 672   // for photons at high p_T
 673   Float_t param = 0;
 674   if (i>3 || i<0) {
 675     AliError(Form("Wrong parameter number: %d\n",i));
 676   } else
 677     param = (*fParameters)(14,i) ;
 678   return param;
 679 }
 680
 681 //____________________________________________________________________________
 682 Float_t  AliPHOSPIDv1::GetParameterPi0Boundary (Int_t i) const
 683 {
 684   // Get the parameter "i" to calculate the boundary on the moment M2x
 685   // for pi0 at high p_T
 686   Float_t param = 0;
 687   if (i>2 || i<0) {
 688     AliError(Form("Wrong parameter number: %d\n",i));
 689   } else
 690     param = (*fParameters)(15,i) ;
 691   return param;
 692 }
 693
 694 //____________________________________________________________________________
 695 Float_t  AliPHOSPIDv1::GetParameterTimeGate(Int_t i) const
 696 {
 697   // Get TimeGate parameter depending on Purity-Efficiency i:
 698   // i=0 - Low purity, i=1 - Medium purity, i=2 - High purity
 699   Float_t param = 0.;
 700   if(i>2 || i<0) {
 701     AliError(Form("Invalid Efficiency-Purity choice %d",i));
 702   } else
 703     param = (*fParameters)(3,i) ;
 704   return param;
 705 }
 706
 707 //_____________________________________________________________________________
 708 Float_t  AliPHOSPIDv1::GetParameterToCalculateEllipse(TString particle, TString param, Int_t i) const
 709 {
 710   // Get the parameter "i" that is needed to calculate the ellipse
 711   // parameter "param" for the particle "particle" ("photon" or "pi0")
 712
 713   particle.ToLower();
 714   param.   ToLower();
 715   Int_t offset = -1;
 716   if      (particle == "photon")
 717     offset=0;
 718   else if (particle == "pi0")
 719     offset=5;
 720   else
 721     AliError(Form("Wrong particle name: %s (choose from pi0/photon)\n",
 722                   particle.Data()));
 723
 724   Int_t p= -1;
 725   Float_t par = 0;
 726
 727   if     (param.Contains("a")) p=4+offset;
 728   else if(param.Contains("b")) p=5+offset;
 729   else if(param.Contains("c")) p=6+offset;
 730   else if(param.Contains("x0"))p=7+offset;
 731   else if(param.Contains("y0"))p=8+offset;
 732
 733   if      (i>4 || i<0) {
 734     AliError(Form("No parameter with index %d", i)) ;
 735   } else if (p==-1) {
 736     AliError(Form("No parameter with name %s", param.Data() )) ;
 737   } else
 738     par = (*fParameters)(p,i) ;
 739
 740   return par;
 741 }
 742
 743
 744 //DP____________________________________________________________________________
 745 //Float_t  AliPHOSPIDv1::GetDistance(AliPHOSEmcRecPoint * emc,AliPHOSCpvRecPoint * cpv, Option_t *  axis)const
 746 //{
 747 //  // Calculates the distance between the EMC RecPoint and the PPSD RecPoint
 748 //
 749 //  AliPHOSGeometry * geom =  AliPHOSGeometry::GetInstance();
 750 //  TVector3 vecEmc ;
 751 //  TVector3 vecCpv ;
 752 //  if(cpv){
 753 //    emc->GetLocalPosition(vecEmc) ;
 754 //    cpv->GetLocalPosition(vecCpv) ;
 755 //
 756 //    if(emc->GetPHOSMod() == cpv->GetPHOSMod()){
 757 //      // Correct to difference in CPV and EMC position due to different distance to center.
 758 //      // we assume, that particle moves from center
 759 //      Float_t dCPV = geom->GetIPtoOuterCoverDistance();
 760 //      Float_t dEMC = geom->GetIPtoCrystalSurface() ;
 761 //      dEMC         = dEMC / dCPV ;
 762 //      vecCpv = dEMC * vecCpv  - vecEmc ;
 763 //      if (axis == "X") return vecCpv.X();
 764 //      if (axis == "Y") return vecCpv.Y();
 765 //      if (axis == "Z") return vecCpv.Z();
 766 //      if (axis == "R") return vecCpv.Mag();
 767 //    }
 768 //    return 100000000 ;
 769 //  }
 770 //  return 100000000 ;
 771 //}
 772 //____________________________________________________________________________
 773 Int_t  AliPHOSPIDv1::GetCPVBit(AliPHOSTrackSegment * ts, Int_t effPur, Float_t e) const
 774 {
 775   //Calculates the pid bit for the CPV selection per each purity.
 776   if(effPur>2 || effPur<0)
 777     AliError(Form("Invalid Efficiency-Purity choice %d",effPur));
 778
 779 //DP  if(ts->GetCpvIndex()<0)
 780 //DP    return 1 ; //no CPV cluster
 781
 782   Float_t sigX = GetCpv2EmcDistanceCut("X",e);
 783   Float_t sigZ = GetCpv2EmcDistanceCut("Z",e);
 784
 785   Float_t deltaX = TMath::Abs(ts->GetCpvDistance("X"));
 786   Float_t deltaZ = TMath::Abs(ts->GetCpvDistance("Z"));
 787 //  Info("GetCPVBit"," xdist %f, sigx %f, zdist %f, sigz %f",deltaX, sigX, deltaZ,sigZ) ;
 788
 789   //if(deltaX>sigX*(effPur+1))
 790   //if((deltaX>sigX*(effPur+1)) || (deltaZ>sigZ*(effPur+1)))
 791   if((deltaX>sigX*(effPur+1)) && (deltaZ>sigZ*(effPur+1)))
 792     return 1;//Neutral
 793   else
 794     return 0;//Charged
 795 }
 796
 797 //____________________________________________________________________________
 798 Int_t  AliPHOSPIDv1::GetPrincipalBit(TString particle, const Double_t* p, Int_t effPur, Float_t e)const
 799 {
 800   //Is the particle inside de PCA ellipse?
 801
 802   particle.ToLower();
 803   Int_t    prinbit  = 0 ;
 804   Float_t a  = GetEllipseParameter(particle,"a" , e);
 805   Float_t b  = GetEllipseParameter(particle,"b" , e);
 806   Float_t c  = GetEllipseParameter(particle,"c" , e);
 807   Float_t x0 = GetEllipseParameter(particle,"x0", e);
 808   Float_t y0 = GetEllipseParameter(particle,"y0", e);
 809
 810   Float_t r = TMath::Power((p[0] - x0)/a,2) +
 811               TMath::Power((p[1] - y0)/b,2) +
 812             c*(p[0] -  x0)*(p[1] - y0)/(a*b) ;
 813   //3 different ellipses defined
 814   if((effPur==2) && (r<1./2.)) prinbit= 1;
 815   if((effPur==1) && (r<2.   )) prinbit= 1;
 816   if((effPur==0) && (r<9./2.)) prinbit= 1;
 817
 818   if(r<0)
 819     AliError("Negative square?") ;
 820
 821   return prinbit;
 822
 823 }
 824 //____________________________________________________________________________
 825 Int_t  AliPHOSPIDv1::GetHardPhotonBit(AliPHOSEmcRecPoint * emc) const
 826 {
 827   // Set bit for identified hard photons (E > 30 GeV)
 828   // if the second moment M2x is below the boundary
 829
 830   Float_t e   = emc->GetEnergy();
 831   if (e < 30.0) return 0;
 832   Float_t m2x = emc->GetM2x();
 833   Float_t m2xBoundary = GetParameterPhotonBoundary(0) *
 834     TMath::Exp(-TMath::Power(e-GetParameterPhotonBoundary(1),2)/2.0/
 835                 TMath::Power(GetParameterPhotonBoundary(2),2)) +
 836     GetParameterPhotonBoundary(3);
 837   AliDebug(1, Form("GetHardPhotonBit","E=%f, m2x=%f, boundary=%f",
 838                        e,m2x,m2xBoundary));
 839   if (m2x < m2xBoundary)
 840     return 1;// A hard photon
 841   else
 842     return 0;// Not a hard photon
 843 }
 844
 845 //____________________________________________________________________________
 846 Int_t  AliPHOSPIDv1::GetHardPi0Bit(AliPHOSEmcRecPoint * emc) const
 847 {
 848   // Set bit for identified hard pi0  (E > 30 GeV)
 849   // if the second moment M2x is above the boundary
 850
 851   Float_t e   = emc->GetEnergy();
 852   if (e < 30.0) return 0;
 853   Float_t m2x = emc->GetM2x();
 854   Float_t m2xBoundary = GetParameterPi0Boundary(0) +
 855                     e * GetParameterPi0Boundary(1);
 856   AliDebug(1,Form("E=%f, m2x=%f, boundary=%f",e,m2x,m2xBoundary));
 857   if (m2x > m2xBoundary)
 858     return 1;// A hard pi0
 859   else
 860     return 0;// Not a hard pi0
 861 }
 862
 863 //____________________________________________________________________________
 864 TVector3 AliPHOSPIDv1::GetMomentumDirection(AliPHOSEmcRecPoint * emc, AliPHOSCpvRecPoint * )const
 865 {
 866   // Calculates the momentum direction:
 867   //   1. if only a EMC RecPoint, direction is given by IP and this RecPoint
 868   //   2. if a EMC RecPoint and CPV RecPoint, direction is given by the line through the 2 recpoints
 869   //  However because of the poor position resolution of PPSD the direction is always taken as if we were
 870   //  in case 1.
 871
 872   TVector3 local ;
 873   emc->GetLocalPosition(local) ;
 874
 875   AliPHOSGeometry * phosgeom = AliPHOSGeometry::GetInstance() ;
 876   //Correct for the non-perpendicular incidence
 877   // Correction for the depth of the shower starting point (TDR p 127)
 878   Float_t para = 0.925 ;
 879   Float_t parb = 6.52 ;
 880
 881   //Remove Old correction (vertex at 0,0,0)
 882   TVector3 vtxOld(0.,0.,0.) ;
 883   TVector3 vInc ;
 884   Float_t x=local.X() ;
 885   Float_t z=local.Z() ;
 886   phosgeom->GetIncidentVector(vtxOld,emc->GetPHOSMod(),x,z,vInc) ;
 887   Float_t depthxOld = 0.;
 888   Float_t depthzOld = 0.;
 889   Float_t energy = emc->GetEnergy() ;
 890   if (energy > 0 && vInc.Y()!=0.) {
 891     depthxOld = ( para * TMath::Log(energy) + parb ) * vInc.X()/TMath::Abs(vInc.Y()) ;
 892     depthzOld = ( para * TMath::Log(energy) + parb ) * vInc.Z()/TMath::Abs(vInc.Y()) ;
 893   }
 894   else{
 895     AliError("Cluster with zero energy \n");
 896   }
 897   //Apply Real vertex
 898   phosgeom->GetIncidentVector(fVtx,emc->GetPHOSMod(),x,z,vInc) ;
 899   Float_t depthx = 0.;
 900   Float_t depthz = 0.;
 901   if (energy > 0 && vInc.Y()!=0.) {
 902     depthx = ( para * TMath::Log(energy) + parb ) * vInc.X()/TMath::Abs(vInc.Y()) ;
 903     depthz = ( para * TMath::Log(energy) + parb ) * vInc.Z()/TMath::Abs(vInc.Y()) ;
 904   }
 905
 906   //Correct for the vertex position and shower depth
 907   Double_t xd=x+(depthxOld-depthx) ;
 908   Double_t zd=z+(depthzOld-depthz) ;
 909   TVector3 dir(0,0,0) ;
 910   phosgeom->Local2Global(emc->GetPHOSMod(),xd,zd,dir) ;
 911
 912   dir-=fVtx ;
 913   dir.SetMag(1.) ;
 914
 915   return dir ;
 916 }
 917
 918 //________________________________________________________________________
 919 Double_t  AliPHOSPIDv1::LandauF(Double_t  x, Double_t y, Double_t * par)
 920 {
 921   //Given the energy x and the parameter y (tof, shower dispersion or cpv-emc distance),
 922   //this method returns a density probability of this parameter, given by a landau
 923   //function whose parameters depend with the energy  with a function: a/(x*x)+b/x+b
 924
 925   if (x > par[9]) x = par[9];
 926
 927   //Double_t cnt    = par[1] / (x*x) + par[2] / x + par[0] ;
 928   Double_t cnt    = par[0] + par[1] * x + par[2] * x * x ;
 929   Double_t mean   = par[4] / (x*x) + par[5] / x + par[3] ;
 930   Double_t sigma  = par[7] / (x*x) + par[8] / x + par[6] ;
 931
 932   if(TMath::Abs(sigma) > 1.e-10){
 933     return cnt*TMath::Landau(y,mean,sigma);
 934   }
 935   else
 936     return 0.;
 937
 938 }
 939 //________________________________________________________________________
 940 Double_t  AliPHOSPIDv1::LandauPol2(Double_t  x, Double_t y, Double_t * par)
 941 {
 942
 943   //Given the energy x and the parameter y (tof, shower dispersion or cpv-emc distance),
 944   //this method returns a density probability of this parameter, given by a landau
 945   //function whose parameters depend with the energy like second order polinomial
 946
 947   Double_t cnt    = par[2] * (x*x) + par[1] * x + par[0] ;
 948   Double_t mean   = par[5] * (x*x) + par[4] * x + par[3] ;
 949   Double_t sigma  = par[8] * (x*x) + par[7] * x + par[6] ;
 950
 951    if(TMath::Abs(sigma) > 1.e-10){
 952     return cnt*TMath::Landau(y,mean,sigma);
 953   }
 954   else
 955     return 0.;
 956
 957
 958 }
 959 // //________________________________________________________________________
 960 // Double_t  AliPHOSPIDv1::ChargedHadronDistProb(Double_t  x, Double_t y, Double_t * parg, Double_t * parl)
 961 // {
 962 //   Double_t cnt   = 0.0 ;
 963 //   Double_t mean  = 0.0 ;
 964 //   Double_t sigma = 0.0 ;
 965 //   Double_t arg   = 0.0 ;
 966 //   if (y < parl[4] / (x*x) + parl[5] / x + parl[3]){
 967 //     cnt    = parg[1] / (x*x) + parg[2] / x + parg[0] ;
 968 //     mean   = parg[4] / (x*x) + parg[5] / x + parg[3] ;
 969 //     sigma  = parg[7] / (x*x) + parg[8] / x + parg[6] ;
 970 //     TF1 * f = new TF1("gaus","gaus",0.,100.);
 971 //     f->SetParameters(cnt,mean,sigma);
 972 //     arg  = f->Eval(y) ;
 973 //   }
 974 //   else{
 975 //     cnt    = parl[1] / (x*x) + parl[2] / x + parl[0] ;
 976 //     mean   = parl[4] / (x*x) + parl[5] / x + parl[3] ;
 977 //     sigma  = parl[7] / (x*x) + parl[8] / x + parl[6] ;
 978 //     TF1 * f = new TF1("landau","landau",0.,100.);
 979 //     f->SetParameters(cnt,mean,sigma);
 980 //     arg  = f->Eval(y) ;
 981 //   }
 982 //   //  Double_t mean   = par[3] + par[4] * x + par[5] * x * x ;
 983 //   //   Double_t sigma  = par[6] + par[7] * x + par[8] * x * x ;
 984
 985 //   //Double_t arg    = -(y-mean)*(y-mean)/(2*sigma*sigma) ;
 986 //   //return cnt * TMath::Exp(arg) ;
 987
 988 //   return arg;
 989
 990 // }
 991 //____________________________________________________________________________
 992 void  AliPHOSPIDv1::MakePID()
 993 {
 994   // construct the PID weight from a Bayesian Method
 995
 996   const Int_t kSPECIES = AliPID::kSPECIESN ;
 997
 998   Int_t nparticles = fRecParticles->GetEntriesFast() ;
 999
1000   if ( !fEMCRecPoints || !fCPVRecPoints || !fTrackSegments ) {
1001     AliFatal("RecPoints or TrackSegments not found !") ;
1002   }
1003
1004   TIter next(fTrackSegments) ;
1005   AliPHOSTrackSegment * ts ;
1006   Int_t index = 0 ;
1007
1008   Double_t * stof[kSPECIES] ;
1009   Double_t * sdp [kSPECIES]  ;
1010   Double_t * scpv[kSPECIES] ;
1011   Double_t * sw  [kSPECIES] ;
1012   //Info("MakePID","Begin MakePID");
1013
1014   for (Int_t i =0; i< kSPECIES; i++){
1015     stof[i] = new Double_t[nparticles] ;
1016     sdp [i] = new Double_t[nparticles] ;
1017     scpv[i] = new Double_t[nparticles] ;
1018     sw  [i] = new Double_t[nparticles] ;
1019   }
1020
1021
1022   while ( (ts = (AliPHOSTrackSegment *)next()) ) {
1023
1024     //cout<<">>>>>> Bayesian Index "<<index<<endl;
1025
1026     AliPHOSEmcRecPoint * emc = 0 ;
1027     if(ts->GetEmcIndex()>=0)
1028       emc = (AliPHOSEmcRecPoint *) fEMCRecPoints->At(ts->GetEmcIndex()) ;
1029
1030 //    AliPHOSCpvRecPoint * cpv = 0 ;
1031 //    if(ts->GetCpvIndex()>=0)
1032 //      cpv = (AliPHOSCpvRecPoint *) cpvRecPoints->At(ts->GetCpvIndex()) ;
1033 //
1034 ////     Int_t track = 0 ;
1035 ////     track = ts->GetTrackIndex() ; //TPC tracks ?
1036
1037     if (!emc) {
1038       AliFatal(Form("-> emc(%d) = %d", ts->GetEmcIndex(), emc )) ;
1039     }
1040
1041
1042     // ############Tof#############################
1043
1044     //    Info("MakePID", "TOF");
1045     Float_t  en   = emc->GetEnergy();
1046     Double_t time = emc->GetTime() ;
1047     //    cout<<">>>>>>>Energy "<<en<<"Time "<<time<<endl;
1048
1049     // now get the signals probability
1050     // s(pid) in the Bayesian formulation
1051
1052     stof[AliPID::kPhoton][index]   = 1.;
1053     stof[AliPID::kElectron][index] = 1.;
1054     stof[AliPID::kEleCon][index]   = 1.;
1055     //We assing the same prob to charged hadrons, sum is 1
1056     stof[AliPID::kPion][index]     = 1./3.;
1057     stof[AliPID::kKaon][index]     = 1./3.;
1058     stof[AliPID::kProton][index]   = 1./3.;
1059     //We assing the same prob to neutral hadrons, sum is 1
1060     stof[AliPID::kNeutron][index]  = 1./2.;
1061     stof[AliPID::kKaon0][index]    = 1./2.;
1062     stof[AliPID::kMuon][index]     = 1.;
1063
1064     if(en <  fTOFEnThreshold) {
1065
1066       Double_t pTofPion = fTFpiong ->Eval(time) ; //gaus distribution
1067       Double_t pTofKaon = 0;
1068
1069       if(time < fTkaonl[1])
1070         pTofKaon = fTFkaong  ->Eval(time) ; //gaus distribution
1071       else
1072         pTofKaon = fTFkaonl  ->Eval(time) ; //landau distribution
1073
1074       Double_t pTofNucleon = 0;
1075
1076       if(time < fThhadronl[1])
1077         pTofNucleon = fTFhhadrong   ->Eval(time) ; //gaus distribution
1078       else
1079         pTofNucleon = fTFhhadronl   ->Eval(time) ; //landau distribution
1080       //We assing the same prob to neutral hadrons, sum is the average prob
1081       Double_t pTofNeHadron =  (pTofKaon + pTofNucleon)/2. ;
1082       //We assing the same prob to charged hadrons, sum is the average prob
1083       Double_t pTofChHadron =  (pTofPion + pTofKaon + pTofNucleon)/3. ;
1084
1085       stof[AliPID::kPhoton][index]   = fTFphoton     ->Eval(time) ;
1086       //gaus distribution
1087       stof[AliPID::kEleCon][index]   = stof[AliPID::kPhoton][index] ;
1088       //a conversion electron has the photon ToF
1089       stof[AliPID::kMuon][index]     = stof[AliPID::kPhoton][index] ;
1090
1091       stof[AliPID::kElectron][index] = pTofPion  ;
1092
1093       stof[AliPID::kPion][index]     =  pTofChHadron ;
1094       stof[AliPID::kKaon][index]     =  pTofChHadron ;
1095       stof[AliPID::kProton][index]   =  pTofChHadron ;
1096
1097       stof[AliPID::kKaon0][index]    =  pTofNeHadron ;
1098       stof[AliPID::kNeutron][index]  =  pTofNeHadron ;
1099     }
1100
1101     //    Info("MakePID", "Dispersion");
1102
1103     // ###########Shower shape: Dispersion####################
1104     Float_t dispersion = emc->GetDispersion();
1105     //DP: Correct for non-perpendicular incidence
1106     //DP: still to be done
1107
1108     //dispersion is not well defined if the cluster is only in few crystals
1109
1110     sdp[AliPID::kPhoton][index]   = 1. ;
1111     sdp[AliPID::kElectron][index] = 1. ;
1112     sdp[AliPID::kPion][index]     = 1. ;
1113     sdp[AliPID::kKaon][index]     = 1. ;
1114     sdp[AliPID::kProton][index]   = 1. ;
1115     sdp[AliPID::kNeutron][index]  = 1. ;
1116     sdp[AliPID::kEleCon][index]   = 1. ;
1117     sdp[AliPID::kKaon0][index]    = 1. ;
1118     sdp[AliPID::kMuon][index]     = 1. ;
1119
1120     if(en > fDispEnThreshold && emc->GetMultiplicity() >  fDispMultThreshold){
1121       sdp[AliPID::kPhoton][index]   = GausF(en , dispersion, fDphoton) ;
1122       sdp[AliPID::kElectron][index] = sdp[AliPID::kPhoton][index] ;
1123       sdp[AliPID::kPion][index]     = LandauF(en , dispersion, fDhadron ) ;
1124       sdp[AliPID::kKaon][index]     = sdp[AliPID::kPion][index]  ;
1125       sdp[AliPID::kProton][index]   = sdp[AliPID::kPion][index]  ;
1126       sdp[AliPID::kNeutron][index]  = sdp[AliPID::kPion][index]  ;
1127       sdp[AliPID::kEleCon][index]   = sdp[AliPID::kPhoton][index];
1128       sdp[AliPID::kKaon0][index]    = sdp[AliPID::kPion][index]  ;
1129       sdp[AliPID::kMuon][index]     = fDFmuon ->Eval(dispersion) ;
1130       //landau distribution
1131     }
1132
1133 //      Info("MakePID","multiplicity %d, dispersion %f", emc->GetMultiplicity(), dispersion);
1134 //      Info("MakePID","ss: photon %f, hadron %f ",  sdp[AliPID::kPhoton][index],  sdp[AliPID::kPion][index]);
1135 //       cout<<">>>>>multiplicity "<<emc->GetMultiplicity()<<", dispersion "<< dispersion<<endl ;
1136 //       cout<<"<<<<<ss: photon   "<<sdp[AliPID::kPhoton][index]<<", hadron    "<<sdp[AliPID::kPion][index]<<endl;
1137
1138     //########## CPV-EMC  Distance#######################
1139     //     Info("MakePID", "Distance");
1140
1141     Float_t x             = TMath::Abs(ts->GetCpvDistance("X")) ;
1142     Float_t z             = ts->GetCpvDistance("Z") ;
1143
1144     Double_t pcpv         = 0 ;
1145     Double_t pcpvneutral  = 0. ;
1146
1147     Double_t elprobx      = GausF(en , x, fXelectron) ;
1148     Double_t elprobz      = GausF(en , z, fZelectron) ;
1149     Double_t chprobx      = GausF(en , x, fXcharged)  ;
1150     Double_t chprobz      = GausF(en , z, fZcharged)  ;
1151     Double_t pcpvelectron = elprobx * elprobz;
1152     Double_t pcpvcharged  = chprobx * chprobz;
1153
1154 //     cout<<">>>>energy "<<en<<endl;
1155 //     cout<<">>>>electron : x "<<x<<" xprob "<<elprobx<<" z "<<z<<" zprob "<<elprobz<<endl;
1156 //     cout<<">>>>hadron   : x "<<x<<" xprob "<<chprobx<<" z "<<z<<" zprob "<<chprobz<<endl;
1157 //     cout<<">>>>electron : px*pz "<<pcpvelectron <<" hadron: px*pz "<<pcpvcharged<<endl;
1158
1159     // Is neutral or charged?
1160     if(pcpvelectron >= pcpvcharged)
1161       pcpv = pcpvelectron ;
1162     else
1163       pcpv = pcpvcharged ;
1164
1165     if(pcpv < fChargedNeutralThreshold)
1166       {
1167         pcpvneutral  = 1. ;
1168         pcpvcharged  = 0. ;
1169         pcpvelectron = 0. ;
1170       }
1171     //    else
1172     //      cout<<">>>>>>>>>>>CHARGED>>>>>>>>>>>"<<endl;
1173
1174     scpv[AliPID::kPion][index]     =  pcpvcharged  ;
1175     scpv[AliPID::kKaon][index]     =  pcpvcharged  ;
1176     scpv[AliPID::kProton][index]   =  pcpvcharged  ;
1177
1178     scpv[AliPID::kMuon][index]     =  pcpvelectron ;
1179     scpv[AliPID::kElectron][index] =  pcpvelectron ;
1180     scpv[AliPID::kEleCon][index]   =  pcpvelectron ;
1181
1182     scpv[AliPID::kPhoton][index]   =  pcpvneutral  ;
1183     scpv[AliPID::kNeutron][index]  =  pcpvneutral  ;
1184     scpv[AliPID::kKaon0][index]    =  pcpvneutral  ;
1185
1186
1187     //   Info("MakePID", "CPV passed");
1188
1189     //############## Pi0 #############################
1190     stof[AliPID::kPi0][index]      = 0. ;
1191     scpv[AliPID::kPi0][index]      = 0. ;
1192     sdp [AliPID::kPi0][index]      = 0. ;
1193
1194     if(en > 30.){
1195       // pi0 are detected via decay photon
1196       stof[AliPID::kPi0][index]  =   stof[AliPID::kPhoton][index];
1197       scpv[AliPID::kPi0][index]  = pcpvneutral  ;
1198       if(emc->GetMultiplicity() >  fDispMultThreshold)
1199         sdp [AliPID::kPi0][index]  = GausF(en , dispersion, fDpi0) ;
1200         //sdp [AliPID::kPi0][index]  = GausPol2(en , dispersion, fDpi0) ;
1201 //       cout<<"E = "<<en<<" GeV; disp = "<<dispersion<<"; mult = "
1202 //        <<emc->GetMultiplicity()<<endl;
1203 //       cout<<"PDF: photon = "<<sdp [AliPID::kPhoton][index]<<"; pi0 = "
1204 //        <<sdp [AliPID::kPi0][index]<<endl;
1205     }
1206
1207
1208
1209
1210     //############## muon #############################
1211
1212     if(en > 0.5){
1213       //Muons deposit few energy
1214       scpv[AliPID::kMuon][index]     =  0 ;
1215       stof[AliPID::kMuon][index]     =  0 ;
1216       sdp [AliPID::kMuon][index]     =  0 ;
1217     }
1218
1219     //Weight to apply to hadrons due to energy reconstruction
1220
1221     Float_t weight = fERecWeight ->Eval(en) ;
1222
1223     sw[AliPID::kPhoton][index]   = 1. ;
1224     sw[AliPID::kElectron][index] = 1. ;
1225     sw[AliPID::kPion][index]     = weight ;
1226     sw[AliPID::kKaon][index]     = weight ;
1227     sw[AliPID::kProton][index]   = weight ;
1228     sw[AliPID::kNeutron][index]  = weight ;
1229     sw[AliPID::kEleCon][index]   = 1. ;
1230     sw[AliPID::kKaon0][index]    = weight ;
1231     sw[AliPID::kMuon][index]     = weight ;
1232     sw[AliPID::kPi0][index]      = 1. ;
1233
1234 //     if(en > 0.5){
1235 //       cout<<"######################################################"<<endl;
1236 //       //cout<<"MakePID: energy "<<en<<", tof "<<time<<", distance "<<distance<<", dispersion "<<dispersion<<endl ;
1237 //       cout<<"MakePID: energy "<<en<<", tof "<<time<<", dispersion "<<dispersion<<", x "<<x<<", z "<<z<<endl ;
1238 //       cout<<">>>>>multiplicity "<<emc->GetMultiplicity()<<endl;
1239 //       cout<<">>>>electron : xprob "<<elprobx<<" zprob "<<elprobz<<endl;
1240 //       cout<<">>>>hadron   : xprob "<<chprobx<<" zprob "<<chprobz<<endl;
1241 //       cout<<">>>>electron : px*pz "<<pcpvelectron <<" hadron: px*pz "<<pcpvcharged<<endl;
1242
1243 //        cout<<"Photon   , pid "<< fInitPID[AliPID::kPhoton]<<" tof "<<stof[AliPID::kPhoton][index]
1244 //        <<", cpv "<<scpv[AliPID::kPhoton][index]<<", ss "<<sdp[AliPID::kPhoton][index]<<endl;
1245 //       cout<<"EleCon   , pid "<< fInitPID[AliPID::kEleCon]<<", tof "<<stof[AliPID::kEleCon][index]
1246 //        <<", cpv "<<scpv[AliPID::kEleCon][index]<<" ss "<<sdp[AliPID::kEleCon][index]<<endl;
1247 //       cout<<"Electron , pid "<< fInitPID[AliPID::kElectron]<<", tof "<<stof[AliPID::kElectron][index]
1248 //        <<", cpv "<<scpv[AliPID::kElectron][index]<<" ss "<<sdp[AliPID::kElectron][index]<<endl;
1249 //       cout<<"Muon     , pid "<< fInitPID[AliPID::kMuon]<<", tof "<<stof[AliPID::kMuon][index]
1250 //        <<", cpv "<<scpv[AliPID::kMuon][index]<<" ss "<<sdp[AliPID::kMuon][index]<<endl;
1251 //        cout<<"Pi0      , pid "<< fInitPID[AliPID::kPi0]<<", tof "<<stof[AliPID::kPi0][index]
1252 //        <<", cpv "<<scpv[AliPID::kPi0][index]<<" ss "<<sdp[AliPID::kPi0][index]<<endl;
1253 //       cout<<"Pion     , pid "<< fInitPID[AliPID::kPion]<<", tof "<<stof[AliPID::kPion][index]
1254 //        <<", cpv "<<scpv[AliPID::kPion][index]<<" ss "<<sdp[AliPID::kPion][index]<<endl;
1255 //       cout<<"Kaon0    , pid "<< fInitPID[AliPID::kKaon0]<<", tof "<<stof[AliPID::kKaon0][index]
1256 //        <<", cpv "<<scpv[AliPID::kKaon0][index]<<" ss "<<sdp[AliPID::kKaon0][index]<<endl;
1257 //       cout<<"Kaon     , pid "<< fInitPID[AliPID::kKaon]<<", tof "<<stof[AliPID::kKaon][index]
1258 //        <<", cpv "<<scpv[AliPID::kKaon][index]<<" ss "<<sdp[AliPID::kKaon][index]<<endl;
1259 //       cout<<"Neutron  , pid "<< fInitPID[AliPID::kNeutron]<<", tof "<<stof[AliPID::kNeutron][index]
1260 //        <<", cpv "<<scpv[AliPID::kNeutron][index]<<" ss "<<sdp[AliPID::kNeutron][index]<<endl;
1261 //       cout<<"Proton   , pid "<< fInitPID[AliPID::kProton]<<", tof "<<stof[AliPID::kProton][index]
1262 //        <<", cpv "<<scpv[AliPID::kProton][index]<<" ss "<<sdp[AliPID::kProton][index]<<endl;
1263 //       cout<<"######################################################"<<endl;
1264 //     }
1265       index++;
1266   }
1267
1268   //for (index = 0 ; index < kSPECIES ; index++)
1269   // pid[index] /= nparticles ;
1270
1271
1272   //  Info("MakePID", "Total Probability calculation");
1273
1274   for(index = 0 ; index < nparticles ; index ++) {
1275
1276     AliPHOSRecParticle * recpar = static_cast<AliPHOSRecParticle *>(fRecParticles->At(index));
1277
1278     //Conversion electron?
1279
1280     if(recpar->IsEleCon()){
1281       fInitPID[AliPID::kEleCon]   = 1. ;
1282       fInitPID[AliPID::kPhoton]   = 0. ;
1283       fInitPID[AliPID::kElectron] = 0. ;
1284     }
1285     else{
1286       fInitPID[AliPID::kEleCon]   = 0. ;
1287       fInitPID[AliPID::kPhoton]   = 1. ;
1288       fInitPID[AliPID::kElectron] = 1. ;
1289     }
1290     //  fInitPID[AliPID::kEleCon]   = 0. ;
1291
1292
1293     // calculates the Bayesian weight
1294
1295     Int_t jndex ;
1296     Double_t wn = 0.0 ;
1297     for (jndex = 0 ; jndex < kSPECIES ; jndex++)
1298       wn += stof[jndex][index] * sdp[jndex][index]  * scpv[jndex][index] *
1299         sw[jndex][index] * fInitPID[jndex] ;
1300
1301     //    cout<<"*************wn "<<wn<<endl;
1302     if (TMath::Abs(wn)>0)
1303       for (jndex = 0 ; jndex < kSPECIES ; jndex++) {
1304         //cout<<"jndex "<<jndex<<" wn "<<wn<<" SetPID * wn"
1305         //<<stof[jndex][index] * sdp[jndex][index] * pid[jndex]  << endl;
1306         //cout<<" tof "<<stof[jndex][index] << " disp " <<sdp[jndex][index] << " pid "<< fInitPID[jndex] << endl;
1307         //      if(jndex ==  AliPID::kPi0 || jndex ==  AliPID::kPhoton){
1308         //        cout<<"Particle "<<jndex<<"  final prob * wn   "
1309         //            <<stof[jndex][index] * sdp[jndex][index] * scpv[jndex][index] *
1310         //          fInitPID[jndex] <<"  wn  "<< wn<<endl;
1311         //        cout<<"pid "<< fInitPID[jndex]<<", tof "<<stof[jndex][index]
1312         //            <<", cpv "<<scpv[jndex][index]<<" ss "<<sdp[jndex][index]<<endl;
1313         //      }
1314         recpar->SetPID(jndex, stof[jndex][index] * sdp[jndex][index] *
1315                        sw[jndex][index] * scpv[jndex][index] *
1316                        fInitPID[jndex] / wn) ;
1317       }
1318   }
1319   //  Info("MakePID", "Delete");
1320
1321   for (Int_t i =0; i< kSPECIES; i++){
1322     delete [] stof[i];
1323     delete [] sdp [i];
1324     delete [] scpv[i];
1325     delete [] sw  [i];
1326   }
1327   //  Info("MakePID","End MakePID");
1328 }
1329
1330 //____________________________________________________________________________
1331 void  AliPHOSPIDv1::MakeRecParticles()
1332 {
1333   // Makes a RecParticle out of a TrackSegment
1334
1335   if ( !fEMCRecPoints || !fCPVRecPoints || !fTrackSegments ) {
1336     AliFatal("RecPoints or TrackSegments not found !") ;
1337   }
1338   fRecParticles->Clear();
1339
1340   TIter next(fTrackSegments) ;
1341   AliPHOSTrackSegment * ts ;
1342   Int_t index = 0 ;
1343   AliPHOSRecParticle * rp ;
1344   while ( (ts = (AliPHOSTrackSegment *)next()) ) {
1345     //  cout<<">>>>>>>>>>>>>>>PCA Index "<<index<<endl;
1346     new( (*fRecParticles)[index] ) AliPHOSRecParticle() ;
1347     rp = (AliPHOSRecParticle *)fRecParticles->At(index) ;
1348     rp->SetTrackSegment(index) ;
1349     rp->SetIndexInList(index) ;
1350
1351     AliPHOSEmcRecPoint * emc = 0 ;
1352     if(ts->GetEmcIndex()>=0)
1353       emc = (AliPHOSEmcRecPoint *) fEMCRecPoints->At(ts->GetEmcIndex()) ;
1354
1355     AliPHOSCpvRecPoint * cpv = 0 ;
1356     if(ts->GetCpvIndex()>=0)
1357       cpv = (AliPHOSCpvRecPoint *) fCPVRecPoints->At(ts->GetCpvIndex()) ;
1358
1359     Int_t track = 0 ;
1360     track = ts->GetTrackIndex() ;
1361
1362     // Now set type (reconstructed) of the particle
1363
1364     // Choose the cluster energy range
1365
1366     if (!emc) {
1367       AliFatal(Form("-> emc(%d) = %d", ts->GetEmcIndex(), emc )) ;
1368     }
1369
1370     Float_t e = emc->GetEnergy() ;
1371
1372     Float_t  lambda[2] ;
1373     emc->GetElipsAxis(lambda) ;
1374
1375     if((lambda[0]>0.01) && (lambda[1]>0.01)){
1376       // Looking PCA. Define and calculate the data (X),
1377       // introduce in the function X2P that gives the components (P).
1378
1379       Float_t  spher = 0. ;
1380       Float_t  emaxdtotal = 0. ;
1381
1382       if((lambda[0]+lambda[1])!=0)
1383         spher=TMath::Abs(lambda[0]-lambda[1])/(lambda[0]+lambda[1]);
1384
1385       emaxdtotal=emc->GetMaximalEnergy()/emc->GetEnergy();
1386
1387       fX[0] = lambda[0] ;
1388       fX[1] = lambda[1] ;
1389       fX[2] = emc->GetDispersion() ;
1390       fX[3] = spher ;
1391       fX[4] = emc->GetMultiplicity() ;
1392       fX[5] = emaxdtotal ;
1393       fX[6] = emc->GetCoreEnergy() ;
1394
1395       fPrincipalPhoton->X2P(fX,fPPhoton);
1396       fPrincipalPi0   ->X2P(fX,fPPi0);
1397
1398     }
1399     else{
1400       fPPhoton[0]=-100.0;  //We do not accept clusters with
1401       fPPhoton[1]=-100.0;  //one cell as a photon-like
1402       fPPi0[0]   =-100.0;
1403       fPPi0[1]   =-100.0;
1404     }
1405
1406     Float_t time = emc->GetTime() ;
1407     rp->SetTof(time) ;
1408
1409     // Loop of Efficiency-Purity (the 3 points of purity or efficiency
1410     // are taken into account to set the particle identification)
1411     for(Int_t effPur = 0; effPur < 3 ; effPur++){
1412
1413       // Looking at the CPV detector. If RCPV greater than CpvEmcDistance,
1414       // 1st,2nd or 3rd bit (depending on the efficiency-purity point )
1415       // is set to 1
1416       if(GetCPVBit(ts, effPur,e) == 1 ){
1417         rp->SetPIDBit(effPur) ;
1418         //cout<<"CPV bit "<<effPur<<endl;
1419       }
1420       // Looking the TOF. If TOF smaller than gate,  4th, 5th or 6th
1421       // bit (depending on the efficiency-purity point )is set to 1
1422       if(time< (*fParameters)(3,effPur))
1423         rp->SetPIDBit(effPur+3) ;
1424
1425       //Photon PCA
1426       //If we are inside the ellipse, 7th, 8th or 9th
1427       // bit (depending on the efficiency-purity point )is set to 1
1428       if(GetPrincipalBit("photon",fPPhoton,effPur,e) == 1)
1429         rp->SetPIDBit(effPur+6) ;
1430
1431       //Pi0 PCA
1432       //If we are inside the ellipse, 10th, 11th or 12th
1433       // bit (depending on the efficiency-purity point )is set to 1
1434       if(GetPrincipalBit("pi0"   ,fPPi0   ,effPur,e) == 1)
1435         rp->SetPIDBit(effPur+9) ;
1436     }
1437     if(GetHardPhotonBit(emc))
1438       rp->SetPIDBit(12) ;
1439     if(GetHardPi0Bit   (emc))
1440       rp->SetPIDBit(13) ;
1441
1442     if(track >= 0)
1443       rp->SetPIDBit(14) ;
1444
1445     //Set momentum, energy and other parameters
1446     Float_t  encal = GetCalibratedEnergy(e);
1447     TVector3 dir   = GetMomentumDirection(emc,cpv) ;
1448     dir.SetMag(encal) ;
1449     rp->SetMomentum(dir.X(),dir.Y(),dir.Z(),encal) ;
1450     rp->SetCalcMass(0);
1451     rp->Name(); //If photon sets the particle pdg name to gamma
1452     rp->SetProductionVertex(fVtx.X(),fVtx.Y(),fVtx.Z(),0);
1453     rp->SetFirstMother(-1);
1454     rp->SetLastMother(-1);
1455     rp->SetFirstDaughter(-1);
1456     rp->SetLastDaughter(-1);
1457     rp->SetPolarisation(0,0,0);
1458     //Set the position in global coordinate system from the RecPoint
1459     AliPHOSTrackSegment * ts1  = static_cast<AliPHOSTrackSegment *>(fTrackSegments->At(rp->GetPHOSTSIndex()));
1460     AliPHOSEmcRecPoint  * erp = static_cast<AliPHOSEmcRecPoint *>(fEMCRecPoints->At(ts1->GetEmcIndex()));
1461     TVector3 pos ;
1462     fGeom->GetGlobalPHOS(erp, pos) ;
1463     rp->SetPos(pos);
1464     index++ ;
1465   }
1466 }
1467
1468 //____________________________________________________________________________
1469 void  AliPHOSPIDv1::Print(const Option_t *) const
1470 {
1471   // Print the parameters used for the particle type identification
1472
1473     AliInfo("=============== AliPHOSPIDv1 ================") ;
1474     printf("Making PID\n") ;
1475     printf("    Pricipal analysis file from 0.5 to 100 %s\n", fFileNamePrincipalPhoton.Data() )   ;
1476     printf("    Name of parameters file     %s\n", fFileNameParameters.Data() )  ;
1477     printf("    Matrix of Parameters: 14x4\n") ;
1478     printf("        Energy Calibration  1x3 [3 parametres to calibrate energy: A + B* E + C * E^2]\n") ;
1479     printf("        RCPV 2x3 rows x and z, columns function cut parameters\n") ;
1480     printf("        TOF  1x3 [High Eff-Low Pur,Medium Eff-Pur, Low Eff-High Pur]\n") ;
1481     printf("        PCA  5x4 [5 ellipse parametres and 4 parametres to calculate them: A/Sqrt(E) + B* E + C * E^2 + D]\n") ;
1482     printf("    Pi0 PCA  5x3 [5 ellipse parametres and 3 parametres to calculate them: A + B* E + C * E^2]\n") ;
1483     fParameters->Print() ;
1484 }
1485
1486
1487
1488 //____________________________________________________________________________
1489 void AliPHOSPIDv1::PrintRecParticles(Option_t * option)
1490 {
1491   // Print table of reconstructed particles
1492
1493   TString message ;
1494   message  = "       found " ;
1495   message += fRecParticles->GetEntriesFast();
1496   message += " RecParticles\n" ;
1497
1498   if(strstr(option,"all")) {  // printing found TS
1499     message += "\n  PARTICLE         Index    \n" ;
1500
1501     Int_t index ;
1502     for (index = 0 ; index < fRecParticles->GetEntries() ; index++) {
1503       AliPHOSRecParticle * rp = (AliPHOSRecParticle * ) fRecParticles->At(index) ;
1504       message += "\n" ;
1505       message += rp->Name().Data() ;
1506       message += " " ;
1507       message += rp->GetIndexInList() ;
1508       message += " " ;
1509       message += rp->GetType()  ;
1510     }
1511   }
1512   AliInfo(message.Data() ) ;
1513 }
1514
1515 //____________________________________________________________________________
1516 void  AliPHOSPIDv1::SetParameters()
1517 {
1518   // PCA : To do the Principal Components Analysis it is necessary
1519   // the Principal file, which is opened here
1520   fX       = new double[7]; // Data for the PCA
1521   fPPhoton = new double[7]; // Eigenvalues of the PCA
1522   fPPi0    = new double[7]; // Eigenvalues of the Pi0 PCA
1523
1524   // Read photon principals from the photon file
1525
1526   fFileNamePrincipalPhoton = "$ALICE_ROOT/PHOS/PCA8pa15_0.5-100.root" ;
1527   TFile f( fFileNamePrincipalPhoton.Data(), "read" ) ;
1528   fPrincipalPhoton = dynamic_cast<TPrincipal*> (f.Get("principal")) ;
1529   f.Close() ;
1530
1531   // Read pi0 principals from the pi0 file
1532
1533   fFileNamePrincipalPi0    = "$ALICE_ROOT/PHOS/PCA_pi0_40-120.root" ;
1534   TFile fPi0( fFileNamePrincipalPi0.Data(), "read" ) ;
1535   fPrincipalPi0    = dynamic_cast<TPrincipal*> (fPi0.Get("principal")) ;
1536   fPi0.Close() ;
1537
1538   // Open parameters file and initialization of the Parameters matrix.
1539   // In the File Parameters.dat are all the parameters. These are introduced
1540   // in a matrix of 16x4
1541   //
1542   // All the parameters defined in this file are, in order of row:
1543   // line   0   : calibration
1544   // lines  1,2 : CPV rectangular cat for X and Z
1545   // line   3   : TOF cut
1546   // lines  4-8 : parameters to calculate photon PCA ellipse
1547   // lines  9-13: parameters to calculate pi0 PCA ellipse
1548   // lines 14-15: parameters to calculate border for high-pt photons and pi0
1549
1550   fFileNameParameters = gSystem->ExpandPathName("$ALICE_ROOT/PHOS/Parameters.dat");
1551   fParameters = new TMatrixF(16,4) ;
1552   const Int_t kMaxLeng=255;
1553   char string[kMaxLeng];
1554
1555   // Open a text file with PID parameters
1556   FILE *fd = fopen(fFileNameParameters.Data(),"r");
1557   if (!fd)
1558     AliFatal(Form("File %s with a PID parameters cannot be opened\n",
1559           fFileNameParameters.Data()));
1560
1561   Int_t i=0;
1562   // Read parameter file line-by-line and skip empty line and comments
1563   while (fgets(string,kMaxLeng,fd) != NULL) {
1564     if (string[0] == '\n' ) continue;
1565     if (string[0] == '!'  ) continue;
1566     sscanf(string, "%f %f %f %f",
1567            &(*fParameters)(i,0), &(*fParameters)(i,1),
1568            &(*fParameters)(i,2), &(*fParameters)(i,3));
1569     i++;
1570     AliDebug(1, Form("SetParameters", "line %d: %s",i,string));
1571   }
1572   fclose(fd);
1573 }
1574
1575 //____________________________________________________________________________
1576 void  AliPHOSPIDv1::SetParameterCalibration(Int_t i,Float_t param)
1577 {
1578   // Set parameter "Calibration" i to a value param
1579   if(i>2 || i<0) {
1580     AliError(Form("Invalid parameter number: %d",i));
1581   } else
1582     (*fParameters)(0,i) = param ;
1583 }
1584
1585 //____________________________________________________________________________
1586 void  AliPHOSPIDv1::SetParameterCpv2Emc(Int_t i, TString axis, Float_t cut)
1587 {
1588   // Set the parameters to calculate Cpv-to-Emc Distance Cut depending on
1589   // Purity-Efficiency point i
1590
1591   if(i>2 || i<0) {
1592     AliError(Form("Invalid parameter number: %d",i));
1593   } else {
1594     axis.ToLower();
1595     if      (axis == "x") (*fParameters)(1,i) = cut;
1596     else if (axis == "z") (*fParameters)(2,i) = cut;
1597     else {
1598       AliError(Form("Invalid axis name: %s",axis.Data()));
1599     }
1600   }
1601 }
1602
1603 //____________________________________________________________________________
1604 void  AliPHOSPIDv1::SetParameterPhotonBoundary(Int_t i,Float_t param)
1605 {
1606   // Set parameter "Hard photon boundary" i to a value param
1607   if(i>4 || i<0) {
1608     AliError(Form("Invalid parameter number: %d",i));
1609   } else
1610     (*fParameters)(14,i) = param ;
1611 }
1612
1613 //____________________________________________________________________________
1614 void  AliPHOSPIDv1::SetParameterPi0Boundary(Int_t i,Float_t param)
1615 {
1616   // Set parameter "Hard pi0 boundary" i to a value param
1617   if(i>1 || i<0) {
1618     AliError(Form("Invalid parameter number: %d",i));
1619   } else
1620     (*fParameters)(15,i) = param ;
1621 }
1622
1623 //_____________________________________________________________________________
1624 void  AliPHOSPIDv1::SetParameterTimeGate(Int_t i, Float_t gate)
1625 {
1626   // Set the parameter TimeGate depending on Purity-Efficiency point i
1627   if (i>2 || i<0) {
1628     AliError(Form("Invalid Efficiency-Purity choice %d",i));
1629   } else
1630     (*fParameters)(3,i)= gate ;
1631 }
1632
1633 //_____________________________________________________________________________
1634 void  AliPHOSPIDv1::SetParameterToCalculateEllipse(TString particle, TString param, Int_t i, Float_t par)
1635 {
1636   // Set the parameter "i" that is needed to calculate the ellipse
1637   // parameter "param" for a particle "particle"
1638
1639   particle.ToLower();
1640   param.   ToLower();
1641   Int_t p= -1;
1642   Int_t offset=0;
1643
1644   if      (particle == "photon") offset=0;
1645   else if (particle == "pi0")    offset=5;
1646   else
1647     AliError(Form("Wrong particle name: %s (choose from pi0/photon)\n",
1648                   particle.Data()));
1649
1650   if     (param.Contains("a")) p=4+offset;
1651   else if(param.Contains("b")) p=5+offset;
1652   else if(param.Contains("c")) p=6+offset;
1653   else if(param.Contains("x0"))p=7+offset;
1654   else if(param.Contains("y0"))p=8+offset;
1655   if((i>4)||(i<0)) {
1656     AliError(Form("No parameter with index %d", i)) ;
1657   } else if(p==-1) {
1658     AliError(Form("No parameter with name %s", param.Data() )) ;
1659   } else
1660     (*fParameters)(p,i) = par ;
1661 }
1662
1663 //____________________________________________________________________________
1664 void AliPHOSPIDv1::GetVertex(void)
1665 { //extract vertex either using ESD or generator
1666
1667   //Try to extract vertex from data
1668   if(fESD){
1669     const AliESDVertex *esdVtx = fESD->GetVertex() ;
1670     if(esdVtx && esdVtx->GetChi2()!=0.){
1671       fVtx.SetXYZ(esdVtx->GetXv(),esdVtx->GetYv(),esdVtx->GetZv()) ;
1672       return ;
1673     }
1674   }
1675
1676   // Use vertex diamond from CDB GRP folder if the one from ESD is missing
1677   // PLEASE FIX IT
1678   AliWarning("Can not read vertex from data, use fixed \n") ;
1679   fVtx.SetXYZ(0.,0.,0.) ;
1680
1681 }
1682 //_______________________________________________________________________
1683 void AliPHOSPIDv1::SetInitPID(const Double_t *p) {
1684   // Sets values for the initial population of each particle type
1685   for (Int_t i=0; i<AliPID::kSPECIESN; i++) fInitPID[i] = p[i];
1686 }
1687 //_______________________________________________________________________
1688 void AliPHOSPIDv1::GetInitPID(Double_t *p) const {
1689   // Gets values for the initial population of each particle type
1690   for (Int_t i=0; i<AliPID::kSPECIESN; i++) p[i] = fInitPID[i];
1691 }