#include "Utilities/Configuration/interface/Architecture.h"
//COBRA
#include "CARF/G3Event/interface/G3EventProxy.h"
#include "CARF/SimEvent/interface/SimEvent.h"
#include "CARF/BaseSimEvent/interface/SimVertex.h"
#include "CARF/Reco/interface/RecQuery.h"
#include "CARF/Reco/interface/RecCollection.h"
#include "CARF/Reco/interface/FilteredRecCollection.h"
#include "CARF/Reco/interface/RecConfig.h"
// MC
#include "GeneratorInterface/Particle/interface/RawParticle.h"
#include "GeneratorInterface/HepEvent/interface/RawHepEvent.h"
#include "GeneratorInterface/RecParticle/interface/RawRecParticle.h"
// muons
#include "Muon/PersistentMuon/interface/RecMuon.h"
#include "MuonReco/MuonIsolation/interface/MuIsoBaseAlgorithm.h"
#include "MuonReco/MuonIsolation/interface/MuIsoByTrackerPt.h"
#include "MuonReco/MuonIsolation/interface/MuIsoByCaloEt.h"
#include "MuonReco/MuonReconstruction/interface/GlobalMuonReconstructor.h"
#include "Muon/MuonInternalReco/interface/MuonStandAloneReconstructor.h"
// electrons
#include "ElectronPhoton/EgammaElectron/interface/ElectronCandidate.h"
#include "ElectronPhoton/EgammaOfflineReco/interface/OfflineElectronReco.h"
// Jets
#include "Jets/JetFramework/interface/JetFinderCARFFacade.h"
#include "Jets/JetTypes/interface/RecJet.h"
#include "Jets/JetInputGenerators/interface/JetFinderGeneratorInput.h"
#include "Jets/JetInputGenerators/interface/JetFinderCaloRecHitInput.h"
#include "Jets/JetInputGenerators/interface/HLTJetFinderEcalPlusHcalTowerInput.h"
#include "Jets/JetInputGenerators/interface/JetFinderEcalPlusHcalTowerInput.h"
#include "Jets/JetAlgos/interface/KtJetAlgorithm.h"
#include "Jets/JetAlgos/interface/IterativeConeAlgorithm.h"
//MET
#include "MET/MetType/interface/RecMET.h"
#include "JetMetAnalysis/JetMetRoot/interface/JetMetRootEvent.h"
#include "JetMetAnalysis/JetMetRoot/interface/JetMetRootOutputSetup.h"
#include "JetMetAnalysis/JetMetParticles/interface/JetMetParticleSetup.h"
#include "JetMetAnalysis/JetMetParticles/interface/JetMetJetFinderSetup.h"
#include "JetMetAnalysis/JetMetRoot/interface/JetMetRootOutputSetup.h"

//Famos
#include "FamosGeneric/FamosManager/interface/FamosEventMgr.h"
#include "FamosGeneric/FamosUtils/interface/FamosHistos.h"
#include "FamosGeneric/FamosSimEvent/interface/FSimEvent.h"
#include "FamosGeneric/FamosSimEvent/interface/FSimTrack.h"
#include "FamosGeneric/FamosSimEvent/interface/FSimVertex.h"
#include "ORCAInterface/EgammaInterface/interface/FastBasicReco.h"
#include "FamosGeneric/FamosManager/interface/FamosSimulator.h"
#include "FamosGeneric/FamosEventReader/interface/FastGenLevelDetector.h"
#include "FamosGeneric/FamosEventReader/interface/FamosPUGenerator.h"
#include "FamosGeneric/FamosSimEvent/interface/FBaseSimEvent.h"
//ORCA
#include "ORCAInterface/FastPersistentJet/interface/FastPersistentJetFinder.h"
#include "ORCAInterface/JetInputGenerators/interface/JetFinderFastEcalPlusHcalTowerInput.h"
#include "ORCAInterface/JetInputGenerators/interface/JetFinderFastCaloRecHitInput.h"
#include "ORCAInterface/JetInputGenerators/interface/JetFinderFastGeneratorInput.h"
#include "ORCAInterface/FastEcalPlusHcalTower/interface/FastEcalPlusHcalTowerBuilder.h"

#include "CLHEP/Units/PhysicalConstants.h"
#include "CLHEP/Vector/LorentzVector.h"
//
#include <iostream>
#include <iomanip>
#include <map>

#include "TROOT.h"
#include "TTree.h"
#include "TFile.h"
#include "TH1F.h"



using namespace std;



/**
 * simple test FamosEvtMgr with PileUP events
 * \author Valery Zhukov
 * \date 20-May-2005
 */


class testFamosWithPU : private Observer<G3EventProxy *> {
public:
testFamosWithPU();
virtual ~testFamosWithPU();
virtual void analyzeEv(G3EventProxy * ev);
void store();
private:
//// don't change the name "upDate" - this method is mandatory.
void upDate(G3EventProxy * ev) { if (ev!=0) analyzeEv(ev);}
int theEvent;
FamosEventMgr* myEventMgr;
  TH1F* h_npu;
  TH1F* h_nmcpart;
  TH1F* h_nmcm;
  TH1F* h_nm;
  TH1F* h_ptm;
  TH1F* h_ne;
  TH1F* h_pte;
  TH1F* h_nj;
  TH1F* h_ej;
  TH1F* h_met;
  TH1F* h_ntrck;
  bool bstore;
};

testFamosWithPU::testFamosWithPU() {
 cout<<">>>>>>testFamosWithPU: construct "<<endl; 
Observer<G3EventProxy*>::init();
myEventMgr = FamosEventMgr::instance();
RecoRegistry::current().debugDict();
theEvent=0;
 bstore=false;
// root
  h_npu = new TH1F("h_npu","h_npu",100,0.,100.);
  h_nmcpart = new TH1F("h_nmcpart","h_nmcpart",1000,0.,2000.);
  h_nmcm = new TH1F("h_nmcm","h_nmcm",100,0.,100.);
  h_nm = new TH1F("h_nm","h_nm",100,0.,100.);
  h_ptm = new TH1F("h_ptm","h_ptm",200,0.,200.);
  h_ne = new TH1F("h_ne","h_ne",100,0.,100.);
  h_pte = new TH1F("h_pte","h_pte",200,0.,200.);
  h_nj = new TH1F("h_nj","h_nj",100,0.,100.);
  h_ej = new TH1F("h_ej","h_ej",1000,0.,1000.);
  h_met = new TH1F("h_met","h_met",100,0.,500.);
  h_ntrck = new TH1F("h_ntrck","h_ntrck",1000,0.,1000.);
  // check PUgenerator
  if(myEventMgr->theFamosPUGenerator()->isActive()) {
 const char* pufile= myEventMgr->theFamosPUGenerator()->pufile();
 cout<<"---PUGenerator is active and using "<<pufile<<endl;
  } else 
    cout<<"---PUGenerator is not activated "<<endl;

cout<<"<<<<<<<<<testFamosWithPU  construct done"<<endl;
}

void testFamosWithPU::store() {
 TFile* theFile=new TFile("out.root","RECREATE");
 theFile->cd();
  h_npu->Write();
  h_nmcpart->Write();
  h_nmcm->Write();
  h_nm->Write();
  h_ptm->Write();
  h_ne->Write();
  h_pte->Write();
  h_nj->Write();
  h_ntrck->Write();
  h_met->Write();
  theFile->Close();
  bstore=true;
}


testFamosWithPU::~testFamosWithPU() {
  if(!bstore) store();
}

// analyze event called by upDate
void testFamosWithPU::analyzeEv(G3EventProxy * ev) {
  cout<<">>>>>start analysis event="<<theEvent<<endl;
  if(theEvent==5000) store();
  if(myEventMgr->theFamosPUGenerator()->isActive()) {
  int npuev= myEventMgr->theFamosPUGenerator()->nPUev();
  cout<<"==N PU events="<<npuev<<endl;
  h_npu->Fill(npuev);
  }
  // MC
FSimEvent& aSimEvent = myEventMgr->simSignal();
RawHepEvent& aHepEvent = myEventMgr->event();
 int nmcpart= aHepEvent.n();
 int nsev=aSimEvent.nGenParts();
 cout<<" ==N MC particles Nsig"<<nmcpart<<" Ntot "<<nsev<<endl;
 h_nmcpart->Fill(nsev);
 int Nm=0;
 genpart_container* agenpart=aSimEvent.genparts();
  for(genpart_iterator it=agenpart->begin();it!=agenpart->end();++it) {
 int part_id = abs((*it).pid());
if(part_id==13) {
  HepLorentzVector p=(*it).fourmomentum();
 cout<<" --MC muon "<<Nm<<" pt="<<p.perp()<<" eta="
 <<  p.eta()<<" sgn="<<(*it).pid()<<endl; 
    Nm++;  
     }
  }
 h_nmcm->Fill(Nm);
cout<<" ==N MC muons "<<Nm<<endl;
  //reconstruction
 //tracks
RecQuery qtrcks("FamosCombinatorialTrackFinder");
RecCollection<TTrack>* thetracks = new RecCollection<TTrack>(qtrcks);
 int ntrck=thetracks->size();
cout<<" == N rec tracks Nrectrk="<<ntrck<<endl;
 h_ntrck->Fill(ntrck);
  //muons
RecQuery qmuons("FamosL3MuonReconstructor");
RecCollection<RecMuon>* themuons=new  RecCollection<RecMuon> (qmuons);
 int nms=themuons->size();
 cout<<" ==N rec muons Nrecm="<<nms<<endl;;
 GlobalVector p;
 int Nrecm=0;
  for (RecCollection<RecMuon>::const_iterator it =
 themuons->begin(); it != themuons->end(); ++it) 
 {
 TrajectoryStateOnSurface traj = ((*it).stateAtVertex().isValid()) ?
                                     (*it).stateAtVertex() : (*it).innermostState();
 p = traj.globalMomentum();
cout<<" --rec muon "<<Nrecm<<" pt="<<p.perp()<<" eta="
<<p.eta()<<" sgn="<< it->charge()<<endl;
 h_ptm->Fill(p.perp());
   Nrecm++;
 }
  h_nm->Fill(Nrecm);
  // electron
RecQuery qele("EGFElectron","1.0");
RecCollection<ElectronCandidate>* theelectrons  = new RecCollection<ElectronCandidate>(qele);
int nes=theelectrons->size();
 cout<<" ==N rec ele Nrecm="<<nes<<endl;
 int Nrece=0;
for (RecCollection<ElectronCandidate>::const_iterator it =theelectrons->begin(); 
it != theelectrons->end();++it) {
 const   TTrack* track=(*it).GetTRElectron();
TrajectoryStateOnSurface traj =(*track).impactPointState();
   p = traj.globalMomentum();
cout<<" --rec ele "<<Nrece<<" pt="<<p.perp()<<" eta="
<<p.eta()<<" sgn="<<traj.charge()<<endl;
 h_pte->Fill(p.perp());
  Nrece++;
}
 h_ne->Fill(Nrece);

 //jets
 RecQuery qjet("FamosPersistentJetFinder");
 qjet.setParameter<int>("JetAlgorithm",2);  // iterative cone
 qjet.setParameter<int>("JetRecom",4);
 qjet.setParameter("ConeCut",0.5);
 qjet.setParameter("ConeSeedEtCut",0.0);
 qjet.setParameter("JetInput","EcalPlusHcalTowerInput");
 qjet.setParameter("EcalPlusHcalTowerEt",0.5);
 qjet.setParameter("JetCalibration","no");
 RecCollection<RecJet>* thejets=new  RecCollection<RecJet>(qjet);
 int njs=thejets->size();
cout<<" ==N rec jets Nrecj="<<njs<<endl;;
 int Nj=0;
for (RecCollection<RecJet>::const_iterator it = 
thejets->begin(); it != thejets->end(); ++it) 
    { 
      cout << " --rec Jet "<<Nj<< " Etot="<<(*it).getEnergy()
	   <<" eta="<<(*it).getEta()<<endl;
      h_ej->Fill((*it).getEnergy());
	Nj++;
}
  h_nj->Fill(Nj);
  // met
RecQuery qmet("FamosMETfromCaloTower");
RecCollection<RecMET>* themet=new RecCollection<RecMET> (qmet);
 for (RecCollection<RecMET>::const_iterator it = themet->begin(); it != themet->end(); ++it) {
   cout<<" ==MET "<<(*it).getEt()<<" sumET="<<(*it).getSumEt()<<endl;
   h_met->Fill((*it).getEt());
 }

 cout<<"<<<analyzed done== "<<theEvent<<endl;
 theEvent++;
}

#include "Utilities/Notification/interface/PackageInitializer.h"
static PKBuilder<testFamosWithPU>  testPUBuilder("testFamosWithPU");