PWGHF/hfe/AliHFEemcalPIDqa.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 // Class AliHFEemcalPIDqa
  17 //
  18 // Monitoring EMCAL PID in the HFE PID montioring framework. The following
  19 // quantities are monitored:
  20 //   TPC Signal distribution
  21 // (Always as function of momentum, particle species and centrality
  22 // before and after cut)
  23 // More information about the PID monitoring framework can be found in
  24 // AliHFEpidQAmanager.cxx and AliHFEdetPIDqa.cxx
  25 //
  26 // Author:
  27 //
  28 //   Shingo Sakai <ssakai@lbl.gov>
  29 #include <TClass.h>
  30 #include <TH2.h>
  31 #include <THnSparse.h>
  32 #include <TString.h>
  33
  34 #include <TMath.h>
  35 #include "AliESDInputHandler.h"
  36 //#include "AliVCluster.h"
  37 //#include "AliVCaloCells.h"
  38 //#include "AliVEvent.h"
  39 #include "AliMagF.h"
  40 #include "AliPIDResponse.h"
  41
  42 #include "AliLog.h"
  43 #include "AliPID.h"
  44 #include "AliVParticle.h"
  45 //#include "AliVTrack.h"
  46 //#include "AliESDtrack.h"
  47 #include "AliHFEcollection.h"
  48 #include "AliHFEpid.h"
  49 #include "AliHFEpidBase.h"
  50 #include "AliHFEpidQAmanager.h"
  51 #include "AliHFEpidTPC.h"
  52 #include "AliHFEpidEMCAL.h"
  53 #include "AliHFEtools.h"
  54 #include "AliHFEemcalPIDqa.h"
  55 #include "AliTrackerBase.h"
  56
  57 ClassImp(AliHFEemcalPIDqa)
  58
  59 //_________________________________________________________
  60 AliHFEemcalPIDqa::AliHFEemcalPIDqa():
  61     AliHFEdetPIDqa()
  62   , fHistos(NULL)
  63 {
  64   //
  65   // Dummy constructor
  66   //
  67 }
  68
  69 //_________________________________________________________
  70 AliHFEemcalPIDqa::AliHFEemcalPIDqa(const char* name):
  71     AliHFEdetPIDqa(name, "QA for EMCAL")
  72   , fHistos(NULL)
  73 {
  74   //
  75   // Default constructor
  76   //
  77 }
  78
  79 //_________________________________________________________
  80 AliHFEemcalPIDqa::AliHFEemcalPIDqa(const AliHFEemcalPIDqa &o):
  81     AliHFEdetPIDqa(o)
  82   , fHistos()
  83 {
  84   //
  85   // Copy constructor
  86   //
  87   o.Copy(*this);
  88 }
  89
  90 //_________________________________________________________
  91 AliHFEemcalPIDqa &AliHFEemcalPIDqa::operator=(const AliHFEemcalPIDqa &o){
  92   //
  93   // Do assignment
  94   //
  95   AliHFEdetPIDqa::operator=(o);
  96   if(&o != this) o.Copy(*this);
  97   return *this;
  98 }
  99
 100 //_________________________________________________________
 101 AliHFEemcalPIDqa::~AliHFEemcalPIDqa(){
 102   //
 103   // Destructor
 104   //
 105   if(fHistos) delete fHistos;
 106 }
 107
 108 //_________________________________________________________
 109 void AliHFEemcalPIDqa::Copy(TObject &o) const {
 110   //
 111   // Make copy
 112   //
 113   AliHFEemcalPIDqa &target = dynamic_cast<AliHFEemcalPIDqa &>(o);
 114   if(target.fHistos){
 115     delete target.fHistos;
 116     target.fHistos = NULL;
 117   }
 118   if(fHistos) target.fHistos = new AliHFEcollection(*fHistos);
 119 }
 120
 121 //_________________________________________________________
 122 Long64_t AliHFEemcalPIDqa::Merge(TCollection *coll){
 123   //
 124   // Merge with other objects
 125   //
 126   if(!coll) return 0;
 127   if(coll->IsEmpty()) return 1;
 128
 129   TIter it(coll);
 130   AliHFEemcalPIDqa *refQA = NULL;
 131   TObject *o = NULL;
 132   Long64_t count = 0;
 133   TList listHistos;
 134   while((o = it())){
 135     refQA = dynamic_cast<AliHFEemcalPIDqa *>(o);
 136     if(!refQA) continue;
 137
 138     listHistos.Add(refQA->fHistos);
 139     count++;
 140   }
 141   fHistos->Merge(&listHistos);
 142   return count + 1;
 143 }
 144
 145 //_________________________________________________________
 146 void AliHFEemcalPIDqa::Initialize(){
 147   //
 148   // Define Histograms
 149   //
 150
 151   fHistos = new AliHFEcollection("emcalqahistos", "Collection of EMCAL QA histograms");
 152
 153   // Make common binning
 154   const Int_t kCentralityBins = 11;
 155   const Double_t kMinP = 0.;
 156   const Double_t kMaxP = 50.;
 157   const Double_t kTPCSigMim = 40.;
 158   const Double_t kTPCSigMax = 140.;
 159
 160   // 1st histogram: TPC dEdx with/without EMCAL (p, pT, TPC Signal, phi, eta,  Sig,  e/p,  ,match, cell, M02, M20, Disp, Centrality, select)
 161   Int_t nBins[16] = {AliPID::kSPECIES + 1, 500, 500,          400, 300,   100,   600,  300, 100,   40,   300, 300, 300, kCentralityBins, 6, 2};
 162   Double_t min[16] = {-1,               kMinP, kMinP,  kTPCSigMim,  0.,  -1.0,  -8.0,    0,   0,    0,   0.0, 0.0, 0.0,   0, -3.0, 0.};
 163   Double_t max[16] = {AliPID::kSPECIES, kMaxP, kMaxP,  kTPCSigMax, 6.0,   1.0,   4.0,  3.0, 0.1,   40,   3.0, 3.0, 3.0,   11., 3.0, 2.};
 164   fHistos->CreateTHnSparse("EMCAL_TPCdedx", "EMCAL signal; species; p [GeV/c]; pT [GeV/c] ; TPC signal [a.u.]; phi ; eta ; nSig ; E/p ; Rmatch ; Ncell ; M02 ; M20 ; Disp ; Centrality;charge ;PID Step; ", 16, nBins, min, max);
 165
 166   //2nd histogram: EMCAL signal - E/p
 167   //Int_t nBins2[7] = {AliPID::kSPECIES + 1, 500, 500, 500, 100, 600, 2};
 168   //Double_t min2[7] = {-1, kMinP, kMinP, 0,  0, -8., 0.};
 169   //Double_t max2[7] = {AliPID::kSPECIES, kMaxP, kMaxP, 5,  0.1, 4., 2.};
 170   //fHistos->CreateTHnSparse("EMCAL_Signal", "EMCAL true signal; species; p [GeV/c]; pT [GeV/c] ; E/p; Rmatch; TPCnsigma; PID Step", 7, nBins2, min2, max2);
 171
 172 }
 173
 174
 175
 176
 177 //_________________________________________________________
 178 void AliHFEemcalPIDqa::ProcessTrack(const AliHFEpidObject *track,AliHFEdetPIDqa::EStep_t step){
 179   //
 180   // Fill TPC histograms
 181   //
 182   //AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
 183   Float_t centrality = track->GetCentrality();
 184
 185   //const AliVTrack *vtrack = dynamic_cast<const AliVTrack *>(track->GetRecTrack());
 186   //const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(vtrack);
 187   const AliESDtrack *esdtrack = static_cast<const AliESDtrack *>(track->GetRecTrack());
 188
 189   Int_t species = track->GetAbInitioPID();
 190   if(species >= AliPID::kSPECIES) species = -1;
 191
 192   // Get TPC nSigma (based on AliHFEtpcPIDqa)
 193   AliHFEpidTPC *tpcpid = dynamic_cast<AliHFEpidTPC *>(fQAmanager->GetDetectorPID(AliHFEpid::kTPCpid));
 194   const AliPIDResponse *pidResponse = tpcpid ? tpcpid->GetPIDResponse() : NULL;
 195   Float_t nSigmatpc = pidResponse ? pidResponse->NumberOfSigmasTPC(static_cast<const AliVTrack *>(track->GetRecTrack()), AliPID::kElectron) : 0.;
 196   //printf("nSigmatpc = %f\n",nSigmatpc);
 197   AliDebug(2, Form("nSigmatpc = %f\n",nSigmatpc));
 198   //
 199
 200   double charge = esdtrack->Charge();
 201   //printf("charge %f\n",charge);
 202
 203
 204   //printf("phi %f, eta %f\n;",phi,eta);
 205
 206   // Get E/p
 207   Double_t showerEMC[4];
 208   TVector3 emcsignal = MomentumEnergyMatchV2(esdtrack,showerEMC);
 209   AliDebug(2, Form("shower shape ; %f %f %f %f \n",showerEMC[0],showerEMC[1],showerEMC[2],showerEMC[3]));
 210
 211   /*
 212   Double_t contentSignal2[7];
 213   contentSignal2[0] = species;
 214   contentSignal2[1] = track->GetRecTrack()->P();
 215   contentSignal2[2] = track->GetRecTrack()->Pt();
 216   contentSignal2[3] = emcsignal(0);//e over p
 217   contentSignal2[4] = emcsignal(1);//residual
 218   contentSignal2[5] = nSigmatpc;
 219   contentSignal2[6] = step == AliHFEdetPIDqa::kBeforePID ? 0. : 1.;
 220   */
 221   // QA array
 222   Double_t contentSignal[16];
 223   contentSignal[0] = species;
 224   contentSignal[1] = track->GetRecTrack()->P();
 225   contentSignal[2] = track->GetRecTrack()->Pt();
 226   contentSignal[3] = esdtrack->GetTPCsignal(); //?
 227   double phi  = track->GetRecTrack()->Phi();
 228   double eta  = track->GetRecTrack()->Eta();
 229   contentSignal[4] = phi;
 230   contentSignal[5] = eta;
 231   contentSignal[6] = nSigmatpc;
 232   contentSignal[7] = emcsignal(0);
 233   contentSignal[8] = emcsignal(2);
 234   contentSignal[9] = showerEMC[0];
 235   contentSignal[10] = showerEMC[1];
 236   contentSignal[11] = showerEMC[2];
 237   contentSignal[12] = showerEMC[3];
 238   contentSignal[13] = centrality;
 239   contentSignal[14] = charge;
 240   contentSignal[15] = step == AliHFEdetPIDqa::kBeforePID ? 0. : 1.;
 241
 242
 243   //printf("ProcessTrack ; Print Content %g; %g; %g; %g \n",contentSignal[0],contentSignal[1],contentSignal[2],contentSignal[3]);
 244   fHistos->Fill("EMCAL_TPCdedx", contentSignal);
 245   //fHistos->Fill("EMCAL_Signal", contentSignal2);
 246 }
 247
 248 //_________________________________________________________
 249 TH1 *AliHFEemcalPIDqa::GetHistogram(const char *name){
 250   //
 251   // Get the histogram
 252   //
 253   if(!fHistos) return NULL;
 254   TObject *histo = fHistos->Get(name);
 255   if(!histo->InheritsFrom("TH1")) return NULL;
 256   return dynamic_cast<TH1 *>(histo);
 257 }
 258
 259
 260 //___________________________________________________________________________
 261 TVector3 AliHFEemcalPIDqa::MomentumEnergyMatchV2(const AliESDtrack *esdtrack, Double_t *shower) const
 262 { // Returns e/p & Rmatch
 263   /*
 264   Float_t  clsPos[3];
 265   Double_t trkPos[3];
 266   Double_t fMass=0.139;
 267   Double_t fStep=10; //This is taken from EMCAL tender, hardcoded!
 268   */
 269   TVector3 refVec(-9999,-9999,-9999);
 270   Double_t matchclsE = 9999.9;
 271   for(int i=0; i<4; i++)shower[i] = -9999;
 272
 273   const AliESDEvent *evt = esdtrack->GetESDEvent();
 274
 275   //Proper trackmatching, added by Tomas
 276   //Get matched cluster
 277   Int_t icl = esdtrack->GetEMCALcluster(); //From tender
 278   AliVCluster *cluster = (AliVCluster*) evt->GetCaloCluster(icl);
 279   if(cluster==NULL) return refVec;
 280   if(!cluster->IsEMCAL()) return refVec;
 281   //printf("EMCal cluster pointer ? %p\n",(void*)cluster);
 282   //printf("EMCal cluster ? %d\n",cluster->IsEMCAL());
 283
 284   // from ESDs
 285   double delphi = cluster->GetTrackDx();
 286   double deleta = cluster->GetTrackDz();
 287   double rmatch1 = sqrt(pow(delphi,2)+pow(deleta,2));
 288
 289   matchclsE = cluster->E();
 290
 291   //double feop = -9999.9;
 292   //if(matchclsE<9999)
 293   double feop = matchclsE/esdtrack->P();
 294
 295         shower[0] = cluster->GetNCells(); // number of cells in cluster
 296         shower[1] = cluster->GetM02(); // long axis
 297         shower[2] = cluster->GetM20(); // short axis
 298         shower[3] = cluster->GetDispersion(); // dispersion
 299
 300   //   if(feop!=-9999.9)printf("%f\n",feop) ;
 301   TVector3 emcsignal(feop,0.0,rmatch1);
 302
 303   return emcsignal;
 304
 305 }