diff --git a/Validation/Generator/src/BBbarAnalyzer.cc b/Validation/Generator/src/BBbarAnalyzer.cc deleted file mode 100644 index ffe72f993ee16..0000000000000 --- a/Validation/Generator/src/BBbarAnalyzer.cc +++ /dev/null @@ -1,118 +0,0 @@ -/* 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.2 2008/07/01 02:55:39 ksmith Exp $ -// -// - - -// system include files -#include -#include - -// user include files -#include "BBbarAnalyzer.h" - - - - - -#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h" -#include "DataFormats/Math/interface/LorentzVector.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() -{ - -} - -// ------------ method called to for each event ------------ -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); - } - } - } -} - - -// ------------ method called once each job just before starting event loop ------------ -void -BBbarAnalyzer::beginJob(const edm::EventSetup&) -{ -} - -// ------------ method called once each job just after ending the event loop ------------ -void -BBbarAnalyzer::endJob() { - // save histograms into file - TFile file(outputFilename.c_str(),"RECREATE"); - Pt_histo->Write(); - invmass_histo->Write(); - file.Close(); - -} - -//define this as a plug-in -DEFINE_FWK_MODULE(BBbarAnalyzer); diff --git a/Validation/Generator/src/H4muAnalyzer.cc b/Validation/Generator/src/H4muAnalyzer.cc deleted file mode 100644 index 1dd5cdcc7e9a3..0000000000000 --- a/Validation/Generator/src/H4muAnalyzer.cc +++ /dev/null @@ -1,112 +0,0 @@ -/* 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.2 2008/07/01 02:55:39 ksmith Exp $ -// -// - - -// system include files -#include -#include - -// user include files -#include "H4muAnalyzer.h" - - - - - -#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h" -#include "DataFormats/Math/interface/LorentzVector.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() -{ - -} - -// ------------ method called to for each event ------------ -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); - -} - - -// ------------ method called once each job just before starting event loop ------------ -void -H4muAnalyzer::beginJob(const edm::EventSetup&) -{ -} - -// ------------ method called once each job just after ending the event loop ------------ -void -H4muAnalyzer::endJob() { - // save histograms into file - TFile file(outputFilename.c_str(),"RECREATE"); - invmass_histo->Write(); - Z_invmass_histo->Write(); - file.Close(); - -} - -//define this as a plug-in -DEFINE_FWK_MODULE(H4muAnalyzer); diff --git a/Validation/Generator/src/H4muAnalyzer.h b/Validation/Generator/src/H4muAnalyzer.h deleted file mode 100644 index ee8b578062860..0000000000000 --- a/Validation/Generator/src/H4muAnalyzer.h +++ /dev/null @@ -1,45 +0,0 @@ -#ifndef H4MU_ANALYZER -#define H4MU_ANALYZER - -#include "FWCore/Framework/interface/Frameworkfwd.h" -#include "FWCore/Framework/interface/EDAnalyzer.h" - -#include "FWCore/Framework/interface/Event.h" -#include "FWCore/Framework/interface/MakerMacros.h" - -#include "FWCore/ParameterSet/interface/ParameterSet.h" - -#include "HepMC/WeightContainer.h" -#include "HepMC/GenEvent.h" -#include "HepMC/GenParticle.h" - -//#include "SimDataFormats/HepMCProduct/interface/HepMCProduct.h" - -#include "TH1D.h" -#include "TFile.h" - -// -// class decleration -// - -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() ; - - // ----------member data --------------------------- - - std::string outputFilename; - TH1D* weight_histo; - TH1D* invmass_histo; - TH1D* Z_invmass_histo; - -}; - -#endif diff --git a/Validation/Generator/src/HZZ4muAnalyzer.cc b/Validation/Generator/src/HZZ4muAnalyzer.cc deleted file mode 100755 index bda3f7042f192..0000000000000 --- a/Validation/Generator/src/HZZ4muAnalyzer.cc +++ /dev/null @@ -1,219 +0,0 @@ -#include - -#include "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 "FWCore/Framework/interface/MakerMacros.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 ; -} - -DEFINE_FWK_MODULE(HZZ4muAnalyzer); diff --git a/Validation/Generator/src/HZZ4muAnalyzer.h b/Validation/Generator/src/HZZ4muAnalyzer.h deleted file mode 100755 index 8cdf5dcb3a2c4..0000000000000 --- a/Validation/Generator/src/HZZ4muAnalyzer.h +++ /dev/null @@ -1,35 +0,0 @@ -#ifndef HZZ4muAnalyzer_H -#define HZZ4muAnalyzer_H - -#include "FWCore/Framework/interface/EDAnalyzer.h" - -// forward declarations -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 ; - -}; - -#endif diff --git a/Validation/Generator/src/MCValidation.cc b/Validation/Generator/src/MCValidation.cc deleted file mode 100644 index 87a4571d29b2c..0000000000000 --- a/Validation/Generator/src/MCValidation.cc +++ /dev/null @@ -1,551 +0,0 @@ -// -*- 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.1 2008/07/01 22:21:07 ksmith 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/MakerMacros.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" - -using namespace std ; -using namespace edm ; -using namespace reco; - -// -// class decleration -// - -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 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 ; - -}; - -// -// constants, enums and typedefs -// - -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 ; - -// -// static data member definitions -// - -// -// constructors and destructor -// -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.) - -} - - -// -// member functions -// - -// ------------ method called to for each event ------------ -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 - -// ------------ method called once each job just before starting event loop ------------ -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 ; - -} - -// ------------ method called once each job just after ending the event loop ------------ -void -MCValidation::endJob() { - -} - -//define this as a plug-in -DEFINE_FWK_MODULE(MCValidation); - diff --git a/Validation/Generator/src/MuMuAnalyzer.cc b/Validation/Generator/src/MuMuAnalyzer.cc deleted file mode 100644 index 4eaf7b2f9db1a..0000000000000 --- a/Validation/Generator/src/MuMuAnalyzer.cc +++ /dev/null @@ -1,285 +0,0 @@ -/* 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.2 2009/04/27 21:15:22 ksmith Exp $ -// -// - - -// system include files -#include -#include - -// user include files -#include "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() -{ - -} - -// ------------ method called to for each event ------------ -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; -} - - - -// ------------ method called once each job just before starting event loop ------------ -void -MuMuAnalyzer::beginJob(const edm::EventSetup&) -{ -} - -// ------------ method called once each job just after ending the event loop ------------ -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(); - -} - -//define this as a plug-in -DEFINE_FWK_MODULE(MuMuAnalyzer); diff --git a/Validation/Generator/src/MuMuAnalyzer.h b/Validation/Generator/src/MuMuAnalyzer.h deleted file mode 100644 index 57ac342bd868b..0000000000000 --- a/Validation/Generator/src/MuMuAnalyzer.h +++ /dev/null @@ -1,46 +0,0 @@ -#ifndef H4MU_ANALYZER -#define H4MU_ANALYZER - -#include "FWCore/Framework/interface/Frameworkfwd.h" -#include "FWCore/Framework/interface/EDAnalyzer.h" - -#include "FWCore/Framework/interface/Event.h" -#include "FWCore/Framework/interface/MakerMacros.h" - -#include "FWCore/ParameterSet/interface/ParameterSet.h" - -#include "HepMC/WeightContainer.h" -#include "HepMC/GenEvent.h" -#include "HepMC/GenParticle.h" - -//#include "SimDataFormats/HepMCProduct/interface/HepMCProduct.h" - -#include "TH1D.h" -#include "TFile.h" - -// -// class decleration -// - -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() ; - - // ----------member data --------------------------- - - 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; -}; - -#endif diff --git a/Validation/Generator/src/NewAnalyzer.cc b/Validation/Generator/src/NewAnalyzer.cc deleted file mode 100644 index 9b5f3a331fc48..0000000000000 --- a/Validation/Generator/src/NewAnalyzer.cc +++ /dev/null @@ -1,586 +0,0 @@ -/* 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.5 2009/04/27 21:15:22 ksmith Exp $ -// -// - - -// system include files -#include -#include - -// user include files -#include "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() -{ - -} - -// ------------ method called to for each event ------------ -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; -} - - - -// ------------ method called once each job just before starting event loop ------------ -void -NewAnalyzer::beginJob(const edm::EventSetup&) -{ -} - -// ------------ method called once each job just after ending the event loop ------------ -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(); - -} - -//define this as a plug-in -DEFINE_FWK_MODULE(NewAnalyzer); diff --git a/Validation/Generator/src/NewAnalyzer.h b/Validation/Generator/src/NewAnalyzer.h deleted file mode 100644 index eede22226e37b..0000000000000 --- a/Validation/Generator/src/NewAnalyzer.h +++ /dev/null @@ -1,53 +0,0 @@ -#ifndef H4MU_ANALYZER -#define H4MU_ANALYZER - -#include "FWCore/Framework/interface/Frameworkfwd.h" -#include "FWCore/Framework/interface/EDAnalyzer.h" - -#include "FWCore/Framework/interface/Event.h" -#include "FWCore/Framework/interface/MakerMacros.h" - -#include "FWCore/ParameterSet/interface/ParameterSet.h" - -#include "HepMC/WeightContainer.h" -#include "HepMC/GenEvent.h" -#include "HepMC/GenParticle.h" - -//#include "SimDataFormats/HepMCProduct/interface/HepMCProduct.h" - -#include "TH1D.h" -#include "TFile.h" - -// -// class decleration -// - -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() ; - - // ----------member data --------------------------- - - 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; -}; - -#endif diff --git a/Validation/Generator/src/PhotonJetAnalyzer.cc b/Validation/Generator/src/PhotonJetAnalyzer.cc deleted file mode 100644 index bf8fae6967506..0000000000000 --- a/Validation/Generator/src/PhotonJetAnalyzer.cc +++ /dev/null @@ -1,341 +0,0 @@ -// -*- 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/23 17:02:43 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" - -using namespace reco; -using namespace std; -// -// class decleration -// - -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; -}; - -// -// constants, enums and typedefs -// - -// -// static data member definitions -// - -// -// constructors and destructor -// -PhotonJetAnalyzer::PhotonJetAnalyzer(const edm::ParameterSet& iConfig) - -{ - //now do what ever initialization is needed - -} - - -PhotonJetAnalyzer::~PhotonJetAnalyzer() -{ - - // do anything here that needs to be done at desctruction time - // (e.g. close files, deallocate resources etc.) - -} - - -// -// member functions -// - -// ------------ method called to for each event ------------ -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 ************************************************* - // *************************************************************************** -} - - -// ------------ method called once each job just before starting event loop ------------ -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); - -} - -// ------------ method called once each job just after ending the event loop ------------ -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 -// ********************** - -//define this as a plug-in -DEFINE_FWK_MODULE(PhotonJetAnalyzer); diff --git a/Validation/Generator/src/TauAnalyzer.cc b/Validation/Generator/src/TauAnalyzer.cc deleted file mode 100644 index de95bdef64216..0000000000000 --- a/Validation/Generator/src/TauAnalyzer.cc +++ /dev/null @@ -1,197 +0,0 @@ -/* Analyize tau decays in various processes using tauola -*/ - -// -// Created: Avto Kharchilava, Feb. 2008 -// Modified: -// - -// system include files -#include -#include - -// user include files -#include "TauAnalyzer.h" - - -#include "FWCore/Framework/interface/Frameworkfwd.h" -#include "FWCore/Framework/interface/EDAnalyzer.h" - -#include "FWCore/Framework/interface/Event.h" -#include "FWCore/Framework/interface/MakerMacros.h" - -#include "FWCore/ParameterSet/interface/ParameterSet.h" - -#include "HepMC/GenEvent.h" -#include "HepMC/GenParticle.h" - -#include "DataFormats/Math/interface/LorentzVector.h" - -//#include "SimDataFormats/HepMCProduct/interface/HepMCProduct.h" -#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h" - -#include "TH1D.h" -#include "TFile.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(); -} - -//define this as a plug-in -#include "FWCore/Framework/interface/MakerMacros.h" -DEFINE_FWK_MODULE(TauAnalyzer); diff --git a/Validation/Generator/src/TauAnalyzer.h b/Validation/Generator/src/TauAnalyzer.h deleted file mode 100644 index e5146e97a6bb1..0000000000000 --- a/Validation/Generator/src/TauAnalyzer.h +++ /dev/null @@ -1,53 +0,0 @@ -#ifndef TAU_ANALYZER -#define TAU_ANALYZER - -#include "FWCore/Framework/interface/Frameworkfwd.h" -#include "FWCore/Framework/interface/EDAnalyzer.h" - -#include "FWCore/Framework/interface/Event.h" -#include "FWCore/Framework/interface/MakerMacros.h" - -#include "FWCore/ParameterSet/interface/ParameterSet.h" - -#include "HepMC/WeightContainer.h" -#include "HepMC/GenEvent.h" -#include "HepMC/GenParticle.h" - -//#include "SimDataFormats/HepMCProduct/interface/HepMCProduct.h" - -#include "TH1D.h" -#include "TFile.h" - -// -// class decleration -// - -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() ; - - // ----------member data --------------------------- - - 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; - -}; - -#endif