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