JSwordfish.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>

#include "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/MultiHead.hh"
#include "km3net-dataformat/offline/Evt.hh" 

#include "JAAnet/JHead.hh"
#include "JAAnet/JHeadToolkit.hh"
#include "JAAnet/JAAnetToolkit.hh"

#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JFileRecorder.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"
#include "JSupport/JSupport.hh"

#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"


/**
 * \file
 *
 * Example program to resize coordinate system of Monte Carlo events.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  JMultipleFileScanner<Evt>    inputFile;
  JLimit_t&                    numberOfEvents = inputFile.getLimit();
  JFileRecorder<JAAnetTypes_t> outputFile;
  double         scale;
  int            debug;

  try { 

    JParser<> zap("Example program to resize coordinate system of Monte Carlo events.");
    
    zap['f'] = make_field(inputFile);
    zap['o'] = make_field(outputFile);
    zap['n'] = make_field(numberOfEvents)      = JLimit::max();
    zap['F'] = make_field(scale)               = 1.0;
    zap['d'] = make_field(debug)               = 2;
    
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }


  Head   header = getHeader(inputFile);
  Vec    center(0,0,0);
  bool   genhen = false;
  bool   mupage = false;
  
  {
    JHead buffer(header);

    center = get<Vec>(buffer);
    genhen = is_genhen(buffer);
    mupage = is_mupage(buffer);
    
    if        (genhen && buffer.is_valid(&JHead::genvol)) {

      buffer.genvol.mul(scale, center.z);

    } else if (mupage && buffer.is_valid(&JHead::livetime)) {

      buffer.livetime.mul(scale);

    } else {
      
      FATAL("Invalid generator." << endl);
    }
    
    copy(buffer, header);
  }


  outputFile.open();

  outputFile.put(header);

  while (inputFile.hasNext()) {

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

    Evt* event = inputFile.next();

    if (genhen) {
    
      for (vector<Trk>::iterator track = event->mc_trks.begin(); track != event->mc_trks.end(); ++track) {

        track->pos -= center;
        track->pos *= scale;
        track->pos += center;
      }
    }

    if (mupage) {

      double E    = 0.0;
      Vec    pos;
      
      for (vector<Trk>::const_iterator track = event->mc_trks.begin(); track != event->mc_trks.end(); ++track) {
        E   += track->E;
        pos += track->pos * track->E;
      }
      
      pos   /=  E;
      pos    = (pos - center) * (scale - 1.0);
      pos.z  =  0.0;
      
      for (vector<Trk>::iterator track = event->mc_trks.begin(); track != event->mc_trks.end(); ++track) {
        track->pos += pos;
      }
    }
    
    outputFile.put(*event);
  }
  STATUS(endl); 

  outputFile.close();
}