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