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

/***********************************************************************
**                                                                    **
**          /----------------------------------------------\          **
**         |  Delphes, a framework for the fast simulation  |         **
**         |       of a  generic collider experiment        |         **
**          \----------------------------------------------/          **
**                                                                    **
**                                                                    **
**   This package uses:                                               **
**   ------------------                                               **
**   FastJet algorithm: Phys. Lett. B641 (2006) [hep-ph/0512210]      ** 
**   Hector: JINST 2:P09005 (2007) [physics.acc-ph:0707.1198v2]       ** 
**   FROG: [hep-ex/0901.2718v1]                                       **
**                                                                    **
** ------------------------------------------------------------------ **
**                                                                    **
**   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)
{
  int Nmax = -1;  // useless for the moment -- should be changed in order
  if(Nmax>0) Nmax = min(Nmax,myhepevt.nhep); // to use Nevt in DataConverter
  else Nmax = myhepevt.nhep;

  TRootC::GenParticle *element;

  Double_t signPz;
  TLorentzVector momentum;
  Int_t number;

  for(number = 0; number < Nmax; ++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;

  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);
    }
  while(1)
    {
      infile >> filename;
      if(!infile.good()) break;
        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(ierr != 0)
                {
                  cerr << left  << setw(49) <<"**   ERROR: Unexpected end of file after"<<""
                       << left  << setw(10) << entry                                    <<""
                       << right << setw(10) <<"entries **"<<""<<endl;
                  break;
                }
              
              // 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;
                    }
                  
                  treeWriter->Clear();
                  
                  AnalyseEvent(branchGenEvent, myhepevt.nevhep);
                  AnalyseParticles(branchGenParticle);
                  
                  treeWriter->Fill();
                  
                  ++recordNumber;
                  
                }
            }
        }
      StdHepXdrEnd(istr);
    }
  treeWriter->Write();
  
  delete treeWriter;

}