photosLCG_pythia_example.cxx

/**
 * Example of use of photos C++ interface. Pythia events are
 * generated first and photos used for FSR.
 *
 * @author LCG & Tomasz Przedzinski
 * @date 6 May 2011
 */

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

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

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

bool ShowersOn=true;
unsigned long NumberOfEvents = 100000;

// Calculates energy ratio between (l + bar_l) and (l + bar_l + X)
double calculate_ratio(HepMC::GenEvent *evt, double *ratio_2)
{
        double ratio = 0.0;
        for(HepMC::GenEvent::particle_const_iterator p=evt->particles_begin();p!=evt->particles_end(); p++)
        {
                // Search for Z
                if( (*p)->pdg_id() != 23 ) continue;

                // Ignore this Z if it does not have at least two daughters
                if( !(*p)->end_vertex() || (*p)->end_vertex()->particles_out_size() < 2 ) continue;

                // Sum all daughters other than photons
                double sum = 0.0;
                for(HepMC::GenVertex::particle_iterator pp = (*p)->end_vertex()->particles_begin(HepMC::children);
                    pp != (*p)->end_vertex()->particles_end(HepMC::children);
                    ++pp)
                {
                   // (*pp)->print();
                   if( (*pp)->pdg_id() != 22 ) sum += (*pp)->momentum().e();
                }

                // Calculate ratio and ratio^2
                ratio    = sum     /   (*p)->momentum().e();
                *ratio_2 = sum*sum / ( (*p)->momentum().e() * (*p)->momentum().e() );

                // Assuming only one Z decay per event
                return ratio;
        }
        return 0.0;
}

int main(int argc,char **argv)
{
        // Initialization of pythia
        HepMC::I_Pythia8 ToHepMC;
        Pythia pythia;
        Event& event = pythia.event;
        //pythia.settings.listAll();

        if(!ShowersOn)
        {
                //pythia.readString("HadronLevel:all = off");
                pythia.readString("HadronLevel:Hadronize = off");
                pythia.readString("SpaceShower:QEDshower = off");
                pythia.readString("SpaceShower:QEDshowerByL = off");
                pythia.readString("SpaceShower:QEDshowerByQ = off");
        }
        pythia.readString("PartonLevel:ISR = on");
        pythia.readString("PartonLevel:FSR = off");

        pythia.readString("WeakSingleBoson:ffbar2gmZ = on");
        pythia.readString("23:onMode = off");
        pythia.readString("23:onIfAny = 13");
        //pythia.readString("23:onIfAny = 11");
        pythia.init( 11, -11, 91.187);                           //e+ e- collisions

        Photos::initialize();

        // Turn on NLO corrections - only for PHOTOS 3.2 or higher

        //Photos::setMeCorrectionWtForZ(zNLO);
        //Photos::maxWtInterference(4.0);
        //Photos::iniInfo();

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

        double ratio_exp   = 0.0, ratio_alpha   = 0.0;
        double ratio_exp_2 = 0.0, ratio_alpha_2 = 0.0;
        double buf = 0.0;

        Photos::setDoubleBrem(true);
        Photos::setExponentiation(true);
        Photos::setInfraredCutOff(0.000001);

        Log::Info() << "---------------- First run - EXP ----------------" <<endl;

        // Begin event loop
        for(unsigned long iEvent = 0; iEvent < NumberOfEvents; ++iEvent)
        {
                if(iEvent%10000==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);
                //HepMCEvt->print();

                //Log::LogPhlupa(1,3);

                // Call photos
                PhotosHepMCEvent evt(HepMCEvt);
                evt.process();

                ratio_exp   += calculate_ratio(HepMCEvt,&buf);
                ratio_exp_2 += buf;

                //HepMCEvt->print();

                // Clean up
                delete HepMCEvt;
        }

        Photos::setDoubleBrem(false);
        Photos::setExponentiation(false);
        Photos::setInfraredCutOff(1./91.187);  // that is 1 GeV for
                                               // pythia.init( 11, -11, 91.187);

        Log::Info() << "---------------- Second run - ALPHA ORDER ----------------" <<endl;
        
        // Begin event loop
        for(unsigned long iEvent = 0; iEvent < NumberOfEvents; ++iEvent)
        {
                if(iEvent%10000==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);
                //HepMCEvt->print();

                //Log::LogPhlupa(1,3);

                // Call photos
                PhotosHepMCEvent evt(HepMCEvt);
                evt.process();

                ratio_alpha   += calculate_ratio(HepMCEvt,&buf);
                ratio_alpha_2 += buf;

                //HepMCEvt->print();

                // Clean up
                delete HepMCEvt;
        }
        
        pythia.statistics();

        ratio_exp   = ratio_exp   / (double)NumberOfEvents;
        ratio_exp_2 = ratio_exp_2 / (double)NumberOfEvents;

        ratio_alpha   = ratio_alpha   / (double)NumberOfEvents;
        ratio_alpha_2 = ratio_alpha_2 / (double)NumberOfEvents;

        double err_exp   = sqrt( (ratio_exp_2   - ratio_exp   * ratio_exp  ) / (double)NumberOfEvents );
        double err_alpha = sqrt( (ratio_alpha_2 - ratio_alpha * ratio_alpha) / (double)NumberOfEvents );
        
        cout.precision(6);
        cout.setf(ios::fixed);
        cout << "********************************************************" << endl;
        cout << "* Z -> l + bar_l + gammas                              *" << endl;
        cout << "* E(l + bar_l) / E(l + bar_l + gammas) ratio           *" << endl;
        cout << "*                                                      *" << endl;
        cout << "* PHOTOS - EXP:          " << ratio_exp   <<" +/- "<<err_exp  <<"      *" << endl;
        cout << "* PHOTOS - ALPHA ORDER:  " << ratio_alpha <<" +/- "<<err_alpha<<"      *" << endl;
        cout << "********************************************************" << endl;

}