PYTHIA8/pythia8175/examples/main72.cc

   1 // main72.cc is a part of the PYTHIA event generator.
   2 // Copyright (C) 2013 Torbjorn Sjostrand.
   3 // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
   4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
   5
   6 // This is a simple test program.
   7 // It compares SlowJet and FastJet, showing that they find the same jets.
   8
   9 #include "Pythia.h"
  10
  11 // The FastJet3.h header enables automatic initialisation of
  12 // fastjet::PseudoJet objects from Pythia8 Particle and Vec4 objects,
  13 // as well as advanced features such as access to (a copy of)
  14 // the original Pythia 8 Particle directly from the PseudoJet,
  15 // and fastjet selectors that make use of the Particle properties.
  16 // See the extensive comments in the header file for further details
  17 // and examples.
  18 #include "FastJet3.h"
  19
  20 using namespace Pythia8;
  21
  22 int main() {
  23
  24   // Number of events, generated and listed ones (for jets).
  25   int nEvent    = 1000;
  26   int nListJets = 3;
  27
  28   // Select common parameters for SlowJet and FastJet analyses.
  29   int    power   = -1;     // -1 = anti-kT; 0 = C/A; 1 = kT.
  30   double R       = 0.7;    // Jet size.
  31   double pTMin   = 5.0;    // Min jet pT.
  32   double etaMax  = 5.0;    // Pseudorapidity range of detector.
  33   int    select  = 2;      // Which particles are included?
  34   int    massSet = 2;      // Which mass are they assumed to have?
  35
  36   // Generator. Shorthand for event.
  37   Pythia pythia;
  38   Event& event = pythia.event;
  39
  40   // Process selection.
  41   pythia.readString("HardQCD:all = on");
  42   pythia.readString("PhaseSpace:pTHatMin = 200.");
  43
  44   // No event record printout.
  45   pythia.readString("Next:numberShowInfo = 0");
  46   pythia.readString("Next:numberShowProcess = 0");
  47   pythia.readString("Next:numberShowEvent = 0");
  48
  49   // LHC initialization.
  50   pythia.readString("Beams:eCM = 14000.");
  51   pythia.init();
  52
  53   // Set up SlowJet jet finder.
  54   SlowJet slowJet( power, R, pTMin, etaMax, select, massSet);
  55
  56   // Set up FastJet jet finder.
  57   //   one can use either explicitly use antikt, cambridge, etc., or
  58   //   just use genkt_algorithm with specification of power
  59   //fastjet::JetAlgorithm algorithm;
  60   //if (power == -1)      algorithm = fastjet::antikt_algorithm;
  61   //if (power ==  0)      algorithm = fastjet::cambridge_algorithm;
  62   //if (power ==  1)      algorithm = fastjet::kt_algorithm;
  63   //fastjet::JetDefinition jetDef(algorithm, R);
  64   // there's no need for a pointer to the jetDef (it's a fairly small object)
  65   fastjet::JetDefinition jetDef(fastjet::genkt_algorithm, R, power);
  66   std::vector <fastjet::PseudoJet> fjInputs;
  67
  68   // Histograms.
  69   Hist nJetsS("number of jets, SlowJet", 100, -0.5, 99.5);
  70   Hist nJetsF("number of jets, FastJet", 100, -0.5, 99.5);
  71   Hist nJetsD("number of jets, SlowJet-FastJet", 99, -49.5, 49.5);
  72   Hist pTjetsS("pT for jets, SlowJet", 100, 0., 500.);
  73   Hist pTjetsF("pT for jets, FastJet", 100, 0., 500.);
  74   Hist pTjetsD("pT for jets, SlowJet - FastJet", 100, 0., 500.);
  75   Hist RdistD("R distance SlowJet to FastJet", 100, 0., 1.);
  76   Hist etaJets("eta for jets", 100, -5., 5.);
  77   Hist phiJets("phi for jets", 100, -M_PI, M_PI);
  78   Hist distJets("R distance between jets", 100, 0., 10.);
  79   Hist pTdiff("pT difference between consecutive jets", 100, -100., 400.);
  80   Hist nAna("multiplicity of analyzed event", 100, -0.5, 999.5);
  81   Hist tGen("generation time as fn of multiplicity", 100, -0.5, 999.5);
  82   Hist tSlow("SlowJet time as fn of multiplicity", 100, -0.5, 999.5);
  83   Hist tFast("FastJet time as fn of multiplicity", 100, -0.5, 999.5);
  84   Hist tSlowGen("SlowJet/generation time as fn of multiplicity",
  85     100, -0.5, 999.5);
  86   Hist tFastGen("FastJet/generation time as fn of multiplicity",
  87     100, -0.5, 999.5);
  88
  89   // Begin event loop. Generate event. Skip if error.
  90   for (int iEvent = 0; iEvent < nEvent; ++iEvent) {
  91     clock_t befGen = clock();
  92     if (!pythia.next()) continue;
  93     clock_t aftGen = clock();
  94
  95     // Begin SlowJet analysis of jet properties. List first few.
  96     clock_t befSlow = clock();
  97     slowJet.analyze( pythia.event );
  98     clock_t aftSlow = clock();
  99     if (iEvent < nListJets) slowJet.list();
 100
 101     // Fill inclusive SlowJet jet distributions.
 102     int nSlow = slowJet.sizeJet();
 103     nJetsS.fill( nSlow );
 104     for (int i = 0; i < nSlow; ++i) {
 105       pTjetsS.fill( slowJet.pT(i) );
 106       etaJets.fill( slowJet.y(i) );
 107       phiJets.fill( slowJet.phi(i) );
 108     }
 109
 110     // Fill distance between SlowJet jets.
 111     for (int i = 0; i < nSlow - 1; ++i)
 112     for (int j = i + 1; j < nSlow; ++j) {
 113       double dY = slowJet.y(i)  - slowJet.y(j);
 114       double dPhi = abs( slowJet.phi(i) - slowJet.phi(j) );
 115       if (dPhi > M_PI) dPhi = 2. * M_PI - dPhi;
 116       double dR = sqrt( pow2(dY) + pow2(dPhi) );
 117       distJets.fill( dR );
 118     }
 119
 120     // Fill pT-difference between SlowJet jets (to check ordering of list).
 121     for (int i = 1; i < nSlow; ++i)
 122       pTdiff.fill( slowJet.pT(i-1)- slowJet.pT(i) );
 123
 124     // Begin FastJet analysis: extract particles from event record.
 125     clock_t befFast = clock();
 126     fjInputs.resize(0);
 127     Vec4   pTemp;
 128     double mTemp;
 129     int nAnalyze = 0;
 130     for (int i = 0; i < event.size(); ++i) if (event[i].isFinal()) {
 131
 132       // Require visible/charged particles inside detector.
 133       if      (select > 2 &&  event[i].isNeutral() ) continue;
 134       else if (select == 2 && !event[i].isVisible() ) continue;
 135       if (etaMax < 20. && abs(event[i].eta()) > etaMax) continue;
 136
 137       // Create a PseudoJet from the complete Pythia particle.
 138       fastjet::PseudoJet particleTemp = event[i];
 139
 140       // Optionally modify mass and energy.
 141       pTemp = event[i].p();
 142       mTemp = event[i].m();
 143       if (massSet < 2) {
 144         mTemp = (massSet == 0 || event[i].id() == 22) ? 0. : 0.13957;
 145         pTemp.e( sqrt(pTemp.pAbs2() + mTemp*mTemp) );
 146         particleTemp.reset_momentum( pTemp.px(), pTemp.py(),
 147            pTemp.pz(), pTemp.e() );
 148       }
 149
 150       // Store acceptable particles as input to Fastjet.
 151       // Conversion to PseudoJet is performed automatically
 152       // with the help of the code in FastJet3.h.
 153       fjInputs.push_back( particleTemp);
 154       ++nAnalyze;
 155     }
 156
 157     // Run Fastjet algorithm and sort jets in pT order.
 158     vector <fastjet::PseudoJet> inclusiveJets, sortedJets;
 159     fastjet::ClusterSequence clustSeq(fjInputs, jetDef);
 160     inclusiveJets = clustSeq.inclusive_jets(pTMin);
 161     sortedJets    = sorted_by_pt(inclusiveJets);
 162     clock_t aftFast = clock();
 163
 164     // List first few FastJet jets and some info about them.
 165     // Note: the final few columns are illustrative of what information
 166     // can be extracted, but does not exhaust the possibilities.
 167     if (iEvent < nListJets) {
 168       cout << "\n --------  FastJet jets, p = " << setw(2) << power
 169            << "  --------------------------------------------------\n\n "
 170            << "  i         pT        y      phi  mult chgmult photons"
 171            << "      hardest  pT in neutral " << endl
 172            << "                                                       "
 173            << "  constituent        hadrons " << endl;
 174       for (int i = 0; i < int(sortedJets.size()); ++i) {
 175         vector<fastjet::PseudoJet> constituents
 176           = sortedJets[i].constituents();
 177         fastjet::PseudoJet hardest
 178           = fastjet::SelectorNHardest(1)(constituents)[0];
 179         vector<fastjet::PseudoJet> neutral_hadrons
 180           = ( fastjet::SelectorIsHadron()
 181            && fastjet::SelectorIsNeutral())(constituents);
 182         double neutral_hadrons_pt = join(neutral_hadrons).perp();
 183         cout << setw(4) << i << fixed << setprecision(3) << setw(11)
 184              << sortedJets[i].perp() << setw(9)  << sortedJets[i].rap()
 185              << setw(9) << sortedJets[i].phi_std()
 186              << setw(6) << constituents.size()
 187              << setw(8) << fastjet::SelectorIsCharged().count(constituents)
 188              << setw(8) << fastjet::SelectorId(22).count(constituents)
 189              << setw(13) << hardest.user_info<Particle>().name()
 190              << "     " << setw(10) << neutral_hadrons_pt << endl;
 191       }
 192       cout << "\n --------  End FastJet Listing  ------------------"
 193            << "---------------------------------" << endl;
 194     }
 195
 196     // Fill inclusive FastJet jet distributions.
 197     int nFast         = sortedJets.size();
 198     nJetsF.fill( nFast );
 199     for (int i = 0; i < int(sortedJets.size()); ++i)
 200       pTjetsF.fill( sortedJets[i].perp() );
 201
 202     // Comparison of SlowJet and FastJet.
 203     nJetsD.fill( nSlow - nFast);
 204     if (nFast == nSlow) {
 205       for (int i = 0; i < nSlow; ++i) {
 206         double dist2 = pow2( slowJet.y(i) - sortedJets[i].rap())
 207           + pow2( slowJet.phi(i) - sortedJets[i].phi_std());
 208         RdistD.fill( sqrt(dist2) );
 209       }
 210     }
 211
 212     // Comparison of time consumption by analyzed multiplicity.
 213     nAna.fill( nAnalyze);
 214     tGen.fill( nAnalyze, aftGen - befGen);
 215     tSlow.fill( nAnalyze, aftSlow - befSlow);
 216     tFast.fill( nAnalyze, aftFast - befFast);
 217
 218   // End of event loop.
 219   }
 220
 221   // Statistics. Histograms.
 222   pythia.stat();
 223   pTjetsD = pTjetsS - pTjetsF;
 224   tSlowGen = tSlow / tGen;
 225   tFastGen = tFast / tGen;
 226   tGen  /= nAna;
 227   tSlow /= nAna;
 228   tFast /= nAna;
 229
 230   cout << nJetsS << nJetsF << nJetsD << pTjetsS << pTjetsF << pTjetsD
 231        << RdistD << etaJets << phiJets << distJets << pTdiff
 232        << nAna << tGen << tSlow << tFast << tSlowGen << tFastGen;
 233
 234   // Done.
 235   return 0;
 236 }