JetAnaPythia< Jet > Class Template Reference

#include <JetAnaPythia.h>

Inheritance diagram for JetAnaPythia< Jet >:

Public Member Functions
	JetAnaPythia (edm::ParameterSet const &cfg)
Public Member Functions inherited from edm::EDAnalyzer
	EDAnalyzer ()
std::string	workerType () const
virtual	~EDAnalyzer ()

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

Additional Inherited Members
Public Types inherited from edm::EDAnalyzer
typedef EDAnalyzer	ModuleType
typedef WorkerT< EDAnalyzer >	WorkerType
Static Public Member Functions inherited from edm::EDAnalyzer
static const std::string &	baseType ()
static void	fillDescriptions (ConfigurationDescriptions &descriptions)
static void	prevalidate (ConfigurationDescriptions &)
Protected Member Functions inherited from edm::EDAnalyzer
CurrentProcessingContext const *	currentContext () const

Detailed Description

template<class Jet>
class JetAnaPythia< Jet >

Definition at line 16 of file JetAnaPythia.h.

Member Typedef Documentation

template<class Jet >

typedef std::vector<Jet> JetAnaPythia< Jet >::JetCollection

private

Definition at line 21 of file JetAnaPythia.h.

Constructor & Destructor Documentation

template<class Jet >

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");

}

Member Function Documentation

template<class Jet >

void JetAnaPythia< Jet >::analyze ( edm::Event const &  e, edm::EventSetup const &  iSetup  )

privatevirtual

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, fwrapper::jets, scaleCards::mass, 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

}

template<class Jet >

void JetAnaPythia< Jet >::beginJob ( void  )

privatevirtual

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");

  
}

template<class Jet >

void JetAnaPythia< Jet >::endJob ( void  )

privatevirtual

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;      
    }
}

template<class Jet >

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);
}

Member Data Documentation

template<class Jet >

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.

template<class Jet >

bool JetAnaPythia< Jet >::debug

private

Definition at line 48 of file JetAnaPythia.h.

template<class Jet >

float JetAnaPythia< Jet >::diJetMass

private

Definition at line 37 of file JetAnaPythia.h.

template<class Jet >

float JetAnaPythia< Jet >::diPartMass

private

Definition at line 38 of file JetAnaPythia.h.

template<class Jet >

float JetAnaPythia< Jet >::etaJet1

private

Definition at line 33 of file JetAnaPythia.h.

template<class Jet >

float JetAnaPythia< Jet >::etaJet2

private

Definition at line 33 of file JetAnaPythia.h.

template<class Jet >

float JetAnaPythia< Jet >::etaPart1

private

Definition at line 35 of file JetAnaPythia.h.

template<class Jet >

float JetAnaPythia< Jet >::etaPart2

private

Definition at line 35 of file JetAnaPythia.h.

template<class Jet >

int JetAnaPythia< Jet >::eventsGen

private

Definition at line 50 of file JetAnaPythia.h.

template<class Jet >

std::string JetAnaPythia< Jet >::HistoFileName

private

Definition at line 44 of file JetAnaPythia.h.

template<class Jet >

std::string JetAnaPythia< Jet >::JetAlgorithm

private

Definition at line 42 of file JetAnaPythia.h.

template<class Jet >

TFile* JetAnaPythia< Jet >::m_file

private

Definition at line 39 of file JetAnaPythia.h.

template<class Jet >

std::map<TString, TH1*> JetAnaPythia< Jet >::m_HistNames1D

private

Definition at line 26 of file JetAnaPythia.h.

template<class Jet >

TTree* JetAnaPythia< Jet >::mcTruthTree_

private

Definition at line 28 of file JetAnaPythia.h.

template<class Jet >

int JetAnaPythia< Jet >::nJets

private

Definition at line 32 of file JetAnaPythia.h.

template<class Jet >

int JetAnaPythia< Jet >::NJets

private

Definition at line 46 of file JetAnaPythia.h.

template<class Jet >

float JetAnaPythia< Jet >::pt_hat

private

Definition at line 31 of file JetAnaPythia.h.

template<class Jet >

std::vector<double> JetAnaPythia< Jet >::ptHatEdges

private

Definition at line 61 of file JetAnaPythia.h.

template<class Jet >

float JetAnaPythia< Jet >::ptJet1

private

Definition at line 34 of file JetAnaPythia.h.

template<class Jet >

float JetAnaPythia< Jet >::ptJet2

private

Definition at line 34 of file JetAnaPythia.h.

template<class Jet >

float JetAnaPythia< Jet >::ptPart1

private

Definition at line 36 of file JetAnaPythia.h.

template<class Jet >

float JetAnaPythia< Jet >::ptPart2

private

Definition at line 36 of file JetAnaPythia.h.

template<class Jet >

float JetAnaPythia< Jet >::weight

private

Definition at line 30 of file JetAnaPythia.h.

Referenced by dirstructure.Directory::calcStats(), dirstructure.Directory::get_fail_rate(), dirstructure.Directory::get_null_rate(), dirstructure.Directory::get_skiped_rate(), dirstructure.Directory::get_success_rate(), and dirstructure.Directory::print_report().

template<class Jet >

float JetAnaPythia< Jet >::xsec

private

Definition at line 29 of file JetAnaPythia.h.

template<class Jet >

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.