CMS 3D CMS Logo
edm::AlpgenSource

edm::AlpgenSource Class Reference

#include <GeneratorInterface/AlpgenInterface/interface/AlpgenSource.h>

Inheritance diagram for edm::AlpgenSource:

edm::ExternalInputSource edm::ConfigurableInputSource edm::InputSource edm::ProductRegistryHelper

List of all members.

Public Member Functions

 AlpgenSource (const ParameterSet &, const InputSourceDescription &)
 Constructor.
virtual ~AlpgenSource ()
 Destructor.

Private Member Functions

void beginRun (Run &r)
bool call_pygive (const std::string &iParm)
 Interface to the PYGIVE/TXGIVE pythia routine, with add'l protections.
bool call_txgive (const std::string &iParm)
void clear ()
virtual bool produce (Event &e)

Private Attributes

bool doubleParticle
double etamax
double etamin
HepMC::GenEvent * evt
std::string fileName_
CLHEP::HepRandomEngine * fRandomEngine
CLHEP::RandFlat * fRandomGenerator
LHERunInfoProduct::Header lheAlpgenUnwParHeader
unsigned int maxEventsToPrint_
 Events to print if verbosity.
int Nev_
int particleID
double phimax
double phimin
double ptmax
double ptmin
bool pythiaHepMCVerbosity_
 HepMC verbosity flag.
unsigned int pythiaPylistVerbosity_
 Pythia PYLIST Verbosity flag.

Detailed Description

Definition at line 35 of file AlpgenSource.h.

Constructor & Destructor Documentation

AlpgenSource::AlpgenSource ( const ParameterSet &  pset,
const InputSourceDescription &  desc 
)

Constructor.

Definition at line 51 of file AlpgenSource.cc.

References LHERunInfoProduct::Header::addLine(), call_pygive(), call_txgive(), edm::errors::Configuration, GenMuonPlsPt100GeV_cfg::cout, lat::endl(), fileName_, edm::ExternalInputSource::fileNames(), fRandomEngine, edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), lheAlpgenUnwParHeader, edm::InputSource::maxEvents(), maxEventsToPrint_, Nev_, pars, pythiaHepMCVerbosity_, pythiaPylistVerbosity_, and randomEngine.

00052                                                                  :
00053   ExternalInputSource(pset, desc, false), evt(0), 
00054   pythiaPylistVerbosity_ (pset.getUntrackedParameter<int>("pythiaPylistVerbosity",0)),
00055   pythiaHepMCVerbosity_ (pset.getUntrackedParameter<bool>("pythiaHepMCVerbosity",false)),
00056   maxEventsToPrint_ (pset.getUntrackedParameter<int>("maxEventsToPrint",1)),
00057   lheAlpgenUnwParHeader("AlpgenUnwParFile")
00058 {
00059   
00060   fileName_ = fileNames()[0];
00061   // strip the file: 
00062   if ( fileName_.find("file:") || fileName_.find("rfio:")){
00063     fileName_.erase(0,5);
00064   }   
00065 
00066   // open the .unw file to store additional 
00067   // informations in the AlpgenInfoProduct
00068 //  unwfile = new ifstream((fileName_+".unw").c_str());
00069   // get the number of input events from  _unw.par files
00070   char buffer[256];
00071   ifstream reader((fileName_+"_unw.par").c_str());
00072   char sNev[80];
00073   //  lheAlpgenUnwParHeader.addLine("\n");
00074   while ( reader.getline (buffer,256) ) {
00075     istringstream is(buffer);
00076     lheAlpgenUnwParHeader.addLine(std::string(buffer) + "\n");
00077     is >> sNev;
00078     Nev_ = atoi(sNev);
00079   }
00080 
00081   //check that N(asked events) <= N(input events)
00082   if(maxEvents()>Nev_) {
00083     edm::LogInfo("Generator|TooLittleData") << "ALPGEN warning: Number of events requested > Number of unweighted events.\n"
00084                                             << "                Execution will stop after processing the last unweighted event";
00085   }
00086 
00087   if(maxEvents() != -1 && maxEvents() < Nev_) // stop at N(asked events) if N(asked events)<N(input events)
00088     Nev_ = maxEvents();
00089   
00090   // PYLIST Verbosity Level
00091   // Valid PYLIST arguments are: 1, 2, 3, 5, 7, 11, 12, 13
00092   pythiaPylistVerbosity_ = pset.getUntrackedParameter<int>("pythiaPylistVerbosity",0);
00093   
00094   // HepMC event verbosity Level
00095   pythiaHepMCVerbosity_ = pset.getUntrackedParameter<bool>("pythiaHepMCVerbosity",false);
00096   
00097   //Max number of events printed on verbosity level 
00098   maxEventsToPrint_ = pset.getUntrackedParameter<int>("maxEventsToPrint",0);
00099   
00100   // Set PYTHIA parameters in a single ParameterSet
00101   {
00102     ParameterSet pythia_params = 
00103       pset.getParameter<ParameterSet>("PythiaParameters") ;
00104     
00105     // Read the PYTHIA parameters for each set of parameters
00106     vector<string> pars = 
00107       pythia_params.getParameter<vector<string> >("pythia");
00108     
00109     // Loop over all parameters and stop in case of mistake
00110     for( vector<string>::const_iterator  
00111            itPar = pars.begin(); itPar != pars.end(); ++itPar ) {
00112       static string sRandomValueSetting("MRPY(1)");
00113       if( 0 == itPar->compare(0,sRandomValueSetting.size(),sRandomValueSetting) ) {
00114         throw edm::Exception(edm::errors::Configuration,"PythiaError")
00115           <<" attempted to set random number seed.";
00116       }
00117       if( ! call_pygive(*itPar) ) {
00118         throw edm::Exception(edm::errors::Configuration,"PythiaError") 
00119           <<" pythia did not accept the following \""<<*itPar<<"\"";
00120       }
00121     }
00122   }
00123 
00124   // Read the Alpgen parameters
00125 
00126   // read External Generator parameters
00127   {   ParameterSet generator_params = 
00128         pset.getParameter<ParameterSet>("GeneratorParameters") ;
00129   vector<string> pars = 
00130     generator_params.getParameter<vector<string> >("generator");
00131   for( vector<string>::const_iterator  
00132          itPar = pars.begin(); itPar != pars.end(); ++itPar ) 
00133     {
00134       call_txgive(*itPar);          
00135     }
00136   // giving to txgive a string with the alpgen filename
00137   string tmpstring = "UNWFILE = " + fileName_;
00138   call_txgive(tmpstring);
00139   }
00140   
00141   
00142   //In the future, we will get the random number seed on each event and tell 
00143   // pythia to use that new seed
00144   edm::Service<RandomNumberGenerator> rng;
00145   randomEngine = fRandomEngine = &(rng->getEngine());
00146   uint32_t seed = rng->mySeed();
00147   ostringstream sRandomSet;
00148   sRandomSet <<"MRPY(1)="<<seed;
00149   call_pygive(sRandomSet.str());
00150   
00151   //  call_pretauola(-1);     // TAUOLA initialization
00152   call_pyinit( "USER", "p", "p", 14000. );
00153   
00154   cout << endl; // Stetically add for the output
00155   //********                                      
00156   
00157   produces<HepMCProduct>();
00158   produces<LHEEventProduct>();
00159 
00160 //  produces<AlpWgtFileInfoProduct, edm::InRun>();
00161   produces<LHERunInfoProduct, edm::InRun>();
00162 }

AlpgenSource::~AlpgenSource (  )  [virtual]

Destructor.

Definition at line 165 of file AlpgenSource.cc.

References alpgen_end, and clear().

00165                            {
00166   call_pystat(1);
00167   //  call_pretauola(1);  // output from TAUOLA 
00168   alpgen_end();
00169   clear(); 
00170 }

Member Function Documentation

void AlpgenSource::beginRun ( Run &  r  )  [private, virtual]

Reimplemented from edm::ConfigurableInputSource.

Definition at line 176 of file AlpgenSource.cc.

References LHERunInfoProduct::Header::addLine(), indexGen::comments, fileName_, lheAlpgenUnwParHeader, edm::Run::put(), lhef::CommonBlocks::readHEPRUP(), and edm::runInfo.

00176                                    {
00177   // get the LHE init information from Fortran code
00178   lhef::HEPRUP heprup;
00179   lhef::CommonBlocks::readHEPRUP(&heprup);
00180   auto_ptr<LHERunInfoProduct> runInfo(new LHERunInfoProduct(heprup));
00181 
00182   // information on weighted events
00183   LHERunInfoProduct::Header lheAlpgenWgtHeader("AlpgenWgtFile");
00184   //  lheAlpgenWgtHeader.addLine("\n");
00185   ifstream wgtascii((fileName_+".wgt").c_str());
00186   char buffer[512];
00187   while(wgtascii.getline(buffer,512)) {
00188     lheAlpgenWgtHeader.addLine(std::string(buffer) + "\n");
00189   }
00190 
00191   LHERunInfoProduct::Header lheAlpgenParHeader("AlpgenParFile");
00192   //  lheAlpgenParHeader.addLine("\n");
00193   ifstream parascii((fileName_+".par").c_str());
00194   while(parascii.getline(buffer,512)) {
00195     lheAlpgenParHeader.addLine(std::string(buffer) + "\n");
00196   }
00197 
00198   // comments on top
00199   LHERunInfoProduct::Header comments;
00200   comments.addLine("\n");
00201   comments.addLine("\tExtracted by AlpgenInterface\n");
00202 
00203   // build the final Run info object
00204   runInfo->addHeader(comments);
00205   runInfo->addHeader(lheAlpgenUnwParHeader);
00206   runInfo->addHeader(lheAlpgenWgtHeader);
00207   runInfo->addHeader(lheAlpgenParHeader);
00208   r.put(runInfo);
00209 }

bool AlpgenSource::call_pygive ( const std::string &  iParm  )  [private]

Interface to the PYGIVE/TXGIVE pythia routine, with add'l protections.

Definition at line 320 of file AlpgenSource.cc.

References pydat1, and PYGIVE.

Referenced by AlpgenSource().

00320                                                  {
00321 
00322   int numWarn = pydat1.mstu[26]; //# warnings
00323   int numErr = pydat1.mstu[22];// # errors
00324   
00325 //call the fortran routine pygive with a fortran string
00326   PYGIVE( iParm.c_str(), iParm.length() );  
00327   //  PYGIVE( iParm );  
00328 //if an error or warning happens it is problem
00329   return pydat1.mstu[26] == numWarn && pydat1.mstu[22] == numErr;   
00330 }

bool AlpgenSource::call_txgive ( const std::string &  iParm  )  [private]

Definition at line 333 of file AlpgenSource.cc.

References TXGIVE.

Referenced by AlpgenSource().

00334    {
00335     //call the fortran routine txgive with a fortran string
00336     TXGIVE( iParm.c_str(), iParm.length() );  
00337     return 1;  
00338    }

void AlpgenSource::clear ( void   )  [private]

Definition at line 172 of file AlpgenSource.cc.

Referenced by ~AlpgenSource().

00172                          {
00173   
00174 }

bool AlpgenSource::produce ( Event &  e  )  [private, virtual]

Implements edm::ConfigurableInputSource.

Definition at line 211 of file AlpgenSource.cc.

References lhef::HEPEUP::AQCDUP, lhef::HEPEUP::AQEDUP, call_pylist(), conv, GenMuonPlsPt100GeV_cfg::cout, lat::endl(), edm::ConfigurableInputSource::event(), evt, i, maxEventsToPrint_, Nev_, edm::ConfigurableInputSource::numberEventsInRun(), lhef::HEPEUP::NUP, edm::Event::put(), pyint1, pypars, pythiaHepMCVerbosity_, pythiaPylistVerbosity_, lhef::CommonBlocks::readHEPEUP(), edm::InputSource::remainingEvents(), and lhef::HEPEUP::SPINUP.

00211                                     {
00212   
00213   // exit if N(events asked) has been exceeded
00214   if(event()> Nev_) {
00215     return false;
00216   } else {
00217     
00218     auto_ptr<HepMCProduct> bare_product(new HepMCProduct());  
00219     
00220 //    // Additional information from unweighted file
00221 //    auto_ptr<AlpgenInfoProduct> alp_product(new AlpgenInfoProduct());
00222 
00223     // Extract from .unw file the info for AlpgenInfoProduct
00224     
00225 //    char buffer[512];
00226 //    if(unwfile->getline(buffer,512)) {
00227 //      alp_product->EventInfo(buffer);
00228 //    }
00229 //    if(unwfile->getline(buffer,512)) 
00230 //      alp_product->InPartonInfo(buffer);
00231 //    if(unwfile->getline(buffer,512)) 
00232 //      alp_product->InPartonInfo(buffer);
00233 //    for(int i_out = 0; i_out <  alp_product->nTot()-2; i_out++) {
00234 //      if(unwfile->getline(buffer,512)) 
00235 //      alp_product->OutPartonInfo(buffer);
00236 //    }
00237 
00238     call_pyevnt();      // generate one event with Pythia
00239     //        call_pretauola(0);  // tau-lepton decays with TAUOLA 
00240 
00241     // fill the parton level LHE event information
00242     lhef::HEPEUP hepeup;
00243     lhef::CommonBlocks::readHEPEUP(&hepeup);
00244     hepeup.AQEDUP = hepeup.AQCDUP = -1.0; // alphas are not saved by Alpgen
00245     for(int i = 0; i < hepeup.NUP; i++)
00246       hepeup.SPINUP[i] = -9;    // Alpgen does not store spin information
00247     auto_ptr<LHEEventProduct> lheEvent(new LHEEventProduct(hepeup));
00248     // Moved to later, first we have to make sure that the event generated by Pythia is OK.
00249     // e.put(lheEvent);
00250 
00251     // The way this is implemented, call_pyevnt() will keep reading events from the .unw until
00252     // one of them is able to pass the matching veto. If it can't read a last event, it will
00253     // simply continue with whatever event was in memory. In this case, there are two possibilities:
00254     // A) The event passed the matching veto. There is meaningful information in HEPEUP. Continue.
00255     // (Will fail next Event with HEPEUP.NUP == 0. Harmless.)
00256     // B) The event didn't pass. HEPEUP.NUP == 0. This should never happen, and we should abort here.
00257     if(hepeup.NUP==0) {
00258       edm::LogInfo("Generator|TooLittleData") << "ALPGEN warning: last unweighted event reached.\n"
00259                                               << "                (hepeup.NUP == 0)\n"
00260                                               << "                The event number " << event() << " will not be written to disk.";  
00261       return false;
00262     }
00263     
00264     call_pyhepc( 1 );
00265     //    HepMC::GenEvent* evt = conv.getGenEventfromHEPEVT();
00266     HepMC::GenEvent* evt = conv.read_next_event();
00267     
00268     evt->set_signal_process_id(pypars.msti[0]);
00269     evt->set_event_number(numberEventsInRun() - remainingEvents() - 1);
00270     
00271     int id1 = pyint1.mint[14];
00272     int id2 = pyint1.mint[15];
00273     if ( id1 == 21 ) id1 = 0;
00274     if ( id2 == 21 ) id2 = 0; 
00275     double x1 = pyint1.vint[40];
00276     double x2 = pyint1.vint[41];  
00277     double Q  = pyint1.vint[50];
00278     double pdf1 = pyint1.vint[38];
00279     pdf1 /= x1 ;
00280     double pdf2 = pyint1.vint[39];
00281     pdf2 /= x2 ;
00282     evt->set_pdf_info( HepMC::PdfInfo(id1,id2,x1,x2,Q,pdf1,pdf2) ) ;
00283     
00284     evt->weights().push_back( pyint1.vint[96] );
00285 
00286     //******** Verbosity ********
00287     
00288     if(event() <= maxEventsToPrint_ &&
00289        (pythiaPylistVerbosity_ || pythiaHepMCVerbosity_)) {
00290       
00291       // Prints PYLIST info
00292       if(pythiaPylistVerbosity_) {
00293         call_pylist(pythiaPylistVerbosity_);
00294       }
00295       
00296       // Prints HepMC event
00297       if(pythiaHepMCVerbosity_) {
00298         cout << "Event process = " << pypars.msti[0] << endl 
00299              << "----------------------" << endl;
00300         evt->print();
00301       }
00302     }
00303 
00304     if(evt)  bare_product->addHepMCData(evt );
00305 
00306     // Last check to make sure there are no empty events. This supersedes the
00307     // AlpgenEmptyEventFilter.
00308     if(!(evt->particles_empty())) {
00309       e.put(lheEvent);
00310       e.put(bare_product);
00311       return true;
00312     }
00313     else
00314       return false;
00315   
00316   }
00317 }

Member Data Documentation

bool edm::AlpgenSource::doubleParticle [private]

Definition at line 72 of file AlpgenSource.h.

double edm::AlpgenSource::etamax [private]

Definition at line 74 of file AlpgenSource.h.

double edm::AlpgenSource::etamin [private]

Definition at line 74 of file AlpgenSource.h.

HepMC::GenEvent* edm::AlpgenSource::evt [private]

Definition at line 59 of file AlpgenSource.h.

Referenced by produce().

std::string edm::AlpgenSource::fileName_ [private]

Definition at line 68 of file AlpgenSource.h.

Referenced by AlpgenSource(), and beginRun().

CLHEP::HepRandomEngine* edm::AlpgenSource::fRandomEngine [private]

Definition at line 77 of file AlpgenSource.h.

Referenced by AlpgenSource().

CLHEP::RandFlat* edm::AlpgenSource::fRandomGenerator [private]

Definition at line 78 of file AlpgenSource.h.

LHERunInfoProduct::Header edm::AlpgenSource::lheAlpgenUnwParHeader [private]

Definition at line 80 of file AlpgenSource.h.

Referenced by AlpgenSource(), and beginRun().

unsigned int edm::AlpgenSource::maxEventsToPrint_ [private]

Events to print if verbosity.

Definition at line 66 of file AlpgenSource.h.

Referenced by AlpgenSource(), and produce().

int edm::AlpgenSource::Nev_ [private]

Definition at line 49 of file AlpgenSource.h.

Referenced by AlpgenSource(), and produce().

int edm::AlpgenSource::particleID [private]

Definition at line 71 of file AlpgenSource.h.

double edm::AlpgenSource::phimax [private]

Definition at line 75 of file AlpgenSource.h.

double edm::AlpgenSource::phimin [private]

Definition at line 75 of file AlpgenSource.h.

double edm::AlpgenSource::ptmax [private]

Definition at line 73 of file AlpgenSource.h.

double edm::AlpgenSource::ptmin [private]

Definition at line 73 of file AlpgenSource.h.

bool edm::AlpgenSource::pythiaHepMCVerbosity_ [private]

HepMC verbosity flag.

Definition at line 64 of file AlpgenSource.h.

Referenced by AlpgenSource(), and produce().

unsigned int edm::AlpgenSource::pythiaPylistVerbosity_ [private]

Pythia PYLIST Verbosity flag.

Definition at line 62 of file AlpgenSource.h.

Referenced by AlpgenSource(), and produce().

The documentation for this class was generated from the following files:
Generated on Tue Jun 9 18:37:44 2009 for CMSSW by  doxygen 1.5.4