source: svn/trunk/src/LHEFConverter.cc@ 889

/***********************************************************************
**                                                                    **
**          /----------------------------------------------\          **
**         |  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 <cstdlib> // for exit()

#include "TLorentzVector.h"

#include "ExRootTreeWriter.h"
#include "BlockClasses.h"
#include "LHEFConverter.h"
#include "SmearUtil.h"
#include "PdgParticle.h"

using namespace std;


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

void LHEFConverter::AnalyseEvent(LHEF::Reader *reader, ExRootTreeBranch *branch, const Long64_t eventNumber)
{
  const LHEF::HEPEUP &hepeup = reader->hepeup;

  TRootLHEFEvent *element;

  element = (TRootLHEFEvent*) branch->NewEntry();

  element->Number = eventNumber;
  element->Nparticles = hepeup.NUP;
  element->ProcessID = hepeup.IDPRUP;
  element->Weight = hepeup.XWGTUP;
  element->ScalePDF = hepeup.SCALUP;
  element->CouplingQED = hepeup.AQEDUP;
  element->CouplingQCD = hepeup.AQCDUP;
}

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

void LHEFConverter::AnalyseParticles(LHEF::Reader *reader, ExRootTreeBranch *branch)
{
  const LHEF::HEPEUP &hepeup = reader->hepeup;
  Double_t signPz;
  TLorentzVector momentum;
  TRootLHEFParticle *element;

  for(Int_t particle = 0; particle < hepeup.NUP; ++particle)
  {
    element = (TRootLHEFParticle*) branch->NewEntry();

    element->PID = hepeup.IDUP[particle];
    element->Status = hepeup.ISTUP[particle];
    element->M1 = hepeup.MOTHUP[particle].first;
    element->M2 = hepeup.MOTHUP[particle].second;
    element->D1 = 0;
    element->D2 = 0;
    element->ColorLine1 = hepeup.ICOLUP[particle].first;
    element->ColorLine2 = hepeup.ICOLUP[particle].second;
    element->Px = hepeup.PUP[particle][0];
    element->Py = hepeup.PUP[particle][1];
    element->Pz = hepeup.PUP[particle][2];
    element->E = hepeup.PUP[particle][3];
    //element->M = hepeup.PUP[particle][4];
    //element->Charge = ChargeVal(element->PID);
    
    PdgParticle pdg_part = pdg_table[element->PID];
    element->Charge = pdg_part.charge();
    element->M = pdg_part.mass();
    
    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->Rapidity = element->PT < 1E-6 ? signPz*999.9 : momentum.Rapidity();
    element->LifeTime = hepeup.VTIMUP[particle];
    element->Spin = hepeup.SPINUP[particle];
  }
}


//---------------------------------------------------------------------------
/*
void LHEFConverter::AnalyseParticles(LHEF::Reader *reader, ExRootTreeBranch *branch)
{
  const LHEF::HEPEUP &hepeup = reader->hepeup;
  Double_t signPz;
  TLorentzVector momentum;
  TRootC::GenParticle *element;

  for(Int_t particle = 0; particle < hepeup.NUP; ++particle)
  {
    element = (TRootC::GenParticle*) branch->NewEntry();

    element->PID = hepeup.IDUP[particle];
    element->Status = hepeup.ISTUP[particle];
    element->M1 = hepeup.MOTHUP[particle].first;
    element->M2 = hepeup.MOTHUP[particle].second;
    //element->ColorLine1 = hepeup.ICOLUP[particle].first;
    //element->ColorLine2 = hepeup.ICOLUP[particle].second;
    element->D1 = 0;
    element->D2 = 0;

    PdgParticle pdg_part = pdg_table[element->PID];
    element->Charge = pdg_part.charge();
    element->M = pdg_part.mass();
    //element->Charge = ChargeVal(element->PID);
    //element->M = hepeup.PUP[particle][4];

    element->E = hepeup.PUP[particle][3];
    element->Px = hepeup.PUP[particle][0];
    element->Py = hepeup.PUP[particle][1];
    element->Pz = hepeup.PUP[particle][2];

    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->Rapidity = element->PT < 1E-6 ? signPz*999.9 : momentum.Rapidity();
//    element->LifeTime = hepeup.VTIMUP[particle];
//    element->Spin = hepeup.SPINUP[particle];
  }
}
*/

LHEFConverter::~LHEFConverter()
{
}

//------------------------------------------------------------------------------
LHEFConverter::LHEFConverter(const string& inputFileList, const string& outputFileName, const PdgTable& pdg, const int& Nevents) : DataConverter(pdg,Nevents) {

  ExRootTreeWriter *treeWriter = new ExRootTreeWriter(outputFileName, "GEN");
  
  // generated event from LHEF
  ExRootTreeBranch *branchEvent = treeWriter->NewBranch("Event", TRootLHEFEvent::Class());
  
  // generated partons from LHEF
  ExRootTreeBranch *branchParticle = treeWriter->NewBranch("Particle", TRootLHEFParticle::Class());
  //ExRootTreeBranch *branchParticle = treeWriter->NewBranch("Particle", TRootC::GenParticle::Class());

  int nevt_already_processed=0;
  
  // 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) { // parsing the list of files
    
      infile >> filename;
      if(!infile.good()) break;
      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(36) << filename                       <<""
                  << right << setw(3) <<" **"<<""<<endl;
              continue;
           }
        else checking_the_file.close();

      // Create the Reader object:
      LHEF::Reader *inputReader = new LHEF::Reader(filename);
      
      Long64_t allEntries = inputReader->getNumberOfEvents();
      allEntries = (Nevt<0)?allEntries:min((int)allEntries,Nevt); // do not miss the "+1"      

      if(allEntries > 0) {
          // Loop over all events
          Long64_t entry = 0;
          while(inputReader->readEvent()) {
              if (Nevt>0 && nevt_already_processed >=Nevt) break; // enough events already processed
              treeWriter->Clear();
              AnalyseEvent(inputReader, branchEvent, entry + 1);
              AnalyseParticles(inputReader, branchParticle);
              treeWriter->Fill();
              ++entry;
              ++nevt_already_processed;
              if(allEntries<entry+1) break;
          }
      }
  }
  treeWriter->Write();
  delete treeWriter;
  //delete inputReader;
}