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

  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,1.0);
  if (finalstateType == 10 )
    fHWWToMuMuSelection->Fill(-1,1.0);
  else if(finalstateType == 11 )
    fHWWToEESelection->Fill(-1,1.0);
  else if(finalstateType == 12 )
    fHWWToEMuSelection->Fill(-1,1.0);

  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,1.0);
      if (finalstateType == 10 )
        fHWWToMuMuSelection->Fill(k,1.0);
      else if(finalstateType == 11)
        fHWWToEESelection->Fill(k,1.0);
      else if(finalstateType == 12)
        fHWWToEMuSelection->Fill(k,1.0);
    }
  }
  
  //*****************************************************************************************
  //Make Preselection Histograms  
  //*****************************************************************************************
  fLeptonEta->Fill(CleanLeptons->At(0)->Eta(),1.0); 
  fLeptonEta->Fill(CleanLeptons->At(1)->Eta(),1.0);
  fLeptonPtMax->Fill(CleanLeptons->At(0)->Pt(),1.0);
  fLeptonPtMin->Fill(CleanLeptons->At(1)->Pt(),1.0);
  fMetPtHist->Fill(caloMet->Pt(),1.0);                             
  fMetPhiHist->Fill(caloMet->Phi(),1.0);                            
  fDeltaPhiLeptons->Fill(deltaPhiLeptons,1.0);
  fDeltaEtaLeptons->Fill(deltaEtaLeptons,1.0);
  fDileptonMass->Fill(dilepton->Mass(),1.0);    

  //*********************************************************************************************
  // 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,1.0);
  }     
  
  // 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(),1.0);  
  }

  // 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(),1.0);
  }
  
  // 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(),1.0);
    fLeptonPtMin_NMinusOne->Fill(CleanLeptons->At(1)->Pt(),1.0);
    fDeltaPhiLeptons_NMinusOne->Fill(deltaPhiLeptons,1.0); 
  }
  
  // NSoftMuons
  pass = true;
  for (int k=0;k<nCuts;k++) {
    if (k != 6)
      pass = (pass && passCut[k]);    
  }
  if (pass) {
    fNSoftMuonsHist_NMinusOne->Fill(SoftMuons->GetEntries(),1.0);
  }

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