Merge pull request #14 from mverwe/prodHybrid

get parton production points in hybrid events

Author: mverwe

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;
     }
   }
 

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<fastjet::PseudoJet> _particles;
   double _hard_event_weight, _pu_event_weight;
+  double _posX, _posY;
 };
 
 #endif  // __EVENTMIXER_HH__

PU14/EventSource.cc

@@ -30,7 +30,8 @@ void EventSource::open_stream(const std::string & filename) {
 
 //----------------------------------------------------------------------
 bool EventSource::append_next_event(std::vector<fastjet::PseudoJet> & 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<fastjet::PseudoJet> & 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<fastjet::PseudoJet> & 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());

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<fastjet::PseudoJet> & particles,
-        double &event_weight, int vertex_number = 0);
+                         double &event_weight, double &prodX, double &prodY,
+                         int vertex_number = 0);
 
 private:
   std::istream * _stream;

runHybridJetAnalysis.cc

@@ -0,0 +1,160 @@
+#include <iostream>
+#include <chrono>
+
+#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<int>("-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<PseudoJet> particlesMerged = mixer.particles();
+
+    vector<double> eventWeight;
+    eventWeight.push_back(mixer.hard_weight());
+    eventWeight.push_back(mixer.pu_weight());
+
+    //get parton production points from hybrid file
+    vector<double> prodX;
+    vector<double> 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<PseudoJet> partons = parton_selector(particlesMerged);
+    
+    // select final state particles from hard event only
+    vector<PseudoJet> 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<double> widthSig; widthSig.reserve(jetCollectionSig.getJet().size());
+    vector<double> 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<string>("-output", "JetToyHIResultHybridJetAnalysis.root").c_str(), "RECREATE");
+  trOut->Write();
+  fout->Write();
+  fout->Close();
+
+  double time_in_seconds = chrono::duration_cast<chrono::milliseconds>
+    (chrono::steady_clock::now() - start_time).count() / 1000.0;
+  cout << "runFromFile: " << time_in_seconds << endl;
+}