CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
myFilter.cc
Go to the documentation of this file.
8 
11 
15 
17 // #include "FWCore/Common/interface/TriggerNames.h"
21 
25 
26 // #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
28 
31 
34 
35 // #include "DataFormats/PhotonReco/interface/PhotonFwd.h"
36 // #include "DataFormats/PhotonReco/interface/Photon.h"
37 
38 
39 // include files
41 
42 
43 using namespace edm;
44 using namespace reco;
45 using namespace std;
46 
47 #define DEBUG false
48 #define INVALID 9999.
49 
50 typedef struct RBX_struct {
51  double et;
52  double hadEnergy;
53  double emEnergy;
54  float hcalTime;
55  float ecalTime;
56  int nTowers;
57 } RBX ;
58 
59 typedef struct HPD_struct {
60  double et;
61  double hadEnergy;
62  double emEnergy;
63  double time;
64  float hcalTime;
65  float ecalTime;
66  int nTowers;
67 } HPD ;
68 
69 
70 
71 
72 //enum HcalSubdetector { HcalEmpty=0, HcalBarrel=1, HcalEndcap=2, HcalOuter=3, HcalForward=4, HcalTriggerTower=5, HcalOther=7 };
73 
74 //enum SectorId { HBplus=1, HBminus=2,
75 // HEplus=3, HEminus=4,
76 // HO2plus=5, HO1plus=6, HOzero=7, HO1minus=8, HO2minus=9,
77 // HFplus=10, HFminus=11 };
78 
79 
81  CaloJetAlgorithm( cfg.getParameter<string>( "CaloJetAlgorithm" ) ),
82  hcalNoiseSummaryTag_(cfg.getParameter<edm::InputTag>("hcalNoiseSummaryTag"))
83 {
84  _nTotal = 0;
85  _nEvent = 0;
86  _acceptedEvt = 0;
87  _passPt = 0;
88  _passNJets = 0;
89  _passDiJet = 0;
90  _passNTrks = 0;
91  _passEMF = 0;
92  _passNTowers = 0;
93  _passMET = 0;
94  _passMETSig = 0;
95  _passHighPtTower = 0;
96  _passNRBX = 0;
97  _passNHPDNoise = 0;
98  _passHLT = 0;
99  _passNPMTHits = 0;
100  _passNMultiPMTHits = 0;
101  _passPKAM = 0;
102  _passHFMET = 0;
103  _passNoiseSummary = 0;
109  _passOERatio = 0;
110  _passTime = 0;
111  _passHBHETime = 0;
112  _passHFTime = 0;
113  _passHFFlagged = 0;
114  _passHFHighEnergy = 0;
115 
116  for (int i=0; i<10; i++) _NoiseResult[i] = 0;
117 
118  theTriggerResultsLabel = cfg.getParameter<edm::InputTag>("TriggerResultsLabel");
119 
120 }
121 
123 }
124 
126 }
127 
129 
130  std::cout << "=============================================================" << std::endl;
131  std::cout << "myFilter: accepted "
132  << _acceptedEvt << " / " << _nEvent << " / " << _nTotal << " events" << std::endl;
133  std::cout << "Pt = " << _passPt << std::endl;
134  std::cout << "NJets = " << _passNJets << std::endl;
135  std::cout << "DiJets = " << _passDiJet << std::endl;
136  std::cout << "NTrks = " << _passNTrks << std::endl;
137  std::cout << "EMF = " << _passEMF << std::endl;
138  std::cout << "NTowers = " << _passNTowers << std::endl;
139  std::cout << "MET = " << _passMET << std::endl;
140  std::cout << "METSig = " << _passMETSig << std::endl;
141  std::cout << "HighPtTower = " << _passHighPtTower << std::endl;
142  std::cout << "NRBX = " << _passNRBX << std::endl;
143  std::cout << "NHPDNoise = " << _passNHPDNoise << std::endl;
144  std::cout << "NPMTHits = " << _passNPMTHits << std::endl;
145  std::cout << "NMultPMTHits = " << _passNMultiPMTHits << std::endl;
146  std::cout << "PKAM = " << _passPKAM << std::endl;
147  std::cout << "HFMET = " << _passHFMET << std::endl;
148  std::cout << "Noise Summary = " << _passNoiseSummary << std::endl;
149  std::cout << "Noise Summary EMF = " << _passNoiseSummaryEMF << std::endl;
150  std::cout << "Noise Summary E2E10 = " << _passNoiseSummaryE2E10 << std::endl;
151  std::cout << "Noise Summary NHITS = " << _passNoiseSummaryNHITS << std::endl;
152  std::cout << "Noise Summary ADC0 = " << _passNoiseSummaryADC0 << std::endl;
153  std::cout << "Noise Summary NoOther = " << _passNoiseSummaryNoOther << std::endl;
154  std::cout << "OERatio = " << _passOERatio << std::endl;
155  std::cout << "Time = " << _passTime << std::endl;
156  std::cout << "HF Time = " << _passHFTime << std::endl;
157  std::cout << "HBHE Time = " << _passHBHETime << std::endl;
158  std::cout << "HF Flagged = " << _passHFFlagged << std::endl;
159  std::cout << "HF High Energy= " << _passHFHighEnergy << std::endl;
160  std::cout << "=============================================================" << std::endl;
161 
162 
163  for (int i=0; i<10; i++) {
164  std::cout << "Noise Results = " << _NoiseResult[i] << std::endl;
165  }
166 
167 }
168 
169 bool
171 
172  double HFRecHit[100][100][2];
173 
174 
175  double HFThreshold = 4.0;
176  //double HOThreshold = 1.0;
177 
178 
179  bool result = false;
180  bool filter_Pt = false;
181  bool filter_DiJet = false;
182  bool filter_NTrks = false;
183  bool filter_EMF = false;
184  bool filter_NJets = false;
185  //bool filter_NTowers = false;
186  bool filter_MET = false;
187  bool filter_METSig = false;
188  bool filter_HighPtTower = false;
189  bool filter_NRBX = false;
190  bool filter_NHPDNoise = false;
191  bool filter_HLT = false;
192  bool filter_NPMTHits = false;
193  bool filter_NMultiPMTHits = false;
194 // bool filter_PKAM = false;
195  bool filter_HFMET = false;
196  bool filter_NoiseSummary = false;
197  bool filter_NoiseSummaryEMF = false;
198  bool filter_NoiseSummaryE2E10 = false;
199  bool filter_NoiseSummaryNHITS = false;
200  bool filter_NoiseSummaryADC0 = false;
201  bool filter_NoiseSummaryNoOther = false;
202  bool filter_OERatio = false;
203  bool filter_Time = false;
204  bool filter_HFTime = false;
205  bool filter_HBHETime = false;
206  bool filter_HFFlagged = false;
207  bool filter_HFHighEnergy = false;
208 
209 
210  bool Pass = false;
211  if (evt.id().run() == 124009) {
212  if ( (evt.bunchCrossing() == 51) ||
213  (evt.bunchCrossing() == 151) ||
214  (evt.bunchCrossing() == 2824) ) {
215  Pass = true;
216  }
217  }
218  if (evt.id().run() == 124020) {
219  if ( (evt.bunchCrossing() == 51) ||
220  (evt.bunchCrossing() == 151) ||
221  (evt.bunchCrossing() == 2824) ) {
222  Pass = true;
223  }
224  }
225  if (evt.id().run() == 124024) {
226  if ( (evt.bunchCrossing() == 51) ||
227  (evt.bunchCrossing() == 151) ||
228  (evt.bunchCrossing() == 2824) ) {
229  Pass = true;
230  }
231  }
232 
233  if ( (evt.bunchCrossing() == 51) ||
234  (evt.bunchCrossing() == 151) ||
235  (evt.bunchCrossing() == 2824) ) {
236  Pass = true;
237  }
238 
239 
240  // ***********************
241  // ***********************
242  // get the Noise summary object
243 
245  evt.getByLabel(hcalNoiseSummaryTag_, summary_h);
246  if(!summary_h.isValid()) {
247  throw edm::Exception(edm::errors::ProductNotFound) << " could not find HcalNoiseSummary.\n";
248  // return true;
249  }
250  const HcalNoiseSummary summary = *summary_h;
251 
252  if(summary.minE2Over10TS()<0.7) {
253  filter_NoiseSummaryE2E10 = true;
254  filter_NoiseSummary = true;
255  _NoiseResult[0]++;
256  }
257  if(summary.maxE2Over10TS()>0.96) {
258  filter_NoiseSummaryE2E10 = true;
259  filter_NoiseSummary = true;
260  _NoiseResult[1]++;
261  }
262  if(summary.maxHPDHits()>=17) {
263  filter_NoiseSummaryNHITS = true;
264  filter_NoiseSummary = true;
265  _NoiseResult[2]++;
266  }
267  if(summary.maxRBXHits()>=999) {
268  filter_NoiseSummary = true;
269  _NoiseResult[3]++;
270  }
271  if(summary.maxHPDNoOtherHits()>=10) {
272  filter_NoiseSummary = true;
273  filter_NoiseSummaryNoOther = true;
274  _NoiseResult[4]++;
275  }
276  if(summary.maxZeros()>=10) {
277  filter_NoiseSummaryADC0 = true;
278  filter_NoiseSummary = true;
279  _NoiseResult[5]++;
280  }
281  if(summary.min25GeVHitTime()<-9999.0) {
282  filter_NoiseSummary = true;
283  _NoiseResult[6]++;
284  }
285  if(summary.max25GeVHitTime()>9999.0) {
286  filter_NoiseSummary = true;
287  _NoiseResult[7]++;
288  }
289  if(summary.minRBXEMF()<0.01) {
290  filter_NoiseSummaryEMF = true;
291  // filter_NoiseSummary = true;
292  _NoiseResult[8]++;
293  }
294 
295  // if (filter_NoiseSummary)
296  // std::cout << ">>> Noise Filter = " << filter_NoiseSummary << std::endl;
297 
298 
299  // summary.passLooseNoiseFilter();
300  // summary.passTightNoiseFilter();
301  // summary.passHighLevelNoiseFilter();
302 
303 
304  // ***********************
305  // ***********************
306 
307  for (int i=0; i<100; i++) {
308  for (int j=0; j<100; j++) {
309  HFRecHit[i][j][0] = -10.;
310  HFRecHit[i][j][1] = -10.;
311  }
312  }
313 
314 
315  double HFM_ETime, HFP_ETime;
316  double HFM_E, HFP_E;
317  double HF_PMM;
318  double MaxRecHitEne;
319 
320  MaxRecHitEne = 0;
321  HFM_ETime = 0.;
322  HFM_E = 0.;
323  HFP_ETime = 0.;
324  HFP_E = 0.;
325  int NPMTHits;
326  int NHFDigiTimeHits;
327  int NHFLongShortHits;
328  int nTime = 0;
329 
330  NPMTHits = 0;
331  NHFDigiTimeHits = 0;
332  NHFLongShortHits = 0;
333 
334  try {
335  std::vector<edm::Handle<HFRecHitCollection> > colls;
336  evt.getManyByType(colls);
337  std::vector<edm::Handle<HFRecHitCollection> >::iterator i;
338  for (i=colls.begin(); i!=colls.end(); i++) {
339  for (HFRecHitCollection::const_iterator j=(*i)->begin(); j!=(*i)->end(); j++) {
340  if (j->id().subdet() == HcalForward) {
341 
342  int myFlag;
343  myFlag= j->flagField(HcalCaloFlagLabels::HFLongShort);
344  if (myFlag==1) {
345  filter_HFFlagged=true;
346  NHFLongShortHits++;
347  }
348  myFlag= j->flagField(HcalCaloFlagLabels::HFDigiTime);
349  if (myFlag==1) {
350  filter_HFFlagged=true;
351  NHFDigiTimeHits++;
352  }
353 
354  if ( ( (j->flagField(HcalCaloFlagLabels::HFLongShort)) == 0) &&
355  ( (j->flagField(HcalCaloFlagLabels::HFDigiTime)) == 0) ) {
356  if (j->energy() > MaxRecHitEne) MaxRecHitEne = j->energy();
357  }
358 
359 
360  // if (filter_HFFlagged) {
361  // std::cout << "HF Flagged = " << _passHFFlagged << std::endl;
362  // }
363 
364 
365  // Long: depth = 1
366  // Short: depth = 2
367  float en = j->energy();
368  float time = j->time();
369  if ( (en > 20.) && (time > 10.)) {
370  nTime++;
371  }
372  int ieta = j->id().ieta();
373  int iphi = j->id().iphi();
374  int depth = j->id().depth();
375  HFRecHit[ieta+41][iphi][depth-1] = en;
376 
377  // Exclude PMTs with crystals
378  if ( (j->id().iphi() == 67) &&
379  ((j->id().ieta() == 29) ||
380  ((j->id().ieta() == 30) && (j->id().depth() == 2)) ||
381  ((j->id().ieta() == 32) && (j->id().depth() == 2)) ||
382  ((j->id().ieta() == 35) && (j->id().depth() == 1)) ||
383  ((j->id().ieta() == 36) && (j->id().depth() == 2)) ||
384  ((j->id().ieta() == 37) && (j->id().depth() == 1)) ||
385  ((j->id().ieta() == 38) && (j->id().depth() == 2)) ) ) {
386  } else {
387  if ( (j->flagField(0) != 0) && (j->energy() > 1.) ) {
388  NPMTHits++;
389  if (NPMTHits > 1) {
390  std::cout << ">>>> MultiHit: Run = " << evt.id().run()
391  << " Event = " << evt.id().event()
392  << " iEta = " << j->id().ieta()
393  << " iPhi = " << j->id().iphi()
394  << " depth = " << j->id().depth()
395  << " flag = " << j->flagField(0)
396  << " hits = " << NPMTHits
397  << " energy = " << j->energy()
398  << std::endl;
399  }
400 
401  }
402  }
403 
404  if (j->id().ieta()<0) {
405  if (j->energy() > HFThreshold) {
406  HFM_ETime += j->energy()*j->time();
407  HFM_E += j->energy();
408  }
409  } else {
410  if (j->energy() > HFThreshold) {
411  HFP_ETime += j->energy()*j->time();
412  HFP_E += j->energy();
413  }
414  }
415 
416 
417  }
418  }
419  break;
420  }
421  } catch (...) {
422  cout << "No HF RecHits." << endl;
423  }
424 
425  if (MaxRecHitEne > 1000.) filter_HFHighEnergy = true;
426 
427  if (nTime > 0) filter_Time = true;
428  if (nTime > 0) filter_HFTime = true;
429 
430 
431  double OER, OddEne, EvenEne;
432  int nOdd, nEven;
433 
434  OER = 0.0;
435  for (int iphi=0; iphi<100; iphi++) {
436  OddEne = EvenEne = 0.;
437  nOdd = 0;
438  nEven = 0;
439  for (int ieta=0; ieta<100; ieta++) {
440  if (HFRecHit[ieta][iphi][0] > 1.0) {
441  if (ieta%2 == 0) {
442  EvenEne += HFRecHit[ieta][iphi][0];
443  nEven++;
444  } else {
445  OddEne += HFRecHit[ieta][iphi][0];
446  nOdd++;
447  }
448  }
449  if (HFRecHit[ieta][iphi][1] > 1.0) {
450  if (ieta%2 == 0) {
451  EvenEne += HFRecHit[ieta][iphi][1];
452  nEven++;
453  } else {
454  OddEne += HFRecHit[ieta][iphi][1];
455  nOdd++;
456  }
457  }
458  }
459  if (((OddEne + EvenEne) > 10.) && (nOdd > 1) && (nEven > 1)) {
460  OER = (OddEne - EvenEne) / (OddEne + EvenEne);
461  }
462  }
463 
464 
465  if (NPMTHits > 0) filter_NPMTHits = true;
466  if (NPMTHits > 1) filter_NMultiPMTHits = true;
467  // cout << "NPMTHits = " << NPMTHits << endl;
468 
469 
470  if ((HFP_E > 0.) && (HFM_E > 0.)) {
471  HF_PMM = (HFP_ETime / HFP_E) - (HFM_ETime / HFM_E);
472  } else {
473  HF_PMM = INVALID;
474  }
475 
476 
477  if (fabs(HF_PMM) < 10.) {
478  Pass = true;
479  } else {
480  Pass = false;
481  }
482 
483  _nTotal++;
484  Pass = true;
485  if (Pass) {
486  /***
487  std::cout << ">>>> FIL: Run = " << evt.id().run()
488  << " Event = " << evt.id().event()
489  << " Bunch Crossing = " << evt.bunchCrossing()
490  << " Orbit Number = " << evt.orbitNumber()
491  << std::endl;
492  ***/
493  // *********************************************************
494  // --- Event Classification
495  // *********************************************************
496 
497  RBX RBXColl[36];
498  HPD HPDColl[144];
499 
500  int evtType = 0;
501 
503  evt.getByLabel( "towerMaker", caloTowers );
504 
505  for (int i=0;i<36;i++) {
506  RBXColl[i].et = 0;
507  RBXColl[i].hadEnergy = 0;
508  RBXColl[i].emEnergy = 0;
509  RBXColl[i].hcalTime = 0;
510  RBXColl[i].ecalTime = 0;
511  RBXColl[i].nTowers = 0;
512  }
513  for (int i=0;i<144;i++) {
514  HPDColl[i].et = 0;
515  HPDColl[i].hadEnergy = 0;
516  HPDColl[i].emEnergy = 0;
517  HPDColl[i].hcalTime = 0;
518  HPDColl[i].ecalTime = 0;
519  HPDColl[i].nTowers = 0;
520  }
521 
522  double HFMET = 0.0;
523  double HFsum_et = 0.0;
524  double HFsum_ex = 0.0;
525  double HFsum_ey = 0.0;
526 
527  for (CaloTowerCollection::const_iterator tower = caloTowers->begin();
528  tower != caloTowers->end(); tower++) {
529 
530  if (tower->hadEnergy() < 0.) {
531  }
532  if (tower->emEnergy() < 0.) {
533  }
534 
535 
536  if (tower->et()>0.5) {
537 
538  int iRBX = tower->iphi();
539  iRBX = iRBX-2;
540  if (iRBX == 0) iRBX = 17;
541  if (iRBX == -1) iRBX = 18;
542  iRBX = (iRBX-1)/4;
543 
544  if (tower->ieta() < 0) iRBX += 18;
545  if (iRBX < 36) {
546  RBXColl[iRBX].et += tower->et();
547  RBXColl[iRBX].hadEnergy += tower->hadEnergy();
548  RBXColl[iRBX].emEnergy += tower->emEnergy();
549  RBXColl[iRBX].hcalTime += tower->hcalTime();
550  RBXColl[iRBX].ecalTime += tower->ecalTime();
551  RBXColl[iRBX].nTowers++;
552  }
553  /***
554  std::cout << "iRBX = " << iRBX << " "
555  << "ieta/iphi = " << tower->ieta() << " / " << tower->iphi()
556  << " et = " << tower->et()
557  << std::endl;
558  ***/
559  int iHPD = tower->iphi();
560  if (tower->ieta() < 0) iHPD = iHPD + 72;
561  if (iHPD < 144) {
562  HPDColl[iHPD].et += tower->et();
563  HPDColl[iHPD].hadEnergy += tower->hadEnergy();
564  HPDColl[iHPD].emEnergy += tower->emEnergy();
565  HPDColl[iHPD].hcalTime += tower->hcalTime();
566  HPDColl[iHPD].ecalTime += tower->ecalTime();
567  HPDColl[iHPD].nTowers++;
568  }
569  /***
570  std::cout << "iHPD = " << iHPD << " "
571  << "ieta/iphi = " << tower->ieta() << " / " << tower->iphi()
572  << " et = " << tower->et()
573  << std::endl;
574  ***/
575 
576  Double_t et = tower->et();
577  Double_t phix = tower->phi();
578 
579  if (fabs(tower->ieta()) > 29) {
580  HFsum_et += et;
581  HFsum_ex += et*cos(phix);
582  HFsum_ey += et*sin(phix);
583  }
584 
585  }
586  }
587 
588  HFMET = sqrt( HFsum_ex*HFsum_ex + HFsum_ey*HFsum_ey);
589  if ( (HFMET > 40.) && (NPMTHits == 0) ) filter_HFMET = true;
590 
591  // Loop over the RBX Collection
592  int nRBX = 0;
593  int nTowers = 0;
594  for (int i=0;i<36;i++) {
595  if (RBXColl[i].hadEnergy > 3.0) {
596  nRBX++;
597  nTowers = RBXColl[i].nTowers;
598  }
599  }
600  if ( (nRBX == 1) && (nTowers > 24) ) {
601  evtType = 1;
602  }
603 
604  // Loop over the HPD Collection
605  int nHPD = 0;
606  for (int i=0;i<144;i++) {
607  if (HPDColl[i].hadEnergy > 3.0) {
608  nHPD++;
609  nTowers = HPDColl[i].nTowers;
610  }
611  }
612  if ( (nHPD == 1) && (nTowers > 6) ) {
613  evtType = 2;
614  // cout << " nHPD = " << nHPD
615  // << " Towers = " << nTowers
616  // << " Type = " << evtType
617  // << endl;
618  }
619 
620 
621 
622  // *********************************************************
623  // --- Access Trigger Info
624  // *********************************************************
625 
626  // **** Get the TriggerResults container
628  evt.getByLabel(theTriggerResultsLabel, triggerResults);
629 
630  Int_t JetLoPass = 0;
631 
632  /****
633  if (triggerResults.isValid()) {
634  if (DEBUG) std::cout << "trigger valid " << std::endl;
635  edm::TriggerNames triggerNames; // TriggerNames class
636  triggerNames.init(*triggerResults);
637  unsigned int n = triggerResults->size();
638  for (unsigned int i=0; i!=n; i++) {
639  if ( triggerNames.triggerName(i) == "HLT_Jet30" ) {
640  JetLoPass = triggerResults->accept(i);
641  if (DEBUG) std::cout << "Found HLT_Jet30" << std::endl;
642  }
643  }
644  }
645  *****/
646 
647  // *********************************************************
648  // --- Vertex Selection
649  // *********************************************************
650 
651  // *********************************************************
652  // --- Pixel Track and Clusters
653  // *********************************************************
654  /*******
655  // -- Tracks
656  edm::Handle<std::vector<reco::Track> > hTrackCollection;
657  try {
658  evt.getByLabel("generalTracks", hTrackCollection);
659  } catch (cms::Exception &ex) {
660  if (fVerbose > 1) cout << "No Track collection with label " << fTrackCollectionLabel << endl;
661  }
662 
663  if (hTrackCollection.isValid()) {
664  const std::vector<reco::Track> trackColl = *(hTrackCollection.product());
665  nTk = trackColl.size();
666  }
667  *******/
668 
669  // *********************************************************
670  // --- Pixel Clusters
671  // *********************************************************
672  // -- Pixel cluster
673  /***
674  edm::Handle<reco::SiPixelCluster> hClusterColl;
675  evt.getByLabel("siPixelClusters", hClusterColl);
676  const reco::SiPixelCluster cC = *(hClusterColl.product());
677  ***/
678 
680  evt.getByLabel("siPixelClusters", hClusterColl);
681  const edmNew::DetSetVector<SiPixelCluster> clustColl = *(hClusterColl.product());
682  // nCl = clustColl.size();
683 
684  /***
685  int nCl = 0;
686  if (hClusterColl.isValid()) {
687  const edmNew::DetSetVector<SiPixelCluster> clustColl = *(hClusterColl.product());
688  nCl = clustColl.size();
689  }
690  ***/
691 
692  // *********************************************************
693  // --- Track Selection
694  // *********************************************************
696  // evt.getByLabel("ctfWithMaterialTracks", trackCollection);
697  evt.getByLabel("generalTracks", trackCollection);
698 
699  const reco::TrackCollection tC = *(trackCollection.product());
700  // std::cout << "FIL: Reconstructed "<< tC.size() << " tracks" << std::endl ;
701 
702  if (tC.size() > 3) filter_NTrks = true;
703 
704  // h_Trk_NTrk->Fill(tC.size());
705 
706  // for (reco::TrackCollection::const_iterator track=tC.begin(); track!=tC.end(); track++){
707  // h_Trk_pt->Fill(track->pt());
708  // }
709 
710  /****
711  std::cout << "Track number "<< i << std::endl ;
712  std::cout << "\tmomentum: " << track->momentum()<< std::endl;
713  std::cout << "\tPT: " << track->pt()<< std::endl;
714  std::cout << "\tvertex: " << track->vertex()<< std::endl;
715  std::cout << "\timpact parameter: " << track->d0()<< std::endl;
716  std::cout << "\tcharge: " << track->charge()<< std::endl;
717  std::cout << "\tnormalizedChi2: " << track->normalizedChi2()<< std::endl;
718 
719  cout<<"\tFrom EXTRA : "<<endl;
720  cout<<"\t\touter PT "<< track->outerPt()<<endl;
721  std::cout << "\t direction: " << track->seedDirection() << std::endl;
722  ****/
723 
724 
725  if ((tC.size() > 100) && (clustColl.size() > 1000)) {
726  _passPKAM++;
727 // filter_PKAM = true;
728  }
729  // std::cout << "N Tracks = " << tC.size()
730  // << " N Cluster = " << clustColl.size() << std::endl ;
731 
732 
733 
734  // *********************************************************
735  // --- RecHits
736  // *********************************************************
737  // Handle<CaloTowerCollection> caloTowers;
738  // evt.getByLabel( "towerMaker", caloTowers );
740 
741  int nHPDNoise = 0;
742  nTime = 0;
743 
744  try {
745  std::vector<edm::Handle<HBHERecHitCollection> > colls;
746  evt.getManyByType(colls);
747  std::vector<edm::Handle<HBHERecHitCollection> >::iterator i;
748  for (i=colls.begin(); i!=colls.end(); i++) {
749  for (HBHERecHitCollection::const_iterator j=(*i)->begin(); j!=(*i)->end(); j++) {
750  // std::cout << *j << std::endl;
751  if (j->id().subdet() == HcalBarrel) {
752  // std::cout << "Barrel : " << j->id() << std::endl;
753  }
754  if (j->id().subdet() == HcalEndcap) {
755  }
756 
757  if (j->flagField(0) != 0) nHPDNoise++;
758 
759  float en = j->energy();
760  float time = j->time();
761 
762  if ( (en > 10.) && (time > 20.)) {
763  nTime++;
764  }
765 
766  /***
767  std::cout << j->id() << " "
768  << j->id().subdet() << " "
769  << j->id().ieta() << " "
770  << j->id().iphi() << " "
771  << j->id().depth() << " "
772  << j->energy() << " "
773  << j->time() << std::endl;
774  ****/
775  }
776  }
777  } catch (...) {
778  cout << "No HB/HE RecHits." << endl;
779  }
780 
781  if (nHPDNoise > 10) filter_NHPDNoise = true;
782  if (nTime > 0) filter_HBHETime = true;
783 
784 
785  try {
786  std::vector<edm::Handle<HFRecHitCollection> > colls;
787  evt.getManyByType(colls);
788  std::vector<edm::Handle<HFRecHitCollection> >::iterator i;
789  for (i=colls.begin(); i!=colls.end(); i++) {
790  for (HFRecHitCollection::const_iterator j=(*i)->begin(); j!=(*i)->end(); j++) {
791  // std::cout << *j << std::endl;
792  }
793  }
794  } catch (...) {
795  cout << "No HF RecHits." << endl;
796  }
797 
798  try {
799  std::vector<edm::Handle<HORecHitCollection> > colls;
800  evt.getManyByType(colls);
801  std::vector<edm::Handle<HORecHitCollection> >::iterator i;
802  for (i=colls.begin(); i!=colls.end(); i++) {
803  for (HORecHitCollection::const_iterator j=(*i)->begin(); j!=(*i)->end(); j++) {
804  // std::cout << *j << std::endl;
805  }
806  }
807  } catch (...) {
808  cout << "No HO RecHits." << endl;
809  }
810 
811 
812 
813  // *********************************************************
814  // --- CaloTower Selection
815  // *********************************************************
816  // Handle<CaloTowerCollection> caloTowers;
817  // evt.getByLabel( "towerMaker", caloTowers );
818 
819  // --- Loop over towers and make a lists of used and unused towers
820  int nTow = 0;
821  for (CaloTowerCollection::const_iterator tower = caloTowers->begin();
822  tower != caloTowers->end(); tower++) {
823  // std::cout << *tower << std::endl;
824  if (tower->et() > 0.5) {
825  nTow++;
826  /****
827  std::cout << "Tower Et = "
828  << tower->et() << " "
829  << tower->emEnergy() << " EmEne = "
830  << tower->hadEnergy() << " HadEne = "
831  << tower->outerEnergy() << " ETime = "
832  << tower->ecalTime() << " HTime = "
833  << tower->hcalTime() << " ieta = "
834  << tower->ieta() << " iphi = "
835  << tower->iphi() << " "
836  << tower->iphi() / 4
837  << endl;
838  ****/
839  }
840 
841  }
842  /****
843  std::cout << "Number of caloTowers = "
844  << caloTowers->size()
845  << " / "
846  << nTow
847  << std::endl;
848  ****/
849 
850  // *********************************************************
851  // --- Jet Selection
852  // *********************************************************
853 
855  evt.getByLabel( CaloJetAlgorithm, jets );
856  int njet = 0;
857  int nDiJet = 0;
858  for ( CaloJetCollection::const_iterator ijet=jets->begin(); ijet!=jets->end(); ijet++) {
859 
860  if ( (ijet->pt() > 100.) && (JetLoPass != 0) ) {
861  filter_HLT = true;
862  }
863 
864 
865  if ((ijet->pt() > 100.) && (evtType == 0)) {
866  filter_HighPtTower = true;
867  }
868 
869  if (ijet->pt() > 50.) nDiJet++;
870  if (ijet->pt() > 50.) filter_Pt = true;
871  if (ijet->pt() > 10.) njet++;
872  if (ijet->pt() > 10.) {
873  if (ijet->emEnergyFraction() > 0.05) filter_EMF = true;
874  }
875 
876  // if (filter_EMF) {
877  // std::cout << "pt = " << ijet->pt()
878  // << " EMF = " << ijet->emEnergyFraction() << std::endl;
879  // }
880 
881  // std::cout << "pt = " << ijet->pt() << std::endl;
882 
883  }
884 
885  if (nDiJet > 1) filter_DiJet = true;
886  if (njet > 1) filter_NJets = true;
887  // if (filter_EMF) {
888  // std::cout << "NJets = " << njet << std::endl;
889  // }
890 
891 
892  // *********************************************************
893  // --- MET Selection
894  // *********************************************************
895 
897  evt.getByLabel("met", calometcoll);
898  double caloMET = 0;
899  if (calometcoll.isValid()) {
900  const CaloMETCollection *calometcol = calometcoll.product();
901  const CaloMET *calomet;
902  calomet = &(calometcol->front());
903  caloMET = calomet->pt();
904  //double caloMETSig = calomet->mEtSig();
905  //double caloSumET = calomet->sumEt();
906  // if ((caloMET > 50.) && (evtType = 0)) filter_MET = true;
907  if (caloMET > 40.) filter_MET = true;
908  }
909  if ((std::abs(OER) > 0.9) && (caloMET > 20.0)) filter_OERatio = true;
910  if (nRBX > 3) filter_NRBX = true;
911 
912  // *********************************************************
913  _nEvent++;
914 
915  // if ( (filter_HLT) || (filter_NJets) ) {
916  // result = true;
917  // _acceptedEvt++;
918  // }
919 
920  /***
921  if ( (filter_Pt) || (filter_NTrks) || (filter_EMF) || (filter_NJets) ||
922  (filter_MET) || (filter_METSig) || (filter_HighPtTower) ) {
923  result = true;
924  _acceptedEvt++;
925  }
926  ***/
927 
928  // if ( (filter_Pt) || (filter_NJets) ) {
929  // result = true;
930  // _acceptedEvt++;
931  // }
932 
933  // if ((filter_PKAM) || (filter_HFMET) ||(filter_NMultiPMTHits) ) {
934  // if ( (filter_DiJet) || (filter_HFMET) ||(filter_NMultiPMTHits) ) {
935  // if ( (filter_NHPDNoise) && ( (filter_Pt) || (filter_MET) ) ) {
936  // result = true;
937  // _acceptedEvt++;
938  // }
939 
940  // if ( (filter_NoiseSummary) && ( (filter_Pt) || (filter_MET) ) ) {
941  // result = true;
942  // _acceptedEvt++;
943  // }
944 
945  // if (filter_NoiseSummaryEMF) {
946  // result = true;
947  // _acceptedEvt++;
948  // }
949 
950  // if (filter_Time) {
951  // if ( (filter_MET) && (!filter_NoiseSummary) ) {
952 
953  // if (filter_NoiseSummary) {
954  // result = true;
955  // _acceptedEvt++;
956  // }
957 
958  // if (filter_HFTime) {
959  // result = true;
960  // _acceptedEvt++;
961  // }
962 
963  // if (filter_HBHETime) {
964  // result = true;
965  // _acceptedEvt++;
966  // }
967 
968 
969  // if (!filter_NoiseSummary && filter_HBHETime) {
970  // result = true;
971  // _acceptedEvt++;
972  // }
973 
974  // if (filter_NoiseSummary && filter_NTrks) {
975  // result = true;
976  // _acceptedEvt++;
977  // }
978 
979  // if (!filter_NoiseSummaryADC0 && !filter_NoiseSummaryNHITS &&
980  // filter_NoiseSummaryE2E10 && !filter_NoiseSummaryNoOther) {
981  // result = true;
982  // _acceptedEvt++;
983  // }
984 
985 
986  // if ((filter_HFFlagged) && ((NHFLongShortHits > 2) || (NHFDigiTimeHits > 2))) {
987  // result = true;
988  // _acceptedEvt++;
989  // }
990 
991  if (filter_HFHighEnergy) {
992  result = true;
993  _acceptedEvt++;
994  }
995 
996  if (filter_Pt) _passPt++;
997  if (filter_NJets) _passNJets++;
998  if (filter_DiJet) _passDiJet++;
999  if (filter_NTrks) _passNTrks++;
1000  if (filter_EMF) _passEMF++;
1001  if (filter_MET) _passMET++;
1002  if (filter_METSig) _passMETSig++;
1003  if (filter_HighPtTower) _passHighPtTower++;
1004  if (filter_NRBX) _passNRBX++;
1005  if (filter_NHPDNoise) _passNHPDNoise++;
1006  if (filter_HLT) _passHLT++;
1007  if (filter_NPMTHits) _passNPMTHits++;
1008  if (filter_NMultiPMTHits) _passNMultiPMTHits++;
1009  if (filter_HFMET) _passHFMET++;
1010  if (filter_NoiseSummary) _passNoiseSummary++;
1011  if (filter_NoiseSummaryEMF) _passNoiseSummaryEMF++;
1012  if (filter_NoiseSummaryE2E10) _passNoiseSummaryE2E10++;
1013  if (filter_NoiseSummaryNHITS) _passNoiseSummaryNHITS++;
1014  if (filter_NoiseSummaryADC0) _passNoiseSummaryADC0++;
1015  if (filter_NoiseSummaryNoOther) _passNoiseSummaryNoOther++;
1016  if (filter_OERatio) _passOERatio++;
1017  if (filter_Time) _passTime++;
1018  if (filter_HFTime) _passHFTime++;
1019  if (filter_HBHETime) _passHBHETime++;
1020  if (filter_HFFlagged) _passHFFlagged++;
1021  if (filter_HFHighEnergy) _passHFHighEnergy++;
1022 
1023  /****
1024  if ((evt.id().run() == 120020) && (evt.id().event() == 453)) {
1025  result = true;
1026  _acceptedEvt++;
1027  } else {
1028  result = false;
1029  }
1030  ****/
1031 
1032  // result = true;
1033 
1034  }
1035 
1036  if (result) {
1037  std::cout << "<<<< Event Passed"
1038  << std::endl;
1039  }
1040  return result;
1041 
1042 }
int _passHBHETime
Definition: myFilter.h:53
RunNumber_t run() const
Definition: EventID.h:39
void getManyByType(std::vector< Handle< PROD > > &results) const
Definition: Event.h:454
T getParameter(std::string const &) const
EventNumber_t event() const
Definition: EventID.h:41
int i
Definition: DBlmapReader.cc:9
double emEnergy
Definition: myFilter.cc:62
int _nTotal
Definition: myFilter.h:24
int _passNHPDNoise
Definition: myFilter.h:38
tuple cfg
Definition: looper.py:293
double et
Definition: myFilter.cc:60
int _passNoiseSummaryADC0
Definition: myFilter.h:48
int _NoiseResult[10]
Definition: myFilter.h:57
Definition: hltDiff.cc:316
virtual ~myFilter()
Definition: myFilter.cc:122
int _acceptedEvt
Definition: myFilter.h:26
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
double emEnergy
Definition: myFilter.cc:53
int _passNRBX
Definition: myFilter.h:37
std::vector< Track > TrackCollection
collection of Tracks
Definition: TrackFwd.h:14
std::vector< HFRecHit >::const_iterator const_iterator
float maxE2Over10TS(void) const
int bunchCrossing() const
Definition: EventBase.h:66
int _passNMultiPMTHits
Definition: myFilter.h:41
double hadEnergy
Definition: myFilter.cc:61
int nTowers
Definition: myFilter.cc:56
HPD HPDColl[144]
Definition: myJetAna.cc:116
int njet
Definition: HydjetWrapper.h:95
float min25GeVHitTime(void) const
virtual void beginJob()
Definition: myFilter.cc:125
virtual double pt() const
transverse momentum
int _passPt
Definition: myFilter.h:28
double et
Definition: myFilter.cc:51
Collection of Calo MET.
int _passOERatio
Definition: myFilter.h:50
struct HPD_struct HPD
int _passNoiseSummaryE2E10
Definition: myFilter.h:46
int _passHighPtTower
Definition: myFilter.h:36
int nTowers
Definition: myFilter.cc:66
RBX RBXColl[36]
Definition: myJetAna.cc:115
int _passHFTime
Definition: myFilter.h:52
float minE2Over10TS(void) const
int _passHFHighEnergy
Definition: myFilter.h:55
int _passHLT
Definition: myFilter.h:39
int _passTime
Definition: myFilter.h:51
virtual bool filter(edm::Event &e, edm::EventSetup const &c)
Definition: myFilter.cc:170
int maxRBXHits(void) const
T sqrt(T t)
Definition: SSEVec.h:48
vector< PseudoJet > jets
#define INVALID
Definition: myFilter.cc:48
tuple result
Definition: query.py:137
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
float ecalTime
Definition: myFilter.cc:65
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
int j
Definition: DBlmapReader.cc:9
static std::string const triggerResults
Definition: EdmProvDump.cc:40
bool isValid() const
Definition: HandleBase.h:75
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:420
int _passHFMET
Definition: myFilter.h:43
virtual void endJob()
Definition: myFilter.cc:128
edm::InputTag theTriggerResultsLabel
Definition: myFilter.h:60
float hcalTime
Definition: myFilter.cc:54
int _passMET
Definition: myFilter.h:34
struct RBX_struct RBX
int _passPKAM
Definition: myFilter.h:42
T const * product() const
Definition: Handle.h:81
int _passMETSig
Definition: myFilter.h:35
int _passNTrks
Definition: myFilter.h:29
int maxZeros(void) const
int _passDiJet
Definition: myFilter.h:32
int _passNPMTHits
Definition: myFilter.h:40
double hadEnergy
Definition: myFilter.cc:52
edm::InputTag hcalNoiseSummaryTag_
Definition: myFilter.h:61
size_type size() const
int _passNoiseSummaryNHITS
Definition: myFilter.h:47
edm::EventID id() const
Definition: EventBase.h:60
int maxHPDNoOtherHits(void) const
float hcalTime
Definition: myFilter.cc:64
int _passHFFlagged
Definition: myFilter.h:54
float ecalTime
Definition: myFilter.cc:55
int _passNoiseSummaryNoOther
Definition: myFilter.h:49
float minRBXEMF(void) const
std::string CaloJetAlgorithm
Definition: myFilter.h:59
int _nEvent
Definition: myFilter.h:23
myFilter(const edm::ParameterSet &)
Definition: myFilter.cc:80
tuple cout
Definition: gather_cfg.py:121
int _passNoiseSummaryEMF
Definition: myFilter.h:45
double time
Definition: myFilter.cc:63
int _passNoiseSummary
Definition: myFilter.h:44
int _passNTowers
Definition: myFilter.h:33
int maxHPDHits(void) const
float max25GeVHitTime(void) const
int _passEMF
Definition: myFilter.h:30
int _passNJets
Definition: myFilter.h:31