PYTHIA8/pythia8170/examples/main16.cc

   1 // main16.cc is a part of the PYTHIA event generator.
   2 // Copyright (C) 2012 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 illustrates (a) how to collect the analysis code in a separate class
   8 // and (b) how to provide the .cmnd filename on the command line
   9
  10 // Once you have linked the main program you can run it with a command line
  11 // ./main16.exe main16.cmnd > out16
  12
  13 #include "Pythia.h"
  14
  15 using namespace Pythia8;
  16
  17 //==========================================================================
  18
  19 // Put all your own analysis code in the myAnalysis class.
  20
  21 class MyAnalysis {
  22
  23 public:
  24
  25   // Constructor can be empty.
  26   MyAnalysis() {}
  27
  28   // Initialization actions.
  29   void init();
  30
  31   // Analysis of each new event.
  32   void analyze(Event& event);
  33
  34   // Show final results.
  35   void finish();
  36
  37 private:
  38
  39   // Declare variables and objects that span init - analyze - finish.
  40   int  nEvt;
  41   Hist yH, etaChg, mult;
  42
  43 };
  44
  45 //--------------------------------------------------------------------------
  46
  47 // The initialization code.
  48
  49 void MyAnalysis::init() {
  50
  51   // Initialize counter for number of events.
  52   nEvt = 0;
  53
  54   // Book histograms.
  55   yH.book("Higgs rapidity", 100, -10., 10.);
  56   etaChg.book("charged pseudorapidity", 100, -10., 10.);
  57   mult.book( "charged multiplicity", 100, -0.5, 799.5);
  58
  59 }
  60
  61 //--------------------------------------------------------------------------
  62
  63 // The event analysis code.
  64
  65 void MyAnalysis::analyze(Event& event) {
  66
  67   // Increase counter.
  68   ++nEvt;
  69
  70   // Find latest copy of Higgs and plot its rapidity.
  71   int iH = 0;
  72   for (int i = 0; i < event.size(); ++i)
  73     if (event[i].id() == 25) iH = i;
  74   yH.fill( event[iH].y() );
  75
  76   // Plot pseudorapidity distribution. Sum up charged multiplicity.
  77   int nChg = 0;
  78   for (int i = 0; i < event.size(); ++i)
  79   if (event[i].isFinal() && event[i].isCharged()) {
  80     etaChg.fill( event[i].eta() );
  81     ++nChg;
  82   }
  83   mult.fill( nChg );
  84
  85 }
  86
  87 //--------------------------------------------------------------------------
  88
  89 // The finishing code.
  90
  91 void MyAnalysis::finish() {
  92
  93   // Normalize histograms.
  94   double binFactor = 5. / nEvt;
  95   yH     *= binFactor;
  96   etaChg *= binFactor;
  97
  98   // Print histograms.
  99   cout << yH << etaChg << mult;
 100
 101 }
 102
 103 //==========================================================================
 104
 105 // You should not need to touch the main program: its actions are
 106 // determined by the .cmnd file and the rest belongs in MyAnalysis.
 107
 108 int main(int argc, char* argv[]) {
 109
 110   // Check that correct number of command-line arguments
 111   if (argc != 2) {
 112     cerr << " Unexpected number of command-line arguments. \n"
 113          << " You are expected to provide a file name and nothing else. \n"
 114          << " Program stopped! " << endl;
 115     return 1;
 116   }
 117
 118   // Check that the provided file name corresponds to an existing file.
 119   ifstream is(argv[1]);
 120   if (!is) {
 121     cerr << " Command-line file " << argv[1] << " was not found. \n"
 122          << " Program stopped! " << endl;
 123     return 1;
 124   }
 125
 126   // Confirm that external file will be used for settings..
 127   cout << " PYTHIA settings will be read from file " << argv[1] << endl;
 128
 129   // Declare generator. Read in commands from external file.
 130   Pythia pythia;
 131   pythia.readFile(argv[1]);
 132
 133   // Initialization.
 134   pythia.init();
 135
 136   // Declare user analysis class. Do initialization part of it.
 137   MyAnalysis myAnalysis;
 138   myAnalysis.init();
 139
 140   // Read in number of event and maximal number of aborts.
 141   int nEvent = pythia.mode("Main:numberOfEvents");
 142   int nAbort = pythia.mode("Main:timesAllowErrors");
 143
 144   // Begin event loop.
 145   int iAbort = 0;
 146   for (int iEvent = 0; iEvent < nEvent; ++iEvent) {
 147
 148     // Generate events. Quit if too many failures.
 149     if (!pythia.next()) {
 150       if (++iAbort < nAbort) continue;
 151       cout << " Event generation aborted prematurely, owing to error!\n";
 152       break;
 153     }
 154
 155     // User Analysis of current event.
 156     myAnalysis.analyze( pythia.event);
 157
 158   // End of event loop.
 159   }
 160
 161   // Final statistics.
 162   pythia.stat();
 163
 164   // User finishing.
 165   myAnalysis.finish();
 166
 167   // Done.
 168   return 0;
 169 }