photos_tauola_test.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"

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

//PHOTOS header files
#include "Photos/Photos.h"
#include "Photos/PhotosHepMCParticle.h"
#include "Photos/PhotosHepMCEvent.h"
#include "Photos/Log.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)
{

        // Program needs at least 4 parameters
        if(argc<5)
        {
                cout<<endl<<"Usage: "<<argv[0]<<" <pythia_conf> <pythia_mode> <no_events> <tauola_mode> [ <alpha_order> <ScalarNLO_mode> ]"<<endl;
                cout<<endl<<"   eg. "<<argv[0]<<" pythia_H.conf 0 10000 4 0 0"<<endl;
                cout<<endl;
                return -1;
        }

        HepMC::I_Pythia8 ToHepMC;

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

        pythia.readString("HadronLevel:Hadronize = off");
        pythia.readString("SpaceShower:QEDshower = off");
        pythia.readString("SpaceShower:QEDshowerByL = off");
        pythia.readString("SpaceShower:QEDshowerByQ = off");
        pythia.readString("PartonLevel:ISR = off");
        pythia.readString("PartonLevel:FSR = off");

        // Tauola is currently set to undecay taus. Otherwise, uncomment this line.
        //pythia.particleData.readString("15:mayDecay = off");

        /********************************************************
          Read input parameters from console. List of parameters:
          1. Pythia configuration filename
          2. Pythia mode - e+e-@200GeV , e+e-@91.187GeV, pp@14TeV or e+e-@500GeV
             (only e+e-@91.187GeV and e+e-@500GeV are used in this example)
          3. Number of events
          4. Tauola decay mode (refer to Tauola documentation)
          5. Photos - use 1-photon mode on/off
          6. Photos - use ScalarNLO mode on/off

          Example where all input parameters are used:

          ./photos_tauola_test.exe pythia_H.conf 1 100000 4 0 0
          - use pythia_H.conf
          - initialize using e+ e- @ 91.187GeV collisions
          - generate 100 000 events
          - fix TAUOLA decay to channel 4 (RHORHO_MODE)
          - Photos is not using 1-photon mode (default option)
          - Photos is not in ScalarNLO mode (default option)
        *********************************************************/

        // 1. Load pythia configuration file (argv[1], from console)
        if(argc>1) pythia.readFile(argv[1]);

        // 2. Initialize pythia to e+e-@91.17GeV or e+e-@500GeV collisions (argv[2], from console)
        if(atoi(argv[2])==1)      pythia.init( 11, -11, 91.187); // e+ e- collisions
        else if(atoi(argv[2])==3) pythia.init( 11, -11, 500.);   // e+ e- collisions
        else
        {
                cout<<"ERROR: Wrong Pythia mode ("<<atoi(argv[4])<<")"<<endl;
                cout<<"       Only modes '1' and '3' are used by this program."<<endl;
                return -1;
        }

        // 3. Get number of events (argv[3], from console)
        if(argc>3) NumberOfEvents=atoi(argv[3]);

        // 4. Set Tauola decay mode (argv[4], from console)
        if(argc>4)
        {
                // argv[4]=3 (tau => pi nu_tau)    for Ztautau
                // argv[4]=4 (tau => pi pi nu_tau) for Htautau
                Tauola::setSameParticleDecayMode(atoi(argv[4]));
                Tauola::setOppositeParticleDecayMode(atoi(argv[4]));
        }

        Tauola::initialize();
        Photos::initialize();

        Photos::setExponentiation(true);
        Photos::setInfraredCutOff(1.e-6);
        Photos::maxWtInterference(3.0);

        // 5. Check if we're using 1-photon mode
        if( argc>5 && atoi(argv[5]) )
        {
                Photos::setDoubleBrem(false);
                Photos::setExponentiation(false);

                // Set infrared cutoff to 10MeV for scale M_Z=91.187GeV or 500 GeV
                if(atoi(argv[2])==1) Photos::setInfraredCutOff(0.01/91.187);
                else                 Photos::setInfraredCutOff(0.01/500.);

                Photos::maxWtInterference(2.0);
        }

        // 6. Check if we're in ScalarNLO mode
        if( argc>6 )
        {
                Tauola::setEtaK0sPi(1,1,0);
    
                // Check if we are using NLO
                if(atoi(argv[6])) Photos::setMeCorrectionWtForScalar(true);
        }

        Log::SummaryAtExit();
        cout.setf(ios::fixed);

        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;

                HepMC::GenEvent * HepMCEvt = new HepMC::GenEvent();
                ToHepMC.fill_next_event(event, HepMCEvt);

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

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

                // Since we let Pythia decay taus, we have to undecay them first.
                t_event->undecayTaus();
                t_event->decayTaus();
                delete t_event;

                // 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);

                //clean up
                delete HepMCEvt;
        }
        pythia.statistics();
        MC_Finalize();
}