1 #include "AliAnalysisTaskSE.h"
2 #include "AliAnalysisManager.h"
3 #include "AliAnalysisDataContainer.h"
4 #include "AliESDEvent.h"
5 #include "AliESDtrackCuts.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 "AliMultiplicity.h"
15 #include <TParticle.h>
23 #include "AliESDInputHandlerRP.h"
24 #include "AliAnalysisTaskCheckHFMCProd.h"
26 /**************************************************************************
27 * Copyright(c) 1998-2012, ALICE Experiment at CERN, All rights reserved. *
29 * Author: The ALICE Off-line Project. *
30 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
43 //*************************************************************************
44 // Implementation of class AliAnalysisTaskCheckHFMCProd
45 // AliAnalysisTask to check MC production at ESD+Kine level
48 // Authors: F. Prino, prino@to.infn.it
50 //*************************************************************************
52 ClassImp(AliAnalysisTaskCheckHFMCProd)
53 //______________________________________________________________________________
54 AliAnalysisTaskCheckHFMCProd::AliAnalysisTaskCheckHFMCProd() : AliAnalysisTaskSE("HFMCChecks"),
61 fHistoTrackletsEta1(0),
76 fHistOriginFeeddown(0),
81 fHistEtaPhiPtGenEle(0),
82 fHistEtaPhiPtGenPi(0),
84 fHistEtaPhiPtGenPro(0),
85 fHistEtaPhiPtRecEle(0),
86 fHistEtaPhiPtRecPi(0),
88 fHistEtaPhiPtRecPro(0),
89 fSearchUpToQuark(kFALSE),
95 DefineInput(0, TChain::Class());
96 DefineOutput(1, TList::Class());
100 //___________________________________________________________________________
101 AliAnalysisTaskCheckHFMCProd::~AliAnalysisTaskCheckHFMCProd(){
103 if (AliAnalysisManager::GetAnalysisManager()->IsProofMode()) return;
108 delete fESDtrackCuts;
112 //___________________________________________________________________________
113 void AliAnalysisTaskCheckHFMCProd::UserCreateOutputObjects() {
114 // create output histos
116 fOutput = new TList();
118 fOutput->SetName("OutputHistos");
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);
125 Double_t maxMult=100.;
126 if(fSystem==1) maxMult=10000.;
127 if(fSystem==2) maxMult=500.;
128 fHistoPhysPrim = new TH1F("hPhysPrim","",100,-0.5,maxMult-0.5);
129 fHistoPhysPrim->Sumw2();
130 fOutput->Add(fHistoPhysPrim);
131 fHistoTracks = new TH1F("hTracks","",100,-0.5,maxMult*2-0.5);
132 fHistoTracks->Sumw2();
133 fOutput->Add(fHistoTracks);
134 fHistoSelTracks = new TH1F("hSelTracks","",100,-0.5,maxMult-0.5);
135 fHistoSelTracks->Sumw2();
136 fOutput->Add(fHistoSelTracks);
137 fHistoTracklets = new TH1F("hTracklets","",100,-0.5,maxMult-0.5);
138 fHistoTracklets->Sumw2();
139 fOutput->Add(fHistoTracklets);
140 fHistoTrackletsEta1 = new TH1F("hTrackletsEta1","",100,-0.5,maxMult-0.5);
141 fHistoTrackletsEta1->Sumw2();
142 fOutput->Add(fHistoTrackletsEta1);
143 fHistoPtPhysPrim = new TH1F("hPtPhysPrim","",100,0.,20.);
144 fHistoPtPhysPrim->Sumw2();
145 fOutput->Add(fHistoPtPhysPrim);
146 fHistoEtaPhysPrim = new TH1F("hEtaPhysPrim","",100,-10.,10.);
147 fHistoEtaPhysPrim->Sumw2();
148 fOutput->Add(fHistoEtaPhysPrim);
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);
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);
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);
182 if(fSystem==1) nBinscb=200;
183 if(fSystem==2) nBinscb=21;
184 Double_t maxncn=nBinscb-0.5;
185 fHistoNcharmed = new TH2F("hncharmed","",100,-0.5,maxMult-0.5,nBinscb,-0.5,maxncn);
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);
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.);
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.);
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.);
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.);
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.);
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.);
229 for(Int_t ih=0; ih<5; ih++){
230 fHistBYPtAllDecay[ih]->Sumw2();
231 fHistBYPtAllDecay[ih]->SetMinimum(0);
232 fOutput->Add(fHistBYPtAllDecay[ih]);
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]);
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]);
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.);
257 for(Int_t ih=0; ih<2; ih++){
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]);
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);
278 fHistMotherID=new TH1F("hMotherID","",1000,-1.5,998.5);
279 fHistMotherID->SetMinimum(0);
280 fOutput->Add(fHistMotherID);
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);
308 fHistNcollHFtype=new TH2F("hNcollHFtype","",5,-1.5,3.5,30,-0.5,29.5);
309 fOutput->Add(fHistNcollHFtype);
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,
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);
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);
346 //______________________________________________________________________________
347 void AliAnalysisTaskCheckHFMCProd::UserExec(Option_t *)
351 AliESDEvent *esd = (AliESDEvent*) (InputEvent());
355 printf("AliAnalysisTaskSDDRP::Exec(): bad ESD\n");
359 fHistoNEvents->Fill(0);
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);
373 Int_t nTracks=esd->GetNumberOfTracks();
374 fHistoTracks->Fill(nTracks);
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;
383 fHistoSelTracks->Fill(nSelTracks);
385 const AliMultiplicity* mult=esd->GetMultiplicity();
386 Int_t nTracklets=mult->GetNumberOfTracklets();
388 for(Int_t it=0; it<nTracklets; it++){
389 Double_t eta=TMath::Abs(mult->GetEta(it));
390 if(eta<1) nTracklets1++;
392 fHistoTracklets->Fill(nTracklets);
393 fHistoTrackletsEta1->Fill(nTracklets1);
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());
401 if(spdv && spdv->IsFromVertexerZ()){
402 fHistoSPDZVtxX->Fill(spdv->GetXv());
403 fHistoSPDZVtxY->Fill(spdv->GetYv());
404 fHistoSPDZVtxZ->Fill(spdv->GetZv());
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());
415 AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
417 Printf("ERROR: Could not retrieve MC event handler");
420 AliMCEvent* mcEvent = eventHandler->MCEvent();
422 Printf("ERROR: Could not retrieve MC event");
425 stack = mcEvent->Stack();
427 Printf("ERROR: stack not available");
430 const AliVVertex* mcVert=mcEvent->GetPrimaryVertex();
432 Printf("ERROR: generated vertex not available");
435 if(TMath::Abs(mcVert->GetZ())>10) return;
437 // const AliHeader* h=(AliHeader*)mcEvent->GetHeader();
439 TString genname=mcEvent->GenEventHeader()->ClassName();
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;
454 nColl=lgen->GetEntries();
455 fHistNcollHFtype->Fill(typeHF,nColl);
457 Int_t nParticles=stack->GetNtrack();
458 Double_t dNchdy = 0.;
462 for (Int_t i=0;i<nParticles;i++){
463 TParticle* part = (TParticle*)stack->Particle(i);
464 Int_t absPdg=TMath::Abs(part->GetPdgCode());
465 Int_t pdg=part->GetPdgCode();
468 if(stack->IsPhysicalPrimary(i)){
469 Double_t eta=part->Eta();
470 fHistoEtaPhysPrim->Fill(eta);
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());
476 if(TMath::Abs(eta)<0.5){
477 dNchdy+=0.6666; // 2/3 for the ratio charged/all
480 if(TMath::Abs(eta)<0.9){
481 fHistoPtPhysPrim->Fill(part->Pt());
485 if (part->Energy() != TMath::Abs(part->Pz())){
486 rapid=0.5*TMath::Log((part->Energy()+part->Pz())/(part->Energy()-part->Pz()));
495 iType=CheckD0Decay(i,stack);
496 if(iType>=0) iPart=0;
498 else if(absPdg==411){
500 iType=CheckDplusDecay(i,stack);
501 if(iType>=0) iPart=1;
503 else if(absPdg==413){
505 iType=CheckDstarDecay(i,stack);
506 if(iType>=0) iPart=2;
508 else if(absPdg==431){
510 iType=CheckDsDecay(i,stack);
511 if(iType==0 || iType==1) iPart=3;
513 else if(absPdg==4122){
515 iType=CheckLcDecay(i,stack);
516 if(iType>=0) iPart=4;
518 if(iSpecies>=0) fHistYPtAllDecay[iSpecies]->Fill(part->Pt(),rapid);
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);
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);
538 if(iSpecies<0) continue; // not a charm meson
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);
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);
555 fHistMotherID->Fill(part->GetFirstMother());
556 TParticle* runningpart=part;
559 if(fSearchUpToQuark){
561 Int_t labmoth=runningpart->GetFirstMother();
562 if(labmoth==-1) break;
563 TParticle *mot=(TParticle*)stack->Particle(labmoth);
564 pdgmoth=TMath::Abs(mot->GetPdgCode());
568 }else if(pdgmoth==4){
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){
585 if(pdgmoth>=5000 && pdgmoth<=5999){
594 fHistYPtPromptAllDecay[iSpecies]->Fill(part->Pt(),rapid);
595 fHistOriginPrompt->Fill(distToVert);
598 fHistYPtFeeddownAllDecay[iSpecies]->Fill(part->Pt(),rapid);
599 fHistOriginFeeddown->Fill(distToVert);
602 if(iPart<0) continue;
603 if(iType<0) continue;
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);
613 if(iFromB==0 && iPart>=0 && iPart<5) fHistYPtPrompt[iPart]->Fill(part->Pt(),rapid);
614 else if(iFromB==1 && iPart>=0 && iPart<5) fHistYPtFeeddown[iPart]->Fill(part->Pt(),rapid);
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());
632 fHistoNcharmed->Fill(dNchdy,nCharmed);
633 fHistoNbVsNc->Fill(nc,nb);
634 fHistoPhysPrim->Fill(nPhysPrim);
640 //______________________________________________________________________________
641 void AliAnalysisTaskCheckHFMCProd::Terminate(Option_t */*option*/)
643 // Terminate analysis
644 fOutput = dynamic_cast<TList*> (GetOutputData(1));
646 printf("ERROR: fOutput not available\n");
656 //______________________________________________________________________________
657 Int_t AliAnalysisTaskCheckHFMCProd::CheckD0Decay(Int_t labD0, AliStack* stack) const{
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();
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;
675 }else if(TMath::Abs(pdgdau)==211){
676 if(pdgdp<0 && pdgdau>0) return -1;
677 if(pdgdp>0 && pdgdau<0) return -1;
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;
691 if(nDau==3 || nDau==4){
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;
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;
712 if(TMath::Abs(pdgresdau)==211){
716 }else if(TMath::Abs(pdgdau)==211){
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;
734 //______________________________________________________________________________
735 Int_t AliAnalysisTaskCheckHFMCProd::CheckDplusDecay(Int_t labDplus, AliStack* stack) const{
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();
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;
753 }else if(TMath::Abs(pdgdau)==211){
754 if(pdgdp<0 && pdgdau>0) return -1;
755 if(pdgdp>0 && pdgdau<0) return -1;
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;
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;
786 if(TMath::Abs(pdgresdau)==211){
787 if(pdgdp<0 && pdgresdau>0) return -1;
788 if(pdgdp>0 && pdgresdau<0) return -1;
792 }else if(TMath::Abs(pdgdau)==211){
793 if(pdgdp<0 && pdgdau>0) return -1;
794 if(pdgdp>0 && pdgdau<0) return -1;
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;
810 //______________________________________________________________________________
811 Int_t AliAnalysisTaskCheckHFMCProd::CheckDsDecay(Int_t labDs, AliStack* stack) const{
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();
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){
827 }else if(TMath::Abs(pdgdau)==211){
828 if(pdgdp<0 && pdgdau>0) return -1;
829 if(pdgdp>0 && pdgdau<0) return -1;
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;
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){
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){
861 if(TMath::Abs(pdgresdau)==211){
862 if(pdgdp<0 && pdgresdau>0) return -1;
863 if(pdgdp>0 && pdgresdau<0) return -1;
867 }else if(TMath::Abs(pdgdau)==333){
869 for(Int_t resDau=dau->GetFirstDaughter(); resDau<=dau->GetLastDaughter(); resDau++){
870 if(resDau<0) return -1;
871 TParticle* resdau=(TParticle*)stack->Particle(resDau);
872 if(!resdau) return -1;
873 Int_t pdgresdau=resdau->GetPdgCode();
874 if(TMath::Abs(pdgresdau)==321){
880 }else if(TMath::Abs(pdgdau)==211){
881 if(pdgdp<0 && pdgdau>0) return -1;
882 if(pdgdp>0 && pdgdau<0) return -1;
884 }else if(TMath::Abs(pdgdau)==321){
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;
896 else if(isPhi) return 0;
902 //______________________________________________________________________________
903 Int_t AliAnalysisTaskCheckHFMCProd::CheckDstarDecay(Int_t labDstar, AliStack* stack) const{
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;
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;
927 if(TMath::Abs(pdgresdau)==211){
928 if(pdgdp<0 && pdgresdau>0) return -1;
929 if(pdgdp>0 && pdgresdau<0) return -1;
933 }else if(TMath::Abs(pdgdau)==211){
934 if(pdgdp<0 && pdgdau>0) return -1;
935 if(pdgdp>0 && pdgdau<0) return -1;
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;
950 //______________________________________________________________________________
951 Int_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();
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;
969 }else if(TMath::Abs(pdgdau)==211){
970 if(pdgdp<0 && pdgdau>0) return -1;
971 if(pdgdp>0 && pdgdau<0) return -1;
973 }else if(TMath::Abs(pdgdau)==2212){
974 if(pdgdp<0 && pdgdau>0) return -1;
975 if(pdgdp>0 && pdgdau<0) return -1;
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;
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;
1011 if(TMath::Abs(pdgresdau)==211){
1012 if(pdgdp<0 && pdgresdau>0) return -1;
1013 if(pdgdp>0 && pdgresdau<0) return -1;
1016 if(TMath::Abs(pdgresdau)==2212){
1017 if(pdgdp<0 && pdgresdau>0) return -1;
1018 if(pdgdp>0 && pdgresdau<0) return -1;
1022 }else if(TMath::Abs(pdgdau)==321){
1023 if(pdgdp>0 && pdgdau>0) return -1;
1024 if(pdgdp<0 && pdgdau<0) return -1;
1026 }else if(TMath::Abs(pdgdau)==211){
1027 if(pdgdp<0 && pdgdau>0) return -1;
1028 if(pdgdp>0 && pdgdau<0) return -1;
1030 }else if(TMath::Abs(pdgdau)==2212){
1031 if(pdgdp<0 && pdgdau>0) return -1;
1032 if(pdgdp>0 && pdgdau<0) return -1;
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;