|
|
|
#include <JetAnaPythia.h>
Public Member Functions | |
| JetAnaPythia (edm::ParameterSet const &cfg) | |
Private Types | |
| typedef std::vector< Jet > | JetCollection |
Private Member Functions | |
| void | analyze (edm::Event const &e, edm::EventSetup const &iSetup) |
| void | beginJob () |
| void | endJob () |
| void | FillHist1D (const TString &histName, const Double_t &x, const Double_t &wt) |
Private Attributes | |
| std::string | anaLevel |
| bool | debug |
| float | diJetMass |
| float | diPartMass |
| float | etaJet1 |
| float | etaJet2 |
| float | etaPart1 |
| float | etaPart2 |
| int | eventsGen |
| std::string | HistoFileName |
| std::string | JetAlgorithm |
| TFile * | m_file |
| std::map< TString, TH1 * > | m_HistNames1D |
| TTree * | mcTruthTree_ |
| int | NJets |
| int | nJets |
| float | pt_hat |
| std::vector< double > | ptHatEdges |
| float | ptJet1 |
| float | ptJet2 |
| float | ptPart1 |
| float | ptPart2 |
| float | weight |
| float | xsec |
| std::vector< double > | xsecGen |
Definition at line 16 of file JetAnaPythia.h.
typedef std::vector<Jet> JetAnaPythia< Jet >::JetCollection [private] |
Definition at line 21 of file JetAnaPythia.h.
| JetAnaPythia< Jet >::JetAnaPythia | ( | edm::ParameterSet const & | cfg | ) |
Definition at line 27 of file JetAnaPythia.cc.
References debug, and edm::ParameterSet::getParameter().
{
JetAlgorithm = cfg.getParameter<std::string> ("JetAlgorithm");
HistoFileName = cfg.getParameter<std::string> ("HistoFileName");
NJets = cfg.getParameter<int> ("NJets");
debug = cfg.getParameter<bool> ("debug");
eventsGen = cfg.getParameter<int> ("eventsGen");
anaLevel = cfg.getParameter<std::string> ("anaLevel");
xsecGen = cfg.getParameter< vector<double> > ("xsecGen");
ptHatEdges = cfg.getParameter< vector<double> > ("ptHatEdges");
}
| void JetAnaPythia< Jet >::analyze | ( | edm::Event const & | e, |
| edm::EventSetup const & | iSetup | ||
| ) | [private, virtual] |
evt.getRun().getByLabel("generator", genInfoProduct );
Histograms for Dijet Mass Analysis ////
Histograms for Dijet Ratio Analysis: Inner region ///
Histograms for Dijet Ratio Analysis: Outer region ////
We are looking at dijet resonances. ////
We are looking at QCD ////
Diparton mass for dijet mass analysis ////
Diparton mass for dijet ratio analysis: inner region ///
Diparton mass for dijet ratio analysis: outer region ///
Implements edm::EDAnalyzer.
Definition at line 140 of file JetAnaPythia.cc.
References abs, gather_cfg::cout, debug, edm::Event::getByLabel(), i, getHLTprescales::index, analyzePatCleaning_cfg::jets, AlCaHLTBitMon_ParallelJobs::p, reco::LeafCandidate::p4(), reco::LeafCandidate::pdgId(), lumiQueryAPI::q, alignCSCRings::r, reco::LeafCandidate::status(), and CommonMethods::weight().
{
int notDone=1;
while(notDone)
{ //while loop to allow us to tailor the analysis level for faster running.
TString hname;
// Process Info
//edm::Handle< double > genEventScale;
//evt.getByLabel("genEventScale", genEventScale );
//pt_hat = *genEventScale;
edm::Handle<edm::HepMCProduct> MCevt;
evt.getByLabel("generator", MCevt);
HepMC::GenEvent * myGenEvent = new HepMC::GenEvent(*(MCevt->GetEvent()));
double pthat = myGenEvent->event_scale();
pt_hat = float(pthat);
delete myGenEvent;
if(anaLevel != "generating"){ //We are not generating events, so xsec is there
//edm::Handle< GenRunInfoProduct > genInfoProduct;
//xsec = (double)genInfoProduct->externalXSecLO();
xsec=0.0;
if( ptHatEdges.size()>xsecGen.size() ){
for(unsigned int i_pthat = 0; i_pthat < xsecGen.size(); ++i_pthat){
if( pthat >= ptHatEdges[i_pthat] && pthat < ptHatEdges[i_pthat+1])xsec=float(xsecGen[i_pthat]);
}
}
else
{
std::cout << "Number of PtHat bin edges too small. Xsec set to zero" << std::endl;
}
}
else
{
xsec = xsecGen[0]; //Generating events, no xsec in event, get xsec from user input
}
if(debug)std::cout << "cross section=" <<xsec << " pb" << std::endl;
weight = xsec/eventsGen;
if(debug)std::cout << "pt_hat=" <<pt_hat << std::endl;
hname = "PtHat";
FillHist1D(hname, pt_hat, 1.0);
hname = "PtHatFine";
FillHist1D(hname, pt_hat, 1.0);
hname = "PtHatWt";
FillHist1D(hname, pt_hat, weight);
hname = "PtHatFineWt";
FillHist1D(hname, pt_hat, weight);
if(anaLevel=="PtHatOnly")break; //ptHatOnly should be very fast
// Jet Info
math::XYZTLorentzVector p4jet[2];
float etajet[2];
Handle<JetCollection> jets;
evt.getByLabel(JetAlgorithm,jets);
typename JetCollection::const_iterator i_jet;
int index = 0;
hname = "NumberOfJets";
nJets = jets->size();
FillHist1D(hname,nJets,1.0);
// Two Leading Jet Info
for(i_jet = jets->begin(); i_jet != jets->end() && index < 2; ++i_jet)
{
hname = "JetPt";
FillHist1D(hname,i_jet->pt(),1.0);
hname = "JetEta";
FillHist1D(hname,i_jet->eta(),1.0);
hname = "JetPhi";
FillHist1D(hname,i_jet->phi(),1.0);
p4jet[index] = i_jet->p4();
etajet[index] = i_jet->eta();
if(debug)std::cout << "jet " << index+1 <<": pt=" <<i_jet->pt() << ", eta=" <<etajet[index] << std::endl;
index++;
}
// TTree variables //
etaJet1 = etajet[0];
etaJet2 = etajet[1];
ptJet1 = p4jet[0].pt();
ptJet2 = p4jet[1].pt();
diJetMass = (p4jet[0]+p4jet[1]).mass();
if(index==2&&abs(etaJet1)<1.3&&abs(etaJet2)<1.3){
hname = "DijetMass";
FillHist1D(hname,diJetMass ,1.0);
hname = "DijetMassWt";
FillHist1D(hname,diJetMass ,weight);
}
if(index==2&&abs(etaJet1)<0.7&&abs(etaJet2)<0.7){
hname = "DijetMassIn";
FillHist1D(hname,diJetMass ,1.0);
hname = "DijetMassInWt";
FillHist1D(hname,diJetMass ,weight);
}
if(index==2 && (abs(etaJet1)>0.7&&abs(etaJet1)<1.3)
&& (abs(etaJet2)>0.7&&abs(etaJet2)<1.3) ){
hname = "DijetMassOut";
FillHist1D(hname, diJetMass ,1.0);
hname = "DijetMassOutWt";
FillHist1D(hname,diJetMass ,weight);
}
if(anaLevel=="Jets")break; //Jets level for samples without genParticles
// Parton Info
edm::Handle<std::vector<reco::GenParticle> > genParticlesHandle_;
evt.getByLabel("genParticles",genParticlesHandle_);
if(debug)for( size_t i = 0; i < genParticlesHandle_->size(); ++ i ) {
const reco::GenParticle & p = (*genParticlesHandle_)[i];
int id = p.pdgId();
int st = p.status();
math::XYZTLorentzVector genP4 = p.p4();
if(i>=2&&i<=8)std::cout << "particle " << i << ": id=" << id << ", status=" << st << ", mass=" << genP4.mass() << ", pt=" << genP4.pt() << ", eta=" << genP4.eta() << std::endl;
}
// Examine the 7th particle in pythia.
// It should be either a resonance (abs(id)>=32) or the first outgoing parton
// for the processes we will consider: dijet resonances, QCD, or QCD +contact interactions.
const reco::GenParticle & p = (*genParticlesHandle_)[6];
int id = p.pdgId();
math::XYZTLorentzVector resonance_p, parton1_p, parton2_p;
if(abs(id)>=32){
resonance_p = p.p4();
hname = "ResonanceMass";
FillHist1D(hname,resonance_p.mass() ,1.0);
const reco::GenParticle & q = (*genParticlesHandle_)[7];
parton1_p = q.p4();
const reco::GenParticle & r = (*genParticlesHandle_)[8];
parton2_p = r.p4();
if(debug)std::cout << "Resonance mass=" << resonance_p.mass() << ", parton 1 pt=" << parton1_p.pt() << ", parton 2 pt=" << parton2_p.pt() << ", diparton mass=" << (parton1_p+parton2_p).mass() << std::endl;
}
else
{
parton1_p = p.p4();
const reco::GenParticle & q = (*genParticlesHandle_)[7];
parton2_p = q.p4();
if(debug)std::cout << "parton 1 pt=" << parton1_p.pt() << ", parton 2 pt=" << parton2_p.pt() << ", diparton mass=" << (parton1_p+parton2_p).mass() << std::endl;
}
etaPart1 = parton1_p.eta();
etaPart2 = parton2_p.eta();
ptPart1 = parton1_p.pt();
ptPart2 = parton2_p.pt();
diPartMass = (parton1_p+parton2_p).mass();
if(abs(etaPart1)<1.3&&abs(etaPart2)<1.3){
hname = "DipartonMass";
FillHist1D(hname,diPartMass ,1.0);
hname = "DipartonMassWt";
FillHist1D(hname,diPartMass ,weight);
}
if(abs(etaPart1)<0.7&&abs(etaPart2)<0.7){
hname = "DipartonMassIn";
FillHist1D(hname,diPartMass ,1.0);
hname = "DipartonMassInWt";
FillHist1D(hname,diPartMass ,weight);
}
if( (abs(etaPart1)>0.7&&abs(etaPart1)<1.3)
&& (abs(etaPart2)>0.7&&abs(etaPart2)<1.3) ){
hname = "DipartonMassOut";
FillHist1D(hname,diPartMass ,1.0);
hname = "DipartonMassOutWt";
FillHist1D(hname,diPartMass ,weight);
}
// Fill the TTree //
mcTruthTree_->Fill();
notDone=0; //We are done, exit the while loop
}//end of while
}
| void JetAnaPythia< Jet >::beginJob | ( | void | ) | [private, virtual] |
Reimplemented from edm::EDAnalyzer.
Definition at line 41 of file JetAnaPythia.cc.
References M_PI, and CommonMethods::weight().
{
TString hname;
m_file = new TFile(HistoFileName.c_str(),"RECREATE");
const int nMassBins = 103;
double massBoundaries[nMassBins+1] = {1, 3, 6, 10, 16, 23, 31, 40, 50, 61, 74, 88, 103, 119, 137, 156, 176, 197, 220, 244, 270, 296, 325, 354, 386, 419, 453, 489, 526, 565, 606, 649, 693, 740, 788, 838, 890, 944, 1000, 1058, 1118, 1181, 1246, 1313, 1383, 1455, 1530, 1607, 1687, 1770, 1856, 1945, 2037, 2132, 2231, 2332, 2438, 2546, 2659, 2775, 2895, 3019, 3147, 3279, 3416, 3558, 3704, 3854, 4010, 4171, 4337, 4509, 4686, 4869, 5058, 5253, 5455, 5663, 5877, 6099, 6328, 6564, 6808, 7060, 7320, 7589, 7866, 8152, 8447, 8752, 9067, 9391, 9726, 10072, 10430, 10798, 11179, 11571, 11977, 12395, 12827, 13272, 13732, 14000};
hname = "JetPt";
m_HistNames1D[hname] = new TH1F(hname,hname,500,0,5000);
hname = "JetEta";
m_HistNames1D[hname] = new TH1F(hname,hname,120,-6,6);
hname = "JetPhi";
m_HistNames1D[hname] = new TH1F(hname,hname,100,-M_PI,M_PI);
hname = "NumberOfJets";
m_HistNames1D[hname] = new TH1F(hname,hname,100,0,100);
hname = "DijetMass";
m_HistNames1D[hname] = new TH1F(hname,hname,nMassBins,massBoundaries);
hname = "DijetMassWt";
m_HistNames1D[hname] = new TH1F(hname,hname,nMassBins,massBoundaries);
m_HistNames1D.find(hname)->second->Sumw2();
hname = "DijetMassIn";
m_HistNames1D[hname] = new TH1F(hname,hname,nMassBins,massBoundaries);
hname = "DijetMassInWt";
m_HistNames1D[hname] = new TH1F(hname,hname,nMassBins,massBoundaries);
m_HistNames1D.find(hname)->second->Sumw2();
hname = "DijetMassOut";
m_HistNames1D[hname] = new TH1F(hname,hname,nMassBins,massBoundaries);
hname = "DijetMassOutWt";
m_HistNames1D[hname] = new TH1F(hname,hname,nMassBins,massBoundaries);
m_HistNames1D.find(hname)->second->Sumw2();
hname = "ResonanceMass";
m_HistNames1D[hname] = new TH1F(hname,hname,nMassBins,massBoundaries);
hname = "DipartonMass";
m_HistNames1D[hname] = new TH1F(hname,hname,nMassBins,massBoundaries);
hname = "DipartonMassWt";
m_HistNames1D[hname] = new TH1F(hname,hname,nMassBins,massBoundaries);
m_HistNames1D.find(hname)->second->Sumw2();
hname = "DipartonMassIn";
m_HistNames1D[hname] = new TH1F(hname,hname,nMassBins,massBoundaries);
hname = "DipartonMassInWt";
m_HistNames1D[hname] = new TH1F(hname,hname,nMassBins,massBoundaries);
m_HistNames1D.find(hname)->second->Sumw2();
hname = "DipartonMassOut";
m_HistNames1D[hname] = new TH1F(hname,hname,nMassBins,massBoundaries);
hname = "DipartonMassOutWt";
m_HistNames1D[hname] = new TH1F(hname,hname,nMassBins,massBoundaries);
m_HistNames1D.find(hname)->second->Sumw2();
hname = "PtHat";
m_HistNames1D[hname] = new TH1F(hname,hname,1000,0,5000);
hname = "PtHatWt";
m_HistNames1D[hname] = new TH1F(hname,hname,1000,0,5000);
m_HistNames1D.find(hname)->second->Sumw2();
hname = "PtHatFine";
m_HistNames1D[hname] = new TH1F(hname,hname,5000,0,5000);
hname = "PtHatFineWt";
m_HistNames1D[hname] = new TH1F(hname,hname,5000,0,5000);
m_HistNames1D.find(hname)->second->Sumw2();
mcTruthTree_ = new TTree("mcTruthTree","mcTruthTree");
mcTruthTree_->Branch("xsec", &xsec, "xsec/F");
mcTruthTree_->Branch("weight", &weight, "weight/F");
mcTruthTree_->Branch("pt_hat", &pt_hat, "pt_hat/F");
mcTruthTree_->Branch("nJets", &nJets, "nJets/I");
mcTruthTree_->Branch("etaJet1", &etaJet1, "etaJet1/F");
mcTruthTree_->Branch("etaJet2", &etaJet2, "etaJet2/F");
mcTruthTree_->Branch("ptJet1", &ptJet1, "ptJet1/F");
mcTruthTree_->Branch("ptJet2", &ptJet2, "ptJet2/F");
mcTruthTree_->Branch("diJetMass", &diJetMass, "diJetMass/F");
mcTruthTree_->Branch("etaPart1", &etaPart1, "etaPart1/F");
mcTruthTree_->Branch("etaPart2", &etaPart2, "etaPart2/F");
mcTruthTree_->Branch("ptPart1", &ptPart1, "ptPart1/F");
mcTruthTree_->Branch("ptPart2", &ptPart2, "ptPart2/F");
mcTruthTree_->Branch("diPartMass", &diPartMass, "diPartMass/F");
}
| void JetAnaPythia< Jet >::endJob | ( | void | ) | [private, virtual] |
Reimplemented from edm::EDAnalyzer.
Definition at line 331 of file JetAnaPythia.cc.
{
if (m_file !=0)
{
m_file->cd();
mcTruthTree_->Write();
for (std::map<TString, TH1*>::iterator hid = m_HistNames1D.begin(); hid != m_HistNames1D.end(); hid++)
hid->second->Write();
delete m_file;
m_file = 0;
}
}
| void JetAnaPythia< Jet >::FillHist1D | ( | const TString & | histName, |
| const Double_t & | x, | ||
| const Double_t & | wt | ||
| ) | [private] |
Definition at line 346 of file JetAnaPythia.cc.
References gather_cfg::cout.
{
std::map<TString, TH1*>::iterator hid=m_HistNames1D.find(histName);
if (hid==m_HistNames1D.end())
std::cout << "%fillHist -- Could not find histogram with name: " << histName << std::endl;
else
hid->second->Fill(value,wt);
}
std::string JetAnaPythia< Jet >::anaLevel [private] |
Analysis level string. Can speed up job by looking at less /// PtHatOnly: only get PtHat and make PtHat histos Jets: do histogram analysis of jets, but not partons all: do analysis of everything and make histos and root tree generating: analysis of everything, make histos and root tree
Definition at line 56 of file JetAnaPythia.h.
bool JetAnaPythia< Jet >::debug [private] |
Definition at line 48 of file JetAnaPythia.h.
float JetAnaPythia< Jet >::diJetMass [private] |
Definition at line 37 of file JetAnaPythia.h.
float JetAnaPythia< Jet >::diPartMass [private] |
Definition at line 38 of file JetAnaPythia.h.
float JetAnaPythia< Jet >::etaJet1 [private] |
Definition at line 33 of file JetAnaPythia.h.
float JetAnaPythia< Jet >::etaJet2 [private] |
Definition at line 33 of file JetAnaPythia.h.
float JetAnaPythia< Jet >::etaPart1 [private] |
Definition at line 35 of file JetAnaPythia.h.
float JetAnaPythia< Jet >::etaPart2 [private] |
Definition at line 35 of file JetAnaPythia.h.
int JetAnaPythia< Jet >::eventsGen [private] |
Definition at line 50 of file JetAnaPythia.h.
std::string JetAnaPythia< Jet >::HistoFileName [private] |
Definition at line 44 of file JetAnaPythia.h.
std::string JetAnaPythia< Jet >::JetAlgorithm [private] |
Definition at line 42 of file JetAnaPythia.h.
TFile* JetAnaPythia< Jet >::m_file [private] |
Definition at line 39 of file JetAnaPythia.h.
std::map<TString, TH1*> JetAnaPythia< Jet >::m_HistNames1D [private] |
Definition at line 26 of file JetAnaPythia.h.
TTree* JetAnaPythia< Jet >::mcTruthTree_ [private] |
Definition at line 28 of file JetAnaPythia.h.
int JetAnaPythia< Jet >::NJets [private] |
Definition at line 46 of file JetAnaPythia.h.
int JetAnaPythia< Jet >::nJets [private] |
Definition at line 32 of file JetAnaPythia.h.
float JetAnaPythia< Jet >::pt_hat [private] |
Definition at line 31 of file JetAnaPythia.h.
std::vector<double> JetAnaPythia< Jet >::ptHatEdges [private] |
Definition at line 61 of file JetAnaPythia.h.
float JetAnaPythia< Jet >::ptJet1 [private] |
Definition at line 34 of file JetAnaPythia.h.
float JetAnaPythia< Jet >::ptJet2 [private] |
Definition at line 34 of file JetAnaPythia.h.
float JetAnaPythia< Jet >::ptPart1 [private] |
Definition at line 36 of file JetAnaPythia.h.
float JetAnaPythia< Jet >::ptPart2 [private] |
Definition at line 36 of file JetAnaPythia.h.
float JetAnaPythia< Jet >::weight [private] |
Definition at line 30 of file JetAnaPythia.h.
float JetAnaPythia< Jet >::xsec [private] |
Definition at line 29 of file JetAnaPythia.h.
std::vector<double> JetAnaPythia< Jet >::xsecGen [private] |
Generator cross section Only 1 entry in case analysis level is "generating" //// Multiple entries when analyzing ///
Definition at line 60 of file JetAnaPythia.h.
1.7.3