CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
IsoTrig.cc
Go to the documentation of this file.
1 // -*- C++ -*-//
2 // Package: IsoTrig
3 // Class: IsoTrig
4 //
12 //
13 // Original Author: Ruchi Gupta
14 // Created: Fri May 25 12:02:48 CDT 2012
15 // $Id$
16 //
17 //
18 #include "IsoTrig.h"
19 
20 
21 
23 
28 
29 IsoTrig::IsoTrig(const edm::ParameterSet& iConfig) : changed(false) {
30  //now do whatever initialization is neededOA
31  trigNames = iConfig.getUntrackedParameter<std::vector<std::string> >("Triggers");
32  PixcandTag_ = iConfig.getParameter<edm::InputTag> ("PixcandTag");
33  L1candTag_ = iConfig.getParameter<edm::InputTag> ("L1candTag");
34  L2candTag_ = iConfig.getParameter<edm::InputTag> ("L2candTag");
35  doL2L3 = iConfig.getUntrackedParameter<bool>("DoL2L3",true);
36  doTiming = iConfig.getUntrackedParameter<bool>("DoTimingTree",true);
37  doMipCutTree = iConfig.getUntrackedParameter<bool>("DoMipCutTree",true);
38  doTrkResTree = iConfig.getUntrackedParameter<bool>("DoTrkResTree",true);
39  doChgIsolTree = iConfig.getUntrackedParameter<bool>("DoChgIsolTree",true);
40  doStudyIsol = iConfig.getUntrackedParameter<bool>("DoStudyIsol",true);
41  verbosity = iConfig.getUntrackedParameter<int>("Verbosity",0);
42  theTrackQuality = iConfig.getUntrackedParameter<std::string>("TrackQuality","highPurity");
43  processName = iConfig.getUntrackedParameter<std::string>("ProcessName","HLT");
45  selectionParameters.minPt = iConfig.getUntrackedParameter<double>("MinTrackPt", 10.0);
46  selectionParameters.minQuality = trackQuality_;
47  selectionParameters.maxDxyPV = iConfig.getUntrackedParameter<double>("MaxDxyPV", 0.2);
48  selectionParameters.maxDzPV = iConfig.getUntrackedParameter<double>("MaxDzPV", 5.0);
49  selectionParameters.maxChi2 = iConfig.getUntrackedParameter<double>("MaxChi2", 5.0);
50  selectionParameters.maxDpOverP = iConfig.getUntrackedParameter<double>("MaxDpOverP", 0.1);
51  selectionParameters.minOuterHit = iConfig.getUntrackedParameter<int>("MinOuterHit", 4);
52  selectionParameters.minLayerCrossed = iConfig.getUntrackedParameter<int>("MinLayerCrossed", 8);
53  selectionParameters.maxInMiss = iConfig.getUntrackedParameter<int>("MaxInMiss", 0);
54  selectionParameters.maxOutMiss = iConfig.getUntrackedParameter<int>("MaxOutMiss", 0);
55  dr_L1 = iConfig.getUntrackedParameter<double>("IsolationL1",1.0);
56  a_coneR = iConfig.getUntrackedParameter<double>("ConeRadius",34.98);
57  a_charIsoR = a_coneR + 28.9;
58  a_neutIsoR = a_charIsoR*0.726;
59  a_mipR = iConfig.getUntrackedParameter<double>("ConeRadiusMIP",14.0);
60  a_neutR1 = iConfig.getUntrackedParameter<double>("ConeRadiusNeut1",21.0);
61  a_neutR2 = iConfig.getUntrackedParameter<double>("ConeRadiusNeut2",29.0);
62  cutMip = iConfig.getUntrackedParameter<double>("MIPCut", 1.0);
63  cutCharge = iConfig.getUntrackedParameter<double>("ChargeIsolation", 2.0);
64  cutNeutral = iConfig.getUntrackedParameter<double>("NeutralIsolation", 2.0);
65  minRunNo = iConfig.getUntrackedParameter<int>("minRun");
66  maxRunNo = iConfig.getUntrackedParameter<int>("maxRun");
67  pixelTracksSources_ = iConfig.getUntrackedParameter<std::vector<edm::InputTag> >("PixelTracksSources");
68  pixelIsolationConeSizeAtEC_ = iConfig.getUntrackedParameter<std::vector<double> >("PixelIsolationConeSizeAtEC");
69  minPTrackValue_ = iConfig.getUntrackedParameter<double>("MinPTrackValue");
70  vtxCutSeed_ = iConfig.getUntrackedParameter<double>("VertexCutSeed");
71  vtxCutIsol_ = iConfig.getUntrackedParameter<double>("VertexCutIsol");
72  tauUnbiasCone_ = iConfig.getUntrackedParameter<double>("TauUnbiasCone");
73  prelimCone_ = iConfig.getUntrackedParameter<double>("PrelimCone");
74  // define tokens for access
75  tok_lumi = consumes<LumiDetails, edm::InLumi>(edm::InputTag("lumiProducer"));
76  edm::InputTag triggerEvent_ ("hltTriggerSummaryAOD","",processName);
77  tok_trigEvt = consumes<trigger::TriggerEvent>(triggerEvent_);
78  edm::InputTag theTriggerResultsLabel ("TriggerResults","",processName);
79  tok_trigRes = consumes<edm::TriggerResults>(theTriggerResultsLabel);
80  tok_genTrack_ = consumes<reco::TrackCollection>(edm::InputTag("generalTracks"));
81  tok_recVtx_ = consumes<reco::VertexCollection>(edm::InputTag("offlinePrimaryVertices"));
82  tok_bs_ = consumes<reco::BeamSpot>(edm::InputTag("offlineBeamSpot"));
83  tok_EB_ = consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit","EcalRecHitsEB"));
84  tok_EE_ = consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit","EcalRecHitsEE"));
85  tok_hbhe_ = consumes<HBHERecHitCollection>(edm::InputTag("hbhereco"));
86  tok_pixtk_ = consumes<reco::IsolatedPixelTrackCandidateCollection>(PixcandTag_);
87  tok_l1cand_ = consumes<trigger::TriggerFilterObjectWithRefs>(L1candTag_);
88  tok_l2cand_ = consumes<reco::IsolatedPixelTrackCandidateCollection>(L2candTag_);
89  if(doTiming){
90  tok_verthb_ = consumes<reco::VertexCollection>(edm::InputTag("hltHITPixelVerticesHB"));
91  tok_verthe_ = consumes<reco::VertexCollection>(edm::InputTag("hltHITPixelVerticesHB"));
92  tok_hlt_ = consumes<trigger::TriggerFilterObjectWithRefs>(edm::InputTag("hltL1sL1SingleJet68"));
93  }
94  if(doChgIsolTree) {
95  for (unsigned int k=0; k<pixelTracksSources_.size(); ++k) {
96  // edm::InputTag pix (pixelTracksSources_[k],"",processName);
97  // tok_pixtks_.push_back(consumes<reco::TrackCollection>(pix));
98  tok_pixtks_.push_back(consumes<reco::TrackCollection>(pixelTracksSources_[k]));
99  }
100  }
101  if (verbosity>=0) {
102  std::cout <<"Parameters read from config file \n"
103  <<"\t minPt " << selectionParameters.minPt
104  <<"\t theTrackQuality " << theTrackQuality
105  <<"\t minQuality " << selectionParameters.minQuality
106  <<"\t maxDxyPV " << selectionParameters.maxDxyPV
107  <<"\t maxDzPV " << selectionParameters.maxDzPV
108  <<"\t maxChi2 " << selectionParameters.maxChi2
109  <<"\t maxDpOverP " << selectionParameters.maxDpOverP
110  <<"\t minOuterHit " << selectionParameters.minOuterHit
111  <<"\t minLayerCrossed " << selectionParameters.minLayerCrossed
112  <<"\t maxInMiss " << selectionParameters.maxInMiss
113  <<"\t maxOutMiss " << selectionParameters.maxOutMiss
114  <<"\t a_coneR " << a_coneR
115  <<"\t a_charIsoR " << a_charIsoR
116  <<"\t a_neutIsoR " << a_neutIsoR
117  <<"\t a_mipR " << a_mipR
118  <<"\t a_neutR " << a_neutR1 << ":" << a_neutR2
119  <<"\t cuts (MIP " << cutMip << " : Charge " << cutCharge
120  <<" : Neutral " << cutNeutral << ")"
121  << std::endl;
122  std::cout <<"Charge Isolation parameters:"
123  <<"\t minPTrackValue " << minPTrackValue_
124  <<"\t vtxCutSeed " << vtxCutSeed_
125  <<"\t vtxCutIsol " << vtxCutIsol_
126  <<"\t tauUnbiasCone " << tauUnbiasCone_
127  <<"\t prelimCone " << prelimCone_
128  <<"\t pixelIsolationConeSizeAtEC";
129  for (unsigned int k=0; k<pixelIsolationConeSizeAtEC_.size(); ++k)
130  std::cout << " " << pixelIsolationConeSizeAtEC_[k];
131  std::cout << std::endl;
132  }
133  double pl[] = {20,30,40,60,80,120};
134  for (int i=0; i<6; ++i) pLimits[i] = pl[i];
135  rEB_ = 123.8;
136  zEE_ = 317.0;
137 }
138 
140  // do anything here that needs to be done at desctruction time
141  // (e.g. close files, deallocate resources etc.)
142 
143 }
144 
146  //The following says we do not know what parameters are allowed so do no validation
147  // Please change this to state exactly what you do use, even if it is no parameters
149  desc.setUnknown();
150  descriptions.addDefault(desc);
151 }
152 
153 void IsoTrig::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
154  if (verbosity%10 > 0) std::cout << "Event starts====================================" << std::endl;
155  int RunNo = iEvent.id().run();
156  int EvtNo = iEvent.id().event();
157  int Lumi = iEvent.luminosityBlock();
158  int Bunch = iEvent.bunchCrossing();
159 
160  iSetup.get<IdealMagneticFieldRecord>().get(bFieldH);
161  iSetup.get<CaloGeometryRecord>().get(pG);
162  const MagneticField *bField = bFieldH.product();
163  GlobalVector BField=bField->inTesla(GlobalPoint(0,0,0));
164  bfVal = BField.mag();
165 
166  trigger::TriggerEvent triggerEvent;
167  edm::Handle<trigger::TriggerEvent> triggerEventHandle;
168  iEvent.getByToken(tok_trigEvt, triggerEventHandle);
169  if (!triggerEventHandle.isValid()) {
170  if (verbosity%10 > 0) std::cout << "Error! Can't get the product hltTriggerSummaryAOD"<< std::endl;
171  } else {
172  triggerEvent = *(triggerEventHandle.product());
173  }
174  const trigger::TriggerObjectCollection& TOC(triggerEvent.getObjects());
177  iEvent.getByToken(tok_trigRes, triggerResults);
178 
180  iEvent.getByToken(tok_genTrack_, trkCollection);
181 
184 
185  iEvent.getByToken(tok_hbhe_, hbhe);
186 
187  iEvent.getByToken(tok_recVtx_, recVtxs);
188  iEvent.getByToken(tok_bs_, beamSpotH);
189  if (recVtxs->size()>0 && !((*recVtxs)[0].isFake())) {
190  leadPV = math::XYZPoint( (*recVtxs)[0].x(),(*recVtxs)[0].y(), (*recVtxs)[0].z() );
191  } else if (beamSpotH.isValid()) {
192  leadPV = beamSpotH->position();
193  }
194  if ((verbosity/100)%10>0) {
195  std::cout << "Primary Vertex " << leadPV;
196  if (beamSpotH.isValid()) std::cout << " Beam Spot "
197  << beamSpotH->position();
198  std::cout << std::endl;
199  }
200  pixelTrackRefsHE.clear(); pixelTrackRefsHB.clear();
201  for (unsigned int iPix=0; iPix<pixelTracksSources_.size(); iPix++) {
203  iEvent.getByToken(tok_pixtks_[iPix],iPixCol);
204  if(iPixCol.isValid()){
205  for (reco::TrackCollection::const_iterator pit=iPixCol->begin(); pit!=iPixCol->end(); pit++) {
206  if(iPix==0)
207  pixelTrackRefsHB.push_back(reco::TrackRef(iPixCol,pit-iPixCol->begin()));
208  pixelTrackRefsHE.push_back(reco::TrackRef(iPixCol,pit-iPixCol->begin()));
209  }
210  }
211  }
212  if(doTiming)
213  studyTiming(iEvent);
214  for (unsigned int ifilter=0; ifilter<triggerEvent.sizeFilters();
215  ++ifilter) {
216  std::string FilterNames[7] = {"hltL1sL1SingleJet68", "hltIsolPixelTrackL2FilterHE", "ecalIsolPixelTrackFilterHE", "hltIsolPixelTrackL3FilterHE",
217  "hltIsolPixelTrackL2FilterHB", "ecalIsolPixelTrackFilterHB", "hltIsolPixelTrackL3FilterHB"};
218  std::string label = triggerEvent.filterTag(ifilter).label();
219  for(int i=0; i<7; i++) {
220  if(label==FilterNames[i]) h_Filters->Fill(i);
221  }
222  }
223  edm::InputTag lumiProducer("LumiProducer", "", "RECO");
225  iEvent.getLuminosityBlock().getByToken(tok_lumi, Lumid);
226  float mybxlumi=-1;
227  if (Lumid.isValid())
228  mybxlumi=Lumid->lumiValue(LumiDetails::kOCC1,iEvent.bunchCrossing())*6.37;
229  if (verbosity%10 > 0)
230  std::cout << "RunNo " << RunNo << " EvtNo " << EvtNo << " Lumi " << Lumi
231  << " Bunch " << Bunch << " mybxlumi " << mybxlumi << std::endl;
232 
233  if (!triggerResults.isValid()) {
234  std::cout << "Error! Can't get the product triggerResults" << std::endl;
235  // boost::shared_ptr<cms::Exception> const & error = triggerResults.whyFailed();
236  // edm::LogWarning(error->category()) << error->what();
237  } else {
238  std::vector<std::string> modules;
239  h_nHLT->Fill(triggerResults->size());
240  const edm::TriggerNames & triggerNames = iEvent.triggerNames(*triggerResults);
241 
242  const std::vector<std::string> & triggerNames_ = triggerNames.triggerNames();
243  // std::cout << "number of HLTs " << triggerNames_.size() << std::endl;
244  int hlt(-1), preL1(-1), preHLT(-1), prescale(-1);
245  for (unsigned int i=0; i<triggerResults->size(); i++) {
246  unsigned int triggerindx = hltConfig_.triggerIndex(triggerNames_[i]);
247  const std::vector<std::string>& moduleLabels(hltConfig_.moduleLabels(triggerindx));
248 
249  for (unsigned int in=0; in<trigNames.size(); ++in) {
250  // if (triggerNames_[i].find(trigNames[in].c_str())!=std::string::npos || triggerNames_[i]==" ") {
251  if (triggerNames_[i].find(trigNames[in].c_str())!=std::string::npos) {
252  if (verbosity%10 > 0) std::cout << "trigger that i want " << triggerNames_[i] << " accept " << triggerResults->accept(i) << std::endl;
253  hlt = triggerResults->accept(i);
254  h_HLT -> Fill(hlt);
255  // if (hlt>0 || triggerNames_[i]==" ") {
256  if (hlt>0) {
258  iEvent.getByToken(tok_pixtk_,Pixcands);
260  iEvent.getByToken(tok_l1cand_, L1cands);
261 
262  const std::pair<int,int> prescales(hltConfig_.prescaleValues(iEvent,iSetup,triggerNames_[i]));
263  preL1 = prescales.first;
264  preHLT = prescales.second;
265  prescale = preL1*preHLT;
266  if (verbosity%10 > 0)
267  std::cout << triggerNames_[i] << " accept " << hlt << " preL1 "
268  << preL1 << " preHLT " << preHLT << std::endl;
269 
270  for (int iv=0; iv<3; ++iv) vec[iv].clear();
271  if (TrigList.find(RunNo) != TrigList.end() ) {
272  TrigList[RunNo] += 1;
273  } else {
274  TrigList.insert(std::pair<unsigned int, unsigned int>(RunNo,1));
275  TrigPreList.insert(std::pair<unsigned int, std::pair<int, int>>(RunNo,prescales));
276  }
277  //loop over all trigger filters in event (i.e. filters passed)
278  for (unsigned int ifilter=0; ifilter<triggerEvent.sizeFilters();
279  ++ifilter) {
280  std::vector<int> Keys;
281  std::string label = triggerEvent.filterTag(ifilter).label();
282  //loop over keys to objects passing this filter
283  for (unsigned int imodule=0; imodule<moduleLabels.size();
284  imodule++) {
285  if (label.find(moduleLabels[imodule]) != std::string::npos) {
286  if (verbosity%10 > 0) std::cout << "FILTERNAME " << label << std::endl;
287  for (unsigned int ifiltrKey=0; ifiltrKey<triggerEvent.filterKeys(ifilter).size(); ++ifiltrKey) {
288  Keys.push_back(triggerEvent.filterKeys(ifilter)[ifiltrKey]);
289  const trigger::TriggerObject& TO(TOC[Keys[ifiltrKey]]);
290  math::XYZTLorentzVector v4(TO.px(), TO.py(), TO.pz(), TO.energy());
291  if (label.find("L2Filter") != std::string::npos) {
292  vec[1].push_back(v4);
293  } else if (label.find("L3Filter") != std::string::npos) {
294  vec[2].push_back(v4);
295  } else {
296  vec[0].push_back(v4);
297  h_L1ObjEnergy->Fill(TO.energy());
298  }
299  if (verbosity%10 > 0)
300  std::cout << "key " << ifiltrKey << " : pt " << TO.pt() << " eta " << TO.eta() << " phi " << TO.phi() << " mass " << TO.mass() << " Id " << TO.id() << std::endl;
301  }
302  }
303  }
304  }
305  std::vector<reco::TrackCollection::const_iterator> goodTks;
306  if (doL2L3) {
307  h_nL3Objs -> Fill(vec[2].size());
308  studyTrigger(trkCollection, goodTks);
309  } else {
310  if (trkCollection.isValid()) {
311  reco::TrackCollection::const_iterator trkItr;
312  for (trkItr=trkCollection->begin();
313  trkItr!=trkCollection->end(); trkItr++)
314  goodTks.push_back(trkItr);
315  }
316  }
317  // Now study isolation etc
318  if(doStudyIsol) studyIsolation(trkCollection, goodTks);
319 
320  if(doTrkResTree) StudyTrkEbyP(trkCollection);
321 
322  std::pair<double,double> etaphi = etaPhiTrigger();
324  iEvent.getByToken(tok_l2cand_,L2cands);
325  if (!L2cands.isValid()) {
326  if (verbosity%10 > 0) std::cout << " trigCand is not valid " << std::endl;
327  } else {
328  if(doMipCutTree) studyMipCut(trkCollection, L2cands);
329  }
330  if (pixelTracksSources_.size()>0)
331  if(doChgIsolTree && pixelTrackRefsHE.size()>0) chgIsolation(etaphi.first, etaphi.second, trkCollection, iEvent);
332  }
333  break;
334  }
335  }
336  }
337  h_PreL1 -> Fill(preL1);
338  h_PreHLT -> Fill(preHLT);
339  h_Pre -> Fill(prescale);
340  h_PreL1wt -> Fill(preL1, mybxlumi);
341  h_PreHLTwt -> Fill(preHLT, mybxlumi);
342 
343  // check if trigger names in (new) config
344  // std::cout << "changed " <<changed << std::endl;
345  if (changed) {
346  changed = false;
347  if ((verbosity/10)%10 > 1) {
348  std::cout<<"New trigger menu found !!!" << std::endl;
349  const unsigned int n(hltConfig_.size());
350  for (unsigned itrig=0; itrig<triggerNames_.size(); itrig++) {
351  unsigned int triggerindx = hltConfig_.triggerIndex(triggerNames_[itrig]);
352  std::cout << triggerNames_[itrig] << " " << triggerindx << " ";
353  if (triggerindx >= n)
354  std::cout << "does not exist in the current menu" << std::endl;
355  else
356  std::cout << "exists" << std::endl;
357  }
358  }
359  }
360  }
361 }
362 
364  t_PixcandP ->clear();
365  t_PixcandPt ->clear();
366  t_PixcandEta ->clear();
367  t_PixcandPhi ->clear();
368  for(unsigned int i=0; i< t_PixcandMaxP->size(); i++)
369  t_PixcandMaxP[i].clear();
370  t_PixcandMaxP ->clear();
371  t_PixTrkcandP ->clear();
372  t_PixTrkcandPt ->clear();
373  t_PixTrkcandEta ->clear();
374  t_PixTrkcandPhi ->clear();
375  t_PixTrkcandMaxP ->clear();
376  t_PixTrkcandselTk ->clear();
377 }
378 
380  t_NFcandP ->clear();
381  t_NFcandPt ->clear();
382  t_NFcandEta ->clear();
383  t_NFcandPhi ->clear();
384  t_NFcandEmip ->clear();
385  t_NFTrkcandP ->clear();
386  t_NFTrkcandPt ->clear();
387  t_NFTrkcandEta ->clear();
388  t_NFTrkcandPhi ->clear();
389  t_NFTrkcandEmip ->clear();
390  t_NFTrkMinDR ->clear();
391  t_NFTrkMinDP1 ->clear();
392  t_NFTrkselTkFlag ->clear();
393  t_NFTrkqltyFlag ->clear();
394  t_NFTrkMissFlag ->clear();
395  t_NFTrkPVFlag ->clear();
396  t_NFTrkPropFlag ->clear();
397  t_NFTrkChgIsoFlag->clear();
398  t_NFTrkNeuIsoFlag->clear();
399  t_NFTrkMipFlag ->clear();
400  t_ECone ->clear();
401 }
402 
404  std::vector<double> &PixMaxP, double &TrkMaxP,
405  bool &selTk) {
406  t_PixcandP ->push_back(Pixcand.r());
407  t_PixcandPt ->push_back(Pixcand.pt());
408  t_PixcandEta ->push_back(Pixcand.eta());
409  t_PixcandPhi ->push_back(Pixcand.phi());
410  t_PixcandMaxP ->push_back(PixMaxP);
411  t_PixTrkcandP ->push_back(Trkcand.r());
412  t_PixTrkcandPt ->push_back(Trkcand.pt());
413  t_PixTrkcandEta ->push_back(Trkcand.eta());
414  t_PixTrkcandPhi ->push_back(Trkcand.phi());
415  t_PixTrkcandMaxP ->push_back(TrkMaxP);
416  t_PixTrkcandselTk ->push_back(selTk);
417 
418 }
419 
421  math::XYZTLorentzVector &Trkcand,
422  double &EmipNFcand, double &EmipTrkcand,
423  double &mindR, double &mindP1,
424  std::vector<bool> &Flags, double hCone) {
425  t_NFcandP ->push_back(NFcand.r());
426  t_NFcandPt ->push_back(NFcand.pt());
427  t_NFcandEta ->push_back(NFcand.eta());
428  t_NFcandPhi ->push_back(NFcand.phi());
429  t_NFcandEmip ->push_back(EmipNFcand);
430  t_NFTrkcandP ->push_back(Trkcand.r());
431  t_NFTrkcandPt ->push_back(Trkcand.pt());
432  t_NFTrkcandEta ->push_back(Trkcand.eta());
433  t_NFTrkcandPhi ->push_back(Trkcand.phi());
434  t_NFTrkcandEmip ->push_back(EmipTrkcand);
435  t_NFTrkMinDR ->push_back(mindR);
436  t_NFTrkMinDP1 ->push_back(mindP1);
437  t_NFTrkselTkFlag ->push_back(Flags[0]);
438  t_NFTrkqltyFlag ->push_back(Flags[1]);
439  t_NFTrkMissFlag ->push_back(Flags[2]);
440  t_NFTrkPVFlag ->push_back(Flags[3]);
441  t_NFTrkPropFlag ->push_back(Flags[4]);
442  t_NFTrkChgIsoFlag->push_back(Flags[5]);
443  t_NFTrkNeuIsoFlag->push_back(Flags[6]);
444  t_NFTrkMipFlag ->push_back(Flags[7]);
445  t_ECone ->push_back(hCone);
446 }
447 
449  char hname[100], htit[100];
450  std::string levels[20] = {"L1", "L2", "L3",
451  "Reco", "RecoMatch", "RecoNoMatch",
452  "L2Match", "L2NoMatch", "L3Match", "L3NoMatch",
453  "HLTTrk", "HLTGoodTrk", "HLTIsoTrk", "HLTMip", "HLTSelect",
454  "nonHLTTrk", "nonHLTGoodTrk", "nonHLTIsoTrk", "nonHLTMip", "nonHLTSelect"};
455  if(doTiming) {
456  TimingTree = fs->make<TTree>("TimingTree", "TimingTree");
457  t_timeL2Prod = new std::vector<double>();
458  t_nPixCand = new std::vector<int>();
459  t_nPixSeed = new std::vector<int>();
460 
461  TimingTree->Branch("t_timeL2Prod", "std::vector<double>", &t_timeL2Prod);
462  TimingTree->Branch("t_nPixCand", "std::vector<int>", &t_nPixCand);
463  TimingTree->Branch("t_nPixSeed", "std::vector<int>", &t_nPixSeed);
464  }
465  if(doTrkResTree) {
466  TrkResTree = fs->make<TTree>("TrkRestree", "TrkRestree");
467  t_TrkhCone = new std::vector<double>();
468  t_TrkP = new std::vector<double>();
469  t_TrkselTkFlag = new std::vector<bool>();
470  t_TrkqltyFlag = new std::vector<bool>();
471  t_TrkMissFlag = new std::vector<bool>();
472  t_TrkPVFlag = new std::vector<bool>();
473  t_TrkNuIsolFlag= new std::vector<bool>();
474 
475  TrkResTree->Branch("t_TrkhCone", "std::vector<double>", &t_TrkhCone);
476  TrkResTree->Branch("t_TrkP", "std::vector<double>", &t_TrkP);
477  TrkResTree->Branch("t_TrkselTkFlag", "std::vector<bool>", &t_TrkselTkFlag);
478  TrkResTree->Branch("t_TrkqltyFlag", "std::vector<bool>", &t_TrkqltyFlag);
479  TrkResTree->Branch("t_TrkMissFlag", "std::vector<bool>", &t_TrkMissFlag);
480  TrkResTree->Branch("t_TrkPVFlag", "std::vector<bool>", &t_TrkPVFlag);
481  TrkResTree->Branch("t_TrkNuIsolFlag","std::vector<bool>", &t_TrkNuIsolFlag);
482  }
483  if(doChgIsolTree) {
484  ChgIsolnTree = fs->make<TTree>("ChgIsolnTree", "ChgIsolntree");
485  t_PixcandP = new std::vector<double>();
486  t_PixcandPt = new std::vector<double>();
487  t_PixcandEta = new std::vector<double>();
488  t_PixcandPhi = new std::vector<double>();
489  t_PixcandMaxP = new std::vector<std::vector<double> >();
490  t_PixTrkcandP = new std::vector<double>();
491  t_PixTrkcandPt = new std::vector<double>();
492  t_PixTrkcandEta = new std::vector<double>();
493  t_PixTrkcandPhi = new std::vector<double>();
494  t_PixTrkcandMaxP = new std::vector<double>();
495  t_PixTrkcandselTk = new std::vector<bool>();
496 
497  ChgIsolnTree->Branch("t_PixcandP", "std::vector<double>", &t_PixcandP);
498  ChgIsolnTree->Branch("t_PixcandPt", "std::vector<double>", &t_PixcandPt);
499  ChgIsolnTree->Branch("t_PixcandEta", "std::vector<double>", &t_PixcandEta);
500  ChgIsolnTree->Branch("t_PixcandPhi", "std::vector<double>", &t_PixcandPhi);
501  ChgIsolnTree->Branch("t_PixcandMaxP", "std::vector<std::vector<double> >", &t_PixcandMaxP);
502  ChgIsolnTree->Branch("t_PixTrkcandP", "std::vector<double>", &t_PixTrkcandP);
503  ChgIsolnTree->Branch("t_PixTrkcandPt", "std::vector<double>", &t_PixTrkcandPt );
504  ChgIsolnTree->Branch("t_PixTrkcandEta", "std::vector<double>", &t_PixTrkcandEta );
505  ChgIsolnTree->Branch("t_PixTrkcandPhi", "std::vector<double>", &t_PixTrkcandPhi );
506  ChgIsolnTree->Branch("t_PixTrkcandMaxP", "std::vector<double>", &t_PixTrkcandMaxP);
507  ChgIsolnTree->Branch("t_PixTrkcandselTk","std::vector<bool>", &t_PixTrkcandselTk);
508  }
509  if(doMipCutTree) {
510  MipCutTree = fs->make<TTree>("MipCutTree", "MipCuttree");
511  t_NFcandP = new std::vector<double>();
512  t_NFcandPt = new std::vector<double>();
513  t_NFcandEta = new std::vector<double>();
514  t_NFcandPhi = new std::vector<double>();
515  t_NFcandEmip = new std::vector<double>();
516  t_NFTrkcandP = new std::vector<double>();
517  t_NFTrkcandPt = new std::vector<double>();
518  t_NFTrkcandEta = new std::vector<double>();
519  t_NFTrkcandPhi = new std::vector<double>();
520  t_NFTrkcandEmip = new std::vector<double>();
521  t_NFTrkMinDR = new std::vector<double>();
522  t_NFTrkMinDP1 = new std::vector<double>();
523  t_NFTrkselTkFlag = new std::vector<bool>();
524  t_NFTrkqltyFlag = new std::vector<bool>();
525  t_NFTrkMissFlag = new std::vector<bool>();
526  t_NFTrkPVFlag = new std::vector<bool>();
527  t_NFTrkPropFlag = new std::vector<bool>();
528  t_NFTrkChgIsoFlag= new std::vector<bool>();
529  t_NFTrkNeuIsoFlag= new std::vector<bool>();
530  t_NFTrkMipFlag = new std::vector<bool>();
531  t_ECone = new std::vector<double>();
532 
533  MipCutTree->Branch("t_NFcandP", "std::vector<double>", &t_NFcandP);
534  MipCutTree->Branch("t_NFcandPt", "std::vector<double>", &t_NFcandPt);
535  MipCutTree->Branch("t_NFcandEta", "std::vector<double>", &t_NFcandEta);
536  MipCutTree->Branch("t_NFcandPhi", "std::vector<double>", &t_NFcandPhi);
537  MipCutTree->Branch("t_NFcandEmip", "std::vector<double>", &t_NFcandEmip);
538  MipCutTree->Branch("t_NFTrkcandP", "std::vector<double>", &t_NFTrkcandP);
539  MipCutTree->Branch("t_NFTrkcandPt", "std::vector<double>", &t_NFTrkcandPt );
540  MipCutTree->Branch("t_NFTrkcandEta", "std::vector<double>", &t_NFTrkcandEta );
541  MipCutTree->Branch("t_NFTrkcandPhi", "std::vector<double>", &t_NFTrkcandPhi );
542  MipCutTree->Branch("t_NFTrkcandEmip", "std::vector<double>", &t_NFTrkcandEmip);
543  MipCutTree->Branch("t_NFTrkMinDR", "std::vector<double>", &t_NFTrkMinDR);
544  MipCutTree->Branch("t_NFTrkMinDP1", "std::vector<double>", &t_NFTrkMinDP1);
545  MipCutTree->Branch("t_NFTrkselTkFlag", "std::vector<bool>", &t_NFTrkselTkFlag);
546  MipCutTree->Branch("t_NFTrkqltyFlag", "std::vector<bool>", &t_NFTrkqltyFlag);
547  MipCutTree->Branch("t_NFTrkMissFlag", "std::vector<bool>", &t_NFTrkMissFlag);
548  MipCutTree->Branch("t_NFTrkPVFlag", "std::vector<bool>", &t_NFTrkPVFlag);
549  MipCutTree->Branch("t_NFTrkPropFlag", "std::vector<bool>", &t_NFTrkPropFlag);
550  MipCutTree->Branch("t_NFTrkChgIsoFlag","std::vector<bool>", &t_NFTrkChgIsoFlag);
551  MipCutTree->Branch("t_NFTrkNeuIsoFlag","std::vector<bool>", &t_NFTrkNeuIsoFlag);
552  MipCutTree->Branch("t_NFTrkMipFlag", "std::vector<bool>", &t_NFTrkMipFlag);
553  MipCutTree->Branch("t_ECone", "std::vector<double>", &t_ECone);
554  }
555  h_Filters = fs->make<TH1I>("h_Filters", "Filter Accepts", 10, 0, 10);
556  std::string FilterNames[7] = {"hltL1sL1SingleJet68", "hltIsolPixelTrackL2FilterHE", "ecalIsolPixelTrackFilterHE", "hltIsolPixelTrackL3FilterHE",
557  "hltIsolPixelTrackL2FilterHB", "ecalIsolPixelTrackFilterHB", "hltIsolPixelTrackL3FilterHB"};
558  for(int i=0; i<7; i++) {
559  h_Filters->GetXaxis()->SetBinLabel(i+1, FilterNames[i].c_str());
560  }
561 
562  h_nHLT = fs->make<TH1I>("h_nHLT" , "Size of trigger Names", 1000, 1, 1000);
563  h_HLT = fs->make<TH1I>("h_HLT" , "HLT accept", 3, -1, 2);
564  h_PreL1 = fs->make<TH1I>("h_PreL1", "L1 Prescale", 500, 0, 500);
565  h_PreHLT = fs->make<TH1I>("h_PreHLT", "HLT Prescale", 50, 0, 50);
566  h_Pre = fs->make<TH1I>("h_Pre", "Prescale", 3000, 0, 3000);
567 
568  h_PreL1wt = fs->make<TH1D>("h_PreL1wt", "Weighted L1 Prescale", 500, 0, 500);
569  h_PreHLTwt = fs->make<TH1D>("h_PreHLTwt", "Weighted HLT Prescale", 500, 0, 100);
570  h_L1ObjEnergy = fs->make<TH1D>("h_L1ObjEnergy", "Energy of L1Object", 500, 0.0, 500.0);
571 
572  h_EnIn = fs->make<TH1D>("h_EnInEcal", "EnergyIn Ecal", 200, 0.0, 20.0);
573  h_EnOut = fs->make<TH1D>("h_EnOutEcal", "EnergyOut Ecal", 200, 0.0, 20.0);
574  h_MipEnMatch = fs->make<TH2D>("h_MipEnMatch", "MipEn: HLT level vs Reco Level (Matched)", 200, 0.0, 20.0, 200, 0.0, 20.0);
575  h_MipEnNoMatch = fs->make<TH2D>("h_MipEnNoMatch", "MipEn: HLT level vs Reco Level (No Match Found)", 200, 0.0, 20.0, 200, 0.0, 20.0);
576 
577  if (doL2L3) {
578  h_nL3Objs = fs->make<TH1I>("h_nL3Objs", "Number of L3 objects", 10, 0, 10);
579 
580  std::string pairs[9] = {"L2L3", "L2L3Match", "L2L3NoMatch", "L3Reco", "L3RecoMatch", "L3RecoNoMatch", "NewFilterReco", "NewFilterRecoMatch", "NewFilterRecoNoMatch"};
581  for (int ipair=0; ipair<9; ipair++) {
582  sprintf(hname, "h_dEta%s", pairs[ipair].c_str());
583  sprintf(htit, "#Delta#eta for %s", pairs[ipair].c_str());
584  h_dEta[ipair] = fs->make<TH1D>(hname, htit, 200, -10.0, 10.0);
585  h_dEta[ipair]->GetXaxis()->SetTitle("d#eta");
586 
587  sprintf(hname, "h_dPhi%s", pairs[ipair].c_str());
588  sprintf(htit, "#Delta#phi for %s", pairs[ipair].c_str());
589  h_dPhi[ipair] = fs->make<TH1D>(hname, htit, 140, -7.0, 7.0);
590  h_dPhi[ipair]->GetXaxis()->SetTitle("d#phi");
591 
592  sprintf(hname, "h_dPt%s", pairs[ipair].c_str());
593  sprintf(htit, "#Delta dp_{T} for %s objects", pairs[ipair].c_str());
594  h_dPt[ipair] = fs->make<TH1D>(hname, htit, 400, -200.0, 200.0);
595  h_dPt[ipair]->GetXaxis()->SetTitle("dp_{T} (GeV)");
596 
597  sprintf(hname, "h_dP%s", pairs[ipair].c_str());
598  sprintf(htit, "#Delta p for %s objects", pairs[ipair].c_str());
599  h_dP[ipair] = fs->make<TH1D>(hname, htit, 400, -200.0, 200.0);
600  h_dP[ipair]->GetXaxis()->SetTitle("dP (GeV)");
601 
602  sprintf(hname, "h_dinvPt%s", pairs[ipair].c_str());
603  sprintf(htit, "#Delta (1/p_{T}) for %s objects", pairs[ipair].c_str());
604  h_dinvPt[ipair] = fs->make<TH1D>(hname, htit, 500, -0.4, 0.1);
605  h_dinvPt[ipair]->GetXaxis()->SetTitle("d(1/p_{T})");
606  sprintf(hname, "h_mindR%s", pairs[ipair].c_str());
607  sprintf(htit, "min(#Delta R) for %s objects", pairs[ipair].c_str());
608  h_mindR[ipair] = fs->make<TH1D>(hname, htit, 500, 0.0, 1.0);
609  h_mindR[ipair]->GetXaxis()->SetTitle("dR");
610  }
611 
612  for (int lvl=0; lvl<2; lvl++) {
613  sprintf(hname, "h_dEtaL1%s", levels[lvl+1].c_str());
614  sprintf(htit, "#Delta#eta for L1 and %s objects", levels[lvl+1].c_str());
615  h_dEtaL1[lvl] = fs->make<TH1D>(hname, htit, 400, -10.0, 10.0);
616 
617  sprintf(hname, "h_dPhiL1%s", levels[lvl+1].c_str());
618  sprintf(htit, "#Delta#phi for L1 and %s objects", levels[lvl+1].c_str());
619  h_dPhiL1[lvl] = fs->make<TH1D>(hname, htit, 280, -7.0, 7.0);
620 
621  sprintf(hname, "h_dRL1%s", levels[lvl+1].c_str());
622  sprintf(htit, "#Delta R for L1 and %s objects", levels[lvl+1].c_str());
623  h_dRL1[lvl] = fs->make<TH1D>(hname, htit, 100, 0.0, 10.0);
624  }
625  }
626 
627  int levmin = (doL2L3 ? 0 : 10);
628  for (int ilevel=levmin; ilevel<20; ilevel++) {
629  sprintf(hname, "h_p%s", levels[ilevel].c_str());
630  sprintf(htit, "p for %s objects", levels[ilevel].c_str());
631  h_p[ilevel] = fs->make<TH1D>(hname, htit, 100, 0.0, 500.0);
632  h_p[ilevel]->GetXaxis()->SetTitle("p (GeV)");
633 
634  sprintf(hname, "h_pt%s", levels[ilevel].c_str());
635  sprintf(htit, "p_{T} for %s objects", levels[ilevel].c_str());
636  h_pt[ilevel] = fs->make<TH1D>(hname, htit, 100, 0.0, 500.0);
637  h_pt[ilevel]->GetXaxis()->SetTitle("p_{T} (GeV)");
638 
639  sprintf(hname, "h_eta%s", levels[ilevel].c_str());
640  sprintf(htit, "#eta for %s objects", levels[ilevel].c_str());
641  h_eta[ilevel] = fs->make<TH1D>(hname, htit, 100, -5.0, 5.0);
642  h_eta[ilevel]->GetXaxis()->SetTitle("#eta");
643 
644  sprintf(hname, "h_phi%s", levels[ilevel].c_str());
645  sprintf(htit, "#phi for %s objects", levels[ilevel].c_str());
646  h_phi[ilevel] = fs->make<TH1D>(hname, htit, 70, -3.5, 3.50);
647  h_phi[ilevel]->GetXaxis()->SetTitle("#phi");
648  }
649 
650  std::string cuts[2] = {"HLTMatched", "HLTNotMatched"};
651  std::string cuts2[2] = {"All", "Away from L1"};
652  for (int icut=0; icut<2; icut++) {
653  sprintf(hname, "h_eMip%s", cuts[icut].c_str());
654  sprintf(htit, "eMip for %s tracks", cuts[icut].c_str());
655  h_eMip[icut] =fs->make<TH1D>(hname, htit, 200, 0.0, 10.0);
656  h_eMip[icut]->GetXaxis()->SetTitle("E_{Mip} (GeV)");
657 
658  sprintf(hname, "h_eMaxNearP%s", cuts[icut].c_str());
659  sprintf(htit, "eMaxNearP for %s tracks", cuts[icut].c_str());
660  h_eMaxNearP[icut]=fs->make<TH1D>(hname, htit, 240, -2.0, 10.0);
661  h_eMaxNearP[icut]->GetXaxis()->SetTitle("E_{MaxNearP} (GeV)");
662 
663  sprintf(hname, "h_eNeutIso%s", cuts[icut].c_str());
664  sprintf(htit, "eNeutIso for %s ", cuts[icut].c_str());
665  h_eNeutIso[icut] =fs->make<TH1D>(hname, htit, 200, 0.0, 10.0);
666  h_eNeutIso[icut]->GetXaxis()->SetTitle("E_{NeutIso} (GeV)");
667 
668  for (int kcut=0; kcut<2; ++kcut) {
669  for (int lim=0; lim<5; ++lim) {
670  sprintf(hname, "h_etaCalibTracks%sCut%dLim%d", cuts[icut].c_str(), kcut, lim);
671  sprintf(htit, "#eta for %s isolated MIP tracks (%4.1f < p < %5.1f Gev/c %s)", cuts[icut].c_str(), pLimits[lim], pLimits[lim+1], cuts2[kcut].c_str());
672  h_etaCalibTracks[lim][icut][kcut]=fs->make<TH1D>(hname, htit, 60, -30.0, 30.0);
673  h_etaCalibTracks[lim][icut][kcut]->GetXaxis()->SetTitle("i#eta");
674 
675  sprintf(hname, "h_etaMipTracks%sCut%dLim%d", cuts[icut].c_str(), kcut, lim);
676  sprintf(htit, "#eta for %s charge isolated MIP tracks (%4.1f < p < %5.1f Gev/c %s)", cuts[icut].c_str(), pLimits[lim], pLimits[lim+1], cuts2[kcut].c_str());
677  h_etaMipTracks[lim][icut][kcut]=fs->make<TH1D>(hname, htit, 60, -30.0, 30.0);
678  h_etaMipTracks[lim][icut][kcut]->GetXaxis()->SetTitle("i#eta");
679  }
680  }
681  }
682 
683  std::string ecut1[3] = {"all","HLTMatched","HLTNotMatched"};
684  std::string ecut2[2] = {"without","with"};
685  int etac[48] = {-1,-2,-3,-4,-5,-6,-7,-8,-9,-10,-11,-12,-13,-14,-15,-16,-17,-18,-19,-20,-21,-22,-23,-24,
686  1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24};
687  for (int icut=0; icut<6; icut++) {
688  // int i1 = (icut>3 ? 1 : 0);
689  int i1 = (icut>2 ? 1 : 0);
690  int i2 = icut - i1*3;
691  for (int kcut=0; kcut<48; kcut++) {
692  for (int lim=0; lim<5; ++lim) {
693  sprintf(hname, "h_eta%dEnHcal%s%s%d", etac[kcut], ecut1[i2].c_str(), ecut2[i1].c_str(), lim);
694  sprintf(htit, "HCAL energy for #eta=%d for %s tracks (p=%4.1f:%5.1f Gev) %s neutral isolation", etac[kcut], ecut1[i2].c_str(), pLimits[lim], pLimits[lim+1], ecut2[i1].c_str());
695  h_eHcal[lim][icut][kcut]=fs->make<TH1D>(hname, htit, 750, 0.0, 150.0);
696  h_eHcal[lim][icut][kcut]->GetXaxis()->SetTitle("Energy (GeV)");
697  sprintf(hname, "h_eta%dEnCalo%s%s%d", etac[kcut], ecut1[i2].c_str(), ecut2[i1].c_str(), lim);
698  sprintf(htit, "Calorimter energy for #eta=%d for %s tracks (p=%4.1f:%5.1f Gev) %s neutral isolation", etac[kcut], ecut1[i2].c_str(), pLimits[lim], pLimits[lim+1], ecut2[i1].c_str());
699  h_eCalo[lim][icut][kcut]=fs->make<TH1D>(hname, htit, 750, 0.0, 150.0);
700  h_eCalo[lim][icut][kcut]->GetXaxis()->SetTitle("Energy (GeV)");
701  }
702  }
703  }
704 }
705 
706 // ------------ method called once each job just after ending the event loop ------------
708  unsigned int preL1, preHLT;
709  std::map<unsigned int, unsigned int>::iterator itr;
710  std::map<unsigned int, const std::pair<int, int>>::iterator itrPre;
711  if (verbosity%10 > 0) std::cout << trigNames.size() << "Triggers were run. RunNo vs HLT accepts for";
712  for (unsigned int i=0; i<trigNames.size(); ++i)
713  std::cout << " " << trigNames[i];
714  std::cout << std::endl;
715  unsigned int n = maxRunNo - minRunNo +1;
716  g_Pre = fs->make<TH1I>("h_PrevsRN", "PreScale Vs Run Number", n, minRunNo, maxRunNo);
717  g_PreL1 = fs->make<TH1I>("h_PreL1vsRN", "L1 PreScale Vs Run Number", n, minRunNo, maxRunNo);
718  g_PreHLT = fs->make<TH1I>("h_PreHLTvsRN", "HLT PreScale Vs Run Number", n, minRunNo, maxRunNo);
719  g_Accepts = fs->make<TH1I>("h_HLTAcceptsvsRN", "HLT Accepts Vs Run Number", n, minRunNo, maxRunNo);
720 
721  for (itr=TrigList.begin(), itrPre=TrigPreList.begin(); itr!=TrigList.end(); itr++, itrPre++) {
722  preL1 = (itrPre->second).first;
723  preHLT = (itrPre->second).second;
724  std::cout << itr->first << " " << itr->second << " " << itrPre->first << " " << preL1 << " " << preHLT << std::endl;
725  g_Accepts->Fill(itr->first, itr->second);
726  g_PreL1->Fill(itr->first, preL1);
727  g_PreHLT->Fill(itr->first, preHLT);
728  g_Pre->Fill(itr->first, preL1*preHLT);
729  }
730 }
731 
732 // ------------ method called when starting to processes a run ------------
733 void IsoTrig::beginRun(edm::Run const& iRun, edm::EventSetup const& iSetup) {
734  std::cout << "Run " << iRun.run() << " hltconfig.init "
735  << hltConfig_.init(iRun,iSetup,processName,changed) << std::endl;;
736 }
737 
738 // ------------ method called when ending the processing of a run ------------
739 void IsoTrig::endRun(edm::Run const&, edm::EventSetup const&) {}
740 
741 // ------------ method called when starting to processes a luminosity block ------------
743 // ------------ method called when ending the processing of a luminosity block ------------
745 // ------------ method fills 'descriptions' with the allowed parameters for the module ------------
746 
748 
749  t_TrkselTkFlag->clear();
750  t_TrkqltyFlag->clear();
751  t_TrkMissFlag->clear();
752  t_TrkPVFlag->clear();
753  t_TrkNuIsolFlag->clear();
754  t_TrkhCone->clear();
755  t_TrkP->clear();
756 
757  if (!trkCollection.isValid()) {
758  std::cout << "trkCollection.isValid is false" << std::endl;
759  } else {
760  std::vector<spr::propagatedTrackDirection>::const_iterator trkDetItr;
761  const MagneticField *bField = bFieldH.product();
762  const CaloGeometry* geo = pG.product();
763  std::vector<spr::propagatedTrackDirection> trkCaloDirections1;
764  spr::propagateCALO(trkCollection, geo, bField, theTrackQuality, trkCaloDirections1, ((verbosity/100)%10>2));
765  unsigned int nTracks=0;
766  int nRH_eMipDR=0, nNearTRKs=0;
767  std::vector<bool> selFlags;
768  for (trkDetItr = trkCaloDirections1.begin(); trkDetItr != trkCaloDirections1.end(); trkDetItr++,nTracks++) {
769  double conehmaxNearP = 0, hCone=0, eMipDR=0.0;
770  const reco::Track* pTrack = &(*(trkDetItr->trkItr));
771  if (verbosity%10>0) std::cout << "track no. " << nTracks << " p(): " << pTrack->p() << std::endl;
772  if (pTrack->p() > 20) {
773  math::XYZTLorentzVector v2(pTrack->px(), pTrack->py(),
774  pTrack->pz(), pTrack->p());
776  trkDetItr->pointHCAL, trkDetItr->pointECAL,
777  a_mipR, trkDetItr->directionECAL, nRH_eMipDR);
778  // std::cout << "mip " << eMipDR << std::endl;
779  bool selectTk = spr::goodTrack(pTrack,leadPV,selectionParameters,((verbosity/100)%10>1));
781  oneCutParameters.maxDxyPV = 10;
782  oneCutParameters.maxDzPV = 100;
783  oneCutParameters.maxInMiss = 2;
784  oneCutParameters.maxOutMiss= 2;
785  bool qltyFlag = spr::goodTrack(pTrack,leadPV,oneCutParameters,((verbosity/100)%10>1));
786  oneCutParameters = selectionParameters;
787  oneCutParameters.maxDxyPV = 10;
788  oneCutParameters.maxDzPV = 100;
789  bool qltyMissFlag = spr::goodTrack(pTrack,leadPV,oneCutParameters,((verbosity/100)%10>1));
790  oneCutParameters = selectionParameters;
791  oneCutParameters.maxInMiss = 2;
792  oneCutParameters.maxOutMiss= 2;
793  bool qltyPVFlag = spr::goodTrack(pTrack,leadPV,oneCutParameters,((verbosity/100)%10>1));
794  /*
795  std::cout << "sel " << selectTk << std::endl;
796  std::cout << "ntracks " << nTracks;
797  std::cout << " a_charIsoR " << a_charIsoR;
798  std::cout << " nNearTRKs " << nNearTRKs << std::endl;
799  */
800  conehmaxNearP = spr::chargeIsolationCone(nTracks, trkCaloDirections1, a_charIsoR, nNearTRKs, ((verbosity/100)%10>1));
801  /*
802  std::cout << "coneh " << conehmaxNearP << std::endl;
803  std::cout << "ok " << trkDetItr->okECAL << " " << trkDetItr->okHCAL << std::endl;
804  */
806  trkDetItr->pointHCAL, trkDetItr->pointECAL,
807  a_neutR1, trkDetItr->directionECAL, nRH_eMipDR);
808 // std::cout << "e1 " << e1 << std::endl;
810  trkDetItr->pointHCAL, trkDetItr->pointECAL,
811  a_neutR2, trkDetItr->directionECAL, nRH_eMipDR);
812  double e_inCone = e2 - e1;
813  bool chgIsolFlag = (conehmaxNearP < cutCharge);
814  bool mipFlag = (eMipDR < cutMip);
815  bool neuIsolFlag = (e_inCone < cutNeutral);
816  bool trkpropFlag = ((trkDetItr->okECAL) && (trkDetItr->okHCAL));
817  selFlags.clear();
818  selFlags.push_back(selectTk); selFlags.push_back(qltyFlag);
819  selFlags.push_back(qltyMissFlag); selFlags.push_back(qltyPVFlag);
820  if (verbosity%10>0) std::cout << "emip: " << eMipDR << "<" << cutMip << "(" << mipFlag << ")"
821  << " ; ok: " << trkDetItr->okECAL << "/" << trkDetItr->okHCAL
822  << " ; chgiso: " << conehmaxNearP << "<" << cutCharge << "(" << chgIsolFlag << ")" << std::endl;
823  // std::cout << selectTk << " " << qltyFlag << " " << qltyMissFlag << " " << qltyPVFlag << " " << trkpropFlag << " " << chgIsolFlag << " " << neuIsolFlag << " " << mipFlag << std::endl;
824 
825  if(chgIsolFlag && mipFlag && trkpropFlag) {
826  double distFromHotCell=-99.0;
827  int nRecHitsCone=-99, ietaHotCell=-99, iphiHotCell=-99;
828  GlobalPoint gposHotCell(0.,0.,0.);
829  std::vector<DetId> coneRecHitDetIds;
830  hCone = spr::eCone_hcal(geo, hbhe, trkDetItr->pointHCAL,
831  trkDetItr->pointECAL,
832  a_coneR, trkDetItr->directionHCAL,
833  nRecHitsCone, coneRecHitDetIds,
834  distFromHotCell, ietaHotCell, iphiHotCell,
835  gposHotCell, -1);
836  // push vectors into the Tree
837  t_TrkselTkFlag ->push_back(selFlags[0]);
838  t_TrkqltyFlag ->push_back(selFlags[1]);
839  t_TrkMissFlag ->push_back(selFlags[2]);
840  t_TrkPVFlag ->push_back(selFlags[3]);
841  t_TrkNuIsolFlag->push_back(neuIsolFlag);
842  t_TrkhCone ->push_back(hCone);
843  t_TrkP ->push_back(pTrack->p());
844  }
845  }
846  }
847  if (verbosity%10>0) std::cout << "Filling " << t_TrkP->size() << " tracks in TrkRestree out of " << nTracks << std::endl;
848  }
849  TrkResTree->Fill();
850 }
851 
852 void IsoTrig::studyTiming(const edm::Event& theEvent) {
853  t_timeL2Prod->clear(); t_nPixCand->clear(); t_nPixSeed->clear();
854 
855  double tHB=0.0, tHE=0.0;
856  int nCandHB=pixelTrackRefsHB.size(), nCandHE=pixelTrackRefsHE.size();
857  int nSeedHB=0, nSeedHE=0;
858 
859  if(nCandHE>0) {
860  edm::Handle<reco::VertexCollection> pVertHB, pVertHE;
861  theEvent.getByToken(tok_verthb_,pVertHB);
862  theEvent.getByToken(tok_verthe_,pVertHE);
864  theEvent.getByToken(tok_hlt_, l1trigobj);
865 
866  std::vector< edm::Ref<l1extra::L1JetParticleCollection> > l1tauobjref;
867  std::vector< edm::Ref<l1extra::L1JetParticleCollection> > l1jetobjref;
868  std::vector< edm::Ref<l1extra::L1JetParticleCollection> > l1forjetobjref;
869 
870  l1trigobj->getObjects(trigger::TriggerL1TauJet, l1tauobjref);
871  l1trigobj->getObjects(trigger::TriggerL1CenJet, l1jetobjref);
872  l1trigobj->getObjects(trigger::TriggerL1ForJet, l1forjetobjref);
873 
874  double ptTriggered = -10;
875  double etaTriggered = -100;
876  double phiTriggered = -100;
877  for (unsigned int p=0; p<l1tauobjref.size(); p++) {
878  if (l1tauobjref[p]->pt()>ptTriggered) {
879  ptTriggered = l1tauobjref[p]->pt();
880  phiTriggered = l1tauobjref[p]->phi();
881  etaTriggered = l1tauobjref[p]->eta();
882  }
883  }
884  for (unsigned int p=0; p<l1jetobjref.size(); p++) {
885  if (l1jetobjref[p]->pt()>ptTriggered) {
886  ptTriggered = l1jetobjref[p]->pt();
887  phiTriggered = l1jetobjref[p]->phi();
888  etaTriggered = l1jetobjref[p]->eta();
889  }
890  }
891  for (unsigned int p=0; p<l1forjetobjref.size(); p++) {
892  if (l1forjetobjref[p]->pt()>ptTriggered) {
893  ptTriggered=l1forjetobjref[p]->pt();
894  phiTriggered=l1forjetobjref[p]->phi();
895  etaTriggered=l1forjetobjref[p]->eta();
896  }
897  }
898  for(unsigned iS=0; iS<pixelTrackRefsHE.size(); iS++) {
899  reco::VertexCollection::const_iterator vitSel;
900  double minDZ = 100;
901  bool vtxMatch;
902  for (reco::VertexCollection::const_iterator vit=pVertHE->begin(); vit!=pVertHE->end(); vit++) {
903  if (fabs(pixelTrackRefsHE[iS]->dz(vit->position()))<minDZ) {
904  minDZ = fabs(pixelTrackRefsHE[iS]->dz(vit->position()));
905  vitSel = vit;
906  }
907  }
908  //cut on dYX:
909  if (minDZ!=100&&fabs(pixelTrackRefsHE[iS]->dxy(vitSel->position()))<vtxCutSeed_) vtxMatch=true;
910  if (minDZ==100) vtxMatch=true;
911 
912  //select tracks not matched to triggered L1 jet
913  double R=deltaR(etaTriggered, phiTriggered, pixelTrackRefsHE[iS]->eta(), pixelTrackRefsHE[iS]->phi());
914  if (R>tauUnbiasCone_ && vtxMatch) nSeedHE++;
915  }
916  for(unsigned iS=0; iS<pixelTrackRefsHB.size(); iS++) {
917  reco::VertexCollection::const_iterator vitSel;
918  double minDZ = 100;
919  bool vtxMatch;
920  for (reco::VertexCollection::const_iterator vit=pVertHB->begin(); vit!=pVertHB->end(); vit++) {
921  if (fabs(pixelTrackRefsHB[iS]->dz(vit->position()))<minDZ) {
922  minDZ = fabs(pixelTrackRefsHB[iS]->dz(vit->position()));
923  vitSel = vit;
924  }
925  }
926  //cut on dYX:
927  if (minDZ!=100&&fabs(pixelTrackRefsHB[iS]->dxy(vitSel->position()))<101.0) vtxMatch=true;
928  if (minDZ==100) vtxMatch=true;
929 
930  //select tracks not matched to triggered L1 jet
931  double R=deltaR(etaTriggered, phiTriggered, pixelTrackRefsHB[iS]->eta(), pixelTrackRefsHB[iS]->phi());
932  if (R>1.2 && vtxMatch) nSeedHB++;
933  }
934 
936  tHE = ft->queryModuleTimeByLabel(theEvent.streamID(),"hltIsolPixelTrackProdHE") ;
937  tHB = ft->queryModuleTimeByLabel(theEvent.streamID(),"hltIsolPixelTrackProdHB");
938 
939  if (verbosity%10>0) std::cout << "(HB/HE) time: " << tHB <<"/" << tHE <<
940  " nCand: " << nCandHB << "/" << nCandHE <<
941  "nSeed: " << nSeedHB << "/" << nSeedHE << std::endl;
942  }
943  t_timeL2Prod->push_back(tHB); t_timeL2Prod->push_back(tHE);
944  t_nPixSeed->push_back(nSeedHB); t_nPixSeed->push_back(nSeedHE);
945  t_nPixCand->push_back(nCandHB); t_nPixCand->push_back(nCandHE);
946 
947  TimingTree->Fill();
948 }
951 
953  if (verbosity%10>0) std::cout << "inside studymipcut" << std::endl;
954  if (!trkCollection.isValid()) {
955  std::cout << "trkCollection.isValid is false" << std::endl;
956  } else {
957  std::vector<spr::propagatedTrackDirection>::const_iterator trkDetItr;
958  const MagneticField *bField = bFieldH.product();
959  const CaloGeometry* geo = pG.product();
960  std::vector<spr::propagatedTrackDirection> trkCaloDirections1;
961  spr::propagateCALO(trkCollection, geo, bField, theTrackQuality, trkCaloDirections1, ((verbosity/100)%10>2));
962  if (verbosity%10>0) std::cout << "Number of L2cands:" << L2cands->size()
963  << " to be matched to something out of "
964  << trkCaloDirections1.size() << " reco tracks" << std::endl;
965  for (unsigned int i=0; i<L2cands->size(); i++) {
968  double enIn = candref->energyIn();
969  h_EnIn->Fill(candref->energyIn());
970  h_EnOut->Fill(candref->energyOut());
971  math::XYZTLorentzVector v1(candref->track()->px(),candref->track()->py(),
972  candref->track()->pz(),candref->track()->p());
973  if (verbosity%10>1) std::cout << "HLT Level Candidate eta/phi/pt/enIn:"
974  << candref->track()->eta() << "/" << candref->track()->phi()
975  << "/" << candref->track()->pt() << "/" << candref->energyIn()
976  << std::endl;
977  math::XYZTLorentzVector mindRvec;
978  double mindR=999.9, mindP1=999.9, eMipDR=0.0;
979  std::vector<bool> selFlags;
980  unsigned int nTracks=0;
981  double conehmaxNearP = 0, hCone=0;
982  for (trkDetItr = trkCaloDirections1.begin(); trkDetItr != trkCaloDirections1.end(); trkDetItr++,nTracks++){
983  const reco::Track* pTrack = &(*(trkDetItr->trkItr));
984  math::XYZTLorentzVector v2(pTrack->px(), pTrack->py(),
985  pTrack->pz(), pTrack->p());
986  double dr = dR(v1,v2);
987  double dp1 = std::abs(1./v1.r() - 1./v2.r());
988  // std::cout <<"This recotrack(eta/phi/pt) " << pTrack->eta() << "/"
989  // << pTrack->phi() << "/" << pTrack->pt() << " has dr/dp= "
990  // << dr << "/" << dp << "/" << dp1 << std::endl;
991  if (dr<mindR) {
992  mindR = dr;
993  mindP1= dp1;
994  mindRvec=v2;
995  int nRH_eMipDR=0, nNearTRKs=0;
996 // std::cout << "debug starts here ";
997 // std::cout << geo;
998 // std::cout << "_" << barrelRecHitsHandle;
999 // std::cout << "_" << endcapRecHitsHandle << "_" << std::endl;
1001  trkDetItr->pointHCAL, trkDetItr->pointECAL,
1002  a_mipR, trkDetItr->directionECAL, nRH_eMipDR);
1003  // std::cout << "mip " << eMipDR << std::endl;
1004  bool selectTk = spr::goodTrack(pTrack,leadPV,selectionParameters,((verbosity/100)%10>1));
1006  oneCutParameters.maxDxyPV = 10;
1007  oneCutParameters.maxDzPV = 100;
1008  oneCutParameters.maxInMiss = 2;
1009  oneCutParameters.maxOutMiss= 2;
1010  bool qltyFlag = spr::goodTrack(pTrack,leadPV,oneCutParameters,((verbosity/100)%10>1));
1011  oneCutParameters = selectionParameters;
1012  oneCutParameters.maxDxyPV = 10;
1013  oneCutParameters.maxDzPV = 100;
1014  bool qltyMissFlag = spr::goodTrack(pTrack,leadPV,oneCutParameters,((verbosity/100)%10>1));
1015  oneCutParameters = selectionParameters;
1016  oneCutParameters.maxInMiss = 2;
1017  oneCutParameters.maxOutMiss= 2;
1018  bool qltyPVFlag = spr::goodTrack(pTrack,leadPV,oneCutParameters,((verbosity/100)%10>1));
1019 // std::cout << "sel " << selectTk << std::endl;
1020 // std::cout << "ntracks " << nTracks;
1021 // std::cout << " a_charIsoR " << a_charIsoR;
1022 // std::cout << " nNearTRKs " << nNearTRKs << std::endl;
1023  conehmaxNearP = spr::chargeIsolationCone(nTracks, trkCaloDirections1, a_charIsoR, nNearTRKs, ((verbosity/100)%10>1));
1024  // std::cout << "coneh " << conehmaxNearP << std::endl;
1025  // std::cout << "ok " << trkDetItr->okECAL << " " << trkDetItr->okHCAL << std::endl;
1027  trkDetItr->pointHCAL, trkDetItr->pointECAL,
1028  a_neutR1, trkDetItr->directionECAL, nRH_eMipDR);
1029  // std::cout << "e1 " << e1 << std::endl;
1031  trkDetItr->pointHCAL, trkDetItr->pointECAL,
1032  a_neutR2, trkDetItr->directionECAL, nRH_eMipDR);
1033  double e_inCone = e2 - e1;
1034  bool chgIsolFlag = (conehmaxNearP < cutCharge);
1035  bool mipFlag = (eMipDR < cutMip);
1036  bool neuIsolFlag = (e_inCone < cutNeutral);
1037  bool trkpropFlag = ((trkDetItr->okECAL) && (trkDetItr->okHCAL));
1038  selFlags.clear();
1039  selFlags.push_back(selectTk); selFlags.push_back(qltyFlag);
1040  selFlags.push_back(qltyMissFlag); selFlags.push_back(qltyPVFlag);
1041  selFlags.push_back(trkpropFlag); selFlags.push_back(chgIsolFlag);
1042  selFlags.push_back(neuIsolFlag); selFlags.push_back(mipFlag);
1043  // std::cout << selectTk << " " << qltyFlag << " " << qltyMissFlag << " " << qltyPVFlag << " " << trkpropFlag << " " << chgIsolFlag << " " << neuIsolFlag << " " << mipFlag << std::endl;
1044  double distFromHotCell=-99.0;
1045  int nRecHitsCone=-99, ietaHotCell=-99, iphiHotCell=-99;
1046  GlobalPoint gposHotCell(0.,0.,0.);
1047  std::vector<DetId> coneRecHitDetIds;
1048  hCone = spr::eCone_hcal(geo, hbhe, trkDetItr->pointHCAL,
1049  trkDetItr->pointECAL,
1050  a_coneR, trkDetItr->directionHCAL,
1051  nRecHitsCone, coneRecHitDetIds,
1052  distFromHotCell, ietaHotCell, iphiHotCell,
1053  gposHotCell, -1);
1054  }
1055  }
1056  pushMipCutTreeVecs(v1, mindRvec, enIn, eMipDR, mindR, mindP1, selFlags, hCone);
1057  fillDifferences(6, v1, mindRvec, (verbosity%10 >0));
1058  if (mindR>=0.05) {
1059  fillDifferences(8, v1, mindRvec, (verbosity%10 >0));
1060  // std::cout << "No Matching found between reco tracks and HLT Level Pixel Tracks: " << mindR << ":" << mindP1 << std::endl;
1061  h_MipEnNoMatch->Fill(candref->energyIn(), eMipDR);
1062  } else {
1063  fillDifferences(7, v1, mindRvec, (verbosity%10 >0));
1064  // std::cout << "Reco Track seems to be macthing eta/phi/pt/eMipDR:"
1065  // << mindRvec.eta() << "/" << mindRvec.phi() << "/"
1066  // << mindRvec.pt() << "/" << eMipDR << " Match " << mindR
1067  // << ":" << mindP1 << std::endl;
1068  h_MipEnMatch->Fill(candref->energyIn(), eMipDR);
1069  }
1070  }
1071  }
1072  MipCutTree->Fill();
1073 }
1075  std::vector<reco::TrackCollection::const_iterator>& goodTks) {
1076 
1077  if (verbosity%10 > 0) std::cout << "Inside StudyTrigger" << std::endl;
1079  for (int j=0; j<3; j++) {
1080  for (unsigned int k=0; k<vec[j].size(); k++) {
1081  if (verbosity%10 > 0) std::cout << "vec[" << j << "][" << k << "] pt " << vec[j][k].pt() << " eta " << vec[j][k].eta() << " phi " << vec[j][k].phi() << std::endl;
1082  fillHist(j, vec[j][k]);
1083  }
1084  }
1085 
1086  double deta, dphi, dr;
1088  for (int lvl=1; lvl<3; lvl++) {
1089  for (unsigned int i=0; i<vec[lvl].size(); i++) {
1090  deta = dEta(vec[0][0],vec[lvl][i]);
1091  dphi = dPhi(vec[0][0],vec[lvl][i]);
1092  dr = dR(vec[0][0],vec[lvl][i]);
1093  if (verbosity%10 > 1) std::cout << "lvl " <<lvl << " i " << i << " deta " << deta << " dphi " << dphi << " dR " << dr << std::endl;
1094  h_dEtaL1[lvl-1] -> Fill(deta);
1095  h_dPhiL1[lvl-1] -> Fill(dphi);
1096  h_dRL1[lvl-1] -> Fill(std::sqrt(dr));
1097  }
1098  }
1099 
1100  math::XYZTLorentzVector mindRvec;
1101  double mindR;
1102  for (unsigned int k=0; k<vec[2].size(); ++k) {
1104  mindR=999.9;
1105  if (verbosity%10 > 1) std::cout << "L3obj: pt " << vec[2][k].pt() << " eta " << vec[2][k].eta() << " phi " << vec[2][k].phi() << std::endl;
1106  for (unsigned int j=0; j<vec[1].size(); j++) {
1107  dr = dR(vec[2][k],vec[1][j]);
1108  if (dr<mindR) {
1109  mindR=dr;
1110  mindRvec=vec[1][j];
1111  }
1112  }
1113  fillDifferences(0, vec[2][k], mindRvec, (verbosity%10 >0));
1114  if (mindR < 0.03) {
1115  fillDifferences(1, vec[2][k], mindRvec, (verbosity%10 >0));
1116  fillHist(6, mindRvec);
1117  fillHist(8, vec[2][k]);
1118  } else {
1119  fillDifferences(2, vec[2][k], mindRvec, (verbosity%10 >0));
1120  fillHist(7, mindRvec);
1121  fillHist(9, vec[2][k]);
1122  }
1123 
1125  mindR=999.9;
1126  if (verbosity%10 > 0)
1127  std::cout << "Now Matching L3 track with reco: L3 Track (eta, phi) "
1128  << vec[2][k].eta() << ":" << vec[2][k].phi() << " L2 Track "
1129  << mindRvec.eta() << ":" << mindRvec.phi() << " dR "
1130  << mindR << std::endl;
1131  reco::TrackCollection::const_iterator goodTk = trkCollection->end();
1132  if (trkCollection.isValid()) {
1133  double mindP(9999.9), mindP1(9999.9);
1134  reco::TrackCollection::const_iterator trkItr;
1135  for (trkItr=trkCollection->begin();
1136  trkItr!=trkCollection->end(); trkItr++) {
1137  math::XYZTLorentzVector v4(trkItr->px(), trkItr->py(),
1138  trkItr->pz(), trkItr->p());
1139  double deltaR = dR(v4, vec[2][k]);
1140  double dp = std::abs(v4.r()/vec[2][k].r()-1.0);
1141  double dp1 = std::abs(1./v4.r()-1./vec[2][k].r());
1142  if (deltaR<mindR) {
1143  mindR = deltaR;
1144  mindP = dp;
1145  mindP1 = dp1;
1146  mindRvec = v4;
1147  goodTk = trkItr;
1148  }
1149  if ((verbosity/10)%10>1 && deltaR<1.0) {
1150  std::cout << "reco track: pt " << v4.pt() << " eta " << v4.eta()
1151  << " phi " << v4.phi() << " DR " << deltaR
1152  << std::endl;
1153  }
1154  }
1155  if (verbosity%10 > 0)
1156  std::cout << "Now Matching at Reco level in step 1 DR: " << mindR
1157  << ":" << mindP << ":" << mindP1 << " eta:phi "
1158  << mindRvec.eta() << ":" << mindRvec.phi() << std::endl;
1159  if (mindR < 0.03 && mindP > 0.1) {
1160  for (trkItr=trkCollection->begin();
1161  trkItr!=trkCollection->end(); trkItr++) {
1162  math::XYZTLorentzVector v4(trkItr->px(), trkItr->py(),
1163  trkItr->pz(), trkItr->p());
1164  double deltaR = dR(v4, vec[2][k]);
1165  double dp = std::abs(v4.r()/vec[2][k].r()-1.0);
1166  double dp1 = std::abs(1./v4.r()-1./vec[2][k].r());
1167  if (dp<mindP && deltaR<0.03) {
1168  mindR = deltaR;
1169  mindP = dp;
1170  mindP1 = dp1;
1171  mindRvec = v4;
1172  goodTk = trkItr;
1173  }
1174  }
1175  if (verbosity%10 > 0)
1176  std::cout << "Now Matching at Reco level in step 2 DR: " << mindR
1177  << ":" << mindP << ":" << mindP1 << " eta:phi "
1178  << mindRvec.eta() << ":" << mindRvec.phi() << std::endl;
1179 
1180  }
1181  fillDifferences(3, vec[2][k], mindRvec, (verbosity%10 >0));
1182  fillHist(3, mindRvec);
1183  if(mindR < 0.03) {
1184  fillDifferences(4, vec[2][k], mindRvec, (verbosity%10 >0));
1185  fillHist(4, mindRvec);
1186  } else {
1187  fillDifferences(5, vec[2][k], mindRvec, (verbosity%10 >0));
1188  fillHist(5, mindRvec);
1189  }
1190  if (goodTk != trkCollection->end()) goodTks.push_back(goodTk);
1191  }
1192  }
1193 }
1194 
1196  std::vector<reco::TrackCollection::const_iterator>& goodTks) {
1197 
1198  if (trkCollection.isValid()) {
1199  // get handles to calogeometry and calotopology
1200  const CaloGeometry* geo = pG.product();
1201  const MagneticField *bField = bFieldH.product();
1202  std::vector<spr::propagatedTrackDirection> trkCaloDirections;
1203  spr::propagateCALO(trkCollection, geo, bField, theTrackQuality, trkCaloDirections, ((verbosity/100)%10>2));
1204 
1205  std::vector<spr::propagatedTrackDirection>::const_iterator trkDetItr;
1206 
1207  if ((verbosity/1000)%10 > 1) {
1208  std::cout << "n of barrelRecHitsHandle " << barrelRecHitsHandle->size() << std::endl;
1210  std::cout << "hit : detid(ieta,iphi) " << (EBDetId)(hit->id()) << " time " << hit->time() << " energy " << hit->energy() << std::endl;
1211  }
1212  std::cout << "n of endcapRecHitsHandle " << endcapRecHitsHandle->size() << std::endl;
1214  std::cout << "hit : detid(ieta,iphi) " << (EEDetId)(hit->id()) << " time " << hit->time() << " energy " << hit->energy() << std::endl;
1215  }
1216  std::cout << "n of hbhe " << hbhe->size() << std::endl;
1217  for (HBHERecHitCollection::const_iterator hit = hbhe->begin(); hit != hbhe->end(); ++hit) {
1218  std::cout << "hit : detid(ieta,iphi) " << hit->id() << " time " << hit->time() << " energy " << hit->energy() << std::endl;
1219  }
1220  }
1221  unsigned int nTracks=0, ngoodTk=0, nselTk=0;
1222  int ieta=999;
1223  for (trkDetItr = trkCaloDirections.begin(); trkDetItr != trkCaloDirections.end(); trkDetItr++,nTracks++){
1224  bool l3Track = (std::find(goodTks.begin(),goodTks.end(),trkDetItr->trkItr) != goodTks.end());
1225  const reco::Track* pTrack = &(*(trkDetItr->trkItr));
1226  math::XYZTLorentzVector v4(pTrack->px(), pTrack->py(),
1227  pTrack->pz(), pTrack->p());
1228  bool selectTk = spr::goodTrack(pTrack,leadPV,selectionParameters,((verbosity/100)%10>1));
1229  double eMipDR=9999., e_inCone=0, conehmaxNearP=0, mindR=999.9, hCone=0;
1230  if (trkDetItr->okHCAL) {
1231  HcalDetId detId = (HcalDetId)(trkDetItr->detIdHCAL);
1232  ieta = detId.ieta();
1233  }
1234  for (unsigned k=0; k<vec[0].size(); ++k) {
1235  double deltaR = dR(v4, vec[0][k]);
1236  if (deltaR<mindR) mindR = deltaR;
1237  }
1238  if ((verbosity/100)%10 > 1) std::cout << "Track ECAL " << trkDetItr->okECAL << " HCAL " << trkDetItr->okHCAL << " Flag " << selectTk << std::endl;
1239  if (selectTk && trkDetItr->okECAL && trkDetItr->okHCAL) {
1240  ngoodTk++;
1241  int nRH_eMipDR=0, nNearTRKs=0;
1243  trkDetItr->pointHCAL, trkDetItr->pointECAL,
1244  a_neutR1, trkDetItr->directionECAL, nRH_eMipDR);
1246  trkDetItr->pointHCAL, trkDetItr->pointECAL,
1247  a_neutR2, trkDetItr->directionECAL, nRH_eMipDR);
1249  trkDetItr->pointHCAL, trkDetItr->pointECAL,
1250  a_mipR, trkDetItr->directionECAL, nRH_eMipDR);
1251  conehmaxNearP = spr::chargeIsolationCone(nTracks, trkCaloDirections, a_charIsoR, nNearTRKs, ((verbosity/100)%10>1));
1252  e_inCone = e2 - e1;
1253  double distFromHotCell=-99.0;
1254  int nRecHitsCone=-99, ietaHotCell=-99, iphiHotCell=-99;
1255  GlobalPoint gposHotCell(0.,0.,0.);
1256  std::vector<DetId> coneRecHitDetIds;
1257  hCone = spr::eCone_hcal(geo, hbhe, trkDetItr->pointHCAL,
1258  trkDetItr->pointECAL,
1259  a_coneR, trkDetItr->directionHCAL,
1260  nRecHitsCone, coneRecHitDetIds,
1261  distFromHotCell, ietaHotCell, iphiHotCell,
1262  gposHotCell, -1);
1263  if (eMipDR<1.0) nselTk++;
1264  }
1265  if (l3Track) {
1266  fillHist(10,v4);
1267  if (selectTk) {
1268  fillHist(11,v4);
1269  fillCuts(0, eMipDR, conehmaxNearP, e_inCone, v4, ieta, (mindR>dr_L1));
1270  if (conehmaxNearP < cutCharge) {
1271  fillHist(12,v4);
1272  if (eMipDR < cutMip) {
1273  fillHist(13,v4);
1274  fillEnergy(1, ieta, hCone, eMipDR, v4);
1275  fillEnergy(0, ieta, hCone, eMipDR, v4);
1276  if (e_inCone < cutNeutral) {
1277  fillHist(14, v4);
1278  fillEnergy(4, ieta, hCone, eMipDR, v4);
1279  fillEnergy(3, ieta, hCone, eMipDR, v4);
1280  }
1281  }
1282  }
1283  }
1284  } else if (doL2L3) {
1285  fillHist(15,v4);
1286  if (selectTk) {
1287  fillHist(16,v4);
1288  fillCuts(1, eMipDR, conehmaxNearP, e_inCone, v4, ieta, (mindR>dr_L1));
1289  if (conehmaxNearP < cutCharge) {
1290  fillHist(17,v4);
1291  if (eMipDR < cutMip) {
1292  fillHist(18,v4);
1293  fillEnergy(2, ieta, hCone, eMipDR, v4);
1294  fillEnergy(0, ieta, hCone, eMipDR, v4);
1295  if (e_inCone < cutNeutral) {
1296  fillHist(19, v4);
1297  fillEnergy(5, ieta, hCone, eMipDR, v4);
1298  fillEnergy(3, ieta, hCone, eMipDR, v4);
1299  }
1300  }
1301  }
1302  }
1303  }
1304  }
1305  // std::cout << "Number of tracks selected offline " << nselTk << std::endl;
1306  }
1307 }
1308 
1309 void IsoTrig::chgIsolation(double& etaTriggered, double& phiTriggered,
1310  edm::Handle<reco::TrackCollection>& trkCollection,
1311  const edm::Event& theEvent) {
1313  if (verbosity%10>0) std::cout << "Inside chgIsolation() with eta/phi Triggered: " << etaTriggered << "/" << phiTriggered << std::endl;
1314  std::vector<double> maxP;
1315 
1316  std::vector<spr::propagatedTrackDirection>::const_iterator trkDetItr;
1317  const MagneticField *bField = bFieldH.product();
1318  const CaloGeometry* geo = pG.product();
1319  std::vector<spr::propagatedTrackDirection> trkCaloDirections1;
1320  spr::propagateCALO(trkCollection, geo, bField, theTrackQuality, trkCaloDirections1, ((verbosity/100)%10>2));
1321  if (verbosity%10>0) std::cout << "Propagated TrkCollection" << std::endl;
1322  for (unsigned int k=0; k<pixelIsolationConeSizeAtEC_.size(); ++k)
1323  maxP.push_back(0);
1324  unsigned i = pixelTrackRefsHE.size();
1325  std::vector<std::pair<unsigned int, std::pair<double, double>>> VecSeedsatEC;
1326  //loop to select isolated tracks
1327  for (unsigned iS=0; iS<pixelTrackRefsHE.size(); iS++) {
1328  if (pixelTrackRefsHE[iS]->p() > minPTrackValue_) {
1329  bool vtxMatch = false;
1330  //associate to vertex (in Z)
1331  unsigned int ivSel = recVtxs->size();
1332  double minDZ = 100;
1333  for (unsigned int iv = 0; iv < recVtxs->size(); ++iv) {
1334  if (fabs(pixelTrackRefsHE[iS]->dz((*recVtxs)[iv].position()))<minDZ) {
1335  minDZ = fabs(pixelTrackRefsHE[iS]->dz((*recVtxs)[iv].position()));
1336  ivSel = iv;
1337  }
1338  }
1339  //cut on dYX:
1340  if (ivSel == recVtxs->size()) {
1341  vtxMatch = true;
1342  } else if (fabs(pixelTrackRefsHE[iS]->dxy((*recVtxs)[ivSel].position()))<vtxCutSeed_){
1343  vtxMatch = true;
1344  }
1345  //select tracks not matched to triggered L1 jet
1346  double R = deltaR(etaTriggered, phiTriggered, pixelTrackRefsHE[iS]->eta(),
1347  pixelTrackRefsHE[iS]->phi());
1348  if (R > tauUnbiasCone_ && vtxMatch) {
1349  //propagate seed track to ECAL surface:
1350  std::pair<double,double> seedCooAtEC;
1351  // in case vertex is found:
1352  if (minDZ!=100) seedCooAtEC=GetEtaPhiAtEcal(pixelTrackRefsHE[iS]->eta(), pixelTrackRefsHE[iS]->phi(), pixelTrackRefsHE[iS]->pt(), pixelTrackRefsHE[iS]->charge(), (*recVtxs)[ivSel].z());
1353  //in case vertex is not found:
1354  else seedCooAtEC=GetEtaPhiAtEcal(pixelTrackRefsHE[iS]->eta(), pixelTrackRefsHE[iS]->phi(), pixelTrackRefsHE[iS]->pt(), pixelTrackRefsHE[iS]->charge(), 0);
1355  VecSeedsatEC.push_back(std::make_pair(iS, seedCooAtEC));
1356  }
1357  }
1358  }
1359 // std::cout << " index of the matched L2cand: " << i << " eta/phi/p: " << pixelTrackRefs[i]->eta()
1360 // << "/" << pixelTrackRefs[i]->phi() << "/" << pixelTrackRefs[i]->p() << std::endl;;
1361 // std::cout << " VecSeedsatEC.size() " << VecSeedsatEC.size() << " pixelTrackRefs.size() " << pixelTrackRefs.size() << std::endl;
1362  for (unsigned int l=0; l<VecSeedsatEC.size(); l++) {
1363  unsigned int iSeed = VecSeedsatEC[l].first;
1364  // std::cout << l << "th Pixcandidate iseed:" << iSeed << std::endl;
1365  math::XYZTLorentzVector v1(pixelTrackRefsHE[iSeed]->px(),pixelTrackRefsHE[iSeed]->py(),
1366  pixelTrackRefsHE[iSeed]->pz(),pixelTrackRefsHE[iSeed]->p());
1367 
1368  for (unsigned int j=0; j<VecSeedsatEC.size(); j++) {
1369  unsigned int iSurr = VecSeedsatEC[j].first;
1370  if (iSeed != iSurr) {
1371  //define preliminary cone around seed track impact point from which tracks will be extrapolated:
1372  // edm::Ref<reco::IsolatedPixelTrackCandidateCollection> cand2ref =
1373  // edm::Ref<reco::IsolatedPixelTrackCandidateCollection>(L2cands, iSurr);
1374  if (deltaR(pixelTrackRefsHE[iSeed]->eta(), pixelTrackRefsHE[iSeed]->phi(), pixelTrackRefsHE[iSurr]->eta(),
1375  pixelTrackRefsHE[iSurr]->phi()) < prelimCone_) {
1376  unsigned int ivSel = recVtxs->size();
1377  double minDZ2=100;
1378  for (unsigned int iv = 0; iv < recVtxs->size(); ++iv) {
1379  if (fabs(pixelTrackRefsHE[iSurr]->dz((*recVtxs)[iv].position()))<minDZ2) {
1380  minDZ2 = fabs(pixelTrackRefsHE[iSurr]->dz((*recVtxs)[iv].position()));
1381  ivSel = iv;
1382  }
1383  }
1384  // std::cout << "minDZ2 " << minDZ2 << " dxy " << pixelTrackRefs[iSurr]->dxy((*recVtxs)[ivSel].position()) << " vtxCutIsol_ " << vtxCutIsol_ <<std::endl;
1385  //cut ot dXY:
1386  if (minDZ2==100 || fabs(pixelTrackRefsHE[iSurr]->dxy((*recVtxs)[ivSel].position()))<vtxCutIsol_) {
1387  //calculate distance at ECAL surface and update isolation:
1388  double dist = getDistInCM(VecSeedsatEC[i].second.first, VecSeedsatEC[i].second.second, VecSeedsatEC[j].second.first, VecSeedsatEC[j].second.second);
1389  // std::cout << "pixelIsolationConeSizeAtEC_.size() " <<pixelIsolationConeSizeAtEC_.size() << std::endl;
1390  for (unsigned int k=0; k<pixelIsolationConeSizeAtEC_.size(); ++k) {
1391  if (dist<pixelIsolationConeSizeAtEC_[k]) {
1392  if (pixelTrackRefsHE[iSurr]->p() > maxP[k])
1393  maxP[k] = pixelTrackRefsHE[iSurr]->p();
1394  }
1395  }
1396  }
1397  }
1398  }
1399  }
1400 
1401  double conehmaxNearP = -1; bool selectTk=false;
1402  double mindR=999.9; int nTracks=0;
1403  math::XYZTLorentzVector mindRvec;
1404  for (trkDetItr = trkCaloDirections1.begin(); trkDetItr != trkCaloDirections1.end(); trkDetItr++, nTracks++){
1405  int nNearTRKs=0;
1406  const reco::Track* pTrack = &(*(trkDetItr->trkItr));
1407  math::XYZTLorentzVector v2(pTrack->px(), pTrack->py(),
1408  pTrack->pz(), pTrack->p());
1409  double dr = dR(v1,v2);
1410  // double dp1 = std::abs(1./v1.r() - 1./v2.r());
1411  // std::cout <<"THis recotrack(eta/phi/pt) " << pTrack->eta() << "/"
1412  // << pTrack->phi() << "/" << pTrack->pt() << " has dr/dp= "
1413  // << dr << "/" << dp << "/" << dp1 << std::endl;
1414  if (dr<mindR) {
1415  selectTk = spr::goodTrack(pTrack,leadPV,selectionParameters,((verbosity/100)%10>1));
1416  conehmaxNearP = spr::chargeIsolationCone(nTracks, trkCaloDirections1, a_charIsoR, nNearTRKs, ((verbosity/100)%10>1));
1417  mindR = dr;
1418  mindRvec = v2;
1419  }
1420  }
1421  pushChgIsolnTreeVecs(v1, mindRvec, maxP, conehmaxNearP, selectTk);
1422  // std::cout << " conehmaxNearP " << conehmaxNearP << std::endl;
1423  // for (unsigned int k=0; k<pixelIsolationConeSizeAtEC_.size(); ++k)
1424  // std::cout << "k/maxP[k]: " << k << "/" << maxP[k] << std::endl;
1425  }
1426  ChgIsolnTree->Fill();
1427 }
1428 
1430  h_p[indx]->Fill(vec.r());
1431  h_pt[indx]->Fill(vec.pt());
1432  h_eta[indx]->Fill(vec.eta());
1433  h_phi[indx]->Fill(vec.phi());
1434 }
1435 
1438  double dr = dR(vec1,vec2);
1439  double deta = dEta(vec1, vec2);
1440  double dphi = dPhi(vec1, vec2);
1441  double dpt = dPt(vec1, vec2);
1442  double dp = dP(vec1, vec2);
1443  double dinvpt = dinvPt(vec1, vec2);
1444  h_dEta[indx] ->Fill(deta);
1445  h_dPhi[indx] ->Fill(dphi);
1446  h_dPt[indx] ->Fill(dpt);
1447  h_dP[indx] ->Fill(dp);
1448  h_dinvPt[indx]->Fill(dinvpt);
1449  h_mindR[indx] ->Fill(dr);
1450  if (debug) std::cout << "mindR for index " << indx << " is " << dr << " deta " << deta << " dphi " << dphi << " dpt " << dpt << " dinvpt " << dinvpt <<std::endl;
1451 }
1452 
1453 void IsoTrig::fillCuts(int indx, double eMipDR, double conehmaxNearP, double e_inCone, math::XYZTLorentzVector& vec, int ieta, bool cut) {
1454  h_eMip[indx] ->Fill(eMipDR);
1455  h_eMaxNearP[indx]->Fill(conehmaxNearP);
1456  h_eNeutIso[indx] ->Fill(e_inCone);
1457  if ((conehmaxNearP < cutCharge) && (eMipDR < cutMip)) {
1458  for (int lim=0; lim<5; ++lim) {
1459  if ((vec.r()>pLimits[lim]) && (vec.r()<=pLimits[lim+1])) {
1460  h_etaMipTracks[lim][indx][0]->Fill((double)(ieta));
1461  if (cut) h_etaMipTracks[lim][indx][1]->Fill((double)(ieta));
1462  if (e_inCone < cutNeutral) {
1463  h_etaCalibTracks[lim][indx][0]->Fill((double)(ieta));
1464  if (cut) h_etaCalibTracks[lim][indx][1]->Fill((double)(ieta));
1465  }
1466  }
1467  }
1468  }
1469 }
1470 
1471 void IsoTrig::fillEnergy(int indx, int ieta, double hCone, double eMipDR, math::XYZTLorentzVector& vec) {
1472  int kk(-1);
1473  if (ieta > 0 && ieta < 25) kk = 23 + ieta;
1474  else if (ieta > -25 && ieta < 0) kk = -(ieta + 1);
1475  if (kk >= 0 && eMipDR > 0.01 && hCone > 1.0) {
1476  for (int lim=0; lim<5; ++lim) {
1477  if ((vec.r()>pLimits[lim]) && (vec.r()<=pLimits[lim+1])) {
1478  h_eHcal[lim][indx][kk] ->Fill(hCone);
1479  h_eCalo[lim][indx][kk] ->Fill(hCone+eMipDR);
1480  }
1481  }
1482  }
1483 }
1484 
1486  return (vec1.eta()-vec2.eta());
1487 }
1488 
1490 
1491  double phi1 = vec1.phi();
1492  if (phi1 < 0) phi1 += 2.0*M_PI;
1493  double phi2 = vec2.phi();
1494  if (phi2 < 0) phi2 += 2.0*M_PI;
1495  double dphi = phi1-phi2;
1496  if (dphi > M_PI) dphi -= 2.*M_PI;
1497  else if (dphi < -M_PI) dphi += 2.*M_PI;
1498  return dphi;
1499 }
1500 
1502  double deta = dEta(vec1,vec2);
1503  double dphi = dPhi(vec1,vec2);
1504  return std::sqrt(deta*deta + dphi*dphi);
1505 }
1506 
1508  return (vec1.pt()-vec2.pt());
1509 }
1510 
1512  return (std::abs(vec1.r()-vec2.r()));
1513 }
1514 
1516  return ((1/vec1.pt())-(1/vec2.pt()));
1517 }
1518 
1519 std::pair<double,double> IsoTrig::etaPhiTrigger() {
1520  double eta(0), phi(0), ptmax(0);
1521  for (unsigned int k=0; k<vec[0].size(); ++k) {
1522  if (k == 0) {
1523  eta = vec[0][k].eta();
1524  phi = vec[0][k].phi();
1525  ptmax = vec[0][k].pt();
1526  } else if (vec[0][k].pt() > ptmax) {
1527  eta = vec[0][k].eta();
1528  phi = vec[0][k].phi();
1529  ptmax = vec[0][k].pt();
1530  }
1531  }
1532  return std::pair<double,double>(eta,phi);
1533 }
1534 
1535 std::pair<double,double> IsoTrig::GetEtaPhiAtEcal(double etaIP, double phiIP,
1536  double pT, int charge,
1537  double vtxZ) {
1538 
1539  double deltaPhi=0;
1540  double etaEC = 100;
1541  double phiEC = 100;
1542 
1543  double Rcurv = 9999999;
1544  if (bfVal!=0) Rcurv=pT*33.3*100/(bfVal*10); //r(m)=pT(GeV)*33.3/B(kG)
1545 
1546  double ecDist = zEE_;
1547  double ecRad = rEB_; //radius of ECAL barrel (cm)
1548  double theta = 2*atan(exp(-etaIP));
1549  double zNew = 0;
1550  if (theta>0.5*M_PI) theta=M_PI-theta;
1551  if (fabs(etaIP)<1.479) {
1552  if ((0.5*ecRad/Rcurv)>1) {
1553  etaEC = 10000;
1554  deltaPhi = 0;
1555  } else {
1556  deltaPhi =-charge*asin(0.5*ecRad/Rcurv);
1557  double alpha1 = 2*asin(0.5*ecRad/Rcurv);
1558  double z = ecRad/tan(theta);
1559  if (etaIP>0) zNew = z*(Rcurv*alpha1)/ecRad+vtxZ; //new z-coordinate of track
1560  else zNew =-z*(Rcurv*alpha1)/ecRad+vtxZ; //new z-coordinate of track
1561  double zAbs=fabs(zNew);
1562  if (zAbs<ecDist) {
1563  etaEC = -log(tan(0.5*atan(ecRad/zAbs)));
1564  deltaPhi = -charge*asin(0.5*ecRad/Rcurv);
1565  }
1566  if (zAbs>ecDist) {
1567  zAbs = (fabs(etaIP)/etaIP)*ecDist;
1568  double Zflight = fabs(zAbs-vtxZ);
1569  double alpha = (Zflight*ecRad)/(z*Rcurv);
1570  double Rec = 2*Rcurv*sin(alpha/2);
1571  deltaPhi =-charge*alpha/2;
1572  etaEC =-log(tan(0.5*atan(Rec/ecDist)));
1573  }
1574  }
1575  } else {
1576  zNew = (fabs(etaIP)/etaIP)*ecDist;
1577  double Zflight = fabs(zNew-vtxZ);
1578  double Rvirt = fabs(Zflight*tan(theta));
1579  double Rec = 2*Rcurv*sin(Rvirt/(2*Rcurv));
1580  deltaPhi =-(charge)*(Rvirt/(2*Rcurv));
1581  etaEC =-log(tan(0.5*atan(Rec/ecDist)));
1582  }
1583 
1584  if (zNew<0) etaEC=-etaEC;
1585  phiEC = phiIP+deltaPhi;
1586 
1587  if (phiEC<-M_PI) phiEC += 2*M_PI;
1588  if (phiEC>M_PI) phiEC -= 2*M_PI;
1589 
1590  std::pair<double,double> retVal(etaEC,phiEC);
1591  return retVal;
1592 }
1593 
1594 double IsoTrig::getDistInCM(double eta1,double phi1, double eta2,double phi2) {
1595  double Rec;
1596  double theta1=2*atan(exp(-eta1));
1597  double theta2=2*atan(exp(-eta2));
1598  if (fabs(eta1)<1.479) Rec=rEB_; //radius of ECAL barrel
1599  else if (fabs(eta1)>1.479&&fabs(eta1)<7.0) Rec=tan(theta1)*zEE_; //distance from IP to ECAL endcap
1600  else return 1000;
1601 
1602  //|vect| times tg of acos(scalar product)
1603  double angle=acos((sin(theta1)*sin(theta2)*(sin(phi1)*sin(phi2)+cos(phi1)*cos(phi2))+cos(theta1)*cos(theta2)));
1604  if (angle<0.5*M_PI) return fabs((Rec/sin(theta1))*tan(angle));
1605  else return 1000;
1606 }
1607 
1608 //define this as a plug-in
RunNumber_t run() const
Definition: EventID.h:39
unsigned int size() const
number of trigger paths in trigger table
double prelimCone_
Definition: IsoTrig.h:130
double p() const
momentum vector magnitude
Definition: TrackBase.h:663
T getParameter(std::string const &) const
EventNumber_t event() const
Definition: EventID.h:41
std::vector< bool > * t_NFTrkPVFlag
Definition: IsoTrig.h:209
std::vector< double > * t_PixTrkcandPhi
Definition: IsoTrig.h:190
std::vector< double > * t_PixTrkcandP
Definition: IsoTrig.h:187
T getUntrackedParameter(std::string const &, T const &) const
double tauUnbiasCone_
Definition: IsoTrig.h:130
int i
Definition: DBlmapReader.cc:9
edm::EDGetTokenT< reco::IsolatedPixelTrackCandidateCollection > tok_l2cand_
Definition: IsoTrig.h:152
float alpha
Definition: AMPTWrapper.h:95
virtual edm::TriggerNames const & triggerNames(edm::TriggerResults const &triggerResults) const
Definition: Event.cc:199
double a_coneR
Definition: IsoTrig.h:133
edm::InputTag L1candTag_
Definition: IsoTrig.h:122
std::vector< spr::propagatedTrackID > propagateCALO(edm::Handle< reco::TrackCollection > &trkCollection, const CaloGeometry *geo, const MagneticField *bField, std::string &theTrackQuality, bool debug=false)
std::vector< double > * t_NFTrkMinDR
Definition: IsoTrig.h:204
int id() const
getters
Definition: TriggerObject.h:55
RunNumber_t run() const
Definition: RunBase.h:42
The single EDProduct to be saved for each event (AOD case)
Definition: TriggerEvent.h:25
trigger::size_type sizeFilters() const
Definition: TriggerEvent.h:135
void chgIsolation(double &etaTriggered, double &phiTriggered, edm::Handle< reco::TrackCollection > &trkCollection, const edm::Event &theEvent)
Definition: IsoTrig.cc:1309
virtual void endLuminosityBlock(edm::LuminosityBlock const &, edm::EventSetup const &)
Definition: IsoTrig.cc:744
std::vector< bool > * t_NFTrkMissFlag
Definition: IsoTrig.h:208
std::vector< reco::TrackRef > pixelTrackRefsHE
Definition: IsoTrig.h:155
double dr_L1
Definition: IsoTrig.h:133
edm::EDGetTokenT< LumiDetails > tok_lumi
Definition: IsoTrig.h:137
TH1D * h_phi[20]
Definition: IsoTrig.h:221
std::vector< bool > * t_PixTrkcandselTk
Definition: IsoTrig.h:192
TH1D * h_EnIn
Definition: IsoTrig.h:216
bool doStudyIsol
Definition: IsoTrig.h:124
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:446
TH1D * h_dEta[9]
Definition: IsoTrig.h:223
float phi() const
Definition: TriggerObject.h:58
bool getByToken(EDGetToken token, Handle< PROD > &result) const
TrackQuality
track quality
Definition: TrackBase.h:133
void studyIsolation(edm::Handle< reco::TrackCollection > &, std::vector< reco::TrackCollection::const_iterator > &)
Definition: IsoTrig.cc:1195
virtual GlobalVector inTesla(const GlobalPoint &gp) const =0
Field value ad specified global point, in Tesla.
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
TH1I * h_Pre
Definition: IsoTrig.h:218
std::vector< double > * t_NFTrkcandEta
Definition: IsoTrig.h:201
std::vector< double > * t_NFcandP
Definition: IsoTrig.h:194
edm::Handle< EcalRecHitCollection > barrelRecHitsHandle
Definition: IsoTrig.h:159
TH1D * h_dEtaL1[2]
Definition: IsoTrig.h:222
bool doTiming
Definition: IsoTrig.h:124
std::string theTrackQuality
Definition: IsoTrig.h:132
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
edm::EDGetTokenT< EcalRecHitCollection > tok_EB_
Definition: IsoTrig.h:144
const Keys & filterKeys(trigger::size_type index) const
Definition: TriggerEvent.h:111
TH1I * g_PreHLT
Definition: IsoTrig.h:228
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
math::XYZPoint leadPV
Definition: IsoTrig.h:163
std::vector< EcalRecHit >::const_iterator const_iterator
Geom::Theta< T > theta() const
std::vector< reco::TrackRef > pixelTrackRefsHB
Definition: IsoTrig.h:155
int bunchCrossing() const
Definition: EventBase.h:62
bool changed
Definition: IsoTrig.h:167
edm::LuminosityBlockNumber_t luminosityBlock() const
Definition: EventBase.h:59
std::vector< math::XYZTLorentzVector > vec[3]
Definition: IsoTrig.h:229
std::vector< int > * t_nPixSeed
Definition: IsoTrig.h:172
std::vector< double > * t_TrkhCone
Definition: IsoTrig.h:174
double rEB_
Definition: IsoTrig.h:127
edm::EDGetTokenT< reco::IsolatedPixelTrackCandidateCollection > tok_pixtk_
Definition: IsoTrig.h:150
std::vector< bool > * t_NFTrkqltyFlag
Definition: IsoTrig.h:207
T * make(const Args &...args) const
make new ROOT object
Definition: TFileService.h:64
bool doTrkResTree
Definition: IsoTrig.h:124
std::pair< double, double > GetEtaPhiAtEcal(double etaIP, double phiIP, double pT, int charge, double vtxZ)
Definition: IsoTrig.cc:1535
std::vector< double > * t_PixTrkcandMaxP
Definition: IsoTrig.h:191
float energy() const
Definition: TriggerObject.h:65
double chargeIsolationCone(unsigned int trkIndex, std::vector< spr::propagatedTrackDirection > &trkDirs, double dR, int &nNearTRKs, bool debug=false)
float eta() const
Definition: TriggerObject.h:57
double px() const
x coordinate of momentum vector
Definition: TrackBase.h:675
T eta() const
std::vector< bool > * t_NFTrkselTkFlag
Definition: IsoTrig.h:206
TH1I * h_HLT
Definition: IsoTrig.h:218
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:7
double dP(math::XYZTLorentzVector &, math::XYZTLorentzVector &)
Definition: IsoTrig.cc:1511
std::vector< double > * t_PixcandP
Definition: IsoTrig.h:182
std::vector< std::string > trigNames
Definition: IsoTrig.h:121
std::pair< double, double > etaPhiTrigger()
Definition: IsoTrig.cc:1519
Strings const & triggerNames() const
Definition: TriggerNames.cc:24
TH1D * h_dP[9]
Definition: IsoTrig.h:223
float float float z
virtual void beginJob()
Definition: IsoTrig.cc:448
HLTConfigProvider hltConfig_
Definition: IsoTrig.h:120
double vtxCutIsol_
Definition: IsoTrig.h:129
std::vector< bool > * t_NFTrkMipFlag
Definition: IsoTrig.h:213
double a_neutIsoR
Definition: IsoTrig.h:133
void pushMipCutTreeVecs(math::XYZTLorentzVector &NFcand, math::XYZTLorentzVector &Trkcand, double &EmipNFcand, double &EmipTrkcand, double &mindR, double &mindP1, std::vector< bool > &Flags, double hCone)
Definition: IsoTrig.cc:420
unsigned int triggerIndex(const std::string &triggerName) const
slot position of trigger path in trigger table (0 to size-1)
Single trigger physics object (e.g., an isolated muon)
Definition: TriggerObject.h:22
U second(std::pair< T, U > const &p)
std::vector< double > * t_NFTrkMinDP1
Definition: IsoTrig.h:205
TH1D * h_EnOut
Definition: IsoTrig.h:216
TH1D * h_dPhiL1[2]
Definition: IsoTrig.h:222
TH1D * h_L1ObjEnergy
Definition: IsoTrig.h:220
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:29
edm::InputTag PixcandTag_
Definition: IsoTrig.h:122
void studyMipCut(edm::Handle< reco::TrackCollection > &trkCollection, edm::Handle< reco::IsolatedPixelTrackCandidateCollection > &L2cands)
Definition: IsoTrig.cc:949
double dinvPt(math::XYZTLorentzVector &, math::XYZTLorentzVector &)
Definition: IsoTrig.cc:1515
int maxRunNo
Definition: IsoTrig.h:135
double bfVal
Definition: IsoTrig.h:127
int verbosity
Definition: IsoTrig.h:126
TH1D * h_dPhi[9]
Definition: IsoTrig.h:223
edm::ESHandle< CaloGeometry > pG
Definition: IsoTrig.h:157
int iEvent
Definition: GenABIO.cc:230
T mag() const
Definition: PV3DBase.h:67
void studyTrigger(edm::Handle< reco::TrackCollection > &, std::vector< reco::TrackCollection::const_iterator > &)
Definition: IsoTrig.cc:1074
void fillHist(int, math::XYZTLorentzVector &)
Definition: IsoTrig.cc:1429
void addDefault(ParameterSetDescription const &psetDescription)
std::vector< double > * t_NFTrkcandP
Definition: IsoTrig.h:199
edm::Handle< HBHERecHitCollection > hbhe
Definition: IsoTrig.h:158
std::vector< bool > * t_TrkNuIsolFlag
Definition: IsoTrig.h:180
TH1I * h_nHLT
Definition: IsoTrig.h:218
double a_mipR
Definition: IsoTrig.h:133
edm::EDGetTokenT< reco::BeamSpot > tok_bs_
Definition: IsoTrig.h:143
bool goodTrack(const reco::Track *pTrack, math::XYZPoint leadPV, trackSelectionParameters parameters, bool debug=false)
std::vector< bool > * t_TrkPVFlag
Definition: IsoTrig.h:179
const TriggerObjectCollection & getObjects() const
Definition: TriggerEvent.h:98
void Fill(HcalDetId &id, double val, std::vector< TH2F > &depth)
std::vector< double > * t_ECone
Definition: IsoTrig.h:214
std::vector< double > vec1
Definition: HCALResponse.h:15
T sqrt(T t)
Definition: SSEVec.h:48
void clear(CLHEP::HepGenMatrix &m)
Helper function: Reset all elements of a matrix to 0.
Definition: matutil.cc:167
edm::EDGetTokenT< reco::TrackCollection > tok_genTrack_
Definition: IsoTrig.h:141
double cutNeutral
Definition: IsoTrig.h:134
TTree * ChgIsolnTree
Definition: IsoTrig.h:169
TH1I * h_PreHLT
Definition: IsoTrig.h:218
TH1D * h_dinvPt[9]
Definition: IsoTrig.h:224
void fillEnergy(int, int, double, double, math::XYZTLorentzVector &)
Definition: IsoTrig.cc:1471
double queryModuleTimeByLabel(edm::StreamID, const std::string &) const
double a_charIsoR
Definition: IsoTrig.h:133
edm::Handle< reco::VertexCollection > recVtxs
Definition: IsoTrig.h:162
double zEE_
Definition: IsoTrig.h:127
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
std::vector< bool > * t_NFTrkChgIsoFlag
Definition: IsoTrig.h:211
int ieta() const
get the cell ieta
Definition: HcalDetId.h:36
std::vector< edm::EDGetTokenT< reco::TrackCollection > > tok_pixtks_
Definition: IsoTrig.h:153
std::vector< double > * t_TrkP
Definition: IsoTrig.h:175
double minPTrackValue_
Definition: IsoTrig.h:129
TH1D * h_eta[20]
Definition: IsoTrig.h:221
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
double cutCharge
Definition: IsoTrig.h:134
std::pair< T, T > etaphi(T x, T y, T z)
Definition: FastMath.h:128
TH1I * h_Filters
Definition: IsoTrig.h:218
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
int j
Definition: DBlmapReader.cc:9
std::vector< double > * t_NFcandEta
Definition: IsoTrig.h:196
std::vector< double > * t_timeL2Prod
Definition: IsoTrig.h:170
double eCone_ecal(const CaloGeometry *geo, edm::Handle< T > &barrelhits, edm::Handle< T > &endcaphits, const GlobalPoint &hpoint1, const GlobalPoint &point1, double dR, const GlobalVector &trackMom, int &nRecHits, double ebThr=-100, double eeThr=-100, double tMin=-500, double tMax=500, bool debug=false)
virtual void analyze(const edm::Event &, const edm::EventSetup &)
Definition: IsoTrig.cc:153
TH1D * h_eHcal[5][6][48]
Definition: IsoTrig.h:227
LuminosityBlock const & getLuminosityBlock() const
Definition: Event.h:77
edm::EDGetTokenT< EcalRecHitCollection > tok_EE_
Definition: IsoTrig.h:145
TTree * MipCutTree
Definition: IsoTrig.h:169
bool doL2L3
Definition: IsoTrig.h:124
static std::string const triggerResults
Definition: EdmProvDump.cc:41
IsoTrig(const edm::ParameterSet &)
Definition: IsoTrig.cc:29
virtual void endJob()
Definition: IsoTrig.cc:707
std::vector< double > * t_PixcandPt
Definition: IsoTrig.h:183
bool first
Definition: L1TdeRCT.cc:75
bool isValid() const
Definition: HandleBase.h:76
double a_neutR2
Definition: IsoTrig.h:134
TH1D * h_PreL1wt
Definition: IsoTrig.h:220
std::map< unsigned int, const std::pair< int, int > > TrigPreList
Definition: IsoTrig.h:166
std::vector< double > * t_NFcandEmip
Definition: IsoTrig.h:198
TH1D * h_eMaxNearP[2]
Definition: IsoTrig.h:225
double dPhi(math::XYZTLorentzVector &, math::XYZTLorentzVector &)
Definition: IsoTrig.cc:1489
void studyTiming(const edm::Event &theEvent)
Definition: IsoTrig.cc:852
#define M_PI
void clearChgIsolnTreeVectors()
Definition: IsoTrig.cc:363
std::vector< bool > * t_NFTrkPropFlag
Definition: IsoTrig.h:210
double pz() const
z coordinate of momentum vector
Definition: TrackBase.h:687
std::vector< double > * t_NFTrkcandPhi
Definition: IsoTrig.h:202
const std::vector< std::string > & moduleLabels(unsigned int trigger) const
label(s) of module(s) on a trigger path
unsigned int id
TH2D * h_MipEnMatch
Definition: IsoTrig.h:217
double deltaR(double eta1, double eta2, double phi1, double phi2)
Definition: TreeUtility.cc:17
TH1D * h_p[20]
Definition: IsoTrig.h:221
TH1D * h_eNeutIso[2]
Definition: IsoTrig.h:225
TH1I * h_PreL1
Definition: IsoTrig.h:218
edm::EDGetTokenT< edm::TriggerResults > tok_trigRes
Definition: IsoTrig.h:139
std::vector< TriggerObject > TriggerObjectCollection
collection of trigger physics objects (e.g., all isolated muons)
Definition: TriggerObject.h:81
const edm::InputTag filterTag(trigger::size_type index) const
Definition: TriggerEvent.h:103
int minRunNo
Definition: IsoTrig.h:135
edm::InputTag L2candTag_
Definition: IsoTrig.h:122
std::vector< bool > * t_TrkMissFlag
Definition: IsoTrig.h:178
double getDistInCM(double eta1, double phi1, double eta2, double phi2)
Definition: IsoTrig.cc:1594
edm::EDGetTokenT< reco::VertexCollection > tok_recVtx_
Definition: IsoTrig.h:142
std::vector< int > * t_nPixCand
Definition: IsoTrig.h:171
static TrackQuality qualityByName(const std::string &name)
Definition: TrackBase.cc:108
#define debug
Definition: HDRShower.cc:19
void StudyTrkEbyP(edm::Handle< reco::TrackCollection > &trkCollection)
Definition: IsoTrig.cc:747
std::vector< double > * t_NFTrkcandPt
Definition: IsoTrig.h:200
TH1D * h_eCalo[5][6][48]
Definition: IsoTrig.h:227
T const * product() const
Definition: Handle.h:81
TH1I * g_Pre
Definition: IsoTrig.h:228
double eCone_hcal(const CaloGeometry *geo, edm::Handle< T > &hits, const GlobalPoint &hpoint1, const GlobalPoint &point1, double dR, const GlobalVector &trackMom, int &nRecHits, double hbThr=-100, double heThr=-100, double hfThr=-100, double hoThr=-100, double tMin=-500, double tMax=500, int detOnly=-1, bool debug=false)
virtual void endRun(edm::Run const &, edm::EventSetup const &)
Definition: IsoTrig.cc:739
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
std::vector< size_type > Keys
edm::Handle< reco::BeamSpot > beamSpotH
Definition: IsoTrig.h:161
std::vector< double > * t_PixTrkcandPt
Definition: IsoTrig.h:188
double dEta(math::XYZTLorentzVector &, math::XYZTLorentzVector &)
Definition: IsoTrig.cc:1485
bool init(const edm::Run &iRun, const edm::EventSetup &iSetup, const std::string &processName, bool &changed)
d&#39;tor
const T & get() const
Definition: EventSetup.h:55
T const * product() const
Definition: ESHandle.h:86
edm::EDGetTokenT< reco::VertexCollection > tok_verthb_
Definition: IsoTrig.h:148
std::pair< int, int > prescaleValues(const edm::Event &iEvent, const edm::EventSetup &iSetup, const std::string &trigger) const
Combined L1T (pair.first) and HLT (pair.second) prescales per HLT path.
bool doChgIsolTree
Definition: IsoTrig.h:124
std::vector< double > * t_NFcandPt
Definition: IsoTrig.h:195
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
Definition: IsoTrig.cc:145
spr::trackSelectionParameters selectionParameters
Definition: IsoTrig.h:131
TH1D * h_pt[20]
Definition: IsoTrig.h:221
std::map< unsigned int, unsigned int > TrigList
Definition: IsoTrig.h:165
std::vector< bool > * t_TrkqltyFlag
Definition: IsoTrig.h:177
TH1D * h_etaMipTracks[5][2][2]
Definition: IsoTrig.h:226
std::string const & label() const
Definition: InputTag.h:42
std::vector< std::vector< double > > * t_PixcandMaxP
Definition: IsoTrig.h:186
std::vector< double > * t_NFTrkcandEmip
Definition: IsoTrig.h:203
edm::EDGetTokenT< trigger::TriggerFilterObjectWithRefs > tok_l1cand_
Definition: IsoTrig.h:151
std::vector< edm::InputTag > pixelTracksSources_
Definition: IsoTrig.h:123
edm::EventID id() const
Definition: EventBase.h:56
double dR(math::XYZTLorentzVector &, math::XYZTLorentzVector &)
Definition: IsoTrig.cc:1501
TH1D * h_mindR[9]
Definition: IsoTrig.h:224
double a_neutR1
Definition: IsoTrig.h:134
void clearMipCutTreeVectors()
Definition: IsoTrig.cc:379
std::vector< double > * t_PixcandEta
Definition: IsoTrig.h:184
edm::EDGetTokenT< HBHERecHitCollection > tok_hbhe_
Definition: IsoTrig.h:146
static int position[264][3]
Definition: ReadPGInfo.cc:509
reco::TrackBase::TrackQuality minQuality
double vtxCutSeed_
Definition: IsoTrig.h:129
StreamID streamID() const
Definition: Event.h:72
void fillDifferences(int, math::XYZTLorentzVector &, math::XYZTLorentzVector &, bool)
Definition: IsoTrig.cc:1436
TH1D * h_dPt[9]
Definition: IsoTrig.h:223
TH2D * h_MipEnNoMatch
Definition: IsoTrig.h:217
double dPt(math::XYZTLorentzVector &, math::XYZTLorentzVector &)
Definition: IsoTrig.cc:1507
std::vector< bool > * t_NFTrkNeuIsoFlag
Definition: IsoTrig.h:212
edm::Handle< EcalRecHitCollection > endcapRecHitsHandle
Definition: IsoTrig.h:160
TTree * TrkResTree
Definition: IsoTrig.h:169
tuple cout
Definition: gather_cfg.py:121
std::vector< double > * t_PixcandPhi
Definition: IsoTrig.h:185
std::vector< double > * t_PixTrkcandEta
Definition: IsoTrig.h:189
bool doMipCutTree
Definition: IsoTrig.h:124
edm::Service< TFileService > fs
Definition: IsoTrig.h:168
TTree * TimingTree
Definition: IsoTrig.h:169
double pLimits[6]
Definition: IsoTrig.h:136
volatile std::atomic< bool > shutdown_flag false
TH1I * g_PreL1
Definition: IsoTrig.h:228
Definition: DDAxes.h:10
edm::ESHandle< MagneticField > bFieldH
Definition: IsoTrig.h:156
float mass() const
Definition: TriggerObject.h:59
edm::EDGetTokenT< reco::VertexCollection > tok_verthe_
Definition: IsoTrig.h:148
list pairs
sort tag files by run number
TH1I * h_nL3Objs
Definition: IsoTrig.h:218
edm::EDGetTokenT< trigger::TriggerFilterObjectWithRefs > tok_hlt_
Definition: IsoTrig.h:147
~IsoTrig()
Definition: IsoTrig.cc:139
std::vector< vec1 > vec2
Definition: HCALResponse.h:16
TH1D * h_dRL1[2]
Definition: IsoTrig.h:222
TH1D * h_etaCalibTracks[5][2][2]
Definition: IsoTrig.h:226
std::vector< bool > * t_TrkselTkFlag
Definition: IsoTrig.h:176
double cutMip
Definition: IsoTrig.h:134
std::string processName
Definition: IsoTrig.h:132
tuple size
Write out results.
void fillCuts(int, double, double, double, math::XYZTLorentzVector &, int, bool)
Definition: IsoTrig.cc:1453
double py() const
y coordinate of momentum vector
Definition: TrackBase.h:681
Definition: Run.h:41
virtual void beginLuminosityBlock(edm::LuminosityBlock const &, edm::EventSetup const &)
Definition: IsoTrig.cc:742
TH1D * h_PreHLTwt
Definition: IsoTrig.h:220
edm::EDGetTokenT< trigger::TriggerEvent > tok_trigEvt
Definition: IsoTrig.h:138
std::vector< double > pixelIsolationConeSizeAtEC_
Definition: IsoTrig.h:128
virtual void beginRun(edm::Run const &, edm::EventSetup const &)
Definition: IsoTrig.cc:733
void pushChgIsolnTreeVecs(math::XYZTLorentzVector &Pixcand, math::XYZTLorentzVector &Trkcand, std::vector< double > &PixMaxP, double &TrkMaxP, bool &selTk)
Definition: IsoTrig.cc:403
TH1I * g_Accepts
Definition: IsoTrig.h:228
T angle(T x1, T y1, T z1, T x2, T y2, T z2)
Definition: angle.h:11
TH1D * h_eMip[2]
Definition: IsoTrig.h:224
Definition: DDAxes.h:10
std::vector< double > * t_NFcandPhi
Definition: IsoTrig.h:197