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