SelMods/src/WBFExampleAnalysisMod.cc

 // $Id $

#include "MitPhysics/SelMods/interface/WBFExampleAnalysisMod.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::WBFExampleAnalysisMod)

//--------------------------------------------------------------------------------------------------
WBFExampleAnalysisMod::WBFExampleAnalysisMod(const char *name, const char *title) : 
  BaseMod(name,title),
  fMetName("NoDefaultNameSet"),
  fCleanJetsName("NoDefaultNameSet"),
  fVertexName(ModNames::gkGoodVertexesName),
  fMet(0),
  fVertices(0),
  fJetPtMax(30),
  fJetPtMin(30),
  fDeltaEtaMin(4),
  fDiJetMassMin(600)
{
  // Constructor.
}

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

//--------------------------------------------------------------------------------------------------
void WBFExampleAnalysisMod::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.

  //Create your histograms here

  //*************************************************************************************************
  // Selection Histograms
  //*************************************************************************************************
  AddTH1(fWBFSelection,"fWBFSelection", ";Cut Number;Number of Events", 7, -1.5, 5.5);
 
  //***********************************************************************************************
  // N-1 Histograms
  //***********************************************************************************************
  //All events
  AddTH1(fWBFPtJetMax_NMinusOne, "fWBFPtJetMax_NMinusOne", ";Pt Jet Max;Number of Events",300,0.,300. );
  AddTH1(fWBFPtJetMin_NMinusOne, "fWBFPtJetMin_NMinusOne", ";Pt Jet Min;Number of Events",300,0.,300. );
  AddTH1(fWBFEta12_NMinusOne,    "fWBFEta12_NMinusOne",    ";eta1*eta2;Number of Events", 2,-1.5,1.5  );
  AddTH1(fWBFdeltaEta_NMinusOne, "fWBFdeltaEta_NMinusOne", ";Delta Eta;Number of Events", 100,0.,10.  );
  AddTH1(fWBFdijetMass_NMinusOne,"fWBFdijetMass_NMinusOne",";DiJet Mass;Number of Events",300,0.,3000.);
  AddTH1(fWBFZVar_NMinusOne,     "fWBFZVar_NMinusOne",     ";Z variable;Number of Events",50,0.,5.    );

  //***********************************************************************************************
  // After all cuts Histograms
  //***********************************************************************************************
  AddTH1(fWBFSSMass_afterCuts,"fWBFSSMass_afterCuts",";SS dilepton Mass;Number of Events",200,0.,400);
  AddTH1(fWBFSSDeltaPhi_afterCuts,"fWBFSSDeltaPhi_afterCuts",";SS DeltaPhi 2l;Number of Events",90,0.,180);
  AddTH1(fWBFOSMass_afterCuts,"fWBFOSMass_afterCuts",";OS dilepton Mass;Number of Events",200,0.,400);
  AddTH1(fWBFOSDeltaPhi_afterCuts,"fWBFOSDeltaPhi_afterCuts",";OS DeltaPhi 2l;Number of Events",90,0.,180);

}

//--------------------------------------------------------------------------------------------------
void WBFExampleAnalysisMod::Process()
{
  //Obtain all the good objects from the event cleaning module
  fVertices = GetObjThisEvt<VertexOArr>(fVertexName);
  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;
  }

  if(CleanJets->GetEntries() < 2) return;

  //*********************************************************************************************
  //Define Cuts
  //*********************************************************************************************
  const int nCuts = 6;
  bool passCut[nCuts] = {false, false, false, false, false, false};
  
  // ptjet max cut
  if(CleanJets->At(0)->Pt() >  fJetPtMax)                  passCut[0] = true;

  // ptjet min cut
  if(CleanJets->At(1)->Pt() >  fJetPtMin)                  passCut[1] = true;
  
  // this cut is always applied
  if(CleanJets->At(0)->Eta()*CleanJets->At(1)->Eta() < 0)  passCut[2] = true;
  
  // deltaEta cut
  double deltaEta = TMath::Abs(CleanJets->At(0)->Eta()-CleanJets->At(1)->Eta());
  if(deltaEta > fDeltaEtaMin)                              passCut[3] = true;
  
  // dijet mass cut
  CompositeParticle dijet;
  dijet.AddDaughter(CleanJets->At(0));
  dijet.AddDaughter(CleanJets->At(1));
  if(dijet.Mass() > fDiJetMassMin)                         passCut[4] = true;
  
  // jet veto cut, use of Zeppenfeld variable
  passCut[5] = true;
  double zVarMin = 30.;
  if(CleanJets->GetEntries() >=2 ){
    for(UInt_t i=2; i<CleanJets->GetEntries(); i++){
      if(CleanJets->At(i)->Pt() <= 30.0) return;
      double zVar = TMath::Abs(CleanJets->At(i)->Eta()-(CleanJets->At(0)->Eta()+CleanJets->At(1)->Eta())/2)/
                    TMath::Abs(CleanJets->At(0)->Eta()-CleanJets->At(1)->Eta());
      if(zVar < zVarMin) zVarMin = zVar;
    }
    if(zVarMin < 1) passCut[5] = false;
  }
  //*********************************************************************************************
  //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
  fWBFSelection->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) {
      fWBFSelection->Fill(k,NNLOWeight->GetVal());
    }
  }
  
  //*********************************************************************************************
  // N-1 Histograms
  //*********************************************************************************************
  // N-1 ptjet max
  bool pass = true;
  for (int k=0;k<nCuts;k++) {
    if (k != 0) {
      pass = (pass && passCut[k]);      
    }
  }
  if (pass) {
    fWBFPtJetMax_NMinusOne->Fill(TMath::Min(CleanJets->At(0)->Pt(),299.999),NNLOWeight->GetVal());
  }     
  
  // N-1 ptjet min
  pass = true;
  for (int k=0;k<nCuts;k++) {
    if (k != 1) {
      pass = (pass && passCut[k]);      
    }
  }
  if (pass) {
    fWBFPtJetMin_NMinusOne->Fill(TMath::Min(CleanJets->At(1)->Pt(),299.999),NNLOWeight->GetVal());
  }     
  
  // N-1 eta1*eta2
  pass = true;
  for (int k=0;k<nCuts;k++) {
    if (k != 2) {
      pass = (pass && passCut[k]);      
    }
  }
  if (pass) {
    fWBFEta12_NMinusOne->Fill(CleanJets->At(0)->Eta()*CleanJets->At(1)->Eta()/
                   TMath::Abs(CleanJets->At(0)->Eta()*CleanJets->At(1)->Eta()),NNLOWeight->GetVal());
  }     
  
  // N-1 deltaEta
  pass = true;
  for (int k=0;k<nCuts;k++) {
    if (k != 3) {
      pass = (pass && passCut[k]);      
    }
  }
  if (pass) {
    fWBFdeltaEta_NMinusOne->Fill(TMath::Min(deltaEta,9.999),NNLOWeight->GetVal());
  }
  
  // N-1 dijet mass
  pass = true;
  for (int k=0;k<nCuts;k++) {
    if (k != 4) {
      pass = (pass && passCut[k]);      
    }
  }
  if (pass) {
    fWBFdijetMass_NMinusOne->Fill(TMath::Min(dijet.Mass(),2999.999),NNLOWeight->GetVal());
  }
  
  // N-1 dijet mass
  pass = true;
  for (int k=0;k<nCuts;k++) {
    if (k != 5) {
      pass = (pass && passCut[k]);      
    }
  }
  if (pass) {
    fWBFZVar_NMinusOne->Fill(TMath::Min(zVarMin,4.999),NNLOWeight->GetVal());
  }

  //*********************************************************************************************
  //Plots after all Cuts
  //*********************************************************************************************
  if (passAllCuts) {
    // Distributions for dileptons events
    if (CleanLeptons->GetEntries() >= 2 && 
        CleanLeptons->At(0)->Pt()  > 20 && CleanLeptons->At(1)->Pt() > 10){

      CompositeParticle dilepton;
      dilepton.AddDaughter(CleanLeptons->At(0));
      dilepton.AddDaughter(CleanLeptons->At(1));
      double deltaPhiLeptons = MathUtils::DeltaPhi(CleanLeptons->At(0)->Phi(), 
                                                   CleanLeptons->At(1)->Phi())* 180.0 / TMath::Pi();
      
      if(CleanLeptons->At(0)->Charge() * CleanLeptons->At(1)->Charge() > 0){  // same-sign
        fWBFSSMass_afterCuts->Fill(TMath::Min(dilepton.Mass(),399.999),NNLOWeight->GetVal());
        fWBFSSDeltaPhi_afterCuts->Fill(deltaPhiLeptons,NNLOWeight->GetVal());
      }
      else { // opposite-sign
        fWBFOSMass_afterCuts->Fill(TMath::Min(dilepton.Mass(),399.999),NNLOWeight->GetVal());
        fWBFOSDeltaPhi_afterCuts->Fill(deltaPhiLeptons,NNLOWeight->GetVal());
      }
    }
  }
  
  return;
}

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

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

}
Revision:	1.3
Committed:	Wed Dec 22 20:23:45 2010 UTC (14 years, 4 months ago) by ceballos
Content type:	text/plain
Branch:	MAIN
Changes since 1.2:	+1 -1 lines
Log Message:	some fixes
#	User	Rev	Content
1	ceballos	1.1	// $Id $
2
3			#include "MitPhysics/SelMods/interface/WBFExampleAnalysisMod.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::WBFExampleAnalysisMod)
18
19			//--------------------------------------------------------------------------------------------------
20			WBFExampleAnalysisMod::WBFExampleAnalysisMod(const char name, const char title) :
21			BaseMod(name,title),
22			fMetName("NoDefaultNameSet"),
23			fCleanJetsName("NoDefaultNameSet"),
24			fVertexName(ModNames::gkGoodVertexesName),
25			fMet(0),
26			fVertices(0),
27			fJetPtMax(30),
28			fJetPtMin(30),
29			fDeltaEtaMin(4),
30			fDiJetMassMin(600)
31			{
32			// Constructor.
33			}
34
35			//--------------------------------------------------------------------------------------------------
36			void WBFExampleAnalysisMod::Begin()
37			{
38			// Run startup code on the client machine. For this module, we dont do
39			// anything here.
40			}
41
42			//--------------------------------------------------------------------------------------------------
43			void WBFExampleAnalysisMod::SlaveBegin()
44			{
45			// Run startup code on the computer (slave) doing the actual analysis. Here,
46			// we typically initialize histograms and other analysis objects and request
47			// branches. For this module, we request a branch of the MitTree.
48
49			//Create your histograms here
50
51			//*************************************************************************************************
52			// Selection Histograms
53			//*************************************************************************************************
54	ceballos	1.2	AddTH1(fWBFSelection,"fWBFSelection", ";Cut Number;Number of Events", 7, -1.5, 5.5);
55	ceballos	1.1
56			//***********************************************************************************************
57			// N-1 Histograms
58			//***********************************************************************************************
59			//All events
60	ceballos	1.2	AddTH1(fWBFPtJetMax_NMinusOne, "fWBFPtJetMax_NMinusOne", ";Pt Jet Max;Number of Events",300,0.,300. );
61			AddTH1(fWBFPtJetMin_NMinusOne, "fWBFPtJetMin_NMinusOne", ";Pt Jet Min;Number of Events",300,0.,300. );
62			AddTH1(fWBFEta12_NMinusOne, "fWBFEta12_NMinusOne", ";eta1*eta2;Number of Events", 2,-1.5,1.5 );
63			AddTH1(fWBFdeltaEta_NMinusOne, "fWBFdeltaEta_NMinusOne", ";Delta Eta;Number of Events", 100,0.,10. );
64			AddTH1(fWBFdijetMass_NMinusOne,"fWBFdijetMass_NMinusOne",";DiJet Mass;Number of Events",300,0.,3000.);
65			AddTH1(fWBFZVar_NMinusOne, "fWBFZVar_NMinusOne", ";Z variable;Number of Events",50,0.,5. );
66	ceballos	1.1
67			//***********************************************************************************************
68			// After all cuts Histograms
69			//***********************************************************************************************
70			AddTH1(fWBFSSMass_afterCuts,"fWBFSSMass_afterCuts",";SS dilepton Mass;Number of Events",200,0.,400);
71	ceballos	1.2	AddTH1(fWBFSSDeltaPhi_afterCuts,"fWBFSSDeltaPhi_afterCuts",";SS DeltaPhi 2l;Number of Events",90,0.,180);
72			AddTH1(fWBFOSMass_afterCuts,"fWBFOSMass_afterCuts",";OS dilepton Mass;Number of Events",200,0.,400);
73			AddTH1(fWBFOSDeltaPhi_afterCuts,"fWBFOSDeltaPhi_afterCuts",";OS DeltaPhi 2l;Number of Events",90,0.,180);
74	ceballos	1.1
75			}
76
77			//--------------------------------------------------------------------------------------------------
78			void WBFExampleAnalysisMod::Process()
79			{
80			//Obtain all the good objects from the event cleaning module
81			fVertices = GetObjThisEvt<VertexOArr>(fVertexName);
82			ParticleOArr CleanLeptons = dynamic_cast<mithep::ParticleOArr>
83			(FindObjThisEvt(ModNames::gkMergedLeptonsName));
84			ObjArray<Jet> CleanJets = dynamic_cast<ObjArray<Jet> >
85			(FindObjThisEvt(fCleanJetsName.Data()));
86			TParameter<Double_t> *NNLOWeight = GetObjThisEvt<TParameter<Double_t> >("NNLOWeight");
87
88			MetCol met = dynamic_cast<ObjArray<Met> >(FindObjThisEvt(fMetName));
89			const Met *caloMet = 0;
90			if (met) {
91			caloMet = met->At(0);
92			} else {
93			cout << "Error: Met Collection " << fMetName << " could not be loaded.\n";
94			return;
95			}
96
97			if(CleanJets->GetEntries() < 2) return;
98
99			//*********************************************************************************************
100			//Define Cuts
101			//*********************************************************************************************
102	ceballos	1.2	const int nCuts = 6;
103			bool passCut[nCuts] = {false, false, false, false, false, false};
104	ceballos	1.1
105			// ptjet max cut
106			if(CleanJets->At(0)->Pt() > fJetPtMax) passCut[0] = true;
107
108			// ptjet min cut
109			if(CleanJets->At(1)->Pt() > fJetPtMin) passCut[1] = true;
110
111			// this cut is always applied
112			if(CleanJets->At(0)->Eta()*CleanJets->At(1)->Eta() < 0) passCut[2] = true;
113
114			// deltaEta cut
115			double deltaEta = TMath::Abs(CleanJets->At(0)->Eta()-CleanJets->At(1)->Eta());
116			if(deltaEta > fDeltaEtaMin) passCut[3] = true;
117
118			// dijet mass cut
119			CompositeParticle dijet;
120			dijet.AddDaughter(CleanJets->At(0));
121			dijet.AddDaughter(CleanJets->At(1));
122			if(dijet.Mass() > fDiJetMassMin) passCut[4] = true;
123	ceballos	1.2
124	ceballos	1.3	// jet veto cut, use of Zeppenfeld variable
125	ceballos	1.2	passCut[5] = true;
126			double zVarMin = 30.;
127			if(CleanJets->GetEntries() >=2 ){
128			for(UInt_t i=2; i<CleanJets->GetEntries(); i++){
129			if(CleanJets->At(i)->Pt() <= 30.0) return;
130			double zVar = TMath::Abs(CleanJets->At(i)->Eta()-(CleanJets->At(0)->Eta()+CleanJets->At(1)->Eta())/2)/
131			TMath::Abs(CleanJets->At(0)->Eta()-CleanJets->At(1)->Eta());
132			if(zVar < zVarMin) zVarMin = zVar;
133			}
134			if(zVarMin < 1) passCut[5] = false;
135			}
136	ceballos	1.1	//*********************************************************************************************
137			//Make Selection Histograms. Number of events passing each level of cut
138			//*********************************************************************************************
139			bool passAllCuts = true;
140			for(int c=0; c<nCuts; c++) passAllCuts = passAllCuts & passCut[c];
141
142			//Cut Selection Histograms
143			fWBFSelection->Fill(-1,NNLOWeight->GetVal());
144			for (int k=0;k<nCuts;k++) {
145			bool pass = true;
146			bool passPreviousCut = true;
147			for (int p=0;p<=k;p++) {
148			pass = (pass && passCut[p]);
149			if (p<k)
150			passPreviousCut = (passPreviousCut&& passCut[p]);
151			}
152
153			if (pass) {
154			fWBFSelection->Fill(k,NNLOWeight->GetVal());
155			}
156			}
157
158			//*********************************************************************************************
159			// N-1 Histograms
160			//*********************************************************************************************
161			// N-1 ptjet max
162			bool pass = true;
163			for (int k=0;k<nCuts;k++) {
164			if (k != 0) {
165			pass = (pass && passCut[k]);
166			}
167			}
168			if (pass) {
169	ceballos	1.2	fWBFPtJetMax_NMinusOne->Fill(TMath::Min(CleanJets->At(0)->Pt(),299.999),NNLOWeight->GetVal());
170	ceballos	1.1	}
171
172			// N-1 ptjet min
173			pass = true;
174			for (int k=0;k<nCuts;k++) {
175			if (k != 1) {
176			pass = (pass && passCut[k]);
177			}
178			}
179			if (pass) {
180	ceballos	1.2	fWBFPtJetMin_NMinusOne->Fill(TMath::Min(CleanJets->At(1)->Pt(),299.999),NNLOWeight->GetVal());
181			}
182
183			// N-1 eta1*eta2
184			pass = true;
185			for (int k=0;k<nCuts;k++) {
186			if (k != 2) {
187			pass = (pass && passCut[k]);
188			}
189			}
190			if (pass) {
191			fWBFEta12_NMinusOne->Fill(CleanJets->At(0)->Eta()*CleanJets->At(1)->Eta()/
192			TMath::Abs(CleanJets->At(0)->Eta()*CleanJets->At(1)->Eta()),NNLOWeight->GetVal());
193	ceballos	1.1	}
194
195			// N-1 deltaEta
196			pass = true;
197			for (int k=0;k<nCuts;k++) {
198			if (k != 3) {
199			pass = (pass && passCut[k]);
200			}
201			}
202			if (pass) {
203			fWBFdeltaEta_NMinusOne->Fill(TMath::Min(deltaEta,9.999),NNLOWeight->GetVal());
204			}
205
206			// N-1 dijet mass
207			pass = true;
208			for (int k=0;k<nCuts;k++) {
209			if (k != 4) {
210			pass = (pass && passCut[k]);
211			}
212			}
213			if (pass) {
214	ceballos	1.2	fWBFdijetMass_NMinusOne->Fill(TMath::Min(dijet.Mass(),2999.999),NNLOWeight->GetVal());
215			}
216
217			// N-1 dijet mass
218			pass = true;
219			for (int k=0;k<nCuts;k++) {
220			if (k != 5) {
221			pass = (pass && passCut[k]);
222			}
223			}
224			if (pass) {
225			fWBFZVar_NMinusOne->Fill(TMath::Min(zVarMin,4.999),NNLOWeight->GetVal());
226	ceballos	1.1	}
227
228			//*********************************************************************************************
229			//Plots after all Cuts
230			//*********************************************************************************************
231			if (passAllCuts) {
232			// Distributions for dileptons events
233			if (CleanLeptons->GetEntries() >= 2 &&
234			CleanLeptons->At(0)->Pt() > 20 && CleanLeptons->At(1)->Pt() > 10){
235
236			CompositeParticle dilepton;
237			dilepton.AddDaughter(CleanLeptons->At(0));
238			dilepton.AddDaughter(CleanLeptons->At(1));
239			double deltaPhiLeptons = MathUtils::DeltaPhi(CleanLeptons->At(0)->Phi(),
240			CleanLeptons->At(1)->Phi())* 180.0 / TMath::Pi();
241
242	ceballos	1.2	if(CleanLeptons->At(0)->Charge() * CleanLeptons->At(1)->Charge() > 0){ // same-sign
243	ceballos	1.1	fWBFSSMass_afterCuts->Fill(TMath::Min(dilepton.Mass(),399.999),NNLOWeight->GetVal());
244			fWBFSSDeltaPhi_afterCuts->Fill(deltaPhiLeptons,NNLOWeight->GetVal());
245			}
246	ceballos	1.2	else { // opposite-sign
247	ceballos	1.1	fWBFOSMass_afterCuts->Fill(TMath::Min(dilepton.Mass(),399.999),NNLOWeight->GetVal());
248			fWBFOSDeltaPhi_afterCuts->Fill(deltaPhiLeptons,NNLOWeight->GetVal());
249			}
250			}
251			}
252
253			return;
254			}
255
256			//--------------------------------------------------------------------------------------------------
257			void WBFExampleAnalysisMod::SlaveTerminate()
258			{
259
260			// Run finishing code on the computer (slave) that did the analysis. For this
261			// module, we dont do anything here.
262
263			}
264			//--------------------------------------------------------------------------------------------------
265			void WBFExampleAnalysisMod::Terminate()
266			{
267			// Run finishing code on the client computer. For this module, we dont do
268			// anything here.
269
270			}