agile is hosted by Hepforge, IPPP Durham

Rapgap.cc

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 #include "AGILe/Rapgap/Rapgap.hh"
00003 #include "AGILe/Rapgap/RapgapWrapper32.hh"
00004 //#include "AGILe/FPythia/FPythia.hh"
00005 #include "AGILe/FPythia/FPythiaWrapper.hh"
00006 #include "AGILe/Utils.hh"
00007 #include "AGILe/HepMCTools.hh"
00008 
00009 #include "HepMC/HEPEVT_Wrapper.h"
00010 
00011 
00012 //-----------------------------------------------------------------------------
00013 // Implementation file for Rapgap helper class.
00014 //
00015 // Authors: Doug Gingrich
00016 //          Andy Buckley
00017 //          Holger Schulz
00018 //-----------------------------------------------------------------------------
00019 
00020 
00021 namespace AGILe {
00022 
00023 
00024   // Standard constructor
00025   Rapgap::Rapgap() {
00026     HepMC::HEPEVT_Wrapper::set_sizeof_real(8);
00027 
00028     _particleNames[ELECTRON]    = "e-";
00029     _particleNames[POSITRON]    = "e+";
00030     _particleNames[PROTON]      = "p+";
00031     _particleNames[ANTIPROTON]  = "pbar-";
00032     _particleNames[NEUTRON]     = "n0";
00033     _particleNames[ANTINEUTRON] = "nbar0";
00034     _particleNames[PHOTON]      = "gamma";
00035     _particleNames[MUON]        = "mu-";
00036     _particleNames[ANTIMUON]    = "mu+";
00037     // TAU, ANTITAU
00038     _particleNames[NU_E]        = "nu_e";
00039     _particleNames[NU_EBAR]     = "nu_ebar";
00040     _particleNames[NU_MU]       = "nu_mu";
00041     _particleNames[NU_MUBAR]    = "nu_mubar";
00042     _particleNames[NU_TAU]      = "nu_tau";
00043     _particleNames[NU_TAUBAR]   = "nu_taubar";
00044     _particleNames[PIPLUS]      = "pi+";
00045     _particleNames[PIMINUS]     = "pi-";
00046 
00047     // Start counting events at 1.
00048     _nevt = 1;
00049   }
00050 
00052   void Rapgap::initialize() {
00053     Generator::initialize();
00054     //Rapgap::setParam("LEPI", "2212");
00055     //Rapgap::setParam("HADI", "2212");
00056       //int ISEED = 213123;
00057       //int N1 = 0;
00058       //int N2 = 0;
00059       //call_rm48in(ISEED, N1, N2);
00060     //call_rapgap();
00061     // Initialise Rapgap parameters
00062     getLog() << Log::INFO << "Now calling GRAINI..." << endl;
00063     call_graini();
00064     getLog() << Log::INFO << "GRAINI is done..." << endl;
00065 
00066     // Switch on QCD processes by default
00067     getLog() << Log::INFO << "Setup basic processes..." << endl;
00068     //Rapgap::setParam("IPRO", 18);
00069     //Rapgap::setParam("IDIR", 0);  // Don't do DIS but diffractive stuff by default
00070     //Rapgap::setParam("INTER", 2);  // Switch to charged current interaction
00074 
00075     getLog() << Log::INFO << "Now calling Steer..." << endl;
00076     getLog() << Log::INFO << "In case nothing happens from this point on," << endl;
00077     getLog() << Log::INFO << "make sure you specified a cascade steering file on CL:" << endl;
00078     getLog() << Log::INFO << "agile-runmc RapgapPythia6:X:Y < steer_pp ... " << endl;
00079     call_steer();
00080     //getLog() << Log::INFO << "Steer is done..." << endl;
00081 
00082     getLog() << Log::INFO << "Change standard RAPGAP parameters..." << endl;
00083     call_rapcha();
00084     //getLog() << Log::INFO << "Change standard HERACLES parameters..." << endl;
00085     call_hercha();
00086     //getLog() << Log::INFO << "Change standard PYTHIA parameters..." << endl;
00087     call_pytcha();
00088 
00089     //getLog() << Log::INFO << "Change standard ARIADNE parameters..." << endl;
00090     //call_aricha();
00091     //getLog() << Log::INFO << "Init ARIADNE..." << endl;
00092     //const char* NAME = "RAPGAP";
00093     //call_arinit(NAME);
00094 
00095     getLog() << Log::INFO << "Calling PTIME..." << endl;
00096     const char* name = "rapgap";
00097     int j = 1;
00098     int k = 0;
00099     call_ptime(name, j, k);
00100 
00101     getLog() << Log::INFO << "Now perfroming basic cross-section integration..." << endl;
00102     call_rapgap();
00103     getLog() << Log::INFO << "RAPGAP is done..." << endl;
00104     call_raend(1);
00105     _initialized = true;
00106   }
00107 
00108   void Rapgap::setGenSpecificInitialState(PdgCode p1, double e1, PdgCode p2, double e2) {
00109     getLog() << Log::INFO << "Setting initial state..." << endl;
00110     // For Pythia, we have to hold on to this state until initialize() is called.
00111     //_particleName1 = _particleNames[p1];
00112     //_particleName2 = _particleNames[p2];
00114     //Rapgap::setParam("KBE1", _particleNames[ParticleName(p1)]);
00115     //Rapgap::setParam("KBE2", _particleNames[ParticleName(p2)]);
00116     Rapgap::setParam("KBE1", p1);
00117     Rapgap::setParam("KBE2", p2);
00118     Rapgap::setParam("PBE1", e1);
00119     Rapgap::setParam("PBE2", e2);
00120     //Rapgap::setParam("PLEP", e1);
00121     //Rapgap::setParam("PPIN", e2);
00122     //Rapgap::setParam("LEPI", p1);
00123     //Rapgap::setParam("HADI", p2);
00124   }
00125 
00126   // Run the generator for one event
00127   void Rapgap::makeEvent(HepMC::GenEvent& evt) {
00128     Generator::makeEvent(evt);
00129     call_event();
00130     call_pyhepc(1);
00131     // fill the event using HepMC
00132     Rapgap::fillEvent(evt);
00133     evt.set_event_number(_nevt);
00134     // Increment an event counter (Pythia does not count for itself).
00135     _nevt++;
00136   }
00137 
00138   // Fill a HepMC event
00139   void Rapgap::fillEvent(HepMC::GenEvent& evt) {
00140     _hepevt.fill_next_event(&evt);
00141     fixHepMCUnitsFromGeVmm(evt);
00142     #ifdef HEPMC_HAS_CROSS_SECTION
00143     const double xsecval = getCrossSection();
00144     const double xsecerr = FC_EFFIC.SD;
00145     getLog() << Log::DEBUG << "Writing cross-section = " << xsecval << " +- " << xsecerr << endl;
00146     HepMC::GenCrossSection xsec;
00147     xsec.set_cross_section(xsecval, xsecerr);
00148     evt.set_cross_section(xsec);
00149     #endif
00150   }
00151 
00152   const double Rapgap::getCrossSection(){
00153     _crossSection = FC_EFFIC.AVGI * 1000.0;
00154     return _crossSection;
00155   }
00156 
00157   // Set string parameters
00158   bool Rapgap::setParam(const string& name, const string& value) {
00159     Generator::setParam(name, value);
00160 
00162     if (name == "PPIN") {
00163       getLog() << Log::INFO << "Setting ... (PPIN) = " << value << endl;
00164       FC_INPU.PIN = -1*asDouble(value);
00165     } else if (name == "PLEP") {
00166       getLog() << Log::INFO << "Setting ... (PLEP) = " << value << endl;
00167       FC_INPU.PLEPIN = asDouble(value);
00168     }
00169     else if (name == "PBE1") {
00170       getLog() << Log::INFO << "Setting ... (PBE1) = " << value << endl;
00171       //FC_BEAM.PBEAM[0] = -1*asDouble(value);
00172       FC_BEAM.PBE1 = -1*asDouble(value);
00173     }
00174     else if (name == "PBE2") {
00175       getLog() << Log::INFO << "Setting ... (PBE2) = " << value << endl;
00176       FC_BEAM.PBEAM[1] = asDouble(value);
00177     }
00178     else if (name == "PT2CUT") {
00179       getLog() << Log::INFO << "Setting ... (PT2CUT) = " << value << endl;
00180       FC_PTCUT.PT2CUT[1] = asDouble(value);
00181     }
00182 
00184     else if (name == "KBE1") {
00185       getLog() << Log::INFO << "Setting ... (KBE1) = " << value << endl;
00186       FC_BEAM.KBEAM[0] = -1*asInt(value);
00187     }
00188     else if (name == "KBE2") {
00189       getLog() << Log::INFO << "Setting ... (KBE2) = " << value << endl;
00190       FC_BEAM.KBEAM[1] = asInt(value);
00191     }
00192     else if (name == "HADI") {
00193       getLog() << Log::INFO << "Setting ... (HADI) = " << value << endl;
00194       FC_INPU.HADI = asInt(value);
00195     } else if (name == "LEPI") {
00196       getLog() << Log::INFO << "Setting ... (LEPI) = " << value << endl;
00197       FC_INPU.LEPI = asInt(value);
00198     } else if (name == "IPRO") {
00199       getLog() << Log::INFO << "Setting ... (IPRO) = " << value << endl;
00200       FC_RAPA.IPRO = asInt(value);
00201     } else if (name == "IDIR") {
00202       getLog() << Log::INFO << "Setting ... (IDIR) = " << value << endl;
00203       FC_DISDIF.IDIR = asInt(value);
00204     } else if (name == "NFLAV") {
00205       getLog() << Log::INFO << "Setting ... (NFLAV) = " << value << endl;
00206       FC_LUCO.NFLAV = asInt(value);
00207     } else if (name == "NFLQCDC") {
00208       getLog() << Log::INFO << "Setting ... (NFLQCDC) = " << value << endl;
00209       FC_LUCO.NFLQCDC = asInt(value);
00210     } else if (name == "INTER") {
00211       getLog() << Log::INFO << "Setting ... (INTER) = " << value << endl;
00212       FC_INPU.INTER = asInt(value);
00213     }
00214 
00215     else {
00216       getLog() << Log::ERROR << "Rapgap doesn't have a parameter called " << name << endl;
00217       //return FAILURE;
00218     }
00219     return SUCCESS;
00220   }
00221 
00222   // Tidy up after ourselves
00223   void Rapgap::finalize() {
00224     getLog() << Log::INFO << "Finalising..." << endl;
00225     // Print out stats from run
00226     call_raend(20);
00227   }
00228 
00229 }
00230 
00231 extern "C" {
00232   AGILe::Generator* create() { return new AGILe::Rapgap(); }
00233   void destroy(AGILe::Generator* gen) { delete gen; }
00234 }
Generated on Tue Mar 6 10:39:39 2012 for AGILe - A Generator Interface Library (+ executable) by  doxygen 1.6.3