CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
AlCaHOCalibProducer.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Feb09 2009
4 // Move the initialisation of SteppingHelixPropagator from ::beginJob() to ::produce()
5 //
6 // Oct3 2008
7 // Difference in tag V00-02-45 with previous code
8 
9 // 1. One new object on data format, which was realised in
10 // CRUZET data analysis.
11 //2. Remove all histogram and cout in the code
12 //3. An upgrade in code, which increases the acceptance of
13 // muon near the edge (this also realised in CRUZET data).
14 // Difference in wrt V00-02-45
15 // 1. initialisation tmpHOCalib.htime = -1000;
16 // 2. By mistake HLT was commented out
17 
18 // Package: AlCaHOCalibProducer
19 // Class: AlCaHOCalibProducer
20 //
53 //
54 // Original Author: Gobinda Majumder
55 // Created: Fri Jul 6 17:17:21 CEST 2007
56 //
57 //
58 
59 
60 // system include files
61 #include <memory>
62 
63 // user include files
66 
73 
78 
85 
91 
92 
95 
100 
101 //08/07/07 #include "CondTools/Hcal/interface/HcalDbPool.h"
102 //#include "CondFormats/HcalObjects/interface/HcalPedestals.h"
103 //#include "CondFormats/HcalObjects/interface/HcalPedestalWidths.h"
104 
105 
106 // #include "TrackingTools/GeomPropagators/interface/HelixArbitraryPlaneCrossing.h"
109 
112 
115 #include "CLHEP/Vector/LorentzVector.h"
116 
119 
124 //#include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixStateInfo.h"
125 
128 
132 
133 
134 #include "TFile.h"
135 #include "TH1F.h"
136 #include "TH2F.h"
137 #include "TProfile.h"
138 #include "TTree.h"
139 /* C++ Headers */
140 #include <string>
141 
142 #include <iostream>
143 #include <fstream>
144 //
145 // class decleration
146 //
147 using namespace std;
148 using namespace edm;
149 using namespace reco;
150 
151 const int netabin= 16;
152 const int nphimx = 72;
153 const int netamx = 32;
154 const int nchnmx = 10;
155 const int ncidmx = 5;
156 
157 //GMA #ifdef DEBUG_OFFLINE_GM
158 const int nsigpk = 7;
159 const int nstrbn = 0;
160 const int ntrgp_gm = 11;
161 
162 
163 const int netahbmx = 60;
164 const int netahb3mx = 32;
165 
166 static const unsigned int nL1trg = 200;
167 
168 static const unsigned int nL1mx=140;
169 static const unsigned int nHLTmx=140;
170 //GMA #endif
171 
173  public:
174  explicit AlCaHOCalibProducer(const edm::ParameterSet&);
176 
180 
181 
182  private:
183  void findHOEtaPhi(int iphsect, int& ietaho, int& iphiho);
184  virtual void beginJob() override ;
185  virtual void produce(edm::Event&, const edm::EventSetup&) override;
186  virtual void endJob() override ;
187 
188  // ----------member data ---------------------------
189 
190  float xhor0; //x-position in ring 0
191  float yhor0; //y-position in ring 0
192  float xhor1; //x-position in ring 1
193  float yhor1; //y-position in ring 1
194  int iring; //Ring number -2,-1,0,1,2
195 
196  float localxhor0; //local x-distance from edege in ring 0
197  float localyhor0; //local y-distance from edege in ring 0
198  float localxhor1; //local x-distance from edege in ring 1
199  float localyhor1; //local y-distance from edege in ring 1
200 
201  float pedestal[netamx][nphimx][ncidmx];
202 
204 
205  bool debug;
207 
208  //GMA #ifdef DEBUG_OFFLINE_GM
209 
210  TH1F* libhoped;
211  TH1F* libhoped1;
212 
213  TH1F* allhotime;
214  TH1F* hotime[ntrgp_gm+1];
215  TH1F* hopedtime;
216 
217  TProfile* hopedpr;
218  TH1F* hopedrms;
219  TH1F* hst_hopedrms;
220 
221  TProfile* hopeak[ntrgp_gm+1];
222  TProfile* horatio;
223 
224  TH1F* Nallhotime;
225  TH1F* Nhotime[ntrgp_gm+1];
226  TH1F* Nhopedtime;
227 
228  TH1F* allhb1;
229  TH1F* allhb2;
230  TH1F* allhb3;
231 
232  TH1F* Nallhb1;
233  TH1F* Nallhb2;
234  TH1F* Nallhb3;
235 
236  TProfile* hb1pedpr;
237  TH1F* hb1pedrms;
239 
240  TH1F* ho_occupency[5];
241 
242  bool m_hotime;
243  //GM #endif
244 
245  edm::InputTag muonTags_; // cosmicMuons or standAloneMuons
246 
253 
254  bool m_digiInput; // digi (true) or rechit (false)
255  bool m_hbinfo;
257  int m_endTS;
258  double m_magscale;
259  double m_sigma;
260 
262  //#ifdef DEBUG_OFFLINE_GM
263  // int Nevents;
264  int Noccu;
265  // int Npass;
266  int nRuns;
267  //#endif
268 
269  int irunold;
270  // SteppingHelixPropagator* stepProp;
271  FreeTrajectoryState getFreeTrajectoryState( const reco::Track& tk, const MagneticField* field, int itag, bool dir);
272 
276 
279 
280  unsigned int Ntp; // # of HLT trigger paths (should be the same for all events!)
281  std::map<std::string, bool> fired;
282 
283 };
284 
285 //
286 // constants, enums and typedefs
287 //
288 
289 //
290 // static data member definitions
291 //
292 
293 //
294 // constructors and destructor
295 //
297 
298 {
299  //register your products
300 
301  theRootFileName = iConfig.getUntrackedParameter<string>("RootFileName","tmp.root");
302  m_digiInput = iConfig.getUntrackedParameter<bool>("digiInput", true);
303  m_hbinfo = iConfig.getUntrackedParameter<bool>("hbinfo", false);
304  m_startTS = iConfig.getUntrackedParameter<int>("firstTS", 4);
305 
306  m_hotime = iConfig.getUntrackedParameter<bool>("hotime", false);
307 
308  if(m_startTS<0) m_startTS=0;
309  m_endTS = iConfig.getUntrackedParameter<int>("lastTS", 7);
310  if (m_endTS < m_startTS) m_endTS = m_startTS + 3;
311  if (m_endTS >9) m_endTS=9;
312  m_magscale = iConfig.getUntrackedParameter<double>("m_scale", 4.0);
313  m_sigma = iConfig.getUntrackedParameter<double>("sigma", 1.0);
314 
315  // keep InputTag muonTags_ since it is used below. - cowden
316  muonTags_ = iConfig.getUntrackedParameter<edm::InputTag>("muons");
317  tok_muons_ = consumes<reco::TrackCollection>(muonTags_);
318  tok_ho_ = consumes<HORecHitCollection>(iConfig.getParameter<edm::InputTag>("hoInput"));
319  tok_hbhe_ = consumes<HBHERecHitCollection>(iConfig.getParameter<edm::InputTag>("hbheInput"));
320  tok_tower_ = consumes<CaloTowerCollection>(iConfig.getParameter<edm::InputTag>("towerInput"));
321 
322  // Since these accesses are currently commented out, I put the registration as "mayConsume". - cowden
323  tok_hlt_ = mayConsume<edm::TriggerResults>(iConfig.getParameter<edm::InputTag>("hltInput"));
324  tok_l1_ = mayConsume<L1GlobalTriggerReadoutRecord>(iConfig.getParameter<edm::InputTag>("l1Input"));
325 
326  produces<HOCalibVariableCollection>("HOCalibVariableCollection").setBranchAlias("HOCalibVariableCollection");
327 
328 
329  if (m_hotime) {
331 
332  char title[200];
333  if ( m_digiInput) {
334  libhoped = fs->make<TH1F>("libhoped", "libhoped", ncidmx*netamx*nphimx, -0.5, ncidmx*netamx*nphimx-0.5);
335  libhoped1 = fs->make<TH1F>("libhoped1", "libhoped1", nchnmx*netamx*nphimx, -0.5, nchnmx*netamx*nphimx-0.5);
336  allhotime = fs->make<TH1F>("allhotime", "allhotime", nchnmx*netamx*nphimx, -0.5, nchnmx*netamx*nphimx-0.5);
337  for (int ij=0; ij<=ntrgp_gm; ij++) {
338  sprintf(title, "hotime_trgp_%i", ij+1);
339  hotime[ij] = fs->make<TH1F>(title, title, nchnmx*netamx*nphimx, -0.5, nchnmx*netamx*nphimx-0.5);
340  sprintf(title, "hopeak_trgp_%i", ij+1);
341  hopeak[ij] = fs->make<TProfile>(title, title,netamx*nphimx, -0.5, netamx*nphimx-0.5);
342  }
343 
344  horatio = fs->make<TProfile>("horatio", "horatio",netamx*nphimx, -0.5, netamx*nphimx-0.5);
345  hopedtime = fs->make<TH1F>("hopedtime", "hopedtime", nchnmx*netamx*nphimx, -0.5, nchnmx*netamx*nphimx-0.5);
346 
347  Nallhotime = fs->make<TH1F>("Nallhotime", "Nallhotime", nchnmx*netamx*nphimx, -0.5, nchnmx*netamx*nphimx-0.5);
348  hopedpr = fs->make<TProfile>("hopedpr", "hopedpr", nchnmx*netamx*nphimx, -0.5, nchnmx*netamx*nphimx-0.5);
349  hopedrms = fs->make<TH1F>("hopedrms", "hopedrms", nchnmx*netamx*nphimx, -0.5, nchnmx*netamx*nphimx-0.5);
350  hst_hopedrms = fs->make<TH1F>("hst_hopedrms", "hst_hopedrms", 100, 0.0, 0.1);
351  for (int ij=0; ij<=ntrgp_gm; ij++) {
352  sprintf(title, "Nhotime_trgp_%i", ij+1);
353  Nhotime[ij] = fs->make<TH1F>(title, title, nchnmx*netamx*nphimx, -0.5, nchnmx*netamx*nphimx-0.5);
354  }
355  Nhopedtime = fs->make<TH1F>("Nhopedtime", "Nhopedtime", nchnmx*netamx*nphimx, -0.5, nchnmx*netamx*nphimx-0.5);
356  allhb1 = fs->make<TH1F>("allhb1", "allhb1", nchnmx*netahbmx*nphimx, -0.5, nchnmx*netahbmx*nphimx-0.5);
357  allhb2 = fs->make<TH1F>("allhb2", "allhb2", nchnmx*netahb3mx*nphimx, -0.5, nchnmx*netahb3mx*nphimx-0.5);
358  allhb3 = fs->make<TH1F>("allhb3", "allhb3", nchnmx*netahb3mx*nphimx, -0.5, nchnmx*netahb3mx*nphimx-0.5);
359  Nallhb1 = fs->make<TH1F>("Nallhb1", "Nallhb1", nchnmx*netahbmx*nphimx, -0.5, nchnmx*netahbmx*nphimx-0.5);
360  Nallhb2 = fs->make<TH1F>("Nallhb2", "Nallhb2", nchnmx*netahb3mx*nphimx, -0.5, nchnmx*netahb3mx*nphimx-0.5);
361  Nallhb3 = fs->make<TH1F>("Nallhb3", "Nallhb3", nchnmx*netahb3mx*nphimx, -0.5, nchnmx*netahb3mx*nphimx-0.5);
362  hb1pedpr = fs->make<TProfile>("hb1pedpr", "hb1pedpr", nchnmx*netahbmx*nphimx, -0.5, nchnmx*netahbmx*nphimx-0.5);
363  hb1pedrms = fs->make<TH1F>("hb1pedrms", "hb1pedrms", nchnmx*netahbmx*nphimx, -0.5, nchnmx*netahbmx*nphimx-0.5);
364  hst_hb1pedrms = fs->make<TH1F>("hst_hb1pedrms", "hst_hb1pedrms", 100, 0., 0.1);
365 
366  }
367  for (int i=0; i<5; i++) {
368  sprintf(title, "ho_occupency (>%i #sigma)", i+2);
369  ho_occupency[i] = fs->make<TH1F>(title, title, netamx*nphimx, -0.5, netamx*nphimx-0.5);
370  }
371  }
372 
373 }
374 
376 {
377 
378  // do anything here that needs to be done at desctruction time
379  // (e.g. close files, deallocate resources etc.)
380 
381  if (m_hotime) {
382  // Write the histos to file
383  if ( m_digiInput) {
384  allhotime->Divide(Nallhotime);
385  for (int ij=0; ij<=ntrgp_gm; ij++) {
386  hotime[ij]->Divide(Nhotime[ij]);
387  }
388  hopedtime->Divide(Nhopedtime);
389  libhoped->Scale(1./max(1,nRuns));
390  libhoped1->Scale(1./max(1,nRuns));
391  for (int i=0; i<nchnmx*netamx*nphimx; i++) {
392  float xx = hopedpr->GetBinError(i+1);
393  if (hopedpr->GetBinEntries(i+1) >0) {
394  hopedrms->Fill(i, xx);
395  hst_hopedrms->Fill(xx);
396  }
397  }
398  allhb1->Divide(Nallhb1);
399  allhb2->Divide(Nallhb2);
400  allhb3->Divide(Nallhb3);
401  for (int i=0; i<nchnmx*netahbmx*nphimx; i++) {
402  float xx = hb1pedpr->GetBinError(i+1);
403  if (hb1pedpr->GetBinEntries(i+1) >0) {
404  hb1pedrms->Fill(i, xx);
405  hst_hb1pedrms->Fill(xx);
406  }
407  }
408  }
409  for (int i=0; i<5; i++) {
410  ho_occupency[i]->Scale(1./max(1,Noccu));
411  }
412  }
413 
414 }
415 
416 
417 //
418 // member functions
419 //
420 
421 // ------------ method called to produce the data ------------
422 void
424 {
425 
426  using namespace edm;
427  int irun = iEvent.id().run();
428  if (m_digiInput) {
429  if (irunold !=irun) {
430  iSetup.get<HcalDbRecord>().get(conditions_);
431 
432  for (int i=0; i<netamx; i++) {
433  for (int j=0; j<nphimx; j++) {
434  for (int k=0; k<ncidmx; k++) {
435  pedestal[i][j][k]=0.0;
436  }
437  }
438  }
439  }
440  }
441 
442  // if (m_hotime && m_digiInput) {
443  if (m_digiInput) {
444  if (irunold !=irun) {
445  nRuns++;
446  for (int i =-netabin+1; i<=netabin-1; i++) {
447  if (i==0) continue;
448  int tmpeta1 = (i>0) ? i -1 : -i +14;
449  if (tmpeta1 <0 || tmpeta1 >netamx) continue;
450  for (int j=0; j<nphimx; j++) {
451 
452  HcalDetId id(HcalOuter, i, j+1, 4);
453  calibped = conditions_->getHcalCalibrations(id);
454 
455  for (int k =0; k<ncidmx-1; k++) {
456  pedestal[tmpeta1][j][k] = calibped.pedestal(k); // pedm->getValue(k);
457  pedestal[tmpeta1][j][ncidmx-1] += (1./(ncidmx-1))*pedestal[tmpeta1][j][k];
458  }
459 
460  if (m_hotime) {
461  for (int k =0; k<ncidmx; k++) {
462  libhoped->Fill(nphimx*ncidmx*tmpeta1 + ncidmx*j + k, pedestal[tmpeta1][j][k]);
463  }
464  for (int k =0; k<nchnmx; k++) {
465  libhoped1->Fill(nphimx*nchnmx*tmpeta1 + nchnmx*j + k, pedestal[tmpeta1][j][min(k,ncidmx-1)]);
466  }
467  }
468 
469  }
470  }
471  }
472  }
473 
474  // Nevents++;
475  irunold = irun;
476 
477  //GMA if (Nevents%500==1)
478  //GMA cout <<"AlCaHOCalibProducer Processing event # "<<Nevents<<" "<<Npass<<" "<<Noccu<<" "<<irun<<" "<<iEvent.id().event()<<endl;
479 
480  std::auto_ptr<HOCalibVariableCollection> hostore (new HOCalibVariableCollection);
481 
483 
485 
486  if (m_digiInput) {
487  iEvent.getByToken(tok_ho_,ho);
488  iEvent.getByToken(tok_hbhe_,hbhe);
489  }
490 
491  if (m_hotime && m_digiInput) {
492  if ((*ho).size()>0) {
493  for (HODigiCollection::const_iterator j=(*ho).begin(); j!=(*ho).end(); j++){
494  HcalDetId id =(*j).id();
495  m_coder = (*conditions_).getHcalCoder(id);
496  m_shape = (*conditions_).getHcalShape(m_coder);
497  int tmpeta= id.ieta();
498  int tmpphi= id.iphi();
499  float tmpdata[nchnmx];
500  int tmpeta1 = (tmpeta>0) ? tmpeta -1 : -tmpeta +14;
501  for (int i=0; i<(*j).size() && i<nchnmx; i++) {
502  tmpdata[i] = m_coder->charge(*m_shape,(*j).sample(i).adc(),(*j).sample(i).capid());
503  allhotime->Fill(nphimx*nchnmx*tmpeta1 + nchnmx*(tmpphi-1) + i, tmpdata[i]);
504  Nallhotime->Fill(nphimx*nchnmx*tmpeta1 + nchnmx*(tmpphi-1) + i, 1.);
505  }
506  }
507  }
508  if ((*hbhe).size()>0) {
509  for (HBHEDigiCollection::const_iterator j=(*hbhe).begin(); j!=(*hbhe).end(); j++){
510  HcalDetId id =(*j).id();
511  m_coder = (*conditions_).getHcalCoder(id);
512  m_shape = (*conditions_).getHcalShape(m_coder);
513  int tmpeta= id.ieta();
514  int tmpphi= id.iphi();
515  int tmpdepth =id.depth();
516  int tmpeta1 = (tmpeta>0) ? tmpeta -15 : -tmpeta + 1;
517  if (tmpdepth==1) tmpeta1 = (tmpeta>0) ? tmpeta -1 : -tmpeta +29;
518  for (int i=0; i<(*j).size() && i<nchnmx; i++) {
519  float signal = m_coder->charge(*m_shape,(*j).sample(i).adc(),(*j).sample(i).capid());
520  if (tmpdepth==1) {
521  allhb1->Fill(nphimx*nchnmx*tmpeta1 + nchnmx*(tmpphi-1) + i, signal);
522  Nallhb1->Fill(nphimx*nchnmx*tmpeta1 + nchnmx*(tmpphi-1) + i, 1);
523  hb1pedpr->Fill(nphimx*nchnmx*tmpeta1 + nchnmx*(tmpphi-1) + i, signal);}
524  if (tmpdepth==2) {
525  allhb2->Fill(nphimx*nchnmx*tmpeta1 + nchnmx*(tmpphi-1) + i, signal);
526  Nallhb2->Fill(nphimx*nchnmx*tmpeta1 + nchnmx*(tmpphi-1) + i, 1);}
527  if (tmpdepth==3) {
528  allhb3->Fill(nphimx*nchnmx*tmpeta1 + nchnmx*(tmpphi-1) + i, signal);
529  Nallhb3->Fill(nphimx*nchnmx*tmpeta1 + nchnmx*(tmpphi-1) + i, 1);}
530  }
531  }
532  }
533  }
534 
535  double pival = acos(-1.);
536 
538  iEvent.getByToken(tok_muons_, cosmicmuon);
539 
540  if (cosmicmuon->size()>0) {
541 
542  int l1trg = 0;
543  int hlttr = 0;
544 
545  int ntrgpas_gm[ntrgp_gm]={0,0,0,0,0,0,0,0,0,0};
546 
547  /*
548  //L1 trigger
549  Handle<L1GlobalTriggerReadoutRecord> L1GTRR;
550  iEvent.getByToken(tok_l1_,L1GTRR); //gtDigis
551 
552  if ( L1GTRR.isValid()) {
553  const unsigned int n(L1GTRR->decisionWord().size());
554  const bool accept(L1GTRR->decision());
555  if (accept) {
556  for (unsigned int i=0; i!=n && i<32; ++i) {
557  // for (unsigned int i=0; i!=n ; ++i) {
558  int il1trg = (L1GTRR->decisionWord()[i]) ? 1 : 0;
559  if (il1trg>0 && i<32) l1trg +=int(std::pow(2., double(i%32))*il1trg);
560  }
561  }
562  }// else { return;}
563 
564  //HLT
565 
566  Handle<edm::TriggerResults> trigRes;
567  iEvent.getByToken(tok_hlt_, trigRes);
568 
569 
570  unsigned int size = trigRes->size();
571  edm::TriggerNames triggerNames(*trigRes);
572 
573  // loop over all paths, get trigger decision
574  for(unsigned i = 0; i != size && i<32; ++i) {
575  std::string name = triggerNames.triggerName(i);
576  fired[name] = trigRes->accept(i);
577  int ihlt = trigRes->accept(i);
578  if (m_hotime){
579  if (ihlt >0 && i < (int)ntrgp_gm) { ntrgpas_gm[i] = 1;}
580  }
581  if (i<32 && ihlt>0) hlttr += int(std::pow(2., double(i%32))*ihlt);
582  }
583 
584  */
585 
586  int Noccu_old = Noccu;
587 
588  for(reco::TrackCollection::const_iterator ncosm = cosmicmuon->begin();
589  ncosm != cosmicmuon->end(); ++ncosm) {
590 
591  if ((*ncosm).ndof() < 15) continue;
592  if ((*ncosm).normalizedChi2() >30.0) continue;
593 
594  HOCalibVariables tmpHOCalib;
595 
596  tmpHOCalib.trig1 = l1trg;
597  tmpHOCalib.trig2 = hlttr;
598 
599  int charge = ncosm->charge();
600 
601  double innerr = (*ncosm).innerPosition().Perp2();
602  double outerr = (*ncosm).outerPosition().Perp2();
603  int iiner = (innerr <outerr) ? 1 : 0;
604 
605  //---------------------------------------------------
606  // in_to_out Dir in_to_out Dir
607  // StandAlone ^ ^ Cosmic ^ |
608  // | | | v
609  //---------------------------------------------------Y=0
610  // StandAlone | | Cosmic ^ |
611  // v v | v
612  //----------------------------------------------------
613 
614  double posx, posy, posz;
615  double momx, momy, momz;
616 
617  if (iiner==1) {
618  posx = (*ncosm).innerPosition().X();
619  posy = (*ncosm).innerPosition().Y();
620  posz = (*ncosm).innerPosition().Z();
621 
622  momx = (*ncosm).innerMomentum().X();
623  momy = (*ncosm).innerMomentum().Y();
624  momz = (*ncosm).innerMomentum().Z();
625 
626  } else {
627  posx = (*ncosm).outerPosition().X();
628  posy = (*ncosm).outerPosition().Y();
629  posz = (*ncosm).outerPosition().Z();
630 
631  momx = (*ncosm).outerMomentum().X();
632  momy = (*ncosm).outerMomentum().Y();
633  momz = (*ncosm).outerMomentum().Z();
634  }
635 
636 
637  PositionType trkpos(posx, posy, posz);
638 
639  CLHEP::Hep3Vector tmpmuon3v(posx, posy, posz);
640  CLHEP::Hep3Vector tmpmuondir(momx, momy, momz);
641 
642  bool samedir = (tmpmuon3v.dot(tmpmuondir) >0) ? true : false;
643  for (int i=0; i<3; i++) {tmpHOCalib.caloen[i] = 0.0;}
644  int inearbymuon = 0;
645  for(reco::TrackCollection::const_iterator ncosmcor = cosmicmuon->begin();
646  ncosmcor != cosmicmuon->end(); ++ncosmcor) {
647  if (ncosmcor==ncosm) continue;
648 
649  CLHEP::Hep3Vector tmpmuon3vcor;
650  CLHEP::Hep3Vector tmpmom3v;
651  if (iiner==1) {
652  tmpmuon3vcor = CLHEP::Hep3Vector((*ncosmcor).innerPosition().X(),(*ncosmcor).innerPosition().Y(),(*ncosmcor).innerPosition().Z());
653  tmpmom3v = CLHEP::Hep3Vector((*ncosmcor).innerMomentum().X(),(*ncosmcor).innerMomentum().Y(),(*ncosmcor).innerMomentum().Z());
654  } else {
655  tmpmuon3vcor = CLHEP::Hep3Vector((*ncosmcor).outerPosition().X(),(*ncosmcor).outerPosition().Y(),(*ncosmcor).outerPosition().Z());
656  tmpmom3v = CLHEP::Hep3Vector((*ncosmcor).outerMomentum().X(),(*ncosmcor).outerMomentum().Y(),(*ncosmcor).outerMomentum().Z());
657 
658  }
659  if (tmpmom3v.mag()<0.2 || (*ncosmcor).ndof()<5) continue;
660 
661  double angle = tmpmuon3v.angle(tmpmuon3vcor);
662  if (angle < 7.5*pival/180.) {inearbymuon=1;} // break;}
663 
664  if (muonTags_.label() =="cosmicMuons") {
665  if (angle <7.5*pival/180.) { tmpHOCalib.caloen[0] +=1.;}
666  if (angle <15.0*pival/180.) { tmpHOCalib.caloen[1] +=1.;}
667  if (angle <35.0*pival/180.) { tmpHOCalib.caloen[2] +=1.;}
668  }
669  }
670 
671  localxhor0 = localyhor0 = 20000; //GM for 22OCT07 data
672 
673  if (muonTags_.label() =="standAloneMuons") {
674 
675  Handle<CaloTowerCollection> calotower;
676  iEvent.getByToken(tok_tower_, calotower);
677 
678  for (CaloTowerCollection::const_iterator calt = calotower->begin();
679  calt !=calotower->end(); calt++) {
680  //CMSSW_2_1_x const math::XYZVector towermom = (*calt).momentum();
681  double ith = (*calt).momentum().theta();
682  double iph = (*calt).momentum().phi();
683 
684  CLHEP::Hep3Vector calo3v(sin(ith)*cos(iph), sin(ith)*sin(iph), cos(ith));
685 
686  double angle = tmpmuon3v.angle(calo3v);
687 
688  if (angle < 7.5*pival/180.) {tmpHOCalib.caloen[0] += calt->emEnergy()+calt->hadEnergy();}
689  if (angle < 15*pival/180.) {tmpHOCalib.caloen[1] += calt->emEnergy()+calt->hadEnergy();}
690  if (angle < 35*pival/180.) {tmpHOCalib.caloen[2] += calt->emEnergy()+calt->hadEnergy();}
691  }
692 
693 
694  }
695  if (tmpHOCalib.caloen[0] >10.0) continue;
696 
697  GlobalPoint glbpt(posx, posy, posz);
698 
699  double mom = sqrt(momx*momx + momy*momy +momz*momz);
700 
701  momx /= mom;
702  momy /= mom;
703  momz /= mom;
704 
705  DirectionType trkdir(momx, momy, momz);
706 
707  tmpHOCalib.trkdr = (*ncosm).d0();
708  tmpHOCalib.trkdz = (*ncosm).dz();
709 
710  tmpHOCalib.nmuon = cosmicmuon->size();
711  tmpHOCalib.trkvx = glbpt.x();
712  tmpHOCalib.trkvy = glbpt.y();
713  tmpHOCalib.trkvz = glbpt.z();
714  tmpHOCalib.trkmm = mom*charge;
715  tmpHOCalib.trkth = trkdir.theta();
716  tmpHOCalib.trkph = trkdir.phi();
717 
718  tmpHOCalib.ndof = (inearbymuon ==0) ? (int)(*ncosm).ndof() : -(int)(*ncosm).ndof();
719  tmpHOCalib.chisq = (*ncosm).normalizedChi2(); // max(1.,tmpHOCalib.ndof);
720  tmpHOCalib.therr = 0.;
721  tmpHOCalib.pherr = 0.;
722 
723  if (iiner==1) {
724  reco::TrackBase::CovarianceMatrix innercov = (*ncosm).innerStateCovariance();
725  tmpHOCalib.therr = innercov(1,1); //thetaError();
726  tmpHOCalib.pherr = innercov(2,2); //phi0Error();
727  } else {
728  reco::TrackBase::CovarianceMatrix outercov = (*ncosm).outerStateCovariance();
729  tmpHOCalib.therr = outercov(1,1); //thetaError();
730  tmpHOCalib.pherr = outercov(2,2); //phi0Error();
731  }
732 
733  ESHandle<MagneticField> theMagField;
734  iSetup.get<IdealMagneticFieldRecord>().get(theMagField );
735 
736  SteppingHelixPropagator myHelix(&*theMagField,anyDirection);
737  myHelix.setMaterialMode(false);
738  myHelix.applyRadX0Correction(true);
739 
740  double phiho = trkpos.phi();
741  if (phiho<0) phiho +=2*pival;
742 
743  int iphisect_dt=int(6*(phiho+pival/18.)/pival); //for u 18/12/06
744  if (iphisect_dt>=12) iphisect_dt=0;
745 
746  int iphisect = -1;
747 
748  bool ipath = false;
749  for (int kl = 0; kl<=2; kl++) {
750 
751  int iphisecttmp = (kl<2) ? iphisect_dt + kl : iphisect_dt - 1;
752  if (iphisecttmp <0) iphisecttmp = 11;
753  if (iphisecttmp >=12) iphisecttmp = 0;
754 
755  double phipos = iphisecttmp*pival/6.;
756  double phirot = phipos;
757 
758  GlobalVector xLocal(-sin(phirot), cos(phirot), 0.);
759 
760  GlobalVector yLocal(0., 0., 1.);
761  GlobalVector zLocal = xLocal.cross(yLocal).unit();
762  // GlobalVector zLocal(cos(phirot), sin(phirot), 0.0);
763 
764 
765  FreeTrajectoryState freetrajectorystate_ = getFreeTrajectoryState(*ncosm,&(*theMagField), iiner, samedir);
766 
767  Surface::RotationType rot(xLocal, yLocal, zLocal);
768 
769  for (int ik=1; ik>=0; ik--) { //propagate track in two HO layers
770 
771  double radial = 407.0;
772  if (ik==0) radial = 382.0;
773 
774  Surface::PositionType pos(radial*cos(phipos), radial*sin(phipos), 0.);
775  PlaneBuilder::ReturnType aPlane = PlaneBuilder().plane(pos,rot);
776 
777  auto aPlane2 = new Plane(pos,rot);
778 
779  SteppingHelixStateInfo steppingHelixstateinfo_ = myHelix.propagate(SteppingHelixStateInfo(freetrajectorystate_), (*aPlane2));
780 
781  if (steppingHelixstateinfo_.isValid()) {
782 
783  GlobalVector hotrkpos2(steppingHelixstateinfo_.position().x(), steppingHelixstateinfo_.position().y(), steppingHelixstateinfo_.position().z());
784  CLHEP::Hep3Vector hotrkdir2(steppingHelixstateinfo_.momentum().x(), steppingHelixstateinfo_.momentum().y(),steppingHelixstateinfo_.momentum().z());
785 
786  LocalVector lclvt0 = (*aPlane).toLocal(hotrkpos2);
787 
788  double xx = lclvt0.x();
789  double yy = lclvt0.y();
790 
791  if (ik ==1) {
792  if ((std::abs(yy) < 130 && xx >-64.7 && xx <138.2)
793  ||(std::abs(yy) > 130 && std::abs(yy) <700 && xx >-76.3 && xx <140.5)) {
794  ipath = true; //Only look for tracks which as hits in layer 1
795  iphisect = iphisecttmp;
796  }
797  }
798 
799  if (iphisect != iphisecttmp) continue; //Look for ring-0 only when ring1 is accepted for that sector
800 
801  switch (ik)
802  {
803  case 0 :
804  xhor0 = xx; //lclvt0.x();
805  yhor0 = yy; //lclvt0.y();
806  break;
807  case 1 :
808  xhor1 = xx; //lclvt0.x();
809  yhor1 = yy; //lclvt0.y();
810 
811  tmpHOCalib.hoang = CLHEP::Hep3Vector(zLocal.x(),zLocal.y(),zLocal.z()).dot(hotrkdir2.unit());
812  break;
813  default : break;
814  }
815  } else {
816  break;
817  }
818  }
819  if (ipath) break;
820  }
821  if (ipath) { //If muon crossed HO laeyrs
822 
823  int ietaho = 50;
824  int iphiho = -1;
825 
826  for (int i=0; i<9; i++) {tmpHOCalib.hosig[i]=-100.0;}
827  for (int i=0; i<18; i++) {tmpHOCalib.hocorsig[i]=-100.0;}
828  for (int i=0; i<9; i++) {tmpHOCalib.hbhesig[i]=-100.0;}
829  tmpHOCalib.hocro = -100;
830  tmpHOCalib.htime = -1000;
831 
832  int isect = 0;
833 
834  findHOEtaPhi(iphisect, ietaho, iphiho);
835 
836  if (ietaho !=0 && iphiho !=0 && std::abs(iring)<=2) { //Muon passed through a tower
837  isect = 100*std::abs(ietaho+30)+std::abs(iphiho);
838  if (std::abs(ietaho) >=netabin || iphiho<0) isect *=-1; //Not extrapolated to any tower
839  if (std::abs(ietaho) >=netabin) isect -=1000000; //not matched with eta
840  if (iphiho<0) isect -=2000000; //not matched with phi
841  tmpHOCalib.isect = isect;
842 
843  tmpHOCalib.hodx = localxhor1;
844  tmpHOCalib.hody = localyhor1;
845 
846  if (iring==0) {
847  tmpHOCalib.hocorsig[8] = localxhor0;
848  tmpHOCalib.hocorsig[9] = localyhor0;
849  }
850 
851  int etamn=-4;
852  int etamx=4;
853  if (iring==1) {etamn=5; etamx = 10;}
854  if (iring==2) {etamn=11; etamx = 16;}
855  if (iring==-1){etamn=-10; etamx = -5;}
856  if (iring==-2){etamn=-16; etamx = -11;}
857 
858  int phimn = 1;
859  int phimx = 2;
860  if (iring ==0) {
861  phimx =2*int((iphiho+1)/2.);
862  phimn = phimx - 1;
863  } else {
864  phimn = 3*int((iphiho+1)/3.) - 1;
865  phimx = phimn + 2;
866  }
867 
868  if (phimn <1) phimn += nphimx;
869  if (phimx >72) phimx -= nphimx;
870 
871  int sigstr = m_startTS; // 5;
872  int sigend = m_endTS; // 8;
873 
874  // if (iphiho <=nphimx/2) { //GMA310508
875  // sigstr -=1; //GMA300608
876  // sigend -=1;
877  // }
878 
879  if (m_hbinfo) {
880  for (int i=0; i<9; i++) {tmpHOCalib.hbhesig[i]=-100.0;}
881 
882  if (m_digiInput) {
883  if ((*hbhe).size() >0) {
884  for (HBHEDigiCollection::const_iterator j=(*hbhe).begin(); j!=(*hbhe).end(); j++){
885  // const HBHEDataFrame digi = (const HBHEDataFrame)(*j);
886  // HcalDetId id =digi.id();
887  HcalDetId id =(*j).id();
888  m_coder = (*conditions_).getHcalCoder(id);
889  m_shape = (*conditions_).getHcalShape(m_coder);
890  int tmpeta= id.ieta();
891  int tmpphi= id.iphi();
892  calibped = conditions_->getHcalCalibrations(id);
893 
894  int deta = tmpeta-ietaho;
895  if (tmpeta==-1 && ietaho== 1) deta = -1;
896  if (tmpeta== 1 && ietaho==-1) deta = 1;
897  int dphi = tmpphi-iphiho;
898  if (phimn >phimx) {
899  if (dphi==71) dphi=-1;
900  if (dphi==-71) dphi=1;
901  }
902 
903  int ipass2 = (std::abs(deta) <=1 && std::abs(dphi)<=1) ? 1 : 0; //NEED correction in full CMS detector
904 
905  if (ipass2 ==0 ) continue;
906 
907  float tmpdata[nchnmx];
908  for (int i=0; i<(*j).size() && i<nchnmx; i++) {
909  tmpdata[i] = m_coder->charge(*m_shape,(*j).sample(i).adc(),(*j).sample(i).capid());
910  }
911 
912  float signal = 0;
913  for (int i=1; i<(*j).size() && i<=8; i++) {
914  signal += tmpdata[i] - calibped.pedestal((*j).sample(i).capid());;
915  }
916 
917  if (ipass2 == 1) {
918  if (3*(deta+1)+dphi+1<9) tmpHOCalib.hbhesig[3*(deta+1)+dphi+1] = signal;
919  }
920  }
921  }
922 
923  } else {
924 
925  edm::Handle<HBHERecHitCollection> hbheht;// iEvent.getByType(hbheht);
926  iEvent.getByToken(tok_hbhe_,hbheht);
927 
928 
929  if ((*hbheht).size()>0) {
930  if(!(*hbheht).size()) throw (int)(*hbheht).size();
931 
932  for (HBHERecHitCollection::const_iterator j=(*hbheht).begin(); j!=(*hbheht).end(); j++){
933  // const HBHERecHit hbhehtrec = (const HBHERecHit)(*j);
934  // HcalDetId id =hbhehtrec.id();
935  HcalDetId id =(*j).id();
936  int tmpeta= id.ieta();
937  int tmpphi= id.iphi();
938 
939  int deta = tmpeta-ietaho;
940  if (tmpeta==-1 && ietaho== 1) deta = -1;
941  if (tmpeta== 1 && ietaho==-1) deta = 1;
942  int dphi = tmpphi-iphiho;
943  if (phimn >phimx) {
944  if (dphi==71) dphi=-1;
945  if (dphi==-71) dphi=1;
946  }
947 
948  int ipass2 = (std::abs(deta) <=1 && std::abs(dphi)<=1) ? 1 : 0; //NEED correction in full CMS detector
949  if ( ipass2 ==0 ) continue;
950 
951  float signal = (*j).energy();
952 
953  if (3*(deta+1)+dphi+1<9) tmpHOCalib.hbhesig[3*(deta+1)+dphi+1] = signal;
954  }
955  }
956 
957  } //else m_digilevel
958 
959  } //m_hbinfo #endif
960 
961  if (m_digiInput) {
962  if ((*ho).size()>0) {
963  int isFilled[netamx*nphimx];
964  for (int j=0; j<netamx*nphimx; j++) {isFilled[j]=0;}
965 
966  double sumEt = 0;
967  double sumE = 0;
968 
969  for (HODigiCollection::const_iterator j=(*ho).begin(); j!=(*ho).end(); j++){
970  // const HODataFrame digi = (const HODataFrame)(*j);
971  // HcalDetId id =digi.id();
972 
973  HcalDetId id =(*j).id();
974  m_coder = (*conditions_).getHcalCoder(id);
975  m_shape = (*conditions_).getHcalShape(m_coder);
976  int tmpeta= id.ieta();
977  int tmpphi= id.iphi();
978 
979  int ipass1 =0;
980  if (tmpeta >=etamn && tmpeta <=etamx) {
981  if (phimn < phimx) {
982  ipass1 = (tmpphi >=phimn && tmpphi <=phimx ) ? 1 : 0;
983  } else {
984  ipass1 = (tmpphi==71 || tmpphi ==72 || tmpphi==1) ? 1 : 0;
985  }
986  }
987 
988  int deta = tmpeta-ietaho;
989  if (tmpeta==-1 && ietaho== 1) deta = -1;
990  if (tmpeta== 1 && ietaho==-1) deta = 1;
991 
992  int dphi = tmpphi -iphiho;
993  if (phimn>phimx) {
994  if (dphi==71) dphi=-1;
995  if (dphi==-71) dphi=1;
996  }
997 
998  int ipass2 = (std::abs(deta) <=1 && std::abs(dphi)<=1) ? 1 : 0;
999 
1000  int tmpeta1 = (tmpeta>0) ? tmpeta -1 : -tmpeta +14;
1001 
1002  float tmpdata[nchnmx]={0,0,0,0,0,0,0,0,0,0};
1003  float sigvall[nsigpk]={0,0,0,0,0,0,0};
1004 
1005  for (int i=0; i<(*j).size() && i<nchnmx; i++) {
1006  tmpdata[i] = m_coder->charge(*m_shape,(*j).sample(i).adc(),(*j).sample(i).capid());
1007  if (deta==0 && dphi==0) {
1008  double tmpE = tmpdata[i] - pedestal[tmpeta1][tmpphi-1][(*j).sample(i).capid()];
1009  if (tmpE >0) {
1010  sumEt += i*tmpE;
1011  sumE += tmpE;
1012  }
1013  if (m_hotime) {
1014  //calculate signals in 4 time slices, 0-3,.. 6-9
1015  if (i>=7-nsigpk) {
1016  for (int ncap=0; ncap<nsigpk; ncap++) {
1017  if (i-ncap >= nstrbn && i-ncap <= nstrbn+3) {
1018  sigvall[ncap] +=tmpdata[i];
1019  }
1020  }
1021  }
1022  if (i==(*j).size()-1) {
1023  float mxled=-1;
1024  int imxled = 0;
1025  for (int ij=0; ij<nsigpk; ij++) {
1026  if (sigvall[ij] > mxled) {mxled = sigvall[ij]; imxled=ij;}
1027  }
1028  double pedx = 0.0;
1029  for (int ij=0; ij<4; ij++) {
1030  pedx +=pedestal[tmpeta1][tmpphi-1][ij];
1031  }
1032  if (mxled-pedx >2 && mxled-pedx <20 ) {
1033  hopeak[ntrgp_gm]->Fill(nphimx*tmpeta1 + tmpphi-1, imxled+nstrbn);
1034  for (int jk=0; jk<ntrgp_gm; jk++) {
1035  if (ntrgpas_gm[jk]>0) {
1036  hopeak[jk]->Fill(nphimx*tmpeta1 + tmpphi-1, imxled+nstrbn);
1037  }
1038  }
1039  if (tmpdata[5]+tmpdata[6] >1) {
1040  horatio->Fill(nphimx*tmpeta1 + tmpphi-1, (tmpdata[5]-tmpdata[6])/(tmpdata[5]+tmpdata[6]));
1041  }
1042  for (int ij=0; ij<(*j).size() && ij<nchnmx; ij++) {
1043  hotime[ntrgp_gm]->Fill(nphimx*nchnmx*tmpeta1 + nchnmx*(tmpphi-1) + ij, tmpdata[ij]);
1044  Nhotime[ntrgp_gm]->Fill(nphimx*nchnmx*tmpeta1 + nchnmx*(tmpphi-1) + ij, 1.);
1045  for (int jk=0; jk<ntrgp_gm; jk++) {
1046  if (ntrgpas_gm[jk]>0) {
1047  hotime[jk]->Fill(nphimx*nchnmx*tmpeta1 + nchnmx*(tmpphi-1) + ij, tmpdata[ij]);
1048  Nhotime[jk]->Fill(nphimx*nchnmx*tmpeta1 + nchnmx*(tmpphi-1) + ij, 1.);
1049  }
1050  }
1051  }
1052  }
1053  }
1054  }
1055  }
1056  }
1057 
1058  if (std::abs(tmpeta) <=15 && deta==0 && dphi ==0) {
1059  float signal = 0;
1060  size_t icnt = 0;
1061  for (int i =0; i<nchnmx && i< (*j).size(); i++) {
1062  if (i >=sigstr && i<=sigend) continue;
1063  signal += tmpdata[i] - pedestal[tmpeta1][tmpphi-1][(*j).sample(i).capid()];
1064  if (++icnt >=4) break;
1065  }
1066  tmpHOCalib.hocro = signal;
1067  }
1068 
1069  if (m_hotime) {
1070  if (ipass1 ==0 && ipass2 ==0 && cosmicmuon->size()<=2) {
1071  if (std::abs(ietaho) <=netabin && iphiho >0) {
1072  if ((iphiho >=1 && iphiho<=nphimx/2 && tmpphi >=1 && tmpphi <=nphimx/2) ||
1073  (iphiho >nphimx/2 && iphiho<=nphimx && tmpphi >nphimx/2 && tmpphi <=nphimx)) {
1074  if (isFilled[nphimx*tmpeta1+tmpphi-1]==0) {
1075  isFilled[nphimx*tmpeta1+tmpphi-1]=1;
1076  for (int i=0; i<(*j).size() && i<nchnmx; i++) {
1077  hopedtime->Fill(nphimx*nchnmx*tmpeta1 + nchnmx*(tmpphi-1) + i, tmpdata[i]);
1078  Nhopedtime->Fill(nphimx*nchnmx*tmpeta1 + nchnmx*(tmpphi-1) + i, 1.);
1079  hopedpr->Fill(nphimx*nchnmx*tmpeta1 + nchnmx*(tmpphi-1) + i, tmpdata[i]);
1080  }
1081  } //isFilled
1082  }
1083  }
1084  }
1085  }
1086 
1087  if (ipass1 ==0 && ipass2 ==0 ) continue;
1088 
1089  float signal = 0;
1090  for (int i=sigstr; i<(*j).size() && i<=sigend; i++) {
1091  signal += tmpdata[i] - pedestal[tmpeta1][tmpphi-1][(*j).sample(i).capid()];
1092  }
1093  if (signal <-100 || signal >100000) signal = -100;
1094  if (m_hotime) {
1095  if (signal >-100 && Noccu == Noccu_old) {
1096  for (int i=0; i<5; i++) {
1097  if (signal >(i+2)*m_sigma) {
1098  ho_occupency[i]->Fill(nphimx*tmpeta1+tmpphi-1);
1099  }
1100  }
1101  }
1102  }
1103 
1104  if (ipass1 ==0 && ipass2 ==0 ) continue;
1105 
1106  if (ipass1 ==1) {
1107  int tmpdph = tmpphi-phimn;
1108  if (tmpdph<0) tmpdph = 2; //only case of iphi==1, where phimn=71
1109 
1110  int ilog = 2*(tmpeta-etamn)+tmpdph;
1111  if (iring !=0) {
1112  if (iring >0) {
1113  ilog = 3*(tmpeta-etamn)+tmpdph; //Again CMS correction
1114  } else {
1115  ilog = 3*(etamx-tmpeta)+tmpdph; //Again CMS correction
1116  }
1117  }
1118  if (ilog>-1 && ilog<18) {
1119  tmpHOCalib.hocorsig[ilog] = signal;
1120  }
1121  }
1122 
1123  if (ipass2 ==1) {
1124  if (3*(deta+1)+dphi+1<9) tmpHOCalib.hosig[3*(deta+1)+dphi+1] = signal; //Again CMS azimuthal near phi 1&72
1125  }
1126 
1127  /*
1128  // Possibility to store pedestal by shifting phi tower by 6
1129  // But, due to missing tower at +-5, we do not have always proper
1130  // statistics and also in pedestal subtracted data, we do not have
1131  // signal in that tower
1132  //
1133  if (deta==0 && dphi ==0) {
1134  int crphi = tmpphi + 6;
1135  if (crphi >72) crphi -=72;
1136 
1137  for (HODigiCollection::const_iterator jcr=(*ho).begin(); jcr!=(*ho).end(); jcr++){
1138  // const HODataFrame (*jcr) = (const HODataFrame)(*jcr);
1139  // HcalDetId idcr =(*jcr).id();
1140  HcalDetId id =(*jcr).id();
1141  m_coder = (*conditions_).getHcalCoder(idcr);
1142  m_shape = (*conditions_).getHcalShape(m_coder);
1143  int etacr= idcr.ieta();
1144  int phicr= idcr.iphi();
1145 
1146  if (tmpeta==etacr && crphi ==phicr) {
1147 
1148  float tmpdatacr[nchnmx];
1149  for (int i=0; i<(*jcr).size() && i<nchnmx; i++) {
1150  tmpdatacr[i] = m_coder->charge(*m_shape,(*jcr).sample(i).adc(),(*jcr).sample(i).capid());
1151  }
1152  }
1153  }
1154  }
1155  */
1156 
1157  }
1158  tmpHOCalib.htime = sumEt/max(sumE,1.e-6);
1159  }
1160  } else {
1162  iEvent.getByToken(tok_ho_,hoht);
1163 
1164 
1165  if ((*hoht).size()>0) {
1166  for (HORecHitCollection::const_iterator j=(*hoht).begin(); j!=(*hoht).end(); j++){
1167  // const HORecHit hohtrec = (const HORecHit)(*j);
1168  // HcalDetId id =hohtrec.id();
1169  HcalDetId id =(*j).id();
1170  int tmpeta= id.ieta();
1171  int tmpphi= id.iphi();
1172 
1173  int ipass1 =0;
1174  if (tmpeta >=etamn && tmpeta <=etamx) {
1175  if (phimn < phimx) {
1176  ipass1 = (tmpphi >=phimn && tmpphi <=phimx ) ? 1 : 0;
1177  } else {
1178  ipass1 = (tmpphi==71 || tmpphi ==72 || tmpphi==1) ? 1 : 0;
1179  }
1180  }
1181 
1182  int deta = tmpeta-ietaho;
1183  if (tmpeta==-1 && ietaho== 1) deta = -1;
1184  if (tmpeta== 1 && ietaho==-1) deta = 1;
1185 
1186  int dphi = tmpphi -iphiho;
1187  if (phimn>phimx) {
1188  if (dphi==71) dphi=-1;
1189  if (dphi==-71) dphi=1;
1190  }
1191 
1192  float signal = (*j).energy();
1193  if (m_hotime) {
1194  int tmpeta1 = (tmpeta>0) ? tmpeta -1 : -tmpeta +14;
1195  if (signal >-100 && Noccu == Noccu_old) {
1196  for (int i=0; i<5; i++) {
1197  if (signal >(i+2)*m_sigma) {
1198  ho_occupency[i]->Fill(nphimx*tmpeta1+tmpphi-1);
1199  }
1200  }
1201  }
1202  }
1203 
1204  int ipass2 = (std::abs(deta) <=1 && std::abs(dphi)<=1) ? 1 : 0;
1205 
1206  if (ipass1 ==0 && ipass2 ==0 ) continue;
1207 
1208  if (ipass1 ==1) {
1209  int tmpdph = tmpphi-phimn;
1210  if (tmpdph<0) tmpdph = 2; //only case of iphi==1, where phimn=71
1211 
1212  int ilog = 2*(tmpeta-etamn)+tmpdph;
1213  if (iring !=0) {
1214  if (iring >0) {
1215  ilog = 3*(tmpeta-etamn)+tmpdph; //Again CMS correction
1216  } else {
1217  ilog = 3*(etamx-tmpeta)+tmpdph; //Again CMS correction
1218  }
1219  }
1220  if (ilog>-1 && ilog<18) {
1221  tmpHOCalib.hocorsig[ilog] = signal;
1222  }
1223  }
1224 
1225  if (ipass2 ==1) {
1226 
1227  if (3*(deta+1)+dphi+1<9) {
1228  tmpHOCalib.hosig[3*(deta+1)+dphi+1] = signal; //Again CMS azimuthal near phi 1&72
1229  }
1230  }
1231 
1232  if (deta==0 && dphi ==0) {
1233  tmpHOCalib.htime = (*j).time();
1234  int crphi = tmpphi + 6;
1235  if (crphi >72) crphi -=72;
1236 
1237  for (HORecHitCollection::const_iterator jcr=(*hoht).begin(); jcr!=(*hoht).end(); jcr++){
1238  const HORecHit reccr = (const HORecHit)(*jcr);
1239  HcalDetId idcr =reccr.id();
1240  int etacr= idcr.ieta();
1241  int phicr= idcr.iphi();
1242  if (tmpeta==etacr && crphi ==phicr) {
1243 
1244  tmpHOCalib.hocro = reccr.energy();
1245 
1246  }
1247  }
1248  }
1249  }
1250  }
1251  }
1252 
1253  //GMA Npass++;
1254  if (Noccu == Noccu_old) Noccu++;
1255  hostore->push_back(tmpHOCalib);
1256 
1257  }
1258  }
1259 
1260  }
1261  }
1262 
1263  iEvent.put(hostore, "HOCalibVariableCollection");
1264 
1265 }
1266 
1267 // ------------ method called once each job just before starting event loop ------------
1268 void
1270 {
1271  //GMA Nevents = 0;
1272  //GMA Npass = 0;
1273  //GMA Noccu = 0;
1274 
1275  irunold = -1;
1276  nRuns = 0;
1277  // edm::ESHandle<MagneticField> bField;
1278  // iSetup.get<IdealMagneticFieldRecord>().get(bField);
1279  // stepProp = new SteppingHelixPropagator(&*bField,anyDirection);
1280  // stepProp->setMaterialMode(false);
1281  // stepProp->applyRadX0Correction(true);
1282 
1283  for (int i=0; i<netamx; i++) {
1284  for (int j=0; j<nphimx; j++) {
1285  for (int k=0; k<ncidmx; k++) {
1286  pedestal[i][j][k]=0.0;
1287  }
1288  }
1289  }
1290 
1291 
1292 }
1293 
1294 // ------------ method called once each job just after ending the event loop ------------
1295 void
1297 
1298 
1299 }
1300 
1301 void AlCaHOCalibProducer::findHOEtaPhi(int iphisect, int& ietaho, int& iphiho) {
1302 
1303  //18/12/06 : use only position, not angle phi
1304 
1305 double etalow[netabin]={ 0.025, 35.195, 70.625, 106.595, 141.565, 180.765, 220.235, 261.385, 304.525, 349.975, 410.025, 452.085, 506.645, 565.025, 627.725, 660.25};
1306 double etahgh[netabin]={ 35.145, 70.575, 106.545, 125.505, 180.715, 220.185, 261.335, 304.475, 349.925, 392.575, 452.035, 506.595, 564.975, 627.675, 661.075, 700.25};
1307 
1308  double philow[6]={-76.27, -35.11, 0.35, 35.81, 71.77, 108.93}; //Ring+/-1 & 2
1309  double phihgh[6]={-35.81, -0.35, 35.11, 71.07, 108.23, 140.49};
1310 
1311  double philow00[6]={-60.27, -32.91, 0.35, 33.61, 67.37, 102.23}; //Ring0 L0
1312  double phihgh00[6]={-33.61, -0.35, 32.91, 66.67, 101.53, 129.49};
1313 
1314  double philow01[6]={-64.67, -34.91, 0.35, 35.61, 71.37, 108.33}; //Ring0 L1
1315  double phihgh01[6]={-35.61, -0.35, 34.91, 70.67, 107.63, 138.19};
1316 
1317 
1318  iring = -10;
1319 
1320  double tmpdy = std::abs(yhor1);
1321  for (int i=0; i<netabin; i++) {
1322  if (tmpdy >etalow[i] && tmpdy <etahgh[i]) {
1323  ietaho = i+1;
1324  float tmp1 = fabs(tmpdy-etalow[i]);
1325  float tmp2 = fabs(tmpdy-etahgh[i]);
1326 
1327  localyhor1 = (tmp1 < tmp2) ? -tmp1 : tmp2;
1328 
1329  if (i<4) iring =0;
1330  if (i>=4 && i<10) iring=1;
1331  if (i>=10 && i<netabin) iring=2;
1332  break;
1333  }
1334  }
1335 
1336  int tmpphi = 0;
1337  int tmpphi0 = 0;
1338 
1339  if (ietaho >4) { //Ring 1 and 2
1340  for (int i=0; i<6; i++) {
1341  if (xhor1 >philow[i] && xhor1 <phihgh[i]) {
1342  tmpphi=i+1;
1343  float tmp1 = fabs(xhor1-philow[i]);
1344  float tmp2 = fabs(xhor1-phihgh[i]);
1345  localxhor1 = (tmp1 < tmp2) ? -tmp1 : tmp2;
1346  break;
1347  }
1348  }
1349  } else { //Ring 0
1350  for (int i=0; i<6; i++) {
1351  if (xhor1 >philow01[i] && xhor1 <phihgh01[i]) {
1352  tmpphi=i+1;
1353  float tmp1 = fabs(xhor1-philow01[i]);
1354  float tmp2 = fabs(xhor1-phihgh01[i]);
1355  localxhor1 = (tmp1 < tmp2) ? -tmp1 : tmp2;
1356  break;
1357  }
1358  }
1359 
1360  for (int i=0; i<6; i++) {
1361  if (xhor0 >philow00[i] && xhor0 <phihgh00[i]) {
1362  tmpphi0=i+1;
1363  float tmp1 = fabs(xhor0-philow00[i]);
1364  float tmp2 = fabs(xhor0-phihgh00[i]);
1365  localxhor0 = (tmp1 < tmp2) ? -tmp1 : tmp2;
1366  if (tmpphi !=tmpphi0) localxhor0 +=10000.;
1367  break;
1368  }
1369  }
1370 
1371  double tmpdy = std::abs(yhor0);
1372  for (int i=0; i<4; i++) {
1373  if (tmpdy >etalow[i] && tmpdy <etahgh[i]) {
1374  float tmp1 = fabs(tmpdy-etalow[i]);
1375  float tmp2 = fabs(tmpdy-etahgh[i]);
1376  localyhor0 = (tmp1 < tmp2) ? -tmp1 : tmp2;
1377  if (i+1 != ietaho) localyhor0 +=10000.;
1378  break;
1379  }
1380  }
1381  }
1382 
1383  if (tmpphi!=0) {
1384  iphiho = 6*iphisect -2 + tmpphi;
1385  if (iphiho <=0) iphiho +=nphimx;
1386  if (iphiho >nphimx) iphiho -=nphimx;
1387  }
1388 
1389  // isect2 = 15*iring+iphisect+1;
1390 
1391  if (yhor1 <0) {
1392  if (std::abs(ietaho) >netabin) { //Initialised with 50
1393  ietaho +=1;
1394  } else {
1395  ietaho *=-1;
1396  }
1397  // isect2 *=-1;
1398  iring *=-1;
1399  }
1400 }
1401 
1403 {
1404 
1405  if (iiner ==0) {
1406  GlobalPoint gpos( tk.outerX(), tk.outerY(), tk.outerZ());
1407  GlobalVector gmom( tk.outerPx(), tk.outerPy(), tk.outerPz());
1408  if (dir) gmom *=-1.;
1409  GlobalTrajectoryParameters par( gpos, gmom, tk.charge(), field);
1410  CurvilinearTrajectoryError err( tk.extra()->outerStateCovariance());
1411  return FreeTrajectoryState( par, err);
1412  } else {
1413  GlobalPoint gpos( tk.innerPosition().X(), tk.innerPosition().Y(), tk.innerPosition().Z());
1414  GlobalVector gmom( tk.innerMomentum().X(), tk.innerMomentum().Y(), tk.innerMomentum().Z());
1415  if (dir) gmom *=-1.;
1416  GlobalTrajectoryParameters par( gpos, -gmom, tk.charge(), field);
1417  CurvilinearTrajectoryError err( tk.extra()->innerStateCovariance());
1418  return FreeTrajectoryState( par, err);
1419  }
1420 
1421 }
1422 
1424 
1425 //define this as a plug-in
1427 
RunNumber_t run() const
Definition: EventID.h:42
T getParameter(std::string const &) const
Basic3DVector< float > DirectionType
T getUntrackedParameter(std::string const &, T const &) const
int i
Definition: DBlmapReader.cc:9
edm::EDGetTokenT< CaloTowerCollection > tok_tower_
FreeTrajectoryState getFreeTrajectoryState(const reco::Track &tk, const MagneticField *field, int itag, bool dir)
edm::EDGetTokenT< reco::TrackCollection > tok_muons_
const int netabin
double outerPy() const
y coordinate of momentum vector at the outermost hit position
Definition: Track.h:72
const TrackExtraRef & extra() const
reference to &quot;extra&quot; object
Definition: Track.h:96
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:434
HcalCalibrations calibped
const int ncidmx
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
Geom::Theta< T > theta() const
ReturnType plane(const PositionType &pos, const RotationType &rot) const
Definition: PlaneBuilder.h:22
std::vector< HODataFrame >::const_iterator const_iterator
edm::ESHandle< HcalDbService > conditions_
T y() const
Definition: PV3DBase.h:63
virtual void produce(edm::Event &, const edm::EventSetup &) override
std::map< std::string, bool > fired
T * make(const Args &...args) const
make new ROOT object
Definition: TFileService.h:64
ErrorD< N >::type type
Definition: Error.h:29
Geom::Phi< T > phi() const
const HcalQIEShape * m_shape
virtual void endJob() override
GlobalVector momentum() const
void findHOEtaPhi(int iphsect, int &ietaho, int &iphiho)
double charge(const std::vector< uint8_t > &Ampls)
Definition: Plane.h:17
tuple field
Definition: statics.py:62
void beginJob()
Definition: Breakpoints.cc:15
double outerZ() const
z coordinate of the outermost hit position
Definition: Track.h:80
static const unsigned int nL1mx
edm::EDGetTokenT< edm::TriggerResults > tok_hlt_
const math::XYZPoint & innerPosition() const
position of the innermost hit
Definition: Track.h:41
const int netahbmx
const HcalQIECoder * m_coder
int iEvent
Definition: GenABIO.cc:243
const T & max(const T &a, const T &b)
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
Definition: Event.h:116
float energy() const
Definition: CaloRecHit.h:17
GlobalPoint position() const
T sqrt(T t)
Definition: SSEVec.h:48
static const unsigned int nL1trg
Basic3DVector< float > RotationType
Vector3DBase< typename PreciseFloatType< T, U >::Type, FrameTag > cross(const Vector3DBase< U, FrameTag > &v) const
Definition: Vector3DBase.h:119
T z() const
Definition: PV3DBase.h:64
const int nsigpk
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
int ieta() const
get the cell ieta
Definition: HcalDetId.h:36
math::Error< 5 >::type CovarianceMatrix
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
int j
Definition: DBlmapReader.cc:9
virtual void beginJob() override
Basic3DVector< float > PositionType
const int nstrbn
double outerX() const
x coordinate of the outermost hit position
Definition: Track.h:76
const int ntrgp_gm
double outerPz() const
z coordinate of momentum vector at the outermost hit position
Definition: Track.h:74
int k[5][pyjets_maxn]
edm::EDGetTokenT< L1GlobalTriggerReadoutRecord > tok_l1_
int iphi() const
get the cell iphi
Definition: HcalDetId.h:38
const T & get() const
Definition: EventSetup.h:55
std::vector< HOCalibVariables > HOCalibVariableCollection
collection of HOcalibration variabale
edm::EDGetTokenT< HORecHitCollection > tok_ho_
edm::EventID id() const
Definition: EventBase.h:56
AlCaHOCalibProducer(const edm::ParameterSet &)
T dot(const Basic3DVector &v) const
Scalar product, or &quot;dot&quot; product, with a vector of same type.
const int nphimx
const math::XYZVector & innerMomentum() const
momentum vector at the innermost hit position
Definition: Track.h:44
const int netamx
double outerY() const
y coordinate of the outermost hit position
Definition: Track.h:78
HcalCalibrationWidths calibwidth
int charge() const
track electric charge
Definition: TrackBase.h:111
const int nchnmx
dbl *** dir
Definition: mlp_gen.cc:35
T x() const
Definition: PV3DBase.h:62
const int netahb3mx
static const unsigned int nHLTmx
HcalDetId id() const
get the id
Definition: HORecHit.h:19
double outerPx() const
x coordinate of momentum vector at the outermost hit position
Definition: Track.h:70
edm::EDGetTokenT< HBHERecHitCollection > tok_hbhe_
math::Error< dimension >::type CovarianceMatrix
5 parameter covariance matrix
Definition: TrackBase.h:70
T angle(T x1, T y1, T z1, T x2, T y2, T z2)
Definition: angle.h:11