ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/OSUT3Analysis/AnaTools/src/OSUAnalysis.cc

Comparing UserCode/OSUT3Analysis/AnaTools/src/OSUAnalysis.cc (file contents):
Revision 1.4 by ahart, Mon Aug 27 18:19:37 2012 UTC vs.
Revision 1.7 by lantonel, Mon Jan 14 23:53:06 2013 UTC

#	Line 3 | Line 3
3		OSUAnalysis::OSUAnalysis (const edm::ParameterSet &cfg) :
4		// Retrieve parameters from the configuration file.
5		muons_ (cfg.getParameter<edm::InputTag> ("muons")),
6	<	muonCfg_ (cfg.getParameter<edm::ParameterSet> ("muonCfg"))
6	>	electrons_ (cfg.getParameter<edm::InputTag> ("electrons")),
7	>	channels_  (cfg.getParameter<vector<edm::ParameterSet> >("channels"))
8	>
9		{
10	<	// Construct CutFlow objects. These store the results of cut decisions and
10	>	TH1::SetDefaultSumw2 ();
11	>
12	>	// Construct Cutflow Objects. These store the results of cut decisions and
13		// handle filling cut flow histograms.
14	<	cuts_ = new CutFlow (hists1D_, fs_);
15	<	muonCuts_ = new CutFlow (hists1D_, fs_, "muon");
14	>	masterCutFlow_ = new CutFlow (fs_);
15	>	std::vector<TFileDirectory> directories;
16	>
17	>
18	>	//loop over all channels (event selections)
19	>	for(uint currentChannel = 0; currentChannel != channels_.size(); currentChannel++){
20	>
21	>	//get name of channel
22	>	string channelName  (channels_.at(currentChannel).getParameter<string>("name"));
23	>	channels.push_back(channelName);
24	>	TString channelLabel = channelName;
25	>
26	>
27	>	//create cutFlow for this channel
28	>	cutFlows_.push_back (new CutFlow (fs_, channelName));
29	>
30	>	//book a directory in the output file with the name of the channel
31	>	TFileDirectory subDir = fs_->mkdir( channelName );
32	>	directories.push_back(subDir);
33	>	std::map<std::string, TH1D*> histoMap;
34	>	oneDHists_.push_back(histoMap);
35	>
36	>
37	>	//muon histograms
38	>	oneDHists_.at(currentChannel)["muonPt"]  = directories.at(currentChannel).make<TH1D> ("muonPt",channelLabel+" channel: Muon Transverse Momentum; p_{T} [GeV]", 100, 0, 500);
39	>	oneDHists_.at(currentChannel)["muonEta"] = directories.at(currentChannel).make<TH1D> ("muonEta",channelLabel+" channel: Muon Eta; #eta", 100, -5, 5);
40	>	oneDHists_.at(currentChannel)["muonD0"] = directories.at(currentChannel).make<TH1D> ("muonD0",channelLabel+" channel: Muon d_{0}; d_{0} [cm] ", 1000, -1, 1);
41	>	oneDHists_.at(currentChannel)["muonAbsD0"] = directories.at(currentChannel).make<TH1D> ("muonAbsD0",channelLabel+" channel: Muon d_{0}; |d_{0}| [cm] ", 1000, 0, 1);
42	>	oneDHists_.at(currentChannel)["muonD0Sig"] = directories.at(currentChannel).make<TH1D> ("muonD0Sig",channelLabel+" channel: Muon d_{0} Significance; d_{0} / #sigma_{d_{0}} ", 1000, -10.0, 10.0);
43	>	oneDHists_.at(currentChannel)["muonAbsD0Sig"] = directories.at(currentChannel).make<TH1D> ("muonAbsD0Sig",channelLabel+" channel: Muon d_{0} Significance; |d_{0}| / #sigma_{d_{0}} ", 1000, 0, 10.0);
44	>	oneDHists_.at(currentChannel)["muonIso"] = directories.at(currentChannel).make<TH1D> ("muonIso",channelLabel+" channel: Muon Combined Relative Isolation; rel. iso. ", 1000, 0, 1);
45	>	oneDHists_.at(currentChannel)["numMuons"] = directories.at(currentChannel).make<TH1D> ("numMuons",channelLabel+" channel: Number of Selected Muons; # muons", 10, 0, 10);
46	>
47	>
48	>	//electron histograms
49	>	oneDHists_.at(currentChannel)["electronPt"]  = directories.at(currentChannel).make<TH1D> ("electronPt",channelLabel+" channel: Electron Transverse Momentum; p_{T} [GeV]", 100, 0, 500);
50	>	oneDHists_.at(currentChannel)["electronEta"] = directories.at(currentChannel).make<TH1D> ("electronEta",channelLabel+" channel: Electron Eta; #eta", 50, -5, 5);
51	>	oneDHists_.at(currentChannel)["electronD0"] = directories.at(currentChannel).make<TH1D> ("electronD0",channelLabel+" channel: Electron d_{0}; d_{0} [cm] ", 1000, -1, 1);
52	>	oneDHists_.at(currentChannel)["electronAbsD0"] = directories.at(currentChannel).make<TH1D> ("electronAbsD0",channelLabel+" channel: Electron d_{0}; |d_{0}| [cm] ", 1000, 0, 1);
53	>	oneDHists_.at(currentChannel)["electronD0Sig"] = directories.at(currentChannel).make<TH1D> ("electronD0Sig",channelLabel+" channel: Electron d_{0} Significance; d_{0} / #sigma_{d_{0}} ", 1000, -10.0, 10.0);
54	>	oneDHists_.at(currentChannel)["electronAbsD0Sig"] = directories.at(currentChannel).make<TH1D> ("electronAbsD0Sig",channelLabel+" channel: Electron d_{0} Significance; |d_{0}| / #sigma_{d_{0}} ", 1000, 0, 10.0);
55	>	oneDHists_.at(currentChannel)["electronIso"] = directories.at(currentChannel).make<TH1D> ("electronIso",channelLabel+" channel: Electron Combined Relative Isolation; rel. iso. ", 1000, 0, 1);
56	>	oneDHists_.at(currentChannel)["electronFbrem"] = directories.at(currentChannel).make<TH1D> ("electronFbrem",channelLabel+" channel: Electron Brem. Energy Fraction; fbrem", 1000, 0, 2);
57	>	oneDHists_.at(currentChannel)["numElectrons"] = directories.at(currentChannel).make<TH1D> ("numElectrons",channelLabel+" channel: Number of Selected Electrons; # electrons", 10, 0, 10);
58	>
59	>
60	>
61	>	//get list of cuts for this channel
62	>	vector<edm::ParameterSet> cuts_  (channels_.at(currentChannel).getParameter<vector<edm::ParameterSet> >("cuts"));
63	>
64	>	//loop over and parse all cuts
65	>	for(uint currentCut = 0; currentCut != cuts_.size(); currentCut++){
66	>
67	>	//store input collection for cut
68	>	string inputCollection = cuts_.at(currentCut).getParameter<string> ("inputCollection");
69	>	inputCollections[channelName].push_back(inputCollection);
70	>
71	>	//split cut string into parts and store them
72	>	string cutString = cuts_.at(currentCut).getParameter<string> ("cutString");
73	>	std::vector<string> cutStringVector = splitString(cutString);
74	>	cutVariables[channelName].push_back(cutStringVector.at(0)); // variable to cut on
75	>	cutComparativeOperators[channelName].push_back(cutStringVector.at(1)); // comparison to make
76	>	cutValues[channelName].push_back(atof(cutStringVector.at(2).c_str())); // threshold value to pass cut
77	>
78	>	//get number of objects required to pass cut for event to pass
79	>	string numberRequiredString = cuts_.at(currentCut).getParameter<string> ("numberRequired");
80	>	std::vector<string> numberRequiredVector = splitString(numberRequiredString);
81	>
82	>	// determine number required if comparison contains "="
83	>	int numberRequiredInt = atoi(numberRequiredVector.at(1).c_str());
84	>	if(numberRequiredVector.at(0) == ">") numberRequiredInt++;
85	>	else if(numberRequiredVector.at(0) == "<") numberRequiredInt--;
86	>
87	>
88	>	numberRequiredInts[channelName].push_back(numberRequiredInt); // number of objects required to pass the cut
89	>	eventComparativeOperators[channelName].push_back(numberRequiredVector.at(0)); // comparison to make
90	>
91	>
92	>	if(cuts_.at(currentCut).exists("alias")){
93	>	cutNames[channelName].push_back(cuts_.at(currentCut).getParameter<string> ("alias"));
94	>	}
95	>	else{
96	>	//construct string for cutflow table
97	>	bool plural = numberRequiredInt != 1;
98	>	string collectionString = plural ? inputCollection : inputCollection.substr(0, inputCollection.size()-1);
99	>	string cutName =  numberRequiredString + " " + collectionString + " with " + cutString;
100	>	cutNames[channelName].push_back(cutName);
101	>	}
102	>
103	>	}//end loop over cuts
104	>
105	>
106	>
107	>	}//end loop over channels
108	>
109	>
110	>
111
13	–	// Create additional histograms.
14	–	TH1::SetDefaultSumw2 ();
15	–	hists1D_["muonPhi"] = fs_->make<TH1D> ("muonPhi", ";#phi", 100, -5.0, 5.0);
16	–	hists1D_["muonEta"] = fs_->make<TH1D> ("muonEta", ";#eta", 100, -3.0, 3.0);
17	–	hists1D_["muonPt"] = fs_->make<TH1D> ("muonPt", ";p_{T} (GeV/c)", 100, 0.0, 300.0);
18	–	hists1D_["dimuonMass"] = fs_->make<TH1D> ("dimuonMass", ";m_{#mu^{+} #mu^{-}} (GeV/c^{2})", 90, 30.0, 120.0);
112		}
113
114		OSUAnalysis::~OSUAnalysis ()
115		{
116	<	// Destroying the CutFlow objects causes the cut flow numbers to be printed
117	<	// to standard output, as well as time information.
118	<	delete cuts_;
119	<	delete muonCuts_;
116	>	// Destroying the CutFlow objects causes the cut flow numbers and time
117	>	// information to be printed to standard output.
118	>	for(uint currentChannel = 0; currentChannel != channels_.size(); currentChannel++){
119	>	delete cutFlows_.at(currentChannel);
120	>	}
121		}
122
123		void
#	Line 32 | Line 126 | OSUAnalysis::analyze (const edm::Event &
126		// Retrieve necessary collections from the event.
127		edm::Handle<BNmuonCollection> muons;
128		event.getByLabel (muons_, muons);
129	+	edm::Handle<BNelectronCollection> electrons;
130	+	event.getByLabel (electrons_, electrons);
131	+
132	+
133	+
134	+	map < string, std::vector < std::vector<bool> > > electronIndividualFlags; //for each channel and each cut, vector of electrons that pass that cut
135	+	map < string, std::vector < std::vector<bool> > > electronCumulativeFlags; //for each channel and each cut, vector of electrons that pass all cuts up to and including that cut
136	+	map < string, vector <int> > passingElectronsCounter;                      //for each channel and each cut, number of electrons that pass all cuts up to and including that cut
137	+
138	+	map < string, std::vector < std::vector<bool> > > muonIndividualFlags; //for each channel and each cut, vector of muons that pass that cut
139	+	map < string, std::vector < std::vector<bool> > > muonCumulativeFlags; //for each channel and each cut, vector of muons that pass all cuts up to and including that cut
140	+	map < string, vector <int> > passingMuonsCounter;                      //for each channel and each cut, number of muons that pass all cuts up to and including that cut
141	+
142	+
143	+	//loop over all channels
144	+	for(uint currentChannel = 0; currentChannel != channels.size(); currentChannel++){
145	+	string currentChannelName = channels.at(currentChannel);
146	+
147	+
148	+	//loop over all cuts
149	+	for(uint currentCut = 0; currentCut != cutNames[currentChannelName].size(); currentCut++){
150	+
151	+	electronIndividualFlags[currentChannelName].push_back (vector<bool> ());
152	+	electronCumulativeFlags[currentChannelName].push_back (vector<bool> ());
153	+	muonIndividualFlags[currentChannelName].push_back (vector<bool> ());
154	+	muonCumulativeFlags[currentChannelName].push_back (vector<bool> ());
155	+
156	+
157	+
158	+
159	+	//loop over all electrons and set flags
160	+	for (BNelectronCollection::const_iterator electron = electrons->begin (); electron != electrons->end (); electron++){
161	+
162	+	int electronIndex = electron-electrons->begin();
163	+	bool decision = true;//electron passes if this cut doesn't cut on electrons
164	+
165	+	if(inputCollections[currentChannelName].at(currentCut) == "electrons"){
166
167	<	// Perform selection on objects within the event. The full collection should
168	<	// only be gone through once.
169	<	GoodMuonCollection goodMuons (muonCfg_, muons, muonCuts_);
170	<
171	<	// Calculate cut decisions and store the results in the CutFlow object, just
172	<	// like a map.
173	<	cuts_->at ("minMuons") = goodMuons.size () > 1;
174	<	cuts_->at ("oppositeSign") = cuts_->at ("minMuons") && goodMuons.leadMuon ()->charge * goodMuons.nextToLeadMuon ()->charge < 0;
175	<
176	<	// This causes the cut flow histograms to be filled.
177	<	cuts_->fillCutFlow ();
178	<
179	<	// Fill additional histograms while applying desired cuts.
180	<	if (cuts_->at ("minMuons") && cuts_->at ("oppositeSign"))
181	<	{
182	<	for (GoodMuonCollection::const_iterator goodMuon = goodMuons.begin (); goodMuon != goodMuons.end (); goodMuon++)
183	<	{
184	<	if (goodMuon->pass ("minPt"))
185	<	hists1D_["muonEta"]->Fill (goodMuon->eta);
186	<	if (goodMuon->pass ("maxEta"))
187	<	hists1D_["muonPt"]->Fill (goodMuon->pt);
188	<	if (goodMuon->pass ())
189	<	hists1D_["muonPhi"]->Fill (goodMuon->phi);
190	<	}
191	<	hists1D_["dimuonMass"]->Fill (goodMuons.dimuonMass ());
167	>	string cutVariable = cutVariables[currentChannelName].at(currentCut);
168	>	double value = valueLookup(electron, cutVariable);
169	>
170	>	decision = evaluateComparison(value,cutComparativeOperators[currentChannelName].at(currentCut),cutValues[currentChannelName].at(currentCut));
171	>
172	>	}
173	>	electronIndividualFlags[currentChannelName].at(currentCut).push_back(decision);
174	>
175	>	//set flags for electrons that pass each cut AND all the previous cuts
176	>	bool previousCumulativeFlag = true;
177	>	for(uint previousCut = 0; previousCut != currentCut; previousCut++){
178	>	if(previousCumulativeFlag && electronIndividualFlags[currentChannelName].at(previousCut).at(electronIndex)) previousCumulativeFlag = true;
179	>	else{ previousCumulativeFlag = false; break;}
180	>	}
181	>	electronCumulativeFlags[currentChannelName].at(currentCut).push_back(previousCumulativeFlag && decision);
182	>	}
183	>
184	>
185	>	//loop over all muons and set flags
186	>	for (BNmuonCollection::const_iterator muon = muons->begin (); muon != muons->end (); muon++){
187	>
188	>	int muonIndex = muon-muons->begin();
189	>	bool decision = true;//electron passes if this cut doesn't cut on electrons
190	>
191	>	//set flag to true if cut doesn't cut on muons
192	>	if(inputCollections[currentChannelName].at(currentCut) == "muons"){
193	>
194	>	string cutVariable = cutVariables[currentChannelName].at(currentCut);
195	>	double value = valueLookup(muon, cutVariable);
196	>
197	>	decision = evaluateComparison(value,cutComparativeOperators[currentChannelName].at(currentCut),cutValues[currentChannelName].at(currentCut));
198	>
199	>	}
200	>	muonIndividualFlags[currentChannelName].at(currentCut).push_back(decision);
201	>
202	>	//set flags for muons that pass each cut AND all the previous cuts
203	>	bool previousCumulativeFlag = true;
204	>	for(uint previousCut = 0; previousCut != currentCut; previousCut++){
205	>	if(previousCumulativeFlag && muonIndividualFlags[currentChannelName].at(previousCut).at(muonIndex)) previousCumulativeFlag = true;
206	>	else{ previousCumulativeFlag = false; break;}
207	>	}
208	>	muonCumulativeFlags[currentChannelName].at(currentCut).push_back(previousCumulativeFlag && decision);
209	>	}
210	>
211	>
212	>
213	>	}//end loop over all cuts
214	>
215	>
216	>	//use cumulative flags to apply cuts at event level
217	>	for(uint currentCut = 0; currentCut != cutNames[currentChannelName].size(); currentCut++){
218	>	//loop over all objects and count how many passed the cumulative selection up to this point
219	>	int numberPassing = 0;
220	>	if(inputCollections[currentChannelName].at(currentCut) == "electrons"){
221	>	for (uint electron = 0; electron != electronCumulativeFlags[currentChannelName].at(currentCut).size() ; electron++){
222	>	if(electronCumulativeFlags[currentChannelName].at(currentCut).at(electron)) numberPassing++;
223	>	}
224	>	}
225	>	else if(inputCollections[currentChannelName].at(currentCut) == "muons"){
226	>	for (uint muon = 0; muon != muonCumulativeFlags[currentChannelName].at(currentCut).size() ; muon++){
227	>	if(muonCumulativeFlags[currentChannelName].at(currentCut).at(muon)) numberPassing++;
228	>	}
229	>	}
230	>	bool cutDecision = evaluateComparison(numberPassing,eventComparativeOperators[currentChannelName].at(currentCut),numberRequiredInts[currentChannelName].at(currentCut));
231	>	cutFlows_.at(currentChannel)->at (cutNames[currentChannelName].at(currentCut)) = cutDecision;
232	>	}
233	>
234	>	cutFlows_.at(currentChannel)->fillCutFlow();
235	>
236	>	string lastCutName = cutNames[currentChannelName].at(cutNames[currentChannelName].size()-1);
237	>	if( !cutFlows_.at(currentChannel)->at (lastCutName) ) continue; //don't fill histograms if event didn't pass
238	>
239	>	vector<bool> electronFlags = electronCumulativeFlags[currentChannelName].at(cutNames[currentChannelName].size()-1);
240	>	int electronCounter = 0;
241	>	for (uint electronIndex = 0; electronIndex != electronFlags.size(); electronIndex++){
242	>	if(!electronFlags.at(electronIndex)) continue;
243	>	electronCounter++;
244	>	oneDHists_.at(currentChannel)["electronPt"]->Fill (electrons->at(electronIndex).pt);
245	>	oneDHists_.at(currentChannel)["electronEta"]->Fill (electrons->at(electronIndex).eta);
246	>	oneDHists_.at(currentChannel)["electronD0"]->Fill (electrons->at(electronIndex).correctedD0Vertex);
247	>	oneDHists_.at(currentChannel)["electronAbsD0"]->Fill (fabs(electrons->at(electronIndex).correctedD0Vertex));
248	>	oneDHists_.at(currentChannel)["electronD0Sig"]->Fill (electrons->at(electronIndex).correctedD0Vertex / electrons->at(electronIndex).tkD0err);
249	>	oneDHists_.at(currentChannel)["electronAbsD0Sig"]->Fill (fabs(electrons->at(electronIndex).correctedD0Vertex) / electrons->at(electronIndex).tkD0err);
250	>	oneDHists_.at(currentChannel)["electronIso"]->Fill (((electrons->at(electronIndex).trackIso + electrons->at(electronIndex).caloIso) / electrons->at(electronIndex).pt));
251	>	oneDHists_.at(currentChannel)["electronFbrem"]->Fill (electrons->at(electronIndex).fbrem);
252	>	}
253	>	oneDHists_.at(currentChannel)["numElectrons"]->Fill (electronCounter);
254	>
255	>	vector<bool> muonFlags = muonCumulativeFlags[currentChannelName].at(cutNames[currentChannelName].size()-1);
256	>	int muonCounter = 0;
257	>	for (uint muonIndex = 0; muonIndex != muonFlags.size(); muonIndex++){
258	>	if(!muonFlags.at(muonIndex)) continue;
259	>	muonCounter++;
260	>	oneDHists_.at(currentChannel)["muonPt"]->Fill (muons->at(muonIndex).pt);
261	>	oneDHists_.at(currentChannel)["muonEta"]->Fill (muons->at(muonIndex).eta);
262	>	oneDHists_.at(currentChannel)["muonD0"]->Fill (muons->at(muonIndex).correctedD0Vertex);
263	>	oneDHists_.at(currentChannel)["muonAbsD0"]->Fill (fabs(muons->at(muonIndex).correctedD0Vertex));
264	>	oneDHists_.at(currentChannel)["muonD0Sig"]->Fill (muons->at(muonIndex).correctedD0Vertex / muons->at(muonIndex).tkD0err);
265	>	oneDHists_.at(currentChannel)["muonAbsD0Sig"]->Fill (fabs(muons->at(muonIndex).correctedD0Vertex) / muons->at(muonIndex).tkD0err);
266	>	oneDHists_.at(currentChannel)["muonIso"]->Fill (((muons->at(muonIndex).trackIso + muons->at(muonIndex).caloIso) / muons->at(muonIndex).pt));
267		}
268	+	oneDHists_.at(currentChannel)["numMuons"]->Fill (muonCounter);
269	+
270	+
271	+
272	+	} //end loop over channel
273	+
274	+	masterCutFlow_->fillCutFlow();
275	+
276		}
277
278	+
279	+	bool
280	+	OSUAnalysis::evaluateComparison (double testValue, string comparison, double cutValue){
281	+
282	+	if(comparison == ">")       return testValue >  cutValue;
283	+	else if(comparison == ">=") return testValue >= cutValue;
284	+	else if(comparison == "<")  return testValue <  cutValue;
285	+	else if(comparison == "<=") return testValue <= cutValue;
286	+	else {std::cout << "WARNING: invalid comparison operator!\n"; return false;}
287	+
288	+	}
289	+
290	+	std::vector<std::string>
291	+	OSUAnalysis::splitString (string inputString){
292	+
293	+	std::stringstream stringStream(inputString);
294	+	std::istream_iterator<std::string> begin(stringStream);
295	+	std::istream_iterator<std::string> end;
296	+	std::vector<std::string> stringVector(begin, end);
297	+	return stringVector;
298	+
299	+	}
300	+
301	+	double
302	+	OSUAnalysis::valueLookup (vector<BNelectron>::const_iterator electron, string variable){
303	+
304	+	double value = 0.0;
305	+	if(variable == "pt") value = electron->pt;
306	+	else if(variable == "eta") value = abs(electron->eta);
307	+	else if(variable == "fbrem") value = electron->fbrem;
308	+
309	+
310	+	else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
311	+
312	+
313	+	return value;
314	+	}
315	+
316	+	double
317	+	OSUAnalysis::valueLookup (vector<BNmuon>::const_iterator muon, string variable){
318	+
319	+	double value = 0.0;
320	+	if(variable == "pt") value = muon->pt;
321	+	else if(variable == "eta") value = abs(muon->eta);
322	+
323	+	else{std::cout << "WARNING: invalid variable '" << variable << "'\n"; value = -999;}
324	+
325	+
326	+	return value;
327	+	}
328	+
329	+
330		DEFINE_FWK_MODULE(OSUAnalysis);

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines