/* File l3Metefficiency.cpp
 *
 * Created       : June 2007
 * Author        : Vesna Cuplov
 * Purpose       : MET studies
 */

#include "tmb_tree/TMBEMCluster.hpp"
#include "cafe/Processor.hpp"
#include "cafe/Collection.hpp"
#include "cafe/Config.hpp"
#include "TBranch.h"
#include "TTree.h"
#include "TLeaf.h"
#include <fstream>
#include <iostream>
#include <vector>
#include "L3MetEff/l3Metefficiency.hpp"
#include "TGraphAsymmErrors.h"

using namespace std;
using namespace cafe;

// Constructor, destructor: 
l3Metefficiency::l3Metefficiency(const char *name) : cafe::Processor(name)
{

}

l3Metefficiency::~l3Metefficiency()
{

}

void l3Metefficiency::begin()
{
  total_events =0.;
  total_events_run_cut = 0. ;
  total_events_offlinecut =0.;
  events_all_cuts =0.;

  // CUTS on Level3 MEt
  Offline_events_L3METcut_10Gev =0. ;
  Offline_events_L3METcut_15Gev =0. ;
  Offline_events_L3METcut_20Gev =0. ;
  Offline_events_L3METcut_25Gev =0. ;

  Offline_events_L3noICDMETcut_10Gev =0. ;
  Offline_events_L3noICDMETcut_15Gev =0. ;
  Offline_events_L3noICDMETcut_20Gev =0. ;
  Offline_events_L3noICDMETcut_25Gev =0. ;

  // Offline events that passed the trigger: E1SHT25
  Offline_events_E1SHT25 =0. ;
  fired_triggers = new fstream("fired_triggers.txt",ios::out);

  // Instantaneous Luminosity
  h_instLumi_E1SHT25 = new TH1F("h_instLumi_E1SHT25"," Instantaneous luminosity after E1SHT25", 30,0.0,5.0);
  h_instLumi_L3MET10 = new TH1F("h_instLumi_L3MET10"," Instantaneous luminosity L3MET10", 30,0.0,5.0);
  h_instLumi_L3MET15 = new TH1F("h_instLumi_L3MET15"," Instantaneous luminosity L3MET15 ", 30,0.0,5.0);
  h_instLumi_L3MET20 = new TH1F("h_instLumi_L3MET20"," Instantaneous luminosity L3MET20 ", 30,0.0,5.0);
  h_instLumi_L3MET25 = new TH1F("h_instLumi_L3MET25"," Instantaneous luminosity L3MET25 ", 30,0.0,5.0);

  // Offline histograms:
  // Transverse W mass and Wgroup MET
  h_Tmass = new TH1F("h_Tmass","T Mass", 150,0.0,150.0);
  h_Offline_MEt = new TH1F("h_W_MEt"," ", 100,0.0,100.0);
  h_OfflineMEt_E1SHT25    = new TH1F("h_OfflineMEt_E1SHT25","Offline MEt passed E1SHT25",100,0.0,100.0);


  //Level 3 MEt histograms:
  h_L3MEt = new TH1F("h_L3MEt","L3 MEt", 100,0.0,100.0);
  h_L3MEt_E1SHT25 = new TH1F("h_L3MEt_E1SHT25","L3 MEt E1SHT25", 100,0.0,100.0);

  h_L3MEt_noicd = new TH1F("h_L3MEt_noicd","L3 MEt noICD", 100,0.0,100.0);
  h_L3MEt_noicd_E1SHT25 = new TH1F("h_L3MEt_noicd_E1SHT25","L3 MEt noICD E1SHT25", 100,0.0,100.0);


  // Level3 MEt cut histograms:
  h_OfflineMEt_E1SHT25_L3MEt10noICD = new TH1F("h_Offline_MEt_E1SHT25_L3MEt10noICD","Offline MEt passed E1SHT25 and L3MEtnoICD  > 10",100,0.0,100.0);
  h_OfflineMEt_E1SHT25_L3MEt15noICD = new TH1F("h_Offline_MEt_E1SHT25_L3MEt15noICD","Offline MEt passed E1SHT25 and L3MEtnoICD > 15",100,0.0,100.0);
  h_OfflineMEt_E1SHT25_L3MEt20noICD = new TH1F("h_Offline_MEt_E1SHT25_L3MEt20noICD","Offline MEt passed E1SHT25 and L3MEtnoICD > 20",100,0.0,100.0);
  h_OfflineMEt_E1SHT25_L3MEt25noICD = new TH1F("h_Offline_MEt_E1SHT25_L3MEt25noICD","Offline MEt passed E1SHT25 and L3MEtnoICD > 25",100,0.0,100.0);

  // Resolution histograms: OfflineMEt - Level3MEt

  h_diffOfflineL3_MEt = new TH1F("h_Offline-L3_MEt","Offline - L3 MEt", 100,-50.0,50.0);
  h_diffOfflineL3_MEt_noicd= new TH1F("h_Offline-L3_MEt_noicd","Offline - L3 MEt noicd", 100,-50.0,50.0);

  h_diffOfflineL3_MEt_E1SHT25= new TH1F("h_Offline-L3_MEt_E1SHT25","Offline - L3 MEt after E1SHT25", 100,-50.0,50.0);
  h_diffOfflineL3_MEt_noicd_E1SHT25= new TH1F("h_Offline-L3_MEt_noicd_E1SHT25","Offline - L3 MEt noicd after E1SHT25", 100,-50.0,50.0);

}

bool l3Metefficiency::processEvent(cafe::Event& event)
{
  bool passed = true;

    ++total_events;

  // global info (run, event, LBN) ============================================
  const TMBGlobal *global = event.getGlobal();
  int run   = global->runno();
  int evt   = global->evtno();
  float InstLumino   = global->instlum();

  bool run_number = false;

  if (run > 231960) {
  ++total_events_run_cut;

// Begin Get Trigger names ******************************************
  TString trigger;
  trigger+="*";
  Collection<TMBTrigger> triggers=event.getTriggers();
  for (int t = 0; t < triggers.size(); t++) {
  trigger+=triggers[t].getTrgName();
    trigger+="***";
  }
  (*fired_triggers)<<trigger<<endl;   // print output in a file
// End Get Trigger names ********************************************

// Get Level3 MET, Offline MET and Extended/Expert Offline MET ****************
  Collection<TMBL3MEt>   L3MEts   = event.getL3MEts();
  const TMBMet *missingET = event.getMet(); 
  const TMBMetEx *missingET_Ex = event.getMetEx();
// ******************************************************************

// Begin Cuts to get sample "a la WMass Group" ***********************************************************

  Collection<TMBEMCluster> EMs = event.getEMscone();
  Collection<TMBEMCluster>::iterator emcl;
  Collection<TMBEMCluster> emobj;

  Collection<TMBTrack> tracks = event.getTracks();
  Collection<TMBTrack>::iterator trk;
  Collection<TMBTrack> trkobj;

  bool quality = false;
  bool is_highpt = false;
  bool good_track = false;
  float Electronpt, Electron_theta, Electron_phi, Electron_Energy, Electronpx, Electronpy, Electronpz  ;

  // Loop over EM objects
  ////////////////////////////
  //is_highpt only set to true if object passes quality.
  //emobj container now only contains quality, high pT objects.
  //loop over tracks only performed if there is one and only one good high pT object.

  for ( emcl=EMs.begin(); emcl!=EMs.end(); ++emcl){ //loop over TMBEMclusters.
    
    //Check for vertex, if vertex not present, use z==0
    const TMBVertex* vertex = emcl->GetVertex();
    float z_vtx=0;
    if( emcl->GetVertex()!= NULL)
      float z_vtx = vertex->vz();
    //Done checking vertex.    

    //quality cuts
    if(
       abs(emcl->id()) < 12
       && emcl->is_in_eta_fiducial()
       && emcl->iso() < 0.15
       && emcl->emfrac() >0.9
       && fabs(emcl->CalDetectorEta())< 1.1
       && fabs(z_vtx) < 40 
       && emcl->track_match_spatialchi2prob()>0.01
       && emcl->flrS1(3) < 14){ 
      quality = true;  
      
      //If and only if object passes quality, calculate pT.
      if (quality){
	Electronpt = emcl->CalE()*sin(kinem::theta(emcl->CalEta()));
	Electron_theta = kinem::theta(emcl->CalEta());
	Electron_phi = emcl->CalPhi();
	Electron_Energy = Electronpt/sin(Electron_theta);
	Electronpx = Electronpt*cos(emcl->CalPhi());
	Electronpy = Electronpt*sin(emcl->CalPhi());
	Electronpz = Electronpt/tan(Electron_theta);
	if ( Electronpt > 25.0 ){ 
	  is_highpt = true;
	  emobj.push_back(*emcl);//put object in container.
	}//if object is high pT
      }//if object quality
    }//quality criterion check.
  }//Loop over EM clusters closes.
   
  //Note that this will veto events with two good central track matched EMs. There should be some
  //Z events in the sample, this will partially remove them.

  if(emobj.size() == 1 && quality && is_highpt){

    for(trk=tracks.begin(); trk!=tracks.end(); ++trk) {
      float dphi = kinem::delta_phi(trk->Phi(), emobj[0].CalPhi());
      float deta = trk->Eta() - emobj[0].CalEta();
      float dr = sqrt(dphi*dphi+deta*deta);
      if (dr<=0.4) {good_track=true; 
      
	trkobj.push_back(*trk);}
      
    }//Looop over tracks.
    
  }//event quality checker

  //Check again, this time impose track isolation.

  // Here impose CUTS:
  if(emobj.size() == 1 && quality && is_highpt && good_track && trkobj.size()==1){
  //  if(emobj.size() == 1 && quality && is_highpt){
        

    // Begin W-mass GROUP MEt calculation **********************************************************************
    float metx=0., mety=0., met=0., scalarEt=0.;
    float metx_noicd=0., mety_noicd=0., met_noicd=0., scalarEt_noicd=0., TMass=0.;
    
    float SumEx_em=0., SumEy_em=0., SumE_em=0., scalarEt_em=0.;
    float SumEx_mg=0., SumEy_mg=0., SumE_mg=0., scalarEt_mg=0.;
    float SumEx_fh=0., SumEy_fh=0., SumE_fh=0., scalarEt_fh=0.;
    float SumEx_icd=0., SumEy_icd=0., SumE_icd=0., scalarEt_icd=0.;
    
    // Energies measured from EM layers(1-7), massless gap layers(8&10), FH layers(11-14) and ICD
    missingET_Ex->getVETEM(scalarEt_em, SumEx_em, SumEy_em, SumE_em);
    missingET_Ex->getVETMG(scalarEt_mg, SumEx_mg, SumEy_mg, SumE_mg);
    missingET_Ex->getVETFH(scalarEt_fh, SumEx_fh, SumEy_fh, SumE_fh);
    missingET_Ex->getVETICD(scalarEt_icd, SumEx_icd, SumEy_icd, SumE_icd);
    
    //This quantity doesn't have ICD in it.
    // Calculate MEtx, MEty, MEt and Scalar Et
    metx = -(SumEx_em + SumEx_mg + SumEx_fh);
    mety = -(SumEy_em + SumEy_mg + SumEy_fh);
    met = sqrt(metx*metx + mety*mety);
    scalarEt = scalarEt_em + scalarEt_mg + scalarEt_fh;
    
    // Apply EM corrections
    double Ex_EMObjs_raw = 0., Ey_EMObjs_raw = 0., scalarEt_EMObjs_raw = 0.;
    double Ex_EMObjs_corr = 0., Ey_EMObjs_corr = 0., scalarEt_EMObjs_corr = 0.;
    
    //Note from A. A. :  W-Mass group is reconsidering this calculation.  New version would
    //compensate missing ET for only the selected EM object that was chosen above.
    //It would be easy to substitute emobj collection rather than the full one here.

    // numbers used below are taken from em_evt/rcp/GoodEMParticle.rcp
    cafe::Collection<TMBEMCluster> em = event.getEMscone();
    for (Collection<TMBEMCluster>::iterator it = em.begin(); it!= em.end(); ++it) {
      
      if ((abs((*it).id())==10 || abs((*it).id())==11) && ((*it).Pt()>5.)) {
	if((*it).emfrac()>0.9 && (*it).iso()<0.2 && fabs((*it).CalDetectorEta())<2.5) {
	  double rawE = (*it).uncorrE();
	  double corrE = (*it).E();
	  // calculated used positive cells, EM corrections not applied
	  double rawScalarEt = (*it).EMScalarEt(); 
	  
	  double theta = (*it).Theta();
	  double phi = (*it).Phi();
	  Ex_EMObjs_corr += corrE * sin(theta) * cos(phi);
	  Ey_EMObjs_corr += corrE * sin(theta) * sin(phi);
	  
	  Ex_EMObjs_raw += rawE * sin(theta) * cos(phi);
	  Ey_EMObjs_raw += rawE * sin(theta) * sin(phi);
	  
	  // because EM corrections are applied on EM cluster, not individual cells
	  // the corrected scalarEt = rawScalarEt * correction factor
	  scalarEt_EMObjs_raw += rawScalarEt;
	  scalarEt_EMObjs_corr += rawScalarEt * corrE/rawE; 
	}
      }
    } // loop over all EM particles
    
    double metx_EMCorr = 0., mety_EMCorr = 0., met_EMCorr = 0., scalarEt_EMCorr = 0.;
    metx_EMCorr = metx + (Ex_EMObjs_raw - Ex_EMObjs_corr);
    mety_EMCorr = mety + (Ey_EMObjs_raw - Ey_EMObjs_corr);
    met_EMCorr = sqrt(metx_EMCorr*metx_EMCorr + mety_EMCorr*mety_EMCorr);
    scalarEt_EMCorr = scalarEt + (scalarEt_EMObjs_raw - scalarEt_EMObjs_corr);
    
    h_Offline_MEt->Fill(met_EMCorr);
    
    TMass = GetTransverseMass(metx_EMCorr, mety_EMCorr, Electron_phi, Electronpt);
    h_Tmass->Fill(TMass);
    ++events_all_cuts ;

 // End W-mass GROUP MEt calculation ***************************************************************

  // Level3 MET ****************************************************
    
    float MET_diffOfflineL3,  MET_diffOfflineL3_E1SHT25;
    float MET_diffOfflineL3_noicd, MET_diffOfflineL3_noicd_E1SHT25  ;

    float L3met, L3met_noicd ;
  
  for(int i=0;i<L3MEts.size();i++){
    
    L3met =L3MEts[i].MEt();
    L3met_noicd = L3MEts[i].MEt_noICD();
    h_L3MEt->Fill(L3met);
    h_L3MEt_noicd->Fill(L3met_noicd);
  }
  
  // Difference between OfflineMEt and Level3MEt: Resolution.

  MET_diffOfflineL3 = met_EMCorr-L3met;
  h_diffOfflineL3_MEt -> Fill(MET_diffOfflineL3);

  MET_diffOfflineL3_noicd = met_EMCorr-L3met_noicd;
  h_diffOfflineL3_MEt_noicd -> Fill(MET_diffOfflineL3_noicd);


  // AFTER TRIGGER SELECTION: E1_SHT25
  TString selected_trigger;
  for (int j=0; j < triggers.size(); j++) {
    selected_trigger=triggers[j].getTrgName();
    if (selected_trigger=="E1_SHT25") {

      h_L3MEt_E1SHT25 ->Fill(L3met);
      h_L3MEt_noicd_E1SHT25 ->Fill(L3met_noicd);
      h_OfflineMEt_E1SHT25 ->Fill(met_EMCorr);

      MET_diffOfflineL3_E1SHT25 = met_EMCorr-L3met;
      h_diffOfflineL3_MEt_E1SHT25 -> Fill(MET_diffOfflineL3_E1SHT25);

      MET_diffOfflineL3_noicd_E1SHT25 = met_EMCorr-L3met_noicd;
      h_diffOfflineL3_MEt_noicd_E1SHT25-> Fill(MET_diffOfflineL3_noicd_E1SHT25);

  // Number of offline events that passed trigger E1SHT25:
      ++Offline_events_E1SHT25;

  // Instantaneous luminosity after trigger E1SHT25
      h_instLumi_E1SHT25 ->Fill(InstLumino);

      // We can use an offlineMEt cut:
      //      if(met_EMCorr>25){
      //	++total_events_offlinecut ;

    for(int i=0;i<L3MEts.size();i++){
      
  // Events passed the trigger. We look at the OfflineMEt 
  // after the Level3 MEt passed a cut of 10, 15, 20 or 25 GeV:
  // Histograms are filled and we keep track of number of events
  // that passed for efficiency studies:

      if (selected_trigger=="E1_SHT25" && L3MEts[i].MEt_noICD()>10){
	h_OfflineMEt_E1SHT25_L3MEt10noICD-> Fill(met_EMCorr);
	++Offline_events_L3noICDMETcut_10Gev;
	h_instLumi_L3MET10 ->Fill(InstLumino);}
	
      if (selected_trigger=="E1_SHT25" && L3MEts[i].MEt_noICD()>15){
	h_OfflineMEt_E1SHT25_L3MEt15noICD-> Fill(met_EMCorr);
	++Offline_events_L3noICDMETcut_15Gev;
	h_instLumi_L3MET15 ->Fill(InstLumino);}
      
      if (selected_trigger=="E1_SHT25" && L3MEts[i].MEt_noICD()>20){
	h_OfflineMEt_E1SHT25_L3MEt20noICD-> Fill(met_EMCorr);
	++Offline_events_L3noICDMETcut_20Gev;
	h_instLumi_L3MET20 ->Fill(InstLumino);}

      if (selected_trigger=="E1_SHT25" && L3MEts[i].MEt_noICD()>25){
	h_OfflineMEt_E1SHT25_L3MEt25noICD-> Fill(met_EMCorr);
	++Offline_events_L3noICDMETcut_25Gev;
	h_instLumi_L3MET25 ->Fill(InstLumino);}

      if (selected_trigger=="E1_SHT25" && L3MEts[i].MEt()>10){
	++Offline_events_with_L3METcut_10Gev; }

      if (selected_trigger=="E1_SHT25" && L3MEts[i].MEt()>15){
	++Offline_events_L3METcut_15Gev; }

      if (selected_trigger=="E1_SHT25" && L3MEts[i].MEt()>20){
	++Offline_events_L3METcut_20Gev; }

      if (selected_trigger=="E1_SHT25" && L3MEts[i].MEt()>25){
	++Offline_events_L3METcut_25Gev; }
    }
    //	   }
    }
  }
  }//Closes condition on object quality. 
  }
  return passed ;

}

void l3Metefficiency::finish() 
{
  cout <<"Total Events: "<<total_events<<endl; 
  cout <<"Total Events after run cut: "<<total_events_after_run_cut<<endl;
  cout <<"Total Events after ALL cuts: "<<events_after_all_cuts<<endl;
  cout <<"Events after trigger E1SHT25: "<<Offline_events_after_E1_SHT25<<endl; 
  //  cout <<"Events after offline cut of 25GeV:"<<total_events_offline25cut <<endl; 

  cout <<"number Offline events after L3 MEt cut of 10GeV: "<<Offline_events_L3METcut_10Gev<< endl;
  cout <<"number Offline events after L3 MEt cut of 15GeV: "<<Offline_events_L3METcut_15Gev<< endl;
  cout <<"number Offline events after L3 MEt cut of 20GeV: "<<Offline_events_L3METcut_20Gev<< endl;
  cout <<"number Offline events after L3 MEt cut of 25GeV: "<<Offline_events_L3METcut_25Gev<< endl;

  cout <<"number Offline events after L3 noICDMEt cut of 10GeV: "<<Offline_events_L3noICDMETcut_10Gev<< endl;
  cout <<"number Offline events after L3 noICDMEt cut of 15GeV: "<<Offline_events_L3noICDMETcut_15Gev<< endl;
  cout <<"number Offline events after L3 noICDMEt cut of 20GeV: "<<Offline_events_L3noICDMETcut_20Gev<< endl;
  cout <<"number Offline events after L3 noICDMEt cut of 25GeV: "<<Offline_events_L3noICDMETcut_25Gev<< endl;
}

void l3Metefficiency::inputFileOpened(TFile *file)
{

}

void l3Metefficiency::inputFileClosing(TFile *file)
{

}

float l3Metefficiency::GetTransverseMass(Float_t Metx, Float_t Mety, Float_t ElecPhi, Float_t ElecPT)
{
  Float_t TransverseMass = 0;
  Float_t MET = sqrt(Metx*Metx + Mety*Mety);
  Float_t PX = ElecPT * cos(ElecPhi);
  Float_t PY = ElecPT * sin(ElecPhi);
  Float_t DotProd = PX * Metx + PY * Mety;
  Float_t CosDeltaPhi = DotProd/(MET*ElecPT);

  //The following does the transverse mass calculation for W->e nu
  TransverseMass =sqrt(2*ElecPT* MET *(1-CosDeltaPhi));
  return TransverseMass;
}


ClassImp(l3Metefficiency);