File indexing completed on 2021-10-04 02:55:06
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066
0067
0068
0069
0070
0071
0072
0073
0074
0075 #include <memory>
0076
0077
0078 #include "FWCore/Framework/interface/Frameworkfwd.h"
0079 #include "FWCore/Framework/interface/one/EDProducer.h"
0080 #include "FWCore/Framework/interface/Event.h"
0081 #include "FWCore/Framework/interface/EventSetup.h"
0082 #include "FWCore/Framework/interface/MakerMacros.h"
0083 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0084 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0085 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0086 #include "FWCore/Utilities/interface/InputTag.h"
0087
0088 #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
0089 #include "DataFormats/HcalCalibObjects/interface/HOCalibVariables.h"
0090 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0091 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0092 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0093 #include "DataFormats/GeometrySurface/interface/PlaneBuilder.h"
0094 #include "DataFormats/Luminosity/interface/LumiDetails.h"
0095 #include "DataFormats/Math/interface/Error.h"
0096 #include "DataFormats/MuonReco/interface/Muon.h"
0097 #include "DataFormats/RecoCandidate/interface/IsoDeposit.h"
0098 #include "DataFormats/Scalers/interface/LumiScalers.h"
0099 #include "DataFormats/TrackReco/interface/Track.h"
0100 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0101 #include "DataFormats/TrajectorySeed/interface/PropagationDirection.h"
0102 #include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
0103 #include "DataFormats/VertexReco/interface/Vertex.h"
0104 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0105
0106 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0107 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0108 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0109 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0110 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0111 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0112
0113 #include "MagneticField/Engine/interface/MagneticField.h"
0114 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0115
0116 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0117 #include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h"
0118
0119 #include "FWCore/ServiceRegistry/interface/Service.h"
0120 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0121
0122
0123
0124 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalSeverityLevelComputer.h"
0125 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalSeverityLevelComputerRcd.h"
0126 #include "CondFormats/HcalObjects/interface/HcalChannelQuality.h"
0127 #include "CondFormats/DataRecord/interface/HcalChannelQualityRcd.h"
0128 #include "CLHEP/Vector/LorentzVector.h"
0129 #include "CLHEP/Units/GlobalPhysicalConstants.h"
0130 #include "CLHEP/Units/GlobalSystemOfUnits.h"
0131
0132 #include "TH2F.h"
0133
0134
0135 #include <string>
0136
0137 #include <iostream>
0138 #include <fstream>
0139
0140
0141
0142
0143 class AlCaHOCalibProducer : public edm::one::EDProducer<edm::one::SharedResources> {
0144 public:
0145 explicit AlCaHOCalibProducer(const edm::ParameterSet&);
0146 ~AlCaHOCalibProducer() override = default;
0147
0148 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0149 typedef Basic3DVector<float> PositionType;
0150 typedef Basic3DVector<float> DirectionType;
0151 typedef Basic3DVector<float> RotationType;
0152
0153 private:
0154 void produce(edm::Event&, const edm::EventSetup&) override;
0155 void beginJob() override;
0156 void endJob() override;
0157
0158 void fillHOStore(const reco::TrackRef& ncosm,
0159 HOCalibVariables& tmpHOCalib,
0160 std::unique_ptr<HOCalibVariableCollection>& hostore,
0161 int Noccu_old,
0162 int indx,
0163 edm::Handle<reco::TrackCollection> cosmicmuon,
0164 edm::View<reco::Muon>::const_iterator muon1,
0165 const edm::Event& iEvent,
0166 const CaloSubdetectorGeometry*,
0167 const MagneticField&);
0168
0169 void findHOEtaPhi(int iphsect, int& ietaho, int& iphiho);
0170
0171
0172
0173 float xhor0;
0174 float yhor0;
0175 float xhor1;
0176 float yhor1;
0177 int iring;
0178
0179 float localxhor0;
0180 float localyhor0;
0181 float localxhor1;
0182 float localyhor1;
0183
0184 TH2F* ho_occupency[5];
0185 bool m_occupancy;
0186 bool m_cosmic;
0187
0188 const int netabin = 16;
0189 const int nphimx = 72;
0190 const int netamx = 32;
0191 const int ncidmx = 5;
0192 const double rHOL0 = 382.0;
0193 const double rHOL1 = 407.0;
0194
0195 edm::InputTag muonTags_;
0196
0197 edm::EDGetTokenT<reco::TrackCollection> tok_muonsCosmic_;
0198 edm::EDGetTokenT<edm::View<reco::Muon> > tok_muons_;
0199 edm::EDGetTokenT<reco::VertexCollection> tok_vertex_;
0200
0201 edm::EDGetTokenT<LumiScalersCollection> tok_lumi_;
0202
0203 edm::EDGetTokenT<HBHERecHitCollection> tok_hbhe_;
0204 edm::EDGetTokenT<HORecHitCollection> tok_ho_;
0205 edm::EDGetTokenT<CaloTowerCollection> tok_tower_;
0206
0207 edm::ESGetToken<HcalChannelQuality, HcalChannelQualityRcd> tok_hcalChStatus_;
0208 edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tok_geom_;
0209 edm::ESGetToken<HcalSeverityLevelComputer, HcalSeverityLevelComputerRcd> tok_hcalSevLvlComputer_;
0210 edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> tok_magField_;
0211
0212 bool m_hbinfo;
0213 int m_startTS;
0214 int m_endTS;
0215 double m_sigma;
0216
0217 typedef math::Error<5>::type CovarianceMatrix;
0218 int Noccu;
0219 int nRuns;
0220
0221
0222 FreeTrajectoryState getFreeTrajectoryState(const reco::Track& tk, const MagneticField* field, int itag, bool dir);
0223
0224 unsigned int Ntp;
0225 std::map<std::string, bool> fired;
0226
0227
0228 const HcalChannelQuality* theHcalChStatus;
0229 const HcalSeverityLevelComputer* theHcalSevLvlComputer;
0230 int Nevents;
0231 };
0232
0233
0234
0235
0236
0237
0238
0239
0240
0241
0242
0243
0244 AlCaHOCalibProducer::AlCaHOCalibProducer(const edm::ParameterSet& iConfig) {
0245 usesResource(TFileService::kSharedResource);
0246
0247
0248 m_hbinfo = iConfig.getUntrackedParameter<bool>("hbinfo", false);
0249 m_sigma = iConfig.getUntrackedParameter<double>("sigma", 0.05);
0250 m_occupancy = iConfig.getUntrackedParameter<bool>("plotOccupancy", false);
0251 m_cosmic = iConfig.getUntrackedParameter<bool>("CosmicData", false);
0252
0253
0254 muonTags_ = iConfig.getUntrackedParameter<edm::InputTag>("muons");
0255 tok_muonsCosmic_ = consumes<reco::TrackCollection>(muonTags_);
0256 tok_muons_ = consumes<edm::View<reco::Muon> >(muonTags_);
0257 tok_vertex_ = consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("vertexTags"));
0258
0259 tok_lumi_ = consumes<LumiScalersCollection>(iConfig.getParameter<edm::InputTag>("lumiTags"));
0260 tok_ho_ = consumes<HORecHitCollection>(iConfig.getParameter<edm::InputTag>("hoInput"));
0261 tok_hbhe_ = consumes<HBHERecHitCollection>(iConfig.getParameter<edm::InputTag>("hbheInput"));
0262 tok_tower_ = consumes<CaloTowerCollection>(iConfig.getParameter<edm::InputTag>("towerInput"));
0263
0264 tok_hcalChStatus_ = esConsumes<HcalChannelQuality, HcalChannelQualityRcd>(edm::ESInputTag("", "withTopo"));
0265 tok_geom_ = esConsumes<CaloGeometry, CaloGeometryRecord>();
0266 tok_hcalSevLvlComputer_ = esConsumes<HcalSeverityLevelComputer, HcalSeverityLevelComputerRcd>();
0267 tok_magField_ = esConsumes<MagneticField, IdealMagneticFieldRecord>();
0268
0269 produces<HOCalibVariableCollection>("HOCalibVariableCollection").setBranchAlias("HOCalibVariableCollection");
0270
0271 if (m_occupancy) {
0272 edm::Service<TFileService> fs;
0273
0274 char title[200];
0275
0276 for (int ij = 0; ij < 5; ij++) {
0277 sprintf(title, "ho_occupency (>%i #sigma)", ij + 2);
0278 ho_occupency[ij] =
0279 fs->make<TH2F>(title, title, netamx + 1, -netamx - 0.5, netamx / 2 + 0.5, nphimx, 0.5, nphimx + 0.5);
0280 }
0281 }
0282 }
0283
0284
0285
0286
0287
0288 void AlCaHOCalibProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0289 edm::ParameterSetDescription desc;
0290 desc.add<edm::InputTag>("hbheInput", edm::InputTag("hbhereco"));
0291 desc.addUntracked<bool>("hotime", false);
0292 desc.addUntracked<bool>("hbinfo", false);
0293 desc.addUntracked<double>("sigma", 1.0);
0294 desc.addUntracked<bool>("plotOccupancy", false);
0295 desc.addUntracked<bool>("CosmicData", false);
0296 desc.add<edm::InputTag>("hoInput", edm::InputTag("horeco"));
0297 desc.add<edm::InputTag>("towerInput", edm::InputTag("towerMaker"));
0298 desc.addUntracked<std::string>("RootFileName", "test.root");
0299 desc.addUntracked<double>("m_scale", 4.0);
0300 desc.addUntracked<bool>("debug", false);
0301 desc.addUntracked<edm::InputTag>("muons", edm::InputTag("muons"));
0302 desc.add<edm::InputTag>("vertexTags", edm::InputTag("offlinePrimaryVertices"));
0303 desc.add<edm::InputTag>("lumiTags", edm::InputTag("scalersRawToDigi"));
0304 descriptions.add("alcaHOCalibProducer", desc);
0305 }
0306
0307
0308 void AlCaHOCalibProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0309 int irun = iEvent.id().run();
0310
0311
0312 Nevents++;
0313
0314 if (Nevents % 5000 == 1)
0315 edm::LogInfo("HOCalib") << "AlCaHOCalibProducer Processing event # " << Nevents << " " << Noccu << " " << irun
0316 << " " << iEvent.id().event();
0317
0318 theHcalChStatus = &iSetup.getData(tok_hcalChStatus_);
0319
0320 auto hostore = std::make_unique<HOCalibVariableCollection>();
0321
0322 edm::Handle<reco::TrackCollection> cosmicmuon;
0323 edm::Handle<edm::View<reco::Muon> > collisionmuon;
0324
0325 bool muonOK(true);
0326 HOCalibVariables tmpHOCalib;
0327 tmpHOCalib.nprim = -1;
0328 tmpHOCalib.inslumi = -1.;
0329
0330 if (m_cosmic) {
0331 cosmicmuon = iEvent.getHandle(tok_muonsCosmic_);
0332 muonOK = (cosmicmuon.isValid() && !cosmicmuon->empty());
0333 } else {
0334 collisionmuon = iEvent.getHandle(tok_muons_);
0335 muonOK = (collisionmuon.isValid() && !collisionmuon->empty());
0336
0337 if (iEvent.isRealData()) {
0338 auto const& primaryVertices = iEvent.getHandle(tok_vertex_);
0339 if (primaryVertices.isValid()) {
0340 tmpHOCalib.nprim = primaryVertices->size();
0341 }
0342
0343 tmpHOCalib.inslumi = 0.;
0344
0345 auto const& lumiScale = iEvent.getHandle(tok_lumi_);
0346
0347 if (lumiScale.isValid()) {
0348 if (lumiScale->empty()) {
0349 edm::LogError("HOCalib") << "lumiScale collection is empty";
0350 } else {
0351 tmpHOCalib.inslumi = lumiScale->begin()->pileup();
0352 }
0353 }
0354 }
0355 }
0356
0357 if (muonOK) {
0358 int Noccu_old = Noccu;
0359 edm::View<reco::Muon>::const_iterator muon1;
0360
0361 theHcalSevLvlComputer = &iSetup.getData(tok_hcalSevLvlComputer_);
0362
0363 MagneticField const& magField = iSetup.getData(tok_magField_);
0364
0365 const CaloGeometry& geo = iSetup.getData(tok_geom_);
0366 const CaloSubdetectorGeometry* gHO = geo.getSubdetectorGeometry(DetId::Hcal, HcalOuter);
0367
0368 if (m_cosmic) {
0369 int indx(0);
0370 for (reco::TrackCollection::const_iterator ncosm = cosmicmuon->begin(); ncosm != cosmicmuon->end();
0371 ++ncosm, ++indx) {
0372 if ((*ncosm).ndof() < 15)
0373 continue;
0374 if ((*ncosm).normalizedChi2() > 30.0)
0375 continue;
0376 reco::TrackRef tRef = reco::TrackRef(cosmicmuon, indx);
0377 fillHOStore(tRef, tmpHOCalib, hostore, Noccu_old, indx, cosmicmuon, muon1, iEvent, gHO, magField);
0378 }
0379 } else {
0380 for (muon1 = collisionmuon->begin(); muon1 < collisionmuon->end(); muon1++) {
0381 if ((!muon1->isGlobalMuon()) || (!muon1->isTrackerMuon()))
0382 continue;
0383 reco::TrackRef ncosm = muon1->innerTrack();
0384 fillHOStore(ncosm, tmpHOCalib, hostore, Noccu_old, 0, cosmicmuon, muon1, iEvent, gHO, magField);
0385 }
0386 }
0387 }
0388
0389 iEvent.put(std::move(hostore), "HOCalibVariableCollection");
0390 }
0391
0392
0393 void AlCaHOCalibProducer::beginJob() {
0394 Nevents = 0;
0395 nRuns = 0;
0396 Noccu = 0;
0397 }
0398
0399
0400 void AlCaHOCalibProducer::endJob() {
0401 if (m_occupancy) {
0402 for (int ij = 0; ij < 5; ij++) {
0403 ho_occupency[ij]->Scale(1. / std::max(1, Noccu));
0404 }
0405 }
0406 edm::LogInfo("HOCalib") << " AlCaHOCalibProducer processed event " << Nevents;
0407 }
0408
0409 void AlCaHOCalibProducer::fillHOStore(const reco::TrackRef& ncosm,
0410 HOCalibVariables& tmpHOCalib,
0411 std::unique_ptr<HOCalibVariableCollection>& hostore,
0412 int Noccu_old,
0413 int indx,
0414 edm::Handle<reco::TrackCollection> cosmicmuon,
0415 edm::View<reco::Muon>::const_iterator muon1,
0416 const edm::Event& iEvent,
0417 const CaloSubdetectorGeometry* gHO,
0418 const MagneticField& magField) {
0419
0420
0421 int charge = ncosm->charge();
0422
0423 double innerr = (*ncosm).innerPosition().Perp2();
0424 double outerr = (*ncosm).outerPosition().Perp2();
0425 int iiner = (innerr < outerr) ? 1 : 0;
0426
0427
0428
0429
0430
0431
0432
0433
0434
0435
0436 double posx, posy, posz;
0437 double momx, momy, momz;
0438
0439 if (iiner == 1) {
0440 posx = (*ncosm).innerPosition().X();
0441 posy = (*ncosm).innerPosition().Y();
0442 posz = (*ncosm).innerPosition().Z();
0443
0444 momx = (*ncosm).innerMomentum().X();
0445 momy = (*ncosm).innerMomentum().Y();
0446 momz = (*ncosm).innerMomentum().Z();
0447
0448 } else {
0449 posx = (*ncosm).outerPosition().X();
0450 posy = (*ncosm).outerPosition().Y();
0451 posz = (*ncosm).outerPosition().Z();
0452
0453 momx = (*ncosm).outerMomentum().X();
0454 momy = (*ncosm).outerMomentum().Y();
0455 momz = (*ncosm).outerMomentum().Z();
0456 }
0457
0458 PositionType trkpos(posx, posy, posz);
0459
0460 CLHEP::Hep3Vector tmpmuon3v(posx, posy, posz);
0461 CLHEP::Hep3Vector tmpmuondir(momx, momy, momz);
0462
0463 bool samedir = (tmpmuon3v.dot(tmpmuondir) > 0) ? true : false;
0464 for (int ij = 0; ij < 3; ij++) {
0465 tmpHOCalib.caloen[ij] = 0.0;
0466 }
0467 int inearbymuon = 0;
0468 localxhor0 = localyhor0 = 20000;
0469
0470 if (m_cosmic) {
0471 int ind(0);
0472 for (reco::TrackCollection::const_iterator ncosmcor = cosmicmuon->begin(); ncosmcor != cosmicmuon->end();
0473 ++ncosmcor, ++ind) {
0474 if (indx == ind)
0475 continue;
0476 CLHEP::Hep3Vector tmpmuon3vcor;
0477 CLHEP::Hep3Vector tmpmom3v;
0478 if (iiner == 1) {
0479 tmpmuon3vcor = CLHEP::Hep3Vector(
0480 (*ncosmcor).innerPosition().X(), (*ncosmcor).innerPosition().Y(), (*ncosmcor).innerPosition().Z());
0481 tmpmom3v = CLHEP::Hep3Vector(
0482 (*ncosmcor).innerMomentum().X(), (*ncosmcor).innerMomentum().Y(), (*ncosmcor).innerMomentum().Z());
0483 } else {
0484 tmpmuon3vcor = CLHEP::Hep3Vector(
0485 (*ncosmcor).outerPosition().X(), (*ncosmcor).outerPosition().Y(), (*ncosmcor).outerPosition().Z());
0486 tmpmom3v = CLHEP::Hep3Vector(
0487 (*ncosmcor).outerMomentum().X(), (*ncosmcor).outerMomentum().Y(), (*ncosmcor).outerMomentum().Z());
0488 }
0489
0490 if (tmpmom3v.mag() < 0.2 || (*ncosmcor).ndof() < 5)
0491 continue;
0492
0493 double angle = tmpmuon3v.angle(tmpmuon3vcor);
0494 if (angle < 7.5 * CLHEP::deg) {
0495 inearbymuon = 1;
0496 }
0497
0498
0499 if (angle < 7.5 * CLHEP::deg) {
0500 tmpHOCalib.caloen[0] += 1.;
0501 }
0502 if (angle < 15.0 * CLHEP::deg) {
0503 tmpHOCalib.caloen[1] += 1.;
0504 }
0505 if (angle < 35.0 * CLHEP::deg) {
0506 tmpHOCalib.caloen[2] += 1.;
0507 }
0508 }
0509 } else {
0510
0511 auto const& calotower = iEvent.getHandle(tok_tower_);
0512
0513 for (CaloTowerCollection::const_iterator calt = calotower->begin(); calt != calotower->end(); calt++) {
0514
0515 double ith = (*calt).momentum().theta();
0516 double iph = (*calt).momentum().phi();
0517
0518 CLHEP::Hep3Vector calo3v(sin(ith) * cos(iph), sin(ith) * sin(iph), cos(ith));
0519
0520 double angle = tmpmuon3v.angle(calo3v);
0521
0522 if (angle < 7.5 * CLHEP::deg) {
0523 tmpHOCalib.caloen[0] += calt->emEnergy() + calt->hadEnergy();
0524 }
0525 if (angle < 15 * CLHEP::deg) {
0526 tmpHOCalib.caloen[1] += calt->emEnergy() + calt->hadEnergy();
0527 }
0528 if (angle < 35 * CLHEP::deg) {
0529 tmpHOCalib.caloen[2] += calt->emEnergy() + calt->hadEnergy();
0530 }
0531 }
0532 }
0533 if ((m_cosmic) || (tmpHOCalib.caloen[0] <= 10.0)) {
0534 GlobalPoint glbpt(posx, posy, posz);
0535
0536 double mom = sqrt(momx * momx + momy * momy + momz * momz);
0537
0538 momx /= mom;
0539 momy /= mom;
0540 momz /= mom;
0541
0542 DirectionType trkdir(momx, momy, momz);
0543
0544 tmpHOCalib.trkdr = (*ncosm).d0();
0545 tmpHOCalib.trkdz = (*ncosm).dz();
0546 tmpHOCalib.nmuon = (m_cosmic) ? cosmicmuon->size() : 1;
0547 tmpHOCalib.trkvx = glbpt.x();
0548 tmpHOCalib.trkvy = glbpt.y();
0549 tmpHOCalib.trkvz = glbpt.z();
0550 tmpHOCalib.trkmm = mom * charge;
0551 tmpHOCalib.trkth = trkdir.theta();
0552 tmpHOCalib.trkph = trkdir.phi();
0553 tmpHOCalib.isect2 = -2;
0554 tmpHOCalib.isect = -2;
0555 tmpHOCalib.hodx = -100;
0556 tmpHOCalib.hody = -100;
0557 tmpHOCalib.hoang = -2.0;
0558 tmpHOCalib.momatho = -2;
0559 tmpHOCalib.ndof = (inearbymuon == 0) ? (int)(*ncosm).ndof() : -(int)(*ncosm).ndof();
0560 tmpHOCalib.chisq = (*ncosm).normalizedChi2();
0561 if (!m_cosmic) {
0562 reco::MuonEnergy muonenr = muon1->calEnergy();
0563 reco::MuonIsolation iso03 = muon1->isolationR03();
0564 reco::MuonIsolation iso05 = muon1->isolationR05();
0565
0566 tmpHOCalib.tkpt03 = iso03.sumPt;
0567 tmpHOCalib.ecal03 = iso05.sumPt;
0568 tmpHOCalib.hcal03 = iso03.hadEt + muonenr.had;
0569 }
0570 tmpHOCalib.therr = 0.;
0571 tmpHOCalib.pherr = 0.;
0572 if (iiner == 1) {
0573 reco::TrackBase::CovarianceMatrix innercov = (*ncosm).innerStateCovariance();
0574 tmpHOCalib.therr = innercov(1, 1);
0575 tmpHOCalib.pherr = innercov(2, 2);
0576 } else {
0577 reco::TrackBase::CovarianceMatrix outercov = (*ncosm).outerStateCovariance();
0578 tmpHOCalib.therr = outercov(1, 1);
0579 tmpHOCalib.pherr = outercov(2, 2);
0580 }
0581
0582 SteppingHelixPropagator myHelix(&magField, anyDirection);
0583 myHelix.setMaterialMode(false);
0584 myHelix.applyRadX0Correction(true);
0585 double phiho = trkpos.phi();
0586 if (phiho < 0)
0587 phiho += CLHEP::twopi;
0588
0589 int iphisect_dt = int(6 * (phiho + 10.0 * CLHEP::deg) / CLHEP::pi);
0590 if (iphisect_dt >= 12)
0591 iphisect_dt = 0;
0592
0593 int iphisect = -1;
0594 bool ipath = false;
0595 for (int kl = 0; kl <= 2; kl++) {
0596 int iphisecttmp = (kl < 2) ? iphisect_dt + kl : iphisect_dt - 1;
0597 if (iphisecttmp < 0)
0598 iphisecttmp = 11;
0599 if (iphisecttmp >= 12)
0600 iphisecttmp = 0;
0601
0602 double phipos = iphisecttmp * CLHEP::pi / 6.;
0603 double phirot = phipos;
0604
0605 GlobalVector xLocal(-sin(phirot), cos(phirot), 0.);
0606 GlobalVector yLocal(0., 0., 1.);
0607 GlobalVector zLocal = xLocal.cross(yLocal).unit();
0608
0609
0610 FreeTrajectoryState freetrajectorystate_ = getFreeTrajectoryState(*ncosm, &(magField), iiner, samedir);
0611
0612 Surface::RotationType rot(xLocal, yLocal, zLocal);
0613
0614 for (int ik = 1; ik >= 0; ik--) {
0615
0616 double radial = rHOL1;
0617 if (ik == 0)
0618 radial = rHOL0;
0619
0620 Surface::PositionType pos(radial * cos(phipos), radial * sin(phipos), 0.);
0621 PlaneBuilder::ReturnType aPlane = PlaneBuilder().plane(pos, rot);
0622
0623 auto aPlane2 = new Plane(pos, rot);
0624
0625 SteppingHelixStateInfo steppingHelixstateinfo_;
0626 myHelix.propagate(SteppingHelixStateInfo(freetrajectorystate_), (*aPlane2), steppingHelixstateinfo_);
0627
0628 if (steppingHelixstateinfo_.isValid()) {
0629 GlobalPoint hotrkpos2xx(steppingHelixstateinfo_.position().x(),
0630 steppingHelixstateinfo_.position().y(),
0631 steppingHelixstateinfo_.position().z());
0632
0633 if (ik == 1) {
0634 HcalDetId ClosestCell = (HcalDetId)gHO->getClosestCell(hotrkpos2xx);
0635 int ixeta = ClosestCell.ieta();
0636 int ixphi = ClosestCell.iphi();
0637 tmpHOCalib.isect2 = 100 * std::abs(ixeta + 50) + std::abs(ixphi);
0638 }
0639
0640 GlobalVector hotrkpos2(steppingHelixstateinfo_.position().x(),
0641 steppingHelixstateinfo_.position().y(),
0642 steppingHelixstateinfo_.position().z());
0643 CLHEP::Hep3Vector hotrkdir2(steppingHelixstateinfo_.momentum().x(),
0644 steppingHelixstateinfo_.momentum().y(),
0645 steppingHelixstateinfo_.momentum().z());
0646
0647 LocalVector lclvt0 = (*aPlane).toLocal(hotrkpos2);
0648
0649 double xx = lclvt0.x();
0650 double yy = lclvt0.y();
0651
0652 if (ik == 1) {
0653 if ((std::abs(yy) < 130 && xx > -64.7 && xx < 138.2)
0654 || (std::abs(yy) > 130 && std::abs(yy) < 700 && xx > -76.3 && xx < 140.5)) {
0655 ipath = true;
0656 iphisect = iphisecttmp;
0657 }
0658 }
0659
0660 if (iphisect != iphisecttmp)
0661 continue;
0662
0663 switch (ik) {
0664 case 0:
0665 xhor0 = xx;
0666 yhor0 = yy;
0667 break;
0668 case 1:
0669 xhor1 = xx;
0670 yhor1 = yy;
0671 tmpHOCalib.momatho = hotrkdir2.mag();
0672 tmpHOCalib.hoang = CLHEP::Hep3Vector(zLocal.x(), zLocal.y(), zLocal.z()).dot(hotrkdir2.unit());
0673 break;
0674 default:
0675 break;
0676 }
0677 } else {
0678 break;
0679 }
0680 }
0681 if (ipath)
0682 break;
0683 }
0684 if (ipath) {
0685
0686 int ietaho = 50;
0687 int iphiho = -1;
0688
0689 for (int ij = 0; ij < 9; ij++) {
0690 tmpHOCalib.hosig[ij] = -100.0;
0691 }
0692 for (int ij = 0; ij < 18; ij++) {
0693 tmpHOCalib.hocorsig[ij] = -100.0;
0694 }
0695 for (int ij = 0; ij < 9; ij++) {
0696 tmpHOCalib.hbhesig[ij] = -100.0;
0697 }
0698 tmpHOCalib.hocro = -100;
0699 tmpHOCalib.htime = -1000;
0700
0701 int isect = 0;
0702
0703 findHOEtaPhi(iphisect, ietaho, iphiho);
0704
0705 if (ietaho != 0 && iphiho != 0 && std::abs(iring) <= 2) {
0706 isect = 100 * std::abs(ietaho + 50) + std::abs(iphiho);
0707 if (std::abs(ietaho) >= netabin || iphiho < 0)
0708 isect *= -1;
0709 if (std::abs(ietaho) >= netabin)
0710 isect -= 1000000;
0711 if (iphiho < 0)
0712 isect -= 2000000;
0713 tmpHOCalib.isect = isect;
0714
0715 tmpHOCalib.hodx = localxhor1;
0716 tmpHOCalib.hody = localyhor1;
0717
0718 if (iring == 0) {
0719 tmpHOCalib.hocorsig[8] = localxhor0;
0720 tmpHOCalib.hocorsig[9] = localyhor0;
0721 }
0722
0723 int etamn = -4;
0724 int etamx = 4;
0725 if (iring == 1) {
0726 etamn = 5;
0727 etamx = 10;
0728 }
0729 if (iring == 2) {
0730 etamn = 11;
0731 etamx = 16;
0732 }
0733 if (iring == -1) {
0734 etamn = -10;
0735 etamx = -5;
0736 }
0737 if (iring == -2) {
0738 etamn = -16;
0739 etamx = -11;
0740 }
0741
0742 int phimn = 1;
0743 int phimx = 2;
0744 if (iring == 0) {
0745 phimx = 2 * int((iphiho + 1) / 2.);
0746 phimn = phimx - 1;
0747 } else {
0748 phimn = 3 * int((iphiho + 1) / 3.) - 1;
0749 phimx = phimn + 2;
0750 }
0751
0752 if (phimn < 1)
0753 phimn += nphimx;
0754 if (phimx > 72)
0755 phimx -= nphimx;
0756
0757 if (m_hbinfo) {
0758 for (int ij = 0; ij < 9; ij++) {
0759 tmpHOCalib.hbhesig[ij] = -100.0;
0760 }
0761
0762 auto const& hbheht = iEvent.getHandle(tok_hbhe_);
0763 if (!(*hbheht).empty()) {
0764 if ((*hbheht).empty())
0765 throw(int)(*hbheht).size();
0766
0767 for (HBHERecHitCollection::const_iterator jk = (*hbheht).begin(); jk != (*hbheht).end(); jk++) {
0768 HcalDetId id = (*jk).id();
0769 int tmpeta = id.ieta();
0770 int tmpphi = id.iphi();
0771
0772 int deta = tmpeta - ietaho;
0773 if (tmpeta < 0 && ietaho > 0)
0774 deta += 1;
0775 if (tmpeta > 0 && ietaho < 0)
0776 deta -= 1;
0777
0778
0779
0780
0781 int dphi = tmpphi - iphiho;
0782 if (dphi > nphimx / 2) {
0783 dphi -= nphimx;
0784 }
0785 if (dphi < -nphimx / 2) {
0786 dphi += nphimx;
0787 }
0788
0789
0790
0791
0792
0793
0794 if (m_occupancy) {
0795 float signal = (*jk).energy();
0796
0797 if (signal > -100 && Noccu == Noccu_old) {
0798 for (int ij = 0; ij < 5; ij++) {
0799 if (signal > (ij + 2) * m_sigma) {
0800 ho_occupency[ij]->Fill(tmpeta, tmpphi);
0801 }
0802 }
0803 }
0804 }
0805
0806 int ipass2 = (std::abs(deta) <= 1 && std::abs(dphi) <= 1) ? 1 : 0;
0807 if (ipass2 == 0)
0808 continue;
0809
0810 float signal = (*jk).energy();
0811
0812 if (3 * (deta + 1) + dphi + 1 < 9)
0813 tmpHOCalib.hbhesig[3 * (deta + 1) + dphi + 1] = signal;
0814 }
0815 }
0816 }
0817
0818 auto const& hoht = iEvent.getHandle(tok_ho_);
0819
0820 if (!(*hoht).empty()) {
0821 for (HORecHitCollection::const_iterator jk = (*hoht).begin(); jk != (*hoht).end(); jk++) {
0822 HcalDetId id = (*jk).id();
0823 int tmpeta = id.ieta();
0824 int tmpphi = id.iphi();
0825
0826 int ipass1 = 0;
0827 if (tmpeta >= etamn && tmpeta <= etamx) {
0828 if (phimn < phimx) {
0829 ipass1 = (tmpphi >= phimn && tmpphi <= phimx) ? 1 : 0;
0830 } else {
0831 ipass1 = (tmpphi == 71 || tmpphi == 72 || tmpphi == 1) ? 1 : 0;
0832 }
0833 }
0834
0835 int deta = tmpeta - ietaho;
0836 int dphi = tmpphi - iphiho;
0837
0838 if (tmpeta < 0 && ietaho > 0)
0839 deta += 1;
0840 if (tmpeta > 0 && ietaho < 0)
0841 deta -= 1;
0842
0843
0844
0845 if (dphi > nphimx / 2) {
0846 dphi -= nphimx;
0847 }
0848 if (dphi < -nphimx / 2) {
0849 dphi += nphimx;
0850 }
0851
0852
0853
0854
0855
0856 float signal = (*jk).energy();
0857
0858 int ipass2 = (std::abs(deta) <= 1 && std::abs(dphi) <= 1) ? 1 : 0;
0859
0860 if (ipass1 == 0 && ipass2 == 0)
0861 continue;
0862
0863 if (ipass1 == 1) {
0864 int tmpdph = tmpphi - phimn;
0865 if (tmpdph < 0)
0866 tmpdph = 2;
0867
0868 int ilog = 2 * (tmpeta - etamn) + tmpdph;
0869 if (iring != 0) {
0870 if (iring > 0) {
0871 ilog = 3 * (tmpeta - etamn) + tmpdph;
0872 } else {
0873 ilog = 3 * (etamx - tmpeta) + tmpdph;
0874 }
0875 }
0876 if (ilog > -1 && ilog < 18) {
0877 tmpHOCalib.hocorsig[ilog] = signal;
0878 }
0879 }
0880
0881 if (ipass2 == 1) {
0882 if (3 * (deta + 1) + dphi + 1 < 9) {
0883 tmpHOCalib.hosig[3 * (deta + 1) + dphi + 1] = signal;
0884 }
0885 }
0886
0887 if (deta == 0 && dphi == 0) {
0888 tmpHOCalib.htime = (*jk).time();
0889 tmpHOCalib.hoflag = (*jk).flags();
0890
0891
0892 unsigned theStatusValue = theHcalChStatus->getValues(id)->getValue();
0893
0894 int hitSeverity = theHcalSevLvlComputer->getSeverityLevel(id, (*jk).flags(), theStatusValue);
0895 tmpHOCalib.hoflag = hitSeverity;
0896 int crphi = tmpphi + 6;
0897 if (crphi > 72)
0898 crphi -= 72;
0899
0900 for (HORecHitCollection::const_iterator jcr = (*hoht).begin(); jcr != (*hoht).end(); jcr++) {
0901 const HORecHit reccr = (const HORecHit)(*jcr);
0902 HcalDetId idcr = reccr.id();
0903 int etacr = idcr.ieta();
0904 int phicr = idcr.iphi();
0905 if (tmpeta == etacr && crphi == phicr) {
0906 tmpHOCalib.hocro = reccr.energy();
0907 }
0908 }
0909 }
0910 }
0911 }
0912 }
0913
0914
0915 if (Noccu == Noccu_old)
0916 Noccu++;
0917 hostore->push_back(tmpHOCalib);
0918 }
0919 }
0920 }
0921
0922 void AlCaHOCalibProducer::findHOEtaPhi(int iphisect, int& ietaho, int& iphiho) {
0923
0924
0925 const double etalow[16] = {0.025,
0926 35.195,
0927 70.625,
0928 106.595,
0929 141.565,
0930 180.765,
0931 220.235,
0932 261.385,
0933 304.525,
0934 349.975,
0935 410.025,
0936 452.085,
0937 506.645,
0938 565.025,
0939 627.725,
0940 660.25};
0941 const double etahgh[16] = {35.145,
0942 70.575,
0943 106.545,
0944 125.505,
0945 180.715,
0946 220.185,
0947 261.335,
0948 304.475,
0949 349.925,
0950 392.575,
0951 452.035,
0952 506.595,
0953 564.975,
0954 627.675,
0955 661.075,
0956 700.25};
0957
0958 const double philow[6] = {-76.27, -35.11, 0.35, 35.81, 71.77, 108.93};
0959 const double phihgh[6] = {-35.81, -0.35, 35.11, 71.07, 108.23, 140.49};
0960
0961 const double philow00[6] = {-60.27, -32.91, 0.35, 33.61, 67.37, 102.23};
0962 const double phihgh00[6] = {-33.61, -0.35, 32.91, 66.67, 101.53, 129.49};
0963
0964 const double philow01[6] = {-64.67, -34.91, 0.35, 35.61, 71.37, 108.33};
0965 const double phihgh01[6] = {-35.61, -0.35, 34.91, 70.67, 107.63, 138.19};
0966
0967 iring = -10;
0968
0969 double tmpdy = std::abs(yhor1);
0970 for (int ij = 0; ij < netabin; ij++) {
0971 if (tmpdy > etalow[ij] && tmpdy < etahgh[ij]) {
0972 ietaho = ij + 1;
0973 float tmp1 = fabs(tmpdy - etalow[ij]);
0974 float tmp2 = fabs(tmpdy - etahgh[ij]);
0975
0976 localyhor1 = (tmp1 < tmp2) ? -tmp1 : tmp2;
0977 if (yhor1 < 0)
0978 localyhor1 *= -1.;
0979
0980 if (ij < 4)
0981 iring = 0;
0982 if (ij >= 4 && ij < 10)
0983 iring = 1;
0984 if (ij >= 10 && ij < netabin)
0985 iring = 2;
0986 break;
0987 }
0988 }
0989
0990 int tmpphi = 0;
0991 int tmpphi0 = 0;
0992
0993 if (ietaho > 4) {
0994 for (int ij = 0; ij < 6; ij++) {
0995 if (xhor1 > philow[ij] && xhor1 < phihgh[ij]) {
0996 tmpphi = ij + 1;
0997 float tmp1 = fabs(xhor1 - philow[ij]);
0998 float tmp2 = fabs(xhor1 - phihgh[ij]);
0999 localxhor1 = (tmp1 < tmp2) ? -tmp1 : tmp2;
1000 break;
1001 }
1002 }
1003 } else {
1004 for (int ij = 0; ij < 6; ij++) {
1005 if (xhor1 > philow01[ij] && xhor1 < phihgh01[ij]) {
1006 tmpphi = ij + 1;
1007 float tmp1 = fabs(xhor1 - philow01[ij]);
1008 float tmp2 = fabs(xhor1 - phihgh01[ij]);
1009 localxhor1 = (tmp1 < tmp2) ? -tmp1 : tmp2;
1010 break;
1011 }
1012 }
1013
1014 for (int ij = 0; ij < 6; ij++) {
1015 if (xhor0 > philow00[ij] && xhor0 < phihgh00[ij]) {
1016 tmpphi0 = ij + 1;
1017 float tmp1 = fabs(xhor0 - philow00[ij]);
1018 float tmp2 = fabs(xhor0 - phihgh00[ij]);
1019 localxhor0 = (tmp1 < tmp2) ? -tmp1 : tmp2;
1020 if (tmpphi != tmpphi0)
1021 localxhor0 += 10000.;
1022 break;
1023 }
1024 }
1025
1026 double tmpdy = std::abs(yhor0);
1027 for (int ij = 0; ij < 4; ij++) {
1028 if (tmpdy > etalow[ij] && tmpdy < etahgh[ij]) {
1029 float tmp1 = fabs(tmpdy - etalow[ij]);
1030 float tmp2 = fabs(tmpdy - etahgh[ij]);
1031 localyhor0 = (tmp1 < tmp2) ? -tmp1 : tmp2;
1032 if (yhor0 < 0)
1033 localyhor0 *= -1.;
1034 if (ij + 1 != ietaho)
1035 localyhor0 += 10000.;
1036 break;
1037 }
1038 }
1039 }
1040
1041 if (tmpphi != 0) {
1042 iphiho = 6 * iphisect - 2 + tmpphi;
1043 if (iphiho <= 0)
1044 iphiho += nphimx;
1045 if (iphiho > nphimx)
1046 iphiho -= nphimx;
1047 }
1048
1049
1050
1051 if (yhor1 < 0) {
1052 if (std::abs(ietaho) > netabin) {
1053 ietaho += 1;
1054 } else {
1055 ietaho *= -1;
1056 }
1057
1058 iring *= -1;
1059 }
1060 }
1061
1062 FreeTrajectoryState AlCaHOCalibProducer::getFreeTrajectoryState(const reco::Track& tk,
1063 const MagneticField* field,
1064 int iiner,
1065 bool dir) {
1066 if (iiner == 0) {
1067 GlobalPoint gpos(tk.outerX(), tk.outerY(), tk.outerZ());
1068 GlobalVector gmom(tk.outerPx(), tk.outerPy(), tk.outerPz());
1069 if (dir)
1070 gmom *= -1.;
1071 GlobalTrajectoryParameters par(gpos, gmom, tk.charge(), field);
1072 CurvilinearTrajectoryError err(tk.extra()->outerStateCovariance());
1073 return FreeTrajectoryState(par, err);
1074 } else {
1075 GlobalPoint gpos(tk.innerPosition().X(), tk.innerPosition().Y(), tk.innerPosition().Z());
1076 GlobalVector gmom(tk.innerMomentum().X(), tk.innerMomentum().Y(), tk.innerMomentum().Z());
1077 if (dir)
1078 gmom *= -1.;
1079 GlobalTrajectoryParameters par(gpos, -gmom, tk.charge(), field);
1080 CurvilinearTrajectoryError err(tk.extra()->innerStateCovariance());
1081 return FreeTrajectoryState(par, err);
1082 }
1083 }
1084
1085 #include "FWCore/Framework/interface/MakerMacros.h"
1086
1087
1088 DEFINE_FWK_MODULE(AlCaHOCalibProducer);