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