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