-/**************************************************************************
- * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- * *
- * Author: The ALICE Off-line Project. *
- * Contributors are mentioned in the code where appropriate. *
- * *
- * Permission to use, copy, modify and distribute this software and its *
- * documentation strictly for non-commercial purposes hereby granted *
- * without fee, provided that the above copyright notice appears in all *
- * copies and that both the copyright notice and this permission notice *
- * appear in the supporting documentation. The authors make no claims *
- * about the suitability of this software for any purpose. It is *
- * provided "as is" without express or implied warranty. *
- **************************************************************************/
-/* $Id: $ */
-
-//_________________________________________________________________________
-//
-// Class for the electron identification.
-// Clusters from EMCAL matched to tracks
-// and kept in the AOD. Few histograms produced.
-//
-// -- Author: J.L. Klay (Cal Poly), M. Heinz (Yale)
-//////////////////////////////////////////////////////////////////////////////
-
-// --- ROOT system ---
-#include <TH2F.h>
-#include <TParticle.h>
-#include <TNtuple.h>
-#include <TClonesArray.h>
-#include <TObjString.h>
-//#include <Riostream.h>
-
-// --- Analysis system ---
-#include "AliAnaElectron.h"
-#include "AliCaloTrackReader.h"
-#include "AliMCAnalysisUtils.h"
-#include "AliAODCaloCluster.h"
-#include "AliFidutialCut.h"
-#include "AliAODTrack.h"
-#include "AliAODPid.h"
-#include "AliCaloPID.h"
-#include "AliAODMCParticle.h"
-#include "AliStack.h"
-#include "AliExternalTrackParam.h"
-#include "AliESDv0.h"
-
-ClassImp(AliAnaElectron)
-
-//____________________________________________________________________________
-AliAnaElectron::AliAnaElectron()
-: AliAnaPartCorrBaseClass(),fCalorimeter(""),
- fpOverEmin(0.),fpOverEmax(0.),fResidualCut(0.),
- fDrCut(0.),fPairDcaCut(0.),fDecayLenCut(0.),fImpactCut(0.),
- fAssocPtCut(0.),fMassCut(0.),fSdcaCut(0.),fITSCut(0),
- fWriteNtuple(kFALSE),
- //matching checks
- fEleNtuple(0),
- fh1pOverE(0),fh1dR(0),fh2EledEdx(0),fh2MatchdEdx(0),fh2dEtadPhi(0),
- fh2dEtadPhiMatched(0),fh2dEtadPhiUnmatched(0),
- fh2TrackPVsClusterE(0),fh2TrackPtVsClusterE(0),fh2TrackPhiVsClusterPhi(0),fh2TrackEtaVsClusterEta(0),
- //Photonic electron checks
- fh1OpeningAngle(0),fh1MinvPhoton(0),
- //reco
- fhPtElectron(0),fhPhiElectron(0),fhEtaElectron(0),
- fhPtNPE(0),fhPhiNPE(0),fhEtaNPE(0),
- fhPtPE(0),fhPhiPE(0),fhEtaPE(0),
- fhPtConversion(0),fhPhiConversion(0),fhEtaConversion(0),
- fhPtBottom(0),fhPhiBottom(0),fhEtaBottom(0),
- fhPtCharm(0),fhPhiCharm(0),fhEtaCharm(0),
- fhPtCFromB(0),fhPhiCFromB(0),fhEtaCFromB(0),
- fhPtDalitz(0),fhPhiDalitz(0),fhEtaDalitz(0),
- fhPtWDecay(0),fhPhiWDecay(0),fhEtaWDecay(0),
- fhPtZDecay(0),fhPhiZDecay(0),fhEtaZDecay(0),
- fhPtPrompt(0),fhPhiPrompt(0),fhEtaPrompt(0),
- fhPtUnknown(0),fhPhiUnknown(0),fhEtaUnknown(0),
- fhPtHadron(0),fhPtEleTrkDet(0),
- //B-tagging
- fhBtagCut1(0),fhBtagCut2(0),fhBtagCut3(0),fhBtagQA1(0),fhBtagQA2(0),
- //MC
- fMCEleNtuple(0),fhPtMCHadron(0)
-{
- //default ctor
-
- //Initialize parameters
- InitParameters();
-
-}
-
-//____________________________________________________________________________
-AliAnaElectron::AliAnaElectron(const AliAnaElectron & g)
- : AliAnaPartCorrBaseClass(g),fCalorimeter(g.fCalorimeter),
- fpOverEmin(g.fpOverEmin),fpOverEmax(g.fpOverEmax),fResidualCut(g.fResidualCut),
- fDrCut(g.fDrCut),fPairDcaCut(g.fPairDcaCut),fDecayLenCut(g.fDecayLenCut),fImpactCut(g.fImpactCut),
- fAssocPtCut(g.fAssocPtCut),fMassCut(g.fMassCut),fSdcaCut(g.fSdcaCut),fITSCut(g.fITSCut),
- fWriteNtuple(g.fWriteNtuple),
- //matching checks
- fEleNtuple(g.fEleNtuple),
- fh1pOverE(g.fh1pOverE),fh1dR(g.fh1dR),
- fh2EledEdx(g.fh2EledEdx),fh2MatchdEdx(g.fh2MatchdEdx),fh2dEtadPhi(g.fh2dEtadPhi),
- fh2dEtadPhiMatched(g.fh2dEtadPhiMatched),fh2dEtadPhiUnmatched(g.fh2dEtadPhiUnmatched),
- fh2TrackPVsClusterE(g.fh2TrackPVsClusterE),fh2TrackPtVsClusterE(g.fh2TrackPtVsClusterE),
- fh2TrackPhiVsClusterPhi(g.fh2TrackPhiVsClusterPhi),fh2TrackEtaVsClusterEta(g.fh2TrackEtaVsClusterEta),
- //Photonic electron checks
- fh1OpeningAngle(g.fh1OpeningAngle),fh1MinvPhoton(g.fh1MinvPhoton),
- //reco
- 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),
- fhPtConversion(g.fhPtConversion),fhPhiConversion(g.fhPhiConversion),fhEtaConversion(g.fhEtaConversion),
- fhPtBottom(g.fhPtBottom),fhPhiBottom(g.fhPhiBottom),fhEtaBottom(g.fhEtaBottom),
- fhPtCharm(g.fhPtCharm),fhPhiCharm(g.fhPhiCharm),fhEtaCharm(g.fhEtaCharm),
- fhPtCFromB(g.fhPtCFromB),fhPhiCFromB(g.fhPhiCFromB),fhEtaCFromB(g.fhEtaCFromB),
- fhPtDalitz(g.fhPtDalitz),fhPhiDalitz(g.fhPhiDalitz),fhEtaDalitz(g.fhEtaDalitz),
- fhPtWDecay(g.fhPtWDecay),fhPhiWDecay(g.fhPhiWDecay),fhEtaWDecay(g.fhEtaWDecay),
- fhPtZDecay(g.fhPtZDecay),fhPhiZDecay(g.fhPhiZDecay),fhEtaZDecay(g.fhEtaZDecay),
- fhPtPrompt(g.fhPtPrompt),fhPhiPrompt(g.fhPhiPrompt),fhEtaPrompt(g.fhEtaPrompt),
- fhPtUnknown(g.fhPtUnknown),fhPhiUnknown(g.fhPhiUnknown),fhEtaUnknown(g.fhEtaUnknown),
- fhPtHadron(g.fhPtHadron),fhPtEleTrkDet(g.fhPtEleTrkDet),
- //B-tagging
- fhBtagCut1(g.fhBtagCut1),fhBtagCut2(g.fhBtagCut2),fhBtagCut3(g.fhBtagCut3),
- fhBtagQA1(g.fhBtagQA1),fhBtagQA2(g.fhBtagQA2),
- //MC
- fMCEleNtuple(g.fMCEleNtuple),fhPtMCHadron(g.fhPtMCHadron)
-{
- // cpy ctor
-
-}
-
-//_________________________________________________________________________
-AliAnaElectron & AliAnaElectron::operator = (const AliAnaElectron & g)
-{
- // assignment operator
-
- if(&g == this) return *this;
- fCalorimeter = g.fCalorimeter;
- fpOverEmin = g.fpOverEmin;
- fpOverEmax = g.fpOverEmax;
- fResidualCut = g.fResidualCut;
- fDrCut = g.fDrCut;
- fPairDcaCut = g.fPairDcaCut;
- fDecayLenCut = g.fDecayLenCut;
- fImpactCut = g.fImpactCut;
- fAssocPtCut = g.fAssocPtCut;
- fMassCut = g.fMassCut;
- fSdcaCut = g.fSdcaCut;
- fITSCut = g.fITSCut;
- fWriteNtuple = g.fWriteNtuple;
- fEleNtuple = g.fEleNtuple;
- fh1pOverE = g.fh1pOverE;
- fh1dR = g.fh1dR;
- fh2EledEdx = g.fh2EledEdx;
- fh2MatchdEdx = g.fh2MatchdEdx;
- fh2dEtadPhi = g.fh2dEtadPhi;
- fh2dEtadPhiMatched = g.fh2dEtadPhiMatched;
- fh2dEtadPhiUnmatched = g.fh2dEtadPhiUnmatched;
- fh2TrackPVsClusterE = g.fh2TrackPVsClusterE;
- fh2TrackPtVsClusterE = g.fh2TrackPtVsClusterE;
- fh2TrackPhiVsClusterPhi = g.fh2TrackPhiVsClusterPhi;
- fh2TrackEtaVsClusterEta = g.fh2TrackEtaVsClusterEta;
- fh1OpeningAngle = g.fh1OpeningAngle;
- fh1MinvPhoton = g.fh1MinvPhoton;
- 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;
- fhPtConversion = g.fhPtConversion;
- fhPhiConversion = g.fhPhiConversion;
- fhEtaConversion = g.fhEtaConversion;
- fhPtBottom = g.fhPtBottom;
- fhPhiBottom = g.fhPhiBottom;
- fhEtaBottom = g.fhEtaBottom;
- fhPtCharm = g.fhPtCharm;
- fhPhiCharm = g.fhPhiCharm;
- fhEtaCharm = g.fhEtaCharm;
- fhPtCFromB = g.fhPtCFromB;
- fhPhiCFromB = g.fhPhiCFromB;
- fhEtaCFromB = g.fhEtaCFromB;
- fhPtDalitz = g.fhPtDalitz;
- fhPhiDalitz = g.fhPhiDalitz;
- fhEtaDalitz = g.fhEtaDalitz;
- fhPtWDecay = g.fhPtWDecay;
- fhPhiWDecay = g.fhPhiWDecay;
- fhEtaWDecay = g.fhEtaWDecay;
- fhPtZDecay = g.fhPtZDecay;
- fhPhiZDecay = g.fhPhiZDecay;
- fhEtaZDecay = g.fhEtaZDecay;
- fhPtPrompt = g.fhPtPrompt;
- fhPhiPrompt = g.fhPhiPrompt;
- fhEtaPrompt = g.fhEtaPrompt;
- fhPtUnknown = g.fhPtUnknown;
- fhPhiUnknown = g.fhPhiUnknown;
- fhEtaUnknown = g.fhEtaUnknown;
-
- fhPtHadron = g.fhPtHadron;
- fhPtEleTrkDet = g.fhPtEleTrkDet;
-
- //B-tagging
- fhBtagCut1 = g.fhBtagCut1;
- fhBtagCut2 = g.fhBtagCut2;
- fhBtagCut3 = g.fhBtagCut3;
- fhBtagQA1 = g.fhBtagQA1;
- fhBtagQA2 = g.fhBtagQA2;
-
- fMCEleNtuple = g.fMCEleNtuple;
- fhPtMCHadron = g.fhPtMCHadron;
-
- return *this;
-
-}
-
-//____________________________________________________________________________
-AliAnaElectron::~AliAnaElectron()
-{
- //dtor
-
-}
-
-
-//________________________________________________________________________
-TList * AliAnaElectron::GetCreateOutputObjects()
-{
- // Create histograms to be saved in output file and
- // store them in outputContainer
- TList * outputContainer = new TList() ;
- outputContainer->SetName("ElectronHistos") ;
-
- //created ele ntuple for further analysis
- if(fWriteNtuple) {
- fEleNtuple = new TNtuple("EleNtuple","Electron Ntuple","tmctag:cmctag:pt:phi:eta:p:E:deta:dphi:nCells:dEdx:pidProb:impXY:impZ");
- outputContainer->Add(fEleNtuple) ;
- }
-
- Int_t nptbins = GetHistoNPtBins();
- Int_t nphibins = GetHistoNPhiBins();
- Int_t netabins = GetHistoNEtaBins();
- Float_t ptmax = GetHistoPtMax();
- Float_t phimax = GetHistoPhiMax();
- Float_t etamax = GetHistoEtaMax();
- Float_t ptmin = GetHistoPtMin();
- Float_t phimin = GetHistoPhiMin();
- Float_t etamin = GetHistoEtaMin();
-
- fh1pOverE = new TH1F("h1pOverE","EMCAL-TRACK matches p/E",100,0.,10.);
- 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.);
- fh2dEtadPhi = new TH2F("h2dEtadPhi","#Delta#eta vs. #Delta#phi for all track-cluster pairs",200,0.,1.4,300,0.,TMath::Pi());
- fh2dEtadPhiMatched = new TH2F("h2dEtadPhiMatched","#Delta#eta vs. #Delta#phi for matched track-cluster pairs",200,0.,1.4,300,0.,TMath::Pi());
- fh2dEtadPhiUnmatched = new TH2F("h2dEtadPhiUnmatched","#Delta#eta vs. #Delta#phi for unmatched track-cluster pairs",200,0.,1.4,300,0.,TMath::Pi());
-
- fh2TrackPVsClusterE = new TH2F("h2TrackPVsClusterE","h2TrackPVsClusterE",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
- fh2TrackPtVsClusterE = new TH2F("h2TrackPtVsClusterE","h2TrackPtVsClusterE",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
- 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(fh1dR) ;
- outputContainer->Add(fh2EledEdx) ;
- outputContainer->Add(fh2MatchdEdx) ;
- outputContainer->Add(fh2dEtadPhi) ;
- outputContainer->Add(fh2dEtadPhiMatched) ;
- outputContainer->Add(fh2dEtadPhiUnmatched) ;
- outputContainer->Add(fh2TrackPVsClusterE) ;
- outputContainer->Add(fh2TrackPtVsClusterE) ;
- outputContainer->Add(fh2TrackPhiVsClusterPhi) ;
- outputContainer->Add(fh2TrackEtaVsClusterEta) ;
-
- //photonic electron checks
- fh1OpeningAngle = new TH1F("hOpeningAngle","Opening angle between e+e- pairs",100,0.,TMath::Pi());
- fh1MinvPhoton = new TH1F("hMinvPhoton","Invariant mass of e+e- pairs",200,0.,2.);
-
- outputContainer->Add(fh1OpeningAngle);
- outputContainer->Add(fh1MinvPhoton);
-
- fhPtElectron = new TH1F("hPtElectron","Electron pT",nptbins,ptmin,ptmax);
- fhPhiElectron = new TH2F("hPhiElectron","Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
- fhEtaElectron = new TH2F("hEtaElectron","Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
- fhPtNPE = new TH1F("hPtNPE","Non-photonic Electron pT",nptbins,ptmin,ptmax);
- fhPhiNPE = new TH2F("hPhiNPE","Non-photonic Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
- fhEtaNPE = new TH2F("hEtaNPE","Non-photonic Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
- fhPtPE = new TH1F("hPtPE","Photonic Electron pT",nptbins,ptmin,ptmax);
- fhPhiPE = new TH2F("hPhiPE","Photonic Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
- fhEtaPE = new TH2F("hEtaPE","Photonic Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
-
- outputContainer->Add(fhPtElectron) ;
- outputContainer->Add(fhPhiElectron) ;
- outputContainer->Add(fhEtaElectron) ;
- outputContainer->Add(fhPtNPE) ;
- outputContainer->Add(fhPhiNPE) ;
- outputContainer->Add(fhEtaNPE) ;
- outputContainer->Add(fhPtPE) ;
- outputContainer->Add(fhPhiPE) ;
- outputContainer->Add(fhEtaPE) ;
-
- fhPtHadron = new TH1F("hPtHadron","Charged hadrons w/in EMCAL acceptance",nptbins,ptmin,ptmax);
- fhPtEleTrkDet = new TH1F("hPtEleTrkDet","Electrons identified by tracking detectors w/in EMCAL acceptance",nptbins,ptmin,ptmax);
-
- outputContainer->Add(fhPtHadron);
- outputContainer->Add(fhPtEleTrkDet);
-
- //B-tagging
- fhBtagCut1 = new TH2F("hbtag_cut1","B-tag result cut1", 10,0,10 ,nptbins,ptmin,ptmax);
- fhBtagCut2 = new TH2F("hbtag_cut2","B-tag result cut2", 10,0,10 ,nptbins,ptmin,ptmax);
- fhBtagCut3 = new TH2F("hbtag_cut3","B-tag result cut3", 10,0,10 ,nptbins,ptmin,ptmax);
- fhBtagQA1 = new TH2F("hbtag_qa1" ,"B-tag QA: pairDCA vs length", 100,0,0.2 ,100,0,1.0);
- fhBtagQA2 = new TH2F("hbtag_qa2" ,"B-tag QA: signDCA vs mass" , 200,-0.5,0.5 ,100,0,10);
-
- outputContainer->Add(fhBtagCut1) ;
- outputContainer->Add(fhBtagCut2) ;
- outputContainer->Add(fhBtagCut3) ;
- outputContainer->Add(fhBtagQA1) ;
- outputContainer->Add(fhBtagQA2) ;
-
- if(IsDataMC()){
-
- fhPtConversion = new TH1F("hPtConversion","Conversion electron pT",nptbins,ptmin,ptmax);
- 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);
- fhPtBottom = new TH1F("hPtBottom","Bottom electron pT",nptbins,ptmin,ptmax);
- fhPhiBottom = new TH2F("hPhiBottom","Bottom Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
- fhEtaBottom = new TH2F("hEtaBottom","Bottom Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
- fhPtCharm = new TH1F("hPtCharm","Charm electron pT",nptbins,ptmin,ptmax);
- fhPhiCharm = new TH2F("hPhiCharm","Charm Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
- fhEtaCharm = new TH2F("hEtaCharm","Charm Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
- fhPtCFromB = new TH1F("hPtCFromB","Charm from Bottom electron pT",nptbins,ptmin,ptmax);
- fhPhiCFromB = new TH2F("hPhiCFromB","Charm from Bottom Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
- fhEtaCFromB = new TH2F("hEtaCFromB","Charm from Bottom Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
- fhPtDalitz = new TH1F("hPtDalitz","Dalitz electron pT",nptbins,ptmin,ptmax);
- fhPhiDalitz = new TH2F("hPhiDalitz","Dalitz Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
- fhEtaDalitz = new TH2F("hEtaDalitz","Dalitz Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
- fhPtWDecay = new TH1F("hPtWDecay","W-boson Electron pT",nptbins,ptmin,ptmax);
- fhPhiWDecay = new TH2F("hPhiWDecay","W-boson electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
- fhEtaWDecay = new TH2F("hEtaWDecay","W-boson Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
- fhPtZDecay = new TH1F("hPtZDecay","Z-boson electron pT",nptbins,ptmin,ptmax);
- fhPhiZDecay = new TH2F("hPhiZDecay","Z-boson Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
- fhEtaZDecay = new TH2F("hEtaZDecay","Z-boson Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
- fhPtPrompt = new TH1F("hPtPrompt","Prompt electron pT",nptbins,ptmin,ptmax);
- fhPhiPrompt = new TH2F("hPhiPrompt","Prompt Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
- fhEtaPrompt = new TH2F("hEtaPrompt","Prompt Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
- fhPtUnknown = new TH1F("hPtUnknown","Unknown electron pT",nptbins,ptmin,ptmax);
- fhPhiUnknown = new TH2F("hPhiUnknown","Unknown Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
- fhEtaUnknown = new TH2F("hEtaUnknown","Unknown Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
-
- outputContainer->Add(fhPtConversion);
- outputContainer->Add(fhPhiConversion);
- outputContainer->Add(fhEtaConversion);
- outputContainer->Add(fhPtBottom);
- outputContainer->Add(fhPhiBottom);
- outputContainer->Add(fhEtaBottom);
- outputContainer->Add(fhPtCharm);
- outputContainer->Add(fhPhiCharm);
- outputContainer->Add(fhEtaCharm);
- outputContainer->Add(fhPtCFromB);
- outputContainer->Add(fhPhiCFromB);
- outputContainer->Add(fhEtaCFromB);
- outputContainer->Add(fhPtDalitz);
- outputContainer->Add(fhPhiDalitz);
- outputContainer->Add(fhEtaDalitz);
- outputContainer->Add(fhPtWDecay);
- outputContainer->Add(fhPhiWDecay);
- outputContainer->Add(fhEtaWDecay);
- outputContainer->Add(fhPtZDecay);
- outputContainer->Add(fhPhiZDecay);
- outputContainer->Add(fhEtaZDecay);
- outputContainer->Add(fhPtPrompt);
- outputContainer->Add(fhPhiPrompt);
- outputContainer->Add(fhEtaPrompt);
- outputContainer->Add(fhPtUnknown);
- outputContainer->Add(fhPhiUnknown);
- outputContainer->Add(fhEtaUnknown);
-
- //created ele ntuple for further analysis
- if(fWriteNtuple) {
- fMCEleNtuple = new TNtuple("MCEleNtuple","MC Electron Ntuple","mctag:pt:phi:eta:x:y:z");
- outputContainer->Add(fMCEleNtuple) ;
- }
-
- fhPtMCHadron = new TH1F("hPtMCHadron","MC Charged hadrons w/in EMCAL acceptance",nptbins,ptmin,ptmax);
-
- outputContainer->Add(fhPtMCHadron);
-
- }//Histos with MC
-
- //Save parameters used for analysis
- TString parList ; //this will be list of parameters used for this analysis.
- char onePar[255] ;
-
- 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,"---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 ;
-
- //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);
-
- return outputContainer ;
-
-}
-
-//____________________________________________________________________________
-void AliAnaElectron::Init()
-{
-
- //do some initialization
- if(fCalorimeter == "PHOS") {
- printf("AliAnaElectron::Init() - !!STOP: You want to use PHOS in analysis but this is not (yet) supported!!\n!!Check the configuration file!!\n");
- fCalorimeter = "EMCAL";
- }
- if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn()){
- printf("AliAnaElectron::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!!\n!!Check the configuration file!!\n");
- abort();
- }
-
-}
-
-
-//____________________________________________________________________________
-void AliAnaElectron::InitParameters()
-{
-
- //Initialize the parameters of the analysis.
- SetOutputAODClassName("AliAODPWG4Particle");
- SetOutputAODName("PWG4Particle");
-
- AddToHistogramsName("AnaElectron_");
-
- fCalorimeter = "EMCAL" ;
- fpOverEmin = 0.5;
- fpOverEmax = 1.5;
- fResidualCut = 0.02;
- //B-tagging
- fDrCut = 1.0;
- fPairDcaCut = 0.02;
- fDecayLenCut = 1.0;
- fImpactCut = 0.5;
- fAssocPtCut = 1.0;
- fMassCut = 1.5;
- fSdcaCut = 0.1;
- fITSCut = 4;
-
-}
-
-//__________________________________________________________________
-void AliAnaElectron::MakeAnalysisFillAOD()
-{
- //
- // Do analysis and fill aods with electron candidates
- // These AODs will be used to do subsequent histogram filling
- //
- // Also fill some QA histograms
- //
-
- TObjArray *cl = new TObjArray();
-
- Double_t bfield = 0.;
- if(GetReader()->GetDataType() != AliCaloTrackReader::kMC) bfield = GetReader()->GetBField();
-
- //Select the calorimeter of the electron
- if(fCalorimeter != "EMCAL") {
- printf("This class not yet implemented for PHOS\n");
- abort();
- }
- cl = GetAODEMCAL();
-
- ////////////////////////////////////////////////
- //Start from tracks and get associated clusters
- ////////////////////////////////////////////////
- if(!GetAODCTS() || GetAODCTS()->GetEntriesFast() == 0) return ;
- Int_t ntracks = GetAODCTS()->GetEntriesFast();
- if(GetDebug() > 0)
- printf("AliAnaElectron::MakeAnalysisFillAOD() - In CTS aod entries %d\n", ntracks);
-
- //Unfortunately, AliAODTracks don't have associated EMCAL clusters.
- //we have to redo track-matching, I guess
- Int_t iCluster = -999;
- Int_t bt = 0; //counter for event b-tags
-
- for (Int_t itrk = 0; itrk < ntracks; itrk++) {////////////// track loop
- iCluster = -999; //start with no match
- AliAODTrack * track = (AliAODTrack*) (GetAODCTS()->At(itrk)) ;
- AliAODPid* pid = (AliAODPid*) track->GetDetPid();
-
- Double_t emcpos[3];
- pid->GetEMCALPosition(emcpos);
- Double_t emcmom[3];
- pid->GetEMCALMomentum(emcmom);
-
- TVector3 pos(emcpos[0],emcpos[1],emcpos[2]);
- TVector3 mom(emcmom[0],emcmom[1],emcmom[2]);
- Double_t tphi = pos.Phi();
- 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);
- if(mom.Pt() > GetMinPt() && in) {
-
- Double_t dEdx = pid->GetTPCsignal();
-
- //NOTE: As of 02-Sep-2009, the XYZAtDCA methods of AOD do not
- //work, but it is possible to get the position of a track at
- //closest approach to the vertex from the GetPosition method
- Double_t xyz[3];
- //track->XYZAtDCA(xyz);
- Bool_t isNotDCA = track->GetPosition(xyz);
- if(isNotDCA) printf("##Problem getting impact parameter!\n");
- //printf("\tTRACK POSITION AT DCA: %2.2f,%2.2f,%2.2f\n",xyz[0],xyz[1],xyz[2]);
- Double_t xy = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
- Double_t z = xyz[2];
-
- Int_t ntot = cl->GetEntriesFast();
- Double_t res = 999.;
- Double_t pOverE = -999.;
-
- Int_t pidProb = track->GetMostProbablePID();
- if(pidProb == AliAODTrack::kPion || pidProb == AliAODTrack::kKaon || pidProb == AliAODTrack::kProton) fhPtHadron->Fill(track->Pt());
- if(pidProb == AliAODTrack::kElectron) fhPtEleTrkDet->Fill(track->Pt());
-
- Bool_t isElectron = kFALSE;
- //For tracks in EMCAL acceptance, pair them with all clusters
- //and fill the dEta vs dPhi for these pairs:
- for(Int_t iclus = 0; iclus < ntot; iclus++) {
- AliAODCaloCluster * clus = (AliAODCaloCluster*) (cl->At(iclus));
- if(!clus) continue;
-
- Double_t x[3];
- clus->GetPosition(x);
- TVector3 cluspos(x[0],x[1],x[2]);
- Double_t deta = teta - cluspos.Eta();
- Double_t dphi = tphi - cluspos.Phi();
- if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi();
- if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi();
- fh2dEtadPhi->Fill(deta,dphi);
- fh2TrackPVsClusterE->Fill(clus->E(),track->P());
- fh2TrackPtVsClusterE->Fill(clus->E(),track->Pt());
- fh2TrackPhiVsClusterPhi->Fill(cluspos.Phi(),mom.Phi());
- fh2TrackEtaVsClusterEta->Fill(cluspos.Eta(),mom.Eta());
-
- res = sqrt(dphi*dphi + deta*deta);
- fh1dR->Fill(res);
-
- /////////////////////////////////
- //Perform electron cut analysis//
- /////////////////////////////////
- //Good match
- if(res < fResidualCut) {
- fh2dEtadPhiMatched->Fill(deta,dphi);
- iCluster = iclus;
-
- Int_t tmctag = -1;
- Int_t cmctag = -1;
-
- if(IsDataMC()) {
- //Input from second AOD?
- Int_t input = 0;
- if(GetReader()->GetAODCTSNormalInputEntries() <= itrk) input = 1;
- tmctag = GetMCAnalysisUtils()->CheckOrigin(track->GetLabel(),GetReader(),input);
-
- //Do you want the cluster or the track label?
- input = 0;
- if(GetReader()->GetAODEMCALNormalInputEntries() <= iclus) input = 1;
- cmctag = GetMCAnalysisUtils()->CheckOrigin(clus->GetLabel(0),GetReader(),input);
- }
-
- if(fWriteNtuple) {
- fEleNtuple->Fill(tmctag,cmctag,track->Pt(),track->Phi(),track->Eta(),track->P(),clus->E(),deta,dphi,clus->GetNCells(),dEdx,pidProb,xy,z);
- }
-
- fh2MatchdEdx->Fill(track->P(),dEdx);
-
- Double_t energy = clus->E();
- if(energy > 0) pOverE = tmom/energy;
- fh1pOverE->Fill(pOverE);
-
- Int_t mult = clus->GetNCells();
- if(mult < 2 && GetDebug() > 0) printf("Single digit cluster.\n");
-
- //////////////////////////////
- //Electron cuts happen here!//
- //////////////////////////////
- if(pOverE > fpOverEmin && pOverE < fpOverEmax) isElectron = kTRUE;
- } else {
- fh2dEtadPhiUnmatched->Fill(deta,dphi);
- }
-
- } //calocluster loop
-
- ///////////////////////////
- //Fill AOD with electrons//
- ///////////////////////////
- if(isElectron) {
-
- //B-tagging
- if(GetDebug() > 1) printf("Found Electron - do b-tagging\n");
- Int_t btag = GetBtag(track); bt += btag;
-
- fh2EledEdx->Fill(track->P(),dEdx);
-
- Double_t eMass = 0.511/1000; //mass in GeV
- Double_t eleE = sqrt(track->P()*track->P() + eMass*eMass);
- AliAODPWG4Particle tr = AliAODPWG4Particle(track->Px(),track->Py(),track->Pz(),eleE);
- 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
- tr.SetDetector(fCalorimeter);
- 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
- tr.SetBtag(btag);
-
- //Play with the MC stack if available
- //Check origin of the candidates
- if(IsDataMC()){
-
- //FIXME: Need to re-think this for track-oriented analysis
- //JLK DO WE WANT TRACK TAG OR CLUSTER TAG?
- tr.SetTag(GetMCAnalysisUtils()->CheckOrigin(tr.GetLabel(),GetReader(),tr.GetInputFileIndex()));
-
- if(GetDebug() > 0) printf("AliAnaElectron::MakeAnalysisFillAOD() - Origin of candidate %d\n",tr.GetTag());
- }//Work with stack also
-
- AddAODParticle(tr);
-
- if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillAOD() - Electron selection cuts passed: pT %3.2f, pdg %d\n",tr.Pt(),tr.GetPdg());
- }//electron
- }//pt, fiducial selection
- }//track loop
-
- //FIXME: Should we also check from the calocluster side, just in
- //case?
-
- if(GetDebug() > 1 && bt > 0) printf("AliAnaElectron::MakeAnalysisFillAOD() *** Event Btagged *** \n");
- if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillAOD() End fill AODs \n");
-
-}
-
-//__________________________________________________________________
-void AliAnaElectron::MakeAnalysisFillHistograms()
-{
- //Do analysis and fill histograms
-
- AliStack * stack = 0x0;
- TParticle * primary = 0x0;
- TClonesArray * mcparticles0 = 0x0;
- TClonesArray * mcparticles1 = 0x0;
- AliAODMCParticle * aodprimary = 0x0;
-
- if(IsDataMC()) {
- if(GetReader()->ReadStack()){
- 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
-
- //Loop on stored AOD electrons
- Int_t naod = GetOutputAODBranch()->GetEntriesFast();
- if(GetDebug() > 0) printf("AliAnaElectron::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
-
- for(Int_t iaod = 0; iaod < naod ; iaod++){
- AliAODPWG4Particle* ele = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
- Int_t pdg = ele->GetPdg();
-
- if(GetDebug() > 3)
- printf("AliAnaElectron::MakeAnalysisFillHistograms() - PDG %d, MC TAG %d, Calorimeter %s\n", ele->GetPdg(),ele->GetTag(), (ele->GetDetector()).Data()) ;
-
- if(TMath::Abs(pdg) != AliCaloPID::kElectron) continue;
- if(ele->GetDetector() != fCalorimeter) continue;
-
- if(GetDebug() > 1)
- printf("AliAnaElectron::MakeAnalysisFillHistograms() - ID Electron: pt %f, phi %f, eta %f\n", ele->Pt(),ele->Phi(),ele->Eta()) ;
-
-
- //Filter for photonic electrons based on opening angle and Minv
- //cuts, also fill histograms
- Bool_t photonic = kFALSE;
- photonic = IsItPhotonic(ele);
-
- //Fill electron histograms
- Float_t ptele = ele->Pt();
- Float_t phiele = ele->Phi();
- Float_t etaele = ele->Eta();
-
- fhPtElectron ->Fill(ptele);
- fhPhiElectron ->Fill(ptele,phiele);
- fhEtaElectron ->Fill(ptele,etaele);
-
- if(photonic) {
- fhPtPE->Fill(ptele);
- fhPhiPE->Fill(ptele,phiele);
- fhEtaPE->Fill(ptele,etaele);
- } else {
- fhPtNPE->Fill(ptele);
- fhPhiNPE->Fill(ptele,phiele);
- fhEtaNPE->Fill(ptele,etaele);
- }
-
- if(IsDataMC()){
- Int_t tag = ele->GetTag();
- if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)){
- fhPtConversion ->Fill(ptele);
- fhPhiConversion ->Fill(ptele,phiele);
- fhEtaConversion ->Fill(ptele,etaele);
- }
- else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEFromB))
- {
- fhPtBottom ->Fill(ptele);
- fhPhiBottom ->Fill(ptele,phiele);
- fhEtaBottom ->Fill(ptele,etaele);
- }
- else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEFromC))
- {
- fhPtCharm ->Fill(ptele);
- fhPhiCharm ->Fill(ptele,phiele);
- fhEtaCharm ->Fill(ptele,etaele);
- }
- else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEFromCFromB))
- {
- fhPtCFromB ->Fill(ptele);
- fhPhiCFromB ->Fill(ptele,phiele);
- fhEtaCFromB ->Fill(ptele,etaele);
- }
- else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay) || GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) || GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))
- {
- fhPtDalitz ->Fill(ptele);
- fhPhiDalitz ->Fill(ptele,phiele);
- fhEtaDalitz ->Fill(ptele,etaele);
- }
- else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCWDecay))
- {
- fhPtWDecay ->Fill(ptele);
- fhPhiWDecay ->Fill(ptele,phiele);
- fhEtaWDecay ->Fill(ptele,etaele);
- }
- else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCZDecay))
- {
- fhPtZDecay ->Fill(ptele);
- fhPhiZDecay ->Fill(ptele,phiele);
- fhEtaZDecay ->Fill(ptele,etaele);
- }
- else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron))
- {
- fhPtPrompt ->Fill(ptele);
- fhPhiPrompt ->Fill(ptele,phiele);
- fhEtaPrompt ->Fill(ptele,etaele);
- }
- else{
- fhPtUnknown ->Fill(ptele);
- fhPhiUnknown ->Fill(ptele,phiele);
- fhEtaUnknown ->Fill(ptele,etaele);
- }
- }//Histograms with MC
-
- }// aod loop
-
- ////////////////////////////////////////////////////////
- //Fill histograms of pure MC kinematics from the stack//
- ////////////////////////////////////////////////////////
- if(IsDataMC()) {
- 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);
-
- //fill ntuple
- if(fWriteNtuple) {
- fMCEleNtuple->Fill(mctag,primary->Pt(),primary->Phi(),primary->Eta(),primary->Vx(),primary->Vy(),primary->Vz());
- }
-
- }
-
- } 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;
- }
-
- 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;
-
- if(aodprimary->Pt() < 0.2) continue;
-
- /*
- if(aodprimary->GetDaughter(0) > 0) {
- Int_t dindex = aodprimary->GetDaughter(0);
- //printf("####AODMCparticle daughter index %d and flag value %d\n",dindex,aodprimary->GetFlag());
- Double_t dxv = 0.;
- Double_t dyv = 0.;
- if(ipart < npart0) {
- dxv = ((AliAODMCParticle*)mcparticles0->At(dindex))->Xv();
- dyv = ((AliAODMCParticle*)mcparticles0->At(dindex))->Yv();
- } else {
- dxv = ((AliAODMCParticle*)mcparticles1->At(dindex))->Xv();
- dyv = ((AliAODMCParticle*)mcparticles1->At(dindex))->Yv();
- }
- Double_t dradius = TMath::Sqrt(dxv*dxv+dyv*dyv);
- //printf("\tDaughter production radius = %2.2f\n",dradius);
- //if you convert/decay within the trackable zone, discard
- //from ntuple ?
- // if(dradius < 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);
-
- //fill ntuple
- if(fWriteNtuple) {
- fMCEleNtuple->Fill(mctag,aodprimary->Pt(),aodprimary->Phi(),aodprimary->Eta(),aodprimary->Xv(),aodprimary->Yv(),aodprimary->Zv());
- }
-
- }
- }
- } //pure MC kine histos
-
-}
-
-//__________________________________________________________________
-Int_t AliAnaElectron::GetBtag(AliAODTrack * tr )
-{
- //This method uses the Displaced Vertex between electron-hadron
- //pairs and the primary vertex to determine whether an electron is
- //likely from a B hadron.
-
- Int_t ncls1 = 0;
- for(Int_t l = 0; l < 6; l++) if(TESTBIT(tr->GetITSClusterMap(),l)) ncls1++;
- if (ncls1 < fITSCut) return 0;
-
- Double_t x[3];
- //Note: 02-Sep-2009, Must use GetPosition, not XYZAtDCA
- //Bool_t gotit = tr->XYZAtDCA(x);
- Bool_t isNotDCA = tr->GetPosition(x);
- if(isNotDCA) { printf("##Problem getting impact parameter!\n"); return 0; }
-
- Double_t d1 = TMath::Sqrt(x[0]*x[0] + x[1]*x[1]);
- if (TMath::Abs(d1) > fImpactCut ) return 0;
- if (TMath::Abs(x[2]) > fImpactCut ) return 0;
- //printf("----- impact parameter: x=%f, y=%f, z=%f -------\n",x[0],x[1], x[2]);
-
- Int_t nvtx1 = 0;
- Int_t nvtx2 = 0;
- Int_t nvtx3 = 0;
-
- for (Int_t k2 =0; k2 < GetAODCTS()->GetEntriesFast() ; k2++) {
- //loop over assoc
- AliAODTrack* track2 = (AliAODTrack*) (GetAODCTS()->At(k2));
- Int_t id1 = tr->GetID();
- Int_t id2 = track2->GetID();
- if(id1 == id2) continue;
-
- Int_t ncls2 = 0;
- for(Int_t l = 0; l < 6; l++) if(TESTBIT(track2->GetITSClusterMap(),l)) ncls2++;
- if (ncls2 < fITSCut) continue;
-
- if(track2->Pt() < fAssocPtCut) continue;
-
- Double_t dphi = tr->Phi() - track2->Phi();
- if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi();
- if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi();
- Double_t deta = tr->Eta() - track2->Eta();
- Double_t dr = sqrt(deta*deta + dphi*dphi);
-
- if(dr > fDrCut) continue;
-
- Double_t sDca1 = ComputeSignDca(tr, track2, 1.0);
- if (sDca1 > fSdcaCut) nvtx1++;
- Double_t sDca2 = ComputeSignDca(tr, track2, 1.5);
- if (sDca2 > fSdcaCut) nvtx2++;
- Double_t sDca3 = ComputeSignDca(tr, track2, 1.8);
- if (sDca3 > fSdcaCut) nvtx3++;
-
- } //loop over hadrons
-
- if(GetDebug() > 0) {
- if (nvtx1>0) printf("result1 of btagging: %d \n",nvtx1);
- if (nvtx2>0) printf("result2 of btagging: %d \n",nvtx2);
- if (nvtx3>0) printf("result3 of btagging: %d \n",nvtx3);
- }
-
- //fill QA histograms
- fhBtagCut1->Fill(nvtx1,tr->Pt());
- fhBtagCut2->Fill(nvtx2,tr->Pt());
- fhBtagCut3->Fill(nvtx3,tr->Pt());
-
- return nvtx2;
-}
-
-//__________________________________________________________________
-Double_t AliAnaElectron::ComputeSignDca(AliAODTrack *tr, AliAODTrack *tr2 , float masscut)
-{
- //Compute the signed dca between two tracks
- //and return the result
-
- Double_t signDca=-999.;
- if(GetDebug() > 2 ) printf(">>ComputeSdca:: track1 %d, track2 %d, masscut %f \n", tr->GetLabel(), tr2->GetLabel(), masscut);
-
- //=====Now calculate DCA between both tracks=======
- 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
- //FIXME: Add a check for whether file 2 is PYTHIA or HIJING
- //If PYTHIA, then set the vertex from file 2, if not, use the
- //vertex from file 1
- if(GetReader()->GetSecondInputAODTree()) GetReader()->GetSecondInputAODVertex(vertex);
- }
-
- TVector3 primV(vertex[0],vertex[1],vertex[2]) ;
-
- if(GetDebug() > 5) printf(">>ComputeSdca:: primary vertex = %2.2f,%2.2f,%2.2f \n",vertex[0],vertex[1],vertex[2]) ;
-
- AliExternalTrackParam *param1 = new AliExternalTrackParam(tr);
- AliExternalTrackParam *param2 = new AliExternalTrackParam(tr2);
-
- Double_t xplane1 = 0.; Double_t xplane2 = 0.;
- Double_t pairdca = param1->GetDCA(param2,bfield,xplane1,xplane2);
-
- Int_t id1 = 0, id2 = 0;
- AliESDv0 bvertex(*param1,id1,*param2,id2);
- Double_t vx,vy,vz;
- bvertex.GetXYZ(vx,vy,vz);
-
- Double_t emom[3];
- Double_t hmom[3];
- param1->PxPyPz(emom);
- param2->PxPyPz(hmom);
- TVector3 emomAtB(emom[0],emom[1],emom[2]);
- TVector3 hmomAtB(hmom[0],hmom[1],hmom[2]);
- TVector3 secvtxpt(vx,vy,vz);
- TVector3 decayvector(0,0,0);
- decayvector = secvtxpt - primV; //decay vector from PrimVtx
- Double_t decaylength = decayvector.Mag();
-
- if(GetDebug() > 0) {
- printf(">>ComputeSdca:: mom1=%f, mom2=%f \n", emomAtB.Perp(), hmomAtB.Perp() );
- printf(">>ComputeSdca:: pairDCA=%f, length=%f \n", pairdca,decaylength );
- }
-
- if (masscut<1.1) fhBtagQA1->Fill(pairdca,decaylength);
-
- if (emomAtB.Mag()>0 && pairdca < fPairDcaCut && decaylength < fDecayLenCut ) {
- TVector3 sumMom = emomAtB+hmomAtB;
- Double_t ener1 = sqrt(pow(emomAtB.Mag(),2) + massE*massE);
- Double_t ener2 = sqrt(pow(hmomAtB.Mag(),2) + massK*massK);
- Double_t ener3 = sqrt(pow(hmomAtB.Mag(),2) + massE*massE);
- Double_t mass = sqrt(pow((ener1+ener2),2) - pow(sumMom.Mag(),2));
- Double_t massPhot = sqrt(pow((ener1+ener3),2) - pow(sumMom.Mag(),2));
- Double_t sDca = decayvector.Dot(emomAtB)/emomAtB.Mag();
-
- if (masscut<1.1) fhBtagQA2->Fill(sDca, mass);
-
- if (mass > masscut && massPhot > 0.1) signDca = sDca;
-
- if(GetDebug() > 0) printf("\t>>ComputeSdca:: mass=%f \n", mass);
- if(GetDebug() > 0) printf("\t>>ComputeSdca:: sec vtx-signdca :%f\n",signDca);
- }
-
- //clean up
- delete param1;
- delete param2;
-
- return signDca;
-}
-
-//__________________________________________________________________
-Bool_t AliAnaElectron::IsItPhotonic(const AliAODPWG4Particle* part)
-{
- //This method checks the opening angle and invariant mass of
- //electron pairs to see if they are likely to be photonic electrons
-
- Bool_t itIS = kFALSE;
-
- Double_t massE = 0.000511;
- Double_t bfield = 5.; //kG
- if(GetReader()->GetDataType() != AliCaloTrackReader::kMC) bfield = GetReader()->GetBField();
-
- Int_t pdg1 = part->GetPdg();
- Int_t trackId = part->GetTrackLabel(0);
- AliAODTrack* track = (AliAODTrack*)GetAODCTS()->At(trackId);
- if(!track) {
- if(GetDebug() > 0) printf("AliAnaElectron::IsItPhotonic - can't get the AOD Track from the particle! Skipping the photonic check");
- return kFALSE; //Don't proceed because we can't get the track
- }
-
- AliExternalTrackParam *param1 = new AliExternalTrackParam(track);
-
- //Loop on stored AOD electrons and compute the angle differences and Minv
- for (Int_t k2 =0; k2 < GetOutputAODBranch()->GetEntriesFast() ; k2++) {
- AliAODPWG4Particle* part2 = (AliAODPWG4Particle*) GetOutputAODBranch()->At(k2);
- Int_t track2Id = part2->GetTrackLabel(0);
- if(trackId == track2Id) continue;
- Int_t pdg2 = part2->GetPdg();
- if(TMath::Abs(pdg2) != AliCaloPID::kElectron) continue;
- if(part2->GetDetector() != fCalorimeter) continue;
-
- //JLK: Check opp. sign pairs only
- if(pdg1*pdg2 > 0) continue; //skip same-sign pairs
-
- //propagate to common vertex and check opening angle
- AliAODTrack* track2 = (AliAODTrack*)GetAODCTS()->At(track2Id);
- if(!track2) {
- if(GetDebug() >0) printf("AliAnaElectron::IsItPhotonic - problem getting the partner track. Continuing on to the next one");
- continue;
- }
- AliExternalTrackParam *param2 = new AliExternalTrackParam(track2);
- Int_t id1 = 0, id2 = 0;
- AliESDv0 photonVtx(*param1,id1,*param2,id2);
- Double_t vx,vy,vz;
- photonVtx.GetXYZ(vx,vy,vz);
-
- Double_t p1mom[3];
- Double_t p2mom[3];
- param1->PxPyPz(p1mom);
- param2->PxPyPz(p2mom);
-
- TVector3 p1momAtB(p1mom[0],p1mom[1],p1mom[2]);
- TVector3 p2momAtB(p2mom[0],p2mom[1],p2mom[2]);
- TVector3 sumMom = p1momAtB+p2momAtB;
-
- Double_t ener1 = sqrt(pow(p1momAtB.Mag(),2) + massE*massE);
- Double_t ener2 = sqrt(pow(p2momAtB.Mag(),2) + massE*massE);
- Double_t mass = sqrt(pow((ener1+ener2),2) - pow(sumMom.Mag(),2));
-
- Double_t dphi = p1momAtB.DeltaPhi(p2momAtB);
- fh1OpeningAngle->Fill(dphi);
- fh1MinvPhoton->Fill(mass);
-
- if(mass < 0.1) {
- if(GetDebug() > 0) printf("######PROBABLY A PHOTON\n");
- itIS = kTRUE;
- }
-
- //clean up
- delete param2;
-
- }
-
- delete param1;
- return itIS;
-
-}
-
-//__________________________________________________________________
-void AliAnaElectron::Print(const Option_t * opt) const
-{
- //Print some relevant parameters set for the analysis
-
- if(! opt)
- return;
-
- printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
- AliAnaPartCorrBaseClass::Print(" ");
-
- printf("Calorimeter = %s\n", fCalorimeter.Data()) ;
- printf("pOverE range = %f - %f\n",fpOverEmin,fpOverEmax);
- printf("residual cut = %f\n",fResidualCut);
- printf("---Btagging\n");
- printf("max IP-cut (e,h) = %f\n",fImpactCut);
- printf("min ITS-hits = %d\n",fITSCut);
- printf("max dR (e,h) = %f\n",fDrCut);
- printf("max pairDCA = %f\n",fPairDcaCut);
- printf("max decaylength = %f\n",fDecayLenCut);
- printf("min Associated Pt = %f\n",fAssocPtCut);
- printf(" \n") ;
-
-}
-
-//________________________________________________________________________
-void AliAnaElectron::ReadHistograms(TList* outputList)
-{
- // Needed when Terminate is executed in distributed environment
- // Refill analysis histograms of this class with corresponding
- // histograms in output list.
-
- // 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"));
-
- //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++);
-
-}
-
-//__________________________________________________________________
-void AliAnaElectron::Terminate(TList* outputList)
-{
-
- //Do some plots to end
- //Recover histograms from output histograms list, needed for
- //distributed analysis.
- //ReadHistograms(outputList);
-
- printf(" AliAnaElectron::Terminate() *** %s Report: %d outputs\n", GetName(), outputList->GetEntries()) ;
-
-}
-
+/**************************************************************************\r
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
+ * *\r
+ * Author: The ALICE Off-line Project. *\r
+ * Contributors are mentioned in the code where appropriate. *\r
+ * *\r
+ * Permission to use, copy, modify and distribute this software and its *\r
+ * documentation strictly for non-commercial purposes hereby granted *\r
+ * without fee, provided that the above copyright notice appears in all *\r
+ * copies and that both the copyright notice and this permission notice *\r
+ * appear in the supporting documentation. The authors make no claims *\r
+ * about the suitability of this software for any purpose. It is *\r
+ * provided "as is" without express or implied warranty. *\r
+ **************************************************************************/\r
+/* $Id: $ */\r
+\r
+//_________________________________________________________________________\r
+//\r
+// Class for the electron identification.\r
+// Clusters from EMCAL matched to tracks\r
+// and kept in the AOD. Few histograms produced.\r
+//\r
+// -- Author: J.L. Klay (Cal Poly), M. Heinz (Yale)\r
+//////////////////////////////////////////////////////////////////////////////\r
+ \r
+// --- ROOT system --- \r
+#include <TH2F.h>\r
+#include <TH3F.h>\r
+#include <TParticle.h>\r
+#include <TNtuple.h>\r
+#include <TClonesArray.h>\r
+#include <TObjString.h>\r
+//#include <Riostream.h>\r
+\r
+// --- Analysis system --- \r
+#include "AliAnaElectron.h" \r
+#include "AliCaloTrackReader.h"\r
+#include "AliMCAnalysisUtils.h"\r
+#include "AliVCluster.h"\r
+#include "AliFiducialCut.h"\r
+#include "AliAODTrack.h"\r
+#include "AliAODPid.h"\r
+#include "AliCaloPID.h"\r
+#include "AliAODMCParticle.h"\r
+#include "AliStack.h"\r
+#include "AliExternalTrackParam.h"\r
+#include "AliESDv0.h"\r
+#include "AliESDtrack.h"\r
+#include "AliAODJet.h"\r
+#include "AliAODEvent.h"\r
+#include "AliGenPythiaEventHeader.h"\r
+\r
+ClassImp(AliAnaElectron)\r
+ \r
+//____________________________________________________________________________\r
+AliAnaElectron::AliAnaElectron() \r
+: AliAnaPartCorrBaseClass(),fCalorimeter(""),\r
+ fpOverEmin(0.),fpOverEmax(0.),fResidualCut(0.),fMinClusEne(0.),\r
+ fDrCut(0.),fPairDcaCut(0.),fDecayLenCut(0.),fImpactCut(0.),\r
+ fAssocPtCut(0.),fMassCut(0.),fSdcaCut(0.),fITSCut(0),\r
+ fNTagTrkCut(0),fIPSigCut(0.),fJetEtaCut(0.3),fJetPhiMin(1.8),fJetPhiMax(2.9),\r
+ fWriteNtuple(kFALSE),\r
+ //event QA histos\r
+ fhImpactXY(0),fhRefMult(0),fhRefMult2(0),\r
+ //matching checks\r
+ fh3pOverE(0),fh3EOverp(0),fh3pOverE2(0),fh3EOverp2(0),fh3pOverE3(0),fh3EOverp3(0),\r
+ fh2pOverE(0),fh2EOverp(0),fh2pOverE2(0),fh2EOverp2(0),\r
+ fh1dR(0),fh2EledEdx(0),fh2MatchdEdx(0),fh2dEtadPhi(0),\r
+ fh2dEtadPhiMatched(0),fh2dEtadPhiUnmatched(0),fh2TrackPVsClusterE(0),\r
+ fh2TrackPtVsClusterE(0),fh2TrackPhiVsClusterPhi(0),fh2TrackEtaVsClusterEta(0),\r
+ //Photonic electron checks\r
+ fh1OpeningAngle(0),fh1MinvPhoton(0),\r
+ //Reconstructed electrons\r
+ fhPtElectron(0),fhPhiElectron(0),fhEtaElectron(0),\r
+ fhPtNPE(0),fhPhiNPE(0),fhEtaNPE(0),\r
+ fhPtPE(0),fhPhiPE(0),fhEtaPE(0),\r
+ //for comparisons with tracking detectors\r
+ fhPtHadron(0),fhPtNPEleTPC(0),fhPtNPEleTPCTRD(0),fhPtNPEleTTE(0),\r
+ fhPtNPEleEMCAL(0),\r
+ //DVM B-tagging\r
+ fhDVMBtagCut1(0),fhDVMBtagCut2(0),fhDVMBtagCut3(0),fhDVMBtagQA1(0),fhDVMBtagQA2(0),\r
+ fhDVMBtagQA3(0),fhDVMBtagQA4(0),fhDVMBtagQA5(0),\r
+ //IPSig B-tagging\r
+ fhIPSigBtagQA1(0),fhIPSigBtagQA2(0),fhTagJetPt1x4(0),fhTagJetPt2x3(0),fhTagJetPt3x2(0),\r
+ fhePlusTagJetPt1x4(0),fhePlusTagJetPt2x3(0),fhePlusTagJetPt3x2(0),\r
+ //B-Jet histograms\r
+ fhJetType(0),fhLeadJetType(0),fhBJetXsiFF(0),fhBJetPtFF(0),fhBJetEtaPhi(0),\r
+ fhNonBJetXsiFF(0),fhNonBJetPtFF(0),fhNonBJetEtaPhi(0),\r
+ /////////////////////////////////////////////////////////////\r
+ //Histograms that rely on MC info (not filled for real data)\r
+ fEleNtuple(0),\r
+ //reco electrons from various sources\r
+ fhPhiConversion(0),fhEtaConversion(0),\r
+ //for comparisons with tracking detectors\r
+ fhPtTrack(0),\r
+ fhPtNPEBHadron(0),\r
+ //for computing efficiency of B-jet tags\r
+ fhBJetPt1x4(0),fhBJetPt2x3(0),fhBJetPt3x2(0),\r
+ fhFakeJetPt1x4(0),fhFakeJetPt2x3(0),fhFakeJetPt3x2(0),fhDVMJet(0),\r
+ //MC rate histograms/ntuple\r
+ fMCEleNtuple(0),fhMCBJetElePt(0),fhMCBHadronElePt(0),fhPtMCHadron(0),fhPtMCElectron(0),\r
+ fhMCXYConversion(0),fhMCRadPtConversion(0)\r
+{\r
+ //default ctor\r
+ \r
+ //Initialize parameters\r
+ InitParameters();\r
+\r
+}\r
+/*\r
+//____________________________________________________________________________\r
+AliAnaElectron::AliAnaElectron(const AliAnaElectron & g) \r
+ : AliAnaPartCorrBaseClass(g),fCalorimeter(g.fCalorimeter),\r
+ fpOverEmin(g.fpOverEmin),fpOverEmax(g.fpOverEmax),\r
+ fResidualCut(g.fResidualCut),fMinClusEne(g.fMinClusEne),\r
+ fDrCut(g.fDrCut),fPairDcaCut(g.fPairDcaCut),fDecayLenCut(g.fDecayLenCut),fImpactCut(g.fImpactCut),\r
+ fAssocPtCut(g.fAssocPtCut),fMassCut(g.fMassCut),fSdcaCut(g.fSdcaCut),fITSCut(g.fITSCut),\r
+ fNTagTrkCut(g.fNTagTrkCut),fIPSigCut(g.fIPSigCut),\r
+ fJetEtaCut(g.fJetEtaCut),fJetPhiMin(g.fJetPhiMin),fJetPhiMax(g.fJetPhiMax),\r
+ fWriteNtuple(g.fWriteNtuple),\r
+ //event QA histos\r
+ fhImpactXY(g.fhImpactXY),fhRefMult(g.fhRefMult),fhRefMult2(g.fhRefMult2),\r
+ //matching checks\r
+ fh3pOverE(g.fh3pOverE),fh3EOverp(g.fh3EOverp),\r
+ fh3pOverE2(g.fh3pOverE2),fh3EOverp2(g.fh3EOverp2),\r
+ fh3pOverE3(g.fh3pOverE3),fh3EOverp3(g.fh3EOverp3),\r
+ fh2pOverE(g.fh2pOverE),fh2EOverp(g.fh2EOverp),\r
+ fh2pOverE2(g.fh2pOverE2),fh2EOverp2(g.fh2EOverp2),\r
+ fh1dR(g.fh1dR),fh2EledEdx(g.fh2EledEdx),\r
+ fh2MatchdEdx(g.fh2MatchdEdx),fh2dEtadPhi(g.fh2dEtadPhi),\r
+ fh2dEtadPhiMatched(g.fh2dEtadPhiMatched),fh2dEtadPhiUnmatched(g.fh2dEtadPhiUnmatched),\r
+ fh2TrackPVsClusterE(g.fh2TrackPVsClusterE),fh2TrackPtVsClusterE(g.fh2TrackPtVsClusterE),\r
+ fh2TrackPhiVsClusterPhi(g.fh2TrackPhiVsClusterPhi),fh2TrackEtaVsClusterEta(g.fh2TrackEtaVsClusterEta),\r
+ //Photonic electron checks\r
+ fh1OpeningAngle(g.fh1OpeningAngle),fh1MinvPhoton(g.fh1MinvPhoton),\r
+ //Reconstructed electrons\r
+ fhPtElectron(g.fhPtElectron),fhPhiElectron(g.fhPhiElectron),fhEtaElectron(g.fhEtaElectron),\r
+ fhPtNPE(g.fhPtNPE),fhPhiNPE(g.fhPhiNPE),fhEtaNPE(g.fhEtaNPE),\r
+ fhPtPE(g.fhPtPE),fhPhiPE(g.fhPhiPE),fhEtaPE(g.fhEtaPE),\r
+ //for comparisons with tracking detectors\r
+ fhPtHadron(g.fhPtHadron),fhPtNPEleTPC(g.fhPtNPEleTPC),\r
+ fhPtNPEleTPCTRD(g.fhPtNPEleTPCTRD),fhPtNPEleTTE(g.fhPtNPEleTTE),\r
+ fhPtNPEleEMCAL(g.fhPtNPEleEMCAL),\r
+ //DVM B-tagging\r
+ fhDVMBtagCut1(g.fhDVMBtagCut1),fhDVMBtagCut2(g.fhDVMBtagCut2),fhDVMBtagCut3(g.fhDVMBtagCut3),\r
+ fhDVMBtagQA1(g.fhDVMBtagQA1),fhDVMBtagQA2(g.fhDVMBtagQA2),\r
+ fhDVMBtagQA3(g.fhDVMBtagQA3),fhDVMBtagQA4(g.fhDVMBtagQA4),fhDVMBtagQA5(g.fhDVMBtagQA5),\r
+ //IPSig B-tagging\r
+ fhIPSigBtagQA1(g.fhIPSigBtagQA1),fhIPSigBtagQA2(g.fhIPSigBtagQA2),\r
+ fhTagJetPt1x4(g.fhTagJetPt1x4),fhTagJetPt2x3(g.fhTagJetPt2x3),fhTagJetPt3x2(g.fhTagJetPt3x2),\r
+ fhePlusTagJetPt1x4(g.fhePlusTagJetPt1x4),fhePlusTagJetPt2x3(g.fhePlusTagJetPt2x3),\r
+ fhePlusTagJetPt3x2(g.fhePlusTagJetPt3x2),\r
+ //B-Jet histograms\r
+ fhJetType(g.fhJetType),fhLeadJetType(g.fhLeadJetType),fhBJetXsiFF(g.fhBJetXsiFF),\r
+ fhBJetPtFF(g.fhBJetPtFF),fhBJetEtaPhi(g.fhBJetEtaPhi),fhNonBJetXsiFF(g.fhNonBJetXsiFF),\r
+ fhNonBJetPtFF(g.fhNonBJetPtFF),fhNonBJetEtaPhi(g.fhNonBJetEtaPhi),\r
+ /////////////////////////////////////////////////////////////\r
+ //Histograms that rely on MC info (not filled for real data)\r
+ fEleNtuple(g.fEleNtuple),\r
+ //reco electrons from various sources\r
+ fhPhiConversion(g.fhPhiConversion),fhEtaConversion(g.fhEtaConversion),\r
+ //for comparisons with tracking detectors\r
+ fhPtTrack(g.fhPtTrack),\r
+ fhPtNPEBHadron(g.fhPtNPEBHadron),\r
+ //for computing efficiency of B-jet tags\r
+ fhBJetPt1x4(g.fhBJetPt1x4),fhBJetPt2x3(g.fhBJetPt2x3),\r
+ fhBJetPt3x2(g.fhBJetPt3x2),\r
+ fhFakeJetPt1x4(g.fhFakeJetPt1x4),fhFakeJetPt2x3(g.fhBJetPt2x3),\r
+ fhFakeJetPt3x2(g.fhFakeJetPt3x2),fhDVMJet(g.fhDVMJet),\r
+ //MC rate histograms/ntuple\r
+ fMCEleNtuple(g.fMCEleNtuple),fhMCBJetElePt(g.fhMCBJetElePt),\r
+ fhMCBHadronElePt(g.fhMCBHadronElePt),\r
+ fhPtMCHadron(g.fhPtMCHadron),fhPtMCElectron(g.fhPtMCElectron),\r
+ fhMCXYConversion(g.fhMCXYConversion),fhMCRadPtConversion(g.fhMCRadPtConversion)\r
+{\r
+ // cpy ctor\r
+ \r
+}\r
+\r
+//_________________________________________________________________________\r
+AliAnaElectron & AliAnaElectron::operator = (const AliAnaElectron & g)\r
+{\r
+ // assignment operator\r
+ \r
+ if(&g == this) return *this;\r
+ fCalorimeter = g.fCalorimeter;\r
+ fpOverEmin = g.fpOverEmin;\r
+ fpOverEmax = g.fpOverEmax;\r
+ fResidualCut = g.fResidualCut;\r
+ fMinClusEne = g.fMinClusEne;\r
+ fDrCut = g.fDrCut;\r
+ fPairDcaCut = g.fPairDcaCut;\r
+ fDecayLenCut = g.fDecayLenCut;\r
+ fImpactCut = g.fImpactCut;\r
+ fAssocPtCut = g.fAssocPtCut;\r
+ fMassCut = g.fMassCut;\r
+ fSdcaCut = g.fSdcaCut;\r
+ fITSCut = g.fITSCut;\r
+ fNTagTrkCut = g.fNTagTrkCut;\r
+ fIPSigCut = g.fIPSigCut;\r
+ fJetEtaCut = g.fJetEtaCut;\r
+ fJetPhiMin = g.fJetPhiMin;\r
+ fJetPhiMax = g.fJetPhiMax;\r
+ fWriteNtuple = g.fWriteNtuple;\r
+ //event QA histos\r
+ fhImpactXY = g.fhImpactXY;\r
+ fhRefMult = g.fhRefMult;\r
+ fhRefMult2 = g.fhRefMult2;\r
+ //matching checks\r
+ fh3pOverE = g.fh3pOverE;\r
+ fh3EOverp = g.fh3EOverp;\r
+ fh3pOverE2 = g.fh3pOverE2;\r
+ fh3EOverp2 = g.fh3EOverp2;\r
+ fh3pOverE3 = g.fh3pOverE3;\r
+ fh3EOverp3 = g.fh3EOverp3;\r
+ fh2pOverE = g.fh2pOverE;\r
+ fh2EOverp = g.fh2EOverp;\r
+ fh2pOverE2 = g.fh2pOverE2;\r
+ fh2EOverp2 = g.fh2EOverp2;\r
+ fh1dR = g.fh1dR;\r
+ fh2EledEdx = g.fh2EledEdx;\r
+ fh2MatchdEdx = g.fh2MatchdEdx;\r
+ fh2dEtadPhi = g.fh2dEtadPhi;\r
+ fh2dEtadPhiMatched = g.fh2dEtadPhiMatched;\r
+ fh2dEtadPhiUnmatched = g.fh2dEtadPhiUnmatched;\r
+ fh2TrackPVsClusterE = g.fh2TrackPVsClusterE;\r
+ fh2TrackPtVsClusterE = g.fh2TrackPtVsClusterE;\r
+ fh2TrackPhiVsClusterPhi = g.fh2TrackPhiVsClusterPhi;\r
+ fh2TrackEtaVsClusterEta = g.fh2TrackEtaVsClusterEta;\r
+ //Photonic electron checks\r
+ fh1OpeningAngle = g.fh1OpeningAngle;\r
+ fh1MinvPhoton = g.fh1MinvPhoton;\r
+ //Reconstructed electrons\r
+ fhPtElectron = g.fhPtElectron; \r
+ fhPhiElectron = g.fhPhiElectron; \r
+ fhEtaElectron = g.fhEtaElectron; \r
+ fhPtNPE = g.fhPtNPE;\r
+ fhPhiNPE = g.fhPhiNPE;\r
+ fhEtaNPE = g.fhEtaNPE; \r
+ fhPtPE = g.fhPtPE;\r
+ fhPhiPE = g.fhPhiPE;\r
+ fhEtaPE = g.fhEtaPE; \r
+ //for comparisons with tracking detectors\r
+ fhPtHadron = g.fhPtHadron; fhPtNPEleTPC = g.fhPtNPEleTPC; \r
+ fhPtNPEleTPCTRD = g.fhPtNPEleTPCTRD; fhPtNPEleTTE = g.fhPtNPEleTTE; \r
+ fhPtNPEleEMCAL = g.fhPtNPEleEMCAL; \r
+ //DVM B-tagging\r
+ fhDVMBtagCut1 = g.fhDVMBtagCut1;\r
+ fhDVMBtagCut2 = g.fhDVMBtagCut2; \r
+ fhDVMBtagCut3 = g.fhDVMBtagCut3; \r
+ fhDVMBtagQA1 = g.fhDVMBtagQA1; \r
+ fhDVMBtagQA2 = g.fhDVMBtagQA2; \r
+ fhDVMBtagQA3 = g.fhDVMBtagQA3; \r
+ fhDVMBtagQA4 = g.fhDVMBtagQA4; \r
+ fhDVMBtagQA5 = g.fhDVMBtagQA5; \r
+ //IPSig B-tagging\r
+ fhIPSigBtagQA1 = g.fhIPSigBtagQA1; \r
+ fhIPSigBtagQA2 = g.fhIPSigBtagQA2; \r
+ fhTagJetPt1x4 = g.fhTagJetPt1x4; \r
+ fhTagJetPt2x3 = g.fhTagJetPt2x3; \r
+ fhTagJetPt3x2 = g.fhTagJetPt3x2; \r
+ fhePlusTagJetPt1x4 = g.fhePlusTagJetPt1x4; \r
+ fhePlusTagJetPt2x3 = g.fhePlusTagJetPt2x3; \r
+ fhePlusTagJetPt3x2 = g.fhePlusTagJetPt3x2; \r
+ //B-Jet histograms\r
+ fhJetType = g.fhJetType; \r
+ fhLeadJetType = g.fhLeadJetType; \r
+ fhBJetXsiFF = g.fhBJetXsiFF; \r
+ fhBJetPtFF = g.fhBJetPtFF; \r
+ fhBJetEtaPhi = g.fhBJetEtaPhi; \r
+ fhNonBJetXsiFF = g.fhNonBJetXsiFF; \r
+ fhNonBJetPtFF = g.fhNonBJetPtFF; \r
+ fhNonBJetEtaPhi = g.fhNonBJetEtaPhi; \r
+ /////////////////////////////////////////////////////////////\r
+ //Histograms that rely on MC info (not filled for real data)\r
+ fEleNtuple = g.fEleNtuple; \r
+ //reco electrons from various sources\r
+ fhPhiConversion = g.fhPhiConversion; \r
+ fhEtaConversion = g.fhEtaConversion;\r
+ //for comparisons with tracking detectors\r
+ fhPtTrack = g.fhPtTrack;\r
+ fhPtNPEBHadron = g.fhPtNPEBHadron;\r
+ //for computing efficiency of B-jet tags\r
+ fhBJetPt1x4 = g.fhBJetPt1x4; fhBJetPt2x3 = g.fhBJetPt2x3; \r
+ fhBJetPt3x2 = g.fhBJetPt3x2;\r
+ fhFakeJetPt1x4 = g.fhFakeJetPt1x4; fhFakeJetPt2x3 = g.fhFakeJetPt2x3; \r
+ fhFakeJetPt3x2 = g.fhFakeJetPt3x2; fhDVMJet = g.fhDVMJet;\r
+ //MC rate histograms/ntuple\r
+ fMCEleNtuple = g.fMCEleNtuple; fhMCBJetElePt = g.fhMCBJetElePt; \r
+ fhMCBHadronElePt = g.fhMCBHadronElePt;\r
+ fhPtMCHadron = g.fhPtMCHadron; fhPtMCElectron = g.fhPtMCElectron; \r
+ fhMCXYConversion = g.fhMCXYConversion;\r
+ fhMCRadPtConversion = g.fhMCRadPtConversion;\r
+\r
+ return *this;\r
+ \r
+}\r
+*/\r
+\r
+//____________________________________________________________________________\r
+AliAnaElectron::~AliAnaElectron() \r
+{\r
+ //dtor\r
+\r
+}\r
+\r
+//________________________________________________________________________\r
+TObjString * AliAnaElectron::GetAnalysisCuts()\r
+{\r
+ //Save parameters used for analysis\r
+ TString parList ; //this will be list of parameters used for this analysis.\r
+ const Int_t buffersize = 255; \r
+ char onePar[buffersize] ;\r
+ \r
+ snprintf(onePar,buffersize,"--- AliAnaElectron ---\n") ;\r
+ parList+=onePar ; \r
+ snprintf(onePar,buffersize,"fCalorimeter: %s\n",fCalorimeter.Data()) ;\r
+ parList+=onePar ; \r
+ snprintf(onePar,buffersize,"fpOverEmin: %f\n",fpOverEmin) ;\r
+ parList+=onePar ; \r
+ snprintf(onePar,buffersize,"fpOverEmax: %f\n",fpOverEmax) ;\r
+ parList+=onePar ; \r
+ snprintf(onePar,buffersize,"fResidualCut: %f\n",fResidualCut) ;\r
+ parList+=onePar ; \r
+ snprintf(onePar,buffersize,"fMinClusEne: %f\n",fMinClusEne) ;\r
+ parList+=onePar ;\r
+ snprintf(onePar,buffersize,"---DVM Btagging\n");\r
+ parList+=onePar ;\r
+ snprintf(onePar,buffersize,"max IP-cut (e,h): %f\n",fImpactCut);\r
+ parList+=onePar ;\r
+ snprintf(onePar,buffersize,"min ITS-hits: %d\n",fITSCut);\r
+ parList+=onePar ;\r
+ snprintf(onePar,buffersize,"max dR (e,h): %f\n",fDrCut);\r
+ parList+=onePar ;\r
+ snprintf(onePar,buffersize,"max pairDCA: %f\n",fPairDcaCut);\r
+ parList+=onePar ;\r
+ snprintf(onePar,buffersize,"max decaylength: %f\n",fDecayLenCut);\r
+ parList+=onePar ;\r
+ snprintf(onePar,buffersize,"min Associated Pt: %f\n",fAssocPtCut);\r
+ parList+=onePar ;\r
+ snprintf(onePar,buffersize,"---IPSig Btagging\n");\r
+ parList+=onePar ;\r
+ snprintf(onePar,buffersize,"min tag track: %d\n",fNTagTrkCut);\r
+ parList+=onePar ;\r
+ snprintf(onePar,buffersize,"min IP significance: %f\n",fIPSigCut);\r
+ parList+=onePar ;\r
+ //\r
+ //Get parameters set in base class.\r
+ parList += GetBaseParametersList() ;\r
+ \r
+ //Get parameters set in FiducialCut class (not available yet)\r
+ //parlist += GetFidCut()->GetFidCutParametersList() \r
+ \r
+ return new TObjString(parList) ;\r
+ \r
+}\r
+\r
+//________________________________________________________________________\r
+TList * AliAnaElectron::GetCreateOutputObjects()\r
+{ \r
+ // Create histograms to be saved in output file and \r
+ // store them in outputContainer\r
+ TList * outputContainer = new TList() ; \r
+ outputContainer->SetName("ElectronHistos") ; \r
+\r
+ Int_t nptbins = GetHistoPtBins();\r
+ Int_t nphibins = GetHistoPhiBins();\r
+ Int_t netabins = GetHistoEtaBins();\r
+ Float_t ptmax = GetHistoPtMax();\r
+ Float_t phimax = GetHistoPhiMax();\r
+ Float_t etamax = GetHistoEtaMax();\r
+ Float_t ptmin = GetHistoPtMin();\r
+ Float_t phimin = GetHistoPhiMin();\r
+ Float_t etamin = GetHistoEtaMin(); \r
+\r
+ //event QA\r
+ fhImpactXY = new TH1F("hImpactXY","Impact parameter for all tracks",200,-10,10.);\r
+ fhRefMult = new TH1F("hRefMult" ,"refmult QA: " ,5000,0,5000);\r
+ fhRefMult2 = new TH1F("hRefMult2" ,"refmult2 QA: " ,5000,0,5000);\r
+\r
+ outputContainer->Add(fhImpactXY);\r
+ outputContainer->Add(fhRefMult);\r
+ outputContainer->Add(fhRefMult2);\r
+ \r
+ //matching checks\r
+ fh3pOverE = new TH3F("h3pOverE" ,"EMCAL-TRACK matches p/E",nptbins,ptmin,ptmax,200,0.,10.,30,0,30);\r
+ fh3EOverp = new TH3F("h3EOverp" ,"EMCAL-TRACK matches E/p",nptbins,ptmin,ptmax,200,0.,5. ,30,0,30);\r
+ fh3pOverE2 = new TH3F("h3pOverE_Trk","EMCAL-TRACK matches p/E",nptbins,ptmin,ptmax,200,0.,10.,30,0,30);\r
+ fh3EOverp2 = new TH3F("h3EOverp_Trk","EMCAL-TRACK matches E/p",nptbins,ptmin,ptmax,200,0.,5. ,30,0,30);\r
+ fh3pOverE3 = new TH3F("h3pOverE_Tpc","EMCAL-TRACK matches p/E",nptbins,ptmin,ptmax,200,0.,10.,30,0,30);\r
+ fh3EOverp3 = new TH3F("h3EOverp_Tpc","EMCAL-TRACK matches E/p",nptbins,ptmin,ptmax,200,0.,5. ,30,0,30);\r
+ fh2pOverE = new TH2F("h2pOverE" ,"EMCAL-TRACK matches p/E",nptbins,ptmin,ptmax,200,0.,10.);\r
+ fh2EOverp = new TH2F("h2EOverp" ,"EMCAL-TRACK matches E/p",nptbins,ptmin,ptmax,200,0.,5. );\r
+ fh2pOverE2 = new TH2F("h2pOverE_Trk","EMCAL-TRACK matches p/E",nptbins,ptmin,ptmax,200,0.,10.);\r
+ fh2EOverp2 = new TH2F("h2EOverp_Trk","EMCAL-TRACK matches E/p",nptbins,ptmin,ptmax,200,0.,5. );\r
+\r
+ fh1dR = new TH1F("h1dR","EMCAL-TRACK matches dR",300, 0.,TMath::Pi());\r
+ fh2EledEdx = new TH2F("h2EledEdx","dE/dx vs. p for electrons",200,0.,50.,200,0.,400.);\r
+ fh2MatchdEdx = new TH2F("h2MatchdEdx","dE/dx vs. p for all matches",200,0.,50.,200,0.,400.);\r
+ fh2dEtadPhi = new TH2F("h2dEtadPhi","#Delta#eta vs. #Delta#phi for all track-cluster pairs",200,0.,1.4,300,0.,TMath::Pi());\r
+ fh2dEtadPhiMatched = new TH2F("h2dEtadPhiMatched","#Delta#eta vs. #Delta#phi for matched track-cluster pairs",200,0.,1.4,300,0.,TMath::Pi());\r
+ fh2dEtadPhiUnmatched = new TH2F("h2dEtadPhiUnmatched","#Delta#eta vs. #Delta#phi for unmatched track-cluster pairs",200,0.,1.4,300,0.,TMath::Pi());\r
+\r
+ fh2TrackPVsClusterE = new TH2F("h2TrackPVsClusterE","h2TrackPVsClusterE",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);\r
+ fh2TrackPtVsClusterE = new TH2F("h2TrackPtVsClusterE","h2TrackPtVsClusterE",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);\r
+ fh2TrackPhiVsClusterPhi = new TH2F("h2TrackPhiVsClusterPhi","h2TrackPhiVsClusterPhi",nphibins,phimin,phimax,nphibins,phimin,phimax);\r
+ fh2TrackEtaVsClusterEta = new TH2F("h2TrackEtaVsClusterEta","h2TrackEtaVsClusterEta",netabins,etamin,etamax,netabins,etamin,etamax);\r
+\r
+ outputContainer->Add(fh3pOverE) ;\r
+ outputContainer->Add(fh3EOverp) ;\r
+ outputContainer->Add(fh3pOverE2) ;\r
+ outputContainer->Add(fh3EOverp2) ;\r
+ outputContainer->Add(fh3pOverE3) ;\r
+ outputContainer->Add(fh3EOverp3) ;\r
+ outputContainer->Add(fh2pOverE) ;\r
+ outputContainer->Add(fh2EOverp) ;\r
+ outputContainer->Add(fh2pOverE2) ;\r
+ outputContainer->Add(fh2EOverp2) ;\r
+ outputContainer->Add(fh1dR) ; \r
+ outputContainer->Add(fh2EledEdx) ;\r
+ outputContainer->Add(fh2MatchdEdx) ;\r
+ outputContainer->Add(fh2dEtadPhi) ;\r
+ outputContainer->Add(fh2dEtadPhiMatched) ;\r
+ outputContainer->Add(fh2dEtadPhiUnmatched) ;\r
+ outputContainer->Add(fh2TrackPVsClusterE) ;\r
+ outputContainer->Add(fh2TrackPtVsClusterE) ;\r
+ outputContainer->Add(fh2TrackPhiVsClusterPhi) ;\r
+ outputContainer->Add(fh2TrackEtaVsClusterEta) ;\r
+ \r
+ //photonic electron checks\r
+ fh1OpeningAngle = new TH1F("hOpeningAngle","Opening angle between e+e- pairs",100,0.,TMath::Pi());\r
+ fh1MinvPhoton = new TH1F("hMinvPhoton","Invariant mass of e+e- pairs",200,0.,2.);\r
+\r
+ outputContainer->Add(fh1OpeningAngle);\r
+ outputContainer->Add(fh1MinvPhoton);\r
+\r
+ //Reconstructed electrons\r
+ fhPtElectron = new TH1F("hPtElectron","Electron pT",nptbins,ptmin,ptmax);\r
+ fhPhiElectron = new TH2F("hPhiElectron","Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax);\r
+ fhEtaElectron = new TH2F("hEtaElectron","Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);\r
+ fhPtNPE = new TH1F("hPtNPE","Non-photonic Electron pT",nptbins,ptmin,ptmax);\r
+ fhPhiNPE = new TH2F("hPhiNPE","Non-photonic Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax);\r
+ fhEtaNPE = new TH2F("hEtaNPE","Non-photonic Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);\r
+ fhPtPE = new TH1F("hPtPE","Photonic Electron pT",nptbins,ptmin,ptmax);\r
+ fhPhiPE = new TH2F("hPhiPE","Photonic Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax);\r
+ fhEtaPE = new TH2F("hEtaPE","Photonic Electron eta vs. eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);\r
+\r
+ outputContainer->Add(fhPtElectron) ; \r
+ outputContainer->Add(fhPhiElectron) ; \r
+ outputContainer->Add(fhEtaElectron) ;\r
+ outputContainer->Add(fhPtNPE) ; \r
+ outputContainer->Add(fhPhiNPE) ; \r
+ outputContainer->Add(fhEtaNPE) ;\r
+ outputContainer->Add(fhPtPE) ; \r
+ outputContainer->Add(fhPhiPE) ; \r
+ outputContainer->Add(fhEtaPE) ;\r
+\r
+ //These histograms are mixed REAL/MC:\r
+ //Bins along y-axis are: \r
+ //0 - unfiltered (filled for both real and MC data) \r
+ //1 - bottom, 2 - charm, 3 - charm from bottom (MC only)\r
+ //4 - conversion, 5 - Dalitz, 6 - W and Z, 7 - junk/unknown (MC only)\r
+ //8 - misidentified (MC only)\r
+\r
+ //histograms for comparison to tracking detectors\r
+ fhPtHadron = new TH2F("hPtHadron","Charged hadrons w/in EMCAL acceptance",nptbins,ptmin,ptmax,10,0,10);\r
+ fhPtNPEleTPC = new TH2F("hPtNPEleTPC","Non-phot. Electrons identified by TPC w/in EMCAL acceptance",nptbins,ptmin,ptmax,10,0,10);\r
+ fhPtNPEleTPCTRD = new TH2F("hPtNPEleTPCTRD","Non-phot. Electrons identified by TPC+TRD w/in EMCAL acceptance",nptbins,ptmin,ptmax,10,0,10);\r
+ fhPtNPEleTTE = new TH2F("hPtNPEleTTE","Non-phot. Electrons identified by TPC+TRD+EMCAL w/in EMCAL acceptance",nptbins,ptmin,ptmax,10,0,10); \r
+ fhPtNPEleEMCAL = new TH2F("hPtNPEleEMCAL","Non-phot. Electrons identified by EMCAL w/in EMCAL acceptance",nptbins,ptmin,ptmax,10,0,10);\r
+ \r
+ outputContainer->Add(fhPtHadron);\r
+ outputContainer->Add(fhPtNPEleTPC);\r
+ outputContainer->Add(fhPtNPEleTPCTRD);\r
+ outputContainer->Add(fhPtNPEleTTE);\r
+ outputContainer->Add(fhPtNPEleEMCAL);\r
+\r
+ //B-tagging\r
+ fhDVMBtagCut1 = new TH2F("hdvmbtag_cut1","DVM B-tag result cut1", 10,0,10 ,nptbins,ptmin,ptmax);\r
+ fhDVMBtagCut2 = new TH2F("hdvmbtag_cut2","DVM B-tag result cut2", 10,0,10 ,nptbins,ptmin,ptmax);\r
+ fhDVMBtagCut3 = new TH2F("hdvmbtag_cut3","DVM B-tag result cut3", 10,0,10 ,nptbins,ptmin,ptmax);\r
+ fhDVMBtagQA1 = new TH2F("hdvmbtag_qa1" ,"DVM B-tag QA: pairDCA vs length", 100,0,0.2 ,100,0,1.0);\r
+ fhDVMBtagQA2 = new TH2F("hdvmbtag_qa2" ,"DVM B-tag QA: signDCA vs mass" , 200,-0.5,0.5 ,100,0,10);\r
+ fhDVMBtagQA3 = new TH1F("hdvmbtag_qa3" ,"DVM B-tag QA: ITS-Hits electron" ,7,0,7);\r
+ fhDVMBtagQA4 = new TH1F("hdvmbtag_qa4" ,"DVM B-tag QA: IP d electron" ,200,-3,3);\r
+ fhDVMBtagQA5 = new TH1F("hdvmbtag_qa5" ,"DVM B-tag QA: IP z electron" ,200,-3,3);\r
+\r
+ outputContainer->Add(fhDVMBtagCut1) ;\r
+ outputContainer->Add(fhDVMBtagCut2) ;\r
+ outputContainer->Add(fhDVMBtagCut3) ;\r
+ outputContainer->Add(fhDVMBtagQA1) ;\r
+ outputContainer->Add(fhDVMBtagQA2) ;\r
+ outputContainer->Add(fhDVMBtagQA3) ;\r
+ outputContainer->Add(fhDVMBtagQA4) ;\r
+ outputContainer->Add(fhDVMBtagQA5) ;\r
+\r
+ //IPSig B-tagging\r
+ fhIPSigBtagQA1 = new TH1F("hipsigbtag_qa1" ,"IPSig B-tag QA: # tag tracks", 20,0,20);\r
+ fhIPSigBtagQA2 = new TH1F("hipsigbtag_qa2" ,"IPSig B-tag QA: IP significance", 200,-10.,10.);\r
+ fhTagJetPt1x4 = new TH1F("hTagJetPt1x4","tagged jet pT (1 track, ipSignif>4);p_{T}",300,0.,300.);\r
+ fhTagJetPt2x3 = new TH1F("hTagJetPt2x3","tagged jet pT (2 track, ipSignif>3);p_{T}",300,0.,300.);\r
+ fhTagJetPt3x2 = new TH1F("hTagJetPt3x2","tagged jet pT (3 track, ipSignif>2);p_{T}",300,0.,300.);\r
+ fhePlusTagJetPt1x4 = new TH1F("hePlusTagJetPt1x4","tagged eJet pT (1 track, ipSignif>4);p_{T}",300,0.,300.);\r
+ fhePlusTagJetPt2x3 = new TH1F("hePlusTagJetPt2x3","tagged eJet pT (2 track, ipSignif>3);p_{T}",300,0.,300.);\r
+ fhePlusTagJetPt3x2 = new TH1F("hePlusTagJetPt3x2","tagged eJet pT (3 track, ipSignif>2);p_{T}",300,0.,300.);\r
+\r
+ outputContainer->Add(fhIPSigBtagQA1) ;\r
+ outputContainer->Add(fhIPSigBtagQA2) ;\r
+ outputContainer->Add(fhTagJetPt1x4);\r
+ outputContainer->Add(fhTagJetPt2x3);\r
+ outputContainer->Add(fhTagJetPt3x2);\r
+ outputContainer->Add(fhePlusTagJetPt1x4);\r
+ outputContainer->Add(fhePlusTagJetPt2x3);\r
+ outputContainer->Add(fhePlusTagJetPt3x2);\r
+\r
+ //B-Jet histograms\r
+ fhJetType = new TH2F("hJetType","# jets passing each tag method vs jet pt",15,0,15,300,0.,300.);\r
+ fhLeadJetType = new TH2F("hLeadJetType","# leading jets passing each tag method vs jet pt",15,0,15,300,0.,300.);\r
+ fhBJetXsiFF = new TH2F("hBJetXsiFF","B-jet #Xsi Frag. Fn.",100,0.,10.,300,0.,300.);\r
+ fhBJetPtFF = new TH2F("hBJetPtFF","B-jet p_{T} Frag. Fn.",nptbins,ptmin,ptmax,300,0.,300.);\r
+ fhBJetEtaPhi = new TH2F("hBJetEtaPhi","B-jet eta-phi distribution",netabins,etamin,etamax,nphibins,phimin,phimax);\r
+ fhNonBJetXsiFF = new TH2F("hNonBJetXsiFF","Non B-jet #Xsi Frag. Fn.",100,0.,10.,300,0.,300.);\r
+ fhNonBJetPtFF = new TH2F("hNonBJetPtFF","Non B-jet p_{T} Frag. Fn.",nptbins,ptmin,ptmax,300,0.,300.);\r
+ fhNonBJetEtaPhi = new TH2F("hNonBJetEtaPhi","Non B-jet eta-phi distribution",netabins,etamin,etamax,nphibins,phimin,phimax);\r
+\r
+ outputContainer->Add(fhJetType);\r
+ outputContainer->Add(fhLeadJetType);\r
+ outputContainer->Add(fhBJetXsiFF);\r
+ outputContainer->Add(fhBJetPtFF);\r
+ outputContainer->Add(fhBJetEtaPhi);\r
+ outputContainer->Add(fhNonBJetXsiFF);\r
+ outputContainer->Add(fhNonBJetPtFF);\r
+ outputContainer->Add(fhNonBJetEtaPhi);\r
+\r
+ //Histograms that use MC information\r
+ if(IsDataMC()){\r
+\r
+ //electron ntuple for further analysis\r
+ if(fWriteNtuple) {\r
+ fEleNtuple = new TNtuple("EleNtuple","Electron Ntuple","tmctag:cmctag:pt:phi:eta:p:E:deta:dphi:nCells:dEdx:pidProb:impXY:impZ");\r
+ outputContainer->Add(fEleNtuple) ;\r
+ }\r
+\r
+ //electrons from various MC sources\r
+ fhPhiConversion = new TH2F("hPhiConversion","Conversion Electron phi vs pT",nptbins,ptmin,ptmax,nphibins,phimin,phimax);\r
+ fhEtaConversion = new TH2F("hEtaConversion","Conversion Electron eta vs. pT",nptbins,ptmin,ptmax,netabins,etamin,etamax);\r
+\r
+ outputContainer->Add(fhPhiConversion);\r
+ outputContainer->Add(fhEtaConversion);\r
+\r
+ //Bins along y-axis are: 0 - unfiltered, 1 - bottom, 2 - charm, 3 - charm from bottom,\r
+ //4 - conversion, 5 - Dalitz, 6 - W and Z, 7 - junk/unknown, 8 - misidentified\r
+\r
+ //histograms for comparison to tracking detectors\r
+ fhPtTrack = new TH2F("hPtTrack","Track w/in EMCAL acceptance",nptbins,ptmin,ptmax,10,0,10);\r
+ 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);\r
+\r
+ outputContainer->Add(fhPtTrack);\r
+ outputContainer->Add(fhPtNPEBHadron);\r
+\r
+ //for computing efficiency of IPSig tag\r
+ fhBJetPt1x4 = new TH1F("hBJetPt1x4","tagged B-jet pT (1 track, ipSignif>4);p_{T}",300,0.,300.);\r
+ fhBJetPt2x3 = new TH1F("hBJetPt2x3","tagged B-jet pT (2 track, ipSignif>3);p_{T}",300,0.,300.);\r
+ fhBJetPt3x2 = new TH1F("hBJetPt3x2","tagged B-jet pT (3 track, ipSignif>2);p_{T}",300,0.,300.);\r
+ fhFakeJetPt1x4 = new TH1F("hFakeJetPt1x4","fake tagged B-jet pT (1 track, ipSignif>4);p_{T}",300,0.,300.);\r
+ fhFakeJetPt2x3 = new TH1F("hFakeJetPt2x3","fake tagged B-jet pT (2 track, ipSignif>3);p_{T}",300,0.,300.);\r
+ fhFakeJetPt3x2 = new TH1F("hFakeJetPt3x2","fake tagged B-jet pT (3 track, ipSignif>2);p_{T}",300,0.,300.);\r
+ fhDVMJet = new TH2F("hDVM_algo","# DVM jets passing vs Mc-Bjet",10,0,10,300,0.,300.);\r
+\r
+ outputContainer->Add(fhBJetPt1x4);\r
+ outputContainer->Add(fhBJetPt2x3);\r
+ outputContainer->Add(fhBJetPt3x2);\r
+ outputContainer->Add(fhFakeJetPt1x4);\r
+ outputContainer->Add(fhFakeJetPt2x3);\r
+ outputContainer->Add(fhFakeJetPt3x2);\r
+ outputContainer->Add(fhDVMJet);\r
+\r
+ //MC Only histograms\r
+ \r
+ //MC ele ntuple for further analysis\r
+ if(fWriteNtuple) {\r
+ fMCEleNtuple = new TNtuple("MCEleNtuple","MC Electron Ntuple","mctag:pt:phi:eta:x:y:z");\r
+ outputContainer->Add(fMCEleNtuple) ;\r
+ }\r
+\r
+ fhMCBJetElePt = new TH2F("hMCBJetElePt","MC B-jet pT vs. electron pT",300,0.,300.,300,0.,300.);\r
+ fhMCBHadronElePt = new TH2F("hMCBHadronElePt","MC B-hadron pT vs. electron pT",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);\r
+ fhPtMCHadron = new TH1F("hPtMCHadron","MC Charged hadrons w/in EMCAL acceptance",nptbins,ptmin,ptmax);\r
+\r
+ //Bins along y-axis are: 0 - unfiltered, 1 - bottom, 2 - charm, 3 - charm from bottom,\r
+ //4 - conversion, 5 - Dalitz, 6 - W and Z, 7 - junk/unknown\r
+ fhPtMCElectron = new TH2F("hPtMCElectron","MC electrons from various sources w/in EMCAL acceptance",nptbins,ptmin,ptmax,10,0,10);\r
+\r
+ fhMCXYConversion = new TH2F("hMCXYConversion","XvsY of conversion electrons",400,-400.,400.,400,-400.,400.);\r
+ fhMCRadPtConversion = new TH2F("hMCRadPtConversion","Radius vs pT of conversion electrons",200,0.,400.,nptbins,ptmin,ptmax);\r
+\r
+ outputContainer->Add(fhMCBJetElePt);\r
+ outputContainer->Add(fhMCBHadronElePt);\r
+ outputContainer->Add(fhPtMCHadron);\r
+ outputContainer->Add(fhPtMCElectron);\r
+ outputContainer->Add(fhMCXYConversion);\r
+ outputContainer->Add(fhMCRadPtConversion);\r
+\r
+ }//Histos with MC\r
+ \r
+\r
+ return outputContainer ;\r
+ \r
+}\r
+\r
+//____________________________________________________________________________\r
+void AliAnaElectron::Init()\r
+{\r
+\r
+ //do some initialization\r
+ if(fCalorimeter == "PHOS") {\r
+ printf("AliAnaElectron::Init() - !!STOP: You want to use PHOS in analysis but this is not (yet) supported!!\n!!Check the configuration file!!\n");\r
+ fCalorimeter = "EMCAL";\r
+ }\r
+ if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn()){\r
+ printf("AliAnaElectron::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!!\n!!Check the configuration file!!\n");\r
+ abort();\r
+ }\r
+\r
+}\r
+\r
+\r
+//____________________________________________________________________________\r
+void AliAnaElectron::InitParameters()\r
+{\r
+ \r
+ //Initialize the parameters of the analysis.\r
+ SetOutputAODClassName("AliAODPWG4Particle");\r
+ SetOutputAODName("PWG4Particle");\r
+\r
+ AddToHistogramsName("AnaElectron_");\r
+\r
+ fCalorimeter = "EMCAL" ;\r
+ fpOverEmin = 0.5;\r
+ fpOverEmax = 1.2;\r
+ fResidualCut = 0.02;\r
+ fMinClusEne = 4.0;\r
+ //DVM B-tagging\r
+ fDrCut = 1.0; \r
+ fPairDcaCut = 0.02;\r
+ fDecayLenCut = 1.0;\r
+ fImpactCut = 0.5;\r
+ fAssocPtCut = 1.0;\r
+ fMassCut = 1.5;\r
+ fSdcaCut = 0.1;\r
+ fITSCut = 4;\r
+ //IPSig B-tagging\r
+ fNTagTrkCut = 2;\r
+ fIPSigCut = 3.0;\r
+ //Jet fiducial cuts\r
+ fJetEtaCut = 0.3;\r
+ fJetPhiMin = 1.8;\r
+ fJetPhiMax = 2.9;\r
+}\r
+\r
+//__________________________________________________________________\r
+void AliAnaElectron::MakeAnalysisFillAOD() \r
+{\r
+ //\r
+ // Do analysis and fill aods with electron candidates\r
+ // These AODs will be used to do subsequent histogram filling\r
+ //\r
+ // Also fill some QA histograms\r
+ //\r
+\r
+ Double_t bfield = 0.;\r
+ if(GetReader()->GetDataType() != AliCaloTrackReader::kMC) bfield = GetReader()->GetBField();\r
+\r
+ //Select the calorimeter of the electron\r
+ if(fCalorimeter != "EMCAL") {\r
+ printf("This class not yet implemented for PHOS\n");\r
+ abort();\r
+ }\r
+ \r
+ TObjArray *cl = GetAODEMCAL();\r
+ \r
+ ////////////////////////////////////////////////\r
+ //Start from tracks and get associated clusters \r
+ ////////////////////////////////////////////////\r
+ if(!GetAODCTS() || GetAODCTS()->GetEntriesFast() == 0) return ;\r
+ Int_t ntracks = GetAODCTS()->GetEntriesFast();\r
+ Int_t refmult = 0; Int_t refmult2 = 0;\r
+ if(GetDebug() > 0)\r
+ printf("AliAnaElectron::MakeAnalysisFillAOD() - In CTS aod entries %d\n", ntracks);\r
+\r
+ //Unfortunately, AliAODTracks don't have associated EMCAL clusters.\r
+ //we have to redo track-matching, I guess\r
+ Int_t iCluster = -999;\r
+ Int_t bt = 0; //counter for event b-tags\r
+\r
+ for (Int_t itrk = 0; itrk < ntracks; itrk++) {////////////// track loop\r
+ iCluster = -999; //start with no match\r
+ AliAODTrack * track = (AliAODTrack*) (GetAODCTS()->At(itrk)) ;\r
+ if (TMath::Abs(track->Eta())< 0.5) refmult++;\r
+ Double_t imp[2] = {-999.,-999.}; Double_t cov[3] = {-999.,-999.,-999.};\r
+ Bool_t dcaOkay = GetDCA(track,imp,cov); //homegrown dca calculation until AOD is fixed\r
+ if(!dcaOkay) printf("AliAnaElectron::Problem computing DCA to primary vertex for track %d. Skipping it...\n",itrk);\r
+ if(TMath::Abs(track->Eta())< 0.5 && TMath::Abs(imp[0])<1.0 && TMath::Abs(imp[1])<1.0) refmult2++;\r
+ fhImpactXY->Fill(imp[0]);\r
+\r
+ //JLK CHECK\r
+ //AliESDtrack esdTrack(track);\r
+ //Double_t tpcpid[AliPID::kSPECIES];\r
+ //esdTrack.GetTPCpid(tpcpid);\r
+ //Double_t eProb = tpcpid[AliPID::kElectron];\r
+ //if(eProb > 0) printf("<%d> ESD eProb = %2.2f\n",itrk,eProb);\r
+\r
+ AliAODPid* pid = (AliAODPid*) track->GetDetPid();\r
+ if(pid == 0) {\r
+ if(GetDebug() > 0) printf("AliAnaElectron::MakeAnalysisFillAOD() - No PID object - skipping track %d",itrk);\r
+ continue;\r
+ } else {\r
+ Double_t emcpos[3];\r
+ pid->GetEMCALPosition(emcpos);\r
+ Double_t emcmom[3];\r
+ pid->GetEMCALMomentum(emcmom);\r
+ \r
+ TVector3 pos(emcpos[0],emcpos[1],emcpos[2]);\r
+ TVector3 mom(emcmom[0],emcmom[1],emcmom[2]);\r
+ Double_t tphi = pos.Phi();\r
+ Double_t teta = pos.Eta();\r
+ Double_t tmom = mom.Mag();\r
+ \r
+ //TLorentzVector mom2(mom,0.);\r
+ Bool_t in = kFALSE;\r
+ if(mom.Phi()*180./TMath::Pi() > 80. && mom.Phi()*180./TMath::Pi() < 190. &&\r
+ mom.Eta() > -0.7 && mom.Eta() < 0.7) in = kTRUE;\r
+ //Also check the track\r
+ if(track->Phi()*180./TMath::Pi() > 80. && track->Phi()*180./TMath::Pi() < 190. &&\r
+ track->Eta() > -0.7 && track->Eta() < 0.7) in = kTRUE;\r
+ ////////////////////////////\r
+ //THIS HAS A MEM LEAK JLK 24-Oct-09\r
+ //Bool_t in = GetFiducialCut()->IsInFiducialCut(mom2,fCalorimeter) ;\r
+ ///////////////////////////\r
+ 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);\r
+\r
+ if(mom.Pt() > GetMinPt() && in) {\r
+ \r
+ Double_t dEdx = pid->GetTPCsignal();\r
+ \r
+ Int_t ntot = cl->GetEntriesFast();\r
+ Double_t res = 999.;\r
+ Double_t pOverE = -999.;\r
+ \r
+ Int_t pidProb = track->GetMostProbablePID();\r
+ Bool_t tpcEle = kFALSE; if(dEdx > 70.) tpcEle = kTRUE;\r
+ Bool_t trkEle = kFALSE; if(pidProb == AliAODTrack::kElectron) trkEle = kTRUE;\r
+ Bool_t trkChgHad = kFALSE; if(pidProb == AliAODTrack::kPion || pidProb == AliAODTrack::kKaon || pidProb == AliAODTrack::kProton) trkChgHad = kTRUE;\r
+\r
+ Int_t tmctag = -1;\r
+\r
+ //Check against V0 for conversion, only if it is flagged as electron\r
+ Bool_t photonic = kFALSE;\r
+ if(tpcEle || trkEle) photonic = PhotonicV0(itrk);\r
+ if(trkEle && !photonic) fhPtNPEleTPCTRD->Fill(track->Pt(),0); //0 = no MC info\r
+ if(tpcEle && !photonic) fhPtNPEleTPC->Fill(track->Pt(),0); //0 = no MC info\r
+\r
+ if(trkChgHad) fhPtHadron->Fill(track->Pt(),0); //0 = no MC info\r
+ if(IsDataMC()) {\r
+ //Input from second AOD?\r
+ Int_t input = 0;\r
+ //if(GetReader()->GetAODCTSNormalInputEntries() <= itrk) input = 1;\r
+ tmctag = GetMCAnalysisUtils()->CheckOrigin(track->GetLabel(),GetReader(),input);\r
+\r
+ if(trkChgHad) fhPtHadron->Fill(track->Pt(),GetMCSource(tmctag));\r
+ if(tpcEle && !photonic) fhPtNPEleTPC->Fill(track->Pt(),GetMCSource(tmctag));\r
+ if(trkEle && !photonic) fhPtNPEleTPCTRD->Fill(track->Pt(),GetMCSource(tmctag));\r
+ fhPtTrack->Fill(track->Pt(),GetMCSource(tmctag));\r
+ }\r
+\r
+ Bool_t emcEle = kFALSE; \r
+ //For tracks in EMCAL acceptance, pair them with all clusters\r
+ //and fill the dEta vs dPhi for these pairs:\r
+\r
+ Double_t minR = 99;\r
+ Double_t minPe =-1;\r
+ Double_t minEp =-1;\r
+ Double_t minMult = -1;\r
+ Double_t minPt = -1;\r
+\r
+ for(Int_t iclus = 0; iclus < ntot; iclus++) {\r
+ AliVCluster * clus = (AliVCluster*) (cl->At(iclus));\r
+ if(!clus) continue;\r
+\r
+ //As of 11-Oct-2009\r
+ //only select "good" clusters \r
+ if (clus->GetNCells() < 2 ) continue;\r
+ if (clus->GetNCells() > 30 ) continue;\r
+ if (clus->E() < fMinClusEne ) continue;\r
+ if (clus->GetDispersion() > 1 ) continue;\r
+ if (clus->GetM20() > 0.4 ) continue;\r
+ if (clus->GetM02() > 0.4 ) continue;\r
+ if (clus->GetM20() < 0.03 ) continue;\r
+ if (clus->GetM02() < 0.03 ) continue;\r
+ \r
+ Float_t x[3];\r
+ clus->GetPosition(x);\r
+ TVector3 cluspos(x[0],x[1],x[2]);\r
+ Double_t deta = teta - cluspos.Eta();\r
+ Double_t dphi = tphi - cluspos.Phi();\r
+ if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi();\r
+ if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi();\r
+ fh2dEtadPhi->Fill(deta,dphi);\r
+ fh2TrackPVsClusterE->Fill(clus->E(),track->P());\r
+ fh2TrackPtVsClusterE->Fill(clus->E(),track->Pt());\r
+ fh2TrackPhiVsClusterPhi->Fill(cluspos.Phi(),mom.Phi());\r
+ fh2TrackEtaVsClusterEta->Fill(cluspos.Eta(),mom.Eta());\r
+ \r
+ res = sqrt(dphi*dphi + deta*deta);\r
+ fh1dR->Fill(res);\r
+ \r
+ /////////////////////////////////\r
+ //Perform electron cut analysis//\r
+ /////////////////////////////////\r
+ //Good match\r
+ if(res < fResidualCut) {\r
+ fh2dEtadPhiMatched->Fill(deta,dphi);\r
+ fh2MatchdEdx->Fill(track->P(),dEdx);\r
+ iCluster = iclus;\r
+\r
+ Double_t energy = clus->E(); \r
+ if(energy > 0) pOverE = tmom/energy;\r
+\r
+ if (res< minR) {\r
+ minR = res;\r
+ minPe = pOverE;\r
+ minEp = energy/tmom;\r
+ minMult = clus->GetNCells() ;\r
+ minPt = track->Pt();\r
+ }\r
+\r
+ Int_t cmctag = -1; \r
+ if(IsDataMC()) { \r
+ //Do you want the cluster or the track label?\r
+ Int_t input = 0;\r
+ //if(GetReader()->GetAODEMCALNormalInputEntries() <= iclus) input = 1;\r
+ cmctag = GetMCAnalysisUtils()->CheckOrigin(clus->GetLabel(),GetReader(),input);\r
+ }\r
+ \r
+ if(fWriteNtuple) {\r
+ fEleNtuple->Fill(tmctag,cmctag,track->Pt(),track->Phi(),track->Eta(),track->P(),clus->E(),deta,dphi,clus->GetNCells(),dEdx,pidProb,imp[0],imp[1]);\r
+ }\r
+\r
+ } else {\r
+ fh2dEtadPhiUnmatched->Fill(deta,dphi);\r
+ }//res cut\r
+ }//calo cluster loop\r
+\r
+ fh3pOverE->Fill(minPt,minPe ,minMult);\r
+ fh3EOverp->Fill(minPt,minEp ,minMult);\r
+ if (trkEle) {\r
+ fh3pOverE2->Fill(minPt,minPe ,minMult);\r
+ fh3EOverp2->Fill(minPt,minEp ,minMult);\r
+ }\r
+ if (tpcEle) {\r
+ fh3pOverE3->Fill(minPt,minPe ,minMult);\r
+ fh3EOverp3->Fill(minPt,minEp ,minMult);\r
+ }\r
+ //new\r
+ if (tmctag>-1 && GetMCSource(tmctag)<8 ) {\r
+ fh2pOverE->Fill(minPt,minPe );\r
+ fh2EOverp->Fill(minPt,minEp );\r
+ if (trkEle) {\r
+ fh2pOverE2->Fill(minPt,minPe );\r
+ fh2EOverp2->Fill(minPt,minEp );\r
+ }\r
+ }\r
+\r
+ //////////////////////////////\r
+ //Electron cuts happen here!//\r
+ //////////////////////////////\r
+ if(minPe > fpOverEmin && minPe < fpOverEmax) emcEle = kTRUE;\r
+\r
+ ///////////////////////////\r
+ //Fill AOD with electrons//\r
+ ///////////////////////////\r
+ //Take all emcal electrons, but the others only if pT < 10 GeV\r
+ if(emcEle || ( (tpcEle || trkEle) && track->Pt() < 10.) ) {\r
+\r
+ //B-tagging\r
+ if(GetDebug() > 1) printf("Found Electron - do b-tagging\n");\r
+ Int_t dvmbtag = GetDVMBtag(track); bt += dvmbtag;\r
+\r
+ fh2EledEdx->Fill(track->P(),dEdx);\r
+ \r
+ Double_t eMass = 0.511/1000; //mass in GeV\r
+ Double_t eleE = sqrt(track->P()*track->P() + eMass*eMass);\r
+ AliAODPWG4Particle tr = AliAODPWG4Particle(track->Px(),track->Py(),track->Pz(),eleE);\r
+ tr.SetLabel(track->GetLabel());\r
+ tr.SetCaloLabel(iCluster,-1); //sets the indices of the original caloclusters\r
+ tr.SetTrackLabel(itrk,-1); //sets the indices of the original tracks\r
+ if(emcEle) {//PID determined by EMCAL\r
+ tr.SetDetector(fCalorimeter);\r
+ } else {\r
+ tr.SetDetector("CTS"); //PID determined by CTS\r
+ }\r
+\r
+ //if(GetReader()->GetAODCTSNormalInputEntries() <= itrk) tr.SetInputFileIndex(1);\r
+ //Make this preserve sign of particle\r
+ if(track->Charge() < 0) tr.SetPdg(11); //electron is 11\r
+ else tr.SetPdg(-11); //positron is -11\r
+ Int_t btag = 0;\r
+ if(dvmbtag > 0) tr.SetBTagBit(btag,tr.kDVMTag0);\r
+ if(dvmbtag > 1) tr.SetBTagBit(btag,tr.kDVMTag1);\r
+ if(dvmbtag > 2) tr.SetBTagBit(btag,tr.kDVMTag2);\r
+ tr.SetBtag(btag);\r
+ \r
+ //Play with the MC stack if available\r
+ //Check origin of the candidates\r
+ if(IsDataMC()){\r
+ \r
+ //FIXME: Need to re-think this for track-oriented analysis\r
+ //JLK DO WE WANT TRACK TAG OR CLUSTER TAG?\r
+ tr.SetTag(GetMCAnalysisUtils()->CheckOrigin(tr.GetLabel(),GetReader(),tr.GetInputFileIndex()));\r
+ \r
+ if(GetDebug() > 0) printf("AliAnaElectron::MakeAnalysisFillAOD() - Origin of candidate %d\n",tr.GetTag());\r
+ }//Work with stack also \r
+ \r
+ AddAODParticle(tr);\r
+ \r
+ if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillAOD() - Electron selection cuts passed: pT %3.2f, pdg %d\n",tr.Pt(),tr.GetPdg()); \r
+ }//electron\r
+ }//pt, fiducial selection\r
+ }//pid check\r
+ }//track loop \r
+ \r
+ fhRefMult->Fill(refmult);\r
+ fhRefMult2->Fill(refmult2);\r
+\r
+ if(GetDebug() > 1 && bt > 0) printf("AliAnaElectron::MakeAnalysisFillAOD() *** Event Btagged *** \n");\r
+ if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillAOD() End fill AODs \n"); \r
+ \r
+}\r
+\r
+//__________________________________________________________________\r
+void AliAnaElectron::MakeAnalysisFillHistograms() \r
+{\r
+ //Do analysis and fill histograms\r
+\r
+ AliStack * stack = 0x0;\r
+ TParticle * primary = 0x0;\r
+ AliAODMCParticle * aodprimary = 0x0;\r
+\r
+ Int_t ph1 = 0; //photonic 1 count\r
+ Int_t ph2 = 0; //photonic 2 count\r
+ Int_t phB = 0; //both count\r
+\r
+ if(IsDataMC()) {\r
+ if(GetReader()->ReadStack()){\r
+ stack = GetMCStack() ; \r
+ if(!stack)\r
+ printf("AliAnaElectron::MakeAnalysisFillHistograms() *** no stack ***: \n");\r
+ \r
+ }\r
+ }// is data and MC\r
+\r
+ ////////////////////////////////////\r
+ //Loop over jets and check for b-tag\r
+ ////////////////////////////////////\r
+ Int_t njets = (GetReader()->GetOutputEvent())->GetNJets();\r
+ if(njets > 0) {\r
+ if(GetDebug() > 0) printf("AliAnaElectron::MakeAnalysisFillHistograms() - Jet AOD branch has %d jets. Performing b-jet tag analysis\n",njets);\r
+\r
+ for(Int_t ijet = 0; ijet < njets ; ijet++) {\r
+ AliAODJet * jet = (AliAODJet*)(GetReader()->GetOutputEvent())->GetJet(ijet) ;\r
+ //Only consider jets with pt > 10 GeV (the rest have to be junk)\r
+ //printf("AODJet<%d> pt = %2.2f\n",ijet,jet->Pt());\r
+ if(jet->Pt() < 10.) continue;\r
+\r
+ if(GetDebug() > 3) {\r
+ printf("AliAODJet ijet = %d\n",ijet);\r
+ jet->Print("");\r
+ }\r
+ //Skip jets not inside a smaller fiducial volume to ensure that\r
+ //they are completely contained in the EMCAL\r
+ if(TMath::Abs(jet->Eta()) > fJetEtaCut) continue;\r
+ if(jet->Phi() < fJetPhiMin || jet->Phi() > fJetPhiMax) continue;\r
+\r
+ //To "tag" the jet, we will look for it to pass our various criteria\r
+ //For e jet tag, we just look to see which ones have NPEs\r
+ //For DVM jet tag, we will look for DVM electrons\r
+ //For IPSig, we compute the IPSig for all tracks and if the\r
+ //number passing is above the cut, it passes\r
+ Bool_t leadJet = kFALSE;\r
+ if (ijet==0) leadJet= kTRUE;\r
+ Bool_t eJet = kFALSE; \r
+ Bool_t eJet2 = kFALSE; //electron triggered\r
+ Bool_t hadJet = kFALSE; //hadron triggered \r
+ Bool_t dvmJet = kFALSE; \r
+ Bool_t ipsigJet = kFALSE;\r
+ TRefArray* rt = jet->GetRefTracks();\r
+ Int_t ntrk = rt->GetEntries();\r
+ Int_t trackCounter[4] = {0,0,0,0}; //for ipsig\r
+ for(Int_t itrk = 0; itrk < ntrk; itrk++) {\r
+ AliAODTrack* jetTrack = (AliAODTrack*)jet->GetTrack(itrk);\r
+ if( GetIPSignificance(jetTrack, jet->Phi()) > fIPSigCut) trackCounter[0]++;\r
+ if( GetIPSignificance(jetTrack, jet->Phi()) > 4.) trackCounter[1]++;\r
+ if( GetIPSignificance(jetTrack, jet->Phi()) > 3.) trackCounter[2]++;\r
+ if( GetIPSignificance(jetTrack, jet->Phi()) > 2.) trackCounter[3]++;\r
+ Bool_t isNPE = CheckTrack(jetTrack,"NPE");\r
+ if(isNPE) eJet = kTRUE;\r
+ if ( isNPE && jetTrack->Pt()>10.0 ) eJet2 =kTRUE;\r
+ if (!isNPE && jetTrack->Pt()>10.0) hadJet =kTRUE;\r
+ Bool_t isDVM = CheckTrack(jetTrack,"DVM");\r
+ if(isDVM) dvmJet = kTRUE;\r
+ }\r
+ fhIPSigBtagQA1->Fill(trackCounter[0]);\r
+ if(trackCounter[1]>0) fhTagJetPt1x4->Fill(jet->Pt());\r
+ if(trackCounter[2]>1) fhTagJetPt2x3->Fill(jet->Pt());\r
+ if(trackCounter[3]>2) fhTagJetPt3x2->Fill(jet->Pt());\r
+\r
+ if(trackCounter[1]>0 && eJet) fhePlusTagJetPt1x4->Fill(jet->Pt());\r
+ if(trackCounter[2]>1 && eJet) fhePlusTagJetPt2x3->Fill(jet->Pt());\r
+ if(trackCounter[3]>2 && eJet) fhePlusTagJetPt3x2->Fill(jet->Pt());\r
+\r
+ if(trackCounter[0] > fNTagTrkCut) ipsigJet = kTRUE;\r
+\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());\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
+\r
+ if(isTrueBjet && dvmJet) fhDVMJet->Fill(0.,jet->Pt()); // good tagged\r
+ if(isTrueBjet && !dvmJet) fhDVMJet->Fill(1.,jet->Pt()); // missed tagged\r
+ if(!isTrueBjet && dvmJet) fhDVMJet->Fill(2.,jet->Pt()); // fake tagged\r
+ if(!isTrueBjet && !dvmJet) fhDVMJet->Fill(3.,jet->Pt()); // others\r
+ if(isTrueDjet && !isTrueBjet && dvmJet) fhDVMJet->Fill(4.,jet->Pt()); // charm-tagged\r
+ if(isTrueDjet && !isTrueBjet && !dvmJet) fhDVMJet->Fill(5.,jet->Pt()); // charm -not tagged\r
+ if(!(isTrueDjet||isTrueBjet ) && dvmJet) fhDVMJet->Fill(6.,jet->Pt()); // light flavor -tagged\r
+ if(!(isTrueDjet||isTrueBjet ) && !dvmJet) fhDVMJet->Fill(7.,jet->Pt()); // light flavor -not tagged\r
+ if(isTrueBjet && eJet && dvmJet) fhDVMJet->Fill(8.,jet->Pt()); // bjet with electron\r
+ if(isTrueBjet && !eJet && dvmJet) fhDVMJet->Fill(9.,jet->Pt()); // needs more thought\r
+\r
+ if(isTrueBjet) {\r
+ if(trackCounter[1]>0) fhBJetPt1x4->Fill(jet->Pt());\r
+ if(trackCounter[2]>1) fhBJetPt2x3->Fill(jet->Pt());\r
+ if(trackCounter[3]>2) fhBJetPt3x2->Fill(jet->Pt());\r
+ } else {\r
+ if(trackCounter[1]>0) fhFakeJetPt1x4->Fill(jet->Pt());\r
+ if(trackCounter[2]>1) fhFakeJetPt2x3->Fill(jet->Pt());\r
+ if(trackCounter[3]>2) fhFakeJetPt3x2->Fill(jet->Pt());\r
+ }\r
+ }\r
+\r
+ //Fill bjet histograms here\r
+ if(!(eJet || ipsigJet || dvmJet)) fhJetType->Fill(0.,jet->Pt()); //none\r
+ if(eJet && !(ipsigJet || dvmJet)) fhJetType->Fill(1.,jet->Pt()); //only ejet\r
+ if(dvmJet && !(eJet || ipsigJet)) fhJetType->Fill(2.,jet->Pt()); //only dvm\r
+ if(ipsigJet && !(eJet || dvmJet)) fhJetType->Fill(3.,jet->Pt()); //only ipsig\r
+ if(eJet && dvmJet && !ipsigJet) fhJetType->Fill(4.,jet->Pt()); //ejet & dvm\r
+ if(eJet && ipsigJet && !dvmJet) fhJetType->Fill(5.,jet->Pt()); //ejet & ipsig\r
+ if(dvmJet && ipsigJet && !eJet) fhJetType->Fill(6.,jet->Pt()); //dvm & ipsig\r
+ if(dvmJet && ipsigJet && eJet) fhJetType->Fill(7.,jet->Pt()); //all\r
+ if(dvmJet || ipsigJet || eJet) fhJetType->Fill(8.,jet->Pt()); //any of them\r
+ if(eJet ) fhJetType->Fill(9.,jet->Pt()); //any of them\r
+ if(dvmJet) fhJetType->Fill(10.,jet->Pt()); //any of them\r
+ if(eJet2 ) fhJetType->Fill(11.,jet->Pt()); //any of them\r
+ if(hadJet) fhJetType->Fill(12.,jet->Pt()); //any of them\r
+\r
+ if(eJet || ipsigJet || dvmJet) fhBJetEtaPhi->Fill(jet->Eta(),jet->Phi());\r
+ else fhNonBJetEtaPhi->Fill(jet->Eta(),jet->Phi());\r
+\r
+ //leading jets\r
+ if (leadJet) {\r
+ fhLeadJetType->Fill(0.,jet->Pt()); //all\r
+ if(eJet ) fhLeadJetType->Fill(1.,jet->Pt());\r
+ if(eJet2 ) fhLeadJetType->Fill(2.,jet->Pt());\r
+ if(hadJet ) fhLeadJetType->Fill(3.,jet->Pt());\r
+ if(eJet && (dvmJet || ipsigJet) ) fhLeadJetType->Fill(4.,jet->Pt());\r
+ if(eJet2 && (dvmJet || ipsigJet) ) fhLeadJetType->Fill(5.,jet->Pt());\r
+ if(hadJet && (dvmJet || ipsigJet) ) fhLeadJetType->Fill(6.,jet->Pt());\r
+ }\r
+\r
+ for(Int_t itrk = 0; itrk < ntrk; itrk++) {\r
+ AliAODTrack* jetTrack = (AliAODTrack*)jet->GetTrack(itrk);\r
+ Double_t xsi = TMath::Log(jet->Pt()/jetTrack->Pt());\r
+ if(eJet || ipsigJet || dvmJet) {\r
+ if(GetDebug() > 0) printf("AliAnaElectron::MakeAnalysisFillHistograms - We have a bjet!\n");\r
+ fhBJetXsiFF->Fill(xsi,jet->Pt());\r
+ fhBJetPtFF->Fill(jetTrack->Pt(),jet->Pt());\r
+ } else {\r
+ //Fill non-bjet histograms here\r
+ fhNonBJetXsiFF->Fill(xsi,jet->Pt());\r
+ fhNonBJetPtFF->Fill(jetTrack->Pt(),jet->Pt());\r
+ }\r
+ }\r
+\r
+ } //jet loop\r
+ } //jets exist\r
+ \r
+ //////////////////////////////\r
+ //Loop on stored AOD electrons\r
+ //////////////////////////////\r
+ Int_t naod = GetOutputAODBranch()->GetEntriesFast();\r
+ if(GetDebug() > 0) printf("AliAnaElectron::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);\r
+ \r
+ for(Int_t iaod = 0; iaod < naod ; iaod++){\r
+ AliAODPWG4Particle* ele = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));\r
+ Int_t pdg = ele->GetPdg();\r
+ \r
+ if(GetDebug() > 3) \r
+ printf("AliAnaElectron::MakeAnalysisFillHistograms() - PDG %d, MC TAG %d, Calorimeter %s\n", ele->GetPdg(),ele->GetTag(), (ele->GetDetector()).Data()) ;\r
+ \r
+ if(TMath::Abs(pdg) != AliCaloPID::kElectron) continue; \r
+ \r
+ if(GetDebug() > 1) \r
+ printf("AliAnaElectron::MakeAnalysisFillHistograms() - ID Electron: pt %f, phi %f, eta %f\n", ele->Pt(),ele->Phi(),ele->Eta()) ;\r
+\r
+ //MC tag of this electron\r
+ Int_t mctag = ele->GetTag();\r
+\r
+ //Filter for photonic electrons based on opening angle and Minv\r
+ //cuts, also fill histograms\r
+ Bool_t photonic = kFALSE;\r
+ Bool_t photonic1 = kFALSE;\r
+ photonic1 = PhotonicPrim(ele); //check against primaries\r
+ if(photonic1) ph1++;\r
+ Bool_t photonic2 = kFALSE;\r
+ photonic2 = PhotonicV0(ele->GetTrackLabel(0)); //check against V0s\r
+ if(photonic2) ph2++;\r
+ if(photonic1 && photonic2) phB++;\r
+ if(photonic1 || photonic2) photonic = kTRUE;\r
+\r
+ //Fill electron histograms \r
+ Float_t ptele = ele->Pt();\r
+ Float_t phiele = ele->Phi();\r
+ Float_t etaele = ele->Eta();\r
+\r
+ //"Best reconstructed electron spectrum" = EMCAL or tracking\r
+ //detectors say it is an electron and it does not form a V0\r
+ //with Minv near a relevant resonance\r
+ if(!photonic) {\r
+ fhPtNPEleTTE->Fill(ptele,0); //0 = no MC info\r
+ if(ele->GetDetector() == fCalorimeter) fhPtNPEleEMCAL->Fill(ptele,0);\r
+ if(IsDataMC()) {\r
+ fhPtNPEleTTE->Fill(ptele,GetMCSource(mctag));\r
+ if(ele->GetDetector() == "EMCAL") fhPtNPEleEMCAL->Fill(ptele,GetMCSource(mctag));\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
+\r
+ //kept for historical reasons?\r
+ fhPtElectron ->Fill(ptele);\r
+ fhPhiElectron ->Fill(ptele,phiele);\r
+ fhEtaElectron ->Fill(ptele,etaele);\r
+\r
+ if(photonic) {\r
+ fhPtPE->Fill(ptele);\r
+ fhPhiPE->Fill(ptele,phiele);\r
+ fhEtaPE->Fill(ptele,etaele);\r
+ } else {\r
+ fhPtNPE->Fill(ptele);\r
+ fhPhiNPE->Fill(ptele,phiele);\r
+ fhEtaNPE->Fill(ptele,etaele);\r
+ }\r
+\r
+ if(IsDataMC()){\r
+ if(GetMCAnalysisUtils()->CheckTagBit(mctag,AliMCAnalysisUtils::kMCConversion)){\r
+ fhPhiConversion ->Fill(ptele,phiele);\r
+ fhEtaConversion ->Fill(ptele,etaele);\r
+ }\r
+ }//Histograms with MC\r
+ \r
+ }// aod loop\r
+\r
+ ////////////////////////////////////////////////////////\r
+ //Fill histograms of pure MC kinematics from the stack//\r
+ ////////////////////////////////////////////////////////\r
+ if(IsDataMC()) {\r
+\r
+ //MC Jets\r
+ TVector3 bjetVect[4];\r
+ Int_t nPythiaGenJets = 0;\r
+ AliGenPythiaEventHeader* pythiaGenHeader = (AliGenPythiaEventHeader*)GetReader()->GetGenEventHeader();\r
+ if(pythiaGenHeader){\r
+ //Get Jets from MC header\r
+ nPythiaGenJets = pythiaGenHeader->NTriggerJets();\r
+ Int_t iCount = 0;\r
+ for(int ip = 0;ip < nPythiaGenJets;++ip){\r
+ if (iCount>3) break;\r
+ Float_t p[4];\r
+ pythiaGenHeader->TriggerJet(ip,p);\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())) {\r
+ bjetVect[iCount].SetXYZ(p[0], p[1], p[2]);\r
+ iCount++;\r
+ }\r
+ }\r
+ }\r
+\r
+ Int_t nPart = GetNumAODMCParticles();\r
+ if(GetReader()->ReadStack()) nPart = stack->GetNtrack();\r
+\r
+ for(Int_t ipart = 0; ipart < nPart; ipart++) {\r
+\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
+ 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
+ } 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
+ TLorentzVector mom(px,py,pz,e);\r
+ TLorentzVector pos(vx,vy,vz,vt);\r
+ Bool_t in = kFALSE;\r
+ if(mom.Phi()*180./TMath::Pi() > 80. && mom.Phi()*180./TMath::Pi() < 190. &&\r
+ mom.Eta() > -0.7 && mom.Eta() < 0.7) in = kTRUE;\r
+ /////////////////////////////////\r
+ //THIS HAS A MEM LEAK JLK 24-Oct-09\r
+ //Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter);\r
+ ////////////////////////////////\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
+ //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 ((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
+ if(GetMCSource(mctag) == 4) {//conversion\r
+ fhMCXYConversion->Fill(vx,vy);\r
+ fhMCRadPtConversion->Fill(TMath::Sqrt(vx*vx+vy*vy),mom.Pt());\r
+ }\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
+ printf("\tAliAnaElectron::Photonic electron counts: ph1 %d, ph2 %d, Both %d\n",ph1,ph2,phB);\r
+}\r
+\r
+//__________________________________________________________________\r
+Int_t AliAnaElectron::GetDVMBtag(AliAODTrack * tr )\r
+{\r
+ //This method uses the Displaced Vertex between electron-hadron\r
+ //pairs and the primary vertex to determine whether an electron is\r
+ //likely from a B hadron.\r
+\r
+ Int_t ncls1 = 0;\r
+ for(Int_t l = 0; l < 6; l++) if(TESTBIT(tr->GetITSClusterMap(),l)) ncls1++;\r
+\r
+ fhDVMBtagQA3->Fill(ncls1);\r
+ if (ncls1 < fITSCut) return 0;\r
+\r
+ Double_t imp[2] = {-999.,-999.}; Double_t cov[3] = {-999.,-999.,-999.};\r
+ Bool_t dcaOkay = GetDCA(tr,imp,cov); //homegrown dca calculation until AOD is fixed \r
+ if(!dcaOkay) {\r
+ printf("AliAnaElectron::Problem computing DCA to primary vertex for track %d",tr->GetID());\r
+ return 0;\r
+ }\r
+\r
+ fhDVMBtagQA4->Fill(imp[0]);\r
+ if (TMath::Abs(imp[0]) > fImpactCut ) return 0;\r
+ fhDVMBtagQA5->Fill(imp[1]);\r
+ if (TMath::Abs(imp[1]) > fImpactCut ) return 0;\r
+\r
+ Int_t nvtx1 = 0;\r
+ Int_t nvtx2 = 0;\r
+ Int_t nvtx3 = 0;\r
+\r
+ for (Int_t k2 =0; k2 < GetAODCTS()->GetEntriesFast() ; k2++) {\r
+ //loop over assoc\r
+ AliAODTrack* track2 = (AliAODTrack*)GetAODCTS()->At(k2);\r
+ Int_t id1 = tr->GetID();\r
+ Int_t id2 = track2->GetID();\r
+ if(id1 == id2) continue;\r
+\r
+ Int_t ncls2 = 0;\r
+ for(Int_t l = 0; l < 6; l++) if(TESTBIT(track2->GetITSClusterMap(),l)) ncls2++;\r
+ if (ncls2 < fITSCut) continue;\r
+\r
+ if(track2->Pt() < fAssocPtCut) continue;\r
+\r
+ Double_t dphi = tr->Phi() - track2->Phi();\r
+ if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi();\r
+ if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi();\r
+ Double_t deta = tr->Eta() - track2->Eta();\r
+ Double_t dr = sqrt(deta*deta + dphi*dphi);\r
+\r
+ if(dr > fDrCut) continue;\r
+ \r
+ Double_t sDca1 = ComputeSignDca(tr, track2, 1.0);\r
+ if (sDca1 > fSdcaCut) nvtx1++;\r
+ Double_t sDca2 = ComputeSignDca(tr, track2, 1.5);\r
+ if (sDca2 > fSdcaCut) nvtx2++;\r
+ Double_t sDca3 = ComputeSignDca(tr, track2, 1.8);\r
+ if (sDca3 > fSdcaCut) nvtx3++;\r
+\r
+ } //loop over hadrons\r
+\r
+ if(GetDebug() > 0) {\r
+ if (nvtx1>0) printf("result1 of btagging: %d \n",nvtx1);\r
+ if (nvtx2>0) printf("result2 of btagging: %d \n",nvtx2);\r
+ if (nvtx3>0) printf("result3 of btagging: %d \n",nvtx3);\r
+ }\r
+\r
+ //fill QA histograms\r
+ fhDVMBtagCut1->Fill(nvtx1,tr->Pt());\r
+ fhDVMBtagCut2->Fill(nvtx2,tr->Pt());\r
+ fhDVMBtagCut3->Fill(nvtx3,tr->Pt());\r
+\r
+ return nvtx2;\r
+\r
+}\r
+\r
+//__________________________________________________________________\r
+Double_t AliAnaElectron::ComputeSignDca(AliAODTrack *tr, AliAODTrack *tr2 , float masscut)\r
+{\r
+ //Compute the signed dca between two tracks\r
+ //and return the result\r
+\r
+ Double_t signDca=-999.;\r
+ if(GetDebug() > 2 ) printf(">>ComputeSdca:: track1 %d, track2 %d, masscut %f \n", tr->GetLabel(), tr2->GetLabel(), masscut);\r
+\r
+ //=====Now calculate DCA between both tracks======= \r
+ Double_t massE = 0.000511;\r
+ Double_t massK = 0.493677;\r
+\r
+ Double_t vertex[3] = {-999.,-999.,-999}; //vertex\r
+ if(GetReader()->GetDataType() != AliCaloTrackReader::kMC) {\r
+ GetVertex(vertex); //If only one file, get the vertex from there\r
+ //FIXME: Add a check for whether file 2 is PYTHIA or HIJING\r
+ //If PYTHIA, then set the vertex from file 2, if not, use the\r
+ //vertex from file 1\r
+ //if(GetReader()->GetSecondInputAODTree()) GetReader()->GetSecondInputAODVertex(vertex);\r
+ }\r
+ \r
+ TVector3 primV(vertex[0],vertex[1],vertex[2]) ;\r
+\r
+ if(GetDebug() > 5) printf(">>ComputeSdca:: primary vertex = %2.2f,%2.2f,%2.2f \n",vertex[0],vertex[1],vertex[2]) ;\r
+\r
+ AliExternalTrackParam *param1 = new AliExternalTrackParam(tr);\r
+ AliExternalTrackParam *param2 = new AliExternalTrackParam(tr2);\r
+\r
+ Double_t bfield[3];\r
+ param1->GetBxByBz(bfield);\r
+ Double_t bz = param1->GetBz();\r
+\r
+ Double_t xplane1 = 0.; Double_t xplane2 = 0.;\r
+ Double_t pairdca = param1->GetDCA(param2,bz,xplane1,xplane2);\r
+\r
+ param1->PropagateToBxByBz(xplane1,bfield);\r
+ param2->PropagateToBxByBz(xplane2,bfield);\r
+\r
+ Int_t id1 = 0, id2 = 0;\r
+ AliESDv0 bvertex(*param1,id1,*param2,id2);\r
+ Double_t vx,vy,vz;\r
+ bvertex.GetXYZ(vx,vy,vz);\r
+\r
+ Double_t emom[3];\r
+ Double_t hmom[3];\r
+ param1->PxPyPz(emom);\r
+ param2->PxPyPz(hmom);\r
+ TVector3 emomAtB(emom[0],emom[1],emom[2]);\r
+ TVector3 hmomAtB(hmom[0],hmom[1],hmom[2]);\r
+ TVector3 secvtxpt(vx,vy,vz);\r
+ TVector3 decayvector(0,0,0);\r
+ decayvector = secvtxpt - primV; //decay vector from PrimVtx\r
+ Double_t decaylength = decayvector.Mag();\r
+\r
+ printf("\t JLK pairDCA = %2.2f\n",pairdca);\r
+\r
+ if(GetDebug() > 0) {\r
+ printf(">>ComputeSdca:: mom1=%f, mom2=%f \n", emomAtB.Perp(), hmomAtB.Perp() );\r
+ printf(">>ComputeSdca:: pairDCA=%f, length=%f \n", pairdca,decaylength );\r
+ }\r
+\r
+ if (masscut<1.1) fhDVMBtagQA1->Fill(pairdca,decaylength);\r
+\r
+ if (emomAtB.Mag()>0 && pairdca < fPairDcaCut && decaylength < fDecayLenCut ) {\r
+ TVector3 sumMom = emomAtB+hmomAtB;\r
+ Double_t ener1 = sqrt(pow(emomAtB.Mag(),2) + massE*massE);\r
+ Double_t ener2 = sqrt(pow(hmomAtB.Mag(),2) + massK*massK);\r
+ Double_t ener3 = sqrt(pow(hmomAtB.Mag(),2) + massE*massE);\r
+ Double_t mass = sqrt(pow((ener1+ener2),2) - pow(sumMom.Mag(),2));\r
+ Double_t massPhot = sqrt(pow((ener1+ener3),2) - pow(sumMom.Mag(),2));\r
+ Double_t sDca = decayvector.Dot(emomAtB)/emomAtB.Mag();\r
+\r
+ if (masscut<1.1) fhDVMBtagQA2->Fill(sDca, mass);\r
+\r
+ if (mass > masscut && massPhot > 0.1) signDca = sDca;\r
+ \r
+ if(GetDebug() > 0) printf("\t>>ComputeSdca:: mass=%f \n", mass);\r
+ if(GetDebug() > 0) printf("\t>>ComputeSdca:: sec vtx-signdca :%f\n",signDca);\r
+ }\r
+\r
+ //clean up\r
+ delete param1;\r
+ delete param2;\r
+\r
+ return signDca;\r
+}\r
+\r
+//__________________________________________________________________\r
+Double_t AliAnaElectron::GetIPSignificance(AliAODTrack *tr, Double_t jetPhi)\r
+{\r
+ //get signed impact parameter significance of the given AOD track\r
+ //for the given jet\r
+\r
+ Int_t trackIndex = 0;\r
+ Int_t ntrk = GetAODCTS()->GetEntriesFast();\r
+ for (Int_t k2 =0; k2 < ntrk ; k2++) {\r
+ //loop over assoc\r
+ AliAODTrack* track2 = (AliAODTrack*)GetAODCTS()->At(k2);\r
+ int id1 = tr->GetID();\r
+ int id2 = track2->GetID();\r
+ if(id1 == id2) {\r
+ trackIndex = k2;//FIXME: check if GetAODCTS stores tracks in the\r
+ //same order of the event\r
+ break;\r
+ }\r
+ }\r
+\r
+ Double_t significance=0;\r
+ Double_t maxD = 10000.;\r
+ Double_t impPar[] = {0,0};\r
+ Double_t ipCov[]={0,0,0};\r
+ Double_t ipVec2D[] = {0,0};\r
+\r
+ AliVEvent* vEvent = (AliVEvent*)GetReader()->GetInputEvent();\r
+ if(!vEvent) return -97;\r
+ AliVVertex* vv = (AliVVertex*)vEvent->GetPrimaryVertex();\r
+ if(!vv) return -98;\r
+ AliVTrack* vTrack = (AliVTrack*)vEvent->GetTrack(trackIndex);\r
+ if(!vTrack) return -99;\r
+ AliESDtrack esdTrack(vTrack);\r
+ Double_t bfield[3];\r
+ esdTrack.GetBxByBz(bfield);\r
+ if(!esdTrack.PropagateToDCABxByBz(vv, bfield, maxD, impPar, ipCov)) return -100;\r
+ if(ipCov[0]<0) return -101;\r
+\r
+ Double_t Pxy[] = {esdTrack.Px(), esdTrack.Py()};\r
+ Double_t Txy[] = {esdTrack.Xv(), esdTrack.Yv()};\r
+ Double_t Vxy[] = {vv->GetX(), vv->GetY()};\r
+ GetImpactParamVect(Pxy, Txy, Vxy, ipVec2D);\r
+ Double_t phiIP = TMath::ATan2(ipVec2D[1], ipVec2D[0]) + (TMath::Abs(ipVec2D[1])-ipVec2D[1])/TMath::Abs(ipVec2D[1])*TMath::Pi();\r
+ Double_t cosTheta = TMath::Cos(jetPhi - phiIP);\r
+ Double_t sign = cosTheta/TMath::Abs(cosTheta);\r
+ significance = TMath::Abs(impPar[0])/TMath::Sqrt(ipCov[0])*sign;\r
+ printf("\t JLK significance = %2.2f\n",significance);\r
+ //ip = fabs(impPar[0]);\r
+ fhIPSigBtagQA2->Fill(significance);\r
+ return significance;\r
+}\r
+\r
+//__________________________________________________________________\r
+void AliAnaElectron::GetImpactParamVect(Double_t Pxy[2], Double_t Txy[2], Double_t Vxy[2], Double_t IPxy[2])\r
+{\r
+ //px,py: momentum components at the origin of the track; tx, ty:\r
+ //origin (x,y) of track; vx, vy: coordinates of primary vertex\r
+ // analytical geometry auxiliary variables\r
+ Double_t mr = Pxy[1]/Pxy[0]; //angular coeficient of the straight\r
+ //line that lies on top of track\r
+ //momentum\r
+ Double_t br = Txy[1] - mr*Txy[0]; //linear coeficient of the straight\r
+ //line that lies on top of track\r
+ //momentum\r
+ Double_t ms = -1./mr; //angular coeficient of the straight line that\r
+ //lies on top of the impact parameter line\r
+ // Double_t bs = Vxy[1] - ms*Vxy[0]; //linear coeficient of the straight\r
+ //line that lies on top of the\r
+ //impact parameter line \r
+ Double_t xIntersection = (mr*Txy[0] - ms*Vxy[0] + Vxy[1] - Txy[1])/(mr - ms);\r
+ Double_t yIntersection = mr*xIntersection + br;\r
+ //if(ceil(10000*yIntersection) - ceil(10000*(ms*xIntersection + bs))\r
+ //!= 0 )cout<<yIntersection<<", "<<ms*xIntersection + bs<<endl;\r
+ IPxy[0] = xIntersection - Vxy[0];\r
+ IPxy[1] = yIntersection - Vxy[1];\r
+ return;\r
+}\r
+\r
+//__________________________________________________________________\r
+Bool_t AliAnaElectron::PhotonicPrim(const AliAODPWG4Particle* part) \r
+{\r
+ //This method checks the opening angle and invariant mass of\r
+ //electron pairs within the AliAODPWG4Particle list to see if \r
+ //they are likely to be photonic electrons\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 massPhi = 1.020;\r
+\r
+ Int_t pdg1 = part->GetPdg();\r
+ Int_t trackId = part->GetTrackLabel(0);\r
+ AliAODTrack* track = (AliAODTrack*)GetAODCTS()->At(trackId);\r
+ if(!track) {\r
+ if(GetDebug() > 0) printf("AliAnaElectron::PhotonicPrim - can't get the AOD Track from the particle! Skipping the photonic check");\r
+ return kFALSE; //Don't proceed because we can't get the track\r
+ }\r
+\r
+ AliExternalTrackParam *param1 = new AliExternalTrackParam(track);\r
+\r
+ //Loop on stored AOD electrons and compute the angle differences and Minv\r
+ for (Int_t k2 =0; k2 < GetOutputAODBranch()->GetEntriesFast() ; k2++) {\r
+ AliAODPWG4Particle* part2 = (AliAODPWG4Particle*) GetOutputAODBranch()->At(k2);\r
+ Int_t track2Id = part2->GetTrackLabel(0);\r
+ if(trackId == track2Id) continue;\r
+ Int_t pdg2 = part2->GetPdg();\r
+ if(TMath::Abs(pdg2) != AliCaloPID::kElectron) continue;\r
+ if(part2->GetDetector() != fCalorimeter) continue;\r
+\r
+ //JLK: Check opp. sign pairs only\r
+ if(pdg1*pdg2 > 0) continue; //skip same-sign pairs\r
+\r
+ //propagate to common vertex and check opening angle\r
+ AliAODTrack* track2 = (AliAODTrack*)GetAODCTS()->At(track2Id);\r
+ if(!track2) {\r
+ if(GetDebug() >0) printf("AliAnaElectron::PhotonicPrim - problem getting the partner track. Continuing on to the next one");\r
+ continue;\r
+ }\r
+ AliExternalTrackParam *param2 = new AliExternalTrackParam(track2);\r
+ Int_t id1 = 0, id2 = 0;\r
+ AliESDv0 photonVtx(*param1,id1,*param2,id2);\r
+ Double_t vx,vy,vz;\r
+ photonVtx.GetXYZ(vx,vy,vz);\r
+\r
+ Double_t p1mom[3];\r
+ Double_t p2mom[3];\r
+ param1->PxPyPz(p1mom);\r
+ param2->PxPyPz(p2mom);\r
+\r
+ TVector3 p1momAtB(p1mom[0],p1mom[1],p1mom[2]);\r
+ TVector3 p2momAtB(p2mom[0],p2mom[1],p2mom[2]);\r
+ TVector3 sumMom = p1momAtB+p2momAtB;\r
+\r
+ Double_t ener1 = sqrt(pow(p1momAtB.Mag(),2) + massE*massE);\r
+ Double_t ener2 = sqrt(pow(p2momAtB.Mag(),2) + massE*massE);\r
+ Double_t mass = sqrt(pow((ener1+ener2),2) - pow(sumMom.Mag(),2));\r
+\r
+ Double_t dphi = p1momAtB.DeltaPhi(p2momAtB);\r
+ fh1OpeningAngle->Fill(dphi);\r
+ fh1MinvPhoton->Fill(mass);\r
+\r
+ if(mass < 0.1 ||\r
+ (mass > massEta-0.05 || mass < massEta+0.05) ||\r
+ (mass > massRho0-0.05 || mass < massRho0+0.05) ||\r
+ (mass > massOmega-0.05 || mass < massOmega+0.05) ||\r
+ (mass > massPhi-0.05 || mass < massPhi+0.05)) \r
+ {\r
+ \r
+ if(GetDebug() > 0) printf("######PROBABLY A PHOTON\n");\r
+ itIS = kTRUE;\r
+ }\r
+ \r
+ //clean up\r
+ delete param2;\r
+ \r
+ }\r
+\r
+ delete param1;\r
+ return itIS;\r
+\r
+}\r
+\r
+//__________________________________________________________________\r
+Bool_t AliAnaElectron::PhotonicV0(Int_t id) \r
+{\r
+ //This method checks to see whether a track that has been flagged as\r
+ //an electron was determined to match to a V0 candidate with\r
+ //invariant mass consistent with photon conversion\r
+\r
+ Bool_t itIS = kFALSE;\r
+\r
+ Double_t massEta = 0.547;\r
+ Double_t massRho0 = 0.770;\r
+ Double_t massOmega = 0.782;\r
+ Double_t massPhi = 1.020;\r
+ \r
+ //---Get V0s---\r
+ AliAODEvent *aod = (AliAODEvent*) GetReader()->GetInputEvent();\r
+ int nv0s = aod->GetNumberOfV0s();\r
+ for (Int_t iV0 = 0; iV0 < nv0s; iV0++) {\r
+ AliAODv0 *v0 = aod->GetV0(iV0);\r
+ if (!v0) continue;\r
+ double radius = v0->RadiusV0();\r
+ double mass = v0->InvMass2Prongs(0,1,11,11);\r
+ if(GetDebug() > 0) {\r
+ printf("## PhotonicV0() :: v0: %d, radius: %f \n", iV0 , radius );\r
+ printf("## PhotonicV0() :: neg-id: %d, pos-id: %d, THIS id: %d\n", v0->GetNegID(), v0->GetPosID(), id);\r
+ printf("## PhotonicV0() :: Minv(e,e): %f \n", v0->InvMass2Prongs(0,1,11,11) );\r
+ }\r
+ if (mass < 0.100 ||\r
+ (mass > massEta-0.05 || mass < massEta+0.05) ||\r
+ (mass > massRho0-0.05 || mass < massRho0+0.05) ||\r
+ (mass > massOmega-0.05 || mass < massOmega+0.05) ||\r
+ (mass > massPhi-0.05 || mass < massPhi+0.05)) {\r
+ if ( id == v0->GetNegID() || id == v0->GetPosID()) {\r
+ itIS=kTRUE;\r
+ if(GetDebug() > 0) printf("## PhotonicV0() :: It's a conversion electron!!! \n" );\r
+ }\r
+ } }\r
+ return itIS;\r
+\r
+}\r
+\r
+//__________________________________________________________________\r
+Bool_t AliAnaElectron::GetDCA(const AliAODTrack* track,Double_t impPar[2], Double_t cov[3]) \r
+{\r
+ //Use the Event vertex and AOD track information to get\r
+ //a real impact parameter for the track\r
+ //Once alice-off gets its act together and fixes the AOD, this\r
+ //should become obsolete.\r
+\r
+ Double_t maxD = 100000.; //max transverse IP\r
+ if(GetReader()->GetDataType() != AliCaloTrackReader::kMC) {\r
+ AliVEvent* ve = (AliVEvent*)GetReader()->GetInputEvent();\r
+ AliVVertex *vv = (AliVVertex*)ve->GetPrimaryVertex();\r
+ AliESDtrack esdTrack(track);\r
+ Double_t bfield[3];\r
+ esdTrack.GetBxByBz(bfield);\r
+ Bool_t gotit = esdTrack.PropagateToDCABxByBz(vv,bfield,maxD,impPar,cov);\r
+ printf("\t JLK impPar = %2.2f\n",impPar[0]);\r
+ return gotit;\r
+ }\r
+\r
+ return kFALSE;\r
+\r
+}\r
+\r
+//__________________________________________________________________\r
+Bool_t AliAnaElectron::CheckTrack(const AliAODTrack* track, const char* type) \r
+{\r
+ //Check this track to see if it is also tagged as an electron in the\r
+ //AliAODPWG4Particle list and if it is non-photonic\r
+\r
+ Bool_t pass = kFALSE;\r
+\r
+ Int_t trackId = track->GetID(); //get the index in the reader\r
+\r
+ Int_t naod = GetOutputAODBranch()->GetEntriesFast();\r
+ if(GetDebug() > 3) printf("AliAnaElectron::CheckTrack() - aod branch entries %d\n", naod);\r
+ for(Int_t iaod = 0; iaod < naod ; iaod++){\r
+ AliAODPWG4Particle* ele = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));\r
+ Int_t label = ele->GetTrackLabel(0);\r
+ if(label != trackId) continue; //skip to the next one if they don't match\r
+\r
+ if(strcmp(type,"DVM")==0) { \r
+ if(ele->CheckBTagBit(ele->GetBtag(),AliAODPWG4Particle::kDVMTag1) ||\r
+ ele->CheckBTagBit(ele->GetBtag(),AliAODPWG4Particle::kDVMTag2))\r
+ pass = kTRUE;\r
+\r
+ } else if (strcmp(type,"NPE")==0) {\r
+\r
+ Bool_t photonic = kFALSE;\r
+ Bool_t photonic1 = kFALSE;\r
+ photonic1 = PhotonicPrim(ele); //check against primaries\r
+ Bool_t photonic2 = kFALSE;\r
+ photonic2 = PhotonicV0(ele->GetTrackLabel(0)); //check against V0s\r
+ if(photonic1 || photonic2) photonic = kTRUE;\r
+ \r
+ if(!photonic) pass = kTRUE;\r
+\r
+ } else {\r
+ return kFALSE;\r
+ }\r
+ }\r
+\r
+ return pass;\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
+ //For determining how to classify electrons using MC info\r
+ //the number returned is the bin along one axis of 2-d histograms in\r
+ //which to fill this electron\r
+\r
+ //Do this first\r
+ if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)) return 4;\r
+\r
+ if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron)) {\r
+ //Bottom\r
+ if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEFromB)) return 1;\r
+ //Charm only\r
+ else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEFromC)\r
+ && !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEFromB)) return 2;\r
+ //Charm from bottom\r
+ else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEFromCFromB)) return 3;\r
+ // //Conversion\r
+ //else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)) return 4;\r
+ //Dalitz\r
+ else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay) \r
+ || GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) \r
+ || GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay)) return 5; \r
+ //W,Z\r
+ else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCWDecay)\r
+ || GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCZDecay)) return 6;\r
+ //Everything else\r
+ else \r
+ return 7;\r
+ } else {\r
+ //Misidentified electron\r
+ return 8;\r
+ }\r
+\r
+}\r
+\r
+//__________________________________________________________________\r
+Int_t AliAnaElectron::GetNumAODMCParticles() \r
+{\r
+ //Get the number of AliAODMCParticles, if any\r
+ Int_t num = 0;\r
+ Int_t npart0 = 0;\r
+ TClonesArray * mcparticles0 = 0x0;\r
+// TClonesArray * mcparticles1 = 0x0;\r
+\r
+ if(GetReader()->ReadAODMCParticles()){\r
+ //Get the list of MC particles\r
+ // \r
+ mcparticles0 = GetReader()->GetAODMCParticles(0);\r
+ if(!mcparticles0) {\r
+ if(GetDebug() > 0) 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
+ else{\r
+ npart0 = mcparticles0->GetEntriesFast();\r
+ }\r
+ //Int_t npart1 = 0;\r
+ //if(mcparticles1) npart1 = mcparticles1->GetEntriesFast();\r
+ //Int_t npart = npart0;//+npart1;\r
+ return npart0;\r
+\r
+ }\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) {\r
+ if (GetDebug() > 0) 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
+ else{\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
+ } else {\r
+ printf("AliAnaElectron::GetMCParticle() - Asked for AliAODMCParticle but we have a stack reader.\n");\r
+ }\r
+ return aodprimary;\r
+ \r
+}\r
+\r
+//__________________________________________________________________\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
+ AliStack* stack = 0x0;\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::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
+ 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
+ bjet=kTRUE;\r
+ //printf("MTH: **** found matching MC-Bjet: PDG=%d, pt=%f,dr=%f \n", pdgcode, primary->Pt(),dr );\r
+ break;\r
+ }\r
+ }\r
+ return bjet;\r
+\r
+}\r
+\r
+//__________________________________________________________________\r
+Bool_t AliAnaElectron::IsMcDJet(Double_t jeta, Double_t jphi)\r
+{\r
+ //Check if this jet is a charm jet\r
+ Bool_t cjet=kFALSE;\r
+\r
+ AliStack* stack = 0x0;\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
+ 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
+ if ( TMath::Abs(pdg) != 4) continue;\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
+\r
+ return cjet;\r
+\r
+}\r
+\r
+//__________________________________________________________________\r
+void AliAnaElectron::Print(const Option_t * opt) const\r
+{\r
+ //Print some relevant parameters set for the analysis\r
+ \r
+ if(! opt)\r
+ return;\r
+ \r
+ printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;\r
+ AliAnaPartCorrBaseClass::Print(" ");\r
+\r
+ printf("Calorimeter = %s\n", fCalorimeter.Data()) ;\r
+ printf("pOverE range = %f - %f\n",fpOverEmin,fpOverEmax);\r
+ printf("residual cut = %f\n",fResidualCut);\r
+ printf("---DVM Btagging\n");\r
+ printf("max IP-cut (e,h) = %f\n",fImpactCut);\r
+ printf("min ITS-hits = %d\n",fITSCut);\r
+ printf("max dR (e,h) = %f\n",fDrCut);\r
+ printf("max pairDCA = %f\n",fPairDcaCut);\r
+ printf("max decaylength = %f\n",fDecayLenCut);\r
+ printf("min Associated Pt = %f\n",fAssocPtCut);\r
+ printf("---IPSig Btagging\n");\r
+ printf("min tag track = %d\n",fNTagTrkCut);\r
+ printf("min IP significance = %f\n",fIPSigCut);\r
+ printf(" \n") ;\r
+ \r
+} \r
+\r
+//________________________________________________________________________\r
+void AliAnaElectron::ReadHistograms(TList* /* outputList */)\r
+{\r
+ // Needed when Terminate is executed in distributed environment \r
+ // Refill analysis histograms of this class with corresponding\r
+ // histograms in output list. \r
+\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
+\r
+ //Read histograms, must be in the same order as in\r
+ //GetCreateOutputObject. \r
+ //fh1pOverE = (TH1F *) outputList->At(index);\r
+ //fh1dR = (TH1F *) outputList->At(index++);\r
+ //fh2EledEdx = (TH2F *) outputList->At(index++);\r
+ //fh2MatchdEdx = (TH2F *) outputList->At(index++);\r
+ \r
+}\r
+\r
+//__________________________________________________________________\r
+void AliAnaElectron::Terminate(TList* outputList)\r
+{\r
+\r
+ //Do some plots to end\r
+ //Recover histograms from output histograms list, needed for\r
+ //distributed analysis. \r
+ //ReadHistograms(outputList);\r
+\r
+ printf(" AliAnaElectron::Terminate() *** %s Report: %d outputs\n", GetName(), outputList->GetEntries()) ;\r
+\r
+}\r
+\r
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff