Multiple Hadronization

GeneratorInterface/Core/interface/BaseHepMCFilter.h

@@ -0,0 +1,26 @@
+#ifndef BaseHepMCFilter_H
+#define BaseHepMCFilter_H
+
+// J.Bendavid
+
+// base class for simple filter to run inside of HadronizerFilter for
+// multiple hadronization attempts on lhe events
+
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
+
+//
+// class declaration
+//
+
+class BaseHepMCFilter  {
+public:
+  BaseHepMCFilter();
+  virtual ~BaseHepMCFilter();
+
+
+  virtual bool filter(const HepMC::GenEvent* evt)  = 0;
+
+};
+
+#endif

GeneratorInterface/Core/interface/GenXSecAnalyzer.h

@@ -46,16 +46,23 @@ class GenXSecAnalyzer : public edm::one::EDAnalyzer<edm::one::WatchLuminosityBlo
   void compute();
 
   edm::EDGetTokenT<GenFilterInfo> genFilterInfoToken_;
+  edm::EDGetTokenT<GenFilterInfo> hepMCFilterInfoToken_;
   edm::EDGetTokenT<GenLumiInfoProduct> genLumiInfoToken_;
 
   // ----------member data --------------------------
 
   int hepidwtup_;
   unsigned int theProcesses_size;
+  bool           hasHepMCFilterInfo_;
+
   // final cross sections
   GenLumiInfoProduct::XSec xsec_;
   // statistics from additional generator filter
   GenFilterInfo  filterOnlyEffStat_;     
+
+  // statistics from HepMC filter
+  GenFilterInfo  hepMCFilterEffStat_;   
+
   // statistics for event level efficiency, the size is the number of processes + 1 
   std::vector<GenFilterInfo>  eventEffStat_; 
   // statistics from jet matching, the size is the number of processes + 1 

GeneratorInterface/Core/interface/GenericDauHepMCFilter.h

@@ -0,0 +1,56 @@
+#ifndef GENERICDAUHEPMCFILTER_h
+#define GENERICDAUHEPMCFILTER_h
+// -*- C++ -*-
+//
+// Package:    GenericDauHepMCFilter
+// Class:      GenericDauHepMCFilter
+//
+/**\class GenericDauHepMCFilter GenericDauHepMCFilter.cc
+
+Description: Filter events using MotherId and ChildrenIds infos
+
+Implementation:
+<Notes on implementation>
+*/
+//
+// Original Author:  Daniele Pedrini
+//         Created:  Apr 29 2008
+// $Id: GenericDauHepMCFilter.h,v 1.2 2010/07/21 04:23:24 wmtan Exp $
+//
+//
+
+
+// system include files
+#include <memory>
+
+// user include files
+#include "GeneratorInterface/Core/interface/BaseHepMCFilter.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+
+
+//
+// class decleration
+//
+
+class GenericDauHepMCFilter : public BaseHepMCFilter {
+  public:
+     GenericDauHepMCFilter(const edm::ParameterSet&);
+    ~GenericDauHepMCFilter();
+
+    virtual bool filter(const HepMC::GenEvent* evt);
+
+  private:
+    // ----------memeber function----------------------
+
+    // ----------member data ---------------------------
+
+    int particleID;
+    bool chargeconju;
+    int ndaughters;
+    std::vector<int> dauIDs;    
+    double minptcut;
+    double maxptcut;
+    double minetacut;
+    double maxetacut;
+    };
+#endif

GeneratorInterface/Core/interface/HadronizerFilter.h

@@ -31,6 +31,8 @@
 
 // #include "GeneratorInterface/ExternalDecays/interface/ExternalDecayDriver.h"
 
+#include "GeneratorInterface/Core/interface/HepMCFilterDriver.h"
+
 // LHE Run
 #include "SimDataFormats/GeneratorProducts/interface/LHERunInfoProduct.h"
 #include "GeneratorInterface/LHEInterface/interface/LHERunInfo.h"
@@ -39,6 +41,7 @@
 #include "SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h"
 #include "GeneratorInterface/LHEInterface/interface/LHEEvent.h"
 
+
 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
 #include "SimDataFormats/GeneratorProducts/interface/GenRunInfoProduct.h"
 #include "SimDataFormats/GeneratorProducts/interface/GenLumiInfoProduct.h"
@@ -75,9 +78,11 @@ namespace edm
     Hadronizer hadronizer_;
     // gen::ExternalDecayDriver* decayer_;
     Decayer* decayer_;
+    HepMCFilterDriver *filter_;
     bool fromSource_;
     InputTag runInfoProductTag_;
     unsigned int counterRunInfoProducts_;
+    unsigned int nAttempts_;
   };
 
   //------------------------------------------------------------------------
@@ -89,9 +94,11 @@ namespace edm
     EDFilter(),
     hadronizer_(ps),
     decayer_(0),
+    filter_(0),
     fromSource_(true),
     runInfoProductTag_(),
-    counterRunInfoProducts_(0)
+    counterRunInfoProducts_(0),
+    nAttempts_(1)
   {
     callWhenNewProductsRegistered([this]( BranchDescription const& iBD) {
       //this is called each time a module registers that it will produce a LHERunInfoProduct
@@ -128,6 +135,17 @@ namespace edm
          usesResource(resource);
        }
     }
+
+    if ( ps.exists("HepMCFilter") ) {
+      ParameterSet psfilter = ps.getParameter<ParameterSet>("HepMCFilter");
+      filter_ = new HepMCFilterDriver(psfilter);
+    }
+
+    //initialize setting for multiple hadronization attempts
+    if (ps.exists("nAttempts")) {
+      nAttempts_ = ps.getParameter<unsigned int>("nAttempts");
+    }
+
     // This handles the case where there are no shared resources, because you
     // have to declare something when the SharedResources template parameter was used.
     if(sharedResources.empty() && (!decayer_ || decayer_->sharedResources().empty())) {
@@ -138,11 +156,16 @@ namespace edm
     produces<GenEventInfoProduct>();
     produces<GenLumiInfoProduct, edm::InLumi>();
     produces<GenRunInfoProduct, edm::InRun>();
+    if(filter_)
+      produces<GenFilterInfo, edm::InLumi>(); 
   }
 
   template <class HAD, class DEC>
   HadronizerFilter<HAD,DEC>::~HadronizerFilter()
-  { if (decayer_) delete decayer_; }
+  {
+    if (decayer_) delete decayer_;
+    if (filter_) delete filter_;
+  }
 
   template <class HAD, class DEC>
   bool
@@ -160,59 +183,93 @@ namespace edm
       ev.getByLabel("source", product);
     else
       ev.getByLabel("externalLHEProducer", product);
-
-
-    lhef::LHEEvent *lheEvent =
-		new lhef::LHEEvent(hadronizer_.getLHERunInfo(), *product);
-    hadronizer_.setLHEEvent( lheEvent );
 
-    // hadronizer_.generatePartons();
-    if ( !hadronizer_.hadronize() ) return false ;
-
-    //  this is "fake" stuff
-    // in principle, decays are done as part of full event generation,
-    // except for particles that are marked as to be kept stable
-    // but we currently keep in it the design, because we might want
-    // to use such feature for other applications
-    //
-    if ( !hadronizer_.decay() ) return false;
+    std::auto_ptr<HepMC::GenEvent> finalEvent;
+    std::auto_ptr<GenEventInfoProduct> finalGenEventInfo;
 
-    std::auto_ptr<HepMC::GenEvent> event (hadronizer_.getGenEvent());
-    if( !event.get() ) return false; 
+    //sum of weights for events passing hadronization
+    double waccept = 0;
+    //number of accepted events
+    unsigned int naccept = 0;
+
+    for (unsigned int itry = 0; itry<nAttempts_; ++itry) {
+
+      lhef::LHEEvent *lheEvent =
+		  new lhef::LHEEvent(hadronizer_.getLHERunInfo(), *product);
+      hadronizer_.setLHEEvent( lheEvent );
+
+      // hadronizer_.generatePartons();
+      if ( !hadronizer_.hadronize() ) continue ;
+
+      //  this is "fake" stuff
+      // in principle, decays are done as part of full event generation,
+      // except for particles that are marked as to be kept stable
+      // but we currently keep in it the design, because we might want
+      // to use such feature for other applications
+      //
+      if ( !hadronizer_.decay() ) continue;
+
+      std::auto_ptr<HepMC::GenEvent> event (hadronizer_.getGenEvent());
+      if( !event.get() ) continue; 
+
+      // The external decay driver is being added to the system,
+      // it should be called here
+      //
+      if ( decayer_ ) 
+      {
+        event.reset( decayer_->decay( event.get(),lheEvent ) );
+      }
 
-    // The external decay driver is being added to the system,
-    // it should be called here
-    //
-    if ( decayer_ ) 
-    {
-      event.reset( decayer_->decay( event.get(),lheEvent ) );
+      if ( !event.get() ) continue;
+
+      // check and perform if there're any unstable particles after 
+      // running external decay packges
+      //
+      hadronizer_.resetEvent( event.release() );
+      if ( !hadronizer_.residualDecay() ) continue;
+
+      hadronizer_.finalizeEvent();
+
+      event.reset( hadronizer_.getGenEvent() );
+      if ( !event.get() ) continue;
+
+      event->set_event_number( ev.id().event() );
+
+      std::auto_ptr<GenEventInfoProduct> genEventInfo(hadronizer_.getGenEventInfo());      
+      if (!genEventInfo.get())
+      { 
+        // create GenEventInfoProduct from HepMC event in case hadronizer didn't provide one
+        genEventInfo.reset(new GenEventInfoProduct(event.get()));
+      }
+
+      //if HepMCFilter was specified, test event
+      if (filter_ && !filter_->filter(event.get(), genEventInfo->weight())) continue;      
+
+      waccept += genEventInfo->weight();
+      ++naccept;
+
+      //keep the LAST accepted event (which is equivalent to choosing randomly from the accepted events)
+      finalEvent.reset(event.release());
+      finalGenEventInfo.reset(genEventInfo.release());
+
     }
+
+    if (!naccept) return false;
+
 
-    if ( !event.get() ) return false;
-
-    // check and perform if there're any unstable particles after 
-    // running external decay packges
-    //
-    hadronizer_.resetEvent( event.release() );
-    if ( !hadronizer_.residualDecay() ) return false;
-
-    hadronizer_.finalizeEvent();
-
-    event.reset( hadronizer_.getGenEvent() );
-    if ( !event.get() ) return false;
-
-    event->set_event_number( ev.id().event() );
-
-    std::auto_ptr<GenEventInfoProduct> genEventInfo(hadronizer_.getGenEventInfo());
-    if (!genEventInfo.get())
-    { 
-      // create GenEventInfoProduct from HepMC event in case hadronizer didn't provide one
-      genEventInfo.reset(new GenEventInfoProduct(event.get()));
+
+    //adjust event weights if necessary (in case input event was accepted multiple times)
+    if (naccept>1) {
+      std::vector<double> genEventInfoWeights = finalGenEventInfo->weights();
+      genEventInfoWeights.push_back(double(naccept)/double(nAttempts_));
+      finalGenEventInfo->setWeights(genEventInfoWeights);
     }
-    ev.put(genEventInfo);
+
+
+    ev.put(finalGenEventInfo);
 
     std::auto_ptr<HepMCProduct> bare_product(new HepMCProduct());
-    bare_product->addHepMCData( event.release() );
+    bare_product->addHepMCData( finalEvent.release() );
     ev.put(bare_product);
 
     return true;
@@ -267,6 +324,7 @@ namespace edm
 
     hadronizer_.statistics();
     if ( decayer_ ) decayer_->statistics();
+    if ( filter_ ) filter_->statistics();
     lheRunInfo->statistics();
 
     std::auto_ptr<GenRunInfoProduct> griproduct( new GenRunInfoProduct(genRunInfo) );
@@ -302,6 +360,10 @@ namespace edm
             << hadronizer_.classname()
             << "\n";
     }
+
+    if (filter_) {
+      filter_->resetStatistics();
+    }
 
     if (! hadronizer_.initializeForExternalPartons())
       throw edm::Exception(errors::Configuration) 
@@ -347,6 +409,23 @@ namespace edm
     lumi.put(genLumiInfo);
 
 
+    // produce GenFilterInfo if HepMCFilter is called
+    if (filter_) {
+
+      std::auto_ptr<GenFilterInfo> thisProduct(new GenFilterInfo(
+								 filter_->numEventsPassPos(),
+								 filter_->numEventsPassNeg(),
+								 filter_->numEventsTotalPos(),
+								 filter_->numEventsTotalNeg(),
+								 filter_->sumpass_w(),
+								 filter_->sumpass_w2(),
+								 filter_->sumtotal_w(),
+								 filter_->sumtotal_w2()
+								 ));
+      lumi.put(thisProduct);
+    }
+
+
   }