analysis/EICTrigger/PHPythia6.C

0001 #include "PHPythia6.h"
0002 #include "PHPy6GenTrigger.h"
0003
0004 #include <phhepmc/PHHepMCGenEvent.h>
0005 #include <phhepmc/PHHepMCGenEventMap.h>
0006
0007 #include <fun4all/Fun4AllReturnCodes.h>
0008
0009 #include <phool/getClass.h>
0010 #include <phool/recoConsts.h>
0011 #include <phool/PHIODataNode.h>
0012 #include <phool/PHDataNode.h>
0013 #include <phool/PHObject.h>
0014 #include <phool/PHCompositeNode.h>
0015 #include <phool/PHNodeIterator.h>
0016 #include <phool/PHNodeReset.h>
0017 #include <phool/PHTimeStamp.h>
0018 #include <phool/PHRandomSeed.h>
0019
0020 #include <HepMC/PythiaWrapper.h>
0021 #include <HepMC/IO_HEPEVT.h>
0022 #include <HepMC/IO_GenEvent.h>
0023 #include <HepMC/IO_AsciiParticles.h>
0024 #include <HepMC/GenEvent.h>
0025
0026 #include <gsl/gsl_randist.h>
0027
0028 #include <iostream>
0029 #include <iomanip>
0030 #include <fstream>
0031 #include <sstream>
0032 #include <cstdlib>
0033 #include <stdlib.h>
0034 #include <ctime>
0035 #include <sys/time.h>
0036 #include <algorithm>
0037 #include <cctype>
0038
0039 #define pytune pytune_
0040 extern "C" int  pytune(int *itune);
0041
0042 using namespace std;
0043
0044 typedef PHIODataNode<PHObject> PHObjectNode_t;
0045
0046 PHPythia6::PHPythia6(const std::string &name):
0047   SubsysReco(name),
0048   _eventcount(0),
0049   _geneventcount(0),
0050   _configFile("phpythia6.cfg"),
0051   _save_ascii( false ),
0052   _filename_ascii("pythia_hepmc.dat"),
0053   _registeredTriggers(),
0054   _triggersOR(true),
0055   _triggersAND(false){
0056
0057   PHHepMCGenHelper::set_embedding_id(1); // default embedding ID to 1
0058 }
0059
0060 PHPythia6::~PHPythia6() {
0061   //gsl_rng_free (RandomGenerator);
0062 }
0063
0064 int PHPythia6::Init(PHCompositeNode *topNode) {
0065
0066   /* Create node tree */
0067   CreateNodeTree(topNode);
0068
0069   /* event numbering will start from 1 */
0070   _eventcount = 0;
0071
0072   /* HepMC/example_MyPythia.cc:
0073    *
0074    * Pythia 6.1 uses HEPEVT with 4000 entries and 8-byte floating point
0075    *  numbers. We need to explicitly pass this information to the
0076    *  HEPEVT_Wrapper.
0077    */
0078   HepMC::HEPEVT_Wrapper::set_max_number_entries(4000);
0079   HepMC::HEPEVT_Wrapper::set_sizeof_real(8);
0080
0081   /* set pythia random number seed (mandatory!) */
0082   int fSeed = PHRandomSeed();
0083   fSeed = abs(fSeed)%900000000;
0084
0085   if ( (fSeed>=0) && (fSeed<=900000000) ) {
0086     pydatr.mrpy[0] = fSeed;                   // set seed
0087     pydatr.mrpy[1] = 0;                       // use new seed
0088   } else {
0089     cout << PHWHERE << " ERROR: seed " << fSeed << " is not valid" << endl;
0090     exit(2);
0091   }
0092
0093   /* read pythia configuration and initialize */
0094   if (!_configFile.empty()) ReadConfig( _configFile );
0095
0096   /* Initialization done */
0097   return Fun4AllReturnCodes::EVENT_OK;
0098 }
0099
0100 int PHPythia6::End(PHCompositeNode *topNode) {
0101
0102   //........................................TERMINATION
0103   // write out some information from Pythia to the screen
0104   call_pystat( 1 );
0105
0106   //match pythia printout
0107   cout << " |                                                                "
0108        << "                                                 | " << endl;
0109   cout << "                         PHPythia6::End - " << _eventcount
0110        << " events passed trigger" << endl;
0111   cout << "                         Fraction passed: " << _eventcount
0112        << "/" << _geneventcount
0113        << " = " << _eventcount/float(_geneventcount) << endl;
0114   cout << " *-------  End PYTHIA Trigger Statistics  ------------------------"
0115        << "-------------------------------------------------* " << endl;
0116
0117   return Fun4AllReturnCodes::EVENT_OK;
0118 }
0119
0120 //__________________________________________________________
0121 int PHPythia6::ReadConfig(const string cfg_file) {
0122
0123   if ( cfg_file != "" ) _configFile = cfg_file;
0124   cout << "PHPythia6::read_config - Reading " << _configFile << endl;
0125
0126   ifstream infile( _configFile );
0127   if (infile.fail ()) {
0128     cout << "PHPythia6::read_config - Failed to open file " << _configFile << endl;
0129     exit(2);
0130   }
0131
0132   // initialize variables
0133   Int_t   _nevents(0);
0134   Float_t _roots(0);
0135   string  _proj;
0136   string  _targ;
0137   string  _frame;
0138
0139   string FullLine;      // a complete line in the config file
0140   string label;         // the label
0141
0142   int index = 999999;
0143   double value = 1e9;
0144
0145   // get one line first
0146   getline(infile, FullLine);
0147   while ( !infile.eof() )
0148   {
0149
0150     // skip lines that begin with #, or "//"
0151     if ( FullLine[0]=='#' || FullLine.substr(0, 2) == "//" )
0152     {
0153       getline(infile,FullLine);
0154       continue;
0155     }
0156
0157     // make FullLine an istringstream
0158     istringstream line( FullLine.c_str() );
0159
0160     // get label
0161     line >> label;
0162
0163     // to lower case
0164     std::transform(label.begin(), label.end(), label.begin(), (int(*)(int)) std::tolower);
0165
0166     // based on label, fill correct item
0167     if ( label == "nevents" )
0168     {
0169       line >> _nevents;
0170       cout << "nevents\t" << _nevents << endl;
0171     }
0172     else if ( label == "roots" )
0173     {
0174       line >> _roots;
0175       cout << "roots\t" << _roots << endl;
0176     }
0177     else if ( label == "proj" )
0178     {
0179       line >> _proj;
0180       cout << "proj\t" << _proj << endl;
0181     }
0182     else if ( label == "targ" )
0183     {
0184       line >> _targ;
0185       cout << "targ\t" << _targ << endl;
0186     }
0187     else if ( label == "frame" )
0188     {
0189       line >> _frame;
0190       cout << "frame\t" << _frame << endl;
0191     }
0192     else if ( label == "p1" || label == "p2")
0193     {
0194       if (label == "p1") //Momentum of Projectile Beam (e- for e-p)
0195     {
0196       line >> index >> value;
0197       //Index Options(3MOM): 1 = x-momentum; 2 = y-momentum; 3 = z-momentum
0198       pyjets.p[index-1][0] = value;
0199       cout << "p1\t" << index << " " << value << endl;
0200     }
0201       if (label == "p2") //Momentum of Target Beam (p for e-p)
0202     {
0203       line >> index >> value;
0204       //Index Options(3MOM): 1 = x-momentum; 2 = y-momentum; 3 = z-momentum
0205       pyjets.p[index-1][1] = value;
0206       cout << "p2\t" << index << " " << value << endl;
0207     }
0208       /*
0209       int entry = 0;
0210       if ( label=="p1") entry = 1;
0211       else if ( label=="p2") entry = 2;
0212
0213       char temp_line[10000];
0214       strcpy(temp_line,FullLine.c_str());
0215       char *sptr = strtok(temp_line," \t"); // skip first word
0216       sptr = strtok(NULL," \t");
0217       cout << label;
0218       int index = 1;
0219       while ( sptr != NULL )
0220         {
0221           double val = atof(sptr);
0222       cout << "Entry: " << entry << endl;
0223       index++;
0224           cout << "\t" << val;
0225           sptr = strtok(NULL," \t");
0226         }
0227     cout << endl;*/
0228     }
0229     else if ( label == "msel" )
0230       {
0231     line >> value;
0232     pysubs.msel= (int) value;
0233     cout << "msel\t" << value << endl;
0234     IntegerTest(value);
0235       }
0236     else if ( label == "msub" )
0237     {
0238       line >> index >> value;
0239       // careful with C/F77 differences: arrays in C start at 0, F77 at 1,
0240       // so we need to subtract 1 from the process #)
0241       pysubs.msub[index-1] = (int) value;
0242       cout << "msub\t" << index << " " << value << endl;
0243       IntegerTest(value);
0244     }
0245     else if ( label == "mstp" )
0246       {
0247     line >> index >> value;
0248     pypars.mstp[index-1] = (int) value;
0249     cout << "mstp\t" << index << " " << value << endl;
0250     IntegerTest(value);
0251       }
0252     else if ( label == "mstj" )
0253       {
0254     line >> index >> value;
0255     pydat1.mstj[index-1] = (int) value;
0256     cout << "mstj\t" << index << " " << value << endl;
0257     IntegerTest(value);
0258       }
0259     else if ( label == "mstu" )
0260       {
0261     line >> index >> value;
0262     pydat1.mstu[index-1] = (int) value;
0263     cout << "mstu\t" << index << " " << value << endl;
0264     IntegerTest(value);
0265       }
0266     else if ( label == "ckin" )
0267       {
0268     line >> index >> value;
0269     pysubs.ckin[index-1] =value;
0270     cout << "ckin\t" << index << " " << value << endl;
0271       }
0272     else if ( label == "parp" )
0273       {
0274     line >> index >> value;
0275     pypars.parp[index-1] = value;
0276     cout << "parp\t" << index << " " << value << endl;
0277       }
0278     else if ( label == "parj" )
0279       {
0280     line >> index >> value;
0281     pydat1.parj[index-1] = value;
0282     cout << "parj\t" << index << " " << value << endl;
0283       }
0284     else if ( label == "paru" )
0285       {
0286     line >> index >> value;
0287     pydat1.paru[index-1] = value;
0288     cout << "paru\t" << index << " " << value << endl;
0289       }
0290     else if ( label == "parf" )
0291       {
0292     line >> index >> value;
0293     pydat2.parf[index-1] = value;
0294     cout << "parf\t" << index << " " << value << endl;
0295       }
0296     else if ( label == "mdme" )
0297     {
0298       int idc = 0;          // decay channel
0299       line >> idc >> index >> value;
0300
0301       // if (ivalue==1/0) turn on/off decay channel idc
0302       pydat3.mdme[index-1][idc-1] = value;
0303       cout << "mdme\t" << idc << " " << index << " " << value << endl;
0304     }
0305     else if ( label == "pmas" )
0306     {
0307       int idc = 0;
0308       line >> idc >> index >> value;
0309
0310       pydat2.pmas[index-1][idc-1] = value;
0311       cout << "pmas\t" << idc << " " << index << " " << value << endl;
0312     }
0313     else if ( label == "pytune" )
0314     {
0315       int ivalue;
0316       line >> ivalue;
0317       pytune(&ivalue);
0318       cout << "pytune\t" << ivalue << endl;
0319     }
0320      else
0321       {
0322     // label was not understood
0323     cout << "************************************************************" << endl;
0324     cout << "PHPythia6::ReadConfig(), ERROR this option is not supported: " << FullLine << endl;
0325     cout << "************************************************************" << endl;
0326       }
0327
0328     // get next line in file
0329     getline( infile, FullLine );
0330   }
0331
0332   // Call pythia initialization
0333   call_pyinit( _frame.c_str(), _proj.c_str(), _targ.c_str(), _roots );
0334
0335   //call_pylist(12);
0336
0337   infile.close();
0338
0339   return _nevents;
0340 }
0341
0342 //-* print pythia config info
0343 void PHPythia6::print_config() const {
0344 }
0345
0346 int PHPythia6::process_event(PHCompositeNode *topNode) {
0347
0348   if (verbosity > 1) cout << "PHPythia6::process_event - event: " << _eventcount << endl;
0349
0350   bool passedTrigger = false;
0351   std::vector<bool> theTriggerResults;
0352   int genCounter = 0;
0353
0354   /* based on HepMC/example_MyPythia.cc
0355    *........................................HepMC INITIALIZATIONS
0356    *
0357    * Instantiate an IO strategy for reading from HEPEVT. */
0358   HepMC::IO_HEPEVT hepevtio;
0359   HepMC::GenEvent* evt;
0360
0361   while (!passedTrigger) {
0362     ++genCounter;
0363
0364     call_pyevnt();      // generate one event with Pythia
0365     _geneventcount++;
0366     // pythia pyhepc routine converts common PYJETS in common HEPEVT
0367     call_pyhepc( 1 );
0368     evt = hepevtio.read_next_event();
0369
0370     // define the units (Pythia uses GeV and mm)
0371     evt->use_units(HepMC::Units::GEV, HepMC::Units::MM);
0372
0373     // add some information to the event
0374     evt->set_event_number(_eventcount);
0375
0376     /* process ID from pythia */
0377     evt->set_signal_process_id(pypars.msti[1-1]);
0378
0379     // set number of multi parton interactions
0380     evt->set_mpi( pypars.msti[31-1] );
0381
0382     // set cross section information
0383     evt->set_cross_section( HepMC::getPythiaCrossSection() );
0384
0385     // Set the PDF information
0386     HepMC::PdfInfo pdfinfo;
0387     pdfinfo.set_x1(pypars.pari[33-1]);
0388     pdfinfo.set_x2(pypars.pari[34-1]);
0389     pdfinfo.set_scalePDF(pypars.pari[22-1]);
0390     pdfinfo.set_id1(pypars.msti[15-1]);
0391     pdfinfo.set_id2(pypars.msti[16-1]);
0392     evt->set_pdf_info(pdfinfo);
0393
0394     // test trigger logic
0395
0396     bool andScoreKeeper = true;
0397     if (verbosity > 2) {
0398       cout << "PHPythia6::process_event - triggersize: " << _registeredTriggers.size() << endl;
0399     }
0400
0401     for (unsigned int tr = 0; tr < _registeredTriggers.size(); tr++) {
0402       bool trigResult = _registeredTriggers[tr]->Apply(evt);
0403
0404       if (verbosity > 2) {
0405     cout << "PHPythia6::process_event trigger: "
0406          << _registeredTriggers[tr]->GetName() << "  " << trigResult << endl;
0407       }
0408
0409       if (_triggersOR && trigResult) {
0410     passedTrigger = true;
0411     break;
0412       } else if (_triggersAND) {
0413     andScoreKeeper &= trigResult;
0414       }
0415
0416       if (verbosity > 2 && !passedTrigger) {
0417     cout << "PHPythia8::process_event - failed trigger: "
0418          << _registeredTriggers[tr]->GetName() <<  endl;
0419       }
0420     }
0421
0422     if ((andScoreKeeper && _triggersAND) || (_registeredTriggers.size() == 0)) {
0423       passedTrigger = true;
0424       genCounter = 0;
0425     }
0426
0427     // delete failed events
0428     if(!passedTrigger) delete evt;
0429
0430   }
0431
0432   /* write the event out to the ascii files */
0433   if ( _save_ascii )
0434   {
0435     HepMC::IO_GenEvent ascii_io(_filename_ascii.c_str(),std::ios::out);
0436     ascii_io << evt;
0437   }
0438
0439   /* pass HepMC to PHNode*/
0440
0441   PHHepMCGenEvent * success = PHHepMCGenHelper::  insert_event(evt);
0442   if (!success) {
0443     cout << "PHPythia6::process_event - Failed to add event to HepMC record!" << endl;
0444     return Fun4AllReturnCodes::ABORTRUN;
0445   }
0446   /* print outs*/
0447   if (verbosity > 2) cout << "PHPythia6::process_event - FINISHED WHOLE EVENT" << endl;
0448
0449   ++_eventcount;
0450   return Fun4AllReturnCodes::EVENT_OK;
0451 }
0452
0453
0454
0455 void PHPythia6::IntegerTest(double number ) {
0456
0457   if (fmod(number, 1.0) != 0) {
0458     cout << "Warning: Value " << number << " is not an integer." << endl;
0459     cout << "This parameter requires an integer value." << endl;
0460     cout << "Value of parameter truncated to " << (int) number  << endl;
0461
0462     //...End simulation if a double value is input for an integer parameter
0463     //    throw Fun4AllReturnCodes::ABORTRUN;
0464   }
0465   return;
0466 }
0467
0468 int PHPythia6::ResetEvent(PHCompositeNode *topNode) {
0469   return Fun4AllReturnCodes::EVENT_OK;
0470 }
0471
0472 void PHPythia6::register_trigger(PHPy6GenTrigger *theTrigger) {
0473   if(verbosity > 1) cout << "PHPythia6::registerTrigger - trigger " << theTrigger->GetName() << " registered" << endl;
0474   _registeredTriggers.push_back(theTrigger);
0475 }