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
  18 //   Compute PID weights for all the clusters that are in AliESDs.root file
  19 //   the AliESDs.root have to be in the same directory as the class
  20 //
  21 //   and do:
  22 //   AliEMCALPID *pid = new AliEMCALPID(kFALSE); // this calls the constructor which avoids the call to recparam
  23 //   pid->SetReconstructor(kFALSE);
  24 //   pid->SetPrintInfo(kTRUE);
  25 //   pid->SetHighFluxParam(); //   pid->SetLowFluxParam();
  26 //
  27 //   then in cluster loop do
  28 //   pid->ComputePID(energy, lambda0);
  29 //
  30 //        Compute PID Weight for all clusters in AliESDs.root file
  31 //        keep this function for the moment for a simple verification, could be removed
  32 //
  33 //   pid->GetPIDFinal(idx) gives the probabilities
  34 //
  35 //   Double_t PIDFinal[AliPID::kSPECIESCN]  is the standard PID for :
  36 //
  37 //      kElectron :  fPIDFinal[0]
  38 //      kMuon     :  fPIDFinal[1]
  39 //      kPion       :  fPIDFinal[2]
  40 //      kKaon       :  fPIDFinal[3]
  41 //      kProton   :  fPIDFinal[4]
  42 //      kPhoton   :  fPIDFinal[5]
  43 //      kPi0        :  fPIDFinal[6]
  44 //      kNeutron  :  fPIDFinal[7]
  45 //      kKaon0    :  fPIDFinal[8]
  46 //      kEleCon   :  fPIDFinal[9]
  47 //      kUnknown  :  fPIDFinal[10]
  48 //
  49 //
  50 //    PID[3] is a simple PID for
  51 //      Electron & Photon  PID[0]
  52 //                    Pi0  PID[1]
  53 //                 Hadron  PID[2]
  54 //
  55 // Author: Genole Bourdaud 2007 (SUBATECH)
  56 //         Marie Germain 07/2009 (SUBATECH), new parametrization for low and high flux environment
  57 //         Gustavo Conesa 08/2009 (LNF), divide class in AliEMCALPID and AliEMCALPIDUtils, PIDUtils belong to library EMCALUtils
  58 // --- standard c ---
  59
  60 // standard C++ includes
  61 //#include <Riostream.h>
  62
  63 // ROOT includes
  64
  65 // STEER includes
  66 #include "AliESDEvent.h"
  67 #include "AliEMCALPID.h"
  68 #include "AliESDCaloCluster.h"
  69 #include "AliEMCALReconstructor.h"
  70
  71
  72 ClassImp(AliEMCALPID)
  73
  74 //______________________________________________
  75   AliEMCALPID::AliEMCALPID()
  76         : AliEMCALPIDUtils(), fReconstructor(kTRUE)
  77 {
  78   //
  79   // Constructor.
  80   // Initialize all constant values which have to be used
  81   // during PID algorithm execution
  82   //
  83
  84   InitParameters();
  85
  86
  87 }
  88
  89 //______________________________________________
  90 AliEMCALPID::AliEMCALPID(Bool_t reconstructor)
  91 : AliEMCALPIDUtils(), fReconstructor(reconstructor)
  92 {
  93   //
  94   // Constructor.
  95   // Initialize all constant values which have to be used
  96   // during PID algorithm execution called when used in standalone mode
  97   //
  98
  99   InitParameters();
 100
 101 }
 102
 103 //______________________________________________
 104 void AliEMCALPID::RunPID(AliESDEvent *esd)
 105 {
 106   //
 107   // Make the PID for all the EMCAL clusters containedin the ESDs File
 108   // but just gamma/PiO/Hadron
 109   //
 110   // trivial check against NULL object passed
 111
 112   if (esd == 0x0) {
 113     AliInfo("NULL ESD object passed !!" );
 114     return ;
 115   }
 116
 117   Int_t nClusters = esd->GetNumberOfCaloClusters();
 118   Int_t firstCluster = 0;
 119   Double_t energy=0., lambda0=0.;
 120   for (Int_t iCluster = firstCluster; iCluster < (nClusters + firstCluster); iCluster++) {
 121
 122     AliESDCaloCluster *clust = esd->GetCaloCluster(iCluster);
 123     if (!clust->IsEMCAL()) continue ;
 124
 125     energy = clust->E();
 126     lambda0 = clust->GetM02();
 127
 128     if (lambda0 != 0  && energy < 1000) {
 129
 130       // reject clusters with lambda0 = 0
 131
 132
 133       ComputePID(energy, lambda0);
 134
 135
 136       if (fPrintInfo) {
 137         AliInfo("___________________________________________________");
 138         AliInfo(Form( "Particle Energy = %f",energy));
 139         AliInfo(Form( "Particle Lambda0 of the particle = %f", lambda0) );
 140         AliInfo("PIDWeight of the particle :" );
 141         AliInfo(Form( " GAMMA  : %f",fPID[0] ));
 142         AliInfo(Form( " PiZero : %f",fPID[1] ));
 143         AliInfo(Form( " HADRON : %f", fPID[2] ));
 144         AliInfo("_________________________________________");
 145         AliInfo(Form( " kElectron : %f", fPIDFinal[0]) );
 146         AliInfo(Form( " kMuon     : %f", fPIDFinal[1] ));
 147         AliInfo(Form( " kPion       : %f", fPIDFinal[2] ));
 148         AliInfo(Form( " kKaon       : %f", fPIDFinal[3] ));
 149         AliInfo(Form( " kProton   : %f", fPIDFinal[4] ));
 150         AliInfo(Form( " kPhoton   : %f", fPIDFinal[5] ));
 151         AliInfo(Form( " kPi0        : %f", fPIDFinal[6] ));
 152         AliInfo(Form( " kNeutron  : %f", fPIDFinal[7] ));
 153         AliInfo(Form( " kKaon0  : %f", fPIDFinal[8] ));
 154         AliInfo(Form( " kEleCon   : %f", fPIDFinal[9] ));
 155         AliInfo(Form( " kUnknown  : %f", fPIDFinal[10] ));
 156         AliInfo("___________________________________________________");
 157       }
 158
 159       if(fReconstructor){ // In case it is called during reconstruction.
 160         //      cout << "############# Fill ESDs with PIDWeight ##########" << endl;
 161         clust->SetPID(fPIDFinal);}
 162     } // end if (lambda0 != 0  && energy < 1000)
 163   } // end for (iCluster...)
 164 }
 165
 166
 167 //_______________________________________________________
 168 void AliEMCALPID::InitParameters()
 169 {
 170   // Initialize PID parameters, depending on the use or not of the reconstructor
 171   // and the kind of event type if the reconstructor is not used.
 172   //  fWeightHadronEnergy=0.;
 173   //  fWeightPiZeroEnergy=0.;
 174   //  fWeightGammaEnergy=0.;
 175
 176   fPIDWeight[0] = -1;
 177   fPIDWeight[1] = -1;
 178   fPIDWeight[2] = -1;
 179
 180   for(Int_t i=0; i<AliPID::kSPECIESCN+1; i++)
 181     fPIDFinal[i]= 0;
 182
 183   const AliEMCALRecParam* recParam = AliEMCALReconstructor::GetRecParam();
 184
 185   if(fReconstructor){
 186
 187     if(!recParam) {
 188       AliFatal("Reconstruction parameters for EMCAL not set!");
 189     }
 190     else {
 191
 192       for(Int_t i=0; i<6; i++){
 193         for(Int_t j=0; j<6; j++){
 194           fGamma[i][j]       = recParam->GetGamma(i,j);
 195           fGamma1to10[i][j]  = recParam->GetGamma1to10(i,j);
 196           fHadron[i][j]      = recParam->GetHadron(i,j);
 197           fHadron1to10[i][j] = recParam->GetHadron1to10(i,j);
 198           fPiZero[i][j]      = recParam->GetPiZero(i,j);
 199
 200
 201           //    AliDebug(1,Form("PID parameters (%d, %d): fGamma=%.3f, fPi=%.3f, fHadron=%.3f",
 202           //                    i,j, fGamma[i][j],fPiZero[i][j],fHadron[i][j] ));
 203           //    cout << "PID parameters (" << i << " ,"<<j<<") fGamma= "<<  fGamma[i][j]<<" fPi0 ="<<  fPiZero[i][j]<< endl;
 204
 205         } // end loop j
 206         fHadronEnergyProb[i] = recParam->GetHadronEnergyProb(i);
 207         fPiZeroEnergyProb[i] = recParam->GetPiZeroEnergyProb(i);
 208         fGammaEnergyProb[i]  = recParam->GetGammaEnergyProb(i);
 209       } //end loop i
 210
 211
 212     } // end if !recparam
 213
 214   }
 215
 216   else{
 217     //   init the parameters here instead of from loading from recparam
 218     //   default parameters are PbPb parameters.
 219     SetHighFluxParam();
 220
 221   }
 222
 223 }
 224