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