source: svn/trunk/src/STDHEPConverter.cc@ 466

/***********************************************************************
**                                                                    **
**          /----------------------------------------------\          **
**         |  Delphes, a framework for the fast simulation  |         **
**         |       of a  generic collider experiment        |         **
**          \-------------  arXiv:0903.2225v1  ------------/          **
**                                                                    **
**                                                                    **
**   This package uses:                                               **
**   ------------------                                               **
**   ROOT: Nucl. Inst. & Meth. in Phys. Res. A389 (1997) 81-86        **
**   FastJet algorithm: Phys. Lett. B641 (2006) [hep-ph/0512210]      **
**   Hector: JINST 2:P09005 (2007) [physics.acc-ph:0707.1198v2]       **
**   FROG: [hep-ex/0901.2718v1]                                       **
**   HepMC: Comput. Phys. Commun.134 (2001) 41                        **
**                                                                    **
** ------------------------------------------------------------------ **
**                                                                    **
**   Main authors:                                                    **
**   -------------                                                    **
**                                                                    **
**                Severine Ovyn                Xavier Rouby           **
**          severine.ovyn@uclouvain.be      xavier.rouby@cern         **
**                                                                    **
**         Center for Particle Physics and Phenomenology (CP3)        **
**                Universite catholique de Louvain (UCL)              **
**                     Louvain-la-Neuve, Belgium                      **
**                                                                    **
**                      Copyright (C) 2008-2009,                      **
**                        All rights reserved.                        **
**                                                                    **
***********************************************************************/

#include <iostream>
#include <fstream>
#include <cmath>
#include "TLorentzVector.h"
#include "BlockClasses.h"
#include "SmearUtil.h"
#include "ExRootTreeWriter.h"
#include "ExRootTreeBranch.h"
#include "stdhep_mcfio.h"
#include "stdhep_declarations.h"
#include "STDHEPConverter.h"
#include "PdgParticle.h"

using namespace std;

//---------------------------------------------------------------------------

void STDHEPConverter::AnalyseEvent(ExRootTreeBranch *branch, const Long64_t eventNumber)
{
  TRootGenEvent *element;

  element = static_cast<TRootGenEvent*>(branch->NewEntry());

  element->Number = eventNumber;
}

//---------------------------------------------------------------------------

void STDHEPConverter::AnalyseParticles(ExRootTreeBranch *branch) 
{
  TRootC::GenParticle *element;

  Double_t signPz;
  TLorentzVector momentum;
  Int_t number;

  for(number = 0; number < myhepevt.nhep; ++number)
  {
    element = static_cast<TRootC::GenParticle*>(branch->NewEntry());

    element->PID = myhepevt.idhep[number];
    element->Status = myhepevt.isthep[number];
    element->M1 = myhepevt.jmohep[number][0] - 1;
    element->M2 = myhepevt.jmohep[number][1] - 1;
    element->D1 = myhepevt.jdahep[number][0] - 1;
    element->D2 = myhepevt.jdahep[number][1] - 1;
    PdgParticle pdg_part = pdg_table[element->PID];
    element->Charge = pdg_part.charge();
    element->M = pdg_part.mass();

    element->E = myhepevt.phep[number][3];
    element->Px = myhepevt.phep[number][0];
    element->Py = myhepevt.phep[number][1];
    element->Pz = myhepevt.phep[number][2];
    //element->M = myhepevt.phep[number][4];

    momentum.SetPxPyPzE(element->Px, element->Py, element->Pz, element->E);
    element->PT = momentum.Perp();
    signPz = (element->Pz >= 0.0) ? 1.0 : -1.0;
    //element->Eta = element->PT == 0.0 ? signPz*999.9 : momentum.Eta(); to avoid a warning from ROOT, replace the "==0" by "< 1e-6"
    element->Eta = element->PT <1e-6 ? signPz*999.9 : momentum.Eta();
    element->Phi = momentum.Phi();

    element->T = myhepevt.vhep[number][3];
    element->X = myhepevt.vhep[number][0];
    element->Y = myhepevt.vhep[number][1];
    element->Z = myhepevt.vhep[number][2];
  }
}


//------------------------------------------------------------------------------

STDHEPConverter::~STDHEPConverter()
{
}

//------------------------------------------------------------------------------

STDHEPConverter::STDHEPConverter(const string& inputFileList, const string& outputFileName, const PdgTable& pdg, const int& Nevents) : DataConverter(pdg, Nevents)
{
  int ierr, entryType;
  int istr = 0;
  int nevt = 0, nevt_already_processed=0;

  ExRootTreeWriter *treeWriter = new ExRootTreeWriter(outputFileName, "GEN");

  // information about generated event
  ExRootTreeBranch *branchGenEvent = treeWriter->NewBranch("Event", TRootGenEvent::Class());
  // generated particles from HEPEVT
  ExRootTreeBranch *branchGenParticle = treeWriter->NewBranch("Particle", TRootC::GenParticle::Class());
  
  // Open a stream connected to an event file:
  ifstream infile(inputFileList.c_str());
  string filename;
  if(!infile.is_open())
    {
      cerr << left  << setw(30) <<"**   ERROR: Can't open  "<<""
                    << left  << setw(20) << inputFileList                    <<""
                    << right << setw(19) <<"for input **"<<""<<endl;
      exit(1);
    }

  // Parsing the input list
  while(1)
    {
      infile >> filename;
      if(!infile.good()) break; // end of the list
      if (Nevt>0 && nevt_already_processed >=Nevt) break; // enough events already processed

      ifstream checking_the_file(filename.c_str());
      if(!checking_the_file.good()) { 
             cerr << left  << setw(30) <<"**   ERROR: Can't find file  "<<""
                  << left  << setw(20) << filename                      <<""
                  << right << setw(19) <<"for input **"<<""<<endl;
             continue;
       } else checking_the_file.close();
      
      ierr = StdHepXdrReadInit(const_cast<char*>(filename.c_str()), &nevt, istr);

      if(ierr != 0) {
          cerr <<"**                ERROR: Can't open file for input                 **"<<endl;
          break;
      }


      Long64_t allEntries = (Nevt<0)?nevt:min(nevt,Nevt+1); // do not miss the "+1"
      if(allEntries > 0)
        {
          // Loop over all objects
          Long64_t entry = 0;
          Long64_t recordNumber = 1;
          for(entry = 0; entry < allEntries; ++entry)
            {
              ierr = StdHepXdrRead(&entryType, istr);
              if (Nevt>0 && nevt_already_processed >=Nevt) break; // enough events already processed
              
              if(ierr != 0) {
                  cerr << left  << setw(49) <<"**   ERROR: Unexpected end of file after"<<""
                       << left  << setw(10) << entry                                    <<""
                       << right << setw(10) <<"entries **"<<""<<endl;
                  break;
              } // if error
              
              // analyse only entries with standard HEPEVT common block
              if(entryType == 1) {
                  // add empty events for missing event numbers
                  while(recordNumber < myhepevt.nevhep) {
                      treeWriter->Clear();
                      AnalyseEvent(branchGenEvent, recordNumber);
                      treeWriter->Fill();
                      ++recordNumber;
                      ++nevt_already_processed;
                  } // while
                  
                  treeWriter->Clear();
                  AnalyseEvent(branchGenEvent, myhepevt.nevhep);
                  AnalyseParticles(branchGenParticle);
                  treeWriter->Fill();
                  ++recordNumber;
                  ++nevt_already_processed;
                } // if
            } // for
        } // if allEntries >0
      StdHepXdrEnd(istr);
    } // while inputlist not empty
  treeWriter->Write();
  
  delete treeWriter;

}