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/Configuration.h"
#include "../../include/JetFinder.h"
#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/EtMissAnalysis.h"
#include "../../include/UserAnalysis.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);
  
  JetFinder FindJet;
  FindJet.Init(currentconfig->max_eta,currentconfig->min_pt);
  FindJet.InitJetFinder(currentconfig->coneRadius,currentconfig->overlapThreshold,currentconfig->jet_ptmin,currentconfig->lepton_ptmin,currentconfig->DeltaR_lepton_track);

  std::vector<baseAnalysis*> analysis;

  // create analysis Objects and initialize them                                                                                          
  if (currentconfig->dijet) {
    std::cout << "Adding module DiJetAnalysis" << std::endl;
    analysis.push_back(new DiJetAnalysis);
  }
  if (currentconfig->wplusjet) {
    std::cout << "Adding module WplusJetAnalysis" << std::endl;
    analysis.push_back(new WplusJetAnalysis);
  }
  if(currentconfig->z){
    std::cout << "Adding module ZAnalysis" << std::endl;
    analysis.push_back(new ZAnalysis);
  }
  if(currentconfig->tau){
    std::cout << "Adding module TauAnalysis" << std::endl;
    analysis.push_back(new TauAnalysis);
  }
  if(currentconfig->top){
    std::cout << "Adding module TopAnalysis" << std::endl;
    analysis.push_back(new TopAnalysis);
  }
  if(currentconfig->ue){
    std::cout << "Adding module UEAnalysis" << std::endl;
    analysis.push_back(new UEAnalysis);
  }  
  if(currentconfig->etmiss){
    std::cout << "Adding module EtMissAnalysis" << std::endl;
    analysis.push_back(new EtMissAnalysis);
  }
  if(currentconfig->user){
    std::cout << "Adding module UserAnalysis" << std::endl;
    analysis.push_back(new UserAnalysis);
  }

  // init requested analysis modules and jet finder
  for (std::vector<baseAnalysis*>::const_iterator i(analysis.begin()); i != analysis.end(); ++i) {
    (*i)->Init(currentconfig->max_eta,currentconfig->min_pt);
    (*i)->InitJetFinder(currentconfig->coneRadius,currentconfig->overlapThreshold,currentconfig->jet_ptmin,currentconfig->lepton_ptmin,currentconfig->DeltaR_lepton_track);
  }

  
  HepMC::I_Pythia8 ToHepMC;

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

  char processline[64];
  // 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());
      }

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

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

    //call the Jetfinder
    vector<fastjet::PseudoJet> inclusive_jets = FindJet.GetJet(hepmcevt);

    //call analysis
    int ret_all = true;
    for (std::vector<baseAnalysis*>::const_iterator i(analysis.begin()); i != analysis.end(); ++i) {
      (*i)->SetJet(&inclusive_jets);
      int ret = (*i)->Process(hepmcevt);
      if (ret == false) ret_all = false;
    }
    
    if (ret_all==false) {
      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
    FindJet.ClearEvent(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);

  TFile* f = new TFile(currentconfig->OutputFileName.c_str(),"RECREATE");
  for (std::vector<baseAnalysis*>::const_iterator i(analysis.begin()); i != analysis.end(); ++i) {
    (*i)->Finalize(f);
  }

  // clean up analysis modules
  for (std::vector<baseAnalysis*>::const_iterator i(analysis.begin()); i != analysis.end(); ++i) {
    delete (*i);
  }
  analysis.resize(0);

  // clean up config object
  delete currentconfig;

  // Tidy up
  std::cout<<"pythia8: program ends"<<std::endl;
  f->Close();

  std::cout<<"Run successfully!"<<std::endl;
  return 0;
}

---