AliPhysics  608b256 (608b256)
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, Int_t &idLcResChan);
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 if(fDDecay==kDplusKKpi){
120  pdgCode=411;
121  nPionDau=1;
122  nKaonDau=2;
123  nProtonDau=0;
124  outFileName.Append("DplusKKpi_");
125  }else{
126  printf("ERROR: Wrong decay selected\n");
127  return;
128  }
129  if(fOptionYFiducial==kFixedY) outFileName.Append(Form("yfid%02d_",(Int_t)(fYMaxFidAccCut*10)));
130  else outFileName.Append("yfidPtDep_");
131  outFileName.Append(Form("etaDau%02d_",(Int_t)(fEtaMaxDau*10)));
132  outFileName.Append(Form("ptDau%d_",(Int_t)(fPtMinDau*1000)));
133  TDatabasePDG* db=TDatabasePDG::Instance();
134  Float_t massD=db->GetParticle(pdgCode)->Mass();
135  TClonesArray *array = new TClonesArray("TParticle",100);
136 
137  TH2D* hPtVsYGen=new TH2D("hPtVsYGen","",400,0.,40.,20.,-1.,1.);
138  hPtVsYGen->GetXaxis()->SetTitle("p_{T} (GeV/c)");
139  hPtVsYGen->GetYaxis()->SetTitle("y");
140  TH2D* hPtVsYGenLimAcc=new TH2D("hPtVsYGenLimAcc","",400,0.,40.,20.,-1.,1.);
141  hPtVsYGenLimAcc->GetXaxis()->SetTitle("p_{T} (GeV/c)");
142  hPtVsYGenLimAcc->GetYaxis()->SetTitle("y");
143  TH2D* hPtVsYGenAcc=new TH2D("hPtVsYGenAcc","",400,0.,40.,20.,-1.,1.);
144  hPtVsYGenAcc->GetXaxis()->SetTitle("p_{T} (GeV/c)");
145  hPtVsYGenAcc->GetYaxis()->SetTitle("y");
146 
147  // Histograms for Lc->pKpi split by resonance
148  TH2D* hPtVsYGenLimAccLcpKpi[4];
149  TH2D* hPtVsYGenAccLcpKpi[4];
150  TString lcChan[4]={"NonRes","L1520","Kstar","Delta"};
151  for(Int_t ich=0; ich<4; ich++){
152  hPtVsYGenLimAccLcpKpi[ich]=new TH2D(Form("hPtVsYGenLimAcc%s",lcChan[ich].Data()),"",400,0.,40.,20.,-1.,1.);
153  hPtVsYGenLimAccLcpKpi[ich]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
154  hPtVsYGenLimAccLcpKpi[ich]->GetYaxis()->SetTitle("y");
155  hPtVsYGenAccLcpKpi[ich]=new TH2D(Form("hPtVsYGenAcc%s",lcChan[ich].Data()),"",400,0.,40.,20.,-1.,1.);
156  hPtVsYGenAccLcpKpi[ich]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
157  hPtVsYGenAccLcpKpi[ich]->GetYaxis()->SetTitle("y");
158  }
159 
160 
161  TF1* funcPt=0x0;
162  TH1D* histPt=0x0;
163  if(fPtShape==kFONLL8TeV){
164  funcPt=new TF1("fFONLL","[0]*x/TMath::Power((1+TMath::Power(x/[1],[3])),[2])",0.,40.);
165  funcPt->SetParameters(0.518046,3.01138,3.38914,1.75899); // Prompt
166  outFileName.Append("FONLL8ptshape.root");
167  }else if(fPtShape==kFONLL8TeVfeeddown){
168  funcPt=new TF1("fFONLL","[0]*x/TMath::Power((1+TMath::Power(x/[1],[3])),[2])",0.,40.);
169  funcPt->SetParameters(0.398252, 3.9603, 3.915, 1.51853); // FeedDown
170  outFileName.Append("FONLL8ptshapeFeedDown.root");
171  }else if(fPtShape==kFONLL7TeV){
172  funcPt=new TF1("fFONLL","[0]*x/TMath::Power((1+TMath::Power(x/[1],[3])),[2])",0.,40.);
173  funcPt->SetParameters(0.322643,2.96275,2.30301,2.5);
174  outFileName.Append("FONLL7ptshape.root");
175  }else if(fPtShape==kFONLL5TeV){
176  funcPt=new TF1("fFONLL","[0]*x/TMath::Power((1+TMath::Power(x/[1],[3])),[2])",0.,40.);
177  funcPt->SetParameters(0.302879,2.9750,3.68139,1.68855);
178  outFileName.Append("FONLL5ptshape.root");
179  }else if(fPtShape==kPythia7TeV){
180  funcPt=new TF1("fFONLL","[0]*x/TMath::Power((1+TMath::Power(x/[1],[3])),[2])",0.,40.);
181  funcPt->SetParameters(0.322643,1.94635,1.40463,2.5);
182  outFileName.Append("PYTHIA7ptshape.root");
183  }else if(fPtShape==kFONLL13TeVprompt){
184  if(fDDecay==kDplusKpipi){
185  histPt = LoadFONLL13TeV_promptDplus();
186  outFileName.Append("promptDplus");
187  }else if (fDDecay==kDstarD0pi){
188  histPt = LoadFONLL13TeV_promptDstar();
189  outFileName.Append("promptDstar");
190  }else{
191  histPt = LoadFONLL13TeV_promptD0();
192  outFileName.Append("promptD0");
193  }
194  outFileName.Append("FONLL13ptshape.root");
195 
196  }else if (fPtShape==kFONLL13TeVfeeddown){
197  if (fDDecay==kDstarD0pi){
198  histPt = LoadFONLL13TeV_feeddownDstar();
199  outFileName.Append("feeddownDstar");
200  }else{
201  histPt = LoadFONLL13TeV_feeddownD();
202  outFileName.Append("feeddownD");
203  }
204  outFileName.Append("FONLL13ptshape.root");
205  }else if(fPtShape==kPythia13TeVprompt){
206  if(fDDecay==kDplusKpipi){
207  histPt = LoadPYTHIA13TeV_promptDplus();
208  outFileName.Append("promptDplus");
209  }else if (fDDecay==kDstarD0pi){
210  histPt = LoadPYTHIA13TeV_promptDstar();
211  outFileName.Append("promptDstar");
212  }else if (fDDecay==kDsKKpi){
213  histPt = LoadPYTHIA13TeV_promptDs();
214  outFileName.Append("promptDs");
215  }else{
216  histPt = LoadPYTHIA13TeV_promptD0();
217  outFileName.Append("promptD0");
218  }
219  outFileName.Append("PYTHIA13ptshape.root");
220  }else if (fPtShape==kPythia13TeVfeeddown){
221  if(fDDecay==kDplusKpipi){
223  outFileName.Append("feeddownDplus");
224  }else if (fDDecay==kDstarD0pi){
226  outFileName.Append("feeddownDstar");
227  }else if (fDDecay==kDsKKpi){
228  histPt = LoadPYTHIA13TeV_feeddownDs();
229  outFileName.Append("feeddownDs");
230  }else{
231  histPt = LoadPYTHIA13TeV_feeddownD0();
232  outFileName.Append("feeddownD0");
233  }
234  outFileName.Append("PYTHIA13ptshape.root");
235  }else{
236  funcPt=new TF1("fFlat","pol0",0.,40.);
237  funcPt->SetParameter(0,1.);
238  outFileName.Append("flatpt.root");
239  }
240 
241  if (funcPt) funcPt->SetNpx(10000);
242 
243  TRandom3* gener=new TRandom3(0);
244  TLorentzVector* vec=new TLorentzVector();
245 
246 
247  for(Int_t itry=0; itry<totTrials; itry++){
248  if(itry%10000==0) printf("Event %d\n",itry);
249  Float_t ptD = funcPt ? funcPt->GetRandom() : histPt->GetRandom();
250  Float_t phiD=gener->Rndm()*2*TMath::Pi();
251  Float_t yD=gener->Rndm()*2.-1.; // flat in -1<y<1
252  Float_t px=ptD*TMath::Cos(phiD);
253  Float_t py=ptD*TMath::Sin(phiD);
254  Float_t mt=TMath::Sqrt(massD*massD+ptD*ptD);
255  Float_t pz=mt*TMath::SinH(yD);
256  Float_t E=TMath::Sqrt(massD*massD+px*px+py*py+pz*pz);
257 
258  // TLorentzVector* vec=new TLorentzVector(px,py,pz,E);
259  vec->SetPxPyPzE(px,py,pz,E);
260  pdec->Decay(pdgCode,vec);
261  array->Clear();
262  Int_t nentries = pdec->ImportParticles(array);
263  TParticle* dmes=(TParticle*)array->At(0);
264  Int_t nDaughters=dmes->GetNDaughters();
265  if(fDDecay==kD0Kpi && nDaughters!=2) continue;
266  if(fDDecay==kLcK0Sp && nentries>6) continue;
267  Int_t nPionsInAcc=0;
268  Int_t nProtonsInAcc=0;
269  Int_t nKaonsInAcc=0;
270  Int_t nPions=0;
271  Int_t nProtons=0;
272  Int_t nKaons=0;
273  Int_t idLcResChan=0; // non resonant by default;
274  Bool_t isOk=CountPKpi(array,nentries,nPions,nKaons,nProtons,nPionsInAcc,nKaonsInAcc,nProtonsInAcc, idLcResChan);
275 
276  if(isOk){
277  if(nPions==nPionDau && nKaons==nKaonDau && nProtons==nProtonDau){
278  hPtVsYGen->Fill(ptD,yD);
279  if(TMath::Abs(yD)<0.5){
280  hPtVsYGenLimAcc->Fill(ptD,yD);
281  hPtVsYGenLimAccLcpKpi[idLcResChan]->Fill(ptD,yD);
282  }
283  if(IsInFiducialAcceptance(ptD,yD)){
284  if(nPionsInAcc==nPionDau && nKaonsInAcc==nKaonDau && nProtonsInAcc==nProtonDau){
285  hPtVsYGenAcc->Fill(ptD,yD);
286  hPtVsYGenAccLcpKpi[idLcResChan]->Fill(ptD,yD);
287  }
288  }
289  }
290  }
291  // delete vec;
292  }
293 
294  TH1D* hPtGenAcc=(TH1D*)hPtVsYGenAcc->ProjectionX("hPtGenAcc");
295  hPtGenAcc->GetYaxis()->SetTitle("Entries");
296  TH1D* hPtGenLimAcc=(TH1D*)hPtVsYGenLimAcc->ProjectionX("hPtGenLimAcc");
297  hPtGenLimAcc->GetYaxis()->SetTitle("Entries");
298  hPtGenAcc->Sumw2();
299  hPtGenLimAcc->Sumw2();
300  TH1D* hAccVsPt=(TH1D*)hPtGenAcc->Clone("hAccVsPt");
301  hAccVsPt->Divide(hPtGenAcc,hPtGenLimAcc,1,1,"B");
302  hAccVsPt->GetYaxis()->SetTitle("Acceptance");
303  hAccVsPt->SetStats(0);
304 
305  TCanvas* c2d=new TCanvas("c2d","Pt vs y",1200,600);
306  c2d->Divide(3,1);
307  c2d->cd(1);
308  hPtVsYGen->Draw("colz");
309  c2d->cd(2);
310  hPtVsYGenLimAcc->Draw("colz");
311  c2d->cd(3);
312  hPtVsYGenAcc->Draw("colz");
313 
314  TCanvas* c1d=new TCanvas("c1d","Acceptance",1200,600);
315  c1d->Divide(2,1);
316  c1d->cd(1);
317  hPtGenLimAcc->Draw("");
318  Double_t ymax=1.2*TMath::Max(hPtGenLimAcc->GetMaximum(),hPtGenAcc->GetMaximum());
319  hPtGenLimAcc->SetMaximum(ymax);
320  gPad->Update();
321  TPaveStats *st1=(TPaveStats*)hPtGenLimAcc->GetListOfFunctions()->FindObject("stats");
322  st1->SetY1NDC(0.71);
323  st1->SetY2NDC(0.9);
324  hPtGenAcc->SetLineColor(kRed+1);
325  hPtGenAcc->Draw("sames");
326  gPad->Update();
327  TPaveStats *st2=(TPaveStats*)hPtGenAcc->GetListOfFunctions()->FindObject("stats");
328  st2->SetY1NDC(0.51);
329  st2->SetY2NDC(0.7);
330  st2->SetTextColor(hPtGenAcc->GetLineColor());
331  gPad->Modified();
332  c1d->cd(2);
333  hAccVsPt->Draw();
334 
335 
336  TH1D* hPtGenAccLcpKpi[4]={0x0,0x0,0x0,0x0};
337  TH1D* hPtGenLimAccLcpKpi[4]={0x0,0x0,0x0,0x0};
338  TH1D* hAccVsPtLcpKpi[4]={0x0,0x0,0x0,0x0};
339  if(fDDecay==kLcpKpi){
340  Int_t nentr[4];
341  Int_t nSum=0;
342  Int_t nTot=hPtGenLimAcc->GetEntries();
343  printf("Lc decay modes fractions: ");
344  for(Int_t ich=0; ich<4; ich++){
345  hPtGenAccLcpKpi[ich]=(TH1D*)hPtVsYGenAccLcpKpi[ich]->ProjectionX(Form("hPtGenAcc%s",lcChan[ich].Data()));
346  hPtGenAccLcpKpi[ich]->GetYaxis()->SetTitle("Entries");
347  hPtGenAccLcpKpi[ich]->Sumw2();
348  hPtGenLimAccLcpKpi[ich]=(TH1D*)hPtVsYGenLimAccLcpKpi[ich]->ProjectionX(Form("hPtGenLimAcc%s",lcChan[ich].Data()));
349  hPtGenLimAccLcpKpi[ich]->GetYaxis()->SetTitle("Entries");
350  hPtGenLimAccLcpKpi[ich]->Sumw2();
351  hAccVsPtLcpKpi[ich]=(TH1D*)hPtGenAccLcpKpi[ich]->Clone(Form("hAccVsPt%s",lcChan[ich].Data()));
352  hAccVsPtLcpKpi[ich]->Divide(hPtGenAccLcpKpi[ich],hPtGenLimAccLcpKpi[ich],1.,1.,"B");
353  hAccVsPtLcpKpi[ich]->GetYaxis()->SetTitle("Acceptance");
354  hAccVsPtLcpKpi[ich]->SetStats(0);
355  nentr[ich]=hPtGenLimAccLcpKpi[ich]->GetEntries();
356  nSum+=nentr[ich];
357  printf(" %s = %.3f, ",lcChan[ich].Data(),(Float_t)nentr[ich]/(Float_t)nTot);
358  }
359  printf("Check = %.3f\n",(Float_t)nSum/(Float_t)nTot);
360  }
361 
362  TFile* outfil=new TFile(outFileName.Data(),"recreate");
363  hPtVsYGen->Write();
364  hPtVsYGenLimAcc->Write();
365  hPtVsYGenAcc->Write();
366  hPtGenLimAcc->Write();
367  hPtGenAcc->Write();
368  hAccVsPt->Write();
369  if(fDDecay==kLcpKpi){
370  for(Int_t ich=0; ich<4; ich++){
371  if(hPtGenAccLcpKpi[ich]) hPtGenAccLcpKpi[ich]->Write();
372  if(hPtGenLimAccLcpKpi[ich]) hPtGenLimAccLcpKpi[ich]->Write();
373  if(hAccVsPtLcpKpi[ich]) hAccVsPtLcpKpi[ich]->Write();
374  }
375  }
376  outfil->Close();
377 
378 }
379 
380 //___________________________________________________
382  // check fiducial acceptance
383 
385  if(TMath::Abs(y) > fYMaxFidAccCut) return kFALSE;
386  else return kTRUE;
387  }
388 
389  if(pt > 5.) {
390  if (TMath::Abs(y) > 0.8) return kFALSE;
391  } else {
392  Double_t maxFiducialY = -0.2/15*pt*pt+1.9/15*pt+0.5;
393  Double_t minFiducialY = 0.2/15*pt*pt-1.9/15*pt-0.5;
394  if (y < minFiducialY || y > maxFiducialY) return kFALSE;
395  }
396 
397  return kTRUE;
398 }
399 //___________________________________________________
400 Bool_t CountKpi(TClonesArray *array, Int_t nentries, Int_t &nPions, Int_t &nKaons, Int_t &nPionsInAcc, Int_t &nKaonsInAcc){
401  // count K and pi in Acc
402 
403  TParticle* dmes=(TParticle*)array->At(0);
404  Double_t sumPx=0;
405  Double_t sumPy=0;
406  Double_t sumPz=0;
407  for(int j=0; j<nentries; j++){
408  TParticle * o = (TParticle*)array->At(j);
409  Int_t pdgdau=TMath::Abs(o->GetPdgCode());
410  if(fDebugLevel>0) printf("%d ",pdgdau);
411  if(pdgdau==130) {
412  if(fDebugLevel>0) printf("K0 dacaying into K0L\n");
413  return kFALSE;
414  }
415  Float_t ptdau=TMath::Sqrt(o->Px()*o->Px()+o->Py()*o->Py());
416  Float_t etadau=o->Eta();
417  if(pdgdau==211){
418  nPions++;
419  sumPx+=o->Px();
420  sumPy+=o->Py();
421  sumPz+=o->Pz();
422  }
423  if(pdgdau==321){
424  nKaons++;
425  sumPx+=o->Px();
426  sumPy+=o->Py();
427  sumPz+=o->Pz();
428  }
429  if(TMath::Abs(etadau)<fEtaMaxDau && ptdau>fPtMinDau){
430  if(pdgdau==211) nPionsInAcc++;
431  if(pdgdau==321) nKaonsInAcc++;
432  }
433  }
434  if(fDebugLevel>0) printf("\n");
435  if(TMath::Abs(sumPx-dmes->Px())>0.001 ||
436  TMath::Abs(sumPy-dmes->Py())>0.001 ||
437  TMath::Abs(sumPz-dmes->Pz())>0.001){
438  printf("Momentum conservation violation\n");
439  return kFALSE;
440  }
441  return kTRUE;
442 }
443 
444 
445 //___________________________________________________
446 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, Int_t &idLcResChan){
447  // count K and pi in Acc
448 
449  TParticle* dmes=(TParticle*)array->At(0);
450  Double_t sumPx=0;
451  Double_t sumPy=0;
452  Double_t sumPz=0;
453 
454  for(int j=0; j<nentries; j++){
455  TParticle * o = (TParticle*)array->At(j);
456  Int_t pdgdau=TMath::Abs(o->GetPdgCode());
457  if(fDebugLevel>0) printf("%d ",pdgdau);
458  if(pdgdau==130) {
459  if(fDebugLevel>0) printf("K0 dacaying into K0L\n");
460  return kFALSE;
461  }
462  Float_t ptdau=TMath::Sqrt(o->Px()*o->Px()+o->Py()*o->Py());
463  Float_t etadau=o->Eta();
464  if(pdgdau==211){
465  nPions++;
466  sumPx+=o->Px();
467  sumPy+=o->Py();
468  sumPz+=o->Pz();
469  }
470  if(pdgdau==321){
471  nKaons++;
472  sumPx+=o->Px();
473  sumPy+=o->Py();
474  sumPz+=o->Pz();
475  }
476  if(pdgdau==2212){
477  nProtons++;
478  sumPx+=o->Px();
479  sumPy+=o->Py();
480  sumPz+=o->Pz();
481  }
482  if(pdgdau==313) idLcResChan=2; //K*0
483  else if(pdgdau==2224) idLcResChan=3; //Delta++
484  else if(pdgdau==3124) idLcResChan=1; //Lambda1520
485  if(TMath::Abs(etadau)<fEtaMaxDau && ptdau>fPtMinDau){
486  if(pdgdau==211) nPionsInAcc++;
487  if(pdgdau==321) nKaonsInAcc++;
488  if(pdgdau==2212) nProtonsInAcc++;
489  }
490  }
491  if(fDebugLevel>0) printf("\n");
492  if(TMath::Abs(sumPx-dmes->Px())>0.001 ||
493  TMath::Abs(sumPy-dmes->Py())>0.001 ||
494  TMath::Abs(sumPz-dmes->Pz())>0.001){
495  printf("Momentum conservation violation\n");
496  return kFALSE;
497  }
498  return kTRUE;
499 }
500 
501 
502 
503 //___________________________________________________
505 {
506  TH1D *hFONLL13 = new TH1D("hFONLL13TeV_D0", "", 80, 0., 40.);
507  Float_t val[80] = {
508  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,
509  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,
510  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,
511  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,
512  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,
513  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,
514  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,
515  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
516  };
517  for (Int_t ibin=0; ibin<80; ++ibin) hFONLL13->SetBinContent(ibin+1, val[ibin]);
518 
519  return hFONLL13;
520 }
521 
522 
523 
524 //___________________________________________________
526 {
527  TH1D *hFONLL13 = new TH1D("hFONLL13TeV_Dplus", "", 80, 0., 40.);
528  Float_t val[80] = {
529  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,
530  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,
531  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,
532  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,
533  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,
534  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,
535  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,
536  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
537  };
538  for (Int_t ibin=0; ibin<80; ++ibin) hFONLL13->SetBinContent(ibin+1, val[ibin]);
539 
540  return hFONLL13;
541 }
542 
543 
544 
545 //___________________________________________________
547 {
548  TH1D *hFONLL13 = new TH1D("hFONLL13TeV_Dstar", "", 80, 0., 40.);
549  Float_t val[80] = {
550  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,
551  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,
552  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,
553  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,
554  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,
555  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,
556  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,
557  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
558  };
559  for (Int_t ibin=0; ibin<80; ++ibin) hFONLL13->SetBinContent(ibin+1, val[ibin]);
560 
561  return hFONLL13;
562 }
563 
564 
565 
566 //___________________________________________________
568 { // B->D with B.R=1
569  TH1D *hFONLL13 = new TH1D("FONLL13TeV_feeddownD", "", 80, 0., 40.);
570  Float_t val[80] = {
571  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,
572  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,
573  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,
574  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,
575  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,
576  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,
577  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,
578  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
579  };
580  for (Int_t ibin=0; ibin<80; ++ibin) hFONLL13->SetBinContent(ibin+1, val[ibin]);
581 
582  return hFONLL13;
583 }
584 
585 
586 
587 //___________________________________________________
589 { // B->D* with B.R=1
590  TH1D *hFONLL13 = new TH1D("FONLL13TeV_feeddownDstar", "", 80, 0., 40.);
591  Float_t val[80] = {
592  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,
593  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,
594  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,
595  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,
596  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,
597  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,
598  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,
599  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
600  };
601  for (Int_t ibin=0; ibin<80; ++ibin) hFONLL13->SetBinContent(ibin+1, val[ibin]);
602 
603  return hFONLL13;
604 }
605 
606 
607 
608 //___________________________________________________
610 {
611  TH1D *hPYTHIA13 = new TH1D("hPYTHIA13TeV_D0", "", 40, 0., 40.);
612  Float_t val[40] = {
613  2617743, 3763836, 2235903, 1140807, 587707, 317881, 180756, 107762, 67299,42882,
614  28550, 19690, 13817, 9754, 7205, 5299, 4100, 3177, 2386,1822,
615  1484, 1248, 933, 798, 639, 534, 435, 405, 304,260,
616  210, 196, 161, 143, 139, 99, 89, 80, 66,63
617  };
618  for (Int_t ibin=0; ibin<40; ++ibin) hPYTHIA13->SetBinContent(ibin+1, val[ibin]);
619 
620  return hPYTHIA13;
621 }
622 
623 
624 
625 //___________________________________________________
627 {
628  TH1D *hPYTHIA13 = new TH1D("hPYTHIA13TeV_Dplus", "", 40, 0., 40.);
629  Float_t val[40] = {
630  1192379, 1736283, 1047485, 540499, 281523, 152153, 87150, 52200, 32709, 20991,
631  13986, 9611, 6832, 4705, 3484, 2650, 2001, 1518, 1197, 924,
632  735, 576, 497, 402, 292, 261, 214, 173, 169, 150,
633  106, 102, 95, 70, 67, 47, 53, 46, 35, 35
634  };
635  for (Int_t ibin=0; ibin<40; ++ibin) hPYTHIA13->SetBinContent(ibin+1, val[ibin]);
636 
637  return hPYTHIA13;
638 }
639 
640 
641 
642 //___________________________________________________
644 {
645  TH1D *hPYTHIA13 = new TH1D("hPYTHIA13TeV_Dstar", "", 40, 0., 40.);
646  Float_t val[40] = {
647  922985, 1419749, 904587, 484407, 257862, 142386, 82177, 49839, 31103, 20422,
648  13493, 9272, 6549, 4746, 3344, 2575, 1959, 1493, 1161, 935,
649  672, 585, 456, 415, 328, 260, 199, 180, 161, 134,
650  109, 97, 103, 72, 63, 56, 40, 46, 38, 27
651  };
652  for (Int_t ibin=0; ibin<40; ++ibin) hPYTHIA13->SetBinContent(ibin+1, val[ibin]);
653 
654  return hPYTHIA13;
655 }
656 
657 
658 
659 //___________________________________________________
661 {
662  TH1D *hPYTHIA13 = new TH1D("hPYTHIA13TeV_Ds", "", 40, 0., 40.);
663  Float_t val[40] = {
664  346381, 519143, 320488, 167944, 87484, 47325, 26932, 16376, 10058, 6527,
665  4347, 3041, 2112, 1521, 1069, 849, 609, 468, 359, 275,
666  226, 197, 152, 124, 111, 90, 63, 60, 52, 46,
667  39, 35, 26, 29, 16, 20, 10, 13, 18, 9
668  };
669  for (Int_t ibin=0; ibin<40; ++ibin) hPYTHIA13->SetBinContent(ibin+1, val[ibin]);
670 
671  return hPYTHIA13;
672 }
673 
674 
675 
676 //___________________________________________________
678 {
679  TH1D *hPYTHIA13 = new TH1D("hPYTHIA13_feeddownD0", "", 40, 0., 40.);
680  Float_t val[40] = {
681  1814030, 2939561, 2059721, 1197596, 684753, 397435, 237432, 146442, 93029, 60222,
682  40940, 28861, 19992, 14344, 10602, 7742, 5781, 4395, 3421, 2732,
683  2216, 1763, 1395, 1208, 888, 735, 599, 509, 418, 374,
684  302, 281, 242, 200, 174, 158, 133, 104, 98, 88
685  };
686  for (Int_t ibin=0; ibin<40; ++ibin) hPYTHIA13->SetBinContent(ibin+1, val[ibin]);
687 
688  return hPYTHIA13;
689 }
690 
691 
692 
693 //___________________________________________________
695 {
696  TH1D *hPYTHIA13 = new TH1D("hPYTHIA13_feeddownDplus", "", 40, 0., 40.);
697  Float_t val[40] = {
698  779389, 1268000, 899979, 527442, 302782, 177610, 105677, 66601, 42008, 27618,
699  18802, 13024, 8975, 6530, 4773, 3495, 2708, 2020, 1634, 1212,
700  1017, 809, 665, 507, 389, 355, 285, 250, 218, 180,
701  143, 137, 106, 85, 71, 64, 57, 43, 50, 47
702  };
703  for (Int_t ibin=0; ibin<40; ++ibin) hPYTHIA13->SetBinContent(ibin+1, val[ibin]);
704 
705  return hPYTHIA13;
706 }
707 
708 
709 
710 //___________________________________________________
712 {
713  TH1D *hPYTHIA13 = new TH1D("hPYTHIA13_feeddownDstar", "", 40, 0., 40.);
714  Float_t val[40] = {
715  665735, 1144214, 856647, 523704, 309061, 183569, 111994, 70187, 45266, 29601,
716  20079, 14042, 10338, 7184, 5234, 4002, 2935, 2301, 1788, 1349,
717  1088, 895, 748, 591, 480, 378, 302, 263, 216, 182,
718  172, 122, 133, 100, 92, 87, 65, 50, 42, 45
719  };
720  for (Int_t ibin=0; ibin<40; ++ibin) hPYTHIA13->SetBinContent(ibin+1, val[ibin]);
721 
722  return hPYTHIA13;
723 }
724 
725 
726 
727 //___________________________________________________
729 {
730  TH1D *hPYTHIA13 = new TH1D("hPYTHIA13_feeddownDs", "", 40, 0., 40.);
731  Float_t val[40] = {
732  377925, 686435, 518802, 318819, 189396, 113651, 68963, 43048, 27776, 18177,
733  12356, 8643, 6100, 4512, 3170, 2479, 1835, 1424, 1086, 825,
734  659, 490, 465, 361, 308, 246, 181, 170, 117, 111,
735  112, 91, 69, 66, 54, 57, 39, 33, 31, 25
736  };
737  for (Int_t ibin=0; ibin<40; ++ibin) hPYTHIA13->SetBinContent(ibin+1, val[ibin]);
738 
739  return hPYTHIA13;
740 }
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
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
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, Int_t &idLcResChan)
Bool_t IsInFiducialAcceptance(Double_t pt, Double_t y)
TH1D * LoadPYTHIA13TeV_promptDplus()