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 }