X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PWG4%2FPartCorrDep%2FAliAnaElectron.cxx;h=9d1e312b5fbd108e5105987e8e34d294315359b9;hb=8cfc6870dc38988186d5104c84476803ab89d750;hp=8399523384a691b07029284e8c19b38d6dc9b976;hpb=f5b89c0e407e9daacb934be482567363820cc25c;p=u%2Fmrichter%2FAliRoot.git diff --git a/PWG4/PartCorrDep/AliAnaElectron.cxx b/PWG4/PartCorrDep/AliAnaElectron.cxx index 8399523384a..9d1e312b5fb 100755 --- a/PWG4/PartCorrDep/AliAnaElectron.cxx +++ b/PWG4/PartCorrDep/AliAnaElectron.cxx @@ -25,10 +25,11 @@ // --- ROOT system --- #include +#include #include #include #include -#include +//#include //#include // --- Analysis system --- @@ -36,7 +37,7 @@ #include "AliCaloTrackReader.h" #include "AliMCAnalysisUtils.h" #include "AliAODCaloCluster.h" -#include "AliFidutialCut.h" +#include "AliFiducialCut.h" #include "AliAODTrack.h" #include "AliAODPid.h" #include "AliCaloPID.h" @@ -54,7 +55,7 @@ ClassImp(AliAnaElectron) //____________________________________________________________________________ AliAnaElectron::AliAnaElectron() : AliAnaPartCorrBaseClass(),fCalorimeter(""), - fpOverEmin(0.),fpOverEmax(0.),fResidualCut(0.), + fpOverEmin(0.),fpOverEmax(0.),fResidualCut(0.),fMinClusEne(0.), fDrCut(0.),fPairDcaCut(0.),fDecayLenCut(0.),fImpactCut(0.), fAssocPtCut(0.),fMassCut(0.),fSdcaCut(0.),fITSCut(0), fNTagTrkCut(0),fIPSigCut(0.),fJetEtaCut(0.3),fJetPhiMin(1.8),fJetPhiMax(2.9), @@ -62,7 +63,9 @@ AliAnaElectron::AliAnaElectron() //event QA histos fhImpactXY(0),fhRefMult(0),fhRefMult2(0), //matching checks - fh1pOverE(0),fh1EOverp(0),fh1dR(0),fh2EledEdx(0),fh2MatchdEdx(0),fh2dEtadPhi(0), + fh3pOverE(0),fh3EOverp(0),fh3pOverE2(0),fh3EOverp2(0),fh3pOverE3(0),fh3EOverp3(0), + fh2pOverE(0),fh2EOverp(0),fh2pOverE2(0),fh2EOverp2(0), + fh1dR(0),fh2EledEdx(0),fh2MatchdEdx(0),fh2dEtadPhi(0), fh2dEtadPhiMatched(0),fh2dEtadPhiUnmatched(0),fh2TrackPVsClusterE(0), fh2TrackPtVsClusterE(0),fh2TrackPhiVsClusterPhi(0),fh2TrackEtaVsClusterEta(0), //Photonic electron checks @@ -71,13 +74,17 @@ AliAnaElectron::AliAnaElectron() fhPtElectron(0),fhPhiElectron(0),fhEtaElectron(0), fhPtNPE(0),fhPhiNPE(0),fhEtaNPE(0), fhPtPE(0),fhPhiPE(0),fhEtaPE(0), + //for comparisons with tracking detectors + fhPtHadron(0),fhPtNPEleTPC(0),fhPtNPEleTPCTRD(0),fhPtNPEleTTE(0), + fhPtNPEleEMCAL(0), //DVM B-tagging fhDVMBtagCut1(0),fhDVMBtagCut2(0),fhDVMBtagCut3(0),fhDVMBtagQA1(0),fhDVMBtagQA2(0), fhDVMBtagQA3(0),fhDVMBtagQA4(0),fhDVMBtagQA5(0), //IPSig B-tagging fhIPSigBtagQA1(0),fhIPSigBtagQA2(0),fhTagJetPt1x4(0),fhTagJetPt2x3(0),fhTagJetPt3x2(0), + fhePlusTagJetPt1x4(0),fhePlusTagJetPt2x3(0),fhePlusTagJetPt3x2(0), //B-Jet histograms - fhJetType(0),fhBJetXsiFF(0),fhBJetPtFF(0),fhBJetEtaPhi(0), + fhJetType(0),fhLeadJetType(0),fhBJetXsiFF(0),fhBJetPtFF(0),fhBJetEtaPhi(0), fhNonBJetXsiFF(0),fhNonBJetPtFF(0),fhNonBJetEtaPhi(0), ///////////////////////////////////////////////////////////// //Histograms that rely on MC info (not filled for real data) @@ -85,11 +92,14 @@ AliAnaElectron::AliAnaElectron() //reco electrons from various sources fhPhiConversion(0),fhEtaConversion(0), //for comparisons with tracking detectors - fhPtHadron(0),fhPtNPEleTPC(0),fhPtNPEleTPCTRD(0),fhPtNPEleTTE(0), + fhPtTrack(0), + fhPtNPEBHadron(0), //for computing efficiency of B-jet tags - fhBJetPt1x4(0),fhBJetPt2x3(0),fhBJetPt3x2(0),fhDVMJet(0), + fhBJetPt1x4(0),fhBJetPt2x3(0),fhBJetPt3x2(0), + fhFakeJetPt1x4(0),fhFakeJetPt2x3(0),fhFakeJetPt3x2(0),fhDVMJet(0), //MC rate histograms/ntuple - fMCEleNtuple(0),fhMCBJetElePt(0),fhPtMCHadron(0),fhPtMCElectron(0) + fMCEleNtuple(0),fhMCBJetElePt(0),fhMCBHadronElePt(0),fhPtMCHadron(0),fhPtMCElectron(0), + fhMCXYConversion(0),fhMCRadPtConversion(0) { //default ctor @@ -101,7 +111,8 @@ AliAnaElectron::AliAnaElectron() //____________________________________________________________________________ AliAnaElectron::AliAnaElectron(const AliAnaElectron & g) : AliAnaPartCorrBaseClass(g),fCalorimeter(g.fCalorimeter), - fpOverEmin(g.fpOverEmin),fpOverEmax(g.fpOverEmax),fResidualCut(g.fResidualCut), + fpOverEmin(g.fpOverEmin),fpOverEmax(g.fpOverEmax), + fResidualCut(g.fResidualCut),fMinClusEne(g.fMinClusEne), fDrCut(g.fDrCut),fPairDcaCut(g.fPairDcaCut),fDecayLenCut(g.fDecayLenCut),fImpactCut(g.fImpactCut), fAssocPtCut(g.fAssocPtCut),fMassCut(g.fMassCut),fSdcaCut(g.fSdcaCut),fITSCut(g.fITSCut), fNTagTrkCut(g.fNTagTrkCut),fIPSigCut(g.fIPSigCut), @@ -110,7 +121,12 @@ AliAnaElectron::AliAnaElectron(const AliAnaElectron & g) //event QA histos fhImpactXY(g.fhImpactXY),fhRefMult(g.fhRefMult),fhRefMult2(g.fhRefMult2), //matching checks - fh1pOverE(g.fh1pOverE),fh1EOverp(g.fh1EOverp),fh1dR(g.fh1dR),fh2EledEdx(g.fh2EledEdx), + fh3pOverE(g.fh3pOverE),fh3EOverp(g.fh3EOverp), + fh3pOverE2(g.fh3pOverE2),fh3EOverp2(g.fh3EOverp2), + fh3pOverE3(g.fh3pOverE3),fh3EOverp3(g.fh3EOverp3), + fh2pOverE(g.fh2pOverE),fh2EOverp(g.fh2EOverp), + fh2pOverE2(g.fh2pOverE2),fh2EOverp2(g.fh2EOverp2), + fh1dR(g.fh1dR),fh2EledEdx(g.fh2EledEdx), fh2MatchdEdx(g.fh2MatchdEdx),fh2dEtadPhi(g.fh2dEtadPhi), fh2dEtadPhiMatched(g.fh2dEtadPhiMatched),fh2dEtadPhiUnmatched(g.fh2dEtadPhiUnmatched), fh2TrackPVsClusterE(g.fh2TrackPVsClusterE),fh2TrackPtVsClusterE(g.fh2TrackPtVsClusterE), @@ -121,6 +137,10 @@ AliAnaElectron::AliAnaElectron(const AliAnaElectron & g) fhPtElectron(g.fhPtElectron),fhPhiElectron(g.fhPhiElectron),fhEtaElectron(g.fhEtaElectron), fhPtNPE(g.fhPtNPE),fhPhiNPE(g.fhPhiNPE),fhEtaNPE(g.fhEtaNPE), fhPtPE(g.fhPtPE),fhPhiPE(g.fhPhiPE),fhEtaPE(g.fhEtaPE), + //for comparisons with tracking detectors + fhPtHadron(g.fhPtHadron),fhPtNPEleTPC(g.fhPtNPEleTPC), + fhPtNPEleTPCTRD(g.fhPtNPEleTPCTRD),fhPtNPEleTTE(g.fhPtNPEleTTE), + fhPtNPEleEMCAL(g.fhPtNPEleEMCAL), //DVM B-tagging fhDVMBtagCut1(g.fhDVMBtagCut1),fhDVMBtagCut2(g.fhDVMBtagCut2),fhDVMBtagCut3(g.fhDVMBtagCut3), fhDVMBtagQA1(g.fhDVMBtagQA1),fhDVMBtagQA2(g.fhDVMBtagQA2), @@ -128,24 +148,30 @@ AliAnaElectron::AliAnaElectron(const AliAnaElectron & g) //IPSig B-tagging fhIPSigBtagQA1(g.fhIPSigBtagQA1),fhIPSigBtagQA2(g.fhIPSigBtagQA2), fhTagJetPt1x4(g.fhTagJetPt1x4),fhTagJetPt2x3(g.fhTagJetPt2x3),fhTagJetPt3x2(g.fhTagJetPt3x2), + fhePlusTagJetPt1x4(g.fhePlusTagJetPt1x4),fhePlusTagJetPt2x3(g.fhePlusTagJetPt2x3), + fhePlusTagJetPt3x2(g.fhePlusTagJetPt3x2), //B-Jet histograms - fhJetType(g.fhJetType),fhBJetXsiFF(g.fhBJetXsiFF),fhBJetPtFF(g.fhBJetPtFF), - fhBJetEtaPhi(g.fhBJetEtaPhi),fhNonBJetXsiFF(g.fhNonBJetXsiFF),fhNonBJetPtFF(g.fhNonBJetPtFF), - fhNonBJetEtaPhi(g.fhNonBJetEtaPhi), + fhJetType(g.fhJetType),fhLeadJetType(g.fhLeadJetType),fhBJetXsiFF(g.fhBJetXsiFF), + fhBJetPtFF(g.fhBJetPtFF),fhBJetEtaPhi(g.fhBJetEtaPhi),fhNonBJetXsiFF(g.fhNonBJetXsiFF), + fhNonBJetPtFF(g.fhNonBJetPtFF),fhNonBJetEtaPhi(g.fhNonBJetEtaPhi), ///////////////////////////////////////////////////////////// //Histograms that rely on MC info (not filled for real data) fEleNtuple(g.fEleNtuple), //reco electrons from various sources fhPhiConversion(g.fhPhiConversion),fhEtaConversion(g.fhEtaConversion), //for comparisons with tracking detectors - fhPtHadron(g.fhPtHadron),fhPtNPEleTPC(g.fhPtNPEleTPC), - fhPtNPEleTPCTRD(g.fhPtNPEleTPCTRD),fhPtNPEleTTE(g.fhPtNPEleTTE), + fhPtTrack(g.fhPtTrack), + fhPtNPEBHadron(g.fhPtNPEBHadron), //for computing efficiency of B-jet tags fhBJetPt1x4(g.fhBJetPt1x4),fhBJetPt2x3(g.fhBJetPt2x3), - fhBJetPt3x2(g.fhBJetPt3x2),fhDVMJet(g.fhDVMJet), + fhBJetPt3x2(g.fhBJetPt3x2), + fhFakeJetPt1x4(g.fhFakeJetPt1x4),fhFakeJetPt2x3(g.fhBJetPt2x3), + fhFakeJetPt3x2(g.fhFakeJetPt3x2),fhDVMJet(g.fhDVMJet), //MC rate histograms/ntuple fMCEleNtuple(g.fMCEleNtuple),fhMCBJetElePt(g.fhMCBJetElePt), - fhPtMCHadron(g.fhPtMCHadron),fhPtMCElectron(g.fhPtMCElectron) + fhMCBHadronElePt(g.fhMCBHadronElePt), + fhPtMCHadron(g.fhPtMCHadron),fhPtMCElectron(g.fhPtMCElectron), + fhMCXYConversion(g.fhMCXYConversion),fhMCRadPtConversion(g.fhMCRadPtConversion) { // cpy ctor @@ -161,6 +187,7 @@ AliAnaElectron & AliAnaElectron::operator = (const AliAnaElectron & g) fpOverEmin = g.fpOverEmin; fpOverEmax = g.fpOverEmax; fResidualCut = g.fResidualCut; + fMinClusEne = g.fMinClusEne; fDrCut = g.fDrCut; fPairDcaCut = g.fPairDcaCut; fDecayLenCut = g.fDecayLenCut; @@ -180,8 +207,16 @@ AliAnaElectron & AliAnaElectron::operator = (const AliAnaElectron & g) fhRefMult = g.fhRefMult; fhRefMult2 = g.fhRefMult2; //matching checks - fh1pOverE = g.fh1pOverE; - fh1EOverp = g.fh1EOverp; + fh3pOverE = g.fh3pOverE; + fh3EOverp = g.fh3EOverp; + fh3pOverE2 = g.fh3pOverE2; + fh3EOverp2 = g.fh3EOverp2; + fh3pOverE3 = g.fh3pOverE3; + fh3EOverp3 = g.fh3EOverp3; + fh2pOverE = g.fh2pOverE; + fh2EOverp = g.fh2EOverp; + fh2pOverE2 = g.fh2pOverE2; + fh2EOverp2 = g.fh2EOverp2; fh1dR = g.fh1dR; fh2EledEdx = g.fh2EledEdx; fh2MatchdEdx = g.fh2MatchdEdx; @@ -205,6 +240,10 @@ AliAnaElectron & AliAnaElectron::operator = (const AliAnaElectron & g) fhPtPE = g.fhPtPE; fhPhiPE = g.fhPhiPE; fhEtaPE = g.fhEtaPE; + //for comparisons with tracking detectors + fhPtHadron = g.fhPtHadron; fhPtNPEleTPC = g.fhPtNPEleTPC; + fhPtNPEleTPCTRD = g.fhPtNPEleTPCTRD; fhPtNPEleTTE = g.fhPtNPEleTTE; + fhPtNPEleEMCAL = g.fhPtNPEleEMCAL; //DVM B-tagging fhDVMBtagCut1 = g.fhDVMBtagCut1; fhDVMBtagCut2 = g.fhDVMBtagCut2; @@ -220,8 +259,12 @@ AliAnaElectron & AliAnaElectron::operator = (const AliAnaElectron & g) fhTagJetPt1x4 = g.fhTagJetPt1x4; fhTagJetPt2x3 = g.fhTagJetPt2x3; fhTagJetPt3x2 = g.fhTagJetPt3x2; + fhePlusTagJetPt1x4 = g.fhePlusTagJetPt1x4; + fhePlusTagJetPt2x3 = g.fhePlusTagJetPt2x3; + fhePlusTagJetPt3x2 = g.fhePlusTagJetPt3x2; //B-Jet histograms fhJetType = g.fhJetType; + fhLeadJetType = g.fhLeadJetType; fhBJetXsiFF = g.fhBJetXsiFF; fhBJetPtFF = g.fhBJetPtFF; fhBJetEtaPhi = g.fhBJetEtaPhi; @@ -235,14 +278,19 @@ AliAnaElectron & AliAnaElectron::operator = (const AliAnaElectron & g) fhPhiConversion = g.fhPhiConversion; fhEtaConversion = g.fhEtaConversion; //for comparisons with tracking detectors - fhPtHadron = g.fhPtHadron; fhPtNPEleTPC = g.fhPtNPEleTPC; - fhPtNPEleTPCTRD = g.fhPtNPEleTPCTRD; fhPtNPEleTTE = g.fhPtNPEleTTE; + fhPtTrack = g.fhPtTrack; + fhPtNPEBHadron = g.fhPtNPEBHadron; //for computing efficiency of B-jet tags fhBJetPt1x4 = g.fhBJetPt1x4; fhBJetPt2x3 = g.fhBJetPt2x3; - fhBJetPt3x2 = g.fhBJetPt3x2; fhDVMJet = g.fhDVMJet; + fhBJetPt3x2 = g.fhBJetPt3x2; + fhFakeJetPt1x4 = g.fhFakeJetPt1x4; fhFakeJetPt2x3 = g.fhFakeJetPt2x3; + fhFakeJetPt3x2 = g.fhFakeJetPt3x2; fhDVMJet = g.fhDVMJet; //MC rate histograms/ntuple fMCEleNtuple = g.fMCEleNtuple; fhMCBJetElePt = g.fhMCBJetElePt; + fhMCBHadronElePt = g.fhMCBHadronElePt; fhPtMCHadron = g.fhPtMCHadron; fhPtMCElectron = g.fhPtMCElectron; + fhMCXYConversion = g.fhMCXYConversion; + fhMCRadPtConversion = g.fhMCRadPtConversion; return *this; @@ -264,9 +312,9 @@ TList * AliAnaElectron::GetCreateOutputObjects() TList * outputContainer = new TList() ; outputContainer->SetName("ElectronHistos") ; - Int_t nptbins = GetHistoNPtBins(); - Int_t nphibins = GetHistoNPhiBins(); - Int_t netabins = GetHistoNEtaBins(); + Int_t nptbins = GetHistoPtBins(); + Int_t nphibins = GetHistoPhiBins(); + Int_t netabins = GetHistoEtaBins(); Float_t ptmax = GetHistoPtMax(); Float_t phimax = GetHistoPhiMax(); Float_t etamax = GetHistoEtaMax(); @@ -276,16 +324,25 @@ TList * AliAnaElectron::GetCreateOutputObjects() //event QA fhImpactXY = new TH1F("hImpactXY","Impact parameter for all tracks",200,-10,10.); - fhRefMult = new TH1F("hRefMult" ,"refmult QA: " ,100,0,200); - fhRefMult2 = new TH1F("hRefMult2" ,"refmult2 QA: " ,100,0,200); + fhRefMult = new TH1F("hRefMult" ,"refmult QA: " ,5000,0,5000); + fhRefMult2 = new TH1F("hRefMult2" ,"refmult2 QA: " ,5000,0,5000); outputContainer->Add(fhImpactXY); outputContainer->Add(fhRefMult); outputContainer->Add(fhRefMult2); //matching checks - fh1pOverE = new TH1F("h1pOverE","EMCAL-TRACK matches p/E",200,0.,10.); - fh1EOverp = new TH1F("h1EOverp","EMCAL-TRACK matches E/p",200,0.,10.); + fh3pOverE = new TH3F("h3pOverE" ,"EMCAL-TRACK matches p/E",nptbins,ptmin,ptmax,200,0.,10.,30,0,30); + fh3EOverp = new TH3F("h3EOverp" ,"EMCAL-TRACK matches E/p",nptbins,ptmin,ptmax,200,0.,5. ,30,0,30); + fh3pOverE2 = new TH3F("h3pOverE_Trk","EMCAL-TRACK matches p/E",nptbins,ptmin,ptmax,200,0.,10.,30,0,30); + fh3EOverp2 = new TH3F("h3EOverp_Trk","EMCAL-TRACK matches E/p",nptbins,ptmin,ptmax,200,0.,5. ,30,0,30); + fh3pOverE3 = new TH3F("h3pOverE_Tpc","EMCAL-TRACK matches p/E",nptbins,ptmin,ptmax,200,0.,10.,30,0,30); + fh3EOverp3 = new TH3F("h3EOverp_Tpc","EMCAL-TRACK matches E/p",nptbins,ptmin,ptmax,200,0.,5. ,30,0,30); + fh2pOverE = new TH2F("h2pOverE" ,"EMCAL-TRACK matches p/E",nptbins,ptmin,ptmax,200,0.,10.); + fh2EOverp = new TH2F("h2EOverp" ,"EMCAL-TRACK matches E/p",nptbins,ptmin,ptmax,200,0.,5. ); + fh2pOverE2 = new TH2F("h2pOverE_Trk","EMCAL-TRACK matches p/E",nptbins,ptmin,ptmax,200,0.,10.); + fh2EOverp2 = new TH2F("h2EOverp_Trk","EMCAL-TRACK matches E/p",nptbins,ptmin,ptmax,200,0.,5. ); + fh1dR = new TH1F("h1dR","EMCAL-TRACK matches dR",300, 0.,TMath::Pi()); fh2EledEdx = new TH2F("h2EledEdx","dE/dx vs. p for electrons",200,0.,50.,200,0.,400.); fh2MatchdEdx = new TH2F("h2MatchdEdx","dE/dx vs. p for all matches",200,0.,50.,200,0.,400.); @@ -298,8 +355,16 @@ TList * AliAnaElectron::GetCreateOutputObjects() fh2TrackPhiVsClusterPhi = new TH2F("h2TrackPhiVsClusterPhi","h2TrackPhiVsClusterPhi",nphibins,phimin,phimax,nphibins,phimin,phimax); fh2TrackEtaVsClusterEta = new TH2F("h2TrackEtaVsClusterEta","h2TrackEtaVsClusterEta",netabins,etamin,etamax,netabins,etamin,etamax); - outputContainer->Add(fh1pOverE) ; - outputContainer->Add(fh1EOverp) ; + outputContainer->Add(fh3pOverE) ; + outputContainer->Add(fh3EOverp) ; + outputContainer->Add(fh3pOverE2) ; + outputContainer->Add(fh3EOverp2) ; + outputContainer->Add(fh3pOverE3) ; + outputContainer->Add(fh3EOverp3) ; + outputContainer->Add(fh2pOverE) ; + outputContainer->Add(fh2EOverp) ; + outputContainer->Add(fh2pOverE2) ; + outputContainer->Add(fh2EOverp2) ; outputContainer->Add(fh1dR) ; outputContainer->Add(fh2EledEdx) ; outputContainer->Add(fh2MatchdEdx) ; @@ -339,6 +404,26 @@ TList * AliAnaElectron::GetCreateOutputObjects() outputContainer->Add(fhPhiPE) ; outputContainer->Add(fhEtaPE) ; + //These histograms are mixed REAL/MC: + //Bins along y-axis are: + //0 - unfiltered (filled for both real and MC data) + //1 - bottom, 2 - charm, 3 - charm from bottom (MC only) + //4 - conversion, 5 - Dalitz, 6 - W and Z, 7 - junk/unknown (MC only) + //8 - misidentified (MC only) + + //histograms for comparison to tracking detectors + fhPtHadron = new TH2F("hPtHadron","Charged hadrons w/in EMCAL acceptance",nptbins,ptmin,ptmax,10,0,10); + fhPtNPEleTPC = new TH2F("hPtNPEleTPC","Non-phot. Electrons identified by TPC w/in EMCAL acceptance",nptbins,ptmin,ptmax,10,0,10); + fhPtNPEleTPCTRD = new TH2F("hPtNPEleTPCTRD","Non-phot. Electrons identified by TPC+TRD w/in EMCAL acceptance",nptbins,ptmin,ptmax,10,0,10); + fhPtNPEleTTE = new TH2F("hPtNPEleTTE","Non-phot. Electrons identified by TPC+TRD+EMCAL w/in EMCAL acceptance",nptbins,ptmin,ptmax,10,0,10); + fhPtNPEleEMCAL = new TH2F("hPtNPEleEMCAL","Non-phot. Electrons identified by EMCAL w/in EMCAL acceptance",nptbins,ptmin,ptmax,10,0,10); + + outputContainer->Add(fhPtHadron); + outputContainer->Add(fhPtNPEleTPC); + outputContainer->Add(fhPtNPEleTPCTRD); + outputContainer->Add(fhPtNPEleTTE); + outputContainer->Add(fhPtNPEleEMCAL); + //B-tagging fhDVMBtagCut1 = new TH2F("hdvmbtag_cut1","DVM B-tag result cut1", 10,0,10 ,nptbins,ptmin,ptmax); fhDVMBtagCut2 = new TH2F("hdvmbtag_cut2","DVM B-tag result cut2", 10,0,10 ,nptbins,ptmin,ptmax); @@ -361,18 +446,25 @@ TList * AliAnaElectron::GetCreateOutputObjects() //IPSig B-tagging fhIPSigBtagQA1 = new TH1F("hipsigbtag_qa1" ,"IPSig B-tag QA: # tag tracks", 20,0,20); fhIPSigBtagQA2 = new TH1F("hipsigbtag_qa2" ,"IPSig B-tag QA: IP significance", 200,-10.,10.); - fhTagJetPt1x4 = new TH1F("hTagJetPt1x4","tagged jet pT (1 track, ipSignif>4);p_{T}",1000,0.,100.); - fhTagJetPt2x3 = new TH1F("hTagJetPt2x3","tagged jet pT (2 track, ipSignif>3);p_{T}",1000,0.,100.); - fhTagJetPt3x2 = new TH1F("hTagJetPt3x2","tagged jet pT (3 track, ipSignif>2);p_{T}",1000,0.,100.); + fhTagJetPt1x4 = new TH1F("hTagJetPt1x4","tagged jet pT (1 track, ipSignif>4);p_{T}",300,0.,300.); + fhTagJetPt2x3 = new TH1F("hTagJetPt2x3","tagged jet pT (2 track, ipSignif>3);p_{T}",300,0.,300.); + fhTagJetPt3x2 = new TH1F("hTagJetPt3x2","tagged jet pT (3 track, ipSignif>2);p_{T}",300,0.,300.); + fhePlusTagJetPt1x4 = new TH1F("hePlusTagJetPt1x4","tagged eJet pT (1 track, ipSignif>4);p_{T}",300,0.,300.); + fhePlusTagJetPt2x3 = new TH1F("hePlusTagJetPt2x3","tagged eJet pT (2 track, ipSignif>3);p_{T}",300,0.,300.); + fhePlusTagJetPt3x2 = new TH1F("hePlusTagJetPt3x2","tagged eJet pT (3 track, ipSignif>2);p_{T}",300,0.,300.); outputContainer->Add(fhIPSigBtagQA1) ; outputContainer->Add(fhIPSigBtagQA2) ; outputContainer->Add(fhTagJetPt1x4); outputContainer->Add(fhTagJetPt2x3); outputContainer->Add(fhTagJetPt3x2); + outputContainer->Add(fhePlusTagJetPt1x4); + outputContainer->Add(fhePlusTagJetPt2x3); + outputContainer->Add(fhePlusTagJetPt3x2); //B-Jet histograms - fhJetType = new TH2F("hJetType","# jets passing each tag method vs jet pt",10,0,10,300,0.,300.); + fhJetType = new TH2F("hJetType","# jets passing each tag method vs jet pt",15,0,15,300,0.,300.); + fhLeadJetType = new TH2F("hLeadJetType","# leading jets passing each tag method vs jet pt",15,0,15,300,0.,300.); fhBJetXsiFF = new TH2F("hBJetXsiFF","B-jet #Xsi Frag. Fn.",100,0.,10.,300,0.,300.); fhBJetPtFF = new TH2F("hBJetPtFF","B-jet p_{T} Frag. Fn.",nptbins,ptmin,ptmax,300,0.,300.); fhBJetEtaPhi = new TH2F("hBJetEtaPhi","B-jet eta-phi distribution",netabins,etamin,etamax,nphibins,phimin,phimax); @@ -381,6 +473,7 @@ TList * AliAnaElectron::GetCreateOutputObjects() fhNonBJetEtaPhi = new TH2F("hNonBJetEtaPhi","Non B-jet eta-phi distribution",netabins,etamin,etamax,nphibins,phimin,phimax); outputContainer->Add(fhJetType); + outputContainer->Add(fhLeadJetType); outputContainer->Add(fhBJetXsiFF); outputContainer->Add(fhBJetPtFF); outputContainer->Add(fhBJetEtaPhi); @@ -399,7 +492,7 @@ TList * AliAnaElectron::GetCreateOutputObjects() //electrons from various MC sources fhPhiConversion = new TH2F("hPhiConversion","Conversion Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax); - fhEtaConversion = new TH2F("hEtaConversion","Conversion Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax); + fhEtaConversion = new TH2F("hEtaConversion","Conversion Electron eta vs. pT",nptbins,ptmin,ptmax,netabins,etamin,etamax); outputContainer->Add(fhPhiConversion); outputContainer->Add(fhEtaConversion); @@ -408,25 +501,27 @@ TList * AliAnaElectron::GetCreateOutputObjects() //4 - conversion, 5 - Dalitz, 6 - W and Z, 7 - junk/unknown, 8 - misidentified //histograms for comparison to tracking detectors - fhPtHadron = new TH2F("hPtHadron","Charged hadrons w/in EMCAL acceptance",nptbins,ptmin,ptmax,10,0,10); - fhPtNPEleTPC = new TH2F("hPtNPEleTPC","Non-phot. Electrons identified by TPC w/in EMCAL acceptance",nptbins,ptmin,ptmax,10,0,10); - fhPtNPEleTPCTRD = new TH2F("hPtNPEleTPCTRD","Non-phot. Electrons identified by TPC+TRD w/in EMCAL acceptance",nptbins,ptmin,ptmax,10,0,10); - fhPtNPEleTTE = new TH2F("hPtNPEleTTE","Non-phot. Electrons identified by TPC+TRD+EMCAL w/in EMCAL acceptance",nptbins,ptmin,ptmax,10,0,10); + fhPtTrack = new TH2F("hPtTrack","Track w/in EMCAL acceptance",nptbins,ptmin,ptmax,10,0,10); + fhPtNPEBHadron = new TH2F("hPtNPEBHadron","Non-phot. b-electrons (TPC+TRD+EMCAL) vs B-hadron pt w/in EMCAL acceptance",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); - outputContainer->Add(fhPtHadron); - outputContainer->Add(fhPtNPEleTPC); - outputContainer->Add(fhPtNPEleTPCTRD); - outputContainer->Add(fhPtNPEleTTE); + outputContainer->Add(fhPtTrack); + outputContainer->Add(fhPtNPEBHadron); //for computing efficiency of IPSig tag - fhBJetPt1x4 = new TH1F("hBJetPt1x4","tagged B-jet pT (1 track, ipSignif>4);p_{T}",1000,0.,100.); - fhBJetPt2x3 = new TH1F("hBJetPt2x3","tagged B-jet pT (2 track, ipSignif>3);p_{T}",1000,0.,100.); - fhBJetPt3x2 = new TH1F("hBJetPt3x2","tagged B-jet pT (3 track, ipSignif>2);p_{T}",1000,0.,100.); + fhBJetPt1x4 = new TH1F("hBJetPt1x4","tagged B-jet pT (1 track, ipSignif>4);p_{T}",300,0.,300.); + fhBJetPt2x3 = new TH1F("hBJetPt2x3","tagged B-jet pT (2 track, ipSignif>3);p_{T}",300,0.,300.); + fhBJetPt3x2 = new TH1F("hBJetPt3x2","tagged B-jet pT (3 track, ipSignif>2);p_{T}",300,0.,300.); + fhFakeJetPt1x4 = new TH1F("hFakeJetPt1x4","fake tagged B-jet pT (1 track, ipSignif>4);p_{T}",300,0.,300.); + fhFakeJetPt2x3 = new TH1F("hFakeJetPt2x3","fake tagged B-jet pT (2 track, ipSignif>3);p_{T}",300,0.,300.); + fhFakeJetPt3x2 = new TH1F("hFakeJetPt3x2","fake tagged B-jet pT (3 track, ipSignif>2);p_{T}",300,0.,300.); fhDVMJet = new TH2F("hDVM_algo","# DVM jets passing vs Mc-Bjet",10,0,10,300,0.,300.); outputContainer->Add(fhBJetPt1x4); outputContainer->Add(fhBJetPt2x3); outputContainer->Add(fhBJetPt3x2); + outputContainer->Add(fhFakeJetPt1x4); + outputContainer->Add(fhFakeJetPt2x3); + outputContainer->Add(fhFakeJetPt3x2); outputContainer->Add(fhDVMJet); //MC Only histograms @@ -438,61 +533,70 @@ TList * AliAnaElectron::GetCreateOutputObjects() } fhMCBJetElePt = new TH2F("hMCBJetElePt","MC B-jet pT vs. electron pT",300,0.,300.,300,0.,300.); + fhMCBHadronElePt = new TH2F("hMCBHadronElePt","MC B-hadron pT vs. electron pT",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); fhPtMCHadron = new TH1F("hPtMCHadron","MC Charged hadrons w/in EMCAL acceptance",nptbins,ptmin,ptmax); //Bins along y-axis are: 0 - unfiltered, 1 - bottom, 2 - charm, 3 - charm from bottom, //4 - conversion, 5 - Dalitz, 6 - W and Z, 7 - junk/unknown fhPtMCElectron = new TH2F("hPtMCElectron","MC electrons from various sources w/in EMCAL acceptance",nptbins,ptmin,ptmax,10,0,10); + fhMCXYConversion = new TH2F("hMCXYConversion","XvsY of conversion electrons",400,-400.,400.,400,-400.,400.); + fhMCRadPtConversion = new TH2F("hMCRadPtConversion","Radius vs pT of conversion electrons",200,0.,400.,nptbins,ptmin,ptmax); + outputContainer->Add(fhMCBJetElePt); + outputContainer->Add(fhMCBHadronElePt); outputContainer->Add(fhPtMCHadron); outputContainer->Add(fhPtMCElectron); + outputContainer->Add(fhMCXYConversion); + outputContainer->Add(fhMCRadPtConversion); }//Histos with MC //Save parameters used for analysis - TString parList ; //this will be list of parameters used for this analysis. - char onePar[500] ; - - sprintf(onePar,"--- AliAnaElectron ---\n") ; - parList+=onePar ; - sprintf(onePar,"fCalorimeter: %s\n",fCalorimeter.Data()) ; - parList+=onePar ; - sprintf(onePar,"fpOverEmin: %f\n",fpOverEmin) ; - parList+=onePar ; - sprintf(onePar,"fpOverEmax: %f\n",fpOverEmax) ; - parList+=onePar ; - sprintf(onePar,"fResidualCut: %f\n",fResidualCut) ; - parList+=onePar ; - sprintf(onePar,"---DVM Btagging\n"); - parList+=onePar ; - sprintf(onePar,"max IP-cut (e,h): %f\n",fImpactCut); - parList+=onePar ; - sprintf(onePar,"min ITS-hits: %d\n",fITSCut); - parList+=onePar ; - sprintf(onePar,"max dR (e,h): %f\n",fDrCut); - parList+=onePar ; - sprintf(onePar,"max pairDCA: %f\n",fPairDcaCut); - parList+=onePar ; - sprintf(onePar,"max decaylength: %f\n",fDecayLenCut); - parList+=onePar ; - sprintf(onePar,"min Associated Pt: %f\n",fAssocPtCut); - parList+=onePar ; - sprintf(onePar,"---IPSig Btagging\n"); - parList+=onePar ; - sprintf(onePar,"min tag track: %d\n",fNTagTrkCut); - parList+=onePar ; - sprintf(onePar,"min IP significance: %f\n",fIPSigCut); - parList+=onePar ; - - //Get parameters set in base class. - parList += GetBaseParametersList() ; - - //Get parameters set in FidutialCut class (not available yet) - //parlist += GetFidCut()->GetFidCutParametersList() - - TObjString *oString= new TObjString(parList) ; - outputContainer->Add(oString); +// TString parList ; //this will be list of parameters used for this analysis. +// char onePar[500] ; +// +// sprintf(onePar,"--- AliAnaElectron ---\n") ; +// parList+=onePar ; +// sprintf(onePar,"fCalorimeter: %s\n",fCalorimeter.Data()) ; +// parList+=onePar ; +// sprintf(onePar,"fpOverEmin: %f\n",fpOverEmin) ; +// parList+=onePar ; +// sprintf(onePar,"fpOverEmax: %f\n",fpOverEmax) ; +// parList+=onePar ; +// sprintf(onePar,"fResidualCut: %f\n",fResidualCut) ; +// parList+=onePar ; +// sprintf(onePar,"fMinClusEne: %f\n",fMinClusEne) ; +// parList+=onePar ; +// sprintf(onePar,"---DVM Btagging\n"); +// parList+=onePar ; +// sprintf(onePar,"max IP-cut (e,h): %f\n",fImpactCut); +// parList+=onePar ; +// sprintf(onePar,"min ITS-hits: %d\n",fITSCut); +// parList+=onePar ; +// sprintf(onePar,"max dR (e,h): %f\n",fDrCut); +// parList+=onePar ; +// sprintf(onePar,"max pairDCA: %f\n",fPairDcaCut); +// parList+=onePar ; +// sprintf(onePar,"max decaylength: %f\n",fDecayLenCut); +// parList+=onePar ; +// sprintf(onePar,"min Associated Pt: %f\n",fAssocPtCut); +// parList+=onePar ; +// sprintf(onePar,"---IPSig Btagging\n"); +// parList+=onePar ; +// sprintf(onePar,"min tag track: %d\n",fNTagTrkCut); +// parList+=onePar ; +// sprintf(onePar,"min IP significance: %f\n",fIPSigCut); +// parList+=onePar ; +// +// //Get parameters set in base class. +// parList += GetBaseParametersList() ; +// +// //Get parameters set in FiducialCut class (not available yet) +// //parlist += GetFidCut()->GetFidCutParametersList() +// +// TObjString *oString= new TObjString(parList) ; +// outputContainer->Add(oString); return outputContainer ; @@ -527,8 +631,9 @@ void AliAnaElectron::InitParameters() fCalorimeter = "EMCAL" ; fpOverEmin = 0.5; - fpOverEmax = 1.5; + fpOverEmax = 1.2; fResidualCut = 0.02; + fMinClusEne = 4.0; //DVM B-tagging fDrCut = 1.0; fPairDcaCut = 0.02; @@ -557,8 +662,6 @@ void AliAnaElectron::MakeAnalysisFillAOD() // Also fill some QA histograms // - TObjArray *cl = new TObjArray(); - Double_t bfield = 0.; if(GetReader()->GetDataType() != AliCaloTrackReader::kMC) bfield = GetReader()->GetBField(); @@ -567,7 +670,8 @@ void AliAnaElectron::MakeAnalysisFillAOD() printf("This class not yet implemented for PHOS\n"); abort(); } - cl = GetAODEMCAL(); + + TObjArray *cl = GetAODEMCAL(); //////////////////////////////////////////////// //Start from tracks and get associated clusters @@ -594,12 +698,11 @@ void AliAnaElectron::MakeAnalysisFillAOD() fhImpactXY->Fill(imp[0]); //JLK CHECK - AliESDtrack esdTrack(track); - Double_t tpcpid[AliPID::kSPECIES]; - esdTrack.GetTPCpid(tpcpid); - Double_t eProb = tpcpid[AliPID::kElectron]; - printf("<%d> ESD eProb = %2.2f\n",itrk,eProb); - + //AliESDtrack esdTrack(track); + //Double_t tpcpid[AliPID::kSPECIES]; + //esdTrack.GetTPCpid(tpcpid); + //Double_t eProb = tpcpid[AliPID::kElectron]; + //if(eProb > 0) printf("<%d> ESD eProb = %2.2f\n",itrk,eProb); AliAODPid* pid = (AliAODPid*) track->GetDetPid(); if(pid == 0) { @@ -617,9 +720,19 @@ void AliAnaElectron::MakeAnalysisFillAOD() Double_t teta = pos.Eta(); Double_t tmom = mom.Mag(); - TLorentzVector mom2(mom,0.); - Bool_t in = GetFidutialCut()->IsInFidutialCut(mom2,fCalorimeter) ; - if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillAOD() - Track pt %2.2f, phi %2.2f, eta %2.2f in fidutial cut %d\n",track->Pt(), track->Phi(), track->Eta(), in); + //TLorentzVector mom2(mom,0.); + Bool_t in = kFALSE; + if(mom.Phi()*180./TMath::Pi() > 80. && mom.Phi()*180./TMath::Pi() < 190. && + mom.Eta() > -0.7 && mom.Eta() < 0.7) in = kTRUE; + //Also check the track + if(track->Phi()*180./TMath::Pi() > 80. && track->Phi()*180./TMath::Pi() < 190. && + track->Eta() > -0.7 && track->Eta() < 0.7) in = kTRUE; + //////////////////////////// + //THIS HAS A MEM LEAK JLK 24-Oct-09 + //Bool_t in = GetFiducialCut()->IsInFiducialCut(mom2,fCalorimeter) ; + /////////////////////////// + if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillAOD() - Track(Extrap) pt %2.2f(%2.2f), phi %2.2f(%2.2f), eta %2.2f(%2.2f) in fiducial cut %d\n",track->Pt(), mom.Pt(), track->Phi(), mom.Phi(), track->Eta(),mom.Eta(), in); + if(mom.Pt() > GetMinPt() && in) { Double_t dEdx = pid->GetTPCsignal(); @@ -651,14 +764,33 @@ void AliAnaElectron::MakeAnalysisFillAOD() if(trkChgHad) fhPtHadron->Fill(track->Pt(),GetMCSource(tmctag)); if(tpcEle && !photonic) fhPtNPEleTPC->Fill(track->Pt(),GetMCSource(tmctag)); if(trkEle && !photonic) fhPtNPEleTPCTRD->Fill(track->Pt(),GetMCSource(tmctag)); + fhPtTrack->Fill(track->Pt(),GetMCSource(tmctag)); } Bool_t emcEle = kFALSE; //For tracks in EMCAL acceptance, pair them with all clusters //and fill the dEta vs dPhi for these pairs: + + Double_t minR = 99; + Double_t minPe =-1; + Double_t minEp =-1; + Double_t minMult = -1; + Double_t minPt = -1; + for(Int_t iclus = 0; iclus < ntot; iclus++) { AliAODCaloCluster * clus = (AliAODCaloCluster*) (cl->At(iclus)); if(!clus) continue; + + //As of 11-Oct-2009 + //only select "good" clusters + if (clus->GetNCells() < 2 ) continue; + if (clus->GetNCells() > 30 ) continue; + if (clus->E() < fMinClusEne ) continue; + if (clus->GetDispersion() > 1 ) continue; + if (clus->GetM20() > 0.4 ) continue; + if (clus->GetM02() > 0.4 ) continue; + if (clus->GetM20() < 0.03 ) continue; + if (clus->GetM02() < 0.03 ) continue; Double_t x[3]; clus->GetPosition(x); @@ -682,10 +814,21 @@ void AliAnaElectron::MakeAnalysisFillAOD() //Good match if(res < fResidualCut) { fh2dEtadPhiMatched->Fill(deta,dphi); + fh2MatchdEdx->Fill(track->P(),dEdx); iCluster = iclus; - - Int_t cmctag = -1; - + + Double_t energy = clus->E(); + if(energy > 0) pOverE = tmom/energy; + + if (res< minR) { + minR = res; + minPe = pOverE; + minEp = energy/tmom; + minMult = clus->GetNCells() ; + minPt = track->Pt(); + } + + Int_t cmctag = -1; if(IsDataMC()) { //Do you want the cluster or the track label? Int_t input = 0; @@ -696,31 +839,42 @@ void AliAnaElectron::MakeAnalysisFillAOD() if(fWriteNtuple) { fEleNtuple->Fill(tmctag,cmctag,track->Pt(),track->Phi(),track->Eta(),track->P(),clus->E(),deta,dphi,clus->GetNCells(),dEdx,pidProb,imp[0],imp[1]); } - - fh2MatchdEdx->Fill(track->P(),dEdx); - - Double_t energy = clus->E(); - if(energy > 0) pOverE = tmom/energy; - fh1pOverE->Fill(pOverE); - fh1EOverp->Fill(energy/tmom); - - Int_t mult = clus->GetNCells(); - if(mult < 2 && GetDebug() > 0) printf("Single digit cluster.\n"); - - ////////////////////////////// - //Electron cuts happen here!// - ////////////////////////////// - if(pOverE > fpOverEmin && pOverE < fpOverEmax) emcEle = kTRUE; + } else { fh2dEtadPhiUnmatched->Fill(deta,dphi); - } - - } //calocluster loop - + }//res cut + }//calo cluster loop + + fh3pOverE->Fill(minPt,minPe ,minMult); + fh3EOverp->Fill(minPt,minEp ,minMult); + if (trkEle) { + fh3pOverE2->Fill(minPt,minPe ,minMult); + fh3EOverp2->Fill(minPt,minEp ,minMult); + } + if (tpcEle) { + fh3pOverE3->Fill(minPt,minPe ,minMult); + fh3EOverp3->Fill(minPt,minEp ,minMult); + } + //new + if (tmctag>-1 && GetMCSource(tmctag)<8 ) { + fh2pOverE->Fill(minPt,minPe ); + fh2EOverp->Fill(minPt,minEp ); + if (trkEle) { + fh2pOverE2->Fill(minPt,minPe ); + fh2EOverp2->Fill(minPt,minEp ); + } + } + + ////////////////////////////// + //Electron cuts happen here!// + ////////////////////////////// + if(minPe > fpOverEmin && minPe < fpOverEmax) emcEle = kTRUE; + /////////////////////////// //Fill AOD with electrons// /////////////////////////// - if(emcEle || trkEle) { + //Take all emcal electrons, but the others only if pT < 10 GeV + if(emcEle || ( (tpcEle || trkEle) && track->Pt() < 10.) ) { //B-tagging if(GetDebug() > 1) printf("Found Electron - do b-tagging\n"); @@ -734,15 +888,17 @@ void AliAnaElectron::MakeAnalysisFillAOD() tr.SetLabel(track->GetLabel()); tr.SetCaloLabel(iCluster,-1); //sets the indices of the original caloclusters tr.SetTrackLabel(itrk,-1); //sets the indices of the original tracks - if(emcEle) //PID determined by EMCAL + if(emcEle) {//PID determined by EMCAL tr.SetDetector(fCalorimeter); - else + } else { tr.SetDetector("CTS"); //PID determined by CTS + } + if(GetReader()->GetAODCTSNormalInputEntries() <= itrk) tr.SetInputFileIndex(1); //Make this preserve sign of particle if(track->Charge() < 0) tr.SetPdg(11); //electron is 11 else tr.SetPdg(-11); //positron is -11 - Int_t btag = -1; + Int_t btag = 0; if(dvmbtag > 0) tr.SetBTagBit(btag,tr.kDVMTag0); if(dvmbtag > 1) tr.SetBTagBit(btag,tr.kDVMTag1); if(dvmbtag > 2) tr.SetBTagBit(btag,tr.kDVMTag2); @@ -782,8 +938,6 @@ void AliAnaElectron::MakeAnalysisFillHistograms() AliStack * stack = 0x0; TParticle * primary = 0x0; - TClonesArray * mcparticles0 = 0x0; - TClonesArray * mcparticles1 = 0x0; AliAODMCParticle * aodprimary = 0x0; Int_t ph1 = 0; //photonic 1 count @@ -792,26 +946,11 @@ void AliAnaElectron::MakeAnalysisFillHistograms() if(IsDataMC()) { if(GetReader()->ReadStack()){ - stack = GetMCStack() ; - + stack = GetMCStack() ; if(!stack) printf("AliAnaElectron::MakeAnalysisFillHistograms() *** no stack ***: \n"); } - else if(GetReader()->ReadAODMCParticles()){ - //Get the list of MC particles - mcparticles0 = GetReader()->GetAODMCParticles(0); - if(!mcparticles0 && GetDebug() > 0) { - printf("AliAnaElectron::MakeAnalysisFillHistograms() - Standard MCParticles not available!\n"); - } - if(GetReader()->GetSecondInputAODTree()){ - mcparticles1 = GetReader()->GetAODMCParticles(1); - if(!mcparticles1 && GetDebug() > 0) { - printf("AliAnaElectron::MakeAnalysisFillHistograms() - Second input MCParticles not available!\n"); - } - } - - } }// is data and MC //////////////////////////////////// @@ -841,7 +980,11 @@ void AliAnaElectron::MakeAnalysisFillHistograms() //For DVM jet tag, we will look for DVM electrons //For IPSig, we compute the IPSig for all tracks and if the //number passing is above the cut, it passes + Bool_t leadJet = kFALSE; + if (ijet==0) leadJet= kTRUE; Bool_t eJet = kFALSE; + Bool_t eJet2 = kFALSE; //electron triggered + Bool_t hadJet = kFALSE; //hadron triggered Bool_t dvmJet = kFALSE; Bool_t ipsigJet = kFALSE; TRefArray* rt = jet->GetRefTracks(); @@ -855,6 +998,8 @@ void AliAnaElectron::MakeAnalysisFillHistograms() if( GetIPSignificance(jetTrack, jet->Phi()) > 2.) trackCounter[3]++; Bool_t isNPE = CheckTrack(jetTrack,"NPE"); if(isNPE) eJet = kTRUE; + if ( isNPE && jetTrack->Pt()>10.0 ) eJet2 =kTRUE; + if (!isNPE && jetTrack->Pt()>10.0) hadJet =kTRUE; Bool_t isDVM = CheckTrack(jetTrack,"DVM"); if(isDVM) dvmJet = kTRUE; } @@ -863,24 +1008,40 @@ void AliAnaElectron::MakeAnalysisFillHistograms() if(trackCounter[2]>1) fhTagJetPt2x3->Fill(jet->Pt()); if(trackCounter[3]>2) fhTagJetPt3x2->Fill(jet->Pt()); + if(trackCounter[1]>0 && eJet) fhePlusTagJetPt1x4->Fill(jet->Pt()); + if(trackCounter[2]>1 && eJet) fhePlusTagJetPt2x3->Fill(jet->Pt()); + if(trackCounter[3]>2 && eJet) fhePlusTagJetPt3x2->Fill(jet->Pt()); + if(trackCounter[0] > fNTagTrkCut) ipsigJet = kTRUE; if(IsDataMC()) { //determine tagging efficiency & mis-tagging rate //using b-quarks from stack - Bool_t isTrueBjet = IsMcBJet(jet->Eta(), jet->Phi() ,stack); + Bool_t isTrueBjet = IsMcBJet(jet->Eta(), jet->Phi()); + Bool_t isTrueDjet = IsMcDJet(jet->Eta(), jet->Phi()); if (isTrueBjet && GetDebug() > 0) printf("== True Bjet==\n"); + if (isTrueDjet && GetDebug() > 0) printf("== True Charm-jet==\n"); if (dvmJet && GetDebug() > 0) printf("== found DVM jet==\n"); if(isTrueBjet && dvmJet) fhDVMJet->Fill(0.,jet->Pt()); // good tagged if(isTrueBjet && !dvmJet) fhDVMJet->Fill(1.,jet->Pt()); // missed tagged if(!isTrueBjet && dvmJet) fhDVMJet->Fill(2.,jet->Pt()); // fake tagged if(!isTrueBjet && !dvmJet) fhDVMJet->Fill(3.,jet->Pt()); // others + if(isTrueDjet && !isTrueBjet && dvmJet) fhDVMJet->Fill(4.,jet->Pt()); // charm-tagged + if(isTrueDjet && !isTrueBjet && !dvmJet) fhDVMJet->Fill(5.,jet->Pt()); // charm -not tagged + if(!(isTrueDjet||isTrueBjet ) && dvmJet) fhDVMJet->Fill(6.,jet->Pt()); // light flavor -tagged + if(!(isTrueDjet||isTrueBjet ) && !dvmJet) fhDVMJet->Fill(7.,jet->Pt()); // light flavor -not tagged + if(isTrueBjet && eJet && dvmJet) fhDVMJet->Fill(8.,jet->Pt()); // bjet with electron + if(isTrueBjet && !eJet && dvmJet) fhDVMJet->Fill(9.,jet->Pt()); // needs more thought if(isTrueBjet) { if(trackCounter[1]>0) fhBJetPt1x4->Fill(jet->Pt()); if(trackCounter[2]>1) fhBJetPt2x3->Fill(jet->Pt()); if(trackCounter[3]>2) fhBJetPt3x2->Fill(jet->Pt()); + } else { + if(trackCounter[1]>0) fhFakeJetPt1x4->Fill(jet->Pt()); + if(trackCounter[2]>1) fhFakeJetPt2x3->Fill(jet->Pt()); + if(trackCounter[3]>2) fhFakeJetPt3x2->Fill(jet->Pt()); } } @@ -894,10 +1055,25 @@ void AliAnaElectron::MakeAnalysisFillHistograms() if(dvmJet && ipsigJet && !eJet) fhJetType->Fill(6.,jet->Pt()); //dvm & ipsig if(dvmJet && ipsigJet && eJet) fhJetType->Fill(7.,jet->Pt()); //all if(dvmJet || ipsigJet || eJet) fhJetType->Fill(8.,jet->Pt()); //any of them + if(eJet ) fhJetType->Fill(9.,jet->Pt()); //any of them + if(dvmJet) fhJetType->Fill(10.,jet->Pt()); //any of them + if(eJet2 ) fhJetType->Fill(11.,jet->Pt()); //any of them + if(hadJet) fhJetType->Fill(12.,jet->Pt()); //any of them if(eJet || ipsigJet || dvmJet) fhBJetEtaPhi->Fill(jet->Eta(),jet->Phi()); else fhNonBJetEtaPhi->Fill(jet->Eta(),jet->Phi()); + //leading jets + if (leadJet) { + fhLeadJetType->Fill(0.,jet->Pt()); //all + if(eJet ) fhLeadJetType->Fill(1.,jet->Pt()); + if(eJet2 ) fhLeadJetType->Fill(2.,jet->Pt()); + if(hadJet ) fhLeadJetType->Fill(3.,jet->Pt()); + if(eJet && (dvmJet || ipsigJet) ) fhLeadJetType->Fill(4.,jet->Pt()); + if(eJet2 && (dvmJet || ipsigJet) ) fhLeadJetType->Fill(5.,jet->Pt()); + if(hadJet && (dvmJet || ipsigJet) ) fhLeadJetType->Fill(6.,jet->Pt()); + } + for(Int_t itrk = 0; itrk < ntrk; itrk++) { AliAODTrack* jetTrack = (AliAODTrack*)jet->GetTrack(itrk); Double_t xsi = TMath::Log(jet->Pt()/jetTrack->Pt()); @@ -958,8 +1134,17 @@ void AliAnaElectron::MakeAnalysisFillHistograms() //with Minv near a relevant resonance if(!photonic) { fhPtNPEleTTE->Fill(ptele,0); //0 = no MC info - if(IsDataMC()) fhPtNPEleTTE->Fill(ptele,GetMCSource(mctag)); - } + if(ele->GetDetector() == fCalorimeter) fhPtNPEleEMCAL->Fill(ptele,0); + if(IsDataMC()) { + fhPtNPEleTTE->Fill(ptele,GetMCSource(mctag)); + if(ele->GetDetector() == "EMCAL") fhPtNPEleEMCAL->Fill(ptele,GetMCSource(mctag)); + if(GetMCSource(mctag) == 1) { //it's a bottom electron, now + //get the parent's pt + Double_t ptbHadron = GetBParentPt(ele->GetLabel()); + fhPtNPEBHadron->Fill(ptele,ptbHadron); + } //mctag + } //isdatamc + } //!photonic //kept for historical reasons? fhPtElectron ->Fill(ptele); @@ -1005,112 +1190,112 @@ void AliAnaElectron::MakeAnalysisFillHistograms() TVector3 tempVect(p[0],p[1],p[2]); if ( TMath::Abs(tempVect.Eta())>fJetEtaCut || tempVect.Phi() < fJetPhiMin || tempVect.Phi() > fJetPhiMax) continue; //Only store it if it has a b-quark within dR < 0.2 of jet axis ? - if(IsMcBJet(tempVect.Eta(),tempVect.Phi(),stack)) { + if(IsMcBJet(tempVect.Eta(),tempVect.Phi())) { bjetVect[iCount].SetXYZ(p[0], p[1], p[2]); iCount++; } } } - - if(GetReader()->ReadStack()) { - for(Int_t ipart = 0; ipart < stack->GetNtrack(); ipart++) { - primary = stack->Particle(ipart); - TLorentzVector mom; - primary->Momentum(mom); - Bool_t in = GetFidutialCut()->IsInFidutialCut(mom,fCalorimeter); - if(primary->Pt() < GetMinPt()) continue; - if(!in) continue; - - Int_t pdgcode = primary->GetPdgCode(); - if(TMath::Abs(pdgcode) == 211 || TMath::Abs(pdgcode) == 321 || TMath::Abs(pdgcode) == 2212) - fhPtMCHadron->Fill(primary->Pt()); - - //we only care about electrons - if(TMath::Abs(pdgcode) != 11) continue; - //we only want TRACKABLE electrons (TPC 85-250cm) - if(primary->R() > 200.) continue; - //Ignore low pt electrons - if(primary->Pt() < 0.2) continue; - - //find out what the ancestry of this electron is - Int_t mctag = -1; - Int_t input = 0; - mctag = GetMCAnalysisUtils()->CheckOrigin(ipart,GetReader(),input); - - if(GetMCSource(mctag)==1) { //bottom electron - //See if it is within dR < 0.4 of a bjet - for(Int_t ij = 0; ij < nPythiaGenJets; ij++) { - Double_t deta = primary->Eta() - bjetVect[ij].Eta(); - Double_t dphi = primary->Phi() - bjetVect[ij].Phi(); - Double_t dR = TMath::Sqrt(deta*deta + dphi*dphi); - if(dR < 0.4) { - fhMCBJetElePt->Fill(primary->Pt(),bjetVect[ij].Pt()); - } - } - } - fhPtMCElectron->Fill(primary->Pt(),0); //0 = unfiltered - fhPtMCElectron->Fill(primary->Pt(),GetMCSource(mctag)); + Int_t nPart = GetNumAODMCParticles(); + if(GetReader()->ReadStack()) nPart = stack->GetNtrack(); - //fill ntuple - if(fWriteNtuple) { - fMCEleNtuple->Fill(mctag,primary->Pt(),primary->Phi(),primary->Eta(),primary->Vx(),primary->Vy(),primary->Vz()); - } + for(Int_t ipart = 0; ipart < nPart; ipart++) { - } //stack loop + //All the variables we want from MC particles + Double_t px = 0.; Double_t py = 0.; Double_t pz = 0.; Double_t e = 0.; + Double_t vx = -999.; Double_t vy = -999.; Double_t vz = -999.; Double_t vt = -999.; + Int_t pdg = 0; Int_t mpdg = 0; Double_t mpt = 0.; - } else if(GetReader()->ReadAODMCParticles()) { - Int_t npart0 = mcparticles0->GetEntriesFast(); - Int_t npart1 = 0; - if(mcparticles1) npart1 = mcparticles1->GetEntriesFast(); - Int_t npart = npart0+npart1; - for(Int_t ipart = 0; ipart < npart; ipart++) { - if(ipart < npart0) aodprimary = (AliAODMCParticle*)mcparticles0->At(ipart); - else aodprimary = (AliAODMCParticle*)mcparticles1->At(ipart-npart0); - if(!aodprimary) { - printf("AliAnaElectron::MakeAnalysisFillHistograms() *** no primary ***: label %d \n", ipart); - continue; + if(GetReader()->ReadStack()) { + primary = stack->Particle(ipart); + pdg = primary->GetPdgCode(); + px = primary->Px(); py = primary->Py(); pz = primary->Pz(); e = primary->Energy(); + vx = primary->Vx(); vy = primary->Vy(); vz = primary->Vz(); vt = primary->T(); + if(primary->GetMother(0)>=0) { + TParticle *parent = stack->Particle(primary->GetMother(0)); + if (parent) { + mpdg = parent->GetPdgCode(); + mpt = parent->Pt(); + } } + } else if(GetReader()->ReadAODMCParticles()) { + aodprimary = (AliAODMCParticle*)GetMCParticle(ipart); + pdg = aodprimary->GetPdgCode(); + px = aodprimary->Px(); py = aodprimary->Py(); pz = aodprimary->Pz(); e = aodprimary->E(); + vx = aodprimary->Xv(); vy = aodprimary->Yv(); vz = aodprimary->Zv(); vt = aodprimary->T(); + Int_t parentId = aodprimary->GetMother(); + if(parentId>=0) { + AliAODMCParticle *parent = (AliAODMCParticle*)GetMCParticle(parentId); + if (parent) { + mpdg = parent->GetPdgCode(); + mpt = parent->Pt(); + } + } + } - Double_t mom[3] = {0.,0.,0.}; - aodprimary->PxPyPz(mom); - TLorentzVector mom2(mom,0.); - Bool_t in = GetFidutialCut()->IsInFidutialCut(mom2,fCalorimeter); - if(aodprimary->Pt() < GetMinPt()) continue; - if(!in) continue; - - Int_t pdgcode = aodprimary->GetPdgCode(); - if(TMath::Abs(pdgcode) == 211 || TMath::Abs(pdgcode) == 321 || TMath::Abs(pdgcode) == 2212) - fhPtMCHadron->Fill(aodprimary->Pt()); - - //we only care about electrons - if(TMath::Abs(pdgcode) != 11) continue; - //we only want TRACKABLE electrons (TPC 85-250cm) - Double_t radius = TMath::Sqrt(aodprimary->Xv()*aodprimary->Xv() + aodprimary->Yv()*aodprimary->Yv()); - if(radius > 200.) continue; - - //find out what the ancestry of this electron is - Int_t mctag = -1; - Int_t input = 0; - Int_t ival = ipart; - if(ipart > npart0) { ival -= npart0; input = 1;} - mctag = GetMCAnalysisUtils()->CheckOrigin(ival,GetReader(),input); - - fhPtMCElectron->Fill(aodprimary->Pt(),0); //0 = unfiltered - fhPtMCElectron->Fill(aodprimary->Pt(),GetMCSource(mctag)); + TLorentzVector mom(px,py,pz,e); + TLorentzVector pos(vx,vy,vz,vt); + Bool_t in = kFALSE; + if(mom.Phi()*180./TMath::Pi() > 80. && mom.Phi()*180./TMath::Pi() < 190. && + mom.Eta() > -0.7 && mom.Eta() < 0.7) in = kTRUE; + ///////////////////////////////// + //THIS HAS A MEM LEAK JLK 24-Oct-09 + //Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter); + //////////////////////////////// + if(mom.Pt() < GetMinPt()) continue; + if(!in) continue; + + if(TMath::Abs(pdg) == 211 || TMath::Abs(pdg) == 321 || TMath::Abs(pdg) == 2212) + fhPtMCHadron->Fill(mom.Pt()); + + //we only care about electrons + if(TMath::Abs(pdg) != 11) continue; + //we only want TRACKABLE electrons (TPC 85-250cm) + if(pos.Rho() > 200.) continue; + //Ignore low pt electrons + if(mom.Pt() < 0.2) continue; - //fill ntuple - if(fWriteNtuple) { - fMCEleNtuple->Fill(mctag,aodprimary->Pt(),aodprimary->Phi(),aodprimary->Eta(), - aodprimary->Xv(),aodprimary->Yv(),aodprimary->Zv()); + //find out what the ancestry of this electron is + Int_t mctag = -1; + Int_t input = 0; + mctag = GetMCAnalysisUtils()->CheckOrigin(ipart,GetReader(),input); + + if(GetMCSource(mctag)==1) { //bottom electron + //See if it is within dR < 0.4 of a bjet + for(Int_t ij = 0; ij < nPythiaGenJets; ij++) { + Double_t deta = primary->Eta() - bjetVect[ij].Eta(); + Double_t dphi = primary->Phi() - bjetVect[ij].Phi(); + Double_t dR = TMath::Sqrt(deta*deta + dphi*dphi); + if(dR < 0.4) { + fhMCBJetElePt->Fill(primary->Pt(),bjetVect[ij].Pt()); + } } + } - } //AODMC particles - } //input type - } //pure MC kine histos - + if ((TMath::Abs(mpdg) >500 && TMath::Abs(mpdg) <600 ) || + (TMath::Abs(mpdg) >5000 && TMath::Abs(mpdg) <6000 ) ) + { + fhMCBHadronElePt->Fill(mom.Pt(), mpt); + } + //CHECK THAT THIS IS CORRECTLY FILLED - SHOULD WE USE MCSOURCE HERE? + fhPtMCElectron->Fill(mom.Pt(),0); //0 = unfiltered + fhPtMCElectron->Fill(mom.Pt(),GetMCSource(mctag)); + + if(GetMCSource(mctag) == 4) {//conversion + fhMCXYConversion->Fill(vx,vy); + fhMCRadPtConversion->Fill(TMath::Sqrt(vx*vx+vy*vy),mom.Pt()); + } + + //fill ntuple + if(fWriteNtuple) { + fMCEleNtuple->Fill(mctag,mom.Pt(),mom.Phi(),mom.Eta(),vx,vy,vz); + } + } + } //MC loop + //if(GetDebug() > 0) - printf("\tAliAnaElectron::Photonic electron counts: ph1 %d, ph2 %d, Both %d\n",ph1,ph2,phB); + printf("\tAliAnaElectron::Photonic electron counts: ph1 %d, ph2 %d, Both %d\n",ph1,ph2,phB); } //__________________________________________________________________ @@ -1200,9 +1385,6 @@ Double_t AliAnaElectron::ComputeSignDca(AliAODTrack *tr, AliAODTrack *tr2 , floa Double_t massE = 0.000511; Double_t massK = 0.493677; - Double_t bfield = 5.; //kG - if(GetReader()->GetDataType() != AliCaloTrackReader::kMC) bfield = GetReader()->GetBField(); - Double_t vertex[3] = {-999.,-999.,-999}; //vertex if(GetReader()->GetDataType() != AliCaloTrackReader::kMC) { GetReader()->GetVertex(vertex); //If only one file, get the vertex from there @@ -1219,11 +1401,15 @@ Double_t AliAnaElectron::ComputeSignDca(AliAODTrack *tr, AliAODTrack *tr2 , floa AliExternalTrackParam *param1 = new AliExternalTrackParam(tr); AliExternalTrackParam *param2 = new AliExternalTrackParam(tr2); + Double_t bfield[3]; + param1->GetBxByBz(bfield); + Double_t bz = param1->GetBz(); + Double_t xplane1 = 0.; Double_t xplane2 = 0.; - Double_t pairdca = param1->GetDCA(param2,bfield,xplane1,xplane2); + Double_t pairdca = param1->GetDCA(param2,bz,xplane1,xplane2); - param1->PropagateTo(xplane1,bfield); - param2->PropagateTo(xplane2,bfield); + param1->PropagateToBxByBz(xplane1,bfield); + param2->PropagateToBxByBz(xplane2,bfield); Int_t id1 = 0, id2 = 0; AliESDv0 bvertex(*param1,id1,*param2,id2); @@ -1241,6 +1427,8 @@ Double_t AliAnaElectron::ComputeSignDca(AliAODTrack *tr, AliAODTrack *tr2 , floa decayvector = secvtxpt - primV; //decay vector from PrimVtx Double_t decaylength = decayvector.Mag(); + printf("\t JLK pairDCA = %2.2f\n",pairdca); + if(GetDebug() > 0) { printf(">>ComputeSdca:: mom1=%f, mom2=%f \n", emomAtB.Perp(), hmomAtB.Perp() ); printf(">>ComputeSdca:: pairDCA=%f, length=%f \n", pairdca,decaylength ); @@ -1293,7 +1481,6 @@ Double_t AliAnaElectron::GetIPSignificance(AliAODTrack *tr, Double_t jetPhi) } Double_t significance=0; - Double_t magField = 0; Double_t maxD = 10000.; Double_t impPar[] = {0,0}; Double_t ipCov[]={0,0,0}; @@ -1306,7 +1493,9 @@ Double_t AliAnaElectron::GetIPSignificance(AliAODTrack *tr, Double_t jetPhi) AliVTrack* vTrack = (AliVTrack*)vEvent->GetTrack(trackIndex); if(!vTrack) return -99; AliESDtrack esdTrack(vTrack); - if(!esdTrack.PropagateToDCA(vv, magField, maxD, impPar, ipCov)) return -100; + Double_t bfield[3]; + esdTrack.GetBxByBz(bfield); + if(!esdTrack.PropagateToDCABxByBz(vv, bfield, maxD, impPar, ipCov)) return -100; if(ipCov[0]<0) return -101; Double_t Pxy[] = {esdTrack.Px(), esdTrack.Py()}; @@ -1317,6 +1506,7 @@ Double_t AliAnaElectron::GetIPSignificance(AliAODTrack *tr, Double_t jetPhi) Double_t cosTheta = TMath::Cos(jetPhi - phiIP); Double_t sign = cosTheta/TMath::Abs(cosTheta); significance = TMath::Abs(impPar[0])/TMath::Sqrt(ipCov[0])*sign; + printf("\t JLK significance = %2.2f\n",significance); //ip = fabs(impPar[0]); fhIPSigBtagQA2->Fill(significance); return significance; @@ -1358,8 +1548,10 @@ Bool_t AliAnaElectron::PhotonicPrim(const AliAODPWG4Particle* part) Bool_t itIS = kFALSE; Double_t massE = 0.000511; - Double_t bfield = 5.; //kG - if(GetReader()->GetDataType() != AliCaloTrackReader::kMC) bfield = GetReader()->GetBField(); + Double_t massEta = 0.547; + Double_t massRho0 = 0.770; + Double_t massOmega = 0.782; + Double_t massPhi = 1.020; Int_t pdg1 = part->GetPdg(); Int_t trackId = part->GetTrackLabel(0); @@ -1412,14 +1604,20 @@ Bool_t AliAnaElectron::PhotonicPrim(const AliAODPWG4Particle* part) fh1OpeningAngle->Fill(dphi); fh1MinvPhoton->Fill(mass); - if(mass < 0.1) { - if(GetDebug() > 0) printf("######PROBABLY A PHOTON\n"); - itIS = kTRUE; - } - + if(mass < 0.1 || + (mass > massEta-0.05 || mass < massEta+0.05) || + (mass > massRho0-0.05 || mass < massRho0+0.05) || + (mass > massOmega-0.05 || mass < massOmega+0.05) || + (mass > massPhi-0.05 || mass < massPhi+0.05)) + { + + if(GetDebug() > 0) printf("######PROBABLY A PHOTON\n"); + itIS = kTRUE; + } + //clean up delete param2; - + } delete param1; @@ -1435,6 +1633,11 @@ Bool_t AliAnaElectron::PhotonicV0(Int_t id) //invariant mass consistent with photon conversion Bool_t itIS = kFALSE; + + Double_t massEta = 0.547; + Double_t massRho0 = 0.770; + Double_t massOmega = 0.782; + Double_t massPhi = 1.020; //---Get V0s--- AliAODEvent *aod = (AliAODEvent*) GetReader()->GetInputEvent(); @@ -1449,7 +1652,11 @@ Bool_t AliAnaElectron::PhotonicV0(Int_t id) printf("## PhotonicV0() :: neg-id: %d, pos-id: %d, THIS id: %d\n", v0->GetNegID(), v0->GetPosID(), id); printf("## PhotonicV0() :: Minv(e,e): %f \n", v0->InvMass2Prongs(0,1,11,11) ); } - if (mass < 0.100) { + if (mass < 0.100 || + (mass > massEta-0.05 || mass < massEta+0.05) || + (mass > massRho0-0.05 || mass < massRho0+0.05) || + (mass > massOmega-0.05 || mass < massOmega+0.05) || + (mass > massPhi-0.05 || mass < massPhi+0.05)) { if ( id == v0->GetNegID() || id == v0->GetPosID()) { itIS=kTRUE; if(GetDebug() > 0) printf("## PhotonicV0() :: It's a conversion electron!!! \n" ); @@ -1467,14 +1674,15 @@ Bool_t AliAnaElectron::GetDCA(const AliAODTrack* track,Double_t impPar[2], Doubl //Once alice-off gets its act together and fixes the AOD, this //should become obsolete. - Double_t bfield = 5.; //kG Double_t maxD = 100000.; //max transverse IP if(GetReader()->GetDataType() != AliCaloTrackReader::kMC) { - bfield = GetReader()->GetBField(); AliVEvent* ve = (AliVEvent*)GetReader()->GetInputEvent(); AliVVertex *vv = (AliVVertex*)ve->GetPrimaryVertex(); AliESDtrack esdTrack(track); - Bool_t gotit = esdTrack.PropagateToDCA(vv,bfield,maxD,impPar,cov); + Double_t bfield[3]; + esdTrack.GetBxByBz(bfield); + Bool_t gotit = esdTrack.PropagateToDCABxByBz(vv,bfield,maxD,impPar,cov); + printf("\t JLK impPar = %2.2f\n",impPar[0]); return gotit; } @@ -1499,13 +1707,12 @@ Bool_t AliAnaElectron::CheckTrack(const AliAODTrack* track, const char* type) Int_t label = ele->GetTrackLabel(0); if(label != trackId) continue; //skip to the next one if they don't match - if(type=="DVM") { - if(ele->CheckBTagBit(ele->GetBtag(),AliAODPWG4Particle::kDVMTag0) || - ele->CheckBTagBit(ele->GetBtag(),AliAODPWG4Particle::kDVMTag1) || + if(strcmp(type,"DVM")==0) { + if(ele->CheckBTagBit(ele->GetBtag(),AliAODPWG4Particle::kDVMTag1) || ele->CheckBTagBit(ele->GetBtag(),AliAODPWG4Particle::kDVMTag2)) pass = kTRUE; - } else if (type=="NPE") { + } else if (strcmp(type,"NPE")==0) { Bool_t photonic = kFALSE; Bool_t photonic1 = kFALSE; @@ -1525,6 +1732,40 @@ Bool_t AliAnaElectron::CheckTrack(const AliAODTrack* track, const char* type) } +//__________________________________________________________________ +Double_t AliAnaElectron::GetBParentPt(Int_t ipart) +{ + //return MC B parent pt + if(GetReader()->ReadStack()) { //only done if we have the stack + AliStack* stack = GetMCStack(); + if(!stack) { + printf("Problem getting stack\n"); + return 0.; + } + TParticle* prim = stack->Particle(ipart); + if(prim->GetMother(0)>=0) { + Int_t mpdg = 0; + TParticle *parent = stack->Particle(prim->GetMother(0)); + if(parent) mpdg = parent->GetPdgCode(); + + if ((TMath::Abs(mpdg) >500 && TMath::Abs(mpdg) <600 ) || + (TMath::Abs(mpdg) >5000 && TMath::Abs(mpdg) <6000 ) ) + return parent->Pt(); + } + } else if(GetReader()->ReadAODMCParticles()){ + AliAODMCParticle* prim = (AliAODMCParticle*)GetMCParticle(ipart); + if(prim->GetMother()>=0) { + Int_t mpdg = 0; + AliAODMCParticle* parent = (AliAODMCParticle*)GetMCParticle(prim->GetMother()); + if(parent) mpdg = parent->GetPdgCode(); + if ((TMath::Abs(mpdg) >500 && TMath::Abs(mpdg) <600 ) || + (TMath::Abs(mpdg) >5000 && TMath::Abs(mpdg) <6000 ) ) + return parent->Pt(); + } + } + return 0.; +} + //__________________________________________________________________ Int_t AliAnaElectron::GetMCSource(Int_t tag) { @@ -1532,6 +1773,9 @@ Int_t AliAnaElectron::GetMCSource(Int_t tag) //the number returned is the bin along one axis of 2-d histograms in //which to fill this electron + //Do this first + if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)) return 4; + if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron)) { //Bottom if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEFromB)) return 1; @@ -1540,8 +1784,8 @@ Int_t AliAnaElectron::GetMCSource(Int_t tag) && !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEFromB)) return 2; //Charm from bottom else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEFromCFromB)) return 3; - //Conversion - else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)) return 4; + // //Conversion + //else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)) return 4; //Dalitz else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay) || GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) @@ -1560,24 +1804,116 @@ Int_t AliAnaElectron::GetMCSource(Int_t tag) } //__________________________________________________________________ -Bool_t AliAnaElectron::IsMcBJet(Double_t eta, Double_t phi, AliStack* stack) +Int_t AliAnaElectron::GetNumAODMCParticles() +{ + //Get the number of AliAODMCParticles, if any + Int_t num = 0; + + TClonesArray * mcparticles0 = 0x0; + TClonesArray * mcparticles1 = 0x0; + + if(GetReader()->ReadAODMCParticles()){ + //Get the list of MC particles + // + mcparticles0 = GetReader()->GetAODMCParticles(0); + if(!mcparticles0 && GetDebug() > 0) { + printf("AliAnaElectron::MakeAnalysisFillHistograms() - Standard MCParticles not available!\n"); + } + if(GetReader()->GetSecondInputAODTree()){ + mcparticles1 = GetReader()->GetAODMCParticles(1); + if(!mcparticles1 && GetDebug() > 0) { + printf("AliAnaElectron::MakeAnalysisFillHistograms() - Second input MCParticles not available!\n"); + } + } + + Int_t npart0 = mcparticles0->GetEntriesFast(); + Int_t npart1 = 0; + if(mcparticles1) npart1 = mcparticles1->GetEntriesFast(); + Int_t npart = npart0+npart1; + return npart; + + } + + return num; +} +//__________________________________________________________________ +AliAODMCParticle* AliAnaElectron::GetMCParticle(Int_t ipart) +{ + //Get the MC particle at position ipart + + AliAODMCParticle* aodprimary = 0x0; + TClonesArray * mcparticles0 = 0x0; + TClonesArray * mcparticles1 = 0x0; + + if(GetReader()->ReadAODMCParticles()){ + //Get the list of MC particles + mcparticles0 = GetReader()->GetAODMCParticles(0); + if(!mcparticles0 && GetDebug() > 0) { + printf("AliAnaElectron::MakeAnalysisFillHistograms() - Standard MCParticles not available!\n"); + } + if(GetReader()->GetSecondInputAODTree()){ + mcparticles1 = GetReader()->GetAODMCParticles(1); + if(!mcparticles1 && GetDebug() > 0) { + printf("AliAnaElectron::MakeAnalysisFillHistograms() - Second input MCParticles not available!\n"); + } + } + + Int_t npart0 = mcparticles0->GetEntriesFast(); + Int_t npart1 = 0; + if(mcparticles1) npart1 = mcparticles1->GetEntriesFast(); + if(ipart < npart0) aodprimary = (AliAODMCParticle*)mcparticles0->At(ipart); + else aodprimary = (AliAODMCParticle*)mcparticles1->At(ipart-npart0); + if(!aodprimary) { + printf("AliAnaElectron::GetMCParticle() *** no primary ***: label %d \n", ipart); + return 0x0; + } + + } else { + printf("AliAnaElectron::GetMCParticle() - Asked for AliAODMCParticle but we have a stack reader.\n"); + } + return aodprimary; + +} + +//__________________________________________________________________ +Bool_t AliAnaElectron::IsMcBJet(Double_t jeta, Double_t jphi) { //Check the jet eta,phi against that of the b-quark //to decide whether it is an MC B-jet Bool_t bjet=kFALSE; // printf("MTH: McStack ,nparticles=%d \n", stack->GetNtrack() ); + + AliStack* stack = 0x0; for(Int_t ipart = 0; ipart < 100; ipart++) { - TParticle* primary = stack->Particle(ipart); - if (!primary) continue; - Int_t pdgcode = primary->GetPdgCode(); - if ( TMath::Abs(pdgcode) != 5) continue; + Double_t pphi = -999.; + Double_t peta = -999.; + Int_t pdg = 0; + if(GetReader()->ReadStack()) { + stack = GetMCStack(); + if(!stack) { + printf("AliAnaElectron::IsMCBJet() *** no stack ***: \n"); + return kFALSE; + } + TParticle* primary = stack->Particle(ipart); + if (!primary) continue; + pdg = primary->GetPdgCode(); + pphi = primary->Phi(); + peta = primary->Eta(); + } else if(GetReader()->ReadAODMCParticles()) { + AliAODMCParticle* aodprimary = GetMCParticle(ipart); + if(!aodprimary) continue; + pdg = aodprimary->GetPdgCode(); + pphi = aodprimary->Phi(); + peta = aodprimary->Eta(); + } + if ( TMath::Abs(pdg) != 5) continue; // printf("MTH: IsMcBJet : %d, pdg=%d : pt=%f \n", ipart, pdgcode, primary->Pt()); - Double_t dphi = phi - primary->Phi(); - Double_t deta = eta - primary->Eta(); + Double_t dphi = jphi - pphi; + Double_t deta = jeta - peta; Double_t dr = sqrt(deta*deta + dphi*dphi); if (dr < 0.2) { @@ -1590,6 +1926,54 @@ Bool_t AliAnaElectron::IsMcBJet(Double_t eta, Double_t phi, AliStack* stack) } +//__________________________________________________________________ +Bool_t AliAnaElectron::IsMcDJet(Double_t jeta, Double_t jphi) +{ + //Check if this jet is a charm jet + Bool_t cjet=kFALSE; + + AliStack* stack = 0x0; + + for(Int_t ipart = 0; ipart < 100; ipart++) { + + Double_t pphi = -999.; + Double_t peta = -999.; + Int_t pdg = 0; + if(GetReader()->ReadStack()) { + stack = GetMCStack(); + if(!stack) { + printf("AliAnaElectron::IsMCDJet() *** no stack ***: \n"); + return kFALSE; + } + TParticle* primary = stack->Particle(ipart); + if (!primary) continue; + pdg = primary->GetPdgCode(); + pphi = primary->Phi(); + peta = primary->Eta(); + } else if(GetReader()->ReadAODMCParticles()) { + AliAODMCParticle* aodprimary = GetMCParticle(ipart); + if(!aodprimary) continue; + pdg = aodprimary->GetPdgCode(); + pphi = aodprimary->Phi(); + peta = aodprimary->Eta(); + } + + if ( TMath::Abs(pdg) != 4) continue; + + Double_t dphi = jphi - pphi; + Double_t deta = jeta - peta; + Double_t dr = sqrt(deta*deta + dphi*dphi); + + if (dr < 0.2) { + cjet=kTRUE; + break; + } + } + + return cjet; + +} + //__________________________________________________________________ void AliAnaElectron::Print(const Option_t * opt) const { @@ -1619,7 +2003,7 @@ void AliAnaElectron::Print(const Option_t * opt) const } //________________________________________________________________________ -void AliAnaElectron::ReadHistograms(TList* outputList) +void AliAnaElectron::ReadHistograms(TList* /* outputList */) { // Needed when Terminate is executed in distributed environment // Refill analysis histograms of this class with corresponding @@ -1628,14 +2012,14 @@ void AliAnaElectron::ReadHistograms(TList* outputList) // Histograms of this analsys are kept in the same list as other // analysis, recover the position of // the first one and then add the next - Int_t index = outputList->IndexOf(outputList->FindObject(GetAddedHistogramsStringToName()+"fh1pOverE")); + //Int_t index = outputList->IndexOf(outputList->FindObject(GetAddedHistogramsStringToName()+"fh1pOverE")); //Read histograms, must be in the same order as in //GetCreateOutputObject. - fh1pOverE = (TH1F *) outputList->At(index); - fh1dR = (TH1F *) outputList->At(index++); - fh2EledEdx = (TH2F *) outputList->At(index++); - fh2MatchdEdx = (TH2F *) outputList->At(index++); + //fh1pOverE = (TH1F *) outputList->At(index); + //fh1dR = (TH1F *) outputList->At(index++); + //fh2EledEdx = (TH2F *) outputList->At(index++); + //fh2MatchdEdx = (TH2F *) outputList->At(index++); }