// main01.cc is a part of the PYTHIA event generator. // Copyright (C) 2008 Torbjorn Sjostrand. // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details. // Please respect the MCnet Guidelines, see GUIDELINES for details. // This is a simple test program. It fits on one slide in a talk. // It studies the charged multiplicity distribution at the LHC. #include #include #include #include #include #include #include #include #include #include #include "Pythia.h" #include "TTree.h" #include "TFile.h" #include "TRandom1.h" using namespace Pythia8; using namespace std; #ifdef _SINGLE_PRECISION_ # define _FLOAT_T_ float #else # define _FLOAT_T_ double #endif #ifndef _AUTHOR_ #define _AUTHOR_ "Peter Loch" #endif #ifndef _EMAIL_ #define _EMAIL_ "" #endif //---------------------------------------------------------------------- /// a function that pretty prints a list of jets int main(int argc, char* argv[]) { // timing tms* startTime = new tms(); clock_t sT = times(startTime); // module name string moduleName("MinBiasProcessor"); // messages #ifdef _SINGLE_PRECISION_ printf("\n\n[%s] - INFO - Single precision version compiled on %s at %s <%s>\n", moduleName.c_str(),__DATE__,__TIME__,__FILE__); #else printf("\n\n[%s] - INFO - Double precision version compiled on %s at %s <%s>\n", moduleName.c_str(),__DATE__,__TIME__,__FILE__); #endif printf("\n[%s] - INFO - Author %s %s\n\n",moduleName.c_str(),_AUTHOR_,_EMAIL_); printf("[%s] - INFO - BeginJob (user time/system time) (%10i,%10i) - return (%10i)\n", moduleName.c_str(),startTime->tms_utime,startTime->tms_stime,sT); string fName("minbias_pythia8_"); int mPart(100000); TFile* rFile((TFile*)0); TTree* rTree((TTree*)0); // ROOT tuple structure int Nentry; int Npartons; int Nparticles; int ID[mPart]; int Stat[mPart]; _FLOAT_T_ Charge[mPart]; _FLOAT_T_ Px[mPart]; _FLOAT_T_ Py[mPart]; _FLOAT_T_ Pz[mPart]; _FLOAT_T_ P0[mPart]; _FLOAT_T_ Pm[mPart]; _FLOAT_T_ Pt[mPart]; _FLOAT_T_ Rap[mPart]; _FLOAT_T_ Phi[mPart]; _FLOAT_T_ Eta[mPart]; // processing directives int nEvts(1000); int nFile(1); _FLOAT_T_ cMass(13000.); // scan argument list if ( argc > 1 ) { // catch all input map dataInput; int i(1); while ( i::iterator fMap(dataInput.begin()); map::iterator lMap(dataInput.end()); for ( ; fMap != lMap; ++fMap ) { if ( (*fMap).first == "-ne" ) { int dEvts(nEvts); nEvts = atoi((*fMap).second.c_str()); printf("[%s] - INFO - number of events/file updated %i (default %i)\n", moduleName.c_str(),nEvts,dEvts); } else if ( (*fMap).first == "-nf" ) { int dFile(nFile); nFile = atoi((*fMap).second.c_str()); printf("[%s] - INFO - number of files updated %i (default %i)\n", moduleName.c_str(),nFile,dFile); } else if ( (*fMap).first == "-s" ) { _FLOAT_T_ dMass(cMass); cMass = atof((*fMap).second.c_str()); printf("[%s] - INFO - center-of-mass energy updated %7.1f GeV (default %7.1f GeV)\n", moduleName.c_str(),cMass,dMass); } else if ( (*fMap).first == "-f" ) { string dName(fName); fName = (*fMap).second; printf("[%s] - INFO - file name base updated to \042%s\042 (default \042%s\042)\n", moduleName.c_str(),fName.c_str(),dName.c_str()); } else { printf("[%s] - WARNING - Unknown option \042%s\042 with value <%s>\n", moduleName.c_str(),(*fMap).first.c_str(),(*fMap).second.c_str()); } } dataInput.clear(); } // random (seed) generator time_t seconds = time(0); // time since epoch TRandom1* randomEng = new TRandom1((unsigned int)seconds); int pySeed = (int)(randomEng->Rndm() * 900000000); if ( pySeed == 0 ) pySeed == 1; if ( pySeed > 900000000 ) pySeed = 900000000; // Generator. Process selection. LHC initialization. Histogram. Pythia pythia; // pythia.readString("PartonLevel:MI = on"); // pythia.readString("PartonLevel:ISR = on"); // pythia.readString("PartonLevel:FSR = on"); // pythia.readString("HadronLevel:Hadronize = on"); pythia.readString("SoftQCD:all = on"); pythia.readString("Tune:pp = 5"); // tune 4C (default) // pythia.readString("Tune:pp = 8"); // tune AM2 (ATLAS) ostringstream ostr; ostr << "Random:seed = " << pySeed; pythia.readString(ostr.str()); printf("[%s] - DEBUG - config string <%s>\n",moduleName.c_str(),ostr.str().c_str()); printf("[%s] - INFO - configured random number generator with seed %i\n", moduleName.c_str(),pySeed); //J0 0 17 //J1 17 35 //J2 35 70 //J3 70 140 //J4 140 280 //J5 280 560 //J6 560 1120 //J7 1120 2240 //J8 2240 inf pythia.init( 2212, 2212, cMass); // Pythia8::settings.listChanged(); // Begin event loop. Generate event. Skip if error. List first one. int totEvts(nEvts*nFile); int iFile(0); printf("[%s] - INFO - process %i events in %i files (total events %i)\n", moduleName.c_str(),nEvts,nFile,totEvts); int iEvt(0); while (iEvt < totEvts) { if ( !pythia.next() ) continue; _FLOAT_T_ pthat(pythia.info.pTHat()); if ( iEvt < 1 ) { pythia.info.list(); pythia.event.list(); } if ( iEvt%nEvts == 0 ) { if (iEvt!=0) { printf("[%s] - INFO - write tuple after %i events\n", moduleName.c_str(),iEvt); rTree->Write(); printf("[%s] - INFO - close file #%i\n", moduleName.c_str(),iFile); rFile->Close(); ++iFile; } string fN = fName; ostringstream xStr; xStr << fName << setfill('0') << setw(5) << iFile << ".root"; printf("[%s] - INFO - opening file #%i \042%s\042\n", moduleName.c_str(),iFile+1,xStr.str().c_str()); rFile = new TFile(xStr.str().c_str(),"RECREATE"); rTree = new TTree("MB_Py8","MinimumBias root tree", 1); printf("[%s] - INFO - setting branches in tree %s\n", moduleName.c_str(),rTree->GetName()); rTree->Branch ("Nentry", &Nentry, "Nentry/I"); rTree->Branch ("Npartons", &Npartons, "Npartons/I"); rTree->Branch ("Nparticles",&Nparticles,"Nparticles/I"); rTree->Branch ("ID", ID, "ID[Nentry]/I"); rTree->Branch ("Stat", Stat, "Stat[Nentry]/I"); #ifdef _SINGLE_PRECISION_ rTree->Branch ("Charge", Charge, "Charge[Nentry]/F"); rTree->Branch ("Px", Px, "Px[Nentry]/F"); rTree->Branch ("Py", Py, "Py[Nentry]/F"); rTree->Branch ("Pz", Pz, "Pz[Nentry]/F"); rTree->Branch ("P0", P0, "P0[Nentry]/F"); rTree->Branch ("Pm", Pm, "Pm[Nentry]/F"); rTree->Branch ("Pt", Pt, "Pt[Nentry]/F"); rTree->Branch ("Rap", Rap, "Rap[Nentry]/F"); rTree->Branch ("Phi", Phi, "Phi[Nentry]/F"); rTree->Branch ("Eta", Eta, "Eta[Nentry]/F"); #else rTree->Branch ("Charge", Charge, "Charge[Nentry]/D"); rTree->Branch ("Px", Px, "Px[Nentry]/D"); rTree->Branch ("Py", Py, "Py[Nentry]/D"); rTree->Branch ("Pz", Pz, "Pz[Nentry]/D"); rTree->Branch ("P0", P0, "P0[Nentry]/D"); rTree->Branch ("Pm", Pm, "Pm[Nentry]/D"); rTree->Branch ("Pt", Pt, "Pt[Nentry]/D"); rTree->Branch ("Rap", Rap, "Rap[Nentry]/D"); rTree->Branch ("Phi", Phi, "Phi[Nentry]/D"); rTree->Branch ("Eta", Eta, "Eta[Nentry]/D"); #endif } // loop partonic event Npartons = 0; for (int i=0;iFill(); ++iEvt; } pythia.statistics(); rTree->Write(); rFile->Close(); // timing tms* endTime = new tms(); clock_t eT = times(endTime); printf("[%s] - INFO - EndJob (user time/system time) (%10i,%10i) - return (%10i)\n", moduleName.c_str(),endTime->tms_utime,endTime->tms_stime,eT); printf("[%s] - INFO - EndJob differences (user time/system time) (%10i,%10i) - return (%10i)\n", moduleName.c_str(),endTime->tms_utime-startTime->tms_utime,endTime->tms_stime-startTime->tms_stime,eT-sT); return 0; }