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