Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-09-07 04:34:59 

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

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!