CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
HeavyFlavorValidation.cc
Go to the documentation of this file.
1 
8 // Original Author: Zoltan Gecse
9 
10 #include <memory>
20 
36 
38 
41 
43 
44 #include "TLorentzVector.h"
45 
46 using namespace std;
47 using namespace edm;
48 using namespace reco;
49 using namespace l1extra;
50 using namespace trigger;
51 
53  public:
56  private:
57  virtual void beginRun(const Run & iRun, const EventSetup & iSetup);
58  virtual void analyze(const edm::Event&, const edm::EventSetup&);
59  virtual void endJob();
60  int getMotherId( const Candidate * p );
61  void match( MonitorElement * me, vector<LeafCandidate> & from, vector<LeafCandidate> & to, double deltaRMatchingCut, vector<int> & map );
62  void myBook2D( TString name, vector<double> &xBins, TString xLabel, vector<double> &yBins, TString yLabel, TString title);
63  void myBook2D( TString name, vector<double> &xBins, TString xLabel, vector<double> &yBins, TString yLabel){
64  myBook2D( name, xBins, xLabel, yBins, yLabel, name);
65  }
66  void myBookProfile2D( TString name, vector<double> &xBins, TString xLabel, vector<double> &yBins, TString yLabel, TString title);
67  void myBookProfile2D( TString name, vector<double> &xBins, TString xLabel, vector<double> &yBins, TString yLabel){
68  myBookProfile2D( name, xBins, xLabel, yBins, yLabel, name);
69  }
70  void myBook1D( TString name, vector<double> &xBins, TString label, TString title );
71  void myBook1D( TString name, vector<double> &xBins, TString label ){
72  myBook1D( name, xBins, label, name );
73  }
74  string dqmFolder;
82  vector<int> motherIDs;
87  vector<double> deltaEtaBins;
88  vector<double> deltaPhiBins;
89  vector<double> muonPtBins;
90  vector<double> muonEtaBins;
91  vector<double> muonPhiBins;
92  vector<double> dimuonPtBins;
93  vector<double> dimuonEtaBins;
94  vector<double> dimuonDRBins;
96  map<TString, MonitorElement *> ME;
97  vector<pair<string,int> > filterNamesLevels;
98  const double muonMass;
99 };
100 
102 //get parameters
103  dqmFolder(pset.getUntrackedParameter<string>("DQMFolder")),
104  triggerProcessName(pset.getUntrackedParameter<string>("TriggerProcessName")),
105  triggerPathName(pset.getUntrackedParameter<string>("TriggerPathName")),
106  triggerSummaryRAWTag(InputTag( pset.getUntrackedParameter<string>("TriggerSummaryRAW"), "", triggerProcessName)),
107  triggerSummaryAODTag(InputTag( pset.getUntrackedParameter<string>("TriggerSummaryAOD"), "", triggerProcessName)),
108  triggerResultsTag(InputTag( pset.getUntrackedParameter<string>("TriggerResults"), "", triggerProcessName)),
109  recoMuonsTag(pset.getParameter<InputTag>("RecoMuons")),
110  genParticlesTag(pset.getParameter<InputTag>("GenParticles")),
111  motherIDs(pset.getUntrackedParameter<vector<int> >("MotherIDs")),
112  genGlobDeltaRMatchingCut(pset.getUntrackedParameter<double>("GenGlobDeltaRMatchingCut")),
113  globL1DeltaRMatchingCut(pset.getUntrackedParameter<double>("GlobL1DeltaRMatchingCut")),
114  globL2DeltaRMatchingCut(pset.getUntrackedParameter<double>("GlobL2DeltaRMatchingCut")),
115  globL3DeltaRMatchingCut(pset.getUntrackedParameter<double>("GlobL3DeltaRMatchingCut")),
116  deltaEtaBins(pset.getUntrackedParameter<vector<double> >("DeltaEtaBins")),
117  deltaPhiBins(pset.getUntrackedParameter<vector<double> >("DeltaPhiBins")),
118  muonPtBins(pset.getUntrackedParameter<vector<double> >("MuonPtBins")),
119  muonEtaBins(pset.getUntrackedParameter<vector<double> >("MuonEtaBins")),
120  muonPhiBins(pset.getUntrackedParameter<vector<double> >("MuonPhiBins")),
121  dimuonPtBins(pset.getUntrackedParameter<vector<double> >("DimuonPtBins")),
122  dimuonEtaBins(pset.getUntrackedParameter<vector<double> >("DimuonEtaBins")),
123  dimuonDRBins(pset.getUntrackedParameter<vector<double> >("DimuonDRBins")),
124  muonMass(0.106)
125 {}
126 
127 void HeavyFlavorValidation::beginRun(const Run & iRun, const EventSetup & iSetup){
128 //discover HLT configuration
130  bool isChanged;
131  if(hltConfig.init(iRun, iSetup, triggerProcessName, isChanged)){
132  LogDebug("HLTriggerOfflineHeavyFlavor") << "Successfully initialized HLTConfigProvider with process name: "<<triggerProcessName<<endl;
133  }else{
134  LogWarning("HLTriggerOfflineHeavyFlavor") << "Could not initialize HLTConfigProvider with process name: "<<triggerProcessName<<endl;
135  return;
136  }
137  stringstream os;
138  vector<string> triggerNames = hltConfig.triggerNames();
139  for( size_t i = 0; i < triggerNames.size(); i++) {
140  TString triggerName = triggerNames[i];
141  if (triggerName.Contains(triggerPathName)){
142  vector<string> moduleNames = hltConfig.moduleLabels( triggerNames[i] );
143  for( size_t j = 0; j < moduleNames.size(); j++) {
144  TString name = moduleNames[j];
145  if(name.Contains("Filtered")){
146  int level = 0;
147  if(name.Contains("L1"))
148  level = 1;
149  else if(name.Contains("L2"))
150  level = 2;
151  else if(name.Contains("L3"))
152  level = 3;
153  filterNamesLevels.push_back( pair<string,int>(moduleNames[j],level) );
154  os<<" "<<moduleNames[j];
155  }
156  }
157  break;
158  }
159  }
160  if(filterNamesLevels.size()==0){
161  LogDebug("HLTriggerOfflineHeavyFlavor")<<"Bad Trigger Path: "<<triggerPathName<<endl;
162  return;
163  }else{
164  LogDebug("HLTriggerOfflineHeavyFlavor")<<"Trigger Path: "<<triggerPathName<<" has filters:"<<os.str();
165  }
166 
167 //create Monitor Elements
169  if( !dqmStore ){
170  LogError("HLTriggerOfflineHeavyFlavor") << "Could not find DQMStore service\n";
171  return;
172  }
173  dqmStore->setVerbose(0);
175 // Eta Pt Single
176  myBook2D( "genMuon_genEtaPt", muonEtaBins, "#mu eta", muonPtBins, " #mu pT (GeV)");
177  myBook2D( "globMuon_genEtaPt", muonEtaBins, "#mu eta", muonPtBins, " #mu pT (GeV)");
178  myBook2D( "globMuon_recoEtaPt", muonEtaBins, "#mu eta", muonPtBins, " #mu pT (GeV)");
179  for(size_t i=0; i<filterNamesLevels.size(); i++){
180  myBook2D( TString::Format("filt%dMuon_recoEtaPt",int(i+1)), muonEtaBins, "#mu eta", muonPtBins, " #mu pT (GeV)", filterNamesLevels[i].first);
181  }
182  myBook2D( "pathMuon_recoEtaPt", muonEtaBins, "#mu eta", muonPtBins, " #mu pT (GeV)", triggerPathName);
183  myBook2D( "resultMuon_recoEtaPt", muonEtaBins, "#mu eta", muonPtBins, " #mu pT (GeV)");
184 // Eta Pt Single Resolution
185  myBookProfile2D( "resGlobGen_genEtaPt", muonEtaBins, "#mu eta", muonPtBins, " #mu pT (GeV)");
186  for(size_t i=0; i<filterNamesLevels.size(); i++){
187  myBookProfile2D( TString::Format("resFilt%dGlob_recoEtaPt",int(i+1)), muonEtaBins, "#mu eta", muonPtBins, " #mu pT (GeV)", filterNamesLevels[i].first);
188  }
189  myBookProfile2D( "resPathGlob_recoEtaPt", muonEtaBins, "#mu eta", muonPtBins, " #mu pT (GeV)", triggerPathName);
190 // Eta Pt Double
191  myBook2D( "genDimuon_genEtaPt", dimuonEtaBins, "#mu#mu eta", dimuonPtBins, " #mu#mu pT (GeV)");
192  myBook2D( "globDimuon_genEtaPt", dimuonEtaBins, "#mu#mu eta", dimuonPtBins, " #mu#mu pT (GeV)");
193  myBook2D( "globDimuon_recoEtaPt", dimuonEtaBins, "#mu#mu eta", dimuonPtBins, " #mu#mu pT (GeV)");
194  for(size_t i=0; i<filterNamesLevels.size(); i++){
195  myBook2D( TString::Format("filt%dDimuon_recoEtaPt",int(i+1)), dimuonEtaBins, "#mu#mu eta", dimuonPtBins, " #mu#mu pT (GeV)", filterNamesLevels[i].first);
196  }
197  myBook2D( "pathDimuon_recoEtaPt", dimuonEtaBins, "#mu#mu eta", dimuonPtBins, " #mu#mu pT (GeV)", triggerPathName);
198  myBook2D( "resultDimuon_recoEtaPt", dimuonEtaBins, "#mu#mu eta", dimuonPtBins, " #mu#mu pT (GeV)");
199  for(size_t i=0; i<filterNamesLevels.size(); i++){
200  myBook2D( TString::Format("diFilt%dDimuon_recoEtaPt",int(i+1)), dimuonEtaBins, "#mu#mu eta", dimuonPtBins, " #mu#mu pT (GeV)", filterNamesLevels[i].first);
201  }
202  myBook2D( "diPathDimuon_recoEtaPt", dimuonEtaBins, "#mu#mu eta", dimuonPtBins, " #mu#mu pT (GeV)", triggerPathName);
203 // Eta Phi Single
204  myBook2D( "genMuon_genEtaPhi", muonEtaBins, "#mu eta", muonPhiBins, "#mu phi");
205  myBook2D( "globMuon_genEtaPhi", muonEtaBins, "#mu eta", muonPhiBins, "#mu phi");
206  myBook2D( "globMuon_recoEtaPhi", muonEtaBins, "#mu eta", muonPhiBins, "#mu phi");
207  for(size_t i=0; i<filterNamesLevels.size(); i++){
208  myBook2D( TString::Format("filt%dMuon_recoEtaPhi",int(i+1)), muonEtaBins, "#mu eta", muonPhiBins, "#mu phi", filterNamesLevels[i].first);
209  }
210  myBook2D( "pathMuon_recoEtaPhi", muonEtaBins, "#mu eta", muonPhiBins, "#mu phi", triggerPathName);
211  myBook2D( "resultMuon_recoEtaPhi", muonEtaBins, "#mu eta", muonPhiBins, "#mu phi");
212 // Rap Pt Double
213  myBook2D( "genDimuon_genRapPt", dimuonEtaBins, "#mu#mu rapidity", dimuonPtBins, " #mu#mu pT (GeV)");
214  myBook2D( "globDimuon_genRapPt", dimuonEtaBins, "#mu#mu rapidity", dimuonPtBins, " #mu#mu pT (GeV)");
215  myBook2D( "globDimuon_recoRapPt", dimuonEtaBins, "#mu#mu rapidity", dimuonPtBins, " #mu#mu pT (GeV)");
216  for(size_t i=0; i<filterNamesLevels.size(); i++){
217  myBook2D( TString::Format("filt%dDimuon_recoRapPt",int(i+1)), dimuonEtaBins, "#mu#mu rapidity", dimuonPtBins, " #mu#mu pT (GeV)", filterNamesLevels[i].first);
218  }
219  myBook2D( "pathDimuon_recoRapPt", dimuonEtaBins, "#mu#mu rapidity", dimuonPtBins, " #mu#mu pT (GeV)", triggerPathName);
220  myBook2D( "resultDimuon_recoRapPt", dimuonEtaBins, "#mu#mu rapidity", dimuonPtBins, " #mu#mu pT (GeV)");
221  for(size_t i=0; i<filterNamesLevels.size(); i++){
222  myBook2D( TString::Format("diFilt%dDimuon_recoRapPt",int(i+1)), dimuonEtaBins, "#mu#mu rapidity", dimuonPtBins, " #mu#mu pT (GeV)", filterNamesLevels[i].first);
223  }
224  myBook2D( "diPathDimuon_recoRapPt", dimuonEtaBins, "#mu#mu rapidity", dimuonPtBins, " #mu#mu pT (GeV)", triggerPathName);
225 // Pt DR Double
226  myBook2D( "genDimuon_genPtDR", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R at IP");
227  myBook2D( "globDimuon_genPtDR", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R at IP");
228  myBook2D( "globDimuon_recoPtDR", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R at IP");
229  for(size_t i=0; i<filterNamesLevels.size(); i++){
230  myBook2D( TString::Format("filt%dDimuon_recoPtDR",int(i+1)), dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R at IP", filterNamesLevels[i].first);
231  }
232  myBook2D( "pathDimuon_recoPtDR", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R at IP", triggerPathName);
233  for(size_t i=0; i<filterNamesLevels.size(); i++){
234  myBook2D( TString::Format("diFilt%dDimuon_recoPtDR",int(i+1)), dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R at IP", filterNamesLevels[i].first);
235  }
236  myBook2D( "diPathDimuon_recoPtDR", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R at IP", triggerPathName);
237  myBook2D( "resultDimuon_recoPtDR", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R at IP");
238 // Pt DRpos Double
239 // myBook2D( "genDimuon_genPtDRpos", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R in MS");
240 // myBook2D( "globDimuon_genPtDRpos", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R in MS");
241  myBook2D( "globDimuon_recoPtDRpos", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R in MS");
242  for(size_t i=0; i<filterNamesLevels.size(); i++){
243  myBook2D( TString::Format("filt%dDimuon_recoPtDRpos",int(i+1)), dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R in MS", filterNamesLevels[i].first);
244  }
245  myBook2D( "pathDimuon_recoPtDRpos", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R in MS", triggerPathName);
246  for(size_t i=0; i<filterNamesLevels.size(); i++){
247  myBook2D( TString::Format("diFilt%dDimuon_recoPtDRpos",int(i+1)), dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R in MS", filterNamesLevels[i].first);
248  }
249  myBook2D( "diPathDimuon_recoPtDRpos", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R in MS", triggerPathName);
250  myBook2D( "resultDimuon_recoPtDRpos", dimuonPtBins, " #mu#mu pT (GeV)", dimuonDRBins, "#mu#mu #Delta R in MS");
251 
252 /* ME["genDimuon_genPt+Pt-"] = myBook2D( "genDimuon_genPt+Pt-", muonPtBins, "#mu+ pt (GeV)", muonPtBins, "#mu- pT (GeV)");
253  ME["genGlobDimuon_genPt+Pt-"] = myBook2D( "genGlobDimuon_genPt+Pt-", muonPtBins, "#mu+ pt (GeV)", muonPtBins, "#mu- pT (GeV)");
254  ME["genGlobDimuon_recoPt+Pt-"] = myBook2D( "genGlobDimuon_recoPt+Pt-", muonPtBins, "#mu+ pt (GeV)", muonPtBins, "#mu- pT (GeV)");
255  ME["genGlobL1Dimuon_recoPt+Pt-"] = myBook2D( "genGlobL1Dimuon_recoPt+Pt-", muonPtBins, "#mu+ pt (GeV)", muonPtBins, "#mu- pT (GeV)");
256  ME["genGlobL1L2Dimuon_recoPt+Pt-"] = myBook2D( "genGlobL1L2Dimuon_recoPt+Pt-", muonPtBins, "#mu+ pt (GeV)", muonPtBins, "#mu- pT (GeV)");
257  ME["genGlobL1L2L3Dimuon_recoPt+Pt-"] = myBook2D( "genGlobL1L2L3Dimuon_recoPt+Pt-", muonPtBins, "#mu+ pt (GeV)", muonPtBins, "#mu- pT (GeV)");*/
258 
259 // Matching
260  myBook2D( "globGen_deltaEtaDeltaPhi", deltaEtaBins, "#Delta eta", deltaPhiBins, "#Delta phi");
261  for(size_t i=0; i<filterNamesLevels.size(); i++){
262  myBook2D( TString::Format("filt%dGlob_deltaEtaDeltaPhi",int(i+1)), deltaEtaBins, "#Delta eta", deltaPhiBins, "#Delta phi", filterNamesLevels[i].first);
263  }
264  myBook2D( "pathGlob_deltaEtaDeltaPhi", deltaEtaBins, "#Delta eta", deltaPhiBins, "#Delta phi", triggerPathName);
265 // Size of containers
266  vector<double> sizeBins; sizeBins.push_back(10); sizeBins.push_back(0); sizeBins.push_back(10);
267  myBook1D( "genMuon_size", sizeBins, "container size" );
268  myBook1D( "globMuon_size", sizeBins, "container size" );
269  for(size_t i=0; i<filterNamesLevels.size(); i++){
270  myBook1D( TString::Format("filt%dMuon_size",int(i+1)), sizeBins, "container size", filterNamesLevels[i].first);
271  }
272  myBook1D( "pathMuon_size", sizeBins, "container size", triggerPathName );
273 }
274 
276  if( !dqmStore ){
277  LogDebug("HLTriggerOfflineHeavyFlavor")<<"Could not access DQM Store service"<<endl;
278  return;
279  }
280  if(filterNamesLevels.size()==0){
281  return;
282  }
283 //access the containers and create LeafCandidate copies
284  vector<LeafCandidate> genMuons;
286  iEvent.getByLabel(genParticlesTag, genParticles);
287  if(genParticles.isValid()){
288  for(GenParticleCollection::const_iterator p=genParticles->begin(); p!= genParticles->end(); ++p){
289  if( p->status() == 1 && std::abs(p->pdgId())==13 &&
290  ( find( motherIDs.begin(), motherIDs.end(), -1 )!=motherIDs.end() || find( motherIDs.begin(), motherIDs.end(), getMotherId(&(*p)) )!=motherIDs.end() ) ){
291  genMuons.push_back( *p );
292  }
293  }
294  }else{
295  LogDebug("HLTriggerOfflineHeavyFlavor")<<"Could not access GenParticleCollection"<<endl;
296  }
297  sort(genMuons.begin(), genMuons.end(), GreaterByPt<LeafCandidate>());
298  ME["genMuon_size"]->Fill(genMuons.size());
299  LogDebug("HLTriggerOfflineHeavyFlavor")<<"GenParticleCollection from "<<genParticlesTag<<" has size: "<<genMuons.size()<<endl;
300 
301  vector<LeafCandidate> globMuons;
302  vector<LeafCandidate> globMuons_position;
303  Handle<MuonCollection> recoMuonsHandle;
304  iEvent.getByLabel(recoMuonsTag, recoMuonsHandle);
305  if(recoMuonsHandle.isValid()){
306  for(MuonCollection::const_iterator p=recoMuonsHandle->begin(); p!= recoMuonsHandle->end(); ++p){
307  if(p->isGlobalMuon()){
308  globMuons.push_back(*p);
309  globMuons_position.push_back( LeafCandidate( p->charge(), math::XYZTLorentzVector(p->outerTrack()->innerPosition().x(), p->outerTrack()->innerPosition().y(), p->outerTrack()->innerPosition().z(), 0.) ) );
310  }
311  }
312  }else{
313  LogDebug("HLTriggerOfflineHeavyFlavor")<<"Could not access reco Muons"<<endl;
314  }
315  ME["globMuon_size"]->Fill(globMuons.size());
316  LogDebug("HLTriggerOfflineHeavyFlavor")<<"Global Muons from "<<recoMuonsTag<<" has size: "<<globMuons.size()<<endl;
317 
318 // access RAW trigger event
319  vector<vector<LeafCandidate> > muonsAtFilter;
320  vector<vector<LeafCandidate> > muonPositionsAtFilter;
321  for(size_t i=0; i<filterNamesLevels.size(); i++){
322  muonsAtFilter.push_back(vector<LeafCandidate>());
323  muonPositionsAtFilter.push_back(vector<LeafCandidate>());
324  }
325  Handle<TriggerEventWithRefs> rawTriggerEvent;
326  iEvent.getByLabel( triggerSummaryRAWTag, rawTriggerEvent );
327  if( rawTriggerEvent.isValid() ){
328  for(size_t i=0; i<filterNamesLevels.size(); i++){
329  size_t index = rawTriggerEvent->filterIndex(InputTag( filterNamesLevels[i].first, "", triggerProcessName ));
330  if ( index < rawTriggerEvent->size() ){
331  if( filterNamesLevels[i].second==1 ){
332  vector<L1MuonParticleRef> l1Cands;
333  rawTriggerEvent->getObjects( index, TriggerL1Mu, l1Cands );
334  for(size_t j=0; j<l1Cands.size(); j++){
335  muonsAtFilter[i].push_back(*l1Cands[j]);
336  }
337  }else{
338  vector<RecoChargedCandidateRef> hltCands;
339  rawTriggerEvent->getObjects( index, TriggerMuon, hltCands );
340  for(size_t j=0; j<hltCands.size(); j++){
341  muonsAtFilter[i].push_back(*hltCands[j]);
342  if( filterNamesLevels[i].second==2 ){
343  muonPositionsAtFilter[i].push_back( LeafCandidate( hltCands[j]->charge(), math::XYZTLorentzVector(hltCands[j]->track()->innerPosition().x(), hltCands[j]->track()->innerPosition().y(), hltCands[j]->track()->innerPosition().z(), 0.) ) );
344  }
345  }
346  }
347  }
348  ME[TString::Format("filt%dMuon_size",int(i+1))]->Fill(muonsAtFilter[i].size());
349  LogDebug("HLTriggerOfflineHeavyFlavor")<<"Filter \""<<filterNamesLevels[i].first<<"\" has "<<muonsAtFilter[i].size()<<" muons"<<endl;
350  }
351  }else{
352  LogDebug("HLTriggerOfflineHeavyFlavor")<<"Could not access RAWTriggerEvent"<<endl;
353  }
354 
355 // access AOD trigger event
356  vector<LeafCandidate> pathMuons;
357  Handle<TriggerEvent> aodTriggerEvent;
358  iEvent.getByLabel(triggerSummaryAODTag,aodTriggerEvent);
359  if(aodTriggerEvent.isValid()){
360  TriggerObjectCollection allObjects = aodTriggerEvent->getObjects();
361  for(int i=0; i<aodTriggerEvent->sizeFilters(); i++){
362  if(aodTriggerEvent->filterTag(i)==InputTag((filterNamesLevels.end()-1)->first,"",triggerProcessName)){
363  Keys keys = aodTriggerEvent->filterKeys(i);
364  for(size_t j=0; j<keys.size(); j++){
365  pathMuons.push_back( LeafCandidate( allObjects[keys[j]].id()>0 ? 1 : -1, math::PtEtaPhiMLorentzVector( allObjects[keys[j]].pt(), allObjects[keys[j]].eta(), allObjects[keys[j]].phi(), muonMass ) ) );
366  }
367  }
368  }
369  ME["pathMuon_size"]->Fill(pathMuons.size());
370  LogDebug("HLTriggerOfflineHeavyFlavor")<<"Path \""<<triggerPathName<<"\" has "<<pathMuons.size()<<" muons at last filter \""<<(filterNamesLevels.end()-1)->first<<"\""<<endl;
371  }else{
372  LogDebug("HLTriggerOfflineHeavyFlavor")<<"Could not access AODTriggerEvent"<<endl;
373  }
374 
375 // access Trigger Results
376  bool triggerFired = false;
377  Handle<TriggerResults> triggerResults;
378  iEvent.getByLabel(triggerResultsTag,triggerResults);
379  if(triggerResults.isValid()){
380  LogDebug("HLTriggerOfflineHeavyFlavor")<<"Successfully initialized "<<triggerResultsTag<<endl;
381  const edm::TriggerNames & triggerNames = iEvent.triggerNames(*triggerResults);
382  size_t index = triggerNames.triggerIndex(triggerPathName);
383  if( index < triggerNames.size() ){
384  triggerFired = triggerResults->accept( index );
385  }else{
386  LogDebug("HLTriggerOfflineHeavyFlavor")<<triggerResultsTag<<" has no trigger: "<<triggerPathName<<endl;
387  }
388  }else{
389  LogDebug("HLTriggerOfflineHeavyFlavor")<<"Could not initialize "<<triggerResultsTag<<endl;
390  }
391 
392 //create matching maps
393  vector<int> glob_gen(genMuons.size(),-1);
394  match( ME["globGen_deltaEtaDeltaPhi"], genMuons, globMuons, genGlobDeltaRMatchingCut, glob_gen );
395  vector<vector<int> > filt_glob;
396  for(size_t i=0; i<filterNamesLevels.size(); i++){
397  filt_glob.push_back( vector<int>(globMuons.size(),-1) );
398  if( filterNamesLevels[i].second == 1 ){
399  match( ME[TString::Format("filt%dGlob_deltaEtaDeltaPhi",int(i+1))], globMuons_position, muonsAtFilter[i] ,globL1DeltaRMatchingCut, filt_glob[i] );
400  }else if( filterNamesLevels[i].second == 2 ){
401  match( ME[TString::Format("filt%dGlob_deltaEtaDeltaPhi",int(i+1))], globMuons_position, muonPositionsAtFilter[i] ,globL2DeltaRMatchingCut, filt_glob[i] );
402  }else if( filterNamesLevels[i].second == 3 ){
403  match( ME[TString::Format("filt%dGlob_deltaEtaDeltaPhi",int(i+1))], globMuons, muonsAtFilter[i] ,globL3DeltaRMatchingCut, filt_glob[i] );
404  }
405  }
406  vector<int> path_glob(globMuons.size(),-1);
407  if( (filterNamesLevels.end()-1)->second == 1 ){
408  match( ME["pathGlob_deltaEtaDeltaPhi"], globMuons_position, pathMuons ,globL1DeltaRMatchingCut, path_glob );
409  }else if( (filterNamesLevels.end()-1)->second == 2 ){
410  match( ME["pathGlob_deltaEtaDeltaPhi"], globMuons, pathMuons ,globL2DeltaRMatchingCut, path_glob );
411  }else if( (filterNamesLevels.end()-1)->second == 3 ){
412  match( ME["pathGlob_deltaEtaDeltaPhi"], globMuons, pathMuons ,globL3DeltaRMatchingCut, path_glob );
413  }
414 
415 //fill histos
416  bool first = true;
417  for(size_t i=0; i<genMuons.size(); i++){
418  ME["genMuon_genEtaPt"]->Fill(genMuons[i].eta(), genMuons[i].pt());
419  ME["genMuon_genEtaPhi"]->Fill(genMuons[i].eta(), genMuons[i].phi());
420  if(glob_gen[i] != -1){
421  ME["resGlobGen_genEtaPt"]->Fill(genMuons[i].eta(), genMuons[i].pt(), (globMuons[glob_gen[i]].pt()-genMuons[i].pt())/genMuons[i].pt() );
422  ME["globMuon_genEtaPt"]->Fill(genMuons[i].eta(), genMuons[i].pt());
423  ME["globMuon_genEtaPhi"]->Fill(genMuons[i].eta(), genMuons[i].phi());
424  ME["globMuon_recoEtaPt"]->Fill(globMuons[glob_gen[i]].eta(), globMuons[glob_gen[i]].pt());
425  ME["globMuon_recoEtaPhi"]->Fill(globMuons[glob_gen[i]].eta(), globMuons[glob_gen[i]].phi());
426  for(size_t f=0; f<filterNamesLevels.size(); f++){
427  if(filt_glob[f][glob_gen[i]] != -1){
428  ME[TString::Format("resFilt%dGlob_recoEtaPt",int(f+1))]->Fill(globMuons[glob_gen[i]].eta(), globMuons[glob_gen[i]].pt(), (muonsAtFilter[f][filt_glob[f][glob_gen[i]]].pt()-globMuons[glob_gen[i]].pt())/globMuons[glob_gen[i]].pt() );
429  ME[TString::Format("filt%dMuon_recoEtaPt",int(f+1))]->Fill(globMuons[glob_gen[i]].eta(), globMuons[glob_gen[i]].pt());
430  ME[TString::Format("filt%dMuon_recoEtaPhi",int(f+1))]->Fill(globMuons[glob_gen[i]].eta(), globMuons[glob_gen[i]].phi());
431  }else{
432  break;
433  }
434  }
435  if(path_glob[glob_gen[i]] != -1){
436  ME["resPathGlob_recoEtaPt"]->Fill(globMuons[glob_gen[i]].eta(), globMuons[glob_gen[i]].pt(), (pathMuons[path_glob[glob_gen[i]]].pt()-globMuons[glob_gen[i]].pt())/globMuons[glob_gen[i]].pt() );
437  ME["pathMuon_recoEtaPt"]->Fill(globMuons[glob_gen[i]].eta(), globMuons[glob_gen[i]].pt());
438  ME["pathMuon_recoEtaPhi"]->Fill(globMuons[glob_gen[i]].eta(), globMuons[glob_gen[i]].phi());
439  }
440 //highest pt muon
441  if( first ){
442  first = false;
443  if( triggerFired ){
444  ME["resultMuon_recoEtaPt"]->Fill(globMuons[glob_gen[i]].eta(), globMuons[glob_gen[i]].pt());
445  ME["resultMuon_recoEtaPhi"]->Fill(globMuons[glob_gen[i]].eta(), globMuons[glob_gen[i]].phi());
446  }
447  }
448  }
449  }
450 
451 //fill dimuon histograms (highest pT, opposite charge)
452  int secondMuon = 0;
453  for(size_t j=1; j<genMuons.size(); j++){
454  if(genMuons[0].charge()*genMuons[j].charge()==-1){
455  secondMuon = j;
456  break;
457  }
458  }
459  if(secondMuon > 0){
460 // int pos = genMuons[0].charge()>0 ? 0 : secondMuon ;
461 // int neg = genMuons[0].charge()<0 ? 0 : secondMuon ;
462 //two generated
463  double genDimuonPt = (genMuons[0].p4()+genMuons[secondMuon].p4()).pt();
464  double genDimuonEta = (genMuons[0].p4()+genMuons[secondMuon].p4()).eta();
465  double genDimuonRap = (genMuons[0].p4()+genMuons[secondMuon].p4()).Rapidity();
466  double genDimuonDR = deltaR<LeafCandidate,LeafCandidate>( genMuons[0], genMuons[secondMuon] );
467  bool highPt = genMuons[0].pt()>7. && genMuons[secondMuon].pt()>7;
468  ME["genDimuon_genEtaPt"]->Fill( genDimuonEta, genDimuonPt );
469  ME["genDimuon_genRapPt"]->Fill( genDimuonRap, genDimuonPt );
470  if(highPt) ME["genDimuon_genPtDR"]->Fill( genDimuonPt, genDimuonDR );
471 //two global
472  if(glob_gen[0]!=-1 && glob_gen[secondMuon]!=-1){
473  ME["globDimuon_genEtaPt"]->Fill( genDimuonEta, genDimuonPt );
474  ME["globDimuon_genRapPt"]->Fill( genDimuonRap, genDimuonPt );
475  if(highPt) ME["globDimuon_genPtDR"]->Fill( genDimuonPt, genDimuonDR );
476  double globDimuonPt = (globMuons[glob_gen[0]].p4()+globMuons[glob_gen[secondMuon]].p4()).pt();
477  double globDimuonEta = (globMuons[glob_gen[0]].p4()+globMuons[glob_gen[secondMuon]].p4()).eta();
478  double globDimuonRap = (globMuons[glob_gen[0]].p4()+globMuons[glob_gen[secondMuon]].p4()).Rapidity();
479  double globDimuonDR = deltaR<LeafCandidate,LeafCandidate>( globMuons[glob_gen[0]], globMuons[glob_gen[secondMuon]] );
480  double globDimuonDRpos = deltaR<LeafCandidate,LeafCandidate>( globMuons_position[glob_gen[0]], globMuons_position[glob_gen[secondMuon]] );
481  ME["globDimuon_recoEtaPt"]->Fill( globDimuonEta, globDimuonPt );
482  ME["globDimuon_recoRapPt"]->Fill( globDimuonRap, globDimuonPt );
483  if(highPt) ME["globDimuon_recoPtDR"]->Fill( globDimuonPt, globDimuonDR );
484  if(highPt) ME["globDimuon_recoPtDRpos"]->Fill( globDimuonPt, globDimuonDRpos );
485 //two filter objects
486  for(size_t f=0; f<filterNamesLevels.size(); f++){
487  if(filt_glob[f][glob_gen[0]] != -1 && filt_glob[f][glob_gen[secondMuon]] != -1){
488  ME[TString::Format("diFilt%dDimuon_recoEtaPt",int(f+1))]->Fill( globDimuonEta, globDimuonPt );
489  ME[TString::Format("diFilt%dDimuon_recoRapPt",int(f+1))]->Fill( globDimuonRap, globDimuonPt );
490  if(highPt) ME[TString::Format("diFilt%dDimuon_recoPtDR",int(f+1))]->Fill( globDimuonPt, globDimuonDR );
491  if(highPt) ME[TString::Format("diFilt%dDimuon_recoPtDRpos",int(f+1))]->Fill( globDimuonPt, globDimuonDRpos );
492  }else{
493  break;
494  }
495  }
496 //one filter object
497  for(size_t f=0; f<filterNamesLevels.size(); f++){
498  if(filt_glob[f][glob_gen[0]] != -1 || filt_glob[f][glob_gen[secondMuon]] != -1){
499  ME[TString::Format("filt%dDimuon_recoEtaPt",int(f+1))]->Fill( globDimuonEta, globDimuonPt );
500  ME[TString::Format("filt%dDimuon_recoRapPt",int(f+1))]->Fill( globDimuonRap, globDimuonPt );
501  if(highPt) ME[TString::Format("filt%dDimuon_recoPtDR",int(f+1))]->Fill( globDimuonPt, globDimuonDR );
502  if(highPt) ME[TString::Format("filt%dDimuon_recoPtDRpos",int(f+1))]->Fill( globDimuonPt, globDimuonDRpos );
503  }else{
504  break;
505  }
506  }
507 //two path objects
508  if(path_glob[glob_gen[0]] != -1 && path_glob[glob_gen[secondMuon]] != -1){
509  ME["diPathDimuon_recoEtaPt"]->Fill( globDimuonEta, globDimuonPt );
510  ME["diPathDimuon_recoRapPt"]->Fill( globDimuonRap, globDimuonPt );
511  if(highPt) ME["diPathDimuon_recoPtDR"]->Fill( globDimuonPt, globDimuonDR );
512  if(highPt) ME["diPathDimuon_recoPtDRpos"]->Fill( globDimuonPt, globDimuonDRpos );
513  }
514 //one path object
515  if(path_glob[glob_gen[0]] != -1 || path_glob[glob_gen[secondMuon]] != -1){
516  ME["pathDimuon_recoEtaPt"]->Fill( globDimuonEta, globDimuonPt );
517  ME["pathDimuon_recoRapPt"]->Fill( globDimuonRap, globDimuonPt );
518  if(highPt) ME["pathDimuon_recoPtDR"]->Fill( globDimuonPt, globDimuonDR );
519  if(highPt) ME["pathDimuon_recoPtDRpos"]->Fill( globDimuonPt, globDimuonDRpos );
520  }
521 //trigger result
522  if( triggerFired ){
523  ME["resultDimuon_recoEtaPt"]->Fill( globDimuonEta, globDimuonPt );
524  ME["resultDimuon_recoRapPt"]->Fill( globDimuonRap, globDimuonPt );
525  if(highPt) ME["resultDimuon_recoPtDR"]->Fill( globDimuonPt, globDimuonDR );
526  if(highPt) ME["resultDimuon_recoPtDRpos"]->Fill( globDimuonPt, globDimuonDRpos );
527  }
528  }
529  }
530 }
531 
533 }
534 
536  const Candidate* mother = p->mother();
537  if( mother ){
538  if( mother->pdgId() == p->pdgId() ){
539  return getMotherId(mother);
540  }else{
541  return mother->pdgId();
542  }
543  }else{
544  return 0;
545  }
546 }
547 
548 void HeavyFlavorValidation::match( MonitorElement * me, vector<LeafCandidate> & from, vector<LeafCandidate> & to, double dRMatchingCut, vector<int> & map ){
549  vector<double> dR(from.size());
550  for(size_t i=0; i<from.size(); i++){
551  map[i] = -1;
552  dR[i] = 10.;
553  //find closest
554  for(size_t j=0; j<to.size(); j++){
555  double dRtmp = deltaR<double>(from[i].eta(), from[i].phi(), to[j].eta(), to[j].phi());
556  if( dRtmp < dR[i] ){
557  dR[i] = dRtmp;
558  map[i] = j;
559  }
560  }
561  //fill matching histo
562  if( map[i] != -1 ){
563  me->Fill( to[map[i]].eta()-from[i].eta(), deltaPhi<double>(to[map[i]].phi(), from[i].phi()) );
564  }
565  //apply matching cut
566  if( dR[i] > dRMatchingCut ){
567  map[i] = -1;
568  }
569  //remove duplication
570  if( map[i] != -1 ){
571  for(size_t k=0; k<i; k++){
572  if( map[k] != -1 && map[i] == map[k] ){
573  if( dR[i] < dR[k] ){
574  map[k] = -1;
575  }else{
576  map[i] = -1;
577  }
578  break;
579  }
580  }
581  }
582  }
583 }
584 
585 void HeavyFlavorValidation::myBook2D( TString name, vector<double> &ptBins, TString ptLabel, vector<double> &etaBins, TString etaLabel, TString title )
586 {
587 // dqmStore->setCurrentFolder(dqmFolder+"/"+folder);
588  int ptN = ptBins.size()==3 ? (int)ptBins[0]+1 : ptBins.size();
589  Double_t *pt = new Double_t[ ptN ];
590  for(int i=0; i<ptN; i++){
591  pt[i] = ptBins.size()==3 ? ptBins[1] + i*(ptBins[2]-ptBins[1])/ptBins[0] : ptBins[i] ;
592  }
593  int etaN = etaBins.size()==3 ? (int)etaBins[0]+1 : etaBins.size();
594  Double_t *eta = new Double_t[ etaN ];
595  for(int i=0; i<etaN; i++){
596  eta[i] = etaBins.size()==3 ? etaBins[1] + i*(etaBins[2]-etaBins[1])/etaBins[0] : etaBins[i] ;
597  }
598  TH2F *h = new TH2F( name, name, ptN-1, pt, etaN-1, eta );
599  h->SetXTitle(ptLabel);
600  h->SetYTitle(etaLabel);
601  h->SetTitle(title);
602  ME[name] = dqmStore->book2D( name.Data(), h );
603  delete h;
604 }
605 
606 void HeavyFlavorValidation::myBookProfile2D( TString name, vector<double> &ptBins, TString ptLabel, vector<double> &etaBins, TString etaLabel, TString title )
607 {
608 // dqmStore->setCurrentFolder(dqmFolder+"/"+folder);
609  int ptN = ptBins.size()==3 ? (int)ptBins[0]+1 : ptBins.size();
610  Double_t *pt = new Double_t[ ptN ];
611  for(int i=0; i<ptN; i++){
612  pt[i] = ptBins.size()==3 ? ptBins[1] + i*(ptBins[2]-ptBins[1])/ptBins[0] : ptBins[i] ;
613  }
614  int etaN = etaBins.size()==3 ? (int)etaBins[0]+1 : etaBins.size();
615  Double_t *eta = new Double_t[ etaN ];
616  for(int i=0; i<etaN; i++){
617  eta[i] = etaBins.size()==3 ? etaBins[1] + i*(etaBins[2]-etaBins[1])/etaBins[0] : etaBins[i] ;
618  }
619  TProfile2D *h = new TProfile2D( name, name, ptN-1, pt, etaN-1, eta );
620  h->SetXTitle(ptLabel);
621  h->SetYTitle(etaLabel);
622  h->SetTitle(title);
623  ME[name] = dqmStore->bookProfile2D( name.Data(), h );
624  delete h;
625 }
626 
627 void HeavyFlavorValidation::myBook1D( TString name, vector<double> &bins, TString label, TString title )
628 {
629 // dqmStore->setCurrentFolder(dqmFolder+"/"+folder);
630  int binsN = bins.size()==3 ? (int)bins[0]+1 : bins.size();
631  Double_t *myBins = new Double_t[ binsN ];
632  for(int i=0; i<binsN; i++){
633  myBins[i] = bins.size()==3 ? bins[1] + i*(bins[2]-bins[1])/bins[0] : bins[i] ;
634  }
635  TH1F *h = new TH1F( name, name, binsN-1, myBins );
636  h->SetXTitle(label);
637  h->SetTitle(title);
638  ME[name] = dqmStore->book1D( name.Data(), h );
639  delete h;
640 }
641 
643 }
644 
645 //define this as a plug-in
#define LogDebug(id)
int i
Definition: DBlmapReader.cc:9
virtual edm::TriggerNames const & triggerNames(edm::TriggerResults const &triggerResults) const
Definition: Event.cc:208
const std::string & label
Definition: MVAComputer.cc:186
virtual const Candidate * mother(size_type i=0) const =0
return pointer to mother
MonitorElement * book1D(const char *name, const char *title, int nchX, double lowX, double highX)
Book 1D histogram.
Definition: DQMStore.cc:514
enum start value shifted to 81 so as to avoid clashes with PDG codes
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
const std::vector< std::string > & triggerNames() const
names of trigger paths
vector< pair< string, int > > filterNamesLevels
#define abs(x)
Definition: mlp_lapack.h:159
void myBook1D(TString name, vector< double > &xBins, TString label, TString title)
Strings::size_type size() const
Definition: TriggerNames.cc:39
T eta() const
vector< double > dimuonEtaBins
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:7
double charge(const std::vector< uint8_t > &Ampls)
double double double z
void myBook2D(TString name, vector< double > &xBins, TString xLabel, vector< double > &yBins, TString yLabel)
void Fill(long long x)
U second(std::pair< T, U > const &p)
PtEtaPhiMLorentzVectorD PtEtaPhiMLorentzVector
Lorentz vector with cartesian internal representation.
Definition: LorentzVector.h:26
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:30
Definition: ME.h:11
int getMotherId(const Candidate *p)
int iEvent
Definition: GenABIO.cc:243
unsigned int triggerIndex(std::string const &name) const
Definition: TriggerNames.cc:32
void match(MonitorElement *me, vector< LeafCandidate > &from, vector< LeafCandidate > &to, double deltaRMatchingCut, vector< int > &map)
int j
Definition: DBlmapReader.cc:9
double f[11][100]
void setVerbose(unsigned level)
Definition: DQMStore.cc:196
bool first
Definition: L1TdeRCT.cc:79
bool isValid() const
Definition: HandleBase.h:76
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:355
HeavyFlavorValidation(const edm::ParameterSet &)
int k[5][pyjets_maxn]
const std::vector< std::string > & moduleLabels(unsigned int trigger) const
label(s) of module(s) on a trigger path
virtual int pdgId() const =0
PDG identifier.
std::vector< TriggerObject > TriggerObjectCollection
collection of trigger physics objects (e.g., all isolated muons)
Definition: TriggerObject.h:83
static std::string from(" from ")
std::vector< size_type > Keys
bool init(const edm::Run &iRun, const edm::EventSetup &iSetup, const std::string &processName, bool &changed)
virtual void beginRun(const Run &iRun, const EventSetup &iSetup)
void myBookProfile2D(TString name, vector< double > &xBins, TString xLabel, vector< double > &yBins, TString yLabel, TString title)
void myBook1D(TString name, vector< double > &xBins, TString label)
void myBook2D(TString name, vector< double > &xBins, TString xLabel, vector< double > &yBins, TString yLabel, TString title)
std::pair< typename Association::data_type::first_type, double > match(Reference key, Association association, bool bestMatchByMaxValue)
Generic matching function.
Definition: Utils.h:6
tuple level
Definition: testEve_cfg.py:34
Definition: DDAxes.h:10
map< TString, MonitorElement * > ME
MonitorElement * book2D(const char *name, const char *title, int nchX, double lowX, double highX, int nchY, double lowY, double highY)
Book 2D histogram.
Definition: DQMStore.cc:642
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4
tuple size
Write out results.
void myBookProfile2D(TString name, vector< double > &xBins, TString xLabel, vector< double > &yBins, TString yLabel)
void setCurrentFolder(const std::string &fullpath)
Definition: DQMStore.cc:232
Definition: Run.h:32
virtual void analyze(const edm::Event &, const edm::EventSetup &)
MonitorElement * bookProfile2D(const char *name, const char *title, int nchX, double lowX, double highX, int nchY, double lowY, double highY, int nchZ, double lowZ, double highZ, const char *option="s")
Definition: DQMStore.cc:972
Definition: DDAxes.h:10