AliPhysics  9fe175b (9fe175b)
 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
27 Double_t fPtMinDau=0.1;
28 Double_t fEtaMaxDau=0.9;
30 Double_t fYMaxFidAccCut=0.8;
32 TString fDecayTableFileName="$ALICE_PHYSICS/../src/PWGHF/vertexingHF/macros/decaytable_acc.dat";
33 Int_t fDebugLevel=0;
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);
38 Bool_t IsInFiducialAcceptance(Double_t pt, Double_t y);
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==kFONLL7TeV){
125  funcPt=new TF1("fFONLL","[0]*x/TMath::Power((1+TMath::Power(x/[1],[3])),[2])",0.,40.);
126  funcPt->SetParameters(0.322643,2.96275,2.30301,2.5);
127  outFileName.Append("FONLL7ptshape.root");
128  }else if(fPtShape==kFONLL5TeV){
129  funcPt=new TF1("fFONLL","[0]*x/TMath::Power((1+TMath::Power(x/[1],[3])),[2])",0.,40.);
130  funcPt->SetParameters(0.302879,2.9750,3.68139,1.68855);
131  outFileName.Append("FONLL5ptshape.root");
132  }else if(fPtShape==kPythia7TeV){
133  funcPt=new TF1("fFONLL","[0]*x/TMath::Power((1+TMath::Power(x/[1],[3])),[2])",0.,40.);
134  funcPt->SetParameters(0.322643,1.94635,1.40463,2.5);
135  outFileName.Append("PYTHIA7ptshape.root");
136  }else{
137  funcPt=new TF1("fFlat","pol0",0.,40.);
138  funcPt->SetParameter(0,1.);
139  outFileName.Append("flatpt.root");
140  }
141  TRandom3* gener=new TRandom3(0);
142 
143  for(Int_t itry=0; itry<totTrials; itry++){
144  if(itry%10000==0) printf("Event %d\n",itry);
145  Float_t ptD=funcPt->GetRandom();
146  Float_t phiD=gener->Rndm()*2*TMath::Pi();
147  Float_t yD=gener->Rndm()*2.-1.; // flat in -1<y<1
148  Float_t px=ptD*TMath::Cos(phiD);
149  Float_t py=ptD*TMath::Sin(phiD);
150  Float_t mt=TMath::Sqrt(massD*massD+ptD*ptD);
151  Float_t pz=mt*TMath::SinH(yD);
152  Float_t E=TMath::Sqrt(massD*massD+px*px+py*py+pz*pz);
153 
154  TLorentzVector* vec=new TLorentzVector(px,py,pz,E);
155  pdec->Decay(pdgCode,vec);
156  array->Clear();
157  Int_t nentries = pdec->ImportParticles(array);
158  TParticle* dmes=(TParticle*)array->At(0);
159  Int_t nDaughters=dmes->GetNDaughters();
160  if(fDDecay==kD0Kpi && nDaughters!=2) return;
161  if(fDDecay==kLcK0Sp && nentries>6) continue;
162  Int_t nPionsInAcc=0;
163  Int_t nProtonsInAcc=0;
164  Int_t nKaonsInAcc=0;
165  Int_t nPions=0;
166  Int_t nProtons=0;
167  Int_t nKaons=0;
168  Bool_t isOk=CountPKpi(array,nentries,nPions,nKaons,nProtons,nPionsInAcc,nKaonsInAcc,nProtonsInAcc);
169 
170  if(isOk){
171  if(nPions==nPionDau && nKaons==nKaonDau && nProtons==nProtonDau){
172  hPtVsYGen->Fill(ptD,yD);
173  if(TMath::Abs(yD)<0.5){
174  hPtVsYGenLimAcc->Fill(ptD,yD);
175  }
176  if(IsInFiducialAcceptance(ptD,yD)){
177  if(nPionsInAcc==nPionDau && nKaonsInAcc==nKaonDau && nProtonsInAcc==nProtonDau){
178  hPtVsYGenAcc->Fill(ptD,yD);
179  }
180  }
181  }
182  }
183  delete vec;
184  }
185 
186  TH1D* hPtGenAcc=(TH1D*)hPtVsYGenAcc->ProjectionX("hPtGenAcc");
187  hPtGenAcc->GetYaxis()->SetTitle("Entries");
188  TH1D* hPtGenLimAcc=(TH1D*)hPtVsYGenLimAcc->ProjectionX("hPtGenLimAcc");
189  hPtGenLimAcc->GetYaxis()->SetTitle("Entries");
190  hPtGenAcc->Sumw2();
191  hPtGenLimAcc->Sumw2();
192  TH1D* hAccVsPt=(TH1D*)hPtGenAcc->Clone("hAccVsPt");
193  hAccVsPt->Divide(hPtGenAcc,hPtGenLimAcc,1,1,"B");
194  hAccVsPt->GetYaxis()->SetTitle("Acceptance");
195  hAccVsPt->SetStats(0);
196 
197  TCanvas* c2d=new TCanvas("c2d","Pt vs y",1200,600);
198  c2d->Divide(3,1);
199  c2d->cd(1);
200  hPtVsYGen->Draw("colz");
201  c2d->cd(2);
202  hPtVsYGenLimAcc->Draw("colz");
203  c2d->cd(3);
204  hPtVsYGenAcc->Draw("colz");
205 
206  TCanvas* c1d=new TCanvas("c1d","Acceptance",1200,600);
207  c1d->Divide(2,1);
208  c1d->cd(1);
209  hPtGenLimAcc->Draw("");
210  Double_t ymax=1.2*TMath::Max(hPtGenLimAcc->GetMaximum(),hPtGenAcc->GetMaximum());
211  hPtGenLimAcc->SetMaximum(ymax);
212  gPad->Update();
213  TPaveStats *st1=(TPaveStats*)hPtGenLimAcc->GetListOfFunctions()->FindObject("stats");
214  st1->SetY1NDC(0.71);
215  st1->SetY2NDC(0.9);
216  hPtGenAcc->SetLineColor(kRed+1);
217  hPtGenAcc->Draw("sames");
218  gPad->Update();
219  TPaveStats *st2=(TPaveStats*)hPtGenAcc->GetListOfFunctions()->FindObject("stats");
220  st2->SetY1NDC(0.51);
221  st2->SetY2NDC(0.7);
222  st2->SetTextColor(hPtGenAcc->GetLineColor());
223  gPad->Modified();
224  c1d->cd(2);
225  hAccVsPt->Draw();
226 
227 
228  TFile* outfil=new TFile(outFileName.Data(),"recreate");
229  hPtVsYGen->Write();
230  hPtVsYGenLimAcc->Write();
231  hPtVsYGenAcc->Write();
232  hPtGenLimAcc->Write();
233  hPtGenAcc->Write();
234  hAccVsPt->Write();
235  outfil->Close();
236 
237 }
238 
239 //___________________________________________________
240 Bool_t IsInFiducialAcceptance(Double_t pt, Double_t y){
241  // check fiducial acceptance
242 
244  if(TMath::Abs(y) > fYMaxFidAccCut) return kFALSE;
245  else return kTRUE;
246  }
247 
248  if(pt > 5.) {
249  if (TMath::Abs(y) > 0.8) return kFALSE;
250  } else {
251  Double_t maxFiducialY = -0.2/15*pt*pt+1.9/15*pt+0.5;
252  Double_t minFiducialY = 0.2/15*pt*pt-1.9/15*pt-0.5;
253  if (y < minFiducialY || y > maxFiducialY) return kFALSE;
254  }
255 
256  return kTRUE;
257 }
258 //___________________________________________________
259 Bool_t CountKpi(TClonesArray *array, Int_t nentries, Int_t &nPions, Int_t &nKaons, Int_t &nPionsInAcc, Int_t &nKaonsInAcc){
260  // count K and pi in Acc
261 
262  TParticle* dmes=(TParticle*)array->At(0);
263  Double_t sumPx=0;
264  Double_t sumPy=0;
265  Double_t sumPz=0;
266  for(int j=0; j<nentries; j++){
267  TParticle * o = (TParticle*)array->At(j);
268  Int_t pdgdau=TMath::Abs(o->GetPdgCode());
269  if(fDebugLevel>0) printf("%d ",pdgdau);
270  if(pdgdau==130) {
271  if(fDebugLevel>0) printf("K0 dacaying into K0L\n");
272  return kFALSE;
273  }
274  Float_t ptdau=TMath::Sqrt(o->Px()*o->Px()+o->Py()*o->Py());
275  Float_t etadau=o->Eta();
276  if(pdgdau==211){
277  nPions++;
278  sumPx+=o->Px();
279  sumPy+=o->Py();
280  sumPz+=o->Pz();
281  }
282  if(pdgdau==321){
283  nKaons++;
284  sumPx+=o->Px();
285  sumPy+=o->Py();
286  sumPz+=o->Pz();
287  }
288  if(TMath::Abs(etadau)<fEtaMaxDau && ptdau>fPtMinDau){
289  if(pdgdau==211) nPionsInAcc++;
290  if(pdgdau==321) nKaonsInAcc++;
291  }
292  }
293  if(fDebugLevel>0) printf("\n");
294  if(TMath::Abs(sumPx-dmes->Px())>0.001 ||
295  TMath::Abs(sumPy-dmes->Py())>0.001 ||
296  TMath::Abs(sumPz-dmes->Pz())>0.001){
297  printf("Momentum conservation violation\n");
298  return kFALSE;
299  }
300  return kTRUE;
301 }
302 
303 
304 //___________________________________________________
305 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){
306  // count K and pi in Acc
307 
308  TParticle* dmes=(TParticle*)array->At(0);
309  Double_t sumPx=0;
310  Double_t sumPy=0;
311  Double_t sumPz=0;
312 
313  for(int j=0; j<nentries; j++){
314  TParticle * o = (TParticle*)array->At(j);
315  Int_t pdgdau=TMath::Abs(o->GetPdgCode());
316  if(fDebugLevel>0) printf("%d ",pdgdau);
317  if(pdgdau==130) {
318  if(fDebugLevel>0) printf("K0 dacaying into K0L\n");
319  return kFALSE;
320  }
321  Float_t ptdau=TMath::Sqrt(o->Px()*o->Px()+o->Py()*o->Py());
322  Float_t etadau=o->Eta();
323  if(pdgdau==211){
324  nPions++;
325  sumPx+=o->Px();
326  sumPy+=o->Py();
327  sumPz+=o->Pz();
328  }
329  if(pdgdau==321){
330  nKaons++;
331  sumPx+=o->Px();
332  sumPy+=o->Py();
333  sumPz+=o->Pz();
334  }
335  if(pdgdau==2212){
336  nProtons++;
337  sumPx+=o->Px();
338  sumPy+=o->Py();
339  sumPz+=o->Pz();
340  }
341  if(TMath::Abs(etadau)<fEtaMaxDau && ptdau>fPtMinDau){
342  if(pdgdau==211) nPionsInAcc++;
343  if(pdgdau==321) nKaonsInAcc++;
344  if(pdgdau==2212) nProtonsInAcc++;
345  }
346  }
347  if(fDebugLevel>0) printf("\n");
348  if(TMath::Abs(sumPx-dmes->Px())>0.001 ||
349  TMath::Abs(sumPy-dmes->Py())>0.001 ||
350  TMath::Abs(sumPz-dmes->Pz())>0.001){
351  printf("Momentum conservation violation\n");
352  return kFALSE;
353  }
354  return kTRUE;
355 }
356 
357 
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
TSystem * gSystem
Int_t fOptionYFiducial
void ComputeAcceptance()
Double_t fPtMinDau
Double_t massD
Int_t totTrials
Int_t fDebugLevel
EPtShape
Int_t fDDecay
TString fDecayTableFileName
Double_t fYMaxFidAccCut
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)