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.16  2007/11/23 13:39:05  gustavo
  21  * Track matching and PID parameters added to AliEMCALRecParam
  22  *
  23  * Revision 1.15  2007/10/09 08:46:10  hristov
  24  * The data members fEMCALClusterCluster and fPHOSCluster are removed from AliESDCaloCluster, the fClusterType is used to select PHOS or EMCAL clusters. Changes, needed to use correctly the new AliESDCaloCluster. (Christian)
  25  *
  26  * Revision 1.14  2007/07/26 16:54:53  morsch
  27  * Changes in AliESDEvent fwd declarartions.
  28  *
  29  * Revision 1.13  2007/07/11 13:43:29  hristov
  30  * New class AliESDEvent, backward compatibility with the old AliESD (Christian)
  31  *
  32  * Revision 1.12  2007/06/11 20:43:06  hristov
  33  * Changes required by the updated AliESDCaloCluster (Gustavo)
  34  *
  35  * Revision 1.11  2007/03/30 13:50:34  gustavo
  36  * PID for particles with E < 5 GeV was not done, temporal solution found (Guenole)
  37  *
  38  * Revision 1.10  2007/03/09 14:34:11  gustavo
  39  * Correct probability calculation, added missing initialization of data members
  40  *
  41  * Revision 1.9  2007/02/20 20:17:43  hristov
  42  * Corrected array size, removed warnings (icc)
  43  *
  44  * Revision 1.8  2006/12/19 08:49:35  gustavo
  45  * New PID class for EMCAL, bayesian analysis done with ESD data, PID information filled when calling AliEMCALPID in AliEMCALReconstructor::FillESD()
  46  *
  47  *
  48  */
  49 //    to compute PID for all the clusters in ESDs.root file
  50 //     the ESDs.root have to be in the same directory as the class
  51 //
  52 //
  53 //
  54 //
  55 //
  56 //      AliEMCALPID::CalculPID(Energy,Lambda0)
  57 //        Calcul PID for all clusters in AliESDs.root file
  58 //        keep this function for the moment for a simple verification, could be removed
  59 //
  60 //
  61 //
  62 //   AliEMCALPID::CalculPID(Energy,Lambda0)
  63 //    calcul PID Weght for a cluster with Energy, Lambda0 .
  64 //    Double_t PIDFinal[AliPID::kSPECIESN]  is the standard PID for :
  65 //
  66 //
  67 //
  68 //      kElectron :  fPIDFinal[0]
  69 //      kMuon     :  fPIDFinal[1]
  70 //      kPion     :  fPIDFinal[2]
  71 //      kKaon     :  fPIDFinal[3]
  72 //      kProton   :  fPIDFinal[4]
  73 //      kPhoton   :  fPIDFinal[5]
  74 //      kPi0      :  fPIDFinal[6]
  75 //      kNeutron  :  fPIDFinal[7]
  76 //      kKaon0    :  fPIDFinal[8]
  77 //      kEleCon   :  fPIDFinal[9]
  78 //      kUnknown  :  fPIDFinal[10]
  79 //
  80 //
  81 //    PID[3] is a simple PID for
  82 //      Electron & Photon  PID[0]
  83 //                    Pi0  PID[1]
  84 //                 Hadron  PID[2]
  85 //
  86 //
  87 //
  88 //
  89 //
  90 // --- ROOT system ---
  91
  92 // standard C++ includes
  93 #include <Riostream.h>
  94
  95 // ROOT includes
  96 #include "TTree.h"
  97 #include "TStyle.h"
  98 #include "TVector3.h"
  99 #include "TBranch.h"
 100 #include "TClonesArray.h"
 101 #include "TCanvas.h"
 102 #include "TLorentzVector.h"
 103 #include "TMath.h"
 104 #include "TFile.h"
 105 #include "TH1.h"
 106 #include "TH2.h"
 107 #include "TParticle.h"
 108
 109 // STEER includes
 110 #include "AliLog.h"
 111 #include "AliEMCALPID.h"
 112 #include "AliESDCaloCluster.h"
 113 #include "AliEMCALRecParam.h"
 114 #include "AliEMCALReconstructor.h"
 115
 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   fPIDWeight[0] = -1;
 130   fPIDWeight[1] = -1;
 131   fPIDWeight[2] = -1;
 132
 133   for(Int_t i=0; i<AliPID::kSPECIESN+1; i++)
 134     fPIDFinal[i]= 0;
 135
 136   const AliEMCALRecParam* recParam = AliEMCALReconstructor::GetRecParam();
 137   if(!recParam) {
 138     AliFatal("Reconstruction parameters for EMCAL not set!");
 139   }
 140   else {
 141     for(Int_t i=0; i<6; i++){
 142       for(Int_t j=0; j<6; j++){
 143         fGamma[i][j] = recParam->GetGamma(i,j);
 144         fHadron[i][j] = recParam->GetHadron(i,j);
 145         fPiZero5to10[i][j] = recParam->GetPiZero5to10(i,j);
 146         fPiZero10to60[i][j] = recParam->GetPiZero10to60(i,j);
 147         AliDebug(1,Form("PID parameters (%d, %d): fGamma=%.3f, fPi=%.3f, fHadron=%.3f",
 148                         i,j, fGamma[i][j],fPiZero5to10[i][j],fHadron[i][j] ));
 149       }
 150     }
 151
 152   }
 153
 154 }
 155
 156 //______________________________________________
 157 void AliEMCALPID::RunPID(AliESDEvent *esd)
 158 {
 159 //
 160 // Make the PID for all the EMCAL clusters containedin the ESDs File
 161 // but just gamma/PiO/Hadron
 162 //
 163         // trivial check against NULL object passed
 164
 165   if (esd == 0x0) {
 166     AliInfo("NULL ESD object passed !!" );
 167     return ;
 168   }
 169
 170   Int_t nClusters = esd->GetNumberOfEMCALClusters();
 171   Int_t firstCluster = esd->GetFirstEMCALCluster();
 172   Double_t energy, lambda0;
 173   for (Int_t iCluster = firstCluster; iCluster < (nClusters + firstCluster); iCluster++) {
 174
 175     AliESDCaloCluster *clust = esd->GetCaloCluster(iCluster);
 176     energy = clust->E();
 177     lambda0 = clust->GetM02();
 178     // verify cluster type
 179     Int_t clusterType= clust->GetClusterType();
 180     if (clusterType == AliESDCaloCluster::kEMCALClusterv1 && lambda0 != 0  && energy < 1000) {
 181
 182
 183       // reject clusters with lambda0 = 0
 184
 185
 186       ComputePID(energy, lambda0);
 187
 188
 189       if (fPrintInfo) {
 190         AliInfo("___________________________________________________");
 191         AliInfo(Form( "Particle Energy = %f",energy));
 192         AliInfo(Form( "Particle Lambda0 of the particle = %f", lambda0) );
 193         AliInfo("PIDWeight of the particle :" );
 194         AliInfo(Form( " GAMMA  : %f",fPID[0] ));
 195         AliInfo(Form( " PiZero : %f",fPID[1] ));
 196         AliInfo(Form( " HADRON : %f", fPID[2] ));
 197         AliInfo("_________________________________________");
 198         AliInfo(Form( " kElectron : %f", fPIDFinal[0]) );
 199         AliInfo(Form( " kMuon     : %f", fPIDFinal[1] ));
 200         AliInfo(Form( " kPion       : %f", fPIDFinal[2] ));
 201         AliInfo(Form( " kKaon       : %f", fPIDFinal[3] ));
 202         AliInfo(Form( " kProton   : %f", fPIDFinal[4] ));
 203         AliInfo(Form( " kPhoton   : %f", fPIDFinal[5] ));
 204         AliInfo(Form( " kPi0        : %f", fPIDFinal[6] ));
 205         AliInfo(Form( " kNeutron  : %f", fPIDFinal[7] ));
 206         AliInfo(Form( " kKaon0  : %f", fPIDFinal[8] ));
 207         AliInfo(Form( " kEleCon   : %f", fPIDFinal[9] ));
 208         AliInfo(Form( " kUnknown  : %f", fPIDFinal[10] ));
 209         AliInfo("___________________________________________________");
 210       }
 211
 212       if(fReconstructor) // In case it is called during reconstruction.
 213         clust->SetPid(fPIDFinal);
 214     } // end if (clusterType...)
 215   } // end for (iCluster...)
 216 }
 217
 218 //__________________________________________________________
 219 void AliEMCALPID::ComputePID(Double_t energy, Double_t lambda0)
 220 {
 221 //
 222 // This is the main command, which uses the distributions computed and parametrised,
 223 // and gives the PID by the bayesian method.
 224 //
 225
 226 if (energy<5){energy =6;}
 227
 228
 229   TArrayD paramDistribGamma  = DistLambda0(energy, 1);
 230   TArrayD paramDistribPiZero = DistLambda0(energy, 2);
 231   TArrayD paramDistribHadron = DistLambda0(energy, 3);
 232
 233   Bool_t norm = kFALSE;
 234
 235   fProbGamma   = TMath::Gaus(lambda0, paramDistribGamma[1], paramDistribGamma[2], norm) * paramDistribGamma[0];
 236   fProbGamma  += TMath::Landau(lambda0, paramDistribGamma[4], paramDistribGamma[5], norm) * paramDistribGamma[3];
 237   fProbPiZero  = TMath::Gaus(lambda0, paramDistribPiZero[1], paramDistribPiZero[2], norm) * paramDistribPiZero[0];
 238   fProbPiZero += TMath::Landau(lambda0, paramDistribPiZero[4], paramDistribPiZero[5], norm) * paramDistribPiZero[3];
 239   fProbHadron  = TMath::Gaus(lambda0, paramDistribHadron[1], paramDistribHadron[2], norm) * paramDistribHadron[0];
 240   fProbHadron += TMath::Landau(lambda0, paramDistribHadron[4], paramDistribHadron[5], norm) * paramDistribHadron[3];
 241
 242   // compute PID Weight
 243   fPIDWeight[0] = fProbGamma / (fProbGamma + fProbPiZero + fProbHadron);
 244   fPIDWeight[1] = fProbPiZero / (fProbGamma+fProbPiZero+fProbHadron);
 245   fPIDWeight[2] = fProbHadron / (fProbGamma+fProbPiZero+fProbHadron);
 246
 247   SetPID(fPIDWeight[0], 0);
 248   SetPID(fPIDWeight[1], 1);
 249   SetPID(fPIDWeight[2], 2);
 250
 251   // sortie ecran pid Weight only for control (= in english ???)
 252   if (fPrintInfo) {
 253     AliInfo(Form( "Energy in loop = %f", energy) );
 254     AliInfo(Form( "Lambda0 in loop = %f", lambda0) );
 255     AliInfo(Form( "fProbGamma in loop = %f", fProbGamma) );
 256     //  AliInfo(Form( "fParametresDistribGamma[2] = %f", fParamDistribGamma[2]) );
 257     AliInfo(Form( "fProbaPiZero = %f", fProbPiZero ));
 258     AliInfo(Form( "fProbaHadron = %f", fProbHadron) );
 259     AliInfo(Form( "PIDWeight in loop = %f ||| %f ||| %f",  fPIDWeight[0] , fPIDWeight[1], fPIDWeight[2]) );
 260     AliInfo(Form( "fGamma[2][2] = %f", fGamma[2][2] ));
 261     AliInfo("********************************************************" );
 262   }
 263
 264   fPIDFinal[0]  = fPIDWeight[0]/2;
 265   fPIDFinal[1]  = fPIDWeight[2]/8;
 266   fPIDFinal[2]  = fPIDWeight[2]/8;
 267   fPIDFinal[3]  = fPIDWeight[2]/8;
 268   fPIDFinal[4]  = fPIDWeight[2]/8;
 269   fPIDFinal[5]  = fPIDWeight[0]/2;
 270   fPIDFinal[6]  = fPIDWeight[1]  ;
 271   fPIDFinal[7]  = fPIDWeight[2]/8;
 272   fPIDFinal[8]  = fPIDWeight[2]/8;
 273   fPIDFinal[9]  = fPIDWeight[2]/8;
 274   fPIDFinal[10] = fPIDWeight[2]/8;
 275 }
 276
 277 //________________________________________________________
 278 TArrayD AliEMCALPID::DistLambda0(Double_t energy, Int_t type)
 279 {
 280   //
 281   // Compute the values of the parametrised distributions using the data initialised before.
 282   //
 283   Double_t constGauss = 0., meanGauss = 0., sigmaGauss = 0.;
 284   Double_t constLandau=0., mpvLandau=0., sigmaLandau=0.;
 285   TArrayD  distributionParam(6);
 286
 287   switch (type) {
 288   case 1:
 289     constGauss  = Polynomial(energy, fGamma[0]);
 290     meanGauss   = Polynomial(energy, fGamma[1]);
 291     sigmaGauss  = Polynomial(energy, fGamma[2]);
 292     constLandau = Polynomial(energy, fGamma[3]);
 293     mpvLandau   = Polynomial(energy, fGamma[4]);
 294     sigmaLandau = Polynomial(energy, fGamma[5]);
 295
 296
 297     break;
 298   case 2:
 299     if (energy < 10) {
 300       constGauss  = Polynomial(energy, fPiZero5to10[0]);
 301       meanGauss   = Polynomial(energy, fPiZero5to10[1]);
 302       sigmaGauss  = Polynomial(energy, fPiZero5to10[2]);
 303       constLandau = Polynomial(energy, fPiZero5to10[3]);
 304       mpvLandau   = Polynomial(energy, fPiZero5to10[4]);
 305       sigmaLandau = Polynomial(energy, fPiZero5to10[5]);
 306
 307     }
 308     else {
 309       constGauss  = Polynomial(energy, fPiZero10to60[0]);
 310       meanGauss   = Polynomial(energy, fPiZero10to60[1]);
 311       sigmaGauss  = Polynomial(energy, fPiZero10to60[2]);
 312       constLandau = Polynomial(energy, fPiZero10to60[3]);
 313       mpvLandau   = Polynomial(energy, fPiZero10to60[4]);
 314       sigmaLandau = Polynomial(energy, fPiZero10to60[5]);
 315
 316     }
 317     break;
 318   case 3:
 319     constGauss  = Polynomial(energy, fHadron[0]);
 320     meanGauss   = Polynomial(energy, fHadron[1]);
 321     sigmaGauss  = Polynomial(energy, fHadron[2]);
 322     constLandau = Polynomial(energy, fHadron[3]);
 323     mpvLandau   = Polynomial(energy, fHadron[4]);
 324     sigmaLandau = Polynomial(energy, fHadron[5]);
 325
 326     break;
 327   }
 328
 329   distributionParam[0] = constGauss;
 330   distributionParam[1] = meanGauss;
 331   distributionParam[2] = sigmaGauss;
 332   distributionParam[3] = constLandau;
 333   distributionParam[4] = mpvLandau;
 334   distributionParam[5] = sigmaLandau;
 335
 336   return distributionParam;
 337 }
 338
 339 //_______________________________________________________
 340 Double_t AliEMCALPID::Polynomial(Double_t x, Double_t *params)
 341 {
 342   //
 343   // Compute a polynomial for a given value of 'x'
 344   // with the array of parameters passed as the second arg
 345   //
 346
 347   Double_t y;
 348   y  = params[0];
 349   y += params[1] * x;
 350   y += params[2] * x * x;
 351   y += params[3] * x * x * x;
 352   y += params[4] * x * x * x * x;
 353   y += params[5] * x * x * x * x * x;
 354
 355   return y;
 356 }