AliPhysics  a0db429 (a0db429)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros
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{
96  printf("ERROR: Wrong decay selected\n");
97  return;
98  }
99  if(fOptionYFiducial==kFixedY) outFileName.Append(Form("yfid%02d_",(Int_t)(fYMaxFidAccCut*10)));
100  else outFileName.Append("yfidPtDep_");
101  outFileName.Append(Form("etaDau%02d_",(Int_t)(fEtaMaxDau*10)));
102  outFileName.Append(Form("ptDau%d_",(Int_t)(fPtMinDau*1000)));
103  TDatabasePDG* db=TDatabasePDG::Instance();
104  Float_t massD=db->GetParticle(pdgCode)->Mass();
105  TClonesArray *array = new TClonesArray("TParticle",100);
106 
107  TH2D* hPtVsYGen=new TH2D("hPtVsYGen","",400,0.,40.,20.,-1.,1.);
108  hPtVsYGen->GetXaxis()->SetTitle("p_{T} (GeV/c)");
109  hPtVsYGen->GetYaxis()->SetTitle("y");
110  TH2D* hPtVsYGenLimAcc=new TH2D("hPtVsYGenLimAcc","",400,0.,40.,20.,-1.,1.);
111  hPtVsYGenLimAcc->GetXaxis()->SetTitle("p_{T} (GeV/c)");
112  hPtVsYGenLimAcc->GetYaxis()->SetTitle("y");
113  TH2D* hPtVsYGenAcc=new TH2D("hPtVsYGenAcc","",400,0.,40.,20.,-1.,1.);
114  hPtVsYGenAcc->GetXaxis()->SetTitle("p_{T} (GeV/c)");
115  hPtVsYGenAcc->GetYaxis()->SetTitle("y");
116 
117  TF1* funcPt=0x0;
118  if(fPtShape==kFONLL7TeV){
119  funcPt=new TF1("fFONLL","[0]*x/TMath::Power((1+TMath::Power(x/[1],[3])),[2])",0.,40.);
120  funcPt->SetParameters(0.322643,2.96275,2.30301,2.5);
121  outFileName.Append("FONLL7ptshape.root");
122  }else if(fPtShape==kFONLL5TeV){
123  funcPt=new TF1("fFONLL","[0]*x/TMath::Power((1+TMath::Power(x/[1],[3])),[2])",0.,40.);
124  funcPt->SetParameters(0.302879,2.9750,3.68139,1.68855);
125  outFileName.Append("FONLL5ptshape.root");
126  }else if(fPtShape==kPythia7TeV){
127  funcPt=new TF1("fFONLL","[0]*x/TMath::Power((1+TMath::Power(x/[1],[3])),[2])",0.,40.);
128  funcPt->SetParameters(0.322643,1.94635,1.40463,2.5);
129  outFileName.Append("PYTHIA7ptshape.root");
130  }else{
131  funcPt=new TF1("fFlat","pol0",0.,40.);
132  funcPt->SetParameter(0,1.);
133  outFileName.Append("flatpt.root");
134  }
135  TRandom3* gener=new TRandom3(0);
136 
137  for(Int_t itry=0; itry<totTrials; itry++){
138  if(itry%10000==0) printf("Event %d\n",itry);
139  Float_t ptD=funcPt->GetRandom();
140  Float_t phiD=gener->Rndm()*2*TMath::Pi();
141  Float_t yD=gener->Rndm()*2.-1.; // flat in -1<y<1
142  Float_t px=ptD*TMath::Cos(phiD);
143  Float_t py=ptD*TMath::Sin(phiD);
144  Float_t mt=TMath::Sqrt(massD*massD+ptD*ptD);
145  Float_t pz=mt*TMath::SinH(yD);
146  Float_t E=TMath::Sqrt(massD*massD+px*px+py*py+pz*pz);
147 
148  TLorentzVector* vec=new TLorentzVector(px,py,pz,E);
149  pdec->Decay(pdgCode,vec);
150  array->Clear();
151  Int_t nentries = pdec->ImportParticles(array);
152  TParticle* dmes=(TParticle*)array->At(0);
153  Int_t nDaughters=dmes->GetNDaughters();
154  if(fDDecay==kD0Kpi && nDaughters!=2) return;
155  Int_t nPionsInAcc=0;
156  Int_t nProtonsInAcc=0;
157  Int_t nKaonsInAcc=0;
158  Int_t nPions=0;
159  Int_t nProtons=0;
160  Int_t nKaons=0;
161  Bool_t isOk=CountPKpi(array,nentries,nPions,nKaons,nProtons,nPionsInAcc,nKaonsInAcc,nProtonsInAcc);
162 
163  if(isOk){
164  if(nPions==nPionDau && nKaons==nKaonDau && nProtons==nProtonDau){
165  hPtVsYGen->Fill(ptD,yD);
166  if(TMath::Abs(yD)<0.5){
167  hPtVsYGenLimAcc->Fill(ptD,yD);
168  }
169  if(IsInFiducialAcceptance(ptD,yD)){
170  if(nPionsInAcc==nPionDau && nKaonsInAcc==nKaonDau && nProtonsInAcc==nProtonDau){
171  hPtVsYGenAcc->Fill(ptD,yD);
172  }
173  }
174  }
175  }
176  delete vec;
177  }
178 
179  TH1D* hPtGenAcc=(TH1D*)hPtVsYGenAcc->ProjectionX("hPtGenAcc");
180  hPtGenAcc->GetYaxis()->SetTitle("Entries");
181  TH1D* hPtGenLimAcc=(TH1D*)hPtVsYGenLimAcc->ProjectionX("hPtGenLimAcc");
182  hPtGenLimAcc->GetYaxis()->SetTitle("Entries");
183  hPtGenAcc->Sumw2();
184  hPtGenLimAcc->Sumw2();
185  TH1D* hAccVsPt=(TH1D*)hPtGenAcc->Clone("hAccVsPt");
186  hAccVsPt->Divide(hPtGenAcc,hPtGenLimAcc,1,1,"B");
187  hAccVsPt->GetYaxis()->SetTitle("Acceptance");
188  hAccVsPt->SetStats(0);
189 
190  TCanvas* c2d=new TCanvas("c2d","Pt vs y",1200,600);
191  c2d->Divide(3,1);
192  c2d->cd(1);
193  hPtVsYGen->Draw("colz");
194  c2d->cd(2);
195  hPtVsYGenLimAcc->Draw("colz");
196  c2d->cd(3);
197  hPtVsYGenAcc->Draw("colz");
198 
199  TCanvas* c1d=new TCanvas("c1d","Acceptance",1200,600);
200  c1d->Divide(2,1);
201  c1d->cd(1);
202  hPtGenLimAcc->Draw("");
203  Double_t ymax=1.2*TMath::Max(hPtGenLimAcc->GetMaximum(),hPtGenAcc->GetMaximum());
204  hPtGenLimAcc->SetMaximum(ymax);
205  gPad->Update();
206  TPaveStats *st1=(TPaveStats*)hPtGenLimAcc->GetListOfFunctions()->FindObject("stats");
207  st1->SetY1NDC(0.71);
208  st1->SetY2NDC(0.9);
209  hPtGenAcc->SetLineColor(kRed+1);
210  hPtGenAcc->Draw("sames");
211  gPad->Update();
212  TPaveStats *st2=(TPaveStats*)hPtGenAcc->GetListOfFunctions()->FindObject("stats");
213  st2->SetY1NDC(0.51);
214  st2->SetY2NDC(0.7);
215  st2->SetTextColor(hPtGenAcc->GetLineColor());
216  gPad->Modified();
217  c1d->cd(2);
218  hAccVsPt->Draw();
219 
220 
221  TFile* outfil=new TFile(outFileName.Data(),"recreate");
222  hPtVsYGen->Write();
223  hPtVsYGenLimAcc->Write();
224  hPtVsYGenAcc->Write();
225  hPtGenLimAcc->Write();
226  hPtGenAcc->Write();
227  hAccVsPt->Write();
228  outfil->Close();
229 
230 }
231 
232 //___________________________________________________
233 Bool_t IsInFiducialAcceptance(Double_t pt, Double_t y){
234  // check fiducial acceptance
235 
237  if(TMath::Abs(y) > fYMaxFidAccCut) return kFALSE;
238  else return kTRUE;
239  }
240 
241  if(pt > 5.) {
242  if (TMath::Abs(y) > 0.8) return kFALSE;
243  } else {
244  Double_t maxFiducialY = -0.2/15*pt*pt+1.9/15*pt+0.5;
245  Double_t minFiducialY = 0.2/15*pt*pt-1.9/15*pt-0.5;
246  if (y < minFiducialY || y > maxFiducialY) return kFALSE;
247  }
248 
249  return kTRUE;
250 }
251 //___________________________________________________
252 Bool_t CountKpi(TClonesArray *array, Int_t nentries, Int_t &nPions, Int_t &nKaons, Int_t &nPionsInAcc, Int_t &nKaonsInAcc){
253  // count K and pi in Acc
254 
255  TParticle* dmes=(TParticle*)array->At(0);
256  Double_t sumPx=0;
257  Double_t sumPy=0;
258  Double_t sumPz=0;
259  for(int j=0; j<nentries; j++){
260  TParticle * o = (TParticle*)array->At(j);
261  Int_t pdgdau=TMath::Abs(o->GetPdgCode());
262  if(fDebugLevel>0) printf("%d ",pdgdau);
263  Float_t ptdau=TMath::Sqrt(o->Px()*o->Px()+o->Py()*o->Py());
264  Float_t etadau=o->Eta();
265  if(pdgdau==211){
266  nPions++;
267  sumPx+=o->Px();
268  sumPy+=o->Py();
269  sumPz+=o->Pz();
270  }
271  if(pdgdau==321){
272  nKaons++;
273  sumPx+=o->Px();
274  sumPy+=o->Py();
275  sumPz+=o->Pz();
276  }
277  if(TMath::Abs(etadau)<fEtaMaxDau && ptdau>fPtMinDau){
278  if(pdgdau==211) nPionsInAcc++;
279  if(pdgdau==321) nKaonsInAcc++;
280  }
281  }
282  if(fDebugLevel>0) printf("\n");
283  if(TMath::Abs(sumPx-dmes->Px())>0.001 ||
284  TMath::Abs(sumPy-dmes->Py())>0.001 ||
285  TMath::Abs(sumPz-dmes->Pz())>0.001){
286  printf("Momentum conservation violation\n");
287  return kFALSE;
288  }
289  return kTRUE;
290 }
291 
292 
293 //___________________________________________________
294 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){
295  // count K and pi in Acc
296 
297  TParticle* dmes=(TParticle*)array->At(0);
298  Double_t sumPx=0;
299  Double_t sumPy=0;
300  Double_t sumPz=0;
301 
302  for(int j=0; j<nentries; j++){
303  TParticle * o = (TParticle*)array->At(j);
304  Int_t pdgdau=TMath::Abs(o->GetPdgCode());
305  if(fDebugLevel>0) printf("%d ",pdgdau);
306  Float_t ptdau=TMath::Sqrt(o->Px()*o->Px()+o->Py()*o->Py());
307  Float_t etadau=o->Eta();
308  if(pdgdau==211){
309  nPions++;
310  sumPx+=o->Px();
311  sumPy+=o->Py();
312  sumPz+=o->Pz();
313  }
314  if(pdgdau==321){
315  nKaons++;
316  sumPx+=o->Px();
317  sumPy+=o->Py();
318  sumPz+=o->Pz();
319  }
320  if(pdgdau==2212){
321  nProtons++;
322  sumPx+=o->Px();
323  sumPy+=o->Py();
324  sumPz+=o->Pz();
325  }
326  if(TMath::Abs(etadau)<fEtaMaxDau && ptdau>fPtMinDau){
327  if(pdgdau==211) nPionsInAcc++;
328  if(pdgdau==321) nKaonsInAcc++;
329  if(pdgdau==2212) nProtonsInAcc++;
330  }
331  }
332  if(fDebugLevel>0) printf("\n");
333  if(TMath::Abs(sumPx-dmes->Px())>0.001 ||
334  TMath::Abs(sumPy-dmes->Py())>0.001 ||
335  TMath::Abs(sumPz-dmes->Pz())>0.001){
336  printf("Momentum conservation violation\n");
337  return kFALSE;
338  }
339  return kTRUE;
340 }
341 
342 
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)