ITS/AliITSdEdxAnalyzer.cxx

   1 /**************************************************************************
   2  * Copyright(c) 2007-2009, 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 ///////////////////////////////////////////////////////////////////
  19 //                                                               //
  20 // Implementation of the base class for dEdx analysis            //
  21 // Origin: F.Prino, Torino, prino@to.infn.it                     //
  22 //                                                               //
  23 ///////////////////////////////////////////////////////////////////
  24
  25 #include <TString.h>
  26 #include <TFile.h>
  27 #include <TParticlePDG.h>
  28 #include <TDatabasePDG.h>
  29 #include "AliStack.h"
  30 #include "AliLog.h"
  31 #include "AliESDEvent.h"
  32 #include "AliStack.h"
  33 #include "AliITSdEdxAnalyzer.h"
  34 #include "AliExternalTrackParam.h"
  35 //#include "AliITSpidESD.h"
  36 #include "AliITSPIDResponse.h"
  37 #include "AliPID.h"
  38
  39 ClassImp(AliITSdEdxAnalyzer)
  40
  41 const Int_t AliITSdEdxAnalyzer::fgkPdgCode[kNParticles]={211,321,2212};
  42 const Int_t AliITSdEdxAnalyzer::fgkLayerCode[kNValuesdEdx]={3,4,5,6,0};
  43
  44 //______________________________________________________________________
  45 AliITSdEdxAnalyzer::AliITSdEdxAnalyzer() :
  46   TObject(),
  47   fNPBins(10),
  48   fPMin(0.1),
  49   fPMax(1.1),
  50   fHistodEdx(0),
  51   fHistoDeltadEdx(0),
  52   fHistodEdxVsP(0),
  53   fThickness(0.03),
  54   fDensity(2.33),
  55   fBBmodel(0),
  56   fMIP(82.),
  57   fTPCpidCut(0.5)
  58 {
  59   // default constructor
  60   BookHistos();
  61 }
  62 //______________________________________________________________________
  63 AliITSdEdxAnalyzer::AliITSdEdxAnalyzer(const Int_t npbins, const Float_t pmin, const Float_t pmax) :
  64   TObject(),
  65   fNPBins(npbins),
  66   fPMin(pmin),
  67   fPMax(pmax),
  68   fHistodEdx(0),
  69   fHistoDeltadEdx(0),
  70   fHistodEdxVsP(0),
  71   fThickness(0.03),
  72   fDensity(2.33),
  73   fBBmodel(0),
  74   fMIP(82.),
  75   fTPCpidCut(0.5)
  76 {
  77   // standard constructor
  78   BookHistos();
  79 }
  80 //______________________________________________________________________
  81 AliITSdEdxAnalyzer::~AliITSdEdxAnalyzer(){
  82   // destructor
  83   DeleteHistos();
  84 }
  85 //______________________________________________________________________
  86 void AliITSdEdxAnalyzer::SetMomentumBins(const Int_t npbins, const Float_t pmin, const Float_t pmax){
  87   // Kill exisiting histos, set new binning, book new histos
  88   DeleteHistos();
  89   fNPBins=npbins;
  90   fPMin=pmin;
  91   fPMax=pmax;
  92   BookHistos();
  93 }
  94 //______________________________________________________________________
  95 void AliITSdEdxAnalyzer::DeleteHistos(){
  96   // deletes the hitograms
  97   if(fHistodEdx){
  98     for(Int_t i=0; i<GetNHistos();i++) delete fHistodEdx[i];
  99     delete [] fHistodEdx;
 100     fHistodEdx=0;
 101   }
 102   if(fHistoDeltadEdx){
 103     for(Int_t i=0; i<GetNHistos();i++) delete fHistoDeltadEdx[i];
 104     delete [] fHistoDeltadEdx;
 105     fHistoDeltadEdx=0;
 106   }
 107   if(fHistodEdxVsP){
 108     for(Int_t i=0; i<GetNHistos2();i++) delete fHistodEdxVsP[i];
 109     delete [] fHistodEdxVsP;
 110     fHistodEdxVsP=0;
 111   }
 112 }
 113 //______________________________________________________________________
 114 void AliITSdEdxAnalyzer::BookHistos(){
 115   // defines the output histograms
 116   fHistodEdx=new TH1F*[GetNHistos()];
 117   fHistoDeltadEdx=new TH1F*[GetNHistos()];
 118   for(Int_t iSpecie=0; iSpecie<kNParticles; iSpecie++){
 119     for(Int_t iPBin=0; iPBin<fNPBins; iPBin++){
 120       for(Int_t iVal=0; iVal<kNValuesdEdx; iVal++){
 121         TString hisnam1=Form("hdEdx%dpbin%dpart%d",fgkLayerCode[iVal],iPBin,fgkPdgCode[iSpecie]);
 122         TString histit1=Form("dEdx layer %d (keV) pbin%d part%d",fgkLayerCode[iVal],iPBin,fgkPdgCode[iSpecie]);
 123         if(iVal==kNValuesdEdx-1) histit1=Form("dEdx trunc. mean (keV) pbin%d part%d",iPBin,fgkPdgCode[iSpecie]);
 124         fHistodEdx[GetIndex(iSpecie,iPBin,iVal)]=new TH1F(hisnam1.Data(),histit1.Data(),100,0.,500.);
 125         fHistodEdx[GetIndex(iSpecie,iPBin,iVal)]->GetXaxis()->SetTitle(histit1.Data());
 126         TString hisnam2=Form("hdeltadEdx%dpbin%dpart%d",fgkLayerCode[iVal],iPBin,fgkPdgCode[iSpecie]);
 127         TString histit2=Form("(dEdx_l%d-BetheBloch)/BetheBloch pbin%d part%d",fgkLayerCode[iVal],iPBin,fgkPdgCode[iSpecie]);
 128         if(iVal==kNValuesdEdx-1) histit1=Form("(dEdx_truncmean-BetheBloch)/BetheBloch pbin%d part%d",iPBin,fgkPdgCode[iSpecie]);
 129         fHistoDeltadEdx[GetIndex(iSpecie,iPBin,iVal)]=new TH1F(hisnam2.Data(),histit2.Data(),100,-0.5,0.5);
 130         fHistoDeltadEdx[GetIndex(iSpecie,iPBin,iVal)]->GetXaxis()->SetTitle(histit2.Data());
 131       }
 132     }
 133   }
 134
 135   fHistodEdxVsP=new TH2F*[GetNHistos2()];
 136   for(Int_t iSpecie=0; iSpecie<kNParticles; iSpecie++){
 137     for(Int_t iVal=0; iVal<kNValuesdEdx; iVal++){
 138       TString hisnam=Form("hdEdx%dVsPpart%d",fgkLayerCode[iVal],fgkPdgCode[iSpecie]);
 139       TString histit=Form("dEdx layer %d (keV) vs P part%d",fgkLayerCode[iVal],fgkPdgCode[iSpecie]);
 140       if(iVal==kNValuesdEdx-1) histit=Form("dEdx trunc. mean (keV) vs P part%d",fgkPdgCode[iSpecie]);
 141       fHistodEdxVsP[GetIndex2(iSpecie,iVal)]=new TH2F(hisnam.Data(),histit.Data(),50,fPMin,fPMax,50,0.,500.);
 142       histit.ReplaceAll(" vs P "," ");
 143       fHistodEdxVsP[GetIndex2(iSpecie,iVal)]->GetYaxis()->SetTitle(histit.Data());
 144       fHistodEdxVsP[GetIndex2(iSpecie,iVal)]->GetXaxis()->SetTitle("Momentum (GeV/c)");
 145     }
 146   }
 147 }
 148 //______________________________________________________________________
 149 void AliITSdEdxAnalyzer:: WriteHistos(TString filename) const {
 150   // stores output histogrma to file
 151   TFile *outfile=new TFile(filename.Data(),"recreate");
 152   outfile->cd();
 153   for(Int_t i=0; i<GetNHistos();i++){
 154     fHistodEdx[i]->Write();
 155     fHistoDeltadEdx[i]->Write();
 156   }
 157   for(Int_t i=0; i<GetNHistos2();i++) fHistodEdxVsP[i]->Write();
 158   outfile->Close();
 159   delete outfile;
 160 }
 161 //______________________________________________________________________
 162 void AliITSdEdxAnalyzer::ReadEvent(const AliESDEvent* esd, AliStack* stack){
 163   // Fill histos
 164   // if stack is !=0 MC truth is used to define particle specie
 165   // else TPC pid is used to define particle specie
 166
 167   if(!esd) AliFatal("Bad ESD event");
 168   for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
 169     AliESDtrack* track = esd->GetTrack(iTrack);
 170     Double_t trmean=track->GetITSsignal();
 171     Double_t sig[4];
 172     track->GetITSdEdxSamples(sig);
 173     Double_t preco=track->P();
 174     Int_t iPBin=GetMomentumBin(preco);
 175     if(iPBin==-1) continue;
 176     Int_t iSpecie=-1;
 177     Double_t dedxbb=0;
 178     if(stack){
 179       Int_t lab=track->GetLabel();
 180       if(lab<0) continue;
 181       TParticle* part=(TParticle*)stack->Particle(lab);
 182       Int_t absPdgCode=TMath::Abs(part->GetPdgCode());
 183       iSpecie=GetSpecieBin(absPdgCode);
 184       if(iSpecie==-1) continue;
 185       dedxbb=BetheBloch(part);
 186     }else{
 187       iSpecie=GetPaticleIdFromTPC(track);
 188       if(iSpecie==-1) continue;
 189       TParticlePDG* p=TDatabasePDG::Instance()->GetParticle(fgkPdgCode[iSpecie]);
 190       dedxbb=BetheBloch(preco,p->Mass());
 191     }
 192     for(Int_t ilay=0;ilay<4;ilay++){
 193       if(sig[ilay]>0.){
 194         fHistodEdx[GetIndex(iSpecie,iPBin,ilay)]->Fill(sig[ilay]);
 195         fHistoDeltadEdx[GetIndex(iSpecie,iPBin,ilay)]->Fill((sig[ilay]-dedxbb)/dedxbb);
 196         fHistodEdxVsP[GetIndex2(iSpecie,ilay)]->Fill(preco,sig[ilay]);
 197       }
 198     }
 199     if(trmean>0.){
 200       fHistodEdx[GetIndex(iSpecie,iPBin,4)]->Fill(trmean);
 201       fHistoDeltadEdx[GetIndex(iSpecie,iPBin,4)]->Fill((trmean-dedxbb)/dedxbb);
 202       fHistodEdxVsP[GetIndex2(iSpecie,4)]->Fill(preco,trmean);
 203     }
 204   }
 205 }
 206 //______________________________________________________________________
 207 Int_t AliITSdEdxAnalyzer::GetPaticleIdFromTPC(const AliESDtrack* track) const {
 208   // Returns the particle specie as identified by TPC
 209   Double_t tpcpid[AliPID::kSPECIES];
 210   track->GetTPCpid(tpcpid);
 211   Int_t maxPos=-1;
 212   Float_t maxProb=0.;
 213   for(Int_t is=0; is<AliPID::kSPECIES; is++){
 214     if(tpcpid[is]>maxProb){
 215       maxProb=tpcpid[is];
 216       maxPos=is;
 217     }
 218   }
 219   Int_t iSpecie=-1;
 220   if(maxProb>fTPCpidCut){
 221     if(maxPos==AliPID::kPion) iSpecie=0;
 222     else if(maxPos==AliPID::kKaon) iSpecie=1;
 223     else if(maxPos==AliPID::kProton) iSpecie=2;
 224   }
 225   return iSpecie;
 226 }
 227 //______________________________________________________________________
 228 Double_t AliITSdEdxAnalyzer::BetheBloch(const Float_t p, const Float_t m) const {
 229   // Computes expected dE/dx value for given particle specie and momentum
 230   static AliITSPIDResponse pidResponse;
 231   Double_t dedxbb=0.;
 232   if(fBBmodel==0){
 233     Double_t betagamma=p/m;
 234     Double_t conv=fDensity*1E6*fThickness/116.24*fMIP;
 235     dedxbb=conv*AliExternalTrackParam::BetheBlochSolid(betagamma);
 236   }else if(fBBmodel==1){
 237     dedxbb=pidResponse.Bethe(p,m);
 238   }
 239   return dedxbb;
 240 }
 241 //______________________________________________________________________
 242 TGraph* AliITSdEdxAnalyzer::GetBetheBlochGraph(const Int_t pdgcode) const {
 243  // Fills a TGraph with Bethe-Bloch expe
 244  TGraph* g=new TGraph(0);
 245  TParticlePDG* part=TDatabasePDG::Instance()->GetParticle(pdgcode);
 246  Float_t mass=part->Mass();
 247  for(Int_t ip=0; ip<100;ip++){
 248    Float_t p=fPMin+(ip+0.5)*(fPMax-fPMin)/100.;
 249    g->SetPoint(ip,p,BetheBloch(p,mass));
 250  }
 251  g->GetXaxis()->SetTitle("Momentum (GeV/c)");
 252  g->GetYaxis()->SetTitle("dE/dx (keV/300 #mum)");
 253  return g;
 254 }