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