PWG1/AliComparisonRes.cxx

   1 //------------------------------------------------------------------------------
   2 // Implementation of AliComparisonRes 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_res.root file.
   8 //
   9 // Author: J.Otwinowski 04/02/2008
  10 //------------------------------------------------------------------------------
  11
  12 /*
  13   //after running analysis, read the file, and get component
  14   gSystem->Load("libPWG1.so");
  15   TFile f("Output.root");
  16   AliComparisonRes * comp = (AliComparisonRes*)f.Get("AliComparisonRes");
  17
  18   // analyse comparison data (output stored in pictures_res.root)
  19   comp->Analyse();
  20
  21   // paramtetrisation of the TPC track length (for information only)
  22   TF1 fl("fl","((min(250./(abs(x+0.000001)),250)-90))",0,2);  // TPC track length function
  23   TF1 fl2("fl2","[0]/((min(250./(abs(x+0.000001)),250)-90))^[1]",0,2);
  24   fl2.SetParameter(1,1);
  25   fl2.SetParameter(0,1);
  26 */
  27
  28 #include <iostream>
  29
  30 #include "TFile.h"
  31 #include "TCint.h"
  32 #include "TH3F.h"
  33 #include "TH2F.h"
  34 #include "TF1.h"
  35 #include "TProfile.h"
  36 #include "TProfile2D.h"
  37 #include "TGraph2D.h"
  38 #include "TCanvas.h"
  39 #include "TGraph.h"
  40
  41 #include "AliESDEvent.h"
  42 #include "AliESD.h"
  43 #include "AliESDfriend.h"
  44 #include "AliESDfriendTrack.h"
  45 #include "AliRecInfoCuts.h"
  46 #include "AliMCInfoCuts.h"
  47 #include "AliLog.h"
  48
  49 #include "AliMathBase.h"
  50 #include "AliTreeDraw.h"
  51
  52 #include "AliMCInfo.h"
  53 #include "AliESDRecInfo.h"
  54 #include "AliComparisonRes.h"
  55
  56 using namespace std;
  57
  58 ClassImp(AliComparisonRes)
  59
  60 //_____________________________________________________________________________
  61 AliComparisonRes::AliComparisonRes():
  62   TNamed("AliComparisonRes","AliComparisonRes"),
  63
  64   // Resolution
  65   fPtResolLPT(0),        // pt resolution - low pt
  66   fPtResolHPT(0),        // pt resolution - high pt
  67   fPtPullLPT(0),         // pt resolution - low pt
  68   fPtPullHPT(0),         // pt resolution - high pt
  69   //
  70   // Resolution constrained param
  71   //
  72   fCPhiResolTan(0),  // angular resolution -  constrained
  73   fCTanResolTan(0),  // angular resolution -  constrained
  74   fCPtResolTan(0),   // pt resolution      -  constrained
  75   fCPhiPullTan(0),   // angular resolution -  constrained
  76   fCTanPullTan(0),   // angular resolution -  constrained
  77   fCPtPullTan(0),    // pt resolution      -  constrained
  78
  79   // Cuts
  80   fCutsRC(0),
  81   fCutsMC(0)
  82 {
  83   InitHisto();
  84   InitCuts();
  85 }
  86
  87 //_____________________________________________________________________________
  88 AliComparisonRes::~AliComparisonRes(){
  89
  90   // Resolution histograms
  91   if(fPtResolLPT) delete  fPtResolLPT; fPtResolLPT=0;
  92   if(fPtResolHPT) delete  fPtResolHPT; fPtResolHPT=0;
  93   if(fPtPullLPT)  delete  fPtPullLPT;  fPtPullLPT=0;
  94   if(fPtPullHPT)  delete  fPtPullHPT;  fPtPullHPT=0;
  95
  96   // Resolution histograms (constrained param)
  97   if(fCPhiResolTan) delete fCPhiResolTan; fCPhiResolTan=0;
  98   if(fCTanResolTan) delete fCTanResolTan; fCTanResolTan=0;
  99   if(fCPtResolTan)  delete fCPtResolTan;  fCPtResolTan=0;
 100   if(fCPhiPullTan)  delete fCPhiPullTan;  fCPhiPullTan=0;
 101   if(fCTanPullTan)  delete fCTanPullTan;  fCTanPullTan=0;
 102   if(fCPtPullTan)   delete fCPtPullTan;   fCPtPullTan=0;
 103 }
 104
 105 //_____________________________________________________________________________
 106 void AliComparisonRes::InitHisto(){
 107
 108   // Init histograms
 109   fCPhiResolTan = new TH2F("CPhiResolTan","CPhiResolTan",50, -2,2,200,-0.025,0.025);
 110   fCPhiResolTan->SetXTitle("tan(#theta)");
 111   fCPhiResolTan->SetYTitle("#Delta#phi");
 112
 113   fCTanResolTan = new TH2F("CTanResolTan","CTanResolTan",50, -2,2,200,-0.025,0.025);
 114   fCTanResolTan->SetXTitle("tan(#theta)");
 115   fCTanResolTan->SetYTitle("#Delta#theta");
 116
 117   fCPtResolTan=new TH2F("CPtResol","CPtResol",50, -2,2,200,-0.2,0.2);
 118   fCPtResolTan->SetXTitle("Tan(#theta)");
 119   fCPtResolTan->SetYTitle("#Deltap_{t}/p_{t}");
 120
 121   fCPhiPullTan = new TH2F("CPhiPullTan","CPhiPullTan",50, -2,2,200,-5,5);
 122   fCPhiPullTan->SetXTitle("Tan(#theta)");
 123   fCPhiPullTan->SetYTitle("#Delta#phi/#Sigma");
 124
 125   fCTanPullTan = new TH2F("CTanPullTan","CTanPullTan",50, -2,2,200,-5,5);
 126   fCTanPullTan->SetXTitle("Tan(#theta)");
 127   fCTanPullTan->SetYTitle("#Delta#theta/#Sigma");
 128
 129   fCPtPullTan=new TH2F("CPtPull","CPtPull",50, -2,2,200,-5,5);
 130   fCPtPullTan->SetXTitle("Tan(#theta)");
 131   fCPtPullTan->SetYTitle("(1/mcp_{t}-1/p_{t})/#Sigma");
 132
 133   fPtResolLPT = new TH2F("Pt_resol_lpt","pt resol",10, 0.1,3,200,-0.2,0.2);
 134   fPtResolLPT->SetXTitle("p_{t}");
 135   fPtResolLPT->SetYTitle("#Deltap_{t}/p_{t}");
 136
 137   fPtResolHPT = new TH2F("Pt_resol_hpt","pt resol",10, 2,100,200,-0.3,0.3);
 138   fPtResolHPT->SetXTitle("p_{t}");
 139   fPtResolHPT->SetYTitle("#Deltap_{t}/p_{t}");
 140
 141   fPtPullLPT = new TH2F("Pt_pull_lpt","pt pool",10, 0.1,3,200,-6,6);
 142   fPtPullLPT->SetXTitle("p_{t}");
 143   fPtPullLPT->SetYTitle("#Deltap_{t}/#Sigma");
 144
 145   fPtPullHPT = new TH2F("Pt_pull_hpt","pt pool",10, 2,100,200,-6,6);
 146   fPtPullHPT->SetXTitle("p_{t}");
 147   fPtPullHPT->SetYTitle("#Deltap_{t}/#Sigma");
 148 }
 149
 150 //_____________________________________________________________________________
 151 void AliComparisonRes::InitCuts()
 152 {
 153
 154  /*
 155     fCutsRC = new AliRecInfoCuts();
 156   if(fCutsRC) {
 157     fCutsRC->SetPtRange(0.15,200.0);
 158     fCutsRC->SetMaxAbsTanTheta(1.0);
 159     fCutsRC->SetMinNClustersTPC(10);
 160     fCutsRC->SetMinTPCsignalN(50);
 161   } else {
 162     AliDebug(AliLog::kError, "ERROR: Cannot create AliRecInfoCuts object");
 163     return;
 164   }
 165
 166   fCutsMC = new AliMCInfoCuts();
 167   if(fCutsMC) {
 168     fCutsMC->SetMinRowsWithDigits(50);
 169     fCutsMC->SetMaxR(0.1);
 170     fCutsMC->SetMaxVz(10.);
 171     fCutsMC->SetRangeTPCSignal(0.5,1.4);
 172   } else {
 173   AliDebug(AliLog::kError, "ERROR: Cannot AliMCInfoCuts object");
 174   return;
 175   }
 176  */
 177   // Init cuts
 178   if(!fCutsMC)
 179     AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");
 180   if(!fCutsRC)
 181     AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");
 182 }
 183
 184 //_____________________________________________________________________________
 185 void AliComparisonRes::Process(AliMCInfo* infoMC, AliESDRecInfo *infoRC)
 186 {
 187   // Fill resolution comparison information
 188
 189   Float_t mcpt = infoMC->GetParticle().Pt();
 190   Bool_t isPrim = infoMC->GetParticle().R() < fCutsMC->GetMaxR() &&
 191                   TMath::Abs(infoMC->GetParticle().Vz()) < fCutsMC->GetMaxVz() ;
 192
 193   // Check selection cuts
 194   if (fCutsMC->IsPdgParticle(TMath::Abs(infoMC->GetParticle().GetPdgCode())) == kFALSE) return;
 195   if (!isPrim) return;
 196   if (infoRC->GetStatus(1)==0) return;
 197   if (!infoRC->GetESDtrack()) return;
 198   if (infoRC->GetESDtrack()->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return;
 199
 200   Float_t deltaPt= (mcpt-infoRC->GetESDtrack()->Pt())/mcpt;
 201   Float_t pullPt= (1/mcpt-infoRC->GetESDtrack()->OneOverPt())/TMath::Sqrt(infoRC->GetESDtrack()->GetSigma1Pt2());
 202
 203   // Fill histograms
 204   fPtResolLPT->Fill(mcpt,deltaPt);
 205   fPtResolHPT->Fill(mcpt,deltaPt);
 206   fPtPullLPT->Fill(mcpt,pullPt);
 207   fPtPullHPT->Fill(mcpt,pullPt);
 208 }
 209
 210 Long64_t AliComparisonRes::Merge(TCollection* list)
 211 {
 212   // Merge list of objects (needed by PROOF)
 213
 214   if (!list)
 215   return 0;
 216
 217   if (list->IsEmpty())
 218   return 1;
 219
 220   TIterator* iter = list->MakeIterator();
 221   TObject* obj = 0;
 222
 223   // collection of generated histograms
 224   Int_t count=0;
 225   while((obj = iter->Next()) != 0)
 226   {
 227   AliComparisonRes* entry = dynamic_cast<AliComparisonRes*>(obj);
 228   if (entry == 0) {
 229     Error("Add","Attempt to add object of class: %s to a %s",
 230         obj->ClassName(),this->ClassName());
 231     return -1;
 232   }
 233
 234   fPtResolLPT->Add(entry->fPtResolLPT);
 235   fPtResolHPT->Add(entry->fPtResolHPT);
 236   fPtPullLPT->Add(entry->fPtPullLPT);
 237   fPtPullHPT->Add(entry->fPtPullHPT);
 238
 239   // Resolution histograms (constrained param)
 240   fCPhiResolTan->Add(entry->fCPhiResolTan);
 241   fCTanResolTan->Add(entry->fCTanResolTan);
 242   fCPtResolTan->Add(entry->fCPtResolTan);
 243   fCPhiPullTan->Add(entry->fCPhiPullTan);
 244   fCTanPullTan->Add(entry->fCTanPullTan);
 245   fCPtPullTan->Add(entry->fCPtPullTan);
 246
 247   count++;
 248   }
 249
 250 return count;
 251 }
 252
 253
 254
 255
 256
 257 //_____________________________________________________________________________
 258 void AliComparisonRes::ProcessConstrained(AliMCInfo* infoMC, AliESDRecInfo *infoRC)
 259 {
 260   // Fill resolution comparison information (constarained parameters)
 261
 262   Float_t mcpt = infoMC->GetParticle().Pt();
 263   Float_t tantheta = TMath::Tan(infoMC->GetParticle().Theta()-TMath::Pi()*0.5);
 264   Bool_t isPrim = infoMC->GetParticle().R()<fCutsMC->GetMaxR() && TMath::Abs(infoMC->GetParticle().Vz())<fCutsMC->GetMaxVz();
 265
 266   // Check selection cuts
 267   if (fCutsMC->IsPdgParticle(TMath::Abs(infoMC->GetParticle().GetPdgCode())) == kFALSE) return;
 268   if (!isPrim) return;
 269   if (infoRC->GetStatus(1)!=3) return;
 270   if (!infoRC->GetESDtrack()) return;
 271   if (infoRC->GetESDtrack()->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return;
 272   if (!infoRC->GetESDtrack()->GetConstrainedParam()) return;
 273
 274   // constrained parameters resolution
 275   const AliExternalTrackParam * cparam = infoRC->GetESDtrack()->GetConstrainedParam();
 276   Float_t deltaCPt= (mcpt-cparam->Pt())/mcpt;
 277   Float_t pullCPt= (1/mcpt-cparam->OneOverPt())/
 278     TMath::Sqrt(cparam->GetSigma1Pt2());
 279   Float_t deltaPhi = TMath::ATan2(cparam->Py(),cparam->Px())-
 280     TMath::ATan2(infoMC->GetParticle().Py(),infoMC->GetParticle().Px());
 281   Float_t pullPhi = deltaPhi/TMath::Sqrt(cparam->GetSigmaSnp2());
 282
 283   Float_t deltaTan = TMath::ATan2(cparam->Pz(),cparam->Pt())-TMath::ATan2(infoMC->GetParticle().Pz(),infoMC->GetParticle().Pt());
 284   Float_t pullTan = deltaPhi/TMath::Sqrt(cparam->GetSigmaSnp2());
 285
 286   // Fill histograms
 287   fCPtResolTan->Fill(tantheta,deltaCPt);
 288   fCPtPullTan->Fill(tantheta,pullCPt);
 289   fCPhiResolTan->Fill(tantheta,deltaPhi);
 290   fCPhiPullTan->Fill(tantheta,pullPhi);
 291   fCTanResolTan->Fill(tantheta,deltaTan);
 292   fCTanPullTan->Fill(tantheta,pullTan);
 293 }
 294
 295 //_____________________________________________________________________________
 296 void AliComparisonRes::Exec(AliMCInfo* infoMC, AliESDRecInfo *infoRC){
 297
 298   // Process comparison information
 299   Process(infoMC,infoRC);
 300   ProcessConstrained(infoMC,infoRC);
 301
 302 }
 303
 304 //_____________________________________________________________________________
 305 TH1F* AliComparisonRes::MakeResol(TH2F * his, Int_t integ, Bool_t type){
 306   // Create resolution histograms
 307
 308   TH1F *hisr, *hism;
 309   if (!gPad) new TCanvas;
 310   hisr = AliTreeDraw::CreateResHistoI(his,&hism,integ);
 311   if (type) return hism;
 312   else
 313     return hisr;
 314 }
 315
 316 //_____________________________________________________________________________
 317 void AliComparisonRes::Analyse(){
 318   // Analyse comparison information and store output histograms
 319   // in the "pictures_res.root" file
 320
 321   AliComparisonRes * comp=this;
 322   TH1F *hiss=0;
 323
 324   TFile *fp = new TFile("pictures_res.root","recreate");
 325   fp->cd();
 326
 327   TCanvas * c = new TCanvas("Phi resol Tan","Phi resol Tan");
 328   c->cd();
 329   //
 330   hiss = comp->MakeResol(comp->fCPtResolTan,1,0);
 331   hiss->SetXTitle("Tan(#theta)");
 332   hiss->SetYTitle("#sigmap_{t}/p_{t}");
 333   hiss->Draw();
 334   hiss->Write("CptResolTan");
 335   //
 336   hiss = comp->MakeResol(comp->fCPhiResolTan,1,0);
 337   hiss->SetXTitle("Tan(#theta)");
 338   hiss->SetYTitle("#sigma#phi (rad)");
 339   hiss->Draw();
 340   hiss->Write("PhiResolTan");
 341   //
 342   hiss = comp->MakeResol(comp->fCTanResolTan,1,0);
 343   hiss->SetXTitle("Tan(#theta)");
 344   hiss->SetYTitle("#sigma#theta (rad)");
 345   hiss->Draw();
 346   hiss->Write("ThetaResolTan");
 347   //
 348   hiss = comp->MakeResol(comp->fCPtPullTan,1,0);
 349   hiss->SetXTitle("Tan(#theta)");
 350   hiss->SetYTitle("1/mcp_{t}-1/p_{t}/#Sigma(1/p_{t})");
 351   hiss->Draw();
 352   hiss->Write("CptPullTan");
 353
 354   fp->Close();
 355 }
 356
 357
 358