diff --git a/PU14/EventMixer.cc b/PU14/EventMixer.cc index e76299b..bf1f763 100644 --- a/PU14/EventMixer.cc +++ b/PU14/EventMixer.cc @@ -49,14 +49,15 @@ bool EventMixer::next_event() { _pu_event_weight = 1; // first get the hard event - if (! _hard->append_next_event(_particles,_hard_event_weight,0)) return false; + if (! _hard->append_next_event(_particles,_hard_event_weight,_posX,_posY,0)) return false; unsigned hard_size = _particles.size(); // add pileup if available if (_pileup.get()){ for (int i = 1; i <= _npu; i++) { - if (! _pileup->append_next_event(_particles,_pu_event_weight,i)) return false; + double dumX, dumY; + if (! _pileup->append_next_event(_particles,_pu_event_weight,dumX,dumY,i)) return false; } } diff --git a/PU14/EventMixer.hh b/PU14/EventMixer.hh index 50b3710..f15b0a1 100644 --- a/PU14/EventMixer.hh +++ b/PU14/EventMixer.hh @@ -42,6 +42,8 @@ public: double weight() {return _hard_event_weight * _pu_event_weight;} double pu_weight() {return _pu_event_weight;} double hard_weight() {return _hard_event_weight;} + double productionX() {return _posX;} + double productionY() {return _posY;} /// returns the number of pileup events generated in the last mixed event int npu() const {return _npu;} @@ -73,6 +75,7 @@ private: std::vector _particles; double _hard_event_weight, _pu_event_weight; + double _posX, _posY; }; #endif // __EVENTMIXER_HH__ diff --git a/PU14/EventSource.cc b/PU14/EventSource.cc index 49834fa..2373d38 100644 --- a/PU14/EventSource.cc +++ b/PU14/EventSource.cc @@ -30,7 +30,8 @@ void EventSource::open_stream(const std::string & filename) { //---------------------------------------------------------------------- bool EventSource::append_next_event(std::vector & particles, - double &event_weight, int vertex_number) { + double &event_weight, double &prodX, double &prodY, + int vertex_number) { PseudoJet particle; string line; double px, py, pz, m, E; @@ -48,6 +49,7 @@ bool EventSource::append_next_event(std::vector & particles, // event (try to make the check as efficient as possible). if (line[0] == 'e' && line.substr(0,3) == "end") break; + bool bStop = false; // if the line says "weight", we multiply current weight by the number that follows if(line[0] == 'w' && line.substr(0, 6) == "weight") { @@ -57,9 +59,25 @@ bool EventSource::append_next_event(std::vector & particles, readline >> dummy >> temp; if(temp > 0) event_weight = event_weight * temp; - continue; + bStop = true; } + //extraction production point info for hybrid event generator + if(line.find("cross") < line.size()) { + std::size_t posX = line.find("X"); + std::string strX = line.substr(posX+2,6); + + std::size_t posY = line.find("Y"); + std::string strY = line.substr(posY+2,6); + + prodX = atof(strX.c_str()); + prodY = atof(strY.c_str()); + + bStop = true; + } + + if(bStop) continue; + // FastIStringStream is not a proper stream, but it's a lot faster // than standard istringstream. FastIStringStream readline(line.c_str()); diff --git a/PU14/EventSource.hh b/PU14/EventSource.hh index 150b4c7..0aa4277 100644 --- a/PU14/EventSource.hh +++ b/PU14/EventSource.hh @@ -25,7 +25,8 @@ public: /// appends the particles from the next event that is read onto the /// particles vector. bool append_next_event(std::vector & particles, - double &event_weight, int vertex_number = 0); + double &event_weight, double &prodX, double &prodY, + int vertex_number = 0); private: std::istream * _stream; diff --git a/runHybridJetAnalysis.cc b/runHybridJetAnalysis.cc new file mode 100644 index 0000000..f7ca685 --- /dev/null +++ b/runHybridJetAnalysis.cc @@ -0,0 +1,160 @@ +#include +#include + +#include "TFile.h" +#include "TTree.h" + +#include "fastjet/PseudoJet.hh" +#include "fastjet/ClusterSequenceArea.hh" + +#include "include/ProgressBar.h" + +#include "PU14/EventMixer.hh" +#include "PU14/CmdLine.hh" +#include "PU14/PU14.hh" + +#include "include/extraInfo.hh" +#include "include/jetCollection.hh" +#include "include/softDropGroomer.hh" +#include "include/treeWriter.hh" +#include "include/jetMatcher.hh" +#include "include/Angularity.hh" + +using namespace std; +using namespace fastjet; + +// ./runHybridJetAnalysis -hard samples/HYBRID_Hadrons_5020_dijet_K20_kappa_0.404_K_20.000.pu14 -nev 10 + +int main (int argc, char ** argv) { + + auto start_time = chrono::steady_clock::now(); + + CmdLine cmdline(argc,argv); + // inputs read from command line + int nEvent = cmdline.value("-nev",1); // first argument: command line option; second argument: default value + //bool verbose = cmdline.present("-verbose"); + + cout << "will run on " << nEvent << " events" << endl; + + // Uncomment to silence fastjet banner + ClusterSequence::set_fastjet_banner_stream(NULL); + + //to write info to root tree + treeWriter trw("jetTree"); + + //Jet definition + double R = 0.4; + double ghostRapMax = 6.0; + double ghost_area = 0.005; + int active_area_repeats = 1; + GhostedAreaSpec ghost_spec(ghostRapMax, active_area_repeats, ghost_area); + AreaDefinition area_def = AreaDefinition(active_area,ghost_spec); + JetDefinition jet_def(antikt_algorithm, R); + + double jetRapMax = 3.0; + Selector jet_selector = SelectorAbsRapMax(jetRapMax); + + Angularity width(1.,1.,R); + Angularity pTD(0.,2.,R); + + ProgressBar Bar(cout, nEvent); + Bar.SetStyle(-1); + + EventMixer mixer(&cmdline); //the mixing machinery from PU14 workshop + + // loop over events + int iev = 0; + unsigned int entryDiv = (nEvent > 200) ? nEvent / 200 : 1; + while ( mixer.next_event() && iev < nEvent ) + { + // increment event number + iev++; + + Bar.Update(iev); + Bar.PrintWithMod(entryDiv); + + vector particlesMerged = mixer.particles(); + + vector eventWeight; + eventWeight.push_back(mixer.hard_weight()); + eventWeight.push_back(mixer.pu_weight()); + + //get parton production points from hybrid file + vector prodX; + vector prodY; + prodX.push_back(mixer.productionX()); + prodY.push_back(mixer.productionY()); + + //std::cout << "X: " << mixer.productionX() << " Y: " << mixer.productionY() << std::endl; + + // extract hard partons that initiated the jets + fastjet::Selector parton_selector = SelectorVertexNumber(-1); + vector partons = parton_selector(particlesMerged); + + // select final state particles from hard event only + vector particlesBkg, particlesSig; + SelectorIsHard().sift(particlesMerged, particlesSig, particlesBkg); // this sifts the full event into two vectors of PseudoJet, one for the hard event, one for the underlying event + + //--------------------------------------------------------------------------- + // jet clustering + //--------------------------------------------------------------------------- + + fastjet::ClusterSequenceArea csSig(particlesSig, jet_def, area_def); + jetCollection jetCollectionSig(sorted_by_pt(jet_selector(csSig.inclusive_jets(10.)))); + + //calculate some angularities + vector widthSig; widthSig.reserve(jetCollectionSig.getJet().size()); + vector pTDSig; pTDSig.reserve(jetCollectionSig.getJet().size()); + for(PseudoJet jet : jetCollectionSig.getJet()) { + widthSig.push_back(width.result(jet)); + pTDSig.push_back(pTD.result(jet)); + } + jetCollectionSig.addVector("widthSig", widthSig); + jetCollectionSig.addVector("pTDSig", pTDSig); + + //--------------------------------------------------------------------------- + // Groom the jets + //--------------------------------------------------------------------------- + + //SoftDrop grooming classic for signal jets (zcut=0.1, beta=0) + softDropGroomer sdgSigBeta00Z01(0.1, 0.0, R); + jetCollection jetCollectionSigSDBeta00Z01(sdgSigBeta00Z01.doGrooming(jetCollectionSig)); + jetCollectionSigSDBeta00Z01.addVector("zgSigSDBeta00Z01", sdgSigBeta00Z01.getZgs()); + jetCollectionSigSDBeta00Z01.addVector("ndropSigSDBeta00Z01", sdgSigBeta00Z01.getNDroppedSubjets()); + jetCollectionSigSDBeta00Z01.addVector("dr12SigSDBeta00Z01", sdgSigBeta00Z01.getDR12()); + + + //--------------------------------------------------------------------------- + // write tree + //--------------------------------------------------------------------------- + + //Give variable we want to write out to treeWriter. + //Only vectors of the types 'jetCollection', and 'double', 'int', 'PseudoJet' are supported + + trw.addCollection("eventWeight", eventWeight); + trw.addCollection("prodX", prodX); + trw.addCollection("prodY", prodY); + trw.addPartonCollection("partons", partons); + + trw.addCollection("sigJet", jetCollectionSig); + trw.addCollection("sigJetSDBeta00Z01", jetCollectionSigSDBeta00Z01); + + trw.fillTree(); + + }//event loop + + Bar.Update(nEvent); + Bar.Print(); + Bar.PrintLine(); + + TTree *trOut = trw.getTree(); + + TFile *fout = new TFile(cmdline.value("-output", "JetToyHIResultHybridJetAnalysis.root").c_str(), "RECREATE"); + trOut->Write(); + fout->Write(); + fout->Close(); + + double time_in_seconds = chrono::duration_cast + (chrono::steady_clock::now() - start_time).count() / 1000.0; + cout << "runFromFile: " << time_in_seconds << endl; +}