When calculating a*a-b*b the form (a-b)*(a+b) is usually more numerically stable.
[u/mrichter/AliRoot.git] / PWG1 / AliComparisonDraw.cxx
1
2
3 //
4 // Comparison draw
5 // Compare the MC information with the reconstructed 
6 //
7
8 /*
9   after running analysis, read the file, and get component
10   gSystem->Load("libPWG1.so");
11   TFile f("Output.root");
12   AliComparisonDraw * comp = (AliComparisonDraw*)f.Get("AliComparisonDraw");
13   TF1 fl("fl","((min(250./(abs(x+0.000001)),250)-90))",0,2);  // length function
14   TF1 fl2("fl2","[0]/((min(250./(abs(x+0.000001)),250)-90))^[1]",0,2);
15   fl2.SetParameter(1,1);
16   fl2.SetParameter(0,1);
17
18 */
19
20
21
22
23 #include "TFile.h"
24 #include "TCint.h"
25 #include "TH3F.h"
26 #include "TH2F.h"
27 #include "TF1.h"
28 #include "TProfile.h"
29 #include "TProfile2D.h"
30 #include "TGraph2D.h"
31 #include "TCanvas.h"
32 #include "TGraph.h"
33 //
34 // 
35 #include "AliESDEvent.h"   // new container
36 #include "AliESD.h"
37 #include "AliESDfriend.h"
38 #include "AliESDfriendTrack.h"
39 //
40 #include "AliMathBase.h"
41 #include "AliTreeDraw.h" 
42
43 #include "AliMCInfo.h" 
44 #include "AliESDRecInfo.h" 
45 #include "AliComparisonDraw.h" 
46
47
48 ClassImp(AliComparisonDraw)
49
50 Bool_t    AliComparisonDraw::fgBDraw=kFALSE;         //option draw temporary results
51
52 AliComparisonDraw::AliComparisonDraw():
53   TNamed("ComparisonDraw","ComparisonDraw"),
54   fEffTPCPt(0),      // TPC efficiency as function of Pt (tan+-1)
55   fEffTPCPtMC(0),    // MC -TPC efficiency as function of Pt (tan+-1)
56   fEffTPCPtF(0),     // efficiency for findable tracks
57   //
58   fEffTPCTan(0),   // TPC efficiency as function of Tan (pt>0.15
59   fEffTPCTanMC(0), // MC -TPC efficiency as function of Tan (pt>0.15)
60   fEffTPCTanF(0),  // efficiency for findable tracks Tan (pt>0.15)
61   //
62   fEffTPCPtTan(0),    // TPC efficiency as function of Pt and tan
63   fEffTPCPtTanMC(0),  // MC -TPC efficiency as function of Pt and tan
64   fEffTPCPtTanF(0),  // TPC efficiency as function of Pt and tan
65   //
66   // dEdx resolution
67   //
68   fTPCSignalNormTan(0), // tpc signal normalized to the mean signal - MC
69   fTPCSignalNormSPhi(0),   // tpc signal normalized to the mean signal - MC
70   fTPCSignalNormTPhi(0),   // tpc signal normalized to the mean signal - MC
71   //
72   fTPCSignalNormTanSPhi(0),   // tpc signal normalized to the mean signal - MC
73   fTPCSignalNormTanTPhi(0),   // tpc signal normalized to the mean signal - MC
74   fTPCSignalNormTanSPt(0),   // tpc signal normalized to the mean signal - MC
75   //
76   //
77   fPtResolLPT(0),        // pt resolution - low pt
78   fPtResolHPT(0),        // pt resolution - high pt 
79   fPtPullLPT(0),         // pt resolution - low pt
80   fPtPullHPT(0),         // pt resolution - high pt 
81   //
82   // Resolution constrained param
83   //
84   fCPhiResolTan(0),   // angular resolution -  constrained
85   fCTanResolTan(0),   // angular resolution -  constrained
86   fCPtResolTan(0),    // pt resolution      -  constrained
87   fCPhiPullTan(0),   // angular resolution -  constrained
88   fCTanPullTan(0),   // angular resolution -  constrained
89   fCPtPullTan(0),    // pt resolution      -  constrained
90   //
91   // DCA resolution
92   //
93   fD0TanSPtB1(0),   // distance to vertex y  
94   fD1TanSPtB1(0),   // distance to vertex z  
95   fD0TanSPtL1(0),   // distance to vertex y  
96   fD1TanSPtL1(0)   // distance to vertex z  
97 {
98   InitHisto();
99 }
100
101 AliComparisonDraw::AliComparisonDraw(const AliComparisonDraw& draw):
102   TNamed(draw.GetName(),draw.GetTitle()),
103   fEffTPCPt(draw.fEffTPCPt),      // TPC efficiency as function of Pt (tan+-1)
104   fEffTPCPtMC(draw.fEffTPCPtMC),    // MC -TPC efficiency as function of Pt (tan+-1)
105   fEffTPCPtF(draw.fEffTPCPtF),     // efficiency for findable tracks
106   //
107   fEffTPCTan(draw.fEffTPCTan),   // TPC efficiency as function of Tan (pt>0.15
108   fEffTPCTanMC(draw.fEffTPCTanMC), // MC -TPC efficiency as function of Tan (pt>0.15)
109   fEffTPCTanF(draw.fEffTPCTanF),  // efficiency for findable tracks Tan (pt>0.15)
110   //
111   fEffTPCPtTan(draw.fEffTPCPtTan),    // TPC efficiency as function of Pt and tan
112   fEffTPCPtTanMC(draw.fEffTPCPtTanMC),  // MC -TPC efficiency as function of Pt and tan
113   fEffTPCPtTanF(draw.fEffTPCPtTanF),  // TPC efficiency as function of Pt and tan
114   //
115   // dEdx resolution
116   //
117   fTPCSignalNormTan(draw.fTPCSignalNormTan), // tpc signal normalized to the mean signal - MC
118   fTPCSignalNormSPhi(draw.fTPCSignalNormSPhi),   // tpc signal normalized to the mean signal - MC
119   fTPCSignalNormTPhi(draw.fTPCSignalNormTPhi),   // tpc signal normalized to the mean signal - MC
120   //
121   fTPCSignalNormTanSPhi(draw.fTPCSignalNormTanSPhi),   // tpc signal normalized to the mean signal - MC
122   fTPCSignalNormTanTPhi(draw.fTPCSignalNormTanTPhi),   // tpc signal normalized to the mean signal - MC
123   fTPCSignalNormTanSPt(draw.fTPCSignalNormTanSPt),   // tpc signal normalized to the mean signal - MC
124   //
125   //
126   fPtResolLPT(draw.fPtResolLPT),        // pt resolution - low pt
127   fPtResolHPT(draw.fPtResolHPT),        // pt resolution - high pt 
128   fPtPullLPT(draw.fPtPullLPT),         // pt resolution - low pt
129   fPtPullHPT(draw.fPtPullHPT),         // pt resolution - high pt 
130   //
131   // Resolution constrained param
132   //
133   fCPhiResolTan(draw.fCPhiResolTan),   // angular resolution -  constrained
134   fCTanResolTan(draw.fCTanResolTan),   // angular resolution -  constrained
135   fCPtResolTan(draw.fCPtResolTan),    // pt resolution      -  constrained
136   fCPhiPullTan(draw.fCPhiPullTan),   // angular resolution -  constrained
137   fCTanPullTan(draw.fCTanPullTan),   // angular resolution -  constrained
138   fCPtPullTan(draw.fCPtPullTan),    // pt resolution      -  constrained
139   //
140   // DCA resolution
141   //
142   fD0TanSPtB1(draw.fD0TanSPtB1),   // distance to vertex y  
143   fD1TanSPtB1(draw.fD1TanSPtB1),   // distance to vertex z  
144   fD0TanSPtL1(draw.fD0TanSPtL1),   // distance to vertex y  
145   fD1TanSPtL1(draw.fD1TanSPtL1)   // distance to vertex z  
146 {
147   //
148   // copy constructor
149   //
150 }
151
152 AliComparisonDraw& AliComparisonDraw::operator=(const AliComparisonDraw& info){
153   //
154   // assignment operator
155   //
156   delete this;
157   new (this) AliComparisonDraw(info);
158   return *this;  
159 }
160
161
162
163
164 AliComparisonDraw::~AliComparisonDraw(){
165   //
166   //
167   //
168   delete  fEffTPCPt;      // TPC efficiency as function of Pt (tan+-1)
169   delete  fEffTPCPtMC;    // MC -TPC efficiency as function of Pt (tan+-1)
170   delete  fEffTPCPtF;     // efficiency for findable tracks
171   //
172   delete  fEffTPCTan;   // TPC efficiency as function of Tan (pt>0.15
173   delete  fEffTPCTanMC; // MC -TPC efficiency as function of Tan (pt>0.15)
174   delete  fEffTPCTanF;  // efficiency for findable tracks Tan (pt>0.15)
175   //
176   delete  fEffTPCPtTan;    // TPC efficiency as function of Pt and tan
177   delete  fEffTPCPtTanMC;  // MC -TPC efficiency as function of Pt and tan
178   delete  fEffTPCPtTanF;  // TPC efficiency as function of Pt and tan
179   //
180   // dEdx resolution
181   //
182   delete  fTPCSignalNormTan; // tpc signal normalized to the mean signal - MC
183   delete  fTPCSignalNormSPhi;   // tpc signal normalized to the mean signal - MC
184   delete  fTPCSignalNormTPhi;   // tpc signal normalized to the mean signal - MC
185   //
186   delete  fTPCSignalNormTanSPhi;   // tpc signal normalized to the mean signal - MC
187   delete  fTPCSignalNormTanTPhi;   // tpc signal normalized to the mean signal - MC
188   delete  fTPCSignalNormTanSPt;   // tpc signal normalized to the mean signal - MC
189   //
190   //
191   delete  fPtResolLPT;        // pt resolution - low pt
192   delete  fPtResolHPT;        // pt resolution - high pt 
193   delete  fPtPullLPT;         // pt resolution - low pt
194   delete  fPtPullHPT;         // pt resolution - high pt 
195   //
196   // Resolution constrained param
197   //
198   delete fCPhiResolTan;   // angular resolution -  constrained
199   delete fCTanResolTan;   // angular resolution -  constrained
200   delete fCPtResolTan;    // pt resolution      -  constrained
201   delete fCPhiPullTan;   // angular resolution -  constrained
202   delete fCTanPullTan;   // angular resolution -  constrained
203   delete fCPtPullTan;    // pt resolution      -  constrained
204   //
205   // DCA resolution
206   //
207   delete fD0TanSPtB1;   // distance to vertex y  
208   delete fD1TanSPtB1;   // distance to vertex z  
209   delete fD0TanSPtL1;   // distance to vertex y  
210   delete fD1TanSPtL1;   // distance to vertex z  
211  
212 }
213
214
215
216
217 void AliComparisonDraw::InitHisto(){
218   //
219   //
220   // EFFICIENCY
221   //  
222   // Efficiency as function of pt
223   fEffTPCPt = new TProfile("Eff_pt","Eff_Pt",50,0.1,3);            // physical
224   fEffTPCPtMC = new TProfile("MC_Eff_pt","MC_Eff_Pt",50,0.1,3);    // MC - particles make more than 50 rowdigits
225   fEffTPCPtF = new TProfile("F_Eff_pt","F_Eff_Pt",50,0.1,3);     // tracking - under condition more than 50 rdigits
226
227   // Efficiency as function of pt
228   fEffTPCTan = new TProfile("Eff_tan","Eff_tan",50,-2.5,2.5);            // physical
229   fEffTPCTanMC = new TProfile("MC_Eff_tan","MC_Eff_tan",50,-2.5,2.5);    // MC - particles make more than 50 rowdigits
230   fEffTPCTanF = new TProfile("F_Eff_tan","F_Eff_tan",50,-2.5,2.5);     // tracking - under condition more than 50 rdigits
231
232   fEffTPCPtTan = new TProfile2D("Eff_pt","Eff_Pt",10,0.1,3,20,-2.,2.);
233   fEffTPCPtTanMC = new TProfile2D("MC_Eff_pt","MC Eff Pt",10,0.1,3,20, -2.,2.);
234   fEffTPCPtTanF = new TProfile2D("MC_Eff_pt","MC Eff Pt",10,0.1,3,20, -2.,2.);
235   
236   //
237   // TPC dEdx
238   //
239   fTPCSignalNormTan = new TH2F("CdEdxTan","CdEdxTan",50, -2,2,  40,30,70); // tpc signal normalized to the MC
240   fTPCSignalNormSPhi   = new TH2F("CdEdxSPhi","CdEdxSPhi",10,0.0,1,40,30,70); // tpc signal normalized to the MC
241   fTPCSignalNormTPhi   = new TH2F("CdEdxTPhi","CdEdxTPhi",10,0.0,2,40,30,70); // tpc signal normalized to the MC
242
243   fTPCSignalNormTanSPhi= new TH3F("CdEdxTanSPhi","CdEdxTanSPhi",20, -2,2, 10,0.0 ,1,  40,30,70);   // tpc signal normalized to the mean signal - MC
244   fTPCSignalNormTanTPhi= new TH3F("CdEdxTanTPhi","CdEdxTanTPhi",20, -2,2, 10,0.0 ,1,  40,30,70);   // tpc signal normalized to the mean signal - MC
245   fTPCSignalNormTanSPt= new TH3F("CdEdxTanSPt","CdEdxTanSPt",20, -2,2, 10,0.3 ,3,  40,30,70);   // tpc signal normalized to the mean signal - MC
246
247
248
249   //
250   // RESOLUTION
251   //
252   fCPtResolTan = new TH2F("C Pt resol","C pt resol",50, -2,2,200,-0.2,0.2);
253   fCPtPullTan = new TH2F("C Pt pull","C pt pull",50, -2,2,200,-5,5);
254   //
255   fCPhiResolTan = new TH2F("CPhiResolTan","CPhiResolTan",50, -2,2,200,-0.025,0.025);   
256   // angular resolution -  constrained
257   fCTanResolTan = new TH2F("CTanResolTan","CTanResolTan",50, -2,2,200,-0.025,0.025);
258   // angular resolution -  constrained
259   fCPtResolTan=new TH2F("CPtResol","CPtResol",50, -2,2,200,-0.2,0.2);;    
260   // pt resolution      -  constrained
261   fCPhiPullTan = new TH2F("CPhiPullTan","CPhiPullTan",50, -2,2,200,-5,5);   
262   // angular resolution -  constrained
263   fCTanPullTan = new TH2F("CTanPullTan","CTanPullTan",50, -2,2,200,-5,5);
264   // angular resolution -  constrained
265   fCPtPullTan=new TH2F("CPtPull","CPtPull",50, -2,2,200,-5,5);    
266   // pt resolution      -  constrained
267   //
268   fPtResolLPT = new TH2F("Pt resol","pt resol",10, 0.1,3,200,-0.2,0.2);
269   fPtResolHPT = new TH2F("Pt resol","pt resol",10, 2,100,200,-0.3,0.3);  
270   fPtPullLPT = new TH2F("Pt pool","pt pool",10, 0.1,3,200,-6,6);
271   fPtPullHPT = new TH2F("Pt pool","pt pool",10, 2,100,200,-6,6);  
272   //
273   fD0TanSPtB1 = new TH3F("DCAyTanSPt","DCAyTanSPt",20,1,2, 10,0.3,2, 100,-4,4);
274   fD1TanSPtB1 = new TH3F("DCAzTanSPt","DCAzTanSPt",20,1,2, 10,0.3,2, 100,-4,4);
275   fD0TanSPtL1 = new TH3F("DCAyTanSPt","DCAyTanSPt",20,0,1, 10,0.3,2, 100,-0.1,0.1);
276   fD1TanSPtL1 = new TH3F("DCAzTanSPt","DCAzTanSPt",20,0,1, 10,0.3,2, 100, -0.1,0.1);
277
278
279
280 }
281
282 void AliComparisonDraw::ProcessEff(AliMCInfo* infoMC, AliESDRecInfo *infoRC){
283   //
284   // make efficiencies histograms
285   //
286   Float_t kptcut = 0.15;
287   Float_t ktancut=1.;
288   Int_t   kmincl =50;
289   Float_t mcpt = infoMC->GetParticle().Pt();
290   Float_t tantheta = TMath::Tan(infoMC->GetParticle().Theta()-TMath::Pi()*0.5);
291   Bool_t isPrim = infoMC->GetParticle().R()<0.1 && TMath::Abs(infoMC->GetParticle().Vz())<10;
292   //z diamond and 
293   
294   if (!isPrim) return;
295
296   //pt
297   if (TMath::Abs(tantheta)<ktancut){
298     fEffTPCPt->Fill(mcpt, infoRC->GetStatus(1)==3);
299     fEffTPCPtMC->Fill(mcpt, infoMC->GetRowsWithDigits()>kmincl);
300     if (infoMC->GetRowsWithDigits()>kmincl){
301       fEffTPCPtF->Fill(mcpt, infoRC->GetStatus(1)==3);
302     }
303   }
304
305   //theta
306   if (TMath::Abs(mcpt)>kptcut){
307     fEffTPCTan->Fill(tantheta, infoRC->GetStatus(1)==3);
308     fEffTPCTanMC->Fill(tantheta, infoMC->GetRowsWithDigits()>kmincl);
309     if (infoMC->GetRowsWithDigits()>kmincl){
310       fEffTPCTanF->Fill(tantheta, infoRC->GetStatus(1)==3);
311     }
312   }
313   // 
314   // pt-theta
315   //
316   fEffTPCPtTan->Fill(mcpt,tantheta,infoRC->GetStatus(1)==3);
317   fEffTPCPtTanMC->Fill(mcpt,tantheta,infoMC->GetRowsWithDigits()>50); 
318   if (infoMC->GetRowsWithDigits()>kmincl){
319     fEffTPCPtTanF->Fill(mcpt,tantheta,infoRC->GetStatus(1)==3); 
320   }
321 }
322
323
324 void AliComparisonDraw::ProcessResolConstrained(AliMCInfo* infoMC, AliESDRecInfo *infoRC){
325   //
326   //
327   //
328   Float_t mcpt = infoMC->GetParticle().Pt();
329   Float_t tantheta = TMath::Tan(infoMC->GetParticle().Theta()-TMath::Pi()*0.5);
330   Bool_t isPrim = infoMC->GetParticle().R()<0.1 && TMath::Abs(infoMC->GetParticle().Vz())<10;
331   //z diamond and 
332   
333   if (!isPrim) return;
334   if (infoRC->GetStatus(1)!=3) return;
335   if (!infoRC->GetESDtrack()) return;  
336   if (infoRC->GetESDtrack()->GetTPCNcls()<10) return;
337   if (!infoRC->GetESDtrack()->GetConstrainedParam()) return;
338   
339   //
340   // constrained parameters resolution
341   //
342   const AliExternalTrackParam * cparam = infoRC->GetESDtrack()->GetConstrainedParam();
343   Float_t deltaCPt= (mcpt-cparam->Pt())/mcpt;  
344   Float_t pullCPt= (1/mcpt-cparam->OneOverPt())/
345     TMath::Sqrt(cparam->GetSigma1Pt2());          
346   Float_t deltaPhi = TMath::ATan2(cparam->Py(),cparam->Px())-
347     TMath::ATan2(infoMC->GetParticle().Py(),infoMC->GetParticle().Px());
348   Float_t pullPhi = deltaPhi/TMath::Sqrt(cparam->GetSigmaSnp2()); 
349
350   Float_t deltaTan = TMath::ATan2(cparam->Pz(),cparam->Pt())-
351     TMath::ATan2(infoMC->GetParticle().Pz(),infoMC->GetParticle().Pt());
352   Float_t pullTan = deltaPhi/TMath::Sqrt(cparam->GetSigmaSnp2()); 
353
354   fCPtResolTan->Fill(tantheta,deltaCPt);
355   fCPtPullTan->Fill(tantheta,pullCPt);
356   fCPhiResolTan->Fill(tantheta,deltaPhi);
357   fCPhiPullTan->Fill(tantheta,pullPhi);
358   fCTanResolTan->Fill(tantheta,deltaTan);
359   fCTanPullTan->Fill(tantheta,pullTan);
360
361 }
362
363
364
365 void  AliComparisonDraw::ProcessTPCdedx(AliMCInfo* infoMC, AliESDRecInfo *infoRC){
366   //
367   //
368   //
369   Float_t mcpt = infoMC->GetParticle().Pt();
370   Float_t tantheta = TMath::Tan(infoMC->GetParticle().Theta()-TMath::Pi()*0.5);
371   Bool_t isPrim = infoMC->GetParticle().R()<0.1 && TMath::Abs(infoMC->GetParticle().Vz())<10;
372   //z diamond and 
373   
374   if (!isPrim) return;
375   if (infoRC->GetStatus(1)!=3) return;
376   if (!infoRC->GetESDtrack()) return;  
377   if (infoRC->GetESDtrack()->GetTPCNcls()<10) return;
378   if (!infoRC->GetESDtrack()->GetConstrainedParam()) return;
379   Float_t mprim = infoMC->GetPrim();
380   if (mprim>1.4) return;
381   if (mprim<0.5) return;
382   if (infoRC->GetESDtrack()->GetTPCsignalN()<50) return;
383   //
384   Float_t ratio = infoRC->GetESDtrack()->GetTPCsignal()/infoMC->GetPrim();
385   Float_t sphi =  infoRC->GetESDtrack()->GetInnerParam()->GetSnp();
386   Float_t tphi =  sphi/TMath::Sqrt((1.-sphi)*(1.+sphi));
387
388
389   if (TMath::Abs(infoMC->GetParticle().GetPdgCode())!=211) return;
390   if (mcpt>0.5){
391     fTPCSignalNormTan->Fill(tantheta,ratio);    //only subset
392   }
393   if (TMath::Abs(tantheta)<0.5){
394     fTPCSignalNormSPhi->Fill(sphi,ratio);        // only subset
395     fTPCSignalNormTPhi->Fill(tphi,ratio);        // only subset
396   }
397   fTPCSignalNormTanSPhi->Fill(tantheta,sphi,ratio);    
398   fTPCSignalNormTanTPhi->Fill(tantheta,tphi,ratio);    
399   fTPCSignalNormTanSPt->Fill(tantheta,TMath::Sqrt(mcpt),ratio);    
400 }
401
402 void      AliComparisonDraw::ProcessDCA(AliMCInfo* infoMC, AliESDRecInfo *infoRC){
403   //
404   //
405   //
406   Float_t mcpt = infoMC->GetParticle().Pt();
407   Float_t tantheta = TMath::Tan(infoMC->GetParticle().Theta()-TMath::Pi()*0.5);
408   Bool_t isPrim = infoMC->GetParticle().R()<0.1 && TMath::Abs(infoMC->GetParticle().Vz())<10;
409   //z diamond and 
410   if (!isPrim) return;
411   if (infoRC->GetStatus(1)!=3) return;
412   if (!infoRC->GetESDtrack()) return;  
413   if (infoRC->GetESDtrack()->GetTPCNcls()<10) return;
414   if (!infoRC->GetESDtrack()->GetConstrainedParam()) return;
415   Float_t spt = TMath::Sqrt(mcpt);
416   Float_t dca[2],cov[3];
417   infoRC->GetESDtrack()->GetImpactParameters(dca,cov);
418   Int_t clusterITS[100];
419   if (infoRC->GetESDtrack()->GetITSclusters(clusterITS)==0){
420     fD0TanSPtB1->Fill(tantheta,spt,dca[0]);
421     fD1TanSPtB1->Fill(tantheta,spt,dca[1]);
422   }
423   fD0TanSPtL1->Fill(tantheta,spt,dca[0]);
424   fD1TanSPtL1->Fill(tantheta,spt,dca[1]);  
425 }
426
427
428
429
430 void AliComparisonDraw::Process(AliMCInfo* infoMC, AliESDRecInfo *infoRC){
431   //
432   // 
433   //
434   ProcessEff(infoMC,infoRC);
435   ProcessResolConstrained(infoMC,infoRC);
436   ProcessTPCdedx(infoMC, infoRC);
437   ProcessDCA(infoMC, infoRC);
438
439   Float_t mcpt = infoMC->GetParticle().Pt();
440   Bool_t isPrim = infoMC->GetParticle().R()<0.1 && TMath::Abs(infoMC->GetParticle().Vz())<10;
441   //z diamond and 
442   
443   if (!isPrim) return;
444   //
445   //
446   if (infoRC->GetStatus(1)==0) return;
447   if (!infoRC->GetESDtrack()) return;  
448   if (infoRC->GetESDtrack()->GetTPCNcls()<10) return;
449   //  printf("Pt\t%f\t%f\n",mcpt, infoRC->GetESDtrack()->Pt());
450   
451   Float_t deltaPt= (mcpt-infoRC->GetESDtrack()->Pt())/mcpt;  
452   Float_t poolPt= (1/mcpt-infoRC->GetESDtrack()->OneOverPt())/
453     TMath::Sqrt(infoRC->GetESDtrack()->GetSigma1Pt2());  
454
455   fPtResolLPT->Fill(mcpt,deltaPt);
456   fPtResolHPT->Fill(mcpt,deltaPt);
457   fPtPullLPT->Fill(mcpt,poolPt);
458   fPtPullHPT->Fill(mcpt,poolPt);  
459 }
460
461
462
463 TH1F* AliComparisonDraw::MakeResol(TH2F * his, Int_t integ, Bool_t type){
464   TH1F *hisr, *hism;
465   if (!gPad) new TCanvas;
466   hisr = AliTreeDraw::CreateResHistoI(his,&hism,integ);
467   if (type) return hism;
468   else 
469     return hisr;
470 }
471
472
473 TGraph2D * AliComparisonDraw::MakeStat2D(TH3 * his, Int_t delta0, Int_t delta1, Int_t type){
474   //
475   //
476   //
477   // delta - number of bins to integrate
478   // type - 0 - mean value
479
480   TAxis * xaxis  = his->GetXaxis();
481   TAxis * yaxis  = his->GetYaxis();
482   //  TAxis * zaxis  = his->GetZaxis();
483   Int_t   nbinx  = xaxis->GetNbins();
484   Int_t   nbiny  = yaxis->GetNbins();
485   char name[1000];
486   Int_t icount=0;
487   TGraph2D  *graph = new TGraph2D(nbinx*nbiny);
488   TF1 f1("f1","gaus");
489   for (Int_t ix=0; ix<nbinx;ix++)
490     for (Int_t iy=0; iy<nbiny;iy++){
491       Float_t xcenter = xaxis->GetBinCenter(ix); 
492       Float_t ycenter = yaxis->GetBinCenter(iy); 
493       sprintf(name,"%s_%d_%d",his->GetName(), ix,iy);
494       TH1 *projection = his->ProjectionZ(name,ix-delta0,ix+delta0,iy-delta1,iy+delta1);
495       Float_t stat= 0;
496       if (type==0) stat = projection->GetMean();
497       if (type==1) stat = projection->GetRMS();
498       if (type==2 || type==3){
499         TVectorD vec(3);
500         AliMathBase::LTM((TH1F*)projection,&vec,0.7);
501         if (type==2) stat= vec[1];
502         if (type==3) stat= vec[0];      
503       }
504       if (type==4|| type==5){
505         projection->Fit(&f1);
506         if (type==4) stat= f1.GetParameter(1);
507         if (type==5) stat= f1.GetParameter(2);
508       }
509       //printf("%d\t%f\t%f\t%f\n", icount,xcenter, ycenter, stat);
510       graph->SetPoint(icount,xcenter, ycenter, stat);
511       icount++;
512     }
513   return graph;
514 }
515
516 TGraph * AliComparisonDraw::MakeStat1D(TH3 * his, Int_t delta1, Int_t type){
517   //
518   //
519   //
520   // delta - number of bins to integrate
521   // type - 0 - mean value
522
523   TAxis * xaxis  = his->GetXaxis();
524   TAxis * yaxis  = his->GetYaxis();
525   //  TAxis * zaxis  = his->GetZaxis();
526   Int_t   nbinx  = xaxis->GetNbins();
527   Int_t   nbiny  = yaxis->GetNbins();
528   char name[1000];
529   Int_t icount=0;
530   TGraph  *graph = new TGraph(nbinx);
531   TF1 f1("f1","gaus");
532   for (Int_t ix=0; ix<nbinx;ix++){
533     Float_t xcenter = xaxis->GetBinCenter(ix); 
534     //    Float_t ycenter = yaxis->GetBinCenter(iy); 
535     sprintf(name,"%s_%d",his->GetName(), ix);
536     TH1 *projection = his->ProjectionZ(name,ix-delta1,ix+delta1,0,nbiny);
537     Float_t stat= 0;
538     if (type==0) stat = projection->GetMean();
539     if (type==1) stat = projection->GetRMS();
540     if (type==2 || type==3){
541       TVectorD vec(3);
542         AliMathBase::LTM((TH1F*)projection,&vec,0.7);
543         if (type==2) stat= vec[1];
544         if (type==3) stat= vec[0];      
545     }
546     if (type==4|| type==5){
547       projection->Fit(&f1);
548       if (type==4) stat= f1.GetParameter(1);
549       if (type==5) stat= f1.GetParameter(2);
550     }
551       //printf("%d\t%f\t%f\t%f\n", icount,xcenter, ycenter, stat);
552     graph->SetPoint(icount,xcenter, stat);
553     icount++;
554   }
555   return graph;
556 }
557
558 //
559 // Make derived plots
560 //
561
562 void AliComparisonDraw::MakePlots(){
563   //
564   //
565   //
566   AliComparisonDraw * comp=this;
567
568   TFile *fp = new TFile("picutures.root","recreate");
569   TH1F *hiss=0;
570   //TH1F *hism=0;
571   TGraph2D * gr=0, gr2=0;
572   TGraph * gr0 = 0;
573   TCanvas * c = new TCanvas("Phi resol Tan","Phi resol Tan");
574   //
575   //
576   //
577   hiss = comp->MakeResol(comp->fCPtResolTan,1,0);
578   hiss->SetXTitle("Tan(#theta)");
579   hiss->SetYTitle("#sigmap_{t}/p_{t}");
580   hiss->Draw(); 
581   hiss->Write("CptResolTan");
582   //
583   //
584   hiss = comp->MakeResol(comp->fCPhiResolTan,1,0);
585   c->cd();
586   hiss->SetXTitle("Tan(#theta)");
587   hiss->SetYTitle("#sigma#phi (rad)");
588   hiss->Draw();
589   fp->cd();
590   hiss->Write("PhiResolTan");
591   //
592   hiss = comp->MakeResol(comp->fCTanResolTan,1,0);
593   c->cd();
594   hiss->SetXTitle("Tan(#theta)");
595   hiss->SetYTitle("#sigma#theta (rad)");
596   hiss->Draw();
597   fp->cd();
598   hiss->Write("ThetaResolTan");
599   //
600   //
601   hiss = comp->MakeResol(comp->fCTanResolTan,1,0);
602   c->cd();
603   hiss->SetXTitle("Tan(#theta)");
604   hiss->SetYTitle("#sigmap_{t}/p_{t} ");
605   hiss->Draw();
606   fp->cd();
607   //
608   //
609   //
610   hiss = comp->MakeResol(comp->fTPCSignalNormTan,4,0);
611   hiss->SetXTitle("Tan(#theta)");
612   hiss->SetYTitle("#sigma_{dEdx}");
613   hiss->Draw();
614   fp->cd();
615   hiss->Write("TPCdEdxResolTan");
616   //
617   //
618   //
619   hiss = comp->MakeResol(comp->fTPCSignalNormTan,4,1); 
620   hiss->SetXTitle("Tan(#theta)");
621   hiss->SetYTitle("<dEdx>");
622   hiss->Draw(); 
623   hiss->Write("TPCdEdxMeanTan");
624   //
625   //
626   gr = comp->MakeStat2D(comp->fTPCSignalNormTanSPt,3,1,4);
627   gr->GetXaxis()->SetTitle("Tan(#theta)");
628   gr->GetYaxis()->SetTitle("#sqrt{p_{t}(GeV)}");
629   gr->GetZaxis()->SetTitle("<dEdx>");
630   gr->Draw("colz"); 
631   gr->GetHistogram()->Write("TPCdEdxMeanTanPt");
632   //
633   //
634   gr = comp->MakeStat2D(comp->fTPCSignalNormTanSPt,3,1,5);
635   gr->GetXaxis()->SetTitle("Tan(#theta)");
636   gr->GetYaxis()->SetTitle("#sqrt{p_{t}(GeV)}");
637   gr->GetZaxis()->SetTitle("#sigma_{dEdx}");
638   gr->Draw("colz"); 
639   gr->GetHistogram()->Write("TPCdEdxMeanTanPt");
640   //
641   //
642   //
643   comp->fEffTPCTanF->SetXTitle("Tan(#theta)");
644   comp->fEffTPCTanF->SetYTitle("eff_{findable}");
645   comp->fEffTPCTanF->Draw();
646   comp->fEffTPCTanF->Write("EffTanFindable");
647   //
648   //
649   comp->fEffTPCTan->SetXTitle("Tan(#theta)");
650   comp->fEffTPCTan->SetYTitle("eff_{all}");
651   comp->fEffTPCTan->Draw();
652   comp->fEffTPCTan->Write("EffTanAll");
653   //
654   //DCA resolution
655   //
656   gr0 = comp->MakeStat1D(comp->fD0TanSPtB1,2,5);
657   gr0->GetXaxis()->SetTitle("Tan(#theta)");
658   gr0->GetYaxis()->SetTitle("#sigmaDCA (cm)");
659   gPad->Clear();
660   gr0->Draw("al*");
661   gr->GetHistogram()->Write("DCAResolTan");
662   //
663   //
664   //
665   gr = comp->MakeStat2D(comp->fD0TanSPtB1,4,2,5); 
666   gr0->GetXaxis()->SetTitle("Tan(#theta)");
667   gr0->GetYaxis()->SetTitle("#sigmaDCA (cm)");
668   gPad->Clear();
669   gr0->Draw("al*");
670   gr->GetHistogram()->Write("DCAResolSPTTan");
671
672   fp->Close();
673
674
675 }
676
677
678