Primary charged tracks added
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / TestAOD / MainAnalysis.C
1 #include <TString.h>
2 #include <TStyle.h>
3 #include <TSystem.h>
4 #include <TFile.h>
5 #include <iostream>
6 #include "TFile.h"
7 #include "TH1I.h"
8 #include "TH2.h"
9 #include "TCanvas.h"
10 #include "TLegend.h"
11 #include "TROOT.h"
12
13 TString Charge[]={"Pos","Neg"};
14 TString Sign[]={"Plus","Minus"};
15 TString Particle[]={"Pion","Kaon","Proton"};
16 Int_t Color[3]={1,2,4};
17 Int_t Marker[6]={20,21,22,24,25,26};
18 Double_t projl[3]={0,20,80};
19 Double_t proju[3]={10,40,90};
20
21
22 void MainAnalysis() {
23
24   //TString fold="5SigmaPID_AOD086";
25   TString fold="5SigmaPID_AOD046";
26   Int_t ibinToCompare=-1;
27   
28   TString sname="Cent0to5_QVec0.0to100.0";ibinToCompare=0;
29   
30   //TString sname="Cent5to10_QVec0.0to100.0";ibinToCompare=1;
31   //TString sname="Cent10to20_QVec0.0to100.0";ibinToCompare=2;
32   //TString sname="Cent20to30_QVec0.0to100.0";ibinToCompare=3;
33   //TString sname="Cent30to40_QVec0.0to100.0";ibinToCompare=4;
34   
35   TString dataFile = Form("output/%s/Pt.AOD.1._data_ptcut_%s.root",fold.Data(),sname.Data());
36   TString mcFile =Form("output/%s/Pt.AOD.1._MC_%s.root",fold.Data(),sname.Data());
37   
38   gSystem->Load("libCore.so");  
39   gSystem->Load("libGeom.so");
40   gSystem->Load("libPhysics.so");
41   gSystem->Load("libVMC");
42   gSystem->Load("libTree");
43   gSystem->Load("libProof");
44   gSystem->Load("libMatrix");
45   gSystem->Load("libSTEERBase");
46   gSystem->Load("libESD");
47   gSystem->Load("libAOD");
48   gSystem->Load("libANALYSIS");
49   gSystem->Load("libOADB");
50   gSystem->Load("libANALYSISalice");
51   gSystem->Load("libTENDER");
52   gSystem->Load("libCORRFW");
53   //gSystem->Load("libPWG0base");
54   gSystem->Load("libMinuit");
55   gSystem->Load("libPWGTools");
56   gSystem->Load("libPWGLFSPECTRA");
57   gSystem->AddIncludePath("-I$ALICE_ROOT/include");
58   gStyle->SetPalette(1);
59   
60   // Open root MC file and get classes
61   cout << "Analysis Macro" << endl;
62   cout << "  > Reading MC data" << endl;
63   TFile *_mc = TFile::Open(mcFile.Data());
64   AliSpectraAODHistoManager* hman_mc = (AliSpectraAODHistoManager*) _mc->Get("SpectraHistos");
65   AliSpectraAODEventCuts* ecuts_mc = (AliSpectraAODEventCuts*) _mc->Get("Event Cuts");
66   AliSpectraAODTrackCuts* tcuts_mc = (AliSpectraAODTrackCuts*) _mc->Get("Track Cuts");
67   // print info about mc track and Event cuts
68   cout << " -- Info about MC -- "<< endl;
69   ecuts_mc->PrintCuts();
70   // tcuts_mc->PrintCuts();
71   // proceed likewise for data
72   TFile *_data = TFile::Open(dataFile.Data());
73   AliSpectraAODHistoManager* hman_data = (AliSpectraAODHistoManager*) _data->Get("SpectraHistos");
74   AliSpectraAODEventCuts* ecuts_data = (AliSpectraAODEventCuts*) _data->Get("Event Cuts");
75   AliSpectraAODTrackCuts* tcuts_data = (AliSpectraAODTrackCuts*) _data->Get("Track Cuts");
76   // print info about track and Event cuts
77   cout << " -- Info about data -- " << endl;
78   ecuts_data->PrintCuts();
79   tcuts_data->PrintCuts();
80   
81   //dedx in data and MC
82   TCanvas *cPIDSig=new TCanvas("cPIDSig","cPIDSig",700,500);
83   cPIDSig->Divide(2,2);
84   cPIDSig->cd(1);
85   TH2F *PIDSig_data = (TH2F*)((TH2F*)hman_data->GetPIDHistogram("histPIDTPC"))->Clone();
86   gPad->SetLogz();
87   gPad->SetGridy();
88   gPad->SetGridx();
89   PIDSig_data->DrawClone("colz");
90   cPIDSig->cd(2);
91   TH2F *PIDSig_mc = (TH2F*)((TH2F*)hman_mc->GetPIDHistogram("histPIDTPC"))->Clone();
92   gPad->SetLogz();
93   gPad->SetGridy();
94   gPad->SetGridx();
95   PIDSig_mc->DrawClone("colz");
96   cPIDSig->cd(3);
97   TH2F *PIDSig_data = (TH2F*)((TH2F*)hman_data->GetPIDHistogram("histPIDTOF"))->Clone();
98   gPad->SetLogz();
99   gPad->SetGridy();
100   gPad->SetGridx();
101   PIDSig_data->DrawClone("colz");
102   cPIDSig->cd(4);
103   TH2F *PIDSig_mc = (TH2F*)((TH2F*)hman_mc->GetPIDHistogram("histPIDTOF"))->Clone();
104   gPad->SetLogz();
105   gPad->SetGridy();
106   gPad->SetGridx();
107   PIDSig_mc->DrawClone("colz");
108
109   //nsig in data and MC
110   for(Int_t ipart=0;ipart<3;ipart++){
111     TCanvas *cnsig=new TCanvas(Form("cnsig%s",Particle[ipart].Data()),Form("cnsig%s",Particle[ipart].Data()),700,500);
112     cnsig->Divide(2,3);
113     cnsig->cd(1);
114     TH2F *nsig_data = (TH2F*)((TH2F*)hman_data->GetNSigHistogram(Form("histNSig%sPtTPC",Particle[ipart].Data())))->Clone();
115     nsig_data->SetXTitle("Pt (GeV/c)");
116     gPad->SetLogz();
117     gPad->SetGridy();
118     gPad->SetGridx();
119     nsig_data->DrawClone("colz");
120     cnsig->cd(2);
121     TH2F *nsig_mc = (TH2F*)((TH2F*)hman_mc->GetNSigHistogram(Form("histNSig%sPtTPC",Particle[ipart].Data())))->Clone();
122     nsig_mc->SetXTitle("Pt (GeV/c)");
123     gPad->SetLogz();
124     gPad->SetGridy();
125     gPad->SetGridx();
126     nsig_mc->DrawClone("colz");
127     cnsig->cd(3);
128     TH2F *nsig_data = (TH2F*)((TH2F*)hman_data->GetNSigHistogram(Form("histNSig%sPtTOF",Particle[ipart].Data())))->Clone();
129     nsig_data->SetXTitle("Pt (GeV/c)");
130     gPad->SetLogz();
131     gPad->SetGridy();
132     gPad->SetGridx();
133     nsig_data->DrawClone("colz");
134     cnsig->cd(4);
135     TH2F *nsig_mc = (TH2F*)((TH2F*)hman_mc->GetNSigHistogram(Form("histNSig%sPtTOF",Particle[ipart].Data())))->Clone();
136     nsig_mc->SetXTitle("Pt (GeV/c)");
137     gPad->SetLogz();
138     gPad->SetGridy();
139     gPad->SetGridx();
140     nsig_mc->DrawClone("colz");
141     cnsig->cd(5);
142     TH2F *nsig_data = (TH2F*)((TH2F*)hman_data->GetNSigHistogram(Form("histNSig%sPtTPCTOF",Particle[ipart].Data())))->Clone();
143     nsig_data->SetXTitle("Pt (GeV/c)");
144     gPad->SetLogz();
145     gPad->SetGridy();
146     gPad->SetGridx();
147     nsig_data->DrawClone("colz");
148     cnsig->cd(6);
149     TH2F *nsig_mc = (TH2F*)((TH2F*)hman_mc->GetNSigHistogram(Form("histNSig%sPtTPCTOF",Particle[ipart].Data())))->Clone();
150     nsig_mc->SetXTitle("Pt (GeV/c)");
151     gPad->SetLogz();
152     gPad->SetGridy();
153     gPad->SetGridx();
154     nsig_mc->DrawClone("colz");
155   }
156   //qvec in data and MC
157   for (Int_t icharge=0;icharge<2;icharge++){
158     TCanvas *cqvec=new TCanvas(Form("cqvec%s",Charge[icharge].Data()),Form("cqvec%s",Charge[icharge].Data()),700,500);
159     cqvec->Divide(2,1);
160     cqvec->cd(1);
161     TString hname=Form("histq%s",Charge[icharge].Data());
162     TH2F *qvec_data = (TH2F*)((TH2F*)hman_data->GetqVecHistogram(hname.Data()))->Clone();
163     gPad->SetLogz();
164     gPad->SetGridy();
165     gPad->SetGridx();
166     qvec_data->DrawClone("colz");
167     cqvec->cd(2);
168     TH2F *qvec_mc = (TH2F*)((TH2F*)hman_mc->GetqVecHistogram(hname.Data()))->Clone();
169     gPad->SetLogz();
170     gPad->SetGridy();
171     gPad->SetGridx();
172     qvec_mc->DrawClone("colz");
173     
174     TCanvas *cProjqvec=new TCanvas(Form("cProjqvec%s",Charge[icharge].Data()),Form("cProjqvec%s",Charge[icharge].Data()),700,500);
175     cProjqvec->Divide(3,2);
176     for(Int_t iproj=0;iproj<3;iproj++){
177       TH1F *proj_data=(TH1F*)qvec_data->ProjectionX("data",qvec_data->GetYaxis()->FindBin(projl[iproj]),qvec_data->GetYaxis()->FindBin(proju[iproj]));
178       if(proj_data->GetEntries()==0)continue;
179       proj_data->Scale(1/proj_data->GetEntries());
180       proj_data->SetTitle(Form("data q%s [%.0f,%.0f]",Charge[icharge].Data(),projl[iproj],proju[iproj]));
181       TH1F *proj_mc=(TH1F*)qvec_mc->ProjectionX("mc",qvec_mc->GetYaxis()->FindBin(projl[iproj]),qvec_mc->GetYaxis()->FindBin(proju[iproj]));
182       proj_mc->Scale(1/proj_mc->GetEntries());
183       proj_mc->SetTitle(Form("mc q%s [%.0f,%.0f]",Charge[icharge].Data(),projl[iproj],proju[iproj]));
184       proj_mc->SetLineColor(2);
185       cProjqvec->cd(iproj+1);
186       gPad->SetGridy();
187       gPad->SetGridx();
188       proj_data->DrawClone();
189       proj_mc->DrawClone("same");
190       gPad->BuildLegend();
191       proj_data->Divide(proj_mc);
192       cProjqvec->cd(iproj+4);
193       proj_data->DrawClone();
194     }
195   }
196   
197   //efficiencies
198   Printf("\n\n-> Calculating MC Correction Factors");
199   TH1F *CorrFact[6];
200   TCanvas *ceff=new TCanvas("ceff","ceff",700,500);
201   Bool_t UseMCDCACorrection=kTRUE;
202   for(Int_t icharge=0;icharge<2;icharge++){
203     for(Int_t ipart=0;ipart<3;ipart++){
204       Int_t index=ipart+3*icharge;
205       //BE CAREFUL! depending on the efficiency you choose, you must change the DCA correction (data or data/mc)
206       
207       TString hname=Form("histPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data()); //MC correction for Prim+Cont+Eff
208       //TString hname=Form("histPtRecTrue%s%s",Particle[ipart].Data(),Sign[icharge].Data());//MC correction for Prim+Eff
209       //TString hname=Form("histPtRecSigmaPrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());UseMCDCACorrection=kFALSE; //MC correction for Cont+Eff
210       //TString hname=Form("histPtRecTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());UseMCDCACorrection=kFALSE; // Pure MC efficiency for Prim. BE CAREFUL WITH MUONS!!!
211       Printf("Getting %s",hname.Data());
212       CorrFact[index]=(TH1F*)((TH1F*) hman_mc->GetPtHistogram1D(hname.Data(),-1,-1))->Clone();
213       CorrFact[index]->SetName(Form("CorrFact_%s%s",Particle[ipart].Data(),Sign[icharge].Data()));
214       CorrFact[index]->SetTitle(Form("CorrFact %s%s",Particle[ipart].Data(),Sign[icharge].Data()));
215       CorrFact[index]->SetMarkerStyle(Marker[index]);
216       CorrFact[index]->SetMarkerColor(Color[ipart]);
217       CorrFact[index]->SetLineColor(Color[ipart]);
218       hname=Form("histPtGenTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());
219       Printf("... and divide it by %s",hname.Data());
220       CorrFact[index]->Divide(CorrFact[index],(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D(hname.Data(),1,1))->Clone(),1,1,"B");//binomial error
221       if(index==0)CorrFact[index]->DrawClone();
222       else CorrFact[index]->DrawClone("same");
223     }
224   } 
225   TFile *fESD=new TFile("EffAlex/pionEffPbPb.root");
226   TH1F *hEffESD=(TH1F*)fESD->Get("effMapPionTpcOnlyNeg0");
227   hEffESD->DrawClone("same");
228   gPad->BuildLegend();
229   
230     //Normalization
231   printf("\n\n-> Spectra Normalization");
232   AliSpectraAODEventCuts* ecuts_data = (AliSpectraAODEventCuts*) _data->Get("Event Cuts");
233   Double_t events_data =  ecuts_data->NumberOfEvents();
234   Printf(": accepted events: %.0f",events_data);
235   
236   //divide RAW for Correction Factor
237   Printf("\n\n-> Using MC correction factor to correct RAW spectra");
238   TH1F *Spectra[6];
239   for(Int_t icharge=0;icharge<2;icharge++){
240     for(Int_t ipart=0;ipart<3;ipart++){
241       Int_t index=ipart+3*icharge;
242       TString hname=Form("histPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data());
243       printf("Getting %s",hname.Data());
244       Spectra[index] =(TH1F*)((TH1F*) hman_data->GetPtHistogram1D(hname.Data(),-1,-1))->Clone();
245       Spectra[index]->SetName(Form("Spectra_%s%s",Particle[ipart].Data(),Sign[icharge].Data()));
246       Spectra[index]->SetTitle(Form("Spectra %s%s",Particle[ipart].Data(),Sign[icharge].Data()));
247       Spectra[index]->SetMarkerStyle(Marker[index]);
248       Spectra[index]->SetMarkerColor(Color[ipart]);
249       Spectra[index]->SetLineColor(Color[ipart]);
250       Printf("... and divide it by %s",hname.Data());
251       Spectra[index]->Divide(CorrFact[index]);//////////////////////////////////////////////////////////////////////////////////////////FIXME
252       // if(index!=3)Spectra[index]->Divide(CorrFact[index]);
253       // else{
254       //        Spectra[index]=AliPWGHistoTools::MyDivideHistosDifferentBins(Spectra[index],hEffESD);
255       // }
256       Spectra[index]->Scale(1./events_data,"width");//NORMALIZATION
257     }
258   } 
259   
260   
261   //Geant/Fluka Correction
262   Printf("\nGF correction for Kaons");
263   TString fnameGeanFlukaK="GFCorrection/correctionForCrossSection.321.root";
264   TFile *fGeanFlukaK= new TFile(fnameGeanFlukaK.Data());
265   TH1F *hGeantFlukaKPos=(TH1F*)fGeanFlukaK->Get("gHistCorrectionForCrossSectionParticles");
266   TH1F *hGeantFlukaKNeg=(TH1F*)fGeanFlukaK->Get("gHistCorrectionForCrossSectionAntiParticles");
267   for(Int_t binK=0;binK<=Spectra[1]->GetNbinsX();binK++){
268     Float_t FlukaCorrKPos=hGeantFlukaKPos->GetBinContent(hGeantFlukaKPos->FindBin(Spectra[1]->GetBinCenter(binK)));
269     Float_t FlukaCorrKNeg=hGeantFlukaKNeg->GetBinContent(hGeantFlukaKNeg->FindBin(Spectra[4]->GetBinCenter(binK)));
270     Spectra[1]->SetBinContent(binK,Spectra[1]->GetBinContent(binK)*FlukaCorrKPos);
271     Spectra[4]->SetBinContent(binK,Spectra[4]->GetBinContent(binK)*FlukaCorrKNeg);
272     Spectra[1]->SetBinError(binK,Spectra[1]->GetBinError(binK)*FlukaCorrKPos);
273     Spectra[4]->SetBinError(binK,Spectra[4]->GetBinError(binK)*FlukaCorrKNeg);
274   }
275   
276   //Geant Fluka for P
277   Printf("\nGF correction for Protons");
278   const Int_t kNCharge=2;
279   Int_t kPos=0;
280   Int_t kNeg=1;
281   TFile* fGFProtons = new TFile ("GFCorrection/correctionForCrossSection.root");
282   TH2D * hCorrFluka[kNCharge];
283   TH2D * hCorrFluka[2];
284   hCorrFluka[kPos] = (TH2D*)fGFProtons->Get("gHistCorrectionForCrossSectionProtons");
285   hCorrFluka[kNeg] = (TH2D*)fGFProtons->Get("gHistCorrectionForCrossSectionAntiProtons");
286   for(Int_t icharge = 0; icharge < kNCharge; icharge++){
287     Int_t nbins = Spectra[2]->GetNbinsX();
288     Int_t nbinsy=hCorrFluka[icharge]->GetNbinsY();
289     for(Int_t ibin = 0; ibin < nbins; ibin++){
290       Float_t pt = Spectra[2]->GetBinCenter(ibin);
291       Float_t minPtCorrection = hCorrFluka[icharge]->GetYaxis()->GetBinLowEdge(1);
292       Float_t maxPtCorrection = hCorrFluka[icharge]->GetYaxis()->GetBinLowEdge(nbinsy+1);
293       if (pt < minPtCorrection) pt = minPtCorrection+0.0001;
294       if (pt > maxPtCorrection) pt = maxPtCorrection;
295       Float_t correction = hCorrFluka[icharge]->GetBinContent(1,hCorrFluka[icharge]->GetYaxis()->FindBin(pt));
296       //cout<<correction<<"     charge "<<icharge<<endl;
297       if(icharge==0){
298         if (correction != 0) {// If the bin is empty this is a  0
299           Spectra[2]->SetBinContent(ibin,Spectra[2]->GetBinContent(ibin)*correction);
300           Spectra[2]->SetBinError(ibin,Spectra[2]->GetBinError  (ibin)*correction);
301         }else if (Spectra[2]->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
302           cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for protons secondaries," << endl;
303           cout << " Bin content: " << Spectra[2]->GetBinContent(ibin)  << endl;
304         }
305       }
306       if(icharge==1){
307         if (correction != 0) {// If the bin is empty this is a  0
308           Spectra[5]->SetBinContent(ibin,Spectra[5]->GetBinContent(ibin)*correction);
309           Spectra[5]->SetBinError(ibin,Spectra[5]->GetBinError  (ibin)*correction);
310         }else if (Spectra[5]->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
311           cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for Antiprotons secondaries," << endl;
312           cout << " Bin content: " << Spectra[5]->GetBinContent(ibin)  << endl;
313         }
314       }
315     }   
316   }
317   
318   //DCA Correction
319   printf("\n\n-> DCA Correction");
320   DCACorrection(Spectra,hman_data,hman_mc,UseMCDCACorrection);
321   
322   
323   TH1F *MCTruth[6];
324   for(Int_t icharge=0;icharge<2;icharge++){
325     for(Int_t ipart=0;ipart<3;ipart++){
326       Int_t index=ipart+3*icharge;
327       hname=Form("histPtGenTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());
328       MCTruth[index]=(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D(hname.Data(),0.9,1.1))->Clone();
329       MCTruth[index]->Scale(1./events_data,"width");//NORMALIZATION
330     }
331   }
332   
333   //Contamination
334   Printf("\n\n-> Contamination from MC");
335   TH1F *Cont[6];
336   TCanvas *ccont=new TCanvas("ccont","ccont",700,500);
337   for(Int_t icharge=0;icharge<2;icharge++){
338     for(Int_t ipart=0;ipart<3;ipart++){
339       Int_t index=ipart+3*icharge;
340       TString hname=Form("histPtRecTrue%s%s",Particle[ipart].Data(),Sign[icharge].Data());
341       Printf("Getting %s",hname.Data());
342       Cont[index] =(TH1F*)((TH1F*) hman_mc->GetPtHistogram1D(hname.Data(),-1,-1))->Clone();
343       Cont[index]->SetName(Form("Cont_%s%s",Particle[ipart].Data(),Sign[icharge].Data()));
344       Cont[index]->SetTitle(Form("RecTrue/RecSigma %s%s",Particle[ipart].Data(),Sign[icharge].Data()));
345       Cont[index]->SetMarkerStyle(Marker[index]);
346       Cont[index]->SetMarkerColor(Color[ipart]);
347       Cont[index]->SetLineColor(Color[ipart]);
348       hname=Form("histPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data());
349       Printf("... and divide it by %s",hname.Data());
350       Cont[index]->Divide((TH1F*) hman_mc->GetPtHistogram1D(hname.Data(),-1,-1));
351       if(index==0)Cont[index]->DrawClone();
352       else Cont[index]->DrawClone("same");
353       //Spectra[index]->Multiply(Cont[index]);
354     }
355   } 
356   gPad->BuildLegend();
357   
358   
359   //Drawing Final Spectra
360   TCanvas *cspectra=new TCanvas("cspectra","cspectra",700,500);
361   for(Int_t icharge=0;icharge<2;icharge++){
362     for(Int_t ipart=0;ipart<3;ipart++){
363       Int_t index=ipart+3*icharge;
364       if(index==0)Spectra[index]->DrawClone();
365       else Spectra[index]->DrawClone("same");
366       //MCTruth[index]->DrawClone("same");
367     }
368   } 
369   gPad->BuildLegend();
370   
371   //saving spectra
372   Printf("\n\n-> Saving spectra in Out%s_%s.root",sname.Data(),fold.Data());
373   TFile * fout=new TFile(Form("results/Out_%s_%s.root",sname.Data(),fold.Data()),"RECREATE");
374   for(Int_t icharge=0;icharge<2;icharge++){
375     for(Int_t ipart=0;ipart<3;ipart++){
376       Int_t index=ipart+3*icharge;
377       Spectra[index]->Write();
378       CorrFact[index]->Write();
379       MCTruth[index]->Write();
380     }
381   }
382   
383   
384   //if Bin 0-5% with no cut ratio with combined analysis
385   if(ibinToCompare!=-1){
386     TCanvas *CratioComb=new TCanvas("CratioComb","CratioComb",700,500);
387     CratioComb->Divide(3,2);
388     TString nameComb[6]={Form("cent%d_pion_plus",ibinToCompare),Form("cent%d_kaon_plus",ibinToCompare),Form("cent%d_proton_plus",ibinToCompare),
389                          Form("cent%d_pion_minus",ibinToCompare),Form("cent%d_kaon_minus",ibinToCompare),Form("cent%d_proton_minus",ibinToCompare)};
390     TFile *fComb=new TFile("Combined05/SPECTRA_COMB_20120323.root");
391     for(Int_t icharge=0;icharge<2;icharge++){
392       for(Int_t ipart=0;ipart<3;ipart++){
393         Int_t index=ipart+3*icharge;
394         Spectra[index]->GetXaxis()->SetRangeUser(0,2.5);
395         TH1F *hcomb=fComb->Get(nameComb[index].Data())->Clone();
396         CratioComb->cd(ipart+1);
397         gPad->SetGridy();
398         gPad->SetGridx();
399         if(icharge==0)Spectra[index]->DrawClone();
400         else Spectra[index]->DrawClone("same");
401         hcomb->DrawClone("same");
402         Spectra[index]->Divide(hcomb);
403         Spectra[index]->SetMaximum(1.3);
404         Spectra[index]->SetMinimum(0.7);
405         CratioComb->cd(ipart+4);
406         gPad->SetGridy();
407         gPad->SetGridx();
408         if(icharge==0)Spectra[index]->DrawClone();
409         else Spectra[index]->DrawClone("same");
410       }
411     }
412   }     
413   
414
415   
416   //comparison with charged hadron
417   Printf("\n\n-> ChargedHadron comparison");
418   TCanvas *cAllCh=new TCanvas("cAllCh","cAllCh",700,500);
419   cAllCh->Divide(1,4);
420   TH1F *hChHad_data=(TH1F*)((TH1F*)hman_data->GetPtHistogram1D("histPtRec",-1,-1))->Clone();
421   //fraction of sec in MC
422   TH1F *hPrimRec_mc=(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D("histPtRecPrimary",-1,-1))->Clone();
423   TH1F *hAllRec_mc=(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D("histPtRec",-1,-1))->Clone();
424   for(Int_t ibin=1;ibin<=hChHad_data->GetNbinsX();ibin++){
425     Double_t en_data=hChHad_data->GetBinContent(ibin);
426     Double_t en_mc=hAllRec_mc->GetBinContent(ibin);
427     Double_t prim_mc=hPrimRec_mc->GetBinContent(ibin);
428     if(en_mc!=0)hChHad_data->SetBinContent(ibin,en_data-(en_data*(en_mc-prim_mc)*1.2/en_mc));
429     Printf("Before: %.0f After: %.0f  fraction: %.1f",en_data,hChHad_data->GetBinContent(ibin),hChHad_data->GetBinContent(ibin)/en_data);
430   }
431   cAllCh->cd(1);
432   gPad->SetGridy();
433   gPad->SetGridx();
434   hPrimRec_mc->Divide(hAllRec_mc);
435   hPrimRec_mc->DrawClone();
436   //efficiency for primaries
437   TH1F *hEff_mc=(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D("histPtRecPrimary",-1,-1))->Clone();
438   hEff_mc->Divide((TH1F*)((TH1F*)hman_mc->GetPtHistogram1D("histPtGen",1,1))->Clone());
439   Printf("+-------------------------------------%f ",((TH1F*)hman_mc->GetPtHistogram1D("histPtGen",1,1))->GetEntries()/1.6/ecuts_mc->NumberOfEvents());
440   hChHad_data->Scale(1./events_data,"width");//NORMALIZATION
441   hChHad_data->Divide(hEff_mc);//Efficiency
442   hChHad_data->Scale(1./(2*tcuts_data->GetEta()));
443   cAllCh->cd(2);
444   gPad->SetGridy();
445   gPad->SetGridx();
446   hEff_mc->DrawClone();
447   cAllCh->cd(3);
448   gPad->SetGridy();
449   gPad->SetGridx();
450   hChHad_data->DrawClone();
451   TFile *fCh=new TFile("ChargedHadron/SPECTRA_UNID_110906.root");
452   TH1F *hCh=fCh->Get("hpt_c0_5");
453   //invariant yield
454   for(Int_t ibin=0;ibin<hCh->GetNbinsX();ibin++){
455     hCh->SetBinContent(ibin,hCh->GetBinContent(ibin)*(2*TMath::Pi()*hCh->GetBinCenter(ibin)));
456     hCh->SetBinError(ibin,hCh->GetBinError(ibin)*(2*TMath::Pi()*hCh->GetBinCenter(ibin)));
457   }
458   hCh->DrawClone("same");
459   TH1F *gRatio=AliPWGHistoTools::MyDivideHistosDifferentBins(hChHad_data,hCh);
460   gRatio->SetMaximum(1.5);
461   gRatio->SetMinimum(.5);
462   cAllCh->cd(4);
463   gPad->SetGridy();
464   gPad->SetGridx();
465   gRatio->DrawClone("");
466   
467   // //Comparison of efficiency with TPCTOF ESD analysis
468   Printf("\n\n-> Calculating Efficiency to be compared with ESD analysis");
469   TH1F *EffTRUEPions;
470   TH1F *EffSIGMAPions;
471   TCanvas *cEffESD=new TCanvas("cEffESD","cEffESD",700,500);
472   cEffESD->Divide(1,2);
473   Int_t icharge=1;
474   Int_t ipart=0;
475   //using MC truth
476   TString hname=Form("histPtRecTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());
477   Printf("Getting %s",hname.Data());
478   EffTRUEPions=(TH1F*)((TH1F*) hman_mc->GetPtHistogram1D(hname.Data(),-1,-1))->Clone();
479   EffTRUEPions->SetName(Form("Eff TRUE_%s%s",Particle[ipart].Data(),Sign[icharge].Data()));
480   EffTRUEPions->SetTitle(Form("Eff TRUE %s%s",Particle[ipart].Data(),Sign[icharge].Data()));
481   EffTRUEPions->SetMarkerStyle(Marker[icharge]);
482   EffTRUEPions->SetMarkerColor(Color[ipart]);
483   EffTRUEPions->SetLineColor(Color[ipart]);
484   hname=Form("histPtGenTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());
485   Printf("... and divide it by %s",hname.Data());
486   EffTRUEPions->Divide(EffTRUEPions,(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D(hname.Data(),1,1))->Clone(),1,1,"B");//binomial error
487   //using NSigma
488   hname=Form("histPtRecSigmaPrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());
489   Printf("Getting %s",hname.Data());
490   EffSIGMAPions=(TH1F*)((TH1F*) hman_mc->GetPtHistogram1D(hname.Data(),-1,-1))->Clone();
491   EffSIGMAPions->SetName(Form("Eff SIGMA_%s%s",Particle[ipart].Data(),Sign[icharge].Data()));
492   EffSIGMAPions->SetTitle(Form("Eff SIGMA %s%s",Particle[ipart].Data(),Sign[icharge].Data()));
493   EffSIGMAPions->SetMarkerStyle(Marker[icharge]);
494   EffSIGMAPions->SetMarkerColor(Color[ipart+1]);
495   EffSIGMAPions->SetLineColor(Color[ipart+1]);
496   hname=Form("histPtGenTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());
497   Printf("... and divide it by %s",hname.Data());
498   EffSIGMAPions->Divide(EffSIGMAPions,(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D(hname.Data(),1,1))->Clone(),1,1,"B");//binomial error
499   cEffESD->cd(1);
500   EffTRUEPions->DrawClone("lhist");
501   EffSIGMAPions->DrawClone("lhistsame");
502   hEffESD->DrawClone("lhistsame");
503   gPad->BuildLegend();
504   cEffESD->cd(2);
505   TH1F *hRatioTRUE=AliPWGHistoTools::MyDivideHistosDifferentBins(EffTRUEPions,hEffESD);
506   hRatioTRUE->DrawClone("lhist");
507   TH1F *hRatioSIGMA=AliPWGHistoTools::MyDivideHistosDifferentBins(EffSIGMAPions,hEffESD);
508   hRatioSIGMA->DrawClone("lhistsame");
509   
510   
511   // //Muon over Pion Ratio
512   // TCanvas *cMu=new TCanvas("cMu","cMu");
513   // TH1F *hMuOverPi[2];
514   // for(Int_t icharge=0;icharge<2;icharge++){
515   //   TString hname=Form("histPtRecTruePrimaryMuon%s",Sign[icharge].Data());
516   //   hMuOverPi[icharge]=(TH1F*)((TH1F*) hman_mc->GetPtHistogram1D(hname.Data(),-1,-1))->Clone();
517   //   hname=Form("histPtRecTruePrimaryPion%s",Sign[icharge].Data());
518   //   hMuOverPi[icharge]->Divide((TH1F*)((TH1F*) hman_mc->GetPtHistogram1D(hname.Data(),-1,-1))->Clone());
519   //   if(icharge==0)hMuOverPi[icharge]->DrawClone();
520   //   else hMuOverPi[icharge]->DrawClone("same");
521   // }
522 }
523
524
525
526 void DCACorrection(TH1F **Spectra, AliSpectraAODHistoManager* hman_data, AliSpectraAODHistoManager* hman_mc,Bool_t UseMCDCACorrection){
527   
528   Double_t FitRange[2]={-3,3};
529   Int_t nrebin=20;
530   Printf("\DCACorr");
531   TH1F *hcorrection[2];
532   TCanvas *ccorrection=new TCanvas("DCAcorrection","DCAcorrection",700,500);
533   ccorrection->Divide(2,1);
534   TString sample[2]={"data","mc"};
535   for(Int_t icharge=0;icharge<2;icharge++){
536     for(Int_t ipart=0;ipart<3;ipart++){
537       Int_t index=ipart+3*icharge;
538       for(Int_t isample=0;isample<2;isample++){
539         TCanvas *cDCA=new TCanvas(Form("cDCA%d%s",index,sample[isample].Data()),Form("cDCA%d%s",index,sample[isample].Data()),700,500);
540         hcorrection[isample]=(TH1F*)Spectra[index]->Clone();
541         hcorrection[isample]->Reset("all");
542         cDCA->Divide(8,4);
543         for(Int_t ibin_data=6;ibin_data<38;ibin_data++){
544           Double_t lowedge=Spectra[index]->GetBinLowEdge(ibin_data);
545           Double_t binwidth=Spectra[index]->GetBinWidth(ibin_data);
546           if(isample==0)TH1F *hToFit =(TH1F*) ((TH1F*)hman_data->GetDCAHistogram1D(Form("histPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge+binwidth))->Clone();
547           if(isample==1)TH1F *hToFit =(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("histPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge+binwidth))->Clone();
548           TH1F *hmc1=(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("histPtRecSigmaPrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge+binwidth))->Clone();
549           TH1F *hmc2=(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("histPtRecSigmaSecondaryWeakDecay%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge+binwidth))->Clone();
550           TH1F *hmc3=(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("histPtRecSigmaSecondaryMaterial%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge+binwidth))->Clone();
551           Double_t minentries=5;
552           if(hToFit->GetEntries()<=minentries || hmc1->GetEntries()<=minentries || hmc2->GetEntries()<=minentries || hmc3->GetEntries()<=minentries)continue;
553           hmc1->Rebin(nrebin);
554           hmc2->Rebin(nrebin);
555           hmc3->Rebin(nrebin);
556           hToFit->Rebin(nrebin);
557           //Data and MC can have different stat
558           hToFit->Sumw2();
559           hmc1->Sumw2();
560           hmc2->Sumw2();
561           hmc3->Sumw2();
562           hToFit->Scale(1./hToFit->GetEntries());
563           Double_t normMC=hmc1->GetEntries()+hmc2->GetEntries()+hmc3->GetEntries();
564           hmc1->Scale(1./normMC);
565           hmc2->Scale(1./normMC);
566           hmc3->Scale(1./normMC);
567           cDCA->cd(ibin_data-5);
568           gPad->SetGridy();
569           gPad->SetGridx();
570           gPad->SetLogy();
571           hToFit->DrawClone("lhist");
572           hmc3->DrawClone("lhist");
573           TObjArray *mc = new TObjArray(3);        // MC histograms are put in this array
574           mc->Add(hmc1);
575           mc->Add(hmc2);
576           mc->Add(hmc3);
577           TFractionFitter* fit = new TFractionFitter(hToFit,mc); // initialise
578           fit->Constrain(0,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
579           fit->Constrain(1,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
580           fit->Constrain(2,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
581           fit->SetRangeX(hToFit->GetXaxis()->FindBin(FitRange[0]),hToFit->GetXaxis()->FindBin(FitRange[1]));
582           hToFit->GetXaxis()->SetRange(hToFit->GetXaxis()->FindBin(FitRange[0]),hToFit->GetXaxis()->FindBin(FitRange[1]));
583           hToFit->SetTitle(Form("DCA distr - %s",sample[isample].Data()));
584           Int_t status = fit->Fit();               // perform the fit
585           cout << "fit status: " << status << endl;
586           if (status == 0) {                       // check on fit status
587             TH1F* result = (TH1F*) fit->GetPlot();
588             TH1F* PrimMCPred=(TH1F*)fit->GetMCPrediction(0);
589             TH1F* secStMCPred=(TH1F*)fit->GetMCPrediction(1);
590             TH1F* secMCPred=(TH1F*)fit->GetMCPrediction(2);
591             //Drawing section
592             PrimMCPred->SetLineColor(2);
593             secStMCPred->SetLineColor(6);
594             secMCPred->SetLineColor(4);
595             hToFit->DrawClone("lhist");
596             result->SetTitle("Fit result");
597             result->DrawClone("lhistsame");
598             PrimMCPred->DrawClone("lhistsame");
599             secStMCPred->DrawClone("lhistsame");
600             secMCPred->DrawClone("lhistsame");
601             Double_t v1=0,v2=0,v3=0;
602             Double_t ev1=0,ev2=0,ev3=0;
603             fit->GetResult(0,v1,ev1);
604             fit->GetResult(1,v2,ev2);
605             fit->GetResult(2,v3,ev3);
606             Printf("\n\n\n\n\nv1:%f  v2:%f  v3:%f   ev1:%f  ev2:%f  ev3:%f   ",v1,v2,v3,ev1,ev2,ev3);
607             hcorrection[isample]->SetBinContent(ibin_data,v1/(v1+v2+v3));
608             //hcorrection[isample]->SetBinError(ibin_data,TMath::Sqrt(ev1*ev1+ev2*ev2+ev3*ev3));
609             hcorrection[isample]->SetBinError(ibin_data,0);
610           }
611           else{
612             hcorrection[isample]->SetBinContent(ibin_data,1);
613             hcorrection[isample]->SetBinError(ibin_data,0);
614           }
615           Printf("deleting");
616           delete hToFit;
617         }
618         
619         ccorrection->cd(icharge+1);
620         gPad->SetGridy();
621         gPad->SetGridx();
622         hcorrection[isample]->SetTitle(Form("DCA corr %s %s %s",Particle[ipart].Data(),Charge[icharge].Data(),sample[isample].Data()));
623         hcorrection[isample]->SetLineWidth(2);
624         hcorrection[isample]->SetLineColor(Color[ipart]);
625         hcorrection[isample]->SetLineStyle(isample+1);
626         hcorrection[isample]->SetMarkerColor(Color[ipart]);
627         hcorrection[isample]->GetXaxis()->SetRangeUser(0.2,2.5);
628         if(ipart==0 && isample==0)hcorrection[isample]->DrawClone("lhist");
629         else hcorrection[isample]->DrawClone("lhistsame");
630         //else hcorrection[isample]->DrawClone("psame");
631         // smooth the DCA correction
632         // TF1 *fitfun = new TF1("fitfun","[0]+[1]/x^[2]",0.2,0.8);
633         // hcorrection[isample]->Fit(fitfun,"WRVMN","",0.35,0.8);
634         // fitfun->SetLineWidth(1.5);
635         // fitfun->SetLineColor(ipart+1);
636         // fitfun->SetLineStyle(ipart+1);
637         // fitfun->DrawClone("same");
638         // for(Int_t ibin=1;ibin<30;ibin++){
639         //      hcorrection[isample]->SetBinContent(ibin,fitfun->Eval(hcorrection[isample]->GetBinCenter(ibin)));
640         // }
641       }
642       Spectra[index]->Multiply(hcorrection[0]);//multiply for data
643       Printf("DCACorrection for DATA used: Spectra[index]->Multiply(hcorrection[0])")
644         if(UseMCDCACorrection){
645           Spectra[index]->Divide(hcorrection[1]); //divide by Monte Carlo
646           Printf("DCACorrection for MC used: Spectra[index]->Divide(hcorrection[1]")
647             }
648     }
649   }
650   ccorrection->cd(1);
651   gPad->BuildLegend();
652   ccorrection->cd(2);
653   gPad->BuildLegend();
654 }
655
656