Adding cuts for comparison studies (Marian)
[u/mrichter/AliRoot.git] / 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