// -*- C++ -*-
//
// Package:    DiscriminationPower
// Class:      DiscriminationPower
// 
/**\class DiscriminationPower DiscriminationPower.cc Discriminator/DiscriminationPower/src/DiscriminationPower.cc

 Description: <one line class summary>

 Implementation:
     <Notes on implementation>
*/
//
// Original Author:  Loic QUERTENMONT
//         Created:  Thu May 22 12:19:07 CEST 2008
// $Id$
//


#include <memory>

#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/EDAnalyzer.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/Utilities/interface/Exception.h"

#include "DataFormats/TrackReco/interface/DeDxData.h"


#include "Geometry/CommonDetUnit/interface/GeomDetUnit.h"
#include "Geometry/CommonDetUnit/interface/GeomDetType.h"
#include "Geometry/CommonTopologies/interface/StripTopology.h"
#include "DataFormats/GeometrySurface/interface/TrapezoidalPlaneBounds.h"
#include "DataFormats/GeometrySurface/interface/RectangularPlaneBounds.h"

#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
#include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
#include "Geometry/TrackerGeometryBuilder/interface/PixelGeomDetUnit.h"
#include "Geometry/TrackerNumberingBuilder/interface/GeometricDet.h"
#include "Geometry/CommonDetUnit/interface/TrackingGeometry.h"

#include "DataFormats/SiStripCluster/interface/SiStripClusterCollection.h"

#include "DataFormats/MuonReco/interface/Muon.h"
#include "DataFormats/MuonReco/interface/MuonFwd.h"

#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/TrackReco/interface/TrackFwd.h"
#include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h"
#include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
#include "DataFormats/SiStripDetId/interface/SiStripSubStructure.h"
#include "DataFormats/DetId/interface/DetId.h"
#include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
#include "DataFormats/SiStripDetId/interface/TECDetId.h"
#include "DataFormats/SiStripDetId/interface/TIBDetId.h"
#include "DataFormats/SiStripDetId/interface/TIDDetId.h"
#include "DataFormats/SiStripDetId/interface/TOBDetId.h"

#include "TrackingTools/PatternTools/interface/Trajectory.h"
#include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"

#include "SimDataFormats/HepMCProduct/interface/HepMCProduct.h"
#include "PhysicsTools/Utilities/interface/deltaR.h"

#include "DataFormats/TrackReco/interface/TrackFwd.h"

#include "FWCore/ServiceRegistry/interface/Service.h"
#include "PhysicsTools/UtilAlgos/interface/TFileService.h"

#include "myTools/HistoricizedEventProducer/interface/HistoricizedEvent.h"
#include "DataFormats/Common/interface/TriggerResults.h"



#include "TFile.h"
#include "TObjString.h"
#include "TString.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TProfile.h"
#include "TF1.h"
#include "TTree.h"
#include "TROOT.h"

#include <ext/hash_map>


using namespace edm;
using namespace reco;
using namespace std;
using namespace __gnu_cxx;


class DiscriminationPower : public edm::EDAnalyzer {
   public:
      explicit DiscriminationPower(const edm::ParameterSet&);
      ~DiscriminationPower();


   private:
      virtual void beginJob(const edm::EventSetup&) ;
      virtual void analyze(const edm::Event&, const edm::EventSetup&);
      virtual void endJob() ;

      const edm::EventSetup* iSetup_;
      const edm::Event*      iEvent_;

      edm::Service<TFileService> tfs;


      TFile*       HistoFile;
      std::string  HistoFileName;

      InputTag     m_dEdxDiscrimTag;
      InputTag     m_tracksTag;


      TH1F*        HNTracks;
      TH1F*        HEstimators;
      TH1F*        HNTracksHits[40];
      TH1F*        HEstimatorsHits[40];


      TH2F*        HP;
      TH2F*        HEta;
      TH2F*        HNHits;

      TH1F*        HCuts;

  TTree*       MyTree;
  float estimator;
  unsigned int NOM;
  int dbx_past;
  int dbx_future;
  int dbx;
  unsigned int eventNumber;
  float momentum;
  float pt;
  float chi2ndof;
  float eta;
  unsigned int isMatchedSA;
  unsigned int isMatchedGL;
  unsigned int isMatchedGLT;

  bool IsMatched(reco::TrackRef, std::vector<reco::MuonRef>, double);

      float        Min;
      float        Max;
      int          NBins;
      unsigned int NTracks;
      unsigned int NTracksHits[40];


      double       MinTrackMomentum;
      double       MaxTrackMomentum;
      double       MinTrackEta;
      double       MaxTrackEta;
      int          MinTrackNHits;
      int          MaxTrackNHits;
      bool         OnlyHSCP;


};

DiscriminationPower::DiscriminationPower(const edm::ParameterSet& iConfig)
{
   m_tracksTag         = iConfig.getParameter<InputTag>("tracks");
   m_dEdxDiscrimTag    = iConfig.getParameter<InputTag>("dEdxDiscrim");  

   HistoFileName       = iConfig.getParameter<string>("HistoFileName");

 
   Min                 = iConfig.getParameter<double>("Min");
   Max                 = iConfig.getParameter<double>("Max");
   NBins               = iConfig.getParameter<int  >("NBins");

   MinTrackMomentum    = iConfig.getUntrackedParameter<double>  ("minTrackMomentum"   ,  1.0);
   MaxTrackMomentum    = iConfig.getUntrackedParameter<double>  ("maxTrackMomentum"   ,  99999.0);
   MinTrackEta         = iConfig.getUntrackedParameter<double>  ("minTrackEta"        , -5.0);
   MaxTrackEta         = iConfig.getUntrackedParameter<double>  ("maxTrackEta"        ,  5.0);
   MinTrackNHits       = iConfig.getUntrackedParameter<int>     ("minTrackNHits"      ,  0  );
   MaxTrackNHits       = iConfig.getUntrackedParameter<int>     ("maxTrackNHits"      ,  50 );

   OnlyHSCP            = iConfig.getUntrackedParameter<bool>    ("OnlyHSCP"           ,  false);
}


DiscriminationPower::~DiscriminationPower()
{
}

void
DiscriminationPower::beginJob(const edm::EventSetup& iSetup)
{
   TH1::AddDirectory(kTRUE);

   iSetup_                  = &iSetup;

   HNTracks    = tfs->make<TH1F>("NTracks"   , "NTracks"    , NBins, Min, Max);
   HEstimators = tfs->make<TH1F>("Estimators", "Estimators" , NBins, Min, Max);

   for(unsigned int i=0;i<40;i++){
      char temp[255];
      sprintf(temp,"NTracksHits%2i",i);
      HNTracksHits[i]    = tfs->make<TH1F> (temp   , temp    , NBins, Min, Max);

      sprintf(temp,"EstimatorsHits%2i",i);
      HEstimatorsHits[i] = tfs->make<TH1F> (temp  , temp , NBins, Min, Max);

      NTracksHits[i]     = 0;
   }

   HP                    = tfs->make<TH2F> ("P"         , "P"     , 1000,0,100,NBins, Min, Max);
   HEta                  = tfs->make<TH2F> ("Eta"       , "Eta"   , 100 ,0,4,  NBins, Min, Max);
   HNHits                = tfs->make<TH2F> ("Nhits"     , "NHits" , 30, 0,30,  NBins, Min, Max);
   HCuts                 = tfs->make<TH1F> ("Cuts"      , "Cuts"  , 50, 0, 20);

   NTracks = 0;

   TTree::SetMaxTreeSize(1000*Long64_t(2000000000)); // authorize Trees up to 2 Terabytes
   MyTree = tfs->make<TTree> ("MyTree","MyTree");
   MyTree->Branch("eventNumber",&eventNumber,"eventNumber/I");
   MyTree->Branch("dbx",&dbx,"dbx/I");
   MyTree->Branch("NOM",&NOM,"NOM/I");
   MyTree->Branch("momentum",&momentum,"momentum/F");
   MyTree->Branch("pt",&pt,"pt/F");
   MyTree->Branch("chi2ndof",&chi2ndof,"chi2ndof/F");
   MyTree->Branch("eta",&eta,"eta/F");
   MyTree->Branch("estimator",&estimator,"estimator/F");
   MyTree->Branch("isMatchedSA",&isMatchedSA,"isMatchedSA/I");
   MyTree->Branch("isMatchedGL",&isMatchedGL,"isMatchedGL/I");
   MyTree->Branch("isMatchedGLT",&isMatchedGLT,"isMatchedGLT/I");

   
}

void
DiscriminationPower::endJob() 
{
/*
   for(int j=0;j<HNTracks->GetXaxis()->GetNbins();j++){
      double eff = HNTracks->GetBinContent(j)/NTracks;
//      double err = (eff*(1-eff))/NTracks;
      double err = sqrt(eff*(1-eff)/NTracks);
      HNTracks->SetBinContent(j,eff);
      HNTracks->SetBinError  (j,err);

      for(unsigned int i=0;i<40;i++){
         double eff = HNTracksHits[i]->GetBinContent(j)/NTracks;
//         double err = (eff*(1-eff))/NTracks;
         double err = sqrt(eff*(1-eff)/NTracks);
         HNTracksHits[i]->SetBinContent(j,eff);
         HNTracksHits[i]->SetBinError  (j,err);
      }
   }
*/
}

bool DiscriminationPower::IsMatched(reco::TrackRef track, std::vector<reco::MuonRef> muons, double drmax)
{
  bool debug = false;
  if (debug) cout << "I'm inside the matching routine." << endl;

  bool matched = true;
  double invpt_tk = 1./track->pt();

  double pt_mu = -1.;
  for (unsigned int i=0; i<muons.size(); i++){
    bool good = true;
    reco::MuonRef muon = muons[i];

    pt_mu = muon->pt();
    if (pt_mu > 30.) {
      // momentum matching
      double invpt_mu = 1./muon->pt();
      double diff = invpt_tk - invpt_mu;
      if (debug) cout << "invpt_tk - invpt_mu = " << diff << endl;
      if (fabs(diff)>0.005) good = false;
      
      // angular matching
      double dr = deltaR(track->momentum(),muon->momentum());
      if (debug) cout << "dr = " << dr << endl;
      if (dr>drmax) good = false;
    }
    matched *= good;
  }

  return matched;
}

void
DiscriminationPower::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
{
  bool debug = false;
  if (debug) cout << "I'm in DiscriminationPower::analyze!" << endl;

   Handle<ValueMap<DeDxData> >          dEdxTrackHandle;
   try {
     iEvent.getByLabel(m_dEdxDiscrimTag, dEdxTrackHandle);
   } catch (...) {;}
   const ValueMap<DeDxData> dEdxTrack = *dEdxTrackHandle.product();


   edm::Handle<reco::TrackCollection> trackCollectionHandle;
   try {
     iEvent.getByLabel(m_tracksTag,trackCollectionHandle);
   } catch (...) {;}
   reco::TrackCollection trackCollection = *trackCollectionHandle.product();


   edm::Handle<reco::MuonCollection> muonCollectionHandle;
   try {
     iEvent.getByLabel("muons",muonCollectionHandle);
   } catch (...) {;}
   reco::MuonCollection muonCollection = *muonCollectionHandle.product();

   Handle<HepMCProduct> HepMC;
   try {
     if(OnlyHSCP) iEvent.getByLabel("source", HepMC);
   } catch (...) {;}

   // HISTORY:

   dbx_past=-9999; dbx_future=-9999;
   dbx=0;

   Handle<HistoricizedEvent> hEvent_past;
   iEvent.getByLabel("consecutiveHEs",hEvent_past);
   dbx_past = (*hEvent_past).DeltaBX();
   Handle<HistoricizedEvent> hEvent_future;
   iEvent.getByLabel("futureHEs",hEvent_future);
   dbx_future = (*hEvent_future).DeltaBX();

   eventNumber = iEvent.id().event();

   std::vector<edm::Handle<edm::TriggerResults> > prod;
   iEvent.getManyByType(prod);

   unsigned int numbits = 2 ;
   if(numbits != prod[0]->size()) {
     std::cerr << "prod[0] should have " << numbits
	       << " entries, got " << prod[0]->size() << " in TriggerResults\n";
     abort();
   }

   bool past = prod[0]->at(0).accept();
   bool future = prod[0]->at(1).accept();
   if (past) {
     cout << "past" << endl;
     cout << "dbx_past = " <<dbx_past << endl;
     cout << "dbx_future = " <<dbx_future << endl;
     dbx=-1*dbx_past;
   }
   if (future) {
     cout << "future" << endl;
     cout << "dbx_past = " <<dbx_past << endl;
     cout << "dbx_future = " <<dbx_future << endl;
     dbx=dbx_future;
   }
   if (past && future) cout << "past & future!" << endl;

   // MUON LOOP:
   std::vector<reco::MuonRef> standaloneMuons;
   std::vector<reco::MuonRef> globalMuons;
   std::vector<reco::MuonRef> globalMuonsTight;
   for(unsigned int i=0; i<muonCollection.size(); i++){
     if (debug) cout << "I'm in the muon loop!" << endl;
     reco::MuonRef muon  = reco::MuonRef( muonCollectionHandle, i );
     bool standaloneMuon = muon->isGood(reco::Muon::AllStandAloneMuons);
     bool globalMuon = muon->isGood(reco::Muon::AllGlobalMuons);
     bool globalMuonTight = muon->isGood(reco::Muon::GlobalMuonPromptTight);
     if (standaloneMuon) {
       if (debug) cout << "This is a stand-alone muon" << endl;
       standaloneMuons.push_back(muon);
     }
     if (globalMuon) {
       if (debug) cout << "This is a global muon" << endl;
       globalMuons.push_back(muon);
     }
     if (globalMuonTight) {
       if (debug) cout << "This is a global muon tight" << endl;
       globalMuonsTight.push_back(muon);
     }
   }
   if (debug) cout << "standaloneMuons.size() = " << standaloneMuons.size() << endl;
   if (debug) cout << "globalMuons.size() = " << globalMuons.size() << endl;
   if (debug) cout << "globalMuonsTight.size() = " << globalMuonsTight.size() << endl;

   ///// muon wiki: https://twiki.cern.ch/twiki/bin/view/CMS/WorkBookMuonAnalysis

   // TRACK LOOP:
   for(unsigned int i=0; i<trackCollection.size(); i++){
     reco::TrackRef track  = reco::TrackRef( trackCollectionHandle, i );
     if (debug) cout << "I'm in the track loop!" << endl;
     estimator = -1.; NOM = 0; momentum = -1.; eta = -99.; pt = -1.; chi2ndof = -1.; isMatchedSA = 0; isMatchedGL = 0; isMatchedGLT = 0;
     momentum = track->p();
     pt = track->pt();
     eta = track->eta();
     chi2ndof = track->chi2()/track->ndof();
     if (debug) cout << "p = " << momentum << endl;
     estimator       = dEdxTrack[track].dEdx();
     if (debug) cout << "estimator = " << estimator << endl;
     NOM      = dEdxTrack[track].numberOfMeasurements();
     if(NOM>=40)NOM=39;
     isMatchedSA = IsMatched(track, standaloneMuons, 0.1);
     isMatchedGL = IsMatched(track, globalMuons, 0.0001);
     isMatchedGLT = IsMatched(track, globalMuonsTight, 0.0001);
     
     if(OnlyHSCP){
         bool GoodTrack = false;
         HepMC::GenEvent* myGenEvent = new  HepMC::GenEvent(*(HepMC->GetEvent()));
         for ( HepMC::GenEvent::particle_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end(); ++p ) {
            if ( abs((*p)->pdg_id()) > 100000 && (*p)->status() == 3 ){
                const Track      tracktmp = *track;
                double Dist = deltaR(tracktmp, (*p)->momentum());
                if(Dist<0.3)GoodTrack = true;

            }
         }
        delete myGenEvent;
        if(!GoodTrack){continue;}
      }

     if(eta<MinTrackEta || eta>MaxTrackEta){HCuts->Fill(4);continue;}
     if(momentum<MinTrackMomentum || momentum>MaxTrackMomentum){HCuts->Fill(6);continue;}
     if((int)NOM<MinTrackNHits ||(int)NOM>MaxTrackNHits){HCuts->Fill(8);continue;}
     if(estimator<0){printf("Negativ Estimator Cut\n"); HCuts->Fill(2); continue;}

     HCuts->Fill(16);

     HEstimators->Fill(estimator);
     HP         ->Fill(momentum,estimator);
     HEta       ->Fill(fabs(eta),estimator);
     HNHits     ->Fill(NOM,estimator);

     if (momentum > 5) MyTree->Fill();

     NTracks++;
     HEstimatorsHits[NOM]->Fill(estimator);
     NTracksHits[NOM]++;


     for(int j=0;j<HNTracks->GetXaxis()->GetNbins();j++){
	 float bincenter = HNTracks->GetXaxis()->GetBinCenter(j);
//         if(estimator>bincenter) HNTracks->Fill(bincenter);
         if(estimator>bincenter){
//             HNTracks->SetBinContent(j,HNTracks->GetBinContent(j)+1);
             HNTracks->Fill(bincenter);


//             HNTracksHits[NOM]->SetBinContent(j,HNTracksHits[NOM]->GetBinContent(j)+1);
		HNTracksHits[NOM]->Fill(bincenter);
         }         
     }

//      printf("%i - %f\n", track.key(),dEdxTrack[i].second);
//      printf("%f\n",dEdxTrack[i].second);
   }

}

//define this as a plug-in
DEFINE_FWK_MODULE(DiscriminationPower);