tauola_photos_pythia_example.cxx

Go to the documentation of this file.

//pythia header files
#include "Pythia.h"
#include "HepMCInterface.h"

//MC-TESTER header files
#include "Generate.h"
#include "HepMCEvent.H"
#include "Setup.H"

//PHOTOS header files
#include "Photos/Photos.h"
#include "Photos/PhotosHepMCEvent.h"
#include "Photos/Log.h"

//TAUOLA header files
#include "Tauola/Tauola.h"
#include "Tauola/TauolaHepMCEvent.h"

using namespace std;
using namespace Pythia8;
using namespace Photospp;
using namespace Tauolapp;

unsigned long NumberOfEvents = 10000;
unsigned int EventsToCheck=20;

// elementary test of HepMC typically executed before
// detector simulation based on http://home.fnal.gov/~mrenna/HCPSS/HCPSShepmc.html
// similar test was performed in Fortran
// we perform it before and after Photos (for the first several events)
void checkMomentumConservationInEvent(HepMC::GenEvent *evt)
{
        //cout<<"List of stable particles: "<<endl;

        double px=0.0,py=0.0,pz=0.0,e=0.0;
        
        for ( HepMC::GenEvent::particle_const_iterator p = evt->particles_begin();
              p != evt->particles_end(); ++p )
        {
                if( (*p)->status() == 1 )
                {
                        HepMC::FourVector m = (*p)->momentum();
                        px+=m.px();
                        py+=m.py();
                        pz+=m.pz();
                        e +=m.e();
                        //(*p)->print();
                }
        }
  cout.precision(6);
  cout.setf(ios_base::floatfield);
        cout<<endl<<"Vector Sum: "<<px<<" "<<py<<" "<<pz<<" "<<e<<endl;
}

int main(int argc,char **argv)
{
        HepMC::I_Pythia8 ToHepMC;

        // Initialization of pythia
        Pythia pythia;
        Event& event = pythia.event;

        pythia.readString("PartonLevel:ISR = off");
        pythia.readString("PartonLevel:FSR = off");

        pythia.readString("WeakSingleBoson:ffbar2gmZ = on");
        pythia.readString("23:onMode = off");
        pythia.readString("23:onIfAny = 15");
        pythia.particleData.readString("15:mayDecay = off"); //<- uncomment for pythia+tauola

        //pythia.init( -2212, -2212, 14000.0);     //proton proton collisions
        pythia.init( 11, -11, 91.187);             //electron positron collisions

        // TAUOLA and PHOTOS initialization
        Tauola::initialize();
        Photos::initialize();

        Photos::setInfraredCutOff(0.01/91.187); // 10MeV for scale to M_Z=91.187
        //Photos::setDoubleBrem(false);
        //Photos::setExponentiation(false);

        Log::SummaryAtExit();
        cout.setf(ios::fixed);
        //Log::LogInfo(false) //To turn printing of last five events and pythia statistics off

        // Example setup - suppress processing of whole Z0 decay,
        // leaving only the Z0 -> tau+ tau- decay and whole branch starting
        // from tau- to be processed
        //Photos::suppressBremForBranch(0,23);
        //Photos::forceBremForDecay (2,23,15,-15);
        //Photos::forceBremForBranch(0,15);

        // Force mass of electron and positron to be 0.000511
        //Photos::forceMass(11,0.000511);

        // Force mass of electron and positron to be taken
        // from event record instead of being calculated from 4-vector
        //Photos::forceMassFromEventRecord(11);

        // Exclude particles with given status code from being processed
        // or taken into account during momentum conservation calculation
        //Photos::ignoreParticlesWithStatus(3);

        // Remove status code from the list of ignored status codes
        //Photos::DeIgnoreParticlesWithStatus(3);

        // Force writing history decay products for vertices
        // modified i.e. with added photons. These particles will
        // have the provided status code. Photos will ignore
        // all particles with this status code.
        //Photos::createHistoryEntries(true,3);

        MC_Initialize();

        // Begin event loop
        for (unsigned long iEvent = 0; iEvent < NumberOfEvents; ++iEvent)
        {
                if(iEvent%1000==0) Log::Info()<<"Event: "<<iEvent<<"\t("<<iEvent*(100./NumberOfEvents)<<"%)"<<endl;
                if(!pythia.next()) continue;

                // Convert event record to HepMC
                HepMC::GenEvent * HepMCEvt = new HepMC::GenEvent();
                ToHepMC.fill_next_event(event, HepMCEvt);

                // Run TAUOLA on the event
                TauolaHepMCEvent * t_event = new TauolaHepMCEvent(HepMCEvt);

                // We may want to undecay previously decayed taus.
                //t_event->undecayTaus();
                t_event->decayTaus();
                delete t_event;

                //Log::LogPhlupa(2,4);

                if(iEvent<EventsToCheck)
                {
                        cout<<"                                          "<<endl;
                        cout<<"Momentum conservation chceck BEFORE/AFTER Photos"<<endl;
                        checkMomentumConservationInEvent(HepMCEvt);
                }

                // Run PHOTOS on the event
                PhotosHepMCEvent evt(HepMCEvt);
                evt.process();

                if(iEvent<EventsToCheck)
                {
                        checkMomentumConservationInEvent(HepMCEvt);
                }

                // Run MC-TESTER on the event
                HepMCEvent temp_event(*HepMCEvt,false);
                MC_Analyze(&temp_event);

                // Print out last 5 events
                if(iEvent>=NumberOfEvents-5)
                {
                        Log::RedirectOutput(Log::Info());
                        //pythia.event.list();
                        HepMCEvt->print();
                        Log::RevertOutput();
                }

                // Clean up
                delete HepMCEvt;
        }

        Log::RedirectOutput(Log::Info());
        pythia.statistics();
        Log::RevertOutput();

        MC_Finalize();
}