Page MenuHomeHEPForge
Files F9501482

No OneTemporary
Actions

View Options

diff --git a/Analysis/HJetsAnalysis2.cc b/Analysis/HJetsAnalysis2.cc
--- a/Analysis/HJetsAnalysis2.cc
+++ b/Analysis/HJetsAnalysis2.cc
@@ -1,754 +1,782 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the HJetsAnalysis2 class.
//
#include "HJetsAnalysis2.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Handlers/EventHandler.h"
#include "ThePEG/Handlers/StandardEventHandler.h"
#include "Herwig++/Sampling/GeneralSampler.h"
#include "Herwig++/Utilities/XML/ElementIO.h"
using namespace HJets;
using namespace Herwig;
using namespace Statistics;
HJetsAnalysis2::HJetsAnalysis2()
- : theFullEvent(false), theBinFactor(1.0), theFuzzy(true), theEnergyCutWidth(1.0*GeV), theRapidityCutWidth(0.1),theRapidityGap(0.0), theopposite_dir(false), thehiggs_ingap(false), themjj_min(0.0*GeV),thehiggsPt_min(0.0*GeV), thehiggsY_max(1000.0) {}
+ : theFullEvent(false), theBinFactor(1.0), theFuzzy(true), theEnergyCutWidth(1.0*GeV), theRapidityCutWidth(0.1),theRapidityGap(0.0),
+ theopposite_dir(false), thenjets_min(2), thehiggs_ingap(false), themjj_min(0.0*GeV),thehiggsPt_min(0.0*GeV), thehiggsY_max(1000.0) {}
HJetsAnalysis2::~HJetsAnalysis2() {}
#ifndef LWH_AIAnalysisFactory_H
#ifndef LWH
#define LWH ThePEGLWH
#endif
#include "ThePEG/Analysis/LWH/AnalysisFactory.h"
#endif
inline double deltaPhi(const LorentzMomentum& a,
const LorentzMomentum& b){
double phi1 = a.phi();
double phi2 = b.phi();
double diff=phi1-phi2;
if(diff<-Constants::pi){
diff+=(2.0*Constants::pi);
}
else if (diff>Constants::pi){
diff-=(2.0*Constants::pi);
}
return diff;
}
inline double deltaY(const LorentzMomentum& a,
const LorentzMomentum& b){
return abs(a.rapidity()-b.rapidity());
}
inline double deltaR(const LorentzMomentum& a,
const LorentzMomentum& b){
return sqrt(sqr(deltaPhi(a,b))+sqr(deltaY(a,b)));
}
inline double mass(const LorentzMomentum& a,
const LorentzMomentum& b){
return (a+b).m()/GeV;
}
inline double mass(const LorentzMomentum& a,
const LorentzMomentum& b, const LorentzMomentum& c){
return (a+b+c).m()/GeV;
}
void HJetsAnalysis2::analyze(tEventPtr event, long ieve, int loop, int state) {
AnalysisHandler::analyze(event, ieve, loop, state);
if ( !theFullEvent ) {
tSubProPtr sub = event->primarySubProcess();
Ptr<SubProcessGroup>::tptr grp =
dynamic_ptr_cast<Ptr<SubProcessGroup>::tptr>(sub);
analyze(sub->outgoing(),ieve,event->weight()*sub->groupWeight());
if ( grp ) {
for ( SubProcessVector::const_iterator s = grp->dependent().begin();
s != grp->dependent().end(); ++s ) {
analyze((**s).outgoing(),ieve,event->weight()*(**s).groupWeight());
}
}
} else {
ParticleVector fs;
event->getFinalState(fs);
analyze(fs,ieve,event->weight());
}
}
struct SortPt {
inline bool operator()(const LorentzMomentum& a,
const LorentzMomentum& b) const {
return a.perp() > b.perp();
}
};
double step(double r) {
if ( r < -0.5 )
return 0.0;
if ( r > 0.5 )
return 1.0;
return r+0.5;
}
bool HJetsAnalysis2::lessThanEnergy(Energy a, Energy b, double& weight) const {
if ( !fuzzy() ) {
if ( a < b ) {
weight = 1.0;
return true;
}
weight = 0.0;
return false;
}
double w = step((b-a)/theEnergyCutWidth);
if ( w == 0.0 ) {
weight = 0.0;
return false;
}
weight *= w;
return true;
}
bool HJetsAnalysis2::lessThanRapidity(double a, double b, double& weight) const {
if ( !fuzzy() ) {
if ( a < b ) {
weight = 1.0;
return true;
}
weight = 0.0;
return false;
}
double w = step((b-a)/theRapidityCutWidth);
if ( w == 0.0 ) {
weight = 0.0;
return false;
}
weight *= w;
return true;
}
bool HJetsAnalysis2::selectionCuts(const LorentzMomentum& h,
const vector<LorentzMomentum>& jets, double& weight){
double theCutWeight(1.0);
+
// higgs pt min and y_max
//std::cout << thehiggsPt_min/GeV << endl;
if ( !(lessThanRapidity(abs(h.rapidity()),thehiggsY_max,theCutWeight))){
theCutWeight = 0.0;
return false;
}
//std::cout << thehiggsPt_min/GeV << endl;
if ( !(lessThanEnergy(thehiggsPt_min,h.perp(),theCutWeight))){
theCutWeight = 0.0;
return false;
}
- if(jets.size() > 1){
- // opposite detectors
- double py = jets[0].rapidity()*jets[1].rapidity();
-
- if (theopposite_dir){
- if( !(lessThanRapidity(py,0.0,theCutWeight)) ){
- theCutWeight=0.0;
- return false;
- }
- }
- // rapidity gap
- if( !lessThanRapidity(theRapidityGap,deltaY(jets[0],jets[1]),theCutWeight)){
+ if(jets.size() > 1){
+ // opposite detectors
+ double py = jets[0].rapidity()*jets[1].rapidity();
+
+ if (theopposite_dir){
+ if( !(lessThanRapidity(py,0.0,theCutWeight)) ){
+ theCutWeight=0.0;
+ return false;
+ }
+ }
+ // rapidity gap
+ if( !lessThanRapidity(theRapidityGap,deltaY(jets[0],jets[1]),theCutWeight)){
theCutWeight = 0.0;
return false;
}
// invariant mass cut
-
+
Energy massj1j2;
massj1j2 = (jets[0]+jets[1]).m();
//std::cout << massj1j2/GeV << endl;
if( !lessThanEnergy(themjj_min,massj1j2,theCutWeight)){
theCutWeight = 0.0;
return false;
}
// higgs in the gap
if(thehiggs_ingap){
double minjrap = min(jets[0].rapidity(),jets[1].rapidity());
double maxjrap = max(jets[0].rapidity(),jets[1].rapidity());
if(!(lessThanRapidity(h.rapidity(),maxjrap,theCutWeight) && lessThanRapidity(minjrap,h.rapidity(),theCutWeight))){
theCutWeight = 0.0;
return false;
}
}
}
// should be a rescaling
weight=theCutWeight*weight;
return true;
}
void HJetsAnalysis2::analyze(const ParticleVector& out, long ieve, double weight) {
if ( weight == 0.0 )
return;
tcPDVector outType;
vector<LorentzMomentum> outMomenta;
for ( ParticleVector::const_iterator p = out.begin();
p != out.end(); ++p ) {
outType.push_back((**p).dataPtr());
outMomenta.push_back((**p).momentum());
}
theJetFinder->cluster(outType, outMomenta,
tcCutsPtr(), tcPDPtr(),tcPDPtr());
for ( vector<Ptr<JetRegion>::ptr>::iterator j =
theJetRegions.begin(); j != theJetRegions.end(); ++j )
(**j).reset();
LorentzMomentum higgs;
vector<LorentzMomentum> allJets;
bool gotHiggs = false;
tcPDVector::const_iterator pd = outType.begin();
vector<LorentzMomentum>::const_iterator p = outMomenta.begin();
for ( ; pd != outType.end(); ++pd, ++p ) {
if ( (**pd).id() == 25 ) {
gotHiggs = true;
higgs = *p;
} else if ( (**pd).coloured() ) {
allJets.push_back(*p);
}
}
assert(gotHiggs);
sort(allJets.begin(),allJets.end(),SortPt());
// @TODO
// perhaps set a PT min threshold and max abs rapidity ABSY_MAX
//
- // vector<LorentzMomentum> selectedJets20;
- // for ( size_t k = 0; k < allJets.size(); ++k ) {
- // if (allJets[k].perp()/GeV > 20.0 ){
- // selectedJets20.push_back(allJets[k]);
- // }
- // }
+ vector<LorentzMomentum> selectedJets20;
+ for ( size_t k = 0; k < allJets.size(); ++k ) {
+ if (allJets[k].perp()/GeV > 20.0 ){
+ selectedJets20.push_back(allJets[k]);
+ }
+ }
- // sort(selectedJets20.begin(),selectedJets20.end(),SortPt());
+ sort(selectedJets20.begin(),selectedJets20.end(),SortPt());
//
//std::cout << "number of jets from selected jets 20 " << selectedJets20.size() << endl;
// Should we only use this for plain NLO runs
size_t matchedJets = 0;
vector<LorentzMomentum> selectedJets;
for ( size_t k = 0; k < allJets.size(); ++k ) {
for ( vector<Ptr<JetRegion>::ptr>::const_iterator r = theJetRegions.begin();
- r != theJetRegions.end(); ++r ) {
+ r != theJetRegions.end(); ++r ) {
if ( (**r).matches(tCutsPtr(),k+1,allJets[k])) {
- ++matchedJets;
- selectedJets.push_back(allJets[k]);
- break;
+ ++matchedJets;
+ selectedJets.push_back(allJets[k]);
+ break;
}
}
+
}
+ // apply jets cuts according to the jet regions
+
+ double cutWeight = 1.0;
+ bool pass = true;
+ bool pass_jetcuts = true;
+
+ for ( vector<Ptr<JetRegion>::ptr>::const_iterator r = theJetRegions.begin();
+ r != theJetRegions.end(); ++r ) {
+ pass &= (**r).didMatch();
+ cutWeight *= (**r).cutWeight();
+ if ( !pass ) {
+ pass_jetcuts = false;
+ cutWeight = 0.0;
+ break;
+ }
+ }
+ weight=cutWeight*weight;
if ( matchedJets < theJetRegions.size() ) {
if ( !theFullEvent )
generator()->log() << "ATTENTION less than number of expected jets!\n";
return;
}
sort(selectedJets.begin(),selectedJets.end(),SortPt());
//
// for plain NLO, we do not want these cuts
// for NLO + PS, we do want these cuts
// For plain NLO, we do not use these cuts for the moment.
// @TODO: Looking for a better suggestion.
//
bool cutsOK = false;
- cutsOK = selectionCuts(higgs,selectedJets,weight);
-
+ // must pass jet cuts
+ if (pass_jetcuts){
+ cutsOK = selectionCuts(higgs,selectedJets,weight);
+ }
if(cutsOK){
- analyze(higgs,selectedJets,ieve,weight);
+ analyze(higgs,allJets,ieve,weight);
}
//
}
void HJetsAnalysis2::analyze(const LorentzMomentum& h,
const vector<LorentzMomentum>& jets,
long xid, double weight) {
const double pi = 3.14159265;
const double massWidth = 4000.0/64.0*theBinFactor;
const double rapWidth1 = 60.0/64.0*theBinFactor;
const double rapWidth = 10.0/32.0*theBinFactor;
const double etaWidth = 10.0/32.0*theBinFactor;
const double ptWidth = 400.0/40.0*theBinFactor;
const double ptHiggsWidth = 400.0/40.0*theBinFactor;
const double phiWidth = 2.*pi/32.0*theBinFactor;
const double ysWidth = 4.0/32.0*theBinFactor;
const double jetMassWidth = 100.0/20.0*theBinFactor;
const double HTWidth = 1000.0/32.0*theBinFactor;
const double higgsDijetDeltaPhiWidth = pi/32.0*theBinFactor;
const double higgsDijetpTWidth=100/32.0*theBinFactor;
const double threeJetMassWidth=2000/32.0*theBinFactor;
const double jetSumDeltaPhiWidth=2.*pi/32.0*theBinFactor;
// const double massWidth = 0.0;
// const double rapWidth = 0.0;
// const double etaWidth = 0.0;
// const double ptWidth = 0.0;
// const double ptHiggsWidth = 0.0;
// const double phiWidth = 0.0;
// const double ysWidth = 0.0;
// const double jetMassWidth = 0.0;
// const double HTWidth = 0.0;
// const double higgsDijetDeltaPhiWidth = 0.0;
// const double higgsDijetpTWidth=0.0;
// const double threeJetMassWidth=0.0;
// const double jetSumDeltaPhiWidth=0.0;
Ptr<StandardEventHandler>::tptr seh =
dynamic_ptr_cast<Ptr<StandardEventHandler>::tptr>(generator()->eventHandler());
Ptr<GeneralSampler>::tptr sampler =
dynamic_ptr_cast<Ptr<GeneralSampler>::tptr>(seh->sampler());
CrossSection maxXSection = sampler->maxXSec();
weight = weight*maxXSection/nanobarn;
//std::cout << "weight " << weight << endl;
//std::cout << "maxXSection " << maxXSection/nanobarn << endl;
size_t njets = jets.size();
size_t id = xid;
//std::cout << "number of jets from jets.size()" << njets << endl;
totalxs.count(EventContribution(0.0,weight,0.0),id);
size_t njets20 = 0;
for ( size_t k = 0; k< njets; ++k ) {
if (jets[k].perp()/GeV > 20.0 ) ++njets20;
}
//std::cout << "number of jets" << njets20 << endl;
njetsex_xs.count(EventContribution(njets20,weight,0.0),id);
for (size_t i = 0; i < njets20; ++i) {
if (njets20 >= i) {
njetsinc_xs.count(EventContribution(i,weight,0.0),id);
}
}
higgsPt.count(EventContribution(h.perp()/GeV,weight,ptHiggsWidth),id);
higgsY.count(EventContribution(h.rapidity(),weight,rapWidth),id);
higgsEta.count(EventContribution(h.eta(),weight,etaWidth),id);
// For containing HT
double HTval=0;
size_t n_max=5; // max number of jets to histogram
size_t njet_max = min(njets,n_max);
for ( size_t k = 0; k < njets; ++k ) {
// Loop over jet momenta to calculate scalar HT
HTval+=jets[k].perp()/GeV;
}
for ( size_t k = 0; k < njet_max; ++k ) {
jetPt[k].count(EventContribution(jets[k].perp()/GeV,weight,ptWidth),id);
jetY[k].count(EventContribution(jets[k].rapidity(),weight,rapWidth),id);
jetHiggsDeltaPhi[k].count(EventContribution(deltaPhi(jets[k],h),weight,phiWidth),id);
jetHiggsDeltaY[k].count(EventContribution(deltaY(jets[k],h),weight,rapWidth),id);
jetHiggsDeltaR[k].count(EventContribution(deltaR(jets[k],h),weight,rapWidth),id);
jetHiggsMass[k].count(EventContribution(mass(jets[k],h),weight,massWidth),id);
jetMass[k].count(EventContribution(jets[k].m()/GeV,weight,jetMassWidth),id);
for ( size_t l = k + 1; l < njet_max; ++l ) {
jetDeltaPhi[k][l-k-1].count(EventContribution(deltaPhi(jets[k],jets[l]),weight,phiWidth),id);
jetDeltaY[k][l-k-1].count(EventContribution(deltaY(jets[k],jets[l]),weight,rapWidth),id);
jetYdotY[k][l-k-1].count(EventContribution(jets[k].rapidity()*jets[l].rapidity(),weight,rapWidth1),id);
jetDeltaR[k][l-k-1].count(EventContribution(deltaR(jets[k],jets[l]),weight,rapWidth),id);
jetPairMass[k][l-k-1].count(EventContribution(mass(jets[k],jets[l]),weight,massWidth),id);
}
}
HT.count(EventContribution(HTval,weight,HTWidth),id);
if ( njets > 2 ) {
double ystar = jets[2].rapidity() - 0.5 * (jets[0].rapidity()+jets[1].rapidity());
ystar /= jets[0].rapidity()-jets[1].rapidity();
yStarNormed.count(EventContribution(ystar,weight,ysWidth),id);
}
if( njets > 1 ){
LorentzMomentum dijetSum=jets[0]+jets[1];
double phiDiff=deltaPhi( dijetSum, h);
higgsDijetDeltaPhi.count(EventContribution( abs( phiDiff ), weight, higgsDijetDeltaPhiWidth),id);
double higgsDijetpTSum=(dijetSum+h).perp()/GeV;
higgsDijetpT.count(EventContribution( higgsDijetpTSum, weight, higgsDijetpTWidth),id );
}
if ( njets > 2 ) {
threeJetMass.count( EventContribution(mass(jets[0],jets[1],jets[2]),weight,threeJetMassWidth),id);
}
// jeppe DeltaPhi, 1001.3822 eq. 9.
if( njets > 1 ){ // need at lest one jet with higher y and one with lower y than the Higgs
// to be changed (probably)
double smallest_rap=jets[0].rapidity();
double largest_rap=jets[0].rapidity();
// Find largest and smallest delta rap w.r.t. Higgs
for ( size_t k = 0; k < njets; ++k ) {
if( smallest_rap > jets[k].rapidity() ) smallest_rap=jets[k].rapidity();
if( largest_rap < jets[k].rapidity() ) largest_rap=jets[k].rapidity();
}
// If one jet has higher rapidity and one lower than Higgs
if( smallest_rap< h.rapidity() and largest_rap> h.rapidity() ){
// Collect jets into those having lower rapidity
// and those having higher rapidity than the Higgs
LorentzMomentum p_a;
LorentzMomentum p_b;
for ( size_t k = 0; k < njets; ++k ) {
if( jets[k].rapidity() < h.rapidity() ) p_a+=jets[k];
else if( jets[k].rapidity() > h.rapidity() ) p_b+=jets[k];
}
double phiDiff=deltaPhi( p_a, p_b );
jetSumDeltaPhi.count( EventContribution(phiDiff,weight,jetSumDeltaPhiWidth),id );
}
}
}
void HJetsAnalysis2::finalizeHistogram(Histogram& h, XML::Element& elem) {
h.finalize();
elem.append(h.toXML());
}
void HJetsAnalysis2::dofinish() {
AnalysisHandler::dofinish();
Ptr<StandardEventHandler>::tptr seh =
dynamic_ptr_cast<Ptr<StandardEventHandler>::tptr>(generator()->eventHandler());
Ptr<GeneralSampler>::tptr sampler =
dynamic_ptr_cast<Ptr<GeneralSampler>::tptr>(seh->sampler());
unsigned long attemptedPoints = sampler->attempts();
double sumOfWeights = sampler->sumWeights();
double sumOfSquaredWeights = sampler->sumWeights2();
CrossSection maxXSection = sampler->maxXSec();
XML::Element elem(XML::ElementTypes::Element,"Run");
elem.appendAttribute("name",generator()->runName());
elem.appendAttribute("attemptedPoints",attemptedPoints);
elem.appendAttribute("sumOfWeights",sumOfWeights*maxXSection/nanobarn);
elem.appendAttribute("sumOfSquaredWeights",sumOfSquaredWeights*sqr(maxXSection/nanobarn));
//elem.appendAttribute("sumOfWeights",sumOfWeights);
//elem.appendAttribute("sumOfSquaredWeights",sumOfSquaredWeights);
XML::Element xhistos(XML::ElementTypes::Element,"Histograms");
finalizeHistogram(totalxs,xhistos);
finalizeHistogram(njetsex_xs,xhistos);
finalizeHistogram(njetsinc_xs,xhistos);
finalizeHistogram(higgsPt,xhistos);
finalizeHistogram(higgsY,xhistos);
finalizeHistogram(higgsEta,xhistos);
finalizeHistogram(HT,xhistos);
finalizeHistogram(higgsDijetpT,xhistos);
finalizeHistogram(threeJetMass,xhistos);
finalizeHistogram(jetSumDeltaPhi,xhistos);
size_t n_min = 5;
for ( size_t k = 0; k < n_min; ++k ) {
finalizeHistogram(jetPt[k],xhistos);
finalizeHistogram(jetY[k],xhistos);
finalizeHistogram(jetHiggsDeltaPhi[k],xhistos);
finalizeHistogram(jetHiggsDeltaY[k],xhistos);
finalizeHistogram(jetHiggsDeltaR[k],xhistos);
finalizeHistogram(jetHiggsMass[k],xhistos);
finalizeHistogram(jetMass[k],xhistos);
for ( size_t l = k + 1; l < n_min; ++l ) {
finalizeHistogram(jetDeltaPhi[k][l-k-1],xhistos);
finalizeHistogram(jetDeltaY[k][l-k-1],xhistos);
finalizeHistogram(jetYdotY[k][l-k-1],xhistos);
finalizeHistogram(jetDeltaR[k][l-k-1],xhistos);
finalizeHistogram(jetPairMass[k][l-k-1],xhistos);
}
}
if ( theJetRegions.size() > 1 ) {
finalizeHistogram(higgsDijetDeltaPhi,xhistos);
}
if ( theJetRegions.size() > 2 ) {
finalizeHistogram(yStarNormed,xhistos);
}
elem.append(xhistos);
string fname = name() + ".xml";
ofstream runXML(fname.c_str());
XML::ElementIO::put(elem,runXML);
}
void HJetsAnalysis2::doinitrun() {
AnalysisHandler::doinitrun();
size_t njets = 5; //theJetRegions.size();
const double pi = 3.14159265;
totalxs = Histogram("sigma_tot",Histogram::regularBinEdges(-0.5,0.5,1));
njetsex_xs = Histogram("sigma_njetex",Histogram::regularBinEdges(-0.5,8.5,9));
njetsinc_xs = Histogram("sigma_njetinc",Histogram::regularBinEdges(-0.5,8.5,9));
higgsPt = Histogram("higgsPt",Histogram::regularBinEdges(0,400,40),true);
higgsY = Histogram("higgsY",Histogram::regularBinEdges(-5,5,32));
higgsEta = Histogram("higgsEta",Histogram::regularBinEdges(-5,5,32));
HT = Histogram("HT",Histogram::regularBinEdges(0,1000,32),true);
higgsDijetDeltaPhi = Histogram("higgsDijetDeltaPhi",Histogram::regularBinEdges(0,pi,32),true,true);
higgsDijetpT = Histogram("higgsDijetpT",Histogram::regularBinEdges(0,100,20),true);
threeJetMass = Histogram("threeJetMass",Histogram::regularBinEdges(0,2000,32),true,false);
jetSumDeltaPhi = Histogram("jetSumDeltaPhi",Histogram::regularBinEdges(-pi,pi,32),true,true);
for ( size_t k = 0; k < njets; ++k ) {
ostringstream name;
name << "jet" << (k+1);
jetPt.push_back(Histogram(name.str() + "Pt",Histogram::regularBinEdges(0,400,40),true));
jetY.push_back(Histogram(name.str() + "Y",Histogram::regularBinEdges(-5,5,32)));
jetHiggsDeltaPhi.push_back(Histogram(name.str() + "HiggsDeltaPhi",Histogram::regularBinEdges(-pi,pi,32),true,true));
jetHiggsDeltaY.push_back(Histogram(name.str() + "HiggsDeltaY",Histogram::regularBinEdges(0,10,32),true));
jetHiggsDeltaR.push_back(Histogram(name.str() + "HiggsDeltaR",Histogram::regularBinEdges(0,10,32),true));
jetHiggsMass.push_back(Histogram(name.str() + "HiggsMass",Histogram::regularBinEdges(0,4000,64),true));
jetMass.push_back(Histogram(name.str() + "JetMass",Histogram::regularBinEdges(0,100,32),true, false));
vector<Histogram> xjetDeltaPhi;
vector<Histogram> xjetDeltaY;
vector<Histogram> xjetYdotY;
vector<Histogram> xjetDeltaR;
vector<Histogram> xjetPairMass;
for ( size_t l = k + 1; l < njets; ++l ) {
ostringstream xname;
xname << "jet" << (k+1) << (l+1);
xjetDeltaPhi.push_back(Histogram(xname.str() + "DeltaPhi",Histogram::regularBinEdges(-pi,pi,32),true,true));
xjetDeltaY.push_back(Histogram(xname.str() + "DeltaY",Histogram::regularBinEdges(0,10,32),true));
xjetYdotY.push_back(Histogram(xname.str() + "YdotY",Histogram::regularBinEdges(-30,30,64),true));
xjetDeltaR.push_back(Histogram(xname.str() + "DeltaR",Histogram::regularBinEdges(0,10,32),true));
xjetPairMass.push_back(Histogram(xname.str() + "Mass",Histogram::regularBinEdges(0,4000,64),true));
}
jetDeltaPhi.push_back(xjetDeltaPhi);
jetDeltaY.push_back(xjetDeltaY);
jetYdotY.push_back(xjetYdotY);
jetDeltaR.push_back(xjetDeltaR);
jetPairMass.push_back(xjetPairMass);
}
yStarNormed = Histogram("yStarNormed",Histogram::regularBinEdges(-2,2,32));
}
IBPtr HJetsAnalysis2::clone() const {
return new_ptr(*this);
}
IBPtr HJetsAnalysis2::fullclone() const {
return new_ptr(*this);
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void HJetsAnalysis2::persistentOutput(PersistentOStream & os) const {
os << theJetFinder << theJetRegions << theFullEvent
<< theBinFactor << theRapidityGap << ounit(themjj_min,GeV)
- << ounit(thehiggsPt_min,GeV) << thehiggsY_max << theopposite_dir
+ << ounit(thehiggsPt_min,GeV) << thehiggsY_max << theopposite_dir << thenjets_min
<< thehiggs_ingap << theFuzzy << ounit(theEnergyCutWidth,GeV) << theRapidityCutWidth;
}
void HJetsAnalysis2::persistentInput(PersistentIStream & is, int) {
is >> theJetFinder >> theJetRegions >> theFullEvent
>> theBinFactor >> theRapidityGap >> iunit(themjj_min,GeV)
- >> iunit(thehiggsPt_min,GeV) >> thehiggsY_max >> theopposite_dir
+ >> iunit(thehiggsPt_min,GeV) >> thehiggsY_max >> theopposite_dir >> thenjets_min
>> thehiggs_ingap >> theFuzzy >> iunit(theEnergyCutWidth,GeV) >> theRapidityCutWidth;
}
// *** Attention *** The following static variable is needed for the type
// description system in ThePEG. Please check that the template arguments
// are correct (the class and its base class), and that the constructor
// arguments are correct (the class name and the name of the dynamically
// loadable library where the class implementation can be found).
DescribeClass<HJetsAnalysis2,AnalysisHandler>
describeHJetsHJetsAnalysis2("HJets::HJetsAnalysis2", "HJetsAnalysis2.so");
void HJetsAnalysis2::Init() {
static ClassDocumentation<HJetsAnalysis2> documentation
("There is no documentation for the HJetsAnalysis2 class");
static Reference<HJetsAnalysis2,JetFinder> interfaceJetFinder
("JetFinder",
"",
&HJetsAnalysis2::theJetFinder, false, false, true, false, false);
static RefVector<HJetsAnalysis2,JetRegion> interfaceJetRegions
("JetRegions",
"",
&HJetsAnalysis2::theJetRegions, -1, false, false, true, false, false);
static Switch<HJetsAnalysis2,bool> interfaceFullEvent
("FullEvent",
"Analyse a fully simulated event",
&HJetsAnalysis2::theFullEvent, false, false, false);
static SwitchOption interfaceFullEventYes
(interfaceFullEvent,
"Yes",
"",
true);
static SwitchOption interfaceFullEventNo
(interfaceFullEvent,
"No",
"",
false);
static Parameter<HJetsAnalysis2,double> interfaceBinFactor
("BinFactor",
"The factor multiplying the bin width",
&HJetsAnalysis2::theBinFactor, 1.0, 0.0, 1.0, true, false, true );
static Parameter<HJetsAnalysis2,double> interfaceRapidityGap
("RapidityGap",
"The distance in pseudorapidity between the two tagging jets",
&HJetsAnalysis2::theRapidityGap, 0.0, 0.0, 20.0, true, false, true);
static Parameter<HJetsAnalysis2,Energy> interfacemjj_min
("mjj_min",
"The minimum dijet invariant mass of the two taggin jets.",
&HJetsAnalysis2::themjj_min, GeV, 1.0*GeV, 0.0*GeV, 0*GeV,
false, false, Interface::lowerlim);
- static Parameter<HJetsAnalysis2,Energy> interfacehiggsPt_min
+ static Parameter<HJetsAnalysis2,Energy> interfacehiggsPt_min
("higgsPt_min",
"Minimun transverse momentum of higgs.",
&HJetsAnalysis2::thehiggsPt_min, GeV, 1.0*GeV, 0.0*GeV, 0*GeV,
false, false, Interface::lowerlim);
static Parameter<HJetsAnalysis2,double> interfacehiggsY_max
("higgsY_max",
"maximum absolute rapidity of higgs",
&HJetsAnalysis2::thehiggsY_max, 0, 0, 14000, true, false, true);
static Switch<HJetsAnalysis2,bool> interfaceopposite_dir
("opposite_dir",
"Should the tagging jets reside in opposite directions (Yes/No)?",
&HJetsAnalysis2::theopposite_dir, false, false, false);
static SwitchOption interfaceopposite_dirYes
(interfaceopposite_dir,
"Yes",
"",
true);
static SwitchOption interfaceopposite_dirNo
(interfaceopposite_dir,
"No",
"",
false);
+ static Parameter<HJetsAnalysis2,unsigned int> interfacenjets_min
+ ("njets_min",
+ "The minimum number of jets",
+ &HJetsAnalysis2::thenjets_min, 2, 1, 10, true, false, true);
+
+
static Switch<HJetsAnalysis2,bool> interfacehiggs_ingap
("higgs_ingap",
"Should the higgs reside in the gap (Yes/No)?",
&HJetsAnalysis2::thehiggs_ingap, false, false, false);
static SwitchOption interfacehiggs_ingapYes
(interfacehiggs_ingap,
"Yes",
"",
true);
static SwitchOption interfacehiggs_ingapNo
(interfacehiggs_ingap,
"No",
"",
false);
static Switch<HJetsAnalysis2,bool> interfaceFuzzy
("Fuzzy",
"Make this jet region a fuzzy cut",
&HJetsAnalysis2::theFuzzy, false, false, false);
static SwitchOption interfaceFuzzyOn
(interfaceFuzzy,
"Yes",
"",
true);
static SwitchOption interfaceFuzzyOff
(interfaceFuzzy,
"No",
"",
false);
static Parameter<HJetsAnalysis2,Energy> interfaceEnergyCutWidth
("EnergyCutWidth",
"The pt cut smearing.",
&HJetsAnalysis2::theEnergyCutWidth, GeV, 1.0*GeV, 0.0*GeV, 0*GeV,
false, false, Interface::lowerlim);
static Parameter<HJetsAnalysis2,double> interfaceRapidityCutWidth
("RapidityCutWidth",
"The rapidity cut smearing.",
&HJetsAnalysis2::theRapidityCutWidth, 0.1, 0.0, 0,
false, false, Interface::lowerlim);
}
diff --git a/Analysis/HJetsAnalysis2.h b/Analysis/HJetsAnalysis2.h
--- a/Analysis/HJetsAnalysis2.h
+++ b/Analysis/HJetsAnalysis2.h
@@ -1,390 +1,394 @@
// -*- C++ -*-
#ifndef HJets_HJetsAnalysis2_H
#define HJets_HJetsAnalysis2_H
//
// This is the declaration of the HJetsAnalysis2 class.
//
#include "ThePEG/Handlers/AnalysisHandler.h"
#include "ThePEG/Cuts/JetFinder.h"
#include "ThePEG/Cuts/JetRegion.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/EventRecord/SubProcess.h"
#include "ThePEG/EventRecord/SubProcessGroup.h"
#include "Herwig++/Utilities/Statistics/Histogram.h"
namespace HJets {
using namespace ThePEG;
/**
* Here is the documentation of the HJetsAnalysis2 class.
*
* @see \ref HJetsAnalysis2Interfaces "The interfaces"
* defined for HJetsAnalysis2.
*/
class HJetsAnalysis2: public AnalysisHandler {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
HJetsAnalysis2();
/**
* The destructor.
*/
virtual ~HJetsAnalysis2();
//@}
public:
/** @name Virtual functions required by the AnalysisHandler class. */
//@{
/**
* Analyze a given Event. Note that a fully generated event
* may be presented several times, if it has been manipulated in
* between. The default version of this function will call transform
* to make a lorentz transformation of the whole event, then extract
* all final state particles and call analyze(tPVector) of this
* analysis object and those of all associated analysis objects. The
* default version will not, however, do anything on events which
* have not been fully generated, or have been manipulated in any
* way.
* @param event pointer to the Event to be analyzed.
* @param ieve the event number.
* @param loop the number of times this event has been presented.
* If negative the event is now fully generated.
* @param state a number different from zero if the event has been
* manipulated in some way since it was last presented.
*/
virtual void analyze(tEventPtr event, long ieve, int loop, int state);
//@}
/**
* Analyze one subprocess, given the event number it belongs to
*/
void analyze(const ParticleVector&, long, double);
/**
* Analyze higgs and jet momenta, given the event number it belongs to
*/
void analyze(const LorentzMomentum& h,
const vector<LorentzMomentum>& jets,
long, double);
protected:
/**
* Initialize this object. Called in the run phase just before
* a run begins.
*/
virtual void doinitrun();
/**
* Finalize this object. Called in the run phase just after a
* run has ended. Used eg. to write out statistics.
*/
virtual void dofinish();
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
protected:
/**
* Return true, if this jet region is fuzzy
*/
bool fuzzy() const { return theFuzzy; }
private:
/**
* Check for cuts
*/
bool selectionCuts(const LorentzMomentum& h, const vector<LorentzMomentum>& jets, double& weight);
/**
* Perform a (potentially) fuzzy check on energy-type quantities
*/
bool lessThanEnergy(Energy a, Energy b, double& weight) const;
/**
* Perform a (potentially) fuzzy check on angular-type quantities
*/
bool lessThanRapidity(double a, double b, double& weight) const;
/**
* The jet finder to use
*/
Ptr<JetFinder>::ptr theJetFinder;
/**
* The jet regions to match.
*/
vector<Ptr<JetRegion>::ptr> theJetRegions;
/**
* True if full evnet is taken
*/
bool theFullEvent;
/**
* Factor multiplying bin width (bin smearing)
*/
double theBinFactor;
/**
* Rapidity gap
*/
double theRapidityGap;
/**
* di-jet mass cut
*/
Energy themjj_min;
/**
* Higgs Pt min
*/
Energy thehiggsPt_min;
/**
* Higgs max rapidity
*/
double thehiggsY_max;
/**
* True if leading jets are opposite in rapidity
*/
bool theopposite_dir;
+ /**
+ * Mininum number of analyzed jets
+ */
+ unsigned int thenjets_min;
/**
* True if higgs is rapidity gap
*/
bool thehiggs_ingap;
/**
* True if Fuzzy cuts
*/
bool theFuzzy;
/**
* The smearing width for the pt or mass cuts, if fuzzy
*/
Energy theEnergyCutWidth;
/**
* The smearing width for the rapidity cut, if fuzzy
*/
double theRapidityCutWidth;
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
HJetsAnalysis2 & operator=(const HJetsAnalysis2 &);
/**
* Histogram indices for total xs
*/
Statistics::Histogram totalxs;
/**
* Histogram indices for exclusive xs
*/
Statistics::Histogram njetsex_xs;
/**
* Histogram indices for inclusive njet xs
*/
Statistics::Histogram njetsinc_xs;
/**
* Histogram indices for higgs pts
*/
Statistics::Histogram higgsPt;
/**
* Statistics::Histogram indices for higgs rapidities
*/
Statistics::Histogram higgsY;
/**
* Statistics::Histogram indices for higgs pseudo rapidities
*/
Statistics::Histogram higgsEta;
/**
* Statistics::Scalar sum of transverse momenta
*/
Statistics::Histogram HT;
/**
* Statistics::Histogram of the transverse momentum
* of higgs plus the 2 hardest jets
*/
Statistics::Histogram higgsDijetpT;
/**
* Statistics::Histogram of delta phi between higgs
* and the two hardest jets
*/
Statistics::Histogram higgsDijetDeltaPhi;
/**
* Statistics::Histogram of the invariant mass of
* the three hardest jets
*/
Statistics::Histogram threeJetMass;
/**
* Statistics::Histogram of the azimuthal angle as defined
* in 1001.3822 eq. 9.
*/
Statistics::Histogram jetSumDeltaPhi;
/**
* Statistics::Histogram indices for jet pts
*/
vector<Statistics::Histogram> jetPt;
/**
* Statistics::Histogram indices for jet rapidities
*/
vector<Statistics::Histogram> jetY;
/**
* Statistics::Histogram indices for jet delta phis
*/
vector<vector<Statistics::Histogram> > jetDeltaPhi;
/**
* Statistics::Histogram indices for jet y dot y
*/
vector<vector<Statistics::Histogram> > jetYdotY;
/**
* Statistics::Histogram indices for jet delta y
*/
vector<vector<Statistics::Histogram> > jetDeltaY;
/**
* Statistics::Histogram indices for jet delta R
*/
vector<vector<Statistics::Histogram> > jetDeltaR;
/**
* Statistics::Histogram indices for jet pair masses
*/
vector<vector<Statistics::Histogram> > jetPairMass;
/**
* Statistics::Histogram indices for jet higgs delta phis
*/
vector<Statistics::Histogram> jetHiggsDeltaPhi;
/**
* Statistics::Histogram indices for jet higgs delta y
*/
vector<Statistics::Histogram> jetHiggsDeltaY;
/**
* Statistics::Histogram indices for jet higgs delta R
*/
vector<Statistics::Histogram> jetHiggsDeltaR;
/**
* Statistics::Histogram indices for jet higgs masses
*/
vector<Statistics::Histogram> jetHiggsMass;
/**
* Statistics::Histogram for invariant mass of jets
*/
vector<Statistics::Histogram> jetMass;
/**
* y*
*/
Statistics::Histogram yStarNormed;
/**
* Finalize histogram and add to XML element
*/
void finalizeHistogram(Statistics::Histogram&,XML::Element&);
};
}
#endif /* HJets_HJetsAnalysis2_H */

File Metadata

Mime Type
text/x-diff
Expires
Sun, Feb 23, 2:30 PM (15 h, 38 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4483767
Default Alt Text
(35 KB)

Event Timeline