Jpp/v12.1.0/JStingray_8cc_source.html

#include <string>

#include <iostream>

#include <iomanip>

#include <map>


#include "TRandom3.h"


#include "evt/Head.hh"

#include "evt/Evt.hh"


#include "JSupport/JMultipleFileScanner.hh"

#include "JSupport/JMonteCarloFileSupportkit.hh"

#include "JSupport/JFileRecorder.hh"

#include "JSupport/JSupport.hh"

#include "JSupport/JMeta.hh"


#include "JPhysics/JPDFToolkit.hh"

#include "JTools/JConstants.hh"


#include "Jeep/JPrint.hh"

#include "Jeep/JParser.hh"

#include "Jeep/JMessage.hh"


/**

 * \file

 *

 * Auxiliary program to manipulate MUPAGE data.

 * \author mdejong

 */

int main(int argc, char **argv)

{

  using namespace std;

  using namespace JPP;


  JMultipleFileScanner<Evt>                                              inputFile;

  JFileRecorder<JTYPELIST<JAAnetTypes_t, JMeta, JRootTypes_t>::typelist> outputFile;

  JLimit_t&      numberOfEvents = inputFile.getLimit();

  double         P;

  double         u;

  double         E;

  UInt_t         seed;

  int            debug;


  try {


    JParser<> zap("Auxiliary program to manipulate MUPAGE data.");


    zap['f'] = make_field(inputFile);

    zap['o'] = make_field(outputFile)      = "mupage.root";

    zap['n'] = make_field(numberOfEvents)  = JLimit::max();

    zap['P'] = make_field(P,  "Keep muons in event with given probability  (<= 1, 1 == all)")  = 1.0;

    zap['u'] = make_field(u,  "Reject events as a function of multiplicity (>= 0, 0 == none)") = 0.0;

    zap['E'] = make_field(E,  "Energy scaling factor")                                         = 1.0;

    zap['S'] = make_field(seed)            = 0;

    zap['d'] = make_field(debug)           = 2;


    zap(argc, argv);

  }

  catch(const exception &error) {

    FATAL(error.what() << endl);

  }


  gRandom->SetSeed(seed);


  outputFile.open();


  if (!outputFile.is_open()) {

    FATAL("Error opening file " << outputFile << endl);

  }


  outputFile.put(JMeta(argc, argv));

  outputFile.put(*gRandom);


  JMultipleFileScanner<Head> io(inputFile);


  io >> outputFile;


  while (inputFile.hasNext()) {


    STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);


    Evt* event = inputFile.next();


    // randomly reject muons in event


    for (vector<Trk>::iterator i = event->mc_trks.begin(); i != event->mc_trks.end(); ) {


      if (gRandom->Rndm() < P)

        ++i;

      else

        i = event->mc_trks.erase(i);

    }


    // apply energy loss


    for (vector<Trk>::iterator i = event->mc_trks.begin(); i != event->mc_trks.end(); ) {


      i->E  = gWater(i->E, i->t * getSpeedOfLight());

      i->E *= E;


      if (i->E > MASS_MUON)

        ++i;

      else

        i = event->mc_trks.erase(i);

    }


    // randonly reject events as a function of muon multiplicity


    if (gRandom->Rndm() < pow(event->mc_trks.size(),-u)) {

      outputFile.put(*event);

    }

  }

  STATUS(endl);


  outputFile.put(*gRandom);

  outputFile.close();

}