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(ModNames::gkGoodVertexesName),
  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);

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