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

#include <iostream>
#include <utility>
#include <deque>
#include <sstream>
#include <fstream>

#include "TROOT.h"
#include "TApplication.h"
#include "TFile.h"
#include "TTree.h"
#include "TBranch.h"
#include "TLeaf.h"
#include "TLorentzVector.h"
#include "TChain.h"

#include "Utilities/ExRootAnalysis/interface/BlockClasses.h"
#include "Utilities/ExRootAnalysis/interface/ExRootTreeReader.h"
#include "Utilities/ExRootAnalysis/interface/ExRootTreeWriter.h"
#include "Utilities/ExRootAnalysis/interface/ExRootTreeBranch.h"
#include "Utilities/ExRootAnalysis/interface/LHEF.h"
#include "Utilities/stdhep/stdhep_mcfio.h"
#include "Utilities/stdhep/stdhep_declarations.h"

#include "interface/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 << "** ERROR: Can't open '" << inputFileList << "' 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()) { cout << filename << ": file not found\n"; 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;
      cout << "** Input file contains " << allEntries << " entries" << endl;
      
      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 << "** ERROR: Unexpected end of file after " << entry << " 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;
                  
                }
            }
        }
    }
  treeWriter->Write();
  cout << "** Exiting..." << endl;
  
  delete treeWriter;
  StdHepXdrEnd(istr);

}