// -*- C++ -*- //Adapted from standard MC_JetAnalysis.cc by R. King #include "MC_JetMatrix.hh" //#include "/afs/cern.ch/sw/lcg/external/MCGenerators/rivet/1.2.1/share/include/Rivet/Tools/Logging.hh" #include "Rivet/Tools/Logging.hh" //#include "/afs/cern.ch/sw/lcg/external/MCGenerators/rivet/1.2.1/share/include/Rivet/Projections/FastJets.hh #include "Rivet/Projections/FastJets.hh" //#include "/afs/cern.ch/sw/lcg/external/MCGenerators/rivet/1.2.1/share/include/Rivet/RivetAIDA.hh" #include "Rivet/RivetAIDA.hh" namespace Rivet { MC_JetMatrix::MC_JetMatrix(const string& name, const size_t& njet, const string& jetpro_name) : Analysis(name), m_njet(njet), m_jetpro_name(jetpro_name), _h_log10_d(njet, NULL), _h_log10_R(njet+1, NULL), _h_pT_jet(njet, NULL), _h_eta_jet(njet, NULL) { setNeedsCrossSection(true); } //Projections: number of jets, jet production? // Book histograms void MC_JetMatrix::init() { //make a filestream JetFile.open("jetcount.txt"); _h_jet_multi_exclusive = bookHistogram1D("jet_multi_exclusive", m_njet+3, -0.5, m_njet+3-0.5); _h_jet_multi_inclusive = bookHistogram1D("jet_multi_inclusive", m_njet+3, -0.5, m_njet+3-0.5); _h_jet_multi_ratio = bookDataPointSet("jet_multi_ratio", m_njet+2, 0.5, m_njet+3-0.5); //a couple of VBF histograms _h_jet1_etaVBF = bookHistogram1D("jet1_etaVBF", 50,-5.0,5.0); //eta _h_jet2_etaVBF = bookHistogram1D("jet2_etaVBF", 50,-5.0,5.0); //eta _h_jet1_ptVBF = bookHistogram1D("jet1_ptVBF", 100, 0, 400); _h_jet2_ptVBF = bookHistogram1D("jet2_ptVBF", 100, 0, 400); _h_jet_mass = bookHistogram1D("jet_mass",100, 0, 1000); _h_jet_massVBF = bookHistogram1D("jet_massVBF",100, 0, 1000); //DR between leading jets. _h_jet_DetaVBF = bookHistogram1D("jet_DetaVBF",20, 0, 10); _h_jet_DRVBF = bookHistogram1D("jet_DRVBF",25,0,8); //now the exclusive hists: _h_jet1_etaVBFex = bookHistogram1D("jet1_etaVBFex", 50,-5.0,5.0); //eta _h_jet2_etaVBFex = bookHistogram1D("jet2_etaVBFex", 50,-5.0,5.0); //eta _h_jet1_ptVBFex = bookHistogram1D("jet1_ptVBFex", 100, 0, 400); _h_jet2_ptVBFex = bookHistogram1D("jet2_ptVBFex", 100, 0, 400); _h_jet_massex = bookHistogram1D("jet_massex",100, 0, 1000); _h_jet_massVBFex = bookHistogram1D("jet_massVBFex",100, 0, 1000); //DR between leading jets. _h_jet_DetaVBFex = bookHistogram1D("jet_DetaVBFex",20, 0, 10); _h_jet_DRVBFex = bookHistogram1D("jet_DRVBFex",25,0,8); } // Do the analysis void MC_JetMatrix::analyze(const Event & e) { double weight = e.weight(); const FastJets& jetpro = applyProjection(e, m_jetpro_name); // jet resolutions and integrated jet rates const fastjet::ClusterSequence* seq = jetpro.clusterSeq(); const Jets& jets = jetpro.jetsByPt(20.0); //does this mean jets>20GeV? _h_jet_multi_exclusive->fill(jets.size(), weight); for (size_t i=0; i=i) { _h_jet_multi_inclusive->fill(i, weight); } } //Write the event number and jet number to the outfile JetFile <<"EventNumber"<< Run::_numEvents<< "Number of Jets "<20GeV //invariant mass>100GeV if(jets.size()>=2){ const double jet1pT = jets[0].momentum().pT(); const double jet2pT = jets[1].momentum().pT(); FourMomentum jet1 = jets[0].momentum(); FourMomentum jet2 = jets[1].momentum(); double m = FourMomentum(jet1+jet2).mass(); _h_jet_mass->fill(m, weight); //DR = sqrt(etadiff + phidiff) double DR = 0; double Deta = jets[0].momentum().eta() - jets[1].momentum().eta(); double Dphi = jets[0].momentum().phi() - jets[1].momentum().phi(); if(Dphi>3.14159) Dphi = Dphi-3.14159; DR = sqrt( Deta*Deta + Dphi*Dphi); if(m>=40){ //use lower cut off for validation purposes //now fill Histograms: _h_jet_massVBF->fill(m, weight); _h_jet1_etaVBF->fill(jets[0].momentum().eta(), weight); _h_jet2_etaVBF->fill(jets[1].momentum().eta(), weight); _h_jet1_ptVBF->fill(jets[0].momentum().pT(), weight); _h_jet2_ptVBF->fill(jets[1].momentum().pT(), weight); _h_jet_DetaVBF->fill(Deta,weight); _h_jet_DRVBF->fill(DR,weight); if(jets.size()==2)//exclusive section { _h_jet_massVBFex->fill(m, weight); _h_jet1_etaVBFex->fill(jets[0].momentum().eta(), weight); _h_jet2_etaVBFex->fill(jets[1].momentum().eta(), weight); _h_jet1_ptVBFex->fill(jets[0].momentum().pT(), weight); _h_jet2_ptVBFex->fill(jets[1].momentum().pT(), weight); _h_jet_DetaVBFex->fill(Deta,weight); _h_jet_DRVBFex->fill(DR,weight); } }//end of the mass cut } }//end of analysis // Finalize void MC_JetMatrix::finalize() { //close the file JetFile.close(); // fill inclusive jet multi ratio int Nbins=_h_jet_multi_inclusive->axis().bins(); std::vector ratio(Nbins-1, 0.0); std::vector err(Nbins-1, 0.0); for (int i=0; ibinHeight(i)>0.0 && _h_jet_multi_inclusive->binHeight(i+1)>0.0) { ratio[i]=_h_jet_multi_inclusive->binHeight(i+1)/_h_jet_multi_inclusive->binHeight(i); double relerr_i=_h_jet_multi_inclusive->binError(i)/_h_jet_multi_inclusive->binHeight(i); double relerr_j=_h_jet_multi_inclusive->binError(i+1)/_h_jet_multi_inclusive->binHeight(i+1); err[i]=ratio[i]*(relerr_i+relerr_j); } } _h_jet_multi_ratio->setCoordinate(1, ratio, err); scale(_h_jet_multi_exclusive, crossSection()/sumOfWeights()); scale(_h_jet_multi_inclusive, crossSection()/sumOfWeights()); scale(_h_jet1_etaVBF, crossSection()/sumOfWeights()); scale(_h_jet2_etaVBF, crossSection()/sumOfWeights()); scale(_h_jet1_ptVBF, crossSection()/sumOfWeights()); scale(_h_jet2_ptVBF, crossSection()/sumOfWeights()); scale(_h_jet_mass, crossSection()/sumOfWeights()); scale(_h_jet_massVBF, crossSection()/sumOfWeights()); scale(_h_jet_DetaVBF, crossSection()/sumOfWeights()); scale(_h_jet_DRVBF, crossSection()/sumOfWeights()); scale(_h_jet1_etaVBFex, crossSection()/sumOfWeights()); scale(_h_jet2_etaVBFex, crossSection()/sumOfWeights()); scale(_h_jet1_ptVBFex, crossSection()/sumOfWeights()); scale(_h_jet2_ptVBFex, crossSection()/sumOfWeights()); scale(_h_jet_massex, crossSection()/sumOfWeights()); scale(_h_jet_massVBFex, crossSection()/sumOfWeights()); scale(_h_jet_DetaVBFex, crossSection()/sumOfWeights()); scale(_h_jet_DRVBFex, crossSection()/sumOfWeights()); } }