diff --git a/Validation/Generator/plugins/BBbarAnalyzer.cc b/Validation/Generator/plugins/BBbarAnalyzer.cc new file mode 100644 index 0000000000000..51419e2b81bcd --- /dev/null +++ b/Validation/Generator/plugins/BBbarAnalyzer.cc @@ -0,0 +1,107 @@ +/* This is en example for an Analyzer of a Herwig HepMCProduct + and looks for muons pairs and fills a histogram + with the invaraint mass of the four. +*/ + +// +// Original Author: Fabian Stoeckli +// Created: Tue Nov 14 13:43:02 CET 2006 +// $Id: BBbarAnalyzer.cc,v 1.3 2009/07/31 23:21:07 kharchil Exp $ +// +// + + +// system include files +#include +#include + +// user include files +#include "Validation/Generator/plugins/BBbarAnalyzer.h" +#include "DataFormats/Math/interface/LorentzVector.h" +#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h" + + +BBbarAnalyzer::BBbarAnalyzer(const edm::ParameterSet& iConfig) +{ + outputFilename=iConfig.getUntrackedParameter("OutputFilename","dummy.root"); + Pt_histo = new TH1F("invmass_histo","invmass_histo",100,0,20); + invmass_histo = new TH1F("mu_invmass_histo","mu_invmass_histo",100,0,20); +} + +BBbarAnalyzer::~BBbarAnalyzer() +{ + +} + +void +BBbarAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) +{ + using namespace edm; + + // get HepMC::GenEvent ... + //Handle evt_h; + //iEvent.getByType(evt_h); + //HepMC::GenEvent * evt = new HepMC::GenEvent(*(evt_h->GetEvent())); + + + // look for stable muons + // std::vector muons; + //muons.resize(0); + //for(HepMC::GenEvent::particle_iterator it = evt->particles_begin(); it != evt->particles_end(); ++it) { + // if(abs((*it)->pdg_id())==13 && (*it)->status()==1) { + // muons.push_back(*it); + // } + // } + + math::XYZTLorentzVector Lvector1(0,0,0,0); + math::XYZTLorentzVector Lvector2(0,0,0,0); + + Handle mcEventHandle; + try{ + iEvent.getByLabel("source",mcEventHandle); + }catch(...) {;} + + if(mcEventHandle.isValid()){ + const HepMC::GenEvent* mcEvent = mcEventHandle->GetEvent() ; + + // if there are at least four muons + // calculate invarant mass of first two and fill it into histogram + math::XYZTLorentzVector tot_momentum; math::XYZTLorentzVector tot_mumomentum; + float inv_mass = 0.0; + //double mu_invmass = 0.0; + float Pt = 0; + + HepMC::GenEvent::particle_const_iterator i; + HepMC::GenEvent::particle_const_iterator j; + for(i = mcEvent->particles_begin(); i!= mcEvent->particles_end(); i++){ + for(j = mcEvent->particles_begin(); j!= mcEvent->particles_end(); j++){ + HepMC::FourVector p41 = (*i)->momentum(); + HepMC::FourVector p42 = (*i)->momentum(); + Lvector1 = math::XYZTLorentzVector(p41.px(),p41.py(),p41.pz(),p41.e()); + Lvector2 = math::XYZTLorentzVector(p42.px(),p42.py(),p42.pz(),p42.e()); + Pt = sqrt(p41.px()*p41.px()+p41.py()*p41.py()); + Pt_histo->Fill(Pt); + inv_mass = sqrt((p41.e()+p42.e())*(p41.e()+p42.e())-((p41.px()+p42.px())*(p41.px()+p42.px())+(p41.py()+p42.py())*(p41.py()+p42.py())+(p41.pz()+p42.pz())*(p41.pz()+p42.pz()))); + invmass_histo->Fill(inv_mass); + } + } + } +} + +void +BBbarAnalyzer::beginJob(const edm::EventSetup&) +{ +} + +void +BBbarAnalyzer::endJob() { + // save histograms into file + TFile file(outputFilename.c_str(),"RECREATE"); + Pt_histo->Write(); + invmass_histo->Write(); + file.Close(); + +} + +#include "FWCore/Framework/interface/MakerMacros.h" +DEFINE_FWK_MODULE(BBbarAnalyzer); diff --git a/Validation/Generator/plugins/BBbarAnalyzer.h b/Validation/Generator/plugins/BBbarAnalyzer.h new file mode 100644 index 0000000000000..cc87a24ff8d1d --- /dev/null +++ b/Validation/Generator/plugins/BBbarAnalyzer.h @@ -0,0 +1,29 @@ +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/EDAnalyzer.h" +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" + +#include "HepMC/GenEvent.h" +#include "HepMC/GenParticle.h" +#include "HepMC/WeightContainer.h" + +#include "TH1D.h" +#include "TFile.h" + + +class BBbarAnalyzer : public edm::EDAnalyzer { + public: + explicit BBbarAnalyzer(const edm::ParameterSet&); + ~BBbarAnalyzer(); + + private: + virtual void beginJob(const edm::EventSetup&) ; + virtual void analyze(const edm::Event&, const edm::EventSetup&); + virtual void endJob() ; + + private: + std::string outputFilename; + TH1D* weight_histo; + TH1F* invmass_histo; + TH1F* Pt_histo; +}; diff --git a/Validation/Generator/plugins/H4muAnalyzer.cc b/Validation/Generator/plugins/H4muAnalyzer.cc new file mode 100644 index 0000000000000..f33ff91d142ec --- /dev/null +++ b/Validation/Generator/plugins/H4muAnalyzer.cc @@ -0,0 +1,100 @@ +/* This is en example for an Analyzer of a Herwig HepMCProduct + and looks for muons pairs and fills a histogram + with the invaraint mass of the four. +*/ + +// +// Original Author: Fabian Stoeckli +// Created: Tue Nov 14 13:43:02 CET 2006 +// $Id: H4muAnalyzer.cc,v 1.3 2009/07/31 23:21:08 kharchil Exp $ +// +// + +// system include files +#include +#include +// user include files +#include "Validation/Generator/plugins/H4muAnalyzer.h" + +#include "DataFormats/Math/interface/LorentzVector.h" +#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h" + + +H4muAnalyzer::H4muAnalyzer(const edm::ParameterSet& iConfig) +{ + outputFilename=iConfig.getUntrackedParameter("OutputFilename","dummy.root"); + invmass_histo = new TH1D("invmass_histo","invmass_histo",60,180,200); + Z_invmass_histo = new TH1D("mu_invmass_histo","mu_invmass_histo",60,70,100); +} + +H4muAnalyzer::~H4muAnalyzer() +{ + +} + +void +H4muAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) +{ + using namespace edm; + + // get HepMC::GenEvent ... + Handle evt_h; + iEvent.getByType(evt_h); + HepMC::GenEvent * evt = new HepMC::GenEvent(*(evt_h->GetEvent())); + + + // look for stable muons + std::vector muons; + muons.resize(0); + for(HepMC::GenEvent::particle_iterator it = evt->particles_begin(); it != evt->particles_end(); ++it) { + if(abs((*it)->pdg_id())==13 && (*it)->status()==1) { + muons.push_back(*it); + } + } + + // if there are at least four muons + // calculate invarant mass of first two and fill it into histogram + math::XYZTLorentzVector tot_momentum; math::XYZTLorentzVector tot_mumomentum; + double inv_mass = 0.0; double mu_invmass = 0.0; + if(muons.size()>1) { + for(unsigned int i=0; imomentum()); + math::XYZTLorentzVector mumom1(muons[j]->momentum()); + tot_mumomentum = mumom + mumom1; + mu_invmass = sqrt(tot_mumomentum.mass2()); + Z_invmass_histo->Fill(mu_invmass); + } + } + //invmass_histo->Fill(inv_mass); + +} + + if(muons.size()>3) { + for(unsigned int i=0; i<4; ++i) { + math::XYZTLorentzVector mom(muons[i]->momentum()); + tot_momentum += mom; + } + inv_mass = sqrt(tot_momentum.mass2()); + } + invmass_histo->Fill(inv_mass); + +} + +void +H4muAnalyzer::beginJob(const edm::EventSetup&) +{ +} + +void +H4muAnalyzer::endJob() { + // save histograms into file + TFile file(outputFilename.c_str(),"RECREATE"); + invmass_histo->Write(); + Z_invmass_histo->Write(); + file.Close(); + +} + +#include "FWCore/Framework/interface/MakerMacros.h" +DEFINE_FWK_MODULE(H4muAnalyzer); diff --git a/Validation/Generator/plugins/H4muAnalyzer.h b/Validation/Generator/plugins/H4muAnalyzer.h new file mode 100644 index 0000000000000..22f9ffd917174 --- /dev/null +++ b/Validation/Generator/plugins/H4muAnalyzer.h @@ -0,0 +1,30 @@ +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/EDAnalyzer.h" +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" + +#include "HepMC/GenEvent.h" +#include "HepMC/GenParticle.h" +#include "HepMC/WeightContainer.h" + +#include "TH1D.h" +#include "TFile.h" + + +class H4muAnalyzer : public edm::EDAnalyzer { + public: + explicit H4muAnalyzer(const edm::ParameterSet&); + ~H4muAnalyzer(); + + private: + virtual void beginJob(const edm::EventSetup&) ; + virtual void analyze(const edm::Event&, const edm::EventSetup&); + virtual void endJob() ; + + private: + std::string outputFilename; + TH1D* weight_histo; + TH1D* invmass_histo; + TH1D* Z_invmass_histo; + +}; diff --git a/Validation/Generator/plugins/HZZ4muAnalyzer.cc b/Validation/Generator/plugins/HZZ4muAnalyzer.cc new file mode 100755 index 0000000000000..01b447c1c8dac --- /dev/null +++ b/Validation/Generator/plugins/HZZ4muAnalyzer.cc @@ -0,0 +1,217 @@ +#include + +#include "Validation/Generator/plugins/HZZ4muAnalyzer.h" +//#include "SimDataFormats/HepMCProduct/interface/HepMCProduct.h" +#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h" +// essentials !!! +#include "FWCore/Framework/interface/Event.h" +#include "DataFormats/Common/interface/Handle.h" + +#include "TFile.h" +#include "TH1.h" + +using namespace edm; +using namespace std; + +HZZ4muAnalyzer::HZZ4muAnalyzer( const ParameterSet& pset ) + : fOutputFileName( pset.getUntrackedParameter("HistOutFile",std::string("TestHiggsMass.root")) ), + fOutputFile(0), fHist2muMass(0), fHist4muMass(0), fHistZZMass(0) +{ +} + +void HZZ4muAnalyzer::beginJob( const EventSetup& ) +{ + + fOutputFile = new TFile( fOutputFileName.c_str(), "RECREATE" ) ; + fHist2muMass = new TH1D( "Hist2muMass", "2-mu inv. mass", 100, 60., 120. ) ; + fHist4muMass = new TH1D( "Hist4muMass", "4-mu inv. mass", 100, 170., 210. ) ; + fHistZZMass = new TH1D( "HistZZMass", "ZZ inv. mass", 100, 170., 210. ) ; + + return ; +} + +void HZZ4muAnalyzer::analyze( const Event& e, const EventSetup& ) +{ + + Handle< HepMCProduct > EvtHandle ; + + // find initial (unsmeared, unfiltered,...) HepMCProduct + // + e.getByLabel( "source", EvtHandle ) ; + + const HepMC::GenEvent* Evt = EvtHandle->GetEvent() ; + + // this a pointer - and not an array/vector/... + // because this example explicitely assumes + // that there one and only Higgs in the record + // + HepMC::GenVertex* HiggsDecVtx = 0 ; + + // find the 1st vertex with outgoing Higgs + // and get Higgs decay vertex from there; + // + // in principal, one can look for the vertex + // with incoming Higgs as well... + // + for ( HepMC::GenEvent::vertex_const_iterator + vit=Evt->vertices_begin(); vit!=Evt->vertices_end(); vit++ ) + { + for ( HepMC::GenVertex::particles_out_const_iterator + pout=(*vit)->particles_out_const_begin(); + pout!=(*vit)->particles_out_const_end(); pout++ ) + { + if ( (*pout)->pdg_id() == 25 && (*pout)->status() == 2 ) + { + if ( (*pout)->end_vertex() != 0 ) + { + HiggsDecVtx = (*pout)->end_vertex() ; + break ; + } + } + } + if ( HiggsDecVtx != 0 ) + { + break ; // break the initial loop over vertices + } + } + + if ( HiggsDecVtx == 0 ) + { + cout << " There is NO Higgs in this event ! " << endl ; + return ; + } + + if ( e.id().event() == 1 ) + { + cout << " " << endl ; + cout << " We do some example printouts in the event 1 " << endl ; + cout << " Higgs decay found at the vertex " << HiggsDecVtx->barcode() <<" (barcode)" << endl ; + + vector HiggsChildren; + + for ( HepMC::GenVertex::particles_out_const_iterator H0in = + HiggsDecVtx->particles_out_const_begin(); + H0in != HiggsDecVtx->particles_out_const_end(); + H0in++ ) + { + HiggsChildren.push_back(*H0in); + } + cout << " Number of Higgs (immediate) children = " << HiggsChildren.size() << endl ; + for (unsigned int ic=0; icprint() ; + } + } + + // select and store stable descendants of the Higgs + // + vector StableHiggsDesc ; + + if ( e.id().event() == 1 ) + cout << " Now let us list all descendents of the Higgs" << endl ; + for ( HepMC::GenVertex::particle_iterator + des=HiggsDecVtx->particles_begin(HepMC::descendants); + des!=HiggsDecVtx->particles_end(HepMC::descendants); des++ ) + { + if ( e.id().event() == 1 ) (*des)->print() ; + if ( (*des)->status() == 1 ) StableHiggsDesc.push_back(*des) ; + } + + HepMC::FourVector Mom2part ; + double XMass2part = 0.; + double XMass4part = 0.; + double XMass2pairs = 0.; + vector< HepMC::FourVector > Mom2partCont ; + + // browse the array of stable descendants + // and do 2-mu inv.mass + // + for ( unsigned int i=0; ipdg_id()) != 13 ) continue ; + + for ( unsigned int j=i+1; jpdg_id()) != 13 ) continue ; + // + // skip same charge combo's + // + if ( (StableHiggsDesc[i]->pdg_id()*StableHiggsDesc[j]->pdg_id()) > 0 ) + continue ; + // + // OK, opposite charges, do the job + // + Mom2part = HepMC::FourVector((StableHiggsDesc[i]->momentum().px()+StableHiggsDesc[j]->momentum().px()), + (StableHiggsDesc[i]->momentum().py()+StableHiggsDesc[j]->momentum().py()), + (StableHiggsDesc[i]->momentum().pz()+StableHiggsDesc[j]->momentum().pz()), + (StableHiggsDesc[i]->momentum().e()+StableHiggsDesc[j]->momentum().e())) ; + + XMass2part = Mom2part.m() ; + fHist2muMass->Fill( XMass2part ) ; + //cout << " counters : " << StableHiggsDesc[i]->barcode() << " " + // << StableHiggsDesc[j]->barcode() + // << " -> 2-part mass = " << XMass2part << endl ; + // + // store if 2-part. inv. mass fits into (roughly) Z-mass interval + // + if ( XMass2part > 80. && XMass2part < 100. ) + { + Mom2partCont.push_back(Mom2part) ; + } + } + } + + // make 4-part inv.mass + // + double px4, py4, pz4, e4; + px4=py4=pz4=e4=0. ; + if ( StableHiggsDesc.size() == 4 ) + { + for ( unsigned int i=0; imomentum().px(); + py4 += StableHiggsDesc[i]->momentum().py(); + pz4 += StableHiggsDesc[i]->momentum().pz(); + e4 += StableHiggsDesc[i]->momentum().e(); + } + XMass4part = HepMC::FourVector(px4,py4,pz4,e4).m() ; + fHist4muMass->Fill( XMass4part ) ; + } + //cout << " 4-part inv. mass = " << XMass4part << endl ; + + // make 2-pairs (ZZ) inv.mass + // + //cout << " selected Z-candidates in this event : " << Mom2partCont.size() << endl ; + for ( unsigned int i=0; iFill( XMass2pairs ) ; + //cout << " 2-pairs (ZZ) inv. mass = " << XMass2pairs << endl ; + } + } + + return ; + +} + +void HZZ4muAnalyzer::endJob() +{ + + fOutputFile->Write() ; + fOutputFile->Close() ; + + return ; +} + +#include "FWCore/Framework/interface/MakerMacros.h" +DEFINE_FWK_MODULE(HZZ4muAnalyzer); diff --git a/Validation/Generator/plugins/HZZ4muAnalyzer.h b/Validation/Generator/plugins/HZZ4muAnalyzer.h new file mode 100755 index 0000000000000..3cbf17ec6330a --- /dev/null +++ b/Validation/Generator/plugins/HZZ4muAnalyzer.h @@ -0,0 +1,26 @@ +#include "FWCore/Framework/interface/EDAnalyzer.h" + +class TFile; +class TH1D; + +class HZZ4muAnalyzer : public edm::EDAnalyzer +{ + public: + // + explicit HZZ4muAnalyzer( const edm::ParameterSet& ) ; + virtual ~HZZ4muAnalyzer() {} // no need to delete ROOT stuff + // as it'll be deleted upon closing TFile + + virtual void analyze( const edm::Event&, const edm::EventSetup& ) ; + virtual void beginJob( const edm::EventSetup& ) ; + virtual void endJob() ; + + private: + // + std::string fOutputFileName ; + TFile* fOutputFile ; + TH1D* fHist2muMass ; + TH1D* fHist4muMass ; + TH1D* fHistZZMass ; + +}; diff --git a/Validation/Generator/plugins/MCValidation.cc b/Validation/Generator/plugins/MCValidation.cc new file mode 100644 index 0000000000000..0f47d37ebe797 --- /dev/null +++ b/Validation/Generator/plugins/MCValidation.cc @@ -0,0 +1,372 @@ +#include "Validation/Generator/plugins/MCValidation.h" + +using namespace std ; +using namespace edm ; +using namespace reco; + + +const double pi = acos(-1.) ; +const double twopi = 2*pi ; +const double eMass = 0.00511 ; +const double eMass2 = eMass*eMass ; +const double muMass = 0.106 ; +const double muMass2 = muMass*muMass ; +const double WMass = 80.403 ; +const double WMass2 = WMass*WMass ; +const double ZMass = 91.1876 ; +const double ZMass2 = ZMass*ZMass ; + +MCValidation::MCValidation(const edm::ParameterSet& iConfig) + +{ + + s_genEvt = iConfig.getUntrackedParameter("genEvt") ; +// l_genEvt = iConfig.getParameter("genEvt") ; + + theOutfileName = iConfig.getUntrackedParameter("OutfileName") ; + +} + +MCValidation::~MCValidation() +{ + + // do anything here that needs to be done at desctruction time + // (e.g. close files, deallocate resources etc.) + +} + +void +MCValidation::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) +{ + using namespace edm; + + nEvt++ ; + + outfile << " anls fn Evt no. " << iEvent.id().event() << endl ; + +// HandlegenPar ; +// iEvent.getByLabel(s_genEvt,genPar) ; +// const GenParticleCollection &genColl = (*(genPar.product())) ; + HandlegenPar ; + iEvent.getByLabel("genParticleCandidates",genPar) ; + const CandidateCollection &genColl = (*(genPar.product())) ; + genparticles(genColl) ; + +} // end analyse + +void MCValidation::genparticles(const CandidateCollection &genParticles){ //GenParticleCollection &genParticles){ + + for(unsigned i = 0; i < genParticles.size(); ++ i) { + +// const GenParticle & part = (genParticles)[i]; + const Candidate & part = (genParticles)[i]; + + int pid = part.pdgId(), pst = part.status(), pch = part.charge(); + unsigned nDau = part.numberOfDaughters() ; + //unsigned nMom = part.numberOfMothers() ; + double ppt = part.pt(), peta = part.eta(), pphi = part.phi(), pmass = part.mass(); + //double pvx = part.vx(), pvy = part.vy(), pvz = part.vz(); + + outfile << " pid " << pid << " pmass " << pmass << " pst " << pst << endl ; + + switch(abs(pid)){ + + case 11 : { // e+ OR e- + + gen_e_pT_->Fill(pch*ppt) ; + + if(pch>0){ + gen_eP_eta_->Fill(peta); gen_eP_phi_->Fill(pphi); + }else if(pch<0){ + gen_eM_eta_->Fill(peta); gen_eM_phi_->Fill(pphi); + } + + icode = 1 ; + break ; + } + + case 13 : { // mu+ OR mu- + + gen_e_pT_->Fill(pch*ppt) ; + + if(pch>0){ + gen_eP_eta_->Fill(peta); gen_eP_phi_->Fill(pphi); + }else if(pch<0){ + gen_eM_eta_->Fill(peta); gen_eM_phi_->Fill(pphi); + } + + icode = 2 ; + break ; + } + + case 22 : { // gamma + + gen_gamma_pT_->Fill(pch*ppt) ; gen_gamma_eta_->Fill(peta); gen_gamma_phi_->Fill(pphi); + + icode = 3 ; + break ; + + } + + case 24 : { // W+ OR W- + + if(pst!=3) break ; + + gen_W_pT_->Fill(ppt) ; + if(pch>0){ + gen_Wp_eta_->Fill(peta) ; gen_Wp_phi_->Fill(pphi) ; + }else if(pch<0){ + gen_Wm_eta_->Fill(peta) ; gen_Wm_phi_->Fill(pphi) ; + } + + if(nDau > 1){ + + for(unsigned i = 0; i < nDau; ++i){ + + const Candidate *dau = part.daughter(i) ; + int did = dau->pdgId(), dst = dau->status(), dch = dau->charge() ; + double dpt = dau->pt(), deta = dau->eta(), dphi = dau->phi();// dmass = dau->mass(); + // double dvx = dau->vx(), dvy = dau->vy(), dvz = dau->vz(); + + if(abs(did) == 11 && dst==1){ + dau_e_pT_->Fill(dpt) ; + if(dch>0){ + dau_eP_eta_->Fill(deta); dau_eP_phi_->Fill(dphi) ; + }else if(dch<0){ + dau_eM_eta_->Fill(deta); dau_eM_phi_->Fill(dphi) ; + } + } + + if(abs(did)== 13 && dst==1){ + dau_mu_pT_->Fill(dpt); + if(dch>0){ + dau_muP_eta_->Fill(deta); dau_muP_phi_->Fill(dphi) ; + }else if(dch<0){ + dau_muM_eta_->Fill(deta); dau_muM_phi_->Fill(dphi) ; + } + } + + } + + } + + icode = 4 ; + break ; + + } + + case 23 : { // Z0 + + gen_Z_pT_->Fill(ppt) ; + gen_Z_eta_->Fill(peta) ; gen_Z_phi_->Fill(pphi) ; + + if(nDau > 1){ + + for(unsigned i = 0; i < nDau; ++i){ + + const Candidate *dau = part.daughter(i) ; + int did = dau->pdgId(), dst = dau->status(), dch = dau->charge() ; + double dpt = dau->pt(), deta = dau->eta(), dphi = dau->phi(); //dmass = dau->mass(); + //double dvx = dau->vx(), dvy = dau->vy(), dvz = dau->vz(); + + if(abs(did) == 11 && dst==1){ + dau_e_pT_->Fill(dpt) ; + if(dch>0){ + dau_eP_eta_->Fill(deta); dau_eP_phi_->Fill(dphi) ; + }else if(dch<0){ + dau_eM_eta_->Fill(deta); dau_eM_phi_->Fill(dphi) ; + } + } + + if(abs(did)== 13 && dst==1){ + dau_mu_pT_->Fill(dpt); + if(dch>0){ + dau_muP_eta_->Fill(deta); dau_muP_phi_->Fill(dphi) ; + }else if(dch<0){ + dau_muM_eta_->Fill(deta); dau_muM_phi_->Fill(dphi) ; + } + } + + } + + } + + icode = 5 ; + break ; + + } + + case 35 : { // H0 + + gen_H0_pT_->Fill(ppt) ; + gen_H0_eta_->Fill(peta) ; gen_H0_phi_->Fill(pphi) ; + + if(nDau > 1){ + + for(unsigned i = 0; i < nDau; ++i){ + + const Candidate *dau = part.daughter(i) ; + int did = dau->pdgId(), dst = dau->status(), dch = dau->charge() ; + double dpt = dau->pt(), deta = dau->eta(), dphi = dau->phi();// dmass = dau->mass(); + //double dvx = dau->vx(), dvy = dau->vy(), dvz = dau->vz(); + + if(abs(did) == 11 && dst==1){ + dau_e_pT_->Fill(dpt) ; + if(dch>0){ + dau_eP_eta_->Fill(deta); dau_eP_phi_->Fill(dphi) ; + }else if(dch<0){ + dau_eM_eta_->Fill(deta); dau_eM_phi_->Fill(dphi) ; + } + } + + if(abs(did)== 13 && dst==1){ + dau_mu_pT_->Fill(dpt); + if(dch>0){ + dau_muP_eta_->Fill(deta); dau_muP_phi_->Fill(dphi) ; + }else if(dch<0){ + dau_muM_eta_->Fill(deta); dau_muM_phi_->Fill(dphi) ; + } + } + + if(did== 22 && dst==1){ + dau_gamma_pT_->Fill(dpt); dau_gamma_eta_->Fill(deta); dau_gamma_phi_->Fill(dphi) ; + } + + } + + } + + icode = 6 ; + break ; + + } + + case 32 : { // Z' + + gen_Zprime_pT_->Fill(ppt) ; + gen_Zprime_eta_->Fill(peta) ; gen_Zprime_phi_->Fill(pphi) ; + + if(nDau > 1){ + + for(unsigned i = 0; i < nDau; ++i){ + + const Candidate *dau = part.daughter(i) ; + int did = dau->pdgId(), dst = dau->status(), dch = dau->charge() ; + double dpt = dau->pt(), deta = dau->eta(), dphi = dau->phi(); //dmass = dau->mass(); + // double dvx = dau->vx(), dvy = dau->vy(), dvz = dau->vz(); + + if(abs(did) == 11 && dst==1){ + dau_e_pT_->Fill(dpt) ; + if(dch>0){ + dau_eP_eta_->Fill(deta); dau_eP_phi_->Fill(dphi) ; + }else if(dch<0){ + dau_eM_eta_->Fill(deta); dau_eM_phi_->Fill(dphi) ; + } + } + + if(abs(did)== 13 && dst==1){ + dau_mu_pT_->Fill(dpt); + if(dch>0){ + dau_muP_eta_->Fill(deta); dau_muP_phi_->Fill(dphi) ; + }else if(dch<0){ + dau_muM_eta_->Fill(deta); dau_muM_phi_->Fill(dphi) ; + } + } + + } + + } + + icode = 7 ; + break ; + + } + + default : icode = 0 ; break ; + + } // switch ends + + } // loop over gen particles ends + +} // end genparticles + +void +MCValidation::beginJob(const edm::EventSetup& iSetup) +{ + + iSetup.get().get(theMagneticField) ; + + outfile.open(theOutfileName.c_str()) ; + + edm::Service fs; + + Tgen = new TTree("Tgen", "GEN") ; + Tgen->Branch("iEvt", &iEvt, "iEvt/I") ; + Tgen->Branch("iRun", &iRun, "iRun/I") ; + + mom_pT_ = fs->make("mom_pT_", "p_{T} mother", 100, -50., 50.) ; + mom_eta_ = fs->make("mom_eta_", "#eta mother", 50, -3.5, 3.5) ; + mom_phi_ = fs->make("mom_phi_", "#phi mother", 50, -3.5, 3.5) ; + + gen_e_pT_ = fs->make("gen_e_pT_", "Gen e p_{T}", 100, -50., 50.) ; + gen_eP_eta_ = fs->make("gen_eP_eta_", "Gen e^{+} #eta spectrum", 50, -3.5, 3.5) ; + gen_eP_phi_ = fs->make("gen_eP_phi_", "Gen e^{+} #phi spectrum", 50, -0.5, 3.5) ; + gen_eM_eta_ = fs->make("gen_eM_eta_", "Gen e^{-} #eta spectrum", 50, -3.5, 3.5) ; + gen_eM_phi_ = fs->make("gen_eM_phi_", "Gen e^{-} #phi spectrum", 50, -0.5, 3.5) ; + + gen_mu_pT_ = fs->make("gen_mu_pT_", "Gen mu p_{T}", 100, -50., 50.) ; + gen_muP_eta_ = fs->make("gen_muP_eta_", "Gen mu^{+} #eta spectrum", 50, -3.5, 3.5) ; + gen_muP_phi_ = fs->make("gen_muP_phi_", "Gen mu^{+} #phi spectrum", 50, -0.5, 3.5) ; + gen_muM_eta_ = fs->make("gen_muM_eta_", "Gen mu^{-} #eta spectrum", 50, -3.5, 3.5) ; + gen_muM_phi_ = fs->make("gen_muM_phi_", "Gen mu^{-} #phi spectrum", 50, -0.5, 3.5) ; + + gen_gamma_pT_ = fs->make("gen_gamma_pT_", "Gen #gamma p_{T}", 50, 0., 50.) ; + gen_gamma_eta_ = fs->make("gen_gamma_eta_", "Gen #gamma #eta spectrum", 50, -3.5, 3.5) ; + gen_gamma_phi_ = fs->make("gen_gamma_phi_", "Gen #gamma #phi spectrum", 50, -0.5, 3.5) ; + + gen_W_pT_ = fs->make("gen_W_pT_", "Gen W p_{T}", 100, -50., 50.) ; + gen_Wp_eta_ = fs->make("gen_Wp_eta_", "Gen W^{+} #eta spectrum", 50, -3.5, 3.5) ; + gen_Wp_phi_ = fs->make("gen_Wp_phi_", "Gen W^{+} #phi spectrum", 50, -0.5, 3.5) ; + gen_Wm_eta_ = fs->make("gen_Wm_eta_", "Gen W^{-} #eta spectrum", 50, -3.5, 3.5) ; + gen_Wm_phi_ = fs->make("gen_Wm_phi_", "Gen W^{-} #phi spectrum", 50, -0.5, 3.5) ; + + gen_Z_pT_ = fs->make("gen_Z_pT_", "Gen Z p_{T}", 50, 0., 50.) ; + gen_Z_eta_ = fs->make("gen_Z_eta_", "Gen Z #eta spectrum", 50, -3.5, 3.5) ; + gen_Z_phi_ = fs->make("gen_Z_phi_", "Gen Z #phi spectrum", 50, -0.5, 3.5) ; + + gen_H0_pT_ = fs->make("gen_H0_pT_", "Gen H^{0} p_{T}", 50, 0., 50.) ; + gen_H0_eta_ = fs->make("gen_H0_eta_", "Gen H^{0} #eta spectrum", 50, -3.5, 3.5) ; + gen_H0_phi_ = fs->make("gen_H0_phi_", "Gen H^{0} #phi spectrum", 50, -0.5, 3.5) ; + + gen_Zprime_pT_ = fs->make("gen_Zprime_pT_", "Gen Z^{'} p_{T}", 50, 0., 50.) ; + gen_Zprime_eta_ = fs->make("gen_Zprime_eta_", "Gen Z^{'} #eta spectrum", 50, -3.5, 3.5) ; + gen_Zprime_phi_ = fs->make("gen_Zprime_phi_", "Gen Z^{'} #phi spectrum", 50, -0.5, 3.5) ; + + dau_e_pT_ = fs->make("dau_e_pT_", "Dau e p_{T}", 100, -50., 50.) ; + dau_eP_eta_ = fs->make("dau_eP_eta_", "Dau e^{+} #eta", 50, -3.5, 3.5) ; + dau_eP_phi_ = fs->make("dau_eP_phi_", "Dau e^{+} #phi", 50, -0.5, 3.5) ; + dau_eM_eta_ = fs->make("dau_eM_eta_", "Dau e^{-} #eta", 50, -3.5, 3.5) ; + dau_eM_phi_ = fs->make("dau_eM_phi_", "Dau e^{-} #phi", 50, -0.5, 3.5) ; + + dau_mu_pT_ = fs->make("dau_mu_pT_", "Dau mu p_{T}", 100, -50., 50.) ; + dau_muP_eta_ = fs->make("dau_muP_eta_", "Dau mu^{+} #eta", 50, -3.5, 3.5) ; + dau_muP_phi_ = fs->make("dau_muP_phi_", "Dau mu^{+} #phi", 50, -0.5, 3.5) ; + dau_muM_eta_ = fs->make("dau_muM_eta_", "Dau mu^{-} #eta", 50, -3.5, 3.5) ; + dau_muM_phi_ = fs->make("dau_muM_phi_", "Dau mu^{-} #phi", 50, -0.5, 3.5) ; + + dau_gamma_pT_ = fs->make("dau_gamma_pT_", "Dau #gamma p_{T}", 50, 0., 50.) ; + dau_gamma_eta_ = fs->make("dau_gamma_eta_", "Dau #gamma #eta", 50, -2.5, 2.5) ; + dau_gamma_phi_ = fs->make("dau_gamma_phi_", "Dau #gamma #phi", 50, -0.5, 3.5) ; + + nEvt = iEvt = iRun = 0 ; + +} + +void +MCValidation::endJob() { + +} + +#include "FWCore/Framework/interface/MakerMacros.h" +DEFINE_FWK_MODULE(MCValidation); + diff --git a/Validation/Generator/plugins/MCValidation.h b/Validation/Generator/plugins/MCValidation.h new file mode 100644 index 0000000000000..dedaea8a4c22a --- /dev/null +++ b/Validation/Generator/plugins/MCValidation.h @@ -0,0 +1,141 @@ +// -*- C++ -*- +// +// Package: MCValidation +// Class: MCValidation +// +/**\class MCValidation MCValidation.cc testAnalysis/MCValidation/src/MCValidation.cc + + Description: + + Implementation: + + + ---------------------------- + | particle | pdgId | icode | + ---------------------------- + | e | 11 | 1 | + | mu | 13 | 2 | + | gamma | 22 | 3 | + | Z | 23 | 4 | + | W | 24 | 5 | + | H0 | 35 | 6 | + | Z' | 32 | 7 | + ---------------------------- + +*/ +// +// Original Author: Devdatta MAJUMDER +// Created: Thu Apr 10 19:55:14 CEST 2008 +// $Id: MCValidation.cc,v 1.2 2009/07/31 23:21:08 kharchil Exp $ +// +// + +// system include files +#include +#include +#include +#include +#include + +// user include files +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/EDAnalyzer.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/TriggerNames.h" + +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ServiceRegistry/interface/Service.h" +#include "PhysicsTools/UtilAlgos/interface/TFileService.h" +#include "MagneticField/Engine/interface/MagneticField.h" +#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h" +#include "CLHEP/Vector/LorentzVector.h" + +///#include "SimDataFormats/HepMCProduct/interface/HepMCProduct.h" +#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h" + +#include "DataFormats/HepMCCandidate/interface/GenParticle.h" +#include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h" + +#include "DataFormats/Candidate/interface/Candidate.h" +#include "DataFormats/Candidate/interface/CandMatchMap.h" + +#include "DataFormats/EgammaCandidates/interface/Photon.h" +#include "DataFormats/EgammaCandidates/interface/PhotonFwd.h" + +#include "DataFormats/MuonReco/interface/Muon.h" +#include "DataFormats/MuonReco/interface/MuonFwd.h" + +#include +#include +#include + +#include "DataFormats/VertexReco/interface/Vertex.h" +#include "DataFormats/VertexReco/interface/VertexFwd.h" + +#include "RecoVertex/KalmanVertexFit/interface/KalmanVertexFitter.h" +#include "RecoVertex/VertexPrimitives/interface/TransientVertex.h" +#include "TrackingTools/PatternTools/interface/TwoTrackMinimumDistance.h" +#include "TrackingTools/Records/interface/TransientTrackRecord.h" +#include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h" +#include "TrackingTools/TransientTrack/interface/TransientTrack.h" +#include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h" + +#include "TFile.h" +#include "TH1F.h" +#include "TTree.h" +#include "TH2F.h" +#include "TCanvas.h" +#include "TPostScript.h" +#include "TStyle.h" +#include "TGraphErrors.h" + +class MCValidation : public edm::EDAnalyzer { +public: + explicit MCValidation(const edm::ParameterSet&); + ~MCValidation(); + +private: + virtual void beginJob(const edm::EventSetup&) ; + virtual void analyze(const edm::Event&, const edm::EventSetup&); + virtual void endJob() ; + + void genparticles(const reco::CandidateCollection&) ; +// void genparticles(const GenParticleCollection&) ; + + // ----------member data --------------------------- + + edm::ESHandle theMagneticField ; + +// edm::InputTag l_genEvt ; + std::string s_genEvt ; + + std::string theOutfileName ; + ofstream outfile ; + + edm::Service fs ; + + TTree* Tgen ; + + TH1F *mom_pT_, *mom_eta_, *mom_phi_ ; + + TH1F *gen_e_pT_, *gen_eP_eta_, *gen_eP_phi_, + *gen_eM_eta_, *gen_eM_phi_ ; + TH1F *gen_mu_pT_, *gen_muP_eta_, *gen_muP_phi_, + *gen_muM_eta_, *gen_muM_phi_ ; + TH1F *gen_gamma_pT_, *gen_gamma_eta_, *gen_gamma_phi_ ; + TH1F *gen_W_pT_, *gen_Wp_eta_, *gen_Wp_phi_, + *gen_Wm_eta_, *gen_Wm_phi_ ; + TH1F *gen_Z_pT_, *gen_Z_eta_, *gen_Z_phi_ ; + TH1F *gen_H0_pT_, *gen_H0_eta_, *gen_H0_phi_ ; + TH1F *gen_Zprime_pT_, *gen_Zprime_eta_, *gen_Zprime_phi_ ; + + TH1F *dau_e_pT_, *dau_eP_eta_, *dau_eP_phi_, + *dau_eM_eta_, *dau_eM_phi_ ; + TH1F *dau_mu_pT_, *dau_muP_eta_, *dau_muP_phi_, + *dau_muM_eta_, *dau_muM_phi_ ; + TH1F *dau_gamma_pT_, *dau_gamma_eta_, *dau_gamma_phi_ ; + + unsigned nEvt, iEvt, iRun ; + short icode ; + +}; diff --git a/Validation/Generator/plugins/MuMuAnalyzer.cc b/Validation/Generator/plugins/MuMuAnalyzer.cc new file mode 100644 index 0000000000000..20026e6f6e51e --- /dev/null +++ b/Validation/Generator/plugins/MuMuAnalyzer.cc @@ -0,0 +1,279 @@ +/* This is en example for an Analyzer of a Herwig HepMCProduct + and looks for muons pairs and fills a histogram + with the invaraint mass of the four. +*/ + +// +// Original Author: Fabian Stoeckli +// Created: Tue Nov 14 13:43:02 CET 2006 +// $Id: MuMuAnalyzer.cc,v 1.3 2009/07/31 23:21:08 kharchil Exp $ +// +// + + +// system include files +#include +#include + +// user include files +#include "Validation/Generator/plugins/MuMuAnalyzer.h" +#include "DataFormats/JetReco/interface/GenJet.h" + + +//#include +#include +#include + +#include "DataFormats/HepMCCandidate/interface/GenParticle.h" +#include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h" +//#include "DataFormats/HepMCCandidate/interface/GenParticleCandidate.h" +#include "DataFormats/Candidate/interface/Candidate.h" +#include "DataFormats/Candidate/interface/CandidateFwd.h" +#include "DataFormats/Candidate/interface/Particle.h" +#include "DataFormats/Math/interface/LorentzVector.h" + +#include "TLorentzVector.h" + +MuMuAnalyzer::MuMuAnalyzer(const edm::ParameterSet& iConfig) +{ + /// Copy plots from Steve exactly + outputFilename=iConfig.getUntrackedParameter("OutputFilename","dummy.root"); + J1Pt_histo = new TH1F("J1pT","J1pT",38,0,220); + J2Pt_histo = new TH1F("J2pT","J2pT",38,0,120); + E1Pt_histo = new TH1F("E1pT","E1pT",38,0,150); + E2Pt_histo = new TH1F("E2pT","E2pT",38,0,150); + ZPz_histo = new TH1F("ZPz","ZPz",100,0,20); + ZPt_histo = new TH1F("ZPt","ZPt",38,0,200); + J1Eta_histo = new TH1F("J1Eta_histo", "J1Eta_histo", 40, -3, 3); + J2Eta_histo = new TH1F("J2Eta_histo", "J2Eta_histo", 40, -3, 3); + JDelR_histo = new TH1F("JDelR_histo", "JDelR_histo", 38, 0, 6); + EJDelR_histo = new TH1F("EJDelR_histo", "EJDelR_histo", 38, 0, 6); + JDelPhi_histo = new TH1F("JDelPhi_histo", "JDelPhi_histo", 38, 0, 3.2); + J1Phi_histo = new TH1F("J1Phi_histo", "J1Phi_histo", 38, 0, 5); + J2Phi_histo = new TH1F("J2Phi_histo", "J2Phi_histo", 38, 0, 5); + MuMu_invmass_histo = new TH1F("MuMU_invmass_histo","MuMu_invmass_histo",100,0,100); + // int event = 0 ; +} + +MuMuAnalyzer::~MuMuAnalyzer() +{ + +} + +void +MuMuAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) +{ + ++event; + using namespace edm; + using namespace std; + using namespace reco; + // get HepMC::GenEvent ... + //Handle evt_h; + //iEvent.getByType(evt_h); + //HepMC::GenEvent * evt = new HepMC::GenEvent(*(evt_h->GetEvent())); + + + // look for stable muons + // std::vector muons; + //muons.resize(0); + //for(HepMC::GenEvent::particle_iterator it = evt->particles_begin(); it != evt->particles_end(); ++it) { + // if(abs((*it)->pdg_id())==13 && (*it)->status()==1) { + // muons.push_back(*it); + // } + // } + + math::XYZTLorentzVector Lvector1(0,0,0,0); + math::XYZTLorentzVector Lvector2(0,0,0,0); + + // Handle mcEventHandle; + //try{ + //iEvent.getByLabel("source",mcEventHandle); + //}catch(...) {;} +typedef std::vector GenJetCollection; + Handle genJets; + iEvent.getByLabel( "iterativeCone7GenJetsNoNuBSM", genJets); + //Handle genJets; + //iEvent.getByLabel( "iterativeCone5GenJetsNoNuBSM", genJets); + int elec = 0; + // Handle genPart; + // iEvent.getByLabel("genParticleCandidates",genPart); + Handle genPart; + iEvent.getByLabel("genParticles",genPart); + std::vector elecEta; + std::vector elecPhi; + std::vector elecPx; + std::vector elecPy; + std::vector elecPz; + float MuMuInvaMass; + elecEta.clear(); + elecPhi.clear(); + elecPx.clear(); + elecPy.clear(); + elecPz.clear(); + math::XYZTLorentzVector MyJets1 (0,0,0,0); + math::XYZTLorentzVector MyJets2 (0,0,0,0); + //std::vector MyJets; + for(size_t q = 0; q < genPart->size(); q++) + { + const Candidate & p = (*genPart)[q]; + int id = p.pdgId(); + size_t NMoth = p.numberOfMothers() ; + int motherID1 = 0 ; + // int motherID2 = 0 ; + + if(abs(id) == 23) + { + ZPt_histo->Fill(sqrt(p.px()*p.px() + p.py()*p.py())); + } + if(abs(id) != 13) continue; + for ( size_t moth1=0; moth1 < NMoth; moth1++ ) + { + motherID1 = (p.mother(moth1))->pdgId(); + if(motherID1 == 23) + { + elec++; + elecEta.push_back(p.eta()); + elecPhi.push_back(p.phi()); + elecPx.push_back(p.px()); + elecPy.push_back(p.py()); + elecPz.push_back(p.pz()); + } + } + //if(abs(id) == 13 && sqrt(p.px()*p.px()+p.py()*p.py()) > 20 ) elec++; + } + if (elec > 1) + { + for(size_t elec1 = 0; elec1 < genPart->size()-1; elec1++) + { + const Candidate & p1 = (*genPart)[elec1]; + int id1 = p1.pdgId(); + //int motherId1 = p1.mother()->pdgId; + if(abs(id1) != 13) continue; + //TLorentzVector momentum = p1.momentum(); + //cout <<" PT = " << sqrt(p1.px()*p1.px() +p1.py()*p1.py()) << endl; + for(size_t elec2 = elec1; elec2 < genPart->size(); elec2++) + { + const Candidate & p2 = (*genPart)[elec2]; + int id2 = p2.pdgId(); + //int motherId2 = p2->mother(0); + if(abs(id2) != 13) continue; + } + } + + + if(genJets.isValid()){ + if(genJets->size() > 1) + { + int nmyJets = 0; + for(int Jets = 0; Jets < int(genJets->size()); Jets++) + { + int incone = 0; + const Candidate & J1 = (*genJets)[Jets]; + for(int elecs = 0; elecs < int(elecPhi.size()); elecs++) + { + float EJDelPhi = fabs(J1.phi()-elecPhi[elecs]); + if(EJDelPhi > 3.1415926) EJDelPhi = 6.2831852 - EJDelPhi; + float EJDelR = sqrt((J1.eta()-elecEta[elecs])*(J1.eta()-elecEta[elecs])+EJDelPhi*EJDelPhi); + + if (EJDelR < .2) { cout << EJDelR << endl; incone++;} + for(int elecs1 = elecs+1; elecs1 < int(elecPhi.size()); elecs1++) + { + if(elecs == int(elecPhi.size())) continue; + if(elecs == elecs1) continue; + MuMuInvaMass = sqrt(sqrt(elecPx[elecs]*elecPx[elecs]+elecPy[elecs]*elecPy[elecs]+elecPz[elecs]*elecPz[elecs])+sqrt(elecPx[elecs1]*elecPx[elecs1]+elecPy[elecs1]*elecPy[elecs1]+elecPz[elecs1]*elecPz[elecs1]) - sqrt((elecPx[elecs] + elecPx[elecs1])* (elecPx[elecs] + elecPx[elecs1]) + (elecPy[elecs] + elecPy[elecs1])* (elecPy[elecs] + elecPy[elecs1]) + (elecPz[elecs] + elecPz[elecs1])* (elecPz[elecs] + elecPz[elecs1]))) ; + MuMu_invmass_histo->Fill(MuMuInvaMass); + + } + } + //if(nmyJets > 1) continue; + if (incone == 0 && nmyJets == 0) MyJets1 = math::XYZTLorentzVector(J1.px(),J1.py(),J1.pz(),0.);///// ***** Filled e with 0 + if (incone == 0 && nmyJets == 1) MyJets2 = math::XYZTLorentzVector(J1.px(),J1.py(),J1.pz(),0.); + nmyJets++; + + + } + if(nmyJets >= 2) + { + + //for(int J1int = 0; J1int < genJets->size()-1 ; J1int++) + //{ + //const Candidate & J1 = (*genJets)[J1int]; + //const Candidate & J1 = (*genJets)[0]; + double PtJ1 = MyJets1.pt(); + //float PtJ1 = sqrt(J1.px()*J1.px() + J1.py()*J1.py()); + float J1Eta = MyJets1.eta(); + float J1Phi = MyJets1.phi(); + if(PtJ1 > 20) + { + //const Candidate & Temp = (*genJets)[J1int + 1 ]; + //if(sqrt(Temp.px()*Temp.px()+Temp.py()*Temp.py())) + // { + J1Eta_histo->Fill(J1Eta); + J1Pt_histo->Fill(PtJ1); + J1Phi_histo->Fill(J1Phi); + // } + //for(int J2int = J1int; J2int < genJets->size() ; J2int++) + // { + //const Candidate & J2 = (*genJets)[1]; + //const Candidate & J2 = (*genJets)[J2int]; + double PtJ2 = MyJets2.pt(); + //float PtJ2 = sqrt(J2.px()*J2.px() + J2.py()*J2.py()); + if(PtJ2 > 20) + { + float J2Eta = MyJets2.eta(); + float J2Phi = MyJets2.phi(); + J2Eta_histo->Fill(J2Eta); + J2Pt_histo->Fill(PtJ2); + float DelPhi = fabs(J1Phi-J2Phi); + if(DelPhi > 3.1415926) DelPhi = 6.2831852 - DelPhi; + JDelPhi_histo->Fill(DelPhi); + float DelR = sqrt((J2Eta - J1Eta)*(J2Eta - J1Eta)+DelPhi*DelPhi); + JDelR_histo->Fill(DelR); + J2Phi_histo->Fill(J2Phi); + } + } + + } + } + + } + + // if there are at least four muons + // calculate invarant mass of first two and fill it into histogram + math::XYZTLorentzVector tot_momentum; math::XYZTLorentzVector tot_mumomentum; + //float inv_mass = 0.0; double mu_invmass = 0.0; float Pt = 0; + for( GenJetCollection::const_iterator gen = genJets->begin(); gen != genJets->end() ; ++gen ) + { + //cout << "gen jet pt " << gen->pt() << endl ; + } + + } + + //cout << elec << event << endl; +} + +void +MuMuAnalyzer::beginJob(const edm::EventSetup&) +{ +} + +void +MuMuAnalyzer::endJob() { + // save histograms into file + TFile file("Test.root","RECREATE"); + J1Pt_histo->Write(); + J2Pt_histo->Write(); + MuMu_invmass_histo->Write(); + ZPt_histo->Write(); + JDelR_histo->Write(); + JDelPhi_histo->Write(); + J1Eta_histo->Write(); + J2Eta_histo->Write(); + J1Phi_histo->Write(); + J2Phi_histo->Write(); + file.Close(); + +} + +#include "FWCore/Framework/interface/MakerMacros.h" +DEFINE_FWK_MODULE(MuMuAnalyzer); diff --git a/Validation/Generator/plugins/MuMuAnalyzer.h b/Validation/Generator/plugins/MuMuAnalyzer.h new file mode 100644 index 0000000000000..7d46e0775a027 --- /dev/null +++ b/Validation/Generator/plugins/MuMuAnalyzer.h @@ -0,0 +1,30 @@ +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/EDAnalyzer.h" +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" + +#include "HepMC/GenEvent.h" +#include "HepMC/GenParticle.h" +#include "HepMC/WeightContainer.h" + +#include "TH1D.h" +#include "TFile.h" + +class MuMuAnalyzer : public edm::EDAnalyzer { + public: + explicit MuMuAnalyzer(const edm::ParameterSet&); + ~MuMuAnalyzer(); + + private: + virtual void beginJob(const edm::EventSetup&) ; + virtual void analyze(const edm::Event&, const edm::EventSetup&); + virtual void endJob() ; + + private: + std::string outputFilename; + TH1D* weight_histo; int event; + TH1F* MuMu_invmass_histo; TH1F* J1Eta_histo; TH1F* J2Eta_histo; + TH1F* Pt_histo; TH1F* J1Pt_histo; TH1F* J2Pt_histo; TH1F* EJDelR_histo; + TH1F* E1Pt_histo; TH1F* E2Pt_histo; TH1F* ZPz_histo; TH1F* ZPt_histo; + TH1F* JDelR_histo; TH1F* JDelPhi_histo; TH1F* J1Phi_histo; TH1F* J2Phi_histo; +}; diff --git a/Validation/Generator/plugins/NewAnalyzer.cc b/Validation/Generator/plugins/NewAnalyzer.cc new file mode 100644 index 0000000000000..cc4c396eaeef7 --- /dev/null +++ b/Validation/Generator/plugins/NewAnalyzer.cc @@ -0,0 +1,578 @@ +/* This is en example for an Analyzer of a Herwig HepMCProduct + and looks for muons pairs and fills a histogram + with the invaraint mass of the four. +*/ + +// +// Original Author: Kenneth Smith +// Created: Tue Nov 14 13:43:02 CET 2006 +// $Id: NewAnalyzer.cc,v 1.6 2009/07/31 23:21:08 kharchil Exp $ +// +// + +// system include files +#include +#include + +// user include files +#include "Validation/Generator/plugins/NewAnalyzer.h" +#include "DataFormats/JetReco/interface/GenJet.h" + + +//#include +//#include +//#include +#include + +#include "DataFormats/HepMCCandidate/interface/GenParticle.h" +#include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h" +//#include "DataFormats/HepMCCandidate/interface/GenParticleCandidate.h" +#include "DataFormats/Candidate/interface/Candidate.h" +#include "DataFormats/Candidate/interface/CandidateFwd.h" +#include "DataFormats/Candidate/interface/Particle.h" +#include "DataFormats/Math/interface/LorentzVector.h" + +#include "TLorentzVector.h" + +NewAnalyzer::NewAnalyzer(const edm::ParameterSet& iConfig) +{ + outputFilename=iConfig.getUntrackedParameter("OutputFilename","dummy.root"); + Jetmult_histo = new TH1F("Jetmult_histo","Jet multiplicity",5,0,5); + ChDauMult_histo = new TH1F("ChDauMult_histo","Charged daughters of tau",7,0,7); + ChDaupT_histo = new TH1F("ChDaupT_histo","Charged daughter pT",38,0,220); + J1Pt_histo = new TH1F("J1pT","J1 pT",38,0,220); + J2Pt_histo = new TH1F("J2pT","J2 pT",38,0,120); + JetPt1J = new TH1F("JetpT1J","Jet pT for 1 jet events",38,0,220); + JetPt2J = new TH1F("JetpT2J","Jet pT for 2 jet events",38,0,220); + JetPt3J = new TH1F("JetpT3J","Jet pT for 3 jet events",38,0,220); + JetPt4J = new TH1F("JetpT4J","Jet pT for 4 jet events",38,0,220); + Z1JJ1Pt_histo = new TH1F("Z1JJ1pT","Jet1 pT for Z+1J",38,0,220); + Z2JJ1Pt_histo = new TH1F("Z2JJ1pT","Jet1 pT for Z+2J",38,0,220); + Z2JJ2Pt_histo = new TH1F("Z2JJ2pT","Jet2 pT for Z+2J",38,0,120); + Z3JJ1Pt_histo = new TH1F("Z3JJ1pT","Jet1 pT for Z+3J",38,0,220); + Z3JJ2Pt_histo = new TH1F("Z3JJ2pT","Jet2 pT for Z+3J",38,0,120); + Z4JJ1Pt_histo = new TH1F("Z4JJ1pT","Jet1 pT for Z+4J",38,0,220); + Z4JJ2Pt_histo = new TH1F("Z4JJ2pT","Jet2 pT for Z+4J",38,0,120); + E1Pt_histo = new TH1F("E1pT","Electron pT",38,0,150); + E2Pt_histo = new TH1F("E2pT","E2 pT",38,0,150); + ZPz_histo = new TH1F("ZPz","ZPz",100,0,20); + ZPt_histo = new TH1F("ZPt_histo","Z Pt",38,0,200); + ZPt1J_histo = new TH1F("ZPt1J","Z Pt Z+1J",38,0,200); + ZPt2J_histo = new TH1F("ZPt2J","Z Pt Z+2J",38,0,200); + ZPt3J_histo = new TH1F("ZPt3J","Z Pt Z+3J",38,0,200); + ZPt4J_histo = new TH1F("ZPt4J","Z Pt Z+4J",38,0,200); + ZPt0J_histo = new TH1F("ZPt0J","Z Pt Z+0J",38,0,200); + J1Eta_histo = new TH1F("J1Eta_histo", "J1Eta_histo", 40, -3, 3); + Z1JJ1Eta_histo = new TH1F("Z1JJ1Eta_histo", "Z1JJ1Eta_histo", 40, -3, 3); + Z2JJ1Eta_histo = new TH1F("Z2JJ1Eta_histo", "Z2JJ1Eta_histo", 40, -3, 3); + Z2JJ2Eta_histo = new TH1F("Z2JJ2Eta_histo", "Z2JJ2Eta_histo", 40, -3, 3); + Z3JJ1Eta_histo = new TH1F("Z3JJ1Eta_histo", "Z3JJ1Eta_histo", 40, -3, 3); + Z3JJ2Eta_histo = new TH1F("Z3JJ2Eta_histo", "Z3JJ2Eta_histo", 40, -3, 3); + Z4JJ1Eta_histo = new TH1F("Z4JJ1Eta_histo", "Z4JJ1Eta_histo", 40, -3, 3); + Z4JJ2Eta_histo = new TH1F("Z4JJ2Eta_histo", "Z4JJ2Eta_histo", 40, -3, 3); + ZEta_histo = new TH1F("ZEta_histo", "Z Eta", 40, -3, 3); + ZRap_histo = new TH1F("ZRap_histo", "Z Rapidity", 40, -3, 3); + JDelR_histo = new TH1F("JDelR_histo", "JDelR_histo", 38, 0, 6); + Z2JJDelR_histo = new TH1F("Z2JJDelR_histo", "Z2JJDelR_histo", 38, 0, 6); + Z3JJDelR_histo = new TH1F("Z3JJDelR_histo", "Z3JJDelR_histo", 38, 0, 6); + Z4JJDelR_histo = new TH1F("Z4JJDelR_histo", "Z4JJDelR_histo", 38, 0, 6); + EJDelR_histo = new TH1F("EJDelR_histo", "EJDelR_histo", 38, 0, 6); + JDelPhi_histo = new TH1F("JDelPhi_histo", "JDelPhi_histo", 38, -3.2, 3.2); + Z2JJDelPhi_histo = new TH1F("Z2JJDelPhi_histo", "Z2JJDelPhi_histo", 38, -3.2, 3.2); + Z3JJDelPhi_histo = new TH1F("Z3JJDelPhi_histo", "Z3JJDelPhi_histo", 38, -3.2, 3.2); + Z4JJDelPhi_histo = new TH1F("Z4JJDelPhi_histo", "Z4JJDelPhi_histo", 38, -3.2, 3.2); + J1Phi_histo = new TH1F("J1Phi_histo", "J1Phi_histo", 38, 0, 5); + Z1JJ1Phi_histo = new TH1F("Z1JJ1Phi_histo", "Z1JJ1Phi_histo", 38, 0, 5); + Z2JJ1Phi_histo = new TH1F("Z2JJ1Phi_histo", "Z2JJ1Phi_histo", 38, 0, 5); + Z2JJ2Phi_histo = new TH1F("Z2JJ2Phi_histo", "Z2JJ2Phi_histo", 38, 0, 5); + Z3JJ1Phi_histo = new TH1F("Z3JJ1Phi_histo", "Z3JJ1Phi_histo", 38, 0, 5); + Z3JJ2Phi_histo = new TH1F("Z3JJ2Phi_histo", "Z3JJ2Phi_histo", 38, 0, 5); + Z4JJ1Phi_histo = new TH1F("Z4JJ1Phi_histo", "Z4JJ1Phi_histo", 38, 0, 5); + Z4JJ2Phi_histo = new TH1F("Z4JJ2Phi_histo", "Z4JJ2Phi_histo", 38, 0, 5); + Z_invmass_histo = new TH1F("Z_invmass_histo","Z_invmass_histo",200,0,200); + Z0J_invmass_histo = new TH1F("Z0J_invmass_histo","Z0J_invmass_histo",200,0,200); + Z1J_invmass_histo = new TH1F("Z1J_invmass_histo","Z1J_invmass_histo",200,0,200); + Z2J_invmass_histo = new TH1F("Z2J_invmass_histo","Z2J_invmass_histo",200,0,200); + Z3J_invmass_histo = new TH1F("Z3J_invmass_histo","Z3J_invmass_histo",200,0,200); + Z4J_invmass_histo = new TH1F("Z4J_invmass_histo","Z4J_invmass_histo",200,0,200); + // int event = 0 ; + tau_evt =0; +} + +NewAnalyzer::~NewAnalyzer() +{ +} + +void +NewAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) +{ + ++event; + using namespace edm; + using namespace std; + using namespace reco; + // get HepMC::GenEvent ... + //Handle evt_h; + //iEvent.getByType(evt_h); + //HepMC::GenEvent * evt = new HepMC::GenEvent(*(evt_h->GetEvent())); + + + // look for stable muons + // std::vector muons; + //muons.resize(0); + //for(HepMC::GenEvent::particle_iterator it = evt->particles_begin(); it != evt->particles_end(); ++it) { + // if(abs((*it)->pdg_id())==13 && (*it)->status()==1) { + // muons.push_back(*it); + // } + // } + + math::XYZTLorentzVector Lvector1(0,0,0,0); + math::XYZTLorentzVector Lvector2(0,0,0,0); + + // Handle mcEventHandle; + //try{ + //iEvent.getByLabel("source",mcEventHandle); + //}catch(...) {;} + typedef std::vector GenJetCollection; + Handle genJets; + //iEvent.getByLabel( "iterativeCone5GenJetsNoNuBSM", genJets); + //Handle genJets; + iEvent.getByLabel( "iterativeCone5GenJets", genJets); + int elec = 0; + //std::cout << "A" << std::endl; +//Handle genPart; +//iEvent.getByLabel("genParticleCandidates",genPart); + +// Handle genPart; + // iEvent.getByLabel("genParticles",genPart); + Handle genPart; + iEvent.getByLabel("genParticles",genPart); + //std::cout << "A" << std::endl; + std::vector elecEta; + std::vector elecPhi; + std::vector elecPx; + std::vector elecPy; + std::vector elecPz; + std::vector elecCh; + std::vector JetpT; + float EEInvaMass; + elecEta.clear(); + elecPhi.clear(); + elecPx.clear(); + elecPy.clear(); + elecPz.clear(); + double ZpT = 0; int taur = 0; + float ptot, etot; + double Jet1Pt, Jet2Pt; + Jet1Pt = 0; + Jet2Pt = 0; + math::XYZTLorentzVector MyJets1 (0,0,0,0); + math::XYZTLorentzVector MyJets2 (0,0,0,0); + for(size_t q = 0; q < genPart->size(); q++) + { + const Candidate & p = (*genPart)[q]; + int id = p.pdgId(); + size_t NMoth = p.numberOfMothers() ; + int motherID1 = 0 ; + //int motherID2 = 0 ; + if(abs(id) == 23) + { + //ZPt_histo->Fill(sqrt(p.px()*p.px() + p.py()*p.py())); + ZEta_histo->Fill(p.eta()); + ZRap_histo->Fill(p.rapidity()); + ZpT = p.pt(); + } + if(abs(id) != 15) continue; + int chd = 0; + taur++; + //cout << "Checking daughters" << endl; + if( p.numberOfDaughters() > 0) + { + for(int dau = 0; dau < int(p.numberOfDaughters()); dau++) + {//cout << "Checking charge" << endl; + if(p.daughter(dau)->charge() != 0) + { + ChDaupT_histo -> Fill(p.daughter(dau)->pt()); + chd++; + } + } + ChDauMult_histo->Fill(chd); + } + for ( size_t moth1=0; moth1pdgId(); + if(fabs(motherID1) == 23) + { + elec++; + elecEta.push_back(p.eta()); + elecPhi.push_back(p.phi()); + elecPx.push_back(p.px()); + elecPy.push_back(p.py()); + elecPz.push_back(p.pz()); + elecCh.push_back(p.charge()); + E1Pt_histo->Fill(sqrt(p.px()*p.px() + p.py() * p.py())); + } + } + + } + if(taur > 2) + { + cout << "More than 2 taus in event. Number of taus: " << taur << endl; + tau_evt++; + } + if (elec > 1) + { + ZPt_histo->Fill(ZpT); + if(genJets.isValid()){ + if(genJets->size() > 1) + { + + for(size_t elec1 = 0; int(elec1) < int(elec-1); elec1++) + { + for(size_t elec2 = elec1 + 1; int(elec2) < int(elec); elec2++) + { + if(elecCh[elec2] == elecCh[elec1]) + continue; + etot = sqrt(elecPx[elec1]*elecPx[elec1]+elecPy[elec1]*elecPy[elec1]+elecPz[elec1]*elecPz[elec1])+sqrt(elecPx[elec2]*elecPx[elec2]+elecPy[elec2]*elecPy[elec2]+elecPz[elec2]*elecPz[elec2]); + ptot = sqrt((elecPx[elec1] + elecPx[elec2])* (elecPx[elec1] + elecPx[elec2]) + (elecPy[elec1] + elecPy[elec2]) * (elecPy[elec2] + elecPy[elec1]) + (elecPz[elec1] + elecPz[elec2]) * (elecPz[elec1] + elecPz[elec2])); + EEInvaMass = sqrt((etot+ptot)*(etot-ptot)); + Z_invmass_histo->Fill(EEInvaMass); + } + } + int nmyJets = 0; + JetpT.clear(); + //Jet1Pt = 0.0; + //Jet2Pt = 0.0; + for(int Jets = 0; Jets < int(genJets->size()); Jets++) + { + int incone = 0; + const Candidate & J1 = (*genJets)[Jets]; + for(int elecs = 0; elecs < elec; elecs++) + { + float EJDelPhi = fabs(J1.phi()-elecPhi[elecs]); + if(EJDelPhi > 3.1415926) EJDelPhi = 6.2831852 - EJDelPhi; + float EJDelR = sqrt((J1.eta()-elecEta[elecs])*(J1.eta()-elecEta[elecs])+EJDelPhi*EJDelPhi); + + if (EJDelR < .2) { //cout << EJDelR << endl; + incone++;} + if(elecs == int(elecPhi.size())) continue; + } + //cout << J1.pt() << " Jet pT " << endl; + + if (incone == 0 && J1.pt() > 12) + { + if(nmyJets == 0 ) + { + MyJets1 = math::XYZTLorentzVector(J1.px(),J1.py(),J1.pz(),J1.energy()); + Jet1Pt = J1.pt(); + + //cout << "Jet 1 found " << endl; + } + if ( nmyJets == 1) + { + MyJets2 = math::XYZTLorentzVector(J1.px(),J1.py(),J1.pz(),J1.energy()); + Jet2Pt = J1.pt(); + //cout << "Jet 2 found " << endl; + } + nmyJets++; + JetpT.push_back(J1.pt()); + } + + } + if(JetpT.size() == 1) + for(int i = 0; i < int(JetpT.size()); i++) + { + JetPt1J->Fill(JetpT[i]); + } + if(JetpT.size() == 2) + for(int i = 0; i < int(JetpT.size()); i++) + { + JetPt2J->Fill(JetpT[i]); + } + if(JetpT.size() == 3) + for(int i = 0; i < int(JetpT.size()); i++) + { + JetPt3J->Fill(JetpT[i]); + } + if(JetpT.size() > 3) + for(int i = 0; i < int(JetpT.size()); i++) + { + JetPt4J->Fill(JetpT[i]); + } + //std::cout << nmyJets << " my Jets" << std::endl; + //std::cout << "Jet 1 pt " << Jet1Pt << " Jet 2 Pt " << Jet2Pt << std::endl; + if(nmyJets == 0 || (Jet1Pt < 12 && Jet2Pt < 12)) + { + for(size_t elec1 = 0; int(elec1) < int(elec-1); elec1++) + { + for(size_t elec2 = elec1 + 1; int(elec2) Fill(EEInvaMass); + ZPt0J_histo->Fill(ZpT); + } + } + Jetmult_histo->Fill(0); + } + if(nmyJets > 0 ) + { + double PtJ1 = MyJets1.pt(); + float J1Eta = MyJets1.eta(); + float J1Phi = MyJets1.phi(); + if(Jet1Pt > 12) + { + J1Pt_histo->Fill(PtJ1); + J1Phi_histo->Fill(J1Phi); + J1Eta_histo->Fill(J1Eta); + } + } + if(nmyJets == 1 ) + { + Jetmult_histo->Fill(1); + double PtJ1 = MyJets1.pt(); + cout << Jet1Pt << " Jet Pt" << endl; + float J1Eta = MyJets1.eta(); + float J1Phi = MyJets1.phi(); + if(Jet1Pt > 12) + { + Z1JJ1Eta_histo->Fill(J1Eta); + Z1JJ1Pt_histo->Fill(PtJ1); + Z1JJ1Phi_histo->Fill(J1Phi); + ZPt1J_histo->Fill(ZpT); + } + for(size_t elec1 = 0; int(elec1) < int(elec-1); elec1++) + { + for(size_t elec2 = elec1 + 1; int(elec2) < int(elec); elec2++) + { + if(elecCh[elec2] == elecCh[elec1]) + continue; + etot = sqrt(elecPx[elec1]*elecPx[elec1]+elecPy[elec1]*elecPy[elec1]+elecPz[elec1]*elecPz[elec1])+sqrt(elecPx[elec2]*elecPx[elec2]+elecPy[elec2]*elecPy[elec2]+elecPz[elec2]*elecPz[elec2]); + ptot = sqrt((elecPx[elec1] + elecPx[elec2])* (elecPx[elec1] + elecPx[elec2]) + (elecPy[elec1] + elecPy[elec2]) * (elecPy[elec2] + elecPy[elec1]) + (elecPz[elec1] + elecPz[elec2]) * (elecPz[elec1] + elecPz[elec2])); + EEInvaMass = sqrt((etot+ptot)*(etot-ptot)); + Z1J_invmass_histo->Fill(EEInvaMass); + } + } + } + if(nmyJets == 2) + { + Jetmult_histo->Fill(2); + double PtJ1 = MyJets1.pt(); + float J1Eta = MyJets1.eta(); + float J1Phi = MyJets1.phi(); + if(Jet1Pt > 12) + { + Z2JJ1Eta_histo->Fill(J1Eta); + Z2JJ1Pt_histo->Fill(PtJ1); + Z2JJ1Phi_histo->Fill(J1Phi); + double PtJ2 = MyJets2.pt(); + if(Jet2Pt > 12) + { + float J2Eta = MyJets2.eta(); + float J2Phi = MyJets2.phi(); + Z2JJ2Eta_histo -> Fill(J2Eta); + Z2JJ2Pt_histo -> Fill(PtJ2); + float DelPhi = fabs(J1Phi-J2Phi); + if(DelPhi > 3.1415926) DelPhi = 6.2831852 - DelPhi; + Z2JJDelPhi_histo -> Fill(DelPhi); + float DelR = sqrt((J2Eta - J1Eta)*(J2Eta - J1Eta)+DelPhi*DelPhi); + Z2JJDelR_histo -> Fill(DelR); + Z2JJ2Phi_histo -> Fill(J2Phi); + ZPt2J_histo->Fill(ZpT); + } + } + for(size_t elec1 = 0; int(elec1) < int(elec-1); elec1++) + { + for(size_t elec2 = elec1 +1 ; int(elec2) < int(elec); elec2++) + { + if(elecCh[elec2] == elecCh[elec1]) + continue; + etot = sqrt(elecPx[elec1]*elecPx[elec1]+elecPy[elec1]*elecPy[elec1]+elecPz[elec1]*elecPz[elec1])+sqrt(elecPx[elec2]*elecPx[elec2]+elecPy[elec2]*elecPy[elec2]+elecPz[elec2]*elecPz[elec2]); + ptot = sqrt((elecPx[elec1] + elecPx[elec2])* (elecPx[elec1] + elecPx[elec2]) + (elecPy[elec1] + elecPy[elec2]) * (elecPy[elec2] + elecPy[elec1]) + (elecPz[elec1] + elecPz[elec2]) * (elecPz[elec1] + elecPz[elec2])); + EEInvaMass = sqrt((etot+ptot)*(etot-ptot)); + Z2J_invmass_histo->Fill(EEInvaMass); + + } + } + } + if(nmyJets == 3) + { + Jetmult_histo->Fill(3); + double PtJ1 = MyJets1.pt(); + float J1Eta = MyJets1.eta(); + float J1Phi = MyJets1.phi(); + if(Jet1Pt > 12) + { + Z3JJ1Eta_histo->Fill(J1Eta); + Z3JJ1Pt_histo->Fill(PtJ1); + Z3JJ1Phi_histo->Fill(J1Phi); + double PtJ2 = MyJets2.pt(); + if(Jet2Pt > 12) + { + float J2Eta = MyJets2.eta(); + float J2Phi = MyJets2.phi(); + Z3JJ2Eta_histo->Fill(J2Eta); + Z3JJ2Pt_histo->Fill(PtJ2); + float DelPhi = fabs(J1Phi-J2Phi); + if(DelPhi > 3.1415926) DelPhi = 6.2831852 - DelPhi; + Z3JJDelPhi_histo->Fill(DelPhi); + float DelR = sqrt((J2Eta - J1Eta)*(J2Eta - J1Eta)+DelPhi*DelPhi); + Z3JJDelR_histo->Fill(DelR); + Z3JJ2Phi_histo->Fill(J2Phi); + ZPt3J_histo->Fill(ZpT); + } + } + for(size_t elec1 = 0; int(elec1) < int(elec-1); elec1++) + { + for(size_t elec2 = elec1 +1 ; int(elec2) < int(elec); elec2++) + { + if(elecCh[elec2] == elecCh[elec1]) + continue; + etot = sqrt(elecPx[elec1]*elecPx[elec1]+elecPy[elec1]*elecPy[elec1]+elecPz[elec1]*elecPz[elec1])+sqrt(elecPx[elec2]*elecPx[elec2]+elecPy[elec2]*elecPy[elec2]+elecPz[elec2]*elecPz[elec2]); + ptot = sqrt((elecPx[elec1] + elecPx[elec2])* (elecPx[elec1] + elecPx[elec2]) + (elecPy[elec1] + elecPy[elec2]) * (elecPy[elec2] + elecPy[elec1]) + (elecPz[elec1] + elecPz[elec2]) * (elecPz[elec1] + elecPz[elec2])); + EEInvaMass = sqrt((etot+ptot)*(etot-ptot)); + Z3J_invmass_histo->Fill(EEInvaMass); + + } + } + } + if(nmyJets > 3) + { + Jetmult_histo->Fill(4); + //std::cout << "C7" << std::endl; + double PtJ1 = MyJets1.pt(); + float J1Eta = MyJets1.eta(); + float J1Phi = MyJets1.phi(); + if(Jet1Pt > 12) + { + //std::cout << "C8" << std::endl; + Z4JJ1Eta_histo->Fill(J1Eta); + Z4JJ1Pt_histo->Fill(PtJ1); + Z4JJ1Phi_histo->Fill(J1Phi); + double PtJ2 = MyJets2.pt(); + //std::cout << "C9" << std::endl; + if(Jet2Pt > 12) + { + float J2Eta = MyJets2.eta(); + float J2Phi = MyJets2.phi(); + //std::cout << "C10" << std::endl; + Z4JJ2Eta_histo->Fill(J2Eta); + Z4JJ2Pt_histo->Fill(PtJ2); + float DelPhi = fabs(J1Phi-J2Phi); + if(DelPhi > 3.1415926) DelPhi = 6.2831852 - DelPhi; + Z4JJDelPhi_histo->Fill(DelPhi); + float DelR = sqrt((J2Eta - J1Eta)*(J2Eta - J1Eta)+DelPhi*DelPhi); + Z4JJDelR_histo->Fill(DelR); + Z4JJ2Phi_histo->Fill(J2Phi); + ZPt4J_histo->Fill(ZpT); + } + } + for(size_t elec1 = 0; int(elec1) < int(elec-1); elec1++) + { + for(size_t elec2 = elec1+1; int(elec2) < int(elec); elec2++) + { + if(elecCh[elec2] == elecCh[elec1]) + continue; + etot = sqrt(elecPx[elec1]*elecPx[elec1]+elecPy[elec1]*elecPy[elec1]+elecPz[elec1]*elecPz[elec1])+sqrt(elecPx[elec2]*elecPx[elec2]+elecPy[elec2]*elecPy[elec2]+elecPz[elec2]*elecPz[elec2]); + ptot = sqrt((elecPx[elec1] + elecPx[elec2])* (elecPx[elec1] + elecPx[elec2]) + (elecPy[elec1] + elecPy[elec2]) * (elecPy[elec2] + elecPy[elec1]) + (elecPz[elec1] + elecPz[elec2]) * (elecPz[elec1] + elecPz[elec2])); + EEInvaMass = sqrt((etot+ptot)*(etot-ptot)); + Z4J_invmass_histo->Fill(EEInvaMass); + + } + } + } + } + else + { + ZPt0J_histo->Fill(ZpT); + } + + std::cout << "E" << std::endl; + } + + // if there are at least four muons + // calculate invarant mass of first two and fill it into histogram + math::XYZTLorentzVector tot_momentum; math::XYZTLorentzVector tot_mumomentum; + //float inv_mass = 0.0; double mu_invmass = 0.0; + //float Pt = 0; + for( GenJetCollection::const_iterator gen = genJets->begin(); gen != genJets->end() ; ++gen ) + { + //cout << "gen jet pt " << gen->pt() << endl ; + } + + } + + //cout << elec << event << endl; +} + +void +NewAnalyzer::beginJob(const edm::EventSetup&) +{ +} + +void +NewAnalyzer::endJob() { + // save histograms into file + TFile file("ZTT_final.root","RECREATE"); + std::cout << tau_evt << " : Number of taus" << std::endl; + ChDauMult_histo->Write(); + ChDaupT_histo->Write(); + J1Pt_histo->Write(); + J2Pt_histo->Write(); + Z_invmass_histo->Write(); + Z1J_invmass_histo->Write(); + Z2J_invmass_histo->Write(); + Z3J_invmass_histo->Write(); + Z4J_invmass_histo->Write(); + ZPt_histo->Write(); + JDelR_histo->Write(); + JDelPhi_histo->Write(); + Z1JJ1Eta_histo->Write(); + ZEta_histo->Write(); + Z1JJ1Phi_histo->Write(); + Z2JJ1Eta_histo->Write(); + Z2JJ2Eta_histo->Write(); + Z2JJ1Phi_histo->Write(); + Z2JJ2Phi_histo->Write(); + Z3JJ1Eta_histo->Write(); + Z3JJ2Eta_histo->Write(); + Z3JJ1Phi_histo->Write(); + Z3JJ2Phi_histo->Write(); + Z4JJ1Eta_histo->Write(); + Z4JJ2Eta_histo->Write(); + Z4JJ1Phi_histo->Write(); + Z4JJ2Phi_histo->Write(); + Jetmult_histo->Write(); + Z2JJDelR_histo->Write(); + Z2JJ2Phi_histo->Write(); + Z3JJDelR_histo->Write(); + Z3JJ2Phi_histo->Write(); + Z4JJDelR_histo->Write(); + Z4JJ2Phi_histo->Write(); + JetPt1J->Write(); + JetPt2J->Write(); + JetPt3J->Write(); + JetPt4J->Write(); + ZRap_histo->Write(); + ZPt0J_histo->Write(); + ZPt1J_histo->Write(); + ZPt2J_histo->Write(); + ZPt3J_histo->Write(); + ZPt4J_histo->Write(); + Z0J_invmass_histo->Write(); + Z1J_invmass_histo->Write(); + Z2J_invmass_histo->Write(); + Z3J_invmass_histo->Write(); + Z4J_invmass_histo->Write(); + E1Pt_histo->Write(); + file.Close(); + +} + +#include "FWCore/Framework/interface/MakerMacros.h" +DEFINE_FWK_MODULE(NewAnalyzer); diff --git a/Validation/Generator/plugins/NewAnalyzer.h b/Validation/Generator/plugins/NewAnalyzer.h new file mode 100644 index 0000000000000..d6ade552745bb --- /dev/null +++ b/Validation/Generator/plugins/NewAnalyzer.h @@ -0,0 +1,39 @@ +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/EDAnalyzer.h" +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" + +#include "HepMC/GenEvent.h" +#include "HepMC/GenParticle.h" +#include "HepMC/WeightContainer.h" + +#include "TH1D.h" +#include "TFile.h" + + +class NewAnalyzer : public edm::EDAnalyzer { + public: + explicit NewAnalyzer(const edm::ParameterSet&); + ~NewAnalyzer(); + + + private: + virtual void beginJob(const edm::EventSetup&) ; + virtual void analyze(const edm::Event&, const edm::EventSetup&); + virtual void endJob() ; + + private: + std::string outputFilename; + TH1D* weight_histo; int event; TH1F* ChDauMult_histo; TH1F* ChDaupT_histo; int tau_evt; + TH1F* Z_invmass_histo; TH1F* Z0J_invmass_histo; TH1F* Z1J_invmass_histo; TH1F* Z2J_invmass_histo; TH1F* Z3J_invmass_histo; TH1F* Z4J_invmass_histo; + TH1F* Z1JJ1Eta_histo; TH1F* Z2JJ1Eta_histo; TH1F* Z3JJ1Eta_histo; TH1F* Z4JJ1Eta_histo; TH1F* ZEta_histo; + TH1F* Z2JJ2Eta_histo; TH1F* Z3JJ2Eta_histo; TH1F* Z4JJ2Eta_histo; TH1F* Jetmult_histo; + TH1F* Pt_histo; TH1F* J1Pt_histo; TH1F* J2Pt_histo; TH1F* EJDelR_histo; + TH1F* E1Pt_histo; TH1F* E2Pt_histo; TH1F* ZPz_histo; TH1F* ZPt_histo; TH1F* Z2JJDelPhi_histo; TH1F* Z3JJDelPhi_histo; TH1F* Z4JJDelPhi_histo; + TH1F* JDelR_histo; TH1F* JDelPhi_histo; TH1F* Z1JJ1Phi_histo; TH1F* Z2JJ1Phi_histo; TH1F* Z3JJ1Phi_histo; TH1F* Z4JJ1Phi_histo; + TH1F* Z2JJ2Phi_histo; TH1F* Z3JJ2Phi_histo; TH1F* Z4JJ2Phi_histo; TH1F* Z2JJDelR_histo; TH1F* Z3JJDelR_histo; TH1F* Z4JJDelR_histo; TH1F* ZRap_histo; + TH1F* JetPt1J; TH1F* JetPt2J; TH1F* JetPt3J; TH1F* JetPt4J; + TH1F* ZPt1J_histo; TH1F* ZPt2J_histo; TH1F* ZPt3J_histo; TH1F* ZPt4J_histo; TH1F* ZPt0J_histo; + TH1F* Z1JJ1Pt_histo; TH1F* Z2JJ1Pt_histo; TH1F* Z2JJ2Pt_histo; TH1F* Z3JJ1Pt_histo; TH1F* Z3JJ2Pt_histo; + TH1F* Z4JJ1Pt_histo; TH1F* Z4JJ2Pt_histo; TH1F* J1Eta_histo; TH1F* J1Phi_histo; +}; diff --git a/Validation/Generator/plugins/PhotonJetAnalyzer.cc b/Validation/Generator/plugins/PhotonJetAnalyzer.cc new file mode 100644 index 0000000000000..ebf63a1606ad1 --- /dev/null +++ b/Validation/Generator/plugins/PhotonJetAnalyzer.cc @@ -0,0 +1,228 @@ +#include "Validation/Generator/plugins/PhotonJetAnalyzer.h" + +using namespace reco; +using namespace std; + +PhotonJetAnalyzer::PhotonJetAnalyzer(const edm::ParameterSet& iConfig) +{ +} + +PhotonJetAnalyzer::~PhotonJetAnalyzer() +{ +} + +void +PhotonJetAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) +{ + using namespace edm; + + // *************************************************************************** + // ***** Collect Generated Particle data ************************************* + + // Get Gen Particle Collection + + + basicStruct mcHardPhoton; + mcHardPhoton.pt = -1.0; + int mcHardPhotonMomID = 0; + bool dtctbleDrctPhtnFnd = false; + + /* + Handle genParticles; + iEvent.getByLabel("genParticles", genParticles); + */ + Handle genParticles; + try { + iEvent.getByLabel( "genParticleCandidates", genParticles ); + } catch (...) { + cout << "No " << "genParticleCandidates" << " found!" << endl; + } + + for(size_t i = 0; (i < genParticles->size()) && (!dtctbleDrctPhtnFnd); ++ i) { + + + const Candidate & p = (*genParticles)[ i ]; + + int id = p.pdgId(); + int status = p.status(); + double et = p.et(); + double eta = p.eta(); + double phi = p.phi(); + + // Find gen photon, with a minimum et of 20 GeV. + // Detectable particles have status of 1. + if ( (id == 22) && (status == 1) && (et > mcHardPhoton.pt) && (et > 20) ) { + + const Candidate * mom = p.mother(); + int motherID = mom->pdgId(); + + if (motherID < 25) { + dtctbleDrctPhtnFnd = true; + mcHardPhoton.et = et; // et and pt are same for photon + mcHardPhoton.phi = phi; + mcHardPhoton.eta = eta; + mcHardPhotonMomID = motherID; + } + + } // End of IF a photon + } + /* + cout << "Photon = " << mcHardPhoton.et << "\t" << mcHardPhoton.eta << "\t" << mcHardPhoton.phi << endl; + */ + // ***** End Collect Generated Pharticle data ******************************** + // *************************************************************************** + + + + // *************************************************************************** + // ***** Collect Gen Jet data ************************************************ + + // get gen jet collection + Handle jetsgen; + try { + iEvent.getByLabel("iterativeCone5GenJetsPt10", jetsgen); + } catch (...) { + cout << "No " << "iterativeCone5GenJetsPt10" << " found!" << endl; + } + + basicStruct hiGenJet; + hiGenJet.pt = -1.0; + + int hiGenJetCount = 0; + + if (jetsgen.isValid() ) { + for (size_t j = 0; j < jetsgen->size(); j++) { + double pt = (*jetsgen)[j].pt(); + double eta = (*jetsgen)[j].eta(); + double phi = (*jetsgen)[j].phi(); + + // Count the number of gen jets above 10 GeV + if (pt > 10) { + hiGenJetCount++; + } + + // Only care about storing jet with reasonable pt, and opposite phi + if ( (pt<10) || (fabs(phi - mcHardPhoton.phi) < 0.8*3.14159) || (fabs(phi - mcHardPhoton.phi) > 1.2*3.14159) ) continue; + + // Record this jet if it's gives the highest pt + if (pt > hiGenJet.pt) { + hiGenJet.pt = pt; + hiGenJet.eta = eta; + hiGenJet.phi = phi; + } + } + } + /* + cout << "Found " << hiGenJetCount << " gen jets." << endl; + cout << "Gen Jet :\t" << hiGenJet.pt << "\t" << hiGenJet.eta << "\t" << hiGenJet.phi << endl; + */ + // ***** END Collect Gen Jet data ******************************************** + // *************************************************************************** + + + + // *************************************************************************** + // ***** Fill Histograms ***************************************************** + if ( (mcHardPhoton.et > 0.0) && (fabs(mcHardPhoton.eta) < 2.5) ) { + cout << "Filling histograms..."; + hGenPhtnHardEt-> Fill(mcHardPhoton.et); + hGenPhtnHardEta->Fill(mcHardPhoton.eta); + hGenPhtnHardPhi->Fill(mcHardPhoton.phi); + hGenPhtnHardMom->Fill(mcHardPhotonMomID); + + // Highest pt gen jet + if (hiGenJet.pt > 0.0) { + hHIGenJetPt-> Fill(hiGenJet.pt); + hHIGenJetEta-> Fill(hiGenJet.eta); + hHIGenJetPhi-> Fill(hiGenJet.phi); + hGenPhtnHardDrJet->Fill(calcDeltaR(mcHardPhoton.eta, mcHardPhoton.phi, hiGenJet.eta, hiGenJet.phi)); + } + hHIGenJetCnt-> Fill(hiGenJetCount); + cout << "...filled." << endl; + } + // ***** END Fill Histograms ************************************************* + // *************************************************************************** +} + +void +PhotonJetAnalyzer::beginJob(const edm::EventSetup&) +{ + + // Plotting variables + const float MAX_E = 800.0; + const float MIN_E = 0.0; + const float MAX_ETA = 3.29867223; + const float MAX_PHI = 3.29867223; + const int HI_NUM_BINS = 80; + const int ME_NUM_BINS = 42; + // const int LO_NUM_BINS = 25; + + // Simple Histograms + hGenPhtnHardEt = new TH1F("genPhtnEt" , "Highest Et Gen Photon: Et " , HI_NUM_BINS, MIN_E, MAX_E); + hGenPhtnHardEta = new TH1F("genPhtnEta", "Highest Et Gen Photon: #eta", ME_NUM_BINS, -MAX_ETA, MAX_ETA); + hGenPhtnHardPhi = new TH1F("genPhtnPhi", "Highest Et Gen Photon: #phi", ME_NUM_BINS, -MAX_PHI, MAX_PHI); + hGenPhtnHardMom = new TH1F("genPhtnMom", "Highest Et Gen Photon: Mother ID", 40, -0.5, 39.5); + hGenPhtnHardDrJet = new TH1F("genPhtnDrJet", "Highest Et Gen Photon: #DeltaR to Gen Jet", ME_NUM_BINS, 0.0, 5.0); + + hHIGenJetPt = new TH1F("genJetPt" , "Highest Pt Gen Jet: Pt " , HI_NUM_BINS, MIN_E, MAX_E); + hHIGenJetEta = new TH1F("genJetEta", "Highest Pt Gen Jet: #eta ", ME_NUM_BINS, -2*MAX_ETA, 2*MAX_ETA); + hHIGenJetPhi = new TH1F("genJetPhi", "Highest Pt Gen Jet: #phi ", ME_NUM_BINS, -MAX_PHI, MAX_PHI); + hHIGenJetCnt = new TH1F("genJetCnt", "Number of Gen Jets (pt > 20 GeV)", 10, -0.5, 9.5); + +} + +void +PhotonJetAnalyzer::endJob() { + + TFile *outputRootFile = new TFile("photonJetoutput.root","RECREATE"); + outputRootFile->cd(); + + hGenPhtnHardEt->Write(); + hGenPhtnHardEta->Write(); + hGenPhtnHardPhi->Write(); + hGenPhtnHardMom->Write(); + hGenPhtnHardDrJet->Write(); + + hHIGenJetPt->Write(); + hHIGenJetEta->Write(); + hHIGenJetPhi->Write(); + hHIGenJetCnt->Write(); + + // Write histograms and ntuples + outputRootFile->Write(); + outputRootFile->Close(); + +} + +// ********************** +// Method for Calculating the delta-r between two things +double PhotonJetAnalyzer::calcDeltaR(double eta1, double phi1, double eta2, double phi2) { + + double deltaEta = eta1 - eta2; + double deltaPhi = calcDeltaPhi(phi1, phi2); + + double deltaR = deltaEta*deltaEta + deltaPhi*deltaPhi; + deltaR = sqrt(deltaR); + + return deltaR; +} // End of calcDeltaR +// ********************** + +// ********************** +// Method for Calculating the delta-phi +double PhotonJetAnalyzer::calcDeltaPhi(double phi1, double phi2) { + + double deltaPhi = phi1 - phi2; + + if (deltaPhi < 0) deltaPhi = -deltaPhi; + + if (deltaPhi > 3.1415926) { + deltaPhi = 2 * 3.1415926 - deltaPhi; + } + + return deltaPhi; +} // End of calcDeltaPhi +// ********************** + +#include "FWCore/Framework/interface/MakerMacros.h" +DEFINE_FWK_MODULE(PhotonJetAnalyzer); diff --git a/Validation/Generator/plugins/PhotonJetAnalyzer.h b/Validation/Generator/plugins/PhotonJetAnalyzer.h new file mode 100644 index 0000000000000..73af18c5853f6 --- /dev/null +++ b/Validation/Generator/plugins/PhotonJetAnalyzer.h @@ -0,0 +1,83 @@ +// -*- C++ -*- +// +// Package: PhotonJetAnalyzer +// Class: PhotonJetAnalyzer +// +/**\class PhotonJetAnalyzer PhotonJetAnalyzer.cc Analysis/PhotonJetAnalyzer/src/PhotonJetAnalyzer.cc + + Description: + + Implementation: + +*/ +// +// Original Author: Mike Anderson +// Created: Tue Apr 22 10:19:02 CDT 2008 +// $Id: PhotonJetAnalyzer.cc,v 1.1 2008/05/27 21:52:15 ksmith Exp $ +// +// + + +// system include files +#include +#include + +// user include files +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/EDAnalyzer.h" +#include "Geometry/Records/interface/IdealGeometryRecord.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" + +#include "FWCore/ParameterSet/interface/ParameterSet.h" + + +#include "DataFormats/HepMCCandidate/interface/GenParticle.h" +#include "DataFormats/JetReco/interface/GenJet.h" + +// ROOT headers +#include "TFile.h" +#include "TH1.h" +#include "TH1F.h" + + +class PhotonJetAnalyzer : public edm::EDAnalyzer { + public: + explicit PhotonJetAnalyzer(const edm::ParameterSet&); + ~PhotonJetAnalyzer(); + + + private: + virtual void beginJob(const edm::EventSetup&) ; + virtual void analyze(const edm::Event&, const edm::EventSetup&); + virtual void endJob() ; + + // ----------member data --------------------------- + + // Basic set of variables to store for physics objects + struct basicStruct { + double energy; + double et; + double pt; + double eta; + double phi; + } ; + + // Functions + double calcDeltaR(double eta1, double phi1, double eta2, double phi2); + double calcDeltaPhi(double phi1, double phi2); + + // ***** Simple Histograms ***** + // Generated Photon + TH1F* hGenPhtnHardEt; + TH1F* hGenPhtnHardEta; + TH1F* hGenPhtnHardPhi; + TH1F* hGenPhtnHardMom; + TH1F* hGenPhtnHardDrJet; + + // Generated Jet + TH1F* hHIGenJetPt; + TH1F* hHIGenJetEta; + TH1F* hHIGenJetPhi; + TH1F* hHIGenJetCnt; +}; diff --git a/Validation/Generator/plugins/TauAnalyzer.cc b/Validation/Generator/plugins/TauAnalyzer.cc new file mode 100644 index 0000000000000..015163d8884e6 --- /dev/null +++ b/Validation/Generator/plugins/TauAnalyzer.cc @@ -0,0 +1,181 @@ +/* Analyize tau decays in various processes using tauola +*/ + +// +// Created: Avto Kharchilava, Feb. 2008 +// Modified: +// + +// system include files +#include +#include + +// user include files +#include "Validation/Generator/plugins/TauAnalyzer.h" + + +#include "DataFormats/Math/interface/LorentzVector.h" +//#include "SimDataFormats/HepMCProduct/interface/HepMCProduct.h" +#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h" + + +TauAnalyzer::TauAnalyzer(const edm::ParameterSet& iConfig) +{ + outputFilename=iConfig.getUntrackedParameter("OutputFilename","dummy.root"); + invmass_histo = new TH1D("invmass_histo","invmass_histo",50,0,500); + pT1_histo = new TH1D("pT1_histo","pT1_histo",50,0,500); + pT2_histo = new TH1D("pT2_histo","pT2_histo",50,0,500); + + h_mcRtau = new TH1D("h_mcRtau","",20,0,1); + h_mcRleptonic = new TH1D("h_mcRleptonic","",20,0,1); + + eventCounter = 0; + mcHadronicTauCounter = 0; + mcVisibleTauCounter = 0; + mcTauPtCutCounter = 0; + +} + + +TauAnalyzer::~TauAnalyzer() +{ + +} + +// ------------ method called to for each event ------------ +void +TauAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) +{ + using namespace edm; + + //bool select = false; + eventCounter++; + + bool lepton = false; + + math::XYZTLorentzVector LvectorTau(0,0,0,0); + math::XYZTLorentzVector LvectorLep(0,0,0,0); + math::XYZTLorentzVector visibleTau(0,0,0,0); + math::XYZTLorentzVector leadingTrack(0,0,0,0); + + Handle mcEventHandle; + try{ + iEvent.getByLabel("source",mcEventHandle); + }catch(...) {;} + + if(mcEventHandle.isValid()){ + const HepMC::GenEvent* mcEvent = mcEventHandle->GetEvent() ; + + HepMC::GenEvent::particle_const_iterator i; + for(i = mcEvent->particles_begin(); i!= mcEvent->particles_end(); i++){ + int id = (*i)->pdg_id(); + + if(abs(id) != 15) continue; + + + int motherId = 0; + if( (*i)->production_vertex() ) { + HepMC::GenVertex::particle_iterator iMother; + for(iMother = (*i)->production_vertex()->particles_begin(HepMC::parents); + iMother!= (*i)->production_vertex()->particles_end(HepMC::parents); iMother++){ + motherId = (*iMother)->pdg_id(); + + std::cout << " tau mother " << motherId << std::endl; + } + } + + if( abs(motherId) != 37 ) continue; + + HepMC::FourVector p4 = (*i)->momentum(); + LvectorTau = math::XYZTLorentzVector(p4.px(),p4.py(),p4.pz(),p4.e()); + visibleTau = math::XYZTLorentzVector(p4.px(),p4.py(),p4.pz(),p4.e()); + + if( (*i)->production_vertex() ) { + HepMC::GenVertex::particle_iterator iChild; + for(iChild = (*i)->production_vertex()->particles_begin(HepMC::descendants); + iChild!= (*i)->production_vertex()->particles_end(HepMC::descendants);iChild++){ + int childId = (*iChild)->pdg_id(); + + std::cout << "tau child id " << childId << std::endl; + + HepMC::FourVector fv = (*iChild)->momentum(); + math::XYZTLorentzVector p(fv.px(),fv.py(),fv.pz(),fv.e()); + + if( abs(childId) == 12 || abs(childId) == 14 || abs(childId) == 16){ + if((*iChild)->status() == 1 && childId*id > 0) { + visibleTau -= p; + } + } + + if( abs(childId) == 11 || abs(childId) == 13 ){ + lepton = true; + LvectorLep = p; + + std::cout << "chiled lepton id " << childId << std::endl; + + } + + if( abs(childId) == 211 ){ // pi+,rho+ + if(p.P() > leadingTrack.P()) leadingTrack = p; + } + + } + } + + + } + + + } + // ??? if(lepton) return select; + if( lepton ) { + + double Rtau_lep = LvectorLep.E()/LvectorTau.E(); + h_mcRleptonic->Fill(Rtau_lep); + } + + mcHadronicTauCounter++; + + // ??? if(visibleTau.Pt() == 0) return select; + mcVisibleTauCounter++; + std::cout << "vis tau px,py,pz,Pt " << visibleTau.Px() << " " << visibleTau.Py() << " " << visibleTau.Pz() << " " << visibleTau.Pt() << std::endl; +//std::cout << " eta,phi " << visibleTau.Eta() << " " << visibleTau.Phi() << std::endl; + + +// ??? if(visibleTau.Pt() < 100) return select; + mcTauPtCutCounter++; + + std::cout << "visible tau pt " << visibleTau.Pt() << std::endl; + + if(!lepton && visibleTau.Pt() > 100.) { + + double Rtau = leadingTrack.P()/visibleTau.E(); + h_mcRtau->Fill(Rtau); + std::cout << "check Rtau " << Rtau << " " << leadingTrack.P() << " " << visibleTau.E() << std::endl; + } + +} + + +// ------------ method called once each job just before starting event loop ------------ +void +TauAnalyzer::beginJob(const edm::EventSetup&) +{ +} + +// ------------ method called once each job just after ending the event loop ------------ +void +TauAnalyzer::endJob() { + // save histograms into file + TFile file(outputFilename.c_str(),"RECREATE"); + invmass_histo->Write(); + pT1_histo->Write(); + pT2_histo->Write(); + h_mcRtau->Write(); + h_mcRleptonic->Write(); + + file.Close(); +} + +#include "FWCore/Framework/interface/MakerMacros.h" +DEFINE_FWK_MODULE(TauAnalyzer); diff --git a/Validation/Generator/plugins/TauAnalyzer.h b/Validation/Generator/plugins/TauAnalyzer.h new file mode 100644 index 0000000000000..bb53f4cef13c3 --- /dev/null +++ b/Validation/Generator/plugins/TauAnalyzer.h @@ -0,0 +1,37 @@ +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/EDAnalyzer.h" +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" + +#include "HepMC/GenEvent.h" +#include "HepMC/GenParticle.h" +#include "HepMC/WeightContainer.h" + +#include "TH1D.h" +#include "TFile.h" + + +class TauAnalyzer : public edm::EDAnalyzer { + public: + explicit TauAnalyzer(const edm::ParameterSet&); + ~TauAnalyzer(); + + private: + virtual void beginJob(const edm::EventSetup&) ; + virtual void analyze(const edm::Event&, const edm::EventSetup&); + virtual void endJob() ; + + private: + std::string outputFilename; + TH1D* weight_histo; + TH1D* invmass_histo; + TH1D* pT1_histo; + TH1D* pT2_histo; + TH1D* h_mcRtau; + TH1D* h_mcRleptonic; + + int eventCounter; + int mcHadronicTauCounter; + int mcVisibleTauCounter; + int mcTauPtCutCounter; +};