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    origin(0,0,0);
  bool   genhen = false;
  bool   mupage = false;
  
  {
    JHead buffer(header);
 
    origin = getOrigin(buffer);
    genhen = is_genhen(buffer);
    mupage = is_mupage(buffer);
    
    if        (genhen && buffer.is_valid(&JHead::genvol)) {
 
      buffer.genvol.mul(scale, origin.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 -= origin;
        track->pos *= scale;
        track->pos += origin;
      }
    }
 
    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 - origin) * (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();
}