CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
HOCalibAnalyzer.cc
Go to the documentation of this file.
1 //April 2015 : Removal of itrg1, itrg2, but addition of isect2, same is true in HOCalibVariables.h
2 // -*- C++ -*-
3 //
4 // Package: HOCalibAnalyzer
5 // Class: HOCalibAnalyzer
6 //
19 //
20 // Original Author: Gobinda Majumder
21 // Created: Sat Jul 7 09:51:31 CEST 2007
22 //
23 //
24 
25 // system include files
26 #include <memory>
27 
28 // user include files
29 #include "FWCore/Framework/interface/Frameworkfwd.h"
30 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
31 
32 #include "FWCore/Framework/interface/Event.h"
33 #include "FWCore/Framework/interface/MakerMacros.h"
34 #include "FWCore/MessageLogger/interface/MessageLogger.h"
35 #include "FWCore/ParameterSet/interface/ParameterSet.h"
36 
37 #include "FWCore/ServiceRegistry/interface/Service.h"
38 #include "FWCore/Utilities/interface/InputTag.h"
39 #include "CommonTools/UtilAlgos/interface/TFileService.h"
40 #include "DataFormats/HcalCalibObjects/interface/HOCalibVariables.h"
41 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
42 
43 #include "TFile.h"
44 #include "TH1F.h"
45 #include "TH2F.h"
46 #include "TTree.h"
47 #include "TProfile.h"
48 
49 #include <string>
50 
51 #include <iostream>
52 #include <fstream>
53 #include <iomanip>
54 
55 using namespace std;
56 using namespace edm;
57 
58 //#define EDM_ML_DEBUG
59 //
60 // class decleration
61 //
62 
63 class HOCalibAnalyzer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
64 public:
65  explicit HOCalibAnalyzer(const edm::ParameterSet&);
66  ~HOCalibAnalyzer() override;
67 
68 private:
69  void beginJob() override {}
70  void analyze(const edm::Event&, const edm::EventSetup&) override;
71  void endJob() override {}
72 
73  static constexpr int netamx = 30;
74  static constexpr int nphimx = 72;
75  static constexpr int ringmx = 5;
76  static constexpr int ncut = 14;
77 
78  const char* varcrit[3] = {"All", "steps", "n-1"}; // or opposite
79 
80  const double elosfact = (14.9 + 0.96 * fabs(log(8 * 2.8)) + 0.033 * 8 * (1.0 - pow(8, -0.33)));
81 
82  int getHOieta(int ij) { return (ij < netamx / 2) ? -netamx / 2 + ij : -netamx / 2 + ij + 1; }
83  int invert_HOieta(int ieta) { return (ieta < 0) ? netamx / 2 + ieta : netamx / 2 + ieta - 1; }
84 
85  int mypow_2[31];
86 
87  bool m_cosmic;
88  bool m_zeroField;
89  int m_bins;
90  double m_low;
91  double m_ahigh;
92  bool m_histFill; //Stored signals of individual HO towers with default selection criteria
93  bool m_treeFill; //Store rootuple without almost any selection criteria except a quality on muon
94 
95  int ipass;
96 
97  TTree* T1;
98 
99  TH1F* ho_indenergy[netamx][nphimx];
100 
101  TH1F* muonnm;
102  TH1F* muonmm;
103  TH1F* muonth;
104  TH1F* muonph;
105  TH1F* muonch;
106 
107  TH1F* sel_muonnm;
108  TH1F* sel_muonmm;
109  TH1F* sel_muonth;
110  TH1F* sel_muonph;
111  TH1F* sel_muonch;
112 
113  // TProfile* sigvsevt[15][ncut];
114  TH2F* sig_eta_evt[3 * netamx][ncut]; //For individual eta
115  TH2F* sigvsevt[3 * netamx][ncut];
116  TH1F* variab[3 * netamx][ncut];
117 
118  TH2F* mu_projection[ncut + 1];
119 
120  unsigned ievt, hoflag;
121  int irun, ilumi, nprim, isect, isect2, ndof, nmuon;
122 
123  float inslumi, trkdr, trkdz, trkvx, trkvy, trkvz, trkmm, trkth, trkph, chisq, therr, pherr, hodx, hody, hoang, htime,
124  hosig[9], hocorsig[18], hocro, hbhesig[9], caloen[3];
125  float momatho, tkpt03, ecal03, hcal03;
126  float tmphoang;
127 
128  int nevents[10];
129 
130  float ncount[ringmx][ncut + 10];
131 
132  edm::InputTag hoCalibVariableCollectionTag;
133  edm::EDGetTokenT<HOCalibVariableCollection> tok_ho_;
134  edm::EDGetTokenT<HORecHitCollection> tok_allho_;
135  // ----------member data ---------------------------
136 };
137 
138 //
139 // constants, enums and typedefs
140 //
141 
142 //
143 // static data member definitions
144 //
145 
146 //
147 // constructors and destructor
148 //
149 
150 HOCalibAnalyzer::HOCalibAnalyzer(const edm::ParameterSet& iConfig) {
151  // It is very likely you want the following in your configuration
152  // hoCalibVariableCollectionTag = cms.InputTag('hoCalibProducer', 'HOCalibVariableCollection')
153 
154  usesResource(TFileService::kSharedResource);
155 
156  tok_ho_ = consumes<HOCalibVariableCollection>(iConfig.getParameter<edm::InputTag>("hoCalibVariableCollectionTag"));
157  tok_allho_ = consumes<HORecHitCollection>(iConfig.getParameter<edm::InputTag>("hoInputTag"));
158  //now do what ever initialization is needed
159  ipass = 0;
160  for (int ij = 0; ij < 10; ij++) {
161  nevents[ij] = 0;
162  }
163 
164  m_cosmic = iConfig.getUntrackedParameter<bool>("cosmic", true);
165  m_zeroField = iConfig.getUntrackedParameter<bool>("zeroField", false);
166 
167  m_bins = iConfig.getUntrackedParameter<int>("HOSignalBins", 120);
168  m_low = iConfig.getUntrackedParameter<double>("lowerRange", -1.0);
169  m_ahigh = iConfig.getUntrackedParameter<double>("upperRange", 29.0);
170 
171  m_histFill = iConfig.getUntrackedParameter<bool>("histFill", true);
172  m_treeFill = iConfig.getUntrackedParameter<bool>("treeFill", false);
173 
174  edm::Service<TFileService> fs;
175 
176  T1 = fs->make<TTree>("T1", "HOSignal");
177 
178  T1->Branch("irun", &irun, "irun/I");
179  T1->Branch("ievt", &ievt, "ievt/i");
180 
181  T1->Branch("isect", &isect, "isect/I");
182  T1->Branch("isect2", &isect2, "isect2/I");
183  T1->Branch("ndof", &ndof, "ndof/I");
184  T1->Branch("nmuon", &nmuon, "nmuon/I");
185 
186  T1->Branch("ilumi", &ilumi, "ilumi/I");
187  if (!m_cosmic) {
188  T1->Branch("inslumi", &inslumi, "inslumi/F");
189  T1->Branch("nprim", &nprim, "nprim/I");
190  T1->Branch("tkpt03", &tkpt03, " tkpt03/F");
191  T1->Branch("ecal03", &ecal03, " ecal03/F");
192  T1->Branch("hcal03", &hcal03, " hcal03/F");
193  }
194 
195  T1->Branch("trkdr", &trkdr, "trkdr/F");
196  T1->Branch("trkdz", &trkdz, "trkdz/F");
197 
198  T1->Branch("trkvx", &trkvx, "trkvx/F");
199  T1->Branch("trkvy", &trkvy, "trkvy/F");
200  T1->Branch("trkvz", &trkvz, "trkvz/F");
201  T1->Branch("trkmm", &trkmm, "trkmm/F");
202  T1->Branch("trkth", &trkth, "trkth/F");
203  T1->Branch("trkph", &trkph, "trkph/F");
204 
205  T1->Branch("chisq", &chisq, "chisq/F");
206  T1->Branch("therr", &therr, "therr/F");
207  T1->Branch("pherr", &pherr, "pherr/F");
208  T1->Branch("hodx", &hodx, "hodx/F");
209  T1->Branch("hody", &hody, "hody/F");
210  T1->Branch("hoang", &hoang, "hoang/F");
211 
212  T1->Branch("momatho", &momatho, "momatho/F");
213  T1->Branch("hoflag", &hoflag, "hoflag/i");
214  T1->Branch("htime", &htime, "htime/F");
215  T1->Branch("hosig", hosig, "hosig[9]/F");
216  T1->Branch("hocro", &hocro, "hocro/F");
217  T1->Branch("hocorsig", hocorsig, "hocorsig[18]/F");
218  T1->Branch("caloen", caloen, "caloen[3]/F");
219 
220  char name[200];
221  char title[200];
222 
223  if (m_histFill) {
224  for (int ij = 0; ij < netamx; ij++) {
225  int ieta = getHOieta(ij);
226  for (int jk = 0; jk < nphimx; jk++) {
227  sprintf(name, "ho_indenergy_%i_%i", ij, jk);
228  sprintf(title, "ho IndEnergy (GeV) i#eta=%i i#phi=%i", ieta, jk + 1);
229  ho_indenergy[ij][jk] = fs->make<TH1F>(name, title, 1200, m_low, m_ahigh);
230  }
231  }
232  }
233 
234  muonnm = fs->make<TH1F>("muonnm", "No of muon", 10, -0.5, 9.5);
235  muonmm = fs->make<TH1F>("muonmm", "P_{mu}", 200, -100., 100.);
236  muonth = fs->make<TH1F>("muonth", "{Theta}_{mu}", 180, 0., 180.);
237  muonph = fs->make<TH1F>("muonph", "{Phi}_{mu}", 180, -180., 180.);
238  muonch = fs->make<TH1F>("muonch", "{chi^2}/ndf", 100, 0., 1000.);
239 
240  sel_muonnm = fs->make<TH1F>("sel_muonnm", "No of muon(sel)", 10, -0.5, 9.5);
241  sel_muonmm = fs->make<TH1F>("sel_muonmm", "P_{mu}(sel)", 200, -100., 100.);
242  sel_muonth = fs->make<TH1F>("sel_muonth", "{Theta}_{mu}(sel)", 180, 0., 180.);
243  sel_muonph = fs->make<TH1F>("sel_muonph", "{Phi}_{mu}(sel)", 180, -180., 180.);
244  sel_muonch = fs->make<TH1F>("sel_muonch", "{chi^2}/ndf(sel)", 100, 0., 1000.);
245 
246  float pival = acos(-1.);
247 
248  //if change order, change in iselect_wotime also and other efficiency numbers
249  const char* varnam[ncut] = {"ndof",
250  "chisq",
251  "th",
252  "ph",
253  "therr",
254  "pherr",
255  "dircos",
256  "trkmm",
257  "nmuon",
258  "calo",
259  "trkiso",
260  "#phi-dir",
261  "#eta-dir",
262  "time"};
263  int nbinxx[ncut] = {25, 60, 60, 60, 60, 60, 60, 120, 6, 60, 60, 120, 120, 60};
264  double alowxx[ncut] = {5.5, 0., 0., -pival, 0.0, 0.0, 0.0, 0., 0.5, 0.0, 0.0, -20., -32., -45.0};
265  double ahghxx[ncut] = {30.5, 40., pival, pival, 0.8, 0.02, 0.5, 300., 6.5, 10.0, 24.0, 20.0, 32.0, 45.0};
266 
267  for (int kl = 0; kl < ncut; kl++) {
268  for (int jk = 0; jk < 3; jk++) {
269  for (int ij = 0; ij < netamx; ij++) {
270  sprintf(name, "sigeta_%i_%i_%i", kl, jk, ij);
271  sprintf(title, "sigeta %s %s i#eta=%i", varnam[kl], varcrit[jk], getHOieta(ij));
272  sig_eta_evt[netamx * jk + ij][kl] =
273  fs->make<TH2F>(name, title, nbinxx[kl], alowxx[kl], ahghxx[kl], m_bins, m_low, m_ahigh);
274  }
275  }
276  }
277 
278  for (int kl = 0; kl < ncut; kl++) {
279  for (int ij = 0; ij < ringmx * 3; ij++) {
280  int iring = ij % ringmx - 2;
281  int iset = ij / ringmx;
282  sprintf(name, "sigring_%i_%i", kl, ij);
283  sprintf(title, "Signal %s %s Ring%i", varnam[kl], varcrit[iset], iring);
284  sigvsevt[ij][kl] = fs->make<TH2F>(name, title, nbinxx[kl], alowxx[kl], ahghxx[kl], m_bins, m_low, m_ahigh);
285  }
286  }
287 
288  for (int kl = 0; kl < ncut; kl++) {
289  for (int ij = 0; ij < ringmx * 3; ij++) {
290  int iring = ij % ringmx - 2;
291  int iset = ij / ringmx;
292  sprintf(name, "varring_%i_%i", kl, ij);
293  sprintf(title, "%s %s Ring%i", varnam[kl], varcrit[iset], iring);
294  variab[ij][kl] = fs->make<TH1F>(name, title, nbinxx[kl], alowxx[kl], ahghxx[kl]);
295  }
296  }
297 
298  for (int ij = 0; ij <= ncut; ij++) {
299  sprintf(name, "mu_projection_%i", ij);
300  if (ij == 0) {
301  sprintf(title, "All projected muon");
302  } else {
303  sprintf(title, "Projected muon with selection %s", varnam[ij - 1]);
304  }
305  mu_projection[ij] =
306  fs->make<TH2F>(name, title, netamx + 1, -netamx / 2 - 0.5, netamx / 2 + 0.5, nphimx, 0.5, nphimx + 0.5);
307  }
308 
309  for (int ij = 0; ij < 31; ij++) {
310  mypow_2[ij] = pow(2, ij);
311  }
312  for (int ij = 0; ij < ringmx; ij++) {
313  for (int jk = 0; jk < ncut + 10; jk++) {
314  ncount[ij][jk] = 0.0;
315  }
316  }
317 }
318 
319 HOCalibAnalyzer::~HOCalibAnalyzer() {
320  // do anything here that needs to be done at desctruction time
321  // (e.g. close files, deallocate resources etc.)
322 
323  edm::LogVerbatim("HOCalib") << " Total events = " << setw(7) << nevents[0] << " " << setw(7) << nevents[1] << " "
324  << setw(7) << nevents[2] << " " << setw(7) << nevents[3] << " " << setw(7) << nevents[4]
325  << " " << setw(7) << nevents[5] << " Selected events # is " << ipass;
326 }
327 
328 //
329 // member functions
330 //
331 
332 // ------------ method called to for each event ------------
333 void HOCalibAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
334  nevents[0]++;
335 
336  using namespace edm;
337 
338  float pival = acos(-1.);
339 
340  ievt = iEvent.id().event();
341  ilumi = iEvent.luminosityBlock();
342 
343  edm::Handle<HOCalibVariableCollection> HOCalib;
344 
345  iEvent.getByToken(tok_ho_, HOCalib);
346 
347  if (nevents[0] % 20000 == 1) {
348  edm::LogVerbatim("HOCalib") << "nmuon event # " << setw(7) << nevents[0] << " " << setw(7) << nevents[1] << " "
349  << setw(7) << nevents[2] << " " << setw(7) << nevents[3] << " " << setw(7) << nevents[4]
350  << " " << setw(7) << nevents[5];
351  edm::LogVerbatim("HOCalib") << " Run # " << iEvent.id().run() << " Evt # " << iEvent.id().event() << " "
352  << int(HOCalib.isValid()) << " " << ipass;
353  }
354 
355  if (HOCalib.isValid()) {
356  nevents[1]++;
357  nmuon = (*HOCalib).size();
358 
359  for (HOCalibVariableCollection::const_iterator hoC = (*HOCalib).begin(); hoC != (*HOCalib).end(); hoC++) {
360  trkdr = (*hoC).trkdr;
361  trkdz = (*hoC).trkdz;
362 
363  trkvx = (*hoC).trkvx;
364  trkvy = (*hoC).trkvy;
365  trkvz = (*hoC).trkvz;
366 
367  trkmm = (*hoC).trkmm;
368  trkth = (*hoC).trkth;
369  trkph = (*hoC).trkph;
370 
371  ndof = (int)(*hoC).ndof;
372  chisq = (*hoC).chisq;
373  momatho = (*hoC).momatho;
374 
375  therr = (*hoC).therr;
376  pherr = (*hoC).pherr;
377  trkph = (*hoC).trkph;
378 
379  if (!m_cosmic) {
380  nprim = (*hoC).nprim;
381  inslumi = (*hoC).inslumi;
382  tkpt03 = (*hoC).tkpt03;
383  ecal03 = (*hoC).ecal03;
384  hcal03 = (*hoC).hcal03;
385  }
386 
387  isect = (*hoC).isect;
388  isect2 = (*hoC).isect2;
389  hodx = (*hoC).hodx;
390  hody = (*hoC).hody;
391  hoang = (*hoC).hoang;
392 
393  tmphoang = sin(trkth) - hoang;
394 
395  htime = (*hoC).htime;
396  hoflag = (*hoC).hoflag;
397  for (int ij = 0; ij < 9; ij++) {
398  hosig[ij] = (*hoC).hosig[ij];
399 #ifdef EDM_ML_DEBUG
400  edm::LogVerbatim("HOCalib") << "hosig " << ij << " " << hosig[ij];
401 #endif
402  }
403  for (int ij = 0; ij < 18; ij++) {
404  hocorsig[ij] = (*hoC).hocorsig[ij];
405 #ifdef EDM_ML_DEBUG
406  edm::LogVerbatim("HOCalib") << "hocorsig " << ij << " " << hocorsig[ij];
407 #endif
408  }
409  hocro = (*hoC).hocro;
410  for (int ij = 0; ij < 3; ij++) {
411  caloen[ij] = (*hoC).caloen[ij];
412  }
413 
414  int ipsall = 0;
415  int ips0 = 0;
416  int ips1 = 0;
417  int ips2 = 0;
418  int ips3 = 0;
419  int ips4 = 0;
420  int ips5 = 0;
421  int ips6 = 0;
422  int ips7 = 0;
423  int ips8 = 0;
424  int ips9 = 0;
425  int ips10 = 0;
426  int ips11 = 0;
427  int ips12 = 0;
428  int ips13 = 0;
429 
430  nevents[2]++;
431  bool isZSps = (hosig[4] < -99.0) ? false : true;
432 
433  if ((!m_cosmic) && fabs(trkmm) < momatho)
434  continue;
435 
436  nevents[3]++;
437  if (fabs(trkth - pival / 2) < 0.000001)
438  continue; //22OCT07
439  nevents[4]++;
440 
441  int ieta = int((abs(isect) % 10000) / 100.) - 50; //an offset to acodate -ve eta values
442  if (abs(ieta) >= 16)
443  continue;
444  nevents[5]++;
445  int iphi = abs(isect) % 100;
446 
447  int tmpsect = int((iphi + 1) / 6.) + 1;
448  if (tmpsect > 12)
449  tmpsect = 1;
450 
451  int iring = 0;
452 
453  int iring2 = iring + 2;
454 
455  double abshoang = (m_cosmic) ? fabs(hoang) : hoang;
456 
457  double elos = 1.0 / TMath::Max(0.1, abs(1.0 * hoang));
458 
459  if (!m_zeroField)
460  elos *= ((14.9 + 0.96 * fabs(log(momatho * 2.8)) + 0.033 * momatho * (1.0 - pow(momatho, -0.33))) / elosfact);
461 
462  if (m_cosmic) {
463  if (abs(ndof) >= 20 && abs(ndof) < 55) {
464  ips0 = mypow_2[0];
465  ipsall += ips0;
466  }
467  if (chisq > 0 && chisq < 12) {
468  ips1 = mypow_2[1];
469  ipsall += ips1;
470  } //18Jan2008
471 
472  if (trkth > 0.3 && trkth < pival - 0.3) {
473  ips2 = mypow_2[2];
474  ipsall += ips2;
475  } //No nead for pp evt
476  if (trkph > -pival + 0.1 && trkph < -0.1) {
477  ips3 = mypow_2[3];
478  ipsall += ips3;
479  } //No nead for pp evt
480 
481  if (therr < 0.02) {
482  ips4 = mypow_2[4];
483  ipsall += ips4;
484  }
485  if (pherr < 0.0002) {
486  ips5 = mypow_2[5];
487  ipsall += ips5;
488  }
489  if (abshoang > 0.60 && abshoang < 1.0) {
490  ips6 = mypow_2[6];
491  ipsall += ips6;
492  }
493 
494  if (m_zeroField || (fabs(momatho) > 5.0 && fabs(momatho) < 2000.0)) {
495  ips7 = mypow_2[7];
496  ipsall += ips7;
497  }
498 
499  if (nmuon >= 1 && nmuon <= 3) {
500  ips8 = mypow_2[8];
501  ipsall += ips8;
502  }
503 
504  // if (hodx>0 && hody>0) { }
505  ips9 = mypow_2[9];
506  ipsall += ips9;
507 
508  ips10 = mypow_2[10];
509  ipsall += ips10;
510 
511  if (iring2 == 2) {
512  if (fabs(hodx) < 100 && fabs(hodx) > 2 && fabs(hocorsig[8]) < 40 && fabs(hocorsig[8]) > 2) {
513  ips11 = mypow_2[11];
514  ipsall += ips11;
515  }
516 
517  if (fabs(hody) < 100 && fabs(hody) > 2 && fabs(hocorsig[9]) < 40 && fabs(hocorsig[9]) > 2) {
518  ips12 = mypow_2[12];
519  ipsall += ips12;
520  }
521 
522  } else {
523  if (fabs(hodx) < 100 && fabs(hodx) > 2) {
524  ips11 = mypow_2[11];
525  ipsall += ips11;
526  }
527 
528  if (fabs(hody) < 100 && fabs(hody) > 2) {
529  ips12 = mypow_2[12];
530  ipsall += ips12;
531  }
532  }
533 
534  if (m_zeroField) {
535  if (iring2 == 0) {
536  if (htime > -60 && htime < 60) {
537  ips13 = mypow_2[13];
538  ipsall += ips13;
539  }
540  }
541  if (iring2 == 1) {
542  if (htime > -60 && htime < 60) {
543  ips13 = mypow_2[13];
544  ipsall += ips13;
545  }
546  }
547  if (iring2 == 2) {
548  if (htime > -60 && htime < 60) {
549  ips13 = mypow_2[13];
550  ipsall += ips13;
551  }
552  }
553  if (iring2 == 3) {
554  if (htime > -60 && htime < 60) {
555  ips13 = mypow_2[13];
556  ipsall += ips13;
557  }
558  }
559  if (iring2 == 4) {
560  if (htime > -60 && htime < 60) {
561  ips13 = mypow_2[13];
562  ipsall += ips13;
563  }
564  }
565  } else {
566  if (htime > -100 && htime < 100) {
567  ips13 = mypow_2[13];
568  ipsall += ips13;
569  }
570  }
571  } else {
572  if (abs(ndof) >= 10 && abs(ndof) < 25) {
573  ips0 = mypow_2[0];
574  ipsall += ips0;
575  }
576  if (chisq > 0 && chisq < 10) {
577  ips1 = mypow_2[1];
578  ipsall += ips1;
579  } //18Jan2008
580 
581  if (fabs(trkth - pival / 2) < 21.5) {
582  ips2 = mypow_2[2];
583  ipsall += ips2;
584  } //No nead for pp evt
585  if (fabs(trkph + pival / 2) < 21.5) {
586  ips3 = mypow_2[3];
587  ipsall += ips3;
588  } //No nead for pp evt
589 
590  if (therr < 0.00002) {
591  ips4 = mypow_2[4];
592  ipsall += ips4;
593  }
594  if (pherr < 0.000002) {
595  ips5 = mypow_2[5];
596  ipsall += ips5;
597  }
598  // if (abshoang >0.40 && abshoang <1.0) {ips6 = mypow_2[6]; ipsall +=ips6;}
599  if (tmphoang < 0.065) {
600  ips6 = mypow_2[6];
601  ipsall += ips6;
602  }
603 
604  if (fabs(momatho) < 250.0 && fabs(momatho) > 15.0) {
605  if (iring2 == 2) {
606  ips7 = mypow_2[7];
607  ipsall += ips7;
608  }
609  if ((iring2 == 1 || iring2 == 3) && fabs(momatho) > 17.0) {
610  ips7 = mypow_2[7];
611  ipsall += ips7;
612  }
613  if ((iring2 == 0 || iring2 == 4) && fabs(momatho) > 20.0) {
614  ips7 = mypow_2[7];
615  ipsall += ips7;
616  }
617  }
618 
619  if (nmuon >= 1 && nmuon <= 3) {
620  ips8 = mypow_2[8];
621  ipsall += ips8;
622  }
623 
624  if (ndof > 0 && caloen[0] < 15.0) {
625  ips9 = mypow_2[9];
626  ipsall += ips9;
627  } //5.0
628  if (tkpt03 < 5.0) {
629  ips10 = mypow_2[10];
630  ipsall += ips10;
631  } //4.0
632 
633  if (iring2 == 2) {
634  if (fabs(hodx) < 100 && fabs(hodx) > 2 && fabs(hocorsig[8]) < 40 && fabs(hocorsig[8]) > 2) {
635  ips11 = mypow_2[11];
636  ipsall += ips11;
637  }
638 
639  if (fabs(hody) < 100 && fabs(hody) > 2 && fabs(hocorsig[9]) < 40 && fabs(hocorsig[9]) > 2) {
640  ips12 = mypow_2[12];
641  ipsall += ips12;
642  }
643 
644  } else {
645  if (fabs(hodx) < 100 && fabs(hodx) > 2) {
646  ips11 = mypow_2[11];
647  ipsall += ips11;
648  }
649 
650  if (fabs(hody) < 100 && fabs(hody) > 2) {
651  ips12 = mypow_2[12];
652  ipsall += ips12;
653  }
654  }
655 
656  if (iring2 == 0) {
657  if (htime > -20 && htime < 20) {
658  ips13 = mypow_2[13];
659  ipsall += ips13;
660  }
661  }
662  if (iring2 == 1) {
663  if (htime > -20 && htime < 20) {
664  ips13 = mypow_2[13];
665  ipsall += ips13;
666  }
667  }
668  if (iring2 == 2) {
669  if (htime > -30 && htime < 20) {
670  ips13 = mypow_2[13];
671  ipsall += ips13;
672  }
673  }
674  if (iring2 == 3) {
675  if (htime > -20 && htime < 20) {
676  ips13 = mypow_2[13];
677  ipsall += ips13;
678  }
679  }
680  if (iring2 == 4) {
681  if (htime > -20 && htime < 20) {
682  ips13 = mypow_2[13];
683  ipsall += ips13;
684  }
685  }
686  }
687  int tmpxet = invert_HOieta(ieta);
688  double nomHOSig = hosig[4] / elos;
689 
690  if (ipsall - ips0 == mypow_2[ncut] - mypow_2[0] - 1) {
691  if (isZSps) {
692  sigvsevt[iring2][0]->Fill(abs(ndof), nomHOSig);
693  sig_eta_evt[tmpxet][0]->Fill(abs(ndof), nomHOSig);
694  }
695  variab[iring2][0]->Fill(abs(ndof));
696  }
697  if (ipsall - ips1 == mypow_2[ncut] - mypow_2[1] - 1) {
698  if (isZSps) {
699  sigvsevt[iring2][1]->Fill(chisq, nomHOSig);
700  sig_eta_evt[tmpxet][1]->Fill(chisq, nomHOSig);
701  }
702  variab[iring2][1]->Fill(chisq);
703  }
704  if (ipsall - ips2 == mypow_2[ncut] - mypow_2[2] - 1) {
705  if (isZSps) {
706  sigvsevt[iring2][2]->Fill(trkth, nomHOSig);
707  sig_eta_evt[tmpxet][2]->Fill(trkth, nomHOSig);
708  }
709  variab[iring2][2]->Fill(trkth);
710  }
711  if (ipsall - ips3 == mypow_2[ncut] - mypow_2[3] - 1) {
712  if (isZSps) {
713  sigvsevt[iring2][3]->Fill(trkph, nomHOSig);
714  sig_eta_evt[tmpxet][3]->Fill(trkph, nomHOSig);
715  }
716  variab[iring2][3]->Fill(trkph);
717  }
718  if (ipsall - ips4 == mypow_2[ncut] - mypow_2[4] - 1) {
719  if (isZSps) {
720  sigvsevt[iring2][4]->Fill(1000 * therr, nomHOSig);
721  sig_eta_evt[tmpxet][4]->Fill(1000 * therr, nomHOSig);
722  }
723  variab[iring2][4]->Fill(1000 * therr);
724  }
725  if (ipsall - ips5 == mypow_2[ncut] - mypow_2[5] - 1) {
726  if (isZSps) {
727  sigvsevt[iring2][5]->Fill(1000 * pherr, nomHOSig);
728  sig_eta_evt[tmpxet][5]->Fill(1000 * pherr, nomHOSig);
729  }
730  variab[iring2][5]->Fill(1000 * pherr);
731  }
732  if (ipsall - ips6 == mypow_2[ncut] - mypow_2[6] - 1) {
733  if (isZSps) {
734  sigvsevt[iring2][6]->Fill(tmphoang, (nomHOSig)*abshoang);
735  sig_eta_evt[tmpxet][6]->Fill(tmphoang, (nomHOSig)*abshoang);
736  }
737  variab[iring2][6]->Fill(tmphoang);
738  }
739  if (ipsall - ips7 == mypow_2[ncut] - mypow_2[7] - 1) {
740  if (isZSps) {
741  sigvsevt[iring2][7]->Fill(fabs(trkmm), nomHOSig);
742  sig_eta_evt[tmpxet][7]->Fill(fabs(trkmm), nomHOSig);
743  }
744  variab[iring2][7]->Fill(fabs(trkmm));
745  }
746  if (ipsall - ips8 == mypow_2[ncut] - mypow_2[8] - 1) {
747  if (isZSps) {
748  sigvsevt[iring2][8]->Fill(nmuon, nomHOSig);
749  sig_eta_evt[tmpxet][8]->Fill(nmuon, nomHOSig);
750  }
751  variab[iring2][8]->Fill(nmuon);
752  }
753  if (!m_cosmic) {
754  if (ipsall - ips9 == mypow_2[ncut] - mypow_2[9] - 1) {
755  if (isZSps) {
756  sigvsevt[iring2][9]->Fill(caloen[0], nomHOSig);
757  sig_eta_evt[tmpxet][9]->Fill(caloen[0], nomHOSig);
758  }
759  variab[iring2][9]->Fill(caloen[0]);
760  }
761  }
762  if (ipsall - ips10 == mypow_2[ncut] - mypow_2[10] - 1) {
763  if (isZSps) {
764  sigvsevt[iring2][10]->Fill(tkpt03, nomHOSig);
765  sig_eta_evt[tmpxet][10]->Fill(tkpt03, nomHOSig);
766  }
767  variab[iring2][10]->Fill(tkpt03);
768  }
769  if (ipsall - ips11 == mypow_2[ncut] - mypow_2[11] - 1) {
770  if (isZSps) {
771  sigvsevt[iring2][11]->Fill(hodx, nomHOSig);
772  sig_eta_evt[tmpxet][11]->Fill(hodx, nomHOSig);
773  }
774  variab[iring2][11]->Fill(hodx);
775  }
776  if (ipsall - ips12 == mypow_2[ncut] - mypow_2[12] - 1) {
777  if (isZSps) {
778  sigvsevt[iring2][12]->Fill(hody, nomHOSig);
779  sig_eta_evt[tmpxet][12]->Fill(hody, nomHOSig);
780  }
781  variab[iring2][12]->Fill(hody);
782  }
783 
784  if (ipsall - ips13 == mypow_2[ncut] - mypow_2[13] - 1) {
785  if (isZSps) {
786  sigvsevt[iring2][13]->Fill(htime, nomHOSig);
787  sig_eta_evt[tmpxet][13]->Fill(htime, nomHOSig);
788  }
789  variab[iring2][13]->Fill(htime);
790  }
791 
792  if (isZSps) {
793  sigvsevt[iring2 + ringmx][0]->Fill(abs(ndof), nomHOSig);
794  sig_eta_evt[netamx + tmpxet][0]->Fill(abs(ndof), nomHOSig);
795  }
796  variab[iring2 + 5][0]->Fill(abs(ndof));
797 
798  ncount[iring2][0]++;
799  if (isZSps) {
800  ncount[iring2][1]++;
801  }
802  if (ips0 > 0) {
803  if (isZSps) {
804  ncount[iring2][10]++;
805  sigvsevt[iring2 + ringmx][1]->Fill(chisq, nomHOSig);
806  sig_eta_evt[netamx + tmpxet][1]->Fill(chisq, nomHOSig);
807  }
808  variab[iring2 + ringmx][1]->Fill(chisq);
809  mu_projection[1]->Fill(ieta, iphi);
810  if (ips1 > 0) {
811  if (isZSps) {
812  ncount[iring2][11]++;
813  sigvsevt[iring2 + ringmx][2]->Fill(trkth, nomHOSig);
814  sig_eta_evt[netamx + tmpxet][2]->Fill(trkth, nomHOSig);
815  }
816  variab[iring2 + ringmx][2]->Fill(trkth);
817  mu_projection[2]->Fill(ieta, iphi);
818  if (ips2 > 0) {
819  if (isZSps) {
820  ncount[iring2][12]++;
821  sigvsevt[iring2 + ringmx][3]->Fill(trkph, nomHOSig);
822  sig_eta_evt[netamx + tmpxet][3]->Fill(trkph, nomHOSig);
823  }
824  variab[iring2 + ringmx][3]->Fill(trkph);
825  mu_projection[3]->Fill(ieta, iphi);
826  if (ips3 > 0) {
827  if (isZSps) {
828  ncount[iring2][13]++;
829  sigvsevt[iring2 + ringmx][4]->Fill(1000 * therr, nomHOSig);
830  sig_eta_evt[netamx + tmpxet][4]->Fill(1000 * therr, nomHOSig);
831  }
832  variab[iring2 + ringmx][4]->Fill(1000 * therr);
833  mu_projection[4]->Fill(ieta, iphi);
834  if (ips4 > 0) {
835  if (isZSps) {
836  ncount[iring2][14]++;
837  sigvsevt[iring2 + ringmx][5]->Fill(1000 * pherr, nomHOSig);
838  sig_eta_evt[netamx + tmpxet][5]->Fill(1000 * pherr, nomHOSig);
839  }
840  variab[iring2 + ringmx][5]->Fill(1000 * pherr);
841  mu_projection[5]->Fill(ieta, iphi);
842  if (ips5 > 0) {
843  if (isZSps) {
844  ncount[iring2][15]++;
845  sigvsevt[iring2 + ringmx][6]->Fill(tmphoang, (nomHOSig)*abshoang);
846  sig_eta_evt[netamx + tmpxet][6]->Fill(tmphoang, (nomHOSig)*abshoang);
847  }
848  variab[iring2 + ringmx][6]->Fill(tmphoang);
849  mu_projection[6]->Fill(ieta, iphi);
850  if (ips6 > 0) {
851  if (isZSps) {
852  ncount[iring2][16]++;
853  sigvsevt[iring2 + ringmx][7]->Fill(fabs(trkmm), nomHOSig);
854  sig_eta_evt[netamx + tmpxet][7]->Fill(fabs(trkmm), nomHOSig);
855  }
856  variab[iring2 + ringmx][7]->Fill(fabs(trkmm));
857  mu_projection[7]->Fill(ieta, iphi);
858  if (ips7 > 0) {
859  ncount[iring2][4]++; //Efficiency of Muon detection
860  if (isZSps) {
861  ncount[iring2][17]++;
862  sigvsevt[iring2 + ringmx][8]->Fill(nmuon, nomHOSig);
863  sig_eta_evt[netamx + tmpxet][8]->Fill(nmuon, nomHOSig);
864  }
865  variab[iring2 + ringmx][8]->Fill(nmuon);
866  mu_projection[8]->Fill(ieta, iphi);
867  if (ips8 > 0) {
868  if (!m_cosmic) {
869  if (isZSps) {
870  ncount[iring2][18]++;
871  sigvsevt[iring2 + ringmx][9]->Fill(caloen[0], nomHOSig);
872  sig_eta_evt[netamx + tmpxet][9]->Fill(caloen[0], nomHOSig);
873  }
874  variab[iring2 + ringmx][9]->Fill(caloen[0]);
875  mu_projection[9]->Fill(ieta, iphi);
876  }
877  if (ips9 > 0) {
878  if (isZSps) {
879  ncount[iring2][19]++;
880  sigvsevt[iring2 + ringmx][10]->Fill(tkpt03, nomHOSig);
881  sig_eta_evt[netamx + tmpxet][10]->Fill(tkpt03, nomHOSig);
882  }
883  variab[iring2 + ringmx][10]->Fill(tkpt03);
884  mu_projection[10]->Fill(ieta, iphi);
885  if (ips10 > 0) {
886  ncount[iring2][3]++; //Efficiency of Muon detection
887  if (isZSps) {
888  ncount[iring2][20]++;
889  sigvsevt[iring2 + ringmx][11]->Fill(hodx, nomHOSig);
890  sig_eta_evt[netamx + tmpxet][11]->Fill(hodx, nomHOSig);
891  }
892  variab[iring2 + ringmx][11]->Fill(hodx);
893  mu_projection[11]->Fill(ieta, iphi);
894 
895  if (ips11 > 0) {
896  if (isZSps) {
897  ncount[iring2][21]++;
898  sigvsevt[iring2 + ringmx][12]->Fill(hody, nomHOSig);
899  sig_eta_evt[netamx + tmpxet][12]->Fill(hody, nomHOSig);
900  }
901  variab[iring2 + ringmx][12]->Fill(hody);
902  mu_projection[12]->Fill(ieta, iphi);
903 
904  if (ips12 > 0) {
905  ncount[iring2][2]++; //Efficiency of Muon detection
906  if (isZSps) {
907  ncount[iring2][22]++;
908  sigvsevt[iring2 + ringmx][13]->Fill(htime, nomHOSig);
909  sig_eta_evt[tmpxet + ringmx][13]->Fill(htime, nomHOSig);
910  }
911  variab[iring2 + ringmx][13]->Fill(htime);
912  mu_projection[13]->Fill(ieta, iphi);
913 
914  if (ips13 > 0) {
915  if (isZSps) {
916  ncount[iring2][23]++;
917  mu_projection[14]->Fill(ieta, iphi);
918  }
919  }
920  }
921  }
922  }
923  }
924  }
925  }
926  }
927  }
928  }
929  }
930  }
931  }
932  }
933  if (isZSps) {
934  sigvsevt[iring2 + 2 * ringmx][0]->Fill(abs(ndof), nomHOSig);
935  sigvsevt[iring2 + 2 * ringmx][1]->Fill(chisq, nomHOSig);
936  sigvsevt[iring2 + 2 * ringmx][2]->Fill(trkth, nomHOSig);
937  sigvsevt[iring2 + 2 * ringmx][3]->Fill(trkph, nomHOSig);
938  sigvsevt[iring2 + 2 * ringmx][4]->Fill(1000 * therr, nomHOSig);
939  sigvsevt[iring2 + 2 * ringmx][5]->Fill(1000 * pherr, nomHOSig);
940  if (abshoang > 0.01) {
941  sigvsevt[iring2 + 2 * ringmx][6]->Fill(tmphoang, (nomHOSig)*abshoang);
942  }
943  sigvsevt[iring2 + 2 * ringmx][7]->Fill(fabs(trkmm), nomHOSig);
944  sigvsevt[iring2 + 2 * ringmx][8]->Fill(nmuon, nomHOSig);
945  if (!m_cosmic)
946  sigvsevt[iring2 + 2 * ringmx][9]->Fill(caloen[0], nomHOSig);
947  sigvsevt[iring2 + 2 * ringmx][10]->Fill(tkpt03, nomHOSig);
948  sigvsevt[iring2 + 2 * ringmx][11]->Fill(hodx, nomHOSig);
949  sigvsevt[iring2 + 2 * ringmx][12]->Fill(hody, nomHOSig);
950  sigvsevt[iring2 + 2 * ringmx][13]->Fill(htime, nomHOSig);
951 
952  sig_eta_evt[2 * netamx + tmpxet][0]->Fill(abs(ndof), nomHOSig);
953  sig_eta_evt[2 * netamx + tmpxet][1]->Fill(chisq, nomHOSig);
954  sig_eta_evt[2 * netamx + tmpxet][2]->Fill(trkth, nomHOSig);
955  sig_eta_evt[2 * netamx + tmpxet][3]->Fill(trkph, nomHOSig);
956  sig_eta_evt[2 * netamx + tmpxet][4]->Fill(1000 * therr, nomHOSig);
957  sig_eta_evt[2 * netamx + tmpxet][5]->Fill(1000 * pherr, nomHOSig);
958  if (abshoang > 0.01) {
959  sig_eta_evt[2 * netamx + tmpxet][6]->Fill(tmphoang, (nomHOSig)*abshoang);
960  }
961  sig_eta_evt[2 * netamx + tmpxet][7]->Fill(fabs(trkmm), nomHOSig);
962  sig_eta_evt[2 * netamx + tmpxet][8]->Fill(nmuon, nomHOSig);
963  if (!m_cosmic)
964  sig_eta_evt[2 * netamx + tmpxet][9]->Fill(caloen[0], nomHOSig);
965  sig_eta_evt[2 * netamx + tmpxet][10]->Fill(tkpt03, nomHOSig);
966  sig_eta_evt[2 * netamx + tmpxet][11]->Fill(hodx, nomHOSig);
967  sig_eta_evt[2 * netamx + tmpxet][12]->Fill(hody, nomHOSig);
968  sig_eta_evt[2 * netamx + tmpxet][13]->Fill(htime, nomHOSig);
969  }
970 
971  variab[iring2 + 2 * ringmx][0]->Fill(abs(ndof));
972  variab[iring2 + 2 * ringmx][1]->Fill(chisq);
973  variab[iring2 + 2 * ringmx][2]->Fill(trkth);
974  variab[iring2 + 2 * ringmx][3]->Fill(trkph);
975  variab[iring2 + 2 * ringmx][4]->Fill(1000 * therr);
976  variab[iring2 + 2 * ringmx][5]->Fill(1000 * pherr);
977  variab[iring2 + 2 * ringmx][6]->Fill(tmphoang);
978  variab[iring2 + 2 * ringmx][7]->Fill(fabs(trkmm));
979  variab[iring2 + 2 * ringmx][8]->Fill(nmuon);
980  if (!m_cosmic)
981  variab[iring2 + 2 * ringmx][9]->Fill(caloen[0]);
982  variab[iring2 + 2 * ringmx][10]->Fill(tkpt03);
983  variab[iring2 + 2 * ringmx][11]->Fill(hodx);
984  variab[iring2 + 2 * ringmx][12]->Fill(hody);
985  variab[iring2 + 2 * ringmx][13]->Fill(htime);
986 
987  muonnm->Fill(nmuon);
988  muonmm->Fill(trkmm);
989  muonth->Fill(trkth * 180 / pival);
990  muonph->Fill(trkph * 180 / pival);
991  muonch->Fill(chisq);
992 
993  int iselect = (ipsall == mypow_2[ncut] - 1) ? 1 : 0;
994 
995  if (iselect == 1) {
996  ipass++;
997  sel_muonnm->Fill(nmuon);
998  sel_muonmm->Fill(trkmm);
999  sel_muonth->Fill(trkth * 180 / pival);
1000  sel_muonph->Fill(trkph * 180 / pival);
1001  sel_muonch->Fill(chisq);
1002  if (m_histFill && tmpxet >= 0 && tmpxet < netamx && iphi >= 0 && iphi < nphimx) {
1003  ho_indenergy[tmpxet][iphi - 1]->Fill(nomHOSig);
1004  }
1005  if (m_treeFill) {
1006  T1->Fill();
1007  }
1008  }
1009  } //for (HOCalibVariableCollection::const_iterator hoC=(*HOCalib).begin(); hoC!=(*HOCalib).end(); hoC++){
1010  } //if (isCosMu)
1011 }
1012 
1013 //define this as a plug-in
1014 DEFINE_FWK_MODULE(HOCalibAnalyzer);
edm::EventID::run
RunNumber_t run() const
Definition: EventID.h:38
TFileService::kSharedResource
static const std::string kSharedResource
Definition: TFileService.h:76
edm::LogVerbatim
Log< level::Info, true > LogVerbatim
Definition: MessageLogger.h:128
edm::EventID::event
EventNumber_t event() const
Definition: EventID.h:40
edm::ParameterSet::getUntrackedParameter
T getUntrackedParameter(std::string const &, T const &) const
log
static std::vector< std::string > checklist log
Definition: LHERunInfoProduct.cc:179
edm::Event::getByToken
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:539
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
HOCalibAnalyzer::HOCalibAnalyzer
HOCalibAnalyzer(const edm::ParameterSet &)
Definition: HOCalibAnalyzer.cc:150
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
HOCalibAnalyzer::beginJob
void beginJob() override
Definition: HOCalibAnalyzer.cc:69
edm::EventBase::luminosityBlock
edm::LuminosityBlockNumber_t luminosityBlock() const
Definition: EventBase.h:61
TFileService::make
T * make(const Args &...args) const
make new ROOT object
Definition: TFileService.h:64
HOCalibAnalyzer::endJob
void endJob() override
Definition: HOCalibAnalyzer.cc:71
HOCalibAnalyzer::analyze
void analyze(const edm::Event &, const edm::EventSetup &) override
Definition: HOCalibAnalyzer.cc:333
HOCalibAnalyzer::hoCalibVariableCollectionTag
edm::InputTag hoCalibVariableCollectionTag
Definition: HOCalibAnalyzer.cc:132
analyze
example_stream void analyze(const edm::Event &, const edm::EventSetup &) override
iEvent
int iEvent
Definition: GenABIO.cc:224
if
if(conf_.getParameter< bool >("UseStripCablingDB"))
Definition: SiStripLAProfileBooker.cc:72
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
HOCalibAnalyzer::getHOieta
int getHOieta(int ij)
Definition: HOCalibAnalyzer.cc:82
HOCalibAnalyzer::invert_HOieta
int invert_HOieta(int ieta)
Definition: HOCalibAnalyzer.cc:83
Max
T Max(T a, T b)
Definition: MathUtil.h:44
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
HOCalibAnalyzer::tok_allho_
edm::EDGetTokenT< HORecHitCollection > tok_allho_
Definition: HOCalibAnalyzer.cc:134
edm::EventBase::id
edm::EventID id() const
Definition: EventBase.h:59
HOCalibAnalyzer::tok_ho_
edm::EDGetTokenT< HOCalibVariableCollection > tok_ho_
Definition: HOCalibAnalyzer.cc:133
HOCalibAnalyzer::~HOCalibAnalyzer
~HOCalibAnalyzer() override
Definition: HOCalibAnalyzer.cc:319
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
Generated for CMSSW Reference Manual by doxygen 1.8.5