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