b5dc2602e52a38fac00a39856b4d8e59e173b8a9
[u/mrichter/AliRoot.git] / PWGDQ / dielectron / AliAnalysisTaskDielectronReadAODBranch.cxx
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
38 ClassImp(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 //________________________________________________________________________
82 AliAnalysisTaskDielectronReadAODBranch::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 //________________________________________________________________________
131 AliAnalysisTaskDielectronReadAODBranch::~AliAnalysisTaskDielectronReadAODBranch()
132 {
133         if(fOutput)   delete fOutput;
134         if(fobj)      delete fobj;
135         if(fobjMC)    delete fobjMC;
136 }
137
138 //________________________________________________________________________
139 void AliAnalysisTaskDielectronReadAODBranch::Init()
140 {
141         // Initialization
142         if(fDebug > 1) printf("AnalysisTaskReadAOD::Init() \n");
143         return;
144 }
145
146 //________________________________________________________________________
147 void 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 //________________________________________________________________________
221 void AliAnalysisTaskDielectronReadAODBranch::UserExec(Option_t */*option*/)
222 {
223         // Execute analysis for current event:
224         AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
225         if (!aod) return;
226         
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 //________________________________________________________________________
334 void AliAnalysisTaskDielectronReadAODBranch::Terminate(Option_t */*option*/)
335 {
336         if(fDebug > 1) printf("AliAnalysisTaskDielectronReadAODBranch::Terminate() \n");
337         return;
338 }
339
340 //___________________________________________________________________
341 Bool_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
350         //
351         if (!mcArray) return kFALSE;
352         
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 }