]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG3/dielectron/AliAnalysisTaskDielectronReadAODBranch.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
fix for coverity defects (Jens)
[u/mrichter/AliRoot.git] / PWG3 / dielectron / AliAnalysisTaskDielectronReadAODBranch.cxx
CommitLineData
ba15fdfb 1/**************************************************************************
2 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16////////////////////////////////////////////////////////////////////////////
17// Class to retrieve the branch of the dielectron candidates stored in //
18// filtered AODs file (AliAOD.Dielectron.root). It is possible to //
19// apply tighter cuts to candidates stored in the branch and do the //
20// matching with MC truth. //
21////////////////////////////////////////////////////////////////////////////
22
23#include "TChain.h"
24#include "TNtuple.h"
25#include "TH1F.h"
26#include "TH2F.h"
27#include "TDatabasePDG.h"
28#include "TF1.h"
29#include "AliAnalysisManager.h"
30#include "AliAODHandler.h"
31#include "AliAODEvent.h"
32#include "AliAODVertex.h"
33#include "AliAODTrack.h"
34#include "AliAODMCParticle.h"
35#include "AliAnalysisTaskDielectronReadAODBranch.h"
36#include "AliDielectronPair.h"
37
38ClassImp(AliAnalysisTaskDielectronReadAODBranch)
39 //________________________________________________________________________
40 AliAnalysisTaskDielectronReadAODBranch::AliAnalysisTaskDielectronReadAODBranch():
41 AliAnalysisTaskSE(),
42 fOutput(0),
43 fNtupleJPSI(0),
44 fNentries(0),
45 fInvMass(0),
46 fInvMassNoCuts(0),
47 fpsproperSignal(0),
48 fpsproperSidebands(0),
49 fpsproperAll(0),
50 fpsproperUnder(0),
51 fpsproperUpper(0),
52 fLxyVsPtleg1(0),
53 fLxyVsPtleg2(0),
54 fLxyVsPt(0),
55 fMeeVsPt(0),
56 fMeeVsLxy(0),
57 fprimvtxZ(0),
58 fsecvtxZ(0),
59 fprimvtxX(0),
60 fsecvtxX(0),
61 fprimvtxY(0),
62 fsecvtxY(0),
63 fPt(0),
64 fPtLeg1(0),
65 fPtLeg2(0),
66 fdEdxP(0),
67 fHasMC(0),
68 fobj(0),
69 fobjMC(0),
70 fPtCut(1.),
71 fSpdFirstRequired(kFALSE),
72 fClsTPC(90),
73 fPairType(1),
74 fPtJpsi(1.3),
75 fInvMassSignalLimits(0),
76 fInvMassSideBandsLimits(0)
77{
78 // Default constructor
79}
80
81//________________________________________________________________________
82AliAnalysisTaskDielectronReadAODBranch::AliAnalysisTaskDielectronReadAODBranch(const char *name):
83 AliAnalysisTaskSE(name),
84 fOutput(0),
85 fNtupleJPSI(0),
86 fNentries(0),
87 fInvMass(0),
88 fInvMassNoCuts(0),
89 fpsproperSignal(0),
90 fpsproperSidebands(0),
91 fpsproperAll(0),
92 fpsproperUnder(0),
93 fpsproperUpper(0),
94 fLxyVsPtleg1(0),
95 fLxyVsPtleg2(0),
96 fLxyVsPt(0),
97 fMeeVsPt(0),
98 fMeeVsLxy(0),
99 fprimvtxZ(0),
100 fsecvtxZ(0),
101 fprimvtxX(0),
102 fsecvtxX(0),
103 fprimvtxY(0),
104 fsecvtxY(0),
105 fPt(0),
106 fPtLeg1(0),
107 fPtLeg2(0),
108 fdEdxP(0),
109 fHasMC(0),
110 fobj(0),
111 fobjMC(0),
112 fPtCut(1.),
113 fSpdFirstRequired(kFALSE),
114 fClsTPC(90),
115 fPairType(1),
116 fPtJpsi(1.3),
117 fInvMassSignalLimits(0),
118 fInvMassSideBandsLimits(0)
119{
120 // Default constructor
121 DefineInput(0,TChain::Class());
122 DefineOutput(1,TList::Class()); //My private output
123
124 fInvMassSignalLimits = new Double_t[2];
125 fInvMassSignalLimits[0] = 0.; fInvMassSignalLimits[1] = 0.;
126 fInvMassSideBandsLimits = new Double_t[2];
127 fInvMassSideBandsLimits[0] = 0.; fInvMassSideBandsLimits[1] = 0.;
128}
129
130//________________________________________________________________________
131AliAnalysisTaskDielectronReadAODBranch::~AliAnalysisTaskDielectronReadAODBranch()
132{
133 if(fOutput) delete fOutput;
134 if(fobj) delete fobj;
135 if(fobjMC) delete fobjMC;
136}
137
138//________________________________________________________________________
139void AliAnalysisTaskDielectronReadAODBranch::Init()
140{
141 // Initialization
142 if(fDebug > 1) printf("AnalysisTaskReadAOD::Init() \n");
143 return;
144}
145
146//________________________________________________________________________
147void AliAnalysisTaskDielectronReadAODBranch::UserCreateOutputObjects()
148{
149 //
150 // Create the output container
151 //
152 if(fDebug > 1) printf("AnalysisTaskReadAOD::UserCreateOutputObjects() \n");
153 fOutput = new TList();
154 fOutput->SetOwner();
155 //// Histogram booking
156
157 // invariant mass
158 fInvMass = new TH1F("fInvMass","fInvMass",300,2.0,2.0+300*.04); // step 40MeV
159 fInvMassNoCuts = new TH1F("fInvMass_no_cuts","fInvMass_no_cuts",125,0,125*0.04); //step 40MeV
160 fNentries = new TH1F("numberOfevent","numbOfEvent",1,-0.5,0.5);
161
162 //pseudoproper
163 fpsproperSignal = new TH1F("psproper_decay_length",Form("psproper_decay_length_distrib(AODcuts+#clustTPC>%d, pT(leg)>%fGeV, pT(J/#psi)>%fGeV/c, %f < M < %f);X [#mu m];Entries/40#mu m",fClsTPC,fPtCut,fPtJpsi,fInvMassSignalLimits[0],fInvMassSignalLimits[1]),150,-3000.,3000.);
164 fpsproperSidebands = new TH1F("psproper_decay_length_sidebands",Form("psproper_decay_length_distrib_sidebands(AODcuts+#clustTPC>%d, pT(e+e-)>%f GeV,pT(J/#psi)>%f GeV/c,M> %f && M < %f ); X [#mu m]; Entries/40#mu m",fClsTPC,fPtCut,fPtJpsi,fInvMassSideBandsLimits[1],fInvMassSideBandsLimits[0]),150,-3000.,3000.);
165 fpsproperAll = new TH1F("psproper_decay_length_all",Form("psproper_decay_length_distrib_all(AODcuts+#clustTPC>%d, pT(e+e-)>%f GeV,pT(J/#psi)>%f GeV/c);X [#mu m];Entries/15#mu m",fClsTPC,fPtCut,fPtJpsi),400,-3000.,3000.);
166 fpsproperUnder = new TH1F("psproper_decay_length_under",Form("psproper_decay_length_distrib_under(AODcuts+#clustTPC>%d,pT(e+e-)>%f GeV,pT(J/#psi)>%f GeV/c);X [#mu m];Entries/15#mu m",fClsTPC,fPtCut,fPtJpsi),400,-3000.,3000.);
167 fpsproperUpper = new TH1F("psproper_decay_length_upper",Form("psproper_decay_length_distrib_upper(AODcuts+#clustTPC>%d,pT(e+e-)>%f GeV,pT(J/#psi)>%f GeV/c);X [#mu m];Entries/15#mu m",fClsTPC,fPtCut,fPtJpsi),400,-3000.,3000.);
168 //
169 fLxyVsPtleg1 = new TH2F("Lxy_vs_Pt_leg1",Form("Lxy_vs_Pt_leg1_distrib(AODcuts+#clustTPC>%d,pT(e+e-)>%f GeV,pT(J/#psi)>%f GeV/c,%f<M<%f);p_{T}[GeV/c]; X",fClsTPC,fPtCut,fPtJpsi,fInvMassSignalLimits[0],fInvMassSignalLimits[1]), 700, 0., 7.,400, -3000.,3000.);
170 fLxyVsPtleg2 = new TH2F("Lxy_vs_Pt_leg2",Form("Lxy_vs_Pt_leg2_distrib(AODcuts+#clustTPC>%d,pT(e+e-)>%f GeV,pT(J/#psi)>%f GeV/c,%f<M<%f);p_{T}[GeV/c]; X",fClsTPC,fPtCut,fPtJpsi,fInvMassSignalLimits[0],fInvMassSignalLimits[1]), 700, 0., 7.,400, -3000.,3000.);
171 fLxyVsPt = new TH2F("Lxy_vs_Pt_jpsi",Form("Lxy_vs_Pt_jpsi_distrib(AODcuts+#clustTPC>%d,pT(e+e-)>%f GeV,pT(J/#psi)>%f GeV/c,%f<M<%f); p_{T}[GeV/c]; X",fClsTPC,fPtCut,fPtJpsi,fInvMassSignalLimits[0],fInvMassSignalLimits[1]), 700, 0., 7.,400, -3000.,3000.);
172 fMeeVsPt = new TH2F("Mee_vs_Pt_jpsi",Form("Mee_vs_Pt_jpsi_distrib(AODcuts+#clustTPC>%d,pT(e+e-)>%f GeV,pT(J/#psi)>%f GeV/c,%f<M<%f); p_{T}[GeV/c]; M_{ee}",fClsTPC,fPtCut,fPtJpsi,fInvMassSignalLimits[0],fInvMassSignalLimits[1]), 700, 0., 7.,200, 1.99,4.1);
173 fMeeVsLxy = new TH2F("Mee_vs_Lxy_jpsi",Form("Mee_vs_Lxy_jpsi_distrib(AODcuts+#clustTPC>%d,pT(e+e-)>%f GeV),pT(J/#psi)>%f GeV/c; X; M_{ee}",fClsTPC,fPtCut,fPtJpsi), 400, -3000, 3000.,200, 1.99,4.1);
174
175 // QA plots
176 fprimvtxZ = new TH1F("prim_vtx_Z",Form("prim_vtx_Z_distrib(AODcuts + #clustTPC>%d, pT(e+e-)>%f GeV, pT(J/#psi)>%f); Z[cm]; Entries/20 mm",fClsTPC,fPtCut,fPtJpsi),1000,-10.,10.);
177 fsecvtxZ = new TH1F("sec_vtx_Z",Form("sec_vtx_Z_distrib(AODcuts + #clustTPC>%d, pT(e+e-)>%f GeV, pT(J/#psi)>%f); Z[cm]; Entries/20 mm",fClsTPC,fPtCut,fPtJpsi),1000,-10.,10.);
178 fprimvtxX = new TH1F("prim_vtx_X",Form("prim_vtx_X_distrib(AODcuts + #clustTPC>%d, pT(e+e-)>%f GeV,pT(J/#psi)>%f); X[cm]; Entries/mm",fClsTPC,fPtCut,fPtJpsi),1000,-0.5,0.5);
179 fsecvtxX = new TH1F("sec_vtx_X",Form("sec_vtx_X_distrib(AODcuts + #clustTPC>%d, pT(e+e-)>%f GeV,pT(J/#psi)>%f); X[cm]; Entries/20 mm",fClsTPC,fPtCut,fPtJpsi),100,-1.,1.);
180 fprimvtxY = new TH1F("prim_vtx_Y",Form("prim_vtx_Y_distrib(AODcuts + #clustTPC>%d, pT(e+e-)>%f GeV,pT(J/#psi)>%f); Y[cm]; Entries/mm",fClsTPC,fPtCut,fPtJpsi),1000,-0.5,0.5);
181 fsecvtxY = new TH1F("sec_vtx_Y",Form("sec_vtx_Y_distrib(AODcuts + #clustTPC>%d, pT(e+e-)>%f GeV,pT(J/#psi)>%f); Y[cm]; Entries/20 mm",fClsTPC,fPtCut,fPtJpsi),100,-1.,1.);
182 //
183 fPt = new TH1F("Pt(J/psi)",Form("Pt_Jpsi_distrib(AODcuts + #clustTPC>%d, pT(e+e-)>%f GeV, pT(J/#psi)>%f); p_{T}(J/#psi)[GeV/c]; Entries/25 MeV",fClsTPC,fPtCut,fPtJpsi),400,0.,10.);
184 fPtLeg1 = new TH1F("Pt_leg1",Form("Pt_leg1_distrib(AODcuts + #clustTPC>%d, pT(e+e-) > %f GeV, pT(J/#psi)>%f); p_{T}[GeV/c]; Entries/25 MeV",fClsTPC,fPtCut,fPtJpsi),400,0.,10.);
185 fPtLeg2 = new TH1F("Pt_leg2",Form("Pt_leg2_distrib(AODcuts + #clustTPC>%d, pT(e+e-) > %f GeV, pT(J/#psi)>%f); p_{T}[GeV/c]; Entries/25 MeV",fClsTPC,fPtCut,fPtJpsi),400,0.,10.);
186 //dE/dx plots
187 fdEdxP = new TH2F("dE/dx_vs_Ptpc",Form("dE/dx_vs_Ptpc(AODcuts + #clustTPC>%d, pT(e+e-) > %f GeV, pT(J/#psi)>%f)",fClsTPC,fPtCut,fPtJpsi),400,0.2,20.,200,0.,200.);
188
189 //J/psi n-tuple (X and M values)
190 fNtupleJPSI = new TNtuple("fNtupleJPSI","J/#psi pseudo-proper decay time & invariant mass","Xdecaytime:Mass");
191
192 // add histograms to list
193 fOutput->Add(fNentries);
194 fOutput->Add(fInvMass);
195 fOutput->Add(fInvMassNoCuts);
196 fOutput->Add(fpsproperSignal);
197 fOutput->Add(fpsproperSidebands);
198 fOutput->Add(fpsproperAll);
199 fOutput->Add(fpsproperUnder);
200 fOutput->Add(fpsproperUpper);
201 fOutput->Add(fLxyVsPtleg1);
202 fOutput->Add(fLxyVsPtleg2);
203 fOutput->Add(fLxyVsPt);
204 fOutput->Add(fMeeVsPt);
205 fOutput->Add(fMeeVsLxy);
206 fOutput->Add(fprimvtxZ);
207 fOutput->Add(fsecvtxZ);
208 fOutput->Add(fprimvtxX);
209 fOutput->Add(fsecvtxX);
210 fOutput->Add(fprimvtxY);
211 fOutput->Add(fsecvtxY);
212 fOutput->Add(fPt);
213 fOutput->Add(fPtLeg1);
214 fOutput->Add(fPtLeg2);
215 fOutput->Add(fNtupleJPSI);
216 fOutput->Add(fdEdxP);
217 return;
218}
219
220//________________________________________________________________________
221void AliAnalysisTaskDielectronReadAODBranch::UserExec(Option_t */*option*/)
222{
223 // Execute analysis for current event:
224 AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
c8f0f810 225 if (!aod) return;
226
ba15fdfb 227 Double_t vtxPrim[3] = {0.,0.,0.};
228
229 AliAODVertex* primvtx = aod->GetPrimaryVertex();
230 vtxPrim[0] = primvtx->GetX();
231 vtxPrim[1] = primvtx->GetY();
232 vtxPrim[2] = primvtx->GetZ();
233
234 AliAODHandler* aodHandler = (AliAODHandler*)((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
235
236 TTree *aodTree = aodHandler->GetTree();
237 if(!fobj) aodTree->SetBranchAddress("dielectrons",&fobj);
238
239 if(fHasMC){ aodTree->SetBranchAddress("mcparticles",&fobjMC);
240 if(!fobjMC) printf("AliAnalysisTaskDielectronReadAODBranch::UserExec: MC particles branch not found!\n"); }
241
242 fNentries->Fill(0);
243 aodTree->GetEvent(Entry());
244 Int_t dgLabels[2] = {0,0};
245 Int_t pdgDgJPSItoEE[2] = {11,11};
246
247 // loop over candidates
248 if(fobj) {
249 TObjArray *objArr = (TObjArray*)fobj->UncheckedAt(fPairType);
250 for(int j=0;j<objArr->GetEntriesFast();j++)
251 {
252 AliDielectronPair *pairObj = (AliDielectronPair*)objArr->UncheckedAt(j);
253 fInvMassNoCuts->Fill(pairObj->M());
254 Double_t vtxSec[3] = {0.,0.,0.};
255 fprimvtxX->Fill(vtxPrim[0]);
256 fprimvtxY->Fill(vtxPrim[1]);
257 fprimvtxZ->Fill(vtxPrim[2]);
258 pairObj->XvYvZv(vtxSec);
259 fsecvtxX->Fill(vtxSec[0]);
260 fsecvtxY->Fill(vtxSec[1]);
261 fsecvtxZ->Fill(vtxSec[2]);
262
263 Double_t lxy = ((vtxSec[0]-vtxPrim[0])*(pairObj->Px()) + (vtxSec[1]-vtxPrim[1])*(pairObj->Py()))/pairObj->Pt();
264 Double_t psProperDecayLength = lxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/pairObj->Pt();
265
266 AliAODTrack *trk = (AliAODTrack*)pairObj->GetFirstDaughter();
267 AliAODTrack *trk1 = (AliAODTrack*)pairObj->GetSecondDaughter();
268 if(!trk || !trk1) {printf("ERROR: daughter tracks not available\n"); continue;}
269 dgLabels[0] = trk->GetLabel();
270 dgLabels[1] = trk1->GetLabel();
271
272 //check in case of MC analysis if candidate is a true J/psi->ee
273 if(fHasMC && !MatchToMC(443,fobjMC,dgLabels,2,2,pdgDgJPSItoEE)) continue;
274 AliAODPid *pid = trk->GetDetPid();
275 AliAODPid *pid1 = trk1->GetDetPid();
276 if(!pid || !pid1){
277 if(fDebug>1) printf("No AliAODPid found\n");
278 continue;
279 }
280 // ptLeg cut
281 if((trk->Pt()<fPtCut) || (trk1->Pt()<fPtCut)) continue;
282
283 // pt jpsi cut
284 if((pairObj->Pt())<fPtJpsi) continue;
285
286 //spd first required
287 if(fSpdFirstRequired)
288 {
289 if((!trk->HasPointOnITSLayer(0)) || (!trk1->HasPointOnITSLayer(0))) continue;
290 }
291 //
292 if((trk->GetTPCNcls()) <= fClsTPC || (trk1->GetTPCNcls()) <= fClsTPC) continue;
293 if(trk->Charge()==-1){
294 fPtLeg1->Fill(trk->Pt());
295 } else {fPtLeg2->Fill(trk->Pt());}
296 if(trk1->Charge()==-1){
297 fPtLeg1->Fill(trk1->Pt());
298 } else {fPtLeg2->Fill(trk1->Pt());}
299 //Fill dE/dx related plots (before PID cuts)
300 fdEdxP->Fill(pid->GetTPCmomentum(),pid->GetTPCsignal());
301 fdEdxP->Fill(pid1->GetTPCmomentum(),pid1->GetTPCsignal());
302 fInvMass->Fill(pairObj->M());
303 fPt->Fill(pairObj->Pt());
304 fMeeVsLxy->Fill(10000*psProperDecayLength,pairObj->M());
305 fMeeVsPt->Fill(pairObj->Pt(),pairObj->M());
306 fpsproperAll->Fill(10000*psProperDecayLength);
307
308 //psproper in the signal region
309 if((pairObj->M())<fInvMassSignalLimits[1] && (pairObj->M())>fInvMassSignalLimits[0]){
310 fpsproperSignal->Fill(10000*psProperDecayLength);
311 fLxyVsPt->Fill(pairObj->Pt(),10000*psProperDecayLength); //jpsi
312 fLxyVsPtleg1->Fill(trk->Pt(),10000*psProperDecayLength); //leg1 (warning: mixture of pos and neg tracks)
313 fLxyVsPtleg2->Fill(trk1->Pt(),10000*psProperDecayLength); //leg2 (warning: mixture of pos and neg tracks)
314 fNtupleJPSI->Fill(10000*psProperDecayLength,pairObj->M()); // fill the N-tuple (imput of the minimization algorithm)
315 }
316
317 // psproper in the sidebands
318 if((pairObj->M())> fInvMassSideBandsLimits[1] || (pairObj->M())< fInvMassSideBandsLimits[0])
319 fpsproperSidebands->Fill(10000*psProperDecayLength);
320
321 // psproper in the lower sideband
322 if(pairObj->M()< fInvMassSideBandsLimits[0]) fpsproperUnder->Fill(10000*psProperDecayLength);
323
324 // psproper in the upper sideband
325 if(pairObj->M()> fInvMassSideBandsLimits[1]) fpsproperUpper->Fill(10000*psProperDecayLength);
326 } // end loop over candidates
327 }// end loop over pair types
328 // Post the data
329 PostData(1,fOutput);
330 return;
331}
332
333//________________________________________________________________________
334void AliAnalysisTaskDielectronReadAODBranch::Terminate(Option_t */*option*/)
335{
336 if(fDebug > 1) printf("AliAnalysisTaskDielectronReadAODBranch::Terminate() \n");
337 return;
338}
339
340//___________________________________________________________________
341Bool_t AliAnalysisTaskDielectronReadAODBranch::MatchToMC(Int_t pdgabs,const TObjArray *mcArray,
342 const Int_t *dgLabels,Int_t ndg,
343 Int_t ndgCk, const Int_t *pdgDg) const
344{
345 //
346 // jpsi candidates association to MC truth
347 // Check if this candidate is matched to a MC signal
348 // If no, return kFALSE
349 // If yes, return kTRUE
c8f0f810 350 //
351 if (!mcArray) return kFALSE;
352
ba15fdfb 353 Int_t labMom[2]={0,0};
354 Int_t i,j,lab,labMother,pdgMother,pdgPart,labJPSIMother,pdgJPSIMother;
355 AliAODMCParticle *part=0;
356 AliAODMCParticle *mother=0;
357 AliAODMCParticle *jPSImother=0;
358 Double_t pxSumDgs=0.,pySumDgs=0.,pzSumDgs=0.;
359 Bool_t pdgUsed[2]={kFALSE,kFALSE};
360
361 // loop on daughter labels
362 for(i=0; i<ndg; i++) {
363 labMom[i]=-1;
364 lab = dgLabels[i];
365 if(lab<0) {
366 printf("daughter with negative label %d\n",lab);
367 return kFALSE;
368 }
369 part = (AliAODMCParticle*)mcArray->At(lab);
370 if(!part) {
371 printf("no MC particle\n");
372 return kFALSE;
373 }
374
375 // check the PDG of the daughter, if requested
376 if(ndgCk>0) {
377 pdgPart=TMath::Abs(part->GetPdgCode());
378 printf("pdg code of daughter %d == %d\n",i,pdgPart);
379 for(j=0; j<ndg; j++) {
380 if(!pdgUsed[j] && pdgPart==pdgDg[j]) {
381 pdgUsed[j]=kTRUE;
382 break;
383 }
384 }
385 }
386
387 // for the J/psi, check that the daughters are electrons
388 if(pdgabs==443) {
389 if(TMath::Abs(part->GetPdgCode())!=11) return kFALSE;
390 }
391
392 mother = part;
393 while(mother->GetMother()>=0) {
394 labMother=mother->GetMother();
395 mother = (AliAODMCParticle*)mcArray->At(labMother);//get particle mother
396 if(!mother) {
397 printf("no MC mother particle\n");
398 break;
399 }
400 pdgMother = TMath::Abs(mother->GetPdgCode());//get particle mother pdg code
401 printf("pdg code of mother track == %d\n",pdgMother);
402 if(pdgMother==pdgabs) {
403 labJPSIMother=mother->GetMother();
404 jPSImother = (AliAODMCParticle*)mcArray->At(labJPSIMother);//get J/psi mother
405 if(jPSImother) {
406 pdgJPSIMother = TMath::Abs(jPSImother->GetPdgCode()); //get J/psi mother pdg code
407 if( pdgJPSIMother==511 || pdgJPSIMother==521 ||
408 pdgJPSIMother==10511 || pdgJPSIMother==10521 ||
409 pdgJPSIMother==513 || pdgJPSIMother==523 ||
410 pdgJPSIMother==10513 || pdgJPSIMother==10523 ||
411 pdgJPSIMother==20513 || pdgJPSIMother==20523 ||
412 pdgJPSIMother==515 || pdgJPSIMother==525 ||
413 pdgJPSIMother==531 || pdgJPSIMother==10531 ||
414 pdgJPSIMother==533 || pdgJPSIMother==10533 ||
415 pdgJPSIMother==20533 || pdgJPSIMother==535 ||
416 pdgJPSIMother==541 || pdgJPSIMother==10541 ||
417 pdgJPSIMother==543 || pdgJPSIMother==10543 ||
418 pdgJPSIMother==20543 || pdgJPSIMother==545) return kFALSE;} //check if J/psi comes from B hadron
419
420 labMom[i]=labMother;
421 // keep sum of daughters' momenta, to check for mom conservation
422 pxSumDgs += part->Px();
423 pySumDgs += part->Py();
424 pzSumDgs += part->Pz();
425 break;
426 } else if(pdgMother>pdgabs || pdgMother<10) {
427 break;
428 }
429 }
430 if(labMom[i]==-1) {printf("mother PDG not ok for this daughter\n"); return kFALSE;} // mother PDG not ok for this daughter
431 } // end loop on daughters
432
433 // check if the candidate is signal
434 labMother=labMom[0];
435 // all labels have to be the same and !=-1
436 for(i=0; i<ndg; i++) {
437 if(labMom[i]==-1) return kFALSE;
438 if(labMom[i]!=labMother) return kFALSE;
439 }
440
441 // check that all daughter PDGs are matched
442 if(ndgCk>0) {
443 for(i=0; i<ndg; i++) {
444 if(pdgUsed[i]==kFALSE) return kFALSE;
445 }
446 }
447
448 mother = (AliAODMCParticle*)mcArray->At(labMother);
449 Double_t pxMother = mother->Px();
450 Double_t pyMother = mother->Py();
451 Double_t pzMother = mother->Pz();
452 // within 0.1%
453 if((TMath::Abs(pxMother-pxSumDgs)/(TMath::Abs(pxMother)+1.e-13)) > 0.00001 &&
454 (TMath::Abs(pyMother-pySumDgs)/(TMath::Abs(pyMother)+1.e-13)) > 0.00001 &&
455 (TMath::Abs(pzMother-pzSumDgs)/(TMath::Abs(pzMother)+1.e-13)) > 0.00001) return kFALSE;
456
457 return kTRUE;
458}