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

/***********************************************************************
**                                                                    **
**          /----------------------------------------------\          **
**         |  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"
#include "SmearUtil.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()) {
           cerr << "ERROR: the extension of the input file is not '.root'. Exiting...\n"; valid_=false; 
        }
        else {
                outputfilename_.replace(r_find,5,"_events.lhco");
                ofstream outfile( outputfilename_.c_str());
                outfile.close();
                cout << left  << setw(20) <<"**        INFO: "<<""
                              << left  << setw(29) << outputfilename_ <<""
                              << right << setw(20) <<"created          **"<<endl;

        }
    }
  }
}



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

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

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

    // 2 muon data
    BranchReaderMuon(branchMuon,branchJet,line);

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

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

    // 6 MET
    BranchReaderETmis(branchETmis,line);

  } // event loop

  outfile.close();
  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;

    // parsing the branch
    for (unsigned int i=0; i< branch_size; i++) { 
        particle = (TRootParticle*) branch->At(i);
        TLorentzVector pmu;
        pmu.SetPtEtaPhiE(particle->PT, particle->Eta, particle->Phi, particle->E);
        float jmass = pmu.M();
//if(jmass < 0) cout << "jmass negatif" << particle->PT << " " <<  particle->Eta << " " <<  particle->Phi << " " <<  particle->E << " " << jmass << endl;
        float ntrk = 0.0;
        float btag =0;
        double ratioE = 0;

        switch(lhcoID) {
          case lhcoPhotonID: {
                TRootPhoton gam(*((TRootPhoton*) branch->At(i))); 
                jmass=0; ratioE = gam.EHoverEE; 
          } break;
          case lhcoElectronID: {
                TRootElectron elec(*((TRootElectron*) branch->At(i))); 
                ntrk = elec.Charge; ratioE = elec.EHoverEE; 
          } break;
          case lhcoTauJetID: {
                TRootTauJet taujet(*((TRootTauJet*) branch->At(i))); 
                ntrk = taujet.Charge;  ratioE = taujet.EHoverEE; 
          } break;
          case lhcoJetID: {
                TRootJet jet(*((TRootJet*) branch->At(i))); 
                ntrk = jet.NTracks; btag = jet.Btag; ratioE = jet.EHoverEE; 
          } break;
          case lhcoMuonID: {
           cout << "ERROR: LHCOConverter: for muons use LHCOConverter::BranchReaderMuon!\n";
          } break;
          case lhcoETmisID: {
           cout << "ERROR: LHCOConverter: for ETmis use LHCOConverter::BranchReaderETmis!\n";
          } break;

          default: break;
        } // switch

        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;
        }
    } // branch parsing

    outfile.close();
    // return outfile;
}

//ostringstream LHCOConverter::BranchReader(const TClonesArray* branch, unsigned int& line, unsigned short int lhcoID) const {
void LHCOConverter::BranchReaderETmis(const TClonesArray* branch, unsigned int& line) const {
    //ostringstream outfile;
    ofstream outfile( outputfilename_.c_str(),ios_base::app);
    TRootETmis * particle = 0;
    particle = (TRootETmis*) branch->At(0);
    outfile << fixed << setprecision(3)
            << setw(3) << line++    // line counter
            << setw(4) << 6 << " "   // object ID in lhco format
            << setw(7) << 0.000 <<" "
            << setw(7) << particle->Phi <<" "
            << setw(7) << particle->ET  <<" "
            << 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;
}


void LHCOConverter::BranchReaderMuon(const TClonesArray* branchMuon, const TClonesArray* branchJet, unsigned int& line) const {

    ofstream outfile( outputfilename_.c_str(),ios_base::app);
    unsigned int branch_size = branchMuon->GetEntries();
    TRootMuon * muon = 0;
    for (unsigned int i=0; i< branch_size; i++) {
        muon = (TRootMuon*) branchMuon->At(i);

        TLorentzVector pmu; pmu.SetPtEtaPhiE(muon->PT, muon->Eta, muon->Phi, muon->E);
        float jmass = pmu.M(); if(jmass <0) jmass =0; // small negative values are possible due to smearing. < 0.1 GeV

        // loop on jets, for muon isolation
        unsigned int closest_jet_ID = 0;
        float closest_jet_dR = 1E99; // initialisation value, should be high. No further impact
        for (unsigned int j=0; j< (unsigned int) branchJet->GetEntries(); j++) {
                TRootJet jet(*((TRootJet*)branchJet->At(j)));
                float dr_ij = DeltaR(muon->Phi, muon->Eta, jet.Phi, jet.Eta);
                if( dr_ij < closest_jet_dR ) {
                        closest_jet_ID = j;
                        closest_jet_dR = dr_ij;
                }
        } // loop on jets

        int ratio = 0; char ratio_string[3];
        if(muon->EtRatio != UNDEFINED) ratio =  static_cast<int>(muon->EtRatio);
        if(ratio<0) sprintf(ratio_string,"00");
        else if(ratio<10) sprintf(ratio_string,"0%d",ratio);
        else if(ratio<100)sprintf(ratio_string,"%d",ratio);
        else { cout << "Error: muon EtRatio > 99" << endl; sprintf(ratio_string,"99"); }
        // output
        outfile << fixed << setprecision(3)
                << setw(3) << line++    // line counter
                << setw(4) << lhcoMuonID << " "   // object ID in lhco format
                << setw(7) << muon->Eta <<" "
                << setw(7) << muon->Phi <<" "
                << setw(7) << muon->PT  <<" "
                << setw(7) << jmass <<" "       // invariant mass
                << setw(7) << (float)muon->Charge  <<" "        // -1 for mu-, +1 for mu+
                << setw(7) << (float)closest_jet_ID << " ";
        outfile << fixed << setprecision(0)
                << setw(4) << floor(muon->IsolPt) << "." << ratio_string<<" ";  // E_had/E_em
        outfile << fixed << setprecision(3)
                << setw(7) << 0.0 << setw(7) << 0.0  // dummy/dummy
                << endl;
    } // muon branch loop

    outfile.close();
}