HerwigToHepMCProduction/HerwigHepMCFilter/RunHerwigHepMCFilter.C

0001 #pragma once
0002 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,00,0)
0003 #include <hepmcjettrigger/HepMCJetTrigger.h>
0004 #include <hepmcjettrigger/HepMCParticleTrigger.h>
0005 #include "sPhenixStyle.h"
0006 #include "sPhenixStyle.C"
0007 #include <fun4all/Fun4AllInputManager.h>
0008 #include <fun4all/Fun4AllServer.h>
0009 #include <phhepmc/Fun4AllHepMCInputManager.h>
0010 #include <fun4all/Fun4AllDstInputManager.h>
0011 //#include <fun4allraw/Fun4AllPrdfInputManager.h>
0012 //#include <fun4allraw/Fun4AllPrdfInputPoolManager.h>
0013 //#include <fun4all/Fun4AllDstOutputManager.h>
0014 #include <phhepmc/Fun4AllHepMCOutputManager.h>
0015 #include <fun4all/SubsysReco.h>
0016 #include <phool/PHRandomSeed.h>
0017 #include <sstream>
0018 #include <string>
0019 #include <algorithm>
0020 #include <fstream>
0021 #include <cmath>
0022
0023 R__LOAD_LIBRARY(libfun4all.so);
0024 //R__LOAD_LIBRARY(libfun4allraw.so);
0025 R__LOAD_LIBRARY(libHepMCJetTrigger.so);
0026 R__LOAD_LIBRARY(libHepMCParticleTrigger.so);
0027 R__LOAD_LIBRARY(libffamodules.so);
0028 R__LOAD_LIBRARY(libffarawmodules.so);
0029 R__LOAD_LIBRARY(libphhepmc.so);
0030
0031 int RunHerwigHepMCFilter(std::string filename="/sphenix/user/sgross/sphenix_herwig/herwig_files/sphenix_10GeV_jetpt.hepmc", std::string trigger_type="jet", std::string trig="10", int goal_event_number=1000, int nGen=1000000, std::string xs_file = "none")
0032 {
0033     float threshold=0.;
0034     try{
0035         threshold=std::stof(trig);
0036     }
0037     catch(std::exception& e){threshold=10.;}
0038     Fun4AllServer* se=Fun4AllServer::instance();
0039     std::string outfile=filename;
0040     int segment=-1;
0041         auto first_break = filename.find("-");
0042     auto second_break = filename.find(".he");
0043     std::string seg = filename.substr(first_break + 1 , second_break - first_break -1 );
0044     try{
0045         segment = std::stoi(seg);
0046     }
0047     catch(std::exception& e){}
0048     outfile.insert(filename.find("-"),"_filtered");
0049     Fun4AllHepMCInputManager *in =new Fun4AllHepMCInputManager("in");
0050     Fun4AllHepMCOutputManager *out=new Fun4AllHepMCOutputManager("out", outfile);
0051     HepMCJetTrigger* hf=new HepMCJetTrigger(threshold, goal_event_number, true);
0052     HepMCParticleTrigger* part = new HepMCParticleTrigger(threshold, goal_event_number, true);
0053     if(trigger_type.find("photon") !=std::string::npos){
0054         part->AddParticle(22);
0055     }
0056     se->registerInputManager(in);
0057     in->fileopen(filename);
0058     se->registerOutputManager(out);
0059     if(trigger_type.find("jet") != std::string::npos) se->registerSubsystem(hf);
0060     if(trigger_type.find("photon") != std::string::npos) se->registerSubsystem(part);
0061     if(trigger_type.find("none") == std::string::npos)
0062     {
0063         se->run();
0064         se->End();
0065     }
0066     int n_good, n_evts;
0067
0068     if(trigger_type.find("jet") != std::string::npos ){
0069         n_evts = hf->getNevts();
0070         n_good=hf->getNgood();
0071     }
0072
0073     else if(trigger_type.find("photon") != std::string::npos ){
0074         n_evts = part->getNevts();
0075         n_good=part->getNgood();
0076     }
0077     else {
0078            n_good=goal_event_number;
0079            n_evts=nGen;
0080     }
0081     std::cout<<"Analyzed events = " <<n_evts <<"\n Good events= " <<n_good <<std::endl;
0082     //get the measured cross section
0083     std::fstream f;
0084     f.open(xs_file);
0085     float xs_value = 0.;
0086     float xs_pow = 0.;
0087     std::string line;
0088     while(std::getline(f, line))
0089     {
0090         if(line.find("generated events") != std::string::npos)
0091         {
0092             std::stringstream linestream (line);
0093             std::string subline;
0094             while(std::getline(linestream, subline, ' '))
0095             {
0096                 if(subline.find(")e")!=std::string::npos)
0097                 {
0098                     auto left_par = subline.find("(");
0099                     xs_value=std::stof(subline.substr(0, left_par-1));
0100                     auto pow_e = subline.find("e");
0101                     xs_pow=std::stof(subline.substr(pow_e+2));
0102                 }
0103             }
0104         }
0105     }
0106     float xs = xs_value * std::pow(10, xs_pow) * ((float)n_good)/((float)n_evts*(float)nGen);
0107     f.close();
0108     std::ofstream of;
0109     std::string xs_of="cross_section_data/Cross_Section_"+trigger_type+trig+"_"+std::to_string(segment)+".txt";
0110     of.open(xs_of.c_str());
0111     if(!of.good()) std::cout<<"failed to open file" <<std::endl;
0112     if(of.good())
0113         of<<"Events analyzed= " << n_evts<<"\n Events passing filter= " <<n_good <<"\n Total Generated events= " <<nGen <<"\n XS/N= "<<xs
0114         <<"\n Read xs = " <<xs_value <<" x 10^ " <<xs_pow <<" nb " <<std::endl;
0115     else
0116         std::cout<<"Events analyzed= " << n_evts<<"\n Events passing filter= " <<n_good <<"\n Total Generated events= " <<nGen <<"\n XS/N= "<<xs
0117             <<"\n Read xs = " <<xs_value <<" x 10^ " <<xs_pow <<" nb " <<std::endl;
0118     of.close();
0119     std::cout<<"Cross section / Evt = " <<xs<<std::endl;
0120
0121     return 0;
0122 }
0123 #endif