ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/OSUT3Analysis/AnaTools/src/OSUAnalysis.cc
Revision: 1.7
Committed: Mon Jan 14 23:53:06 2013 UTC (12 years, 4 months ago) by lantonel
Content type: text/plain
Branch: MAIN
CVS Tags: V01-00-01, V00-00-01
Changes since 1.6: +303 -38 lines
Log Message:
new and improved!

File Contents

# User Rev Content
1 ahart 1.1 #include "OSUT3Analysis/AnaTools/interface/OSUAnalysis.h"
2    
3     OSUAnalysis::OSUAnalysis (const edm::ParameterSet &cfg) :
4 ahart 1.2 // Retrieve parameters from the configuration file.
5 ahart 1.1 muons_ (cfg.getParameter<edm::InputTag> ("muons")),
6 lantonel 1.7 electrons_ (cfg.getParameter<edm::InputTag> ("electrons")),
7     channels_ (cfg.getParameter<vector<edm::ParameterSet> >("channels"))
8    
9 ahart 1.1 {
10 lantonel 1.7 TH1::SetDefaultSumw2 ();
11    
12     // Construct Cutflow Objects. These store the results of cut decisions and
13 ahart 1.2 // handle filling cut flow histograms.
14 lantonel 1.7 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 ahart 1.1
112     }
113    
114 ahart 1.2 OSUAnalysis::~OSUAnalysis ()
115 ahart 1.1 {
116 ahart 1.5 // Destroying the CutFlow objects causes the cut flow numbers and time
117     // information to be printed to standard output.
118 lantonel 1.7 for(uint currentChannel = 0; currentChannel != channels_.size(); currentChannel++){
119     delete cutFlows_.at(currentChannel);
120     }
121 ahart 1.1 }
122    
123     void
124     OSUAnalysis::analyze (const edm::Event &event, const edm::EventSetup &setup)
125     {
126 ahart 1.2 // Retrieve necessary collections from the event.
127 ahart 1.1 edm::Handle<BNmuonCollection> muons;
128     event.getByLabel (muons_, muons);
129 lantonel 1.7 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 ahart 1.1
167 lantonel 1.7 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 ahart 1.1 }
268 lantonel 1.7 oneDHists_.at(currentChannel)["numMuons"]->Fill (muonCounter);
269    
270    
271    
272     } //end loop over channel
273    
274     masterCutFlow_->fillCutFlow();
275    
276 ahart 1.1 }
277    
278 lantonel 1.7
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 ahart 1.1 DEFINE_FWK_MODULE(OSUAnalysis);