HerwigToHepMCProduction/HepMCTrigger/HepMCJetTrigger.cc

0001 //____________________________________________________________________________..
0002 //
0003 // This is a template for a Fun4All SubsysReco module with all methods from the
0004 // $OFFLINE_MAIN/include/fun4all/SubsysReco.h baseclass
0005 // You do not have to implement all of them, you can just remove unused methods
0006 // here and in HepMCJetTrigger.h.
0007 //
0008 // HepMCJetTrigger(const std::string &name = "HepMCJetTrigger")
0009 // everything is keyed to HepMCJetTrigger, duplicate names do work but it makes
0010 // e.g. finding culprits in logs difficult or getting a pointer to the module
0011 // from the command line
0012 //
0013 // HepMCJetTrigger::~HepMCJetTrigger()
0014 // this is called when the Fun4AllServer is deleted at the end of running. Be
0015 // mindful what you delete - you do loose ownership of object you put on the node tree
0016 //
0017 // int HepMCJetTrigger::Init(PHCompositeNode *topNode)
0018 // This method is called when the module is registered with the Fun4AllServer. You
0019 // can create historgrams here or put objects on the node tree but be aware that
0020 // modules which haven't been registered yet did not put antyhing on the node tree
0021 //
0022 // int HepMCJetTrigger::InitRun(PHCompositeNode *topNode)
0023 // This method is called when the first event is read (or generated). At
0024 // this point the run number is known (which is mainly interesting for raw data
0025 // processing). Also all objects are on the node tree in case your module's action
0026 // depends on what else is around. Last chance to put nodes under the DST Node
0027 // We mix events during readback if branches are added after the first event
0028 //
0029 // int HepMCJetTrigger::process_event(PHCompositeNode *topNode)
0030 // called for every event. Return codes trigger actions, you find them in
0031 // $OFFLINE_MAIN/include/fun4all/Fun4AllReturnCodes.h
0032 //   everything is good:
0033 //     return Fun4AllReturnCodes::EVENT_OK
0034 //   abort event reconstruction, clear everything and process next event:
0035 //     return Fun4AllReturnCodes::ABORT_EVENT;
0036 //   proceed but do not save this event in output (needs output manager setting):
0037 //     return Fun4AllReturnCodes::DISCARD_EVENT;
0038 //   abort processing:
0039 //     return Fun4AllReturnCodes::ABORT_RUN
0040 // all other integers will lead to an error and abort of processing
0041 //
0042 // int HepMCJetTrigger::ResetEvent(PHCompositeNode *topNode)
0043 // If you have internal data structures (arrays, stl containers) which needs clearing
0044 // after each event, this is the place to do that. The nodes under the DST node are cleared
0045 // by the framework
0046 //
0047 // int HepMCJetTrigger::EndRun(const int runnumber)
0048 // This method is called at the end of a run when an event from a new run is
0049 // encountered. Useful when analyzing multiple runs (raw data). Also called at
0050 // the end of processing (before the End() method)
0051 //
0052 // int HepMCJetTrigger::End(PHCompositeNode *topNode)
0053 // This is called at the end of processing. It needs to be called by the macro
0054 // by Fun4AllServer::End(), so do not forget this in your macro
0055 //
0056 // int HepMCJetTrigger::Reset(PHCompositeNode *topNode)
0057 // not really used - it is called before the dtor is called
0058 //
0059 // void HepMCJetTrigger::Print(const std::string &what) const
0060 // Called from the command line - useful to print information when you need it
0061 //
0062 //____________________________________________________________________________..
0063
0064 #include "HepMCJetTrigger.h"
0065
0066 #include <fun4all/Fun4AllReturnCodes.h>
0067
0068 #include <phool/PHCompositeNode.h>
0069
0070 //____________________________________________________________________________..
0071 HepMCJetTrigger::HepMCJetTrigger(float trigger_thresh, int n_incom, bool up_lim, const std::string &name):
0072  SubsysReco(name)
0073 {
0074   std::cout << "HepMCJetTrigger::HepMCJetTrigger(const std::string &name) Calling ctor" << std::endl;
0075     this->threshold=trigger_thresh;
0076     this->goal_event_number=n_incom;
0077     this->set_event_limit=up_lim;
0078 }
0079
0080
0081 //____________________________________________________________________________..
0082 HepMCJetTrigger::~HepMCJetTrigger()
0083 {
0084   std::cout << "HepMCJetTrigger::~HepMCJetTrigger() Calling dtor" << std::endl;
0085 }
0086
0087 //____________________________________________________________________________..
0088 int HepMCJetTrigger::Init(PHCompositeNode *topNode)
0089 {
0090   std::cout << "HepMCJetTrigger::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0091   return Fun4AllReturnCodes::EVENT_OK;
0092 }
0093
0094 //____________________________________________________________________________..
0095 int HepMCJetTrigger::InitRun(PHCompositeNode *topNode)
0096 {
0097   std::cout << "HepMCJetTrigger::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
0098   return Fun4AllReturnCodes::EVENT_OK;
0099 }
0100
0101 //____________________________________________________________________________..
0102 int HepMCJetTrigger::process_event(PHCompositeNode *topNode)
0103 {
0104   //std::cout << "HepMCJetTrigger::process_event(PHCompositeNode *topNode) Processing Event" << std::endl;
0105     n_evts++;
0106     if(this->set_event_limit == true){ //needed to keep all HepMC output at the same number of events
0107         if(n_good >= this->goal_event_number) return Fun4AllReturnCodes::ABORTEVENT;
0108     }
0109     PHHepMCGenEventMap* phg=findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
0110     if(!phg){
0111         return Fun4AllReturnCodes::ABORTEVENT;
0112     }
0113     for( PHHepMCGenEventMap::ConstIter eventIter=phg->begin(); eventIter != phg->end(); ++eventIter){
0114         PHHepMCGenEvent* hepev=eventIter->second;
0115         if(!hepev) return Fun4AllReturnCodes::ABORTEVENT;
0116         HepMC::GenEvent* ev=hepev->getEvent();
0117         if(!ev) return Fun4AllReturnCodes::ABORTEVENT;
0118         bool good_event=isGoodEvent(ev);
0119         if (good_event) n_good++;
0120         if(!good_event) return Fun4AllReturnCodes::ABORTEVENT;
0121     }
0122     return Fun4AllReturnCodes::EVENT_OK;
0123 }
0124
0125 //____________________________________________________________________________..
0126 int HepMCJetTrigger::ResetEvent(PHCompositeNode *topNode)
0127 {
0128   std::cout << "HepMCJetTrigger::ResetEvent(PHCompositeNode *topNode) Resetting internal structures, prepare for next event" << std::endl;
0129   return Fun4AllReturnCodes::EVENT_OK;
0130 }
0131
0132 //____________________________________________________________________________..
0133 int HepMCJetTrigger::EndRun(const int runnumber)
0134 {
0135   std::cout << "HepMCJetTrigger::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
0136   return Fun4AllReturnCodes::EVENT_OK;
0137 }
0138
0139 //____________________________________________________________________________..
0140 int HepMCJetTrigger::End(PHCompositeNode *topNode)
0141 {
0142   std::cout << "HepMCJetTrigger::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0143   return Fun4AllReturnCodes::EVENT_OK;
0144 }
0145
0146 //____________________________________________________________________________..
0147 int HepMCJetTrigger::Reset(PHCompositeNode *topNode)
0148 {
0149  std::cout << "HepMCJetTrigger::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
0150   return Fun4AllReturnCodes::EVENT_OK;
0151 }
0152
0153 //____________________________________________________________________________..
0154 void HepMCJetTrigger::Print(const std::string &what) const
0155 {
0156   std::cout << "HepMCJetTrigger::Print(const std::string &what) const Printing info for " << what << std::endl;
0157 }
0158 bool HepMCJetTrigger::isGoodEvent( HepMC::GenEvent* e1)
0159 {
0160     //this is really just the call to actually evaluate and return the filter
0161     if(this->threshold == 0 ) return true;
0162     std::vector<fastjet::PseudoJet> jets=findAllJets(e1);
0163     int njetsabove=jetsAboveThreshold(jets);
0164     if(njetsabove > 0 ) return true;
0165     else return false;
0166 }
0167
0168 std::vector<fastjet::PseudoJet> HepMCJetTrigger::findAllJets(HepMC::GenEvent* e1)
0169 {
0170     //do the fast jet clustering, antikt r=-0.4
0171     fastjet::JetDefinition jetdef(fastjet::antikt_algorithm, 0.4);
0172     std::vector<fastjet::PseudoJet> input, output;
0173     for(HepMC::GenEvent::particle_const_iterator iter = e1->particles_begin(); iter != e1->particles_end(); ++iter)
0174     {
0175         if(!(*iter)->end_vertex() && (*iter)->status() == 1){
0176             auto p=(*iter)->momentum();
0177             fastjet::PseudoJet pj( p.px(), p.py(), p.pz(), p.e());
0178             pj.set_user_index((*iter)->barcode());
0179             input.push_back(pj);
0180         }
0181     }
0182     if(input.size() == 0 ) return input;
0183     fastjet::ClusterSequence js (input, jetdef);
0184     auto j = js.inclusive_jets();
0185     for(auto j1:j){
0186         output.push_back(j1); //just keep in the corect format
0187     }
0188     return output;
0189 }
0190 int HepMCJetTrigger::jetsAboveThreshold(std::vector<fastjet::PseudoJet> jets)
0191 {
0192     //search through for the number of identified jets above the threshold
0193     int n_good_jets=0;
0194     for(auto j:jets)
0195     {
0196         float pt=j.pt();
0197         if(pt > this->threshold) n_good_jets++;
0198     }
0199     return n_good_jets;
0200 }
0201