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

// 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->SumPt) << "." << 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();
}