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,Bool_t fullicorr=kTRUE)
   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                                 if(!hman_data->GetNSigHistogram(Form("hHistNSig%sPt%s",Particle[ipart].Data(),pidmethods[imethod].Data())))
  18                                         continue;
  19                                  TH2F *nsig_data = (TH2F*)((TH2F*)hman_data->GetNSigHistogram(Form("hHistNSig%sPt%s",Particle[ipart].Data(),pidmethods[imethod].Data())))->Clone();
  20                                 // nsig_data->RebinX(20);
  21                                 // nsig_data->RebinY(4);
  22                                 // nsig_data->Sumw2();
  23                                 if(!nsig_data)
  24                                         continue;
  25
  26                                 if(!hman_mc->GetNSigHistogram(Form("hHistNSig%sPt%s",Particle[ipart].Data(),pidmethods[imethod].Data())))
  27                                         continue;
  28                                  TH2F *nsig_mc = (TH2F*)((TH2F*)hman_mc->GetNSigHistogram(Form("hHistNSig%sPt%s",Particle[ipart].Data(),pidmethods[imethod].Data())))->Clone();
  29                                  //nsig_mc->RebinX(20);
  30                                 // nsig_mc->RebinY(4);
  31                                 // nsig_mc->Sumw2();
  32                                 if(!nsig_mc)
  33                                         continue;
  34
  35
  36                                 Int_t ibin=1;
  37                                 Float_t binsize=nsig_mc->GetXaxis()->GetBinWidth(1);
  38
  39                                 TH1F* maxposdata=(TH1F*)nsig_data->ProjectionX(Form("%s%sdatamaxpos",Particle[ipart].Data(),pidmethods[imethod].Data()),-1,-1);
  40                                 maxposdata->Reset();
  41                                 maxposdata->SetTitle(";p_{T} (GeV/c);max in (-2,2)");
  42                                 TH1F* maxposmc=(TH1F*)nsig_data->ProjectionX(Form("%s%smcmaxpos",Particle[ipart].Data(),pidmethods[imethod].Data()),-1,-1);
  43                                 maxposmc->Reset();
  44                                 maxposmc->SetTitle(";p_{T} (GeV/c);max in (-2,2)");
  45
  46
  47                                  while (ibin*binsize<3.0)
  48                                  {
  49                                         // 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);
  50
  51                                          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);
  52                                          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);
  53                                          nsig_data_Proj1->GetXaxis()->SetRangeUser(-3,3);
  54                                          nsig_data_Proj1->SetLineColor(kRed);
  55                                          if(nsig_data_Proj1->Integral()<1&&nsig_mc_Proj1->Integral()<1)
  56                                          {
  57                                                 ibin++;
  58                                                 continue;
  59                                          }
  60                                                                         //      nsig_data_Proj1->Sumw2();
  61                                 //      nsig_mc_Proj1->Sumw2();
  62                                           //nsig_data_Proj1->GetXaxis()->SetRangeUser(-5,5);
  63
  64                                          //c->cd()->SetLogy();
  65                                          //nsig_data_Proj1->Draw();
  66                                          //nsig_mc_Proj1->Draw("same");
  67                                         nsig_data_Proj1->GetXaxis()->SetRange(nsig_data_Proj1->GetXaxis()->FindBin(-2.0),nsig_data_Proj1->GetXaxis()->FindBin(2.0));
  68                                         if(nsig_data_Proj1->GetMaximumBin()<=nsig_data_Proj1->GetXaxis()->FindBin(2.0)&&nsig_data_Proj1->GetMaximumBin()>=nsig_data_Proj1->GetXaxis()->FindBin(-2.0))
  69                                         {
  70                                                 maxposdata->SetBinContent(ibin,nsig_data_Proj1->GetXaxis()->GetBinCenter(nsig_data_Proj1->GetMaximumBin()));
  71                                                 maxposdata->SetBinError(ibin,nsig_data_Proj1->GetXaxis()->GetBinWidth(nsig_data_Proj1->GetMaximumBin())/2.0);
  72                                         }
  73                                 //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;
  74
  75                                         nsig_data_Proj1->GetXaxis()->SetRange(0,nsig_data_Proj1->GetXaxis()->GetNbins());
  76
  77                                         nsig_mc_Proj1->GetXaxis()->SetRange(nsig_mc_Proj1->GetXaxis()->FindBin(-2.0),nsig_mc_Proj1->GetXaxis()->FindBin(2.0));
  78                                         if(nsig_mc_Proj1->GetMaximumBin()<=nsig_mc_Proj1->GetXaxis()->FindBin(2.0)&&nsig_mc_Proj1->GetMaximumBin()>=nsig_mc_Proj1->GetXaxis()->FindBin(-2.0))
  79                                         {
  80                                                 maxposmc->SetBinContent(ibin,nsig_mc_Proj1->GetXaxis()->GetBinCenter(nsig_mc_Proj1->GetMaximumBin()));
  81                                                 maxposmc->SetBinError(ibin,nsig_mc_Proj1->GetXaxis()->GetBinWidth(nsig_mc_Proj1->GetMaximumBin())/2.0);
  82                                         }
  83                                 //              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;
  84
  85                                         nsig_mc_Proj1->GetXaxis()->SetRange(0,nsig_mc_Proj1->GetXaxis()->GetNbins());
  86
  87                                         if(neventsdata>0.0&&neventsmc>0.0)
  88                                         {
  89                                                 nsig_data_Proj1->Sumw2();
  90                                                 nsig_mc_Proj1->Sumw2();
  91                                                 nsig_data_Proj1->Scale(1.0/neventsdata);
  92                                                 nsig_mc_Proj1->Scale(1.0/neventsmc);
  93                                         }
  94
  95                                         flistqa->Add(nsig_data_Proj1);
  96                                         flistqa->Add(nsig_mc_Proj1);
  97                                         ibin++;
  98                                  }
  99                                 flistqa->Add(maxposmc);
 100                                 flistqa->Add(maxposdata);
 101                         }
 102         }
 103         TH1F* fHistoVtxAftSeldata=(TH1F*)ecuts_data->GetHistoVtxAftSel();
 104         TH1F* fHistoVtxAftSelmc=(TH1F*)ecuts_mc->GetHistoVtxAftSel();
 105         Int_t binstartx=1;
 106          while(fHistoVtxAftSeldata->GetBinContent(binstartx)<1&&binstartx<fHistoVtxAftSeldata->GetXaxis()->GetNbins())
 107                 binstartx++;
 108         binstartx++;
 109         flistcanvas->Add(plot_on_canvas("vertex",fHistoVtxAftSeldata,fHistoVtxAftSelmc));
 110         /*
 111
 112         TF1* fdata=new TF1("dataveretxfit","gausn");
 113         TF1* fmc=new TF1("mcveretxfit","gausn");
 114         //we strat fit a second not empty bin
 115         Float_t minfit=fHistoVtxAftSeldata->GetXaxis()->GetBinCenter(binstartx);
 116         cout<<"fit starts "<<minfit<<endl;
 117         fHistoVtxAftSeldata->Fit("dataveretxfit","0","",minfit,-1.0*minfit);
 118         fHistoVtxAftSelmc->Fit("mcveretxfit","0","",minfit,-1.0*minfit);
 119         Float_t datavertexratio=fHistoVtxAftSeldata->Integral(-1,-1,"width")/fdata->GetParameter(0);
 120         Float_t mcvertexratio=fHistoVtxAftSelmc->Integral(-1,-1,"width")/fmc->GetParameter(0);
 121         */
 122         //Event cut histo
 123         //TH1I* histodata=ecuts_data->GetHistoVtxAftSel();
 124         //TH1I* histomc=ecuts_mc->GetHistoVtxAftSel();
 125
 126         Int_t events_data=ecuts_data->GetHistoVtxAftSel()->GetEntries();
 127         Int_t events_mc=ecuts_mc->GetHistoVtxAftSel()->GetEntries();
 128
 129         if(events_data==0&&events_mc==0)
 130                 return 0;
 131
 132
 133         Float_t datavertexratio=((Float_t)(events_data))/((Float_t)ecuts_data->GetHistoVtxAftSelwithoutZvertexCut()->GetEntries());
 134          Float_t mcvertexratio=((Float_t)(events_mc))/((Float_t)ecuts_mc->GetHistoVtxAftSelwithoutZvertexCut()->GetEntries());
 135          TH1F* fHistoEtaAftSeldata=(TH1F*)ecuts_data->GetHistoEtaAftSel();
 136          TH1F* fHistoEtaAftSelmc=(TH1F*)ecuts_mc->GetHistoEtaAftSel();
 137         flistcanvas->Add(plot_on_canvas("ETA",fHistoEtaAftSeldata,fHistoEtaAftSelmc));
 138
 139
 140          TH1F* fITSclustershistdata=(TH1F*)tcuts_data->GetHistoNclustersITS();
 141           TH1F* fITSclustershistmc=(TH1F*)tcuts_mc->GetHistoNclustersITS();
 142
 143         flistcanvas->Add(plot_on_canvas("NITS",fITSclustershistdata,fITSclustershistmc));
 144         cout<<" data "<<datavertexratio<<" mc "<<mcvertexratio<<endl;
 145
 146         TH2F* hmul=(TH2F*)hman_mc->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul");
 147         hmul->Sumw2();
 148         TCanvas* cbc=new TCanvas("broken chunks","broken chunks",1200,600);
 149         cbc->Divide(2,1);
 150         cbc->cd(1);
 151         hmul->Draw();
 152         cbc->cd(2);
 153         TH1F* nonzero=(TH1F*)hmul->ProjectionX("nonzerotracks",2,-1);
 154         nonzero->SetMarkerColor(kRed);
 155         nonzero->SetMarkerStyle(21);
 156         TH1F* binzero=(TH1F*)hmul->ProjectionX("binzerotracks",1,1);
 157         binzero->SetMarkerColor(kBlack);
 158         binzero->SetMarkerStyle(22);
 159         binzero->Sumw2();
 160         nonzero->Sumw2();
 161         binzero->Divide(nonzero);
 162         TF1* badchunk=new TF1("badchunkfit","pol0",10,40);
 163         binzero->Fit("badchunkfit","R");
 164         Float_t badchunksfraction=badchunk->GetParameter(0);
 165         cout<<"Bad chunks "<<badchunksfraction<<endl;
 166         if(badchunksfraction<0.005)
 167         {
 168                 badchunksfraction=0.0; //only applied if higer than 0.5%
 169                 cout<<"reset bad chunks correction"<<badchunksfraction<<endl;
 170
 171         }
 172
 173
 174         binzero->Draw("E1");
 175         flistcanvas->Add(cbc);
 176         if(TMath::Abs(hmul->GetEntries()/events_mc-1.0)>0.001)
 177                 cout<<"MC merging problem"<<endl;
 178         if(TMath::Abs(hman_data->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul")->GetEntries()/events_data-1.0)>0.001)
 179                 cout<<"Data merging problem"<<endl;
 180
 181         if(fullicorr)
 182                 return (1.0-badchunksfraction)*mcvertexratio/datavertexratio;
 183         else
 184                 return (1.0-badchunksfraction)*mcvertexratio;
 185 }
 186 TCanvas* plot_on_canvas(TString name, TH1* h1,TH1* h2)
 187 {
 188         TCanvas* cvrt=new TCanvas(name.Data(),name.Data(),600,600);
 189         cvrt->cd();
 190         h1->SetLineColor(kRed);
 191         h2->SetLineColor(kBlue);
 192         h1->Sumw2();
 193         h2->Sumw2();
 194         TLegend *lvtr=new TLegend(0.2,0.2,0.5,0.3,"","NDC");
 195         lvtr->SetLineColor(kWhite);
 196         lvtr->AddEntry(h1,"data","l");
 197         lvtr->AddEntry(h2,"MC","l");
 198         if(h1->GetBinContent(h1->GetXaxis()->FindBin(0.0))>0.0&&h2->GetBinContent(h2->GetXaxis()->FindBin(0.0))>0.0)
 199         {
 200                 h1->Scale(1.0/h1->GetBinContent(h1->GetXaxis()->FindBin(0.0)));
 201                 h2->Scale(1.0/h2->GetBinContent(h2->GetXaxis()->FindBin(0.0)));
 202         }
 203         h1->DrawCopy("L");
 204         h2->DrawCopy("Lsame");
 205         lvtr->Draw();
 206         return cvrt;
 207 }
 208