source: svn/trunk/src/LHCOConverter.cc@ 296

/***********************************************************************
**                                                                    **
**          /----------------------------------------------\          **
**         |  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.                        **  
**                                                                    **
***********************************************************************/

// local includes
#include "LHCOConverter.h"
#include "ExRootTreeReader.h"
#include "BlockClasses.h"

// root includes
#include "TFile.h"
#include "TChain.h"
#include "TClonesArray.h"

// c++ includes
#include <iostream>
#include <string>
#include <fstream>
#include <sstream>
#include <iomanip>
using namespace std;


//------------------------------------------------------------------------------
extern const float UNDEFINED;

LHCOConverter::LHCOConverter(const string& inputroot, const string& inputlog):
  inputfilename_(inputroot), inputlogname_(inputlog), outputfilename_(inputroot), valid_(true) {
  if (inputfilename_ == "") {valid_ = false;}
  else { 
    TFile ftest(inputfilename_.c_str(),"READ");
    if (!ftest.IsOpen()) { cout << "ERROR: file " << inputfilename_ << " could not be opened\n"; valid_=false; 
    }
    else {
        ftest.Close();
        const unsigned int r_find = outputfilename_.rfind(".root");
        if(r_find > outputfilename_.size()) {
           cout << "ERROR: the extension of the input file is not '.root'. Exiting...\n"; valid_=false; 
        }
        else {
                outputfilename_.replace(r_find,5,".lhco");
                ofstream outfile( outputfilename_.c_str());
                outfile.close();
                cout << "INFO: " << outputfilename_ << " has been created\n";
        }
    }
  }
}



void LHCOConverter::CopyRunLogFile() {
  if(!valid_) {/*cout << "Nothing done\n";*/ return; }
  if(inputlogname_ == "") {
        // re-creates the logrun file name from the rootfile name
        inputlogname_ = inputfilename_;
        const unsigned int r_find = inputlogname_.rfind(".root");
        inputlogname_.replace(r_find,5,"_run.log");
  }
  ifstream infile( inputlogname_.c_str() );
  if (!infile.good()) { cout << "Warning: no input logfile found\n"; return; } 

  ofstream outfile( outputfilename_.c_str(),ios_base::app);
  
  // else, if find is found
  string linereader;
  while ( getline(infile,linereader) ) {
    outfile << "  # " << linereader << endl;
  }
  outfile << "  #" << endl;
  infile.close();
  outfile.close();
}



void LHCOConverter::ConvertExRootAnalysisToLHCO() {
  if(!valid_) {/*cout << "Nothing done\n";*/ return; }

  //output file
  ofstream outfile( outputfilename_.c_str(),ios_base::app);
  //#  typ      eta    phi      pt    jmas  ntrk  btag   had/em  dum1  dum2
  outfile << " ## More info on LHCO files: http://cp3wks05.fynu.ucl.ac.be/Manual/lhco.html\n  #" << endl;
  outfile << "  #  typ   eta     phi     pt     jmas    ntrk    btag   had/em    dum1   dum2" << endl;

  // input files
  TChain analysisChain("Analysis");
  TChain triggerChain("Trigger");
  analysisChain.Add(inputfilename_.c_str());
  triggerChain.Add(inputfilename_.c_str());
  ExRootTreeReader *analysisTree = new ExRootTreeReader(&analysisChain);
  ExRootTreeReader *triggerTree  = new ExRootTreeReader(&triggerChain);

  TClonesArray *branchPhoton = analysisTree->UseBranch("Photon");
  TClonesArray *branchElectron = analysisTree->UseBranch("Electron");
  TClonesArray *branchMuon = analysisTree->UseBranch("Muon");
  TClonesArray *branchTauJet = analysisTree->UseBranch("TauJet");
  TClonesArray *branchJet = analysisTree->UseBranch("Jet");
  TClonesArray *branchETmis = analysisTree->UseBranch("ETmis");

  Long64_t Nevents = analysisTree->GetEntries();
  cout << "** TTree 'Analysis' contains " << Nevents << " events" << endl;
  Nevents = min(Nevents,triggerTree->GetEntries());
  if (Nevents != analysisTree->GetEntries()) 
     cout << "** WARNING: not the same number of entries in 'Analysis' and      **\n**          'Trigger' trees\n";
  for(Long64_t event = 0; event < Nevents; ++event) {
    analysisTree->ReadEntry(event);
    triggerTree->ReadEntry(event);
    unsigned int triginfo = 0;
    unsigned int line =0;

    outfile << setw(3) << line++  << setw(4) << " " << setw(7) << event << setw(7) << triginfo << endl;
    outfile.close();

    // 0 photon data
    //outfile << BranchReader(branchPhoton,line,lhcoPhotonID);
    BranchReader(branchPhoton,line,lhcoPhotonID);

    // 1 electron/positron data
    BranchReader(branchElectron,line,lhcoElectronID);

    // 2 muon data
    BranchReader(branchMuon,line,lhcoMuonID);

    // 3 tau-jets
    BranchReader(branchTauJet,line,lhcoTauJetID);

    // 4 jets
    BranchReader(branchJet,line,lhcoJetID);

    // 6 MET
    BranchReader(branchETmis,line,lhcoETmisID);

  } // event loop

  delete triggerTree;
  delete analysisTree;
}


//ostringstream LHCOConverter::BranchReader(const TClonesArray* branch, unsigned int& line, unsigned short int lhcoID) const {
void LHCOConverter::BranchReader(const TClonesArray* branch, unsigned int& line, unsigned short int lhcoID) const {
    //ostringstream outfile;
    ofstream outfile( outputfilename_.c_str(),ios_base::app);
    unsigned int branch_size = branch->GetEntries();
    TRootParticle * particle = 0;
    for (unsigned int i=0; i< branch_size; i++) {
        particle = (TRootParticle*) branch->At(i);
        double temp =  sqrt(particle->PT*particle->PT + particle->Pz * particle->Pz);
        double jmass = (particle->E - temp) * (particle->E + temp); // prevents some numerical sensitivity
        TLorentzVector pmu;
        pmu.SetPtEtaPhiE(particle->PT, particle->Eta, particle->Phi, particle->E);
        jmass = pmu.M();
        float ntrk = 0.0;
        if(lhcoID == lhcoElectronID) { TRootElectron elec(*((TRootElectron*) branch->At(i))); ntrk = elec.Charge; }
        else if (lhcoID == lhcoMuonID) { TRootMuon muon(*((TRootMuon*) branch->At(i))); ntrk = muon.Charge; }
        else if (lhcoID == lhcoTauJetID) { TRootTauJet taujet(*((TRootTauJet*) branch->At(i))); ntrk = taujet.Charge; }
        float btag =0;
        double ratioE = 0;
        if(lhcoID != lhcoETmisID) {
        outfile << fixed << setprecision(3)
                << setw(3) << line++    // line counter
                << setw(4) << lhcoID << " "   // object ID in lhco format
                << setw(7) << particle->Eta <<" "
                << setw(7) << particle->Phi <<" "
                << setw(7) << particle->PT  <<" "
                << setw(7) << jmass <<" "       // invariant mass
                << setw(7) << ntrk  <<" "       // number of tracks
                << setw(7) << btag  <<" "       
                << setw(7) << ratioE <<" "      // E_had/E_em
                << setw(7) << 0.0 << setw(7) << 0.0  // dummy/dummy
                << endl;
        } else {
        outfile << fixed << setprecision(3)
                << setw(3) << line++    // line counter
                << setw(4) << lhcoID << " "   // object ID in lhco format
                << setw(7) << 0.000 <<" "
                << setw(7) << particle->Phi <<" "
                << setw(7) << particle->PT  <<" "
                << setw(7) << 0. <<" "       // invariant mass
                << setw(7) << 0  <<" "       // number of tracks
                << setw(7) << 0.  <<" "
                << setw(7) << 0. <<" "      // E_had/E_em
                << setw(7) << 0.0 << setw(7) << 0.0  // dummy/dummy
                << endl;
        }
    }
    outfile.close();
    // return outfile;
}