PhotonConversion/pythia/generator.cc

0001 #include "Pythia8/Pythia.h"
0002 #include "Pythia8Plugins/HepMC2.h" //added plugin for HepMC, think we will need some new library in pythia for this
0003 #include <TTree.h>
0004 #include <TFile.h>
0005 using namespace Pythia8;
0006 using namespace std;
0007
0008 void generator(std::string filename, long nEvents, bool signalOnly=false){
0009 /*  using namespace HepMC;
0010   string hepName = filename+".dat";    //filenames
0011   HepMC::Pythia8ToHepMC ToHepMC;    // Interface for conversion from Pythia8::Event to HepMC event.
0012   HepMC::IO_GenEvent ascii_io(hepName, std::ios::out);*/ //file where HepMC events will be stored.
0013
0014   /*pythia set up*/
0015   Pythia pythiaengine;
0016   pythiaengine.readString("Beams:eCM = 200.");
0017   pythiaengine.readString("SoftQCD:nonDiffractive = on");
0018   pythiaengine.readString("SoftQCD:singleDiffractive = on");
0019   pythiaengine.readString("SoftQCD:doubleDiffractive = on");
0020   pythiaengine.readString("PhaseSpace:pTHatMin = 0");
0021   pythiaengine.readString("Random::setSeed = on");
0022   pythiaengine.readString("Random::seed =0");
0023   //pythiaengine.readString("111:onMode = off"); ///pi0 won't decay
0024   pythiaengine.init();
0025
0026   string tfilename = filename+"_analysis.root";
0027   TFile *outFile = new TFile(tfilename.c_str(),"RECREATE");
0028   TTree *photonTree = new TTree("photonTree","phat phirn tree");
0029   photonTree->SetAutoSave(300);
0030   vector<float> photon_pT;
0031   photonTree->Branch("photon_pT",&photon_pT);
0032   TTree *nophotonTree = new TTree("nophotonTree","phat phirn tree");
0033   unsigned noPhotonEvents=0;
0034   nophotonTree->Branch("nNoPhoton",&noPhotonEvents);
0035
0036   for (int iEvent = 0; iEvent < nEvents; ++iEvent)
0037   {
0038     if (!pythiaengine.next()){
0039       cout<<"pythia.next() failed"<<"\n";
0040       iEvent--;
0041       continue;
0042     }
0043     photon_pT.clear();
0044     for(unsigned ipart=0; ipart!=pythiaengine.event.size(); ipart++){
0045       if(pythiaengine.event[ipart].id()==22&&pythiaengine.event[ipart].isFinal()&&pythiaengine.event[ipart].pT()>5
0046           &&TMath::Abs(pythiaengine.event[ipart].eta()))photon_pT.push_back(pythiaengine.event[ipart].pT());
0047     }
0048     if (photon_pT.size()>0)photonTree->Fill();
0049     else{
0050       noPhotonEvents++;
0051     }
0052 /*    if(!signalOnly||photon_pT.size()>0){
0053       HepMC::GenEvent* hepmcevtfrag = new HepMC::GenEvent(); //create HepMC event
0054       ToHepMC.fill_next_event( pythiaengine, hepmcevtfrag ); //convert event from pythia to HepMC
0055       ascii_io << hepmcevtfrag;//write event to file
0056       delete hepmcevtfrag; //delete event so it can be redeclared next time
0057     }*/
0058   }
0059   nophotonTree->Fill();
0060   outFile->Write();
0061   outFile->Close();
0062   pythiaengine.stat();
0063 }
0064
0065 int main(int argc, char const *argv[] )
0066 {
0067   string fileOut = string(argv[1]);
0068   bool signalOnly=false;
0069   if(argv[2][0]=='1'||argv[2][0]=='t')signalOnly=true;
0070   long nEvents =strtol(argv[2],NULL,10);  // 5000000;
0071   generator(fileOut,nEvents,signalOnly);
0072   cout<<"All done"<<endl;
0073   return 0;
0074 }