/*
  ---- Delphes ----
  A Fast Simulator for general purpose LHC detector
  S. Ovyn ~~~~ severine.ovyn@uclouvain.be

  Center for Particle Physics and Phenomenology (CP3)
  Universite Catholique de Louvain (UCL)
  Louvain-la-Neuve, Belgium
*/

/// \file Delphes.cpp
/// \brief executable for the Delphes

#include "TChain.h"
#include "TApplication.h"

#include "Utilities/ExRootAnalysis/interface/ExRootTreeReader.h"
#include "Utilities/ExRootAnalysis/interface/ExRootTreeWriter.h"
#include "Utilities/ExRootAnalysis/interface/ExRootTreeBranch.h"

#include "H_BeamParticle.h"
#include "H_BeamLine.h"
#include "H_RomanPot.h"

#include "interface/DataConverter.h"
#include "interface/HEPEVTConverter.h"
#include "interface/LHEFConverter.h"
#include "interface/STDHEPConverter.h"

#include "interface/SmearUtil.h"
#include "Utilities/Fastjet/include/fastjet/PseudoJet.hh"
#include "Utilities/Fastjet/include/fastjet/ClusterSequence.hh"

// get info on how fastjet was configured
#include "Utilities/Fastjet/include/fastjet/config.h"

// make sure we have what is needed
#ifdef ENABLE_PLUGIN_SISCONE
#  include "Utilities/Fastjet/plugins/SISCone/SISConePlugin.hh"
#endif
#ifdef ENABLE_PLUGIN_CDFCONES
#  include "Utilities/Fastjet/plugins/CDFCones/CDFMidPointPlugin.hh"
#  include "Utilities/Fastjet/plugins/CDFCones/CDFJetCluPlugin.hh"
#endif

#include<vector>
#include<iostream>



using namespace std;

//------------------------------------------------------------------------------
void todo(string filename) {
  ifstream infile(filename.c_str());
  cout << "** TODO list ..." << endl;
  while(infile.good()) {
	string temp;
	getline(infile,temp);
	cout << "*" << temp << endl;
  }
  cout << "** done...\n"; 
}

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

int main(int argc, char *argv[])
{
  int appargc = 2;
  char *appName = "JetsSim";
  char *appargv[] = {appName, "-b"};
  TApplication app(appName, &appargc, appargv);
  
  if(argc != 4 && argc != 3) {
      cout << " Usage: " << argv[0] << " input_file" << " output_file" << " data_card " << endl;
      cout << " input_list - list of files in Ntpl, StdHep of LHEF format," << endl;
      cout << " output_file - output file." << endl;
      cout << " data_card - Datacard containing resolution variables for the detector simulation (optional) "<<endl;
      exit(1);
  }

  srand (time (NULL));         /* Initialisation du générateur */
  
  //read the input TROOT file
  string inputFileList(argv[1]), outputfilename(argv[2]);
  if(outputfilename.find(".root") > outputfilename.length() ) {
	cout << "output_file should be a .root file!\n"; 
  	exit(1);
  }

  TFile *outputFile = TFile::Open(outputfilename.c_str(), "RECREATE"); // Creates the file, but should be closed just after
  outputFile->Close();

  string line;
  ifstream infile(inputFileList.c_str());
  infile >> line; // the first line determines the type of input files
 
  DataConverter *converter=0;
   
  if(strstr(line.c_str(),".hep"))
    {                           
      cout<<"*************************************************************************"<<endl;
      cout<<"************         StdHEP file format detected           **************"<<endl;
      cout<<"************     Starting convertion to TRoot format       **************"<<endl;
      cout<<"*************************************************************************"<<endl;
      converter = new STDHEPConverter(inputFileList,outputfilename);//case ntpl file in input list
    }
  else if(strstr(line.c_str(),".lhe"))
    {
      cout<<"*************************************************************************"<<endl;
      cout<<"************          LHEF file format detected            **************"<<endl;
      cout<<"************     Starting convertion to TRoot format       **************"<<endl;
      cout<<"*************************************************************************"<<endl;
      converter = new LHEFConverter(inputFileList,outputfilename);//case ntpl file in input list
    }
  else if(strstr(line.c_str(),".root"))
    {
      cout<<"*************************************************************************"<<endl;
      cout<<"************         h2root file format detected           **************"<<endl;
      cout<<"************     Starting convertion to TRoot format       **************"<<endl;
      cout<<"*************************************************************************"<<endl;
      converter = new HEPEVTConverter(inputFileList,outputfilename);//case ntpl file in input list
    }
  else { cout << "***  " << line.c_str() << "\n***  file format not identified\n***  Exiting\n"; return -1;};

  TChain chain("GEN");
  chain.Add(outputfilename.c_str());
  ExRootTreeReader *treeReader = new ExRootTreeReader(&chain);
  const TClonesArray *branchGen = treeReader->UseBranch("Particle");
  TIter itGen((TCollection*)branchGen);
  
  //write the output root file
  ExRootTreeWriter *treeWriter = new ExRootTreeWriter(outputfilename, "Analysis");
  ExRootTreeBranch *branchJet = treeWriter->NewBranch("Jet", TRootJet::Class());
  ExRootTreeBranch *branchTauJet = treeWriter->NewBranch("TauJet", TRootTauJet::Class());
  ExRootTreeBranch *branchElectron = treeWriter->NewBranch("Electron", TRootElectron::Class());
  ExRootTreeBranch *branchMuon = treeWriter->NewBranch("Muon", TRootMuon::Class());
  ExRootTreeBranch *branchPhoton = treeWriter->NewBranch("Photon", TRootPhoton::Class());
  ExRootTreeBranch *branchTracks = treeWriter->NewBranch("Tracks", TRootTracks::Class());
  ExRootTreeBranch *branchETmis = treeWriter->NewBranch("ETmis", TRootETmis::Class());
  ExRootTreeBranch *branchCalo = treeWriter->NewBranch("CaloTower", TRootCalo::Class());
  ExRootTreeBranch *branchZDC = treeWriter->NewBranch("ZDChits", TRootZdcHits::Class());
  ExRootTreeBranch *branchRP220 = treeWriter->NewBranch("RP220hits", TRootRomanPotHits::Class());
  ExRootTreeBranch *branchFP420 = treeWriter->NewBranch("FP420hits", TRootRomanPotHits::Class());

 
  TRootGenParticle *particle; 
  TRootETmis *elementEtmis;
  TRootElectron *elementElec;
  TRootMuon *elementMu;
  TRootPhoton *elementPhoton;
  TRootJet *elementJet;
  TRootTauJet *elementTauJet;
  TRootTracks *elementTracks;
  TRootCalo *elementCalo;
  TRootZdcHits *elementZdc;
  TRootRomanPotHits *elementRP220, *elementFP420;
  
  //read the datacard input file
  string DetDatacard("");
  if(argc==4)  DetDatacard =argv[3];
  RESOLution *DET = new RESOLution();   
  DET->ReadDataCard(DetDatacard);
  
  
  TLorentzVector genMomentum(0,0,0,0);
  LorentzVector jetMomentum;
  vector<TLorentzVector> TrackCentral;  
  vector<PhysicsTower> towers;

  vector<fastjet::PseudoJet> input_particles;//for FastJet algorithm
  vector<fastjet::PseudoJet> inclusive_jets;
  vector<fastjet::PseudoJet> sorted_jets;
  
  //Initialisation of Hector
  extern bool relative_energy;
  relative_energy = true; // should always be true
  extern int kickers_on;
  kickers_on = 1;	  // should always be 1
  
  // user should provide : (1) optics file for each beamline, and IPname, 
  // and offset data (s,x) for optical elements
  H_BeamLine* beamline1 = new H_BeamLine(1,500.);
  beamline1->fill("data/LHCB1IR5_v6.500.tfs",1,"IP5");
  beamline1->offsetElements(120,-0.097);
  H_RomanPot * rp220_1 = new H_RomanPot("rp220_1",220,2000); // RP 220m, 2mm, beam 1
  H_RomanPot * rp420_1 = new H_RomanPot("rp420_1",420,4000); // RP 420m, 4mm, beam 1
  beamline1->add(rp220_1);
  beamline1->add(rp420_1);
  
  H_BeamLine* beamline2 = new H_BeamLine(1,500.);
  beamline2->fill("data/LHCB1IR5_v6.500.tfs",-1,"IP5");
  beamline2->offsetElements(120,+0.097);
  H_RomanPot * rp220_2 = new H_RomanPot("rp220_2",220,2000);// RP 220m, 2mm, beam 2
  H_RomanPot * rp420_2 = new H_RomanPot("rp420_2",420,4000);// RP 420m, 4mm, beam 2
  beamline2->add(rp220_2);
  beamline2->add(rp420_2);

  // we will have four jet definitions, and the first three will be
  // plugins
  fastjet::JetDefinition jet_def;
  fastjet::JetDefinition::Plugin * plugins;

  switch(DET->JETALGO) {
  default:
  case 1: {

   // set up a CDF midpoint jet definition
   #ifdef ENABLE_PLUGIN_CDFCONES
   plugins = new fastjet::CDFJetCluPlugin(DET->C_SEEDTHRESHOLD,DET->CONERADIUS,DET->C_ADJACENCYCUT,DET->C_MAXITERATIONS,DET->C_IRATCH,DET->C_OVERLAPTHRESHOLD);
   jet_def = fastjet::JetDefinition(plugins);
   #else
   plugins = NULL;
   #endif
   }
   break;

   case 2: {

   // set up a CDF midpoint jet definition
   #ifdef ENABLE_PLUGIN_CDFCONES
   plugins = new fastjet::CDFMidPointPlugin (DET->M_SEEDTHRESHOLD,DET->CONERADIUS,DET->M_CONEAREAFRACTION,DET->M_MAXPAIRSIZE,DET->M_MAXPAIRSIZE,DET->C_OVERLAPTHRESHOLD);
   jet_def = fastjet::JetDefinition(plugins);
   #else
   plugins = NULL;
   #endif
   }
   break;
   case 3: {
   // set up a siscone jet definition
   #ifdef ENABLE_PLUGIN_SISCONE
   int npass = 0;               // do infinite number of passes
   double protojet_ptmin = 0.0; // use all protojets
   plugins = new fastjet::SISConePlugin (DET->CONERADIUS,DET->C_OVERLAPTHRESHOLD,npass, protojet_ptmin);
   jet_def = fastjet::JetDefinition(plugins);
   #else
   plugins = NULL;
   #endif
   }
   break;

  case 4: {
  jet_def = fastjet::JetDefinition(fastjet::kt_algorithm, DET->CONERADIUS);
  //jet_defs[4] = fastjet::JetDefinition(fastjet::cambridge_algorithm,jet_radius);
  //jet_defs[5] = fastjet::JetDefinition(fastjet::antikt_algorithm,jet_radius);
  }
  break;
  }
 
   // Loop over all events
  Long64_t entry, allEntries = treeReader->GetEntries();
  cout << "** Chain contains " << allEntries << " events" << endl;
  for(entry = 0; entry < allEntries; ++entry)
    {
      TLorentzVector PTmis(0,0,0,0);
      treeReader->ReadEntry(entry);
      treeWriter->Clear();
      
      if((entry % 100) == 0 && entry > 0 )  cout << "** Processing element # " << entry << endl;
      
      TSimpleArray<TRootGenParticle> bGen;
      itGen.Reset();
      TrackCentral.clear();
      towers.clear();
      input_particles.clear();
      TSimpleArray<TRootGenParticle> NFCentralQ;
      inclusive_jets.clear();
      sorted_jets.clear();

      // Loop over all particles in event
      while( (particle = (TRootGenParticle*) itGen.Next()) )
	{
          genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);

	  int pid = abs(particle->PID);
	  float eta = fabs(particle->Eta);

          //subarray of the quarks (i.e. not final particles, i.e status not equal to 1) 
          // in the tracker (i.e. for those we know the tracks)
          // with enough p_T
          //// This subarray is needed for the B-jet algorithm
	  // optimization for speed : put first PID condition, then ETA condition, then either pt or status
	  if( (pid <= pB || pid == pGLUON) &&// is it a light quark or a gluon, i.e. is it one of these : u,d,c,s,b,g ?
              eta < DET->MAX_TRACKER &&
              particle->Status != 1 &&
              particle->PT > DET->PT_QUARKS_MIN ) {
            NFCentralQ.Add(particle);
          }


	  // keeps only final particles, visible by the central detector, including the fiducial volume
	  // the ordering of conditions have been optimised for speed : put first the STATUS condition
          if( (particle->Status == 1)    &&
	      ( 
		(pid == pMU  && eta < DET->MAX_MU) ||
                (pid != pMU  && (pid != pNU1) && (pid != pNU2) && (pid != pNU3) && eta < DET->MAX_CALO_FWD) 
              )
            ) { 
		switch(pid) {
		    
		    case pE: // all electrons with eta < DET->MAX_CALO_FWD
		      DET->SmearElectron(genMomentum);
		      if(genMomentum.E()!=0 && eta < DET->MAX_TRACKER) {
			  elementElec = (TRootElectron*) branchElectron->NewEntry();
			  elementElec->Set(genMomentum);
			  elementElec->Charge = sign(particle->PID);
		      }
		      break; // case pE

		    case pGAMMA: // all photons with eta < DET->MAX_CALO_FWD
		      DET->SmearElectron(genMomentum);
		      if(genMomentum.E()!=0 && eta < DET->MAX_TRACKER) {
			  elementPhoton = (TRootPhoton*) branchPhoton->NewEntry();
			  elementPhoton->Set(genMomentum);
		      }
		      break; // case pGAMMA

		    case pMU: // all muons with eta < DET->MAX_MU
		      DET->SmearMu(genMomentum);
		      if(genMomentum.E() !=0 ) {
			  elementMu = (TRootMuon*) branchMuon->NewEntry();
			  elementMu->Charge = sign(particle->PID);
			  elementMu->Set(genMomentum);
		      }	
		      break; // case pMU

		    case pLAMBDA: // all lambdas with eta < DET->MAX_CALO_FWD
		    case pK0S:    // all K0s with eta < DET->MAX_CALO_FWD
		      DET->SmearHadron(genMomentum, 0.7);
		      break; // case hadron
		      
		    default:   // all other final particles with eta < DET->MAX_CALO_FWD
		      DET->SmearHadron(genMomentum, 1.0);
		      break;
		} // switch (pid)
		
		// all final particles but muons and neutrinos	
		// for calorimetric towers and mission PT
		if(genMomentum.E()!=0) {
		  PTmis = PTmis + genMomentum;//ptmis
		  if(pid !=pMU) {
		    towers.push_back(PhysicsTower(LorentzVector(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E())));
		    // create a fastjet::PseudoJet with these components and put it onto
		    // back of the input_particles vector
		    input_particles.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E())); 
		    elementCalo = (TRootCalo*) branchCalo->NewEntry();
		    elementCalo->Set(genMomentum);
		  }
		}
		
		
		// all final charged particles
		if(
		   ((rand()%100) < DET->TRACKING_EFF)  &&
		   (genMomentum.E()!=0) &&
		   (fabs(particle->Eta) < DET->MAX_TRACKER) && 
		   (genMomentum.Pt() > DET->PT_TRACKS_MIN ) &&     // pt too small to be taken into account
		   (pid != pGAMMA)    && 
		   (pid != pPI0)      &&
		   (pid != pK0L)      &&
		   (pid != pN)        &&
		   (pid != pSIGMA0)   &&
		   (pid != pDELTA0)   &&
		   (pid != pK0S)   // not charged particles : invisible by tracker
		   )
		  {
		    elementTracks = (TRootTracks*) branchTracks->NewEntry();
		    elementTracks->Set(genMomentum);
		    TrackCentral.push_back(genMomentum);
		  }
	  } // switch
	  
	  // Forward particles in CASTOR ?
	  /*	  if (particle->Status == 1 && (fabs(particle->Eta) > DET->MIN_CALO_VFWD) 
		  && (fabs(particle->Eta) < DET->MAX_CALO_VFWD)) {
		  
		  
		  } // CASTOR
	  */
	  // Zero degree calorimeter, for forward neutrons and photons
	  if (particle->Status ==1 && (pid == pN || pid == pGAMMA ) && eta > DET->MIN_ZDC ) {
	    // !!!!!!!!! vérifier que particle->Z est bien en micromètres!!!
	    // !!!!!!!!! vérifier que particle->T est bien en secondes!!!
	    // !!!!!!!!! pas de smearing ! on garde trop d'info !
	    elementZdc = (TRootZdcHits*) branchZDC->NewEntry();
	    elementZdc->Set(genMomentum);
	    
	    // time of flight t is t = T + d/[ cos(theta) v ]
	    //double tx = acos(particle->Px/particle->Pz);
	    //double ty = acos(particle->Py/particle->Pz);
	    //double theta = (1E-6)*sqrt( pow(tx,2) + pow(ty,2) );
	    //double flight_distance = (DET->ZDC_S - particle->Z*(1E-6))/cos(theta) ; // assumes that Z is in micrometers
	    double flight_distance = (DET->ZDC_S - particle->Z*(1E-6)); 
	    // assumes also that the emission angle is so small that 1/(cos theta) = 1
	    elementZdc->T = 0*particle->T + flight_distance/speed_of_light; // assumes highly relativistic particles 
	    elementZdc->side = sign(particle->Eta);
	    
	  }
	  
	  // if forward proton
          if( (pid == pP)  && (particle->Status == 1) &&  (fabs(particle->Eta) > DET->MAX_CALO_FWD) ) 
            {
	      // !!!!!!!! put here particle->CHARGE and particle->MASS
	      H_BeamParticle p1; /// put here particle->CHARGE and particle->MASS
	      p1.smearAng();
	      p1.smearPos();
	      p1.setPosition(p1.getX()-500.,p1.getY(),p1.getTX()-1*kickers_on*CRANG,p1.getTY(),0);
	      p1.set4Momentum(particle->Px,particle->Py,particle->Pz,particle->E);
	      
	      H_BeamLine *beamline;
	      if(particle->Eta >0) beamline = beamline1;
	      else beamline = beamline2;
	      
	      p1.computePath(beamline,1);
	      
	      if(p1.stopped(beamline)) {
		if (p1.getStoppingElement()->getName()=="rp220_1" || p1.getStoppingElement()->getName()=="rp220_2") {
		  p1.propagate(DET->RP220_S);
		  elementRP220 = (TRootRomanPotHits*) branchRP220->NewEntry();
		  elementRP220->X  = (1E-6)*p1.getX();  // [m]
		  elementRP220->Y  = (1E-6)*p1.getY();  // [m]
		  elementRP220->Tx = (1E-6)*p1.getTX(); // [rad]
		  elementRP220->Ty = (1E-6)*p1.getTY(); // [rad]
		  elementRP220->S = p1.getS();          // [m]
		  elementRP220->T = -1; 		// not yet implemented
		  elementRP220->E = p1.getE(); 	// not yet implemented
		  elementRP220->q2 = -1; 		// not yet implemented
		  elementRP220->side = sign(particle->Eta); 
		  
		} else if (p1.getStoppingElement()->getName()=="rp420_1" || p1.getStoppingElement()->getName()=="rp420_2") {
		  p1.propagate(DET->FP420_S);
		  elementFP420 = (TRootRomanPotHits*) branchFP420->NewEntry();
		  elementFP420->X  = (1E-6)*p1.getX();  // [m]
		  elementFP420->Y  = (1E-6)*p1.getY();  // [m]
		  elementFP420->Tx = (1E-6)*p1.getTX(); // [rad]
		  elementFP420->Ty = (1E-6)*p1.getTY(); // [rad]
		  elementFP420->S = p1.getS();	      // [m]
		  elementFP420->T = -1; 		      // not yet implemented
		  elementFP420->E = p1.getE();          // not yet implemented
		  elementFP420->q2 = -1;                // not yet implemented
		  elementFP420->side = sign(particle->Eta);
		}
	      }
	      
	      //		if(p1.stopped(beamline) && (p1.getStoppingElement()->getS() > 100))
	      //			cout << "Eloss ="  << 7000.-p1.getE() << " ; " << p1.getStoppingElement()->getName() << endl;
	    } // if forward proton
	  
	} // while
      
      // computes the Missing Transverse Momentum
      elementEtmis = (TRootETmis*) branchETmis->NewEntry();
      elementEtmis->ET = (-PTmis).Pt();
      elementEtmis->Phi = (-PTmis).Phi();
      elementEtmis->Px = (-PTmis).Px();
      elementEtmis->Py = (-PTmis).Py();
      
      //*****************************

      // run the jet clustering with the above jet definition
      if(input_particles.size()!=0)
	{
	  fastjet::ClusterSequence clust_seq(input_particles, jet_def);
	  
	  
	  // extract the inclusive jets with pt > 5 GeV
	  double ptmin = 5.0;
	  inclusive_jets = clust_seq.inclusive_jets(ptmin);
	  
	  // sort jets into increasing pt
	  sorted_jets = sorted_by_pt(inclusive_jets);
	}
      for (unsigned int i = 0; i < sorted_jets.size(); i++) {
	elementJet = (TRootJet*) branchJet->NewEntry();
	TLorentzVector JET;
	JET.SetPxPyPzE(sorted_jets[i].px(),sorted_jets[i].py(),sorted_jets[i].pz(),sorted_jets[i].E());
        //cout<<"Jet.Pt() "<<JET.Pt()<<endl;
	elementJet->Set(JET);
	// b-jets
	bool btag=false;
	if((fabs(JET.Eta()) < DET->MAX_TRACKER && DET->Btaggedjet(JET, NFCentralQ)))btag=true;
	elementJet->Btag = btag;
	
	// Tau jet identification : 1! track and electromagnetic collimation
	if(fabs(JET.Eta()) < (DET->MAX_TRACKER - DET->TAU_CONE_TRACKS)) {
	  double Energie_tau_central = DET->EnergySmallCone(towers,JET.Eta(),JET.Phi());
	  if(
	     ( Energie_tau_central/JET.E() > DET->TAU_EM_COLLIMATION ) &&
	     ( DET->NumTracks(TrackCentral,DET->PT_TRACK_TAU,JET.Eta(),JET.Phi()) == 1 )
	     ) {
	    elementTauJet = (TRootTauJet*) branchTauJet->NewEntry();
	    elementTauJet->Set(JET);
	  } // if tau jet
	} // if JET.eta < tracker - tau_cone : Tau jet identification
      } // for itJet : loop on all jets
      
      treeWriter->Fill();
      // Add here the trigger
      // Should test all the trigger table on the event, based on reconstructed objects
    } // Loop over all events
  treeWriter->Write();
  
  cout << "** Exiting..." << endl;
  
  delete treeWriter;
  delete treeReader;
  delete DET;
  if(converter) delete converter;

  todo("TODO");
}