diff --git a/Cards/pythia8.cfi b/Cards/pythia8.cfi index f6c7267..c70632a 100644 --- a/Cards/pythia8.cfi +++ b/Cards/pythia8.cfi @@ -1,4 +1,28 @@ hadroniser = { name = "pythia8"; - decayWto = ( 11, 13 ); + seed = 0L; + pythiaConfiguration = ( + "Random:setSeed = on", // random engine + "Tune:preferLHAPDF = 2", + "Main:timesAllowErrors = 10000", + "Check:epTolErr = 0.01", + "Beams:setProductionScalesFromLHEF = off", + "SLHA:keepSM = on", + "SLHA:minMassSM = 1000.", + "ParticleDecays:limitTau0 = on", + "ParticleDecays:tau0Max = 10", + "ParticleDecays:allowPhotonRadiation = on", + //"ParticleDecays:allowPhotonRadiation = off", //FIXME + // CUEP8M1 tuning + "Tune:pp 14", + "Tune:ee 7", + "MultipartonInteractions:pT0Ref = 2.4024", + "MultipartonInteractions:ecmPow = 0.25208", + "MultipartonInteractions:expPow = 1.6", + // process-specific + "ProcessLevel:all = off", + "ProcessLevel:resonanceDecays = on", + "24:onMode = off", // disable all W decays, but... + "24:onIfAny = 11 13" // enable e-nue + mu-numu final states + ); }; diff --git a/CepGen/Cards/ConfigHandler.cpp b/CepGen/Cards/ConfigHandler.cpp index 20cdf6a..3d704bf 100644 --- a/CepGen/Cards/ConfigHandler.cpp +++ b/CepGen/Cards/ConfigHandler.cpp @@ -1,319 +1,304 @@ #include "ConfigHandler.h" #include "CepGen/Core/Exception.h" #include "CepGen/Processes/GamGamLL.h" #include "CepGen/Processes/PPtoLL.h" #include "CepGen/Processes/PPtoWW.h" #include "CepGen/Hadronisers/Pythia8Hadroniser.h" #ifdef LIBCONFIG namespace CepGen { namespace Cards { //----- specialization for CepGen input cards ConfigHandler::ConfigHandler( const char* file ) { #ifdef LIBCONFIG libconfig::Config cfg; try { cfg.readFile( file ); } catch ( const libconfig::ParseException& pe ) { FatalError( Form( "Failed to parse the configuration card \"%s\".\n\tParser error: %s (at line %d)", file, pe.getError(), pe.getLine() ) ); } try { const libconfig::Setting& root = cfg.getRoot(); const libconfig::Setting& proc = root["process"]; //--- type of process to consider const char* proc_name = proc["name"]; const std::string str_proc_name = proc_name; if ( str_proc_name == "lpair" ) params_.setProcess( new Process::GamGamLL ); else if ( str_proc_name == "pptoll" ) params_.setProcess( new Process::PPtoLL ); else if ( str_proc_name == "pptoww" ) params_.setProcess( new Process::PPtoWW ); else FatalError( Form( "Unrecognised process: %s", proc_name ) ); //--- process mode int int_mode; const char* mode; if ( proc.lookupValue( "mode", int_mode ) ) { params_.kinematics.mode = (Kinematics::ProcessMode)int_mode; } else if ( proc.lookupValue( "mode", mode ) ) { const std::string str_mode = mode; if ( str_mode == "elastic/elastic" ) params_.kinematics.mode = Kinematics::ProcessMode::ElasticElastic; else if ( str_mode == "elastic/inelastic" ) params_.kinematics.mode = Kinematics::ProcessMode::ElasticInelastic; else if ( str_mode == "inelastic/elastic" ) params_.kinematics.mode = Kinematics::ProcessMode::InelasticElastic; else if ( str_mode == "inelastic/inelastic" ) params_.kinematics.mode = Kinematics::ProcessMode::InelasticInelastic; else FatalError( Form( "Unrecognised interaction mode: %s", mode ) ); } //--- process kinematics if ( proc.exists( "in_kinematics" ) ) parseIncomingKinematics( proc["in_kinematics"] ); if ( proc.exists( "out_kinematics" ) ) parseOutgoingKinematics( proc["out_kinematics"] ); //--- hadroniser parameters if ( root.exists( "hadroniser" ) ) parseHadroniser( root["hadroniser"] ); //--- generation parameters if ( root.exists( "integrator" ) ) parseIntegrator( root["integrator"] ); if ( root.exists( "vegas" ) ) parseIntegrator( root["vegas"] ); //FIXME for backward compatibility if ( root.exists( "generator" ) ) parseGenerator( root["generator"] ); //--- taming functions if ( proc.exists( "taming_functions" ) ) parseTamingFunctions( proc["taming_functions"] ); } catch ( const libconfig::SettingNotFoundException& nfe ) { FatalError( Form( "Failed to retrieve the field \"%s\".", nfe.getPath() ) ); } catch ( const libconfig::SettingTypeException& te ) { FatalError( Form( "Field \"%s\" has wrong type.", te.getPath() ) ); } #else InWarning( "libconfig++ is not present on this machine" ); #endif } #ifdef LIBCONFIG void ConfigHandler::parseIncomingKinematics( const libconfig::Setting& kin ) { try { if ( kin.exists( "beam1_pz" ) ) params_.kinematics.inp.first = (double)kin["beam1_pz"]; if ( kin.exists( "beam2_pz" ) ) params_.kinematics.inp.second = (double)kin["beam2_pz"]; if ( kin.exists( "structure_functions" ) ) { const char* sf = kin["structure_functions" ]; const std::string sf_str = sf; if ( sf_str == "electron" ) params_.kinematics.structure_functions = StructureFunctions::Electron; else if ( sf_str == "elastic proton" ) params_.kinematics.structure_functions = StructureFunctions::ElasticProton; else if ( sf_str == "Suri-Yennie" ) params_.kinematics.structure_functions = StructureFunctions::SuriYennie; else if ( sf_str == "Szczurek-Uleshchenko" ) params_.kinematics.structure_functions = StructureFunctions::SzczurekUleshchenko; else if ( sf_str == "Fiore" ) params_.kinematics.structure_functions = StructureFunctions::FioreBrasse; else if ( sf_str == "ALLM" ) params_.kinematics.structure_functions = StructureFunctions::ALLM97; else if ( sf_str == "ALLM;91" ) params_.kinematics.structure_functions = StructureFunctions::ALLM91; else if ( sf_str == "ALLM;97" ) params_.kinematics.structure_functions = StructureFunctions::ALLM97; //else if ( sf_str == "ALLM;HHT" ) params_.kinematics.structure_functions = StructureFunctions::ALLM_HHT; //else if ( sf_str == "ALLM;HHT-FT" ) params_.kinematics.structure_functions = StructureFunctions::ALLM_HHT_FT; else if ( sf_str == "ALLM;GD07p" ) params_.kinematics.structure_functions = StructureFunctions::GD07p; else if ( sf_str == "ALLM;GD11p" ) params_.kinematics.structure_functions = StructureFunctions::GD11p; else FatalError( Form( "Invalid structure functions mode: %s", sf ) ); } } catch ( const libconfig::SettingNotFoundException& nfe ) { FatalError( Form( "Failed to retrieve the field \"%s\".", nfe.getPath() ) ); } catch ( const libconfig::SettingTypeException& te ) { FatalError( Form( "Field \"%s\" has wrong type.", te.getPath() ) ); } } void ConfigHandler::parseOutgoingKinematics( const libconfig::Setting& kin ) { try { if ( kin.exists( "pair" ) ) { Particle::ParticleCode pair = (Particle::ParticleCode)(int)kin["pair"]; params_.kinematics.central_system = { pair, pair }; } if ( kin.exists( "min_pt" ) ) params_.kinematics.cuts.central[Cuts::pt_single].min() = (double)kin["min_pt"]; if ( kin.exists( "max_pt" ) ) params_.kinematics.cuts.central[Cuts::pt_single].max() = (double)kin["max_pt"]; if ( kin.exists( "min_ptdiff" ) ) params_.kinematics.cuts.central[Cuts::pt_diff].min() = (double)kin["min_ptdiff"]; if ( kin.exists( "max_ptdiff" ) ) params_.kinematics.cuts.central[Cuts::pt_diff].max() = (double)kin["max_ptdiff"]; if ( kin.exists( "min_rapiditydiff" ) ) params_.kinematics.cuts.central[Cuts::rapidity_diff].min() = (double)kin["min_rapiditydiff"]; if ( kin.exists( "max_rapiditydiff" ) ) params_.kinematics.cuts.central[Cuts::rapidity_diff].max() = (double)kin["max_rapiditydiff"]; if ( kin.exists( "min_energy" ) ) params_.kinematics.cuts.central[Cuts::energy_single].min() = (double)kin["min_energy"]; if ( kin.exists( "max_energy" ) ) params_.kinematics.cuts.central[Cuts::energy_single].max() = (double)kin["max_energy"]; if ( kin.exists( "min_eta" ) ) params_.kinematics.cuts.central[Cuts::eta_single].min() = (double)kin["min_eta"]; if ( kin.exists( "max_eta" ) ) params_.kinematics.cuts.central[Cuts::eta_single].max() = (double)kin["max_eta"]; if ( kin.exists( "min_rapidity" ) ) params_.kinematics.cuts.central[Cuts::rapidity_single].min() = (double)kin["min_rapidity"]; if ( kin.exists( "max_rapidity" ) ) params_.kinematics.cuts.central[Cuts::rapidity_single].max() = (double)kin["max_rapidity"]; if ( kin.exists( "min_mx" ) ) params_.kinematics.cuts.remnants[Cuts::mass].min() = (double)kin["min_mx"]; if ( kin.exists( "max_mx" ) ) params_.kinematics.cuts.remnants[Cuts::mass].max() = (double)kin["max_mx"]; } catch ( const libconfig::SettingNotFoundException& nfe ) { FatalError( Form( "Failed to retrieve the field \"%s\".", nfe.getPath() ) ); } catch ( const libconfig::SettingTypeException& te ) { FatalError( Form( "Field \"%s\" has wrong type.", te.getPath() ) ); } } void ConfigHandler::parseIntegrator( const libconfig::Setting& integr ) { try { if ( integr.exists( "algorithm" ) ) { const std::string algo = integr["algorithm"]; if ( algo == "Vegas" ) params_.integrator.type = Integrator::Vegas; else if ( algo == "MISER" ) params_.integrator.type = Integrator::MISER; } if ( integr.exists( "num_points" ) ) params_.integrator.npoints = (int)integr["num_points"]; if ( integr.exists( "num_integration_calls" ) ) params_.integrator.ncvg = (int)integr["num_integration_calls"]; if ( integr.exists( "num_integration_iterations" ) ) params_.integrator.itvg = (int)integr["num_integration_iterations"]; if ( integr.exists( "seed" ) ) params_.integrator.seed = (unsigned long)integr["seed"]; } catch ( const libconfig::SettingNotFoundException& nfe ) { FatalError( Form( "Failed to retrieve the field \"%s\".", nfe.getPath() ) ); } } void ConfigHandler::parseGenerator( const libconfig::Setting& gen ) { params_.generation.enabled = true; try { if ( gen.exists( "num_events" ) ) params_.generation.maxgen = (int)gen["num_events"]; if ( gen.exists( "print_every" ) ) params_.generation.gen_print_every = (int)gen["print_every"]; } catch ( const libconfig::SettingNotFoundException& nfe ) { FatalError( Form( "Failed to retrieve the field \"%s\".", nfe.getPath() ) ); } } void ConfigHandler::parseTamingFunctions( const libconfig::Setting& tf ) { if ( !tf.isList() ) FatalError( "The taming functions definition must be wrapped within a list!" ); for ( unsigned short i = 0; i < tf.getLength(); ++i ) { params_.taming_functions.add( tf[i]["variable"], tf[i]["expression"] ); } } void ConfigHandler::parseHadroniser( const libconfig::Setting& hadr ) { try { std::string hadr_name = hadr["name"]; if ( hadr_name == "pythia8" ) { #ifdef PYTHIA8 Hadroniser::Pythia8Hadroniser* pythia8 = new Hadroniser::Pythia8Hadroniser; + long long seed = ( hadr.exists( "seed" ) ) ? hadr["seed"] : -1ll; + pythia8->setSeed( seed ); pythia8->readString( Form( "Beams:eCM = %.2f", params_.kinematics.sqrtS() ) ); - if ( hadr.exists( "decayWto" ) ) { - libconfig::Setting& pdg_ids = hadr["decayWto"]; - if ( !pdg_ids.isList() || pdg_ids.getLength() != 2 ) throw libconfig::SettingTypeException( pdg_ids ); - std::ostringstream os; - pythia8->readString( "24:onMode = off" ); // disable all W decays, but... - for ( unsigned short i = 0; i < pdg_ids.getLength(); ++i ) { - int pdg_id = pdg_ids[i]; - switch ( abs( pdg_id ) ) { - case 1: case 2: // u/d quark - os << " 1 2"; break; - case 3: case 4: // c/s quark - os << " 3 4"; break; - case 5: case 6: // t/b quark - os << " 5 6"; break; - case 11: // e/nu_e - os << " 11 12"; break; - case 13: // mu/nu_mu - os << " 13 14"; break; - case 15: // tau/nu_tau - os << " 15 16"; break; - default: { FatalError( Form( "Unsupported W decay mode: %d", pdg_id ) ); } break; - } + if ( hadr.exists( "pythiaConfiguration" ) ) { + libconfig::Setting& configs = hadr["pythiaConfiguration"]; + if ( !configs.isList() ) throw libconfig::SettingTypeException( configs ); + for ( unsigned short i = 0; i < configs.getLength(); ++i ) { + std::string config = configs[i]; + pythia8->readString( config ); } - pythia8->readString( Form( "24:onIfAny = %s", os.str().c_str() ) ); } pythia8->init(); params_.setHadroniser( pythia8 ); #endif } } catch ( const libconfig::SettingTypeException& nfe ) { FatalError( Form( "Wrong type to the field \"%s\".", nfe.getPath() ) ); } catch ( const libconfig::SettingNotFoundException& nfe ) { FatalError( Form( "Failed to retrieve the field \"%s\".", nfe.getPath() ) ); } } void ConfigHandler::writeProcess( const Parameters* params, libconfig::Setting& root ) { libconfig::Setting& proc = root.add( "process", libconfig::Setting::TypeGroup ); proc.add( "name", libconfig::Setting::TypeString ) = params->processName(); std::ostringstream os; os << params->kinematics.mode; proc.add( "mode", libconfig::Setting::TypeString ) = os.str(); } void ConfigHandler::writeIncomingKinematics( const Parameters* params, libconfig::Setting& root ) { libconfig::Setting& kin = root.add( "in_kinematics", libconfig::Setting::TypeGroup ); kin.add( "beam1_pz", libconfig::Setting::TypeFloat ) = params->kinematics.inp.first; kin.add( "beam2_pz", libconfig::Setting::TypeFloat ) = params->kinematics.inp.second; std::ostringstream os; os << params->kinematics.structure_functions; kin.add( "structure_function", libconfig::Setting::TypeString ) = os.str(); } void ConfigHandler::writeOutgoingKinematics( const Parameters* params, libconfig::Setting& root ) { libconfig::Setting& kin = root.add( "out_kinematics", libconfig::Setting::TypeGroup ); if ( params->kinematics.central_system.size() > 0 ) kin.add( "pair", libconfig::Setting::TypeInt ) = (int)params->kinematics.central_system[0]; if ( params->kinematics.cuts.central.count( Cuts::pt_single ) ) { kin.add( "min_pt", libconfig::Setting::TypeFloat ) = params->kinematics.cuts.central.at( Cuts::pt_single ).min(); kin.add( "max_pt", libconfig::Setting::TypeFloat ) = params->kinematics.cuts.central.at( Cuts::pt_single ).max(); } if ( params->kinematics.cuts.central.count( Cuts::pt_diff ) ) { kin.add( "min_ptdiff", libconfig::Setting::TypeFloat ) = params->kinematics.cuts.central.at( Cuts::pt_diff ).min(); kin.add( "max_ptdiff", libconfig::Setting::TypeFloat ) = params->kinematics.cuts.central.at( Cuts::pt_diff ).max(); } if ( params->kinematics.cuts.central.count( Cuts::rapidity_diff ) ) { kin.add( "min_dely", libconfig::Setting::TypeFloat ) = params->kinematics.cuts.central.at( Cuts::rapidity_diff ).min(); kin.add( "max_dely", libconfig::Setting::TypeFloat ) = params->kinematics.cuts.central.at( Cuts::rapidity_diff ).max(); } if ( params->kinematics.cuts.central.count( Cuts::energy_single ) ) { kin.add( "min_energy", libconfig::Setting::TypeFloat ) = params->kinematics.cuts.central.at( Cuts::energy_single ).min(); kin.add( "max_energy", libconfig::Setting::TypeFloat ) = params->kinematics.cuts.central.at( Cuts::energy_single ).max(); } if ( params->kinematics.cuts.central.count( Cuts::eta_single ) ) { kin.add( "min_eta", libconfig::Setting::TypeFloat ) = params->kinematics.cuts.central.at( Cuts::eta_single ).min(); kin.add( "max_eta", libconfig::Setting::TypeFloat ) = params->kinematics.cuts.central.at( Cuts::eta_single ).max(); } if ( params->kinematics.cuts.remnants.count( Cuts::mass ) ) { kin.add( "min_mx", libconfig::Setting::TypeFloat ) = params->kinematics.cuts.remnants.at( Cuts::mass ).min(); kin.add( "max_mx", libconfig::Setting::TypeFloat ) = params->kinematics.cuts.remnants.at( Cuts::mass ).max(); } } void ConfigHandler::writeTamingFunctions( const Parameters* params, libconfig::Setting& root ) { libconfig::Setting& tf = root.add( "taming_functions", libconfig::Setting::TypeList ); for ( std::map<std::string,TamingFunction>::const_iterator it = params->taming_functions.begin(); it != params->taming_functions.end(); ++it ) { libconfig::Setting& fun = tf.add( libconfig::Setting::TypeGroup ); fun.add( "variable", libconfig::Setting::TypeString ) = it->first; fun.add( "expression", libconfig::Setting::TypeString ) = it->second.expression; } } void ConfigHandler::writeIntegrator( const Parameters* params, libconfig::Setting& root ) { libconfig::Setting& integr = root.add( "integrator", libconfig::Setting::TypeGroup ); std::ostringstream os; os << params->integrator.type; integr.add( "algorithm", libconfig::Setting::TypeString ) = os.str(); integr.add( "num_points", libconfig::Setting::TypeInt ) = (int)params->integrator.npoints; integr.add( "num_integration_calls", libconfig::Setting::TypeInt ) = (int)params->integrator.ncvg; integr.add( "num_integration_iterations", libconfig::Setting::TypeInt ) = (int)params->integrator.itvg; integr.add( "seed", libconfig::Setting::TypeInt64 ) = (long)params->integrator.seed; } void ConfigHandler::writeGenerator( const Parameters* params, libconfig::Setting& root ) { if ( !params->generation.enabled ) return; libconfig::Setting& gen = root.add( "generator", libconfig::Setting::TypeGroup ); gen.add( "num_events", libconfig::Setting::TypeInt ) = (int)params->generation.maxgen; gen.add( "print_every", libconfig::Setting::TypeInt ) = (int)params->generation.gen_print_every; } #endif void ConfigHandler::store( const Parameters* params, const char* file ) { #ifdef LIBCONFIG libconfig::Config cfg; libconfig::Setting& root = cfg.getRoot(); writeProcess( params, root ); writeIncomingKinematics( params, root["process"] ); writeOutgoingKinematics( params, root["process"] ); writeTamingFunctions( params, root["process"] ); writeIntegrator( params, root ); writeGenerator( params, root ); cfg.writeFile( file ); #endif } } } #endif diff --git a/CepGen/Core/Integrand.cpp b/CepGen/Core/Integrand.cpp index 97a9f0f..b223c01 100644 --- a/CepGen/Core/Integrand.cpp +++ b/CepGen/Core/Integrand.cpp @@ -1,141 +1,143 @@ #include "CepGen/Core/Timer.h" #include "CepGen/Core/Exception.h" #include "CepGen/Event/Event.h" #include "CepGen/Event/Particle.h" #include "CepGen/Physics/Kinematics.h" #include "CepGen/Processes/GenericProcess.h" #include "CepGen/Parameters.h" #include <sstream> namespace CepGen { double f( double* x, size_t ndim, void* params ) { std::ostringstream os; Parameters* p = static_cast<Parameters*>( params ); std::shared_ptr<Event> ev = p->process()->event(); if ( p->process()->hasEvent() ) { p->process()->clearEvent(); + std::cout << "event cleared!!! " << ev->getByRole( Particle::CentralSystem ).size() << std::endl; + const Particle::Momentum p1( 0., 0., p->kinematics.inp.first ), p2( 0., 0., -p->kinematics.inp.second ); p->process()->setIncomingKinematics( p1, p2 ); // at some point introduce non head-on colliding beams? DebuggingInsideLoop( Form( "Function f called -- some parameters:\n\t" " pz(p1) = %5.2f pz(p2) = %5.2f\n\t" " remnant mode: %d", p->kinematics.inp.first, p->kinematics.inp.second, p->kinematics.structure_functions ) ); if ( p->integrator.first_run ) { if ( Logger::get().level >= Logger::Debug ) { std::ostringstream oss; oss << p->kinematics.mode; Debugging( Form( "Process mode considered: %s", oss.str().c_str() ) ); } Particle& op1 = ev->getOneByRole( Particle::OutgoingBeam1 ), &op2 = ev->getOneByRole( Particle::OutgoingBeam2 ); //--- add outgoing protons or remnants switch ( p->kinematics.mode ) { case Kinematics::ElectronProton: default: { InError( Form( "Process mode %d not yet handled!", (int)p->kinematics.mode ) ); } case Kinematics::ElasticElastic: break; // nothing to change in the event case Kinematics::ElasticInelastic: op2.setPdgId( Particle::uQuark, 1 ); break; case Kinematics::InelasticElastic: // set one of the outgoing protons to be fragmented op1.setPdgId( Particle::uQuark, 1 ); break; case Kinematics::InelasticInelastic: // set both the outgoing protons to be fragmented op1.setPdgId( Particle::uQuark, 1 ); op2.setPdgId( Particle::uQuark, 1 ); break; } //--- prepare the function to be integrated p->process()->prepareKinematics(); //--- add central system Particles& central_system = ev->getByRole( Particle::CentralSystem ); if ( central_system.size() == p->kinematics.central_system.size() ) { unsigned short i = 0; for ( Particles::iterator part = central_system.begin(); part != central_system.end(); ++part ) { if ( p->kinematics.central_system[i] == Particle::invalidParticle ) continue; part->setPdgId( p->kinematics.central_system[i] ); part->computeMass(); i++; } } p->process()->clearRun(); p->integrator.first_run = false; } } // event is not empty p->process()->setPoint( ndim, x ); if ( Logger::get().level >= Logger::DebugInsideLoop ) { std::ostringstream oss; for ( unsigned int i=0; i<ndim; i++ ) { oss << x[i] << " "; } DebuggingInsideLoop( Form( "Computing dim-%d point ( %s)", ndim, oss.str().c_str() ) ); } //--- from this step on, the phase space point is supposed to be set by 'GenericProcess::setPoint()' p->process()->beforeComputeWeight(); Timer tmr; // start the timer double integrand = p->process()->computeWeight(); if ( integrand < 0. ) return 0.; //--- only fill in the process' Event object if storage is requested // or if taming functions are to be applied p->process()->fillKinematics(); //--- once the kinematics variables have been populated, // can apply the collection of taming functions double taming = 1.0; if ( p->taming_functions.has( "m_central" ) || p->taming_functions.has( "pt_central" ) ) { const Particle::Momentum central_system( ev->getByRole( Particle::CentralSystem )[0].momentum() + ev->getByRole( Particle::CentralSystem )[1].momentum() ); taming *= p->taming_functions.eval( "m_central", central_system.mass() ); taming *= p->taming_functions.eval( "pt_central", central_system.pt() ); } if ( p->taming_functions.has( "q2" ) ) { taming *= p->taming_functions.eval( "q2", ev->getOneByRole( Particle::Parton1 ).momentum().mass() ); taming *= p->taming_functions.eval( "q2", ev->getOneByRole( Particle::Parton2 ).momentum().mass() ); } integrand *= taming; //--- full event content (+ hadronisation) if generating events if ( p->hadroniser() ) { p->hadroniser()->hadronise( ev.get() ); } if ( p->storage() ) { ev->time_generation = tmr.elapsed(); ev->time_total = tmr.elapsed(); p->process()->addGenerationTime( ev->time_total ); Debugging( Form( "Generation time: %5.6f sec\n\t" "Total time (gen+hadr): %5.6f sec", ev->time_generation, ev->time_total ) ); p->generation.last_event = ev; } // generating events if ( Logger::get().level>=Logger::DebugInsideLoop ) { os.str( "" ); for ( unsigned int i=0; i<ndim; i++ ) { os << Form( "%10.8f ", x[i] ); } Debugging( Form( "f value for dim-%d point ( %s): %4.4e", ndim, os.str().c_str(), integrand ) ); } return integrand; } } diff --git a/CepGen/Event/Event.cpp b/CepGen/Event/Event.cpp index f732c84..89c2da1 100644 --- a/CepGen/Event/Event.cpp +++ b/CepGen/Event/Event.cpp @@ -1,264 +1,294 @@ #include "Event.h" #include "CepGen/Core/Exception.h" #include "CepGen/Core/utils.h" #include <algorithm> #include <math.h> namespace CepGen { Event::Event() : num_hadronisation_trials( 0 ), time_generation( -1. ), time_total( -1. ) {} Event::~Event() {} Event& Event::operator=( const Event &ev_ ) { particles_ = ev_.particles_; time_generation = ev_.time_generation; time_total = ev_.time_total; num_hadronisation_trials = ev_.num_hadronisation_trials; return *this; } void Event::clear() { particles_.clear(); time_generation = -1.; time_total = -1.; } void - Event::init() + Event::freeze() { - last_particle_ = particles_.end(); + //--- store a snapshot of the primordial event block + if ( particles_.count( Particle::CentralSystem ) != 0 ) + last_cs_particle_ = particles_[Particle::CentralSystem].end(); + if ( particles_.count( Particle::OutgoingBeam1 ) != 0 ) + last_op1_particle_ = particles_[Particle::OutgoingBeam1].end(); + if ( particles_.count( Particle::OutgoingBeam2 ) != 0 ) + last_op2_particle_ = particles_[Particle::OutgoingBeam2].end(); } void Event::restore() { //--- remove all particles after the primordial event block - particles_.erase( last_particle_, particles_.end() ); + if ( particles_.count( Particle::CentralSystem ) != 0 ) + particles_[Particle::CentralSystem].erase( last_cs_particle_, particles_[Particle::CentralSystem].end() ); + if ( particles_.count( Particle::OutgoingBeam1 ) != 0 ) + particles_[Particle::OutgoingBeam1].erase( last_op1_particle_, particles_[Particle::OutgoingBeam1].end() ); + if ( particles_.count( Particle::OutgoingBeam2 ) != 0 ) + particles_[Particle::OutgoingBeam2].erase( last_op2_particle_, particles_[Particle::OutgoingBeam2].end() ); } Particles& Event::getByRole( const Particle::Role& role ) { //--- retrieve all particles with a given role return particles_[role]; } + ParticlesIds + Event::getIdsByRole( const Particle::Role& role ) const + { + //--- retrieve all particles ids with a given role + ParticlesIds out; + Particles parts; + try { + parts = particles_.at( role ); + } catch ( std::out_of_range ) { return out; } + for ( Particles::const_iterator it = parts.begin(); it != parts.end(); ++it ) { + out.insert( it->id() ); + } + return out; + } + Particle& Event::getOneByRole( const Particle::Role& role ) { //--- retrieve the first particle a the given role Particles& parts_by_role = getByRole( role ); if ( parts_by_role.size() > 1 ) { FatalError( Form( "More than one particle with role %d: %d particles", (int)role, parts_by_role.size() ) ); } return *parts_by_role.begin(); } Particle& Event::getById( int id ) { for ( ParticlesMap::iterator out=particles_.begin(); out!=particles_.end(); out++ ) { for ( Particles::iterator part=out->second.begin(); part!=out->second.end(); part++ ) { if ( part->id() == id ) return *part; } } throw Exception( __PRETTY_FUNCTION__, Form( "Failed to retrieve the particle with id=%d", id ), FatalError ); } const Particle& Event::getConstById( int id ) const { for ( ParticlesMap::const_iterator out=particles_.begin(); out!=particles_.end(); out++ ) { for ( Particles::const_iterator part=out->second.begin(); part!=out->second.end(); part++ ) { if ( part->id() == id ) return *part; } } throw Exception( __PRETTY_FUNCTION__, Form( "Failed to retrieve the particle with id=%d", id ), FatalError ); } Particles Event::getByIds( const ParticlesIds& ids ) const { Particles out; for ( ParticlesIds::const_iterator id=ids.begin(); id!=ids.end(); id++ ) { out.emplace_back( getConstById( *id ) ); } return out; } Particles Event::mothers( const Particle& part ) { return getByIds( part.mothers() ); } Particles Event::daughters( const Particle& part ) { return getByIds( part.daughters() ); } ParticleRoles Event::roles() const { ParticleRoles out; ParticlesMap::const_iterator it, end; for ( it=particles_.begin(), end=particles_.end(); it!=end; it=particles_.upper_bound( it->first ) ) { out.emplace_back( it->first ); } return out; } - void - Event::addParticle( Particle part, bool replace ) + Particle& + Event::addParticle( Particle& part, bool replace ) { DebuggingInsideLoop( Form( "Particle with PDGid = %d has role %d", part.pdgId(), part.role() ) ); if ( part.role() <= 0 ) FatalError( Form( "Trying to add a particle with role=%d", (int)part.role() ) ); //--- retrieve the list of particles with the same role Particles& part_with_same_role = getByRole( part.role() ); //--- specify the id if ( part_with_same_role.empty() && part.id() < 0 ) part.setId( numParticles() ); // set the id if previously invalid/inexistent if ( !part_with_same_role.empty() ) { if ( replace ) part.setId( part_with_same_role[0].id() ); // set the previous id if replacing a particle else part.setId( numParticles() ); } //--- add the particle to the collection if ( replace ) part_with_same_role = Particles( 1, part ); // generate a vector containing only this particle else part_with_same_role.emplace_back( part ); + + return part_with_same_role.back(); } - void + Particle& Event::addParticle( const Particle::Role& role, bool replace ) { - addParticle( Particle( role ), replace ); + Particle np( role ); + return addParticle( np, replace ); } size_t Event::numParticles() const { size_t out = 0; for ( ParticlesMap::const_iterator it=particles_.begin(); it!=particles_.end(); it++ ) { out += it->second.size(); } return out; } const Particles Event::particles() const { Particles out; for ( ParticlesMap::const_iterator it=particles_.begin(); it!=particles_.end(); it++ ) { out.insert( out.end(), it->second.begin(), it->second.end() ); } std::sort( out.begin(), out.end() ); return out; } const Particles Event::stableParticles() const { Particles out; for ( ParticlesMap::const_iterator it=particles_.begin(); it!=particles_.end(); it++ ) { for ( Particles::const_iterator part=it->second.begin(); part!=it->second.end(); part++ ) { if ( part->status() == Particle::Undefined || part->status() == Particle::FinalState ) { out.emplace_back( *part ); } } } std::sort( out.begin(), out.end() ); return out; } void Event::checkKinematics() const { try { const Particles& parts = particles(); // check the kinematics through parentage for ( Particles::const_iterator p = parts.begin(); p != parts.end(); ++p ) { ParticlesIds daughters = p->daughters(); if ( daughters.empty() ) continue; Particle::Momentum ptot; for ( ParticlesIds::const_iterator daugh = daughters.begin(); daugh != daughters.end(); ++daugh ) { const Particle& d = getConstById( *daugh ); const ParticlesIds mothers = d.mothers(); if ( mothers.size() > 1 ) { for ( ParticlesIds::const_iterator moth = mothers.begin(); moth != mothers.end(); ++moth ) { if ( *moth == p->id() ) continue; ptot -= getConstById( *moth ).momentum(); } } ptot += d.momentum(); } const double mass_diff = ( ptot-p->momentum() ).mass(); if ( fabs( mass_diff ) > 1.e-10 ) { dump(); throw Exception( __PRETTY_FUNCTION__, Form( "Error in momentum balance for particle %d: mdiff = %.5e", p->id(), mass_diff ), FatalError ); } } } catch ( const Exception& e ) { throw Exception( __PRETTY_FUNCTION__, Form( "Event kinematics check failed:\n\t%s", e.what() ), FatalError ); } } void Event::dump( std::ostream& out, bool stable ) const { Particles parts = ( stable ) ? stableParticles() : particles(); std::ostringstream os; double pxtot = 0., pytot = 0., pztot = 0., etot = 0.; for ( Particles::const_iterator part_ref=parts.begin(); part_ref!=parts.end(); part_ref++ ) { const Particle& part = *part_ref; { std::ostringstream oss; oss << part.pdgId(); os << Form( "\n %2d\t%+6d%8s", part.id(), part.integerPdgId(), oss.str().c_str() ); } os << "\t"; if ( part.charge() != 999. ) os << Form( "%6.2f ", part.charge() ); else os << "\t"; { std::ostringstream oss; oss << part.role(); os << Form( "%8s\t%6d\t", oss.str().c_str(), part.status() ); } const ParticlesIds mothers = part.mothers(); if ( !mothers.empty() ) { std::ostringstream oss; for ( ParticlesIds::const_iterator moth = mothers.begin(); moth != mothers.end(); ++moth ) { oss << ( ( moth != mothers.begin() ) ? "+" : "" ) << *moth; } os << Form( "%6s ", oss.str().c_str() ); } else os << " "; const Particle::Momentum mom = part.momentum(); - os << Form( "% 9.6e % 9.6e % 9.6e % 9.6e % 9.5e", mom.px(), mom.py(), mom.pz(), part.energy(), part.mass() ); + os << Form( "% 9.6e % 9.6e % 9.6e % 9.6e % 12.7f", mom.px(), mom.py(), mom.pz(), part.energy(), part.mass() ); if ( part.status() == Particle::Undefined || part.status() == Particle::Undecayed + || part.status() == Particle::Unfragmented || part.status() == Particle::FinalState ) { const int sign = ( part.status() == Particle::Undefined ) ? -1 : 1; pxtot += sign*mom.px(); pytot += sign*mom.py(); pztot += sign*mom.pz(); etot += sign*part.energy(); } } //--- set a threshold to the computation precision if ( fabs( pxtot ) < 1.e-10 ) pxtot = 0.; if ( fabs( pytot ) < 1.e-10 ) pytot = 0.; if ( fabs( pztot ) < 1.e-10 ) pztot = 0.; if ( fabs( etot ) < 1.e-10 ) etot = 0.; // Information( Form( "Dump of event content:\n" "Part.\tPDG id\t\tCharge\t Role\tStatus\tMother\t\t\t\t4-Momentum (GeV)\t\tMass (GeV)\n" "----\t------\t\t------\t ----\t------\t------\t------------------------------------------------------ -----------" "%s\n" "---------------------------------------------------------------------------------------------------------------------------\n" "\t\t\t\t\t\tTotal: % 9.6e % 9.6e % 9.6e % 9.6e", os.str().c_str(), pxtot, pytot, pztot, etot ) ); } } diff --git a/CepGen/Event/Event.h b/CepGen/Event/Event.h index 9fef87c..5ec67de 100644 --- a/CepGen/Event/Event.h +++ b/CepGen/Event/Event.h @@ -1,119 +1,120 @@ #ifndef CepGen_Event_Event_h #define CepGen_Event_Event_h #include "Particle.h" #include "CepGen/Core/Logger.h" namespace CepGen { /** * Class containing all the information on the in- and outgoing particles' kinematics * \brief Kinematic information on the particles in the event */ class Event { public: Event(); ~Event(); /** * \brief Copies all the relevant quantities from one Event object to another */ Event& operator=( const Event& ); /** * \brief Empties the whole event content */ void clear(); /// Initialize an "empty" event collection - void init(); + void freeze(); /// Restore the event to its "empty" state void restore(); /// Dump all the known information on every Particle object contained in this Event container in the output stream /// \param[out] os Output stream where to dump the information /// \param[in] stable_ Do we only show the stable particles in this event? void dump( std::ostream& os = Logger::get().outputStream, bool stable_=false ) const; //----- particles adders /// Set the information on one particle in the process /// \param[in] part The Particle object to insert or modify in the event /// \param[in] replace Do we replace the particle if already present in the event or do we append another particle with the same role ? - void addParticle( Particle part, bool replace=false ); + Particle& addParticle( Particle& part, bool replace=false ); /// \brief Create a new particle in the event, with no kinematic information but the role it has to play in the process /// \param[in] role The role the particle will play in the process /// \param[in] replace Do we replace the particle if already present in the event or do we append another particle with the same role ? - void addParticle( const Particle::Role& role, bool replace=false ); + Particle& addParticle( const Particle::Role& role, bool replace=false ); //----- particles retrievers /// Number of particles in the event size_t numParticles() const; /// \brief Vector of all particles in the event const Particles particles() const; /// \brief Vector of all stable particles in the event const Particles stableParticles() const; /** * Returns the list of Particle objects corresponding to a certain role in the process kinematics * \brief Gets a list of particles by their role in the event * \param[in] role The role the particles have to play in the process * \return A vector of references to the requested Particle objects */ Particles& getByRole( const Particle::Role& role ); + ParticlesIds getIdsByRole( const Particle::Role& role ) const; /** * Returns the first Particle object in the particles list whose role corresponds to the given argument * \param[in] role The role the particle has to play in the event * \return A Particle object corresponding to the first particle with the role */ Particle& getOneByRole( const Particle::Role& role ); /** * Returns the reference to the Particle object corresponding to a unique identifier in the event * \brief Gets one particle by its unique identifier in the event * \param[in] id_ The unique identifier to this particle in the event * \return A reference to the requested Particle object */ Particle& getById( int id_ ); /// Get a const Particle object using its unique identifier /// \param[in] id_ Unique identifier of the particle in the event /// \return Constant object to be retrieved const Particle& getConstById( int id_ ) const; /** * Returns the references to the Particle objects corresponding to the unique identifiers in the event * \brief Gets a vector of particles by their unique identifier in the event * \param[in] ids_ The unique identifiers to the particles to be selected in the event * \return A vector of references to the requested Particle objects */ Particles getByIds( const ParticlesIds& ids_ ) const; /** * Returns the list of mother particles of any given Particle object in this event * \param[in] part The reference to the Particle object from which we want to extract the mother particles * \return A list of parenting Particle object */ //----- general particles information retriever Particles mothers( const Particle& part ); /// Get a vector containing all the daughters from a particle /// \param[in] part The particle for which the daughter particles have to be retrieved /// \return Vector of Particle objects containing all the daughters' kinematic information Particles daughters( const Particle& part ); /// Get a list of roles for the given event (really process-dependant for the central system) /// \return Vector of integers corresponding to all the roles the particles can play in the event ParticleRoles roles() const; /// Number of trials before the event was "correctly" hadronised int num_hadronisation_trials; /// Time needed to generate the event at parton level (in seconds) float time_generation; /// Time needed to generate the hadronised (if needed) event (in seconds) float time_total; private: /// Check if the event kinematics is properly defined void checkKinematics() const; /// List of particles in the event, mapped to their role in the process ParticlesMap particles_; /// Last particle in an "empty" event - ParticlesMap::iterator last_particle_; + Particles::iterator last_cs_particle_, last_op1_particle_, last_op2_particle_; }; } #endif diff --git a/CepGen/Event/Momentum.cpp b/CepGen/Event/Momentum.cpp index 235d6a6..ce55b33 100644 --- a/CepGen/Event/Momentum.cpp +++ b/CepGen/Event/Momentum.cpp @@ -1,315 +1,321 @@ #include "Particle.h" #include "CepGen/Core/Exception.h" #include "CepGen/Core/utils.h" #include <math.h> namespace CepGen { using Momentum = Particle::Momentum; //----- Particle momentum methods Momentum::Momentum() : px_( 0. ), py_( 0. ), pz_( 0. ), p_( 0. ), energy_( 0. ) {} Momentum::Momentum( double x, double y, double z, double t ) : px_( x ), py_( y ), pz_( z ), energy_( t ) { computeP(); } //--- static constructors Momentum Momentum::fromPtEtaPhi( double pt, double eta, double phi, double e ) { const double px = pt*cos( phi ), py = pt*sin( phi ), pz = pt*sinh( eta ); return Momentum( px, py, pz, e ); } Momentum Momentum::fromPThetaPhi( double p, double theta, double phi, double e ) { const double px = p*sin( theta )*cos( phi ), py = p*sin( theta )*sin( phi ), pz = p*cos( theta ); return Momentum( px, py, pz, e ); } Momentum Momentum::fromPxPyPzE( double px, double py, double pz, double e ) { return Momentum( px, py, pz, e ); } //--- arithmetic operators Momentum& Momentum::operator+=( const Momentum& mom ) { px_ += mom.px_; py_ += mom.py_; pz_ += mom.pz_; energy_ += mom.energy_; //FIXME not supposed to be this way! computeP(); return *this; } Momentum& Momentum::operator-=( const Momentum& mom ) { px_ -= mom.px_; py_ -= mom.py_; pz_ -= mom.pz_; energy_ -= mom.energy_; //FIXME not supposed to be this way! computeP(); return *this; } + bool + Momentum::operator==( const Momentum& mom ) const + { + return ( px_ == mom.px_ && py_ == mom.py_ && pz_ == mom.pz_ && energy_ == mom.energy_ ); + } + double Momentum::threeProduct( const Momentum& mom ) const { DebuggingInsideLoop( Form( " (%f, %f, %f, %f)\n\t* (%f, %f, %f, %f)\n\t= %f", px_, py_, pz_, energy_, mom.px_, mom.py_, mom.pz_, mom.energy_, px_*mom.px_+py_*mom.py_+pz_*mom.pz_ ) ); return px_*mom.px_+py_*mom.py_+pz_*mom.pz_; } double Momentum::fourProduct( const Momentum& mom ) const { DebuggingInsideLoop( Form( " (%f, %f, %f, %f)\n\t* (%f, %f, %f, %f)\n\t= %f", px_, py_, pz_, energy_, mom.px_, mom.py_, mom.pz_, mom.energy_, px_*mom.px_+py_*mom.py_+pz_*mom.pz_ ) ); return energy_*mom.energy_-threeProduct(mom); } double Momentum::operator*=( const Momentum& mom ) { return threeProduct( mom ); } Momentum& Momentum::operator*=( double c ) { px_ *= c; py_ *= c; pz_ *= c; computeP(); return *this; } Momentum operator*( const Momentum& mom, double c ) { Momentum out = mom; out *= c; return out; } Momentum operator*( double c, const Momentum& mom ) { Momentum out = mom; out *= c; return out; } Momentum operator+( const Momentum& mom1, const Momentum& mom2 ) { Momentum out = mom1; out += mom2; return out; } Momentum operator-( const Momentum& mom1, const Momentum& mom2 ) { Momentum out = mom1; out -= mom2; return out; } double operator*( const Momentum& mom1, const Momentum& mom2 ) { Momentum tmp = mom1; return tmp.threeProduct( mom2 ); } //--- various setters void Momentum::setMass2( double m2 ) { energy_ = sqrt( p2()+m2 ); } void Momentum::setP( double px, double py, double pz, double e ) { setP( px, py, pz ); setEnergy( e ); } void Momentum::setP( double px, double py, double pz ) { px_ = px; py_ = py; pz_ = pz; computeP(); } void Momentum::setP( unsigned int i, double p ) { switch ( i ) { case 0: px_ = p; break; case 1: py_ = p; break; case 2: pz_ = p; break; case 3: energy_ = p; break; default: return; } computeP(); } void Momentum::computeP() { p_ = sqrt( px_*px_ + py_*py_ + pz_*pz_ ); } //--- various getters double Momentum::operator[]( const unsigned int i ) const { switch ( i ) { case 0: return px_; case 1: return py_; case 2: return pz_; case 3: return energy_; default: return -1.; } } const std::vector<double> Momentum::pVector() const { return std::vector<double>( { px(), py(), pz(), energy(), mass() } ); } double Momentum::mass() const { if ( mass2() >= 0. ) return sqrt( mass2() ); return -sqrt( -mass2() ); } double Momentum::theta() const { return atan2( pt(), pz() ); } double Momentum::phi() const { return atan2( py(), px() ); } double Momentum::pt() const { return sqrt( pt2() ); } double Momentum::eta() const { const int sign = ( pz()/fabs( pz() ) ); return ( pt() != 0. ) ? log( ( p()+fabs( pz() ) )/pt() )*sign : 9999.*sign; } double Momentum::rapidity() const { const int sign = ( pz()/fabs( pz() ) ); return ( energy() >= 0. ) ? log( ( energy()+pz() )/( energy()-pz() ) )*0.5 : 999.*sign; } //--- boosts/rotations Momentum& Momentum::betaGammaBoost( double gamma, double betagamma ) { if ( gamma == 1. && betagamma == 0. ) return *this; // trivial case const double pz = pz_, e = energy_; pz_ = gamma*pz+betagamma*e; energy_ = gamma*e +betagamma*pz; computeP(); return *this; } Momentum& Momentum::lorentzBoost( const Momentum& p ) { const double m = p.mass(); if ( m == p.energy() ) return *this; const double pf4 = ( ( *this )*p ) / m, fn = ( pf4+energy() )/( p.energy()+m ); *this -= p*fn; setEnergy( pf4 ); return *this; } Momentum& Momentum::rotatePhi( double phi, double sign ) { const double px = px_*cos( phi )+py_*sin( phi )*sign, py =-px_*sin( phi )+py_*cos( phi )*sign; px_ = px; py_ = py; return *this; } Momentum& Momentum::rotateThetaPhi( double theta, double phi ) { double rotmtx[3][3], mom[3]; //FIXME check this! cos(phi)->-sin(phi) & sin(phi)->cos(phi) --> phi->phi+pi/2 ? rotmtx[0][0] = -sin( phi ); rotmtx[0][1] = -cos( theta )*cos( phi ); rotmtx[0][2] = sin( theta )*cos( phi ); rotmtx[1][0] = cos( phi ); rotmtx[1][1] = -cos( theta )*sin( phi ); rotmtx[1][2] = sin( theta )*sin( phi ); rotmtx[2][0] = 0.; rotmtx[2][1] = sin( theta ); rotmtx[2][2] = cos( theta ); for (int i=0; i<3; i++) { mom[i] = 0.; for (int j=0; j<3; j++) { mom[i] += rotmtx[i][j]*operator[]( j ); } } setP( mom[0], mom[1], mom[2] ); return *this; } //--- printout std::ostream& operator<<( std::ostream& os, const Momentum& mom ) { return os << "(E ; p) = (" << mom.energy_ << " ; " << mom.px_ << ", " << mom.py_ << ", " << mom.pz_ << ")"; } } diff --git a/CepGen/Event/Particle.h b/CepGen/Event/Particle.h index b561a3e..724948a 100644 --- a/CepGen/Event/Particle.h +++ b/CepGen/Event/Particle.h @@ -1,402 +1,404 @@ #ifndef CepGen_Event_Particle_h #define CepGen_Event_Particle_h #include <set> #include <map> #include <vector> namespace CepGen { /// A set of integer-type particle identifiers typedef std::set<int> ParticlesIds; /// Kinematic information for one particle class Particle { public: /** Unique identifier for a particle type. From \cite Beringer:1900zz : * `The Monte Carlo particle numbering scheme [...] is intended to facilitate interfacing between event generators, detector simulators, and analysis packages used in particle physics.` * \brief PDG ids of all known particles */ enum ParticleCode { invalidParticle = 0, //--- fundamental particles dQuark = 1, uQuark = 2, Electron = 11, ElectronNeutrino = 12, Muon = 13, MuonNeutrino = 14, Tau = 15, TauNeutrino = 16, Gluon = 21, Photon = 22, Z = 23, W = 24, //--- composite particles PiPlus = 211, PiZero = 111, Rho770_0 = 113, Rho1450_0 = 100113, Rho1700_0 = 30113, Eta = 221, Omega782 = 223, h1380_1 = 10333, JPsi= 443, Phi1680 = 100333, Upsilon1S = 553, Upsilon2S = 100553, Upsilon3S = 200553, ud0Diquark = 2101, ud1Diquark = 2103, uu1Diquark = 2203, Proton = 2212, Neutron = 2112, Pomeron = 990, Reggeon = 110 }; /// Internal status code for a particle enum Status { PrimordialIncoming = -9, - Undecayed = -3, + Unfragmented = -4, Undecayed = -3, Propagator = -2, Incoming = -1, Undefined = 0, FinalState = 1, Resonance = 2, DebugResonance = 3, PythiaHIncoming = 21 }; /// Role of the particle in the process enum Role { UnknownRole = -1, IncomingBeam1 = 1, IncomingBeam2 = 2, OutgoingBeam1 = 3, OutgoingBeam2 = 5, CentralSystem = 6, Intermediate = 4, Parton1 = 41, Parton2 = 42, Parton3 = 43 }; /** * Container for a particle's 4-momentum, along with useful methods to ease the development of any matrix element level generator * \brief 4-momentum for a particle * \date Dec 2015 * \author Laurent Forthomme <laurent.forthomme@cern.ch> */ class Momentum { public: /// Build a 4-momentum at rest with an invalid energy (no mass information known) Momentum(); /// Build a 4-momentum using its 3-momentum coordinates and its energy Momentum( double x_, double y_, double z_, double t_=-1. ); inline ~Momentum() {} // --- static definitions /// Build a 3-momentum from its three pseudo-cylindric coordinates static Momentum fromPtEtaPhi( double pt, double eta, double phi, double e=-1. ); /// Build a 4-momentum from its scalar momentum, and its polar and azimuthal angles static Momentum fromPThetaPhi( double p, double theta, double phi, double e=-1. ); /// Build a 4-momentum from its four momentum and energy coordinates static Momentum fromPxPyPzE( double px, double py, double pz, double e ); // --- vector and scalar operators /// Scalar product of the 3-momentum with another 3-momentum double threeProduct( const Momentum& ) const; /// Scalar product of the 4-momentum with another 4-momentum double fourProduct( const Momentum& ) const; /// Add a 4-momentum through a 4-vector sum Momentum& operator+=( const Momentum& ); /// Subtract a 4-momentum through a 4-vector sum Momentum& operator-=( const Momentum& ); /// Scalar product of the 3-momentum with another 3-momentum double operator*=( const Momentum& ); /// Multiply all 4-momentum coordinates by a scalar Momentum& operator*=( double c ); + /// Equality operator + bool operator==( const Momentum& ) const; /// Human-readable format for a particle's momentum friend std::ostream& operator<<( std::ostream& os, const Particle::Momentum& mom ); Momentum& betaGammaBoost( double gamma, double betagamma ); /// Forward Lorentz boost Momentum& lorentzBoost( const Particle::Momentum& p ); // --- setters and getters /// Set all the components of the 4-momentum (in GeV) void setP( double px, double py, double pz, double e ); /// Set all the components of the 3-momentum (in GeV) void setP( double px, double py, double pz ); /// Set an individual component of the 4-momentum (in GeV) void setP( unsigned int i, double p ); /// Set the energy (in GeV) inline void setEnergy( double e ) { energy_ = e; } /// Compute the energy from the mass inline void setMass( double m ) { setMass2( m*m ); } /// Compute the energy from the mass void setMass2( double m2 ); /// Get one component of the 4-momentum (in GeV) double operator[]( const unsigned int i ) const; /// Momentum along the \f$x\f$-axis (in GeV) inline double px() const { return px_; } /// Momentum along the \f$y\f$-axis (in GeV) inline double py() const { return py_; } /// Longitudinal momentum (in GeV) inline double pz() const { return pz_; } /// Transverse momentum (in GeV) double pt() const; /// Squared transverse momentum (in GeV\f$^\textrm{2}\f$) inline double pt2() const { return ( px()*px()+py()*py() ); } /// 4-vector of double precision floats (in GeV) const std::vector<double> pVector() const; /// 3-momentum norm (in GeV) inline double p() const { return p_; } /// Squared 3-momentum norm (in GeV\f$^\textrm{2}\f$) inline double p2() const { return p_*p_; } /// Energy (in GeV) inline double energy() const { return energy_; } /// Squared energy (in GeV^2) inline double energy2() const { return energy_*energy_; } /// Squared mass (in GeV^2) as computed from its energy and momentum inline double mass2() const { return energy2()-p2(); } /// Mass (in GeV) as computed from its energy and momentum /// \note Returns \f$-\sqrt{|E^2-\mathbf{p}^2|}<0\f$ if \f$\mathbf{p}^2>E^2\f$ double mass() const; /// Polar angle (angle with respect to the longitudinal direction) double theta() const; /// Azimutal angle (angle in the transverse plane) double phi() const; /// Pseudo-rapidity double eta() const; /// Rapidity double rapidity() const; /// Rotate the transverse components by an angle phi (and reflect the y coordinate) Momentum& rotatePhi( double phi, double sign ); /// Rotate the particle's momentum by a polar/azimuthal angle Momentum& rotateThetaPhi( double theta_, double phi_ ); /// Apply a \f$ z\rightarrow -z\f$ transformation inline Momentum& mirrorZ() { pz_ = -pz_; return *this; } private: /// Compute the 3-momentum's norm void computeP(); /// Momentum along the \f$x\f$-axis double px_; /// Momentum along the \f$y\f$-axis double py_; /// Momentum along the \f$z\f$-axis double pz_; /// 3-momentum's norm (in GeV/c) double p_; /// Energy (in GeV) double energy_; }; /// Human-readable format for a particle's PDG code friend std::ostream& operator<<( std::ostream& os, const Particle::ParticleCode& pc ); /// Human-readable format for a particle's role in the event friend std::ostream& operator<<( std::ostream& os, const Particle::Role& rl ); /// Compute the 4-vector sum of two 4-momenta friend Particle::Momentum operator+( const Particle::Momentum& mom1, const Particle::Momentum& mom2 ); /// Compute the 4-vector difference of two 4-momenta friend Particle::Momentum operator-( const Particle::Momentum& mom1, const Particle::Momentum& mom2 ); /// Scalar product of two 3-momenta friend double operator*( const Particle::Momentum& mom1, const Particle::Momentum& mom2 ); /// Multiply all components of a 4-momentum by a scalar friend Particle::Momentum operator*( const Particle::Momentum& mom, double c ); /// Multiply all components of a 4-momentum by a scalar friend Particle::Momentum operator*( double c, const Particle::Momentum& mom ); //----- static getters /** * \brief Gets the mass of a particle * \param pdgId ParticleCode (PDG ID) * \return Mass of the particle in \f$\textrm{GeV}/c^2\f$ */ static double massFromPDGId( const Particle::ParticleCode& pdgId ); /** * \brief Gets the electric charge of a particle * \param pdgId ParticleCode (PDG ID) * \return Charge of the particle in \f$e\f$ */ static double chargeFromPDGId( const Particle::ParticleCode& pdgId ); /** * \brief Total decay width of one unstable particle * \param[in] pdgId ParticleCode (PDG ID) * \return Decay width in GeV */ static double widthFromPDGId( const Particle::ParticleCode& pdgId ); /// Convert a polar angle to a pseudo-rapidity static double thetaToEta( double theta ); /// Convert a pseudo-rapidity to a polar angle static double etaToTheta( double eta ); /// Convert a pseudo-rapidity to a rapidity static double etaToY( double eta_, double m_, double pt_ ); Particle(); /// Build using the role of the particle in the process and its PDG id /// \param[in] pdgId ParticleCode (PDG ID) /// \param[in] role Role of the particle in the process /// \param[in] st Current status Particle( Role role, ParticleCode pdgId=invalidParticle, Status st=Undefined ); /// Copy constructor Particle( const Particle& ); inline ~Particle() {} /// Comparison operator (from unique identifier) bool operator<( const Particle& rhs ) const; /// Comparison operator (from their reference's unique identifier) //bool operator<( Particle *rhs ) const { return ( id < rhs->id ); } Particle& lorentzBoost( double m_, const Momentum& mom_ ); /// Lorentz boost (shamelessly stolen from ROOT) std::vector<double> lorentzBoost( const Momentum& mom_ ); // --- general particle properties /// Unique identifier (in a Event object context) int id() const { return id_; } //void setId( int id ) { id_ = id; } /// Set the particle unique identifier in an event void setId( int id ) { id_ = id; } /// Electric charge (given as a float number, for the quarks and bound states) float charge() const { return charge_sign_ * chargeFromPDGId( pdg_id_ ); } /// Set the electric charge sign (+-1 for charged or 0 for neutral particles) void setChargeSign( int sign ) { charge_sign_ = sign; } /// Role in the considered process Role role() const { return role_; } /// Set the particle role in the process void setRole( const Role& role ) { role_ = role; } /** * Codes 1-10 correspond to currently existing partons/particles, and larger codes contain partons/particles which no longer exist, or other kinds of event information * \brief Particle status */ Status status() const { return status_; } /// Set the particle decay/stability status - void setStatus( const Status& status ) { status_ = status; } + void setStatus( Status status ) { status_ = status; } /// Set the PDG identifier (along with the particle's electric charge) /// \param[in] pdg ParticleCode (PDG ID) /// \param[in] ch Electric charge (0, 1, or -1) void setPdgId( const ParticleCode& pdg, short ch=0 ); /// Retrieve the objectified PDG identifier inline ParticleCode pdgId() const { return pdg_id_; } /// Retrieve the integer value of the PDG identifier int integerPdgId() const; /// Particle's helicity float helicity() const { return helicity_; } /// Set the helicity of the particle void setHelicity( float heli ) { helicity_ = heli; } /** * Gets the particle's mass in \f$\textrm{GeV}/c^{2}\f$. * \brief Gets the particle's mass * \return The particle's mass */ inline double mass() const { return mass_; }; /** * Set the mass of the particle in \f$\textrm{GeV}/c^{2}\f$ while ensuring that the kinematics is properly set (the mass is set according to the energy and the momentum in priority) * \brief Compute the particle's mass in \f$\textrm{GeV}/c^{2}\f$ */ void computeMass( bool off_shell=false ); /** * Set the mass of the particle in \f$\textrm{GeV}/c^{2}\f$ according to a value given as an argument. This method ensures that the kinematics is properly set (the mass is set according to the energy and the momentum in priority) * \param m The mass in \f$\textrm{GeV}/c^{2}\f$ to set * \brief Set the particle's mass in \f$\textrm{GeV}/c^{2}\f$ */ void setMass( double m=-1. ); /// Get the particle's squared mass (in \f$\textrm{GeV}^\textrm{2}\f$) inline double mass2() const { return mass_*mass_; }; /// Retrieve the momentum object associated with this particle inline Momentum& momentum() { return momentum_; } /// Retrieve the momentum object associated with this particle inline Momentum momentum() const { return momentum_; } /// Associate a momentum object to this particle void setMomentum( const Momentum& mom, bool offshell=false ); /** * \brief Set the 3-momentum associated to the particle * \param[in] px Momentum along the \f$x\f$-axis, in \f$\textrm{GeV}/c\f$ * \param[in] py Momentum along the \f$y\f$-axis, in \f$\textrm{GeV}/c\f$ * \param[in] pz Momentum along the \f$z\f$-axis, in \f$\textrm{GeV}/c\f$ */ void setMomentum( double px, double py, double pz ); /** * \brief Set the 4-momentum associated to the particle * \param[in] px Momentum along the \f$x\f$-axis, in \f$\textrm{GeV}/c\f$ * \param[in] py Momentum along the \f$y\f$-axis, in \f$\textrm{GeV}/c\f$ * \param[in] pz Momentum along the \f$z\f$-axis, in \f$\textrm{GeV}/c\f$ * \param[in] e Energy, in GeV */ void setMomentum( double px, double py, double pz, double e ); /** * \brief Set the 4-momentum associated to the particle * \param[in] p 4-momentum */ inline void setMomentum( double p[4] ) { setMomentum( p[0], p[1], p[2], p[3] ); } /** * \brief Set the particle's energy * \param[in] e Energy, in GeV */ void setEnergy( double e=-1. ); /// Get the particle's energy (in GeV) double energy() const; /// Get the particle's squared energy (in \f$\textrm{GeV}^\textrm{2}\f$) inline double energy2() const { return energy()*energy(); }; /// Is this particle a valid particle which can be used for kinematic computations? bool valid(); // --- particle relations /// Is this particle a primary particle? inline bool primary() const { return mothers_.empty(); } /** * \brief Set the mother particle * \param[in] part A Particle object containing all the information on the mother particle */ void addMother( Particle& part ); /** * \brief Gets the unique identifier to the mother particle from which this particle arises * \return An integer representing the unique identifier to the mother of this particle in the event */ inline ParticlesIds mothers() const { return mothers_; } /** * \brief Add a decay product * \param[in] part The Particle object in which this particle will desintegrate or convert * \return A boolean stating if the particle has been added to the daughters list or if it was already present before */ void addDaughter( Particle& part ); /// Gets the number of daughter particles inline unsigned int numDaughters() const { return daughters_.size(); }; /** * \brief Get an identifiers list all daughter particles * \return An integer vector containing all the daughters' unique identifier in the event */ inline ParticlesIds daughters() const { return daughters_; } // --- global particle information extraction /// Dump all the information on this particle into the standard output stream void dump() const; private: /// Unique identifier in an event int id_; /// Electric charge (+-1 or 0) short charge_sign_; /// Momentum properties handler Momentum momentum_; /// Mass in \f$\textrm{GeV}/c^2\f$ double mass_; /// Helicity float helicity_; /// Role in the process Role role_; /// Decay/stability status Status status_; /// List of mother particles ParticlesIds mothers_; /// List of daughter particles ParticlesIds daughters_; /// PDG id ParticleCode pdg_id_; }; /// Compute the centre of mass energy of two particles (incoming or outgoing states) double CMEnergy( const Particle& p1, const Particle& p2 ); /// Compute the centre of mass energy of two particles (incoming or outgoing states) double CMEnergy( const Particle::Momentum& m1, const Particle::Momentum& m2 ); //bool operator<( const Particle& a, const Particle& b ) { return a.id<b.id; } // --- particle containers /// List of Particle objects typedef std::vector<Particle> Particles; /// List of particles' roles typedef std::vector<Particle::Role> ParticleRoles; /// Map between a particle's role and its associated Particle object typedef std::map<Particle::Role,Particles> ParticlesMap; } #endif diff --git a/CepGen/Hadronisers/Pythia8Hadroniser.cpp b/CepGen/Hadronisers/Pythia8Hadroniser.cpp index 1d1414f..0d06ff5 100644 --- a/CepGen/Hadronisers/Pythia8Hadroniser.cpp +++ b/CepGen/Hadronisers/Pythia8Hadroniser.cpp @@ -1,55 +1,116 @@ #ifdef PYTHIA8 #include "Pythia8Hadroniser.h" namespace CepGen { namespace Hadroniser { Pythia8Hadroniser::Pythia8Hadroniser(): GenericHadroniser( "pythia8" ), pythia_( std::unique_ptr<Pythia8::Pythia>( new Pythia8::Pythia ) ) {} Pythia8Hadroniser::~Pythia8Hadroniser() {} + void + Pythia8Hadroniser::setSeed( long long seed ) + { + if ( seed == -1ll ) { + pythia_->readString( "Random:setSeed = off" ); + return; + } + pythia_->readString( "Random:setSeed = on" ); + pythia_->readString( Form( "Random:seed = %llu", seed ) ); + } + bool Pythia8Hadroniser::hadronise( const Particle* part ) { return true; } bool Pythia8Hadroniser::hadronise( Event* ev ) { + //--- first start by cleaning up the previous runs leftovers + + pythia_->event.reset(); + std::cout << ">>>>> before hadronisation" << std::endl; + ev->dump(); + //--- outgoing proton remnants fragmentation Particles& op1 = ev->getByRole( Particle::OutgoingBeam1 ), &op2 = ev->getByRole( Particle::OutgoingBeam2 ); for ( Particles::const_iterator p_it = op1.begin(); p_it != op1.end(); ++p_it ) { //... excited proton fragmentation } for ( Particles::const_iterator p_it = op2.begin(); p_it != op2.end(); ++p_it ) { //... excited proton fragmentation } //--- central system hadronisation/decay/... - Particles& cs = ev->getByRole( Particle::CentralSystem ); - for ( Particles::const_iterator p_it = cs.begin(); p_it != cs.end(); ++p_it ) { - if ( p_it->status() != Particle::Undecayed ) continue; - p_it->dump(); + const ParticlesIds cs = ev->getIdsByRole( Particle::CentralSystem ); + for ( ParticlesIds::iterator p_it = cs.begin(); p_it != cs.end(); ++p_it ) { + const Particle& part = ev->getById( *p_it ); + + // check if the particle can be decayed + if ( part.status() != Particle::Undecayed ) continue; + if ( !pythia_->particleData.canDecay( part.pdgId() ) ) continue; + + // check if the particle has a valid kinematics + const Particle::Momentum mom = part.momentum(); + if ( mom == Particle::Momentum() ) continue; // FIXME probably during warm-up... + + pythia_->event.reset(); + Pythia8::Particle py8part( part.integerPdgId(), 93, 0, 0, 0, 0, 0, 0, mom.px(), mom.py(), mom.pz(), mom.energy(), part.mass() ); + py8part.tau( pythia_->particleData.tau0( part.pdgId() ) ); + pythia_->event.append( py8part ); + + const unsigned short num_before = pythia_->event.size(); // number of particles before the decay + pythia_->next(); // launch the decay + + // map { "pythia id" -> "cepgen id" } + std::map<unsigned short, unsigned short> py_cg_corresp = { { 1, part.id() } }; + + const unsigned short num_after = pythia_->event.size(); // number of particles after the decay + if ( num_before == num_after ) continue; // nothing happened... so no decay! + + for ( unsigned short i = 2; i < pythia_->event.size(); ++i ) { // skip the initial system + const Pythia8::Particle p = pythia_->event[i]; + const std::vector<int> mothers = p.motherList(); + if ( mothers.size() == 0 ) continue; + Particle& op = ev->addParticle( Particle::CentralSystem ); + py_cg_corresp[i] = op.id(); + +//std::cout << "dumping ids correspondence (" << *p_it << ")" << std::endl; for ( const auto& p : py_cg_corresp ) { std::cout << p.first << "--->" << p.second << std::endl; } std::cout << "-------enddump-------" << std::endl; + op.setPdgId( ( Particle::ParticleCode )abs( p.id() ), p.charge() ); + if ( p.isFinal() ) op.setStatus( Particle::FinalState ); + else op.setStatus( Particle::Propagator ); + + op.setMomentum( Particle::Momentum( p.px(), p.py(), p.pz(), p.e() ) ); + for ( std::vector<int>::const_iterator moth = mothers.begin(); moth != mothers.end(); ++moth ) { + if ( py_cg_corresp.count( *moth ) == 0 ) { FatalError( Form( "Particle with id=%d was not found in the event content!", *moth ) ); } + op.addMother( ev->getById( py_cg_corresp[*moth] ) ); + } + } + ev->getById( *p_it ).setStatus( Particle::Resonance ); + //std::cout << part.id() << "->" << part.integerPdgId() << "->" << part.status() << std::endl; } + ev->dump(); + return true; } bool Pythia8Hadroniser::prepareHadronisation( Event *ev ) { - pythia_->init(); + init(); return true; } } } #endif diff --git a/CepGen/Hadronisers/Pythia8Hadroniser.h b/CepGen/Hadronisers/Pythia8Hadroniser.h index ac6a22f..6f76797 100644 --- a/CepGen/Hadronisers/Pythia8Hadroniser.h +++ b/CepGen/Hadronisers/Pythia8Hadroniser.h @@ -1,36 +1,37 @@ #ifndef Pythia8Hadroniser_h #define Pythia8Hadroniser_h #include <Pythia8/Pythia.h> #include <memory> #include "GenericHadroniser.h" namespace CepGen { namespace Hadroniser { /** * Full interface to the Pythia8 hadronisation algorithm. It can be used in a single particle decay mode as well as a full event hadronisation using the string model, as in Jetset. * \brief Pythia8 hadronisation algorithm */ class Pythia8Hadroniser : public GenericHadroniser { public: Pythia8Hadroniser(); ~Pythia8Hadroniser(); bool hadronise( const Particle* part ) override; bool hadronise( Event* ev ) override; + void setSeed( long long seed ); void readString( const char* param ) { pythia_->readString( param ); } void readString( const std::string& param ) { pythia_->readString( param ); } void init() { pythia_->init(); } private: bool prepareHadronisation( Event *ev ); std::unique_ptr<Pythia8::Pythia> pythia_; }; } } #endif diff --git a/CepGen/Processes/GamGamLL.cpp b/CepGen/Processes/GamGamLL.cpp index dab5dfa..4d1fc92 100644 --- a/CepGen/Processes/GamGamLL.cpp +++ b/CepGen/Processes/GamGamLL.cpp @@ -1,877 +1,877 @@ #include "GamGamLL.h" #include "CepGen/Core/Exception.h" using namespace CepGen::Process; GamGamLL::GamGamLL( int nopt ) : GenericProcess( "lpair", "pp -> p(*) (gamma gamma -> l+ l-) p(*)" ), n_opt_( nopt ), MX2_( 0. ), MY2_( 0. ), Ml2_( 0. ), ep1_( 0. ), ep2_( 0. ), p_cm_( 0. ), w12_( 0. ), w31_( 0. ), dw31_( 0. ), w52_( 0. ), dw52_( 0. ), ec4_( 0. ), pc4_( 0. ), mc4_( 0. ), w4_( 0. ), p12_( 0. ), p1k2_( 0. ), p2k1_( 0. ), p13_( 0. ), p14_( 0. ), p25_( 0. ), q1dq_( 0. ), q1dq2_( 0. ), s1_( 0. ), s2_( 0. ), epsi_( 0. ), g5_( 0. ), g6_( 0. ), a5_( 0. ), a6_( 0. ), bb_( 0. ), gram_( 0. ), dd1_( 0. ), dd2_( 0. ), dd3_( 0. ), dd4_( 0. ), dd5_( 0. ), delta_( 0. ), g4_( 0. ), sa1_( 0. ), sa2_( 0. ), sl1_( 0. ), cos_theta4_( 0. ), sin_theta4_( 0. ), al4_( 0. ), be4_( 0. ), de3_( 0. ), de5_( 0. ), pt4_( 0. ), jacobian_( 0. ) {} void GamGamLL::addEventContent() { GenericProcess::setEventContent( { { Particle::IncomingBeam1, Particle::Proton }, { Particle::IncomingBeam2, Particle::Proton }, { Particle::Parton1, Particle::Photon }, { Particle::Parton2, Particle::Photon } }, { { Particle::OutgoingBeam1, { Particle::Proton } }, { Particle::OutgoingBeam2, { Particle::Proton } }, { Particle::CentralSystem, { Particle::Muon, Particle::Muon } } } ); } unsigned int GamGamLL::numDimensions( const Kinematics::ProcessMode& process_mode ) const { switch ( process_mode ) { case Kinematics::ElectronProton: { InError( "Not supported yet!" ); } case Kinematics::ElasticElastic: default: return 7; case Kinematics::ElasticInelastic: case Kinematics::InelasticElastic: return 8; case Kinematics::InelasticInelastic: return 9; } } bool GamGamLL::pickin() { DebuggingInsideLoop( Form( "Optimised mode? %i", n_opt_ ) ); jacobian_ = 0.; w4_ = mc4_*mc4_; // sig1 = sigma and sig2 = sigma' in [1] const double sig = mc4_+MY_; double sig1 = sig*sig, sig2 = sig1; DebuggingInsideLoop( Form( "mc4 = %f\n\t" "sig1 = %f\n\t" "sig2 = %f", mc4_, sig1, sig2 ) ); // Mass difference between the first outgoing particle and the first incoming // particle w31_ = MX2_-w1_; // Mass difference between the second outgoing particle and the second // incoming particle w52_ = MY2_-w2_; // Mass difference between the two incoming particles w12_ = w1_-w2_; // Mass difference between the central two-photons system and the second // outgoing particle const double d6 = w4_-MY2_; DebuggingInsideLoop( Form( "w1 = %f\n\t" "w2 = %f\n\t" "w3 = %f\n\t" "w4 = %f\n\t" "w5 = %f", w1_, w2_, MX2_, w4_, MY2_ ) ); DebuggingInsideLoop( Form( "w31 = %f\n\tw52 = %f\n\tw12 = %f", w31_, w52_, w12_ ) ); const double ss = s_+w12_; const double rl1 = ss*ss-4.*w1_*s_; // lambda(s, m1**2, m2**2) if ( rl1 <= 0. ) { InWarning( Form( "rl1 = %f <= 0", rl1 ) ); return false; } sl1_ = sqrt( rl1 ); s2_ = 0.; double ds2 = 0.; if ( n_opt_ == 0 ) { const double smax = s_+MX2_-2.*MX_*sqs_; Map( x(2), sig1, smax, s2_, ds2, "s2" ); sig1 = s2_; //FIXME!!!!!!!!!!!!!!!!!!!! } DebuggingInsideLoop( Form( "s2 = %f", s2_ ) ); //std::cout << "s=" << _s << ", w3=" << MX2_ << ", sig1=" << sig1 << std::endl; const double sp = s_+MX2_-sig1, d3 = sig1-w2_; const double rl2 = sp*sp-4.*s_*MX2_; // lambda(s, m3**2, sigma) if ( rl2 <= 0. ) { InWarning( Form( "rl2 = %f <= 0", rl2 ) ); return false; } const double sl2 = sqrt( rl2 ); //std::cout << "ss=" << ss << ", sp=" << sp << ", sl1=" << sl1_ << ", sl2=" << sl2 << std::endl; double t1_max = w1_+MX2_-( ss*sp+sl1_*sl2 )/( 2.*s_ ), // definition from eq. (A.4) in [1] t1_min = ( w31_*d3+( d3-w31_ )*( d3*w1_-w31_*w2_ )/s_ )/t1_max; // definition from eq. (A.5) in [1] // FIXME dropped in CDF version const Kinematics::Limits q2_limits = cuts_.cuts.initial[Cuts::q2]; if ( t1_max > -q2_limits.min() ) { InWarning( Form( "t1max = %f > -q2min = %f", t1_max, -q2_limits.min() ) ); return false; } if ( t1_min < -q2_limits.max() && q2_limits.hasMax() ) { Debugging( Form( "t1min = %f < -q2max = %f", t1_min, -q2_limits.max() ) ); return false; } if ( t1_max < -q2_limits.max() && q2_limits.hasMax() ) t1_max = -q2_limits.max(); if ( t1_min > -q2_limits.min() ) t1_min = -q2_limits.min(); ///// // t1, the first photon propagator, is defined here t1_ = 0.; double dt1 = 0.; Map( x(0), t1_min, t1_max, t1_, dt1, "t1" ); // changes wrt mapt1 : dx->-dx dt1 = -dt1; DebuggingInsideLoop( Form( "Definition of t1 = %f according to\n\t" "(t1min, t1max) = (%f, %f)", t1_, t1_min, t1_max ) ); dd4_ = w4_-t1_; const double d8 = t1_-w2_, t13 = t1_-w1_-MX2_; sa1_ = -pow( t1_-w31_, 2 )/4.+w1_*t1_; if ( sa1_ >= 0. ) { InWarning( Form( "sa1_ = %f >= 0", sa1_ ) ); return false; } const double sl3 = sqrt( -sa1_ ); // one computes splus and (s2x=s2max) double splus, s2max; if ( w1_ != 0. ) { const double sb =( s_*( t1_-w31_ )+w12_*t13 )/( 2.*w1_ )+MX2_, sd = sl1_*sl3/w1_, se =( s_*( t1_*( s_+t13-w2_ )-w2_*w31_ )+MX2_*( w12_*d8+w2_*MX2_ ) )/w1_; if ( fabs( ( sb-sd )/sd )>=1. ) { splus = sb-sd; s2max = se/splus; } else { s2max = sb+sd; splus = se/s2max; } } else { // 3 s2max = ( s_*(t1_*(s_+d8-MX2_)-w2_*MX2_ )+w2_*MX2_*( w2_+MX2_-t1_ ) )/( ss*t13 ); splus = sig2; } // 4 double s2x = s2max; DebuggingInsideLoop( Form( "s2x = s2max = %f", s2x ) ); if ( n_opt_ < 0 ) { // 5 if ( splus > sig2 ) { sig2 = splus; DebuggingInsideLoop( Form( "sig2 truncated to splus = %f", splus ) ); } if ( n_opt_ < -1 ) { Map( x(2), sig2, s2max, s2_, ds2, "s2" ); } else { Mapla( t1_, w2_, x(2), sig2, s2max, s2_, ds2 ); } // n_opt_==-1 s2x = s2_; } else if ( n_opt_ == 0 ) { s2x = s2_; } // 6 DebuggingInsideLoop( Form( "s2x = %f", s2x ) ); // 7 const double r1 = s2x-d8, r2 = s2x-d6; const double rl4 = ( r1*r1-4.*w2_*s2x )*( r2*r2-4.*MY2_*s2x ); if ( rl4 <= 0. ) { DebuggingInsideLoop( Form( "rl4 = %f <= 0", rl4 ) ); return false; } const double sl4 = sqrt( rl4 ); // t2max, t2min definitions from eq. (A.12) and (A.13) in [1] const double t2_max = w2_+MY2_-( r1*r2+sl4 )/s2x * 0.5, t2_min = ( w52_*dd4_+( dd4_-w52_ )*( dd4_*w2_-w52_*t1_ )/s2x )/t2_max; // t2, the second photon propagator, is defined here t2_ = 0.; double dt2 = 0.; Map( x(1), t2_min, t2_max, t2_, dt2, "t2" ); // changes wrt mapt2 : dx->-dx dt2 = -dt2; // \f$\delta_6=m_4^2-m_5^2\f$ as defined in Vermaseren's paper const double tau = t1_-t2_, r3 = dd4_-t2_, r4 = w52_-t2_; DebuggingInsideLoop( Form( "r1 = %f\n\tr2 = %f\n\tr3 = %f\n\tr4 = %f", r1, r2, r3, r4 ) ); const double b = r3*r4-2.*( t1_+w2_ )*t2_, c = t2_*d6*d8+( d6-d8 )*( d6*w2_-d8*MY2_ ); const double t25 = t2_-w2_-MY2_; sa2_ = -r4*r4/4.+w2_*t2_; if ( sa2_ >= 0. ) { InWarning( Form( "sa2_ = %f >= 0", sa2_ ) ); return false; } const double sl6 = 2.*sqrt( -sa2_ ); g4_ = -r3*r3/4.+t1_*t2_; if ( g4_ >= 0. ) { InWarning( Form( "g4_ = %f >= 0", g4_ ) ); return false; } const double sl7 = 2.*sqrt( -g4_ ), sl5 = sl6*sl7; double s2p, s2min; if ( fabs( ( sl5-b )/sl5 ) >= 1. ) { s2p = ( sl5-b )/t2_ * 0.5; s2min = c/( t2_*s2p ); } else { // 8 s2min = ( -sl5-b )/t2_ * 0.5; s2p = c/( t2_*s2min ); } // 9 if ( n_opt_ > 1 ) Map( x( 2 ), s2min, s2max, s2_, ds2, "s2" ); else if ( n_opt_ == 1 ) Mapla( t1_, w2_, x( 2 ), s2min, s2max, s2_, ds2 ); const double ap = -0.25*pow( s2_+d8, 2 )+s2_*t1_; DebuggingInsideLoop( Form( "s2 = %f, s2max = %f, splus = %f", s2_, s2max, splus ) ); if ( w1_!=0. ) dd1_ = -0.25 * ( s2_-s2max ) * ( s2_-splus ) * w1_; // 10 else dd1_ = 0.25 * ( s2_-s2max ) * ss * t13; // 11 dd2_ = -t2_*( s2_-s2p )*( s2_-s2min ) * 0.25; DebuggingInsideLoop( Form( "t2 =%f\n\ts2 =%f\n\ts2p=%f\n\ts2min=%f\n\tdd2=%f", t2_, s2_, s2p, s2min, dd2_ ) ); const double yy4 = cos( M_PI*x( 3 ) ); const double dd = dd1_*dd2_; p12_ = ( s_-w1_-w2_ )*0.5; const double st = s2_-t1_-w2_; const double delb = ( 2.*w2_*r3+r4*st )*( 4.*p12_*t1_-( t1_-w31_ )*st )/( 16.*ap ); if ( dd <= 0. ) { DebuggingInsideLoop( Form( "dd = %e <= 0\n\tdd1 = %e\tdd2 = %e", dd, dd1_, dd2_ ) ); return false; } delta_ = delb - yy4*st*sqrt( dd )/ ap * 0.5; s1_ = t2_+w1_+( 2.*p12_*r3-4.*delta_ )/st; if ( ap >= 0. ) { DebuggingInsideLoop( Form( "ap = %f >= 0", ap ) ); return false; } jacobian_ = ds2 * dt1 * dt2 * M_PI*M_PI/( 8.*sl1_*sqrt( -ap ) ); DebuggingInsideLoop( Form( "Jacobian = %e", jacobian_ ) ); gram_ = ( 1.-yy4*yy4 )*dd/ap; p13_ = -t13 * 0.5; p14_ = (tau+s1_-MX2_) * 0.5; p25_ = -t25 * 0.5; /*const double //p15 = (s_+t2_-s1_-w2_)/2., //p23 = (s_+t1_-s2_-w1_)/2., //p24 = (s2_-tau-MY2_)/2., //p34 = (s1_-MX2_-w4_)/2., //p35 = (s_+w4_-s1_-s2_)/2., //p45 = (s2_-w4_-MY2_)/2.;*/ p1k2_ = ( s1_-t2_-w1_ ) * 0.5; p2k1_ = st * 0.5; double s1p, s1m; if ( w2_ != 0. ) { const double sbb = ( s_*( t2_-w52_ )-w12_*t25 )/w2_ * 0.5 + MY2_, sdd = sl1_*sl6/w2_ * 0.5, see = ( s_*( t2_*( s_+t25-w1_ )-w1_*w52_ )+MY2_*( w1_*MY2_-w12_*( t2_-w1_ ) ) )/w2_; if ( sbb/sdd >= 0. ) { s1p = sbb+sdd; s1m = see/s1p; } else { s1m = sbb-sdd; s1p = see/s1m; } // 12 dd3_ = -w2_*( s1p-s1_ )*( s1m-s1_ ) * 0.25; // 13 } else { // 14 s1p = ( s_*( t2_*( s_-MY2_+t2_-w1_ )-w1_*MY2_ )+w1_*MY2_*( w1_+MY2_-t2_ ) )/( t25*( s_-w12_ ) ); dd3_ = -t25*( s_-w12_ )*( s1p-s1_ ) * 0.25; } // 15 //const double acc3 = (s1p-s1_)/(s1p+s1_); const double ssb = t2_+w1_-r3*( w31_-t1_ )/t1_ * 0.5, ssd = sl3*sl7/t1_, sse = ( t2_-w1_ )*( w4_-MX2_ )+( t2_-w4_+w31_ )*( ( t2_-w1_)*MX2_-(w4_-MX2_)*w1_)/t1_; double s1pp, s1pm; if (ssb/ssd>=0.) { s1pp = ssb+ssd; s1pm = sse/s1pp; } else { s1pm = ssb-ssd; s1pp = sse/s1pm; } // 16 // 17 dd4_ = -t1_*( s1_-s1pp )*( s1_-s1pm ) * 0.25; //const double acc4 = ( s1_-s1pm )/( s1_+s1pm ); dd5_ = dd1_+dd3_+( ( p12_*( t1_-w31_ )*0.5-w1_*p2k1_ )*( p2k1_*( t2_-w52_ )-w2_*r3 )-delta_*( 2.*p12_*p2k1_-w2_*( t1_-w31_ ) ) ) / p2k1_; return true; } bool GamGamLL::orient() { if ( !pickin() or jacobian_ == 0. ) { DebuggingInsideLoop( Form( "Pickin failed! Jacobian = %f", jacobian_ ) ); return false; } const double re = 0.5 / sqs_; ep1_ = re*( s_+w12_ ); ep2_ = re*( s_-w12_ ); DebuggingInsideLoop( Form( " re = %e\n\tw12_ = %e", re, w12_ ) ); DebuggingInsideLoop( Form( "Incoming particles' energy = %f, %f", ep1_, ep2_ ) ); p_cm_ = re*sl1_; de3_ = re*(s2_-MX2_+w12_); de5_ = re*(s1_-MY2_-w12_); // Final state energies const double ep3 = ep1_-de3_, ep5 = ep2_-de5_; ec4_ = de3_+de5_; if ( ec4_ < mc4_ ) { InWarning( Form( "ec4_ = %f < mc4_ = %f\n\t==> de3 = %f, de5 = %f", ec4_, mc4_, de3_, de5_ ) ); return false; } // What if the protons' momenta are not along the z-axis? pc4_ = sqrt( ec4_*ec4_-mc4_*mc4_ ); if ( pc4_==0. ) { InWarning( "pzc4==0" ); return false; } const double pp3 = sqrt( ep3*ep3-MX2_ ), pt3 = sqrt( dd1_/s_ )/p_cm_, pp5 = sqrt( ep5*ep5-MY2_ ), pt5 = sqrt( dd3_/s_ )/p_cm_; DebuggingInsideLoop( Form( "Central system's energy: E4 = %f\n\t" " momentum: p4 = %f\n\t" " invariant mass: m4 = %f\n\t" "Outgoing particles' energy: E3 = %f\n\t" " E5 = %f", ec4_, pc4_, mc4_, ep3, ep5 ) ); const double sin_theta3 = pt3/pp3, sin_theta5 = pt5/pp5; DebuggingInsideLoop( Form( "sin_theta3 = %e\n\tsin_theta5 = %e", sin_theta3, sin_theta5 ) ); if (sin_theta3>1.) { InWarning( Form( "sin_theta3 = %e > 1", sin_theta3 ) ); return false; } if (sin_theta5>1.) { InWarning( Form( "sin_theta5 = %e > 1", sin_theta5 ) ); return false; } double ct3 = sqrt( 1.-sin_theta3*sin_theta3 ), ct5 = sqrt( 1.-sin_theta5*sin_theta5 ); if ( ep1_*ep3 < p13_ ) ct3 *= -1.; if ( ep2_*ep5 > p25_ ) ct5 *= -1.; DebuggingInsideLoop( Form( "ct3 = %e\n\tct5 = %e", ct3, ct5 ) ); if ( dd5_ < 0. ) { InWarning( Form( "dd5 = %f < 0", dd5_ ) ); return false; } // Centre of mass system kinematics (theta4 and phi4) pt4_ = sqrt( dd5_/s_ )/p_cm_; sin_theta4_ = pt4_/pc4_; if ( sin_theta4_ > 1. ) { InWarning( Form( "st4 = %f > 1", sin_theta4_ ) ); return false; } cos_theta4_ = sqrt( 1.-sin_theta4_*sin_theta4_ ); if ( ep1_*ec4_ < p14_ ) cos_theta4_ *= -1.; al4_ = 1.-cos_theta4_; be4_ = 1.+cos_theta4_; if (cos_theta4_<0.) be4_ = sin_theta4_*sin_theta4_/al4_; else al4_ = sin_theta4_*sin_theta4_/be4_; DebuggingInsideLoop( Form( "ct4 = %f\n\tal4 = %f, be4 = %f", cos_theta4_, al4_, be4_ ) ); const double rr = sqrt( -gram_/s_ )/( p_cm_*pt4_ ); const double sin_phi3 = rr / pt3, sin_phi5 = -rr / pt5; if ( fabs( sin_phi3 ) > 1. ) { InWarning( Form( "sin(phi_3) = %e while it must be in [-1 ; 1]", sin_phi3 ) ); return false; } if ( fabs( sin_phi5 ) > 1. ) { InWarning( Form( "sin(phi_5) = %e while it must be in [-1 ; 1]", sin_phi5 ) ); return false; } const double cos_phi3 = -sqrt( 1.-sin_phi3*sin_phi3 ), cos_phi5 = -sqrt( 1.-sin_phi5*sin_phi5 ); p3_lab_ = Particle::Momentum( pp3*sin_theta3*cos_phi3, pp3*sin_theta3*sin_phi3, pp3*ct3, ep3 ); p5_lab_ = Particle::Momentum( pp5*sin_theta5*cos_phi5, pp5*sin_theta5*sin_phi5, pp5*ct5, ep5 ); const double a1 = p3_lab_.px()-p5_lab_.px(); DebuggingInsideLoop( Form( "Kinematic quantities\n\t" "cos(theta3) = %1.4f\tsin(theta3) = %1.4f\tcos( phi3 ) = %1.4f\tsin( phi3 ) = %1.4f\n\t" "cos(theta4) = %1.4f\tsin(theta4) = %1.4f\n\t" "cos(theta5) = %1.4f\tsin(theta5) = %1.4f\tcos( phi5 ) = %1.4f\tsin( phi5 ) = %1.4f\n\t" "a1 = %f", ct3, sin_theta3, cos_phi3, sin_phi3, cos_theta4_, sin_theta4_, ct5, sin_theta5, cos_phi5, sin_phi5, a1 ) ); if ( fabs( pt4_+p3_lab_.px()+p5_lab_.px() ) < fabs( fabs( a1 )-pt4_ ) ) { DebuggingInsideLoop( Form( "|pp4+pp3*cos(phi3)+pp5*cos(phi5)| < | |a1|-pp4 |\n\t" "pp4 = %f\tpp5 = %f\n\t" "cos(phi3) = %f\tcos(phi5) = %f" "a1 = %f", pt4_, pt5, cos_phi3, cos_phi5, a1 ) ); return true; } if ( a1 < 0. ) p5_lab_.setP( 0, -p5_lab_.px() ); else p3_lab_.setP( 0, -p3_lab_.px() ); return true; } double GamGamLL::computeOutgoingPrimaryParticlesMasses( double x, double outmass, double lepmass, double& dw ) { const double mx0 = Particle::massFromPDGId( Particle::Proton )+Particle::massFromPDGId( Particle::PiPlus ); // 1.07 const Kinematics::Limits mx_limits = cuts_.cuts.remnants[Cuts::mass]; const double wx2min = pow( std::max( mx0, mx_limits.min() ), 2 ), wx2max = pow( std::min( sqs_-outmass-2.*lepmass, mx_limits.max() ), 2 ); double mx2 = 0., dmx2 = 0.; Map( x, wx2min, wx2max, mx2, dmx2, "mx2" ); DebuggingInsideLoop( Form( "mX^2 in range (%f, %f), x = %f\n\t" "mX^2 = %f, d(mX^2) = %f\n\t" "mX = %f, d(mX) = %f", wx2min, wx2max, x, mx2, dmx2, sqrt( mx2 ), sqrt( dmx2 ) ) ); dw = sqrt( dmx2 ); return sqrt( mx2 ); } void GamGamLL::beforeComputeWeight() { if ( !GenericProcess::is_point_set_ ) return; const Particle& p1 = event_->getOneByRole( Particle::IncomingBeam1 ), &p2 = event_->getOneByRole( Particle::IncomingBeam2 ); ep1_ = p1.energy(); ep2_ = p2.energy(); Ml2_ = event_->getByRole( Particle::CentralSystem )[0].mass2(); switch ( cuts_.mode ) { case Kinematics::ElectronProton: default: { InError( "Case not yet supported!" ); } break; case Kinematics::ElasticElastic: dw31_ = dw52_ = 0.; break; case Kinematics::InelasticElastic: { const double m = computeOutgoingPrimaryParticlesMasses( x( 7 ), p1.mass(), sqrt( Ml2_ ), dw31_ ); event_->getOneByRole( Particle::OutgoingBeam1 ).setMass( m ); event_->getOneByRole( Particle::OutgoingBeam2 ).setMass( Particle::massFromPDGId( p2.pdgId() ) ); } break; case Kinematics::ElasticInelastic: { const double m = computeOutgoingPrimaryParticlesMasses( x( 7 ), p2.mass(), sqrt( Ml2_ ), dw52_ ); event_->getOneByRole( Particle::OutgoingBeam1 ).setMass( Particle::massFromPDGId( p1.pdgId() ) ); event_->getOneByRole( Particle::OutgoingBeam2 ).setMass( m ); } break; case Kinematics::InelasticInelastic: { const double mx = computeOutgoingPrimaryParticlesMasses( x( 7 ), p2.mass(), sqrt( Ml2_ ), dw31_ ); event_->getOneByRole( Particle::OutgoingBeam1 ).setMass( mx ); const double my = computeOutgoingPrimaryParticlesMasses( x( 8 ), p1.mass(), sqrt( Ml2_ ), dw52_ ); event_->getOneByRole( Particle::OutgoingBeam2 ).setMass( my ); } break; } MX_ = event_->getOneByRole( Particle::OutgoingBeam1 ).mass(); MY_ = event_->getOneByRole( Particle::OutgoingBeam2 ).mass(); MX2_ = MX_*MX_; MY2_ = MY_*MY_; } double GamGamLL::computeWeight() { if ( !is_outgoing_state_set_ ) { InWarning( "Output state not set!" ); return 0.; } DebuggingInsideLoop( Form( "sqrt(s)=%f\n\tm(X1)=%f\tm(X2)=%f", sqs_, MX_, MY_ ) ); Kinematics::Limits& w_limits = cuts_.cuts.initial[Cuts::w]; if ( !w_limits.hasMax() ) w_limits.max() = s_; // The minimal energy for the central system is its outgoing leptons' mass energy (or wmin_ if specified) const double wmin = std::max( 4.*Ml2_, w_limits.min() ); // The maximal energy for the central system is its CM energy with the outgoing particles' mass energy substracted (or _wmax if specified) const double wmax = std::min( pow( sqs_-MX_-MY_, 2 ), w_limits.max() ); DebuggingInsideLoop( Form( "wmin = %f\n\twmax = %f\n\twmax/wmin = %f", wmin, wmax, wmax/wmin ) ); // compute the two-photon energy for this point w4_ = 0.; double dw4 = 0.; Map( x( 4 ), wmin, wmax, w4_, dw4, "w4" ); Map( x( 4 ), wmin, wmax, w4_, dw4, "w4" ); mc4_ = sqrt( w4_ ); DebuggingInsideLoop( Form( "Computed value for w4 = %f -> mc4 = %f", w4_, mc4_ ) ); if ( !orient() ) return 0.; if ( jacobian_ == 0. ) { InWarning( Form( "dj = %f", jacobian_ ) ); return 0.; } if ( t1_>0. ) { InWarning( Form( "t1 = %f > 0", t1_ ) ); return 0.; } if ( t2_>0. ) { InWarning( Form( "t2 = %f > 0", t2_ ) ); return 0.; } const double ecm6 = w4_ / ( 2.*mc4_ ), pp6cm = sqrt( ecm6*ecm6-Ml2_ ); jacobian_ *= dw4*pp6cm/( mc4_*Constants::sconstb*s_ ); // Let the most obscure part of this code begin... const double e1mp1 = w1_ / ( ep1_+p_cm_ ), e3mp3 = MX2_ / ( p3_lab_.energy()+p3_lab_.p() ); const double al3 = pow( sin( p3_lab_.theta() ), 2 )/( 1.+( p3_lab_.theta() ) ); // 2-photon system kinematics ?! const double eg = ( w4_+t1_-t2_ )/( 2.*mc4_ ); double pg = sqrt( eg*eg-t1_ ); const double pgx = -p3_lab_.px()*cos_theta4_-sin_theta4_*( de3_-e1mp1 + e3mp3 + p3_lab_.p()*al3 ), pgy = -p3_lab_.py(), pgz = mc4_*de3_/( ec4_+pc4_ )-ec4_*de3_*al4_/mc4_-p3_lab_.px()*ec4_*sin_theta4_/mc4_+ec4_*cos_theta4_/mc4_*( p3_lab_.p()*al3+e3mp3-e1mp1 ); DebuggingInsideLoop( Form( "pg3 = (%f, %f, %f)\n\t" "pg3^2 = %f", pgx, pgy, pgz, sqrt( pgx*pgx+pgy*pgy+pgz*pgz ) ) ); const double pgp = sqrt( pgx*pgx + pgy*pgy ), // outgoing proton (3)'s transverse momentum pgg = sqrt( pgp*pgp + pgz*pgz ); // outgoing proton (3)'s momentum if ( pgg > pgp*0.9 && pgg > pg ) { pg = pgg; } //FIXME ??? // Phi angle for the 2-photon system ?! const double cpg = pgx/pgp, spg = pgy/pgp; // Theta angle for the 2-photon system ?! const double stg = pgp/pg; const int theta_sign = ( pgz>0. ) ? 1 : -1; const double ctg = theta_sign*sqrt( 1.-stg*stg ); double xx6 = x( 5 ); const double amap = 0.5 * (w4_-t1_-t2_), bmap = 0.5 * sqrt( ( pow( w4_-t1_-t2_, 2 )-4.*t1_*t2_ )*( 1.-4.*Ml2_/w4_ ) ), ymap = ( amap+bmap )/( amap-bmap ), beta = pow( ymap, 2.*xx6-1. ); xx6 = 0.5 * ( 1. + amap/bmap*( beta-1. )/( beta+1. ) ); xx6 = std::max( 0., std::min( xx6, 1. ) ); // xx6 in [0., 1.] DebuggingInsideLoop( Form( "amap = %f\n\tbmap = %f\n\tymap = %f\n\tbeta = %f", amap, bmap, ymap, beta ) ); // 3D rotation of the first outgoing lepton wrt the CM system const double theta6cm = acos( 1.-2.*xx6 ); // match the Jacobian jacobian_ *= ( amap+bmap*cos( theta6cm ) ); jacobian_ *= ( amap-bmap*cos( theta6cm ) ); jacobian_ /= amap; jacobian_ /= bmap; jacobian_ *= log( ymap ); jacobian_ *= 0.5; DebuggingInsideLoop( Form( "Jacobian = %e", jacobian_ ) ); DebuggingInsideLoop( Form( "ctcm6 = %f\n\tstcm6 = %f", cos( theta6cm ), sin( theta6cm ) ) ); const double phi6cm = 2.*M_PI*x( 6 ); // First outgoing lepton's 3-momentum in the centre of mass system Particle::Momentum p6cm = Particle::Momentum::fromPThetaPhi( pp6cm, theta6cm, phi6cm ); DebuggingInsideLoop( Form( "p3cm6 = (%f, %f, %f)", p6cm.px(), p6cm.py(), p6cm.pz() ) ); const double h1 = stg*p6cm.pz()+ctg*p6cm.px(); const double pc6z = ctg*p6cm.pz()-stg*p6cm.px(), pc6x = cpg*h1-spg*p6cm.py(); const double qcx = 2.*pc6x, qcz = 2.*pc6z; // qcy == QCY is never defined const double el6 = ( ec4_*ecm6+pc4_*pc6z ) / mc4_, h2 = ( ec4_*pc6z+pc4_*ecm6 ) / mc4_; DebuggingInsideLoop( Form( "h1 = %f\n\th2 = %f", h1, h2 ) ); // first outgoing lepton's 3-momentum const double p6x = cos_theta4_*pc6x+sin_theta4_*h2, p6y = cpg*p6cm.py()+spg*h1, p6z = cos_theta4_*h2-sin_theta4_*pc6x; // first outgoing lepton's kinematics p6_cm_ = Particle::Momentum( p6x, p6y, p6z, el6 ); DebuggingInsideLoop( Form( "E6(cm) = %f\n\tP6(cm) = (%f, %f, %f)", el6, p6x, p6y, p6z ) ); const double hq = ec4_*qcz/mc4_; const Particle::Momentum qve( cos_theta4_*qcx+sin_theta4_*hq, 2.*p6y, cos_theta4_*hq-sin_theta4_*qcx, pc4_*qcz/mc4_ // energy ); // Available energy for the second lepton is the 2-photon system's energy with the first lepton's energy removed const double el7 = ec4_-el6; DebuggingInsideLoop( Form( "Outgoing kinematics\n\t" " first outgoing lepton: p = %f, E = %f\n\t" "second outgoing lepton: p = %f, E = %f", p6_cm_.p(), p6_cm_.energy(), p7_cm_.p(), p6_cm_.energy() ) ); // Second outgoing lepton's 3-momentum const double p7x = -p6x + pt4_, p7y = -p6y, p7z = -p6z + pc4_*cos_theta4_; // second outgoing lepton's kinematics p7_cm_ = Particle::Momentum( p7x, p7y, p7z, el7 ); //p6_cm_ = Particle::Momentum(pl6*st6*cp6, pl6*st6*sp6, pl6*ct6, el6); //p7_cm_ = Particle::Momentum(pl7*st7*cp7, pl7*st7*sp7, pl7*ct7, el7); q1dq_ = eg*( 2.*ecm6-mc4_ )-2.*pg*p6cm.pz(); q1dq2_ = ( w4_-t1_-t2_ ) * 0.5; const double phi3 = p3_lab_.phi(), cos_phi3 = cos( phi3 ), sin_phi3 = sin( phi3 ), phi5 = p5_lab_.phi(), cos_phi5 = cos( phi5 ), sin_phi5 = sin( phi5 ); bb_ = t1_*t2_+( w4_*pow( sin( theta6cm ), 2 ) + 4.*Ml2_*pow( cos( theta6cm ), 2 ) )*pg*pg; const double c1 = p3_lab_.pt() * ( qve.px()*sin_phi3 - qve.py()*cos_phi3 ), c2 = p3_lab_.pt() * ( qve.pz()*ep1_ - qve.energy() *p_cm_ ), c3 = ( w31_*ep1_*ep1_ + 2.*w1_*de3_*ep1_ - w1_*de3_*de3_ + p3_lab_.pt2()*ep1_*ep1_ ) / ( p3_lab_.energy()*p_cm_ + p3_lab_.pz()*ep1_ ); const double b1 = p5_lab_.pt() * ( qve.px()*sin_phi5 - qve.py()*cos_phi5 ), b2 = p5_lab_.pt() * ( qve.pz()*ep2_ + qve.energy() *p_cm_ ), b3 = ( w52_*ep2_*ep2_ + 2.*w2_*de5_*ep2_ - w2_*de5_*de5_ + p5_lab_.pt2()*ep2_*ep2_ ) / ( ep2_*p5_lab_.pz() - p5_lab_.energy()*p_cm_ ); const double r12 = c2*sin_phi3 + qve.py()*c3, r13 = -c2*cos_phi3 - qve.px()*c3; const double r22 = b2*sin_phi5 + qve.py()*b3, r23 = -b2*cos_phi5 - qve.px()*b3; epsi_ = p12_*c1*b1 + r12*r22 + r13*r23; g5_ = w1_*c1*c1 + r12*r12 + r13*r13; g6_ = w2_*b1*b1 + r22*r22 + r23*r23; a5_ = -( qve.px()*cos_phi3 + qve.py()*sin_phi3 )*p3_lab_.pt()*p1k2_-( ep1_*qve.energy()-p_cm_*qve.pz() )*( cos_phi3*cos_phi5 + sin_phi3*sin_phi5 )*p3_lab_.pt()*p5_lab_.pt() + ( de5_*qve.pz()+qve.energy()*( p_cm_+p5_lab_.pz() ) )*c3; a6_ = -( qve.px()*cos_phi5 + qve.py()*sin_phi5 )*p5_lab_.pt()*p2k1_-( ep2_*qve.energy()+p_cm_*qve.pz() )*( cos_phi3*cos_phi5 + sin_phi3*sin_phi5 )*p3_lab_.pt()*p5_lab_.pt() + ( de3_*qve.pz()-qve.energy()*( p_cm_-p3_lab_.pz() ) )*b3; DebuggingInsideLoop( Form( "a5 = %f\n\ta6 = %f", a5_, a6_ ) ); //////////////////////////////////////////////////////////////// // END of GAMGAMLL subroutine in the FORTRAN version //////////////////////////////////////////////////////////////// const Particle::Momentum cm = event_->getOneByRole( Particle::IncomingBeam1 ).momentum() + event_->getOneByRole( Particle::IncomingBeam2 ).momentum(); //////////////////////////////////////////////////////////////// // INFO from f.f //////////////////////////////////////////////////////////////// const double gamma = cm.energy() / sqs_, betgam = cm.pz() / sqs_; //--- kinematics computation for both leptons p6_cm_.betaGammaBoost( gamma, betgam ); p7_cm_.betaGammaBoost( gamma, betgam ); //--- cut on mass of final hadronic system (MX/Y) const Kinematics::Limits mx_limits = cuts_.cuts.remnants[Cuts::mass]; if ( cuts_.mode == Kinematics::InelasticElastic || cuts_.mode == Kinematics::InelasticInelastic ) { if ( mx_limits.hasMin() && MX_ < mx_limits.min() ) return 0.; if ( mx_limits.hasMax() && MX_ > mx_limits.max() ) return 0.; } if ( cuts_.mode == Kinematics::ElasticInelastic || cuts_.mode == Kinematics::InelasticInelastic ) { if ( mx_limits.hasMin() && MY_ < mx_limits.min() ) return 0.; if ( mx_limits.hasMax() && MY_ > mx_limits.max() ) return 0.; } //--- cut on the proton's Q2 (first photon propagator T1) const Kinematics::Limits q2_limits = cuts_.cuts.initial[Cuts::q2]; if ( q2_limits.hasMax() && t1_<-q2_limits.max() ) return 0; if ( q2_limits.hasMin() && t1_>-q2_limits.min() ) return 0.; //--- cuts on outgoing leptons' kinematics const Kinematics::Limits m_limits = cuts_.cuts.central[Cuts::mass_sum]; if ( m_limits.hasMin() && ( p6_cm_+p7_cm_ ).mass() < m_limits.min() ) return 0.; if ( m_limits.hasMax() && ( p6_cm_+p7_cm_ ).mass() > m_limits.max() ) return 0.; //----- cuts on the individual leptons const Kinematics::Limits pt_limits = cuts_.cuts.central[Cuts::pt_single]; if ( pt_limits.hasMin() && ( p6_cm_.pt() < pt_limits.min() || p7_cm_.pt() < pt_limits.min() ) ) return 0.; if ( pt_limits.hasMax() && ( p6_cm_.pt() > pt_limits.max() || p7_cm_.pt() > pt_limits.max() ) ) return 0.; const Kinematics::Limits energy_limits = cuts_.cuts.central[Cuts::energy_single]; if ( energy_limits.hasMin() && ( p6_cm_.energy() < energy_limits.min() || p7_cm_.energy() < energy_limits.min() ) ) return 0.; if ( energy_limits.hasMax() && ( p6_cm_.energy() > energy_limits.max() || p7_cm_.energy() > energy_limits.max() ) ) return 0.; const Kinematics::Limits eta_limits = cuts_.cuts.central[Cuts::eta_single]; if ( eta_limits.hasMin() && ( p6_cm_.eta() < eta_limits.min() || p7_cm_.eta() < eta_limits.min() ) ) return 0.; if ( eta_limits.hasMax() && ( p6_cm_.eta() > eta_limits.max() || p7_cm_.eta() > eta_limits.max() ) ) return 0.; //--- compute the structure functions factors switch ( cuts_.mode ) { // inherited from CDF version case Kinematics::ElectronProton: default: jacobian_ *= periPP( 1, 2 ); break; // ep case case Kinematics::ElasticElastic: jacobian_ *= periPP( 2, 2 ); break; // elastic case case Kinematics::InelasticElastic: jacobian_ *= periPP( 3, 2 )*( dw31_*dw31_ ); break; case Kinematics::ElasticInelastic: jacobian_ *= periPP( 3, 2 )*( dw52_*dw52_ ); break; // single-dissociative case case Kinematics::InelasticInelastic: jacobian_ *= periPP( 3, 3 )*( dw31_*dw31_ )*( dw52_*dw52_ ); break; // double-dissociative case } //--- compute the event weight using the Jacobian return Constants::GeV2toBarn*jacobian_; } void GamGamLL::fillKinematics( bool ) { const Particle::Momentum cm = event_->getOneByRole( Particle::IncomingBeam1 ).momentum() + event_->getOneByRole( Particle::IncomingBeam2 ).momentum(); const double gamma = cm.energy()/sqs_, betgam = cm.pz()/sqs_; Particle::Momentum plab_ip1 = Particle::Momentum( 0., 0., p_cm_, ep1_ ).betaGammaBoost( gamma, betgam ); Particle::Momentum plab_ip2 = Particle::Momentum( 0., 0., -p_cm_, ep2_ ).betaGammaBoost( gamma, betgam ); p3_lab_.betaGammaBoost( gamma, betgam ); p5_lab_.betaGammaBoost( gamma, betgam ); //----- needed to parametrise a random rotation around z-axis const int rany = ( rand() >= RAND_MAX*0.5 ) ? 1 : -1, ransign = ( rand() >= RAND_MAX*0.5 ) ? 1 : -1; const double ranphi = ( rand()/RAND_MAX*0.5 )*2.*M_PI; Particle::Momentum plab_ph1 = ( plab_ip1-p3_lab_ ).rotatePhi( ranphi, rany ); Particle::Momentum plab_ph2 = ( plab_ip2-p5_lab_ ).rotatePhi( ranphi, rany ); p3_lab_.rotatePhi( ranphi, rany ); p5_lab_.rotatePhi( ranphi, rany ); p6_cm_.rotatePhi( ranphi, rany ); p7_cm_.rotatePhi( ranphi, rany ); /*if ( symmetrise_ && rand() >= .5*RAND_MAX ) { p6_cm_.mirrorZ(); p7_cm_.mirrorZ(); }*/ //----- incoming protons event_->getOneByRole( Particle::IncomingBeam1 ).setMomentum( plab_ip1 ); event_->getOneByRole( Particle::IncomingBeam2 ).setMomentum( plab_ip2 ); //----- first outgoing proton Particle& op1 = event_->getOneByRole( Particle::OutgoingBeam1 ); op1.setMomentum( p3_lab_ ); switch ( cuts_.mode ) { case Kinematics::ElasticElastic: case Kinematics::ElasticInelastic: default: op1.setStatus( Particle::FinalState ); // stable proton break; case Kinematics::InelasticElastic: case Kinematics::InelasticInelastic: - op1.setStatus( Particle::Undecayed ); // fragmenting remnants + op1.setStatus( Particle::Unfragmented ); // fragmenting remnants op1.setMass( MX_ ); break; } //----- second outgoing proton Particle& op2 = event_->getOneByRole( Particle::OutgoingBeam2 ); op2.setMomentum( p5_lab_ ); switch ( cuts_.mode ) { case Kinematics::ElasticElastic: case Kinematics::InelasticElastic: default: op2.setStatus( Particle::FinalState ); // stable proton break; case Kinematics::ElasticInelastic: case Kinematics::InelasticInelastic: - op2.setStatus( Particle::Undecayed ); // fragmenting remnants + op2.setStatus( Particle::Unfragmented ); // fragmenting remnants op2.setMass( MY_ ); break; } //----- first incoming photon Particle& ph1 = event_->getOneByRole( Particle::Parton1 ); ph1.setMomentum( plab_ph1 ); ph1.setStatus( Particle::Incoming ); //----- second incoming photon Particle& ph2 = event_->getOneByRole( Particle::Parton2 ); ph2.setMomentum( plab_ph2 ); ph2.setStatus( Particle::Incoming ); Particles& central_system = event_->getByRole( Particle::CentralSystem ); //----- first outgoing lepton Particle& ol1 = central_system[0]; ol1.setPdgId( ol1.pdgId(), ransign ); ol1.setMomentum( p6_cm_ ); ol1.setStatus( Particle::FinalState ); //----- second outgoing lepton Particle& ol2 = central_system[1]; ol2.setPdgId( ol2.pdgId(), -ransign ); ol2.setMomentum( p7_cm_ ); ol2.setStatus( Particle::FinalState ); //----- intermediate two-lepton system event_->getOneByRole( Particle::Intermediate ).setMomentum( p6_cm_+p7_cm_ ); } double GamGamLL::periPP( int nup_, int ndown_ ) { DebuggingInsideLoop( Form( " Nup = %d\n\tNdown = %d", nup_, ndown_ ) ); FormFactors fp1, fp2; GenericProcess::formFactors( -t1_, -t2_, fp1, fp2 ); DebuggingInsideLoop( Form( "u1 = %f\n\tu2 = %f\n\tv1 = %f\n\tv2 = %f", fp1.FM, fp1.FE, fp2.FM, fp2.FE ) ); const double qqq = q1dq_*q1dq_, qdq = 4.*Ml2_-w4_; const double t11 = 64. *( bb_*( qqq-g4_-qdq*( t1_+t2_+2.*Ml2_ ) )-2.*( t1_+2.*Ml2_ )*( t2_+2.*Ml2_ )*qqq ) * t1_*t2_, // magnetic-magnetic t12 = 128.*( -bb_*( dd2_+g6_ )-2.*( t1_+2.*Ml2_ )*( sa2_*qqq+a6_*a6_ ) ) * t1_, // electric-magnetic t21 = 128.*( -bb_*( dd4_+g5_ )-2.*( t2_+2.*Ml2_ )*( sa1_*qqq+a5_*a5_ ) ) * t2_, // magnetic-electric t22 = 512.*( bb_*( delta_*delta_-gram_ )-pow( epsi_-delta_*( qdq+q1dq2_ ), 2 )-sa1_*a6_*a6_-sa2_*a5_*a5_-sa1_*sa2_*qqq ); // electric-electric const double peripp = ( fp1.FM*fp2.FM*t11 + fp1.FE*fp2.FM*t21 + fp1.FM*fp2.FE*t12 + fp1.FE*fp2.FE*t22 ) / pow( 2.*t1_*t2_*bb_, 2 ); DebuggingInsideLoop( Form( "t11 = %5.2f\tt12 = %5.2f\n\t" "t21 = %5.2f\tt22 = %5.2f\n\t" "=> PeriPP = %e", t11, t12, t21, t22, peripp ) ); return peripp; } diff --git a/CepGen/Processes/GenericKTProcess.cpp b/CepGen/Processes/GenericKTProcess.cpp index f186176..82e35a2 100644 --- a/CepGen/Processes/GenericKTProcess.cpp +++ b/CepGen/Processes/GenericKTProcess.cpp @@ -1,269 +1,269 @@ #include "GenericKTProcess.h" #include "CepGen/StructureFunctions/StructureFunctionsBuilder.h" #include "CepGen/Core/Exception.h" namespace CepGen { namespace Process { GenericKTProcess::GenericKTProcess( const std::string& name, const std::string& description, const unsigned int& num_user_dimensions, const std::array<Particle::ParticleCode,2>& partons, const std::vector<Particle::ParticleCode>& central ) : GenericProcess( name, description+" (kT-factorisation approach)" ), kNumUserDimensions( num_user_dimensions ), kIntermediateParts( partons ), kProducedParts( central ) {} GenericKTProcess::~GenericKTProcess() {} void GenericKTProcess::prepareKTKinematics() { DebuggingInsideLoop("Dummy kinematics prepared!"); } void GenericKTProcess::addEventContent() { GenericProcess::setEventContent( { // incoming state { Particle::IncomingBeam1, Particle::Proton }, { Particle::IncomingBeam2, Particle::Proton }, { Particle::Parton1, kIntermediateParts[0] }, { Particle::Parton2, kIntermediateParts[1] } }, { // outgoing state { Particle::OutgoingBeam1, { Particle::Proton } }, { Particle::OutgoingBeam2, { Particle::Proton } }, { Particle::CentralSystem, kProducedParts } } ); setExtraContent(); } unsigned int GenericKTProcess::numDimensions( const Kinematics::ProcessMode& process_mode ) const { switch ( process_mode ) { default: case Kinematics::ElasticElastic: return kNumRequiredDimensions+kNumUserDimensions; case Kinematics::ElasticInelastic: case Kinematics::InelasticElastic: return kNumRequiredDimensions+kNumUserDimensions+1; case Kinematics::InelasticInelastic: return kNumRequiredDimensions+kNumUserDimensions+2; } } void GenericKTProcess::addPartonContent() { // Incoming partons qt1_ = exp( log_qmin_+( log_qmax_-log_qmin_ )*x( 0 ) ); qt2_ = exp( log_qmin_+( log_qmax_-log_qmin_ )*x( 1 ) ); phi_qt1_ = 2.*M_PI*x( 2 ); phi_qt2_ = 2.*M_PI*x( 3 ); DebuggingInsideLoop( Form( "photons transverse virtualities (qt):\n\t" " mag = %f / %f (%.2f < log(qt) < %.2f)\n\t" " phi = %f / %f", qt1_, qt2_, log_qmin_, log_qmax_, phi_qt1_, phi_qt2_ ) ); } void GenericKTProcess::setKinematics( const Kinematics& kin ) { cuts_ = kin; log_qmin_ = -10.; // FIXME //log_qmin_ = std::log( std::sqrt( cuts_.q2min ) ); log_qmax_ = log( cuts_.cuts.initial[Cuts::qt].max() ); } double GenericKTProcess::computeWeight() { addPartonContent(); prepareKTKinematics(); computeOutgoingPrimaryParticlesMasses(); const double jac = computeJacobian(), integrand = computeKTFactorisedMatrixElement(), weight = jac*integrand; DebuggingInsideLoop( Form( "Jacobian = %f\n\tIntegrand = %f\n\tdW = %f", jac, integrand, weight ) ); return weight; } void GenericKTProcess::computeOutgoingPrimaryParticlesMasses() { const unsigned int op_index = kNumRequiredDimensions+kNumUserDimensions; const Kinematics::Limits remn_mx_cuts = cuts_.cuts.remnants[Cuts::mass]; switch ( cuts_.mode ) { case Kinematics::ElectronProton: default: { InError( "This kT factorisation process is intended for p-on-p collisions! Aborting!" ); exit( 0 ); } break; case Kinematics::ElasticElastic: { MX_ = event_->getOneByRole( Particle::IncomingBeam1 ).mass(); MY_ = event_->getOneByRole( Particle::IncomingBeam2 ).mass(); } break; case Kinematics::ElasticInelastic: { const double mx_min = remn_mx_cuts.min(), mx_range = remn_mx_cuts.range(); MX_ = event_->getOneByRole( Particle::IncomingBeam1 ).mass(); MY_ = mx_min + mx_range*x( op_index ); } break; case Kinematics::InelasticElastic: { const double mx_min = remn_mx_cuts.min(), mx_range = remn_mx_cuts.range(); MX_ = mx_min + mx_range*x( op_index ); MY_ = event_->getOneByRole( Particle::IncomingBeam2 ).mass(); } break; case Kinematics::InelasticInelastic: { const double mx_min = remn_mx_cuts.min(), mx_range = remn_mx_cuts.range(); MX_ = mx_min + mx_range*x( op_index ); MY_ = mx_min + mx_range*x( op_index+1 ); } break; } DebuggingInsideLoop( Form( "outgoing remnants invariant mass: %f / %f (%.2f < M(X/Y) < %.2f)", MX_, MY_, remn_mx_cuts.min(), remn_mx_cuts.max() ) ); } void GenericKTProcess::computeIncomingFluxes( double x1, double q1t2, double x2, double q2t2 ) { flux1_ = flux2_ = 0.; switch ( cuts_.mode ) { case Kinematics::ElasticElastic: flux1_ = elasticFlux( x1, q1t2 ); flux2_ = elasticFlux( x2, q2t2 ); break; case Kinematics::ElasticInelastic: flux1_ = elasticFlux( x1, q1t2 ); flux2_ = inelasticFlux( x2, q2t2, MY_ ); break; case Kinematics::InelasticElastic: flux1_ = inelasticFlux( x1, q1t2, MX_ ); flux2_ = elasticFlux( x2, q2t2 ); break; case Kinematics::InelasticInelastic: flux1_ = inelasticFlux( x1, q1t2, MX_ ); flux2_ = inelasticFlux( x2, q2t2, MY_ ); break; default: return; } flux1_ = std::max( flux1_, 1.e-20 ); flux2_ = std::max( flux2_, 1.e-20 ); DebuggingInsideLoop( Form( "Form factors: %e / %e", flux1_, flux2_ ) ); } void GenericKTProcess::fillKinematics( bool ) { fillPrimaryParticlesKinematics(); fillCentralParticlesKinematics(); } void GenericKTProcess::fillPrimaryParticlesKinematics() { //================================================================= // outgoing protons //================================================================= Particle& op1 = event_->getOneByRole( Particle::OutgoingBeam1 ), &op2 = event_->getOneByRole( Particle::OutgoingBeam2 ); op1.setMomentum( PX_ ); op2.setMomentum( PY_ ); switch ( cuts_.mode ) { case Kinematics::ElasticElastic: op1.setStatus( Particle::FinalState ); op2.setStatus( Particle::FinalState ); break; case Kinematics::ElasticInelastic: op1.setStatus( Particle::FinalState ); - op2.setStatus( Particle::Undecayed ); op2.setMass( MY_ ); + op2.setStatus( Particle::Unfragmented ); op2.setMass( MY_ ); break; case Kinematics::InelasticElastic: - op1.setStatus( Particle::Undecayed ); op1.setMass( MX_ ); + op1.setStatus( Particle::Unfragmented ); op1.setMass( MX_ ); op2.setStatus( Particle::FinalState ); break; case Kinematics::InelasticInelastic: - op1.setStatus( Particle::Undecayed ); op1.setMass( MX_ ); - op2.setStatus( Particle::Undecayed ); op2.setMass( MY_ ); + op1.setStatus( Particle::Unfragmented ); op1.setMass( MX_ ); + op2.setStatus( Particle::Unfragmented ); op2.setMass( MY_ ); break; default: { FatalError( "This kT factorisation process is intended for p-on-p collisions! Aborting!" ); } break; } //================================================================= // incoming partons (photons, pomerons, ...) //================================================================= //FIXME ensure the validity of this approach Particle& g1 = event_->getOneByRole( Particle::Parton1 ), &g2 = event_->getOneByRole( Particle::Parton2 ); g1.setMomentum( event_->getOneByRole( Particle::IncomingBeam1 ).momentum()-PX_, true ); g1.setStatus( Particle::Incoming ); g2.setMomentum( event_->getOneByRole( Particle::IncomingBeam2 ).momentum()-PY_, true ); g2.setStatus( Particle::Incoming ); } double GenericKTProcess::minimalJacobian() const { double jac = 1.; jac *= ( log_qmax_-log_qmin_ )*qt1_; // d(q1t) . q1t jac *= ( log_qmax_-log_qmin_ )*qt2_; // d(q2t) . q2t jac *= 2.*M_PI; // d(phi1) jac *= 2.*M_PI; // d(phi2) const double mx_range = cuts_.cuts.remnants.at( Cuts::mass ).range(); switch ( cuts_.mode ) { case Kinematics::ElasticElastic: default: break; case Kinematics::ElasticInelastic: jac *= 2.* mx_range * MY_; break; case Kinematics::InelasticElastic: jac *= 2.* mx_range * MX_; break; case Kinematics::InelasticInelastic: jac *= 2.* mx_range * MX_; jac *= 2.* mx_range * MY_; break; } // d(mx/y**2) return jac; } double GenericKTProcess::elasticFlux( double x, double kt2 ) { const double mp = Particle::massFromPDGId( Particle::Proton ), mp2 = mp*mp; const double Q2_min = x*x*mp2/( 1.-x ), Q2_ela = Q2_min + kt2/( 1.-x ); const FormFactors ela = FormFactors::ProtonElastic( Q2_ela ); const double ela1 = ( 1.-x )*( 1.-Q2_min/Q2_ela ); const double ela2 = ela.FE; //const double ela3 = 1.-( Q2_ela-kt2 )/Q2_ela; const double f_ela = Constants::alphaEM/M_PI/kt2*( ela1*ela2 + 0.5*x*x*ela.FM ); return f_ela * ( 1.-x ) * kt2 / Q2_ela; } double GenericKTProcess::inelasticFlux( double x, double kt2, double mx ) { const double mx2 = mx*mx, mp = Particle::massFromPDGId( Particle::Proton ), mp2 = mp*mp; // F2 structure function const double Q2min = 1. / ( 1.-x )*( x*( mx2-mp2 ) + x*x*mp2 ), Q2 = Q2min + kt2 / ( 1.-x ); float xbj = Q2 / ( Q2 + mx2 - mp2 ); const StructureFunctions sf = StructureFunctionsBuilder::get( cuts_.structure_functions, Q2, xbj ); // Longitudinal/transverse virtual photon cross section R // from Sibirtsev and Blunden (Phys Rev C 88,065202 (2013)) const double R = 0.014 * Q2 * ( exp( -0.07*Q2 ) + 41.*exp( -0.8*Q2 ) ); const double F1 = sf.F2 * ( 1.+4*xbj*xbj*mp2/Q2 )/( 1.+R )/( 2.*xbj ); const double ine1 = ( 1.-x )*( 1.-Q2min / Q2 ); const double f_D = sf.F2/( Q2 + mx2 - mp2 ) * ine1, f_C= 2.*F1/Q2; const double f_ine = Constants::alphaEM/M_PI/kt2*( f_D + 0.5*x*x*f_C ); return f_ine * ( 1.-x ) * kt2 / Q2; } } } diff --git a/CepGen/Processes/GenericProcess.cpp b/CepGen/Processes/GenericProcess.cpp index 686f190..b3a988c 100644 --- a/CepGen/Processes/GenericProcess.cpp +++ b/CepGen/Processes/GenericProcess.cpp @@ -1,181 +1,187 @@ #include "GenericProcess.h" #include "CepGen/Core/Exception.h" namespace CepGen { namespace Process { GenericProcess::GenericProcess( const std::string& name, const std::string& description, bool has_event ) : event_( std::shared_ptr<Event>( new Event ) ), is_point_set_( false ), is_incoming_state_set_( false ), is_outgoing_state_set_( false ), is_kinematics_set_( false ), name_( name ), description_( description ), total_gen_time_( 0. ), num_gen_events_( 0 ), has_event_( has_event ) {} GenericProcess::~GenericProcess() {} void GenericProcess::setPoint( const unsigned int ndim, double* x ) { if ( ndim != x_.size() ) x_.resize( ndim ); x_ = std::vector<double>( x, x+ndim ); is_point_set_ = true; if ( Logger::get().level>=Logger::DebugInsideLoop ) { dumpPoint(); } } void GenericProcess::prepareKinematics() { if ( !isKinematicsDefined() ) return; // FIXME dump some information... const Particle ib1 = event_->getOneByRole( Particle::IncomingBeam1 ), ib2 = event_->getOneByRole( Particle::IncomingBeam2 ); sqs_ = CMEnergy( ib1, ib2 ); s_ = sqs_*sqs_; w1_ = ib1.mass2(); w2_ = ib2.mass2(); Debugging( Form( "Kinematics successfully prepared! sqrt(s) = %.2f", sqs_ ) ); } void GenericProcess::dumpPoint() { std::ostringstream os; for ( unsigned int i=0; i<x_.size(); i++ ) { os << Form( " x(%2d) = %8.6f\n\t", i, x_[i] ); } Information( Form( "Number of integration parameters: %d\n\t" "%s", x_.size(), os.str().c_str() ) ); } void GenericProcess::setEventContent( const IncomingState& is, const OutgoingState& os ) { event_->clear(); //----- add the particles in the event //--- incoming state for ( IncomingState::const_iterator ip=is.begin(); ip!=is.end(); ip++ ) { - event_->addParticle( Particle( ip->first, ip->second ) ); + std::cout << ip->first << "\t" << ip->second << std::endl; + Particle& p = event_->addParticle( ip->first ); + p.setPdgId( ip->second ); + p.dump(); } //--- central system (if not already there) IncomingState::const_iterator central_system = is.find( Particle::CentralSystem ); if ( central_system == is.end() ) { - event_->addParticle( Particle( Particle::Intermediate, Particle::invalidParticle, Particle::Propagator ) ); + Particle& p = event_->addParticle( Particle::Intermediate ); + p.setPdgId( Particle::invalidParticle ); + p.setStatus( Particle::Propagator ); } //--- outgoing state for ( OutgoingState::const_iterator op = os.begin(); op != os.end(); ++op ) { for ( std::vector<Particle::ParticleCode>::const_iterator it = op->second.begin(); it != op->second.end(); ++it ) { - event_->addParticle( Particle( op->first, *it ) ); + Particle& p = event_->addParticle( op->first ); + p.setPdgId( *it ); } } //----- define the particles parentage const Particles parts = event_->particles(); for ( Particles::const_iterator p = parts.begin(); p != parts.end(); ++p ) { Particle& part = event_->getById( p->id() ); switch ( part.role() ) { case Particle::OutgoingBeam1: case Particle::Parton1: part.addMother( event_->getOneByRole( Particle::IncomingBeam1 ) ); break; case Particle::OutgoingBeam2: case Particle::Parton2: part.addMother( event_->getOneByRole( Particle::IncomingBeam2 ) ); break; case Particle::Intermediate: part.addMother( event_->getOneByRole( Particle::Parton1 ) ); part.addMother( event_->getOneByRole( Particle::Parton2 ) ); break; case Particle::CentralSystem: part.addMother( event_->getOneByRole( Particle::Intermediate ) ); break; default: break; } } //----- freeze the event as it is - event_->init(); + event_->freeze(); } void GenericProcess::setIncomingKinematics( const Particle::Momentum& p1, const Particle::Momentum& p2 ) { event_->getOneByRole( Particle::IncomingBeam1 ).setMomentum( p1 ); event_->getOneByRole( Particle::IncomingBeam2 ).setMomentum( p2 ); } void GenericProcess::formFactors( double q1, double q2, FormFactors& fp1, FormFactors& fp2 ) const { const double mx2 = MX_*MX_, my2 = MY_*MY_; switch ( cuts_.mode ) { case Kinematics::ElectronElectron: { fp1 = FormFactors::Trivial(); // electron (trivial) form factor fp2 = FormFactors::Trivial(); // electron (trivial) form factor } break; case Kinematics::ProtonElectron: { fp1 = FormFactors::ProtonElastic( -t1_ ); // proton elastic form factor fp2 = FormFactors::Trivial(); // electron (trivial) form factor } break; case Kinematics::ElectronProton: { fp1 = FormFactors::Trivial(); // electron (trivial) form factor fp2 = FormFactors::ProtonElastic( -t2_ ); // proton elastic form factor } break; case Kinematics::ElasticElastic: { fp1 = FormFactors::ProtonElastic( -t1_ ); // proton elastic form factor fp2 = FormFactors::ProtonElastic( -t2_ ); // proton elastic form factor } break; case Kinematics::ElasticInelastic: { fp1 = FormFactors::ProtonElastic( -t1_ ); fp2 = FormFactors::ProtonInelastic( cuts_.structure_functions, -t2_, w2_, my2 ); } break; case Kinematics::InelasticElastic: { fp1 = FormFactors::ProtonInelastic( cuts_.structure_functions, -t1_, w1_, mx2 ); fp2 = FormFactors::ProtonElastic( -t2_ ); } break; case Kinematics::InelasticInelastic: { fp1 = FormFactors::ProtonInelastic( cuts_.structure_functions, -t1_, w1_, mx2 ); fp2 = FormFactors::ProtonInelastic( cuts_.structure_functions, -t2_, w2_, my2 ); } break; } } inline bool GenericProcess::isKinematicsDefined() { // check the incoming state if ( !particles( Particle::IncomingBeam1 ).empty() && !particles( Particle::IncomingBeam2 ).empty() ) { is_incoming_state_set_ = true; } // check the outgoing state if ( !particles( Particle::OutgoingBeam1 ).empty() && !particles( Particle::OutgoingBeam2 ).empty() && !particles( Particle::CentralSystem ).empty() ) { is_outgoing_state_set_ = true; } // combine both states is_kinematics_set_ = is_incoming_state_set_ && is_outgoing_state_set_; return is_kinematics_set_; } std::ostream& operator<<( std::ostream& os, const GenericProcess& proc ) { return os << proc.name().c_str(); } std::ostream& operator<<( std::ostream& os, const GenericProcess* proc ) { return os << proc->name().c_str(); } } }