pythia8.cc

Go to the documentation of this file.

/***********************************************************************
*                                                                      *
* File: pythia8_jet.cc                                                 *
* This is a simple test program.                                       *
* Copyright \ufffd 2005 Torbj\ufffdrn Sj\ufffdstrand                   *
*                                                                      *
* Modified by sj and ca                                                *
*                                                                      *
***********************************************************************/
#include "Pythia.h"

#include "LHAFortran.h"

#include "HepMCInterface.h"

#include "HepMC/GenEvent.h"

// this is needed to generate some random numbers                                                                                                                                                                    
#include <stdlib.h>
#include <time.h>

#include <iostream>
#include <stdio.h>
#include "HepMC/IO_GenEvent.h"
#include "HepMC/GenEvent.h"
#include "HepMC/IO_AsciiParticles.h"
#include "HepMC/SimpleVector.h"

#include "TH1.h"
#include "TFile.h"
#include "TMath.h"
#include <TCanvas.h>

#include "TF1.h"
#include <TH2.h>
#include <TStyle.h>
#include "TLegend.h"
#include <TTree.h>
#include <TString.h>

using namespace Pythia8; 
using namespace std;

#include "../../include/DiJetAnalysis.h"
#include "../../include/WplusJetAnalysis.h"
#include "../../include/ZAnalysis.h"
#include "../../include/TauAnalysis.h"
#include "../../include/TopAnalysis.h"
#include "../../include/UEAnalysis.h"
#include "../../include/Configuration.h"

/**************************************************************************/

int main(int argc, const char* argv[]) {

    if(argc <= 1 ) {
        std::cout << "Usage: " << argv[0] << " <Configuration File>" << std::endl;
        return 0;
    }


    Configuration* currentconfig = new Configuration();
    if(currentconfig->read(argv[1]))
        {
            std::cout << " ... skip !" << std::endl;
            return 0;
        }

    //argv[2] = (Pythia8DATA);

    DiJetAnalysis dijet;
    WplusJetAnalysis wplusjet;
    ZAnalysis z0;
    TauAnalysis tau;
    TopAnalysis top;
    UEAnalysis ue;

    // create analysis Objects and initialize them                                                                                          
  if(currentconfig->dijet){
    dijet.Init(currentconfig->max_eta,currentconfig->min_pt);
    dijet.InitJetFinder(currentconfig->coneRadius,currentconfig->overlapThreshold,currentconfig->jet_ptmin);
  }

  if(currentconfig->wplusjet){
    wplusjet.Init(currentconfig->max_eta,currentconfig->min_pt);
    wplusjet.InitJetFinder(currentconfig->coneRadius,currentconfig->overlapThreshold,currentconfig->jet_ptmin);
  }

  if(currentconfig->z){
    z0.Init(currentconfig->max_eta,currentconfig->min_pt);
    z0.InitJetFinder(currentconfig->coneRadius,currentconfig->overlapThreshold,currentconfig->jet_ptmin);
  }
  if(currentconfig->tau){
    tau.Init(currentconfig->max_eta,currentconfig->min_pt);
    tau.InitJetFinder(currentconfig->coneRadius,currentconfig->overlapThreshold,currentconfig->jet_ptmin);
  }
  if(currentconfig->top){
    top.Init(currentconfig->max_eta,currentconfig->min_pt);
    top.InitJetFinder(currentconfig->coneRadius,currentconfig->overlapThreshold,currentconfig->jet_ptmin);
  }
  if(currentconfig->ue){
    ue.Init(currentconfig->max_eta,currentconfig->min_pt);
    ue.InitJetFinder(currentconfig->coneRadius,currentconfig->overlapThreshold,currentconfig->jet_ptmin);
  }
  
  HepMC::I_Pythia8 ToHepMC;

  // Generator. Process selection. LHC initialization.
  //Pythia pythia;
  Pythia pythia(std::getenv("Pythia8DATA"));

  // as done for the random seed, if it is specified
  char processline[64];
  pythia.readString("Random:setSeed = on");
  sprintf(processline,"Random:seed = %ld",currentconfig->rseed);
  pythia.readString(processline);

  // Read in commands from external file.
  for (unsigned int files=0; files<currentconfig->ConfigFileNames.size(); files++)
      {
          std::cout << "pythia8: Will use settings from file " << currentconfig->ConfigFileNames[files].c_str() << std::endl;
          pythia.readFile(currentconfig->ConfigFileNames[files].c_str());
      }

  // master switches for FSR, ISR and MI
  sprintf(processline,"PartonLevel:FSR = %d",int(currentconfig->FSR));
  pythia.readString(processline);

  sprintf(processline,"PartonLevel:ISR = %d",int(currentconfig->ISR));
  pythia.readString(processline);

  sprintf(processline,"PartonLevel:MI = %d",int(currentconfig->MI));
  pythia.readString(processline);

  // Settings can also be changed with this command
  //pythia.readString("WeakBosonAndParton:all = on");

  pythia.init( 2212, 2212, 14000.);

  // Shorthand for the event and for settings.
  Event& event = pythia.event;
  Settings& settings = pythia.settings;

  // Extract settings to be used in the main program.
  int nEvent = currentconfig->nevents;
  //sj//  int nEvent = settings.mode("Main:numberOfEvents");
  int nList = settings.mode("Main:numberToList");
  int nShow = settings.mode("Main:timesToShow");
  int nAbort = settings.mode("Main:timesAllowErrors");
  bool showChangedSettings = settings.flag("Main:showChangedSettings");
  bool showAllSettings = settings.flag("Main:showAllSettings");
  bool showAllParticleData = settings.flag("Main:showAllParticleData");

  // List changed data.
  if (showChangedSettings) settings.listChanged();
  if (showAllSettings) settings.listAll();

  // List particle data.  
  if (showAllParticleData) ParticleDataTable::listAll();

  // Begin event loop.
  int nShowPace = max(1,nEvent/nShow); 
  int iAbort = 0; 
  for (int iEvent = 0; iEvent < nEvent; ++iEvent) {

        if (iEvent%nShowPace == 0) 
            cout << "pythia8: Now begin event "  << iEvent << endl;

    // Generate events. Quit if too many failures.
    if (!pythia.next()) {
      if (++iAbort < nAbort) continue;
      cout << "pythia8:  Event generation aborted prematurely, owing to error!\n"; 
      break;
    }
 
    // List first few events, both hard process and complete events.
    //if (iEvent < nList) { 
    //if (iEvent < nEvent) { 
    //  pythia.process.list();
      //event.list();
    // }

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

    //hepmcevt->print();

    int ret_dijet = SUCCESS;
    int ret_wplusjet = SUCCESS;
    int ret_z0 = SUCCESS;
    int ret_tau = SUCCESS;
    int ret_top = SUCCESS;
    int ret_ue = SUCCESS;

    //call analysis
    if(currentconfig->dijet)       ret_dijet = dijet.Process(hepmcevt);
    if(currentconfig->wplusjet)    ret_wplusjet = wplusjet.Process(hepmcevt);
    if(currentconfig->z)           ret_z0 = z0.Process(hepmcevt);
    if(currentconfig->tau)  ret_tau = tau.Process(hepmcevt);
    if(currentconfig->top)  ret_top = top.Process(hepmcevt);
    if(currentconfig->ue)   ret_ue = ue.Process(hepmcevt);

    if (ret_dijet==FAILURE || ret_wplusjet==FAILURE || ret_z0==FAILURE 
        || ret_tau==FAILURE || ret_top==FAILURE || ret_ue==FAILURE)
      {
        std::cout << "Processing event number " << iEvent << std::endl;    
        std::cout << "throw FAILURE"<<std::endl;
        //std::cout<<std::endl;
        //std::cout<<"HepMC Output"<<std::endl;
        //hepmcevt->print();
        std::cout<<std::endl;
      }

    //delete all the jet objects to be free for the next event
    if(currentconfig->dijet)       dijet.DeleteJetObject(hepmcevt);
    if(currentconfig->wplusjet)    wplusjet.DeleteJetObject(hepmcevt);
    if(currentconfig->z)         z0.DeleteJetObject(hepmcevt);
    if(currentconfig->tau)       tau.DeleteJetObject(hepmcevt);
    if(currentconfig->top)       top.DeleteJetObject(hepmcevt);
    if(currentconfig->ue)        ue.DeleteJetObject(hepmcevt);

    // clear memory from this event
    delete hepmcevt;

    // End of event loop.
  }

  // Final statistics.
  pythia.statistics();
  //  pythia.particleData.listAll();
  pythia.settings.writeFile("allSettings.cmnd",true);

  //JetpAnalysis jet;
  TFile* f = new TFile(currentconfig->OutputFileName.c_str(),"RECREATE");
    if(currentconfig->dijet)       dijet.Finalize(f);
    if(currentconfig->wplusjet)    wplusjet.Finalize(f);
    if(currentconfig->z)         z0.Finalize(f);
    if(currentconfig->tau)       tau.Finalize(f);
    if(currentconfig->top)       top.Finalize(f);
    if(currentconfig->ue)       ue.finalize(f);

  std::cout<<"pythia8: program ends"<<std::endl;
  f->Close();
  // Done.
  return 0;
}

---