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

/// \file Resolution_ATLAS.cpp
/// \brief Resolution for Atlas

#include "TChain.h"
#include "TApplication.h"
#include "TFile.h"
#include "TStopwatch.h"

#include "ExRootTreeReader.h"
#include "ExRootTreeWriter.h"
#include "ExRootTreeBranch.h"
#include "ExRootProgressBar.h"
#include "TreeClasses.h"

#include "SmearUtil.h"
#include "JetsUtil.h"
#include "BFieldProp.h"
#include "PdgParticle.h"

#include<vector>
#include<iostream>

using namespace std;

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

// //********************************** PYTHIA INFORMATION*********************************

TSimpleArray<TRootGenParticle> TauHadr(const TClonesArray *GEN)
 {
   TIter it((TCollection*)GEN);
   it.Reset();
   TRootGenParticle *gen1;
   TSimpleArray<TRootGenParticle> array,array2;

   while((gen1 = (TRootGenParticle*) it.Next()))
      {
        array.Add(gen1);
      }
   it.Reset();
   bool tauhad;
   while((gen1 = (TRootGenParticle*) it.Next()))
    {
       tauhad=false;
       if(abs(gen1->PID)==15)
         {
            int d1=gen1->D1;
            int d2=gen1->D2;
            if((d1 < array.GetEntries()) && (d1 > 0) && (d2 < array.GetEntries()) && (d2 > 0))
               {
                 //if(fabs(gen1->Eta)<2.5) tauhad=true;
                 tauhad=true;
                 for(int d=d1; d < d2+1; d++)
                    {
                      if(abs(array[d]->PID)== pE || abs(array[d]->PID)== pMU)tauhad=false;
                    }
               }
         }
        if(tauhad)array2.Add(gen1);
     }
   return array2;
 }

double EnergySmallCone(const vector<TLorentzVector> &towers, const float eta, const float phi,float energy_scone,float JET_seed) {
  double Energie=0;
  for(unsigned int i=0; i < towers.size(); i++) {
      if(towers[i].Pt() < JET_seed) continue;
      if((DeltaR(phi,eta,towers[i].Phi(),towers[i].Eta()) < energy_scone)) {
          Energie += towers[i].E();
      }
  }
  return Energie;
}


void PairingJet(TLorentzVector &JETSm, const TLorentzVector &JET, const TClonesArray *branchJet, const float dR=0.25)
{
  JETSm.SetPxPyPzE(0,0,0,0);
  float deltaRtest=5000;
  TIter itJet((TCollection*)branchJet);
  TRootJet *jet;
  itJet.Reset();
  while( (jet = (TRootJet*) itJet.Next()) )
   {
      TLorentzVector Att;
      Att.SetPtEtaPhiE(jet->PT,jet->Eta,jet->Phi,jet->E);
      if(DeltaR(JET.Phi(),JET.Eta(),Att.Phi(),Att.Eta()) < deltaRtest)
        {
          deltaRtest = DeltaR(JET.Phi(),JET.Eta(),Att.Phi(),Att.Eta());
          if(deltaRtest < dR)
             {
               JETSm = Att;
             }
        }
    }
}

void PairingElec(TLorentzVector &ELECSm, const TLorentzVector &ELEC, const TClonesArray *branchElec)
{
  ELECSm.SetPxPyPzE(0,0,0,0);
  float deltaRtest=5000;
  TIter itElec((TCollection*)branchElec);
  TRootElectron *elec;
  itElec.Reset();
  while( (elec = (TRootElectron*) itElec.Next()) )
   {
      TLorentzVector Att;
      Att.SetPtEtaPhiE(elec->PT,elec->Eta,elec->Phi,elec->E);
      if(DeltaR(ELEC.Phi(),ELEC.Eta(),Att.Phi(),Att.Eta()) < deltaRtest)
        {
          deltaRtest = DeltaR(ELEC.Phi(),ELEC.Eta(),Att.Phi(),Att.Eta());
          if(deltaRtest < 0.025)
             {
               ELECSm = Att;
             }
        }
    }
}

void PairingMuon(TLorentzVector &MUONSm, const TLorentzVector &MUON, const TClonesArray *branchMuon)
{
  MUONSm.SetPxPyPzE(0,0,0,0);
  float deltaRtest=5000;
  TIter itMuon((TCollection*)branchMuon);
  TRootMuon *muon;
  itMuon.Reset();
  while( (muon = (TRootMuon*) itMuon.Next()) )
   {
      TLorentzVector Att;
      Att.SetPxPyPzE(muon->Px,muon->Py,muon->Pz,muon->E);
      if(DeltaR(MUON.Phi(),MUON.Eta(),Att.Phi(),Att.Eta()) < deltaRtest)
        {
          deltaRtest = DeltaR(MUON.Phi(),MUON.Eta(),Att.Phi(),Att.Eta());
          if(deltaRtest < 0.025)
             {
               MUONSm = Att;
             }
        }
    }
}

unsigned int NumTracks(const TClonesArray *branchTracks, const float pt_track, const float eta, const float phi,float track_scone) {
  unsigned int numtrack=0;
  TIter itTrack((TCollection*)branchTracks);
  TRootTracks *track;
  itTrack.Reset();
  while( (track = (TRootTracks*) itTrack.Next()) )
   {
      if((track->PT < pt_track )||
         (DeltaR(phi,eta,track->Phi,track->Eta) > track_scone)
        )continue;
      numtrack++;
  }
  return numtrack;
}



int main(int argc, char *argv[])
{
  int appargc = 2;
  char *appName= new char[20];
  char *appOpt= new char[20];
  sprintf(appName,"Resolution_ATLAS");
  sprintf(appOpt,"-b");
  char *appargv[] = {appName,appOpt};
  TApplication app(appName, &appargc, appargv);
  delete [] appName;
  delete [] appOpt;

  if(argc != 3) {
      cout << " Usage: " << argv[0] << " input_file" << " output_file" << endl;
      cout << " input_file - input file in Delphes-root format," << endl;
      cout << " output_file - output file." << endl;
      exit(1);
  }


  // 1. ********** initialisation ***********

  srand (time (NULL));         
  TStopwatch globalwatch;
  globalwatch.Start();
  
  string inputfilename(argv[1]), outputfilename(argv[2]);

  // 1.1 checks the input file
  if(inputfilename.find(".root") > inputfilename.length()   ) {
    cout << "input_file should be a .root file from Delphes!\n";
    return -1;
  }
  TFile f(inputfilename.c_str());
  if (!f.FindKey("GEN") || !f.FindKey("Analysis")           ) {
    cout << "input_file should be a .root file from Delphes!\n";
    cout << "it should contain both \"GEN\" and \"Analysis\" trees at least.\n";
    return -1;
  }
  f.Close();

  // 1.2 checks the output file
  if(outputfilename.find(".root") > outputfilename.length() ) {
    cout << "output_file should be a .root file!\n";
    return -1;
  }
  TFile *outputFile = TFile::Open(outputfilename.c_str(), "RECREATE");// Creates the file, but should be closed just after
  outputFile->Close();

  // 1.3 Reads the trees in  input file 
  TChain chainGEN("GEN");
  chainGEN.Add(inputfilename.c_str());
  //chainGEN.Add("all_dijets_atlas_kt.root");
  //chainGEN.Add("all_dijets_atlas_kt_2.root");
  ExRootTreeReader *treeReaderGEN = new ExRootTreeReader(&chainGEN);
  TChain chain("Analysis");
  chain.Add(inputfilename.c_str());
  //chain.Add("all_dijets_atlas_kt.root");
  //chain.Add("all_dijets_atlas_kt_2.root");
  ExRootTreeReader *treeReader = new ExRootTreeReader(&chain);
  const TClonesArray *branchJet = treeReader->UseBranch("Jet");
  const TClonesArray *branchElec = treeReader->UseBranch("Electron");
  const TClonesArray *branchMuon = treeReader->UseBranch("Muon");
  const TClonesArray *branchTracks = treeReader->UseBranch("Tracks");
  const TClonesArray *branchTowers = treeReader->UseBranch("CaloTower");
  const TClonesArray *branchGen = treeReaderGEN->UseBranch("Particle");
  TIter itGen((TCollection*)branchGen);
 
 
  // 1.4 Prepares the output root file
  ExRootTreeWriter *treeWriter = new ExRootTreeWriter(outputfilename, "Analysis");
  ExRootTreeBranch *branchjet = treeWriter->NewBranch("JetPTResol", RESOLJET::Class());
  ExRootTreeBranch *branchelec = treeWriter->NewBranch("ElecEResol", RESOLELEC::Class());
  ExRootTreeBranch *branchmuon = treeWriter->NewBranch("MuonPTResol", RESOLMUON::Class());
  ExRootTreeBranch *branchtaujet = treeWriter->NewBranch("TauJetPTResol", TAUHAD::Class());
  ExRootTreeBranch *branchetmis = treeWriter->NewBranch("ETmisResol",ETMIS::Class());
 
  // 1.5 other initialisations 
  TRootGenParticle *particle;
  RESOLELEC * elementElec;
  RESOLMUON *elementMuon;
  RESOLJET *elementJet;
  TAUHAD *elementTaujet;
  ETMIS *elementEtmis;
  
  int numTau=0, numTauRec=0;
  
  RESOLution *DET = new RESOLution();
  /*
  string detectorcard = "data/DetectorCard_CMS.dat";
  const float dR_jetpairing = 0.25;
  const float jet_pt_cut = 1;
  */
  string detectorcard = "data/DetectorCard_ATLAS.dat";
  const float dR_jetpairing = 0.2;
  const float jet_pt_cut = 7;
  DET->ReadDataCard(detectorcard);

  JetsUtil *JETRUN = new JetsUtil(detectorcard);

  TLorentzVector genMomentum(0,0,0,0); // generator level information
  TLorentzVector recoMomentum(0,0,0,0);// reconstruction level information
  LorentzVector jetMomentum;
  
  vector<fastjet::PseudoJet> input_particlesGEN;//for FastJet algorithm
  vector<fastjet::PseudoJet> sorted_jetsGEN;
  vector<int> NTrackJet;
  vector<TLorentzVector> towers;
 
  // 2. Loop over all events
  Long64_t entry, allEntries = treeReader->GetEntries();
  cout << endl << endl;
  cout <<"*********************************************************************"<< endl;
  cout <<"**                                                                 **"<< endl;
  cout <<"**        ####### RESOLUTION STUDIES FOR ATLAS   ########          **"<< endl;
  cout <<"**                                                                 **"<< endl;
  cout <<"**        ####### Start resolution processing    ########          **"<< endl;
  cout << left  << setw(52) <<"**              Total number of events to run: "<<""
       << left  << setw(15) << allEntries <<""
       << right << setw(2) <<"**"<<endl;

  ExRootProgressBar *Progress = new ExRootProgressBar(allEntries);

  for(entry = 0; entry < allEntries; ++entry)
    {
      Progress->Update(entry);
      TLorentzVector PTmisReco(0,0,0,0);
      TLorentzVector PTmisGEN(0,0,0,0);
      treeReader->ReadEntry(entry);
      treeReaderGEN->ReadEntry(entry);
      treeWriter->Clear();
      
      TSimpleArray<TRootGenParticle> bGen;
      itGen.Reset();
      TSimpleArray<TRootGenParticle> NFCentralQ;
      
      input_particlesGEN.clear();
      towers.clear();
      
      // Loop over all particles in event
      while( (particle = (TRootGenParticle*) itGen.Next()) )
	{
          genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
	  PdgParticle pdg_part = DET->PDGtable[particle->PID];
	  
	  int pid = abs(particle->PID);
	  float eta = fabs(particle->Eta);
	  
          //input generator level particle for jet algorithm	  
          if(particle->Status == 1 && eta < DET->CEN_max_calo_fwd)
	    {
	      input_particlesGEN.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E()));
	    }
	  
	  //Calculate ETMIS from generated particles
          //if((pid == pNU1) || (pid == pNU2) || (pid == pNU3))PTmisGEN = PTmisGEN + genMomentum;
          if( (particle->Status == 1) && pdg_part.invisible() )PTmisGEN = PTmisGEN + genMomentum;

	  //Electrons and muons	 
          if( (particle->Status == 1)    &&
            //((pid != pNU1) && (pid != pNU2) && (pid != pNU3)) &&
             (! pdg_part.invisible() ) &&
             (fabs(particle->Eta) < DET->CEN_max_calo_fwd)
            )
            {
              eta=fabs(genMomentum.Eta());

              switch(pid) {

                 case pE: // all electrons with eta < DET->MAX_CALO_FWD
                    PairingElec(recoMomentum,genMomentum,branchElec);
                    if(recoMomentum.E()!=0){
                     elementElec=(RESOLELEC*) branchelec->NewEntry();
                     elementElec->E = genMomentum.E();
                     elementElec->SmearedE = recoMomentum.E();}
                   break; // case pE
                 case pMU: // all muons with eta < DET->MAX_MU
                    PairingMuon(recoMomentum,genMomentum,branchMuon);
                    if(recoMomentum.E()!=0){
                     elementMuon = (RESOLMUON*) branchmuon->NewEntry();
                     elementMuon->OverPT = 1./genMomentum.Pt();
                     elementMuon->OverSmearedPT = 1./recoMomentum.Pt();}
                   break; // case pMU
                 default:
                   break;
              } // switch (pid)
           }

	} // while
     
      //compute missing transverse energy from calo towers    
      TIter itCalo((TCollection*)branchTowers);
      TRootCalo *calo;
      itCalo.Reset();
      TLorentzVector Att(0.,0.,0.,0.);
      float ScalarEt=0;
      while( (calo = (TRootCalo*) itCalo.Next()) )
           {
             if(calo->E !=0){
             Att.SetPtEtaPhiE(calo->getET(),calo->Eta,calo->Phi,calo->E);
             towers.push_back(Att);
             if(fabs(Att.Eta()) < DET->CEN_max_calo_fwd)
               {
                  ScalarEt = ScalarEt + calo->getET();
                  PTmisReco = PTmisReco + Att;
               }
           }  
      }
      elementEtmis= (ETMIS*) branchetmis->NewEntry();
      elementEtmis->Et = (PTmisGEN).Pt();
      elementEtmis->Ex = (-PTmisGEN).Px();
      elementEtmis->SEt = ScalarEt;
      elementEtmis->EtSmeare = (PTmisReco).Pt()-(PTmisGEN).Pt();
      elementEtmis->ExSmeare = (-PTmisReco).Px()-(PTmisGEN).Px();

      //*****************************
      sorted_jetsGEN=JETRUN->RunJetsResol(input_particlesGEN);

      TSimpleArray<TRootGenParticle> TausHadr = TauHadr(branchGen);
      TLorentzVector JETreco(0,0,0,0);
      for (unsigned int i = 0; i < sorted_jetsGEN.size(); i++) {
        TLorentzVector JETgen(0,0,0,0);
        JETgen.SetPxPyPzE(sorted_jetsGEN[i].px(),sorted_jetsGEN[i].py(),sorted_jetsGEN[i].pz(),sorted_jetsGEN[i].E());
        PairingJet(JETreco,JETgen,branchJet,dR_jetpairing);
        if(JETreco.Pt()>jet_pt_cut)
          {
            elementJet= (RESOLJET*) branchjet->NewEntry();
            elementJet->PT = JETgen.Et();
            elementJet->SmearedPT = JETreco.Et()/JETgen.Et();
            elementJet->E = JETgen.E();
            elementJet->dE = (JETreco.E()-JETgen.E())/JETgen.E() ;
            elementJet->dE2 = pow( (JETreco.E()-JETgen.E())/JETgen.E() , 2.) ;
	    elementJet->Eta = JETgen.Eta();
            elementJet->dE_reco = (JETreco.E()-JETgen.E())/JETreco.E() ;
            elementJet->dE2_reco = pow( (JETreco.E()-JETgen.E())/JETreco.E() , 2.) ;
          }
      }
      numTau += TausHadr.GetEntries();

      TIter itJet((TCollection*)branchJet);
      TRootJet *jet;
      itJet.Reset();
      while( (jet = (TRootJet*) itJet.Next()) )
        {
          TLorentzVector JETT(0,0,0,0);
          JETT.SetPxPyPzE(jet->Px,jet->Py,jet->Pz,jet->E);
          if(fabs(JETT.Eta()) < (DET->CEN_max_tracker - DET->TAU_track_scone)) 
           {
            for(Int_t i=0; i<TausHadr.GetEntries();i++)
              {
                if(DeltaR(TausHadr[i]->Phi,TausHadr[i]->Eta,JETT.Phi(),JETT.Eta())<0.1)
                  {
                    elementTaujet= (TAUHAD*) branchtaujet->NewEntry();
                    elementTaujet->EnergieCen = EnergySmallCone(towers,JETT.Eta(),JETT.Phi(),DET->TAU_energy_scone,DET->JET_seed)/JETT.E();
                    elementTaujet->NumTrack = NumTracks(branchTracks,DET->TAU_track_pt,JETT.Eta(),JETT.Phi(),DET->TAU_track_scone);
                    if( (EnergySmallCone(towers,JETT.Eta(),JETT.Phi(),DET->TAU_energy_scone,DET->JET_seed)/JETT.E()) > 0.95 
                     &&  (NumTracks(branchTracks,DET->TAU_track_pt,JETT.Eta(),JETT.Phi(),DET->TAU_track_scone))==1) if(TausHadr[i]->PT>30 && TausHadr[i]->PT<60) numTauRec++;
                 }
              }
           }
         
	
      } // for itJet : loop on all jets
      
      treeWriter->Fill();
    } // Loop over all events
  treeWriter->Write();
  globalwatch.Stop();

  cout << "This sample contains " << numTau << " taus that decayed hadronically.\n";
  cout << numTauRec << " have been reconstructed.\n";

  // Screen output
  cout <<"**                                                                 **"<< endl;
  cout <<"**                                                                 **"<< endl;
  cout <<"**        ################## Time report #################         **"<< endl;
  cout << left  << setw(32) <<"**              Time report for "<<""
       << left  << setw(15) << allEntries <<""
       << right << setw(22) <<"events         **"<<endl;
  cout <<"**                                                                 **"<< endl;
  cout << left  << setw(10) <<"**"<<""
       << left  << setw(15) <<"Time (s):"<<""
       << right << setw(15) <<"CPU"<<""
       << right << setw(15) <<"Real"<<""
       << right << setw(14)  <<"**"<<endl;
  cout << left  << setw(10) <<"**"<<""
       << left  << setw(15) <<" +  Global:"<<""
       << right << setw(15) <<globalwatch.CpuTime()<<""
       << right << setw(15) <<globalwatch.RealTime()<<""
       << right << setw(14)  <<"**"<<endl;


  string tempstring = detectorcard + " has been used.";
  cout <<"**        ################## Configuration ###############         **"<< endl;
  cout << left << setw(16) << "**        " << ""
       << left << setw(51) <<  tempstring << "**" << endl;

  tempstring =  outputfilename + " has been created.";
  cout << left << setw(16) << "**        " << ""
       << left << setw(51) <<  tempstring << "**" << endl;

  cout << left << setw(16) << "**        " << ""
       << left << setw(51) <<  get_time_date() << "**" << endl;

  cout <<"**                                                                 **"<< endl;
  cout <<"**                        Exiting ...                              **"<< endl;
  cout <<"**                                                                 **"<< endl;
  cout <<"**                                                                 **"<< endl;
  cout <<"*********************************************************************"<< endl;

  
  delete treeWriter;
  delete treeReader;
  delete DET;
}