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