AliPhysics  6bc8652 (6bc8652)
 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/../src/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 
42  // main function
43 
44  gSystem->Load("liblhapdf.so"); // Parton density functions
45  gSystem->Load("libEGPythia6.so"); // TGenerator interface
46  gSystem->Load("libpythia6.so"); // Pythia
47  gSystem->Load("libAliPythia6"); // ALICE specific implementations
48 
49  TPythia6Decayer* pdec=TPythia6Decayer::Instance();
50  if(fDecayTableFileName.CompareTo("")!=0){
51  if(fDecayTableFileName.Contains("ALICE_PHYSICS")){
52  gSystem->Exec(Form("cp %s .",fDecayTableFileName.Data()));
53  fDecayTableFileName.ReplaceAll("$ALICE_PHYSICS/../src/PWGHF/vertexingHF/macros/","./");
54  }
55  pdec->SetDecayTableFile(fDecayTableFileName.Data());
56  pdec->ReadDecayTable();
57  }
58  pdec->Init();
59 
60  Int_t pdgCode=0;
61  Int_t nPionDau=-1;
62  Int_t nProtonDau=-1;
63  Int_t nKaonDau=-1;
64  TString outFileName="Acceptance_Toy_";
65  if(fDDecay==kD0Kpi){
66  pdgCode=421;
67  nPionDau=1;
68  nKaonDau=1;
69  nProtonDau=0;
70  outFileName.Append("D0Kpi_");
71  }else if(fDDecay==kDplusKpipi){
72  pdgCode=411;
73  nPionDau=2;
74  nKaonDau=1;
75  nProtonDau=0;
76  outFileName.Append("DplusKpipi_");
77  }else if(fDDecay==kDstarD0pi){
78  pdgCode=413;
79  nPionDau=2;
80  nKaonDau=1;
81  nProtonDau=0;
82  outFileName.Append("DStarD0pi_");
83  }else if(fDDecay==kDsKKpi){
84  pdgCode=431;
85  nPionDau=1;
86  nKaonDau=2;
87  nProtonDau=0;
88  outFileName.Append("DsKKpi_");
89  }else if(fDDecay==kLcpKpi){
90  pdgCode=4122;
91  nPionDau=1;
92  nKaonDau=1;
93  nProtonDau=1;
94  outFileName.Append("LcpKpi_");
95  }else if(fDDecay==kLcK0Sp){
96  pdgCode=4122;
97  nPionDau=2;
98  nKaonDau=0;
99  nProtonDau=1;
100  outFileName.Append("LcK0Sp_");
101  }else{
102  printf("ERROR: Wrong decay selected\n");
103  return;
104  }
105  if(fOptionYFiducial==kFixedY) outFileName.Append(Form("yfid%02d_",(Int_t)(fYMaxFidAccCut*10)));
106  else outFileName.Append("yfidPtDep_");
107  outFileName.Append(Form("etaDau%02d_",(Int_t)(fEtaMaxDau*10)));
108  outFileName.Append(Form("ptDau%d_",(Int_t)(fPtMinDau*1000)));
109  TDatabasePDG* db=TDatabasePDG::Instance();
110  Float_t massD=db->GetParticle(pdgCode)->Mass();
111  TClonesArray *array = new TClonesArray("TParticle",100);
112 
113  TH2D* hPtVsYGen=new TH2D("hPtVsYGen","",400,0.,40.,20.,-1.,1.);
114  hPtVsYGen->GetXaxis()->SetTitle("p_{T} (GeV/c)");
115  hPtVsYGen->GetYaxis()->SetTitle("y");
116  TH2D* hPtVsYGenLimAcc=new TH2D("hPtVsYGenLimAcc","",400,0.,40.,20.,-1.,1.);
117  hPtVsYGenLimAcc->GetXaxis()->SetTitle("p_{T} (GeV/c)");
118  hPtVsYGenLimAcc->GetYaxis()->SetTitle("y");
119  TH2D* hPtVsYGenAcc=new TH2D("hPtVsYGenAcc","",400,0.,40.,20.,-1.,1.);
120  hPtVsYGenAcc->GetXaxis()->SetTitle("p_{T} (GeV/c)");
121  hPtVsYGenAcc->GetYaxis()->SetTitle("y");
122 
123  TF1* funcPt=0x0;
124  if(fPtShape==kFONLL8TeV){
125  funcPt=new TF1("fFONLL","[0]*x/TMath::Power((1+TMath::Power(x/[1],[3])),[2])",0.,40.);
126  funcPt->SetParameters(0.518046,3.01138,3.38914,1.75899); // Prompt
127  outFileName.Append("FONLL8ptshape.root");
128  }else if(fPtShape==kFONLL8TeVfeeddown){
129  funcPt=new TF1("fFONLL","[0]*x/TMath::Power((1+TMath::Power(x/[1],[3])),[2])",0.,40.);
130  funcPt->SetParameters(0.398252, 3.9603, 3.915, 1.51853); // FeedDown
131  outFileName.Append("FONLL8ptshapeFeedDown.root");
132  }else if(fPtShape==kFONLL7TeV){
133  funcPt=new TF1("fFONLL","[0]*x/TMath::Power((1+TMath::Power(x/[1],[3])),[2])",0.,40.);
134  funcPt->SetParameters(0.322643,2.96275,2.30301,2.5);
135  outFileName.Append("FONLL7ptshape.root");
136  }else if(fPtShape==kFONLL5TeV){
137  funcPt=new TF1("fFONLL","[0]*x/TMath::Power((1+TMath::Power(x/[1],[3])),[2])",0.,40.);
138  funcPt->SetParameters(0.302879,2.9750,3.68139,1.68855);
139  outFileName.Append("FONLL5ptshape.root");
140  }else if(fPtShape==kPythia7TeV){
141  funcPt=new TF1("fFONLL","[0]*x/TMath::Power((1+TMath::Power(x/[1],[3])),[2])",0.,40.);
142  funcPt->SetParameters(0.322643,1.94635,1.40463,2.5);
143  outFileName.Append("PYTHIA7ptshape.root");
144  }else{
145  funcPt=new TF1("fFlat","pol0",0.,40.);
146  funcPt->SetParameter(0,1.);
147  outFileName.Append("flatpt.root");
148  }
149  TRandom3* gener=new TRandom3(0);
150 
151  for(Int_t itry=0; itry<totTrials; itry++){
152  if(itry%10000==0) printf("Event %d\n",itry);
153  Float_t ptD=funcPt->GetRandom();
154  Float_t phiD=gener->Rndm()*2*TMath::Pi();
155  Float_t yD=gener->Rndm()*2.-1.; // flat in -1<y<1
156  Float_t px=ptD*TMath::Cos(phiD);
157  Float_t py=ptD*TMath::Sin(phiD);
158  Float_t mt=TMath::Sqrt(massD*massD+ptD*ptD);
159  Float_t pz=mt*TMath::SinH(yD);
160  Float_t E=TMath::Sqrt(massD*massD+px*px+py*py+pz*pz);
161 
162  TLorentzVector* vec=new TLorentzVector(px,py,pz,E);
163  pdec->Decay(pdgCode,vec);
164  array->Clear();
165  Int_t nentries = pdec->ImportParticles(array);
166  TParticle* dmes=(TParticle*)array->At(0);
167  Int_t nDaughters=dmes->GetNDaughters();
168  if(fDDecay==kD0Kpi && nDaughters!=2) return;
169  if(fDDecay==kLcK0Sp && nentries>6) continue;
170  Int_t nPionsInAcc=0;
171  Int_t nProtonsInAcc=0;
172  Int_t nKaonsInAcc=0;
173  Int_t nPions=0;
174  Int_t nProtons=0;
175  Int_t nKaons=0;
176  Bool_t isOk=CountPKpi(array,nentries,nPions,nKaons,nProtons,nPionsInAcc,nKaonsInAcc,nProtonsInAcc);
177 
178  if(isOk){
179  if(nPions==nPionDau && nKaons==nKaonDau && nProtons==nProtonDau){
180  hPtVsYGen->Fill(ptD,yD);
181  if(TMath::Abs(yD)<0.5){
182  hPtVsYGenLimAcc->Fill(ptD,yD);
183  }
184  if(IsInFiducialAcceptance(ptD,yD)){
185  if(nPionsInAcc==nPionDau && nKaonsInAcc==nKaonDau && nProtonsInAcc==nProtonDau){
186  hPtVsYGenAcc->Fill(ptD,yD);
187  }
188  }
189  }
190  }
191  delete vec;
192  }
193 
194  TH1D* hPtGenAcc=(TH1D*)hPtVsYGenAcc->ProjectionX("hPtGenAcc");
195  hPtGenAcc->GetYaxis()->SetTitle("Entries");
196  TH1D* hPtGenLimAcc=(TH1D*)hPtVsYGenLimAcc->ProjectionX("hPtGenLimAcc");
197  hPtGenLimAcc->GetYaxis()->SetTitle("Entries");
198  hPtGenAcc->Sumw2();
199  hPtGenLimAcc->Sumw2();
200  TH1D* hAccVsPt=(TH1D*)hPtGenAcc->Clone("hAccVsPt");
201  hAccVsPt->Divide(hPtGenAcc,hPtGenLimAcc,1,1,"B");
202  hAccVsPt->GetYaxis()->SetTitle("Acceptance");
203  hAccVsPt->SetStats(0);
204 
205  TCanvas* c2d=new TCanvas("c2d","Pt vs y",1200,600);
206  c2d->Divide(3,1);
207  c2d->cd(1);
208  hPtVsYGen->Draw("colz");
209  c2d->cd(2);
210  hPtVsYGenLimAcc->Draw("colz");
211  c2d->cd(3);
212  hPtVsYGenAcc->Draw("colz");
213 
214  TCanvas* c1d=new TCanvas("c1d","Acceptance",1200,600);
215  c1d->Divide(2,1);
216  c1d->cd(1);
217  hPtGenLimAcc->Draw("");
218  Double_t ymax=1.2*TMath::Max(hPtGenLimAcc->GetMaximum(),hPtGenAcc->GetMaximum());
219  hPtGenLimAcc->SetMaximum(ymax);
220  gPad->Update();
221  TPaveStats *st1=(TPaveStats*)hPtGenLimAcc->GetListOfFunctions()->FindObject("stats");
222  st1->SetY1NDC(0.71);
223  st1->SetY2NDC(0.9);
224  hPtGenAcc->SetLineColor(kRed+1);
225  hPtGenAcc->Draw("sames");
226  gPad->Update();
227  TPaveStats *st2=(TPaveStats*)hPtGenAcc->GetListOfFunctions()->FindObject("stats");
228  st2->SetY1NDC(0.51);
229  st2->SetY2NDC(0.7);
230  st2->SetTextColor(hPtGenAcc->GetLineColor());
231  gPad->Modified();
232  c1d->cd(2);
233  hAccVsPt->Draw();
234 
235 
236  TFile* outfil=new TFile(outFileName.Data(),"recreate");
237  hPtVsYGen->Write();
238  hPtVsYGenLimAcc->Write();
239  hPtVsYGenAcc->Write();
240  hPtGenLimAcc->Write();
241  hPtGenAcc->Write();
242  hAccVsPt->Write();
243  outfil->Close();
244 
245 }
246 
247 //___________________________________________________
249  // check fiducial acceptance
250 
252  if(TMath::Abs(y) > fYMaxFidAccCut) return kFALSE;
253  else return kTRUE;
254  }
255 
256  if(pt > 5.) {
257  if (TMath::Abs(y) > 0.8) return kFALSE;
258  } else {
259  Double_t maxFiducialY = -0.2/15*pt*pt+1.9/15*pt+0.5;
260  Double_t minFiducialY = 0.2/15*pt*pt-1.9/15*pt-0.5;
261  if (y < minFiducialY || y > maxFiducialY) return kFALSE;
262  }
263 
264  return kTRUE;
265 }
266 //___________________________________________________
267 Bool_t CountKpi(TClonesArray *array, Int_t nentries, Int_t &nPions, Int_t &nKaons, Int_t &nPionsInAcc, Int_t &nKaonsInAcc){
268  // count K and pi in Acc
269 
270  TParticle* dmes=(TParticle*)array->At(0);
271  Double_t sumPx=0;
272  Double_t sumPy=0;
273  Double_t sumPz=0;
274  for(int j=0; j<nentries; j++){
275  TParticle * o = (TParticle*)array->At(j);
276  Int_t pdgdau=TMath::Abs(o->GetPdgCode());
277  if(fDebugLevel>0) printf("%d ",pdgdau);
278  if(pdgdau==130) {
279  if(fDebugLevel>0) printf("K0 dacaying into K0L\n");
280  return kFALSE;
281  }
282  Float_t ptdau=TMath::Sqrt(o->Px()*o->Px()+o->Py()*o->Py());
283  Float_t etadau=o->Eta();
284  if(pdgdau==211){
285  nPions++;
286  sumPx+=o->Px();
287  sumPy+=o->Py();
288  sumPz+=o->Pz();
289  }
290  if(pdgdau==321){
291  nKaons++;
292  sumPx+=o->Px();
293  sumPy+=o->Py();
294  sumPz+=o->Pz();
295  }
296  if(TMath::Abs(etadau)<fEtaMaxDau && ptdau>fPtMinDau){
297  if(pdgdau==211) nPionsInAcc++;
298  if(pdgdau==321) nKaonsInAcc++;
299  }
300  }
301  if(fDebugLevel>0) printf("\n");
302  if(TMath::Abs(sumPx-dmes->Px())>0.001 ||
303  TMath::Abs(sumPy-dmes->Py())>0.001 ||
304  TMath::Abs(sumPz-dmes->Pz())>0.001){
305  printf("Momentum conservation violation\n");
306  return kFALSE;
307  }
308  return kTRUE;
309 }
310 
311 
312 //___________________________________________________
313 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){
314  // count K and pi in Acc
315 
316  TParticle* dmes=(TParticle*)array->At(0);
317  Double_t sumPx=0;
318  Double_t sumPy=0;
319  Double_t sumPz=0;
320 
321  for(int j=0; j<nentries; j++){
322  TParticle * o = (TParticle*)array->At(j);
323  Int_t pdgdau=TMath::Abs(o->GetPdgCode());
324  if(fDebugLevel>0) printf("%d ",pdgdau);
325  if(pdgdau==130) {
326  if(fDebugLevel>0) printf("K0 dacaying into K0L\n");
327  return kFALSE;
328  }
329  Float_t ptdau=TMath::Sqrt(o->Px()*o->Px()+o->Py()*o->Py());
330  Float_t etadau=o->Eta();
331  if(pdgdau==211){
332  nPions++;
333  sumPx+=o->Px();
334  sumPy+=o->Py();
335  sumPz+=o->Pz();
336  }
337  if(pdgdau==321){
338  nKaons++;
339  sumPx+=o->Px();
340  sumPy+=o->Py();
341  sumPz+=o->Pz();
342  }
343  if(pdgdau==2212){
344  nProtons++;
345  sumPx+=o->Px();
346  sumPy+=o->Py();
347  sumPz+=o->Pz();
348  }
349  if(TMath::Abs(etadau)<fEtaMaxDau && ptdau>fPtMinDau){
350  if(pdgdau==211) nPionsInAcc++;
351  if(pdgdau==321) nKaonsInAcc++;
352  if(pdgdau==2212) nProtonsInAcc++;
353  }
354  }
355  if(fDebugLevel>0) printf("\n");
356  if(TMath::Abs(sumPx-dmes->Px())>0.001 ||
357  TMath::Abs(sumPy-dmes->Py())>0.001 ||
358  TMath::Abs(sumPz-dmes->Pz())>0.001){
359  printf("Momentum conservation violation\n");
360  return kFALSE;
361  }
362  return kTRUE;
363 }
364 
365 
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
TSystem * gSystem
Int_t fOptionYFiducial
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
Int_t fDebugLevel
EPtShape
Int_t fDDecay
TString fDecayTableFileName
Double_t fYMaxFidAccCut
bool Bool_t
Definition: External.C:53
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)