CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
IsolatedTracksNxN.cc
Go to the documentation of this file.
1 // -*- C++ -*
2 //
3 // Package: IsolatedTracksNxN
4 // Class: IsolatedTracksNxN
5 //
13 //
14 // Original Author: Seema Sharma
15 // Created: Mon Aug 10 15:30:40 CST 2009
16 //
17 //
18 
21 
34 
36 
37  //now do what ever initialization is needed
38  doMC = iConfig.getUntrackedParameter<bool> ("DoMC", false);
39  writeAllTracks = iConfig.getUntrackedParameter<bool> ("WriteAllTracks", false );
40  myverbose_ = iConfig.getUntrackedParameter<int> ("Verbosity", 5 );
41  pvTracksPtMin_ = iConfig.getUntrackedParameter<double>("PVTracksPtMin", 1.0 );
42  debugTrks_ = iConfig.getUntrackedParameter<int> ("DebugTracks" );
43  printTrkHitPattern_ = iConfig.getUntrackedParameter<bool> ("PrintTrkHitPattern" );
44  minTrackP_ = iConfig.getUntrackedParameter<double>("minTrackP", 1.0 );
45  maxTrackEta_ = iConfig.getUntrackedParameter<double>("maxTrackEta", 5.0 );
46  debugL1Info_ = iConfig.getUntrackedParameter<bool> ("DebugL1Info", false );
47  L1TriggerAlgoInfo_ = iConfig.getUntrackedParameter<bool> ("L1TriggerAlgoInfo");
48  L1extraTauJetSource_ = iConfig.getParameter<edm::InputTag> ("L1extraTauJetSource" );
49  L1extraCenJetSource_ = iConfig.getParameter<edm::InputTag> ("L1extraCenJetSource" );
50  L1extraFwdJetSource_ = iConfig.getParameter<edm::InputTag> ("L1extraFwdJetSource" );
51  L1extraMuonSource_ = iConfig.getParameter<edm::InputTag> ("L1extraMuonSource" );
52  L1extraIsoEmSource_ = iConfig.getParameter<edm::InputTag> ("L1extraIsoEmSource" );
53  L1extraNonIsoEmSource_ = iConfig.getParameter<edm::InputTag> ("L1extraNonIsoEmSource" );
54  L1GTReadoutRcdSource_ = iConfig.getParameter<edm::InputTag> ("L1GTReadoutRcdSource" );
55  L1GTObjectMapRcdSource_= iConfig.getParameter<edm::InputTag> ("L1GTObjectMapRcdSource");
56  JetSrc_ = iConfig.getParameter<edm::InputTag> ("JetSource");
57  JetExtender_ = iConfig.getParameter<edm::InputTag> ("JetExtender");
58  HBHERecHitSource_ = iConfig.getParameter<edm::InputTag> ("HBHERecHitSource");
59  tMinE_ = iConfig.getUntrackedParameter<double>("TimeMinCutECAL", -500.);
60  tMaxE_ = iConfig.getUntrackedParameter<double>("TimeMaxCutECAL", 500.);
61  tMinH_ = iConfig.getUntrackedParameter<double>("TimeMinCutHCAL", -500.);
62  tMaxH_ = iConfig.getUntrackedParameter<double>("TimeMaxCutHCAL", 500.);
63  nbad = 0;
64 
65  // define tokens for access
66  tok_L1extTauJet_ = consumes<l1extra::L1JetParticleCollection>(L1extraTauJetSource_);
67  tok_L1extCenJet_ = consumes<l1extra::L1JetParticleCollection>(L1extraCenJetSource_);
68  tok_L1extFwdJet_ = consumes<l1extra::L1JetParticleCollection>(L1extraFwdJetSource_);
69  tok_L1extMu_ = consumes<l1extra::L1MuonParticleCollection>(L1extraMuonSource_);
70  tok_L1extIsoEm_ = consumes<l1extra::L1EmParticleCollection>(L1extraIsoEmSource_);
71  tok_L1extNoIsoEm_ = consumes<l1extra::L1EmParticleCollection>(L1extraNonIsoEmSource_);
72  tok_jets_ = consumes<reco::CaloJetCollection>(JetSrc_);
73  tok_hbhe_ = consumes<HBHERecHitCollection>(HBHERecHitSource_);
74 
75 
76  tok_genTrack_ = consumes<reco::TrackCollection>(edm::InputTag("generalTracks"));
77  tok_recVtx_ = consumes<reco::VertexCollection>(edm::InputTag("offlinePrimaryVertices"));
78  tok_bs_ = consumes<reco::BeamSpot>(edm::InputTag("offlineBeamSpot"));
79  tok_EB_ = consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit","EcalRecHitsEB"));
80  tok_EE_ = consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit","EcalRecHitsEE"));
81  tok_simTk_ = consumes<edm::SimTrackContainer>(edm::InputTag("g4SimHits"));
82  tok_simVtx_ = consumes<edm::SimVertexContainer>(edm::InputTag("g4SimHits"));
83  tok_caloEB_ = consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", "EcalHitsEB"));
84  tok_caloEE_ = consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", "EcalHitsEE"));
85  tok_caloHH_ = consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", "HcalHits"));
86 
87  if(myverbose_>=0) {
88  std::cout <<"Parameters read from config file \n"
89  <<" doMC " << doMC
90  <<"\t myverbose_ " << myverbose_
91  <<"\t minTrackP_ " << minTrackP_
92  << "\t maxTrackEta_ " << maxTrackEta_
93  << "\t tMinE_ " << tMinE_
94  << "\t tMaxE_ " << tMaxE_
95  << "\t tMinH_ " << tMinH_
96  << "\t tMaxH_ " << tMaxH_
97  << "\n debugL1Info_ " << debugL1Info_
98  << "\t L1TriggerAlgoInfo_ " << L1TriggerAlgoInfo_
99  << "\t L1extraTauJetSource_ " << L1extraTauJetSource_
100  << "\t L1extraCenJetSource_ " << L1extraCenJetSource_
101  << "\t L1extraFwdJetSource_ " << L1extraFwdJetSource_
102  << "\t L1extraMuonSource_ " << L1extraMuonSource_
103  << "\t L1extraIsoEmSource_ " << L1extraIsoEmSource_
104  << "\t L1extraNonIsoEmSource_ " << L1extraNonIsoEmSource_
105  << "\t L1GTReadoutRcdSource_ " << L1GTReadoutRcdSource_
106  << "\t L1GTObjectMapRcdSource_ " << L1GTObjectMapRcdSource_
107  << "\t JetSrc_ " << JetSrc_
108  << "\t JetExtender_ " << JetExtender_
109  << "\t HBHERecHitSource_ " << HBHERecHitSource_
110  << "\n" << std::endl;
111  }
112 
113  initL1 = false;
114 
115 }
116 
118 }
119 
121 
122  bool haveIsoTrack=false;
123 
125  iSetup.get<IdealMagneticFieldRecord>().get(bFieldH);
126  bField = bFieldH.product();
127 
129 
130  t_RunNo = iEvent.id().run();
131  t_EvtNo = iEvent.id().event();
132  t_Lumi = iEvent.luminosityBlock();
133  t_Bunch = iEvent.bunchCrossing();
134 
135  nEventProc++;
136 
138  iEvent.getByToken(tok_genTrack_, trkCollection);
139  reco::TrackCollection::const_iterator trkItr;
140  if(debugTrks_>1){
141  std::cout << "Track Collection: " << std::endl;
142  std::cout << "Number of Tracks " << trkCollection->size() << std::endl;
143  }
144  std::string theTrackQuality = "highPurity";
145  reco::TrackBase::TrackQuality trackQuality_=reco::TrackBase::qualityByName(theTrackQuality);
146 
147  //===================== save L1 Trigger information =======================
148  if( L1TriggerAlgoInfo_ ) {
149 
150  bool useL1EventSetup = true;
151  bool useL1GtTriggerMenuLite = true;
152 
153  m_l1GtUtils.getL1GtRunCache(iEvent, iSetup, useL1EventSetup, useL1GtTriggerMenuLite);
154 
155  int iErrorCode = -1;
156  int l1ConfCode = -1;
157  const bool l1Conf = m_l1GtUtils.availableL1Configuration(iErrorCode, l1ConfCode);
158  if( !l1Conf) {
159  std::cout << "\nL1 configuration code:" << l1ConfCode
160  << "\nNo valid L1 trigger configuration available."
161  << "\nSee text above for error code interpretation"
162  << "\nNo return here, in order to test each method, protected against configuration error."
163  << std::endl;
164  }
165 
166  const L1GtTriggerMenu* m_l1GtMenu = m_l1GtUtils.ptrL1TriggerMenuEventSetup(iErrorCode);
167  const AlgorithmMap& algorithmMap = m_l1GtMenu->gtAlgorithmMap();
168  const std::string& menuName = m_l1GtMenu->gtTriggerMenuName();
169 
170  if (!initL1) {
171  initL1=true;
172  std::cout << "menuName " << menuName << std::endl;
173  for (CItAlgo itAlgo = algorithmMap.begin(); itAlgo != algorithmMap.end(); itAlgo++) {
174  std::string algName = itAlgo->first;
175  int algBitNumber = ( itAlgo->second ).algoBitNumber();
176  l1AlgoMap.insert( std::pair<std::pair<unsigned int,std::string>,int>( std::pair<unsigned int,std::string>(algBitNumber, algName) , 0) ) ;
177  }
178  std::map< std::pair<unsigned int,std::string>, int>::iterator itr;
179  for(itr=l1AlgoMap.begin(); itr!=l1AlgoMap.end(); itr++) {
180  std::cout << " ********** " << (itr->first).first <<" "<<(itr->first).second <<" "<<itr->second << std::endl;
181  }
182  }
183 
184  std::vector<int> algbits;
185  for (CItAlgo itAlgo = algorithmMap.begin(); itAlgo != algorithmMap.end(); itAlgo++) {
186  std::string algName = itAlgo->first;
187  int algBitNumber = ( itAlgo->second ).algoBitNumber();
188  //bool algResultBeforeMask = m_l1GtUtils.decisionBeforeMask(iEvent, itAlgo->first, iErrorCode);
189  //bool algResultAfterMask = m_l1GtUtils.decisionAfterMask (iEvent, itAlgo->first, iErrorCode);
190  //int triggerMask = m_l1GtUtils.triggerMask (iEvent, itAlgo->first, iErrorCode);
191  bool decision = m_l1GtUtils.decision (iEvent, itAlgo->first, iErrorCode);
192  int preScale = m_l1GtUtils.prescaleFactor (iEvent, itAlgo->first, iErrorCode);
193 
194  // save the algo names which fired
195  if( decision ){
196  l1AlgoMap[std::pair<unsigned int,std::string>(algBitNumber, algName)] += 1;
197  h_L1AlgoNames->Fill( algBitNumber );
198  t_L1AlgoNames->push_back(itAlgo->first);
199  t_L1PreScale->push_back(preScale);
200  t_L1Decision[algBitNumber] = 1;
201  algbits.push_back(algBitNumber);
202  }
203  }
204 
205  if(debugL1Info_) {
206  for(unsigned int ii=0; ii<t_L1AlgoNames->size(); ii++){
207  std::cout<<ii<<" "<<(*t_L1AlgoNames)[ii]<<" "<<(*t_L1PreScale)[ii]<<" "<<algbits[ii]<<std::endl;
208  }
209  std::cout << "L1Decision: ";
210  for (int i=0; i<128; ++i) {
211  std::cout << " " << i << ":" << t_L1Decision[i];
212  if ((i+1)%32 == 0) std::cout << std::endl;
213  }
214  }
215 
216  // L1Taus
218  iEvent.getByToken(tok_L1extTauJet_,l1TauHandle);
219  l1extra::L1JetParticleCollection::const_iterator itr;
220  int iL1Obj=0;
221  for(itr = l1TauHandle->begin(),iL1Obj=0; itr != l1TauHandle->end(); ++itr,iL1Obj++) {
222  if(iL1Obj<1) {
223  t_L1TauJetPt ->push_back( itr->pt() );
224  t_L1TauJetEta ->push_back( itr->eta() );
225  t_L1TauJetPhi ->push_back( itr->phi() );
226  }
227  if(debugL1Info_) {
228  std::cout << "tauJ p/pt " << itr->momentum() << " " << itr->pt()
229  << " eta/phi " << itr->eta() << " " << itr->phi()
230  << std::endl;
231  }
232  }
233 
234  // L1 Central Jets
236  iEvent.getByToken(tok_L1extCenJet_,l1CenJetHandle);
237  for( itr = l1CenJetHandle->begin(),iL1Obj=0; itr != l1CenJetHandle->end(); ++itr,iL1Obj++ ) {
238  if(iL1Obj<1) {
239  t_L1CenJetPt ->push_back( itr->pt() );
240  t_L1CenJetEta ->push_back( itr->eta() );
241  t_L1CenJetPhi ->push_back( itr->phi() );
242  }
243  if(debugL1Info_) {
244  std::cout << "cenJ p/pt " << itr->momentum() << " " << itr->pt()
245  << " eta/phi " << itr->eta() << " " << itr->phi()
246  << std::endl;
247  }
248  }
249 
250  // L1 Forward Jets
252  iEvent.getByToken(tok_L1extFwdJet_,l1FwdJetHandle);
253  for( itr = l1FwdJetHandle->begin(),iL1Obj=0; itr != l1FwdJetHandle->end(); ++itr,iL1Obj++ ) {
254  if(iL1Obj<1) {
255  t_L1FwdJetPt ->push_back( itr->pt() );
256  t_L1FwdJetEta ->push_back( itr->eta() );
257  t_L1FwdJetPhi ->push_back( itr->phi() );
258  }
259  if(debugL1Info_) {
260  std::cout << "fwdJ p/pt " << itr->momentum() << " " << itr->pt()
261  << " eta/phi " << itr->eta() << " " << itr->phi()
262  << std::endl;
263  }
264  }
265 
266  // L1 Isolated EM onjects
267  l1extra::L1EmParticleCollection::const_iterator itrEm;
269  iEvent.getByToken(tok_L1extIsoEm_, l1IsoEmHandle);
270  for( itrEm = l1IsoEmHandle->begin(),iL1Obj=0; itrEm != l1IsoEmHandle->end(); ++itrEm,iL1Obj++ ) {
271  if(iL1Obj<1) {
272  t_L1IsoEMPt ->push_back( itrEm->pt() );
273  t_L1IsoEMEta ->push_back( itrEm->eta() );
274  t_L1IsoEMPhi ->push_back( itrEm->phi() );
275  }
276  if(debugL1Info_) {
277  std::cout << "isoEm p/pt " << itrEm->momentum() << " " << itrEm->pt()
278  << " eta/phi " << itrEm->eta() << " " << itrEm->phi()
279  << std::endl;
280  }
281  }
282 
283  // L1 Non-Isolated EM onjects
285  iEvent.getByToken(tok_L1extNoIsoEm_, l1NonIsoEmHandle);
286  for( itrEm = l1NonIsoEmHandle->begin(),iL1Obj=0; itrEm != l1NonIsoEmHandle->end(); ++itrEm,iL1Obj++ ) {
287  if(iL1Obj<1) {
288  t_L1NonIsoEMPt ->push_back( itrEm->pt() );
289  t_L1NonIsoEMEta ->push_back( itrEm->eta() );
290  t_L1NonIsoEMPhi ->push_back( itrEm->phi() );
291  }
292  if(debugL1Info_) {
293  std::cout << "nonIsoEm p/pt " << itrEm->momentum() << " " << itrEm->pt()
294  << " eta/phi " << itrEm->eta() << " " << itrEm->phi()
295  << std::endl;
296  }
297  }
298 
299  // L1 Muons
300  l1extra::L1MuonParticleCollection::const_iterator itrMu;
302  iEvent.getByToken(tok_L1extMu_, l1MuHandle);
303  for( itrMu = l1MuHandle->begin(),iL1Obj=0; itrMu != l1MuHandle->end(); ++itrMu,iL1Obj++ ) {
304  if(iL1Obj<1) {
305  t_L1MuonPt ->push_back( itrMu->pt() );
306  t_L1MuonEta ->push_back( itrMu->eta() );
307  t_L1MuonPhi ->push_back( itrMu->phi() );
308  }
309  if(debugL1Info_) {
310  std::cout << "l1muon p/pt " << itrMu->momentum() << " " << itrMu->pt()
311  << " eta/phi " << itrMu->eta() << " " << itrMu->phi()
312  << std::endl;
313  }
314  }
315  }
316 
317  //============== store the information about all the Non-Fake vertices ===============
318 
320  iEvent.getByToken(tok_recVtx_,recVtxs);
321 
322  std::vector<reco::Track> svTracks;
323  math::XYZPoint leadPV(0,0,0);
324  double sumPtMax = -1.0;
325  for(unsigned int ind=0;ind<recVtxs->size();ind++) {
326 
327  if ( !((*recVtxs)[ind].isFake()) ) {
328  double vtxTrkSumPt=0.0, vtxTrkSumPtWt=0.0;
329  int vtxTrkNWt =0;
330  double vtxTrkSumPtHP=0.0, vtxTrkSumPtHPWt=0.0;
331  int vtxTrkNHP =0, vtxTrkNHPWt =0;
332 
333  reco::Vertex::trackRef_iterator vtxTrack = (*recVtxs)[ind].tracks_begin();
334 
335  for (vtxTrack = (*recVtxs)[ind].tracks_begin(); vtxTrack!=(*recVtxs)[ind].tracks_end(); vtxTrack++) {
336 
337  if((*vtxTrack)->pt()<pvTracksPtMin_) continue;
338 
339  bool trkQuality = (*vtxTrack)->quality(trackQuality_);
340 
341  vtxTrkSumPt += (*vtxTrack)->pt();
342 
343  vtxTrkSumPt += (*vtxTrack)->pt();
344  if( trkQuality ) {
345  vtxTrkSumPtHP += (*vtxTrack)->pt();
346  vtxTrkNHP++;
347  }
348 
349  double weight = (*recVtxs)[ind].trackWeight(*vtxTrack);
350  h_PVTracksWt ->Fill(weight);
351  if( weight>0.5) {
352  vtxTrkSumPtWt += (*vtxTrack)->pt();
353  vtxTrkNWt++;
354  if( trkQuality ) {
355  vtxTrkSumPtHPWt += (*vtxTrack)->pt();
356  vtxTrkNHPWt++;
357  }
358  }
359  }
360 
361  if(vtxTrkSumPt>sumPtMax) {
362  sumPtMax = vtxTrkSumPt;
363  leadPV = math::XYZPoint( (*recVtxs)[ind].x(),(*recVtxs)[ind].y(), (*recVtxs)[ind].z() );
364  }
365 
366  t_PVx ->push_back( (*recVtxs)[ind].x() );
367  t_PVy ->push_back( (*recVtxs)[ind].y() );
368  t_PVz ->push_back( (*recVtxs)[ind].z() );
369  t_PVisValid ->push_back( (*recVtxs)[ind].isValid() );
370  t_PVNTracks ->push_back( (*recVtxs)[ind].tracksSize() );
371  t_PVndof ->push_back( (*recVtxs)[ind].ndof() );
372  t_PVNTracksWt ->push_back( vtxTrkNWt );
373  t_PVTracksSumPt ->push_back( vtxTrkSumPt );
374  t_PVTracksSumPtWt->push_back( vtxTrkSumPtWt );
375 
376  t_PVNTracksHP ->push_back( vtxTrkNHP );
377  t_PVNTracksHPWt ->push_back( vtxTrkNHPWt );
378  t_PVTracksSumPtHP ->push_back( vtxTrkSumPtHP );
379  t_PVTracksSumPtHPWt->push_back( vtxTrkSumPtHPWt );
380 
381  if(myverbose_==4) {
382  std::cout<<"PV "<<ind<<" isValid "<<(*recVtxs)[ind].isValid()<<" isFake "<<(*recVtxs)[ind].isFake()
383  <<" hasRefittedTracks() "<<ind<<" "<<(*recVtxs)[ind].hasRefittedTracks()
384  <<" refittedTrksSize "<<(*recVtxs)[ind].refittedTracks().size()
385  <<" tracksSize() "<<(*recVtxs)[ind].tracksSize()<<" sumPt "<<vtxTrkSumPt
386  <<std::endl;
387  }
388  } // if vtx is not Fake
389  } // loop over PVs
390  //===================================================================================
391 
392  // Get the beamspot
393  edm::Handle<reco::BeamSpot> beamSpotH;
394  iEvent.getByToken(tok_bs_, beamSpotH);
395  math::XYZPoint bspot;
396  bspot = ( beamSpotH.isValid() ) ? beamSpotH->position() : math::XYZPoint(0, 0, 0);
397 
398  //=====================================================================
399 
401  iEvent.getByToken(tok_jets_,jets);
402  // edm::Handle<reco::JetExtendedAssociation::Container> jetExtender;
403  // iEvent.getByLabel(JetExtender_,jetExtender);
404 
405  for(unsigned int ijet=0;ijet<(*jets).size();ijet++) {
406  t_jetPt ->push_back( (*jets)[ijet].pt() );
407  t_jetEta ->push_back( (*jets)[ijet].eta() );
408  t_jetPhi ->push_back( (*jets)[ijet].phi() );
409  t_nTrksJetVtx ->push_back( -1.0);
410  t_nTrksJetCalo->push_back( -1.0);
411  }
412 
413  //===================================================================================
414 
415  // get handles to calogeometry and calotopology
417  iSetup.get<CaloGeometryRecord>().get(pG);
418  const CaloGeometry* geo = pG.product();
419 
420  edm::ESHandle<CaloTopology> theCaloTopology;
421  iSetup.get<CaloTopologyRecord>().get(theCaloTopology);
422  const CaloTopology *caloTopology = theCaloTopology.product();
423 
425  iSetup.get<IdealGeometryRecord>().get(htopo);
426  const HcalTopology* theHBHETopology = htopo.product();
427 
428  edm::Handle<EcalRecHitCollection> barrelRecHitsHandle;
429  edm::Handle<EcalRecHitCollection> endcapRecHitsHandle;
430  iEvent.getByToken(tok_EB_,barrelRecHitsHandle);
431  iEvent.getByToken(tok_EE_,endcapRecHitsHandle);
432 
433  // Retrieve the good/bad ECAL channels from the DB
435  iSetup.get<EcalChannelStatusRcd>().get(ecalChStatus);
436  const EcalChannelStatus* theEcalChStatus = ecalChStatus.product();
437 
438  // Retrieve trigger tower map
440  iSetup.get<IdealGeometryRecord>().get(hTtmap);
441  const EcalTrigTowerConstituentsMap& ttMap = *hTtmap;
442 
444  iEvent.getByToken(tok_hbhe_, hbhe);
445  if (!hbhe.isValid()) {
446  nbad++;
447  if (nbad < 10) std::cout << "No HBHE rechit collection\n";
448  return;
449  }
450  const HBHERecHitCollection Hithbhe = *(hbhe.product());
451 
452  //get Handles to SimTracks and SimHits
454  if (doMC) iEvent.getByToken(tok_simTk_,SimTk);
455  edm::SimTrackContainer::const_iterator simTrkItr;
456 
458  if (doMC) iEvent.getByToken(tok_simVtx_,SimVtx);
459 
460  //get Handles to PCaloHitContainers of eb/ee/hbhe
462  if (doMC) iEvent.getByToken(tok_caloEB_, pcaloeb);
463 
465  if (doMC) iEvent.getByToken(tok_caloEE_, pcaloee);
466 
468  if (doMC) iEvent.getByToken(tok_caloHH_, pcalohh);
469 
470  //associates tracker rechits/simhits to a track
471  TrackerHitAssociator* associate=0;
472  if (doMC) associate = new TrackerHitAssociator(iEvent);
473 
474  //===================================================================================
475 
476  h_nTracks->Fill(trkCollection->size());
477 
478  int nTracks = 0;
479 
480  t_nTracks = trkCollection->size();
481 
482  // get the list of DetIds closest to the impact point of track on surface calorimeters
483  std::vector<spr::propagatedTrackID> trkCaloDets;
484  spr::propagateCALO(trkCollection, geo, bField, theTrackQuality, trkCaloDets, false);
485  std::vector<spr::propagatedTrackID>::const_iterator trkDetItr;
486 
487  if(myverbose_>2) {
488  for(trkDetItr = trkCaloDets.begin(); trkDetItr != trkCaloDets.end(); trkDetItr++){
489  std::cout<<trkDetItr->trkItr->p()<<" "<<trkDetItr->trkItr->eta()<<" "<<trkDetItr->okECAL<<" ";
490  if(trkDetItr->detIdECAL.subdetId() == EcalBarrel) std::cout << (EBDetId)trkDetItr->detIdECAL <<" ";
491  else std::cout << (EEDetId)trkDetItr->detIdECAL <<" ";
492  std::cout<<trkDetItr->okHCAL<<" ";
493  if(trkDetItr->okHCAL) std::cout<<(HcalDetId)trkDetItr->detIdHCAL;
494  std::cout << std::endl;
495  }
496  }
497 
498  int nvtxTracks=0;
499  for(trkDetItr = trkCaloDets.begin(),nTracks=0; trkDetItr != trkCaloDets.end(); trkDetItr++,nTracks++){
500 
501  const reco::Track* pTrack = &(*(trkDetItr->trkItr));
502 
503  // find vertex index the track is associated with
504  int pVtxTkId = -1;
505  for(unsigned int ind=0; ind<recVtxs->size(); ind++) {
506  if (!((*recVtxs)[ind].isFake())) {
507  reco::Vertex::trackRef_iterator vtxTrack = (*recVtxs)[ind].tracks_begin();
508  for (vtxTrack = (*recVtxs)[ind].tracks_begin(); vtxTrack!=(*recVtxs)[ind].tracks_end(); vtxTrack++) {
509 
510  const edm::RefToBase<reco::Track> pvtxTrack = (*vtxTrack);
511  if ( pTrack == pvtxTrack.get() ) {
512  pVtxTkId = ind;
513  break;
514  if(myverbose_>2) {
515  if( pTrack->pt()>1.0) {
516  std::cout<<"Debug the track association with vertex "<<std::endl;
517  std::cout<< pTrack << " "<< pvtxTrack.get() << std::endl;
518  std::cout<<" trkVtxIndex "<<nvtxTracks<<" vtx "<<ind<<" pt "<<pTrack->pt()
519  <<" eta "<<pTrack->eta()<<" "<<pTrack->pt()-pvtxTrack->pt()
520  <<" "<< pTrack->eta()-pvtxTrack->eta()
521  << std::endl;
522  nvtxTracks++;
523  }
524  }
525  }
526  }
527  }
528  }
529 
530  const reco::HitPattern& hitp = pTrack->hitPattern();
531  const reco::HitPattern& hitpIn = pTrack->trackerExpectedHitsInner();
532  const reco::HitPattern& hitpOut = pTrack->trackerExpectedHitsOuter();
533 
534  int nLayersCrossed = hitp.trackerLayersWithMeasurement() ;
535  int nOuterHits = hitp.stripTOBLayersWithMeasurement()+hitp.stripTECLayersWithMeasurement() ;
536 
537  bool ifGood = pTrack->quality(trackQuality_);
538  double pt1 = pTrack->pt();
539  double p1 = pTrack->p();
540  double eta1 = pTrack->momentum().eta();
541  double phi1 = pTrack->momentum().phi();
542  double etaEcal1 = trkDetItr->etaECAL;
543  double phiEcal1 = trkDetItr->phiECAL;
544  double etaHcal1 = trkDetItr->etaHCAL;
545  double phiHcal1 = trkDetItr->phiHCAL;
546  double dxy1 = pTrack->dxy();
547  double dz1 = pTrack->dz();
548  double chisq1 = pTrack->normalizedChi2();
549  double dxybs1 = beamSpotH.isValid() ? pTrack->dxy(bspot) : pTrack->dxy();
550  double dzbs1 = beamSpotH.isValid() ? pTrack->dz(bspot) : pTrack->dz();
551  double dxypv1 = pTrack->dxy();
552  double dzpv1 = pTrack->dz();
553  if(pVtxTkId>=0) {
554  math::XYZPoint thisTkPV = math::XYZPoint( (*recVtxs)[pVtxTkId].x(),(*recVtxs)[pVtxTkId].y(), (*recVtxs)[pVtxTkId].z() );
555  dxypv1 = pTrack->dxy(thisTkPV);
556  dzpv1 = pTrack->dz (thisTkPV);
557  }
558 
559  h_recEtaPt_0->Fill(eta1, pt1);
560  h_recEtaP_0 ->Fill(eta1, p1);
561  h_recPt_0 ->Fill(pt1);
562  h_recP_0 ->Fill(p1);
563  h_recEta_0 ->Fill(eta1);
564  h_recPhi_0 ->Fill(phi1);
565 
566  if(ifGood && nLayersCrossed>7 ) {
567  h_recEtaPt_1->Fill(eta1, pt1);
568  h_recEtaP_1 ->Fill(eta1, p1);
569  h_recPt_1 ->Fill(pt1);
570  h_recP_1 ->Fill(p1);
571  h_recEta_1 ->Fill(eta1);
572  h_recPhi_1 ->Fill(phi1);
573  }
574 
575  if( ! ifGood ) continue;
576 
577  if( writeAllTracks && p1>2.0 && nLayersCrossed>7) {
578  t_trackPAll ->push_back( p1 );
579  t_trackEtaAll ->push_back( eta1 );
580  t_trackPhiAll ->push_back( phi1 );
581  t_trackPtAll ->push_back( pt1 );
582  t_trackDxyAll ->push_back( dxy1 );
583  t_trackDzAll ->push_back( dz1 );
584  t_trackDxyPVAll ->push_back( dxypv1 );
585  t_trackDzPVAll ->push_back( dzpv1 );
586  t_trackChiSqAll ->push_back( chisq1 );
587  }
588  if (doMC) {
589  edm::SimTrackContainer::const_iterator matchedSimTrkAll = spr::matchedSimTrack(iEvent, SimTk, SimVtx, pTrack, *associate, false);
590  if( writeAllTracks && matchedSimTrkAll != SimTk->end()) t_trackPdgIdAll->push_back( matchedSimTrkAll->type() );
591  }
592 
593  if( pt1>minTrackP_ && std::abs(eta1)<maxTrackEta_ && trkDetItr->okECAL) {
594 
595  double maxNearP31x31=999.0, maxNearP25x25=999.0, maxNearP21x21=999.0, maxNearP15x15=999.0;
596  maxNearP31x31 = spr::chargeIsolationEcal(nTracks, trkCaloDets, geo, caloTopology, 15,15);
597  maxNearP25x25 = spr::chargeIsolationEcal(nTracks, trkCaloDets, geo, caloTopology, 12,12);
598  maxNearP21x21 = spr::chargeIsolationEcal(nTracks, trkCaloDets, geo, caloTopology, 10,10);
599  maxNearP15x15 = spr::chargeIsolationEcal(nTracks, trkCaloDets, geo, caloTopology, 7, 7);
600 
601  int iTrkEtaBin=-1, iTrkMomBin=-1;
602  for(unsigned int ieta=0; ieta<NEtaBins; ieta++) {
603  if(std::abs(eta1)>genPartEtaBins[ieta] && std::abs(eta1)<genPartEtaBins[ieta+1] ) iTrkEtaBin = ieta;
604  }
605  for(unsigned int ipt=0; ipt<NPBins; ipt++) {
606  if( p1>genPartPBins[ipt] && p1<genPartPBins[ipt+1] ) iTrkMomBin = ipt;
607  }
608  if( iTrkMomBin>=0 && iTrkEtaBin>=0 ) {
609  h_maxNearP31x31[iTrkMomBin][iTrkEtaBin]->Fill( maxNearP31x31 );
610  h_maxNearP25x25[iTrkMomBin][iTrkEtaBin]->Fill( maxNearP25x25 );
611  h_maxNearP21x21[iTrkMomBin][iTrkEtaBin]->Fill( maxNearP21x21 );
612  h_maxNearP15x15[iTrkMomBin][iTrkEtaBin]->Fill( maxNearP15x15 );
613  }
614  if( maxNearP31x31<0.0 && nLayersCrossed>7 && nOuterHits>4) {
615  h_recEtaPt_2->Fill(eta1, pt1);
616  h_recEtaP_2 ->Fill(eta1, p1);
617  h_recPt_2 ->Fill(pt1);
618  h_recP_2 ->Fill(p1);
619  h_recEta_2 ->Fill(eta1);
620  h_recPhi_2 ->Fill(phi1);
621  }
622 
623  // if charge isolated in ECAL, store the further quantities
624  if( maxNearP31x31<1.0) {
625 
626  haveIsoTrack = true;
627 
628  // get the matching simTrack
629  double simTrackP = -1;
630  if (doMC) {
631  edm::SimTrackContainer::const_iterator matchedSimTrk = spr::matchedSimTrack(iEvent, SimTk, SimVtx, pTrack, *associate, false);
632  if( matchedSimTrk != SimTk->end() )simTrackP = matchedSimTrk->momentum().P();
633  }
634 
636  iSetup.get<EcalSeverityLevelAlgoRcd>().get(sevlv);
637 
638  // get ECal Tranverse Profile
639  std::pair<double, bool> e7x7P, e9x9P, e11x11P, e15x15P;
640  std::pair<double, bool> e7x7_10SigP, e9x9_10SigP, e11x11_10SigP, e15x15_10SigP;
641  std::pair<double, bool> e7x7_15SigP, e9x9_15SigP, e11x11_15SigP, e15x15_15SigP;
642  std::pair<double, bool> e7x7_20SigP, e9x9_20SigP, e11x11_20SigP, e15x15_20SigP;
643  std::pair<double, bool> e7x7_25SigP, e9x9_25SigP, e11x11_25SigP, e15x15_25SigP;
644  std::pair<double, bool> e7x7_30SigP, e9x9_30SigP, e11x11_30SigP, e15x15_30SigP;
645 
646  spr::caloSimInfo simInfo3x3, simInfo5x5, simInfo7x7, simInfo9x9;
647  spr::caloSimInfo simInfo11x11, simInfo13x13, simInfo15x15, simInfo21x21, simInfo25x25, simInfo31x31;
648  double trkEcalEne=0;
649 
650  const DetId isoCell = trkDetItr->detIdECAL;
651  e7x7P = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology,sevlv.product(),3,3, -100.0, -100.0, tMinE_,tMaxE_);
652  e9x9P = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology,sevlv.product(),4,4, -100.0, -100.0, tMinE_,tMaxE_);
653  e11x11P = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology,sevlv.product(),5,5, -100.0, -100.0, tMinE_,tMaxE_);
654  e15x15P = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology,sevlv.product(),7,7, -100.0, -100.0, tMinE_,tMaxE_);
655 
656  e7x7_10SigP = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology,sevlv.product(),3,3, 0.030, 0.150, tMinE_,tMaxE_);
657  e9x9_10SigP = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology,sevlv.product(),4,4, 0.030, 0.150, tMinE_,tMaxE_);
658  e11x11_10SigP = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology,sevlv.product(),5,5, 0.030, 0.150, tMinE_,tMaxE_);
659  e15x15_10SigP = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology,sevlv.product(),7,7, 0.030, 0.150, tMinE_,tMaxE_);
660 
661  e7x7_15SigP = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology, sevlv.product(),ttMap, 3,3, 0.20,0.45, tMinE_,tMaxE_);
662  e9x9_15SigP = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology, sevlv.product(),ttMap, 4,4, 0.20,0.45, tMinE_,tMaxE_);
663  e11x11_15SigP = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology,sevlv.product(), ttMap, 5,5, 0.20,0.45, tMinE_,tMaxE_);
664  e15x15_15SigP = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology, sevlv.product(),ttMap, 7,7, 0.20,0.45, tMinE_,tMaxE_, false);
665 
666  e7x7_20SigP = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology,sevlv.product(),3,3, 0.060, 0.300, tMinE_,tMaxE_);
667  e9x9_20SigP = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology,sevlv.product(),4,4, 0.060, 0.300, tMinE_,tMaxE_);
668  e11x11_20SigP = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology,sevlv.product(),5,5, 0.060, 0.300, tMinE_,tMaxE_);
669  e15x15_20SigP = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology,sevlv.product(),7,7, 0.060, 0.300, tMinE_,tMaxE_);
670 
671  e7x7_25SigP = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology,sevlv.product(),3,3, 0.075, 0.375, tMinE_,tMaxE_);
672  e9x9_25SigP = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology,sevlv.product(),4,4, 0.075, 0.375, tMinE_,tMaxE_);
673  e11x11_25SigP = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology,sevlv.product(),5,5, 0.075, 0.375, tMinE_,tMaxE_);
674  e15x15_25SigP = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology,sevlv.product(),7,7, 0.075, 0.375, tMinE_,tMaxE_);
675 
676  e7x7_30SigP = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology,sevlv.product(),3,3, 0.090, 0.450, tMinE_,tMaxE_);
677  e9x9_30SigP = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology,sevlv.product(),4,4, 0.090, 0.450, tMinE_,tMaxE_);
678  e11x11_30SigP = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology,sevlv.product(),5,5, 0.090, 0.450, tMinE_,tMaxE_);
679  e15x15_30SigP = spr::eECALmatrix(isoCell,barrelRecHitsHandle,endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology,sevlv.product(),7,7, 0.090, 0.450, tMinE_,tMaxE_);
680  if(myverbose_ == 2) {
681  std::cout << "clean ecal rechit " << std::endl;
682  std::cout<<"e7x7 "<<e7x7P.first<<" e9x9 "<<e9x9P.first<<" e11x11 " << e11x11P.first << " e15x15 "<<e15x15P.first<<std::endl;
683  std::cout<<"e7x7_10Sig "<<e7x7_10SigP.first<<" e11x11_10Sig "<<e11x11_10SigP.first<<" e15x15_10Sig "<<e15x15_10SigP.first<<std::endl;
684  }
685 
686  if (doMC) {
687  // check the energy from SimHits
688  spr::eECALSimInfo(iEvent,isoCell,geo,caloTopology,pcaloeb,pcaloee,SimTk,SimVtx,pTrack, *associate, 1,1, simInfo3x3);
689  spr::eECALSimInfo(iEvent,isoCell,geo,caloTopology,pcaloeb,pcaloee,SimTk,SimVtx,pTrack, *associate, 2,2, simInfo5x5);
690  spr::eECALSimInfo(iEvent,isoCell,geo,caloTopology,pcaloeb,pcaloee,SimTk,SimVtx,pTrack, *associate, 3,3, simInfo7x7);
691  spr::eECALSimInfo(iEvent,isoCell,geo,caloTopology,pcaloeb,pcaloee,SimTk,SimVtx,pTrack, *associate, 4,4, simInfo9x9);
692  spr::eECALSimInfo(iEvent,isoCell,geo,caloTopology,pcaloeb,pcaloee,SimTk,SimVtx,pTrack, *associate, 5,5, simInfo11x11);
693  spr::eECALSimInfo(iEvent,isoCell,geo,caloTopology,pcaloeb,pcaloee,SimTk,SimVtx,pTrack, *associate, 6,6, simInfo13x13);
694  spr::eECALSimInfo(iEvent,isoCell,geo,caloTopology,pcaloeb,pcaloee,SimTk,SimVtx,pTrack, *associate, 7,7, simInfo15x15, 150.0,false);
695  spr::eECALSimInfo(iEvent,isoCell,geo,caloTopology,pcaloeb,pcaloee,SimTk,SimVtx,pTrack, *associate, 10,10, simInfo21x21);
696  spr::eECALSimInfo(iEvent,isoCell,geo,caloTopology,pcaloeb,pcaloee,SimTk,SimVtx,pTrack, *associate, 12,12, simInfo25x25);
697  spr::eECALSimInfo(iEvent,isoCell,geo,caloTopology,pcaloeb,pcaloee,SimTk,SimVtx,pTrack, *associate, 15,15, simInfo31x31);
698 
699  trkEcalEne = spr::eCaloSimInfo(iEvent, geo, pcaloeb,pcaloee, SimTk, SimVtx, pTrack, *associate);
700  if(myverbose_ == 1) {
701  std::cout << "Track momentum " << pt1 << std::endl;
702 
703  std::cout << "ecal siminfo " << std::endl;
704  std::cout << "simInfo3x3: " << "eTotal " << simInfo3x3.eTotal << " eMatched " << simInfo3x3.eMatched << " eRest " << simInfo3x3.eRest << " eGamma "<<simInfo3x3.eGamma << " eNeutralHad " << simInfo3x3.eNeutralHad << " eChargedHad " << simInfo3x3.eChargedHad << std::endl;
705  std::cout << "simInfo5x5: " << "eTotal " << simInfo5x5.eTotal << " eMatched " << simInfo5x5.eMatched << " eRest " << simInfo5x5.eRest << " eGamma "<<simInfo5x5.eGamma << " eNeutralHad " << simInfo5x5.eNeutralHad << " eChargedHad " << simInfo5x5.eChargedHad << std::endl;
706  std::cout << "simInfo7x7: " << "eTotal " << simInfo7x7.eTotal << " eMatched " << simInfo7x7.eMatched << " eRest " << simInfo7x7.eRest << " eGamma "<<simInfo7x7.eGamma << " eNeutralHad " << simInfo7x7.eNeutralHad << " eChargedHad " << simInfo7x7.eChargedHad << std::endl;
707  std::cout << "simInfo9x9: " << "eTotal " << simInfo9x9.eTotal << " eMatched " << simInfo9x9.eMatched << " eRest " << simInfo9x9.eRest << " eGamma "<<simInfo9x9.eGamma << " eNeutralHad " << simInfo9x9.eNeutralHad << " eChargedHad " << simInfo9x9.eChargedHad << std::endl;
708  std::cout << "simInfo11x11: " << "eTotal " << simInfo11x11.eTotal << " eMatched " << simInfo11x11.eMatched << " eRest " << simInfo11x11.eRest << " eGamma "<<simInfo11x11.eGamma << " eNeutralHad " << simInfo11x11.eNeutralHad << " eChargedHad " << simInfo11x11.eChargedHad << std::endl;
709  std::cout << "simInfo15x15: " << "eTotal " << simInfo15x15.eTotal << " eMatched " << simInfo15x15.eMatched << " eRest " << simInfo15x15.eRest << " eGamma "<<simInfo15x15.eGamma << " eNeutralHad " << simInfo15x15.eNeutralHad << " eChargedHad " << simInfo15x15.eChargedHad << std::endl;
710  std::cout << "simInfo31x31: " << "eTotal " << simInfo31x31.eTotal << " eMatched " << simInfo31x31.eMatched << " eRest " << simInfo31x31.eRest << " eGamma "<<simInfo31x31.eGamma << " eNeutralHad " << simInfo31x31.eNeutralHad << " eChargedHad " << simInfo31x31.eChargedHad << std::endl;
711  std::cout << "trkEcalEne" << trkEcalEne << std::endl;
712  }
713  }
714 
715  // ======= Get HCAL information
716  double hcalScale=1.0;
717  if( std::abs(pTrack->eta())<1.4 ) {
718  hcalScale=120.0;
719  } else {
720  hcalScale=135.0;
721  }
722 
723  double maxNearHcalP3x3=-1, maxNearHcalP5x5=-1, maxNearHcalP7x7=-1;
724  maxNearHcalP3x3 = spr::chargeIsolationHcal(nTracks, trkCaloDets, theHBHETopology, 1,1);
725  maxNearHcalP5x5 = spr::chargeIsolationHcal(nTracks, trkCaloDets, theHBHETopology, 2,2);
726  maxNearHcalP7x7 = spr::chargeIsolationHcal(nTracks, trkCaloDets, theHBHETopology, 3,3);
727 
728  double h3x3=0, h5x5=0, h7x7=0;
729  double h3x3Sig=0, h5x5Sig=0, h7x7Sig=0;
730  double trkHcalEne = 0;
731  spr::caloSimInfo hsimInfo3x3, hsimInfo5x5, hsimInfo7x7;
732 
733  if(trkDetItr->okHCAL) {
734  const DetId ClosestCell(trkDetItr->detIdHCAL);
735  // bool includeHO=false, bool algoNew=true, bool debug=false
736  h3x3 = spr::eHCALmatrix(theHBHETopology, ClosestCell, hbhe,1,1, false, true, -100.0, -100.0, -100.0, -100.0, tMinH_,tMaxH_);
737  h5x5 = spr::eHCALmatrix(theHBHETopology, ClosestCell, hbhe,2,2, false, true, -100.0, -100.0, -100.0, -100.0, tMinH_,tMaxH_);
738  h7x7 = spr::eHCALmatrix(theHBHETopology, ClosestCell, hbhe,3,3, false, true, -100.0, -100.0, -100.0, -100.0, tMinH_,tMaxH_);
739  h3x3Sig = spr::eHCALmatrix(theHBHETopology, ClosestCell, hbhe,1,1, false, true, 0.7, 0.8, -100.0, -100.0, tMinH_,tMaxH_);
740  h5x5Sig = spr::eHCALmatrix(theHBHETopology, ClosestCell, hbhe,2,2, false, true, 0.7, 0.8, -100.0, -100.0, tMinH_,tMaxH_);
741  h7x7Sig = spr::eHCALmatrix(theHBHETopology, ClosestCell, hbhe,3,3, false, true, 0.7, 0.8, -100.0, -100.0, tMinH_,tMaxH_);
742 
743  if(myverbose_==2) {
744  std::cout << "HCAL 3x3 " << h3x3 << " " << h3x3Sig << " 5x5 " << h5x5 << " " << h5x5Sig << " 7x7 " << h7x7 << " " << h7x7Sig << std::endl;
745  }
746 
747  if (doMC) {
748  spr::eHCALSimInfo(iEvent, theHBHETopology, ClosestCell, geo,pcalohh, SimTk, SimVtx, pTrack, *associate, 1,1, hsimInfo3x3);
749  spr::eHCALSimInfo(iEvent, theHBHETopology, ClosestCell, geo,pcalohh, SimTk, SimVtx, pTrack, *associate, 2,2, hsimInfo5x5);
750  spr::eHCALSimInfo(iEvent, theHBHETopology, ClosestCell, geo,pcalohh, SimTk, SimVtx, pTrack, *associate, 3,3, hsimInfo7x7, 150.0, false,false);
751  trkHcalEne = spr::eCaloSimInfo(iEvent, geo,pcalohh, SimTk, SimVtx, pTrack, *associate);
752  if(myverbose_ == 1) {
753  std::cout << "Hcal siminfo " << std::endl;
754  std::cout << "hsimInfo3x3: " << "eTotal " << hsimInfo3x3.eTotal << " eMatched " << hsimInfo3x3.eMatched << " eRest " << hsimInfo3x3.eRest << " eGamma "<<hsimInfo3x3.eGamma << " eNeutralHad " << hsimInfo3x3.eNeutralHad << " eChargedHad " << hsimInfo3x3.eChargedHad << std::endl;
755  std::cout << "hsimInfo5x5: " << "eTotal " << hsimInfo5x5.eTotal << " eMatched " << hsimInfo5x5.eMatched << " eRest " << hsimInfo5x5.eRest << " eGamma "<<hsimInfo5x5.eGamma << " eNeutralHad " << hsimInfo5x5.eNeutralHad << " eChargedHad " << hsimInfo5x5.eChargedHad << std::endl;
756  std::cout << "hsimInfo7x7: " << "eTotal " << hsimInfo7x7.eTotal << " eMatched " << hsimInfo7x7.eMatched << " eRest " << hsimInfo7x7.eRest << " eGamma "<<hsimInfo7x7.eGamma << " eNeutralHad " << hsimInfo7x7.eNeutralHad << " eChargedHad " << hsimInfo7x7.eChargedHad << std::endl;
757  std::cout << "trkHcalEne " << trkHcalEne << std::endl;
758  }
759  }
760 
761  // debug the ecal and hcal matrix
762  if(myverbose_==4) {
763  std::cout<<"Run "<<iEvent.id().run()<<" Event "<<iEvent.id().event()<<std::endl;
764  std::vector<std::pair<DetId,double> > v7x7 = spr::eHCALmatrixCell(theHBHETopology, ClosestCell, hbhe,3,3, false, false);
765  double sumv=0.0;
766 
767  for(unsigned int iv=0; iv<v7x7.size(); iv++) {
768  sumv += v7x7[iv].second;
769  }
770  std::cout<<"h7x7 "<<h7x7<<" v7x7 "<<sumv << " in " << v7x7.size() <<std::endl;
771  for(unsigned int iv=0; iv<v7x7.size(); iv++) {
772  HcalDetId id = v7x7[iv].first;
773  std::cout << " Cell " << iv << " 0x" << std::hex << v7x7[iv].first() << std::dec << " " << id << " Energy " << v7x7[iv].second << std::endl;
774  }
775  }
776 
777  }
778  if (doMC) {
779  trkHcalEne = spr::eCaloSimInfo(iEvent, geo,pcalohh, SimTk, SimVtx, pTrack, *associate);
780  }
781 
782  // ====================================================================================================
783  // get diff between track outermost hit position and the propagation point at outermost surface of tracker
784  std::pair<math::XYZPoint,double> point2_TK0 = spr::propagateTrackerEnd( pTrack, bField, false);
785  math::XYZPoint diff(pTrack->outerPosition().X()-point2_TK0.first.X(),
786  pTrack->outerPosition().Y()-point2_TK0.first.Y(),
787  pTrack->outerPosition().Z()-point2_TK0.first.Z() );
788  double trackOutPosOutHitDr = diff.R();
789  double trackL = point2_TK0.second;
790  if (myverbose_==5) {
791  std::cout<<" propagted "<<point2_TK0.first<<" "<< point2_TK0.first.eta()<<" "<<point2_TK0.first.phi()<<std::endl;
792  std::cout<<" outerPosition() "<< pTrack->outerPosition() << " "<< pTrack->outerPosition().eta()<< " " << pTrack->outerPosition().phi()<< std::endl;
793  std::cout<<"diff " << diff << " diffR " <<diff.R()<<" diffR/L "<<diff.R()/point2_TK0.second <<std::endl;
794  }
795 
796  for(unsigned int ind=0;ind<recVtxs->size();ind++) {
797  if (!((*recVtxs)[ind].isFake())) {
798  reco::Vertex::trackRef_iterator vtxTrack = (*recVtxs)[ind].tracks_begin();
799  if( DeltaR(eta1,phi1, (*vtxTrack)->eta(),(*vtxTrack)->phi()) < 0.01 ) t_trackPVIdx ->push_back( ind );
800  else t_trackPVIdx ->push_back( -1 );
801  }
802  }
803 
804  // Fill the tree Branches here
805  t_trackPVIdx ->push_back( pVtxTkId );
806  t_trackP ->push_back( p1 );
807  t_trackPt ->push_back( pt1 );
808  t_trackEta ->push_back( eta1 );
809  t_trackPhi ->push_back( phi1 );
810  t_trackEcalEta ->push_back( etaEcal1 );
811  t_trackEcalPhi ->push_back( phiEcal1 );
812  t_trackHcalEta ->push_back( etaHcal1 );
813  t_trackHcalPhi ->push_back( phiHcal1 );
814  t_trackDxy ->push_back( dxy1 );
815  t_trackDz ->push_back( dz1 );
816  t_trackDxyBS ->push_back( dxybs1 );
817  t_trackDzBS ->push_back( dzbs1 );
818  t_trackDxyPV ->push_back( dxypv1 );
819  t_trackDzPV ->push_back( dzpv1 );
820  t_trackChiSq ->push_back( chisq1 );
821  t_trackNOuterHits ->push_back( nOuterHits );
822  t_NLayersCrossed ->push_back( nLayersCrossed );
823 
824  t_trackHitsTOB ->push_back( hitp.stripTOBLayersWithMeasurement() );
825  t_trackHitsTEC ->push_back( hitp.stripTECLayersWithMeasurement() );
826  t_trackHitInMissTOB ->push_back( hitpIn.stripTOBLayersWithoutMeasurement() );
827  t_trackHitInMissTEC ->push_back( hitpIn.stripTECLayersWithoutMeasurement() );
828  t_trackHitInMissTIB ->push_back( hitpIn.stripTIBLayersWithoutMeasurement() );
831 
833  t_trackHitOutMissTEC ->push_back( hitpOut.stripTECLayersWithoutMeasurement() );
837 
838  t_trackHitInMeasTOB ->push_back( hitpIn.stripTOBLayersWithMeasurement() );
839  t_trackHitInMeasTEC ->push_back( hitpIn.stripTECLayersWithMeasurement() );
840  t_trackHitInMeasTIB ->push_back( hitpIn.stripTIBLayersWithMeasurement() );
841  t_trackHitInMeasTID ->push_back( hitpIn.stripTIDLayersWithMeasurement() );
842  t_trackHitOutMeasTOB ->push_back( hitpOut.stripTOBLayersWithMeasurement() );
843  t_trackHitOutMeasTEC ->push_back( hitpOut.stripTECLayersWithMeasurement() );
844  t_trackHitOutMeasTIB ->push_back( hitpOut.stripTIBLayersWithMeasurement() );
845  t_trackHitOutMeasTID ->push_back( hitpOut.stripTIDLayersWithMeasurement() );
846  t_trackOutPosOutHitDr ->push_back( trackOutPosOutHitDr );
847  t_trackL ->push_back( trackL );
848 
849  t_maxNearP31x31 ->push_back( maxNearP31x31 );
850  t_maxNearP21x21 ->push_back( maxNearP21x21 );
851 
852  t_ecalSpike11x11 ->push_back( e11x11P.second );
853  t_e7x7 ->push_back( e7x7P.first );
854  t_e9x9 ->push_back( e9x9P.first );
855  t_e11x11 ->push_back( e11x11P.first );
856  t_e15x15 ->push_back( e15x15P.first );
857 
858  t_e7x7_10Sig ->push_back( e7x7_10SigP.first );
859  t_e9x9_10Sig ->push_back( e9x9_10SigP.first );
860  t_e11x11_10Sig ->push_back( e11x11_10SigP.first );
861  t_e15x15_10Sig ->push_back( e15x15_10SigP.first );
862  t_e7x7_15Sig ->push_back( e7x7_15SigP.first );
863  t_e9x9_15Sig ->push_back( e9x9_15SigP.first );
864  t_e11x11_15Sig ->push_back( e11x11_15SigP.first );
865  t_e15x15_15Sig ->push_back( e15x15_15SigP.first );
866  t_e7x7_20Sig ->push_back( e7x7_20SigP.first );
867  t_e9x9_20Sig ->push_back( e9x9_20SigP.first );
868  t_e11x11_20Sig ->push_back( e11x11_20SigP.first );
869  t_e15x15_20Sig ->push_back( e15x15_20SigP.first );
870  t_e7x7_25Sig ->push_back( e7x7_25SigP.first );
871  t_e9x9_25Sig ->push_back( e9x9_25SigP.first );
872  t_e11x11_25Sig ->push_back( e11x11_25SigP.first );
873  t_e15x15_25Sig ->push_back( e15x15_25SigP.first );
874  t_e7x7_30Sig ->push_back( e7x7_30SigP.first );
875  t_e9x9_30Sig ->push_back( e9x9_30SigP.first );
876  t_e11x11_30Sig ->push_back( e11x11_30SigP.first );
877  t_e15x15_30Sig ->push_back( e15x15_30SigP.first );
878 
879  if (doMC) {
880  t_esim7x7 ->push_back( simInfo7x7.eTotal );
881  t_esim9x9 ->push_back( simInfo9x9.eTotal );
882  t_esim11x11 ->push_back( simInfo11x11.eTotal );
883  t_esim15x15 ->push_back( simInfo15x15.eTotal );
884 
885  t_esim7x7Matched ->push_back( simInfo7x7.eMatched );
886  t_esim9x9Matched ->push_back( simInfo9x9.eMatched );
887  t_esim11x11Matched ->push_back( simInfo11x11.eMatched );
888  t_esim15x15Matched ->push_back( simInfo15x15.eMatched );
889 
890  t_esim7x7Rest ->push_back( simInfo7x7.eRest );
891  t_esim9x9Rest ->push_back( simInfo9x9.eRest );
892  t_esim11x11Rest ->push_back( simInfo11x11.eRest );
893  t_esim15x15Rest ->push_back( simInfo15x15.eRest );
894 
895  t_esim7x7Photon ->push_back( simInfo7x7.eGamma );
896  t_esim9x9Photon ->push_back( simInfo9x9.eGamma );
897  t_esim11x11Photon ->push_back( simInfo11x11.eGamma );
898  t_esim15x15Photon ->push_back( simInfo15x15.eGamma );
899 
900  t_esim7x7NeutHad ->push_back( simInfo7x7.eNeutralHad );
901  t_esim9x9NeutHad ->push_back( simInfo9x9.eNeutralHad );
902  t_esim11x11NeutHad ->push_back( simInfo11x11.eNeutralHad );
903  t_esim15x15NeutHad ->push_back( simInfo15x15.eNeutralHad );
904 
905  t_esim7x7CharHad ->push_back( simInfo7x7.eChargedHad );
906  t_esim9x9CharHad ->push_back( simInfo9x9.eChargedHad );
907  t_esim11x11CharHad ->push_back( simInfo11x11.eChargedHad );
908  t_esim15x15CharHad ->push_back( simInfo15x15.eChargedHad );
909 
910  t_trkEcalEne ->push_back( trkEcalEne );
911  t_simTrackP ->push_back( simTrackP );
912  t_esimPdgId ->push_back( simInfo11x11.pdgMatched );
913  }
914 
915  t_maxNearHcalP3x3 ->push_back( maxNearHcalP3x3 );
916  t_maxNearHcalP5x5 ->push_back( maxNearHcalP5x5 );
917  t_maxNearHcalP7x7 ->push_back( maxNearHcalP7x7 );
918 
919  t_h3x3 ->push_back( h3x3 );
920  t_h5x5 ->push_back( h5x5 );
921  t_h7x7 ->push_back( h7x7 );
922  t_h3x3Sig ->push_back( h3x3Sig );
923  t_h5x5Sig ->push_back( h5x5Sig );
924  t_h7x7Sig ->push_back( h7x7Sig );
925 
926  t_infoHcal ->push_back( trkDetItr->okHCAL );
927  if (doMC) {
928  t_trkHcalEne ->push_back( hcalScale*trkHcalEne );
929 
930  t_hsim3x3 ->push_back( hcalScale*hsimInfo3x3.eTotal );
931  t_hsim5x5 ->push_back( hcalScale*hsimInfo5x5.eTotal );
932  t_hsim7x7 ->push_back( hcalScale*hsimInfo7x7.eTotal );
933 
934  t_hsim3x3Matched ->push_back( hcalScale*hsimInfo3x3.eMatched );
935  t_hsim5x5Matched ->push_back( hcalScale*hsimInfo5x5.eMatched );
936  t_hsim7x7Matched ->push_back( hcalScale*hsimInfo7x7.eMatched );
937 
938  t_hsim3x3Rest ->push_back( hcalScale*hsimInfo3x3.eRest );
939  t_hsim5x5Rest ->push_back( hcalScale*hsimInfo5x5.eRest );
940  t_hsim7x7Rest ->push_back( hcalScale*hsimInfo7x7.eRest );
941 
942  t_hsim3x3Photon ->push_back( hcalScale*hsimInfo3x3.eGamma );
943  t_hsim5x5Photon ->push_back( hcalScale*hsimInfo5x5.eGamma );
944  t_hsim7x7Photon ->push_back( hcalScale*hsimInfo7x7.eGamma );
945 
946  t_hsim3x3NeutHad ->push_back( hcalScale*hsimInfo3x3.eNeutralHad );
947  t_hsim5x5NeutHad ->push_back( hcalScale*hsimInfo5x5.eNeutralHad );
948  t_hsim7x7NeutHad ->push_back( hcalScale*hsimInfo7x7.eNeutralHad );
949 
950  t_hsim3x3CharHad ->push_back( hcalScale*hsimInfo3x3.eChargedHad );
951  t_hsim5x5CharHad ->push_back( hcalScale*hsimInfo5x5.eChargedHad );
952  t_hsim7x7CharHad ->push_back( hcalScale*hsimInfo7x7.eChargedHad );
953  }
954  /*
955  if(hcalScale*hsimInfo3x3["eTotal"] > 50.0) {
956 
957  std::cout << "Loosely Iso Track : eta " << eta1 << " Rec Mom " << p1 << " SimMom " << simTrackP << " h3x3 " << h3x3 << std::endl;
958 
959  std::cout <<"Closest cell Hcal (atHCAL) " << (HcalDetId)ClosestCell << std::endl;
960 
961  std::cout <<"trkHcalEne, etotal, matched, rest " <<hcalScale*trkHcalEne<<std::setw(15)<<hcalScale*hsimInfo3x3["eTotal"]
962  <<std::setw(15)<<hcalScale*hsimInfo3x3["eMatched"]<<std::setw(15)<<hcalScale*hsimInfo3x3["eRest"]
963  <<std::endl;
964  unsigned int nn = t_trkHcalEne->size();
965  std::cout <<"in Tree " << (*t_trkHcalEne)[nn-1] <<std::setw(15)<< (*t_hsim3x3)[nn-1]
966  <<std::setw(15)<< (*t_hsim3x3Matched)[nn-1] <<std::setw(15)<< (*t_hsim3x3Rest)[nn-1]
967  << std::endl;
968 
969  std::cout << "debug output \n" << std::endl;
970  std::map<std::string, double> hsimInfo3x3_debug = spr::eHCALSimInfo(iEvent, theHBHETopology, ClosestCell, geo,pcalohh, SimTk, SimVtx, pTrack, *associate, 1,1, 150.0, true);
971 
972  }
973  */
974 
975 
976  } // if loosely isolated track
977  } // check p1/eta
978  } // loop over track collection
979 
980  if(haveIsoTrack) tree->Fill();
981 }
982 
983 // ----- method called once each job just before starting event loop ----
985 
986  nEventProc=0;
987 
988  // double tempgen_TH[21] = { 1.0, 2.0, 3.0, 4.0, 5.0,
989  double tempgen_TH[16] = { 0.0, 1.0, 2.0, 3.0, 4.0,
990  5.0, 6.0, 7.0, 9.0, 11.0,
991  15.0, 20.0, 30.0, 50.0, 75.0, 100.0};
992 
993  for(int i=0; i<16; i++) genPartPBins[i] = tempgen_TH[i];
994 
995  double tempgen_Eta[4] = {0.0, 1.131, 1.653, 2.172};
996 
997  for(int i=0; i<4; i++) genPartEtaBins[i] = tempgen_Eta[i];
998 
999  BookHistograms();
1000 }
1001 
1002 // ----- method called once each job just after ending the event loop ----
1004 
1005  if(L1TriggerAlgoInfo_) {
1006  std::map< std::pair<unsigned int,std::string>, int>::iterator itr;
1007  for(itr=l1AlgoMap.begin(); itr!=l1AlgoMap.end(); itr++) {
1008  std::cout << " ****endjob**** " << (itr->first).first <<" "
1009  <<(itr->first).second <<" "<<itr->second
1010  << std::endl;
1011  int ibin = (itr->first).first;
1012  TString name( (itr->first).second );
1013  h_L1AlgoNames->GetXaxis()->SetBinLabel(ibin+1,name);
1014  }
1015  std::cout << "Number of Events Processed " << nEventProc << std::endl;
1016  }
1017 
1018 }
1019 
1020 //========================================================================================================
1021 
1023 
1024  t_PVx ->clear();
1025  t_PVy ->clear();
1026  t_PVz ->clear();
1027  t_PVisValid ->clear();
1028  t_PVndof ->clear();
1029  t_PVNTracks ->clear();
1030  t_PVNTracksWt ->clear();
1031  t_PVTracksSumPt ->clear();
1032  t_PVTracksSumPtWt ->clear();
1033  t_PVNTracksHP ->clear();
1034  t_PVNTracksHPWt ->clear();
1035  t_PVTracksSumPtHP ->clear();
1036  t_PVTracksSumPtHPWt ->clear();
1037 
1038  for(int i=0; i<128; i++) t_L1Decision[i]=0;
1039  t_L1AlgoNames ->clear();
1040  t_L1PreScale ->clear();
1041 
1042  t_L1CenJetPt ->clear();
1043  t_L1CenJetEta ->clear();
1044  t_L1CenJetPhi ->clear();
1045  t_L1FwdJetPt ->clear();
1046  t_L1FwdJetEta ->clear();
1047  t_L1FwdJetPhi ->clear();
1048  t_L1TauJetPt ->clear();
1049  t_L1TauJetEta ->clear();
1050  t_L1TauJetPhi ->clear();
1051  t_L1MuonPt ->clear();
1052  t_L1MuonEta ->clear();
1053  t_L1MuonPhi ->clear();
1054  t_L1IsoEMPt ->clear();
1055  t_L1IsoEMEta ->clear();
1056  t_L1IsoEMPhi ->clear();
1057  t_L1NonIsoEMPt ->clear();
1058  t_L1NonIsoEMEta ->clear();
1059  t_L1NonIsoEMPhi ->clear();
1060  t_L1METPt ->clear();
1061  t_L1METEta ->clear();
1062  t_L1METPhi ->clear();
1063 
1064  t_jetPt ->clear();
1065  t_jetEta ->clear();
1066  t_jetPhi ->clear();
1067  t_nTrksJetCalo ->clear();
1068  t_nTrksJetVtx ->clear();
1069 
1070  t_trackPAll ->clear();
1071  t_trackEtaAll ->clear();
1072  t_trackPhiAll ->clear();
1073  t_trackPdgIdAll ->clear();
1074  t_trackPtAll ->clear();
1075  t_trackDxyAll ->clear();
1076  t_trackDzAll ->clear();
1077  t_trackDxyPVAll ->clear();
1078  t_trackDzPVAll ->clear();
1079  t_trackChiSqAll ->clear();
1080 
1081  t_trackP ->clear();
1082  t_trackPt ->clear();
1083  t_trackEta ->clear();
1084  t_trackPhi ->clear();
1085  t_trackEcalEta ->clear();
1086  t_trackEcalPhi ->clear();
1087  t_trackHcalEta ->clear();
1088  t_trackHcalPhi ->clear();
1089  t_NLayersCrossed ->clear();
1090  t_trackNOuterHits ->clear();
1091  t_trackDxy ->clear();
1092  t_trackDxyBS ->clear();
1093  t_trackDz ->clear();
1094  t_trackDzBS ->clear();
1095  t_trackDxyPV ->clear();
1096  t_trackDzPV ->clear();
1097  t_trackChiSq ->clear();
1098  t_trackPVIdx ->clear();
1099  t_trackHitsTOB ->clear();
1100  t_trackHitsTEC ->clear();
1101  t_trackHitInMissTOB ->clear();
1102  t_trackHitInMissTEC ->clear();
1103  t_trackHitInMissTIB ->clear();
1104  t_trackHitInMissTID ->clear();
1105  t_trackHitInMissTIBTID ->clear();
1106  t_trackHitOutMissTOB ->clear();
1107  t_trackHitOutMissTEC ->clear();
1108  t_trackHitOutMissTIB ->clear();
1109  t_trackHitOutMissTID ->clear();
1110  t_trackHitOutMissTOBTEC ->clear();
1111  t_trackHitInMeasTOB ->clear();
1112  t_trackHitInMeasTEC ->clear();
1113  t_trackHitInMeasTIB ->clear();
1114  t_trackHitInMeasTID ->clear();
1115  t_trackHitOutMeasTOB ->clear();
1116  t_trackHitOutMeasTEC ->clear();
1117  t_trackHitOutMeasTIB ->clear();
1118  t_trackHitOutMeasTID ->clear();
1119  t_trackOutPosOutHitDr ->clear();
1120  t_trackL ->clear();
1121 
1122  t_maxNearP31x31 ->clear();
1123  t_maxNearP21x21 ->clear();
1124 
1125  t_ecalSpike11x11 ->clear();
1126  t_e7x7 ->clear();
1127  t_e9x9 ->clear();
1128  t_e11x11 ->clear();
1129  t_e15x15 ->clear();
1130 
1131  t_e7x7_10Sig ->clear();
1132  t_e9x9_10Sig ->clear();
1133  t_e11x11_10Sig ->clear();
1134  t_e15x15_10Sig ->clear();
1135  t_e7x7_15Sig ->clear();
1136  t_e9x9_15Sig ->clear();
1137  t_e11x11_15Sig ->clear();
1138  t_e15x15_15Sig ->clear();
1139  t_e7x7_20Sig ->clear();
1140  t_e9x9_20Sig ->clear();
1141  t_e11x11_20Sig ->clear();
1142  t_e15x15_20Sig ->clear();
1143  t_e7x7_25Sig ->clear();
1144  t_e9x9_25Sig ->clear();
1145  t_e11x11_25Sig ->clear();
1146  t_e15x15_25Sig ->clear();
1147  t_e7x7_30Sig ->clear();
1148  t_e9x9_30Sig ->clear();
1149  t_e11x11_30Sig ->clear();
1150  t_e15x15_30Sig ->clear();
1151 
1152  if (doMC) {
1153  t_simTrackP ->clear();
1154  t_esimPdgId ->clear();
1155  t_trkEcalEne ->clear();
1156 
1157  t_esim7x7 ->clear();
1158  t_esim9x9 ->clear();
1159  t_esim11x11 ->clear();
1160  t_esim15x15 ->clear();
1161 
1162  t_esim7x7Matched ->clear();
1163  t_esim9x9Matched ->clear();
1164  t_esim11x11Matched ->clear();
1165  t_esim15x15Matched ->clear();
1166 
1167  t_esim7x7Rest ->clear();
1168  t_esim9x9Rest ->clear();
1169  t_esim11x11Rest ->clear();
1170  t_esim15x15Rest ->clear();
1171 
1172  t_esim7x7Photon ->clear();
1173  t_esim9x9Photon ->clear();
1174  t_esim11x11Photon ->clear();
1175  t_esim15x15Photon ->clear();
1176 
1177  t_esim7x7NeutHad ->clear();
1178  t_esim9x9NeutHad ->clear();
1179  t_esim11x11NeutHad ->clear();
1180  t_esim15x15NeutHad ->clear();
1181 
1182  t_esim7x7CharHad ->clear();
1183  t_esim9x9CharHad ->clear();
1184  t_esim11x11CharHad ->clear();
1185  t_esim15x15CharHad ->clear();
1186  }
1187 
1188  t_maxNearHcalP3x3 ->clear();
1189  t_maxNearHcalP5x5 ->clear();
1190  t_maxNearHcalP7x7 ->clear();
1191 
1192  t_h3x3 ->clear();
1193  t_h5x5 ->clear();
1194  t_h7x7 ->clear();
1195  t_h3x3Sig ->clear();
1196  t_h5x5Sig ->clear();
1197  t_h7x7Sig ->clear();
1198 
1199  t_infoHcal ->clear();
1200 
1201  if (doMC) {
1202  t_trkHcalEne ->clear();
1203 
1204  t_hsim3x3 ->clear();
1205  t_hsim5x5 ->clear();
1206  t_hsim7x7 ->clear();
1207  t_hsim3x3Matched ->clear();
1208  t_hsim5x5Matched ->clear();
1209  t_hsim7x7Matched ->clear();
1210  t_hsim3x3Rest ->clear();
1211  t_hsim5x5Rest ->clear();
1212  t_hsim7x7Rest ->clear();
1213  t_hsim3x3Photon ->clear();
1214  t_hsim5x5Photon ->clear();
1215  t_hsim7x7Photon ->clear();
1216  t_hsim3x3NeutHad ->clear();
1217  t_hsim5x5NeutHad ->clear();
1218  t_hsim7x7NeutHad ->clear();
1219  t_hsim3x3CharHad ->clear();
1220  t_hsim5x5CharHad ->clear();
1221  t_hsim7x7CharHad ->clear();
1222  }
1223 }
1224 
1226 
1227  char hname[100], htit[100];
1228 
1229  TFileDirectory dir = fs->mkdir("nearMaxTrackP");
1230 
1231  for(unsigned int ieta=0; ieta<NEtaBins; ieta++) {
1232  double lowEta=-5.0, highEta= 5.0;
1233  lowEta = genPartEtaBins[ieta];
1234  highEta = genPartEtaBins[ieta+1];
1235 
1236  for(unsigned int ipt=0; ipt<NPBins; ipt++) {
1237  double lowP=0.0, highP=300.0;
1238  lowP = genPartPBins[ipt];
1239  highP = genPartPBins[ipt+1];
1240  sprintf(hname, "h_maxNearP31x31_ptBin%i_etaBin%i",ipt, ieta);
1241  sprintf(htit, "maxNearP in 31x31 (%3.2f<|#eta|<%3.2f), (%2.0f<trkP<%3.0f)", lowEta, highEta, lowP, highP );
1242  h_maxNearP31x31[ipt][ieta] = dir.make<TH1F>(hname, htit, 220, -2.0, 20.0);
1243  h_maxNearP31x31[ipt][ieta] ->Sumw2();
1244  sprintf(hname, "h_maxNearP25x25_ptBin%i_etaBin%i",ipt, ieta);
1245  sprintf(htit, "maxNearP in 25x25 (%3.2f<|#eta|<%3.2f), (%2.0f<trkP<%3.0f)", lowEta, highEta, lowP, highP );
1246  h_maxNearP25x25[ipt][ieta] = dir.make<TH1F>(hname, htit, 220, -2.0, 20.0);
1247  h_maxNearP25x25[ipt][ieta] ->Sumw2();
1248  sprintf(hname, "h_maxNearP21x21_ptBin%i_etaBin%i",ipt, ieta);
1249  sprintf(htit, "maxNearP in 21x21 (%3.2f<|#eta|<%3.2f), (%2.0f<trkP<%3.0f)", lowEta, highEta, lowP, highP );
1250  h_maxNearP21x21[ipt][ieta] = dir.make<TH1F>(hname, htit, 220, -2.0, 20.0);
1251  h_maxNearP21x21[ipt][ieta] ->Sumw2();
1252  sprintf(hname, "h_maxNearP15x15_ptBin%i_etaBin%i",ipt, ieta);
1253  sprintf(htit, "maxNearP in 15x15 (%3.2f<|#eta|<%3.2f), (%2.0f<trkP<%3.0f)", lowEta, highEta, lowP, highP );
1254  h_maxNearP15x15[ipt][ieta] = dir.make<TH1F>(hname, htit, 220, -2.0, 20.0);
1255  h_maxNearP15x15[ipt][ieta] ->Sumw2();
1256  }
1257  }
1258 
1259  h_L1AlgoNames = fs->make<TH1I>("h_L1AlgoNames", "h_L1AlgoNames:Bin Labels", 128, -0.5, 127.5);
1260 
1261  // Reconstructed Tracks
1262 
1263  h_PVTracksWt = fs->make<TH1F>("h_PVTracksWt", "h_PVTracksWt", 600, -0.1, 1.1);
1264 
1265  h_nTracks = fs->make<TH1F>("h_nTracks", "h_nTracks", 1000, -0.5, 999.5);
1266 
1267  sprintf(hname, "h_recEtaPt_0");
1268  sprintf(htit, "h_recEtaPt (all tracks Eta vs pT)");
1269  h_recEtaPt_0 = fs->make<TH2F>(hname, htit, 30, -3.0,3.0, 15, genPartPBins);
1270 
1271  sprintf(hname, "h_recEtaP_0");
1272  sprintf(htit, "h_recEtaP (all tracks Eta vs pT)");
1273  h_recEtaP_0 = fs->make<TH2F>(hname, htit, 30, -3.0,3.0, 15, genPartPBins);
1274 
1275  h_recPt_0 = fs->make<TH1F>("h_recPt_0", "Pt (all tracks)", 15, genPartPBins);
1276  h_recP_0 = fs->make<TH1F>("h_recP_0", "P (all tracks)", 15, genPartPBins);
1277  h_recEta_0 = fs->make<TH1F>("h_recEta_0", "Eta (all tracks)", 60, -3.0, 3.0);
1278  h_recPhi_0 = fs->make<TH1F>("h_recPhi_0", "Phi (all tracks)", 100, -3.2, 3.2);
1279  //-------------------------
1280  sprintf(hname, "h_recEtaPt_1");
1281  sprintf(htit, "h_recEtaPt (all good tracks Eta vs pT)");
1282  h_recEtaPt_1 = fs->make<TH2F>(hname, htit, 30, -3.0,3.0, 15, genPartPBins);
1283 
1284  sprintf(hname, "h_recEtaP_1");
1285  sprintf(htit, "h_recEtaP (all good tracks Eta vs pT)");
1286  h_recEtaP_1 = fs->make<TH2F>(hname, htit, 30, -3.0, 3.0, 15, genPartPBins);
1287 
1288  h_recPt_1 = fs->make<TH1F>("h_recPt_1", "Pt (all good tracks)", 15, genPartPBins);
1289  h_recP_1 = fs->make<TH1F>("h_recP_1", "P (all good tracks)", 15, genPartPBins);
1290  h_recEta_1 = fs->make<TH1F>("h_recEta_1", "Eta (all good tracks)", 60, -3.0, 3.0);
1291  h_recPhi_1 = fs->make<TH1F>("h_recPhi_1", "Phi (all good tracks)", 100, -3.2, 3.2);
1292  //-------------------------
1293  sprintf(hname, "h_recEtaPt_2");
1294  sprintf(htit, "h_recEtaPt (charge isolation Eta vs pT)");
1295  h_recEtaPt_2 = fs->make<TH2F>(hname, htit, 30, -3.0, 3.0, 15, genPartPBins);
1296 
1297  sprintf(hname, "h_recEtaP_2");
1298  sprintf(htit, "h_recEtaP (charge isolation Eta vs pT)");
1299  h_recEtaP_2 = fs->make<TH2F>(hname, htit, 30, -3.0, 3.0, 15, genPartPBins);
1300 
1301  h_recPt_2 = fs->make<TH1F>("h_recPt_2", "Pt (charge isolation)", 15, genPartPBins);
1302  h_recP_2 = fs->make<TH1F>("h_recP_2", "P (charge isolation)", 15, genPartPBins);
1303  h_recEta_2 = fs->make<TH1F>("h_recEta_2", "Eta (charge isolation)", 60, -3.0, 3.0);
1304  h_recPhi_2 = fs->make<TH1F>("h_recPhi_2", "Phi (charge isolation)", 100, -3.2, 3.2);
1305 
1306 
1307  tree = fs->make<TTree>("tree", "tree");
1308  tree->SetAutoSave(10000);
1309 
1310 
1311  tree->Branch("t_EvtNo" ,&t_EvtNo ,"t_EvtNo/I");
1312  tree->Branch("t_RunNo" ,&t_RunNo ,"t_RunNo/I");
1313  tree->Branch("t_Lumi" ,&t_Lumi ,"t_Lumi/I");
1314  tree->Branch("t_Bunch" ,&t_Bunch ,"t_Bunch/I");
1315 
1316 
1317  t_PVx = new std::vector<double>();
1318  t_PVy = new std::vector<double>();
1319  t_PVz = new std::vector<double>();
1320  t_PVisValid = new std::vector<int>();
1321  t_PVndof = new std::vector<int>();
1322  t_PVNTracks = new std::vector<int>();
1323  t_PVNTracksWt = new std::vector<int>();
1324  t_PVTracksSumPt = new std::vector<double>();
1325  t_PVTracksSumPtWt = new std::vector<double>();
1326  t_PVNTracksHP = new std::vector<int>();
1327  t_PVNTracksHPWt = new std::vector<int>();
1328  t_PVTracksSumPtHP = new std::vector<double>();
1329  t_PVTracksSumPtHPWt = new std::vector<double>();
1330 
1331  tree->Branch("PVx" ,"vector<double>" ,&t_PVx);
1332  tree->Branch("PVy" ,"vector<double>" ,&t_PVy);
1333  tree->Branch("PVz" ,"vector<double>" ,&t_PVz);
1334  tree->Branch("PVisValid" ,"vector<int>" ,&t_PVisValid);
1335  tree->Branch("PVndof" ,"vector<int>" ,&t_PVndof);
1336  tree->Branch("PVNTracks" ,"vector<int>" ,&t_PVNTracks);
1337  tree->Branch("PVNTracksWt" ,"vector<int>" ,&t_PVNTracksWt);
1338  tree->Branch("t_PVTracksSumPt" ,"vector<double>" ,&t_PVTracksSumPt);
1339  tree->Branch("t_PVTracksSumPtWt" ,"vector<double>" ,&t_PVTracksSumPtWt);
1340  tree->Branch("PVNTracksHP" ,"vector<int>" ,&t_PVNTracksHP);
1341  tree->Branch("PVNTracksHPWt" ,"vector<int>" ,&t_PVNTracksHPWt);
1342  tree->Branch("t_PVTracksSumPtHP" ,"vector<double>" ,&t_PVTracksSumPtHP);
1343  tree->Branch("t_PVTracksSumPtHPWt" ,"vector<double>" ,&t_PVTracksSumPtHPWt);
1344 
1345  //----- L1Trigger information
1346  for(int i=0; i<128; i++) t_L1Decision[i]=0;
1347  t_L1AlgoNames = new std::vector<std::string>();
1348  t_L1PreScale = new std::vector<int>();
1349  t_L1CenJetPt = new std::vector<double>();
1350  t_L1CenJetEta = new std::vector<double>();
1351  t_L1CenJetPhi = new std::vector<double>();
1352  t_L1FwdJetPt = new std::vector<double>();
1353  t_L1FwdJetEta = new std::vector<double>();
1354  t_L1FwdJetPhi = new std::vector<double>();
1355  t_L1TauJetPt = new std::vector<double>();
1356  t_L1TauJetEta = new std::vector<double>();
1357  t_L1TauJetPhi = new std::vector<double>();
1358  t_L1MuonPt = new std::vector<double>();
1359  t_L1MuonEta = new std::vector<double>();
1360  t_L1MuonPhi = new std::vector<double>();
1361  t_L1IsoEMPt = new std::vector<double>();
1362  t_L1IsoEMEta = new std::vector<double>();
1363  t_L1IsoEMPhi = new std::vector<double>();
1364  t_L1NonIsoEMPt = new std::vector<double>();
1365  t_L1NonIsoEMEta = new std::vector<double>();
1366  t_L1NonIsoEMPhi = new std::vector<double>();
1367  t_L1METPt = new std::vector<double>();
1368  t_L1METEta = new std::vector<double>();
1369  t_L1METPhi = new std::vector<double>();
1370 
1371  tree->Branch("t_L1Decision", t_L1Decision, "t_L1Decision[128]/I");
1372  tree->Branch("t_L1AlgoNames", "vector<string>", &t_L1AlgoNames);
1373  tree->Branch("t_L1PreScale", "vector<int>", &t_L1PreScale);
1374  tree->Branch("t_L1CenJetPt", "vector<double>", &t_L1CenJetPt);
1375  tree->Branch("t_L1CenJetEta", "vector<double>", &t_L1CenJetEta);
1376  tree->Branch("t_L1CenJetPhi", "vector<double>", &t_L1CenJetPhi);
1377  tree->Branch("t_L1FwdJetPt", "vector<double>", &t_L1FwdJetPt);
1378  tree->Branch("t_L1FwdJetEta", "vector<double>", &t_L1FwdJetEta);
1379  tree->Branch("t_L1FwdJetPhi", "vector<double>", &t_L1FwdJetPhi);
1380  tree->Branch("t_L1TauJetPt", "vector<double>", &t_L1TauJetPt);
1381  tree->Branch("t_L1TauJetEta", "vector<double>", &t_L1TauJetEta);
1382  tree->Branch("t_L1TauJetPhi", "vector<double>", &t_L1TauJetPhi);
1383  tree->Branch("t_L1MuonPt", "vector<double>", &t_L1MuonPt);
1384  tree->Branch("t_L1MuonEta", "vector<double>", &t_L1MuonEta);
1385  tree->Branch("t_L1MuonPhi", "vector<double>", &t_L1MuonPhi);
1386  tree->Branch("t_L1IsoEMPt", "vector<double>", &t_L1IsoEMPt);
1387  tree->Branch("t_L1IsoEMEta", "vector<double>", &t_L1IsoEMEta);
1388  tree->Branch("t_L1IsoEMPhi", "vector<double>", &t_L1IsoEMPhi);
1389  tree->Branch("t_L1NonIsoEMPt", "vector<double>", &t_L1NonIsoEMPt);
1390  tree->Branch("t_L1NonIsoEMEta", "vector<double>", &t_L1NonIsoEMEta);
1391  tree->Branch("t_L1NonIsoEMPhi", "vector<double>", &t_L1NonIsoEMPhi);
1392  tree->Branch("t_L1METPt", "vector<double>", &t_L1METPt);
1393  tree->Branch("t_L1METEta", "vector<double>", &t_L1METEta);
1394  tree->Branch("t_L1METPhi", "vector<double>", &t_L1METPhi);
1395 
1396  t_jetPt = new std::vector<double>();
1397  t_jetEta = new std::vector<double>();
1398  t_jetPhi = new std::vector<double>();
1399  t_nTrksJetCalo = new std::vector<double>();
1400  t_nTrksJetVtx = new std::vector<double>();
1401  tree->Branch("t_jetPt", "vector<double>",&t_jetPt);
1402  tree->Branch("t_jetEta", "vector<double>",&t_jetEta);
1403  tree->Branch("t_jetPhi", "vector<double>",&t_jetPhi);
1404  tree->Branch("t_nTrksJetCalo", "vector<double>",&t_nTrksJetCalo);
1405  tree->Branch("t_nTrksJetVtx", "vector<double>",&t_nTrksJetVtx);
1406 
1407  t_trackPAll = new std::vector<double>();
1408  t_trackEtaAll = new std::vector<double>();
1409  t_trackPhiAll = new std::vector<double>();
1410  t_trackPdgIdAll = new std::vector<double>();
1411  t_trackPtAll = new std::vector<double>();
1412  t_trackDxyAll = new std::vector<double>();
1413  t_trackDzAll = new std::vector<double>();
1414  t_trackDxyPVAll = new std::vector<double>();
1415  t_trackDzPVAll = new std::vector<double>();
1416  t_trackChiSqAll = new std::vector<double>();
1417  tree->Branch("t_trackPAll", "vector<double>", &t_trackPAll );
1418  tree->Branch("t_trackPhiAll", "vector<double>", &t_trackPhiAll );
1419  tree->Branch("t_trackEtaAll", "vector<double>", &t_trackEtaAll );
1420  tree->Branch("t_trackPtAll", "vector<double>", &t_trackPtAll );
1421  tree->Branch("t_trackDxyAll", "vector<double>", &t_trackDxyAll );
1422  tree->Branch("t_trackDzAll", "vector<double>", &t_trackDzAll );
1423  tree->Branch("t_trackDxyPVAll", "vector<double>", &t_trackDxyPVAll );
1424  tree->Branch("t_trackDzPVAll", "vector<double>", &t_trackDzPVAll );
1425  tree->Branch("t_trackChiSqAll", "vector<double>", &t_trackChiSqAll );
1426  //tree->Branch("t_trackPdgIdAll", "vector<double>", &t_trackPdgIdAll);
1427 
1428  t_trackP = new std::vector<double>();
1429  t_trackPt = new std::vector<double>();
1430  t_trackEta = new std::vector<double>();
1431  t_trackPhi = new std::vector<double>();
1432  t_trackEcalEta = new std::vector<double>();
1433  t_trackEcalPhi = new std::vector<double>();
1434  t_trackHcalEta = new std::vector<double>();
1435  t_trackHcalPhi = new std::vector<double>();
1436  t_trackNOuterHits = new std::vector<int>();
1437  t_NLayersCrossed = new std::vector<int>();
1438  t_trackDxy = new std::vector<double>();
1439  t_trackDxyBS = new std::vector<double>();
1440  t_trackDz = new std::vector<double>();
1441  t_trackDzBS = new std::vector<double>();
1442  t_trackDxyPV = new std::vector<double>();
1443  t_trackDzPV = new std::vector<double>();
1444  t_trackPVIdx = new std::vector<int>();
1445  t_trackChiSq = new std::vector<double>();
1446  t_trackHitsTOB = new std::vector<int>();
1447  t_trackHitsTEC = new std::vector<int>();
1448  t_trackHitInMissTOB = new std::vector<int>();
1449  t_trackHitInMissTEC = new std::vector<int>();
1450  t_trackHitInMissTIB = new std::vector<int>();
1451  t_trackHitInMissTID = new std::vector<int>();
1452  t_trackHitOutMissTOB = new std::vector<int>();
1453  t_trackHitOutMissTEC = new std::vector<int>();
1454  t_trackHitOutMissTIB = new std::vector<int>();
1455  t_trackHitOutMissTID = new std::vector<int>();
1456  t_trackHitInMissTIBTID = new std::vector<int>();
1457  t_trackHitOutMissTOB = new std::vector<int>();
1458  t_trackHitOutMissTEC = new std::vector<int>();
1459  t_trackHitOutMissTIB = new std::vector<int>();
1460  t_trackHitOutMissTID = new std::vector<int>();
1461  t_trackHitOutMissTOBTEC= new std::vector<int>();
1462  t_trackHitInMeasTOB = new std::vector<int>();
1463  t_trackHitInMeasTEC = new std::vector<int>();
1464  t_trackHitInMeasTIB = new std::vector<int>();
1465  t_trackHitInMeasTID = new std::vector<int>();
1466  t_trackHitOutMeasTOB = new std::vector<int>();
1467  t_trackHitOutMeasTEC = new std::vector<int>();
1468  t_trackHitOutMeasTIB = new std::vector<int>();
1469  t_trackHitOutMeasTID = new std::vector<int>();
1470  t_trackOutPosOutHitDr = new std::vector<double>();
1471  t_trackL = new std::vector<double>();
1472 
1473  tree->Branch("t_trackP", "vector<double>", &t_trackP );
1474  tree->Branch("t_trackPt", "vector<double>", &t_trackPt );
1475  tree->Branch("t_trackEta", "vector<double>", &t_trackEta );
1476  tree->Branch("t_trackPhi", "vector<double>", &t_trackPhi );
1477  tree->Branch("t_trackEcalEta", "vector<double>", &t_trackEcalEta );
1478  tree->Branch("t_trackEcalPhi", "vector<double>", &t_trackEcalPhi );
1479  tree->Branch("t_trackHcalEta", "vector<double>", &t_trackHcalEta );
1480  tree->Branch("t_trackHcalPhi", "vector<double>", &t_trackHcalPhi );
1481 
1482  tree->Branch("t_trackNOuterHits", "vector<int>", &t_trackNOuterHits );
1483  tree->Branch("t_NLayersCrossed", "vector<int>", &t_NLayersCrossed );
1484  tree->Branch("t_trackHitsTOB", "vector<int>", &t_trackHitsTOB );
1485  tree->Branch("t_trackHitsTEC", "vector<int>", &t_trackHitsTEC );
1486  tree->Branch("t_trackHitInMissTOB", "vector<int>", &t_trackHitInMissTOB );
1487  tree->Branch("t_trackHitInMissTEC", "vector<int>", &t_trackHitInMissTEC );
1488  tree->Branch("t_trackHitInMissTIB", "vector<int>", &t_trackHitInMissTIB );
1489  tree->Branch("t_trackHitInMissTID", "vector<int>", &t_trackHitInMissTID );
1490  tree->Branch("t_trackHitInMissTIBTID", "vector<int>", &t_trackHitInMissTIBTID );
1491  tree->Branch("t_trackHitOutMissTOB", "vector<int>", &t_trackHitOutMissTOB);
1492  tree->Branch("t_trackHitOutMissTEC", "vector<int>", &t_trackHitOutMissTEC);
1493  tree->Branch("t_trackHitOutMissTIB", "vector<int>", &t_trackHitOutMissTIB);
1494  tree->Branch("t_trackHitOutMissTID", "vector<int>", &t_trackHitOutMissTID);
1495  tree->Branch("t_trackHitOutMissTOBTEC","vector<int>", &t_trackHitOutMissTOBTEC);
1496  tree->Branch("t_trackHitInMeasTOB", "vector<int>", &t_trackHitInMeasTOB );
1497  tree->Branch("t_trackHitInMeasTEC", "vector<int>", &t_trackHitInMeasTEC );
1498  tree->Branch("t_trackHitInMeasTIB", "vector<int>", &t_trackHitInMeasTIB );
1499  tree->Branch("t_trackHitInMeasTID", "vector<int>", &t_trackHitInMeasTID );
1500  tree->Branch("t_trackHitOutMeasTOB", "vector<int>", &t_trackHitOutMeasTOB);
1501  tree->Branch("t_trackHitOutMeasTEC", "vector<int>", &t_trackHitOutMeasTEC);
1502  tree->Branch("t_trackHitOutMeasTIB", "vector<int>", &t_trackHitOutMeasTIB);
1503  tree->Branch("t_trackHitOutMeasTID", "vector<int>", &t_trackHitOutMeasTID);
1504  tree->Branch("t_trackOutPosOutHitDr", "vector<double>", &t_trackOutPosOutHitDr);
1505  tree->Branch("t_trackL", "vector<double>", &t_trackL);
1506 
1507  tree->Branch("t_trackDxy", "vector<double>", &t_trackDxy );
1508  tree->Branch("t_trackDxyBS", "vector<double>", &t_trackDxyBS );
1509  tree->Branch("t_trackDz", "vector<double>", &t_trackDz );
1510  tree->Branch("t_trackDzBS", "vector<double>", &t_trackDzBS );
1511  tree->Branch("t_trackDxyPV", "vector<double>", &t_trackDxyPV );
1512  tree->Branch("t_trackDzPV", "vector<double>", &t_trackDzPV );
1513  tree->Branch("t_trackChiSq", "vector<double>", &t_trackChiSq );
1514  tree->Branch("t_trackPVIdx", "vector<int>", &t_trackPVIdx );
1515 
1516  t_maxNearP31x31 = new std::vector<double>();
1517  t_maxNearP21x21 = new std::vector<double>();
1518 
1519  tree->Branch("t_maxNearP31x31", "vector<double>", &t_maxNearP31x31);
1520  tree->Branch("t_maxNearP21x21", "vector<double>", &t_maxNearP21x21);
1521 
1522  t_ecalSpike11x11 = new std::vector<int>();
1523  t_e7x7 = new std::vector<double>();
1524  t_e9x9 = new std::vector<double>();
1525  t_e11x11 = new std::vector<double>();
1526  t_e15x15 = new std::vector<double>();
1527 
1528  tree->Branch("t_ecalSpike11x11", "vector<int>", &t_ecalSpike11x11);
1529  tree->Branch("t_e7x7", "vector<double>", &t_e7x7);
1530  tree->Branch("t_e9x9", "vector<double>", &t_e9x9);
1531  tree->Branch("t_e11x11", "vector<double>", &t_e11x11);
1532  tree->Branch("t_e15x15", "vector<double>", &t_e15x15);
1533 
1534  t_e7x7_10Sig = new std::vector<double>();
1535  t_e9x9_10Sig = new std::vector<double>();
1536  t_e11x11_10Sig = new std::vector<double>();
1537  t_e15x15_10Sig = new std::vector<double>();
1538  t_e7x7_15Sig = new std::vector<double>();
1539  t_e9x9_15Sig = new std::vector<double>();
1540  t_e11x11_15Sig = new std::vector<double>();
1541  t_e15x15_15Sig = new std::vector<double>();
1542  t_e7x7_20Sig = new std::vector<double>();
1543  t_e9x9_20Sig = new std::vector<double>();
1544  t_e11x11_20Sig = new std::vector<double>();
1545  t_e15x15_20Sig = new std::vector<double>();
1546  t_e7x7_25Sig = new std::vector<double>();
1547  t_e9x9_25Sig = new std::vector<double>();
1548  t_e11x11_25Sig = new std::vector<double>();
1549  t_e15x15_25Sig = new std::vector<double>();
1550  t_e7x7_30Sig = new std::vector<double>();
1551  t_e9x9_30Sig = new std::vector<double>();
1552  t_e11x11_30Sig = new std::vector<double>();
1553  t_e15x15_30Sig = new std::vector<double>();
1554 
1555  tree->Branch("t_e7x7_10Sig" ,"vector<double>", &t_e7x7_10Sig);
1556  tree->Branch("t_e9x9_10Sig" ,"vector<double>", &t_e9x9_10Sig);
1557  tree->Branch("t_e11x11_10Sig" ,"vector<double>", &t_e11x11_10Sig);
1558  tree->Branch("t_e15x15_10Sig" ,"vector<double>", &t_e15x15_10Sig);
1559  tree->Branch("t_e7x7_15Sig" ,"vector<double>", &t_e7x7_15Sig);
1560  tree->Branch("t_e9x9_15Sig" ,"vector<double>", &t_e9x9_15Sig);
1561  tree->Branch("t_e11x11_15Sig" ,"vector<double>", &t_e11x11_15Sig);
1562  tree->Branch("t_e15x15_15Sig" ,"vector<double>", &t_e15x15_15Sig);
1563  tree->Branch("t_e7x7_20Sig" ,"vector<double>", &t_e7x7_20Sig);
1564  tree->Branch("t_e9x9_20Sig" ,"vector<double>", &t_e9x9_20Sig);
1565  tree->Branch("t_e11x11_20Sig" ,"vector<double>", &t_e11x11_20Sig);
1566  tree->Branch("t_e15x15_20Sig" ,"vector<double>", &t_e15x15_20Sig);
1567  tree->Branch("t_e7x7_25Sig" ,"vector<double>", &t_e7x7_25Sig);
1568  tree->Branch("t_e9x9_25Sig" ,"vector<double>", &t_e9x9_25Sig);
1569  tree->Branch("t_e11x11_25Sig" ,"vector<double>", &t_e11x11_25Sig);
1570  tree->Branch("t_e15x15_25Sig" ,"vector<double>", &t_e15x15_25Sig);
1571  tree->Branch("t_e7x7_30Sig" ,"vector<double>", &t_e7x7_30Sig);
1572  tree->Branch("t_e9x9_30Sig" ,"vector<double>", &t_e9x9_30Sig);
1573  tree->Branch("t_e11x11_30Sig" ,"vector<double>", &t_e11x11_30Sig);
1574  tree->Branch("t_e15x15_30Sig" ,"vector<double>", &t_e15x15_30Sig);
1575 
1576  if (doMC) {
1577  t_esim7x7 = new std::vector<double>();
1578  t_esim9x9 = new std::vector<double>();
1579  t_esim11x11 = new std::vector<double>();
1580  t_esim15x15 = new std::vector<double>();
1581 
1582  t_esim7x7Matched = new std::vector<double>();
1583  t_esim9x9Matched = new std::vector<double>();
1584  t_esim11x11Matched = new std::vector<double>();
1585  t_esim15x15Matched = new std::vector<double>();
1586 
1587  t_esim7x7Rest = new std::vector<double>();
1588  t_esim9x9Rest = new std::vector<double>();
1589  t_esim11x11Rest = new std::vector<double>();
1590  t_esim15x15Rest = new std::vector<double>();
1591 
1592  t_esim7x7Photon = new std::vector<double>();
1593  t_esim9x9Photon = new std::vector<double>();
1594  t_esim11x11Photon = new std::vector<double>();
1595  t_esim15x15Photon = new std::vector<double>();
1596 
1597  t_esim7x7NeutHad = new std::vector<double>();
1598  t_esim9x9NeutHad = new std::vector<double>();
1599  t_esim11x11NeutHad = new std::vector<double>();
1600  t_esim15x15NeutHad = new std::vector<double>();
1601 
1602  t_esim7x7CharHad = new std::vector<double>();
1603  t_esim9x9CharHad = new std::vector<double>();
1604  t_esim11x11CharHad = new std::vector<double>();
1605  t_esim15x15CharHad = new std::vector<double>();
1606 
1607  t_trkEcalEne = new std::vector<double>();
1608  t_simTrackP = new std::vector<double>();
1609  t_esimPdgId = new std::vector<double>();
1610 
1611  tree->Branch("t_esim7x7", "vector<double>", &t_esim7x7);
1612  tree->Branch("t_esim9x9", "vector<double>", &t_esim9x9);
1613  tree->Branch("t_esim11x11", "vector<double>", &t_esim11x11);
1614  tree->Branch("t_esim15x15", "vector<double>", &t_esim15x15);
1615 
1616  tree->Branch("t_esim7x7Matched", "vector<double>", &t_esim7x7Matched);
1617  tree->Branch("t_esim9x9Matched", "vector<double>", &t_esim9x9Matched);
1618  tree->Branch("t_esim11x11Matched", "vector<double>", &t_esim11x11Matched);
1619  tree->Branch("t_esim15x15Matched", "vector<double>", &t_esim15x15Matched);
1620 
1621  tree->Branch("t_esim7x7Rest", "vector<double>", &t_esim7x7Rest);
1622  tree->Branch("t_esim9x9Rest", "vector<double>", &t_esim9x9Rest);
1623  tree->Branch("t_esim11x11Rest", "vector<double>", &t_esim11x11Rest);
1624  tree->Branch("t_esim15x15Rest", "vector<double>", &t_esim15x15Rest);
1625 
1626  tree->Branch("t_esim7x7Photon", "vector<double>", &t_esim7x7Photon);
1627  tree->Branch("t_esim9x9Photon", "vector<double>", &t_esim9x9Photon);
1628  tree->Branch("t_esim11x11Photon", "vector<double>", &t_esim11x11Photon);
1629  tree->Branch("t_esim15x15Photon", "vector<double>", &t_esim15x15Photon);
1630 
1631 ;
1632  tree->Branch("t_esim7x7NeutHad", "vector<double>", &t_esim7x7NeutHad);
1633  tree->Branch("t_esim9x9NeutHad", "vector<double>", &t_esim9x9NeutHad);
1634  tree->Branch("t_esim11x11NeutHad", "vector<double>", &t_esim11x11NeutHad);
1635  tree->Branch("t_esim15x15NeutHad", "vector<double>", &t_esim15x15NeutHad);
1636 
1637  tree->Branch("t_esim7x7CharHad", "vector<double>", &t_esim7x7CharHad);
1638  tree->Branch("t_esim9x9CharHad", "vector<double>", &t_esim9x9CharHad);
1639  tree->Branch("t_esim11x11CharHad", "vector<double>", &t_esim11x11CharHad);
1640  tree->Branch("t_esim15x15CharHad", "vector<double>", &t_esim15x15CharHad);
1641 
1642  tree->Branch("t_trkEcalEne", "vector<double>", &t_trkEcalEne);
1643  tree->Branch("t_simTrackP", "vector<double>", &t_simTrackP);
1644  tree->Branch("t_esimPdgId", "vector<double>", &t_esimPdgId);
1645  }
1646 
1647  t_maxNearHcalP3x3 = new std::vector<double>();
1648  t_maxNearHcalP5x5 = new std::vector<double>();
1649  t_maxNearHcalP7x7 = new std::vector<double>();
1650  t_h3x3 = new std::vector<double>();
1651  t_h5x5 = new std::vector<double>();
1652  t_h7x7 = new std::vector<double>();
1653  t_h3x3Sig = new std::vector<double>();
1654  t_h5x5Sig = new std::vector<double>();
1655  t_h7x7Sig = new std::vector<double>();
1656  t_infoHcal = new std::vector<int>();
1657 
1658  if (doMC) {
1659  t_trkHcalEne = new std::vector<double>();
1660  t_hsim3x3 = new std::vector<double>();
1661  t_hsim5x5 = new std::vector<double>();
1662  t_hsim7x7 = new std::vector<double>();
1663  t_hsim3x3Matched = new std::vector<double>();
1664  t_hsim5x5Matched = new std::vector<double>();
1665  t_hsim7x7Matched = new std::vector<double>();
1666  t_hsim3x3Rest = new std::vector<double>();
1667  t_hsim5x5Rest = new std::vector<double>();
1668  t_hsim7x7Rest = new std::vector<double>();
1669  t_hsim3x3Photon = new std::vector<double>();
1670  t_hsim5x5Photon = new std::vector<double>();
1671  t_hsim7x7Photon = new std::vector<double>();
1672  t_hsim3x3NeutHad = new std::vector<double>();
1673  t_hsim5x5NeutHad = new std::vector<double>();
1674  t_hsim7x7NeutHad = new std::vector<double>();
1675  t_hsim3x3CharHad = new std::vector<double>();
1676  t_hsim5x5CharHad = new std::vector<double>();
1677  t_hsim7x7CharHad = new std::vector<double>();
1678  }
1679 
1680  tree->Branch("t_maxNearHcalP3x3", "vector<double>", &t_maxNearHcalP3x3);
1681  tree->Branch("t_maxNearHcalP5x5", "vector<double>", &t_maxNearHcalP5x5);
1682  tree->Branch("t_maxNearHcalP7x7", "vector<double>", &t_maxNearHcalP7x7);
1683  tree->Branch("t_h3x3", "vector<double>", &t_h3x3);
1684  tree->Branch("t_h5x5", "vector<double>", &t_h5x5);
1685  tree->Branch("t_h7x7", "vector<double>", &t_h7x7);
1686  tree->Branch("t_h3x3Sig", "vector<double>", &t_h3x3Sig);
1687  tree->Branch("t_h5x5Sig", "vector<double>", &t_h5x5Sig);
1688  tree->Branch("t_h7x7Sig", "vector<double>", &t_h7x7Sig);
1689  tree->Branch("t_infoHcal", "vector<int>", &t_infoHcal);
1690 
1691  if (doMC) {
1692  tree->Branch("t_trkHcalEne", "vector<double>", &t_trkHcalEne);
1693  tree->Branch("t_hsim3x3", "vector<double>", &t_hsim3x3);
1694  tree->Branch("t_hsim5x5", "vector<double>", &t_hsim5x5);
1695  tree->Branch("t_hsim7x7", "vector<double>", &t_hsim7x7);
1696  tree->Branch("t_hsim3x3Matched", "vector<double>", &t_hsim3x3Matched);
1697  tree->Branch("t_hsim5x5Matched", "vector<double>", &t_hsim5x5Matched);
1698  tree->Branch("t_hsim7x7Matched", "vector<double>", &t_hsim7x7Matched);
1699  tree->Branch("t_hsim3x3Rest", "vector<double>", &t_hsim3x3Rest);
1700  tree->Branch("t_hsim5x5Rest", "vector<double>", &t_hsim5x5Rest);
1701  tree->Branch("t_hsim7x7Rest", "vector<double>", &t_hsim7x7Rest);
1702  tree->Branch("t_hsim3x3Photon", "vector<double>", &t_hsim3x3Photon);
1703  tree->Branch("t_hsim5x5Photon", "vector<double>", &t_hsim5x5Photon);
1704  tree->Branch("t_hsim7x7Photon", "vector<double>", &t_hsim7x7Photon);
1705  tree->Branch("t_hsim3x3NeutHad", "vector<double>", &t_hsim3x3NeutHad);
1706  tree->Branch("t_hsim5x5NeutHad", "vector<double>", &t_hsim5x5NeutHad);
1707  tree->Branch("t_hsim7x7NeutHad", "vector<double>", &t_hsim7x7NeutHad);
1708  tree->Branch("t_hsim3x3CharHad", "vector<double>", &t_hsim3x3CharHad);
1709  tree->Branch("t_hsim5x5CharHad", "vector<double>", &t_hsim5x5CharHad);
1710  tree->Branch("t_hsim7x7CharHad", "vector<double>", &t_hsim7x7CharHad);
1711  }
1712  tree->Branch("t_nTracks", &t_nTracks, "t_nTracks/I");
1713 
1714 }
1715 
1716 
1717 double IsolatedTracksNxN::DeltaPhi(double v1, double v2) {
1718  // Computes the correctly normalized phi difference
1719  // v1, v2 = phi of object 1 and 2
1720 
1721  double pi = 3.141592654;
1722  double twopi = 6.283185307;
1723 
1724  double diff = std::abs(v2 - v1);
1725  double corr = twopi - diff;
1726  if (diff < pi){ return diff;} else { return corr;}
1727 }
1728 
1729 double IsolatedTracksNxN::DeltaR(double eta1, double phi1, double eta2, double phi2) {
1730  double deta = eta1 - eta2;
1731  double dphi = DeltaPhi(phi1, phi2);
1732  return std::sqrt(deta*deta + dphi*dphi);
1733 }
1734 
1736 
1737  std::string theTrackQuality = "highPurity";
1738  reco::TrackBase::TrackQuality trackQuality_ = reco::TrackBase::qualityByName(theTrackQuality);
1739 
1740  std::cout << " Reference Point " << pTrack->referencePoint() <<"\n"
1741  << " TrackMmentum " << pTrack->momentum()
1742  << " (pt,eta,phi)(" << pTrack->pt()<<","<<pTrack->eta()<<","<<pTrack->phi()<<")"
1743  << " p " << pTrack->p() << "\n"
1744  << " Normalized chi2 " << pTrack->normalizedChi2() <<" charge " << pTrack->charge()
1745  << " qoverp() " << pTrack->qoverp() <<"+-" << pTrack->qoverpError()
1746  << " d0 " << pTrack->d0() << "\n"
1747  << " NValidHits " << pTrack->numberOfValidHits() << " NLostHits " << pTrack->numberOfLostHits()
1748  << " TrackQuality " << pTrack->qualityName(trackQuality_) << " " << pTrack->quality(trackQuality_)
1749  << std::endl;
1750 
1751  if( printTrkHitPattern_ ) {
1752  const reco::HitPattern& p = pTrack->hitPattern();
1753  const reco::HitPattern& p1 = pTrack->trackerExpectedHitsInner();
1754  const reco::HitPattern& p2 = pTrack->trackerExpectedHitsOuter();
1755 
1756  std::cout<<"default " << std::endl;
1757  for (int i=0; i<p.numberOfHits(); i++) {
1759  }
1760  std::cout<<"trackerExpectedHitsInner() " << std::endl;
1761  for (int i=0; i<p1.numberOfHits(); i++) {
1763  }
1764  std::cout<<"trackerExpectedHitsOuter() " << std::endl;
1765  for (int i=0; i<p2.numberOfHits(); i++) {
1767  }
1768 
1769 
1770  std::cout << "\n \t trackerLayersWithMeasurement() " << p.trackerLayersWithMeasurement()
1771  << "\n \t pixelLayersWithMeasurement() " << p.pixelLayersWithMeasurement()
1772  << "\n \t stripLayersWithMeasurement() " << p.stripLayersWithMeasurement()
1773  << "\n \t pixelBarrelLayersWithMeasurement() " << p.pixelBarrelLayersWithMeasurement()
1774  << "\n \t pixelEndcapLayersWithMeasurement() " << p.pixelEndcapLayersWithMeasurement()
1775  << "\n \t stripTIBLayersWithMeasurement() " << p.stripTIBLayersWithMeasurement()
1776  << "\n \t stripTIDLayersWithMeasurement() " << p.stripTIDLayersWithMeasurement()
1777  << "\n \t stripTOBLayersWithMeasurement() " << p.stripTOBLayersWithMeasurement()
1778  << "\n \t stripTECLayersWithMeasurement() " << p.stripTECLayersWithMeasurement()
1779  << std::endl;
1780 
1781  }
1782 }
1783 
1784 //define this as a plug-in
RunNumber_t run() const
Definition: EventID.h:42
std::vector< double > * t_hsim7x7CharHad
double qoverp() const
q/p
Definition: TrackBase.h:113
double DeltaPhi(double v1, double v2)
double p() const
momentum vector magnitude
Definition: TrackBase.h:127
std::vector< double > * t_hsim3x3Matched
int stripTOBLayersWithMeasurement() const
Definition: HitPattern.cc:374
std::vector< double > * t_esim7x7CharHad
edm::EDGetTokenT< l1extra::L1JetParticleCollection > tok_L1extCenJet_
const std::string & gtTriggerMenuName() const
std::vector< double > * t_e15x15
std::vector< double > * t_nTrksJetCalo
T getParameter(std::string const &) const
EventNumber_t event() const
Definition: EventID.h:44
std::vector< double > * t_maxNearHcalP7x7
T getUntrackedParameter(std::string const &, T const &) const
std::vector< double > * t_trackPt
int i
Definition: DBlmapReader.cc:9
std::vector< double > * t_trackPAll
std::vector< double > * t_e7x7
std::vector< double > * t_trackHcalPhi
std::vector< double > * t_trackDxyPVAll
static std::string qualityName(TrackQuality)
Definition: TrackBase.h:404
std::vector< int > * t_trackHitOutMeasTEC
std::vector< double > * t_trackPdgIdAll
edm::InputTag JetExtender_
std::vector< spr::propagatedTrackID > propagateCALO(edm::Handle< reco::TrackCollection > &trkCollection, const CaloGeometry *geo, const MagneticField *bField, std::string &theTrackQuality, bool debug=false)
void eCaloSimInfo(std::vector< DetId > vdets, const CaloGeometry *geo, edm::Handle< T > &hitsEB, edm::Handle< T > &hitsEE, edm::Handle< edm::SimTrackContainer > &SimTk, edm::Handle< edm::SimVertexContainer > &SimVtx, edm::SimTrackContainer::const_iterator trkInfo, caloSimInfo &info, double timeCut=150, bool debug=false)
edm::InputTag L1extraCenJetSource_
std::vector< double > * t_trackOutPosOutHitDr
static const bool useL1GtTriggerMenuLite(false)
const Vector & momentum() const
track momentum vector
Definition: TrackBase.h:148
std::vector< double > * t_maxNearHcalP5x5
std::vector< double > * t_L1NonIsoEMPt
std::vector< int > * t_trackHitOutMeasTIB
std::vector< double > * t_esim7x7
std::vector< int > * t_trackHitInMissTOB
std::vector< double > * t_trackDz
double d0() const
dxy parameter in perigee convention (d0 = - dxy)
Definition: TrackBase.h:121
std::map< std::string, double > eHCALSimInfo(const edm::Event &, const HcalTopology *topology, const DetId &det, const CaloGeometry *geo, edm::Handle< T > &hits, edm::Handle< edm::SimTrackContainer > &SimTk, edm::Handle< edm::SimVertexContainer > &SimVtx, const reco::Track *pTrack, TrackerHitAssociator &associate, int ieta, int iphi, double timeCut=150, bool includeHO=false, bool debug=false)
std::vector< double > * t_trackEcalEta
int stripTIBLayersWithMeasurement() const
Definition: HitPattern.cc:356
double normalizedChi2() const
chi-squared divided by n.d.o.f. (or chi-squared * 1e6 if n.d.o.f. is zero)
Definition: TrackBase.h:109
std::vector< double > * t_trackEtaAll
std::vector< double > * t_esim11x11Matched
std::vector< double > * t_L1IsoEMEta
std::vector< double > * t_PVTracksSumPt
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:434
std::vector< double > * t_e9x9_10Sig
TH1F * h_maxNearP25x25[NPBins][NEtaBins]
TrackQuality
track quality
Definition: TrackBase.h:93
std::vector< double > * t_hsim3x3CharHad
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
std::vector< int > * t_PVndof
std::map< std::string, L1GtAlgorithm > AlgorithmMap
map containing the algorithms
std::vector< double > * t_hsim7x7Rest
std::vector< double > * t_L1TauJetPhi
std::vector< double > * t_PVTracksSumPtWt
const bool availableL1Configuration(int &errorCode, int &l1ConfCode) const
Definition: L1GtUtils.cc:2293
std::vector< int > * t_infoHcal
std::vector< double > * t_e11x11_15Sig
static const size_t NEtaBins
std::vector< double > * t_L1TauJetPt
std::vector< int > * t_PVNTracksWt
std::vector< double > * t_trackEcalPhi
std::vector< double > * t_esim11x11NeutHad
std::vector< double > * t_L1TauJetEta
std::vector< double > * t_e15x15_15Sig
std::vector< double > * t_L1MuonEta
std::vector< double > * t_L1MuonPt
int pixelLayersWithMeasurement() const
Definition: HitPattern.h:745
int bunchCrossing() const
Definition: EventBase.h:62
std::vector< double > * t_h3x3
edm::LuminosityBlockNumber_t luminosityBlock() const
Definition: EventBase.h:59
std::vector< double > * t_L1MuonPhi
double phi() const
azimuthal angle of momentum vector
Definition: TrackBase.h:137
unsigned short numberOfLostHits() const
number of cases where track crossed a layer without getting a hit.
Definition: TrackBase.h:232
int stripTIDLayersWithoutMeasurement() const
Definition: HitPattern.cc:422
std::vector< double > * t_trackDxyAll
T * make(const Args &...args) const
make new ROOT object
Definition: TFileService.h:64
std::vector< double > * t_maxNearHcalP3x3
edm::InputTag L1extraIsoEmSource_
IsolatedTracksNxN(const edm::ParameterSet &)
edm::EDGetTokenT< edm::SimVertexContainer > tok_simVtx_
std::vector< double > * t_hsim3x3Photon
edm::EDGetTokenT< reco::BeamSpot > tok_bs_
T eta() const
int ii
Definition: cuy.py:588
std::vector< double > * t_PVTracksSumPtHP
const math::XYZPoint & outerPosition() const
position of the outermost hit
Definition: Track.h:47
const Point & referencePoint() const
Reference point on the track.
Definition: TrackBase.h:151
std::vector< double > * t_esim11x11Rest
int pixelEndcapLayersWithMeasurement() const
Definition: HitPattern.cc:346
edm::InputTag L1GTReadoutRcdSource_
float float float z
std::vector< int > * t_trackNOuterHits
std::vector< double > * t_h7x7Sig
std::vector< double > * t_e15x15_30Sig
std::vector< double > * t_trackPhiAll
std::vector< double > * t_e7x7_25Sig
std::vector< double > * t_trkHcalEne
edm::InputTag L1extraFwdJetSource_
std::vector< double > * t_e7x7_20Sig
std::vector< double > * t_e11x11_10Sig
const Double_t pi
AlgorithmMap::const_iterator CItAlgo
iterators through map containing the algorithms
std::vector< double > * t_trackPtAll
std::vector< double > * t_h3x3Sig
edm::Service< TFileService > fs
double eHCALmatrix(const HcalTopology *topology, const DetId &det, edm::Handle< T > &hits, int ieta, int iphi, bool includeHO=false, bool algoNew=true, double hbThr=-100, double heThr=-100, double hfThr=-100, double hoThr=-100, double tMin=-500, double tMax=500, bool debug=false)
std::vector< double > * t_esimPdgId
double eNeutralHad
Definition: CaloSimInfo.h:47
double eChargedHad
Definition: CaloSimInfo.h:48
std::vector< double > * t_e11x11_20Sig
std::vector< double > * t_jetEta
std::map< std::pair< unsigned int, std::string >, int > l1AlgoMap
std::vector< double > * t_maxNearP31x31
std::vector< int > * t_trackHitInMissTEC
double chargeIsolationEcal(unsigned int trkIndex, std::vector< spr::propagatedTrackID > &vdetIds, const CaloGeometry *geo, const CaloTopology *caloTopology, int ieta, int iphi, bool debug=false)
int iEvent
Definition: GenABIO.cc:230
std::vector< double > * t_trackEta
const HitPattern & trackerExpectedHitsOuter() const
Access the hit pattern counting (in the Tracker) the number of expected crossed layers after the last...
Definition: TrackBase.h:225
const bool decision(const edm::Event &iEvent, const edm::InputTag &l1GtRecordInputTag, const edm::InputTag &l1GtReadoutRecordInputTag, const std::string &nameAlgoTechTrig, int &errorCode) const
Definition: L1GtUtils.cc:1448
std::vector< double > * t_esim9x9Photon
std::vector< int > * t_trackHitInMeasTID
double eta() const
pseudorapidity of momentum vector
Definition: TrackBase.h:139
std::vector< double > * t_hsim5x5Rest
std::vector< double > * t_e11x11
std::vector< double > * t_esim15x15
edm::EDGetTokenT< edm::SimTrackContainer > tok_simTk_
std::vector< double > * t_L1CenJetPt
std::vector< double > * t_trackChiSqAll
edm::InputTag L1extraTauJetSource_
std::vector< double > * t_e15x15_25Sig
std::vector< int > * t_trackHitOutMissTIB
edm::EDGetTokenT< edm::PCaloHitContainer > tok_caloHH_
std::vector< double > * t_esim15x15Rest
std::vector< double > * t_esim11x11
std::vector< double > * t_hsim7x7
std::vector< int > * t_PVNTracksHPWt
int stripTIDLayersWithMeasurement() const
Definition: HitPattern.cc:365
int trackerLayersWithMeasurement() const
Definition: HitPattern.h:740
T sqrt(T t)
Definition: SSEVec.h:48
std::vector< int > * t_trackHitOutMissTID
std::vector< int > * t_trackHitOutMissTOBTEC
vector< PseudoJet > jets
double pt() const
track transverse momentum
Definition: TrackBase.h:129
edm::InputTag HBHERecHitSource_
std::vector< int > * t_trackHitsTEC
std::vector< double > * t_e9x9_30Sig
const MagneticField * bField
edm::EDGetTokenT< edm::PCaloHitContainer > tok_caloEB_
std::vector< double > * t_L1CenJetPhi
std::vector< double > * t_L1FwdJetEta
std::vector< double > * t_h7x7
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
std::vector< int > * t_ecalSpike11x11
int numberOfHits() const
Definition: HitPattern.cc:211
std::vector< double > * t_esim15x15NeutHad
int stripTIBLayersWithoutMeasurement() const
Definition: HitPattern.cc:413
std::vector< double > * t_e9x9_20Sig
std::vector< int > * t_trackHitOutMeasTOB
std::vector< int > * t_trackHitInMeasTIB
const HitPattern & trackerExpectedHitsInner() const
Access the hit pattern counting (in the Tracker) the number of expected crossed layers before the fir...
Definition: TrackBase.h:223
edm::SimTrackContainer::const_iterator matchedSimTrack(const edm::Event &iEvent, edm::Handle< edm::SimTrackContainer > &SimTk, edm::Handle< edm::SimVertexContainer > &SimVtx, const reco::Track *pTrack, TrackerHitAssociator &associate, bool debug=false)
unsigned short numberOfValidHits() const
number of valid hits found
Definition: TrackBase.h:230
edm::EDGetTokenT< reco::TrackCollection > tok_genTrack_
T * make(const Args &...args) const
make new ROOT object
std::vector< double > * t_trackDzPV
std::vector< int > * t_PVNTracksHP
std::vector< double > * t_L1FwdJetPhi
edm::EDGetTokenT< l1extra::L1JetParticleCollection > tok_L1extFwdJet_
edm::EDGetTokenT< l1extra::L1EmParticleCollection > tok_L1extIsoEm_
const HitPattern & hitPattern() const
Access the hit pattern, indicating in which Tracker layers the track has hits.
Definition: TrackBase.h:221
static const bool useL1EventSetup(true)
std::vector< double > * t_nTrksJetVtx
bool first
Definition: L1TdeRCT.cc:75
bool isValid() const
Definition: HandleBase.h:76
double p2[4]
Definition: TauolaWrapper.h:90
std::vector< double > * t_esim7x7NeutHad
std::vector< double > * t_hsim3x3NeutHad
std::vector< double > * t_L1NonIsoEMPhi
std::vector< double > * t_trackDzBS
void analyze(const edm::Event &, const edm::EventSetup &)
std::vector< int > * t_trackHitOutMissTOB
std::vector< double > * t_L1CenJetEta
std::vector< double > * t_trackDzPVAll
JetCorrectorParameters corr
Definition: classes.h:5
edm::EDGetTokenT< l1extra::L1MuonParticleCollection > tok_L1extMu_
std::vector< double > * t_esim7x7Matched
std::vector< double > * t_esim9x9
std::vector< double > * t_hsim5x5
std::vector< double > * t_h5x5
double qoverpError() const
error on signed transverse curvature
Definition: TrackBase.h:190
std::vector< double > * t_trackDxyPV
std::vector< double > * t_hsim5x5NeutHad
std::vector< double > * t_trackChiSq
double dz() const
dz parameter (= dsz/cos(lambda)). This is the track z0 w.r.t (0,0,0) only if the refPoint is close to...
Definition: TrackBase.h:125
std::vector< double > * t_hsim5x5Photon
std::vector< double > * t_PVy
std::vector< int > * t_L1PreScale
void getL1GtRunCache(const edm::Run &, const edm::EventSetup &, const bool, const bool, const edm::InputTag &)
get all the run-constant quantities for L1 trigger and cache them
Definition: L1GtUtils.cc:311
std::vector< double > * t_hsim3x3Rest
Definition: DetId.h:18
edm::InputTag L1extraNonIsoEmSource_
std::vector< int > * t_PVisValid
static TrackQuality qualityByName(const std::string &name)
Definition: TrackBase.cc:46
std::vector< double > * t_L1METPhi
std::vector< double > * t_esim9x9Rest
std::vector< int > * t_trackHitInMissTIB
TFileDirectory mkdir(const std::string &dir, const std::string &descr="")
create a new subdirectory
Definition: TFileService.h:69
std::vector< double > * t_e9x9_15Sig
std::vector< double > * t_e7x7_15Sig
std::vector< double > * t_trackPhi
std::vector< double > * t_e7x7_10Sig
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
std::vector< double > * t_trkEcalEne
edm::EDGetTokenT< l1extra::L1EmParticleCollection > tok_L1extNoIsoEm_
std::vector< double > * t_trackDzAll
edm::EDGetTokenT< reco::VertexCollection > tok_recVtx_
TH1F * h_maxNearP31x31[NPBins][NEtaBins]
std::vector< double > * t_esim11x11CharHad
const T & get() const
Definition: EventSetup.h:55
int stripLayersWithMeasurement() const
Definition: HitPattern.h:750
int pixelBarrelLayersWithMeasurement() const
Definition: HitPattern.cc:336
T const * product() const
Definition: ESHandle.h:62
edm::EDGetTokenT< edm::PCaloHitContainer > tok_caloEE_
edm::EDGetTokenT< EcalRecHitCollection > tok_EB_
TH1F * h_maxNearP15x15[NPBins][NEtaBins]
std::vector< double > * t_e9x9_25Sig
void eECALSimInfo(const edm::Event &, const DetId &det, const CaloGeometry *geo, const CaloTopology *caloTopology, edm::Handle< T > &hitsEB, edm::Handle< T > &hitsEE, edm::Handle< edm::SimTrackContainer > &SimTk, edm::Handle< edm::SimVertexContainer > &SimVtx, const reco::Track *pTrack, TrackerHitAssociator &associate, int ieta, int iphi, caloSimInfo &info, double timeCut=150, bool debug=false)
std::vector< double > * t_esim7x7Photon
std::vector< double > * t_esim9x9NeutHad
std::vector< double > * t_L1FwdJetPt
std::vector< double > * t_simTrackP
std::vector< double > * t_hsim7x7Matched
T const * product() const
Definition: Handle.h:81
std::vector< int > * t_trackPVIdx
bool quality(const TrackQuality) const
Track quality.
Definition: TrackBase.h:375
std::vector< int > * t_trackHitsTOB
std::vector< double > * t_esim9x9Matched
std::vector< double > * t_hsim7x7NeutHad
std::vector< double > * t_L1IsoEMPhi
std::vector< int > * t_PVNTracks
edm::EDGetTokenT< HBHERecHitCollection > tok_hbhe_
std::vector< double > * t_maxNearP21x21
edm::EventID id() const
Definition: EventBase.h:56
std::vector< double > * t_L1NonIsoEMEta
double p1[4]
Definition: TauolaWrapper.h:89
std::vector< double > * t_PVz
std::pair< math::XYZPoint, double > propagateTrackerEnd(const reco::Track *, const MagneticField *, bool debug=false)
void printTrack(const reco::Track *pTrack)
int stripTECLayersWithoutMeasurement() const
Definition: HitPattern.cc:440
std::vector< double > * t_e15x15_20Sig
std::vector< double > * t_PVTracksSumPtHPWt
std::vector< double > * t_trackDxy
std::vector< std::string > * t_L1AlgoNames
std::vector< double > * t_hsim5x5CharHad
std::vector< double > * t_e15x15_10Sig
const AlgorithmMap & gtAlgorithmMap() const
get / set the algorithm map (by name)
int stripTECLayersWithMeasurement() const
Definition: HitPattern.cc:383
std::vector< TrackBaseRef >::const_iterator trackRef_iterator
The iteratator for the vector&lt;TrackRef&gt;
Definition: Vertex.h:37
TH1F * h_maxNearP21x21[NPBins][NEtaBins]
int stripTOBLayersWithoutMeasurement() const
Definition: HitPattern.cc:431
std::vector< double > * t_jetPhi
std::vector< double > * t_trackL
std::vector< double > * t_esim7x7Rest
std::vector< double > * t_trackDxyBS
tuple cout
Definition: gather_cfg.py:121
int charge() const
track electric charge
Definition: TrackBase.h:111
std::vector< double > * t_trackHcalEta
std::vector< double > * t_L1METPt
void printHitPattern(int position, std::ostream &stream) const
Definition: HitPattern.cc:563
edm::InputTag L1extraMuonSource_
dbl *** dir
Definition: mlp_gen.cc:35
volatile std::atomic< bool > shutdown_flag false
std::vector< int > * t_NLayersCrossed
edm::EDGetTokenT< reco::CaloJetCollection > tok_jets_
Definition: DDAxes.h:10
double DeltaR(double eta1, double phi1, double eta2, double phi2)
edm::InputTag JetSrc_
std::vector< std::pair< DetId, double > > eHCALmatrixCell(const HcalTopology *topology, const DetId &det, edm::Handle< T > &hits, int ieta, int iphi, bool includeHO=false, double hbThr=-100, double heThr=-100, double hfThr=-100, double hoThr=-100, bool debug=false)
int weight
Definition: histoStyle.py:50
std::vector< int > * t_trackHitInMeasTOB
std::vector< double > * t_PVx
std::vector< int > * t_trackHitInMeasTEC
std::vector< double > * t_esim15x15Photon
std::vector< double > * t_e9x9
std::vector< double > * t_hsim5x5Matched
const L1GtTriggerMenu * ptrL1TriggerMenuEventSetup(int &errorCode)
return a pointer to the L1 trigger menu from event setup
Definition: L1GtUtils.cc:2228
const int prescaleFactor(const edm::Event &iEvent, const edm::InputTag &l1GtRecordInputTag, const edm::InputTag &l1GtReadoutRecordInputTag, const std::string &nameAlgoTechTrig, int &errorCode) const
return prescale factor for a given algorithm or technical trigger
Definition: L1GtUtils.cc:1485
double dxy() const
dxy parameter. (This is the transverse impact parameter w.r.t. to (0,0,0) ONLY if refPoint is close t...
Definition: TrackBase.h:119
std::vector< double > * t_esim15x15Matched
double chargeIsolationHcal(unsigned int trkIndex, std::vector< spr::propagatedTrackID > &vdetIds, const HcalTopology *topology, int ieta, int iphi, bool debug=false)
std::vector< int > * t_trackHitOutMeasTID
std::vector< double > * t_L1METEta
std::vector< double > * t_e7x7_30Sig
std::vector< int > * t_trackHitOutMissTEC
value_type const * get() const
Definition: RefToBase.h:212
std::vector< double > * t_e11x11_25Sig
std::vector< double > * t_hsim3x3
std::vector< double > * t_trackP
static const size_t NPBins
std::vector< double > * t_hsim7x7Photon
std::vector< double > * t_jetPt
std::vector< double > * t_L1IsoEMPt
std::vector< double > * t_esim11x11Photon
edm::InputTag L1GTObjectMapRcdSource_
std::vector< double > * t_h5x5Sig
std::vector< int > * t_trackHitInMissTID
std::vector< double > * t_e11x11_30Sig
std::vector< int > * t_trackHitInMissTIBTID
std::vector< double > * t_esim9x9CharHad
double eECALmatrix(const DetId &detId, edm::Handle< T > &hitsEB, edm::Handle< T > &hitsEE, const CaloGeometry *geo, const CaloTopology *caloTopology, int ieta, int iphi, double ebThr=-100, double eeThr=-100, double tMin=-500, double tMax=500, bool debug=false)
edm::EDGetTokenT< EcalRecHitCollection > tok_EE_
edm::EDGetTokenT< l1extra::L1JetParticleCollection > tok_L1extTauJet_
Definition: DDAxes.h:10
std::vector< double > * t_esim15x15CharHad