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.);

}

//--------------------------------------------------------------------------------------------------
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() <= 10) 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.4
Committed:	Tue Oct 12 00:44:14 2010 UTC (14 years, 6 months ago) by ceballos
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Mit_015
Changes since 1.3:	+1 -1 lines
Log Message:	typo fixed
#	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
118			//--------------------------------------------------------------------------------------------------
119			void HwwExampleAnalysisMod::Process()
120			{
121			// Process entries of the tree. For this module, we just load the branches and
122			LoadBranch(fMuonBranchName);
123			LoadBranch(fVertexName);
124
125			//Obtain all the good objects from the event cleaning module
126			ObjArray<Muon> CleanMuons = dynamic_cast<ObjArray<Muon> >(FindObjThisEvt(ModNames::gkCleanMuonsName));
127			ParticleOArr CleanLeptons = dynamic_cast<mithep::ParticleOArr>
128			(FindObjThisEvt(ModNames::gkMergedLeptonsName));
129			ObjArray<Jet> CleanJets = dynamic_cast<ObjArray<Jet> >
130			(FindObjThisEvt(fCleanJetsName.Data()));
131	ceballos	1.2	TParameter<Double_t> *NNLOWeight = GetObjThisEvt<TParameter<Double_t> >("NNLOWeight");
132	ceballos	1.1
133			MetCol met = dynamic_cast<ObjArray<Met> >(FindObjThisEvt(fMetName));
134			const Met *caloMet = 0;
135			if (met) {
136			caloMet = met->At(0);
137			} else {
138			cout << "Error: Met Collection " << fMetName << " could not be loaded.\n";
139			return;
140			}
141
142			//***********************************************************************************************
143			//Kinematic PreSelection
144			//***********************************************************************************************
145			// At least two leptons in the event
146			if (CleanLeptons->GetEntries() < 2) return;
147			// Pt1 > 20 && Pt2 > 10
148	ceballos	1.4	if(CleanLeptons->At(0)->Pt() <= 20 \|\| CleanLeptons->At(1)->Pt() <= 10) return;
149	ceballos	1.1	// opposite charge leptons
150			if(CleanLeptons->At(0)->Charge() * CleanLeptons->At(1)->Charge() > 0) return;
151
152			CompositeParticle *dilepton = new CompositeParticle();
153			dilepton->AddDaughter(CleanLeptons->At(0));
154			dilepton->AddDaughter(CleanLeptons->At(1));
155
156			//***********************************************************************************************
157			//Get Dirty Muons: Non-isolated Muons (exclude the clean muons)
158			//***********************************************************************************************
159			ObjArray<Muon> *SoftMuons = new ObjArray<Muon>;
160			for (UInt_t i=0; i<fMuons->GetEntries(); ++i) {
161			const Muon *mu = fMuons->At(i);
162			if(!MuonTools::PassSoftMuonCut(mu, fVertices)) continue;
163
164			bool isCleanMuon = kFALSE;
165			for (UInt_t j=0; j<CleanMuons->GetEntries(); j++) {
166			if(fMuons->At(i) == CleanMuons->At(j) &&
167			CleanMuons->At(j)->Pt() > 10) isCleanMuon = kTRUE;
168			}
169			if(isCleanMuon == kFALSE) SoftMuons->Add(mu);
170			}
171
172			//***********************************************************************************************
173			//Define Event Variables
174			//***********************************************************************************************
175			//delta phi between the 2 leptons in degrees
176			double deltaPhiLeptons = MathUtils::DeltaPhi(CleanLeptons->At(0)->Phi(),
177			CleanLeptons->At(1)->Phi())* 180.0 / TMath::Pi();
178
179			double deltaEtaLeptons = abs(CleanLeptons->At(0)->Eta() - CleanLeptons->At(1)->Eta()) * 180.0 / TMath::Pi();
180
181			double deltaPhiDileptonMet = MathUtils::DeltaPhi(caloMet->Phi(),
182			dilepton->Phi())*180.0 / TMath::Pi();
183
184			double mtHiggs = TMath::Sqrt(2.0dilepton->Pt() caloMet->Pt()*
185			(1.0 - cos(deltaPhiDileptonMet * TMath::Pi() / 180.0)));
186
187			//angle between MET and closest lepton
188			double deltaPhiMetLepton[2] = {MathUtils::DeltaPhi(caloMet->Phi(), CleanLeptons->At(0)->Phi()),
189			MathUtils::DeltaPhi(caloMet->Phi(), CleanLeptons->At(1)->Phi())};
190
191			double mTW[2] = {TMath::Sqrt(2.0CleanLeptons->At(0)->Pt()caloMet->Pt()*
192			(1.0 - cos(deltaPhiMetLepton[0]))),
193			TMath::Sqrt(2.0CleanLeptons->At(1)->Pt()caloMet->Pt()*
194			(1.0 - cos(deltaPhiMetLepton[1])))};
195
196			double minDeltaPhiMetLepton = (deltaPhiMetLepton[0] < deltaPhiMetLepton[1])?
197			deltaPhiMetLepton[0]:deltaPhiMetLepton[1];
198
199			double METdeltaPhilEt = caloMet->Pt();
200			if(minDeltaPhiMetLepton < TMath::Pi()/2.)
201			METdeltaPhilEt = METdeltaPhilEt * sin(minDeltaPhiMetLepton);
202
203			//count the number of central Jets for vetoing
204			int nCentralJets = 0;
205			for(UInt_t i=0; i<CleanJets->GetEntries(); i++){
206			if(TMath::Abs(CleanJets->At(i)->Eta()) < 5.0 &&
207			CleanJets->At(i)->Pt() > 25.0){
208			nCentralJets++;
209			}
210			}
211
212			//Lepton Type
213			int finalstateType = -1;
214			if (CleanLeptons->At(0)->ObjType() == kMuon && CleanLeptons->At(1)->ObjType() == kMuon ){ // mumu
215			finalstateType = 10;
216			} else if(CleanLeptons->At(0)->ObjType() == kElectron && CleanLeptons->At(1)->ObjType() == kElectron ){ // ee
217			finalstateType = 11;
218			} else if((CleanLeptons->At(0)->ObjType() == kElectron && CleanLeptons->At(1)->ObjType() == kMuon) \|\|
219			(CleanLeptons->At(1)->ObjType() == kElectron && CleanLeptons->At(0)->ObjType() == kMuon)) {
220			finalstateType = 12;
221			} else {
222			cerr << "Error: finalstate lepton type not supported\n";
223			}
224
225			//*********************************************************************************************
226			//Define Cuts
227			//*********************************************************************************************
228			const int nCuts = 7;
229			bool passCut[nCuts] = {false, false, false, false, false, false, false};
230
231			if(CleanLeptons->At(0)->Pt() > 20.0 &&
232			CleanLeptons->At(1)->Pt() >= 20.0) passCut[0] = true;
233
234			if(caloMet->Pt() > 20.0) passCut[1] = true;
235
236			if(dilepton->Mass() > 12.0) passCut[2] = true;
237
238			if(nCentralJets < 1) passCut[5] = true;
239
240			if(CleanLeptons->GetEntries() == 2 &&
241			SoftMuons->GetEntries() == 0) passCut[6] = true;
242
243			if (finalstateType == 10 \|\| finalstateType == 11){ // mumu/ee
244			if(fabs(dilepton->Mass()-91.1876) > 15.0) passCut[3] = true;
245			if(METdeltaPhilEt > 35) passCut[4] = true;
246			}
247			else if(finalstateType == 12) { // emu
248			passCut[3] = true;
249			if(METdeltaPhilEt > 20) passCut[4] = true;
250			}
251
252			//*********************************************************************************************
253			//Make Selection Histograms. Number of events passing each level of cut
254			//*********************************************************************************************
255			bool passAllCuts = true;
256			for(int c=0; c<nCuts; c++) passAllCuts = passAllCuts & passCut[c];
257
258			//Cut Selection Histograms
259	ceballos	1.2	fHWWSelection->Fill(-1,NNLOWeight->GetVal());
260	ceballos	1.1	if (finalstateType == 10 )
261	ceballos	1.2	fHWWToMuMuSelection->Fill(-1,NNLOWeight->GetVal());
262	ceballos	1.1	else if(finalstateType == 11 )
263	ceballos	1.2	fHWWToEESelection->Fill(-1,NNLOWeight->GetVal());
264	ceballos	1.1	else if(finalstateType == 12 )
265	ceballos	1.2	fHWWToEMuSelection->Fill(-1,NNLOWeight->GetVal());
266	ceballos	1.1
267			for (int k=0;k<nCuts;k++) {
268			bool pass = true;
269			bool passPreviousCut = true;
270			for (int p=0;p<=k;p++) {
271			pass = (pass && passCut[p]);
272			if (p<k)
273			passPreviousCut = (passPreviousCut&& passCut[p]);
274			}
275
276			if (pass) {
277	ceballos	1.2	fHWWSelection->Fill(k,NNLOWeight->GetVal());
278	ceballos	1.1	if (finalstateType == 10 )
279	ceballos	1.2	fHWWToMuMuSelection->Fill(k,NNLOWeight->GetVal());
280	ceballos	1.1	else if(finalstateType == 11)
281	ceballos	1.2	fHWWToEESelection->Fill(k,NNLOWeight->GetVal());
282	ceballos	1.1	else if(finalstateType == 12)
283	ceballos	1.2	fHWWToEMuSelection->Fill(k,NNLOWeight->GetVal());
284	ceballos	1.1	}
285			}
286
287			//*****************************************************************************************
288			//Make Preselection Histograms
289			//*****************************************************************************************
290	ceballos	1.2	fLeptonEta->Fill(CleanLeptons->At(0)->Eta(),NNLOWeight->GetVal());
291			fLeptonEta->Fill(CleanLeptons->At(1)->Eta(),NNLOWeight->GetVal());
292			fLeptonPtMax->Fill(CleanLeptons->At(0)->Pt(),NNLOWeight->GetVal());
293			fLeptonPtMin->Fill(CleanLeptons->At(1)->Pt(),NNLOWeight->GetVal());
294			fMetPtHist->Fill(caloMet->Pt(),NNLOWeight->GetVal());
295			fMetPhiHist->Fill(caloMet->Phi(),NNLOWeight->GetVal());
296			fDeltaPhiLeptons->Fill(deltaPhiLeptons,NNLOWeight->GetVal());
297			fDeltaEtaLeptons->Fill(deltaEtaLeptons,NNLOWeight->GetVal());
298			fDileptonMass->Fill(dilepton->Mass(),NNLOWeight->GetVal());
299	ceballos	1.1
300			//*********************************************************************************************
301			// N-1 Histograms
302			//*********************************************************************************************
303			bool pass;;
304
305			//N Jet Veto
306			pass = true;
307			for (int k=0;k<nCuts;k++) {
308			if (k != 5) {
309			pass = (pass && passCut[k]);
310			}
311			}
312			if (pass) {
313	ceballos	1.2	fNCentralJets_NMinusOne->Fill(nCentralJets,NNLOWeight->GetVal());
314	ceballos	1.1	}
315
316			// Final Met Cut
317			pass = true;
318			for (int k=0;k<nCuts;k++) {
319			if (k != 4) {
320			pass = (pass && passCut[k]);
321			}
322			}
323			if (pass) {
324	ceballos	1.2	fMetPtHist_NMinusOne->Fill(caloMet->Pt(),NNLOWeight->GetVal());
325	ceballos	1.1	}
326
327			// dilepton mass
328			pass = true;
329			for (int k=0;k<nCuts;k++) {
330			if (k != 2 && k != 3)
331			pass = (pass && passCut[k]);
332			}
333			if (pass) {
334	ceballos	1.2	fDileptonMass_NMinusOne->Fill(dilepton->Mass(),NNLOWeight->GetVal());
335	ceballos	1.1	}
336
337			// Lepton Pt Max, Lepton Pt Min, DeltaPhiLeptons
338			pass = true;
339			for (int k=0;k<nCuts;k++) {
340			pass = (pass && passCut[k]);
341			}
342			if (pass) {
343	ceballos	1.2	fLeptonPtMax_NMinusOne->Fill(CleanLeptons->At(0)->Pt(),NNLOWeight->GetVal());
344			fLeptonPtMin_NMinusOne->Fill(CleanLeptons->At(1)->Pt(),NNLOWeight->GetVal());
345			fDeltaPhiLeptons_NMinusOne->Fill(deltaPhiLeptons,NNLOWeight->GetVal());
346	ceballos	1.1	}
347
348			// NSoftMuons
349			pass = true;
350			for (int k=0;k<nCuts;k++) {
351			if (k != 6)
352			pass = (pass && passCut[k]);
353			}
354			if (pass) {
355	ceballos	1.2	fNSoftMuonsHist_NMinusOne->Fill(SoftMuons->GetEntries(),NNLOWeight->GetVal());
356	ceballos	1.1	}
357
358			//*********************************************************************************************
359			//Plots after all Cuts
360			//*********************************************************************************************
361			if (passAllCuts) {
362	ceballos	1.2	fMinDeltaPhiLeptonMet_afterCuts->Fill(minDeltaPhiMetLepton,NNLOWeight->GetVal());
363			fMtLepton1_afterCuts->Fill(mTW[0],NNLOWeight->GetVal());
364			fMtLepton2_afterCuts->Fill(mTW[1],NNLOWeight->GetVal());
365			fMtHiggs_afterCuts->Fill(mtHiggs,NNLOWeight->GetVal());
366			fLeptonPtPlusMet_afterCuts->Fill(CleanLeptons->At(0)->Pt()+CleanLeptons->At(1)->Pt()+caloMet->Pt(),NNLOWeight->GetVal());
367	ceballos	1.1	}
368
369			delete dilepton;
370			delete SoftMuons;
371			return;
372			}
373
374			//--------------------------------------------------------------------------------------------------
375			void HwwExampleAnalysisMod::SlaveTerminate()
376			{
377
378			// Run finishing code on the computer (slave) that did the analysis. For this
379			// module, we dont do anything here.
380
381			}
382			//--------------------------------------------------------------------------------------------------
383			void HwwExampleAnalysisMod::Terminate()
384			{
385			// Run finishing code on the client computer. For this module, we dont do
386			// anything here.
387
388			}