PhotonConversion/pythia/photonJetGenerator.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("promptphoton:all = on");
0018   pythiaengine.readString("HardQCD:all = on");
0019   pythiaengine.readString("PhaseSpace:pTHatMin = 10.");
0020   pythiaengine.readString("Random::setSeed = on");
0021   pythiaengine.readString("Random::seed =0");
0022   //pythiaengine.readString("111:onMode = off"); ///pi0 won't decay
0023   pythiaengine.init();
0024
0025   string tfilename = filename+"_analysis.root";
0026   TFile *outFile = new TFile(tfilename.c_str(),"RECREATE");
0027   TTree *photonTree = new TTree("photonTree","phat phirn tree");
0028   photonTree->SetAutoSave(300);
0029   vector<float> photon_pT;
0030   photonTree->Branch("photon_pT",&photon_pT);
0031   TTree *nophotonTree = new TTree("nophotonTree","phat phirn tree");
0032   unsigned long noPhotonEvents=0;
0033   nophotonTree->Branch("n",&noPhotonEvents);
0034
0035   for (int iEvent = 0; iEvent < nEvents; ++iEvent)
0036   {
0037     if (!pythiaengine.next()){
0038       cout<<"pythia.next() failed"<<"\n";
0039       iEvent--;
0040       continue;
0041     }
0042     photon_pT.clear();
0043     for(unsigned ipart=0; ipart!=pythiaengine.event.size(); ipart++){
0044       if(pythiaengine.event[ipart].id()==22&&pythiaengine.event[ipart].isFinal()&&pythiaengine.event[ipart].pT()>5
0045           &&TMath::Abs(pythiaengine.event[ipart].eta()))photon_pT.push_back(pythiaengine.event[ipart].pT());
0046     }
0047     if (photon_pT.size()>0)photonTree->Fill();
0048     else{
0049       noPhotonEvents++;
0050     }
0051     if(!signalOnly||photon_pT.size()>0){
0052       HepMC::GenEvent* hepmcevtfrag = new HepMC::GenEvent(); //create HepMC event
0053       ToHepMC.fill_next_event( pythiaengine, hepmcevtfrag ); //convert event from pythia to HepMC
0054       ascii_io << hepmcevtfrag;//write event to file
0055       delete hepmcevtfrag; //delete event so it can be redeclared next time
0056     }
0057   }
0058   nophotonTree->Fill();
0059   outFile->Write();
0060   outFile->Close();
0061   pythiaengine.stat();
0062 }
0063
0064 int main(int argc, char const *argv[] )
0065 {
0066   string fileOut = string(argv[1]);
0067   bool signalOnly=false;
0068   if(argv[2][0]=='1'||argv[2][0]=='t')signalOnly=true;
0069   long nEvents =strtol(argv[2],NULL,10);  // 5000000;
0070   generator(fileOut,nEvents,signalOnly);
0071   cout<<"All done"<<endl;
0072   return 0;
0073 }