Merge pull request #751 from mkirsano/decay_tauolapp_pi0

Decay tauolapp pi0

GeneratorInterface/Pythia8Interface/interface/Py8GunBase.h

@@ -1,20 +1,10 @@
-// -*- C++ -*-
-//
-//
-
-//
-// This class is a "Hadronizer" template (see GeneratorInterface/Core)
-// 
-
 #ifndef gen_Py8GunBase_h
 #define gen_Py8GunBase_h
 
 #include <memory>
 
 #include <boost/shared_ptr.hpp>
 
-#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
-
 #include "SimDataFormats/GeneratorProducts/interface/GenRunInfoProduct.h"
 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
 
@@ -77,4 +67,4 @@ namespace gen {
 
 } // namespace gen
 
-#endif // gen_BaseHadronizer_h
+#endif // gen_Py8GunBase_h

GeneratorInterface/Pythia8Interface/interface/Py8InterfaceBase.h

@@ -20,7 +20,6 @@ namespace gen {
 
          virtual bool generatePartonsAndHadronize() = 0;
          bool decay() { return true; } // NOT used - let's call it "design imperfection"
-         virtual bool residualDecay();
          bool readSettings( int ); // common func
          virtual bool initializeForInternalPartons() = 0;
          bool declareStableParticles( const std::vector<int>& ); // common func
@@ -44,4 +43,4 @@ namespace gen {
 
 }
 
-#endif
+#endif // gen_Py8InterfaceBase_h

GeneratorInterface/Pythia8Interface/plugins/Pythia8Hadronizer.cc

@@ -57,6 +57,9 @@ class Pythia8Hadronizer : public BaseHadronizer, public Py8InterfaceBase {
 
     bool generatePartonsAndHadronize();
     bool hadronize();
+
+    virtual bool residualDecay();
+
     void finalizeEvent();
 
     void statistics();
@@ -295,6 +298,12 @@ bool Pythia8Hadronizer::initializeForInternalPartons()
            fMasterGen->particleData.listAll();
   }
 
+  // init decayer
+  fDecayer->readString("ProcessLevel:all = off"); // trick
+  fDecayer->readString("Standalone:allowResDec=on");
+  // pythia->readString("ProcessLevel::resonanceDecays=on");
+  fDecayer->init();
+
   return true;
 }
 
@@ -335,6 +344,12 @@ bool Pythia8Hadronizer::initializeForExternalPartons()
            fMasterGen->particleData.listAll();
   }
 
+  // init decayer
+  fDecayer->readString("ProcessLevel:all = off"); // trick
+  fDecayer->readString("Standalone:allowResDec=on");
+  // pythia->readString("ProcessLevel::resonanceDecays=on");
+  fDecayer->init();
+
   return true;
 }
 
@@ -403,6 +418,82 @@ bool Pythia8Hadronizer::hadronize()
 }
 
 
+bool Pythia8Hadronizer::residualDecay()
+{
+
+  Event* pythiaEvent = &(fMasterGen->event);
+
+  int NPartsBeforeDecays = pythiaEvent->size();
+  int NPartsAfterDecays = event().get()->particles_size();
+  int NewBarcode = NPartsAfterDecays;
+
+  for ( int ipart=NPartsAfterDecays; ipart>NPartsBeforeDecays; ipart-- )
+  {
+
+    HepMC::GenParticle* part = event().get()->barcode_to_particle( ipart );
+
+    if ( part->status() == 1 )
+    {
+      fDecayer->event.reset();
+      Particle py8part(  part->pdg_id(), 93, 0, 0, 0, 0, 0, 0,
+                         part->momentum().x(),
+                         part->momentum().y(), 
+                         part->momentum().z(),
+                         part->momentum().t(),
+                         part->generated_mass() );
+      HepMC::GenVertex* ProdVtx = part->production_vertex();
+      py8part.vProd( ProdVtx->position().x(), ProdVtx->position().y(),
+                     ProdVtx->position().z(), ProdVtx->position().t() );
+      py8part.tau( (fDecayer->particleData).tau0( part->pdg_id() ) );
+      fDecayer->event.append( py8part );
+      int nentries = fDecayer->event.size();
+      if ( !fDecayer->event[nentries-1].mayDecay() ) continue;
+      fDecayer->next();
+      int nentries1 = fDecayer->event.size();
+      // fDecayer->event.list(std::cout);
+      if ( nentries1 <= nentries ) continue; //same number of particles, no decays...
+
+      part->set_status(2);
+
+      Particle& py8daughter = fDecayer->event[nentries]; // the 1st daughter
+      HepMC::GenVertex* DecVtx = new HepMC::GenVertex( HepMC::FourVector(py8daughter.xProd(),
+                                                       py8daughter.yProd(),
+                                                       py8daughter.zProd(),
+                                                       py8daughter.tProd()) );
+
+      DecVtx->add_particle_in( part ); // this will cleanup end_vertex if exists, replace with the new one
+                                       // I presume (vtx) barcode will be given automatically
+
+      HepMC::FourVector pmom( py8daughter.px(), py8daughter.py(), py8daughter.pz(), py8daughter.e() );
+
+      HepMC::GenParticle* daughter =
+                        new HepMC::GenParticle( pmom, py8daughter.id(), 1 );
+
+      NewBarcode++;
+      daughter->suggest_barcode( NewBarcode );
+      DecVtx->add_particle_out( daughter );
+
+      for ( int ipart1=nentries+1; ipart1<nentries1; ipart1++ )
+      {
+        py8daughter = fDecayer->event[ipart1];
+        HepMC::FourVector pmomN( py8daughter.px(), py8daughter.py(), py8daughter.pz(), py8daughter.e() );
+        HepMC::GenParticle* daughterN =
+                        new HepMC::GenParticle( pmomN, py8daughter.id(), 1 );
+        NewBarcode++;
+        daughterN->suggest_barcode( NewBarcode );
+        DecVtx->add_particle_out( daughterN );
+      }
+
+      event().get()->add_vertex( DecVtx );
+      // fCurrentEventState->add_vertex( DecVtx );
+
+    } 
+ }   
+ return true;
+
+}     
+
+
 void Pythia8Hadronizer::finalizeEvent()
 {
   bool lhe = lheEvent() != 0;

GeneratorInterface/Pythia8Interface/src/Py8GunBase.cc

@@ -1,4 +1,3 @@
-
 #include "GeneratorInterface/Pythia8Interface/interface/Py8GunBase.h"
 // #include "GeneratorInterface/Pythia8Interface/interface/RandomP8.h"
 // #include "GeneratorInterface/Core/interface/RNDMEngineAccess.h"

GeneratorInterface/Pythia8Interface/src/Py8InterfaceBase.cc

@@ -1,9 +1,3 @@
-//#include <iostream>
-//#include <sstream>
-//#include <string>
-//#include <memory>
-//#include <stdint.h>
-
 #include "GeneratorInterface/Pythia8Interface/interface/Py8InterfaceBase.h"
 #include "GeneratorInterface/Pythia8Interface/interface/RandomP8.h"
 #include "GeneratorInterface/Core/interface/RNDMEngineAccess.h"
@@ -94,88 +88,6 @@ bool Py8InterfaceBase:: declareSpecialSettings( const std::vector<std::string>&
    return true;
 }
 
-bool Py8InterfaceBase::residualDecay() 
-{
-
-/*
-  Event* pythiaEvent = &(fMasterPtr->event);
-  
-  assert(fCurrentEventState);
-  
-  int NPartsBeforeDecays = pythiaEvent->size();
-  // int NPartsAfterDecays = event().get()->particles_size();
-  int NPartsAfterDecays = fCurrentEventState->particles_size();
-  int NewBarcode = NPartsAfterDecays;
-   
-  for ( int ipart=NPartsAfterDecays; ipart>NPartsBeforeDecays; ipart-- )
-  {
-
-    // HepMC::GenParticle* part = event().get()->barcode_to_particle( ipart );
-    HepMC::GenParticle* part = fCurrentEventState->barcode_to_particle( ipart );
-
-    if ( part->status() == 1 )
-    {
-      fDecayerPtr->event.reset();
-      Particle py8part(  part->pdg_id(), 93, 0, 0, 0, 0, 0, 0,
-                         part->momentum().x(),
-                         part->momentum().y(),
-                         part->momentum().z(),
-                         part->momentum().t(),
-                         part->generated_mass() );
-      HepMC::GenVertex* ProdVtx = part->production_vertex();
-      py8part.vProd( ProdVtx->position().x(), ProdVtx->position().y(), 
-                     ProdVtx->position().z(), ProdVtx->position().t() );
-      py8part.tau( (fDecayerPtr->particleData).tau0( part->pdg_id() ) );
-      fDecayerPtr->event.append( py8part );
-      int nentries = fDecayerPtr->event.size();
-      if ( !fDecayerPtr->event[nentries-1].mayDecay() ) continue;
-      fDecayerPtr->next();
-      int nentries1 = fDecayerPtr->event.size();
-      // --> fDecayerPtr->event.list(std::cout);
-      if ( nentries1 <= nentries ) continue; //same number of particles, no decays...
-	    
-      part->set_status(2);
-	    
-      Particle& py8daughter = fDecayerPtr->event[nentries]; // the 1st daughter
-      HepMC::GenVertex* DecVtx = new HepMC::GenVertex( HepMC::FourVector(py8daughter.xProd(),
-                                                       py8daughter.yProd(),
-                                                       py8daughter.zProd(),
-                                                       py8daughter.tProd()) );
-
-      DecVtx->add_particle_in( part ); // this will cleanup end_vertex if exists, replace with the new one
-                                       // I presume (vtx) barcode will be given automatically
-	    
-      HepMC::FourVector pmom( py8daughter.px(), py8daughter.py(), py8daughter.pz(), py8daughter.e() );
-	    
-      HepMC::GenParticle* daughter =
-                        new HepMC::GenParticle( pmom, py8daughter.id(), 1 );
-	    
-      NewBarcode++;
-      daughter->suggest_barcode( NewBarcode );
-      DecVtx->add_particle_out( daughter );
-	    	    
-      for ( ipart=nentries+1; ipart<nentries1; ipart++ )
-      {
-        py8daughter = fDecayerPtr->event[ipart];
-        HepMC::FourVector pmomN( py8daughter.px(), py8daughter.py(), py8daughter.pz(), py8daughter.e() );	    
-        HepMC::GenParticle* daughterN =
-                        new HepMC::GenParticle( pmomN, py8daughter.id(), 1 );
-        NewBarcode++;
-        daughterN->suggest_barcode( NewBarcode );
-        DecVtx->add_particle_out( daughterN );
-      }
-	    
-      // event().get()->add_vertex( DecVtx );
-      fCurrentEventState->add_vertex( DecVtx );
-
-    }
- } 
-*/   
- return true;
-
-}
-
-
 void Py8InterfaceBase::statistics()
 {
 

GeneratorInterface/Pythia8Interface/test/BuildFile.xml

@@ -0,0 +1,14 @@
+<use name="FWCore/Framework"/>
+<use name="FWCore/ParameterSet"/>
+<use name="root"/>
+<use name="SimDataFormats/GeneratorProducts"/> 
+<use name="GeneratorInterface/Core"/>
+<use name="CommonTools/UtilAlgos"/>
+
+#<library file="analyserhepmc/LeptonAnalyserHepMC.cc" name="AnalyserHepMC">
+#  <flags EDM_PLUGIN="1"/>
+#</library>
+
+<library file="ZJetsAnalyzer.cc,analyserhepmc/LeptonAnalyserHepMC.cc,analyserhepmc/JetInputHepMC.cc" name="ZJetsTestAnalyzer">
+  <flags EDM_PLUGIN="1"/>
+</library>