Updates in the ITSsa track QA plotting (Leonardo)
[u/mrichter/AliRoot.git] / ITS / PlotOutputQATaskITSsa.C
1 #if !defined(__CINT__) || defined(__MAKECINT__)
2 #include <TCanvas.h>
3 #include <TGrid.h>
4 #include <TFile.h>
5 #include <TList.h>
6 #include <TPaveStats.h>
7 #include <TGraph.h>
8 #include <TGraphErrors.h>
9 #include <TH1.h>
10 #include <TF1.h>
11 #include <TH2.h>
12 #include <TLegend.h>
13 #include <TLegendEntry.h>
14 #include <TLatex.h>
15 #include <TStyle.h>
16 #include <TROOT.h>
17 #endif
18
19 /*  $Id$    */
20
21 //-------------------------------------------------------
22 //
23 // Macro do plot the output histograms of the QA task for ITS standalone tracks
24 // dE/dx for ITS pure standalone tracks with parameterization
25 // dE/dx resolution for ITS pure standalone tracks with parameterization
26 // d0 resolution and bias with gaussian fit
27 // d0 resolution and bias with gaussian+tail fit
28 // Comparison of the two methods for d0 resolution
29 // Pt resolution (matching ITSèureSA with TPC+ITS tracks)
30 // General Plots: ratios between ITSsa, ITSpureSA and TPC+ITS tracks 
31 //                eta phi distributions of tracks
32 //                number of clusters per track
33 //
34 // Authors: Leonardo Milano, Francesco Prino
35 //
36 //-------------------------------------------------------
37
38 Int_t color[3]={1,2,4};
39 TString particle[3]={"Pion","Kaon","Proton"};
40
41
42 void PlotOutputQATaskITSsa(TString filename="LHC11h_QA95_MB_4cls.root",TString foutName="OutputQAITSstandaloneTracks_LHC11h_QA95_MB_4cls.root",Bool_t fMC=0){
43   
44   gROOT->SetStyle("Plain");
45   gStyle->SetPalette(1); 
46
47   if(filename.Contains("alien")) TGrid::Connect("alien:");
48   TFile* f=TFile::Open(filename.Data());
49   TDirectory *dirFile=(TDirectory*)f->Get("ITSsaTracks");
50   TString lname="clistITSsaTracks";
51   TList *li = (TList*)dirFile->Get(lname.Data());
52   
53   TFile *fout=new TFile(foutName.Data(),"recreate");
54   fout->Close();
55   delete fout;
56   
57   Bool_t useHybridITSsaParameterization=1; //BetheITSsaHybrid committed in Rev. 59830
58   Bool_t PlotPIDBands=0; //Bands used for PID in the spectra analysis
59   dedxPlot(li,fMC,foutName,useHybridITSsaParameterization,PlotPIDBands);
60   dedxRes(li,foutName);
61   ImpParRes(li,foutName);    
62   ImpParResGausTail(li,foutName);    
63   CompareImpParRes(foutName);
64   Bool_t optFromMC=0;
65   PlotPtResol(li,foutName,optFromMC);
66   PlotITSsa(li,foutName);
67 }    
68
69
70
71 //----------------------------------------------------------------------------
72
73 void dedxPlot(TList *li,Bool_t fMC,TString foutName,Bool_t useHybridITSsaParameterization,Bool_t PlotPIDBands)    
74 {
75   // Plot the dedx vs momentum from the ITSsaTracks QA output
76   // Bethe Bloch parameterization is taken from AliITSPIDResponse
77   // If useHybridITSsaParameterization the Hybrid ITS parameterization (Phobos + polinomial at low beta*gamma) is used. This requires AliRoot >  Rev. 59830
78   
79   TH2F *fHistDEDX=(TH2F*)li->FindObject("hdedxvsP2clsITSpureSA");
80   fHistDEDX->Add((TH2F*)li->FindObject("hdedxvsP3clsITSpureSA"));
81   fHistDEDX->Add((TH2F*)li->FindObject("hdedxvsP4clsITSpureSA"));
82   
83   TCanvas *c=new TCanvas("dEdx_vs_p","dEdx_vs_p",1);
84   c->SetLogz();
85   c->SetLogx();
86   fHistDEDX->SetXTitle("#it{p} (GeV/#it{c})");
87   fHistDEDX->SetYTitle("d#it{E}/dx (KeV/300#mum)");
88   fHistDEDX->GetXaxis()->SetRangeUser(0.08,3);
89   fHistDEDX->GetYaxis()->SetRangeUser(0.,700);
90   fHistDEDX->GetYaxis()->SetTitleSize(0.05);
91   fHistDEDX->GetXaxis()->SetTitleSize(0.05);
92   //fHistDEDX->SetMinimum(2000);
93   //fHistDEDX->SetMaximum(18000);
94   fHistDEDX->DrawClone("col");
95   
96   Float_t pdgmass[5]={0.13957,0.493677,0.938272,1.875612762,0.00996}; //mass for pi, K, P, d (Gev/c^2)
97   Double_t dedxvalue[5];
98   Double_t funcvalue[7];
99   TGraph *fdedx[7];
100   AliITSPIDResponse *pidresp=new AliITSPIDResponse(fMC);
101   
102   for(Int_t i7=0;i7<7;i7++){
103     fdedx[i7]=new TGraph();
104     fdedx[i7]->SetName(Form("%d",i7));
105     fdedx[i7]->SetTitle(Form("%d",i7));
106     fdedx[i7]->SetLineWidth(2);
107     fdedx[i7]->SetLineColor(1);
108   }
109   
110   Float_t mom=0.07.;//low pt range
111   Float_t step=(3-mom)/100;
112   Double_t lr[7]={0.09,0.09,0.09,0.1,0.2,.25,.3};
113   
114   for(Int_t i=0;i<100;i++){
115     if(useHybridITSsaParameterization)for(Int_t i2=0;i2<5;i2++)dedxvalue[i2]=pidresp->BetheITSsaHybrid(mom,pdgmass[i2]);
116     else for(Int_t i2=0;i2<5;i2++)dedxvalue[i2]=pidresp->Bethe(mom,pdgmass[i2],1);
117     funcvalue[0]=dedxvalue[0]-0.2*dedxvalue[0];
118     funcvalue[1]=dedxvalue[0];
119     funcvalue[2]=0.5*(dedxvalue[0]+dedxvalue[1]);
120     funcvalue[3]=dedxvalue[1];
121     funcvalue[4]=0.5*(dedxvalue[2]+dedxvalue[1]);
122     funcvalue[5]=dedxvalue[2];
123     funcvalue[6]=0.5*(dedxvalue[2]+dedxvalue[3]);
124     
125     for(Int_t i2=0;i2<7;i2++){
126       fdedx[i2]->SetPoint(i,mom,funcvalue[i2]);
127     }
128     mom+=step;
129   }
130   fdedx[0]->SetLineStyle(2);
131   fdedx[2]->SetLineStyle(2);
132   fdedx[4]->SetLineStyle(2);
133   fdedx[6]->SetLineStyle(2);
134   for(Int_t i2=0;i2<7;i2++){
135     fdedx[i2]->RemovePoint(0);
136     fdedx[i2]->GetXaxis()->SetRangeUser(lr[i2],3);
137     if(!PlotPIDBands)if(i2==0 || i2==2 || i2==4 || i2==6)continue;
138     fdedx[i2]->DrawClone("Lsame");
139   }
140   TLatex *tex=new TLatex(0.55,630,"ITS standalone tracking");
141   tex->SetTextSize(0.04);
142   tex->DrawClone();
143   
144   TLatex *pi=new TLatex(0.1223724,176.0205,"#pi");
145   TLatex *k=new TLatex(0.263512,291.1394,"K");
146   TLatex *p=new TLatex(0.517675,296.7096,"p");
147   TLatex *e=new TLatex(0.09856606,101.7502,"e");
148   pi->DrawClone();
149   k->DrawClone();
150   p->DrawClone();
151   e->DrawClone();
152   
153   //add electron line
154   Float_t mom=0.09.;//low pt range
155   Float_t step=(0.16-mom)/10;
156   TGraph* fElectr=new TGraph();
157   fElectr->SetName("fElectr");
158   fElectr->SetTitle("fElectr");
159   AliITSPIDResponse *pidresp=new AliITSPIDResponse(0);
160   for(Int_t i=0;i<10;i++){
161     if(useHybridITSsaParameterization)fElectr->SetPoint(i,mom,pidresp->BetheITSsaHybrid(mom,0.00051));
162     else  fElectr->SetPoint(i,mom,pidresp->Bethe(mom,0.00051,1));
163     mom+=step;
164   }
165   fElectr->SetLineWidth(2);
166   fElectr->SetLineColor(1);
167   fElectr->DrawClone("Lsame");
168   
169   TFile *fout=new TFile(foutName.Data(),"update");
170   c->Write();
171   fout->Close();
172   delete fout;
173   }
174
175 //----------------------------------------------------------------------------
176
177 void dedxRes(TList *li,TString foutName)    
178 {
179   // plot the dedx resolution as a function of momuntum from the ITSsaTracks QA
180   // Resolution is calculated from Gaussian fit to the pion peak for tracks with 2, 3 and 4 clusters in the ITS
181   
182   TCanvas *cdedxRes=new TCanvas("dedxRes_vs_p","dedxRes_vs_p",1000,800);
183   cdedxRes->SetGridx();
184   cdedxRes->SetLogx();
185   
186   for(Int_t ncls=2;ncls<5;ncls++){
187     TH2F *hdedxvsPITSpureSA=(TH2F*)li->FindObject(Form("hdedxvsP%iclsITSpureSA",ncls));
188     //binning
189     const Int_t nbins = 29;
190     Double_t xbins[nbins+1]={0.06,0.08,0.10,0.12,0.14,0.16,0.18,0.20,0.25,0.30,
191                              0.35,0.40,0.45,0.50,0.55,0.60,0.65,0.70,0.75,0.80,
192                              0.85,0.90,0.95,1.00,1.20,1.40,1.60,1.80,1.90,2.00};
193     TH1F *fHistdedxRes = new TH1F(Form("fHistdedxRes%icls",ncls),Form("fHistdedxRes%icls",ncls),nbins,xbins);
194     fHistdedxRes->GetXaxis()->SetTitle("#font[52]{p} (GeV/#font[52]{c})");
195     fHistdedxRes->GetYaxis()->SetTitle("de/dx relative resolution");
196     fHistdedxRes->GetYaxis()->SetTitleOffset(1.5);
197     TH1F *fHistDCA[nbins];
198     TF1 *fPar = new TF1("fPar","gaus",0,1000);
199     for(Int_t m=0;m<nbins;m++){
200       fHistDCA[m]= (TH1F*)hdedxvsPITSpureSA->ProjectionY(Form("%i",m),hdedxvsPITSpureSA->GetXaxis()->FindBin(xbins[m]+0.0000001),hdedxvsPITSpureSA->GetXaxis()->FindBin(xbins[m+1]-0.0000001));
201     }
202     TCanvas *cgaus=new TCanvas(Form("cgausFitdEdxRes_%icls",ncls),Form("cgausFitdEdxRes_%icls",ncls),1000,800);
203     cgaus->Divide(8,4,0.001,0.001);
204     for(Int_t i=0; i<nbins; i++){
205       cgaus->cd(i+1);
206       fHistDCA[i]->SetFillColor(16);
207       Float_t MaxPosition=fHistDCA[i]->GetBinCenter(fHistDCA[i]->GetMaximumBin());
208       Float_t minFit=MaxPosition-0.1*MaxPosition;
209       Float_t maxFit=MaxPosition+0.1*MaxPosition;
210       fPar->SetParameter(1,MaxPosition);
211       fHistDCA[i]->Fit(fPar,"NM","",minFit,maxFit);
212       fHistDCA[i]->GetXaxis()->SetRangeUser(minFit-50,maxFit+50);
213       fPar->SetLineColor(1);
214       fHistDCA[i]->DrawClone();
215       fPar->DrawClone("same");
216       fHistdedxRes->Fill((xbins[i]+xbins[i+1])/2,fPar->GetParameter(2)/fPar->GetParameter(1));
217       fHistdedxRes->SetBinError(fHistdedxRes->FindBin((xbins[i]+xbins[i+1])/2),fPar->GetParError(2)/fPar->GetParameter(1));
218     }
219     cdedxRes->cd();
220     setDrawAtt(20,ncls,1,ncls,1,fHistdedxRes);
221     fHistdedxRes->SetMinimum(0.0);
222     fHistdedxRes->SetMaximum(0.2);
223     fHistdedxRes->GetXaxis()->SetRangeUser(0.2,2);  
224     if(ncls==2)fHistdedxRes->DrawCopy("p");
225     else fHistdedxRes->DrawCopy("psame");
226     
227     TF1 *fFit = new TF1(Form("fFitdEdxRes_%d",ncls),"pol0",0.6,2);
228     fFit->SetLineColor(ncls);
229     fHistdedxRes->Fit(fFit,"NM","",0.6,2);
230     fFit->DrawCopy("same");
231     TPaveText *tpave1=new TPaveText(0.3,0.2-0.15*ncls,0.7,0.89,"brNDC");
232     tpave1->SetBorderSize(0);
233     tpave1->SetFillStyle(0);
234     tpave1->SetFillColor(0);
235     tpave1->SetTextColor(ncls);
236     TText *txt2=tpave1->AddText(Form("pol 0: %.3f",fFit->GetParameter(0)));
237     txt2->SetTextFont(62);
238     tpave1->DrawClone();
239     
240     TFile *fout=new TFile(foutName.Data(),"update");
241     fHistdedxRes->Write();
242     fFit->Write();
243     fout->Close();
244     delete fout;
245   }
246   
247   cdedxRes->BuildLegend()->SetFillColor(0);
248   TPaveText *tpave=new TPaveText(0.5,0.9,0.8,0.99,"brNDC");
249   tpave->SetBorderSize(0);
250   tpave->SetFillStyle(0);
251   tpave->SetFillColor(0);
252   tpave->SetTextColor(2);
253   TText *txt1=tpave->AddText("ITS standalone");
254   txt1->SetTextFont(62);
255   tpave->DrawClone();
256   
257   TFile *fout=new TFile(foutName.Data(),"update");
258   cdedxRes->Write();
259   fout->Close();
260   delete fout;
261 }
262
263 //----------------------------------------------------------------------------
264
265 void ImpParRes(TList *li,TString foutName)    
266 {
267   //Plot the resolution on the tranverse component of impact parmeter from the ITSsaTracks QA output as a function of pt
268   // Gaussian fit is used
269
270   TH2F *hd0rphiITSpureSA[3];
271   //binning
272   const Int_t nbins = 29;
273   Double_t xbins[nbins+1]={0.06,0.08,0.10,0.12,0.14,0.16,0.18,0.20,0.25,0.30,
274                            0.35,0.40,0.45,0.50,0.55,0.60,0.65,0.70,0.75,0.80,
275                            0.85,0.90,0.95,1.00,1.20,1.40,1.60,1.80,1.90,2.00};
276   Int_t MinBin[3]={2,7,9};
277   
278   TH1F *fHistDCA[nbins];  
279   TString particle[3]={"Pion","Kaon","Proton"};
280   
281   TCanvas *cImpPar=new TCanvas("ImpParRes_vs_pt_Gaus","ImpParRes_vs_pt_Gaus",1000,800);
282   cImpPar->SetGridx();
283   cImpPar->SetLogx();
284   TCanvas *cImpParMean=new TCanvas("ImpParMean_vs_pt_Gaus","ImpParMean_vs_pt_Gaus",1000,800);
285   cImpParMean->SetGridx();
286   cImpParMean->SetLogx();
287   
288   for(Int_t iparticle=0;iparticle<3;iparticle++){
289     hd0rphiITSpureSA[iparticle]=(TH2F*)li->FindObject(Form("hd0rphiITSpureSA%s",particle[iparticle].Data()));
290     
291     TH1F *fHistImpParRes = new TH1F(Form("fHistImpParResGaus%s",particle[iparticle].Data()),Form("%s",particle[iparticle].Data()),nbins,xbins);
292     fHistImpParRes->GetXaxis()->SetTitle("#font[52]{p}_{T} (GeV/#font[52]{c})");
293     fHistImpParRes->GetYaxis()->SetTitle("d0 r#phi resolution [#mum]");
294     fHistImpParRes->GetYaxis()->SetTitleOffset(1.5);
295     TH1F *fHistImpParMean = new TH1F(Form("fHistImpParMeanGaus%s",particle[iparticle].Data()),Form("%s",particle[iparticle].Data()),nbins,xbins);
296     fHistImpParMean->GetXaxis()->SetTitle("#font[52]{p}_{T} (GeV/#font[52]{c})");
297     fHistImpParMean->GetYaxis()->SetTitle("d0 r#phi mean [#mum]");
298     fHistImpParMean->GetYaxis()->SetTitleOffset(1.5);
299     
300     TF1 *fPar = new TF1("fPar","gaus",-1,1);
301     
302     for(Int_t m=0;m<nbins;m++){
303       fHistDCA[m]= (TH1F*)hd0rphiITSpureSA[iparticle]->ProjectionY(Form("%s%i",particle[iparticle].Data(),m),hd0rphiITSpureSA[iparticle]->GetXaxis()->FindBin(xbins[m]+0.000001),hd0rphiITSpureSA[iparticle]->GetXaxis()->FindBin(xbins[m+1]-0.000001));
304     }
305     
306     TCanvas *cgaus=new TCanvas(Form("cgausFit_ImpParRes_vs_pt_%s",particle[iparticle].Data()),Form("cgausFit_ImpParRes_vs_pt_%s",particle[iparticle].Data()),1000,800);
307     cgaus->Divide(8,4,0.001,0.001);
308     Float_t sigma=0;
309     
310     for(Int_t i=0; i<nbins; i++){
311       if(i<MinBin[iparticle])continue;
312       
313       cgaus->cd(i+1);
314       gPad->SetLogy();
315       fHistDCA[i]->SetLineColor(color[iparticle]);
316       fHistDCA[i]->SetMarkerColor(color[iparticle]);
317       fPar->SetLineColor(color[iparticle]);
318       fHistDCA[i]->SetFillColor(16);
319       Printf("first fit step [-1,1]");
320       fHistDCA[i]->Fit(fPar,"NM","",-1,1);
321       Printf("second fit step [-2sigma,2sigma]");
322       Float_t nsigmas=2.;
323       sigma=fPar->GetParameter(2);
324       fHistDCA[i]->Fit(fPar,"NM","",fPar->GetParameter(1)-nsigmas*sigma,fPar->GetParameter(1)+nsigmas*sigma);
325       Printf("Third fit step [-0.5 sigma,0.5 sigma]");
326       nsigmas=.5; // narrow range, just want to fit the primaries
327       fHistDCA[i]->Fit(fPar,"NM","",fPar->GetParameter(1)-nsigmas*sigma,fPar->GetParameter(1)+nsigmas*sigma);
328       
329       //set range to 3 sigmas (for the plot)
330       sigma=fPar->GetParameter(2);
331       fHistDCA[i]->GetXaxis()->SetRangeUser(fPar->GetParameter(1)-3*sigma,fPar->GetParameter(1)+3*sigma);
332       fHistDCA[i]->DrawClone();
333       fPar->DrawClone("same");
334       fHistImpParRes->Fill((xbins[i]+xbins[i+1])/2,10000*fPar->GetParameter(2));
335       fHistImpParRes->SetBinError(fHistImpParRes->FindBin((xbins[i]+xbins[i+1])/2),10000*fPar->GetParError(2));
336       fHistImpParMean->Fill((xbins[i]+xbins[i+1])/2,10000*fPar->GetParameter(1));
337       fHistImpParMean->SetBinError(fHistImpParMean->FindBin((xbins[i]+xbins[i+1])/2),10000*fPar->GetParError(1));
338     }
339     
340     fHistImpParRes->SetMaximum(1000);
341     fHistImpParRes->SetMinimum(0);
342     fHistImpParMean->SetMaximum(80);
343     fHistImpParMean->SetMinimum(-80);
344     
345     setDrawAtt(iparticle+20,color[iparticle],1,color[iparticle],1,fHistImpParRes);
346     fHistImpParRes->SetTitle(Form("%s - Gaus Fit",particle[iparticle].Data()));
347     setDrawAtt(iparticle+20,color[iparticle],1,color[iparticle],1,fHistImpParMean);
348     fHistImpParMean->SetTitle(Form("%s - Gaus Fit",particle[iparticle].Data()));
349     
350     cImpPar->cd();
351     if(iparticle==0)fHistImpParRes->DrawCopy("p");
352     else fHistImpParRes->DrawCopy("psame");
353     
354     cImpParMean->cd();
355     if(iparticle==0)fHistImpParMean->DrawCopy("p");
356     else fHistImpParMean->DrawCopy("psame");
357     
358     TFile *fout=new TFile(foutName.Data(),"update");
359     cgaus->Write();
360     fHistImpParRes->Write();
361     fHistImpParMean->Write();
362     fout->Close();
363     delete fout;
364   }
365   
366   cImpPar->BuildLegend()->SetFillColor(0);
367   cImpParMean->BuildLegend()->SetFillColor(0);
368   TFile *fout=new TFile(foutName.Data(),"update");
369   cImpPar->Write();
370   cImpParMean->Write();
371   fout->Close();
372   delete fout;
373   
374 }
375
376 //----------------------------------------------------------------------------
377
378 void ImpParResGausTail(TList *li,TString foutName)    
379 {
380   //Plot the resolution on the tranverse component of impact parmeter from the ITSsaTracks QA output as a function of pt
381   // Gaussian + tail fit is used
382  
383   TH2F *hd0rphiITSpureSA[3];
384   //binning
385   const Int_t nbins = 29;
386   Double_t xbins[nbins+1]={0.06,0.08,0.10,0.12,0.14,0.16,0.18,0.20,0.25,0.30,
387                            0.35,0.40,0.45,0.50,0.55,0.60,0.65,0.70,0.75,0.80,
388                            0.85,0.90,0.95,1.00,1.20,1.40,1.60,1.80,1.90,2.00};
389   TH1F *d0AllpointrphiSkip1_[nbins];  
390   
391   TCanvas *cImpPar=new TCanvas("ImpParResGausTail_vs_pt","ImpParResGausTail_vs_pt",1000,800);
392   cImpPar->SetGridx();
393   cImpPar->SetLogx();
394   TCanvas *cImpParMean=new TCanvas("ImpParMeanGausTail_vs_pt","ImpParMeanGausTail_vs_pt",1000,800);
395   cImpParMean->SetGridx();
396   cImpParMean->SetLogx();
397   
398   
399   //define the histgram
400   TH1F **d0AllpointrphiSkipTail1_=new TH1F*[nbins]; 
401   TH1F **d0AllpointrphiSkipGaus1_=new TH1F*[nbins];  
402   TH1F **d0Pt_=new TH1F*[nbins];
403   Float_t sigma=0;
404   Double_t j =3.;  
405   
406   for(Int_t iparticle=0;iparticle<3;iparticle++){
407     hd0rphiITSpureSA[iparticle]=(TH2F*)li->FindObject(Form("hd0rphiITSpureSA%s",particle[iparticle].Data()));
408     
409     TH1F *fHistImpParRes = new TH1F(Form("fHistImpParResGausTail%s",particle[iparticle].Data()),Form("%s",particle[iparticle].Data()),nbins,xbins);
410     fHistImpParRes->GetXaxis()->SetTitle("#font[52]{p}_{T} (GeV/#font[52]{c})");
411     fHistImpParRes->GetYaxis()->SetTitle("d0 r#phi resolution [#mum]");
412     fHistImpParRes->GetYaxis()->SetTitleOffset(1.5);
413     TH1F *fHistImpParMean = new TH1F(Form("fHistImpParMeanGausTail%s",particle[iparticle].Data()),Form("%s",particle[iparticle].Data()),nbins,xbins);
414     fHistImpParMean->GetXaxis()->SetTitle("#font[52]{p}_{T} (GeV/#font[52]{c})");
415     fHistImpParMean->GetYaxis()->SetTitle("d0 r#phi mean [#mum]");
416     fHistImpParMean->GetYaxis()->SetTitleOffset(1.5);
417     
418     for(Int_t m=0;m<nbins;m++){
419       d0AllpointrphiSkip1_[m]= (TH1F*)hd0rphiITSpureSA[iparticle]->ProjectionY(Form("%s%i",particle[iparticle].Data(),m),hd0rphiITSpureSA[iparticle]->GetXaxis()->FindBin(xbins[m]+0.000001),hd0rphiITSpureSA[iparticle]->GetXaxis()->FindBin(xbins[m+1]-0.000001));
420     }
421     
422     TCanvas *cgaus=new TCanvas(Form("cgausTailFit_ImpParRes_vs_pt_%s",particle[iparticle].Data()),Form("cgausTailFit_ImpParRes_vs_pt_%s",particle[iparticle].Data()),1000,800);
423     cgaus->Divide(8,4,0.001,0.001);
424     
425     for(Int_t i=0; i<nbins; i++){
426       cgaus->cd(i+1);
427       gPad->SetLogy();
428       d0AllpointrphiSkip1_[i]->SetLineColor(1);
429       d0AllpointrphiSkip1_[i]->SetMarkerColor(1);
430       TF1 *h = new TF1("h","gaus",-1,1);
431       d0AllpointrphiSkip1_[i]->Fit(h,"NM","",-1,1);
432       Double_t d0rphirange_allpointSkip1 = h->GetParameter(2);   
433       Double_t d0rphiMean_allpointSkip1 = h->GetParameter(1);   
434       Double_t cutleft1= -j*(d0rphirange_allpointSkip1);
435       Double_t cutright1 =j*d0rphirange_allpointSkip1;
436       
437       //fitting
438       d0AllpointrphiSkipTail1_[i] = new TH1F(*d0AllpointrphiSkip1_[i]); 
439       d0AllpointrphiSkipGaus1_[i] = new TH1F(*d0AllpointrphiSkip1_[i]); 
440       d0AllpointrphiSkipTail1_[i]->Reset(0);
441       d0AllpointrphiSkipGaus1_[i]->Reset(0);
442       
443       //Filling only with tail and Gaus
444       for (Int_t bin=1;bin<d0AllpointrphiSkip1_[i]->GetNbinsX();bin++){
445         Float_t bincenter = d0AllpointrphiSkip1_[i]->GetBinCenter(bin);
446         if(bincenter<cutleft1 || bincenter>cutright1) {
447           d0AllpointrphiSkipTail1_[i]->SetBinContent(bin,d0AllpointrphiSkip1_[i]->GetBinContent(bin));
448           d0AllpointrphiSkipTail1_[i]->SetBinError(bin,d0AllpointrphiSkip1_[i]->GetBinError(bin));
449           d0AllpointrphiSkipGaus1_[i]->SetBinContent(bin,0.);  
450           //This sentence is very important,otherwise we will get the information the data is empty when we fit it .
451         }
452         else if(bincenter>=cutleft1 && bincenter<=cutright1){
453           d0AllpointrphiSkipTail1_[i]->SetBinContent(bin,0);
454           d0AllpointrphiSkipGaus1_[i]->SetBinContent(bin,d0AllpointrphiSkip1_[i]->GetBinContent(bin));
455           d0AllpointrphiSkipGaus1_[i]->SetBinError(bin,d0AllpointrphiSkip1_[i]->GetBinError(bin));  
456         }
457       }
458       d0AllpointrphiSkipGaus1_[i]->SetLineColor(2);
459       d0AllpointrphiSkipGaus1_[i]->SetMarkerColor(2);
460       d0AllpointrphiSkipTail1_[i]->SetLineColor(4);
461       d0AllpointrphiSkipTail1_[i]->SetMarkerColor(4);
462       
463       TF1 *hh;
464       hh =CreateFuncTail(d0AllpointrphiSkipTail1_[i],"hh");
465       hh->SetLineColor(d0AllpointrphiSkipTail1_[i]->GetLineColor());
466       d0AllpointrphiSkipTail1_[i]->Fit(hh,"NM","",-1,1);  
467       Double_t Sigmatail_allpointSkip1 = hh->GetParameter(2);
468       d0AllpointrphiSkipGaus1_[i]->Fit(h,"NM","",d0rphiMean_allpointSkip1-d0rphirange_allpointSkip1,d0rphiMean_allpointSkip1+d0rphirange_allpointSkip1);//narrow around the peak
469       //d0AllpointrphiSkipGaus1_[i]->Fit(h,"NM","",-1,1);
470       h->SetLineColor(d0AllpointrphiSkipGaus1_[i]->GetLineColor());
471       Double_t Sigmagaus_allpointSkip1 = h->GetParameter(2);
472       Double_t Meangaus_allpointSkip1 = h->GetParameter(1);
473       Double_t Constgaus_allpointSkip1 = h->GetParameter(0);
474       d0AllpointrphiSkipGaus1_[i]->DrawClone("");
475       d0AllpointrphiSkipTail1_[i]->DrawClone("same");
476       hh->DrawClone("same");
477       h->DrawClone("same");
478       
479       TF1 * fPar=CreateFuncGaussTail(d0AllpointrphiSkip1_[i],"allpointSkip1",Constgaus_allpointSkip1,Meangaus_allpointSkip1,Sigmagaus_allpointSkip1,Sigmatail_allpointSkip1);
480       d0AllpointrphiSkip1_[i]->Fit(fPar,"NM","",-1,1);
481       fPar->SetLineColor(d0AllpointrphiSkip1_[i]->GetLineColor());
482       fPar->DrawClone("same");
483       
484       sigma=fPar->GetParameter(2);
485       fHistImpParRes->Fill((xbins[i]+xbins[i+1])/2,10000*fPar->GetParameter(2));
486       fHistImpParRes->SetBinError(fHistImpParRes->FindBin((xbins[i]+xbins[i+1])/2),10000*fPar->GetParError(2));
487       fHistImpParMean->Fill((xbins[i]+xbins[i+1])/2,10000*fPar->GetParameter(1));
488       fHistImpParMean->SetBinError(fHistImpParMean->FindBin((xbins[i]+xbins[i+1])/2),10000*fPar->GetParError(1));
489     }
490     
491     fHistImpParRes->SetMaximum(1000);
492     fHistImpParRes->SetMinimum(0);
493     fHistImpParMean->SetMaximum(80);
494     fHistImpParMean->SetMinimum(-80);
495     
496     setDrawAtt(iparticle+20,color[iparticle],1,color[iparticle],1,fHistImpParRes);
497     fHistImpParRes->SetTitle(Form("%s - GausTail Fit",particle[iparticle].Data()));
498     setDrawAtt(iparticle+20,color[iparticle],1,color[iparticle],1,fHistImpParMean);
499     fHistImpParMean->SetTitle(Form("%s - GausTail Fit",particle[iparticle].Data()));
500    
501     cImpPar->cd();
502     if(iparticle==0)fHistImpParRes->DrawCopy("p");
503     else fHistImpParRes->DrawCopy("psame");
504     
505     cImpParMean->cd();
506     if(iparticle==0)fHistImpParMean->DrawCopy("p");
507     else fHistImpParMean->DrawCopy("psame");
508     
509     
510     TFile *fout=new TFile(foutName.Data(),"update");
511     cgaus->Write();
512     fHistImpParRes->Write();
513     fHistImpParMean->Write();
514     fout->Close();
515     delete fout;
516   }
517   
518   cImpPar->BuildLegend()->SetFillColor(0);
519   cImpParMean->BuildLegend()->SetFillColor(0);
520   TFile *fout=new TFile(foutName.Data(),"update");
521   cImpPar->Write();
522   cImpParMean->Write();
523   fout->Close();
524   delete fout;
525   
526 }
527 //----------------------------------------------------------------------------
528
529 TF1 *CreateFuncTail(TH1F *hh,TString funcname,Double_t wholeRangeInt=-1.)
530 {
531   TF1 *tail=new TF1(funcname.Data(),"[0]*(1./(2.*[2])*TMath::Exp(-TMath::Abs(x-[1])/[2]))",-1,1);
532   Double_t binwidth=hh->GetBinWidth(10);
533   Double_t integral=hh->Integral();
534   if(wholeRangeInt!=-1.)tail->SetParLimits(0,(0.2)*wholeRangeInt*binwidth,(0.5)*wholeRangeInt*binwidth);
535   Double_t RMS1=TMath::Abs(hh->GetRMS());
536   Double_t firstvalue1=binwidth*integral;
537   tail->SetParameters(1.,0,100.);//Set the initial value of parameter
538   return tail;
539
540
541 //----------------------------------------------------------------------------
542 TF1 *CreateFuncGaussTail(TH1F *h,TString funcname,Double_t Norm,Double_t parzero,Double_t parone,Double_t partwo)
543 {
544   TF1 *gaustail=new TF1(funcname.Data(),"[0]*([4]/(TMath::Sqrt(2.*TMath::Pi())*[2])*TMath::Exp(-1.*(x-[1])*(x-[1])/(2.*[2]*[2]))+(1.-[4])/(2.*[3])*TMath::Exp(-TMath::Abs(x-[1])/[3]))",-1,1);
545   gaustail->SetParameters(Norm,parzero,parone,partwo,0.82);//Set the initial value of parameter
546   return gaustail;
547
548
549
550 //----------------------------------------------------------------------------
551
552 void CompareImpParRes(TString foutName){
553   //Copare the resolution obtained from Gaus and GausTail fit!
554   // BE CAREFUL!! look at the fits "by eye", GausTail sometimes fail, especially pions at low pt
555   
556   TFile *fout=new TFile(foutName.Data(),"read");
557   
558   TCanvas *c=new TCanvas("CompareImpParRes","CompareImpParRes",1);
559   c->SetLogx();
560   
561   for(Int_t iparticle=0;iparticle<3;iparticle++){
562     TH1F *fHistImpParResGaus = (TH1F*)fout->Get(Form("fHistImpParResGaus%s",particle[iparticle].Data()));
563     TH1F *fHistImpParResGausTail = (TH1F*)fout->Get(Form("fHistImpParResGausTail%s",particle[iparticle].Data()));
564     setDrawAtt(iparticle+24,color[iparticle],1,color[iparticle],1,fHistImpParResGaus);
565     if(iparticle==0)fHistImpParResGaus->DrawClone();
566     else fHistImpParResGaus->DrawClone("same");
567     fHistImpParResGausTail->DrawClone("same");
568   }
569   gPad->BuildLegend()->SetFillColor(0);
570   
571   TFile *fout=new TFile(foutName.Data(),"update");
572   c->Write();
573   fout->Close();
574   delete fout;
575   
576 }
577
578
579 //-----------------------------------------------------
580
581 void PlotPtResol(TList* l,TString foutName, Bool_t optFromMC){
582   TString hNameR,hNameA;
583   TString partName[3]={"Pion","Kaon","Proton"};
584   TString prefix;
585   if(optFromMC) prefix="hMC";
586   else prefix="h";
587   
588   TCanvas* c2d[3];
589   TCanvas* c1dA[3];
590   TCanvas* c1dR[3];
591   
592   TH2F* h2DA[3];
593   TH2F* h2DR[3];
594   TH1F* hptres[3][40];
595   TH1F* h1ptrelres[3][40];
596   TH1F* hptreco[3][40];
597   
598   TGraphErrors* gbias[3];
599   TGraphErrors* grelresol[3];
600   
601   gStyle->SetPalette(1);
602   
603   for(Int_t iSpec=0; iSpec<3; iSpec++){
604     hNameA=Form("%sPtResid%s",prefix.Data(),partName[iSpec].Data());
605     hNameR=Form("%sInvPtRelResid%s",prefix.Data(),partName[iSpec].Data());
606     printf("%s %s\n",hNameA.Data(),hNameR.Data());
607     h2DA[iSpec]=(TH2F*)l->FindObject(hNameA.Data());
608     h2DR[iSpec]=(TH2F*)l->FindObject(hNameR.Data());
609     c2d[iSpec]=new TCanvas(Form("c2d%s",partName[iSpec].Data()),Form("c2d%s",partName[iSpec].Data()));
610     c2d[iSpec]->Divide(2,1);
611     c2d[iSpec]->cd(1);
612     h2DA[iSpec]->Draw("colz"); 
613     c2d[iSpec]->cd(2);
614     h2DR[iSpec]->Draw("colz");
615
616     Int_t nptbins=h2DR[iSpec]->GetNbinsX();
617
618     Int_t nybinsA=h2DA[iSpec]->GetNbinsY();
619     Float_t minyA=h2DA[iSpec]->GetYaxis()->GetBinLowEdge(1);
620     Float_t maxyA=h2DA[iSpec]->GetYaxis()->GetBinUpEdge(nybinsA);
621
622     Int_t nybinsR=h2DR[iSpec]->GetNbinsY();
623     Float_t minyR=h2DR[iSpec]->GetYaxis()->GetBinLowEdge(1);
624     Float_t maxyR=h2DR[iSpec]->GetYaxis()->GetBinUpEdge(nybinsR);
625     printf("%d   %d %f %f   %d %f %f\n",nptbins,nybinsA,minyA,maxyA,nybinsR,minyR,maxyR);
626
627     c1dA[iSpec]=new TCanvas(Form("c1dA%s",partName[iSpec].Data()),Form("c1dA%s",partName[iSpec].Data()));
628     c1dA[iSpec]->Divide(6,5);
629     c1dR[iSpec]=new TCanvas(Form("c1dR%s",partName[iSpec].Data()),Form("c1dR%s",partName[iSpec].Data()));
630     c1dR[iSpec]->Divide(6,5);
631
632
633     gbias[iSpec]=new TGraphErrors(0);
634     grelresol[iSpec]=new TGraphErrors(0);
635     gbias[iSpec]->SetTitle("");
636     grelresol[iSpec]->SetTitle("");
637
638     for(Int_t iptbin=0; iptbin<nptbins;iptbin++){
639       Float_t avept=h2DA[iSpec]->GetXaxis()->GetBinCenter(iptbin+1);
640       Float_t widpt=0.5*h2DA[iSpec]->GetXaxis()->GetBinWidth(iptbin+1);
641       Int_t minptbinmev=(Int_t)(h2DA[iSpec]->GetXaxis()->GetBinLowEdge(iptbin+1)*1000.+0.5);
642
643       hptres[iSpec][iptbin]=new TH1F(Form("hptres%s_%d",partName[iSpec].Data(),minptbinmev),
644                                      Form("hptres%s_%d",partName[iSpec].Data(),minptbinmev),
645                                      nybinsA,minyA,maxyA);
646       h1ptrelres[iSpec][iptbin]=new TH1F(Form("h1ptrelres%s_%d",partName[iSpec].Data(),minptbinmev),
647                                          Form("h1ptrelres%s_%d",partName[iSpec].Data(),minptbinmev),
648                                          nybinsR,minyR,maxyR);
649       hptreco[iSpec][iptbin]=new TH1F(Form("hptreco%s_%d",partName[iSpec].Data(),minptbinmev),
650                                       Form("hptreco%s_%d",partName[iSpec].Data(),minptbinmev),
651                                       400,0.,2.);
652       for(Int_t iBin=1; iBin<=nybinsA; iBin++){
653         hptres[iSpec][iptbin]->SetBinContent(iBin,h2DA[iSpec]->GetBinContent(iptbin+1,iBin));
654         hptres[iSpec][iptbin]->SetBinError(iBin,h2DA[iSpec]->GetBinError(iptbin+1,iBin));
655       }
656       for(Int_t iBin=1; iBin<=nybinsR; iBin++){
657         h1ptrelres[iSpec][iptbin]->SetBinContent(iBin,h2DR[iSpec]->GetBinContent(iptbin+1,iBin));
658         h1ptrelres[iSpec][iptbin]->SetBinError(iBin,h2DR[iSpec]->GetBinError(iptbin+1,iBin));
659       }
660
661       c1dA[iSpec]->cd(iptbin+1);
662       hptres[iSpec][iptbin]->Draw();
663       if(hptres[iSpec][iptbin]->Integral()>50){
664         hptres[iSpec][iptbin]->Fit("gaus");
665         hptres[iSpec][iptbin]->GetXaxis()->SetTitle("Pt residuals (GeV/c)");
666         hptres[iSpec][iptbin]->GetXaxis()->CenterTitle();
667         TF1* fgaus= (TF1*)hptres[iSpec][iptbin]->GetListOfFunctions()->FindObject("gaus");
668         Int_t nPoint=gbias[iSpec]->GetN();
669         gbias[iSpec]->SetPoint(nPoint, avept, fgaus->GetParameter(1));//hptres[iSpec][iptbin]->GetMean());
670         gbias[iSpec]->SetPointError(nPoint, widpt, fgaus->GetParError(1)); //hptres[iSpec][iptbin]->GetMeanError());
671       }
672       c1dR[iSpec]->cd(iptbin+1);
673       h1ptrelres[iSpec][iptbin]->Draw();
674       if(h1ptrelres[iSpec][iptbin]->Integral()>50){
675         h1ptrelres[iSpec][iptbin]->Fit("gaus","","",-0.1,0.1);//,"L");
676         h1ptrelres[iSpec][iptbin]->GetXaxis()->SetTitle("1/Pt relative residuals");
677         h1ptrelres[iSpec][iptbin]->GetXaxis()->CenterTitle();
678         TF1* fgaus= (TF1*)h1ptrelres[iSpec][iptbin]->GetListOfFunctions()->FindObject("gaus");    
679         Int_t nPoint=grelresol[iSpec]->GetN();
680         grelresol[iSpec]->SetPoint(nPoint, avept, fgaus->GetParameter(2));
681         grelresol[iSpec]->SetPointError(nPoint, widpt, fgaus->GetParError(2));    
682       }
683
684     }
685   }
686   
687
688   TCanvas* cb=new TCanvas("cb","Bias");
689   gbias[2]->SetMarkerStyle(22);
690   gbias[2]->SetMarkerColor(4);
691   gbias[2]->SetLineColor(4);
692   gbias[2]->Draw("PA");
693   gbias[0]->SetMarkerStyle(20);
694   gbias[0]->SetMarkerColor(1);
695   gbias[0]->SetLineColor(1);
696   gbias[0]->Draw("PSAME");
697   gbias[1]->SetMarkerStyle(25);
698   gbias[1]->SetMarkerColor(2);
699   gbias[1]->SetLineColor(2);
700   gbias[1]->Draw("PSAME");
701   gbias[2]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
702   if(optFromMC) gbias[2]->GetYaxis()->SetTitle("<p_{T}(ITSsa)-p_{T}(MC)> (GeV/c)");
703   else gbias[2]->GetYaxis()->SetTitle("<p_{T}(ITSsa)-p_{T}(TPCITS)> (GeV/c)");
704   gbias[2]->GetYaxis()->SetTitleOffset(1.2);
705   cb->Update();
706
707   TCanvas* cr=new TCanvas("cr","Resol");
708   grelresol[2]->SetMinimum(0.);
709   grelresol[2]->SetMaximum(0.2);
710   grelresol[2]->SetMarkerStyle(22);
711   grelresol[2]->SetMarkerColor(4);
712   grelresol[2]->SetLineColor(4);
713   grelresol[2]->Draw("PA");
714   grelresol[0]->SetMarkerStyle(20);
715   grelresol[0]->SetMarkerColor(1);
716   grelresol[0]->SetLineColor(1);
717   grelresol[0]->Draw("PSAME");
718   grelresol[1]->SetMarkerStyle(25);
719   grelresol[1]->SetMarkerColor(2);
720   grelresol[1]->SetLineColor(2);
721   grelresol[1]->Draw("PSAME");
722   grelresol[2]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
723   grelresol[2]->GetYaxis()->SetTitle("1/Pt relative resolution (%)");
724   grelresol[2]->GetYaxis()->SetTitleOffset(1.2);
725   cr->Update();
726  
727   TFile *fout=new TFile(foutName.Data(),"update");
728   for(Int_t iSpec=0; iSpec<3; iSpec++){
729     c2d[iSpec]->Write();
730     c1dA[iSpec]->Write();
731     c1dR[iSpec]->Write();
732   }
733   cb->Write();
734   cr->Write();
735   fout->Close();
736   delete fout;
737 }
738 //----------------------------------------------------------------------------
739
740 void PlotITSsa(TList* li,TString foutName){
741
742   TH1F* hPtTPCITS=(TH1F*)li->FindObject("hPtTPCITS");
743   TH1F* hPtITSsa=(TH1F*)li->FindObject("hPtITSsa");
744   TH1F* hPtITSpureSA=(TH1F*)li->FindObject("hPtITSpureSA");
745
746   TH2F* hEtaPhiTPCITS=(TH2F*)li->FindObject("hEtaPhiTPCITS");
747   TH2F* hEtaPhiITSsa=(TH2F*)li->FindObject("hEtaPhiITSsa");
748   TH2F* hEtaPhiITSpureSA=(TH2F*)li->FindObject("hEtaPhiITSpureSA");
749
750   TH1F* hChi2TPCITS=(TH1F*)li->FindObject("hChi2TPCITS");
751   TH1F* hChi2ITSsa=(TH1F*)li->FindObject("hChi2ITSsa");
752   TH1F* hChi2ITSpureSA=(TH1F*)li->FindObject("hChi2ITSpureSA");
753
754   TH1F* hNcluTPCITS=(TH1F*)li->FindObject("hNcluTPCITS");
755   TH1F* hNcluITSsa=(TH1F*)li->FindObject("hNcluITSsa");
756   TH1F* hNcluITSpureSA=(TH1F*)li->FindObject("hNcluITSpureSA");
757
758
759   TH1F* hRatio=(TH1F*)hPtTPCITS->Clone("hRatio");
760   hRatio->Add(hPtITSsa);
761   hRatio->Divide(hPtITSpureSA);
762   hRatio->SetStats(0);
763
764   TCanvas* c1=new TCanvas("c1","Pt",800,1000);
765   c1->Divide(1,2);
766   c1->cd(1);
767   hPtITSpureSA->Draw();
768   hPtITSpureSA->GetXaxis()->SetTitle("Pt (GeV/c)");
769   gPad->Update();
770   TPaveStats *st1=(TPaveStats*)hPtITSpureSA->GetListOfFunctions()->FindObject("stats");
771   st1->SetY1NDC(0.71);
772   st1->SetY2NDC(0.9);
773   hPtTPCITS->SetLineColor(2);
774   hPtTPCITS->Draw("sames");
775   gPad->Update();
776   TPaveStats *st2=(TPaveStats*)hPtTPCITS->GetListOfFunctions()->FindObject("stats");
777   st2->SetY1NDC(0.51);
778   st2->SetY2NDC(0.7);
779   st2->SetTextColor(2);
780
781   hPtITSsa->SetLineColor(4);
782   hPtITSsa->Draw("sames");
783   gPad->Update();
784   TPaveStats *st3=(TPaveStats*)hPtITSsa->GetListOfFunctions()->FindObject("stats");
785   st3->SetY1NDC(0.31);
786   st3->SetY2NDC(0.5);
787   st3->SetTextColor(4);
788   TLegend* leg=new TLegend(0.5,0.5,0.69,0.69);
789   leg->SetFillColor(0);
790   TLegendEntry* ent=leg->AddEntry(hPtITSpureSA,"ITS pureSA","L");
791   ent->SetTextColor(hPtITSpureSA->GetLineColor());
792   ent=leg->AddEntry(hPtTPCITS,"TPC+ITS","L");
793   ent->SetTextColor(hPtTPCITS->GetLineColor());
794   ent=leg->AddEntry(hPtITSsa,"ITSsa","L");
795   ent->SetTextColor(hPtITSsa->GetLineColor());
796   leg->Draw();
797   c1->cd(2);
798   gPad->SetGridx();
799   gPad->SetGridy();
800   hRatio->Draw();
801   hRatio->GetXaxis()->SetTitle("Pt (GeV/c)");
802   hRatio->GetYaxis()->SetTitle("(TPCITS+ITSsa)/ITSpureSA");
803
804   hChi2ITSpureSA->Scale(1./hChi2ITSpureSA->GetEntries());
805   hChi2ITSsa->Scale(1./hChi2ITSsa->GetEntries());
806   hChi2TPCITS->Scale(1./hChi2TPCITS->GetEntries());
807
808   TCanvas* c2=new TCanvas("c2","Chi2");
809   hChi2ITSpureSA->Draw();
810   hChi2ITSpureSA->GetXaxis()->SetTitle("Chi2");
811   gPad->Update();
812   TPaveStats *stc1=(TPaveStats*)hChi2ITSpureSA->GetListOfFunctions()->FindObject("stats");
813   stc1->SetY1NDC(0.71);
814   stc1->SetY2NDC(0.9);
815   hChi2TPCITS->SetLineColor(2);
816   hChi2TPCITS->Draw("sames");
817   gPad->Update();
818   TPaveStats *stc2=(TPaveStats*)hChi2TPCITS->GetListOfFunctions()->FindObject("stats");
819   stc2->SetY1NDC(0.51);
820   stc2->SetY2NDC(0.7);
821   stc2->SetTextColor(2);
822   c2->Update();
823   hChi2ITSsa->SetLineColor(4);
824   hChi2ITSsa->Draw("sames");
825   gPad->Update();
826   TPaveStats *stc3=(TPaveStats*)hChi2ITSsa->GetListOfFunctions()->FindObject("stats");
827   stc3->SetY1NDC(0.31);
828   stc3->SetY2NDC(0.5);
829   stc3->SetTextColor(4);
830   leg->Draw();
831
832   hNcluITSpureSA->Scale(1./hNcluITSpureSA->GetEntries());
833   hNcluITSsa->Scale(1./hNcluITSsa->GetEntries());
834   hNcluTPCITS->Scale(1./hNcluTPCITS->GetEntries());
835
836   TCanvas* c3=new TCanvas("c3","Nclu");
837   c3->SetRightMargin(0.22);
838   hNcluITSpureSA->Draw();
839   hNcluITSpureSA->GetXaxis()->SetTitle("n. ITS clusters");
840   gPad->Update();
841   TPaveStats *stn1=(TPaveStats*)hNcluITSpureSA->GetListOfFunctions()->FindObject("stats");
842   stn1->SetY1NDC(0.71);
843   stn1->SetY2NDC(0.9);
844   hNcluTPCITS->SetLineColor(2);
845   hNcluTPCITS->Draw("sames");
846   gPad->Update();
847   TPaveStats *stn2=(TPaveStats*)hNcluTPCITS->GetListOfFunctions()->FindObject("stats");
848   stn2->SetY1NDC(0.51);
849   stn2->SetY2NDC(0.7);
850   stn2->SetTextColor(2);
851
852   hNcluITSsa->SetLineColor(4);
853   hNcluITSsa->Draw("sames");
854   gPad->Update();
855   TPaveStats *stn3=(TPaveStats*)hNcluITSsa->GetListOfFunctions()->FindObject("stats");
856   stn3->SetY1NDC(0.31);
857   stn3->SetY2NDC(0.5);
858   stn3->SetTextColor(4);
859   leg->Draw();
860
861   gStyle->SetPalette(1);
862   hEtaPhiITSpureSA->SetStats(0);
863   hEtaPhiITSpureSA->SetTitle("ITS pureSA");
864   hEtaPhiITSsa->SetStats(0);
865   hEtaPhiITSsa->SetTitle("ITSsa");
866   hEtaPhiTPCITS->SetStats(0);
867   hEtaPhiTPCITS->SetTitle("TPC+ITS");
868   TCanvas* c4=new TCanvas("c4","EtaPhi",1000,700);
869   c4->Divide(3,1);
870   c4->cd(1);
871   hEtaPhiITSpureSA->Draw("colz");
872   hEtaPhiITSpureSA->GetXaxis()->SetTitle("Eta");
873   hEtaPhiITSpureSA->GetYaxis()->SetTitle("Phi");
874   c4->cd(2);
875   hEtaPhiITSsa->Draw("colz");
876   hEtaPhiITSsa->GetXaxis()->SetTitle("Eta");
877   hEtaPhiITSsa->GetYaxis()->SetTitle("Phi");
878   c4->cd(3);
879   hEtaPhiTPCITS->Draw("colz");
880   hEtaPhiTPCITS->GetXaxis()->SetTitle("Eta");
881   hEtaPhiTPCITS->GetYaxis()->SetTitle("Phi");
882   
883   TFile *fout=new TFile(foutName.Data(),"update");
884   c1->Write();
885   c2->Write();
886   c3->Write();
887   c4->Write();
888   fout->Close();
889   delete fout;
890  
891 }
892
893
894 //----------------------------------------------------------------------------
895
896 void setDrawAtt(Int_t markerstyle,Int_t markercolor,Int_t markersize,Int_t linecolor,Int_t linewidth,TH1 *h1)
897
898   h1->SetMarkerStyle(markerstyle);
899   h1->SetMarkerColor(markercolor);
900   h1->SetMarkerSize(markersize);
901   h1->SetLineColor(linecolor);
902   h1->SetLineWidth(linewidth);
903 }