SelMods/src/HwwExampleAnalysisMod.cc

 // $Id $

#include "MitPhysics/SelMods/interface/HwwExampleAnalysisMod.h"
#include <TH1D.h>
#include <TH2D.h>
#include <TParameter.h>
#include "MitAna/DataUtil/interface/Debug.h"
#include "MitAna/DataTree/interface/Names.h"
#include "MitPhysics/Init/interface/ModNames.h"
#include "MitAna/DataCont/interface/ObjArray.h"
#include "MitCommon/MathTools/interface/MathUtils.h"
#include "MitAna/DataTree/interface/ParticleCol.h"
#include "TFile.h"
#include "TTree.h"

using namespace mithep;
ClassImp(mithep::HwwExampleAnalysisMod)

//--------------------------------------------------------------------------------------------------
HwwExampleAnalysisMod::HwwExampleAnalysisMod(const char *name, const char *title) : 
  BaseMod(name,title),
  fMuonBranchName(Names::gkMuonBrn),
  fMetName("NoDefaultNameSet"),
  fCleanJetsName("NoDefaultNameSet"),
  fVertexName(string("PrimaryVertexes").c_str()),
  fMuons(0),
  fMet(0),
  fVertices(0)
{
  // Constructor.
}

//--------------------------------------------------------------------------------------------------
void HwwExampleAnalysisMod::Begin()
{
  // Run startup code on the client machine. For this module, we dont do
  // anything here.
}

//--------------------------------------------------------------------------------------------------
void HwwExampleAnalysisMod::SlaveBegin()
{
  // Run startup code on the computer (slave) doing the actual analysis. Here,
  // we typically initialize histograms and other analysis objects and request
  // branches. For this module, we request a branch of the MitTree.

  // Load Branches
  ReqBranch(fMuonBranchName,   fMuons);
  ReqBranch(fVertexName,    fVertices);

  //Create your histograms here

  //*************************************************************************************************
  // Selection Histograms
  //*************************************************************************************************
  AddTH1(fHWWSelection,"hHWWSelection", ";Cut Number;Number of Events", 8, -1.5, 6.5);
  AddTH1(fHWWToEESelection,"hHWWToEESelection", ";Cut Number;Number of Events", 8, -1.5, 6.5);
  AddTH1(fHWWToMuMuSelection,"hHWWToMuMuSelection", ";Cut Number;Number of Events", 8, -1.5, 6.5);
  AddTH1(fHWWToEMuSelection,"hHWWToEMuSelection", ";Cut Number;Number of Events", 8, -1.5, 6.5);

  //***********************************************************************************************
  // Histograms after preselection
  //***********************************************************************************************
  AddTH1(fLeptonEta          ,"hLeptonEta",";LeptonEta;Number of Events",100,-5.,5.0);
  AddTH1(fLeptonPtMax        ,"hLeptonPtMax",";Lepton P_t Max;Number of Events",150,0.,150.);
  AddTH1(fLeptonPtMin        ,"hLeptonPtMin",";Lepton P_t Min;Number of Events",150,0.,150.);
  AddTH1(fMetPtHist          ,"hMetPtHist",";Met;Number of Events",150,0.,300.);  
  AddTH1(fMetPhiHist         ,"hMetPhiHist",";#phi;Number of Events",28,-3.5,3.5);
  AddTH1(fUncorrMetPtHist    ,"hUncorrMetPtHist",";Met;Number of Events",150,0.,300.);  
  AddTH1(fUncorrMetPhiHist   ,"hUncorrMetPhiHist",";#phi;Number of Events",28,-3.5,3.5);
  AddTH1(fDeltaPhiLeptons    ,"hDeltaPhiLeptons",";#Delta#phi_{ll};Number of Events",90,0,180);
  AddTH1(fDeltaEtaLeptons    ,"hDeltaEtaLeptons",";#Delta#eta_{ll};Number of Events",100,-50.,5.0);
  AddTH1(fDileptonMass       ,"hDileptonMass",";Mass_{ll};Number of Events",150,0.,300.);
 
  //***********************************************************************************************
  // N-1 Histograms
  //***********************************************************************************************
  //All events
  AddTH1(fLeptonPtMax_NMinusOne            ,"hLeptonPtMax_NMinusOne",
                                            ";Lepton P_t Max;Number of Events",150,0.,150.);
  AddTH1(fLeptonPtMin_NMinusOne            ,"hLeptonPtMin_NMinusOne",
                                            ";Lepton P_t Min;Number of Events",150,0.,150.);
  AddTH1(fMetPtHist_NMinusOne              ,"hMetPtHist_NMinusOne",
                                            ";Met;Number of Events",150,0.,300.);  
  AddTH1(fMetPhiHist_NMinusOne             ,"hMetPhiHist_NMinusOne",
                                            ";#phi;Number of Events",28,-3.5,3.5);
  AddTH1(fMETdeltaPhilEtHist_NMinusOne     ,"hMETdeltaPhilEtHist_NMinusOne",
                                            ";METdeltaPhilEtHist;Number of Events",150,0.,300.);
  AddTH1(fNCentralJets_NMinusOne           ,"hNCentralJets_NMinusOne",
                                            ";Number of Central Jets;Number of Events",6,-0.5,5.5);
  AddTH1(fNSoftMuonsHist_NMinusOne        ,"hNSoftMuonsHist_NMinusOne",
                                            ";Number of Dirty Muons; Number of Events",6,-0.5,5.5);
  AddTH1(fDeltaPhiLeptons_NMinusOne        ,"hDeltaPhiLeptons_NMinusOne",
                                            ";#Delta#phi_{ll};Number of Events",90,0,180);
  AddTH1(fDeltaEtaLeptons_NMinusOne        ,"hDeltaEtaLeptons_NMinusOne",
                                            ";#Delta#eta_{ll};Number of Events",100,-5.0,5.0);
  AddTH1(fDileptonMass_NMinusOne           ,"hDileptonMass_NMinusOne",
                                            ";Mass_{ll};Number of Events",150,0.,300.);
  AddTH1(fMinDeltaPhiLeptonMet_NMinusOne   ,"hMinDeltaPhiLeptonMet_NMinusOne", 
                                            ";Min #Delta#phi_{l,Met};Number of Events",90,0.,180);

  //***********************************************************************************************
  // After all cuts Histograms
  //***********************************************************************************************
  AddTH1(fMinDeltaPhiLeptonMet_afterCuts    ,"hMinDeltaPhiLeptonMet_afterCuts", 
                                             ";Min #Delta#phi_{l,Met};Number of Events",90,0.,180);
  AddTH1(fMtLepton1_afterCuts               ,"hMtLepton1_afterCuts",
                                             ";M_t (Lepton1,Met);Number of Events",100,0.,200.);
  AddTH1(fMtLepton2_afterCuts               ,"hMtLepton2_afterCuts",
                                             ";M_t (Lepton2,Met);Number of Events",100,0.,200.);
  AddTH1(fMtHiggs_afterCuts                 ,"hMtHiggs_afterCuts",
                                             ";M_t (l1+l2+Met);Number of Events",150,0.,300.);
  AddTH1(fLeptonPtPlusMet_afterCuts         ,"hLeptonPtPlusMet_afterCuts",
                                             ";LeptonPtPlusMet;Number of Events",150,0., 300.);

  //*********************
  //Added Test Histograms
  //********************
}

//--------------------------------------------------------------------------------------------------
void HwwExampleAnalysisMod::Process()
{
  // Process entries of the tree. For this module, we just load the branches and  
  LoadBranch(fMuonBranchName);
  LoadBranch(fVertexName);
  
  //Obtain all the good objects from the event cleaning module
  ObjArray<Muon> *CleanMuons = dynamic_cast<ObjArray<Muon>* >(FindObjThisEvt(ModNames::gkCleanMuonsName));
  ParticleOArr *CleanLeptons = dynamic_cast<mithep::ParticleOArr*>
    (FindObjThisEvt(ModNames::gkMergedLeptonsName));
  ObjArray<Jet> *CleanJets = dynamic_cast<ObjArray<Jet>* >
    (FindObjThisEvt(fCleanJetsName.Data()));
  TParameter<Double_t> *NNLOWeight = GetObjThisEvt<TParameter<Double_t> >("NNLOWeight");

  MetCol *met = dynamic_cast<ObjArray<Met>* >(FindObjThisEvt(fMetName));
  const Met *caloMet = 0;
  if (met) {
    caloMet = met->At(0);
  } else {
    cout << "Error: Met Collection " << fMetName << " could not be loaded.\n";
    return;
  }

  //***********************************************************************************************
  //Kinematic PreSelection
  //***********************************************************************************************
  // At least two leptons in the event
  if (CleanLeptons->GetEntries() < 2) return;
  // Pt1 > 20 && Pt2 > 10
  if(CleanLeptons->At(0)->Pt() <= 20 || CleanLeptons->At(1)->Pt() <= 20) return;
  // opposite charge leptons
  if(CleanLeptons->At(0)->Charge() * CleanLeptons->At(1)->Charge() > 0) return;
    
  CompositeParticle *dilepton = new CompositeParticle();
  dilepton->AddDaughter(CleanLeptons->At(0));
  dilepton->AddDaughter(CleanLeptons->At(1));
   
  //***********************************************************************************************
  //Get Dirty Muons: Non-isolated Muons (exclude the clean muons)
  //***********************************************************************************************
  ObjArray<Muon> *SoftMuons = new ObjArray<Muon>;
  for (UInt_t i=0; i<fMuons->GetEntries(); ++i) {
    const Muon *mu = fMuons->At(i);
    if(!MuonTools::PassSoftMuonCut(mu, fVertices)) continue;
    
    bool isCleanMuon = kFALSE;
    for (UInt_t j=0; j<CleanMuons->GetEntries(); j++) {
      if(fMuons->At(i) == CleanMuons->At(j) &&
         CleanMuons->At(j)->Pt() > 10) isCleanMuon = kTRUE;
    }
    if(isCleanMuon == kFALSE) SoftMuons->Add(mu);
  }

  //***********************************************************************************************
  //Define Event Variables
  //***********************************************************************************************
  //delta phi between the 2 leptons in degrees
  double deltaPhiLeptons = MathUtils::DeltaPhi(CleanLeptons->At(0)->Phi(), 
                                               CleanLeptons->At(1)->Phi())* 180.0 / TMath::Pi();

  double deltaEtaLeptons = abs(CleanLeptons->At(0)->Eta() - CleanLeptons->At(1)->Eta()) * 180.0 / TMath::Pi();

  double deltaPhiDileptonMet = MathUtils::DeltaPhi(caloMet->Phi(), 
                                                   dilepton->Phi())*180.0 / TMath::Pi();

  double mtHiggs = TMath::Sqrt(2.0*dilepton->Pt() * caloMet->Pt()*
                               (1.0 - cos(deltaPhiDileptonMet * TMath::Pi() / 180.0)));

  //angle between MET and closest lepton
  double deltaPhiMetLepton[2] = {MathUtils::DeltaPhi(caloMet->Phi(), CleanLeptons->At(0)->Phi()),
                                 MathUtils::DeltaPhi(caloMet->Phi(), CleanLeptons->At(1)->Phi())};
  
  double mTW[2] = {TMath::Sqrt(2.0*CleanLeptons->At(0)->Pt()*caloMet->Pt()*
                               (1.0 - cos(deltaPhiMetLepton[0]))),
                   TMath::Sqrt(2.0*CleanLeptons->At(1)->Pt()*caloMet->Pt()*
                               (1.0 - cos(deltaPhiMetLepton[1])))};

  double minDeltaPhiMetLepton = (deltaPhiMetLepton[0] < deltaPhiMetLepton[1])?
    deltaPhiMetLepton[0]:deltaPhiMetLepton[1];

  double METdeltaPhilEt = caloMet->Pt();
  if(minDeltaPhiMetLepton < TMath::Pi()/2.)
      METdeltaPhilEt = METdeltaPhilEt * sin(minDeltaPhiMetLepton);

  //count the number of central Jets for vetoing
  int nCentralJets = 0;
  for(UInt_t i=0; i<CleanJets->GetEntries(); i++){
    if(TMath::Abs(CleanJets->At(i)->Eta()) < 5.0 &&
       CleanJets->At(i)->Pt() > 25.0){
      nCentralJets++;
    }
  }

  //Lepton Type
  int finalstateType = -1;
  if (CleanLeptons->At(0)->ObjType() == kMuon && CleanLeptons->At(1)->ObjType() == kMuon ){ // mumu
    finalstateType = 10;
  } else if(CleanLeptons->At(0)->ObjType() == kElectron && CleanLeptons->At(1)->ObjType() == kElectron ){ // ee
    finalstateType = 11;
  } else if((CleanLeptons->At(0)->ObjType() == kElectron && CleanLeptons->At(1)->ObjType() == kMuon) ||
            (CleanLeptons->At(1)->ObjType() == kElectron && CleanLeptons->At(0)->ObjType() == kMuon)) {
    finalstateType = 12;
  } else {
    cerr << "Error: finalstate lepton type not supported\n";
  }
                        
  //*********************************************************************************************
  //Define Cuts
  //*********************************************************************************************
  const int nCuts = 7;
  bool passCut[nCuts] = {false, false, false, false, false, false, false};
  
  if(CleanLeptons->At(0)->Pt() >  20.0 &&
     CleanLeptons->At(1)->Pt() >= 20.0) passCut[0] = true;
  
  if(caloMet->Pt()    > 20.0)           passCut[1] = true;
  
  if(dilepton->Mass() > 12.0)           passCut[2] = true;
  
  if(nCentralJets     < 1)              passCut[5] = true;

  if(CleanLeptons->GetEntries() == 2 &&
     SoftMuons->GetEntries() == 0)      passCut[6] = true;

  if (finalstateType == 10 || finalstateType == 11){ // mumu/ee
    if(fabs(dilepton->Mass()-91.1876)   > 15.0)   passCut[3] = true;
    if(METdeltaPhilEt > 35) passCut[4] = true;
  }
  else if(finalstateType == 12) { // emu
    passCut[3] = true;
    if(METdeltaPhilEt > 20) passCut[4] = true;
  }
  
  //*********************************************************************************************
  //Make Selection Histograms. Number of events passing each level of cut
  //*********************************************************************************************  
  bool passAllCuts = true;
  for(int c=0; c<nCuts; c++) passAllCuts = passAllCuts & passCut[c];
    
  //Cut Selection Histograms
  fHWWSelection->Fill(-1,NNLOWeight->GetVal());
  if (finalstateType == 10 )
    fHWWToMuMuSelection->Fill(-1,NNLOWeight->GetVal());
  else if(finalstateType == 11 )
    fHWWToEESelection->Fill(-1,NNLOWeight->GetVal());
  else if(finalstateType == 12 )
    fHWWToEMuSelection->Fill(-1,NNLOWeight->GetVal());

  for (int k=0;k<nCuts;k++) {
    bool pass = true;
    bool passPreviousCut = true;
    for (int p=0;p<=k;p++) {
      pass = (pass && passCut[p]);
      if (p<k)
        passPreviousCut = (passPreviousCut&& passCut[p]);
    }
    
    if (pass) {
      fHWWSelection->Fill(k,NNLOWeight->GetVal());
      if (finalstateType == 10 )
        fHWWToMuMuSelection->Fill(k,NNLOWeight->GetVal());
      else if(finalstateType == 11)
        fHWWToEESelection->Fill(k,NNLOWeight->GetVal());
      else if(finalstateType == 12)
        fHWWToEMuSelection->Fill(k,NNLOWeight->GetVal());
    }
  }
  
  //*****************************************************************************************
  //Make Preselection Histograms  
  //*****************************************************************************************
  fLeptonEta->Fill(CleanLeptons->At(0)->Eta(),NNLOWeight->GetVal()); 
  fLeptonEta->Fill(CleanLeptons->At(1)->Eta(),NNLOWeight->GetVal());
  fLeptonPtMax->Fill(CleanLeptons->At(0)->Pt(),NNLOWeight->GetVal());
  fLeptonPtMin->Fill(CleanLeptons->At(1)->Pt(),NNLOWeight->GetVal());
  fMetPtHist->Fill(caloMet->Pt(),NNLOWeight->GetVal());                             
  fMetPhiHist->Fill(caloMet->Phi(),NNLOWeight->GetVal());                            
  fDeltaPhiLeptons->Fill(deltaPhiLeptons,NNLOWeight->GetVal());
  fDeltaEtaLeptons->Fill(deltaEtaLeptons,NNLOWeight->GetVal());
  fDileptonMass->Fill(dilepton->Mass(),NNLOWeight->GetVal());    

  //*********************************************************************************************
  // N-1 Histograms
  //*********************************************************************************************
  bool pass;;
  
  //N Jet Veto  
  pass = true;
  for (int k=0;k<nCuts;k++) {
    if (k != 5) {
      pass = (pass && passCut[k]);      
    }
  }
  if (pass) {
    fNCentralJets_NMinusOne->Fill(nCentralJets,NNLOWeight->GetVal());
  }     
  
  // Final Met Cut
  pass = true;
  for (int k=0;k<nCuts;k++) {
    if (k != 4) {
      pass = (pass && passCut[k]);      
    }
  }
  if (pass) {
    fMetPtHist_NMinusOne->Fill(caloMet->Pt(),NNLOWeight->GetVal());  
  }

  // dilepton mass
  pass = true;
  for (int k=0;k<nCuts;k++) {
    if (k != 2 && k !=  3)
      pass = (pass && passCut[k]);    
  }
  if (pass) {
    fDileptonMass_NMinusOne->Fill(dilepton->Mass(),NNLOWeight->GetVal());
  }
  
  // Lepton Pt Max, Lepton Pt Min, DeltaPhiLeptons
  pass = true;
  for (int k=0;k<nCuts;k++) {
    pass = (pass && passCut[k]);      
  }
  if (pass) {
    fLeptonPtMax_NMinusOne->Fill(CleanLeptons->At(0)->Pt(),NNLOWeight->GetVal());
    fLeptonPtMin_NMinusOne->Fill(CleanLeptons->At(1)->Pt(),NNLOWeight->GetVal());
    fDeltaPhiLeptons_NMinusOne->Fill(deltaPhiLeptons,NNLOWeight->GetVal()); 
  }
  
  // NSoftMuons
  pass = true;
  for (int k=0;k<nCuts;k++) {
    if (k != 6)
      pass = (pass && passCut[k]);    
  }
  if (pass) {
    fNSoftMuonsHist_NMinusOne->Fill(SoftMuons->GetEntries(),NNLOWeight->GetVal());
  }

  //*********************************************************************************************
  //Plots after all Cuts
  //*********************************************************************************************
  if (passAllCuts) {
    fMinDeltaPhiLeptonMet_afterCuts->Fill(minDeltaPhiMetLepton,NNLOWeight->GetVal());
    fMtLepton1_afterCuts->Fill(mTW[0],NNLOWeight->GetVal());
    fMtLepton2_afterCuts->Fill(mTW[1],NNLOWeight->GetVal());
    fMtHiggs_afterCuts->Fill(mtHiggs,NNLOWeight->GetVal());
    fLeptonPtPlusMet_afterCuts->Fill(CleanLeptons->At(0)->Pt()+CleanLeptons->At(1)->Pt()+caloMet->Pt(),NNLOWeight->GetVal());
  }
  
  delete dilepton;
  delete SoftMuons;
  return;
}

//--------------------------------------------------------------------------------------------------
void HwwExampleAnalysisMod::SlaveTerminate()
{
  
  // Run finishing code on the computer (slave) that did the analysis. For this
  // module, we dont do anything here.

} 
//--------------------------------------------------------------------------------------------------
void HwwExampleAnalysisMod::Terminate()
{
  // Run finishing code on the client computer. For this module, we dont do
  // anything here.

}
Revision:	1.2
Committed:	Tue Oct 5 06:37:11 2010 UTC (14 years, 7 months ago) by ceballos
Content type:	text/plain
Branch:	MAIN
Changes since 1.1:	+30 -29 lines
Log Message:	small update
#	User	Rev	Content
1	ceballos	1.1	// $Id $
2
3			#include "MitPhysics/SelMods/interface/HwwExampleAnalysisMod.h"
4			#include <TH1D.h>
5			#include <TH2D.h>
6			#include <TParameter.h>
7			#include "MitAna/DataUtil/interface/Debug.h"
8			#include "MitAna/DataTree/interface/Names.h"
9			#include "MitPhysics/Init/interface/ModNames.h"
10			#include "MitAna/DataCont/interface/ObjArray.h"
11			#include "MitCommon/MathTools/interface/MathUtils.h"
12			#include "MitAna/DataTree/interface/ParticleCol.h"
13			#include "TFile.h"
14			#include "TTree.h"
15
16			using namespace mithep;
17			ClassImp(mithep::HwwExampleAnalysisMod)
18
19			//--------------------------------------------------------------------------------------------------
20			HwwExampleAnalysisMod::HwwExampleAnalysisMod(const char name, const char title) :
21			BaseMod(name,title),
22			fMuonBranchName(Names::gkMuonBrn),
23			fMetName("NoDefaultNameSet"),
24			fCleanJetsName("NoDefaultNameSet"),
25			fVertexName(string("PrimaryVertexes").c_str()),
26			fMuons(0),
27			fMet(0),
28			fVertices(0)
29			{
30			// Constructor.
31			}
32
33			//--------------------------------------------------------------------------------------------------
34			void HwwExampleAnalysisMod::Begin()
35			{
36			// Run startup code on the client machine. For this module, we dont do
37			// anything here.
38			}
39
40			//--------------------------------------------------------------------------------------------------
41			void HwwExampleAnalysisMod::SlaveBegin()
42			{
43			// Run startup code on the computer (slave) doing the actual analysis. Here,
44			// we typically initialize histograms and other analysis objects and request
45			// branches. For this module, we request a branch of the MitTree.
46
47			// Load Branches
48			ReqBranch(fMuonBranchName, fMuons);
49			ReqBranch(fVertexName, fVertices);
50
51			//Create your histograms here
52
53			//*************************************************************************************************
54			// Selection Histograms
55			//*************************************************************************************************
56			AddTH1(fHWWSelection,"hHWWSelection", ";Cut Number;Number of Events", 8, -1.5, 6.5);
57			AddTH1(fHWWToEESelection,"hHWWToEESelection", ";Cut Number;Number of Events", 8, -1.5, 6.5);
58			AddTH1(fHWWToMuMuSelection,"hHWWToMuMuSelection", ";Cut Number;Number of Events", 8, -1.5, 6.5);
59			AddTH1(fHWWToEMuSelection,"hHWWToEMuSelection", ";Cut Number;Number of Events", 8, -1.5, 6.5);
60
61			//***********************************************************************************************
62			// Histograms after preselection
63			//***********************************************************************************************
64			AddTH1(fLeptonEta ,"hLeptonEta",";LeptonEta;Number of Events",100,-5.,5.0);
65			AddTH1(fLeptonPtMax ,"hLeptonPtMax",";Lepton P_t Max;Number of Events",150,0.,150.);
66			AddTH1(fLeptonPtMin ,"hLeptonPtMin",";Lepton P_t Min;Number of Events",150,0.,150.);
67			AddTH1(fMetPtHist ,"hMetPtHist",";Met;Number of Events",150,0.,300.);
68			AddTH1(fMetPhiHist ,"hMetPhiHist",";#phi;Number of Events",28,-3.5,3.5);
69			AddTH1(fUncorrMetPtHist ,"hUncorrMetPtHist",";Met;Number of Events",150,0.,300.);
70			AddTH1(fUncorrMetPhiHist ,"hUncorrMetPhiHist",";#phi;Number of Events",28,-3.5,3.5);
71			AddTH1(fDeltaPhiLeptons ,"hDeltaPhiLeptons",";#Delta#phi_{ll};Number of Events",90,0,180);
72			AddTH1(fDeltaEtaLeptons ,"hDeltaEtaLeptons",";#Delta#eta_{ll};Number of Events",100,-50.,5.0);
73			AddTH1(fDileptonMass ,"hDileptonMass",";Mass_{ll};Number of Events",150,0.,300.);
74
75			//***********************************************************************************************
76			// N-1 Histograms
77			//***********************************************************************************************
78			//All events
79			AddTH1(fLeptonPtMax_NMinusOne ,"hLeptonPtMax_NMinusOne",
80			";Lepton P_t Max;Number of Events",150,0.,150.);
81			AddTH1(fLeptonPtMin_NMinusOne ,"hLeptonPtMin_NMinusOne",
82			";Lepton P_t Min;Number of Events",150,0.,150.);
83			AddTH1(fMetPtHist_NMinusOne ,"hMetPtHist_NMinusOne",
84			";Met;Number of Events",150,0.,300.);
85			AddTH1(fMetPhiHist_NMinusOne ,"hMetPhiHist_NMinusOne",
86			";#phi;Number of Events",28,-3.5,3.5);
87			AddTH1(fMETdeltaPhilEtHist_NMinusOne ,"hMETdeltaPhilEtHist_NMinusOne",
88			";METdeltaPhilEtHist;Number of Events",150,0.,300.);
89			AddTH1(fNCentralJets_NMinusOne ,"hNCentralJets_NMinusOne",
90			";Number of Central Jets;Number of Events",6,-0.5,5.5);
91			AddTH1(fNSoftMuonsHist_NMinusOne ,"hNSoftMuonsHist_NMinusOne",
92			";Number of Dirty Muons; Number of Events",6,-0.5,5.5);
93			AddTH1(fDeltaPhiLeptons_NMinusOne ,"hDeltaPhiLeptons_NMinusOne",
94			";#Delta#phi_{ll};Number of Events",90,0,180);
95			AddTH1(fDeltaEtaLeptons_NMinusOne ,"hDeltaEtaLeptons_NMinusOne",
96			";#Delta#eta_{ll};Number of Events",100,-5.0,5.0);
97			AddTH1(fDileptonMass_NMinusOne ,"hDileptonMass_NMinusOne",
98			";Mass_{ll};Number of Events",150,0.,300.);
99			AddTH1(fMinDeltaPhiLeptonMet_NMinusOne ,"hMinDeltaPhiLeptonMet_NMinusOne",
100			";Min #Delta#phi_{l,Met};Number of Events",90,0.,180);
101
102			//***********************************************************************************************
103			// After all cuts Histograms
104			//***********************************************************************************************
105			AddTH1(fMinDeltaPhiLeptonMet_afterCuts ,"hMinDeltaPhiLeptonMet_afterCuts",
106			";Min #Delta#phi_{l,Met};Number of Events",90,0.,180);
107			AddTH1(fMtLepton1_afterCuts ,"hMtLepton1_afterCuts",
108			";M_t (Lepton1,Met);Number of Events",100,0.,200.);
109			AddTH1(fMtLepton2_afterCuts ,"hMtLepton2_afterCuts",
110			";M_t (Lepton2,Met);Number of Events",100,0.,200.);
111			AddTH1(fMtHiggs_afterCuts ,"hMtHiggs_afterCuts",
112			";M_t (l1+l2+Met);Number of Events",150,0.,300.);
113			AddTH1(fLeptonPtPlusMet_afterCuts ,"hLeptonPtPlusMet_afterCuts",
114			";LeptonPtPlusMet;Number of Events",150,0., 300.);
115
116			//*********************
117			//Added Test Histograms
118			//********************
119			}
120
121			//--------------------------------------------------------------------------------------------------
122			void HwwExampleAnalysisMod::Process()
123			{
124			// Process entries of the tree. For this module, we just load the branches and
125			LoadBranch(fMuonBranchName);
126			LoadBranch(fVertexName);
127
128			//Obtain all the good objects from the event cleaning module
129			ObjArray<Muon> CleanMuons = dynamic_cast<ObjArray<Muon> >(FindObjThisEvt(ModNames::gkCleanMuonsName));
130			ParticleOArr CleanLeptons = dynamic_cast<mithep::ParticleOArr>
131			(FindObjThisEvt(ModNames::gkMergedLeptonsName));
132			ObjArray<Jet> CleanJets = dynamic_cast<ObjArray<Jet> >
133			(FindObjThisEvt(fCleanJetsName.Data()));
134	ceballos	1.2	TParameter<Double_t> *NNLOWeight = GetObjThisEvt<TParameter<Double_t> >("NNLOWeight");
135	ceballos	1.1
136			MetCol met = dynamic_cast<ObjArray<Met> >(FindObjThisEvt(fMetName));
137			const Met *caloMet = 0;
138			if (met) {
139			caloMet = met->At(0);
140			} else {
141			cout << "Error: Met Collection " << fMetName << " could not be loaded.\n";
142			return;
143			}
144
145			//***********************************************************************************************
146			//Kinematic PreSelection
147			//***********************************************************************************************
148			// At least two leptons in the event
149			if (CleanLeptons->GetEntries() < 2) return;
150			// Pt1 > 20 && Pt2 > 10
151			if(CleanLeptons->At(0)->Pt() <= 20 \|\| CleanLeptons->At(1)->Pt() <= 20) return;
152			// opposite charge leptons
153			if(CleanLeptons->At(0)->Charge() * CleanLeptons->At(1)->Charge() > 0) return;
154
155			CompositeParticle *dilepton = new CompositeParticle();
156			dilepton->AddDaughter(CleanLeptons->At(0));
157			dilepton->AddDaughter(CleanLeptons->At(1));
158
159			//***********************************************************************************************
160			//Get Dirty Muons: Non-isolated Muons (exclude the clean muons)
161			//***********************************************************************************************
162			ObjArray<Muon> *SoftMuons = new ObjArray<Muon>;
163			for (UInt_t i=0; i<fMuons->GetEntries(); ++i) {
164			const Muon *mu = fMuons->At(i);
165			if(!MuonTools::PassSoftMuonCut(mu, fVertices)) continue;
166
167			bool isCleanMuon = kFALSE;
168			for (UInt_t j=0; j<CleanMuons->GetEntries(); j++) {
169			if(fMuons->At(i) == CleanMuons->At(j) &&
170			CleanMuons->At(j)->Pt() > 10) isCleanMuon = kTRUE;
171			}
172			if(isCleanMuon == kFALSE) SoftMuons->Add(mu);
173			}
174
175			//***********************************************************************************************
176			//Define Event Variables
177			//***********************************************************************************************
178			//delta phi between the 2 leptons in degrees
179			double deltaPhiLeptons = MathUtils::DeltaPhi(CleanLeptons->At(0)->Phi(),
180			CleanLeptons->At(1)->Phi())* 180.0 / TMath::Pi();
181
182			double deltaEtaLeptons = abs(CleanLeptons->At(0)->Eta() - CleanLeptons->At(1)->Eta()) * 180.0 / TMath::Pi();
183
184			double deltaPhiDileptonMet = MathUtils::DeltaPhi(caloMet->Phi(),
185			dilepton->Phi())*180.0 / TMath::Pi();
186
187			double mtHiggs = TMath::Sqrt(2.0dilepton->Pt() caloMet->Pt()*
188			(1.0 - cos(deltaPhiDileptonMet * TMath::Pi() / 180.0)));
189
190			//angle between MET and closest lepton
191			double deltaPhiMetLepton[2] = {MathUtils::DeltaPhi(caloMet->Phi(), CleanLeptons->At(0)->Phi()),
192			MathUtils::DeltaPhi(caloMet->Phi(), CleanLeptons->At(1)->Phi())};
193
194			double mTW[2] = {TMath::Sqrt(2.0CleanLeptons->At(0)->Pt()caloMet->Pt()*
195			(1.0 - cos(deltaPhiMetLepton[0]))),
196			TMath::Sqrt(2.0CleanLeptons->At(1)->Pt()caloMet->Pt()*
197			(1.0 - cos(deltaPhiMetLepton[1])))};
198
199			double minDeltaPhiMetLepton = (deltaPhiMetLepton[0] < deltaPhiMetLepton[1])?
200			deltaPhiMetLepton[0]:deltaPhiMetLepton[1];
201
202			double METdeltaPhilEt = caloMet->Pt();
203			if(minDeltaPhiMetLepton < TMath::Pi()/2.)
204			METdeltaPhilEt = METdeltaPhilEt * sin(minDeltaPhiMetLepton);
205
206			//count the number of central Jets for vetoing
207			int nCentralJets = 0;
208			for(UInt_t i=0; i<CleanJets->GetEntries(); i++){
209			if(TMath::Abs(CleanJets->At(i)->Eta()) < 5.0 &&
210			CleanJets->At(i)->Pt() > 25.0){
211			nCentralJets++;
212			}
213			}
214
215			//Lepton Type
216			int finalstateType = -1;
217			if (CleanLeptons->At(0)->ObjType() == kMuon && CleanLeptons->At(1)->ObjType() == kMuon ){ // mumu
218			finalstateType = 10;
219			} else if(CleanLeptons->At(0)->ObjType() == kElectron && CleanLeptons->At(1)->ObjType() == kElectron ){ // ee
220			finalstateType = 11;
221			} else if((CleanLeptons->At(0)->ObjType() == kElectron && CleanLeptons->At(1)->ObjType() == kMuon) \|\|
222			(CleanLeptons->At(1)->ObjType() == kElectron && CleanLeptons->At(0)->ObjType() == kMuon)) {
223			finalstateType = 12;
224			} else {
225			cerr << "Error: finalstate lepton type not supported\n";
226			}
227
228			//*********************************************************************************************
229			//Define Cuts
230			//*********************************************************************************************
231			const int nCuts = 7;
232			bool passCut[nCuts] = {false, false, false, false, false, false, false};
233
234			if(CleanLeptons->At(0)->Pt() > 20.0 &&
235			CleanLeptons->At(1)->Pt() >= 20.0) passCut[0] = true;
236
237			if(caloMet->Pt() > 20.0) passCut[1] = true;
238
239			if(dilepton->Mass() > 12.0) passCut[2] = true;
240
241			if(nCentralJets < 1) passCut[5] = true;
242
243			if(CleanLeptons->GetEntries() == 2 &&
244			SoftMuons->GetEntries() == 0) passCut[6] = true;
245
246			if (finalstateType == 10 \|\| finalstateType == 11){ // mumu/ee
247			if(fabs(dilepton->Mass()-91.1876) > 15.0) passCut[3] = true;
248			if(METdeltaPhilEt > 35) passCut[4] = true;
249			}
250			else if(finalstateType == 12) { // emu
251			passCut[3] = true;
252			if(METdeltaPhilEt > 20) passCut[4] = true;
253			}
254
255			//*********************************************************************************************
256			//Make Selection Histograms. Number of events passing each level of cut
257			//*********************************************************************************************
258			bool passAllCuts = true;
259			for(int c=0; c<nCuts; c++) passAllCuts = passAllCuts & passCut[c];
260
261			//Cut Selection Histograms
262	ceballos	1.2	fHWWSelection->Fill(-1,NNLOWeight->GetVal());
263	ceballos	1.1	if (finalstateType == 10 )
264	ceballos	1.2	fHWWToMuMuSelection->Fill(-1,NNLOWeight->GetVal());
265	ceballos	1.1	else if(finalstateType == 11 )
266	ceballos	1.2	fHWWToEESelection->Fill(-1,NNLOWeight->GetVal());
267	ceballos	1.1	else if(finalstateType == 12 )
268	ceballos	1.2	fHWWToEMuSelection->Fill(-1,NNLOWeight->GetVal());
269	ceballos	1.1
270			for (int k=0;k<nCuts;k++) {
271			bool pass = true;
272			bool passPreviousCut = true;
273			for (int p=0;p<=k;p++) {
274			pass = (pass && passCut[p]);
275			if (p<k)
276			passPreviousCut = (passPreviousCut&& passCut[p]);
277			}
278
279			if (pass) {
280	ceballos	1.2	fHWWSelection->Fill(k,NNLOWeight->GetVal());
281	ceballos	1.1	if (finalstateType == 10 )
282	ceballos	1.2	fHWWToMuMuSelection->Fill(k,NNLOWeight->GetVal());
283	ceballos	1.1	else if(finalstateType == 11)
284	ceballos	1.2	fHWWToEESelection->Fill(k,NNLOWeight->GetVal());
285	ceballos	1.1	else if(finalstateType == 12)
286	ceballos	1.2	fHWWToEMuSelection->Fill(k,NNLOWeight->GetVal());
287	ceballos	1.1	}
288			}
289
290			//*****************************************************************************************
291			//Make Preselection Histograms
292			//*****************************************************************************************
293	ceballos	1.2	fLeptonEta->Fill(CleanLeptons->At(0)->Eta(),NNLOWeight->GetVal());
294			fLeptonEta->Fill(CleanLeptons->At(1)->Eta(),NNLOWeight->GetVal());
295			fLeptonPtMax->Fill(CleanLeptons->At(0)->Pt(),NNLOWeight->GetVal());
296			fLeptonPtMin->Fill(CleanLeptons->At(1)->Pt(),NNLOWeight->GetVal());
297			fMetPtHist->Fill(caloMet->Pt(),NNLOWeight->GetVal());
298			fMetPhiHist->Fill(caloMet->Phi(),NNLOWeight->GetVal());
299			fDeltaPhiLeptons->Fill(deltaPhiLeptons,NNLOWeight->GetVal());
300			fDeltaEtaLeptons->Fill(deltaEtaLeptons,NNLOWeight->GetVal());
301			fDileptonMass->Fill(dilepton->Mass(),NNLOWeight->GetVal());
302	ceballos	1.1
303			//*********************************************************************************************
304			// N-1 Histograms
305			//*********************************************************************************************
306			bool pass;;
307
308			//N Jet Veto
309			pass = true;
310			for (int k=0;k<nCuts;k++) {
311			if (k != 5) {
312			pass = (pass && passCut[k]);
313			}
314			}
315			if (pass) {
316	ceballos	1.2	fNCentralJets_NMinusOne->Fill(nCentralJets,NNLOWeight->GetVal());
317	ceballos	1.1	}
318
319			// Final Met Cut
320			pass = true;
321			for (int k=0;k<nCuts;k++) {
322			if (k != 4) {
323			pass = (pass && passCut[k]);
324			}
325			}
326			if (pass) {
327	ceballos	1.2	fMetPtHist_NMinusOne->Fill(caloMet->Pt(),NNLOWeight->GetVal());
328	ceballos	1.1	}
329
330			// dilepton mass
331			pass = true;
332			for (int k=0;k<nCuts;k++) {
333			if (k != 2 && k != 3)
334			pass = (pass && passCut[k]);
335			}
336			if (pass) {
337	ceballos	1.2	fDileptonMass_NMinusOne->Fill(dilepton->Mass(),NNLOWeight->GetVal());
338	ceballos	1.1	}
339
340			// Lepton Pt Max, Lepton Pt Min, DeltaPhiLeptons
341			pass = true;
342			for (int k=0;k<nCuts;k++) {
343			pass = (pass && passCut[k]);
344			}
345			if (pass) {
346	ceballos	1.2	fLeptonPtMax_NMinusOne->Fill(CleanLeptons->At(0)->Pt(),NNLOWeight->GetVal());
347			fLeptonPtMin_NMinusOne->Fill(CleanLeptons->At(1)->Pt(),NNLOWeight->GetVal());
348			fDeltaPhiLeptons_NMinusOne->Fill(deltaPhiLeptons,NNLOWeight->GetVal());
349	ceballos	1.1	}
350
351			// NSoftMuons
352			pass = true;
353			for (int k=0;k<nCuts;k++) {
354			if (k != 6)
355			pass = (pass && passCut[k]);
356			}
357			if (pass) {
358	ceballos	1.2	fNSoftMuonsHist_NMinusOne->Fill(SoftMuons->GetEntries(),NNLOWeight->GetVal());
359	ceballos	1.1	}
360
361			//*********************************************************************************************
362			//Plots after all Cuts
363			//*********************************************************************************************
364			if (passAllCuts) {
365	ceballos	1.2	fMinDeltaPhiLeptonMet_afterCuts->Fill(minDeltaPhiMetLepton,NNLOWeight->GetVal());
366			fMtLepton1_afterCuts->Fill(mTW[0],NNLOWeight->GetVal());
367			fMtLepton2_afterCuts->Fill(mTW[1],NNLOWeight->GetVal());
368			fMtHiggs_afterCuts->Fill(mtHiggs,NNLOWeight->GetVal());
369			fLeptonPtPlusMet_afterCuts->Fill(CleanLeptons->At(0)->Pt()+CleanLeptons->At(1)->Pt()+caloMet->Pt(),NNLOWeight->GetVal());
370	ceballos	1.1	}
371
372			delete dilepton;
373			delete SoftMuons;
374			return;
375			}
376
377			//--------------------------------------------------------------------------------------------------
378			void HwwExampleAnalysisMod::SlaveTerminate()
379			{
380
381			// Run finishing code on the computer (slave) that did the analysis. For this
382			// module, we dont do anything here.
383
384			}
385			//--------------------------------------------------------------------------------------------------
386			void HwwExampleAnalysisMod::Terminate()
387			{
388			// Run finishing code on the client computer. For this module, we dont do
389			// anything here.
390
391			}