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