]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGHF/vertexingHF/AliAnalysisTaskCheckHFMCProd.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
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"
7b6a4dcd 14#include "AliMultiplicity.h"
15#include <TParticle.h>
16#include <TSystem.h>
17#include <TTree.h>
18#include <TNtuple.h>
19#include <TH1F.h>
20#include <TH2F.h>
7c569cb6 21#include <TH3F.h>
7b6a4dcd 22#include <TChain.h>
23#include "AliESDInputHandlerRP.h"
24#include "AliAnalysisTaskCheckHFMCProd.h"
25
26/**************************************************************************
27 * Copyright(c) 1998-2012, ALICE Experiment at CERN, All rights reserved. *
28 * *
29 * Author: The ALICE Off-line Project. *
30 * Contributors are mentioned in the code where appropriate. *
31 * *
32 * Permission to use, copy, modify and distribute this software and its *
33 * documentation strictly for non-commercial purposes is hereby granted *
34 * without fee, provided that the above copyright notice appears in all *
35 * copies and that both the copyright notice and this permission notice *
36 * appear in the supporting documentation. The authors make no claims *
37 * about the suitability of this software for any purpose. It is *
38 * provided "as is" without express or implied warranty. *
39 **************************************************************************/
40
41/* $Id$ */
42
43//*************************************************************************
44// Implementation of class AliAnalysisTaskCheckHFMCProd
45// AliAnalysisTask to check MC production at ESD+Kine level
46//
47//
48// Authors: F. Prino, prino@to.infn.it
49//
50//*************************************************************************
51
52ClassImp(AliAnalysisTaskCheckHFMCProd)
53//______________________________________________________________________________
54AliAnalysisTaskCheckHFMCProd::AliAnalysisTaskCheckHFMCProd() : AliAnalysisTaskSE("HFMCChecks"),
55 fOutput(0),
56 fHistoNEvents(0),
a35b52e8 57 fHistoPhysPrim(0),
7b6a4dcd 58 fHistoTracks(0),
59 fHistoSelTracks(0),
60 fHistoTracklets(0),
7c3336d1 61 fHistoTrackletsEta1(0),
a916dfdb 62 fHistoPtPhysPrim(0),
c7112a2a 63 fHistoEtaPhysPrim(0),
7b6a4dcd 64 fHistoSPD3DVtxX(0),
65 fHistoSPD3DVtxY(0),
66 fHistoSPD3DVtxZ(0),
67 fHistoSPDZVtxX(0),
68 fHistoSPDZVtxY(0),
69 fHistoSPDZVtxZ(0),
70 fHistoTRKVtxX(0),
71 fHistoTRKVtxY(0),
72 fHistoTRKVtxZ(0),
73 fHistoNcharmed(0),
74 fHistoNbVsNc(0),
fcedd2b1 75 fHistOriginPrompt(0),
76 fHistOriginFeeddown(0),
c833b704 77 fHistMotherID(0),
a916dfdb 78 fHistDSpecies(0),
79 fHistBSpecies(0),
c7112a2a 80 fHistNcollHFtype(0),
7c569cb6 81 fHistEtaPhiPtGenEle(0),
82 fHistEtaPhiPtGenPi(0),
83 fHistEtaPhiPtGenK(0),
84 fHistEtaPhiPtGenPro(0),
85 fHistEtaPhiPtRecEle(0),
86 fHistEtaPhiPtRecPi(0),
87 fHistEtaPhiPtRecK(0),
88 fHistEtaPhiPtRecPro(0),
fcedd2b1 89 fSearchUpToQuark(kFALSE),
a35b52e8 90 fSystem(0),
7c569cb6 91 fESDtrackCuts(0x0),
7b6a4dcd 92 fReadMC(kTRUE)
93{
94 //
95 DefineInput(0, TChain::Class());
96 DefineOutput(1, TList::Class());
97}
98
99
100//___________________________________________________________________________
101AliAnalysisTaskCheckHFMCProd::~AliAnalysisTaskCheckHFMCProd(){
102 //
103 if (AliAnalysisManager::GetAnalysisManager()->IsProofMode()) return;
104 if (fOutput) {
105 delete fOutput;
106 fOutput = 0;
107 }
7c569cb6 108 delete fESDtrackCuts;
109
7b6a4dcd 110}
111
112//___________________________________________________________________________
113void AliAnalysisTaskCheckHFMCProd::UserCreateOutputObjects() {
114 // create output histos
115
116 fOutput = new TList();
117 fOutput->SetOwner();
118 fOutput->SetName("OutputHistos");
119
120 fHistoNEvents = new TH1F("hNEvents", "Number of processed events",3,-0.5,2.5);
121 fHistoNEvents->Sumw2();
122 fHistoNEvents->SetMinimum(0);
123 fOutput->Add(fHistoNEvents);
124
125 Double_t maxMult=100.;
a35b52e8 126 if(fSystem==1) maxMult=10000.;
127 if(fSystem==2) maxMult=500.;
7c3336d1 128 fHistoPhysPrim = new TH1F("hPhysPrim","",100,-0.5,maxMult-0.5);
a35b52e8 129 fHistoPhysPrim->Sumw2();
130 fOutput->Add(fHistoPhysPrim);
7c3336d1 131 fHistoTracks = new TH1F("hTracks","",100,-0.5,maxMult*2-0.5);
7b6a4dcd 132 fHistoTracks->Sumw2();
133 fOutput->Add(fHistoTracks);
7c3336d1 134 fHistoSelTracks = new TH1F("hSelTracks","",100,-0.5,maxMult-0.5);
7b6a4dcd 135 fHistoSelTracks->Sumw2();
136 fOutput->Add(fHistoSelTracks);
7c3336d1 137 fHistoTracklets = new TH1F("hTracklets","",100,-0.5,maxMult-0.5);
7b6a4dcd 138 fHistoTracklets->Sumw2();
139 fOutput->Add(fHistoTracklets);
7c3336d1 140 fHistoTrackletsEta1 = new TH1F("hTrackletsEta1","",100,-0.5,maxMult-0.5);
141 fHistoTrackletsEta1->Sumw2();
142 fOutput->Add(fHistoTrackletsEta1);
a916dfdb 143 fHistoPtPhysPrim = new TH1F("hPtPhysPrim","",100,0.,20.);
144 fHistoPtPhysPrim->Sumw2();
145 fOutput->Add(fHistoPtPhysPrim);
c7112a2a 146 fHistoEtaPhysPrim = new TH1F("hEtaPhysPrim","",100,-10.,10.);
147 fHistoEtaPhysPrim->Sumw2();
148 fOutput->Add(fHistoEtaPhysPrim);
7b6a4dcd 149
150 fHistoSPD3DVtxX = new TH1F("hSPD3DvX","",100,-1.,1.);
151 fHistoSPD3DVtxX->Sumw2();
152 fOutput->Add(fHistoSPD3DVtxX);
153 fHistoSPD3DVtxY = new TH1F("hSPD3DvY","",100,-1.,1.);
154 fHistoSPD3DVtxY->Sumw2();
155 fOutput->Add(fHistoSPD3DVtxY);
156 fHistoSPD3DVtxZ = new TH1F("hSPD3DvZ","",100,-15.,15.);
157 fHistoSPD3DVtxZ->Sumw2();
158 fOutput->Add(fHistoSPD3DVtxZ);
159
160 fHistoSPDZVtxX = new TH1F("hSPDZvX","",100,-1.,1.);
161 fHistoSPDZVtxX->Sumw2();
162 fOutput->Add(fHistoSPDZVtxX);
163 fHistoSPDZVtxY = new TH1F("hSPDZvY","",100,-1.,1.);
164 fHistoSPDZVtxY->Sumw2();
165 fOutput->Add(fHistoSPDZVtxY);
166 fHistoSPDZVtxZ = new TH1F("hSPDZvZ","",100,-15.,15.);
167 fHistoSPDZVtxZ->Sumw2();
168 fOutput->Add(fHistoSPDZVtxZ);
169
170
171 fHistoTRKVtxX = new TH1F("hTRKvX","",100,-1.,1.);
172 fHistoTRKVtxX->Sumw2();
173 fOutput->Add(fHistoTRKVtxX);
174 fHistoTRKVtxY = new TH1F("hTRKvY","",100,-1.,1.);
175 fHistoTRKVtxY->Sumw2();
176 fOutput->Add(fHistoTRKVtxY);
177 fHistoTRKVtxZ = new TH1F("hTRKvZ","",100,-15.,15.);
178 fHistoTRKVtxZ->Sumw2();
179 fOutput->Add(fHistoTRKVtxZ);
180
181 Int_t nBinscb=11;
a916dfdb 182 if(fSystem==1) nBinscb=200;
a35b52e8 183 if(fSystem==2) nBinscb=21;
7b6a4dcd 184 Double_t maxncn=nBinscb-0.5;
7c3336d1 185 fHistoNcharmed = new TH2F("hncharmed","",100,-0.5,maxMult-0.5,nBinscb,-0.5,maxncn);
7b6a4dcd 186 fHistoNcharmed->Sumw2();
187 fOutput->Add(fHistoNcharmed);
188 fHistoNbVsNc = new TH2F("hnbvsnc","",nBinscb,-0.5,maxncn,nBinscb,-0.5,maxncn);
189 fHistoNbVsNc->Sumw2();
190 fOutput->Add(fHistoNbVsNc);
191
a35b52e8 192 fHistYPtPrompt[0] = new TH2F("hyptD0prompt","D0 - Prompt",40,0.,40.,20,-2.,2.);
193 fHistYPtPrompt[1] = new TH2F("hyptDplusprompt","Dplus - Prompt",40,0.,40.,20,-2.,2.);
194 fHistYPtPrompt[2] = new TH2F("hyptDstarprompt","Dstar - Prompt",40,0.,40.,20,-2.,2.);
195 fHistYPtPrompt[3] = new TH2F("hyptDsprompt","Ds - Prompt",40,0.,40.,20,-2.,2.);
196 fHistYPtPrompt[4] = new TH2F("hyptLcprompt","Lc - Prompt",40,0.,40.,20,-2.,2.);
197
198 fHistBYPtAllDecay[0] = new TH2F("hyptB0AllDecay","B0 - All",40,0.,40.,40,-2.,2.);
199 fHistBYPtAllDecay[1] = new TH2F("hyptBplusAllDecay","Bplus - All",40,0.,40.,40,-2.,2.);
200 fHistBYPtAllDecay[2] = new TH2F("hyptBstarAllDecay","Bstar - All",40,0.,40.,40,-2.,2.);
201 fHistBYPtAllDecay[3] = new TH2F("hyptBsAllDecay","Bs - All",40,0.,40.,40,-2.,2.);
202 fHistBYPtAllDecay[4] = new TH2F("hyptLbAllDecay","LB - All",40,0.,40.,40,-2.,2.);
203
204 fHistYPtAllDecay[0] = new TH2F("hyptD0AllDecay","D0 - All",40,0.,40.,40,-2.,2.);
205 fHistYPtAllDecay[1] = new TH2F("hyptDplusAllDecay","Dplus - All",40,0.,40.,40,-2.,2.);
206 fHistYPtAllDecay[2] = new TH2F("hyptDstarAllDecay","Dstar - All",40,0.,40.,40,-2.,2.);
207 fHistYPtAllDecay[3] = new TH2F("hyptDsAllDecay","Ds - All",40,0.,40.,40,-2.,2.);
208 fHistYPtAllDecay[4] = new TH2F("hyptLcAllDecay","Lc - All",40,0.,40.,40,-2.,2.);
209
210 fHistYPtPromptAllDecay[0] = new TH2F("hyptD0promptAllDecay","D0 - Prompt",40,0.,40.,40,-2.,2.);
211 fHistYPtPromptAllDecay[1] = new TH2F("hyptDpluspromptAllDecay","Dplus - Prompt",40,0.,40.,40,-2.,2.);
212 fHistYPtPromptAllDecay[2] = new TH2F("hyptDstarpromptAllDecay","Dstar - Prompt",40,0.,40.,40,-2.,2.);
213 fHistYPtPromptAllDecay[3] = new TH2F("hyptDspromptAllDecay","Ds - Prompt",40,0.,40.,40,-2.,2.);
214 fHistYPtPromptAllDecay[4] = new TH2F("hyptLcpromptAllDecay","Lc - Prompt",40,0.,40.,40,-2.,2.);
215
216 fHistYPtFeeddownAllDecay[0] = new TH2F("hyptD0feeddownAllDecay","D0 - FromB",40,0.,40.,40,-2.,2.);
217 fHistYPtFeeddownAllDecay[1] = new TH2F("hyptDplusfeeddownAllDecay","Dplus - FromB",40,0.,40.,40,-2.,2.);
218 fHistYPtFeeddownAllDecay[2] = new TH2F("hyptDstarfeeddownAllDecay","Dstar - FromB",40,0.,40.,40,-2.,2.);
219 fHistYPtFeeddownAllDecay[3] = new TH2F("hyptDsfeeddownAllDecay","Ds - FromB",40,0.,40.,40,-2.,2.);
220 fHistYPtFeeddownAllDecay[4] = new TH2F("hyptLcfeeddownAllDecay","Lc - FromB",40,0.,40.,40,-2.,2.);
221
222
223 fHistYPtFeeddown[0] = new TH2F("hyptD0feeddown","D0 - Feeddown",40,0.,40.,20,-2.,2.);
224 fHistYPtFeeddown[1] = new TH2F("hyptDplusfeeddown","Dplus - Feeddown",40,0.,40.,20,-2.,2.);
225 fHistYPtFeeddown[2] = new TH2F("hyptDstarfeedown","Dstar - Feeddown",40,0.,40.,20,-2.,2.);
226 fHistYPtFeeddown[3] = new TH2F("hyptDsfeedown","Ds - Feeddown",40,0.,40.,20,-2.,2.);
227 fHistYPtFeeddown[4] = new TH2F("hyptLcfeedown","Lc - Feeddown",40,0.,40.,20,-2.,2.);
7b6a4dcd 228
229 for(Int_t ih=0; ih<5; ih++){
a35b52e8 230 fHistBYPtAllDecay[ih]->Sumw2();
231 fHistBYPtAllDecay[ih]->SetMinimum(0);
232 fOutput->Add(fHistBYPtAllDecay[ih]);
5472ae5b 233 fHistYPtAllDecay[ih]->Sumw2();
234 fHistYPtAllDecay[ih]->SetMinimum(0);
235 fOutput->Add(fHistYPtAllDecay[ih]);
236 fHistYPtPromptAllDecay[ih]->Sumw2();
237 fHistYPtPromptAllDecay[ih]->SetMinimum(0);
238 fOutput->Add(fHistYPtPromptAllDecay[ih]);
239 fHistYPtFeeddownAllDecay[ih]->Sumw2();
240 fHistYPtFeeddownAllDecay[ih]->SetMinimum(0);
241 fOutput->Add(fHistYPtFeeddownAllDecay[ih]);
7b6a4dcd 242 fHistYPtPrompt[ih]->Sumw2();
243 fHistYPtPrompt[ih]->SetMinimum(0);
244 fOutput->Add(fHistYPtPrompt[ih]);
245 fHistYPtFeeddown[ih]->Sumw2();
246 fHistYPtFeeddown[ih]->SetMinimum(0);
247 fOutput->Add(fHistYPtFeeddown[ih]);
248 }
249
a35b52e8 250 fHistYPtD0byDecChannel[0] = new TH2F("hyptD02","D0 - 2prong",40,0.,40.,20,-2.,2.);
251 fHistYPtD0byDecChannel[1] = new TH2F("hyptD04","D0 - 4prong",40,0.,40.,20,-2.,2.);
252 fHistYPtDplusbyDecChannel[0] = new TH2F("hyptDplusnonreson","Dplus - non reson",40,0.,40.,20,-2.,2.);
253 fHistYPtDplusbyDecChannel[1] = new TH2F("hyptDplusreson","Dplus - reson via K0*",40,0.,40.,20,-2.,2.);
254 fHistYPtDsbyDecChannel[0] = new TH2F("hyptDsphi","Ds - vis Phi",40,0.,40.,20,-2.,2.);
255 fHistYPtDsbyDecChannel[1] = new TH2F("hyptDsk0st","Ds - via k0*",40,0.,40.,20,-2.,2.);
7b6a4dcd 256
257 for(Int_t ih=0; ih<2; ih++){
258
259 fHistYPtD0byDecChannel[ih]->Sumw2();
260 fHistYPtD0byDecChannel[ih]->SetMinimum(0);
261 fOutput->Add(fHistYPtD0byDecChannel[ih]);
262 fHistYPtDplusbyDecChannel[ih]->Sumw2();
263 fHistYPtDplusbyDecChannel[ih]->SetMinimum(0);
264 fOutput->Add(fHistYPtDplusbyDecChannel[ih]);
265 fHistYPtDsbyDecChannel[ih]->Sumw2();
266 fHistYPtDsbyDecChannel[ih]->SetMinimum(0);
267 fOutput->Add(fHistYPtDsbyDecChannel[ih]);
268 }
269
fcedd2b1 270 fHistOriginPrompt=new TH1F("hOriginPrompt","",100,0.,0.5);
271 fHistOriginPrompt->Sumw2();
272 fHistOriginPrompt->SetMinimum(0);
273 fOutput->Add(fHistOriginPrompt);
274 fHistOriginFeeddown=new TH1F("hOriginFeeddown","",100,0.,0.5);
275 fHistOriginFeeddown->Sumw2();
276 fHistOriginFeeddown->SetMinimum(0);
277 fOutput->Add(fHistOriginFeeddown);
c833b704 278 fHistMotherID=new TH1F("hMotherID","",1000,-1.5,998.5);
279 fHistMotherID->SetMinimum(0);
280 fOutput->Add(fHistMotherID);
a916dfdb 281 fHistDSpecies=new TH1F("hDSpecies","",10,-0.5,9.5);
282 fHistDSpecies->GetXaxis()->SetBinLabel(1,"D0");
283 fHistDSpecies->GetXaxis()->SetBinLabel(2,"D0bar");
284 fHistDSpecies->GetXaxis()->SetBinLabel(3,"D+");
285 fHistDSpecies->GetXaxis()->SetBinLabel(4,"D-");
286 fHistDSpecies->GetXaxis()->SetBinLabel(5,"D*+");
287 fHistDSpecies->GetXaxis()->SetBinLabel(6,"D*-");
288 fHistDSpecies->GetXaxis()->SetBinLabel(7,"Ds+");
289 fHistDSpecies->GetXaxis()->SetBinLabel(8,"Ds-");
290 fHistDSpecies->GetXaxis()->SetBinLabel(9,"Lc+");
291 fHistDSpecies->GetXaxis()->SetBinLabel(10,"Lc-");
292 fHistDSpecies->SetMinimum(0);
293 fOutput->Add(fHistDSpecies);
294 fHistBSpecies=new TH1F("hBSpecies","",10,-0.5,9.5);
295 fHistBSpecies->GetXaxis()->SetBinLabel(1,"B0");
296 fHistBSpecies->GetXaxis()->SetBinLabel(2,"B0bar");
297 fHistBSpecies->GetXaxis()->SetBinLabel(3,"B+");
298 fHistBSpecies->GetXaxis()->SetBinLabel(4,"B-");
299 fHistBSpecies->GetXaxis()->SetBinLabel(5,"B*+");
300 fHistBSpecies->GetXaxis()->SetBinLabel(6,"B*-");
301 fHistBSpecies->GetXaxis()->SetBinLabel(7,"Bs+");
302 fHistBSpecies->GetXaxis()->SetBinLabel(8,"Bs-");
303 fHistBSpecies->GetXaxis()->SetBinLabel(9,"Lb+");
304 fHistBSpecies->GetXaxis()->SetBinLabel(10,"Lb-");
305 fHistBSpecies->SetMinimum(0);
306 fOutput->Add(fHistBSpecies);
c7112a2a 307
308 fHistNcollHFtype=new TH2F("hNcollHFtype","",5,-1.5,3.5,30,-0.5,29.5);
309 fOutput->Add(fHistNcollHFtype);
7c569cb6 310
311 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};
312 const Int_t nBinsPhi=40;
313 Double_t binsphi[nBinsPhi+1];
314 for(Int_t ib=0; ib<=nBinsPhi; ib++) binsphi[ib]=ib*TMath::Pi()/20.;
315 const Int_t nBinsPt=24;
316 Double_t binspt[nBinsPt+1]={0.,0.10,0.15,0.2,0.25,
317 0.3,0.4,0.5,0.6,0.7,
318 0.8,0.9,1.,1.25,1.5,
319 1.75,2.,2.5,3.,4.,
320 5.,7.5,10.,15.,20.};
321
322 fHistEtaPhiPtGenEle=new TH3F("hEtaPhiPtGenEle","",10,binseta,nBinsPhi,binsphi,nBinsPt,binspt);
323 fOutput->Add(fHistEtaPhiPtGenEle);
324 fHistEtaPhiPtGenPi=new TH3F("hEtaPhiPtGenPi","",10,binseta,nBinsPhi,binsphi,nBinsPt,binspt);
325 fOutput->Add(fHistEtaPhiPtGenPi);
326 fHistEtaPhiPtGenK=new TH3F("hEtaPhiPtGenK","",10,binseta,nBinsPhi,binsphi,nBinsPt,binspt);
327 fOutput->Add(fHistEtaPhiPtGenK);
328 fHistEtaPhiPtGenPro=new TH3F("hEtaPhiPtGenPro","",10,binseta,nBinsPhi,binsphi,nBinsPt,binspt);
329 fOutput->Add(fHistEtaPhiPtGenPro);
330
331
332 fHistEtaPhiPtRecEle=new TH3F("hEtaPhiPtRecEle","",10,binseta,nBinsPhi,binsphi,nBinsPt,binspt);
333 fOutput->Add(fHistEtaPhiPtRecEle);
334 fHistEtaPhiPtRecPi=new TH3F("hEtaPhiPtRecPi","",10,binseta,nBinsPhi,binsphi,nBinsPt,binspt);
335 fOutput->Add(fHistEtaPhiPtRecPi);
336 fHistEtaPhiPtRecK=new TH3F("hEtaPhiPtRecK","",10,binseta,nBinsPhi,binsphi,nBinsPt,binspt);
337 fOutput->Add(fHistEtaPhiPtRecK);
338 fHistEtaPhiPtRecPro=new TH3F("hEtaPhiPtRecPro","",10,binseta,nBinsPhi,binsphi,nBinsPt,binspt);
339 fOutput->Add(fHistEtaPhiPtRecPro);
340
341
342
7b6a4dcd 343 PostData(1,fOutput);
344
345}
346//______________________________________________________________________________
347void AliAnalysisTaskCheckHFMCProd::UserExec(Option_t *)
348{
349 //
350
351 AliESDEvent *esd = (AliESDEvent*) (InputEvent());
352
353
354 if(!esd) {
355 printf("AliAnalysisTaskSDDRP::Exec(): bad ESD\n");
356 return;
357 }
358
359 fHistoNEvents->Fill(0);
360
7c569cb6 361 if(!fESDtrackCuts){
362 Int_t year=2011;
363 if(esd->GetRunNumber()<=139517) year=2010;
364 if(year==2010) fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kFALSE);
365 else fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kFALSE);
366 fESDtrackCuts->SetMaxDCAToVertexXY(2.4);
367 fESDtrackCuts->SetMaxDCAToVertexZ(3.2);
368 fESDtrackCuts->SetDCAToVertex2D(kTRUE);
369 fESDtrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
370 AliESDtrackCuts::kAny);
371 }
372
7b6a4dcd 373 Int_t nTracks=esd->GetNumberOfTracks();
374 fHistoTracks->Fill(nTracks);
375 Int_t nSelTracks=0;
376 for(Int_t it=0; it<nTracks; it++){
377 AliESDtrack* tr=esd->GetTrack(it);
378 UInt_t status=tr->GetStatus();
379 if(!(status&AliESDtrack::kITSrefit)) continue;
380 if(!(status&AliESDtrack::kTPCin)) continue;
381 nSelTracks++;
382 }
383 fHistoSelTracks->Fill(nSelTracks);
384
385 const AliMultiplicity* mult=esd->GetMultiplicity();
386 Int_t nTracklets=mult->GetNumberOfTracklets();
7c3336d1 387 Int_t nTracklets1=0;
388 for(Int_t it=0; it<nTracklets; it++){
389 Double_t eta=TMath::Abs(mult->GetEta(it));
390 if(eta<1) nTracklets1++;
391 }
7b6a4dcd 392 fHistoTracklets->Fill(nTracklets);
7c3336d1 393 fHistoTrackletsEta1->Fill(nTracklets1);
394
7b6a4dcd 395 const AliESDVertex *spdv=esd->GetVertex();
396 if(spdv && spdv->IsFromVertexer3D()){
397 fHistoSPD3DVtxX->Fill(spdv->GetXv());
398 fHistoSPD3DVtxY->Fill(spdv->GetYv());
399 fHistoSPD3DVtxZ->Fill(spdv->GetZv());
400 }
401 if(spdv && spdv->IsFromVertexerZ()){
402 fHistoSPDZVtxX->Fill(spdv->GetXv());
403 fHistoSPDZVtxY->Fill(spdv->GetYv());
404 fHistoSPDZVtxZ->Fill(spdv->GetZv());
405 }
406 const AliESDVertex *trkv=esd->GetPrimaryVertex();
407 if(trkv && trkv->GetNContributors()>1){
408 fHistoTRKVtxX->Fill(trkv->GetXv());
409 fHistoTRKVtxY->Fill(trkv->GetYv());
410 fHistoTRKVtxZ->Fill(trkv->GetZv());
411 }
412
413 AliStack* stack=0;
7b6a4dcd 414 if(fReadMC){
415 AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
416 if (!eventHandler) {
417 Printf("ERROR: Could not retrieve MC event handler");
418 return;
419 }
420 AliMCEvent* mcEvent = eventHandler->MCEvent();
421 if (!mcEvent) {
422 Printf("ERROR: Could not retrieve MC event");
423 return;
424 }
425 stack = mcEvent->Stack();
426 if (!stack) {
427 Printf("ERROR: stack not available");
428 return;
429 }
fcedd2b1 430 const AliVVertex* mcVert=mcEvent->GetPrimaryVertex();
431 if(!mcVert){
432 Printf("ERROR: generated vertex not available");
433 return;
434 }
7c569cb6 435 if(TMath::Abs(mcVert->GetZ())>10) return;
436
c7112a2a 437 // const AliHeader* h=(AliHeader*)mcEvent->GetHeader();
438 // cout<<h<<endl;
439 TString genname=mcEvent->GenEventHeader()->ClassName();
440 Int_t nColl=-1;
441 Int_t typeHF=-1;
442 TList* lgen=0x0;
443 if(genname.Contains("CocktailEventHeader")){
444 AliGenCocktailEventHeader *cockhead=(AliGenCocktailEventHeader*)mcEvent->GenEventHeader();
445 lgen=cockhead->GetHeaders();
446 for(Int_t ig=0; ig<lgen->GetEntries(); ig++){
447 AliGenerator* gen=(AliGenerator*)lgen->At(ig);
448 TString title=gen->GetName();
449 if(title.Contains("bchadr")) typeHF=1;
450 else if(title.Contains("chadr")) typeHF=0;
451 else if(title.Contains("bele")) typeHF=3;
452 else if(title.Contains("cele")) typeHF=2;
453 }
454 nColl=lgen->GetEntries();
c7112a2a 455 fHistNcollHFtype->Fill(typeHF,nColl);
456 }
7b6a4dcd 457 Int_t nParticles=stack->GetNtrack();
458 Double_t dNchdy = 0.;
459 Int_t nb = 0, nc=0;
a0dabf3d 460 Int_t nCharmed=0;
a35b52e8 461 Int_t nPhysPrim=0;
7b6a4dcd 462 for (Int_t i=0;i<nParticles;i++){
463 TParticle* part = (TParticle*)stack->Particle(i);
464 Int_t absPdg=TMath::Abs(part->GetPdgCode());
a916dfdb 465 Int_t pdg=part->GetPdgCode();
7b6a4dcd 466 if(absPdg==4) nc++;
467 if(absPdg==5) nb++;
468 if(stack->IsPhysicalPrimary(i)){
469 Double_t eta=part->Eta();
c7112a2a 470 fHistoEtaPhysPrim->Fill(eta);
7c569cb6 471 if(absPdg==11) fHistEtaPhiPtGenEle->Fill(eta,part->Phi(),part->Pt());
472 else if(absPdg==211) fHistEtaPhiPtGenPi->Fill(eta,part->Phi(),part->Pt());
473 else if(absPdg==321) fHistEtaPhiPtGenK->Fill(eta,part->Phi(),part->Pt());
474 else if(absPdg==2212) fHistEtaPhiPtGenPro->Fill(eta,part->Phi(),part->Pt());
475
a35b52e8 476 if(TMath::Abs(eta)<0.5){
477 dNchdy+=0.6666; // 2/3 for the ratio charged/all
478 nPhysPrim++;
479 }
a916dfdb 480 if(TMath::Abs(eta)<0.9){
481 fHistoPtPhysPrim->Fill(part->Pt());
482 }
7b6a4dcd 483 }
5f7f126e 484 Float_t rapid=-999.;
485 if (part->Energy() != TMath::Abs(part->Pz())){
486 rapid=0.5*TMath::Log((part->Energy()+part->Pz())/(part->Energy()-part->Pz()));
487 }
5472ae5b 488
489
7b6a4dcd 490 Int_t iPart=-1;
491 Int_t iType=0;
5472ae5b 492 Int_t iSpecies=-1;
493 if(absPdg==421){
494 iSpecies=0;
7b6a4dcd 495 iType=CheckD0Decay(i,stack);
496 if(iType>=0) iPart=0;
497 }
498 else if(absPdg==411){
5472ae5b 499 iSpecies=1;
7b6a4dcd 500 iType=CheckDplusDecay(i,stack);
501 if(iType>=0) iPart=1;
502 }
503 else if(absPdg==413){
5472ae5b 504 iSpecies=2;
7b6a4dcd 505 iType=CheckDstarDecay(i,stack);
506 if(iType>=0) iPart=2;
507 }
508 else if(absPdg==431){
5472ae5b 509 iSpecies=3;
7b6a4dcd 510 iType=CheckDsDecay(i,stack);
511 if(iType==0 || iType==1) iPart=3;
512 }
513 else if(absPdg==4122){
5472ae5b 514 iSpecies=4;
7b6a4dcd 515 iType=CheckLcDecay(i,stack);
516 if(iType>=0) iPart=4;
517 }
a35b52e8 518 if(iSpecies>=0) fHistYPtAllDecay[iSpecies]->Fill(part->Pt(),rapid);
519
520 // check beauty mesons
521 if(absPdg==511) fHistBYPtAllDecay[0]->Fill(part->Pt(),rapid);
522 else if(absPdg==521) fHistBYPtAllDecay[1]->Fill(part->Pt(),rapid);
523 else if(absPdg==513) fHistBYPtAllDecay[2]->Fill(part->Pt(),rapid);
524 else if(absPdg==531) fHistBYPtAllDecay[3]->Fill(part->Pt(),rapid);
525 else if(absPdg==5122) fHistBYPtAllDecay[4]->Fill(part->Pt(),rapid);
526
a916dfdb 527 if(pdg==511) fHistBSpecies->Fill(0);
528 else if(pdg==-511) fHistBSpecies->Fill(1);
529 else if(pdg==521) fHistBSpecies->Fill(2);
530 else if(pdg==-521) fHistBSpecies->Fill(3);
531 else if(pdg==513) fHistBSpecies->Fill(4);
532 else if(pdg==-513) fHistBSpecies->Fill(5);
533 else if(pdg==531) fHistBSpecies->Fill(6);
534 else if(pdg==-531) fHistBSpecies->Fill(7);
535 else if(pdg==5122) fHistBSpecies->Fill(8);
536 else if(pdg==-5122) fHistBSpecies->Fill(9);
537
538 if(iSpecies<0) continue; // not a charm meson
539
540 if(pdg==421) fHistDSpecies->Fill(0);
541 else if(pdg==-421) fHistDSpecies->Fill(1);
542 else if(pdg==411) fHistDSpecies->Fill(2);
543 else if(pdg==-411) fHistDSpecies->Fill(3);
544 else if(pdg==413) fHistDSpecies->Fill(4);
545 else if(pdg==-413) fHistDSpecies->Fill(5);
546 else if(pdg==431) fHistDSpecies->Fill(6);
547 else if(pdg==-431) fHistDSpecies->Fill(7);
548 else if(pdg==4122) fHistDSpecies->Fill(8);
549 else if(pdg==-4122) fHistDSpecies->Fill(9);
5472ae5b 550
fcedd2b1 551 Double_t distx=part->Vx()-mcVert->GetX();
552 Double_t disty=part->Vy()-mcVert->GetY();
553 Double_t distz=part->Vz()-mcVert->GetZ();
554 Double_t distToVert=TMath::Sqrt(distx*distx+disty*disty+distz*distz);
c833b704 555 fHistMotherID->Fill(part->GetFirstMother());
7b6a4dcd 556 TParticle* runningpart=part;
557 Int_t iFromB=-1;
5472ae5b 558 Int_t pdgmoth=-1;
fcedd2b1 559 if(fSearchUpToQuark){
560 while(1){
561 Int_t labmoth=runningpart->GetFirstMother();
562 if(labmoth==-1) break;
563 TParticle *mot=(TParticle*)stack->Particle(labmoth);
564 pdgmoth=TMath::Abs(mot->GetPdgCode());
565 if(pdgmoth==5){
566 iFromB=1;
567 break;
568 }else if(pdgmoth==4){
569 iFromB=0;
570 break;
571 }
572 runningpart=mot;
573 }
574 }else{
575 iFromB=0;
576 while(1){
577 Int_t labmoth=runningpart->GetFirstMother();
578 if(labmoth==-1) break;
579 TParticle *mot=(TParticle*)stack->Particle(labmoth);
580 pdgmoth=TMath::Abs(mot->GetPdgCode());
581 if(pdgmoth>=500 && pdgmoth<=599){
582 iFromB=1;
583 break;
584 }
585 if(pdgmoth>=5000 && pdgmoth<=5999){
586 iFromB=1;
587 break;
588 }
589 runningpart=mot;
7b6a4dcd 590 }
7b6a4dcd 591 }
fcedd2b1 592
593 if(iFromB==0){
594 fHistYPtPromptAllDecay[iSpecies]->Fill(part->Pt(),rapid);
595 fHistOriginPrompt->Fill(distToVert);
596 }
597 else if(iFromB==1){
598 fHistYPtFeeddownAllDecay[iSpecies]->Fill(part->Pt(),rapid);
599 fHistOriginFeeddown->Fill(distToVert);
600 }
5472ae5b 601
602 if(iPart<0) continue;
603 if(iType<0) continue;
604 nCharmed++;
605 if(iPart==0 && iType<=1){
606 fHistYPtD0byDecChannel[iType]->Fill(part->Pt(),rapid);
607 }else if(iPart==1 && iType<=1){
608 fHistYPtDplusbyDecChannel[iType]->Fill(part->Pt(),rapid);
609 }else if(iPart==3 && iType<=1){
610 fHistYPtDsbyDecChannel[iType]->Fill(part->Pt(),rapid);
611 }
612
7b6a4dcd 613 if(iFromB==0 && iPart>=0 && iPart<5) fHistYPtPrompt[iPart]->Fill(part->Pt(),rapid);
fcedd2b1 614 else if(iFromB==1 && iPart>=0 && iPart<5) fHistYPtFeeddown[iPart]->Fill(part->Pt(),rapid);
7b6a4dcd 615 }
fcedd2b1 616
7c569cb6 617 for(Int_t i=0; i<nTracks; i++){
618 AliESDtrack* track=esd->GetTrack(i);
619 if(fESDtrackCuts->AcceptTrack(track)){
620 Int_t label=TMath::Abs(track->GetLabel());
621 if(stack->IsPhysicalPrimary(label)){
622 TParticle* part = (TParticle*)stack->Particle(label);
623 Int_t absPdg=TMath::Abs(part->GetPdgCode());
624 Double_t eta=part->Eta();
625 if(absPdg==11) fHistEtaPhiPtRecEle->Fill(eta,part->Phi(),part->Pt());
626 else if(absPdg==211) fHistEtaPhiPtRecPi->Fill(eta,part->Phi(),part->Pt());
627 else if(absPdg==321) fHistEtaPhiPtRecK->Fill(eta,part->Phi(),part->Pt());
628 else if(absPdg==2212) fHistEtaPhiPtRecPro->Fill(eta,part->Phi(),part->Pt());
629 }
630 }
631 }
7b6a4dcd 632 fHistoNcharmed->Fill(dNchdy,nCharmed);
633 fHistoNbVsNc->Fill(nc,nb);
a35b52e8 634 fHistoPhysPrim->Fill(nPhysPrim);
7b6a4dcd 635 }
636
637 PostData(1,fOutput);
638
639}
640//______________________________________________________________________________
641void AliAnalysisTaskCheckHFMCProd::Terminate(Option_t */*option*/)
642{
643 // Terminate analysis
644 fOutput = dynamic_cast<TList*> (GetOutputData(1));
645 if (!fOutput) {
646 printf("ERROR: fOutput not available\n");
647 return;
648 }
649
650 return;
651}
652
653
654
655
656//______________________________________________________________________________
657Int_t AliAnalysisTaskCheckHFMCProd::CheckD0Decay(Int_t labD0, AliStack* stack) const{
658
659 if(labD0<0) return -1;
660 TParticle* dp = (TParticle*)stack->Particle(labD0);
661 Int_t pdgdp=dp->GetPdgCode();
662 Int_t nDau=dp->GetNDaughters();
7b6a4dcd 663
664 if(nDau==2){
665 Int_t nKaons=0;
666 Int_t nPions=0;
667 for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
668 if(iDau<0) return -1;
669 TParticle* dau=(TParticle*)stack->Particle(iDau);
670 Int_t pdgdau=dau->GetPdgCode();
671 if(TMath::Abs(pdgdau)==321){
672 if(pdgdp>0 && pdgdau>0) return -1;
673 if(pdgdp<0 && pdgdau<0) return -1;
674 nKaons++;
675 }else if(TMath::Abs(pdgdau)==211){
676 if(pdgdp<0 && pdgdau>0) return -1;
677 if(pdgdp>0 && pdgdau<0) return -1;
678 nPions++;
679 }else{
680 return -1;
681 }
682 }
683 if(nPions!=1) return -1;
684 if(nKaons!=1) return -1;
685 for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
686 if(iDau<0) return -1;
687 }
688 return 0;
689 }
690
691 if(nDau==3 || nDau==4){
692 Int_t nKaons=0;
693 Int_t nPions=0;
694 for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
695 if(iDau<0) return -1;
696 TParticle* dau=(TParticle*)stack->Particle(iDau);
697 Int_t pdgdau=dau->GetPdgCode();
698 if(TMath::Abs(pdgdau)==321){
699 if(pdgdp>0 && pdgdau>0) return -1;
700 if(pdgdp<0 && pdgdau<0) return -1;
701 nKaons++;
702 }else if(TMath::Abs(pdgdau)==113 || TMath::Abs(pdgdau)==313){
703 for(Int_t resDau=dau->GetFirstDaughter(); resDau<=dau->GetLastDaughter(); resDau++){
704 if(resDau<0) return -1;
705 TParticle* resdau=(TParticle*)stack->Particle(resDau);
706 Int_t pdgresdau=resdau->GetPdgCode();
707 if(TMath::Abs(pdgresdau)==321){
708 if(pdgdp>0 && pdgresdau>0) return -1;
709 if(pdgdp<0 && pdgresdau<0) return -1;
710 nKaons++;
711 }
712 if(TMath::Abs(pdgresdau)==211){
713 nPions++;
714 }
715 }
716 }else if(TMath::Abs(pdgdau)==211){
717 nPions++;
718 }else{
719 return -1;
720 }
721 }
722 if(nPions!=3) return -1;
723 if(nKaons!=1) return -1;
724 for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
725 if(iDau<0) return -1;
726 }
727 return 1;
728 }
729
730 return -1;
731}
732
733
734//______________________________________________________________________________
735Int_t AliAnalysisTaskCheckHFMCProd::CheckDplusDecay(Int_t labDplus, AliStack* stack) const{
736
737 if(labDplus<0) return -1;
738 TParticle* dp = (TParticle*)stack->Particle(labDplus);
739 Int_t pdgdp=dp->GetPdgCode();
740 Int_t nDau=dp->GetNDaughters();
7b6a4dcd 741
742 if(nDau==3){
743 Int_t nKaons=0;
744 Int_t nPions=0;
745 for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
746 if(iDau<0) return -1;
747 TParticle* dau=(TParticle*)stack->Particle(iDau);
748 Int_t pdgdau=dau->GetPdgCode();
749 if(TMath::Abs(pdgdau)==321){
750 if(pdgdp>0 && pdgdau>0) return -1;
751 if(pdgdp<0 && pdgdau<0) return -1;
752 nKaons++;
753 }else if(TMath::Abs(pdgdau)==211){
754 if(pdgdp<0 && pdgdau>0) return -1;
755 if(pdgdp>0 && pdgdau<0) return -1;
756 nPions++;
757 }else{
758 return -1;
759 }
760 }
761 if(nPions!=2) return -1;
762 if(nKaons!=1) return -1;
763 for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
764 if(iDau<0) return -1;
765 }
766 return 0;
767 }
768
769 if(nDau==2){
770 Int_t nKaons=0;
771 Int_t nPions=0;
772 for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
773 if(iDau<0) return -1;
774 TParticle* dau=(TParticle*)stack->Particle(iDau);
775 Int_t pdgdau=dau->GetPdgCode();
776 if(TMath::Abs(pdgdau)==313){
777 for(Int_t resDau=dau->GetFirstDaughter(); resDau<=dau->GetLastDaughter(); resDau++){
778 if(resDau<0) return -1;
779 TParticle* resdau=(TParticle*)stack->Particle(resDau);
780 Int_t pdgresdau=resdau->GetPdgCode();
781 if(TMath::Abs(pdgresdau)==321){
782 if(pdgdp>0 && pdgresdau>0) return -1;
783 if(pdgdp<0 && pdgresdau<0) return -1;
784 nKaons++;
785 }
786 if(TMath::Abs(pdgresdau)==211){
787 if(pdgdp<0 && pdgresdau>0) return -1;
788 if(pdgdp>0 && pdgresdau<0) return -1;
789 nPions++;
790 }
791 }
792 }else if(TMath::Abs(pdgdau)==211){
793 if(pdgdp<0 && pdgdau>0) return -1;
794 if(pdgdp>0 && pdgdau<0) return -1;
795 nPions++;
796 }else{
797 return -1;
798 }
799 }
800 if(nPions!=2) return -1;
801 if(nKaons!=1) return -1;
802 for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
803 if(iDau<0) return -1;
804 }
805 return 1;
806 }
807 return -1;
808}
809
810//______________________________________________________________________________
811Int_t AliAnalysisTaskCheckHFMCProd::CheckDsDecay(Int_t labDs, AliStack* stack) const{
812 // Ds decay
813 if(labDs<0) return -1;
814 TParticle* dp = (TParticle*)stack->Particle(labDs);
815 Int_t pdgdp=dp->GetPdgCode();
816 Int_t nDau=dp->GetNDaughters();
7b6a4dcd 817
818 if(nDau==3){
819 Int_t nKaons=0;
820 Int_t nPions=0;
821 for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
822 if(iDau<0) return -1;
823 TParticle* dau=(TParticle*)stack->Particle(iDau);
824 Int_t pdgdau=dau->GetPdgCode();
825 if(TMath::Abs(pdgdau)==321){
826 nKaons++;
827 }else if(TMath::Abs(pdgdau)==211){
828 if(pdgdp<0 && pdgdau>0) return -1;
829 if(pdgdp>0 && pdgdau<0) return -1;
830 nPions++;
831 }else{
832 return -1;
833 }
834 }
835 if(nPions!=1) return -1;
836 if(nKaons!=2) return -1;
837 for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
838 if(iDau<0) return -1;
839 }
840 return 2;
841 }
842
843 if(nDau==2){
844 Int_t nKaons=0;
845 Int_t nPions=0;
846 Bool_t isPhi=kFALSE;
847 Bool_t isk0st=kFALSE;
848 for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
849 if(iDau<0) return -1;
850 TParticle* dau=(TParticle*)stack->Particle(iDau);
851 Int_t pdgdau=dau->GetPdgCode();
852 if(TMath::Abs(pdgdau)==313){
853 isk0st=kTRUE;
854 for(Int_t resDau=dau->GetFirstDaughter(); resDau<=dau->GetLastDaughter(); resDau++){
855 if(resDau<0) return -1;
856 TParticle* resdau=(TParticle*)stack->Particle(resDau);
857 Int_t pdgresdau=resdau->GetPdgCode();
858 if(TMath::Abs(pdgresdau)==321){
859 nKaons++;
860 }
861 if(TMath::Abs(pdgresdau)==211){
862 if(pdgdp<0 && pdgresdau>0) return -1;
863 if(pdgdp>0 && pdgresdau<0) return -1;
864 nPions++;
865 }
866 }
867 }else if(TMath::Abs(pdgdau)==333){
868 isPhi=kTRUE;
869 for(Int_t resDau=dau->GetFirstDaughter(); resDau<=dau->GetLastDaughter(); resDau++){
870 if(resDau<0) return -1;
871 TParticle* resdau=(TParticle*)stack->Particle(resDau);
48709ec3 872 if(!resdau) return -1;
7b6a4dcd 873 Int_t pdgresdau=resdau->GetPdgCode();
874 if(TMath::Abs(pdgresdau)==321){
875 nKaons++;
876 }else{
877 return -1;
878 }
879 }
880 }else if(TMath::Abs(pdgdau)==211){
881 if(pdgdp<0 && pdgdau>0) return -1;
882 if(pdgdp>0 && pdgdau<0) return -1;
883 nPions++;
884 }else if(TMath::Abs(pdgdau)==321){
885 nKaons++;
886 }else{
887 return -1;
888 }
889 }
890 if(nPions!=1) return -1;
891 if(nKaons!=2) return -1;
892 for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
893 if(iDau<0) return -1;
894 }
895 if(isk0st) return 1;
896 else if(isPhi) return 0;
897 else return 3;
898 }
899 return -1;
900}
901
902//______________________________________________________________________________
903Int_t AliAnalysisTaskCheckHFMCProd::CheckDstarDecay(Int_t labDstar, AliStack* stack) const{
904
905 if(labDstar<0) return -1;
906 TParticle* dp = (TParticle*)stack->Particle(labDstar);
907 Int_t pdgdp=dp->GetPdgCode();
908 Int_t nDau=dp->GetNDaughters();
909 if(nDau!=2) return -1;
910
911 Int_t nKaons=0;
912 Int_t nPions=0;
913 for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
914 if(iDau<0) return -1;
915 TParticle* dau=(TParticle*)stack->Particle(iDau);
916 Int_t pdgdau=dau->GetPdgCode();
917 if(TMath::Abs(pdgdau)==421){
918 for(Int_t resDau=dau->GetFirstDaughter(); resDau<=dau->GetLastDaughter(); resDau++){
919 if(resDau<0) return -1;
920 TParticle* resdau=(TParticle*)stack->Particle(resDau);
921 Int_t pdgresdau=resdau->GetPdgCode();
922 if(TMath::Abs(pdgresdau)==321){
923 if(pdgdp>0 && pdgresdau>0) return -1;
924 if(pdgdp<0 && pdgresdau<0) return -1;
925 nKaons++;
926 }
927 if(TMath::Abs(pdgresdau)==211){
928 if(pdgdp<0 && pdgresdau>0) return -1;
929 if(pdgdp>0 && pdgresdau<0) return -1;
930 nPions++;
931 }
932 }
933 }else if(TMath::Abs(pdgdau)==211){
934 if(pdgdp<0 && pdgdau>0) return -1;
935 if(pdgdp>0 && pdgdau<0) return -1;
936 nPions++;
937 }else{
938 return -1;
939 }
940 }
941 if(nPions!=2) return -1;
942 if(nKaons!=1) return -1;
943 for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
944 if(iDau<0) return -1;
945 }
946 return 0;
947
948}
949
950//______________________________________________________________________________
951Int_t AliAnalysisTaskCheckHFMCProd::CheckLcDecay(Int_t labLc, AliStack* stack) const{
952 if(labLc<0) return -1;
953 TParticle* dp = (TParticle*)stack->Particle(labLc);
954 Int_t pdgdp=dp->GetPdgCode();
955 Int_t nDau=dp->GetNDaughters();
7b6a4dcd 956
957 if(nDau==3){
958 Int_t nKaons=0;
959 Int_t nPions=0;
960 Int_t nProtons=0;
961 for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
962 if(iDau<0) return -1;
963 TParticle* dau=(TParticle*)stack->Particle(iDau);
964 Int_t pdgdau=dau->GetPdgCode();
965 if(TMath::Abs(pdgdau)==321){
966 if(pdgdp>0 && pdgdau>0) return -1;
967 if(pdgdp<0 && pdgdau<0) return -1;
968 nKaons++;
969 }else if(TMath::Abs(pdgdau)==211){
970 if(pdgdp<0 && pdgdau>0) return -1;
971 if(pdgdp>0 && pdgdau<0) return -1;
972 nPions++;
973 }else if(TMath::Abs(pdgdau)==2212){
974 if(pdgdp<0 && pdgdau>0) return -1;
975 if(pdgdp>0 && pdgdau<0) return -1;
976 nProtons++;
977 }else{
978 return -1;
979 }
980 }
981 if(nPions!=1) return -1;
982 if(nKaons!=1) return -1;
983 if(nProtons!=1) return -1;
984 for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
985 if(iDau<0) return -1;
986 }
987 return 0;
988 }
989
990 if(nDau==2){
991 Int_t nKaons=0;
992 Int_t nPions=0;
993 Int_t nProtons=0;
994 for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
995 if(iDau<0) return -1;
996 TParticle* dau=(TParticle*)stack->Particle(iDau);
997 Int_t pdgdau=dau->GetPdgCode();
998 if(TMath::Abs(pdgdau)==313 || TMath::Abs(pdgdau)==3124 || TMath::Abs(pdgdau)==2224
999 || TMath::Abs(pdgdau)==3122 || TMath::Abs(pdgdau)==311){
1000 Int_t nDauRes=dau->GetNDaughters();
1001 if(nDauRes!=2) return -1;
1002 for(Int_t resDau=dau->GetFirstDaughter(); resDau<=dau->GetLastDaughter(); resDau++){
1003 if(resDau<0) return -1;
1004 TParticle* resdau=(TParticle*)stack->Particle(resDau);
1005 Int_t pdgresdau=resdau->GetPdgCode();
1006 if(TMath::Abs(pdgresdau)==321){
1007 if(pdgdp>0 && pdgresdau>0) return -1;
1008 if(pdgdp<0 && pdgresdau<0) return -1;
1009 nKaons++;
1010 }
1011 if(TMath::Abs(pdgresdau)==211){
1012 if(pdgdp<0 && pdgresdau>0) return -1;
1013 if(pdgdp>0 && pdgresdau<0) return -1;
1014 nPions++;
1015 }
1016 if(TMath::Abs(pdgresdau)==2212){
1017 if(pdgdp<0 && pdgresdau>0) return -1;
1018 if(pdgdp>0 && pdgresdau<0) return -1;
1019 nProtons++;
1020 }
1021 }
1022 }else if(TMath::Abs(pdgdau)==321){
1023 if(pdgdp>0 && pdgdau>0) return -1;
1024 if(pdgdp<0 && pdgdau<0) return -1;
1025 nKaons++;
1026 }else if(TMath::Abs(pdgdau)==211){
1027 if(pdgdp<0 && pdgdau>0) return -1;
1028 if(pdgdp>0 && pdgdau<0) return -1;
1029 nPions++;
1030 }else if(TMath::Abs(pdgdau)==2212){
1031 if(pdgdp<0 && pdgdau>0) return -1;
1032 if(pdgdp>0 && pdgdau<0) return -1;
1033 nProtons++;
1034 }else{
1035 return -1;
1036 }
1037 }
1038 if(nPions!=1) return -1;
1039 if(nKaons!=1) return -1;
1040 if(nProtons!=1) return -1;
1041 for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
1042 if(iDau<0) return -1;
1043 }
1044 return 1;
1045 }
1046 return -1;
1047}
1048