Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-18 10:06:36 

0001 // aMCatNLOHooks.h is a part of the PYTHIA event generator.
0002 // Copyright (C) 2024 Torbjorn Sjostrand.
0003 // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
0004 // Please respect the MCnet Guidelines, see GUIDELINES for details.
0005 
0006 // This program is written by Stefan Prestel.
0007 // It illustrates how to do run PYTHIA with LHEF input, allowing a
0008 // sample-by-sample generation of
0009 // a) Non-matched/non-merged events
0010 // b) MLM jet-matched events (kT-MLM, shower-kT, FxFx)
0011 // c) CKKW-L and UMEPS-merged events
0012 // d) UNLOPS NLO merged events
0013 // see the respective sections in the online manual for details.
0014 
0015 #ifndef Pythia8_aMCatNLOHooks_H
0016 #define Pythia8_aMCatNLOHooks_H
0017 
0018 #include "Pythia8/Pythia.h"
0019 
0020 namespace Pythia8 {
0021 
0022 //==========================================================================
0023 
0024 // Use userhooks to set the number of requested partons dynamically, as
0025 // needed when running CKKW-L or UMEPS on a single input file that contains
0026 // all parton multiplicities.
0027 
0028 class amcnlo_unitarised_interface : public UserHooks {
0029 
0030 public:
0031 
0032   // Constructor and destructor.
0033   amcnlo_unitarised_interface() : mergingScheme(0), normFactor(1.) {}
0034   amcnlo_unitarised_interface(int mergingSchemeIn)
0035     : mergingScheme(mergingSchemeIn), normFactor(1.) {}
0036  ~amcnlo_unitarised_interface() {}
0037 
0038   double getNormFactor(){return normFactor;}
0039 
0040   // Allow to set the number of partons.
0041   bool canVetoProcessLevel() { return true;}
0042   // Set the number of partons.
0043   bool doVetoProcessLevel( Event& process) {
0044 
0045     int nPartons = 0;
0046     normFactor = 1.;
0047 
0048     // Do not include resonance decay products in the counting.
0049     omitResonanceDecays(process);
0050 
0051     // Get the maximal quark flavour counted as "additional" parton.
0052     int nQuarksMerge = settingsPtr->mode("Merging:nQuarksMerge");
0053 
0054     // Save information on the process string.
0055     string procSave = settingsPtr->word("Merging:Process");
0056     bool hasIN = (procSave.find(">", 0) != string::npos);
0057     string incoming = procSave.substr(0,procSave.find(">", 0));
0058     string outgoing = procSave.substr(procSave.find(">", 0)+1,procSave.size());
0059     string outgoingSave = outgoing;
0060 
0061     // Count number of user particles.
0062     int nCommas = 0;
0063     for(int n = outgoing.find(",", 0); n != int(string::npos);
0064         n = outgoing.find(",", n)) { nCommas++; n++; }
0065     vector <string> out;
0066     for(int i =0; i < nCommas;++i) {
0067       int n = outgoing.find(",", 0);
0068       out.push_back(outgoing.substr(0,n));
0069       outgoing = outgoing.substr(n+1,outgoing.size());
0070     }
0071     if (outgoing.size()>0) out.push_back(outgoing);
0072 
0073     // Dynamically set the process string.
0074     string proc;
0075     if ( settingsPtr->word("Merging:Process") == "guess" ) {
0076       string processString = "";
0077       // Set incoming particles.
0078       if (hasIN) {
0079         processString += incoming + ">";
0080       } else {
0081         int beamAid = beamAPtr->id();
0082         int beamBid = beamBPtr->id();
0083         if (abs(beamAid) == 2212) processString += "p";
0084         if (beamAid == 11)        processString += "e-";
0085         if (beamAid ==-11)        processString += "e+";
0086         if (abs(beamBid) == 2212) processString += "p";
0087         if (beamBid == 11)        processString += "e-";
0088         if (beamBid ==-11)        processString += "e+";
0089         processString += ">";
0090       }
0091       // Set outgoing particles.
0092       bool foundOutgoing = false;
0093       for(int i=0; i < int(workEvent.size()); ++i) {
0094         if ( workEvent[i].isFinal()
0095           && ( workEvent[i].colType() == 0
0096             || workEvent[i].idAbs() > 21
0097             || ( workEvent[i].id() != 21
0098               && workEvent[i].idAbs() > nQuarksMerge) ) ) {
0099           foundOutgoing = true;
0100           ostringstream procOSS;
0101           procOSS << "{" << workEvent[i].name() << "," << workEvent[i].id()
0102                   << "}";
0103           processString += procOSS.str();
0104         }
0105       }
0106 
0107       for (int i =0; i < int(out.size()); ++i) {
0108         if (out[i].find("guess") != string::npos) continue;
0109         processString += "," + out[i];
0110       }
0111 
0112       if (foundOutgoing) proc = processString;
0113     }
0114 
0115     if (proc.size()>0) settingsPtr->word("Merging:Process", proc);
0116     else proc = procSave;
0117 
0118     // Loop through event and count.
0119     for(int i=0; i < int(workEvent.size()); ++i)
0120       if ( workEvent[i].isFinal()
0121         && workEvent[i].colType()!= 0
0122         && ( workEvent[i].id() == 21 || workEvent[i].idAbs() <= nQuarksMerge))
0123         nPartons++;
0124 
0125     // Store merging scheme.
0126     bool isumeps  = (mergingScheme == 1);
0127     bool isunlops = (mergingScheme == 2);
0128 
0129     // Get number of requested partons.
0130     string nps_nlo = infoPtr->getEventAttribute("npNLO",true);
0131     int np_nlo     = (nps_nlo != "") ? atoi((char*)nps_nlo.c_str()) : -1;
0132     string nps_lo  = infoPtr->getEventAttribute("npLO",true);
0133     int np_lo      = (nps_lo != "") ? atoi((char*)nps_lo.c_str()) : -1;
0134 
0135     if ( (settingsPtr->flag("Merging:doUNLOPSTree")
0136        || settingsPtr->flag("Merging:doUNLOPSSubt")) && np_lo == 0)
0137        return true;
0138 
0139     int nj = 0;
0140     for(int n = proc.find("j", 0); n != int(string::npos);
0141         n = proc.find("j", n)) { nj++; n++; }
0142     nPartons -= nj;
0143     if (settingsPtr->word("Merging:process").compare("e+e->jj") == 0) {
0144       np_lo -= 2;
0145       np_nlo -= 2;
0146     }
0147 
0148     // Set number of requested partons.
0149     if (np_nlo > -1){
0150       settingsPtr->mode("Merging:nRequested", np_nlo);
0151       np_lo = -1;
0152     } else if (np_lo > -1){
0153       settingsPtr->mode("Merging:nRequested", np_lo);
0154       np_nlo = -1;
0155     } else {
0156       settingsPtr->mode("Merging:nRequested", nPartons);
0157       np_nlo = -1;
0158       np_lo = nPartons;
0159     }
0160 
0161     // Reset the event weight to incorporate corrective factor.
0162     bool updateWgt = true;
0163 
0164     // Choose randomly if this event should be treated as subtraction or
0165     // as regular event. Put the correct settings accordingly.
0166     if (isunlops && np_nlo == 0 && np_lo == -1) {
0167       settingsPtr->flag("Merging:doUNLOPSTree", false);
0168       settingsPtr->flag("Merging:doUNLOPSSubt", false);
0169       settingsPtr->flag("Merging:doUNLOPSLoop", true);
0170       settingsPtr->flag("Merging:doUNLOPSSubtNLO", false);
0171       settingsPtr->mode("Merging:nRecluster",0);
0172       normFactor *= 1.;
0173     } else if (isunlops && np_nlo > 0 && np_lo == -1) {
0174       normFactor *= 2.;
0175       if (rndmPtr->flat() < 0.5) {
0176         normFactor *= -1.;
0177         settingsPtr->flag("Merging:doUNLOPSTree", false);
0178         settingsPtr->flag("Merging:doUNLOPSSubt", false);
0179         settingsPtr->flag("Merging:doUNLOPSLoop", false);
0180         settingsPtr->flag("Merging:doUNLOPSSubtNLO", true);
0181         settingsPtr->mode("Merging:nRecluster",1);
0182       } else {
0183         settingsPtr->flag("Merging:doUNLOPSTree", false);
0184         settingsPtr->flag("Merging:doUNLOPSSubt", false);
0185         settingsPtr->flag("Merging:doUNLOPSLoop", true);
0186         settingsPtr->flag("Merging:doUNLOPSSubtNLO", false);
0187         settingsPtr->mode("Merging:nRecluster",0);
0188       }
0189     } else if (isunlops && np_nlo == -1 && np_lo > -1) {
0190       normFactor *= 2.;
0191       if (rndmPtr->flat() < 0.5) {
0192         normFactor *= -1.;
0193         settingsPtr->flag("Merging:doUNLOPSTree", false);
0194         settingsPtr->flag("Merging:doUNLOPSSubt", true);
0195         settingsPtr->flag("Merging:doUNLOPSLoop", false);
0196         settingsPtr->flag("Merging:doUNLOPSSubtNLO", false);
0197         settingsPtr->mode("Merging:nRecluster",1);
0198 
0199         // Double reclustering for exclusive NLO cross sections.
0200         bool isnlotilde = settingsPtr->flag("Merging:doUNLOPSTilde");
0201         int nmaxNLO = settingsPtr->mode("Merging:nJetMaxNLO");
0202         if ( isnlotilde
0203           && nmaxNLO > 0
0204           && np_lo <= nmaxNLO+1
0205           && np_lo > 1 ){
0206           normFactor *= 2.;
0207           if (rndmPtr->flat() < 0.5)
0208             settingsPtr->mode("Merging:nRecluster",2);
0209         }
0210       } else {
0211         settingsPtr->flag("Merging:doUNLOPSTree", true);
0212         settingsPtr->flag("Merging:doUNLOPSSubt", false);
0213         settingsPtr->flag("Merging:doUNLOPSLoop", false);
0214         settingsPtr->flag("Merging:doUNLOPSSubtNLO", false);
0215         settingsPtr->mode("Merging:nRecluster",0);
0216       }
0217     } else if (isumeps && np_lo == 0) {
0218       settingsPtr->flag("Merging:doUMEPSTree", true);
0219       settingsPtr->flag("Merging:doUMEPSSubt", false);
0220       settingsPtr->mode("Merging:nRecluster",0);
0221     } else if (isumeps && np_lo > 0) {
0222       normFactor *= 2.;
0223       if (rndmPtr->flat() < 0.5) {
0224         normFactor *= -1.;
0225         settingsPtr->flag("Merging:doUMEPSTree", false);
0226         settingsPtr->flag("Merging:doUMEPSSubt", true);
0227         settingsPtr->mode("Merging:nRecluster",1);
0228       } else {
0229         settingsPtr->flag("Merging:doUMEPSTree", true);
0230         settingsPtr->flag("Merging:doUMEPSSubt", false);
0231         settingsPtr->mode("Merging:nRecluster",0);
0232       }
0233     }
0234     // Reset the event weight to incorporate corrective factor.
0235     if ( updateWgt) {
0236       infoPtr->weightContainerPtr->weightNominal *= normFactor;
0237       normFactor = 1.;
0238     }
0239 
0240     // Done
0241     return false;
0242   }
0243 
0244 private:
0245 
0246   int mergingScheme;
0247   double normFactor;
0248 };
0249 
0250 //==========================================================================
0251 
0252 } // end namespace Pythia8
0253 
0254 #endif // end Pythia8_aMCatNLOHooks_H
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!