Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-10-04 02:55:06 

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.inslumi=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), inslumi (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/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 // Necessary includes for identify severity of flagged problems in HO rechits
0123 //#include "RecoLocalCalo/HcalRecAlgos/interface/HcalCaloFlagLabels.h"
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 /* C++ Headers */
0135 #include <string>
0136 
0137 #include <iostream>
0138 #include <fstream>
0139 //
0140 // class decleration
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   //  virtual void endRun(edm::Run const &, edm::EventSetup const &) override;
0171   // ----------member data ---------------------------
0172 
0173   float xhor0;  //x-position in ring 0
0174   float yhor0;  //y-position in ring 0
0175   float xhor1;  //x-position in ring 1
0176   float yhor1;  //y-position in ring 1
0177   int iring;    //Ring number -2,-1,0,1,2
0178 
0179   float localxhor0;  //local x-distance from edege in ring 0
0180   float localyhor0;  //local y-distance from edege in ring 0
0181   float localxhor1;  //local x-distance from edege in ring 1
0182   float localyhor1;  //local y-distance from edege in ring 1
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_;  // cosmicMuons (for cosmic run) or muons (for collision run)
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   //  edm::EDGetTokenT<LumiDetails> tok_lumi_;
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   //  SteppingHelixPropagator* stepProp;
0222   FreeTrajectoryState getFreeTrajectoryState(const reco::Track& tk, const MagneticField* field, int itag, bool dir);
0223 
0224   unsigned int Ntp;  // # of HLT trigger paths (should be the same for all events!)
0225   std::map<std::string, bool> fired;
0226 
0227   //hcal severity ES
0228   const HcalChannelQuality* theHcalChStatus;
0229   const HcalSeverityLevelComputer* theHcalSevLvlComputer;
0230   int Nevents;
0231 };
0232 
0233 //
0234 // constants, enums and typedefs
0235 //
0236 
0237 //
0238 // static data member definitions
0239 //
0240 
0241 //
0242 // constructors and destructor
0243 //
0244 AlCaHOCalibProducer::AlCaHOCalibProducer(const edm::ParameterSet& iConfig) {
0245   usesResource(TFileService::kSharedResource);
0246   //register your products
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   // keep InputTag muonTags_ since it is used below. - cowden
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   //  tok_lumi_ = consumes<LumiDetails ,edm::InLumi>(iConfig.getParameter<edm::InputTag>("lumiTags"));
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 // member functions
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 // ------------ method called to produce the data  ------------
0308 void AlCaHOCalibProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0309   int irun = iEvent.id().run();
0310   //  int ilumi = iEvent.luminosityBlock();
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 // ------------ method called once each job just before starting event loop  ------------
0393 void AlCaHOCalibProducer::beginJob() {
0394   Nevents = 0;
0395   nRuns = 0;
0396   Noccu = 0;
0397 }
0398 
0399 // ------------ method called once each job just after ending the event loop  ------------
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   // Get Hcal Severity Level Computer, so that the severity of each rechit flag/status may be determined
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   //             in_to_out  Dir         in_to_out  Dir
0429   //   StandAlone ^         ^     Cosmic    ^    |
0430   //              |         |               |    v
0431   //---------------------------------------------------Y=0
0432   //   StandAlone |         |     Cosmic    ^    |
0433   //              v         v               |    v
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;  //GM for 22OCT07 data
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       }  //  break;}
0497 
0498       //    if (muonTagsi_.label() =="cosmicMuons") {
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     //            if (muonTags_.label() =="muons") {
0511     auto const& calotower = iEvent.getHandle(tok_tower_);
0512 
0513     for (CaloTowerCollection::const_iterator calt = calotower->begin(); calt != calotower->end(); calt++) {
0514       //CMSSW_2_1_x const math::XYZVector towermom = (*calt).momentum();
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();  // max(1.,tmpHOCalib.ndof);
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;  // iso03.emEt+muonenr.em;
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);  //thetaError();
0575       tmpHOCalib.pherr = innercov(2, 2);  //phi0Error();
0576     } else {
0577       reco::TrackBase::CovarianceMatrix outercov = (*ncosm).outerStateCovariance();
0578       tmpHOCalib.therr = outercov(1, 1);  //thetaError();
0579       tmpHOCalib.pherr = outercov(2, 2);  //phi0Error();
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);  //for u 18/12/06
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       //    GlobalVector zLocal(cos(phirot), sin(phirot), 0.0);
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--) {  //propagate track in two HO layers
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)                              //Ring-0
0654                 || (std::abs(yy) > 130 && std::abs(yy) < 700 && xx > -76.3 && xx < 140.5)) {  //Ring +-1,2
0655               ipath = true;  //Only look for tracks which as hits in layer 1
0656               iphisect = iphisecttmp;
0657             }
0658           }
0659 
0660           if (iphisect != iphisecttmp)
0661             continue;  //Look for ring-0 only when ring1 is accepted for that sector
0662 
0663           switch (ik) {
0664             case 0:
0665               xhor0 = xx;  //lclvt0.x();
0666               yhor0 = yy;  //lclvt0.y();
0667               break;
0668             case 1:
0669               xhor1 = xx;  //lclvt0.x();
0670               yhor1 = yy;  //lclvt0.y();
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) {  //If muon crossed HO laeyrs
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) {  //Muon passed through a tower
0706         isect = 100 * std::abs(ietaho + 50) + std::abs(iphiho);
0707         if (std::abs(ietaho) >= netabin || iphiho < 0)
0708           isect *= -1;  //Not extrapolated to any tower
0709         if (std::abs(ietaho) >= netabin)
0710           isect -= 1000000;  //not matched with eta
0711         if (iphiho < 0)
0712           isect -= 2000000;  //not matched with phi
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_);  // iEvent.getByType(hbheht);
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               //        if (tmpeta==-1 && ietaho== 1) deta = -1;
0779               //        if (tmpeta== 1 && ietaho==-1) deta =  1;
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               //        if (phimn >phimx) {
0790               //          if (dphi==71) dphi=-1;
0791               //          if (dphi==-71) dphi=1;
0792               //        }
0793 
0794               if (m_occupancy) {
0795                 float signal = (*jk).energy();
0796                 //      int tmpeta1 = (tmpeta>0) ? tmpeta -1 : -tmpeta +14;
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;  //NEED correction in full CMS detector
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         }  //m_hbinfo #endif
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             //        if (tmpeta==-1 && ietaho== 1) deta = -1;
0843             //        if (tmpeta== 1 && ietaho==-1) deta =  1;
0844 
0845             if (dphi > nphimx / 2) {
0846               dphi -= nphimx;
0847             }
0848             if (dphi < -nphimx / 2) {
0849               dphi += nphimx;
0850             }
0851             //        if (phimn>phimx) {
0852             //      if (dphi==71) dphi=-1;
0853             //      if (dphi==-71) dphi=1;
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;  //only case of iphi==1, where phimn=71
0867 
0868               int ilog = 2 * (tmpeta - etamn) + tmpdph;
0869               if (iring != 0) {
0870                 if (iring > 0) {
0871                   ilog = 3 * (tmpeta - etamn) + tmpdph;  //Again CMS correction
0872                 } else {
0873                   ilog = 3 * (etamx - tmpeta) + tmpdph;  //Again CMS correction
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;  //Again CMS azimuthal near phi 1&72
0884               }
0885             }
0886 
0887             if (deta == 0 && dphi == 0) {
0888               tmpHOCalib.htime = (*jk).time();
0889               tmpHOCalib.hoflag = (*jk).flags();
0890 
0891               // Get Channel Quality information for the given detID
0892               unsigned theStatusValue = theHcalChStatus->getValues(id)->getValue();
0893               // Now get severity of problems for the given detID, based on the rechit flag word and the channel quality status value
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       //GMA   Npass++;
0915       if (Noccu == Noccu_old)
0916         Noccu++;
0917       hostore->push_back(tmpHOCalib);
0918     }  // if (ipath)
0919   }    // Cut on calo energy
0920 }
0921 
0922 void AlCaHOCalibProducer::findHOEtaPhi(int iphisect, int& ietaho, int& iphiho) {
0923   //18/12/06 : use only position, not angle phi
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};  //Ring+/-1 & 2
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};  //Ring0 L0
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};  //Ring0 L1
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) {  //Ring 1 and 2
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 {  //Ring 0
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   //  isect2 = 15*iring+iphisect+1;
1050 
1051   if (yhor1 < 0) {
1052     if (std::abs(ietaho) > netabin) {  //Initialised with 50
1053       ietaho += 1;
1054     } else {
1055       ietaho *= -1;
1056     }
1057     //    isect2 *=-1;
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 //define this as a plug-in
1088 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!