PWG1/AliComparisonDEdx.cxx

   1 //------------------------------------------------------------------------------
   2 // Implementation of AliComparisonDEdx class. It keeps information from
   3 // comparison of reconstructed and MC particle tracks. In addtion,
   4 // it keeps selection cuts used during comparison. The comparison
   5 // information is stored in the ROOT histograms. Analysis of these
   6 // histograms can be done by using Analyse() class function. The result of
   7 // the analysis (histograms) are stored in the output picture_dedx.root file.
   8 //
   9 // Author: J.Otwinowski 04/02/2008
  10 //------------------------------------------------------------------------------
  11
  12
  13 /*
  14   //after running analysis, read the file, and get component
  15   gSystem->Load("libPWG1.so");
  16   TFile f("Output.root");
  17   AliComparisonDEdx * comp = (AliComparisonDEdx*)f.Get("AliComparisonDEdx");
  18
  19   // analyse comparison data (output stored in pictures_dedx.root)
  20   comp->Analyse();
  21
  22   // TPC track length parameterisation
  23   TF1 fl("fl","((min(250./(abs(x+0.000001)),250)-90))",0,2);  // length function
  24   TF1 fl2("fl2","[0]/((min(250./(abs(x+0.000001)),250)-90))^[1]",0,2);
  25   fl2.SetParameter(1,1);
  26   fl2.SetParameter(0,1);
  27 */
  28
  29 #include <iostream>
  30
  31 #include "TFile.h"
  32 #include "TCint.h"
  33 #include "TH3F.h"
  34 #include "TH2F.h"
  35 #include "TF1.h"
  36 #include "TProfile.h"
  37 #include "TProfile2D.h"
  38 #include "TGraph2D.h"
  39 #include "TCanvas.h"
  40 #include "TGraph.h"
  41 //
  42 #include "AliESDEvent.h"
  43 #include "AliESD.h"
  44 #include "AliESDfriend.h"
  45 #include "AliESDfriendTrack.h"
  46 #include "AliRecInfoCuts.h"
  47 #include "AliMCInfoCuts.h"
  48 #include "AliLog.h"
  49 //
  50 #include "AliMathBase.h"
  51 #include "AliTreeDraw.h"
  52 //#include "TStatToolkit.h"
  53
  54 #include "AliMCInfo.h"
  55 #include "AliESDRecInfo.h"
  56 #include "AliComparisonDEdx.h"
  57
  58 using namespace std;
  59
  60 ClassImp(AliComparisonDEdx)
  61
  62 //_____________________________________________________________________________
  63 AliComparisonDEdx::AliComparisonDEdx():
  64   TNamed("AliComparisonDEdx","AliComparisonDEdx"),
  65
  66   // dEdx
  67   fTPCSignalNormTan(0),
  68   fTPCSignalNormSPhi(0),
  69   fTPCSignalNormTPhi(0),
  70   //
  71   fTPCSignalNormTanSPhi(0),
  72   fTPCSignalNormTanTPhi(0),
  73   fTPCSignalNormTanSPt(0),
  74
  75   // Cuts
  76   fCutsRC(0),
  77   fCutsMC(0),
  78   fMCPtMin(0),
  79   fMCAbsTanThetaMax(0),
  80   fMCPdgCode(0)
  81 {
  82   InitHisto();
  83   InitCuts();
  84 }
  85
  86 //_____________________________________________________________________________
  87 AliComparisonDEdx::~AliComparisonDEdx(){
  88
  89   if(fTPCSignalNormTan)  delete fTPCSignalNormTan; fTPCSignalNormTan=0;
  90   if(fTPCSignalNormSPhi) delete fTPCSignalNormSPhi; fTPCSignalNormSPhi=0;
  91   if(fTPCSignalNormTPhi) delete fTPCSignalNormTPhi; fTPCSignalNormTPhi=0;
  92   //
  93   if(fTPCSignalNormTanSPhi) delete fTPCSignalNormTanSPhi; fTPCSignalNormTanSPhi=0;
  94   if(fTPCSignalNormTanTPhi) delete fTPCSignalNormTanTPhi; fTPCSignalNormTanTPhi=0;
  95   if(fTPCSignalNormTanSPt)  delete fTPCSignalNormTanSPt; fTPCSignalNormTanSPt=0;
  96 }
  97
  98 //_____________________________________________________________________________
  99 void AliComparisonDEdx::InitHisto()
 100 {
 101   // Init histograms
 102
 103   // TPC dEdx
 104   fTPCSignalNormTan = new TH2F("CdEdxTan","CdEdxTan",50, -2,2,  40,30,70);
 105   fTPCSignalNormTan->SetXTitle("tan(#theta)");
 106   fTPCSignalNormTan->SetYTitle("rec. dE/dx / calc. dE/dx");
 107
 108   fTPCSignalNormSPhi   = new TH2F("CdEdxSPhi","CdEdxSPhi",10,0.0,1,40,30,70);
 109   fTPCSignalNormSPhi->SetXTitle("sin(#phi)");
 110   fTPCSignalNormSPhi->SetYTitle("rec. dE/dx / calc. dE/dx");
 111
 112   fTPCSignalNormTPhi   = new TH2F("CdEdxTPhi","CdEdxTPhi",10,0.0,2,40,30,70);
 113   fTPCSignalNormTPhi->SetXTitle("tan(#phi)");
 114   fTPCSignalNormTPhi->SetYTitle("rec. dE/dx / calc. dE/dx");
 115
 116   fTPCSignalNormTanSPhi= new TH3F("CdEdxTanSPhi","CdEdxTanSPhi",20, -2,2, 10,0.0 ,1,  40,30,70);
 117   fTPCSignalNormTanSPhi->SetXTitle("tan(#theta)");
 118   fTPCSignalNormTanSPhi->SetYTitle("sin(#phi)");
 119   fTPCSignalNormTanSPhi->SetZTitle("rec. dE/dx / calc. dE/dx");
 120
 121   fTPCSignalNormTanTPhi= new TH3F("CdEdxTanTPhi","CdEdxTanTPhi",20, -2,2, 10,0.0 ,1,  40,30,70);
 122   fTPCSignalNormTanTPhi->SetXTitle("tan(#theta)");
 123   fTPCSignalNormTanTPhi->SetYTitle("tan(#phi)");
 124   fTPCSignalNormTanTPhi->SetZTitle("rec. dE/dx / calc. dE/dx");
 125
 126   fTPCSignalNormTanSPt= new TH3F("CdEdxTanSPt","CdEdxTanSPt",20, -2,2, 10,0.3 ,3, 40,30,70);
 127   fTPCSignalNormTanSPt->SetXTitle("tan(#theta)");
 128   fTPCSignalNormTanSPt->SetYTitle("#sqrt(p_{t})");
 129   fTPCSignalNormTanSPt->SetZTitle("rec. dE/dx / calc. dE/dx");
 130 }
 131
 132 void AliComparisonDEdx::InitCuts()
 133 {
 134   // Init cuts
 135   if(!fCutsMC)
 136     AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");
 137   if(!fCutsRC)
 138     AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");
 139 }
 140
 141 //_____________________________________________________________________________
 142 void AliComparisonDEdx::Process(AliMCInfo* infoMC, AliESDRecInfo *infoRC){
 143
 144   // Fill dE/dx  comparison information
 145
 146   Float_t mcpt = infoMC->GetParticle().Pt();
 147   Float_t tantheta = TMath::Tan(infoMC->GetParticle().Theta()-TMath::Pi()*0.5);
 148   Bool_t isPrim = infoMC->GetParticle().R()<fCutsMC->GetMaxR() && TMath::Abs(infoMC->GetParticle().Vz())<fCutsMC->GetMaxVz() ;
 149
 150   Float_t mprim = infoMC->GetPrim();
 151
 152   // Check selection cuts
 153   if (fCutsMC->IsPdgParticle(TMath::Abs(infoMC->GetParticle().GetPdgCode())) == kFALSE) return;
 154   if (!isPrim) return;
 155   if (infoRC->GetStatus(1)!=3) return;
 156   if (!infoRC->GetESDtrack()) return;
 157   if (infoRC->GetESDtrack()->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return;
 158   if (!infoRC->GetESDtrack()->GetConstrainedParam()) return;
 159   //if (mprim>1.4) return;
 160   //if (mprim<0.5) return;
 161   if (mprim > fCutsMC->GetMaxTPCSignal()) return;
 162   if (mprim < fCutsMC->GetMinTPCSignal()) return;
 163   if (infoRC->GetESDtrack()->GetTPCsignalN()<fCutsRC->GetMinTPCsignalN()) return;
 164   //
 165   Float_t ratio = infoRC->GetESDtrack()->GetTPCsignal()/infoMC->GetPrim();
 166   Float_t sphi =  infoRC->GetESDtrack()->GetInnerParam()->GetSnp();
 167   Float_t tphi =  sphi/TMath::Sqrt(1-sphi*sphi);
 168
 169   if (TMath::Abs(infoMC->GetParticle().GetPdgCode()) != GetMCPdgCode()) return;
 170   //if (mcpt>0.5) {
 171   if (mcpt > GetMCPtMin()) {
 172     fTPCSignalNormTan->Fill(tantheta,ratio);    // only subset
 173   }
 174
 175   //if (TMath::Abs(tantheta)<0.5){
 176   if (TMath::Abs(tantheta) < GetMCAbsTanThetaMax()){
 177     fTPCSignalNormSPhi->Fill(sphi,ratio);       // only subset
 178     fTPCSignalNormTPhi->Fill(tphi,ratio);       // only subset
 179   }
 180   fTPCSignalNormTanSPhi->Fill(tantheta,sphi,ratio);
 181   fTPCSignalNormTanTPhi->Fill(tantheta,tphi,ratio);
 182   fTPCSignalNormTanSPt->Fill(tantheta,TMath::Sqrt(mcpt),ratio);
 183 }
 184
 185 //_____________________________________________________________________________
 186 Long64_t AliComparisonDEdx::Merge(TCollection* list)
 187 {
 188   // Merge list of objects (needed by PROOF)
 189
 190   if (!list)
 191   return 0;
 192
 193   if (list->IsEmpty())
 194   return 1;
 195
 196   TIterator* iter = list->MakeIterator();
 197   TObject* obj = 0;
 198
 199   // collection of generated histograms
 200   Int_t count=0;
 201   while((obj = iter->Next()) != 0)
 202   {
 203     AliComparisonDEdx* entry = dynamic_cast<AliComparisonDEdx*>(obj);
 204     if (entry == 0) {
 205       Error("Add","Attempt to add object of class: %s to a %s",
 206           obj->ClassName(),this->ClassName());
 207       return -1;
 208   }
 209
 210   fTPCSignalNormTan->Add(entry->fTPCSignalNormTan);
 211   fTPCSignalNormSPhi->Add(entry->fTPCSignalNormSPhi);
 212   fTPCSignalNormTPhi->Add(entry->fTPCSignalNormTPhi);
 213   //
 214   fTPCSignalNormTanSPhi->Add(entry->fTPCSignalNormTanSPhi);
 215   fTPCSignalNormTanTPhi->Add(entry->fTPCSignalNormTanTPhi);
 216   fTPCSignalNormTanSPt->Add(entry->fTPCSignalNormTanSPt);
 217
 218   count++;
 219   }
 220
 221 return count;
 222 }
 223
 224 //_____________________________________________________________________________
 225 void AliComparisonDEdx::Exec(AliMCInfo* infoMC, AliESDRecInfo *infoRC)
 226 {
 227   // Process comparison information
 228   Process(infoMC,infoRC);
 229 }
 230
 231 //_____________________________________________________________________________
 232 TH1F* AliComparisonDEdx::MakeResol(TH2F * his, Int_t integ, Bool_t type)
 233 {
 234   // Make resolution histograms
 235   TH1F *hisr, *hism;
 236   if (!gPad) new TCanvas;
 237   hisr = AliTreeDraw::CreateResHistoI(his,&hism,integ);
 238   if (type) return hism;
 239   else
 240     return hisr;
 241 }
 242
 243 //_____________________________________________________________________________
 244 void AliComparisonDEdx::Analyse()
 245 {
 246   // Analyse output histograms
 247
 248   AliComparisonDEdx * comp=this;
 249   TH1F *hiss=0;
 250   TGraph2D * gr=0;
 251
 252   TFile *fp = new TFile("pictures_dedx.root","recreate");
 253   fp->cd();
 254
 255   TCanvas * c = new TCanvas("TPCdedx","TPC dedx");
 256   c->cd();
 257
 258   hiss = comp->MakeResol(comp->fTPCSignalNormTan,4,0);
 259   hiss->SetXTitle("Tan(#theta)");
 260   hiss->SetYTitle("#sigma_{dEdx}");
 261   hiss->Draw();
 262   hiss->Write("TPCdEdxResolTan");
 263   //
 264   hiss = comp->MakeResol(comp->fTPCSignalNormTan,4,1);
 265   hiss->SetXTitle("Tan(#theta)");
 266   hiss->SetYTitle("<dEdx>");
 267   hiss->Draw();
 268   hiss->Write("TPCdEdxMeanTan");
 269   //
 270   gr = AliMathBase::MakeStat2D(comp->fTPCSignalNormTanSPt,3,1,4);
 271   gr->GetXaxis()->SetTitle("Tan(#theta)");
 272   gr->GetYaxis()->SetTitle("#sqrt{p_{t}(GeV)}");
 273   gr->GetZaxis()->SetTitle("<dEdx>");
 274   gr->Draw("colz");
 275   gr->GetHistogram()->Write("TPCdEdxMeanTanPt_1");
 276   //
 277   gr = AliMathBase::MakeStat2D(comp->fTPCSignalNormTanSPt,3,1,5);
 278   gr->GetXaxis()->SetTitle("Tan(#theta)");
 279   gr->GetYaxis()->SetTitle("#sqrt{p_{t}(GeV)}");
 280   gr->GetZaxis()->SetTitle("#sigma_{dEdx}");
 281   gr->Draw("colz");
 282   gr->GetHistogram()->Write("TPCdEdxMeanTanPt_2");
 283
 284   fp->Close();
 285 }