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

/*
  ---- Delphes ----
  A Fast Simulator for general purpose LHC detector
  S. Ovyn ~~~~ severine.ovyn@uclouvain.be

  Center for Particle Physics and Phenomenology (CP3)
  Universite Catholique de Louvain (UCL)
  Louvain-la-Neuve, Belgium
*/

#include <iostream>
#include <fstream>
#include "TLorentzVector.h"
#include "BlockClasses.h"
#include "ExRootTreeWriter.h"
#include "ExRootTreeBranch.h"
#include "stdhep_mcfio.h"
#include "stdhep_declarations.h"
#include "STDHEPConverter.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)
{
  TRootGenParticle *element;

  Double_t signPz;
  TLorentzVector momentum;
  Int_t number;

  for(number = 0; number < myhepevt.nhep; ++number)
  {
    element = static_cast<TRootGenParticle*>(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;
//    element->Charge = myhepevt.hepchg(element->PID);

    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();
    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)
{

  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", TRootGenParticle::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;
      
      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;

}