///////////////////////////////////////////////////////////////////////////////
//
// Summary of fPrintLevel values:
// ----------------------
// fPrintLevel =   1: print Electron block for all the events
// fPrintLevel = 101: print line per event with good electron
// fPrintLevel = 102: found matching non-trigger jet
///////////////////////////////////////////////////////////////////////////////

#include "TF1.h" 
#include "TCanvas.h"
#include "TPad.h"
#include "TText.h"

#include "Stntuple/obj/TStnHeaderBlock.hh"
#include "Stntuple/obj/TStnTrackBlock.hh"
#include "Stntuple/obj/TStnDBManager.hh"
#include "Stntuple/alg/TStntuple.hh"
#include "Stntuple/alg/TStnElectronID.hh"
#include "Stntuple/obj/TStnTriggerTable.hh"
#include "Stntuple/obj/TStnTrigger.hh"

#include "Stntuple/loop/TStnAna.hh"

#include "Stntuple/geom/TTrajectoryPoint.hh"
#include "Stntuple/geom/TSimpleExtrapolator.hh"

#include "TLumiMonModule.hh"

#include <iostream>

ClassImp(TLumiMonModule)
//_____________________________________________________________________________
TLumiMonModule::TLumiMonModule(const char* name, const char* title): 
  TStnModule(name,title)
{
  fListOfGoodElectrons  = new TObjArray(20);
  fTightID              = new TStnElectronID();
//-----------------------------------------------------------------------------
// modes:
//-----------------------------------------------------------------------------
  fCutMode           = 0; 
}

//_____________________________________________________________________________
TLumiMonModule::~TLumiMonModule() {

  fListOfGoodElectrons->Clear();
  delete fListOfGoodElectrons;

  delete fTightID;
}


//_____________________________________________________________________________
int TLumiMonModule::BeginJob() {
  // register the data blocks and book histograms

  RegisterDataBlock("ElectronBlock", "TStnElectronBlock", &fElectronBlock);
  RegisterDataBlock("TriggerBlock" , "TStnTriggerBlock" , &fTriggerBlock);
  RegisterDataBlock("TrackBlock"   , "TStnTrackBlock"   , &fTrackBlock);
  RegisterDataBlock("ZVertexBlock" , "TStnVertexBlock"  , &fVertexBlock);
  RegisterDataBlock("MetBlock"     , "TStnMetBlock"     , &fMetBlock);

  BookHistograms();

  return 0;
}


//_____________________________________________________________________________
int TLumiMonModule::BeginRun() {

  int rn = GetHeaderBlock()->RunNumber();
  TStntuple::Init(rn);

  TStntuple::InitListOfL3Triggers(this);

  return 0;
}

//_____________________________________________________________________________
void TLumiMonModule::FillEventHistograms(EventHist_t&     Hist) {

  int nemo;
  

  TStnHeaderBlock* h = GetHeaderBlock();

  Hist.fRunNumber->Fill(h->RunNumber());
  Hist.fFilterResult->Fill(GetPassed());

  float lumi = h->InstLum()/1.e30;

  Hist.fInstLumi  ->Fill(lumi);

  Hist.fMet0      ->Fill(fMetBlock->Met(0));
  Hist.fMet1      ->Fill(fMetBlock->Met(1));
  Hist.fMet0VsMet1->Fill(fMetBlock->Met(1),fMetBlock->Met(0),1);

  Hist.fNTracks   ->Fill(fNTracks[0]);

  Hist.fNTracksVsLumi[0]->Fill(lumi,fNTracks[0],1);
  Hist.fNTracksVsLumi[1]->Fill(lumi,fNTracks[1],1);
  Hist.fNTracksVsLumi[2]->Fill(lumi,fNTracks[2],1);

  Hist.fNV00VsLumi->Fill(lumi,fNV00,1);
  Hist.fNV12VsLumi->Fill(lumi,fNV12,1);

  Hist.fNGoodEle->Fill(fNGoodElectrons);
//-----------------------------------------------------------------------------
//  transverse mass histograms
//-----------------------------------------------------------------------------
  double metm[2], metx[2], mety[2], mt[2], q[2];

  for (int i=0; i<2; i++) {
    metm[i] = fMetBlock->Met(i);
    metx[i] = metm[i]*cos(fMetBlock->MetPhi(i));
    mety[i] = metm[i]*sin(fMetBlock->MetPhi(i));
  }

  if (fNGoodElectrons > 0) {
    nemo = fElectronBlock->NElectrons();
    for (int i1=0; i1<nemo; i1++) {
      TStnElectron*   e1 = fElectronBlock->Electron(i1);
      const TLorentzVector* m1 = e1->Momentum();
      if (fIdWord[i1] == 0) {
	
       	mt[0] = sqrt(2*(e1->Et()*metm[0] -m1->Px()*metx[0]-m1->Py()*mety[0]));
       	mt[1] = sqrt(2*(e1->Et()*metm[1] -m1->Px()*metx[1]-m1->Py()*mety[1]));
	
	Hist.fMt[0]->Fill(mt[0]);
	Hist.fMt[1]->Fill(mt[1]);
//-----------------------------------------------------------------------------
// distributions for isolation vs MET
//-----------------------------------------------------------------------------
	Hist.fIsoVsMet->Fill(fMetBlock->Met(1),e1->Iso());
	Hist.fIso1VsMet->Fill(fMetBlock->Met(1),e1->Iso()/e1->Et());
//-----------------------------------------------------------------------------
// distributions for the component of MET perpendicular to the electron 
// momentum
//-----------------------------------------------------------------------------
	q[0] = (metx[0]*m1->Py()-mety[0]*m1->Px())/m1->Pt();
	q[1] = (metx[1]*m1->Py()-mety[1]*m1->Px())/m1->Pt();

	Hist.fMetPerp[0]->Fill(q[0]);
	Hist.fMetPerp[1]->Fill(q[1]);
      }
    }
  }
}


//_____________________________________________________________________________
void TLumiMonModule::FillElectronHistograms(EleHist_t& Hist, TStnElectron* Ele) {

  TStnHeaderBlock* h = GetHeaderBlock();
  float lumi = h->InstLum()/1.e30;

  const TLorentzVector* mom = Ele->Momentum();

  Hist.fEt->Fill(Ele->Et());

  Hist.fEtVsEta->Fill(mom->Eta(),Ele->Et());
  Hist.fEOverP->Fill(Ele->EOverP());

  if (Ele->Charge() > 0) {
    Hist.fWedgeVsDelXPos->Fill(Ele->DelX(),Ele->SeedIPhi()+24*Ele->Side());
    Hist.fEOverPPos->Fill(Ele->EOverP());
  }
  else {
    Hist.fWedgeVsDelXNeg->Fill(Ele->DelX(),Ele->SeedIPhi()+24*Ele->Side());
    Hist.fEOverPNeg->Fill(Ele->EOverP());
  }

  Hist.fWedgeVsZCes->Fill(Ele->ZCes(),Ele->WedgeNumber());
  Hist.fWedgeVsDelZ->Fill(Ele->DelZ(),Ele->WedgeNumber());

  Hist.fEmfrVsDelZ->Fill(Ele->DelZ(),Ele->Emfr(),1.);

  Hist.fPhi->Fill(TVector2::Phi_0_2pi(mom->Phi()));
  Hist.fIso->Fill(Ele->Iso());
  Hist.fIso1->Fill(Ele->Iso()/Ele->Et());
  Hist.fHadEm->Fill(Ele->HadEm());
  Hist.fTrackPt->Fill(Ele->TrackPt());
  Hist.fTrackXCes->Fill(Ele->TrackPt());
  Hist.fTrackZCes->Fill(Ele->TrackPt());
  
  Hist.fXCes->Fill(Ele->XCes());
  Hist.fZCes->Fill(Ele->ZCes());

  Hist.fLshr->Fill (Ele->Lshr());
  Hist.fLshr2->Fill(Ele->Lshr2());

  Hist.fDelX->Fill(Ele->DelX());
  Hist.fDelZ->Fill(Ele->DelZ());

  Hist.fNCesStripClusters->Fill(Ele->NCesClusters(0));
  Hist.fNCesWireClusters->Fill(Ele->NCesClusters(1));
  Hist.fCesStripEnergy->Fill(Ele->CesEnergy(0));
  Hist.fCesWireEnergy->Fill(Ele->CesEnergy(1));
  Hist.fCesStripChi2->Fill(Ele->Chi2Strip());
  Hist.fCesWireChi2->Fill(Ele->Chi2Wire());
  Hist.fCesEStripVsEWire->Fill(Ele->CesEnergy(1),Ele->CesEnergy(0),1);
  Hist.fCesEStripVsEta[0]->Fill(Ele->CesEnergy(1),mom->Eta(),1);
  Hist.fCesEStripVsEt [0]->Fill(Ele->CesEnergy(1),mom->E(),1);

  double estrip_corr = Ele->CesEnergy(1)*Ele->Et()/mom->E();
  Hist.fCesEStripVsEta[1]->Fill(estrip_corr,mom->Eta(),1);
  Hist.fCesEStripVsEt [1]->Fill(estrip_corr,mom->E(),1);
//-----------------------------------------------------------------------------
//  fill track histograms
//-----------------------------------------------------------------------------
  //  TLorentzVector* tmom;
  TStnTrack*      trk;
  int it = Ele->TrackNumber(); 
  if (it >= 0) {
    trk   = fTrackBlock->Track(it);
    //    tmom  = trk->Momentum();
    Hist.fTIso->Fill(trk->Iso4());
    Hist.fTrackD0->Fill(trk->D0());
    Hist.fTrackD0Corr->Fill(TStntuple::CorrectedD0(trk));
    Hist.fTrackZ0->Fill(trk->Z0());
    Hist.fTrackChi2->Fill(trk->Chi2()/(trk->NCotHitsTot()+trk->NSvxHits()-4.9999));
    Hist.fTrackCotChi2->Fill(trk->Chi2Cot()/(trk->NCotHitsTot()-4.9999));
    Hist.fTrackNCotHitsTot->Fill(trk->NCotHitsTot());
    Hist.fTrackNCotHitsAx->Fill(trk->NCotHitsAx());
    Hist.fTrackNCotHitsSt->Fill(trk->NCotHitsSt());

    Hist.fNSvxHits->Fill(trk->NSvxHits());
    Hist.fNSvxRPhiHits->Fill(trk->NSvxRPhiHits());
    Hist.fNSvxSASHits->Fill(trk->NSvxSASHits());
    Hist.fNSvxZHits->Fill(trk->NSvxZHits());
  }
  else {
    trk = 0;
    Hist.fTIso->Fill(1e6);
    Hist.fTrackD0->Fill(1e6);
    Hist.fTrackD0Corr->Fill(1.e6);
    Hist.fTrackZ0->Fill(1.e6);
    Hist.fTrackChi2->Fill(1.e6);
    Hist.fTrackCotChi2->Fill(1.e6);

    Hist.fTrackNCotHitsTot->Fill(-1);
    Hist.fTrackNCotHitsAx ->Fill(-1);
    Hist.fTrackNCotHitsSt ->Fill(-1);

    Hist.fNSvxHits    ->Fill(-1);
    Hist.fNSvxRPhiHits->Fill(-1);
    Hist.fNSvxSASHits ->Fill(-1);
    Hist.fNSvxZHits   ->Fill(-1);
  }
//-----------------------------------------------------------------------------
//  lumi histograms
//-----------------------------------------------------------------------------
  if (trk) {
    Hist.fNCotHitsTotVsLumi->Fill(lumi,trk->NCotHitsTot());
    Hist.fNCotHitsAxVsLumi ->Fill(lumi,trk->NCotHitsAx ());
    Hist.fNCotHitsStVsLumi ->Fill(lumi,trk->NCotHitsSt ());
    Hist.fIsoVsLumi        ->Fill(lumi,Ele->Iso        ());
  }
//-----------------------------------------------------------------------------
//  ID bits
//-----------------------------------------------------------------------------
  int id_word = Ele->IDWord();
  for (int bit=0; bit<32; bit++) {
    if (((id_word >> bit) & 0x1) == 1) {
      Hist.fFailedBits->Fill(bit);
    }
  }
}


//_____________________________________________________________________________
int TLumiMonModule::Event(int ientry) {
  
  //  TLorentzVector* mom;
  Double_t        met1, pt, lam0;
  Int_t           nemo;
  TStnTrack       *seed;
//-----------------------------------------------------------------------------
//  check the trigger path
//-----------------------------------------------------------------------------
  fTriggerBlock->GetEntry(ientry);
  int passed = TStntuple::CheckL3TriggerPath(fTriggerBlock,fListOfL3Triggers);

  if (! passed) return 0;
//-----------------------------------------------------------------------------
//  read the necessary data
//-----------------------------------------------------------------------------
  fElectronBlock->GetEntry(ientry);
  fTrackBlock->GetEntry(ientry);
  fVertexBlock->GetEntry(ientry);
  fMetBlock->GetEntry(ientry);

  nemo = fElectronBlock->NElectrons();
  met1 = fMetBlock->Met(1);
//-----------------------------------------------------------------------------
//  dont forget to clear the lists
//-----------------------------------------------------------------------------
  fListOfGoodElectrons->Clear();
  fNGoodElectrons  = 0;

  for (int i=0; i<nemo; i++) {
    TStnElectron* ele = fElectronBlock->Electron(i);
    seed   = TStntuple::ElectronTrack(ele,fTrackBlock);
    ele->SetTrack(seed);
//-----------------------------------------------------------------------------
// check tight ID cuts
//-----------------------------------------------------------------------------
    fIdWord[i] = fTightID->IDWord(ele);
    ele->SetIdWord(fIdWord[i]);
//-----------------------------------------------------------------------------
// set 0: histograms for ALL EM objects
//-----------------------------------------------------------------------------
    FillElectronHistograms(fHist.fEle[0],ele);

    if (fIdWord[i] == 0) { 
//-----------------------------------------------------------------------------
// set 1: all TIGHT electrons
//-----------------------------------------------------------------------------
      fNGoodElectrons++;
      fListOfGoodElectrons->Add(ele);
      FillElectronHistograms(fHist.fEle[1],ele);

      if (met1 > 25.) {
	FillElectronHistograms(fHist.fEle[2],ele);
      }
    }
  }
//-----------------------------------------------------------------------------
//  track and vertex counting
//-----------------------------------------------------------------------------
  fNV00 = fVertexBlock->NVertices(0);
  fNV12 = fVertexBlock->NVertices(12);

  fNTracks[0] = fTrackBlock->NTracks();
  TStnTrack* trk;

  fNTracks[1] = 0;
  fNTracks[2] = 0;

  for (int i=0; i<fNTracks[0]; i++) {
    trk = fTrackBlock->Track(i);
    lam0 = trk->Lam0();
    pt   = trk->Momentum()->Pt();
    if ((pt > 0.4) && (TMath::Abs(lam0) < 1)) {
      fNTracks[1] ++;
      if ((trk->NStSeg() > 2) && (trk->NAxSeg() > 2)) {
	fNTracks[2]++;
      }
    }
  }
//-----------------------------------------------------------------------------
//  event histograms
//-----------------------------------------------------------------------------
  if (fNGoodElectrons > 0) {
    FillEventHistograms(fHist.fEvent[0]);

    if (fMetBlock->Met(1) > 25.) {
//-----------------------------------------------------------------------------
//  new W->enu candidate: look at it
//-----------------------------------------------------------------------------
      FillEventHistograms(fHist.fEvent[1]);
    }
  }

  return 0;			       
}


//_____________________________________________________________________________
int TLumiMonModule::EndJob() {
  return 0;
}