PWGHF/vertexingHF/macros/PlotOutputMCCheck.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2 #include <TCanvas.h>
   3 #include <TH1.h>
   4 #include <TH2.h>
   5 #include <TStyle.h>
   6 #include <TFile.h>
   7 #endif
   8
   9 /* $Id$ */
  10
  11 // Macro to plot the output of AliAnalysisTaskCheckHFMCProd
  12 // Author: F. Prino, prino@to.infn.it
  13
  14
  15 void PlotOutputMCCheck(){
  16
  17   TFile *fil=new TFile("AnalysisResults.root");
  18   TDirectoryFile* df=(TDirectoryFile*)fil->Get("HFMCCheck");
  19   TList* l=(TList*)df->Get("clistHFMCCheck");
  20
  21   TH1F* hNEvents=(TH1F*)l->FindObject("hNEvents");
  22   Int_t nAnalEv=hNEvents->GetBinContent(1);
  23   printf("Number of events= %d\n",nAnalEv);
  24
  25   TCanvas* cv=new TCanvas("cv","Vertex");
  26   cv->Divide(3,3);
  27   cv->cd(1);
  28   TH1F* hSPD3DvX=(TH1F*)l->FindObject("hSPD3DvX");
  29   hSPD3DvX->Draw();
  30   cv->cd(2);
  31   TH1F* hSPD3DvY=(TH1F*)l->FindObject("hSPD3DvY");
  32   hSPD3DvY->Draw();
  33   cv->cd(3);
  34   TH1F* hSPD3DvZ=(TH1F*)l->FindObject("hSPD3DvZ");
  35   hSPD3DvZ->Draw();
  36   cv->cd(4);
  37   TH1F* hSPDZvX=(TH1F*)l->FindObject("hSPDZvX");
  38   hSPDZvX->Draw();
  39   cv->cd(5);
  40   TH1F* hSPDZvY=(TH1F*)l->FindObject("hSPDZvY");
  41   hSPDZvY->Draw();
  42   cv->cd(6);
  43   TH1F* hSPDZvZ=(TH1F*)l->FindObject("hSPDZvZ");
  44   hSPDZvZ->Draw();
  45   cv->cd(7);
  46   TH1F* hTRKvX=(TH1F*)l->FindObject("hTRKvX");
  47   hTRKvX->Draw();
  48   cv->cd(8);
  49   TH1F* hTRKvY=(TH1F*)l->FindObject("hTRKvY");
  50   hTRKvY->Draw();
  51   cv->cd(9);
  52   TH1F* hTRKvZ=(TH1F*)l->FindObject("hTRKvZ");
  53   hTRKvZ->Draw();
  54
  55   TCanvas* c1=new TCanvas("c1","Multipl");
  56   c1->Divide(2,2);
  57   c1->cd(1);
  58   gPad->SetLogy();
  59   TH1F* hPhysPrim=(TH1F*)l->FindObject("hPhysPrim");
  60   hPhysPrim->Draw();
  61   c1->cd(2);
  62   gPad->SetLogy();
  63   TH1F* hTracklets=(TH1F*)l->FindObject("hTracklets");
  64   hTracklets->Draw();
  65   c1->cd(3);
  66   gPad->SetLogy();
  67   TH1F* hTracks=(TH1F*)l->FindObject("hTracks");
  68   hTracks->Draw();
  69   c1->cd(4);
  70   gPad->SetLogy();
  71   TH1F* hSelTracks=(TH1F*)l->FindObject("hSelTracks");
  72   hSelTracks->Draw();
  73
  74   // tracking efficiency
  75   Double_t minEta=-0.8;
  76   Double_t maxEta=0.8;
  77   Double_t minPt=0.5;
  78   Double_t maxPt=2.;
  79   TH3F* hEtaPhiPtGenPi=(TH3F*)l->FindObject("hEtaPhiPtGenPi");
  80   Int_t minEtaBin=hEtaPhiPtGenPi->GetXaxis()->FindBin(minEta+0.00001);
  81   Int_t maxEtaBin=hEtaPhiPtGenPi->GetXaxis()->FindBin(maxEta-0.00001);
  82   cout<<minEtaBin<<"    "<<maxEtaBin<<endl;
  83   Int_t minPtBin=hEtaPhiPtGenPi->GetZaxis()->FindBin(minPt+0.00001);
  84   Int_t maxPtBin=hEtaPhiPtGenPi->GetZaxis()->FindBin(maxPt-0.00001);
  85   cout<<minPtBin<<"    "<<maxPtBin<<endl;
  86   TH1D* hEtaGenPi=hEtaPhiPtGenPi->ProjectionX("hEtaGenPi",0,-1,minPtBin,maxPtBin);
  87   TH1D* hPhiGenPi=hEtaPhiPtGenPi->ProjectionY("hPhiGenPi",minEtaBin,maxEtaBin,minPtBin,maxPtBin);
  88   TH1D* hPtGenPi=hEtaPhiPtGenPi->ProjectionZ("hPtGenPi",minEtaBin,maxEtaBin);
  89   hEtaGenPi->Sumw2();
  90   hPhiGenPi->Sumw2();
  91   hPtGenPi->Sumw2();
  92   hEtaGenPi->GetXaxis()->SetTitle("#eta");
  93   hPhiGenPi->GetXaxis()->SetTitle("#varphi (rad)");
  94   hPtGenPi->GetXaxis()->SetTitle("p_{T} (GeV/c)");
  95   TH3F* hEtaPhiPtRecPi=(TH3F*)l->FindObject("hEtaPhiPtRecPi");
  96   TH1D* hEtaRecPi=hEtaPhiPtRecPi->ProjectionX("hEtaRecPi",0,-1,minPtBin,maxPtBin);
  97   TH1D* hPhiRecPi=hEtaPhiPtRecPi->ProjectionY("hPhiRecPi",minEtaBin,maxEtaBin,minPtBin,maxPtBin);
  98   TH1D* hPtRecPi=hEtaPhiPtRecPi->ProjectionZ("hPtRecPi",minEtaBin,maxEtaBin);
  99   hEtaRecPi->Sumw2();
 100   hPhiRecPi->Sumw2();
 101   hPtRecPi->Sumw2();
 102   hEtaRecPi->GetXaxis()->SetTitle("#eta");
 103   hPhiRecPi->GetXaxis()->SetTitle("#varphi (rad)");
 104   hPtRecPi->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 105   TH1D* hEtaEffPi=(TH1D*)hEtaRecPi->Clone("hEtaEffPi");
 106   hEtaEffPi->Divide(hEtaRecPi,hEtaGenPi,1.,1.,"B");
 107   TH1D* hPhiEffPi=(TH1D*)hPhiRecPi->Clone("hPhiEffPi");
 108   hPhiEffPi->Divide(hPhiRecPi,hPhiGenPi,1.,1.,"B");
 109   TH1D* hPtEffPi=(TH1D*)hPtRecPi->Clone("hPtEffPi");
 110   hPtEffPi->Divide(hPtRecPi,hPtGenPi,1.,1.,"B");
 111
 112   TH3F* hEtaPhiPtGenK=(TH3F*)l->FindObject("hEtaPhiPtGenK");
 113   TH1D* hEtaGenK=hEtaPhiPtGenK->ProjectionX("hEtaGenK",0,-1,minPtBin,maxPtBin);
 114   TH1D* hPhiGenK=hEtaPhiPtGenK->ProjectionY("hPhiGenK",minEtaBin,maxEtaBin,minPtBin,maxPtBin);
 115   TH1D* hPtGenK=hEtaPhiPtGenK->ProjectionZ("hPtGenK",minEtaBin,maxEtaBin);
 116   hEtaGenK->Sumw2();
 117   hPhiGenK->Sumw2();
 118   hPtGenK->Sumw2();
 119   hEtaGenK->GetXaxis()->SetTitle("#eta");
 120   hPhiGenK->GetXaxis()->SetTitle("#varphi (rad)");
 121   hPtGenK->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 122   TH3F* hEtaPhiPtRecK=(TH3F*)l->FindObject("hEtaPhiPtRecK");
 123   TH1D* hEtaRecK=hEtaPhiPtRecK->ProjectionX("hEtaRecK",0,-1,minPtBin,maxPtBin);
 124   TH1D* hPhiRecK=hEtaPhiPtRecK->ProjectionY("hPhiRecK",minEtaBin,maxEtaBin,minPtBin,maxPtBin);
 125   TH1D* hPtRecK=hEtaPhiPtRecK->ProjectionZ("hPtRecK",minEtaBin,maxEtaBin);
 126   hEtaRecK->Sumw2();
 127   hPhiRecK->Sumw2();
 128   hPtRecK->Sumw2();
 129   hEtaRecK->GetXaxis()->SetTitle("#eta");
 130   hPhiRecK->GetXaxis()->SetTitle("#varphi (rad)");
 131   hPtRecK->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 132   TH1D* hEtaEffK=(TH1D*)hEtaRecK->Clone("hEtaEffK");
 133   hEtaEffK->Divide(hEtaRecK,hEtaGenK,1.,1.,"B");
 134   TH1D* hPhiEffK=(TH1D*)hPhiRecK->Clone("hPhiEffK");
 135   hPhiEffK->Divide(hPhiRecK,hPhiGenK,1.,1.,"B");
 136   TH1D* hPtEffK=(TH1D*)hPtRecK->Clone("hPtEffK");
 137   hPtEffK->Divide(hPtRecK,hPtGenK,1.,1.,"B");
 138
 139   TH3F* hEtaPhiPtGenPro=(TH3F*)l->FindObject("hEtaPhiPtGenPro");
 140   TH1D* hEtaGenPro=hEtaPhiPtGenPro->ProjectionX("hEtaGenPro",0,-1,minPtBin,maxPtBin);
 141   TH1D* hPhiGenPro=hEtaPhiPtGenPro->ProjectionY("hPhiGenPro",minEtaBin,maxEtaBin,minPtBin,maxPtBin);
 142   TH1D* hPtGenPro=hEtaPhiPtGenPro->ProjectionZ("hPtGenPro",minEtaBin,maxEtaBin);
 143   hEtaGenPro->Sumw2();
 144   hPhiGenPro->Sumw2();
 145   hPtGenPro->Sumw2();
 146   hEtaGenPro->GetXaxis()->SetTitle("#eta");
 147   hPhiGenPro->GetXaxis()->SetTitle("#varphi (rad)");
 148   hPtGenPro->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 149   TH3F* hEtaPhiPtRecPro=(TH3F*)l->FindObject("hEtaPhiPtRecPro");
 150   TH1D* hEtaRecPro=hEtaPhiPtRecPro->ProjectionX("hEtaRecPro",0,-1,minPtBin,maxPtBin);
 151   TH1D* hPhiRecPro=hEtaPhiPtRecPro->ProjectionY("hPhiRecPro",minEtaBin,maxEtaBin,minPtBin,maxPtBin);
 152   TH1D* hPtRecPro=hEtaPhiPtRecPro->ProjectionZ("hPtRecPro",minEtaBin,maxEtaBin);
 153   hEtaRecPro->Sumw2();
 154   hPhiRecPro->Sumw2();
 155   hPtRecPro->Sumw2();
 156   hEtaRecPro->GetXaxis()->SetTitle("#eta");
 157   hPhiRecPro->GetXaxis()->SetTitle("#varphi (rad)");
 158   hPtRecPro->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 159   TH1D* hEtaEffPro=(TH1D*)hEtaRecPro->Clone("hEtaEffPro");
 160   hEtaEffPro->Divide(hEtaRecPro,hEtaGenPro,1.,1.,"B");
 161   TH1D* hPhiEffPro=(TH1D*)hPhiRecPro->Clone("hPhiEffPro");
 162   hPhiEffPro->Divide(hPhiRecPro,hPhiGenPro,1.,1.,"B");
 163   TH1D* hPtEffPro=(TH1D*)hPtRecPro->Clone("hPtEffPro");
 164   hPtEffPro->Divide(hPtRecPro,hPtGenPro,1.,1.,"B");
 165
 166   TH3F* hEtaPhiPtGenEle=(TH3F*)l->FindObject("hEtaPhiPtGenEle");
 167   TH1D* hEtaGenEle=hEtaPhiPtGenEle->ProjectionX("hEtaGenEle",0,-1,minPtBin,maxPtBin);
 168   TH1D* hPhiGenEle=hEtaPhiPtGenEle->ProjectionY("hPhiGenEle",minEtaBin,maxEtaBin,minPtBin,maxPtBin);
 169   TH1D* hPtGenEle=hEtaPhiPtGenEle->ProjectionZ("hPtGenEle",minEtaBin,maxEtaBin);
 170   hEtaGenEle->Sumw2();
 171   hPhiGenEle->Sumw2();
 172   hPtGenEle->Sumw2();
 173   hEtaGenEle->GetXaxis()->SetTitle("#eta");
 174   hPhiGenEle->GetXaxis()->SetTitle("#varphi (rad)");
 175   hPtGenEle->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 176   TH3F* hEtaPhiPtRecEle=(TH3F*)l->FindObject("hEtaPhiPtRecEle");
 177   TH1D* hEtaRecEle=hEtaPhiPtRecEle->ProjectionX("hEtaRecEle",0,-1,minPtBin,maxPtBin);
 178   TH1D* hPhiRecEle=hEtaPhiPtRecEle->ProjectionY("hPhiRecEle",minEtaBin,maxEtaBin,minPtBin,maxPtBin);
 179   TH1D* hPtRecEle=hEtaPhiPtRecEle->ProjectionZ("hPtRecEle",minEtaBin,maxEtaBin);
 180   hEtaRecEle->Sumw2();
 181   hPhiRecEle->Sumw2();
 182   hPtRecEle->Sumw2();
 183   hEtaRecEle->GetXaxis()->SetTitle("#eta");
 184   hPhiRecEle->GetXaxis()->SetTitle("#varphi (rad)");
 185   hPtRecEle->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 186   TH1D* hEtaEffEle=(TH1D*)hEtaRecEle->Clone("hEtaEffEle");
 187   hEtaEffEle->Divide(hEtaRecEle,hEtaGenEle,1.,1.,"B");
 188   TH1D* hPhiEffEle=(TH1D*)hPhiRecEle->Clone("hPhiEffEle");
 189   hPhiEffEle->Divide(hPhiRecEle,hPhiGenEle,1.,1.,"B");
 190   TH1D* hPtEffEle=(TH1D*)hPtRecEle->Clone("hPtEffEle");
 191   hPtEffEle->Divide(hPtRecEle,hPtGenEle,1.,1.,"B");
 192
 193   for(Int_t iBin=1; iBin<=hPtGenPi->GetNbinsX(); iBin++){
 194     Double_t w=hPtGenPi->GetBinWidth(iBin);
 195     Double_t c=hPtGenPi->GetBinContent(iBin);
 196     hPtGenPi->SetBinContent(iBin,c/w);
 197     c=hPtGenK->GetBinContent(iBin);
 198     hPtGenK->SetBinContent(iBin,c/w);
 199     c=hPtGenPro->GetBinContent(iBin);
 200     hPtGenPro->SetBinContent(iBin,c/w);
 201     c=hPtGenEle->GetBinContent(iBin);
 202     hPtGenEle->SetBinContent(iBin,c/w);
 203
 204   }
 205
 206   hEtaEffPi->SetStats(0);
 207   hPtEffPi->SetStats(0);
 208   hPhiEffPi->SetStats(0);
 209   hEtaEffPi->SetMinimum(0.);
 210   hPtEffPi->SetMinimum(0.);
 211   hPhiEffPi->SetMinimum(0.);
 212   hEtaEffPi->SetMaximum(1.05);
 213   hPtEffPi->SetMaximum(1.05);
 214   hPhiEffPi->SetMaximum(1.05);
 215
 216   hPtGenPi->SetMarkerStyle(27);
 217   hPtEffPi->SetMarkerStyle(27);
 218   hEtaEffPi->SetMarkerStyle(27);
 219   hPhiEffPi->SetMarkerStyle(27);
 220   hPtGenK->SetMarkerStyle(22);
 221   hPtEffK->SetMarkerStyle(22);
 222   hEtaEffK->SetMarkerStyle(22);
 223   hPhiEffK->SetMarkerStyle(22);
 224   hPtGenK->SetMarkerColor(2);
 225   hPtEffK->SetMarkerColor(2);
 226   hEtaEffK->SetMarkerColor(2);
 227   hPhiEffK->SetMarkerColor(2);
 228   hPtGenPro->SetMarkerStyle(25);
 229   hPtEffPro->SetMarkerStyle(25);
 230   hEtaEffPro->SetMarkerStyle(25);
 231   hPhiEffPro->SetMarkerStyle(25);
 232   hPtGenPro->SetMarkerColor(4);
 233   hPtEffPro->SetMarkerColor(4);
 234   hEtaEffPro->SetMarkerColor(4);
 235   hPhiEffPro->SetMarkerColor(4);
 236   hPtGenEle->SetMarkerStyle(20);
 237   hPtEffEle->SetMarkerStyle(20);
 238   hEtaEffEle->SetMarkerStyle(20);
 239   hPhiEffEle->SetMarkerStyle(20);
 240   hPtGenEle->SetMarkerColor(kGreen+1);
 241   hPtEffEle->SetMarkerColor(kGreen+1);
 242   hEtaEffEle->SetMarkerColor(kGreen+1);
 243   hPhiEffEle->SetMarkerColor(kGreen+1);
 244
 245   TCanvas* ctref=new TCanvas("ctref","Track eff",1200,800);
 246   ctref->Divide(2,2);
 247   ctref->cd(1);
 248   gPad->SetLogy();
 249   hPtGenPi->Draw();
 250   ctref->Update();
 251   TPaveStats* st1=(TPaveStats*)hPtGenPi->GetListOfFunctions()->FindObject("stats");
 252   st1->SetY2NDC(0.91);
 253   st1->SetY1NDC(0.76);
 254   st1->SetTextColor(1);
 255   hPtGenK->Draw("sames");
 256   ctref->Update();
 257   TPaveStats* st2=(TPaveStats*)hPtGenK->GetListOfFunctions()->FindObject("stats");
 258   st2->SetY2NDC(0.75);
 259   st2->SetY1NDC(0.60);
 260   st2->SetTextColor(2);
 261   hPtGenPro->Draw("sames");
 262   ctref->Update();
 263   TPaveStats* st3=(TPaveStats*)hPtGenPro->GetListOfFunctions()->FindObject("stats");
 264   st3->SetY2NDC(0.59);
 265   st3->SetY1NDC(0.44);
 266   st3->SetTextColor(4);
 267   hPtGenEle->Draw("sames");
 268   ctref->Update();
 269   TPaveStats* st4=(TPaveStats*)hPtGenEle->GetListOfFunctions()->FindObject("stats");
 270   st4->SetY2NDC(0.43);
 271   st4->SetY1NDC(0.28);
 272   st4->SetTextColor(kGreen+1);
 273   gPad->Modified();
 274   TLegend* leg=new TLegend(0.5,0.5,0.7,0.8);
 275   leg->SetFillColor(0);
 276   leg->SetBorderSize(0);
 277   leg->AddEntry(hPtGenPi,"Pions","P")->SetTextColor(hPtGenPi->GetMarkerColor());
 278   leg->AddEntry(hPtGenK,"Kaons","P")->SetTextColor(hPtGenK->GetMarkerColor());
 279   leg->AddEntry(hPtGenPro,"Protons","P")->SetTextColor(hPtGenPro->GetMarkerColor());
 280   leg->AddEntry(hPtGenEle,"Electrons","P")->SetTextColor(hPtGenEle->GetMarkerColor());
 281   leg->Draw();
 282
 283   gPad->Update();
 284   ctref->cd(2);
 285   hPtEffPi->Draw();
 286   hPtEffPi->GetYaxis()->SetTitle("Efficiency");
 287   hPtEffK->Draw("same");
 288   hPtEffEle->Draw("same");
 289   hPtEffPro->Draw("same");
 290   TLatex* t1=new TLatex(0.62,0.18,Form("%.2f < #eta < %.2f",minEta,maxEta));
 291   t1->SetNDC();
 292   t1->Draw();
 293   ctref->cd(3);
 294   hPhiEffPi->Draw();
 295   hPhiEffPi->GetYaxis()->SetTitle("Efficiency");
 296   hPhiEffK->Draw("same");
 297   hPhiEffEle->Draw("same");
 298   hPhiEffPro->Draw("same");
 299   TLatex* t2=new TLatex(0.18,0.18,Form("%.1f < p_{T} < %.1f GeV/c",minPt,maxPt));
 300   t2->SetNDC();
 301   t2->Draw();
 302   t1->Draw();
 303   ctref->cd(4);
 304   hEtaEffPi->Draw();
 305   hEtaEffPi->GetYaxis()->SetTitle("Efficiency");
 306   hEtaEffK->Draw("same");
 307   hEtaEffEle->Draw("same");
 308   hEtaEffPro->Draw("same");
 309   t2->Draw();
 310
 311   TH1F* hncharmed=(TH1F*)l->FindObject("hncharmed");
 312   TCanvas* cn=new TCanvas("cn","ncharm");
 313   hncharmed->Draw("box");
 314   hncharmed->GetXaxis()->SetTitle("dNch/dy");
 315   hncharmed->GetYaxis()->SetTitle("N Charm hadrons in golden channels");
 316   cn->Update();
 317
 318   TH1F* hnbvsnc=(TH1F*)l->FindObject("hnbvsnc");
 319   TCanvas* cnhf=new TCanvas("cnhf","nb/c");
 320   hnbvsnc->Draw("colztext");
 321   hnbvsnc->GetXaxis()->SetTitle("Nc");
 322   hnbvsnc->GetYaxis()->SetTitle("Nb");
 323   cnhf->Update();
 324   TH2F* hyptD0all=(TH2F*)l->FindObject("hyptD0AllDecay");
 325   TH2F* hyptD0promptall=(TH2F*)l->FindObject("hyptD0promptAllDecay");
 326   TH2F* hyptD0feeddownall=(TH2F*)l->FindObject("hyptD0feeddownAllDecay");
 327   TH2F* hyptD0prompt=(TH2F*)l->FindObject("hyptD0prompt");
 328   TH2F* hyptD0feeddown=(TH2F*)l->FindObject("hyptD0feeddown");
 329   TH2F* hyptD02=(TH2F*)l->FindObject("hyptD02");
 330   TH2F* hyptD04=(TH2F*)l->FindObject("hyptD04");
 331   TH1D* hptD0all=hyptD0all->ProjectionX();
 332   TH1D* hptD0promptall=hyptD0promptall->ProjectionX();
 333   TH1D* hptD0feeddownall=hyptD0feeddownall->ProjectionX();
 334   TH2F* hyptDplusall=(TH2F*)l->FindObject("hyptDplusAllDecay");
 335   TH2F* hyptDpluspromptall=(TH2F*)l->FindObject("hyptDpluspromptAllDecay");
 336   TH2F* hyptDplusfeeddownall=(TH2F*)l->FindObject("hyptDplusfeeddownAllDecay");
 337   TH2F* hyptDplusprompt=(TH2F*)l->FindObject("hyptDplusprompt");
 338   TH2F* hyptDplusfeeddown=(TH2F*)l->FindObject("hyptDplusfeeddown");
 339   TH2F* hyptDplusnonreson=(TH2F*)l->FindObject("hyptDplusnonreson");
 340   TH2F* hyptDplusreson=(TH2F*)l->FindObject("hyptDplusreson");
 341   TH2F* hyptDsall=(TH2F*)l->FindObject("hyptDsAllDecay");
 342   TH2F* hyptDsprompt=(TH2F*)l->FindObject("hyptDsprompt");
 343   TH2F* hyptDsfeeddown=(TH2F*)l->FindObject("hyptDsfeedown");
 344   TH2F* hyptDsphi=(TH2F*)l->FindObject("hyptDsphi");
 345   TH2F* hyptDsK0st=(TH2F*)l->FindObject("hyptDsk0st");
 346   TH2F* hyptDstarall=(TH2F*)l->FindObject("hyptDstarAllDecay");
 347   TH2F* hyptDstarprompt=(TH2F*)l->FindObject("hyptDstarprompt");
 348   TH2F* hyptDstarfeedown=(TH2F*)l->FindObject("hyptDstarfeedown");
 349   TH2F* hyptLcprompt=(TH2F*)l->FindObject("hyptLcprompt");
 350   TH2F* hyptLcfeedown=(TH2F*)l->FindObject("hyptLcfeedown");
 351   TH2F* hyptLcall=(TH2F*)l->FindObject("hyptLcAllDecay");
 352
 353   TH2F* hyptB0all=(TH2F*)l->FindObject("hyptB0AllDecay");
 354   TH2F* hyptBplusall=(TH2F*)l->FindObject("hyptBplusAllDecay");
 355   TH2F* hyptBsall=(TH2F*)l->FindObject("hyptBsAllDecay");
 356   TH2F* hyptBstarall=(TH2F*)l->FindObject("hyptBstarAllDecay");
 357   TH2F* hyptLball=(TH2F*)l->FindObject("hyptLbAllDecay");
 358
 359
 360
 361   TH1F* hD0fonll7=HistoFONLL7TeV();
 362   hD0fonll7->Scale(hptD0all->GetMaximum()/hD0fonll7->GetMaximum());
 363   hD0fonll7->SetLineColor(kGreen+1);
 364   hD0fonll7->SetLineWidth(2);
 365
 366   TH1F* hD0fonll2=HistoFONLL2_76TeV();
 367   hD0fonll2->Scale(hptD0all->GetMaximum()/hD0fonll2->GetMaximum());
 368   hD0fonll2->SetLineColor(kBlue+1);
 369   hD0fonll2->SetLineWidth(2);
 370
 371   TH1F* hptD0pythia=HistoPYTHIA7();
 372   hptD0pythia->Scale(hptD0all->GetMaximum()/hptD0pythia->GetMaximum());
 373   hptD0pythia->SetLineColor(kRed+1);
 374   hptD0pythia->SetLineWidth(2);
 375
 376   hptD0all->SetLineWidth(2);
 377   hptD0all->SetMarkerStyle(20);
 378   hptD0all->SetTitle("");
 379   hptD0all->GetYaxis()->SetTitle("D^{0} dN/dp_{T} (a.u.)");
 380   hptD0all->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 381   hptD0all->SetStats(0);
 382   hptD0all->GetYaxis()->SetTitleOffset(1.2);
 383   hptD0all->GetXaxis()->SetTitleOffset(1.2);
 384
 385   TCanvas* cd0a=new TCanvas("cd0a","D0 spectra",700,700);
 386   cd0a->SetLogy();
 387   cd0a->SetLeftMargin(0.13);
 388   cd0a->SetRightMargin(0.07);
 389   hptD0all->Draw();
 390   hD0fonll7->Draw("lsame");
 391   hD0fonll2->Draw("lsame");
 392   hptD0pythia->Draw("lsame");
 393   TLegend* leg=new TLegend(0.45,0.6,0.89,0.85);
 394   leg->SetFillStyle(0);
 395   leg->SetBorderSize(0);
 396   leg->AddEntry(hptD0all,"MC production","LP");
 397   leg->AddEntry(hD0fonll7,"FONLL, #sqrt{s}=7 TeV","L");
 398   leg->AddEntry(hD0fonll2,"FONLL, #sqrt{s}=2.76 TeV","L");
 399   leg->AddEntry(hptD0pythia,"PYTHIA Perugia0, #sqrt{s}=7 TeV","L");
 400   leg->Draw();
 401
 402   // Prompt and Feeddown
 403   TH1F* hOriginPrompt=(TH1F*)l->FindObject("hOriginPrompt");
 404   TH1F* hOriginFeeddown=(TH1F*)l->FindObject("hOriginFeeddown");
 405
 406   hptD0promptall->SetLineColor(4);
 407   hptD0promptall->SetMarkerColor(4);
 408   hptD0promptall->SetMarkerStyle(26);
 409   hptD0feeddownall->SetLineColor(2);
 410   hptD0feeddownall->SetMarkerColor(2);
 411   hptD0feeddownall->SetMarkerStyle(23)
 412 ;
 413   TCanvas* cprf1=new TCanvas("cprf1","Prompt/Feeddown",700,700);
 414   gPad->SetLogy();
 415   gPad->SetLeftMargin(0.13);
 416   gPad->SetRightMargin(0.07);
 417   hptD0all->Draw();
 418   hptD0promptall->Draw("same");
 419   hptD0feeddownall->Draw("same");
 420   TLegend* leg2=new TLegend(0.4,0.5,0.89,0.85);
 421   leg2->SetFillStyle(0);
 422   leg2->SetBorderSize(0);
 423   leg2->AddEntry(hptD0all,"All D0","LP");
 424   leg2->AddEntry("",Form("Entries=%.0f  <pt>=%.2f GeV/c",hptD0all->Integral(),hptD0all->GetMean()),"");
 425   leg2->AddEntry(hptD0promptall,"Prompt D0","LP");
 426   leg2->AddEntry("",Form("Entries=%.0f  <pt>=%.2f GeV/c",hptD0promptall->Integral(),hptD0promptall->GetMean()),"");
 427   leg2->AddEntry(hptD0feeddownall,"Feeddown D0","LP");
 428   leg2->AddEntry("",Form("Entries=%.0f  <pt>=%.2f GeV/c",hptD0feeddownall->Integral(),hptD0feeddownall->GetMean()),"");
 429   leg2->Draw();
 430   hOriginPrompt->GetXaxis()->SetTitle("Distance of prompt D meson origin to vertex (cm)");
 431   hOriginFeeddown->GetXaxis()->SetTitle("Distance of feed-down D meson origin to vertex (cm)");
 432
 433   TCanvas* cprf2=new TCanvas("cprf2","Origin",1000,600);
 434   cprf2->Divide(2,1);
 435   cprf2->cd(1);
 436   hOriginPrompt->Draw();
 437   cprf2->cd(2);
 438   hOriginFeeddown->Draw();
 439
 440   // Hadrons
 441   Double_t nev=hNEvents->GetBinContent(1);
 442   Double_t nD0=hyptD0all->GetEntries();
 443   Double_t nDp=hyptDplusall->GetEntries();
 444   Double_t nDs=hyptDsall->GetEntries();
 445   Double_t nDst=hyptDstarall->GetEntries();
 446   Double_t nLc=hyptLcall->GetEntries();
 447   Double_t nB0=hyptB0all->GetEntries();
 448   Double_t nBp=hyptBplusall->GetEntries();
 449   Double_t nBs=hyptBsall->GetEntries();
 450   Double_t nBst=hyptBstarall->GetEntries();
 451   Double_t nLb=hyptLball->GetEntries();
 452
 453
 454   Double_t nccbar=0;
 455   Double_t nbbbar=0;
 456   for(Int_t iBinx=1; iBinx<=hnbvsnc->GetNbinsX(); iBinx++){
 457     for(Int_t iBiny=1; iBiny<=hnbvsnc->GetNbinsY(); iBiny++){
 458       Double_t bincentx=hnbvsnc->GetXaxis()->GetBinCenter(iBinx);
 459       Double_t bincenty=hnbvsnc->GetYaxis()->GetBinCenter(iBiny);
 460       Double_t bincont=hnbvsnc->GetBinContent(iBinx,iBiny);
 461       nccbar+=(bincentx*bincont);
 462       nbbbar+=(bincenty*bincont);
 463     }
 464   }
 465
 466
 467   printf("Events =%f\n",nev);
 468   printf("c+cbar =%f\n",nccbar);
 469   printf("D0 =%f\n",nD0);
 470   printf("D+ =%f\n",nDp);
 471   printf("D*+ =%f\n",nDst);
 472   printf("Ds =%f\n",nDs);
 473   printf("Lc =%f\n",nLc);
 474   printf("----------\n");
 475   printf("b+bbar =%f\n",nbbbar);
 476   printf("B0 =%f\n",nB0);
 477   printf("B+ =%f\n",nBp);
 478   printf("B*0 =%f\n",nBst);
 479   printf("Bs =%f\n",nBs);
 480   printf("Lb =%f\n",nLb);
 481
 482   if(nccbar==0) nccbar=nD0+nDp+nDs+nLc;
 483   if(nbbbar==0) nbbbar=nB0+nBp+nBs+nLb;
 484
 485
 486   TH1F* hCharmHad=new TH1F("hCharmHad","",5,-0.5,4.5);
 487   hCharmHad->GetXaxis()->SetBinLabel(1,"D0");
 488   hCharmHad->SetBinContent(1,nD0/nccbar);
 489   hCharmHad->GetXaxis()->SetBinLabel(2,"D+");
 490   hCharmHad->SetBinContent(2,nDp/nccbar);
 491   hCharmHad->GetXaxis()->SetBinLabel(3,"Ds");
 492   hCharmHad->SetBinContent(3,nDs/nccbar);
 493   hCharmHad->GetXaxis()->SetBinLabel(4,"Lc");
 494   hCharmHad->SetBinContent(4,nLc/nccbar);
 495   hCharmHad->GetXaxis()->SetBinLabel(5,"D*+");
 496   hCharmHad->SetBinContent(5,nDst/nccbar);
 497   hCharmHad->SetMinimum(0);
 498   hCharmHad->SetStats(0);
 499   hCharmHad->GetYaxis()->SetTitle("N(species)/N(c+#bar{c})");
 500   hCharmHad->GetYaxis()->SetTitleOffset(1.3);
 501
 502   TH1F* hBeautyHad=new TH1F("hBeautyHad","",5,-0.5,4.5);
 503   hBeautyHad->GetXaxis()->SetBinLabel(1,"B0");
 504   hBeautyHad->SetBinContent(1,nB0/nbbbar);
 505   hBeautyHad->GetXaxis()->SetBinLabel(2,"B+");
 506   hBeautyHad->SetBinContent(2,nBp/nbbbar);
 507   hBeautyHad->GetXaxis()->SetBinLabel(3,"Bs");
 508   hBeautyHad->SetBinContent(3,nBs/nbbbar);
 509   hBeautyHad->GetXaxis()->SetBinLabel(4,"Lb");
 510   hBeautyHad->SetBinContent(4,nLb/nbbbar);
 511   hBeautyHad->GetXaxis()->SetBinLabel(5,"B*0");
 512   hBeautyHad->SetBinContent(5,nBst/nbbbar);
 513   hBeautyHad->SetMinimum(0);
 514   hBeautyHad->SetStats(0);
 515   hBeautyHad->GetYaxis()->SetTitle("N(species)/N(b+#bar{b})");
 516   hBeautyHad->GetYaxis()->SetTitleOffset(1.3);
 517
 518
 519
 520   TCanvas* chad=new TCanvas("chad","Hadrons",1200,600);
 521   chad->Divide(2,1);
 522   chad->cd(1);
 523   hCharmHad->Draw();
 524   chad->cd(2);
 525   hBeautyHad->Draw();
 526
 527
 528
 529   TCanvas* cd0=new TCanvas("cd0","D0");
 530   cd0->Divide(2,2);
 531   cd0->cd(1);
 532   hyptD0prompt->Draw("colz");
 533   cd0->cd(2);
 534   hyptD0feeddown->Draw("colz");
 535   cd0->cd(3);
 536   hyptD02->Draw("colz");
 537   cd0->cd(4);
 538   hyptD04->Draw("colz");
 539
 540
 541   TCanvas* cdplus=new TCanvas("cdplus","Dplus");
 542   cdplus->Divide(2,2);
 543   cdplus->cd(1);
 544   hyptDplusprompt->Draw("colz");
 545   cdplus->cd(2);
 546   hyptDplusfeeddown->Draw("colz");
 547   cdplus->cd(3);
 548   hyptDplusnonreson->Draw("colz");
 549   cdplus->cd(4);
 550   hyptDplusreson->Draw("colz");
 551
 552
 553   TCanvas* cds=new TCanvas("cds","Ds");
 554   cds->Divide(2,2);
 555   cds->cd(1);
 556   hyptDsprompt->Draw("colz");
 557   cds->cd(2);
 558   hyptDsfeeddown->Draw("colz");
 559   cds->cd(3);
 560   hyptDsphi->Draw("colz");
 561   cds->cd(4);
 562   hyptDsK0st->Draw("colz");
 563
 564
 565   TCanvas* cdstlc=new TCanvas("cdstls","Dstar LambdaC");
 566   cdstlc->Divide(2,2);
 567   cdstlc->cd(1);
 568   hyptDstarprompt->Draw("colz");
 569   cdstlc->cd(2);
 570   hyptDstarfeedown->Draw("colz");
 571   cdstlc->cd(3);
 572   hyptLcprompt->Draw("colz");
 573   cdstlc->cd(4);
 574   hyptLcfeedown->Draw("colz");
 575
 576 }
 577
 578
 579
 580 TH1F* HistoFONLL7TeV(){
 581
 582   TH1F* hFONLL7=new TH1F("hFONLLD07TeV","",61,0.,30.5);
 583   Float_t val[61]={
 584     1390542.31,3512269.33,4572233.65,4116353.65,3104057.40,
 585     2185147.21,1507632.40,1038687.03,721889.43,509311.55,
 586     365094.01,265684.43,196609.08,147415.64,112019.94,
 587     86170.85,66997.46,52651.54,41787.55,33458.64,
 588     27012.62,21981.48,18020.69,14873.73,12354.91,
 589     10324.90,8677.32,7331.56,6225.67,5311.74,
 590     4552.30,3918.10,3385.73,2936.79,2556.52,
 591     2233.25,1957.23,1720.64,1517.12,1341.44,
 592     1189.32,1057.17,942.05,841.47,753.33,
 593     675.89,607.68,547.45,494.14,446.87,
 594     404.83,367.39,333.96,304.08,277.30,
 595     253.25,231.62,212.13,194.55,178.66,
 596     164.27};
 597   for(Int_t i=0; i<61; i++) hFONLL7->SetBinContent(i+1,val[i]);
 598   return hFONLL7;
 599 }
 600
 601 TH1F* HistoFONLL2_76TeV(){
 602   TH1F* hFONLL2=new TH1F("hFONLLD02_76TeV","",61,0.,30.5);
 603   Float_t val[61]={
 604     1154043.73,2596175.89,2995937.57,2442988.08,1701598.07,
 605     1122452.81,732896.42,482314.10,322062.75,218493.05,
 606     151653.77,107053.80,76999.08,56239.64,41687.89,
 607     31322.25,23818.86,18326.64,14252.10,11192.03,
 608     8869.29,7089.11,5711.50,4635.50,3788.31,
 609     3116.19,2578.84,2146.51,1796.18,1510.67,
 610     1276.72,1083.93,924.20,791.18,679.88,
 611     586.36,507.49,440.73,383.99,335.54,
 612     294.07,258.43,227.68,201.11,178.07,
 613     158.04,140.58,125.32,111.95,100.20,
 614     89.86,80.73,72.66,65.51,59.16,
 615     53.51,48.48,43.98,39.96,36.36,
 616     33.12};
 617   for(Int_t i=0; i<61; i++) hFONLL2->SetBinContent(i+1,val[i]);
 618   return hFONLL2;
 619 }
 620
 621 TH1F* HistoPYTHIA7(){
 622   TH1F* hPYTHIA7=new TH1F("hPYTHIAD07TeV","",40,0.,20.);
 623   Float_t val[40]={
 624     826307.00,1753264.00,1877464.00,1634664.00,1278586.00,
 625     959137.00,713389.00,535745.00,410250.00,316284.00,
 626     249107.00,198235.00,158832.00,129936.00,106957.00,
 627     89098.00,73690.00,62752.00,53247.00,45004.00,
 628     38475.00,32691.00,28510.00,24516.00,21204.00,
 629     18276.00,15890.00,13702.00,12326.00,10612.00,
 630     9184.00,8028.00,7194.00,6384.00,5767.00,
 631     5102.00,4505.00,3939.00,3578.00,3288.00
 632   };
 633   for(Int_t i=0; i<40; i++) hPYTHIA7->SetBinContent(i+1,val[i]);
 634   return hPYTHIA7;
 635 }