AliPhysics  ef3b16e (ef3b16e)
AliAnalysisTaskCheckHFMCProd.cxx
Go to the documentation of this file.
1 #include "AliAnalysisTaskSE.h"
2 #include "AliAnalysisManager.h"
3 #include "AliAnalysisDataContainer.h"
4 #include "AliESDEvent.h"
5 #include "AliESDtrackCuts.h"
6 #include "AliCentrality.h"
7 #include "AliMCEventHandler.h"
8 #include "AliMCEvent.h"
9 #include "AliHeader.h"
10 #include "AliGenCocktailEventHeader.h"
11 #include "AliGenHijingEventHeader.h"
12 #include "AliGenEventHeader.h"
13 #include "AliCollisionGeometry.h"
14 #include "AliGenerator.h"
15 #include "AliVertexingHFUtils.h"
16 #include "AliMultiplicity.h"
17 #include <TParticle.h>
18 #include <TSystem.h>
19 #include <TTree.h>
20 #include <TNtuple.h>
21 #include <TH1F.h>
22 #include <TH2F.h>
23 #include <TH3F.h>
24 #include <TChain.h>
25 #include "AliESDInputHandlerRP.h"
27 
28 /**************************************************************************
29  * Copyright(c) 1998-2012, ALICE Experiment at CERN, All rights reserved. *
30  * *
31  * Author: The ALICE Off-line Project. *
32  * Contributors are mentioned in the code where appropriate. *
33  * *
34  * Permission to use, copy, modify and distribute this software and its *
35  * documentation strictly for non-commercial purposes is hereby granted *
36  * without fee, provided that the above copyright notice appears in all *
37  * copies and that both the copyright notice and this permission notice *
38  * appear in the supporting documentation. The authors make no claims *
39  * about the suitability of this software for any purpose. It is *
40  * provided "as is" without express or implied warranty. *
41  **************************************************************************/
42 
43 /* $Id$ */
44 
45 //*************************************************************************
46 // Implementation of class AliAnalysisTaskCheckHFMCProd
47 // AliAnalysisTask to check MC production at ESD+Kine level
48 //
49 //
50 // Authors: F. Prino, prino@to.infn.it
51 //
52 //*************************************************************************
53 
57 
58 //______________________________________________________________________________
60  AliAnalysisTaskSE("HFMCChecks"),
61  fOutput(0),
62  fHistoNEvents(0),
63  fHistoPhysPrim(0),
64  fHistoPhysPrimPiKPi09(0),
65  fHistoPhysPrimPiKPi09vsb(0),
66  fHistoTracks(0),
67  fHistoSelTracks(0),
68  fHistoTracklets(0),
69  fHistoTrackletsEta1(0),
70  fHistoPtPhysPrim(0),
71  fHistoEtaPhysPrim(0),
72  fHistoSPD3DVtxX(0),
73  fHistoSPD3DVtxY(0),
74  fHistoSPD3DVtxZ(0),
75  fHistoSPDZVtxX(0),
76  fHistoSPDZVtxY(0),
77  fHistoSPDZVtxZ(0),
78  fHistoTRKVtxX(0),
79  fHistoTRKVtxY(0),
80  fHistoTRKVtxZ(0),
81  fHistoNcharmed(0),
82  fHistoNbVsNc(0),
83  fHistOriginPrompt(0),
84  fHistOriginFeeddown(0),
85  fHistMotherID(0),
86  fHistDSpecies(0),
87  fHistBSpecies(0),
88  fHistLcDecayChan(0),
89  fHistNcollHFtype(0),
90  fHistNinjectedvsb(0),
91  fHistEtaPhiPtGenEle(0),
92  fHistEtaPhiPtGenPi(0),
93  fHistEtaPhiPtGenK(0),
94  fHistEtaPhiPtGenPro(0),
95  fHistEtaPhiPtRecEle(0),
96  fHistEtaPhiPtRecPi(0),
97  fHistEtaPhiPtRecK(0),
98  fHistEtaPhiPtRecPro(0),
99  fHistPtRecVsPtGen(0),
100  fHistPhiRecVsPhiGen(0),
101  fHistEtaRecVsEtaGen(0),
102  fHistPtRecGood(0),
103  fHistPtRecFake(0),
104  fSearchUpToQuark(kFALSE),
105  fSystem(0),
106  fESDtrackCuts(0x0),
107  fReadMC(kTRUE)
108 {
109  //
110  for(Int_t i=0; i<5; i++){
111  fHistBYPtAllDecay[i]=0x0;
112  fHistYPtAllDecay[i]=0x0;
113  fHistYPtPromptAllDecay[i]=0x0;
115  fHistYPtPrompt[i]=0x0;
116  fHistYPtFeeddown[i]=0x0;
117  }
118  for(Int_t i=0; i<2; i++){
119  fHistYPtD0byDecChannel[i]=0x0;
121  fHistYPtDsbyDecChannel[i]=0x0;
122  }
124 
125  DefineInput(0, TChain::Class());
126  DefineOutput(1, TList::Class());
127 }
128 
129 
130 //___________________________________________________________________________
132  //
133  if (AliAnalysisManager::GetAnalysisManager()->IsProofMode()) return;
134  if (fOutput) {
135  delete fOutput;
136  fOutput = 0;
137  }
138  delete fESDtrackCuts;
139 
140 }
141 
142 //___________________________________________________________________________
145 
146  fOutput = new TList();
147  fOutput->SetOwner();
148  fOutput->SetName("OutputHistos");
149 
150  fHistoNEvents = new TH1F("hNEvents", "Number of processed events",3,-0.5,2.5);
151  fHistoNEvents->SetMinimum(0);
152  fOutput->Add(fHistoNEvents);
153 
154  Double_t maxMult=100.;
155  if(fSystem==1) maxMult=10000.;
156  if(fSystem==2) maxMult=500.;
157  fHistoPhysPrim = new TH1F("hPhysPrim","",100,-0.5,maxMult-0.5);
158  fOutput->Add(fHistoPhysPrim);
159  fHistoPhysPrimPiKPi09 = new TH1F("hPhysPrimPiKPi09","",100,-0.5,maxMult-0.5);
161  fHistoPhysPrimPiKPi09vsb = new TH2F("hPhysPrimPiKPi09vsb","",50,0.,20.,100,-0.5,maxMult-0.5);
163 
164  fHistoTracks = new TH1F("hTracks","",100,-0.5,maxMult*2-0.5);
165  fOutput->Add(fHistoTracks);
166  fHistoSelTracks = new TH1F("hSelTracks","",100,-0.5,maxMult-0.5);
167  fOutput->Add(fHistoSelTracks);
168  fHistoTracklets = new TH1F("hTracklets","",100,-0.5,maxMult-0.5);
169  fOutput->Add(fHistoTracklets);
170  fHistoTrackletsEta1 = new TH1F("hTrackletsEta1","",100,-0.5,maxMult-0.5);
172  fHistoPtPhysPrim = new TH1F("hPtPhysPrim","",100,0.,20.);
174  fHistoEtaPhysPrim = new TH1F("hEtaPhysPrim","",100,-10.,10.);
176 
177  fHistoSPD3DVtxX = new TH1F("hSPD3DvX","",100,-1.,1.);
178  fOutput->Add(fHistoSPD3DVtxX);
179  fHistoSPD3DVtxY = new TH1F("hSPD3DvY","",100,-1.,1.);
180  fOutput->Add(fHistoSPD3DVtxY);
181  fHistoSPD3DVtxZ = new TH1F("hSPD3DvZ","",100,-15.,15.);
182  fOutput->Add(fHistoSPD3DVtxZ);
183 
184  fHistoSPDZVtxX = new TH1F("hSPDZvX","",100,-1.,1.);
185  fOutput->Add(fHistoSPDZVtxX);
186  fHistoSPDZVtxY = new TH1F("hSPDZvY","",100,-1.,1.);
187  fOutput->Add(fHistoSPDZVtxY);
188  fHistoSPDZVtxZ = new TH1F("hSPDZvZ","",100,-15.,15.);
189  fOutput->Add(fHistoSPDZVtxZ);
190 
191 
192  fHistoTRKVtxX = new TH1F("hTRKvX","",100,-1.,1.);
193  fOutput->Add(fHistoTRKVtxX);
194  fHistoTRKVtxY = new TH1F("hTRKvY","",100,-1.,1.);
195  fOutput->Add(fHistoTRKVtxY);
196  fHistoTRKVtxZ = new TH1F("hTRKvZ","",100,-15.,15.);
197  fOutput->Add(fHistoTRKVtxZ);
198 
199  Int_t nBinscb=11;
200  if(fSystem==1) nBinscb=200;
201  if(fSystem==2) nBinscb=21;
202  Double_t maxncn=nBinscb-0.5;
203  fHistoNcharmed = new TH2F("hncharmed","",100,-0.5,maxMult-0.5,nBinscb,-0.5,maxncn);
204  fOutput->Add(fHistoNcharmed);
205  fHistoNbVsNc = new TH2F("hnbvsnc","",nBinscb,-0.5,maxncn,nBinscb,-0.5,maxncn);
206  fOutput->Add(fHistoNbVsNc);
207 
208  fHistYPtPrompt[0] = new TH2F("hyptD0prompt","D0 - Prompt",40,0.,40.,20,-2.,2.);
209  fHistYPtPrompt[1] = new TH2F("hyptDplusprompt","Dplus - Prompt",40,0.,40.,20,-2.,2.);
210  fHistYPtPrompt[2] = new TH2F("hyptDstarprompt","Dstar - Prompt",40,0.,40.,20,-2.,2.);
211  fHistYPtPrompt[3] = new TH2F("hyptDsprompt","Ds - Prompt",40,0.,40.,20,-2.,2.);
212  fHistYPtPrompt[4] = new TH2F("hyptLcprompt","Lc - Prompt",40,0.,40.,20,-2.,2.);
213 
214  fHistBYPtAllDecay[0] = new TH2F("hyptB0AllDecay","B0 - All",40,0.,40.,40,-2.,2.);
215  fHistBYPtAllDecay[1] = new TH2F("hyptBplusAllDecay","Bplus - All",40,0.,40.,40,-2.,2.);
216  fHistBYPtAllDecay[2] = new TH2F("hyptBstarAllDecay","Bstar - All",40,0.,40.,40,-2.,2.);
217  fHistBYPtAllDecay[3] = new TH2F("hyptBsAllDecay","Bs - All",40,0.,40.,40,-2.,2.);
218  fHistBYPtAllDecay[4] = new TH2F("hyptLbAllDecay","LB - All",40,0.,40.,40,-2.,2.);
219 
220  fHistYPtAllDecay[0] = new TH2F("hyptD0AllDecay","D0 - All",40,0.,40.,40,-2.,2.);
221  fHistYPtAllDecay[1] = new TH2F("hyptDplusAllDecay","Dplus - All",40,0.,40.,40,-2.,2.);
222  fHistYPtAllDecay[2] = new TH2F("hyptDstarAllDecay","Dstar - All",40,0.,40.,40,-2.,2.);
223  fHistYPtAllDecay[3] = new TH2F("hyptDsAllDecay","Ds - All",40,0.,40.,40,-2.,2.);
224  fHistYPtAllDecay[4] = new TH2F("hyptLcAllDecay","Lc - All",40,0.,40.,40,-2.,2.);
225 
226  fHistYPtPromptAllDecay[0] = new TH2F("hyptD0promptAllDecay","D0 - Prompt",40,0.,40.,40,-2.,2.);
227  fHistYPtPromptAllDecay[1] = new TH2F("hyptDpluspromptAllDecay","Dplus - Prompt",40,0.,40.,40,-2.,2.);
228  fHistYPtPromptAllDecay[2] = new TH2F("hyptDstarpromptAllDecay","Dstar - Prompt",40,0.,40.,40,-2.,2.);
229  fHistYPtPromptAllDecay[3] = new TH2F("hyptDspromptAllDecay","Ds - Prompt",40,0.,40.,40,-2.,2.);
230  fHistYPtPromptAllDecay[4] = new TH2F("hyptLcpromptAllDecay","Lc - Prompt",40,0.,40.,40,-2.,2.);
231 
232  fHistYPtFeeddownAllDecay[0] = new TH2F("hyptD0feeddownAllDecay","D0 - FromB",40,0.,40.,40,-2.,2.);
233  fHistYPtFeeddownAllDecay[1] = new TH2F("hyptDplusfeeddownAllDecay","Dplus - FromB",40,0.,40.,40,-2.,2.);
234  fHistYPtFeeddownAllDecay[2] = new TH2F("hyptDstarfeeddownAllDecay","Dstar - FromB",40,0.,40.,40,-2.,2.);
235  fHistYPtFeeddownAllDecay[3] = new TH2F("hyptDsfeeddownAllDecay","Ds - FromB",40,0.,40.,40,-2.,2.);
236  fHistYPtFeeddownAllDecay[4] = new TH2F("hyptLcfeeddownAllDecay","Lc - FromB",40,0.,40.,40,-2.,2.);
237 
238 
239  fHistYPtFeeddown[0] = new TH2F("hyptD0feeddown","D0 - Feeddown",40,0.,40.,20,-2.,2.);
240  fHistYPtFeeddown[1] = new TH2F("hyptDplusfeeddown","Dplus - Feeddown",40,0.,40.,20,-2.,2.);
241  fHistYPtFeeddown[2] = new TH2F("hyptDstarfeedown","Dstar - Feeddown",40,0.,40.,20,-2.,2.);
242  fHistYPtFeeddown[3] = new TH2F("hyptDsfeedown","Ds - Feeddown",40,0.,40.,20,-2.,2.);
243  fHistYPtFeeddown[4] = new TH2F("hyptLcfeedown","Lc - Feeddown",40,0.,40.,20,-2.,2.);
244 
245  for(Int_t ih=0; ih<5; ih++){
246  fHistBYPtAllDecay[ih]->SetMinimum(0);
247  fOutput->Add(fHistBYPtAllDecay[ih]);
248  fHistYPtAllDecay[ih]->SetMinimum(0);
249  fOutput->Add(fHistYPtAllDecay[ih]);
250  fHistYPtPromptAllDecay[ih]->SetMinimum(0);
252  fHistYPtFeeddownAllDecay[ih]->SetMinimum(0);
254  fHistYPtPrompt[ih]->SetMinimum(0);
255  fOutput->Add(fHistYPtPrompt[ih]);
256  fHistYPtFeeddown[ih]->SetMinimum(0);
257  fOutput->Add(fHistYPtFeeddown[ih]);
258  }
259 
260  fHistYPtD0byDecChannel[0] = new TH2F("hyptD02","D0 - 2prong",40,0.,40.,20,-2.,2.);
261  fHistYPtD0byDecChannel[1] = new TH2F("hyptD04","D0 - 4prong",40,0.,40.,20,-2.,2.);
262  fHistYPtDplusbyDecChannel[0] = new TH2F("hyptDplusnonreson","Dplus - non reson",40,0.,40.,20,-2.,2.);
263  fHistYPtDplusbyDecChannel[1] = new TH2F("hyptDplusreson","Dplus - reson via K0*",40,0.,40.,20,-2.,2.);
264  fHistYPtDplusbyDecChannel[2] = new TH2F("hyptDplusKKpi","Dplus -> KKpi",40,0.,40.,20,-2.,2.);
265  fHistYPtDsbyDecChannel[0] = new TH2F("hyptDsphi","Ds - vis Phi",40,0.,40.,20,-2.,2.);
266  fHistYPtDsbyDecChannel[1] = new TH2F("hyptDsk0st","Ds - via k0*",40,0.,40.,20,-2.,2.);
267 
268  for(Int_t ih=0; ih<2; ih++){
269 
270  fHistYPtD0byDecChannel[ih]->SetMinimum(0);
272  fHistYPtDplusbyDecChannel[ih]->SetMinimum(0);
274  fHistYPtDsbyDecChannel[ih]->SetMinimum(0);
276  }
277  fHistYPtDplusbyDecChannel[2]->SetMinimum(0);
279 
280  fHistOriginPrompt=new TH1F("hOriginPrompt","",100,0.,0.5);
281  fHistOriginPrompt->SetMinimum(0);
283  fHistOriginFeeddown=new TH1F("hOriginFeeddown","",100,0.,0.5);
284  fHistOriginFeeddown->SetMinimum(0);
286  fHistMotherID=new TH1F("hMotherID","",1000,-1.5,998.5);
287  fHistMotherID->SetMinimum(0);
288  fOutput->Add(fHistMotherID);
289  fHistDSpecies=new TH1F("hDSpecies","",10,-0.5,9.5);
290  fHistDSpecies->GetXaxis()->SetBinLabel(1,"D0");
291  fHistDSpecies->GetXaxis()->SetBinLabel(2,"D0bar");
292  fHistDSpecies->GetXaxis()->SetBinLabel(3,"D+");
293  fHistDSpecies->GetXaxis()->SetBinLabel(4,"D-");
294  fHistDSpecies->GetXaxis()->SetBinLabel(5,"D*+");
295  fHistDSpecies->GetXaxis()->SetBinLabel(6,"D*-");
296  fHistDSpecies->GetXaxis()->SetBinLabel(7,"Ds+");
297  fHistDSpecies->GetXaxis()->SetBinLabel(8,"Ds-");
298  fHistDSpecies->GetXaxis()->SetBinLabel(9,"Lc+");
299  fHistDSpecies->GetXaxis()->SetBinLabel(10,"Lc-");
300  fHistDSpecies->SetMinimum(0);
301  fOutput->Add(fHistDSpecies);
302  fHistBSpecies=new TH1F("hBSpecies","",10,-0.5,9.5);
303  fHistBSpecies->GetXaxis()->SetBinLabel(1,"B0");
304  fHistBSpecies->GetXaxis()->SetBinLabel(2,"B0bar");
305  fHistBSpecies->GetXaxis()->SetBinLabel(3,"B+");
306  fHistBSpecies->GetXaxis()->SetBinLabel(4,"B-");
307  fHistBSpecies->GetXaxis()->SetBinLabel(5,"B*+");
308  fHistBSpecies->GetXaxis()->SetBinLabel(6,"B*-");
309  fHistBSpecies->GetXaxis()->SetBinLabel(7,"Bs+");
310  fHistBSpecies->GetXaxis()->SetBinLabel(8,"Bs-");
311  fHistBSpecies->GetXaxis()->SetBinLabel(9,"Lb+");
312  fHistBSpecies->GetXaxis()->SetBinLabel(10,"Lb-");
313  fHistBSpecies->SetMinimum(0);
314  fOutput->Add(fHistBSpecies);
315  fHistLcDecayChan=new TH1F("hLcDecayChan","",9,-2.5,6.5);
316  fHistLcDecayChan->GetXaxis()->SetBinLabel(1,"Violates p cons");
317  fHistLcDecayChan->GetXaxis()->SetBinLabel(2,"Other decay");
318  fHistLcDecayChan->GetXaxis()->SetBinLabel(3,"Error");
319  fHistLcDecayChan->GetXaxis()->SetBinLabel(4,"pK#pi non res");
320  fHistLcDecayChan->GetXaxis()->SetBinLabel(5,"pK#pi via K*0");
321  fHistLcDecayChan->GetXaxis()->SetBinLabel(6,"pK#pi via #Delta++");
322  fHistLcDecayChan->GetXaxis()->SetBinLabel(7,"pK#pi via #Lambda1520");
323  fHistLcDecayChan->GetXaxis()->SetBinLabel(8,"pK0s#rightarrowp#pi#pi");
324  fHistLcDecayChan->GetXaxis()->SetBinLabel(9,"#pi#Lambda#rightarrowp#pi#pi");
325  fHistLcDecayChan->SetMinimum(0);
327 
328  fHistNcollHFtype=new TH2F("hNcollHFtype","",5,-1.5,3.5,30,-0.5,29.5);
330  fHistNinjectedvsb=new TH2F("hNinjectedvsb","",50,0.,20.,140,-0.5,139.5);
332 
333  Double_t binseta[11]={-1.0,-0.8,-0.6,-0.4,-0.2,0.,0.2,0.4,0.6,0.8,1.0};
334  const Int_t nBinsPhi=40;
335  Double_t binsphi[nBinsPhi+1];
336  for(Int_t ib=0; ib<=nBinsPhi; ib++) binsphi[ib]=ib*TMath::Pi()/20.;
337  const Int_t nBinsPt=24;
338  Double_t binspt[nBinsPt+1]={0.,0.10,0.15,0.2,0.25,
339  0.3,0.4,0.5,0.6,0.7,
340  0.8,0.9,1.,1.25,1.5,
341  1.75,2.,2.5,3.,4.,
342  5.,7.5,10.,15.,20.};
343 
344  fHistEtaPhiPtGenEle=new TH3F("hEtaPhiPtGenEle","",10,binseta,nBinsPhi,binsphi,nBinsPt,binspt);
346  fHistEtaPhiPtGenPi=new TH3F("hEtaPhiPtGenPi","",10,binseta,nBinsPhi,binsphi,nBinsPt,binspt);
348  fHistEtaPhiPtGenK=new TH3F("hEtaPhiPtGenK","",10,binseta,nBinsPhi,binsphi,nBinsPt,binspt);
350  fHistEtaPhiPtGenPro=new TH3F("hEtaPhiPtGenPro","",10,binseta,nBinsPhi,binsphi,nBinsPt,binspt);
352 
353 
354  fHistEtaPhiPtRecEle=new TH3F("hEtaPhiPtRecEle","",10,binseta,nBinsPhi,binsphi,nBinsPt,binspt);
356  fHistEtaPhiPtRecPi=new TH3F("hEtaPhiPtRecPi","",10,binseta,nBinsPhi,binsphi,nBinsPt,binspt);
358  fHistEtaPhiPtRecK=new TH3F("hEtaPhiPtRecK","",10,binseta,nBinsPhi,binsphi,nBinsPt,binspt);
360  fHistEtaPhiPtRecPro=new TH3F("hEtaPhiPtRecPro","",10,binseta,nBinsPhi,binsphi,nBinsPt,binspt);
362 
363  fHistPtRecVsPtGen=new TH2F("hPtRecVsPtGen"," ; particle p_{T} (GeV/c); track p_{T} (GeV/c)",100,0.,10.,100,0.,10.);
365  fHistPhiRecVsPhiGen=new TH2F("hPhiRecVsPhiGen"," ; particle #varphi ; track #varphi",100,0.,2.*TMath::Pi(),100,0.,2.*TMath::Pi());
367  fHistEtaRecVsEtaGen=new TH2F("hEtaRecVsEtaGen"," ; particle #eta ; track #eta",100,0.,2.*TMath::Pi(),100,0.,2.*TMath::Pi());
369 
370  fHistPtRecGood=new TH1F("hPtRecGood"," ; track p_{T} (GeV/c)",100,0.,10.);
371  fOutput->Add(fHistPtRecGood);
372  fHistPtRecFake=new TH1F("hPtRecFake"," ; track p_{T} (GeV/c)",100,0.,10.);
373  fOutput->Add(fHistPtRecFake);
374 
375  PostData(1,fOutput);
376 
377 }
378 //______________________________________________________________________________
380 {
381  //
382 
383  AliESDEvent *esd = (AliESDEvent*) (InputEvent());
384 
385 
386  if(!esd) {
387  printf("AliAnalysisTaskSDDRP::Exec(): bad ESD\n");
388  return;
389  }
390 
391  fHistoNEvents->Fill(0);
392 
393  if(!fESDtrackCuts){
394  Int_t year=2011;
395  if(esd->GetRunNumber()<=139517) year=2010;
396  if(year==2010) fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kFALSE);
397  else fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kFALSE);
398  fESDtrackCuts->SetMaxDCAToVertexXY(2.4);
399  fESDtrackCuts->SetMaxDCAToVertexZ(3.2);
400  fESDtrackCuts->SetDCAToVertex2D(kTRUE);
401  fESDtrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
402  AliESDtrackCuts::kAny);
403  }
404 
405  Int_t nTracks=esd->GetNumberOfTracks();
406  fHistoTracks->Fill(nTracks);
407  Int_t nSelTracks=0;
408  for(Int_t it=0; it<nTracks; it++){
409  AliESDtrack* tr=esd->GetTrack(it);
410  UInt_t status=tr->GetStatus();
411  if(!(status&AliESDtrack::kITSrefit)) continue;
412  if(!(status&AliESDtrack::kTPCin)) continue;
413  nSelTracks++;
414  }
415  fHistoSelTracks->Fill(nSelTracks);
416 
417  const AliMultiplicity* mult=esd->GetMultiplicity();
418  Int_t nTracklets=mult->GetNumberOfTracklets();
419  Int_t nTracklets1=0;
420  for(Int_t it=0; it<nTracklets; it++){
421  Double_t eta=TMath::Abs(mult->GetEta(it));
422  if(eta<1) nTracklets1++;
423  }
424  fHistoTracklets->Fill(nTracklets);
425  fHistoTrackletsEta1->Fill(nTracklets1);
426 
427  const AliESDVertex *spdv=esd->GetVertex();
428  if(spdv && spdv->IsFromVertexer3D()){
429  fHistoSPD3DVtxX->Fill(spdv->GetX());
430  fHistoSPD3DVtxY->Fill(spdv->GetY());
431  fHistoSPD3DVtxZ->Fill(spdv->GetZ());
432  }
433  if(spdv && spdv->IsFromVertexerZ()){
434  fHistoSPDZVtxX->Fill(spdv->GetX());
435  fHistoSPDZVtxY->Fill(spdv->GetY());
436  fHistoSPDZVtxZ->Fill(spdv->GetZ());
437  }
438  const AliESDVertex *trkv=esd->GetPrimaryVertex();
439  if(trkv && trkv->GetNContributors()>1){
440  fHistoTRKVtxX->Fill(trkv->GetX());
441  fHistoTRKVtxY->Fill(trkv->GetY());
442  fHistoTRKVtxZ->Fill(trkv->GetZ());
443  }
444 
445  AliMCEvent* mcEvent = 0x0;
446 
447  if(fReadMC){
448  AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
449  if (!eventHandler) {
450  Printf("ERROR: Could not retrieve MC event handler");
451  return;
452  }
453  mcEvent = eventHandler->MCEvent();
454  if (!mcEvent) {
455  Printf("ERROR: Could not retrieve MC event");
456  return;
457  }
458  const AliVVertex* mcVert=mcEvent->GetPrimaryVertex();
459  if(!mcVert){
460  Printf("ERROR: generated vertex not available");
461  return;
462  }
463  if(TMath::Abs(mcVert->GetZ())>10) return;
464 
465  // const AliHeader* h=(AliHeader*)mcEvent->GetHeader();
466  // cout<<h<<endl;
467  TString genname=mcEvent->GenEventHeader()->ClassName();
468  Int_t nColl=-1;
469  Double_t imppar=-999.;
470  Int_t nInjected=0;
471  Int_t typeHF=-1;
472  TList* lgen=0x0;
473  if(genname.Contains("CocktailEventHeader")){
474  AliGenCocktailEventHeader *cockhead=(AliGenCocktailEventHeader*)mcEvent->GenEventHeader();
475  lgen=cockhead->GetHeaders();
476  for(Int_t ig=0; ig<lgen->GetEntries(); ig++){
477  AliGenerator* gen=(AliGenerator*)lgen->At(ig);
478  TString title=gen->GetName();
479  if(title.Contains("bchadr")){
480  typeHF=1;
481  nInjected++;
482  }else if(title.Contains("chadr")) {
483  typeHF=0;
484  nInjected++;
485  }else if(title.Contains("bele")){
486  typeHF=3;
487  nInjected++;
488  }else if(title.Contains("cele")){
489  typeHF=2;
490  nInjected++;
491  }else if(title.Contains("pythiaHF")){
492  nInjected++;
493  }else if(title.Contains("hijing") || title.Contains("Hijing")){
494  AliGenHijingEventHeader* hijh=(AliGenHijingEventHeader*)lgen->At(ig);
495  imppar=hijh->ImpactParameter();
496  }
497  }
498  nColl=lgen->GetEntries();
499  fHistNcollHFtype->Fill(typeHF,nColl);
500  fHistNinjectedvsb->Fill(imppar,nInjected);
501  }else if(genname.Contains("HijingEventHeader")){
502  AliGenHijingEventHeader* hijh=(AliGenHijingEventHeader*)mcEvent->GenEventHeader();
503  imppar=hijh->ImpactParameter();
504  }else{
505  TString genTitle=mcEvent->GenEventHeader()->GetTitle();
506  if(genTitle.Contains("bchadr")) typeHF=1;
507  else if(genTitle.Contains("chadr")) typeHF=0;
508  else if(genTitle.Contains("bele")) typeHF=3;
509  else if(genTitle.Contains("cele")) typeHF=2;
510  fHistNcollHFtype->Fill(typeHF,1.);
511  }
512  Int_t nParticles=mcEvent->GetNumberOfTracks();
513  Double_t dNchdy = 0.;
514  Int_t nb = 0, nc=0;
515  Int_t nCharmed=0;
516  Int_t nPhysPrim=0;
517  Int_t nPiKPeta09=0;
518  for (Int_t i=0;i<nParticles;i++){
519  TParticle* part = (TParticle*)mcEvent->Particle(i);
520  Int_t absPdg=TMath::Abs(part->GetPdgCode());
521  Int_t pdg=part->GetPdgCode();
522  if(absPdg==4) nc++;
523  if(absPdg==5) nb++;
524  if(mcEvent->IsPhysicalPrimary(i)){
525  Double_t eta=part->Eta();
526  fHistoEtaPhysPrim->Fill(eta);
527  if(absPdg==11) fHistEtaPhiPtGenEle->Fill(eta,part->Phi(),part->Pt());
528  else if(absPdg==211) fHistEtaPhiPtGenPi->Fill(eta,part->Phi(),part->Pt());
529  else if(absPdg==321) fHistEtaPhiPtGenK->Fill(eta,part->Phi(),part->Pt());
530  else if(absPdg==2212) fHistEtaPhiPtGenPro->Fill(eta,part->Phi(),part->Pt());
531 
532  if(TMath::Abs(eta)<0.5){
533  dNchdy+=0.6666; // 2/3 for the ratio charged/all
534  nPhysPrim++;
535  }
536  if(TMath::Abs(eta)<0.9){
537  fHistoPtPhysPrim->Fill(part->Pt());
538  if(absPdg==211 || absPdg==321 || absPdg==2212){
539  nPiKPeta09++;
540  }
541  }
542  }
543  Float_t rapid=-999.;
544  if (part->Energy() != TMath::Abs(part->Pz())){
545  rapid=0.5*TMath::Log((part->Energy()+part->Pz())/(part->Energy()-part->Pz()));
546  }
547  Int_t iPart=-1;
548  Int_t iType=0;
549  Int_t iSpecies=-1;
550  Int_t dummy[4];
551  if(absPdg==421){
552  iSpecies=0;
553  iType=AliVertexingHFUtils::CheckD0Decay(mcEvent,i,dummy);
554  if(iType>0) iPart=0;
555  }
556  else if(absPdg==411){
557  iSpecies=1;
558  iType=AliVertexingHFUtils::CheckDplusDecay(mcEvent,i,dummy);
559  if(iType<0){
560  Int_t iTypeKKpi=AliVertexingHFUtils::CheckDplusKKpiDecay(mcEvent,i,dummy);
561  if(iTypeKKpi>0) iType=3;
562  }
563  if(iType>0) iPart=1;
564  }
565  else if(absPdg==413){
566  iSpecies=2;
567  iType=AliVertexingHFUtils::CheckDstarDecay(mcEvent,i,dummy);
568  if(iType>0) iPart=2;
569  }
570  else if(absPdg==431){
571  iSpecies=3;
572  iType=AliVertexingHFUtils::CheckDsDecay(mcEvent,i,dummy);
573  if(iType==1 || iType==2) iPart=3;
574  }
575  else if(absPdg==4122){
576  iSpecies=4;
577  iType=AliVertexingHFUtils::CheckLcpKpiDecay(mcEvent,i,dummy);
578  if(iType<0){
579  Int_t iTypeV0=AliVertexingHFUtils::CheckLcV0bachelorDecay(mcEvent,i,dummy);
580  if(iTypeV0==1) iType=5;
581  if(iTypeV0==2) iType=6;
582  }
583  fHistLcDecayChan->Fill(iType);
584  if(iType>=0) iPart=4;
585  }
586  if(iSpecies>=0) fHistYPtAllDecay[iSpecies]->Fill(part->Pt(),rapid);
587 
588  // check beauty mesons
589  if(absPdg==511) fHistBYPtAllDecay[0]->Fill(part->Pt(),rapid);
590  else if(absPdg==521) fHistBYPtAllDecay[1]->Fill(part->Pt(),rapid);
591  else if(absPdg==513) fHistBYPtAllDecay[2]->Fill(part->Pt(),rapid);
592  else if(absPdg==531) fHistBYPtAllDecay[3]->Fill(part->Pt(),rapid);
593  else if(absPdg==5122) fHistBYPtAllDecay[4]->Fill(part->Pt(),rapid);
594 
595  if(pdg==511) fHistBSpecies->Fill(0);
596  else if(pdg==-511) fHistBSpecies->Fill(1);
597  else if(pdg==521) fHistBSpecies->Fill(2);
598  else if(pdg==-521) fHistBSpecies->Fill(3);
599  else if(pdg==513) fHistBSpecies->Fill(4);
600  else if(pdg==-513) fHistBSpecies->Fill(5);
601  else if(pdg==531) fHistBSpecies->Fill(6);
602  else if(pdg==-531) fHistBSpecies->Fill(7);
603  else if(pdg==5122) fHistBSpecies->Fill(8);
604  else if(pdg==-5122) fHistBSpecies->Fill(9);
605 
606  if(iSpecies<0) continue; // not a charm meson
607 
608  if(pdg==421) fHistDSpecies->Fill(0);
609  else if(pdg==-421) fHistDSpecies->Fill(1);
610  else if(pdg==411) fHistDSpecies->Fill(2);
611  else if(pdg==-411) fHistDSpecies->Fill(3);
612  else if(pdg==413) fHistDSpecies->Fill(4);
613  else if(pdg==-413) fHistDSpecies->Fill(5);
614  else if(pdg==431) fHistDSpecies->Fill(6);
615  else if(pdg==-431) fHistDSpecies->Fill(7);
616  else if(pdg==4122) fHistDSpecies->Fill(8);
617  else if(pdg==-4122) fHistDSpecies->Fill(9);
618 
619  Double_t distx=part->Vx()-mcVert->GetX();
620  Double_t disty=part->Vy()-mcVert->GetY();
621  Double_t distz=part->Vz()-mcVert->GetZ();
622  Double_t distToVert=TMath::Sqrt(distx*distx+disty*disty+distz*distz);
623  fHistMotherID->Fill(part->GetFirstMother());
625  if(iFromB==4){
626  fHistYPtPromptAllDecay[iSpecies]->Fill(part->Pt(),rapid);
627  fHistOriginPrompt->Fill(distToVert);
628  }
629  else if(iFromB==5){
630  fHistYPtFeeddownAllDecay[iSpecies]->Fill(part->Pt(),rapid);
631  fHistOriginFeeddown->Fill(distToVert);
632  }
633 
634  if(iPart<0) continue;
635  if(iType<0) continue;
636  nCharmed++;
637  if(iPart==0 && iType>0 && iType<=2){
638  fHistYPtD0byDecChannel[iType-1]->Fill(part->Pt(),rapid);
639  }else if(iPart==1 && iType>0 && iType<=3){
640  fHistYPtDplusbyDecChannel[iType-1]->Fill(part->Pt(),rapid);
641  }else if(iPart==3 && iType>0 && iType<=2){
642  fHistYPtDsbyDecChannel[iType-1]->Fill(part->Pt(),rapid);
643  }
644 
645  if(iFromB==4 && iPart>=0 && iPart<5) fHistYPtPrompt[iPart]->Fill(part->Pt(),rapid);
646  else if(iFromB==5 && iPart>=0 && iPart<5) fHistYPtFeeddown[iPart]->Fill(part->Pt(),rapid);
647  }
648 
649  for(Int_t i=0; i<nTracks; i++){
650  AliESDtrack* track=esd->GetTrack(i);
651  if(fESDtrackCuts->AcceptTrack(track)){
652  if(track->GetLabel()>0) fHistPtRecGood->Fill(track->Pt());
653  else fHistPtRecFake->Fill(track->Pt());
654  Int_t label=TMath::Abs(track->GetLabel());
655 
656  if(mcEvent->IsPhysicalPrimary(label)){
657  TParticle* part = (TParticle*)mcEvent->Particle(label);
658  Int_t absPdg=TMath::Abs(part->GetPdgCode());
659  if(absPdg==11) fHistEtaPhiPtRecEle->Fill(part->Eta(),part->Phi(),part->Pt());
660  else if(absPdg==211) fHistEtaPhiPtRecPi->Fill(part->Eta(),part->Phi(),part->Pt());
661  else if(absPdg==321) fHistEtaPhiPtRecK->Fill(part->Eta(),part->Phi(),part->Pt());
662  else if(absPdg==2212) fHistEtaPhiPtRecPro->Fill(part->Eta(),part->Phi(),part->Pt());
663  fHistPtRecVsPtGen->Fill(part->Pt(),track->Pt());
664  fHistPhiRecVsPhiGen->Fill(part->Phi(),track->Phi());
665  fHistEtaRecVsEtaGen->Fill(part->Eta(),track->Eta());
666  }
667  }
668  }
669  fHistoNcharmed->Fill(dNchdy,nCharmed);
670  fHistoNbVsNc->Fill(nc,nb);
671  fHistoPhysPrim->Fill(nPhysPrim);
672  fHistoPhysPrimPiKPi09->Fill(nPiKPeta09);
673  fHistoPhysPrimPiKPi09vsb->Fill(imppar,nPiKPeta09);
674  }
675 
676  PostData(1,fOutput);
677 
678 }
679 //______________________________________________________________________________
681 {
683  fOutput = dynamic_cast<TList*> (GetOutputData(1));
684  if (!fOutput) {
685  printf("ERROR: fOutput not available\n");
686  return;
687  }
688 
689  return;
690 }
691 
692 
693 
694 
TH3F * fHistEtaPhiPtRecK
! histo of generated kaons
Int_t pdg
static Int_t CheckD0Decay(AliMCEvent *mcEvent, Int_t label, Int_t *arrayDauLab)
AliESDtrackCuts * fESDtrackCuts
0=pp, 1=PbPb, 2=pPb
TH1F * fHistPtRecGood
! pt distribution of "good" tracks
TH1F * fHistBSpecies
! histo of B hadron species
TH1F * fHistoSPDZVtxY
! histo with distr. of y coord. of SPDZ vertex
double Double_t
Definition: External.C:58
Definition: External.C:260
TH2F * fHistYPtPrompt[5]
! histo of y vs. pt from prompt D0, D+, D*, Ds, Lc
Int_t fSystem
c/b separation using quarks
Definition: External.C:236
TH3F * fHistEtaPhiPtGenEle
! histo of generated electrons
const char * title
Definition: MakeQAPdf.C:27
static Int_t CheckDplusDecay(AliMCEvent *mcEvent, Int_t label, Int_t *arrayDauLab)
TH1F * fHistoSPDZVtxX
! histo with distr. of x coord. of SPDZ vertex
TH2F * fHistYPtDplusbyDecChannel[3]
! histo of y vs. pt for D+->Kpipi and D+->K0*pi
TH1F * fHistMotherID
! histo of mother ID
TH3F * fHistEtaPhiPtGenK
! histo of generated kaons
TH2F * fHistBYPtAllDecay[5]
! histo of y vs. pt from prompt B0, B+, B*, Bs, Lb
TH1F * fHistOriginFeeddown
! histo of D production point (feeddown)
static Int_t CheckDsDecay(AliMCEvent *mcEvent, Int_t label, Int_t *arrayDauLab)
static Int_t CheckLcpKpiDecay(AliMCEvent *mcEvent, Int_t label, Int_t *arrayDauLab)
TH1F * fHistoPhysPrimPiKPi09
! histo of n. of primary pi, K, p
TH2F * fHistEtaRecVsEtaGen
! correlation between rec and gen pt
TH1F * fHistoTracklets
! histo with number of SPD tracklets
static Int_t CheckOrigin(TClonesArray *arrayMC, AliAODMCParticle *mcPart, Bool_t searchUpToQuark=kTRUE)
Functions to check the decay tree.
virtual void UserExec(Option_t *option)
TH2F * fHistYPtDsbyDecChannel[2]
! histo of y vs. pt for Ds->phipi and Ds->K0*K
TH1F * fHistoSPD3DVtxX
! histo with distr. of x coord. of SPD3D vertex
TH2F * fHistoNbVsNc
! histo of n. b quarks vs. n c. quarks
TH1F * fHistoTRKVtxX
! histo with distr. of x coord. of TRK vertex
TH2F * fHistNcollHFtype
! histo of number of injected events vs. type
TH2F * fHistYPtFeeddown[5]
! histo of y vs. pt from feeddown D0, D+, D*, Ds, Lc
static Int_t CheckLcV0bachelorDecay(AliMCEvent *mcEvent, Int_t label, Int_t *arrayDauLab)
TH3F * fHistEtaPhiPtRecPro
! histo of generated protons
TH1F * fHistoSelTracks
! histo with number of TPC+ITS refit ESD tracks
TH1F * fHistoPhysPrim
! histo of n. of physical primaries in |eta|<0.5
TH1F * fHistoNEvents
! histo with N of events
int Int_t
Definition: External.C:63
TH3F * fHistEtaPhiPtGenPi
! histo of generated pions
unsigned int UInt_t
Definition: External.C:33
TH3F * fHistEtaPhiPtGenPro
! histo of generated protons
float Float_t
Definition: External.C:68
TH1F * fHistPtRecFake
! pt distribution of fake tracks
static Int_t CheckDplusKKpiDecay(AliMCEvent *mcEvent, Int_t label, Int_t *arrayDauLab)
TH1F * fHistOriginPrompt
! histo of D production point (prompt)
TH2F * fHistYPtPromptAllDecay[5]
! histo of y vs. pt from prompt D0, D+, D*, Ds, Lc, no selection on decay channel ...
TH2F * fHistoNcharmed
! histo of D mesons vs. dN/dy
TH2F * fHistoPhysPrimPiKPi09vsb
! histo of n. of primary pi, K, p vs. b
TH2F * fHistPhiRecVsPhiGen
! correlation between rec and gen pt
TH1F * fHistoTRKVtxZ
! histo with distr. of z coord. of TRK vertex
TH1F * fHistDSpecies
! histo of D hadron species
static Int_t CheckDstarDecay(AliMCEvent *mcEvent, Int_t label, Int_t *arrayDauLab)
TH2F * fHistYPtAllDecay[5]
! histo of y vs. pt from prompt D0, D+, D*, Ds, Lc, no selection on decay channel ...
TH1F * fHistoEtaPhysPrim
! histo of eta distribution of physical primaries
TH1F * fHistoSPD3DVtxY
! histo with distr. of y coord. of SPD3D vertex
TH1F * fHistoTrackletsEta1
! histo with number of SPD tracklets in |eta|<1
TH2F * fHistYPtD0byDecChannel[2]
! histo of y vs. pt for D0->Kpi and D0->Kpipipi
TH1F * fHistoTracks
! histo with number of ESD tracks
TH1F * fHistoPtPhysPrim
! histo of pt distribution of phys primaries
TH1F * fHistoSPDZVtxZ
! histo with distr. of z coord. of SPDZ vertex
TH1F * fHistoTRKVtxY
! histo with distr. of y coord. of TRK vertex
TH3F * fHistEtaPhiPtRecEle
! histo of generated electrons
TH2F * fHistYPtFeeddownAllDecay[5]
! histo of y vs. pt from prompt D0, D+, D*, Ds, Lc, no selection on decay channel ...
TList * fOutput
! list of output histos
const char Option_t
Definition: External.C:48
TH1F * fHistLcDecayChan
! histo of Lc decay modes
TH2F * fHistNinjectedvsb
! histo of number of injected events vs. b
TH1F * fHistoSPD3DVtxZ
! histo with distr. of z coord. of SPD3D vertex
TH3F * fHistEtaPhiPtRecPi
! histo of generated pions
Double_t maxMult
TH2F * fHistPtRecVsPtGen
! correlation between rec and gen pt
virtual void Terminate(Option_t *option)