AliPhysics  1168478 (1168478)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ComputeAcceptance.C
Go to the documentation of this file.
1 #if !defined(__CINT__) || defined(__MAKECINT__)
2 #include <TDatabasePDG.h>
3 #include <TRandom3.h>
4 #include <TMath.h>
5 #include <TLorentzVector.h>
6 #include <TParticle.h>
7 #include <TParticle.h>
8 #include <TClonesArray.h>
9 #include <TH1D.h>
10 #include <TH2D.h>
11 #include <TF1.h>
12 #include <TFile.h>
13 #include <TCanvas.h>
14 #include <TLatex.h>
15 #include <TSystem.h>
16 #include <TStyle.h>
17 #include <TPythia6Decayer.h>
18 #include <TPaveStats.h>
19 #endif
20 
24 
25 // Configuration
32 TString fDecayTableFileName="$ALICE_PHYSICS/PWGHF/vertexingHF/macros/decaytable_acc.dat";
34 Int_t totTrials=1000000;
35 
36 
37 Bool_t CountKpi(TClonesArray *array, Int_t nentries, Int_t &nPions, Int_t &nKaons, Int_t &nPionsInAcc, Int_t &nKaonsInAcc);
39 Bool_t CountPKpi(TClonesArray *array, Int_t nentries, Int_t &nPions, Int_t &nKaons, Int_t &nProtons, Int_t &nPionsInAcc, Int_t &nKaonsInAcc, Int_t &nProtonsInAcc);
40 
41 
42 // Pt-shape histograms
56 
57 
58 
60  // main function
61 
62  gSystem->Load("liblhapdf.so"); // Parton density functions
63  gSystem->Load("libEGPythia6.so"); // TGenerator interface
64  gSystem->Load("libpythia6.so"); // Pythia
65  gSystem->Load("libAliPythia6"); // ALICE specific implementations
66 
67  TPythia6Decayer* pdec=TPythia6Decayer::Instance();
68  if(fDecayTableFileName.CompareTo("")!=0){
69  if(fDecayTableFileName.Contains("ALICE_PHYSICS")){
70  gSystem->Exec(Form("cp %s .",fDecayTableFileName.Data()));
71  fDecayTableFileName.ReplaceAll("$ALICE_PHYSICS/PWGHF/vertexingHF/macros/","./");
72  }
73  pdec->SetDecayTableFile(fDecayTableFileName.Data());
74  pdec->ReadDecayTable();
75  }
76  pdec->Init();
77 
78  Int_t pdgCode=0;
79  Int_t nPionDau=-1;
80  Int_t nProtonDau=-1;
81  Int_t nKaonDau=-1;
82  TString outFileName="Acceptance_Toy_";
83  if(fDDecay==kD0Kpi){
84  pdgCode=421;
85  nPionDau=1;
86  nKaonDau=1;
87  nProtonDau=0;
88  outFileName.Append("D0Kpi_");
89  }else if(fDDecay==kDplusKpipi){
90  pdgCode=411;
91  nPionDau=2;
92  nKaonDau=1;
93  nProtonDau=0;
94  outFileName.Append("DplusKpipi_");
95  }else if(fDDecay==kDstarD0pi){
96  pdgCode=413;
97  nPionDau=2;
98  nKaonDau=1;
99  nProtonDau=0;
100  outFileName.Append("DStarD0pi_");
101  }else if(fDDecay==kDsKKpi){
102  pdgCode=431;
103  nPionDau=1;
104  nKaonDau=2;
105  nProtonDau=0;
106  outFileName.Append("DsKKpi_");
107  }else if(fDDecay==kLcpKpi){
108  pdgCode=4122;
109  nPionDau=1;
110  nKaonDau=1;
111  nProtonDau=1;
112  outFileName.Append("LcpKpi_");
113  }else if(fDDecay==kLcK0Sp){
114  pdgCode=4122;
115  nPionDau=2;
116  nKaonDau=0;
117  nProtonDau=1;
118  outFileName.Append("LcK0Sp_");
119  }else{
120  printf("ERROR: Wrong decay selected\n");
121  return;
122  }
123  if(fOptionYFiducial==kFixedY) outFileName.Append(Form("yfid%02d_",(Int_t)(fYMaxFidAccCut*10)));
124  else outFileName.Append("yfidPtDep_");
125  outFileName.Append(Form("etaDau%02d_",(Int_t)(fEtaMaxDau*10)));
126  outFileName.Append(Form("ptDau%d_",(Int_t)(fPtMinDau*1000)));
127  TDatabasePDG* db=TDatabasePDG::Instance();
128  Float_t massD=db->GetParticle(pdgCode)->Mass();
129  TClonesArray *array = new TClonesArray("TParticle",100);
130 
131  TH2D* hPtVsYGen=new TH2D("hPtVsYGen","",400,0.,40.,20.,-1.,1.);
132  hPtVsYGen->GetXaxis()->SetTitle("p_{T} (GeV/c)");
133  hPtVsYGen->GetYaxis()->SetTitle("y");
134  TH2D* hPtVsYGenLimAcc=new TH2D("hPtVsYGenLimAcc","",400,0.,40.,20.,-1.,1.);
135  hPtVsYGenLimAcc->GetXaxis()->SetTitle("p_{T} (GeV/c)");
136  hPtVsYGenLimAcc->GetYaxis()->SetTitle("y");
137  TH2D* hPtVsYGenAcc=new TH2D("hPtVsYGenAcc","",400,0.,40.,20.,-1.,1.);
138  hPtVsYGenAcc->GetXaxis()->SetTitle("p_{T} (GeV/c)");
139  hPtVsYGenAcc->GetYaxis()->SetTitle("y");
140 
141  TF1* funcPt=0x0;
142  TH1D* histPt=0x0;
143  if(fPtShape==kFONLL8TeV){
144  funcPt=new TF1("fFONLL","[0]*x/TMath::Power((1+TMath::Power(x/[1],[3])),[2])",0.,40.);
145  funcPt->SetParameters(0.518046,3.01138,3.38914,1.75899); // Prompt
146  outFileName.Append("FONLL8ptshape.root");
147  }else if(fPtShape==kFONLL8TeVfeeddown){
148  funcPt=new TF1("fFONLL","[0]*x/TMath::Power((1+TMath::Power(x/[1],[3])),[2])",0.,40.);
149  funcPt->SetParameters(0.398252, 3.9603, 3.915, 1.51853); // FeedDown
150  outFileName.Append("FONLL8ptshapeFeedDown.root");
151  }else if(fPtShape==kFONLL7TeV){
152  funcPt=new TF1("fFONLL","[0]*x/TMath::Power((1+TMath::Power(x/[1],[3])),[2])",0.,40.);
153  funcPt->SetParameters(0.322643,2.96275,2.30301,2.5);
154  outFileName.Append("FONLL7ptshape.root");
155  }else if(fPtShape==kFONLL5TeV){
156  funcPt=new TF1("fFONLL","[0]*x/TMath::Power((1+TMath::Power(x/[1],[3])),[2])",0.,40.);
157  funcPt->SetParameters(0.302879,2.9750,3.68139,1.68855);
158  outFileName.Append("FONLL5ptshape.root");
159  }else if(fPtShape==kPythia7TeV){
160  funcPt=new TF1("fFONLL","[0]*x/TMath::Power((1+TMath::Power(x/[1],[3])),[2])",0.,40.);
161  funcPt->SetParameters(0.322643,1.94635,1.40463,2.5);
162  outFileName.Append("PYTHIA7ptshape.root");
163  }else if(fPtShape==kFONLL13TeVprompt){
164  if(fDDecay==kDplusKpipi){
165  histPt = LoadFONLL13TeV_promptDplus();
166  outFileName.Append("promptDplus");
167  }else if (fDDecay==kDstarD0pi){
168  histPt = LoadFONLL13TeV_promptDstar();
169  outFileName.Append("promptDstar");
170  }else{
171  histPt = LoadFONLL13TeV_promptD0();
172  outFileName.Append("promptD0");
173  }
174  outFileName.Append("FONLL13ptshape.root");
175 
176  }else if (fPtShape==kFONLL13TeVfeeddown){
177  if (fDDecay==kDstarD0pi){
178  histPt = LoadFONLL13TeV_feeddownDstar();
179  outFileName.Append("feeddownDstar");
180  }else{
181  histPt = LoadFONLL13TeV_feeddownD();
182  outFileName.Append("feeddownD");
183  }
184  outFileName.Append("FONLL13ptshape.root");
185  }else if(fPtShape==kPythia13TeVprompt){
186  if(fDDecay==kDplusKpipi){
187  histPt = LoadPYTHIA13TeV_promptDplus();
188  outFileName.Append("promptDplus");
189  }else if (fDDecay==kDstarD0pi){
190  histPt = LoadPYTHIA13TeV_promptDstar();
191  outFileName.Append("promptDstar");
192  }else if (fDDecay==kDsKKpi){
193  histPt = LoadPYTHIA13TeV_promptDs();
194  outFileName.Append("promptDs");
195  }else{
196  histPt = LoadPYTHIA13TeV_promptD0();
197  outFileName.Append("promptD0");
198  }
199  outFileName.Append("PYTHIA13ptshape.root");
200  }else if (fPtShape==kPythia13TeVfeeddown){
201  if(fDDecay==kDplusKpipi){
203  outFileName.Append("feeddownDplus");
204  }else if (fDDecay==kDstarD0pi){
206  outFileName.Append("feeddownDstar");
207  }else if (fDDecay==kDsKKpi){
208  histPt = LoadPYTHIA13TeV_feeddownDs();
209  outFileName.Append("feeddownDs");
210  }else{
211  histPt = LoadPYTHIA13TeV_feeddownD0();
212  outFileName.Append("feeddownD0");
213  }
214  outFileName.Append("PYTHIA13ptshape.root");
215  }else{
216  funcPt=new TF1("fFlat","pol0",0.,40.);
217  funcPt->SetParameter(0,1.);
218  outFileName.Append("flatpt.root");
219  }
220 
221  if (funcPt) funcPt->SetNpx(10000);
222 
223  TRandom3* gener=new TRandom3(0);
224  TLorentzVector* vec=new TLorentzVector();
225 
226 
227  for(Int_t itry=0; itry<totTrials; itry++){
228  if(itry%10000==0) printf("Event %d\n",itry);
229  Float_t ptD = funcPt ? funcPt->GetRandom() : histPt->GetRandom();
230  Float_t phiD=gener->Rndm()*2*TMath::Pi();
231  Float_t yD=gener->Rndm()*2.-1.; // flat in -1<y<1
232  Float_t px=ptD*TMath::Cos(phiD);
233  Float_t py=ptD*TMath::Sin(phiD);
234  Float_t mt=TMath::Sqrt(massD*massD+ptD*ptD);
235  Float_t pz=mt*TMath::SinH(yD);
236  Float_t E=TMath::Sqrt(massD*massD+px*px+py*py+pz*pz);
237 
238  // TLorentzVector* vec=new TLorentzVector(px,py,pz,E);
239  vec->SetPxPyPzE(px,py,pz,E);
240  pdec->Decay(pdgCode,vec);
241  array->Clear();
242  Int_t nentries = pdec->ImportParticles(array);
243  TParticle* dmes=(TParticle*)array->At(0);
244  Int_t nDaughters=dmes->GetNDaughters();
245  if(fDDecay==kD0Kpi && nDaughters!=2) continue;
246  if(fDDecay==kLcK0Sp && nentries>6) continue;
247  Int_t nPionsInAcc=0;
248  Int_t nProtonsInAcc=0;
249  Int_t nKaonsInAcc=0;
250  Int_t nPions=0;
251  Int_t nProtons=0;
252  Int_t nKaons=0;
253  Bool_t isOk=CountPKpi(array,nentries,nPions,nKaons,nProtons,nPionsInAcc,nKaonsInAcc,nProtonsInAcc);
254 
255  if(isOk){
256  if(nPions==nPionDau && nKaons==nKaonDau && nProtons==nProtonDau){
257  hPtVsYGen->Fill(ptD,yD);
258  if(TMath::Abs(yD)<0.5){
259  hPtVsYGenLimAcc->Fill(ptD,yD);
260  }
261  if(IsInFiducialAcceptance(ptD,yD)){
262  if(nPionsInAcc==nPionDau && nKaonsInAcc==nKaonDau && nProtonsInAcc==nProtonDau){
263  hPtVsYGenAcc->Fill(ptD,yD);
264  }
265  }
266  }
267  }
268  // delete vec;
269  }
270 
271  TH1D* hPtGenAcc=(TH1D*)hPtVsYGenAcc->ProjectionX("hPtGenAcc");
272  hPtGenAcc->GetYaxis()->SetTitle("Entries");
273  TH1D* hPtGenLimAcc=(TH1D*)hPtVsYGenLimAcc->ProjectionX("hPtGenLimAcc");
274  hPtGenLimAcc->GetYaxis()->SetTitle("Entries");
275  hPtGenAcc->Sumw2();
276  hPtGenLimAcc->Sumw2();
277  TH1D* hAccVsPt=(TH1D*)hPtGenAcc->Clone("hAccVsPt");
278  hAccVsPt->Divide(hPtGenAcc,hPtGenLimAcc,1,1,"B");
279  hAccVsPt->GetYaxis()->SetTitle("Acceptance");
280  hAccVsPt->SetStats(0);
281 
282  TCanvas* c2d=new TCanvas("c2d","Pt vs y",1200,600);
283  c2d->Divide(3,1);
284  c2d->cd(1);
285  hPtVsYGen->Draw("colz");
286  c2d->cd(2);
287  hPtVsYGenLimAcc->Draw("colz");
288  c2d->cd(3);
289  hPtVsYGenAcc->Draw("colz");
290 
291  TCanvas* c1d=new TCanvas("c1d","Acceptance",1200,600);
292  c1d->Divide(2,1);
293  c1d->cd(1);
294  hPtGenLimAcc->Draw("");
295  Double_t ymax=1.2*TMath::Max(hPtGenLimAcc->GetMaximum(),hPtGenAcc->GetMaximum());
296  hPtGenLimAcc->SetMaximum(ymax);
297  gPad->Update();
298  TPaveStats *st1=(TPaveStats*)hPtGenLimAcc->GetListOfFunctions()->FindObject("stats");
299  st1->SetY1NDC(0.71);
300  st1->SetY2NDC(0.9);
301  hPtGenAcc->SetLineColor(kRed+1);
302  hPtGenAcc->Draw("sames");
303  gPad->Update();
304  TPaveStats *st2=(TPaveStats*)hPtGenAcc->GetListOfFunctions()->FindObject("stats");
305  st2->SetY1NDC(0.51);
306  st2->SetY2NDC(0.7);
307  st2->SetTextColor(hPtGenAcc->GetLineColor());
308  gPad->Modified();
309  c1d->cd(2);
310  hAccVsPt->Draw();
311 
312 
313  TFile* outfil=new TFile(outFileName.Data(),"recreate");
314  hPtVsYGen->Write();
315  hPtVsYGenLimAcc->Write();
316  hPtVsYGenAcc->Write();
317  hPtGenLimAcc->Write();
318  hPtGenAcc->Write();
319  hAccVsPt->Write();
320  outfil->Close();
321 
322 }
323 
324 //___________________________________________________
326  // check fiducial acceptance
327 
329  if(TMath::Abs(y) > fYMaxFidAccCut) return kFALSE;
330  else return kTRUE;
331  }
332 
333  if(pt > 5.) {
334  if (TMath::Abs(y) > 0.8) return kFALSE;
335  } else {
336  Double_t maxFiducialY = -0.2/15*pt*pt+1.9/15*pt+0.5;
337  Double_t minFiducialY = 0.2/15*pt*pt-1.9/15*pt-0.5;
338  if (y < minFiducialY || y > maxFiducialY) return kFALSE;
339  }
340 
341  return kTRUE;
342 }
343 //___________________________________________________
344 Bool_t CountKpi(TClonesArray *array, Int_t nentries, Int_t &nPions, Int_t &nKaons, Int_t &nPionsInAcc, Int_t &nKaonsInAcc){
345  // count K and pi in Acc
346 
347  TParticle* dmes=(TParticle*)array->At(0);
348  Double_t sumPx=0;
349  Double_t sumPy=0;
350  Double_t sumPz=0;
351  for(int j=0; j<nentries; j++){
352  TParticle * o = (TParticle*)array->At(j);
353  Int_t pdgdau=TMath::Abs(o->GetPdgCode());
354  if(fDebugLevel>0) printf("%d ",pdgdau);
355  if(pdgdau==130) {
356  if(fDebugLevel>0) printf("K0 dacaying into K0L\n");
357  return kFALSE;
358  }
359  Float_t ptdau=TMath::Sqrt(o->Px()*o->Px()+o->Py()*o->Py());
360  Float_t etadau=o->Eta();
361  if(pdgdau==211){
362  nPions++;
363  sumPx+=o->Px();
364  sumPy+=o->Py();
365  sumPz+=o->Pz();
366  }
367  if(pdgdau==321){
368  nKaons++;
369  sumPx+=o->Px();
370  sumPy+=o->Py();
371  sumPz+=o->Pz();
372  }
373  if(TMath::Abs(etadau)<fEtaMaxDau && ptdau>fPtMinDau){
374  if(pdgdau==211) nPionsInAcc++;
375  if(pdgdau==321) nKaonsInAcc++;
376  }
377  }
378  if(fDebugLevel>0) printf("\n");
379  if(TMath::Abs(sumPx-dmes->Px())>0.001 ||
380  TMath::Abs(sumPy-dmes->Py())>0.001 ||
381  TMath::Abs(sumPz-dmes->Pz())>0.001){
382  printf("Momentum conservation violation\n");
383  return kFALSE;
384  }
385  return kTRUE;
386 }
387 
388 
389 //___________________________________________________
390 Bool_t CountPKpi(TClonesArray *array, Int_t nentries, Int_t &nPions, Int_t &nKaons, Int_t &nProtons, Int_t &nPionsInAcc, Int_t &nKaonsInAcc, Int_t &nProtonsInAcc){
391  // count K and pi in Acc
392 
393  TParticle* dmes=(TParticle*)array->At(0);
394  Double_t sumPx=0;
395  Double_t sumPy=0;
396  Double_t sumPz=0;
397 
398  for(int j=0; j<nentries; j++){
399  TParticle * o = (TParticle*)array->At(j);
400  Int_t pdgdau=TMath::Abs(o->GetPdgCode());
401  if(fDebugLevel>0) printf("%d ",pdgdau);
402  if(pdgdau==130) {
403  if(fDebugLevel>0) printf("K0 dacaying into K0L\n");
404  return kFALSE;
405  }
406  Float_t ptdau=TMath::Sqrt(o->Px()*o->Px()+o->Py()*o->Py());
407  Float_t etadau=o->Eta();
408  if(pdgdau==211){
409  nPions++;
410  sumPx+=o->Px();
411  sumPy+=o->Py();
412  sumPz+=o->Pz();
413  }
414  if(pdgdau==321){
415  nKaons++;
416  sumPx+=o->Px();
417  sumPy+=o->Py();
418  sumPz+=o->Pz();
419  }
420  if(pdgdau==2212){
421  nProtons++;
422  sumPx+=o->Px();
423  sumPy+=o->Py();
424  sumPz+=o->Pz();
425  }
426  if(TMath::Abs(etadau)<fEtaMaxDau && ptdau>fPtMinDau){
427  if(pdgdau==211) nPionsInAcc++;
428  if(pdgdau==321) nKaonsInAcc++;
429  if(pdgdau==2212) nProtonsInAcc++;
430  }
431  }
432  if(fDebugLevel>0) printf("\n");
433  if(TMath::Abs(sumPx-dmes->Px())>0.001 ||
434  TMath::Abs(sumPy-dmes->Py())>0.001 ||
435  TMath::Abs(sumPz-dmes->Pz())>0.001){
436  printf("Momentum conservation violation\n");
437  return kFALSE;
438  }
439  return kTRUE;
440 }
441 
442 
443 
444 //___________________________________________________
446 {
447  TH1D *hFONLL13 = new TH1D("hFONLL13TeV_D0", "", 80, 0., 40.);
448  Float_t val[80] = {
449  1.4686e+08, 3.9542e+08, 5.4901e+08, 5.2166e+08, 4.1083e+08, 2.9968e+08, 2.1299e+08, 1.5057e+08, 1.0701e+08, 7.6919e+07,
450  5.6121e+07, 4.1546e+07, 3.1184e+07, 2.3715e+07, 1.8253e+07, 1.4206e+07, 1.1175e+07, 8.8774e+06, 7.1169e+06, 5.7544e+06,
451  4.6899e+06, 3.8509e+06, 3.1841e+06, 2.6499e+06, 2.2189e+06, 1.8687e+06, 1.5823e+06, 1.3467e+06, 1.1516e+06, 9.8933e+05,
452  8.5356e+05, 7.3942e+05, 6.4302e+05, 5.6122e+05, 4.9154e+05, 4.3193e+05, 3.8074e+05, 3.3662e+05, 2.9846e+05, 2.6535e+05,
453  2.3653e+05, 2.1136e+05, 1.8932e+05, 1.6997e+05, 1.5293e+05, 1.3789e+05, 1.2458e+05, 1.1277e+05, 1.0227e+05, 9.2913e+04,
454  8.4560e+04, 7.7087e+04, 7.0388e+04, 6.4371e+04, 5.8956e+04, 5.4074e+04, 4.9667e+04, 4.5680e+04, 4.2068e+04, 3.8791e+04,
455  3.5813e+04, 3.3103e+04, 3.0633e+04, 2.8380e+04, 2.6321e+04, 2.4436e+04, 2.2710e+04, 2.1127e+04, 1.9673e+04, 1.8336e+04,
456  1.7106e+04, 1.5972e+04, 1.4926e+04, 1.3961e+04, 1.3069e+04, 1.2243e+04, 1.1479e+04, 1.0771e+04, 1.0113e+04, 9.5031e+03
457  };
458  for (Int_t ibin=0; ibin<80; ++ibin) hFONLL13->SetBinContent(ibin+1, val[ibin]);
459 
460  return hFONLL13;
461 }
462 
463 
464 
465 //___________________________________________________
467 {
468  TH1D *hFONLL13 = new TH1D("hFONLL13TeV_Dplus", "", 80, 0., 40.);
469  Float_t val[80] = {
470  1.5242e+08, 3.9396e+08, 5.3767e+08, 5.1263e+08, 4.0706e+08, 2.9926e+08, 2.1414e+08, 1.5230e+08, 1.0879e+08, 7.8507e+07,
471  5.7460e+07, 4.2651e+07, 3.2093e+07, 2.4462e+07, 1.8865e+07, 1.4708e+07, 1.1588e+07, 9.2185e+06, 7.3998e+06, 5.9900e+06,
472  4.8871e+06, 4.0167e+06, 3.3240e+06, 2.7686e+06, 2.3200e+06, 1.9552e+06, 1.6566e+06, 1.4107e+06, 1.2070e+06, 1.0374e+06,
473  8.9550e+05, 7.7610e+05, 6.7519e+05, 5.8953e+05, 5.1652e+05, 4.5404e+05, 4.0037e+05, 3.5409e+05, 3.1405e+05, 2.7929e+05,
474  2.4902e+05, 2.2258e+05, 1.9943e+05, 1.7909e+05, 1.6117e+05, 1.4535e+05, 1.3135e+05, 1.1892e+05, 1.0787e+05, 9.8024e+04,
475  8.9230e+04, 8.1359e+04, 7.4302e+04, 6.7963e+04, 6.2256e+04, 5.7112e+04, 5.2465e+04, 4.8261e+04, 4.4452e+04, 4.0996e+04,
476  3.7854e+04, 3.4995e+04, 3.2390e+04, 3.0012e+04, 2.7838e+04, 2.5849e+04, 2.4027e+04, 2.2355e+04, 2.0820e+04, 1.9408e+04,
477  1.8108e+04, 1.6910e+04, 1.5805e+04, 1.4785e+04, 1.3842e+04, 1.2969e+04, 1.2161e+04, 1.1411e+04, 1.0716e+04, 1.0070e+04
478  };
479  for (Int_t ibin=0; ibin<80; ++ibin) hFONLL13->SetBinContent(ibin+1, val[ibin]);
480 
481  return hFONLL13;
482 }
483 
484 
485 
486 //___________________________________________________
488 {
489  TH1D *hFONLL13 = new TH1D("hFONLL13TeV_Dstar", "", 80, 0., 40.);
490  Float_t val[80] = {
491  1.2433e+08, 3.4512e+08, 5.0662e+08, 5.1020e+08, 4.2016e+08, 3.1661e+08, 2.3064e+08, 1.6632e+08, 1.2007e+08, 8.7340e+07,
492  6.4329e+07, 4.8000e+07, 3.6287e+07, 2.7772e+07, 2.1492e+07, 1.6808e+07, 1.3278e+07, 1.0589e+07, 8.5178e+06, 6.9084e+06,
493  5.6462e+06, 4.6479e+06, 3.8519e+06, 3.2124e+06, 2.6951e+06, 2.2738e+06, 1.9285e+06, 1.6437e+06, 1.4076e+06, 1.2108e+06,
494  1.0459e+06, 9.0711e+05, 7.8968e+05, 6.8993e+05, 6.0484e+05, 5.3197e+05, 4.6933e+05, 4.1529e+05, 3.6850e+05, 3.2786e+05,
495  2.9246e+05, 2.6152e+05, 2.3441e+05, 2.1058e+05, 1.8958e+05, 1.7104e+05, 1.5461e+05, 1.4003e+05, 1.2706e+05, 1.1550e+05,
496  1.0517e+05, 9.5920e+04, 8.7626e+04, 8.0172e+04, 7.3461e+04, 6.7409e+04, 6.1940e+04, 5.6992e+04, 5.2508e+04, 4.8437e+04,
497  4.4736e+04, 4.1367e+04, 3.8296e+04, 3.5492e+04, 3.2929e+04, 3.0583e+04, 2.8433e+04, 2.6460e+04, 2.4648e+04, 2.2981e+04,
498  2.1446e+04, 2.0032e+04, 1.8727e+04, 1.7522e+04, 1.6407e+04, 1.5376e+04, 1.4420e+04, 1.3535e+04, 1.2713e+04, 1.1950e+04
499  };
500  for (Int_t ibin=0; ibin<80; ++ibin) hFONLL13->SetBinContent(ibin+1, val[ibin]);
501 
502  return hFONLL13;
503 }
504 
505 
506 
507 //___________________________________________________
509 { // B->D with B.R=1
510  TH1D *hFONLL13 = new TH1D("FONLL13TeV_feeddownD", "", 80, 0., 40.);
511  Float_t val[80] = {
512  1.0310e+07, 2.6790e+07, 3.4480e+07, 3.4430e+07, 3.0200e+07, 2.4740e+07, 1.9600e+07, 1.5300e+07, 1.1880e+07, 9.2260e+06,
513  7.1950e+06, 5.6470e+06, 4.4650e+06, 3.5580e+06, 2.8570e+06, 2.3120e+06, 1.8860e+06, 1.5490e+06, 1.2810e+06, 1.0660e+06,
514  8.9260e+05, 7.5170e+05, 6.3650e+05, 5.4160e+05, 4.6310e+05, 3.9770e+05, 3.4300e+05, 2.9700e+05, 2.5820e+05, 2.2520e+05,
515  1.9710e+05, 1.7310e+05, 1.5240e+05, 1.3460e+05, 1.1920e+05, 1.0590e+05, 9.4290e+04, 8.4160e+04, 7.5290e+04, 6.7500e+04,
516  6.0650e+04, 5.4610e+04, 4.9270e+04, 4.4530e+04, 4.0320e+04, 3.6580e+04, 3.3240e+04, 3.0250e+04, 2.7570e+04, 2.5170e+04,
517  2.3020e+04, 2.1070e+04, 1.9320e+04, 1.7740e+04, 1.6310e+04, 1.5010e+04, 1.3830e+04, 1.2760e+04, 1.1790e+04, 1.0900e+04,
518  1.0090e+04, 9.3520e+03, 8.6760e+03, 8.0560e+03, 7.4890e+03, 6.9670e+03, 6.4880e+03, 6.0470e+03, 5.6410e+03, 5.2670e+03,
519  4.9220e+03, 4.6030e+03, 4.3080e+03, 4.0350e+03, 3.7830e+03, 3.5480e+03, 3.3310e+03, 3.1290e+03, 2.9410e+03, 2.7670e+03
520  };
521  for (Int_t ibin=0; ibin<80; ++ibin) hFONLL13->SetBinContent(ibin+1, val[ibin]);
522 
523  return hFONLL13;
524 }
525 
526 
527 
528 //___________________________________________________
530 { // B->D* with B.R=1
531  TH1D *hFONLL13 = new TH1D("FONLL13TeV_feeddownDstar", "", 80, 0., 40.);
532  Float_t val[80] = {
533  9.5260e+06, 2.5070e+07, 3.2890e+07, 3.3540e+07, 3.0010e+07, 2.5000e+07, 2.0090e+07, 1.5860e+07, 1.2430e+07, 9.7280e+06,
534  7.6360e+06, 6.0240e+06, 4.7820e+06, 3.8240e+06, 3.0800e+06, 2.5000e+06, 2.0430e+06, 1.6810e+06, 1.3930e+06, 1.1610e+06,
535  9.7390e+05, 8.2120e+05, 6.9610e+05, 5.9300e+05, 5.0750e+05, 4.3620e+05, 3.7650e+05, 3.2620e+05, 2.8370e+05, 2.4760e+05,
536  2.1680e+05, 1.9050e+05, 1.6780e+05, 1.4830e+05, 1.3140e+05, 1.1680e+05, 1.0400e+05, 9.2870e+04, 8.3100e+04, 7.4530e+04,
537  6.6990e+04, 6.0330e+04, 5.4440e+04, 4.9220e+04, 4.4580e+04, 4.0440e+04, 3.6760e+04, 3.3460e+04, 3.0510e+04, 2.7860e+04,
538  2.5470e+04, 2.3330e+04, 2.1390e+04, 1.9640e+04, 1.8060e+04, 1.6630e+04, 1.5320e+04, 1.4140e+04, 1.3060e+04, 1.2080e+04,
539  1.1190e+04, 1.0370e+04, 9.6190e+03, 8.9330e+03, 8.3050e+03, 7.7270e+03, 7.1970e+03, 6.7080e+03, 6.2580e+03, 5.8440e+03,
540  5.4610e+03, 5.1080e+03, 4.7810e+03, 4.4790e+03, 4.1980e+03, 3.9390e+03, 3.6970e+03, 3.4740e+03, 3.2650e+03, 3.0720e+03
541  };
542  for (Int_t ibin=0; ibin<80; ++ibin) hFONLL13->SetBinContent(ibin+1, val[ibin]);
543 
544  return hFONLL13;
545 }
546 
547 
548 
549 //___________________________________________________
551 {
552  TH1D *hPYTHIA13 = new TH1D("hPYTHIA13TeV_D0", "", 40, 0., 40.);
553  Float_t val[40] = {
554  2617743, 3763836, 2235903, 1140807, 587707, 317881, 180756, 107762, 67299,42882,
555  28550, 19690, 13817, 9754, 7205, 5299, 4100, 3177, 2386,1822,
556  1484, 1248, 933, 798, 639, 534, 435, 405, 304,260,
557  210, 196, 161, 143, 139, 99, 89, 80, 66,63
558  };
559  for (Int_t ibin=0; ibin<40; ++ibin) hPYTHIA13->SetBinContent(ibin+1, val[ibin]);
560 
561  return hPYTHIA13;
562 }
563 
564 
565 
566 //___________________________________________________
568 {
569  TH1D *hPYTHIA13 = new TH1D("hPYTHIA13TeV_Dplus", "", 40, 0., 40.);
570  Float_t val[40] = {
571  1192379, 1736283, 1047485, 540499, 281523, 152153, 87150, 52200, 32709, 20991,
572  13986, 9611, 6832, 4705, 3484, 2650, 2001, 1518, 1197, 924,
573  735, 576, 497, 402, 292, 261, 214, 173, 169, 150,
574  106, 102, 95, 70, 67, 47, 53, 46, 35, 35
575  };
576  for (Int_t ibin=0; ibin<40; ++ibin) hPYTHIA13->SetBinContent(ibin+1, val[ibin]);
577 
578  return hPYTHIA13;
579 }
580 
581 
582 
583 //___________________________________________________
585 {
586  TH1D *hPYTHIA13 = new TH1D("hPYTHIA13TeV_Dstar", "", 40, 0., 40.);
587  Float_t val[40] = {
588  922985, 1419749, 904587, 484407, 257862, 142386, 82177, 49839, 31103, 20422,
589  13493, 9272, 6549, 4746, 3344, 2575, 1959, 1493, 1161, 935,
590  672, 585, 456, 415, 328, 260, 199, 180, 161, 134,
591  109, 97, 103, 72, 63, 56, 40, 46, 38, 27
592  };
593  for (Int_t ibin=0; ibin<40; ++ibin) hPYTHIA13->SetBinContent(ibin+1, val[ibin]);
594 
595  return hPYTHIA13;
596 }
597 
598 
599 
600 //___________________________________________________
602 {
603  TH1D *hPYTHIA13 = new TH1D("hPYTHIA13TeV_Ds", "", 40, 0., 40.);
604  Float_t val[40] = {
605  346381, 519143, 320488, 167944, 87484, 47325, 26932, 16376, 10058, 6527,
606  4347, 3041, 2112, 1521, 1069, 849, 609, 468, 359, 275,
607  226, 197, 152, 124, 111, 90, 63, 60, 52, 46,
608  39, 35, 26, 29, 16, 20, 10, 13, 18, 9
609  };
610  for (Int_t ibin=0; ibin<40; ++ibin) hPYTHIA13->SetBinContent(ibin+1, val[ibin]);
611 
612  return hPYTHIA13;
613 }
614 
615 
616 
617 //___________________________________________________
619 {
620  TH1D *hPYTHIA13 = new TH1D("hPYTHIA13_feeddownD0", "", 40, 0., 40.);
621  Float_t val[40] = {
622  1814030, 2939561, 2059721, 1197596, 684753, 397435, 237432, 146442, 93029, 60222,
623  40940, 28861, 19992, 14344, 10602, 7742, 5781, 4395, 3421, 2732,
624  2216, 1763, 1395, 1208, 888, 735, 599, 509, 418, 374,
625  302, 281, 242, 200, 174, 158, 133, 104, 98, 88
626  };
627  for (Int_t ibin=0; ibin<40; ++ibin) hPYTHIA13->SetBinContent(ibin+1, val[ibin]);
628 
629  return hPYTHIA13;
630 }
631 
632 
633 
634 //___________________________________________________
636 {
637  TH1D *hPYTHIA13 = new TH1D("hPYTHIA13_feeddownDplus", "", 40, 0., 40.);
638  Float_t val[40] = {
639  779389, 1268000, 899979, 527442, 302782, 177610, 105677, 66601, 42008, 27618,
640  18802, 13024, 8975, 6530, 4773, 3495, 2708, 2020, 1634, 1212,
641  1017, 809, 665, 507, 389, 355, 285, 250, 218, 180,
642  143, 137, 106, 85, 71, 64, 57, 43, 50, 47
643  };
644  for (Int_t ibin=0; ibin<40; ++ibin) hPYTHIA13->SetBinContent(ibin+1, val[ibin]);
645 
646  return hPYTHIA13;
647 }
648 
649 
650 
651 //___________________________________________________
653 {
654  TH1D *hPYTHIA13 = new TH1D("hPYTHIA13_feeddownDstar", "", 40, 0., 40.);
655  Float_t val[40] = {
656  665735, 1144214, 856647, 523704, 309061, 183569, 111994, 70187, 45266, 29601,
657  20079, 14042, 10338, 7184, 5234, 4002, 2935, 2301, 1788, 1349,
658  1088, 895, 748, 591, 480, 378, 302, 263, 216, 182,
659  172, 122, 133, 100, 92, 87, 65, 50, 42, 45
660  };
661  for (Int_t ibin=0; ibin<40; ++ibin) hPYTHIA13->SetBinContent(ibin+1, val[ibin]);
662 
663  return hPYTHIA13;
664 }
665 
666 
667 
668 //___________________________________________________
670 {
671  TH1D *hPYTHIA13 = new TH1D("hPYTHIA13_feeddownDs", "", 40, 0., 40.);
672  Float_t val[40] = {
673  377925, 686435, 518802, 318819, 189396, 113651, 68963, 43048, 27776, 18177,
674  12356, 8643, 6100, 4512, 3170, 2479, 1835, 1424, 1086, 825,
675  659, 490, 465, 361, 308, 246, 181, 170, 117, 111,
676  112, 91, 69, 66, 54, 57, 39, 33, 31, 25
677  };
678  for (Int_t ibin=0; ibin<40; ++ibin) hPYTHIA13->SetBinContent(ibin+1, val[ibin]);
679 
680  return hPYTHIA13;
681 }
Int_t fPtShape
Double_t fEtaMaxDau
Bool_t CountKpi(TClonesArray *array, Int_t nentries, Int_t &nPions, Int_t &nKaons, Int_t &nPionsInAcc, Int_t &nKaonsInAcc)
const Double_t ymax
Definition: AddTaskCFDStar.C:7
double Double_t
Definition: External.C:58
TH1D * LoadPYTHIA13TeV_promptD0()
TSystem * gSystem
Int_t fOptionYFiducial
TH1D * LoadFONLL13TeV_feeddownDstar()
void ComputeAcceptance()
Double_t fPtMinDau
Double_t massD
int Int_t
Definition: External.C:63
float Float_t
Definition: External.C:68
Definition: External.C:228
Definition: External.C:212
Int_t totTrials
TH1D * LoadPYTHIA13TeV_feeddownDs()
TH1D * LoadPYTHIA13TeV_feeddownDplus()
TH1D * LoadPYTHIA13TeV_feeddownD0()
Int_t fDebugLevel
TH1D * LoadPYTHIA13TeV_promptDs()
EPtShape
Int_t fDDecay
TH1D * LoadFONLL13TeV_promptD0()
TH1D * LoadFONLL13TeV_promptDstar()
TString fDecayTableFileName
Double_t fYMaxFidAccCut
TH1D * LoadFONLL13TeV_promptDplus()
TH1D * LoadFONLL13TeV_feeddownD()
TH1D * LoadPYTHIA13TeV_promptDstar()
bool Bool_t
Definition: External.C:53
TH1D * LoadPYTHIA13TeV_feeddownDstar()
Bool_t CountPKpi(TClonesArray *array, Int_t nentries, Int_t &nPions, Int_t &nKaons, Int_t &nProtons, Int_t &nPionsInAcc, Int_t &nKaonsInAcc, Int_t &nProtonsInAcc)
Bool_t IsInFiducialAcceptance(Double_t pt, Double_t y)
TH1D * LoadPYTHIA13TeV_promptDplus()