photos_pythia_example.cxx File Reference

#include "Pythia.h"
#include "HepMCInterface.h"
#include "Generate.h"
#include "HepMCEvent.H"
#include "Setup.H"
#include "Photos/Photos.h"
#include "Photos/PhotosHepMCEvent.h"
#include "Photos/Log.h"

Go to the source code of this file.

Functions
void	checkMomentumConservationInEvent (HepMC::GenEvent *evt)
void	switch_history_entries_status (HepMC::GenEvent *evt)
int	main (int argc, char **argv)

Variables
bool	ShowersOn =true
unsigned long	NumberOfEvents = 10000
unsigned int	EventsToCheck =20

Function Documentation

void checkMomentumConservationInEvent

(

HepMC::GenEvent *

evt

)

Definition at line 35 of file photos_pythia_example.cxx.

References m, p, and pz.

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

Definition at line 145 of file photos_pythia_example.cxx.

References checkMomentumConservationInEvent(), EventsToCheck, evt, chigen::initialize(), NumberOfEvents, Photospp::PhotosEvent::process(), and chigen::pythia::pythia.

{
        // Initialization of pythia
        HepMC::I_Pythia8 ToHepMC;
        Pythia pythia;
        Event& event = pythia.event;
        //pythia.settings.listAll();

        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.init( 11, -11, 91.187);                           //e+ e- collisions

        MC_Initialize();

        Photos::initialize();
        //Photos::setDoubleBrem(false);
        //Photos::setExponentiation(false);

        Photos::setInfraredCutOff(0.01/91.187); // 10MeV for scale to M_Z=91.187
        Photos::maxWtInterference(3.0);
  //Photos::createHistoryEntries(true,3);

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

        // 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);
                //HepMCEvt->print();

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

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

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

    // Uncomment to turn on switching of the status code of history entries
    // with the regular entries for stable particles
    //switch_history_entries_status(HepMCEvt);

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

                //HepMCEvt->print();

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

                // Print out last 5 events
                if(iEvent>=NumberOfEvents-5) HepMCEvt->print();

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

void switch_history_entries_status

(

HepMC::GenEvent *

evt

)

Definition at line 75 of file photos_pythia_example.cxx.

References exit(), p, p2, and v.

{
  for ( HepMC::GenEvent::particle_const_iterator p = evt->particles_begin();
              p != evt->particles_end(); ++p )
        {
    if((*p)->status()==3)
    {
      if((*p)->pdg_id()==22) continue;

      int barcode = (*p)->barcode();

      HepMC::GenVertex *v = (*p)->production_vertex();
      
      // History entries are added after photons, so we check what is the
      // position of current particle relative to photons.
      int position = 0;
      int last_photon_position = -1;
      
      for(HepMC::GenVertex::particles_out_const_iterator p2 = v->particles_out_const_begin();
          p2 != v->particles_out_const_end(); ++p2)
      {
        position++;

        if((*p2)->barcode()==barcode) break;

        if((*p2)->pdg_id()==22) { last_photon_position=position; }
      }
      
      // If particle is found prior to photons - it was already processed, so skip it
      if(last_photon_position<0) continue;
      
      position -= last_photon_position;
      HepMC::GenParticle *part = NULL;
      
      // Now, find the particle that corresponds to this history entry
      for(HepMC::GenVertex::particles_out_const_iterator p2 = v->particles_out_const_begin();
          p2 != v->particles_out_const_end(); ++p2)
      {
        --position;
        
        if     (position >  0) continue;
        else if(position == 0) part = *p2;
        else
        {
          // Give all remaining photons status 3
          if((*p2)->pdg_id()==22 ) (*p2)->set_status(3);

          // Finish if there are no more photons
          else break;
        }
      }

      // Check if this is the particle we are looking for
      if( part->pdg_id() != (*p)->pdg_id())
      {
        cout<<"switch_history_entries_status: mismatch in pdg_id of history entry"<<endl;
        cout<<"and its corresponding particle. The algorithm does not work correctly."<<endl;
        exit(-1);
      }

      // Skip this particle if its status is not 1
      if(part->status()!=1) continue;
      
      // Switch status codes of these particles
      part->set_status(3);
      (*p)->set_status(1);
    }
  }
}

Variable Documentation

unsigned int EventsToCheck =20

Definition at line 29 of file photos_pythia_example.cxx.

unsigned long NumberOfEvents = 10000

Definition at line 28 of file photos_pythia_example.cxx.

bool ShowersOn =true

Definition at line 27 of file photos_pythia_example.cxx.

Referenced by main().