filling of mc info from stack or AliAODMCParticles updated
authorjklay <jklay@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sun, 11 Oct 2009 07:33:14 +0000 (07:33 +0000)
committerjklay <jklay@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sun, 11 Oct 2009 07:33:14 +0000 (07:33 +0000)
PWG4/PartCorrDep/AliAnaElectron.cxx
PWG4/PartCorrDep/AliAnaElectron.h

index e0a9dda91a098d0905c7195ea5fde4c7bff32cb0..287ae3d05a4f3af31431cb7fee1fce9dece57b77 100755 (executable)
@@ -841,8 +841,6 @@ void  AliAnaElectron::MakeAnalysisFillHistograms()
 \r
   AliStack * stack = 0x0;\r
   TParticle * primary = 0x0;\r
 \r
   AliStack * stack = 0x0;\r
   TParticle * primary = 0x0;\r
-  TClonesArray * mcparticles0 = 0x0;\r
-  TClonesArray * mcparticles1 = 0x0;\r
   AliAODMCParticle * aodprimary = 0x0;\r
 \r
   Int_t ph1 = 0;  //photonic 1 count\r
   AliAODMCParticle * aodprimary = 0x0;\r
 \r
   Int_t ph1 = 0;  //photonic 1 count\r
@@ -851,26 +849,11 @@ void  AliAnaElectron::MakeAnalysisFillHistograms()
 \r
   if(IsDataMC()) {\r
     if(GetReader()->ReadStack()){\r
 \r
   if(IsDataMC()) {\r
     if(GetReader()->ReadStack()){\r
-      stack =  GetMCStack() ;\r
-      \r
+      stack =  GetMCStack() ;      \r
       if(!stack)\r
        printf("AliAnaElectron::MakeAnalysisFillHistograms() *** no stack ***: \n");\r
       \r
     }\r
       if(!stack)\r
        printf("AliAnaElectron::MakeAnalysisFillHistograms() *** no stack ***: \n");\r
       \r
     }\r
-    else if(GetReader()->ReadAODMCParticles()){\r
-      //Get the list of MC particles\r
-      mcparticles0 = GetReader()->GetAODMCParticles(0);\r
-      if(!mcparticles0 && GetDebug() > 0)     {\r
-       printf("AliAnaElectron::MakeAnalysisFillHistograms() -  Standard MCParticles not available!\n");\r
-      }\r
-      if(GetReader()->GetSecondInputAODTree()){\r
-       mcparticles1 = GetReader()->GetAODMCParticles(1);\r
-       if(!mcparticles1 && GetDebug() > 0)     {\r
-         printf("AliAnaElectron::MakeAnalysisFillHistograms() -  Second input MCParticles not available!\n");\r
-       }\r
-      }\r
-      \r
-    }\r
   }// is data and MC\r
 \r
   ////////////////////////////////////\r
   }// is data and MC\r
 \r
   ////////////////////////////////////\r
@@ -931,8 +914,8 @@ void  AliAnaElectron::MakeAnalysisFillHistograms()
       if(IsDataMC()) {\r
        //determine tagging efficiency & mis-tagging rate\r
        //using b-quarks from stack\r
       if(IsDataMC()) {\r
        //determine tagging efficiency & mis-tagging rate\r
        //using b-quarks from stack\r
-       Bool_t isTrueBjet = IsMcBJet(jet->Eta(), jet->Phi() ,stack);\r
-       Bool_t isTrueDjet = IsMcDJet(jet->Eta(), jet->Phi() ,stack);\r
+       Bool_t isTrueBjet = IsMcBJet(jet->Eta(), jet->Phi());\r
+       Bool_t isTrueDjet = IsMcDJet(jet->Eta(), jet->Phi());\r
        if (isTrueBjet && GetDebug() > 0) printf("== True Bjet==\n");\r
        if (isTrueDjet && GetDebug() > 0) printf("== True Charm-jet==\n");\r
        if (dvmJet && GetDebug() > 0)     printf("== found DVM jet==\n");\r
        if (isTrueBjet && GetDebug() > 0) printf("== True Bjet==\n");\r
        if (isTrueDjet && GetDebug() > 0) printf("== True Charm-jet==\n");\r
        if (dvmJet && GetDebug() > 0)     printf("== found DVM jet==\n");\r
@@ -1035,19 +1018,10 @@ void  AliAnaElectron::MakeAnalysisFillHistograms()
       if(IsDataMC()) {\r
        fhPtNPEleTTE->Fill(ptele,GetMCSource(mctag));\r
        if(ele->GetDetector() == "EMCAL") fhPtNPEleEMCAL->Fill(ptele,GetMCSource(mctag));\r
       if(IsDataMC()) {\r
        fhPtNPEleTTE->Fill(ptele,GetMCSource(mctag));\r
        if(ele->GetDetector() == "EMCAL") fhPtNPEleEMCAL->Fill(ptele,GetMCSource(mctag));\r
-       if(GetMCSource(mctag) == 1) {\r
-         if(GetReader()->ReadStack()) { //only done if we have the stack\r
-           TParticle* prim = stack->Particle(ele->GetLabel());\r
-           if(prim->GetMother(0)>=0) {\r
-             Int_t mpdg = 0;\r
-             TParticle *parent = stack->Particle(prim->GetMother(0));\r
-             if(parent) mpdg = parent->GetPdgCode();\r
-             \r
-             if ((TMath::Abs(mpdg) >500  && TMath::Abs(mpdg) <600 ) ||\r
-                 (TMath::Abs(mpdg) >5000 && TMath::Abs(mpdg) <6000 ) )\r
-               fhPtNPEBHadron->Fill(ptele,parent->Pt());\r
-           } //mother\r
-         } //stack\r
+       if(GetMCSource(mctag) == 1) { //it's a bottom electron, now\r
+                                     //get the parent's pt\r
+         Double_t ptbHadron = GetBParentPt(ele->GetLabel());\r
+         fhPtNPEBHadron->Fill(ptele,ptbHadron);\r
        } //mctag\r
       } //isdatamc\r
     } //!photonic\r
        } //mctag\r
       } //isdatamc\r
     } //!photonic\r
@@ -1096,125 +1070,101 @@ void  AliAnaElectron::MakeAnalysisFillHistograms()
        TVector3 tempVect(p[0],p[1],p[2]);\r
        if ( TMath::Abs(tempVect.Eta())>fJetEtaCut || tempVect.Phi() < fJetPhiMin || tempVect.Phi() > fJetPhiMax) continue;\r
        //Only store it if it has a b-quark within dR < 0.2 of jet axis ?\r
        TVector3 tempVect(p[0],p[1],p[2]);\r
        if ( TMath::Abs(tempVect.Eta())>fJetEtaCut || tempVect.Phi() < fJetPhiMin || tempVect.Phi() > fJetPhiMax) continue;\r
        //Only store it if it has a b-quark within dR < 0.2 of jet axis ?\r
-       if(IsMcBJet(tempVect.Eta(),tempVect.Phi(),stack)) {\r
+       if(IsMcBJet(tempVect.Eta(),tempVect.Phi())) {\r
          bjetVect[iCount].SetXYZ(p[0], p[1], p[2]);\r
          iCount++;\r
        }\r
       }\r
     }\r
          bjetVect[iCount].SetXYZ(p[0], p[1], p[2]);\r
          iCount++;\r
        }\r
       }\r
     }\r
-        \r
-    if(GetReader()->ReadStack()) {\r
-      for(Int_t ipart = 0; ipart < stack->GetNtrack(); ipart++) {\r
-       primary = stack->Particle(ipart);\r
-       TLorentzVector mom;\r
-       primary->Momentum(mom);\r
-       Bool_t in = GetFidutialCut()->IsInFidutialCut(mom,fCalorimeter);\r
-       if(primary->Pt() < GetMinPt()) continue;\r
-       if(!in) continue;\r
-\r
-       Int_t pdgcode = primary->GetPdgCode();\r
-       if(TMath::Abs(pdgcode) == 211 || TMath::Abs(pdgcode) == 321 || TMath::Abs(pdgcode) == 2212)\r
-         fhPtMCHadron->Fill(primary->Pt());\r
-\r
-       //we only care about electrons\r
-       if(TMath::Abs(pdgcode) != 11) continue;\r
-       //we only want TRACKABLE electrons (TPC 85-250cm)\r
-       if(primary->R() > 200.) continue;\r
-       //Ignore low pt electrons\r
-       if(primary->Pt() < 0.2) continue;\r
-       \r
-       //find out what the ancestry of this electron is\r
-       Int_t mctag = -1;\r
-       Int_t input = 0;\r
-       mctag = GetMCAnalysisUtils()->CheckOrigin(ipart,GetReader(),input);\r
-\r
-       if(GetMCSource(mctag)==1) { //bottom electron\r
-         //See if it is within dR < 0.4 of a bjet\r
-         for(Int_t ij = 0; ij < nPythiaGenJets; ij++) {\r
-           Double_t deta = primary->Eta() - bjetVect[ij].Eta();\r
-           Double_t dphi = primary->Phi() - bjetVect[ij].Phi();\r
-           Double_t dR = TMath::Sqrt(deta*deta + dphi*dphi);\r
-           if(dR < 0.4) {\r
-             fhMCBJetElePt->Fill(primary->Pt(),bjetVect[ij].Pt());\r
-           }\r
-         }\r
-       }\r
-\r
-       if(primary->GetMother(0)>=0) {\r
-         Int_t mpdg = 0;\r
-         TParticle *parent = stack->Particle(primary->GetMother(0));\r
-         if (parent) mpdg = parent->GetPdgCode();\r
-         \r
-         if ((TMath::Abs(mpdg) >500  && TMath::Abs(mpdg) <600 ) ||\r
-             (TMath::Abs(mpdg) >5000 && TMath::Abs(mpdg) <6000 ) )\r
-           {\r
-             fhMCBHadronElePt->Fill(primary->Pt(), parent->Pt()); \r
-           }\r
-       }\r
 \r
 \r
-       //CHECK THAT THIS IS CORRECTLY FILLED - SHOULD WE USE MCSOURCE HERE?\r
-       fhPtMCElectron->Fill(primary->Pt(),0);  //0 = unfiltered\r
-       fhPtMCElectron->Fill(primary->Pt(),GetMCSource(mctag));\r
+    Int_t nPart = GetNumAODMCParticles();\r
+    if(GetReader()->ReadStack()) nPart = stack->GetNtrack();\r
 \r
 \r
-       //fill ntuple\r
-       if(fWriteNtuple) {\r
-         fMCEleNtuple->Fill(mctag,primary->Pt(),primary->Phi(),primary->Eta(),primary->Vx(),primary->Vy(),primary->Vz());\r
-       }\r
+    for(Int_t ipart = 0; ipart < nPart; ipart++) {\r
 \r
 \r
-      } //stack loop\r
+      //All the variables we want from MC particles\r
+      Double_t px = 0.; Double_t py = 0.; Double_t pz = 0.; Double_t e = 0.;\r
+      Double_t vx = -999.; Double_t vy = -999.; Double_t vz = -999.; Double_t vt = -999.;\r
+      Int_t pdg = 0; Int_t mpdg = 0; Double_t mpt = 0.;\r
 \r
 \r
-    } else if(GetReader()->ReadAODMCParticles()) {\r
-      Int_t npart0 = mcparticles0->GetEntriesFast();\r
-      Int_t npart1 = 0;\r
-      if(mcparticles1) npart1 = mcparticles1->GetEntriesFast();\r
-      Int_t npart = npart0+npart1;\r
-      for(Int_t ipart = 0; ipart < npart; ipart++) {\r
-       if(ipart < npart0) aodprimary = (AliAODMCParticle*)mcparticles0->At(ipart);\r
-       else aodprimary = (AliAODMCParticle*)mcparticles1->At(ipart-npart0);\r
-       if(!aodprimary) {\r
-         printf("AliAnaElectron::MakeAnalysisFillHistograms() *** no primary ***:  label %d \n", ipart);\r
-         continue;\r
+      if(GetReader()->ReadStack()) {\r
+       primary = stack->Particle(ipart);\r
+       pdg = primary->GetPdgCode();\r
+       px = primary->Px(); py = primary->Py(); pz = primary->Pz(); e = primary->Energy();\r
+       vx = primary->Vx(); vy = primary->Vy(); vz = primary->Vz(); vt = primary->T();\r
+       if(primary->GetMother(0)>=0) {\r
+         TParticle *parent = stack->Particle(primary->GetMother(0));\r
+         if (parent) {\r
+           mpdg = parent->GetPdgCode();\r
+           mpt = parent->Pt();\r
+         }\r
        }\r
        }\r
+      } else if(GetReader()->ReadAODMCParticles()) {\r
+       aodprimary =  (AliAODMCParticle*)GetMCParticle(ipart);\r
+       pdg = aodprimary->GetPdgCode();\r
+       px = aodprimary->Px(); py = aodprimary->Py(); pz = aodprimary->Pz(); e = aodprimary->E();\r
+       vx = aodprimary->Xv(); vy = aodprimary->Yv(); vz = aodprimary->Zv(); vt = aodprimary->T();\r
+       Int_t parentId = aodprimary->GetMother();\r
+       if(parentId>=0) {\r
+         AliAODMCParticle *parent = (AliAODMCParticle*)GetMCParticle(parentId);\r
+         if (parent) {\r
+           mpdg = parent->GetPdgCode();\r
+           mpt = parent->Pt();\r
+         }\r
+       }       \r
+      }\r
 \r
 \r
-       Double_t mom[3] = {0.,0.,0.};\r
-       aodprimary->PxPyPz(mom);\r
-       TLorentzVector mom2(mom,0.);    \r
-       Bool_t in = GetFidutialCut()->IsInFidutialCut(mom2,fCalorimeter);\r
-       if(aodprimary->Pt() < GetMinPt()) continue;\r
-       if(!in) continue;\r
-\r
-       Int_t pdgcode = aodprimary->GetPdgCode();\r
-       if(TMath::Abs(pdgcode) == 211 || TMath::Abs(pdgcode) == 321 || TMath::Abs(pdgcode) == 2212)\r
-         fhPtMCHadron->Fill(aodprimary->Pt());\r
-\r
-       //we only care about electrons\r
-       if(TMath::Abs(pdgcode) != 11) continue;\r
-       //we only want TRACKABLE electrons (TPC 85-250cm)\r
-       Double_t radius = TMath::Sqrt(aodprimary->Xv()*aodprimary->Xv() + aodprimary->Yv()*aodprimary->Yv());\r
-       if(radius > 200.) continue;\r
-\r
-       //find out what the ancestry of this electron is\r
-       Int_t mctag = -1;\r
-       Int_t input = 0;\r
-       Int_t ival = ipart;\r
-       if(ipart > npart0) { ival -= npart0; input = 1;}\r
-       mctag = GetMCAnalysisUtils()->CheckOrigin(ival,GetReader(),input);\r
-\r
-       fhPtMCElectron->Fill(aodprimary->Pt(),0); //0 = unfiltered\r
-       fhPtMCElectron->Fill(aodprimary->Pt(),GetMCSource(mctag));\r
+      TLorentzVector mom(px,py,pz,e);\r
+      TLorentzVector pos(vx,vy,vz,vt);\r
+      Bool_t in = GetFidutialCut()->IsInFidutialCut(mom,fCalorimeter);\r
+      if(mom.Pt() < GetMinPt()) continue;\r
+      if(!in) continue;\r
+\r
+      if(TMath::Abs(pdg) == 211 || TMath::Abs(pdg) == 321 || TMath::Abs(pdg) == 2212)\r
+       fhPtMCHadron->Fill(mom.Pt());\r
+      \r
+      //we only care about electrons\r
+      if(TMath::Abs(pdg) != 11) continue;\r
+      //we only want TRACKABLE electrons (TPC 85-250cm)\r
+      if(pos.Rho() > 200.) continue;\r
+      //Ignore low pt electrons\r
+      if(mom.Pt() < 0.2) continue;\r
        \r
        \r
-       //fill ntuple\r
-       if(fWriteNtuple) {\r
-         fMCEleNtuple->Fill(mctag,aodprimary->Pt(),aodprimary->Phi(),aodprimary->Eta(),\r
-                            aodprimary->Xv(),aodprimary->Yv(),aodprimary->Zv());\r
+      //find out what the ancestry of this electron is\r
+      Int_t mctag = -1;\r
+      Int_t input = 0;\r
+      mctag = GetMCAnalysisUtils()->CheckOrigin(ipart,GetReader(),input);\r
+      \r
+      if(GetMCSource(mctag)==1) { //bottom electron\r
+       //See if it is within dR < 0.4 of a bjet\r
+       for(Int_t ij = 0; ij < nPythiaGenJets; ij++) {\r
+         Double_t deta = primary->Eta() - bjetVect[ij].Eta();\r
+         Double_t dphi = primary->Phi() - bjetVect[ij].Phi();\r
+         Double_t dR = TMath::Sqrt(deta*deta + dphi*dphi);\r
+         if(dR < 0.4) {\r
+           fhMCBJetElePt->Fill(primary->Pt(),bjetVect[ij].Pt());\r
+         }\r
        }\r
        }\r
+      }\r
 \r
 \r
-      } //AODMC particles\r
-    } //input type\r
-  } //pure MC kine histos\r
-    \r
+      if ((TMath::Abs(mpdg) >500  && TMath::Abs(mpdg) <600 ) ||\r
+         (TMath::Abs(mpdg) >5000 && TMath::Abs(mpdg) <6000 ) )\r
+       {\r
+         fhMCBHadronElePt->Fill(mom.Pt(), mpt); \r
+       }\r
+      //CHECK THAT THIS IS CORRECTLY FILLED - SHOULD WE USE MCSOURCE HERE?\r
+      fhPtMCElectron->Fill(mom.Pt(),0);  //0 = unfiltered\r
+      fhPtMCElectron->Fill(mom.Pt(),GetMCSource(mctag));\r
+      \r
+      //fill ntuple\r
+      if(fWriteNtuple) {\r
+       fMCEleNtuple->Fill(mctag,mom.Pt(),mom.Phi(),mom.Eta(),vx,vy,vz);\r
+      }\r
+    }\r
+  } //MC loop\r
+  \r
   //if(GetDebug() > 0) \r
   //if(GetDebug() > 0) \r
-    printf("\tAliAnaElectron::Photonic electron counts: ph1 %d, ph2 %d, Both %d\n",ph1,ph2,phB);\r
+  printf("\tAliAnaElectron::Photonic electron counts: ph1 %d, ph2 %d, Both %d\n",ph1,ph2,phB);\r
 }\r
 \r
 //__________________________________________________________________\r
 }\r
 \r
 //__________________________________________________________________\r
@@ -1551,7 +1501,6 @@ Bool_t AliAnaElectron::PhotonicV0(Int_t id)
 \r
   Bool_t itIS = kFALSE;\r
 \r
 \r
   Bool_t itIS = kFALSE;\r
 \r
-  Double_t massE = 0.000511;\r
   Double_t massEta = 0.547;\r
   Double_t massRho0 = 0.770;\r
   Double_t massOmega = 0.782;\r
   Double_t massEta = 0.547;\r
   Double_t massRho0 = 0.770;\r
   Double_t massOmega = 0.782;\r
@@ -1649,6 +1598,40 @@ Bool_t AliAnaElectron::CheckTrack(const AliAODTrack* track, const char* type)
 \r
 }\r
 \r
 \r
 }\r
 \r
+//__________________________________________________________________\r
+Double_t AliAnaElectron::GetBParentPt(Int_t ipart)\r
+{\r
+  //return MC B parent pt\r
+  if(GetReader()->ReadStack()) { //only done if we have the stack                                                                                               \r
+    AliStack* stack = GetMCStack();\r
+    if(!stack) {\r
+      printf("Problem getting stack\n");\r
+      return 0.;\r
+    }\r
+    TParticle* prim = stack->Particle(ipart);\r
+    if(prim->GetMother(0)>=0) {\r
+      Int_t mpdg = 0;\r
+      TParticle *parent = stack->Particle(prim->GetMother(0));\r
+      if(parent) mpdg = parent->GetPdgCode();\r
+\r
+      if ((TMath::Abs(mpdg) >500  && TMath::Abs(mpdg) <600 ) ||\r
+         (TMath::Abs(mpdg) >5000 && TMath::Abs(mpdg) <6000 ) )\r
+       return parent->Pt();\r
+    }\r
+  } else if(GetReader()->ReadAODMCParticles()){\r
+    AliAODMCParticle* prim = (AliAODMCParticle*)GetMCParticle(ipart);\r
+    if(prim->GetMother()>=0) {\r
+      Int_t mpdg = 0;\r
+      AliAODMCParticle* parent = (AliAODMCParticle*)GetMCParticle(prim->GetMother());\r
+      if(parent) mpdg = parent->GetPdgCode();\r
+      if ((TMath::Abs(mpdg) >500  && TMath::Abs(mpdg) <600 ) ||\r
+          (TMath::Abs(mpdg) >5000 && TMath::Abs(mpdg) <6000 ) )\r
+        return parent->Pt();\r
+    }\r
+  }\r
+  return 0.;\r
+}\r
+\r
 //__________________________________________________________________\r
 Int_t AliAnaElectron::GetMCSource(Int_t tag)\r
 {\r
 //__________________________________________________________________\r
 Int_t AliAnaElectron::GetMCSource(Int_t tag)\r
 {\r
@@ -1684,26 +1667,48 @@ Int_t AliAnaElectron::GetMCSource(Int_t tag)
 }\r
 \r
 //__________________________________________________________________\r
 }\r
 \r
 //__________________________________________________________________\r
-AliAODMCParticle* AliAnaElectron::GetMCParticle(Int_t part\r
+Int_t AliAnaElectron::GetNumAODMCParticles(\r
 {\r
 {\r
-  //Use the appropriate MC source to get the MC particle at position\r
-  //part\r
-\r
-  AliAODMCParticle* mcp = 0x0;\r
-  TParticle * primary = 0x0;\r
+  //Get the number of AliAODMCParticles, if any\r
+  Int_t num = 0;\r
 \r
 \r
-  AliStack * stack = 0x0;\r
   TClonesArray * mcparticles0 = 0x0;\r
   TClonesArray * mcparticles1 = 0x0;\r
 \r
   TClonesArray * mcparticles0 = 0x0;\r
   TClonesArray * mcparticles1 = 0x0;\r
 \r
-  if(GetReader()->ReadStack()){\r
-    stack =  GetMCStack() ;\r
-    \r
-    if(!stack)\r
-      printf("AliAnaElectron::MakeAnalysisFillHistograms() *** no stack ***: \n");\r
-    \r
+  if(GetReader()->ReadAODMCParticles()){\r
+    //Get the list of MC particles\r
+    //                                                                                                 \r
+    mcparticles0 = GetReader()->GetAODMCParticles(0);\r
+    if(!mcparticles0 && GetDebug() > 0) {\r
+      printf("AliAnaElectron::MakeAnalysisFillHistograms() -  Standard MCParticles not available!\n");\r
+    }\r
+    if(GetReader()->GetSecondInputAODTree()){\r
+      mcparticles1 = GetReader()->GetAODMCParticles(1);\r
+      if(!mcparticles1 && GetDebug() > 0) {\r
+        printf("AliAnaElectron::MakeAnalysisFillHistograms() -  Second input MCParticles not available!\n");\r
+      }\r
+    }\r
+\r
+    Int_t npart0 = mcparticles0->GetEntriesFast();\r
+    Int_t npart1 = 0;\r
+    if(mcparticles1) npart1 = mcparticles1->GetEntriesFast();\r
+    Int_t npart = npart0+npart1;\r
+    return npart;\r
+\r
   }\r
   }\r
-  else if(GetReader()->ReadAODMCParticles()){\r
+\r
+  return num;\r
+}\r
+//__________________________________________________________________\r
+AliAODMCParticle* AliAnaElectron::GetMCParticle(Int_t ipart) \r
+{\r
+  //Get the MC particle at position ipart\r
+\r
+  AliAODMCParticle* aodprimary = 0x0;\r
+  TClonesArray * mcparticles0 = 0x0;\r
+  TClonesArray * mcparticles1 = 0x0;\r
+\r
+  if(GetReader()->ReadAODMCParticles()){\r
     //Get the list of MC particles                                                                                                                           \r
     mcparticles0 = GetReader()->GetAODMCParticles(0);\r
     if(!mcparticles0 && GetDebug() > 0) {\r
     //Get the list of MC particles                                                                                                                           \r
     mcparticles0 = GetReader()->GetAODMCParticles(0);\r
     if(!mcparticles0 && GetDebug() > 0) {\r
@@ -1716,31 +1721,62 @@ AliAODMCParticle* AliAnaElectron::GetMCParticle(Int_t part)
       }\r
     }\r
 \r
       }\r
     }\r
 \r
-  }\r
+    Int_t npart0 = mcparticles0->GetEntriesFast();\r
+    Int_t npart1 = 0;\r
+    if(mcparticles1) npart1 = mcparticles1->GetEntriesFast();\r
+    if(ipart < npart0) aodprimary = (AliAODMCParticle*)mcparticles0->At(ipart);\r
+    else aodprimary = (AliAODMCParticle*)mcparticles1->At(ipart-npart0);\r
+    if(!aodprimary) {\r
+      printf("AliAnaElectron::GetMCParticle() *** no primary ***:  label %d \n", ipart);\r
+      return 0x0;\r
+    }\r
 \r
 \r
-  return mcp;\r
+  } else {\r
+    printf("AliAnaElectron::GetMCParticle() - Asked for AliAODMCParticle but we have a stack reader.\n");\r
+  }\r
+  return aodprimary;\r
 \r
 }\r
 \r
 //__________________________________________________________________\r
 \r
 }\r
 \r
 //__________________________________________________________________\r
-Bool_t  AliAnaElectron::IsMcBJet(Double_t eta, Double_t phi, AliStack* stack)\r
+Bool_t  AliAnaElectron::IsMcBJet(Double_t jeta, Double_t jphi)\r
 {\r
   //Check the jet eta,phi against that of the b-quark\r
   //to decide whether it is an MC B-jet\r
   Bool_t bjet=kFALSE;\r
 \r
   //      printf("MTH: McStack ,nparticles=%d \n", stack->GetNtrack() );\r
 {\r
   //Check the jet eta,phi against that of the b-quark\r
   //to decide whether it is an MC B-jet\r
   Bool_t bjet=kFALSE;\r
 \r
   //      printf("MTH: McStack ,nparticles=%d \n", stack->GetNtrack() );\r
+\r
+  AliStack* stack = 0x0;\r
   \r
   for(Int_t ipart = 0; ipart < 100; ipart++) {\r
 \r
   \r
   for(Int_t ipart = 0; ipart < 100; ipart++) {\r
 \r
-    TParticle* primary = stack->Particle(ipart);\r
-    if (!primary) continue;\r
-    Int_t pdgcode = primary->GetPdgCode();\r
-    if ( TMath::Abs(pdgcode) != 5) continue;\r
+    Double_t pphi = -999.;\r
+    Double_t peta = -999.;\r
+    Int_t pdg = 0;\r
+    if(GetReader()->ReadStack()) {\r
+      stack = GetMCStack();\r
+      if(!stack) {\r
+       printf("AliAnaElectron::IsMCBJet() *** no stack ***: \n");\r
+       return kFALSE;\r
+      }\r
+      TParticle* primary = stack->Particle(ipart);\r
+      if (!primary) continue;\r
+      pdg = primary->GetPdgCode();\r
+      pphi = primary->Phi();\r
+      peta = primary->Eta();\r
+    } else if(GetReader()->ReadAODMCParticles()) {\r
+      AliAODMCParticle* aodprimary = GetMCParticle(ipart);\r
+      if(!aodprimary) continue;\r
+      pdg = aodprimary->GetPdgCode();\r
+      pphi = aodprimary->Phi();\r
+      peta = aodprimary->Eta();\r
+    }\r
+    if ( TMath::Abs(pdg) != 5) continue;\r
     \r
     //      printf("MTH: IsMcBJet : %d, pdg=%d : pt=%f \n", ipart, pdgcode, primary->Pt());\r
     \r
     //      printf("MTH: IsMcBJet : %d, pdg=%d : pt=%f \n", ipart, pdgcode, primary->Pt());\r
-    Double_t dphi = phi - primary->Phi();\r
-    Double_t deta = eta - primary->Eta();\r
+    Double_t dphi = jphi - pphi;\r
+    Double_t deta = jeta - peta;\r
     Double_t dr = sqrt(deta*deta + dphi*dphi);\r
     \r
     if (dr < 0.2) {\r
     Double_t dr = sqrt(deta*deta + dphi*dphi);\r
     \r
     if (dr < 0.2) {\r
@@ -1754,30 +1790,48 @@ Bool_t  AliAnaElectron::IsMcBJet(Double_t eta, Double_t phi, AliStack* stack)
 }\r
 \r
 //__________________________________________________________________\r
 }\r
 \r
 //__________________________________________________________________\r
-Bool_t  AliAnaElectron::IsMcDJet(Double_t eta, Double_t phi, AliStack* stack)\r
+Bool_t  AliAnaElectron::IsMcDJet(Double_t jeta, Double_t jphi)\r
 {\r
 {\r
-\r
+  //Check if this jet is a charm jet\r
   Bool_t cjet=kFALSE;\r
 \r
   Bool_t cjet=kFALSE;\r
 \r
-  if(IsDataMC()) {\r
-\r
-    for(Int_t ipart = 0; ipart < 100; ipart++) {\r
+  AliStack* stack = 0x0;\r
 \r
 \r
+  for(Int_t ipart = 0; ipart < 100; ipart++) {\r
+    \r
+    Double_t pphi = -999.;\r
+    Double_t peta = -999.;\r
+    Int_t pdg = 0;\r
+    if(GetReader()->ReadStack()) {\r
+      stack = GetMCStack();\r
+      if(!stack) {\r
+       printf("AliAnaElectron::IsMCDJet() *** no stack ***: \n");\r
+       return kFALSE;\r
+      }\r
       TParticle* primary = stack->Particle(ipart);\r
       if (!primary) continue;\r
       TParticle* primary = stack->Particle(ipart);\r
       if (!primary) continue;\r
-      Int_t pdgcode = primary->GetPdgCode();\r
-      if ( TMath::Abs(pdgcode) != 4) continue;\r
+      pdg = primary->GetPdgCode();\r
+      pphi = primary->Phi();\r
+      peta = primary->Eta();\r
+    } else if(GetReader()->ReadAODMCParticles()) {\r
+      AliAODMCParticle* aodprimary = GetMCParticle(ipart);\r
+      if(!aodprimary) continue;\r
+      pdg = aodprimary->GetPdgCode();\r
+      pphi = aodprimary->Phi();\r
+      peta = aodprimary->Eta();\r
+    }\r
 \r
 \r
-      Double_t dphi = phi - primary->Phi();\r
-      Double_t deta = eta - primary->Eta();\r
-      Double_t dr = sqrt(deta*deta + dphi*dphi);\r
+    if ( TMath::Abs(pdg) != 4) continue;\r
 \r
 \r
-      if (dr < 0.2) {\r
-       cjet=kTRUE;\r
-       break;\r
-      }\r
+    Double_t dphi = jphi - pphi;\r
+    Double_t deta = jeta - peta;\r
+    Double_t dr = sqrt(deta*deta + dphi*dphi);\r
+    \r
+    if (dr < 0.2) {\r
+      cjet=kTRUE;\r
+      break;\r
     }\r
     }\r
-  }//mc\r
+  }\r
 \r
   return cjet;\r
 \r
 \r
   return cjet;\r
 \r
@@ -1821,7 +1875,7 @@ void AliAnaElectron::ReadHistograms(TList* outputList)
   // Histograms of this analsys are kept in the same list as other\r
   // analysis, recover the position of\r
   // the first one and then add the next                                                      \r
   // Histograms of this analsys are kept in the same list as other\r
   // analysis, recover the position of\r
   // the first one and then add the next                                                      \r
-  Int_t index = outputList->IndexOf(outputList->FindObject(GetAddedHistogramsStringToName()+"fh1pOverE"));\r
+  //Int_t index = outputList->IndexOf(outputList->FindObject(GetAddedHistogramsStringToName()+"fh1pOverE"));\r
 \r
   //Read histograms, must be in the same order as in\r
   //GetCreateOutputObject.                   \r
 \r
   //Read histograms, must be in the same order as in\r
   //GetCreateOutputObject.                   \r
index 4c4e3366caec4e45e8f220946b2f91b7cb79fc1a..705a403e4c1cd6e1c9be4ecf7a4a549004a9f945 100755 (executable)
@@ -56,8 +56,8 @@ public:
   //check if track has been flagged as a non-photonic or DVM electron\r
   //used with the jet tracks to tag bjets\r
   Bool_t CheckTrack(const AliAODTrack* track,const char* type);  \r
   //check if track has been flagged as a non-photonic or DVM electron\r
   //used with the jet tracks to tag bjets\r
   Bool_t CheckTrack(const AliAODTrack* track,const char* type);  \r
-  Bool_t IsMcBJet(Double_t x, Double_t y, AliStack* st);\r
-  Bool_t IsMcDJet(Double_t x, Double_t y, AliStack* st);\r
+  Bool_t IsMcBJet(Double_t x, Double_t y);\r
+  Bool_t IsMcDJet(Double_t x, Double_t y);\r
 \r
   void Print(const Option_t * opt)const;\r
   \r
 \r
   void Print(const Option_t * opt)const;\r
   \r
@@ -113,6 +113,10 @@ public:
   //Need a clean way to get the MC info.  An AliAODMCParticle object\r
   //is returned from whichever source we are operating on\r
   AliAODMCParticle* GetMCParticle(Int_t part);\r
   //Need a clean way to get the MC info.  An AliAODMCParticle object\r
   //is returned from whichever source we are operating on\r
   AliAODMCParticle* GetMCParticle(Int_t part);\r
+  //Get MC B Parent pt\r
+  Double_t GetBParentPt(Int_t label);\r
+  //Get Number of particles in AliAODMCParticle array, if it exists\r
+  Int_t GetNumAODMCParticles();\r
 \r
   private:\r
   TString  fCalorimeter;  //! Which detector? EMCAL or PHOS\r
 \r
   private:\r
   TString  fCalorimeter;  //! Which detector? EMCAL or PHOS\r