EMCAL/AliEMCALPID.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 /* History of cvs commits:
  18  *
  19  * $Log$
  20  * Revision 1.14  2007/07/26 16:54:53  morsch
  21  * Changes in AliESDEvent fwd declarartions.
  22  *
  23  * Revision 1.13  2007/07/11 13:43:29  hristov
  24  * New class AliESDEvent, backward compatibility with the old AliESD (Christian)
  25  *
  26  * Revision 1.12  2007/06/11 20:43:06  hristov
  27  * Changes required by the updated AliESDCaloCluster (Gustavo)
  28  *
  29  * Revision 1.11  2007/03/30 13:50:34  gustavo
  30  * PID for particles with E < 5 GeV was not done, temporal solution found (Guenole)
  31  *
  32  * Revision 1.10  2007/03/09 14:34:11  gustavo
  33  * Correct probability calculation, added missing initialization of data members
  34  *
  35  * Revision 1.9  2007/02/20 20:17:43  hristov
  36  * Corrected array size, removed warnings (icc)
  37  *
  38  * Revision 1.8  2006/12/19 08:49:35  gustavo
  39  * New PID class for EMCAL, bayesian analysis done with ESD data, PID information filled when calling AliEMCALPID in AliEMCALReconstructor::FillESD()
  40  *
  41  *
  42  */
  43 //    to compute PID for all the clusters in ESDs.root file
  44 //     the ESDs.root have to be in the same directory as the class
  45 //
  46 //
  47 //
  48 //
  49 //
  50 //      AliEMCALPID::CalculPID(Energy,Lambda0)
  51 //        Calcul PID for all clusters in AliESDs.root file
  52 //        keep this function for the moment for a simple verification, could be removed
  53 //
  54 //
  55 //
  56 //   AliEMCALPID::CalculPID(Energy,Lambda0)
  57 //    calcul PID Weght for a cluster with Energy, Lambda0 .
  58 //    Double_t PIDFinal[AliPID::kSPECIESN]  is the standard PID for :
  59 //
  60 //
  61 //
  62 //      kElectron :  fPIDFinal[0]
  63 //      kMuon     :  fPIDFinal[1]
  64 //      kPion     :  fPIDFinal[2]
  65 //      kKaon     :  fPIDFinal[3]
  66 //      kProton   :  fPIDFinal[4]
  67 //      kPhoton   :  fPIDFinal[5]
  68 //      kPi0      :  fPIDFinal[6]
  69 //      kNeutron  :  fPIDFinal[7]
  70 //      kKaon0    :  fPIDFinal[8]
  71 //      kEleCon   :  fPIDFinal[9]
  72 //      kUnknown  :  fPIDFinal[10]
  73 //
  74 //
  75 //    PID[3] is a simple PID for
  76 //      Electron & Photon  PID[0]
  77 //                    Pi0  PID[1]
  78 //                 Hadron  PID[2]
  79 //
  80 //
  81 //
  82 //
  83 //
  84 // --- ROOT system ---
  85
  86 // standard C++ includes
  87 #include <Riostream.h>
  88 // #include <cstdlib>
  89 // #include <cmath>
  90
  91 // ROOT includes
  92 #include "TTree.h"
  93 #include "TStyle.h"
  94 #include "TVector3.h"
  95 #include "TBranch.h"
  96 #include "TClonesArray.h"
  97 #include "TCanvas.h"
  98 #include "TLorentzVector.h"
  99 #include "TMath.h"
 100 #include "TFile.h"
 101 #include "TH1.h"
 102 #include "TH2.h"
 103 #include "TParticle.h"
 104
 105 // // STEER includes
 106 // #include "AliRun.h"
 107 // #include "AliRunLoader.h"
 108 // #include "AliHeader.h"
 109 // #include "AliLoader.h"
 110 // #include "AliStack.h"
 111 // #include "AliESDtrack.h"
 112 // #include "AliESDEvent.h"
 113 #include "AliLog.h"
 114 #include "AliEMCALPID.h"
 115 #include "AliESDCaloCluster.h"
 116
 117 ClassImp(AliEMCALPID)
 118
 119 //______________________________________________
 120   AliEMCALPID::AliEMCALPID():
 121     fPrintInfo(kFALSE), fProbGamma(0.),fProbPiZero(0.),fProbHadron(0.),fReconstructor(kFALSE)
 122 {
 123   //
 124   // Constructor.
 125   // Initialize all constant values which have to be used
 126   // during PID algorithm execution
 127   //
 128
 129   // set flag for printing to FALSE by default
 130   fPrintInfo = kFALSE;
 131
 132   // as a first step, all array elements are initialized to 0.0
 133   Int_t i, j;
 134   for (i = 0; i < 6; i++) {
 135     for (j = 0; j < 6; j++) {
 136       fGamma[i][j] = fHadron[i][j] = fPiZero5to10[i][j] = fPiZero10to60[i][j] = 0.;
 137     }
 138   }
 139
 140   // then, only the ones which must be not zero are initialized
 141   // while the others will remain to the value 0.0
 142
 143   fGamma[0][0] =  0.038022;
 144   fGamma[0][1] = -0.0001883;
 145   fGamma[0][2] =  5.449e-06;
 146
 147   fGamma[1][0] =  0.207313;
 148   fGamma[1][1] = -0.000978;
 149   fGamma[1][2] =  0.00001634;
 150
 151   fGamma[2][0] =  0.043364;
 152   fGamma[2][1] = -0.0002048;
 153   fGamma[2][2] =  8.661e-06;
 154   fGamma[2][3] = -1.353e-07;
 155
 156   fGamma[3][0] =  0.265004;
 157   fGamma[3][1] =  0.061298;
 158   fGamma[3][2] = -0.003203;
 159   fGamma[3][3] =  4.73e-05;
 160
 161   fGamma[4][0] =  0.243579;
 162   fGamma[4][1] = -1.614e-05;
 163
 164   fGamma[5][0] =  0.002942;
 165   fGamma[5][1] = -3.976e-05;
 166
 167   fHadron[0][0] =  0.011945 / 3.;
 168   fHadron[0][1] =  0.000386 / 3.;
 169   fHadron[0][2] = -0.000014 / 3.;
 170   fHadron[0][3] =  1.336e-07 / 3.;
 171
 172   fHadron[1][0] =  0.496544;
 173   fHadron[1][1] = -0.003226;
 174   fHadron[1][2] =  0.00001678;
 175
 176   fHadron[2][0] =  0.144838;
 177   fHadron[2][1] = -0.002954;
 178   fHadron[2][2] =  0.00008754;
 179   fHadron[2][3] = -7.587e-07;
 180
 181   fHadron[3][0] =  1.264461 / 7.;
 182   fHadron[3][1] =  0.002097 / 7.;
 183
 184   fHadron[4][0] =  0.261950;
 185   fHadron[4][1] = -0.001078;
 186   fHadron[4][2] =  0.00003237;
 187   fHadron[4][3] = -3.241e-07;
 188   fHadron[4][4] =  0.;
 189   fHadron[4][5] =  0.;
 190   fHadron[5][0] =  0.010317;
 191   fHadron[5][1] =  0.;
 192   fHadron[5][2] =  0.;
 193   fHadron[5][3] =  0.;
 194   fHadron[5][4] =  0.;
 195   fHadron[5][5] =  0.;
 196
 197   fPiZero5to10[0][0] = 0.009138;
 198   fPiZero5to10[0][1] = 0.0006377;
 199
 200   fPiZero5to10[1][0] = 0.08;
 201
 202   fPiZero5to10[2][0] = -0.061119;
 203   fPiZero5to10[2][1] =  0.019013;
 204
 205   fPiZero5to10[3][0] =  0.2;
 206
 207   fPiZero5to10[4][0] =  0.252044;
 208   fPiZero5to10[4][1] = -0.002315;
 209
 210   fPiZero5to10[5][0] =  0.002942;
 211   fPiZero5to10[5][1] = -3.976e-05;
 212
 213   fPiZero10to60[0][0] =  0.009138;
 214   fPiZero10to60[0][1] =  0.0006377;
 215
 216   fPiZero10to60[1][0] =  1.272837;
 217   fPiZero10to60[1][1] = -0.069708;
 218   fPiZero10to60[1][2] =  0.001568;
 219   fPiZero10to60[1][3] = -1.162e-05;
 220
 221   fPiZero10to60[2][0] =  0.139703;
 222   fPiZero10to60[2][1] =  0.003687;
 223   fPiZero10to60[2][2] = -0.000568;
 224   fPiZero10to60[2][3] =  1.498e-05;
 225   fPiZero10to60[2][4] = -1.174e-07;
 226
 227   fPiZero10to60[3][0] = -0.826367;
 228   fPiZero10to60[3][1] =  0.096951;
 229   fPiZero10to60[3][2] = -0.002215;
 230   fPiZero10to60[3][3] =  2.523e-05;
 231
 232   fPiZero10to60[4][0] =  0.249890;
 233   fPiZero10to60[4][1] = -0.000063;
 234
 235   fPiZero10to60[5][0] =  0.002942;
 236   fPiZero10to60[5][1] = -3.976e-05;
 237
 238   fPIDWeight[0] = -1;
 239   fPIDWeight[1] = -1;
 240   fPIDWeight[2] = -1;
 241   fReconstructor = kFALSE;
 242 }
 243
 244 //______________________________________________
 245 void AliEMCALPID::RunPID(AliESDEvent *esd)
 246 {
 247 //
 248 // Make the PID for all the EMCAL clusters containedin the ESDs File
 249 // but just gamma/PiO/Hadron
 250 //
 251         // trivial check against NULL object passed
 252
 253   if (esd == 0x0) {
 254     AliInfo("NULL ESD object passed !!" );
 255     return ;
 256   }
 257
 258   Int_t nClusters = esd->GetNumberOfEMCALClusters();
 259   Int_t firstCluster = esd->GetFirstEMCALCluster();
 260   Double_t energy, lambda0;
 261   for (Int_t iCluster = firstCluster; iCluster < (nClusters + firstCluster); iCluster++) {
 262
 263     AliESDCaloCluster *clust = esd->GetCaloCluster(iCluster);
 264     energy = clust->E();
 265     lambda0 = clust->GetM02();
 266     // verify cluster type
 267     Int_t clusterType= clust->GetClusterType();
 268     if (clusterType == AliESDCaloCluster::kEMCALClusterv1 && lambda0 != 0  && energy < 1000) {
 269
 270
 271       // reject clusters with lambda0 = 0
 272
 273
 274       ComputePID(energy, lambda0);
 275
 276
 277       if (fPrintInfo) {
 278         AliInfo("___________________________________________________");
 279         AliInfo(Form( "Particle Energy = %f",energy));
 280         AliInfo(Form( "Particle Lambda0 of the particle = %f", lambda0) );
 281         AliInfo("PIDWeight of the particle :" );
 282         AliInfo(Form( " GAMMA  : %f",fPID[0] ));
 283         AliInfo(Form( " PiZero : %f",fPID[1] ));
 284         AliInfo(Form( " HADRON : %f", fPID[2] ));
 285         AliInfo("_________________________________________");
 286         AliInfo(Form( " kElectron : %f", fPIDFinal[0]) );
 287         AliInfo(Form( " kMuon     : %f", fPIDFinal[1] ));
 288         AliInfo(Form( " kPion       : %f", fPIDFinal[2] ));
 289         AliInfo(Form( " kKaon       : %f", fPIDFinal[3] ));
 290         AliInfo(Form( " kProton   : %f", fPIDFinal[4] ));
 291         AliInfo(Form( " kPhoton   : %f", fPIDFinal[5] ));
 292         AliInfo(Form( " kPi0        : %f", fPIDFinal[6] ));
 293         AliInfo(Form( " kNeutron  : %f", fPIDFinal[7] ));
 294         AliInfo(Form( " kKaon0  : %f", fPIDFinal[8] ));
 295         AliInfo(Form( " kEleCon   : %f", fPIDFinal[9] ));
 296         AliInfo(Form( " kUnknown  : %f", fPIDFinal[10] ));
 297         AliInfo("___________________________________________________");
 298       }
 299
 300       if(fReconstructor) // In case it is called during reconstruction.
 301         clust->SetPid(fPIDFinal);
 302     } // end if (clusterType...)
 303   } // end for (iCluster...)
 304 }
 305
 306 //__________________________________________________________
 307 void AliEMCALPID::ComputePID(Double_t energy, Double_t lambda0)
 308 {
 309 //
 310 // This is the main command, which uses the distributions computed and parametrised,
 311 // and gives the PID by the bayesian method.
 312 //
 313
 314 if (energy<5){energy =6;}
 315
 316
 317   TArrayD paramDistribGamma  = DistLambda0(energy, 1);
 318   TArrayD paramDistribPiZero = DistLambda0(energy, 2);
 319   TArrayD paramDistribHadron = DistLambda0(energy, 3);
 320
 321   Bool_t norm = kFALSE;
 322
 323   fProbGamma   = TMath::Gaus(lambda0, paramDistribGamma[1], paramDistribGamma[2], norm) * paramDistribGamma[0];
 324   fProbGamma  += TMath::Landau(lambda0, paramDistribGamma[4], paramDistribGamma[5], norm) * paramDistribGamma[3];
 325   fProbPiZero  = TMath::Gaus(lambda0, paramDistribPiZero[1], paramDistribPiZero[2], norm) * paramDistribPiZero[0];
 326   fProbPiZero += TMath::Landau(lambda0, paramDistribPiZero[4], paramDistribPiZero[5], norm) * paramDistribPiZero[3];
 327   fProbHadron  = TMath::Gaus(lambda0, paramDistribHadron[1], paramDistribHadron[2], norm) * paramDistribHadron[0];
 328   fProbHadron += TMath::Landau(lambda0, paramDistribHadron[4], paramDistribHadron[5], norm) * paramDistribHadron[3];
 329
 330   // compute PID Weight
 331   fPIDWeight[0] = fProbGamma / (fProbGamma + fProbPiZero + fProbHadron);
 332   fPIDWeight[1] = fProbPiZero / (fProbGamma+fProbPiZero+fProbHadron);
 333   fPIDWeight[2] = fProbHadron / (fProbGamma+fProbPiZero+fProbHadron);
 334
 335   SetPID(fPIDWeight[0], 0);
 336   SetPID(fPIDWeight[1], 1);
 337   SetPID(fPIDWeight[2], 2);
 338
 339   // sortie ecran pid Weight only for control (= in english ???)
 340   if (fPrintInfo) {
 341     AliInfo(Form( "Energy in loop = %f", energy) );
 342     AliInfo(Form( "Lambda0 in loop = %f", lambda0) );
 343     AliInfo(Form( "fProbGamma in loop = %f", fProbGamma) );
 344     //  AliInfo(Form( "fParametresDistribGamma[2] = %f", fParamDistribGamma[2]) );
 345     AliInfo(Form( "fProbaPiZero = %f", fProbPiZero ));
 346     AliInfo(Form( "fProbaHadron = %f", fProbHadron) );
 347     AliInfo(Form( "PIDWeight in loop = %f ||| %f ||| %f",  fPIDWeight[0] , fPIDWeight[1], fPIDWeight[2]) );
 348     AliInfo(Form( "fGamma[2][2] = %f", fGamma[2][2] ));
 349     AliInfo("********************************************************" );
 350   }
 351
 352   fPIDFinal[0]  = fPIDWeight[0]/2;
 353   fPIDFinal[1]  = fPIDWeight[2]/8;
 354   fPIDFinal[2]  = fPIDWeight[2]/8;
 355   fPIDFinal[3]  = fPIDWeight[2]/8;
 356   fPIDFinal[4]  = fPIDWeight[2]/8;
 357   fPIDFinal[5]  = fPIDWeight[0]/2;
 358   fPIDFinal[6]  = fPIDWeight[1]  ;
 359   fPIDFinal[7]  = fPIDWeight[2]/8;
 360   fPIDFinal[8]  = fPIDWeight[2]/8;
 361   fPIDFinal[9]  = fPIDWeight[2]/8;
 362   fPIDFinal[10] = fPIDWeight[2]/8;
 363 }
 364
 365 //________________________________________________________
 366 TArrayD AliEMCALPID::DistLambda0(Double_t energy, Int_t type)
 367 {
 368   //
 369   // Compute the values of the parametrised distributions using the data initialised before.
 370   //
 371   Double_t constGauss = 0., meanGauss = 0., sigmaGauss = 0.;
 372   Double_t constLandau=0., mpvLandau=0., sigmaLandau=0.;
 373   TArrayD  distributionParam(6);
 374
 375   switch (type) {
 376   case 1:
 377     constGauss  = Polynomial(energy, fGamma[0]);
 378     meanGauss   = Polynomial(energy, fGamma[1]);
 379     sigmaGauss  = Polynomial(energy, fGamma[2]);
 380     constLandau = Polynomial(energy, fGamma[3]);
 381     mpvLandau   = Polynomial(energy, fGamma[4]);
 382     sigmaLandau = Polynomial(energy, fGamma[5]);
 383     break;
 384   case 2:
 385     if (energy < 10) {
 386       constGauss  = Polynomial(energy, fPiZero5to10[0]);
 387       meanGauss   = Polynomial(energy, fPiZero5to10[1]);
 388       sigmaGauss  = Polynomial(energy, fPiZero5to10[2]);
 389       constLandau = Polynomial(energy, fPiZero5to10[3]);
 390       mpvLandau   = Polynomial(energy, fPiZero5to10[4]);
 391       sigmaLandau = Polynomial(energy, fPiZero5to10[5]);
 392     }
 393     else {
 394       constGauss  = Polynomial(energy, fPiZero10to60[0]);
 395       meanGauss   = Polynomial(energy, fPiZero10to60[1]);
 396       sigmaGauss  = Polynomial(energy, fPiZero10to60[2]);
 397       constLandau = Polynomial(energy, fPiZero10to60[3]);
 398       mpvLandau   = Polynomial(energy, fPiZero10to60[4]);
 399       sigmaLandau = Polynomial(energy, fPiZero10to60[5]);
 400     }
 401     break;
 402   case 3:
 403     constGauss  = Polynomial(energy, fHadron[0]);
 404     meanGauss   = Polynomial(energy, fHadron[1]);
 405     sigmaGauss  = Polynomial(energy, fHadron[2]);
 406     constLandau = Polynomial(energy, fHadron[3]);
 407     mpvLandau   = Polynomial(energy, fHadron[4]);
 408     sigmaLandau = Polynomial(energy, fHadron[5]);
 409     break;
 410   }
 411
 412   distributionParam[0] = constGauss;
 413   distributionParam[1] = meanGauss;
 414   distributionParam[2] = sigmaGauss;
 415   distributionParam[3] = constLandau;
 416   distributionParam[4] = mpvLandau;
 417   distributionParam[5] = sigmaLandau;
 418
 419   return distributionParam;
 420 }
 421
 422 //_______________________________________________________
 423 Double_t AliEMCALPID::Polynomial(Double_t x, Double_t *params)
 424 {
 425   //
 426   // Compute a polynomial for a given value of 'x'
 427   // with the array of parameters passed as the second arg
 428   //
 429   Double_t y;
 430   y  = params[0];
 431   y += params[1] * x;
 432   y += params[2] * x * x;
 433   y += params[3] * x * x * x;
 434   y += params[4] * x * x * x * x;
 435   y += params[5] * x * x * x * x * x;
 436
 437   return y;
 438 }