Adding some plots to the QAPlotsBoth.C macro.
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / TestAOD / QAPlotsBoth.C
1 Float_t QAPlotsBoth( AliSpectraBothHistoManager* hman_data, AliSpectraBothHistoManager* hman_mc,
2               AliSpectraBothEventCuts* ecuts_data, AliSpectraBothEventCuts* ecuts_mc,
3               AliSpectraBothTrackCuts* tcuts_data, AliSpectraBothTrackCuts* tcuts_mc,
4               TList * flistqa,TList * flistcanvas)
5 {
6 TString pidmethods[3]={"TPC","TOF","TPCTOF"};   
7         Double_t neventsdata =  ecutsdata->NumberOfPhysSelEvents();
8         Double_t neventsmc =  ecutsmc->NumberOfPhysSelEvents();
9         
10         
11         
12         for(Int_t ipart=0;ipart<3;ipart++)
13         {
14                         
15                         for(Int_t imethod=0;imethod<3;imethod++)
16                         {
17                                  TH2F *nsig_data = (TH2F*)((TH2F*)hman_data->GetNSigHistogram(Form("hHistNSig%sPt%s",Particle[ipart].Data(),pidmethods[imethod].Data())))->Clone();
18                                 // nsig_data->RebinX(20);                        
19                                 // nsig_data->RebinY(4);
20                                 // nsig_data->Sumw2();
21
22                                  TH2F *nsig_mc = (TH2F*)((TH2F*)hman_mc->GetNSigHistogram(Form("hHistNSig%sPt%s",Particle[ipart].Data(),pidmethods[imethod].Data())))->Clone();
23                                  //nsig_mc->RebinX(20);                  
24                                 // nsig_mc->RebinY(4);
25                                 // nsig_mc->Sumw2();
26                                  
27                                 Int_t ibin=1;
28                                 Float_t binsize=nsig_mc->GetXaxis()->GetBinWidth(1);
29
30                                 TH1F* maxposdata=(TH1F*)nsig_data->ProjectionX(Form("%s%sdatamaxpos",Particle[ipart].Data(),pidmethods[imethod].Data()),-1,-1));
31                                 maxposdata->Reset();
32                                 maxposdata->SetTitle(";p_{T} (GeV/c);max in (-2,2)");
33                                 TH1F* maxposmc=(TH1F*)nsig_data->ProjectionX(Form("%s%smcmaxpos",Particle[ipart].Data(),pidmethods[imethod].Data()),-1,-1));
34                                 maxposmc->Reset();
35                                 maxposmc->SetTitle(";p_{T} (GeV/c);max in (-2,2)");
36
37         
38                                  while (ibin*binsize<3.0)
39                                  {
40                                         // TCanvas* c=new TCanvas(Form("canvas%s%s%d",Particle[ipart].Data(),pidmethods[imethod].Data(),ibin),Form("canvas%s%s%d",Particle[ipart].Data(),pidmethods[imethod].Data(),ibin),700,500);
41                                         
42                                          TH1F *nsig_data_Proj1=(TH1F*)nsig_data->ProjectionY(Form("%s%sdata[%.2f,%.2f]",Particle[ipart].Data(),pidmethods[imethod].Data(),nsig_data->GetXaxis()->GetBinLowEdge(ibin),nsig_data->GetXaxis()->GetBinUpEdge(ibin)),ibin,ibin));
43                                          TH1F *nsig_mc_Proj1=(TH1F*)nsig_mc->ProjectionY(Form("%s%smc[%.2f,%.2f]",Particle[ipart].Data(),pidmethods[imethod].Data(),nsig_mc->GetXaxis()->GetBinLowEdge(ibin),nsig_mc->GetXaxis()->GetBinUpEdge(ibin)),ibin,ibin));
44                                          nsig_data_Proj1->GetXaxis()->SetRangeUser(-3,3);
45                                          nsig_data_Proj1->SetLineColor(kRed);
46                                          if(nsig_data_Proj1->Integral()<1&&nsig_mc_Proj1->Integral()<1)
47                                          {      
48                                                 ibin++; 
49                                                 continue;
50                                          } 
51                                                                         //      nsig_data_Proj1->Sumw2();
52                                 //      nsig_mc_Proj1->Sumw2();
53                                           //nsig_data_Proj1->GetXaxis()->SetRangeUser(-5,5);
54                                          
55                                          //c->cd()->SetLogy();
56                                          //nsig_data_Proj1->Draw();
57                                          //nsig_mc_Proj1->Draw("same");
58                                         nsig_data_Proj1->GetXaxis()->SetRange(nsig_data_Proj1->GetXaxis()->FindBin(-2.0),nsig_data_Proj1->GetXaxis()->FindBin(2.0));
59                                         if(nsig_data_Proj1->GetMaximumBin()<=nsig_data_Proj1->GetXaxis()->FindBin(2.0)&&nsig_data_Proj1->GetMaximumBin()>=nsig_data_Proj1->GetXaxis()->FindBin(-2.0))
60                                         {
61                                                 maxposdata->SetBinContent(ibin,nsig_data_Proj1->GetXaxis()->GetBinCenter(nsig_data_Proj1->GetMaximumBin()));    
62                                                 maxposdata->SetBinError(ibin,nsig_data_Proj1->GetXaxis()->GetBinWidth(nsig_data_Proj1->GetMaximumBin())/2.0);   
63                                         }
64                                 cout<<Form("%s%sdatamaxpos",Particle[ipart].Data(),pidmethods[imethod].Data())<<" "<<nsig_data_Proj1->GetMaximumBin()<<" "<<nsig_data_Proj1->GetXaxis()->FindBin(2.0)<<" "<<nsig_data_Proj1->GetXaxis()->FindBin(-2.0)<<" "<<nsig_data_Proj1->GetXaxis()->GetBinCenter(nsig_data_Proj1->GetMaximumBin())<<" "<<ibin<<endl;
65
66                                         nsig_data_Proj1->GetXaxis()->SetRange(0,nsig_data_Proj1->GetXaxis()->GetNbins());
67                                                                         
68                                         nsig_mc_Proj1->GetXaxis()->SetRange(nsig_mc_Proj1->GetXaxis()->FindBin(-2.0),nsig_mc_Proj1->GetXaxis()->FindBin(2.0));
69                                         if(nsig_mc_Proj1->GetMaximumBin()<=nsig_mc_Proj1->GetXaxis()->FindBin(2.0)&&nsig_mc_Proj1->GetMaximumBin()>=nsig_mc_Proj1->GetXaxis()->FindBin(-2.0))
70                                         {
71                                                 maxposmc->SetBinContent(ibin,nsig_mc_Proj1->GetXaxis()->GetBinCenter(nsig_mc_Proj1->GetMaximumBin()));  
72                                                 maxposmc->SetBinError(ibin,nsig_mc_Proj1->GetXaxis()->GetBinWidth(nsig_mc_Proj1->GetMaximumBin())/2.0); 
73                                         }
74                                                 cout<<Form("%s%smcmaxpos",Particle[ipart].Data(),pidmethods[imethod].Data())<<" "<<nsig_mc_Proj1->GetMaximumBin()<<" "<<nsig_mc_Proj1->GetXaxis()->FindBin(2.0)<<" "<<nsig_mc_Proj1->GetXaxis()->FindBin(-2.0)<<" "<<ibin<<" "<<nsig_mc_Proj1->GetXaxis()->GetBinCenter(nsig_mc_Proj1->GetMaximumBin())<<endl;
75
76                                         nsig_mc_Proj1->GetXaxis()->SetRange(0,nsig_mc_Proj1->GetXaxis()->GetNbins());
77
78
79                                         nsig_data_Proj1->Scale(1.0/neventsdata);
80                                          nsig_mc_Proj1->Scale(1.0/neventsmc);
81
82                                         
83                                         flistqa->Add(nsig_data_Proj1);
84                                         flistqa->Add(nsig_mc_Proj1);
85                                         ibin++;
86                                  }
87                                 flistqa->Add(maxposmc);
88                                 flistqa->Add(maxposdata);
89                         }
90         }
91         TH1F* fHistoVtxAftSeldata=(TH1F*)ecuts_data->GetHistoVtxAftSel();
92         TH1F* fHistoVtxAftSelmc=(TH1F*)ecuts_mc->GetHistoVtxAftSel();
93         flistcanvas->Add(plot_on_canvas("vertex",fHistoVtxAftSeldata,fHistoVtxAftSelmc));
94         TF1* fdata=new TF1("dataveretxfit","gausn");
95         TF1* fmc=new TF1("mcveretxfit","gausn");
96         fHistoVtxAftSeldata->Fit("dataveretxfit","0");
97         fHistoVtxAftSelmc->Fit("mcveretxfit","0");
98         Float_t datavertexratio=fHistoVtxAftSeldata->Integral(-1,-1,"width")/fdata->GetParameter(0);
99         Float_t mcvertexratio=fHistoVtxAftSelmc->Integral(-1,-1,"width")/fmc->GetParameter(0);
100         
101
102          TH1F* fHistoEtaAftSeldata=(TH1F*)ecuts_data->GetHistoEtaAftSel();
103          TH1F* fHistoEtaAftSelmc=(TH1F*)ecuts_mc->GetHistoEtaAftSel();
104         flistcanvas->Add(plot_on_canvas("ETA",fHistoEtaAftSeldata,fHistoEtaAftSelmc));
105
106
107          TH1F* fITSclustershistdata=(TH1F*)tcuts_data->GetHistoNclustersITS();
108           TH1F* fITSclustershistmc=(TH1F*)tcuts_mc->GetHistoNclustersITS();
109
110         flistcanvas->Add(plot_on_canvas("NITS",fITSclustershistdata,fITSclustershistmc));
111         cout<<" data "<<datavertexratio<<" mc "<<mcvertexratio<<endl;
112         
113         TH2F* hmul=(TH2F*)hman_mc->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul");   
114         hmul->Sumw2();  
115         TCanvas* cbc=new TCanvas("broken chunks","broken chunks",1200,600);
116         cbc->Divide(2,1);
117         cbc->cd(1);
118         hmul->Draw();
119         cbc->cd(2);
120         TH1F* nonzero=(TH1F*)hmul->ProjectionX("nonzerotracks",2,-1);
121         nonzero->SetMarkerColor(kRed);
122         nonzero->SetMarkerStyle(21);
123         TH1F* binzero=(TH1F*)hmul->ProjectionX("binzerotracks",1,1);
124         binzero->SetMarkerColor(kBlack);
125         binzero->SetMarkerStyle(22);
126         binzero->Sumw2();
127         nonzero->Sumw2();
128         binzero->Divide(nonzero);
129         TF1* badchunk=new TF1("badchunkfit","pol0",10,40);
130         binzero->Fit("badchunkfit","R");
131         Float_t badchunksfraction=badchunk->GetParameter(0);
132         binzero->Draw("E1");
133         flistcanvas->Add(cbc);
134         
135         return (1.0-badchunksfraction)*mcvertexratio/datavertexratio;
136
137 }
138 TCanvas* plot_on_canvas(TString name, TH1* h1,TH1* h2)
139 {
140         TCanvas* cvrt=new TCanvas(name.Data(),name.Data(),600,600);
141         cvrt->cd();
142         h1->SetLineColor(kRed);
143         h2->SetLineColor(kBlue);
144         h1->Sumw2();
145         h2->Sumw2();
146         TLegend *lvtr=new TLegend(0.2,0.2,0.5,0.3,"","NDC");
147         lvtr->SetLineColor(kWhite);
148         lvtr->AddEntry(h1,"data","l");
149         lvtr->AddEntry(h2,"MC","l");
150         h1->Scale(1.0/h1->GetBinContent(h1->GetXaxis()->FindBin(0.0)));
151         h2->Scale(1.0/h2->GetBinContent(h2->GetXaxis()->FindBin(0.0)));
152         h1->DrawCopy("L");
153         h2->DrawCopy("Lsame");
154         lvtr->Draw();
155         return cvrt;
156 }
157