| 1 |
bbetchar |
1.1 |
#ifndef TUPLE_GEN
|
| 2 |
|
|
#define TUPLE_GEN
|
| 3 |
|
|
|
| 4 |
|
|
#include "FWCore/Framework/interface/EDProducer.h"
|
| 5 |
|
|
#include "FWCore/Framework/interface/Frameworkfwd.h"
|
| 6 |
|
|
#include "FWCore/Utilities/interface/InputTag.h"
|
| 7 |
|
|
#include "FWCore/Framework/interface/Event.h"
|
| 8 |
|
|
#include "FWCore/Framework/interface/ESHandle.h"
|
| 9 |
|
|
#include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
|
| 10 |
|
|
#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
|
| 11 |
|
|
|
| 12 |
|
|
// LHE Event headers
|
| 13 |
|
|
#include "SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h"
|
| 14 |
|
|
#include "GeneratorInterface/LHEInterface/interface/LHEEvent.h"
|
| 15 |
|
|
|
| 16 |
|
|
#include <map>
|
| 17 |
|
|
|
| 18 |
|
|
template< typename T >
|
| 19 |
|
|
class Tuple_Gen : public edm::EDProducer {
|
| 20 |
|
|
public:
|
| 21 |
|
|
explicit Tuple_Gen(const edm::ParameterSet&);
|
| 22 |
|
|
private:
|
| 23 |
|
|
void produce(edm::Event &, const edm::EventSetup & );
|
| 24 |
|
|
void produceGenJets(edm::Event &);
|
| 25 |
|
|
static int index(const reco::Candidate*, const std::vector<const T*>&);
|
| 26 |
|
|
typedef reco::Candidate::LorentzVector LorentzVector;
|
| 27 |
|
|
const edm::InputTag inputTag;
|
| 28 |
|
|
const std::vector<edm::InputTag> jetCollections;
|
| 29 |
|
|
const std::string Prefix,Suffix;
|
| 30 |
|
|
const double GenStatus1PtCut;
|
| 31 |
|
|
const double GenJetPtCut;
|
| 32 |
|
|
};
|
| 33 |
|
|
|
| 34 |
|
|
template< typename T > Tuple_Gen<T>::
|
| 35 |
|
|
Tuple_Gen(const edm::ParameterSet& iConfig) :
|
| 36 |
|
|
inputTag(iConfig.getParameter<edm::InputTag>("InputTag")),
|
| 37 |
|
|
jetCollections(iConfig.getParameter<std::vector<edm::InputTag> >("JetCollections")),
|
| 38 |
|
|
Prefix(iConfig.getParameter<std::string>("Prefix")),
|
| 39 |
|
|
Suffix(iConfig.getParameter<std::string>("Suffix")),
|
| 40 |
|
|
GenStatus1PtCut(iConfig.getParameter<double>("GenStatus1PtCut")),
|
| 41 |
|
|
GenJetPtCut(iConfig.getParameter<double>("GenJetPtCut")) {
|
| 42 |
|
|
produces <unsigned int> (Prefix + "signalProcessID" + Suffix);
|
| 43 |
|
|
produces <bool> (Prefix + "GenInfoHandleValid" + Suffix);
|
| 44 |
|
|
produces <bool > (Prefix + "HandleValid" + Suffix);
|
| 45 |
|
|
produces <double> (Prefix + "pthat" + Suffix);
|
| 46 |
|
|
produces <int> (Prefix + "id1" + Suffix);
|
| 47 |
|
|
produces <int> (Prefix + "id2" + Suffix);
|
| 48 |
|
|
produces <double> (Prefix + "x1" + Suffix);
|
| 49 |
|
|
produces <double> (Prefix + "x2" + Suffix);
|
| 50 |
|
|
produces <double> (Prefix + "pdf1" + Suffix);
|
| 51 |
|
|
produces <double> (Prefix + "pdf2" + Suffix);
|
| 52 |
|
|
produces <double> (Prefix + "PartonHT" + Suffix);
|
| 53 |
|
|
produces <std::vector<double> > (Prefix + "cteq66" + Suffix);
|
| 54 |
|
|
produces <std::vector<double> > (Prefix + "NNPDF10" + Suffix);
|
| 55 |
|
|
produces <std::vector<double> > (Prefix + "MRST2006nnlo" + Suffix);
|
| 56 |
|
|
produces <std::vector<double> > (Prefix + "BinningValues" + Suffix);
|
| 57 |
|
|
produces <float> (Prefix + "Q" + Suffix);
|
| 58 |
|
|
produces <std::vector<LorentzVector> > ( Prefix + "P4" + Suffix );
|
| 59 |
|
|
produces <std::vector<int> > (Prefix + "PdgId" + Suffix);
|
| 60 |
|
|
produces <std::vector<int> > (Prefix + "Status" + Suffix);
|
| 61 |
|
|
produces <std::vector<int> > (Prefix + "MotherIndex" + Suffix);
|
| 62 |
|
|
produces <std::vector<int> > (Prefix + "MotherPdgId" + Suffix);
|
| 63 |
|
|
|
| 64 |
|
|
|
| 65 |
|
|
for(unsigned i=0; i<jetCollections.size(); ++i)
|
| 66 |
|
|
produces<std::vector<LorentzVector> >(Prefix + jetCollections[i].label() +"P4" + Suffix);
|
| 67 |
|
|
}
|
| 68 |
|
|
|
| 69 |
|
|
template< typename T > int Tuple_Gen<T>::
|
| 70 |
|
|
index(const reco::Candidate* item, const typename std::vector<const T*>& collection) {
|
| 71 |
|
|
typename std::vector<const T*>::const_iterator it(collection.begin()), begin(collection.begin()), end(collection.end());
|
| 72 |
|
|
for(; it!=end; it++) if ((*it)==item) return it-begin; //Compare addresses
|
| 73 |
|
|
return -1;
|
| 74 |
|
|
}
|
| 75 |
|
|
|
| 76 |
|
|
template< typename T > void Tuple_Gen<T>::
|
| 77 |
|
|
produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
|
| 78 |
|
|
produceGenJets(iEvent);
|
| 79 |
|
|
|
| 80 |
|
|
edm::Handle<std::vector<T> > collection; iEvent.getByLabel(inputTag,collection);
|
| 81 |
|
|
edm::Handle<GenEventInfoProduct> geninfo; iEvent.getByLabel("generator",geninfo);
|
| 82 |
|
|
edm::Handle<std::vector<double> > cteqHandle;
|
| 83 |
|
|
edm::Handle<std::vector<double> > nnpdfHandle;
|
| 84 |
|
|
edm::Handle<std::vector<double> > mrstHandle;
|
| 85 |
|
|
iEvent.getByLabel("pdfWeights", "cteq66", cteqHandle);
|
| 86 |
|
|
iEvent.getByLabel("pdfWeights", "NNPDF10", nnpdfHandle);
|
| 87 |
|
|
iEvent.getByLabel("pdfWeights", "MRST2006nnlo", mrstHandle);
|
| 88 |
|
|
|
| 89 |
|
|
//add handle for LHE event
|
| 90 |
|
|
edm::Handle<LHEEventProduct> product;
|
| 91 |
|
|
iEvent.getByLabel("source", product);
|
| 92 |
|
|
|
| 93 |
|
|
std::auto_ptr<unsigned int> signalProcessID(new unsigned int(geninfo->signalProcessID()));
|
| 94 |
|
|
std::auto_ptr<float> Q (new float(geninfo->pdf()->scalePDF));
|
| 95 |
|
|
std::auto_ptr<int> id1 (new int( geninfo->pdf()->id.first));
|
| 96 |
|
|
std::auto_ptr<int> id2 (new int( geninfo->pdf()->id.second));
|
| 97 |
|
|
std::auto_ptr<double> x1 (new double( geninfo->pdf()->x.first));
|
| 98 |
|
|
std::auto_ptr<double> x2 (new double( geninfo->pdf()->x.second));
|
| 99 |
|
|
std::auto_ptr<double> pdf1 (new double( geninfo->pdf()->xPDF.first));
|
| 100 |
|
|
std::auto_ptr<double> pdf2 (new double( geninfo->pdf()->xPDF.second));
|
| 101 |
|
|
|
| 102 |
|
|
std::auto_ptr<bool> handleValid ( new bool(collection.isValid()) );
|
| 103 |
|
|
std::auto_ptr<bool> genInfoValid ( new bool( geninfo.isValid() && !geninfo->binningValues().empty()));
|
| 104 |
|
|
std::auto_ptr<double> pthat (new double(*genInfoValid ? geninfo->binningValues()[0] : -1.));
|
| 105 |
|
|
std::auto_ptr<std::vector<double> > binningValues (*genInfoValid ? new std::vector<double>(geninfo->binningValues()) : new std::vector<double>());
|
| 106 |
|
|
std::auto_ptr<std::vector<LorentzVector> > p4 ( new std::vector<LorentzVector>() ) ;
|
| 107 |
|
|
std::auto_ptr<std::vector<int> > status ( new std::vector<int>() ) ;
|
| 108 |
|
|
std::auto_ptr<std::vector<int> > pdgId ( new std::vector<int>() ) ;
|
| 109 |
|
|
std::auto_ptr<std::vector<int> > motherIndex ( new std::vector<int>() ) ;
|
| 110 |
|
|
std::auto_ptr<std::vector<int> > motherPdgId ( new std::vector<int>() ) ;
|
| 111 |
|
|
std::auto_ptr<std::vector<double> > cteq ( new std::vector<double>() ) ;
|
| 112 |
|
|
std::auto_ptr<std::vector<double> > nnpdf ( new std::vector<double>() ) ;
|
| 113 |
|
|
std::auto_ptr<std::vector<double> > mrst ( new std::vector<double>() ) ;
|
| 114 |
|
|
std::auto_ptr<double> gHT ( new double(0) );
|
| 115 |
|
|
|
| 116 |
|
|
std::vector<const T*> self;
|
| 117 |
|
|
std::vector<const reco::Candidate*> mom;
|
| 118 |
|
|
|
| 119 |
|
|
if(cteqHandle.isValid()){
|
| 120 |
|
|
for(std::vector<double>::const_iterator it = cteqHandle->begin(); it != cteqHandle->end(); ++it) {
|
| 121 |
|
|
cteq->push_back(*it);
|
| 122 |
|
|
}}
|
| 123 |
|
|
if(nnpdfHandle.isValid()){
|
| 124 |
|
|
for(std::vector<double>::const_iterator it = nnpdfHandle->begin(); it != nnpdfHandle->end(); ++it) {
|
| 125 |
|
|
nnpdf->push_back(*it);
|
| 126 |
|
|
}}
|
| 127 |
|
|
if(mrstHandle.isValid()){
|
| 128 |
|
|
for(std::vector<double>::const_iterator it = mrstHandle->begin(); it != mrstHandle->end(); ++it) {
|
| 129 |
|
|
mrst->push_back(*it);
|
| 130 |
|
|
}}
|
| 131 |
|
|
|
| 132 |
|
|
|
| 133 |
|
|
|
| 134 |
|
|
if(collection.isValid()){
|
| 135 |
|
|
for(typename std::vector<T>::const_iterator it = collection->begin(); it != collection->end(); ++it) {
|
| 136 |
|
|
if ( it->status() == 3 // any status 3 genParticle
|
| 137 |
|
|
|| abs(it->pdgId()) == 11 // any electron
|
| 138 |
|
|
|| abs(it->pdgId()) == 13 // any muon
|
| 139 |
|
|
|| abs(it->pdgId()) == 15 // any tau
|
| 140 |
|
|
|| ( it->status() == 1 // status 1 particles
|
| 141 |
|
|
&& it->pt() > GenStatus1PtCut) // above threshold
|
| 142 |
|
|
) {
|
| 143 |
|
|
p4->push_back(it->p4());
|
| 144 |
|
|
status->push_back(it->status());
|
| 145 |
|
|
pdgId->push_back(it->pdgId());
|
| 146 |
|
|
motherPdgId->push_back( it->numberOfMothers() ? it->mother()->pdgId() : 0 );
|
| 147 |
|
|
self.push_back(&*it);
|
| 148 |
|
|
mom.push_back( it->numberOfMothers() ? it->mother(): 0);
|
| 149 |
|
|
}
|
| 150 |
|
|
}
|
| 151 |
|
|
} //collection
|
| 152 |
|
|
|
| 153 |
|
|
// taken from the discussion here: https://hypernews.cern.ch/HyperNews/CMS/get/generators/1234/1/1/1.html
|
| 154 |
|
|
|
| 155 |
|
|
if(product.isValid()){
|
| 156 |
|
|
//gen Parton HT
|
| 157 |
|
|
|
| 158 |
|
|
//make some new LHE based variables
|
| 159 |
|
|
const lhef::HEPEUP hepeup_ = product->hepeup();
|
| 160 |
|
|
const std::vector<lhef::HEPEUP::FiveVector> pup_ = hepeup_.PUP;
|
| 161 |
|
|
|
| 162 |
|
|
double htEvent=0.;
|
| 163 |
|
|
size_t iMax = hepeup_.NUP;
|
| 164 |
|
|
for(size_t i = 2; i < iMax; ++i) {
|
| 165 |
|
|
if( hepeup_.ISTUP[i] != 1 ) continue;
|
| 166 |
|
|
int idabs = abs( hepeup_.IDUP[i] );
|
| 167 |
|
|
if( idabs != 21 && (idabs<1 || idabs>6) ) continue;
|
| 168 |
|
|
double ptPart = sqrt( pow(hepeup_.PUP[i][0],2) + pow(hepeup_.PUP[i][1],2) );
|
| 169 |
|
|
//std::cout << ">>>>>>>> Pt Parton: " << ptPart << std::endl;
|
| 170 |
|
|
htEvent += ptPart;
|
| 171 |
|
|
} // iMax
|
| 172 |
|
|
|
| 173 |
|
|
*gHT=htEvent;
|
| 174 |
|
|
} //product
|
| 175 |
|
|
|
| 176 |
|
|
for(typename std::vector<const reco::Candidate*>::const_iterator it = mom.begin(); it!=mom.end(); ++it)
|
| 177 |
|
|
motherIndex->push_back( index(*it,self) );
|
| 178 |
|
|
|
| 179 |
|
|
iEvent.put( handleValid, Prefix + "HandleValid" + Suffix);
|
| 180 |
|
|
iEvent.put( genInfoValid, Prefix + "GenInfoHandleValid" + Suffix);
|
| 181 |
|
|
iEvent.put( pthat, Prefix + "pthat" + Suffix);
|
| 182 |
|
|
iEvent.put( binningValues,Prefix + "BinningValues" + Suffix);
|
| 183 |
|
|
iEvent.put( p4, Prefix + "P4" + Suffix );
|
| 184 |
|
|
iEvent.put( status, Prefix + "Status" + Suffix );
|
| 185 |
|
|
iEvent.put( pdgId, Prefix + "PdgId" + Suffix );
|
| 186 |
|
|
iEvent.put( motherIndex, Prefix + "MotherIndex" + Suffix );
|
| 187 |
|
|
iEvent.put( motherPdgId, Prefix + "MotherPdgId" + Suffix );
|
| 188 |
|
|
iEvent.put( signalProcessID, Prefix + "signalProcessID" + Suffix );
|
| 189 |
|
|
iEvent.put( Q, Prefix + "Q" + Suffix );
|
| 190 |
|
|
iEvent.put( x1, Prefix + "x1" + Suffix );
|
| 191 |
|
|
iEvent.put( x2, Prefix + "x2" + Suffix );
|
| 192 |
|
|
iEvent.put( pdf1, Prefix + "pdf1" + Suffix );
|
| 193 |
|
|
iEvent.put( pdf2, Prefix + "pdf2" + Suffix );
|
| 194 |
|
|
iEvent.put( id1, Prefix + "id1" + Suffix );
|
| 195 |
|
|
iEvent.put( id2, Prefix + "id2" + Suffix );
|
| 196 |
|
|
iEvent.put( gHT, Prefix + "PartonHT" + Suffix );
|
| 197 |
|
|
iEvent.put( cteq, Prefix + "cteq66" + Suffix);
|
| 198 |
|
|
iEvent.put( nnpdf, Prefix + "NNPDF10" + Suffix);
|
| 199 |
|
|
iEvent.put( mrst, Prefix + "MRST2006nnlo" + Suffix);
|
| 200 |
|
|
|
| 201 |
|
|
}
|
| 202 |
|
|
|
| 203 |
|
|
template< typename T > void Tuple_Gen<T>::
|
| 204 |
|
|
produceGenJets(edm::Event& iEvent) {
|
| 205 |
|
|
for(unsigned i=0; i<jetCollections.size(); ++i) {
|
| 206 |
|
|
std::auto_ptr<std::vector<LorentzVector> > p4(new std::vector<LorentzVector>());
|
| 207 |
|
|
edm::Handle<edm::View<reco::GenJet> > genjets;
|
| 208 |
|
|
iEvent.getByLabel(jetCollections[i], genjets);
|
| 209 |
|
|
if(genjets.isValid())
|
| 210 |
|
|
for(edm::View<reco::GenJet>::const_iterator it(genjets->begin()), end(genjets->end()); it!=end; ++it) {
|
| 211 |
|
|
if (it->pt() >= GenJetPtCut) p4->push_back(it->p4());
|
| 212 |
|
|
}
|
| 213 |
|
|
iEvent.put(p4, Prefix + jetCollections[i].label() + "P4" + Suffix);
|
| 214 |
|
|
}
|
| 215 |
|
|
}
|
| 216 |
|
|
|
| 217 |
|
|
#endif
|
