Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 11:59:06 

0001 #include "FWCore/Framework/interface/Frameworkfwd.h"
0002 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0003 #include "FWCore/Framework/interface/EventSetup.h"
0004 #include "FWCore/Framework/interface/MakerMacros.h"
0005 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0006 #include "FWCore/Utilities/interface/InputTag.h"
0007 #include "FWCore/ServiceRegistry/interface/Service.h"
0008 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0009 
0010 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0011 #include "TrackingTools/TrackAssociator/interface/TrackDetectorAssociator.h"
0012 
0013 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0014 #include "DataFormats/HcalRecHit/interface/HcalSourcePositionData.h"
0015 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0016 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0017 
0018 #include "DataFormats/HcalDetId/interface/HcalElectronicsId.h"
0019 #include "DataFormats/HcalDigi/interface/HcalDigiCollections.h"
0020 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0021 #include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
0022 
0023 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0024 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0025 #include "CondFormats/HcalObjects/interface/HcalRespCorrs.h"
0026 #include "CondFormats/DataRecord/interface/HcalRespCorrsRcd.h"
0027 #include "CondFormats/HcalObjects/interface/HcalPFCorrs.h"
0028 #include "CondFormats/DataRecord/interface/HcalPFCorrsRcd.h"
0029 #include "MagneticField/Engine/interface/MagneticField.h"
0030 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0031 #include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h"
0032 #include "DataFormats/GeometrySurface/interface/PlaneBuilder.h"
0033 #include "DataFormats/GeometrySurface/interface/Cylinder.h"
0034 #include "DataFormats/GeometrySurface/interface/Plane.h"
0035 
0036 #include "Calibration/HcalCalibAlgos/interface/CommonUsefulStuff.h"
0037 
0038 #include <iostream>
0039 #include "TProfile.h"
0040 #include "TTree.h"
0041 
0042 //#define EDM_ML_DEBUG
0043 
0044 class HcalCorrPFCalculation : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0045 public:
0046   HcalCorrPFCalculation(edm::ParameterSet const& conf);
0047   void analyze(edm::Event const& ev, edm::EventSetup const& c) override;
0048   void beginJob() override;
0049 
0050 private:
0051   double RecalibFactor(HcalDetId id);
0052 
0053   const bool Respcorr_;
0054   const bool PFcorr_;
0055   const bool Conecorr_;
0056   // double radius_;
0057 
0058   const double clusterConeSize_, associationConeSize_, ecalCone_, neutralIsolationCone_, trackIsolationCone_;
0059   float eECAL, eECAL09cm, eECAL40cm;
0060   // double energyECALmip;
0061 
0062   float eCentHit, eTrack, eParticle;
0063   float CentHitFactor;
0064   int iEta, iPhi;
0065   //int iEtaTr, iPhiTr;
0066   float iDr, delR;
0067   float e3x3, e5x5;
0068   float phiTrack, etaTrack;
0069   float phiGPoint, etaGPoint;
0070 
0071   Bool_t doHF;
0072   Bool_t AddRecalib;
0073   int nevtot;
0074 
0075   const HcalRespCorrs* respRecalib;
0076   const HcalPFCorrs* pfRecalib;
0077 
0078   SteppingHelixPropagator* stepPropF;
0079   MagneticField* theMagField;
0080 
0081   TrackDetectorAssociator trackAssociator_;
0082   TrackAssociatorParameters parameters_;
0083 
0084   const CaloGeometry* geo;
0085   const HcalGeometry* gHcal;
0086 
0087   Float_t xTrkHcal, yTrkHcal, zTrkHcal;
0088   Float_t xTrkEcal, yTrkEcal, zTrkEcal;
0089   Float_t xAtHcal, yAtHcal, zAtHcal;
0090 
0091   float eEcalCone, eHcalCone, eHcalConeNoise;
0092   int UsedCells, UsedCellsNoise;
0093   float phiParticle, etaParticle;
0094 
0095   Int_t numValidTrkHits[50], numValidTrkStrips[50], numLayers[50];
0096 
0097   Bool_t trkQual[50];
0098 
0099   TProfile *nCells, *nCellsNoise, *enHcal, *enHcalNoise;
0100   //  TH1F *enEcalB, *enEcalE;
0101   TTree *pfTree, *tracksTree;
0102   //  TFile *rootFile;
0103   edm::Service<TFileService> fs;
0104   Int_t nTracks;
0105   Float_t genEta, genPhi, trackEta[50], trackPhi[50], trackP[50], delRmc[50];
0106 
0107   const edm::EDGetTokenT<HBHERecHitCollection> tok_hbhe_;
0108   const edm::EDGetTokenT<HORecHitCollection> tok_ho_;
0109   const edm::EDGetTokenT<HFRecHitCollection> tok_hf_;
0110 
0111   const edm::EDGetTokenT<EcalRecHitCollection> tok_EE_;
0112   const edm::EDGetTokenT<EcalRecHitCollection> tok_EB_;
0113   const edm::EDGetTokenT<reco::TrackCollection> tok_tracks_;
0114   const edm::EDGetTokenT<edm::HepMCProduct> tok_gen_;
0115 
0116   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> tok_bFieldH_;
0117   const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tok_geom_;
0118   const edm::ESGetToken<HcalRespCorrs, HcalRespCorrsRcd> tok_resp_;
0119   const edm::ESGetToken<HcalPFCorrs, HcalPFCorrsRcd> tok_pfcorr_;
0120 };
0121 
0122 HcalCorrPFCalculation::HcalCorrPFCalculation(edm::ParameterSet const& iConfig)
0123     : Respcorr_(iConfig.getUntrackedParameter<bool>("RespcorrAdd", false)),
0124       PFcorr_(iConfig.getUntrackedParameter<bool>("PFcorrAdd", false)),
0125       Conecorr_(iConfig.getUntrackedParameter<bool>("ConeCorrAdd", true)),
0126       clusterConeSize_(iConfig.getParameter<double>("clusterConeSize")),
0127       associationConeSize_(iConfig.getParameter<double>("associationConeSize")),
0128       ecalCone_(iConfig.getParameter<double>("ecalCone")),
0129       neutralIsolationCone_(iConfig.getParameter<double>("neutralIsolationCone")),
0130       trackIsolationCone_(iConfig.getParameter<double>("trackIsolationCone")),
0131       tok_hbhe_(consumes<HBHERecHitCollection>(iConfig.getParameter<edm::InputTag>("hbheRecHitCollectionTag"))),
0132       tok_ho_(consumes<HORecHitCollection>(iConfig.getParameter<edm::InputTag>("hoRecHitCollectionTag"))),
0133       tok_hf_(consumes<HFRecHitCollection>(iConfig.getParameter<edm::InputTag>("hfRecHitCollectionTag"))),
0134       tok_EE_(consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit", "EcalRecHitsEE"))),
0135       tok_EB_(consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit", "EcalRecHitsEB"))),
0136       tok_tracks_(consumes<reco::TrackCollection>(edm::InputTag("generalTracks"))),
0137       tok_gen_(consumes<edm::HepMCProduct>(edm::InputTag("generatorSmeared"))),
0138       tok_bFieldH_(esConsumes<MagneticField, IdealMagneticFieldRecord>()),
0139       tok_geom_(esConsumes<CaloGeometry, CaloGeometryRecord>()),
0140       tok_resp_(esConsumes<HcalRespCorrs, HcalRespCorrsRcd>()),
0141       tok_pfcorr_(esConsumes<HcalPFCorrs, HcalPFCorrsRcd>()) {
0142   usesResource(TFileService::kSharedResource);
0143 
0144   // should maybe add these options to configuration - cowden
0145 
0146   //  outputFile_ = iConfig.getUntrackedParameter<std::string>("outputFile", "myfile.root");
0147 
0148   //radius_       = iConfig.getUntrackedParameter<double>("ConeRadiusCm", 40.);
0149   //energyECALmip = iConfig.getParameter<double>("energyECALmip");
0150 
0151   edm::ParameterSet parameters = iConfig.getParameter<edm::ParameterSet>("TrackAssociatorParameters");
0152   edm::ConsumesCollector iC = consumesCollector();
0153   parameters_.loadParameters(parameters, iC);
0154   trackAssociator_.useDefaultPropagator();
0155 
0156   // AxB_=iConfig.getParameter<std::string>("AxB");
0157 }
0158 
0159 double HcalCorrPFCalculation::RecalibFactor(HcalDetId id) {
0160   Float_t resprecal = 1.;
0161   Float_t pfrecal = 1.;
0162   if (AddRecalib) {
0163     if (Respcorr_)
0164       resprecal = respRecalib->getValues(id)->getValue();
0165     if (PFcorr_)
0166       pfrecal = pfRecalib->getValues(id)->getValue();
0167   }
0168   Float_t factor = resprecal * pfrecal;
0169   return factor;
0170 }
0171 
0172 void HcalCorrPFCalculation::analyze(edm::Event const& ev, edm::EventSetup const& c) {
0173   AddRecalib = kFALSE;
0174 
0175   try {
0176     respRecalib = &c.getData(tok_resp_);
0177     pfRecalib = &c.getData(tok_pfcorr_);
0178 
0179     AddRecalib = kTRUE;
0180 #ifdef EDM_ML_DEBUG
0181     edm::LogVerbatim("CalibConstants") << "   OK ";
0182 #endif
0183 
0184   } catch (const cms::Exception& e) {
0185     edm::LogWarning("CalibConstants") << "   Not Found!! ";
0186   }
0187 
0188   const HBHERecHitCollection& Hithbhe = ev.get(tok_hbhe_);
0189   const HFRecHitCollection& Hithf = ev.get(tok_hf_);
0190   //  const HORecHitCollection & Hitho = ev.get(tok_ho_);
0191 
0192   const edm::Handle<EERecHitCollection>& ecalEE = ev.getHandle(tok_EE_);
0193   const EERecHitCollection HitecalEE = *(ecalEE.product());
0194 
0195   const edm::Handle<EBRecHitCollection>& ecalEB = ev.getHandle(tok_EB_);
0196   const EBRecHitCollection HitecalEB = *(ecalEB.product());
0197 
0198   // temporary collection of EB+EE recHits
0199   auto tmpEcalRecHitCollection = std::make_unique<EcalRecHitCollection>();
0200   for (EcalRecHitCollection::const_iterator recHit = (*ecalEB).begin(); recHit != (*ecalEB).end(); ++recHit) {
0201     tmpEcalRecHitCollection->push_back(*recHit);
0202   }
0203   for (EcalRecHitCollection::const_iterator recHit = (*ecalEE).begin(); recHit != (*ecalEE).end(); ++recHit) {
0204     tmpEcalRecHitCollection->push_back(*recHit);
0205   }
0206   const EcalRecHitCollection Hitecal = *tmpEcalRecHitCollection;
0207 
0208   const edm::Handle<reco::TrackCollection>& generalTracks = ev.getHandle(tok_tracks_);
0209 
0210   geo = &c.getData(tok_geom_);
0211 
0212   gHcal = static_cast<const HcalGeometry*>(geo->getSubdetectorGeometry(DetId::Hcal, HcalBarrel));
0213 
0214   parameters_.useEcal = true;
0215   parameters_.useHcal = true;
0216   parameters_.useCalo = false;
0217   parameters_.useMuon = false;
0218   parameters_.dREcal = 0.5;
0219   parameters_.dRHcal = 0.6;
0220 
0221   //parameters_.dREcal = taECALCone_;
0222   //parameters_.dRHcal = taHCALCone_;
0223 
0224   const MagneticField* bField = &c.getData(tok_bFieldH_);
0225   stepPropF = new SteppingHelixPropagator(bField, alongMomentum);
0226   stepPropF->setMaterialMode(false);
0227   stepPropF->applyRadX0Correction(true);
0228 
0229   //  double eta_bin[42]={0.,.087,.174,.261,.348,.435,.522,.609,.696,.783,
0230   //.870,.957,1.044,1.131,1.218,1.305,1.392,1.479,1.566,1.653,1.740,1.830,1.930,2.043,2.172,
0231   //2.322,2.500,2.650,2.853,3.000,3.139,3.314,3.489,3.664,3.839,4.013,4.191,4.363,4.538,4.716,4.889,5.191};
0232 
0233   // MC info
0234   //double phi_MC = -999999.;  // phi of initial particle from HepMC
0235   //double eta_MC = -999999.;  // eta of initial particle from HepMC
0236   double mom_MC = 50.;  // P of initial particle from HepMC
0237   //bool MC = false;
0238 
0239   // MC information
0240 
0241   const edm::Handle<edm::HepMCProduct>& evtMC = ev.getHandle(tok_gen_);
0242   if (!evtMC.isValid()) {
0243     edm::LogVerbatim("HcalCalib") << "no HepMCProduct found";
0244   } else {
0245     //MC=true;
0246 #ifdef EDM_ML_DEBUG
0247     edm::LogVerbatim("HcalCalib") << "*** source HepMCProduct found";
0248 #endif
0249   }
0250 
0251   // MC particle with highest pt is taken as a direction reference
0252   double maxPt = -99999.;
0253 
0254   GlobalPoint initpos(0, 0, 0);
0255 
0256   HepMC::GenEvent* myGenEvent = new HepMC::GenEvent(*(evtMC->GetEvent()));
0257   for (HepMC::GenEvent::particle_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end(); ++p) {
0258     phiParticle = (*p)->momentum().phi();
0259     etaParticle = (*p)->momentum().eta();
0260     double pt = (*p)->momentum().perp();
0261     mom_MC = (*p)->momentum().rho();
0262     if (pt > maxPt) {
0263       maxPt = pt; /*phi_MC = phiParticle; eta_MC = etaParticle;*/
0264     }
0265     GlobalVector mom((*p)->momentum().x(), (*p)->momentum().y(), (*p)->momentum().z());
0266     int charge = -1;
0267 
0268     if (abs((*p)->pdg_id()) == 211)
0269       charge = (*p)->pdg_id() / abs((*p)->pdg_id());  // pions only !!!
0270     else
0271       continue;
0272 
0273     /*  Propogate the partoicle to Hcal */
0274 
0275     const FreeTrajectoryState* freetrajectorystate_ = new FreeTrajectoryState(initpos, mom, charge, &(*theMagField));
0276 
0277     TrackDetMatchInfo info = trackAssociator_.associate(ev, c, *freetrajectorystate_, parameters_);
0278 
0279     GlobalPoint barrelMC(0, 0, 0), endcapMC(0, 0, 0), forwardMC1(0, 0, 0), forwardMC2(0, 0, 0);
0280 
0281     GlobalPoint gPointHcal(0., 0., 0.);
0282 
0283     /*   
0284       xTrkHcal=info.trkGlobPosAtHcal.x();
0285       yTrkHcal=info.trkGlobPosAtHcal.y();
0286       zTrkHcal=info.trkGlobPosAtHcal.z();
0287       //GlobalPoint gPointHcal(xTrkHcal,yTrkHcal,zTrkHcal);
0288       
0289       GlobalVector trackMomAtHcal = info.trkMomAtHcal;
0290  
0291       if (xTrkHcal==0 && yTrkHcal==0 && zTrkHcal==0) continue;
0292       */
0293 
0294     if (fabs(etaParticle) < 1.392) {
0295       Cylinder* cylinder = new Cylinder(181.1, Surface::PositionType(0, 0, 0), Surface::RotationType());
0296 
0297       TrajectoryStateOnSurface steppingHelixstateinfo_ = stepPropF->propagate(*freetrajectorystate_, (*cylinder));
0298       if (steppingHelixstateinfo_.isValid()) {
0299         barrelMC = steppingHelixstateinfo_.freeState()->position();
0300       }
0301     }
0302 
0303     doHF = kFALSE;
0304     if (fabs(etaParticle) >= 1.392 && fabs(etaParticle) < 5.191) {
0305       Surface::RotationType rot(GlobalVector(1, 0, 0), GlobalVector(0, 1, 0));
0306 
0307       Surface::PositionType pos(0., 0., 400.5 * fabs(etaParticle) / etaParticle);
0308       PlaneBuilder::ReturnType aPlane = PlaneBuilder().plane(pos, rot);
0309 
0310       TrajectoryStateOnSurface steppingHelixstateinfo_ = stepPropF->propagate(*freetrajectorystate_, (*aPlane));
0311       if (steppingHelixstateinfo_.isValid() && fabs(steppingHelixstateinfo_.freeState()->position().eta()) < 3.0)
0312         endcapMC = steppingHelixstateinfo_.freeState()->position();
0313       if (steppingHelixstateinfo_.isValid() && fabs(steppingHelixstateinfo_.freeState()->position().eta()) > 3.0)
0314         doHF = kTRUE;
0315     }
0316     if (doHF) {
0317       if (abs(etaParticle) > 5.191)
0318         continue;
0319 
0320       if (abs(etaParticle) > 2.99) {
0321         Surface::RotationType rot(GlobalVector(1, 0, 0), GlobalVector(0, 1, 0));
0322 
0323         Surface::PositionType pos1(0., 0., 1125 * fabs(etaParticle) / etaParticle);
0324         //  Surface::PositionType pos1(0., 0.,1115*fabs(etaParticle)/etaParticle);
0325         Surface::PositionType pos2(0., 0., 1137 * fabs(etaParticle) / etaParticle);
0326         PlaneBuilder::ReturnType aPlane1 = PlaneBuilder().plane(pos1, rot);
0327         PlaneBuilder::ReturnType aPlane2 = PlaneBuilder().plane(pos2, rot);
0328 
0329         TrajectoryStateOnSurface steppingHelixstateinfo_ = stepPropF->propagate(*freetrajectorystate_, (*aPlane1));
0330         if (steppingHelixstateinfo_.isValid())
0331           forwardMC1 = steppingHelixstateinfo_.freeState()->position();
0332 
0333         steppingHelixstateinfo_ = stepPropF->propagate(*freetrajectorystate_, (*aPlane2));
0334         if (steppingHelixstateinfo_.isValid())
0335           forwardMC2 = steppingHelixstateinfo_.freeState()->position();
0336       }
0337     }
0338     /*   ------------       ------------    -----------------------------  */
0339 
0340     /*  Finding the closest cell at Hcal  */
0341     Int_t iphitrue = -10;
0342     Int_t ietatrue = 100;
0343     HcalDetId tempId;
0344 
0345     if (abs(etaParticle) < 1.392) {
0346       gPointHcal = barrelMC;
0347       tempId = gHcal->getClosestCell(gPointHcal);
0348     }
0349     if (abs(etaParticle) >= 1.392 && abs(etaParticle) < 3.0) {
0350       gPointHcal = endcapMC;
0351       tempId = gHcal->getClosestCell(gPointHcal);
0352     }
0353     if (abs(etaParticle) >= 3.0 && abs(etaParticle) < 5.191) {
0354       gPointHcal = forwardMC1;
0355       tempId = gHcal->getClosestCell(gPointHcal);
0356       //tempId = gHcal->CaloSubdetectorGeometry::getClosestCell(gPointHcal);
0357     }
0358 
0359     tempId = gHcal->getClosestCell(gPointHcal);
0360 
0361     ietatrue = tempId.ieta();
0362     iphitrue = tempId.iphi();
0363 
0364     etaGPoint = gPointHcal.eta();
0365     phiGPoint = gPointHcal.phi();
0366 
0367     //xAtHcal = gPointHcal.x();
0368     //yAtHcal = gPointHcal.y();
0369     //zAtHcal = gPointHcal.z();
0370     /*       -----------------   ------------------------      */
0371 
0372     if (gPointHcal.x() == 0 && gPointHcal.y() == 0 && gPointHcal.z() == 0) {
0373 #ifdef EDM_ML_DEBUG
0374       edm::LogVerbatim("HcalCalib") << "gPointHcal is Zero!";
0375 #endif
0376       continue;
0377     }
0378 
0379     float etahcal = gPointHcal.eta();
0380     if (abs(etahcal) > 5.192)
0381       continue;
0382 #ifdef EDM_ML_DEBUG
0383     if (std::abs(etahcal) > 3.0 && std::abs(etahcal) < 5.191) {
0384       edm::LogVerbatim("HcalCalib") << gPointHcal.x() << "   " << gPointHcal.y() << "   " << gPointHcal.z() << "    "
0385                                     << gPointHcal.eta() << "  " << gPointHcal.phi() << "   " << ietatrue << "   "
0386                                     << iphitrue;
0387       if (ietatrue == 100 || iphitrue == -10) {
0388         edm::LogVerbatim("HcalCalib") << "ietatrue: " << ietatrue << "   iphitrue: " << iphitrue
0389                                       << "  etahcal: " << etahcal << "  phihcal: " << gPointHcal.phi();
0390       }
0391     }
0392 #endif
0393     /*   -------------   Calculate Ecal Energy using TrackAssociator  ---------------------- */
0394 
0395     //float etaecal=info.trkGlobPosAtEcal.eta();
0396 
0397     xTrkEcal = info.trkGlobPosAtEcal.x();
0398     yTrkEcal = info.trkGlobPosAtEcal.y();
0399     zTrkEcal = info.trkGlobPosAtEcal.z();
0400 
0401     GlobalPoint gPointEcal(xTrkEcal, yTrkEcal, zTrkEcal);
0402 
0403     eECAL = ecalEnergyInCone(gPointEcal, ecalCone_, Hitecal, geo);
0404     eECAL09cm = ecalEnergyInCone(gPointEcal, 9, Hitecal, geo);
0405     eECAL40cm = ecalEnergyInCone(gPointEcal, 40, Hitecal, geo);
0406 
0407     eEcalCone = eECAL;
0408     //if(abs(etaecal)<1.5) enEcalB -> Fill(eEcalCone);
0409     //if(abs(etaecal)>1.5 && abs(etaecal)<3.1) enEcalE -> Fill(eEcalCone);
0410 
0411     /*    ------------------------------              --------------------------       ----------------- */
0412 
0413     /*    ----------------- Find the Hottest Hcal RecHit   --------------------------- */
0414     MaxHit_struct MaxHit;
0415 
0416     MaxHit.hitenergy = -100.;
0417     Float_t recal = 1.0;
0418 
0419     //Hcal:
0420     eHcalCone = 0.;
0421     eHcalConeNoise = 0.;
0422     UsedCells = 0;
0423     UsedCellsNoise = 0;
0424     e3x3 = 0.;
0425     e5x5 = 0.;
0426 
0427     for (HBHERecHitCollection::const_iterator hhit = Hithbhe.begin(); hhit != Hithbhe.end(); hhit++)
0428     //for (HcalRecHitCollection::const_iterator hhit=Hithcal.begin(); hhit!=Hithcal.end(); hhit++)
0429     {
0430       recal = RecalibFactor(hhit->detid());
0431       GlobalPoint pos = gHcal->getPosition(hhit->detid());
0432 
0433       int iphihit = (hhit->id()).iphi();
0434       int ietahit = (hhit->id()).ieta();
0435       int depthhit = (hhit->id()).depth();
0436       float enehit = hhit->energy() * recal;
0437 
0438       if (depthhit != 1)
0439         continue;
0440 
0441       double distAtHcal = getDistInPlaneSimple(gPointHcal, pos);
0442 
0443       if (distAtHcal < clusterConeSize_) {
0444         for (HBHERecHitCollection::const_iterator hhit2 = Hithbhe.begin(); hhit2 != Hithbhe.end(); hhit2++)
0445         //for (HcalRecHitCollection::const_iterator hhit2=Hithcal.begin(); hhit2!=Hithcal.end(); hhit2++)
0446         {
0447           int iphihit2 = (hhit2->id()).iphi();
0448           int ietahit2 = (hhit2->id()).ieta();
0449           int depthhit2 = (hhit2->id()).depth();
0450           float enehit2 = hhit2->energy() * recal;
0451 
0452           if (iphihit == iphihit2 && ietahit == ietahit2 && depthhit != depthhit2)
0453             enehit = enehit + enehit2;
0454         }
0455 
0456         if (enehit > MaxHit.hitenergy) {
0457           MaxHit.hitenergy = enehit;
0458           MaxHit.ietahitm = (hhit->id()).ieta();
0459           MaxHit.iphihitm = (hhit->id()).iphi();
0460           MaxHit.dr = distAtHcal;
0461           //MaxHit.depthhit  = (hhit->id()).depth();
0462           MaxHit.depthhit = 1;
0463         }
0464       }
0465     }
0466 
0467     for (HFRecHitCollection::const_iterator hhit = Hithf.begin(); hhit != Hithf.end(); hhit++) {
0468       recal = RecalibFactor(hhit->detid());
0469 
0470       GlobalPoint pos = gHcal->getPosition(hhit->detid());
0471 
0472       int iphihit = (hhit->id()).iphi();
0473       int ietahit = (hhit->id()).ieta();
0474       int depthhit = (hhit->id()).depth();
0475       float enehit = hhit->energy() * recal;
0476 
0477       double distAtHcal = getDistInPlaneSimple(gPointHcal, pos);
0478 
0479       if (distAtHcal < associationConeSize_) {
0480         for (HFRecHitCollection::const_iterator hhit2 = Hithf.begin(); hhit2 != Hithf.end(); hhit2++) {
0481           int iphihit2 = (hhit2->id()).iphi();
0482           int ietahit2 = (hhit2->id()).ieta();
0483           int depthhit2 = (hhit2->id()).depth();
0484           float enehit2 = hhit2->energy() * recal;
0485 
0486           if (iphihit == iphihit2 && ietahit == ietahit2 && depthhit != depthhit2)
0487             enehit = enehit + enehit2;
0488         }
0489 
0490         if (enehit > MaxHit.hitenergy) {
0491           MaxHit.hitenergy = enehit;
0492           MaxHit.ietahitm = (hhit->id()).ieta();
0493           MaxHit.iphihitm = (hhit->id()).iphi();
0494           MaxHit.dr = distAtHcal;
0495           MaxHit.depthhit = 1;
0496         }
0497       }
0498     }
0499 
0500     /*    ----------------------       ----------------              --------------------------------------------   */
0501 
0502     /*    ----------- Collect Hcal Energy in a Cone (and also 3x3 and 5x5 around the hottest cell)*/
0503 
0504     for (HBHERecHitCollection::const_iterator hhit = Hithbhe.begin(); hhit != Hithbhe.end(); hhit++)
0505     //    for (HcalRecHitCollection::const_iterator hhit=Hithcal.begin(); hhit!=Hithcal.end(); hhit++)
0506     {
0507       recal = RecalibFactor(hhit->detid());
0508 #ifdef EDM_ML_DEBUG
0509       edm::LogVerbatim("HcalCalib") << "recal: " << recal;
0510 #endif
0511       GlobalPoint pos = gHcal->getPosition(hhit->detid());
0512 
0513       int iphihit = (hhit->id()).iphi();
0514       int ietahit = (hhit->id()).ieta();
0515       int depthhit = (hhit->id()).depth();
0516       float enehit = hhit->energy() * recal;
0517 
0518       //if (depthhit!=1) continue;
0519 
0520       //Set noise RecHit opposite to track hits
0521       int iphihitNoise = iphihit > 36 ? iphihit - 36 : iphihit + 36;
0522       int ietahitNoise = -ietahit;
0523       int depthhitNoise = depthhit;
0524 
0525       double distAtHcal = getDistInPlaneSimple(gPointHcal, pos);
0526       if (distAtHcal < clusterConeSize_ && MaxHit.hitenergy > 0.) {
0527         eHcalCone += enehit;
0528         UsedCells++;
0529 
0530         int DIETA = 100;
0531         if (MaxHit.ietahitm * (hhit->id()).ieta() > 0) {
0532           DIETA = MaxHit.ietahitm - (hhit->id()).ieta();
0533         }
0534         if (MaxHit.ietahitm * (hhit->id()).ieta() < 0) {
0535           DIETA = MaxHit.ietahitm - (hhit->id()).ieta();
0536           DIETA = DIETA > 0 ? DIETA - 1 : DIETA + 1;
0537         }
0538         int DIPHI = abs(MaxHit.iphihitm - (hhit->id()).iphi());
0539         DIPHI = DIPHI > 36 ? 72 - DIPHI : DIPHI;
0540 
0541         if (abs(DIETA) <= 2 && (abs(DIPHI) <= 2 || ((abs(MaxHit.ietahitm) > 20 && abs(DIPHI) <= 4) &&
0542                                                     !((abs(MaxHit.ietahitm) == 21 || abs(MaxHit.ietahitm) == 22) &&
0543                                                       abs((hhit->id()).ieta()) <= 20 && abs(DIPHI) > 2)))) {
0544           e5x5 += hhit->energy();
0545         }
0546         if (abs(DIETA) <= 1 &&
0547             (abs(DIPHI) <= 1 || ((abs(MaxHit.ietahitm) > 20 && abs(DIPHI) <= 2) &&
0548                                  !(abs(MaxHit.ietahitm) == 21 && abs((hhit->id()).ieta()) <= 20 && abs(DIPHI) > 1)))) {
0549           e3x3 += hhit->energy();
0550         }
0551 #ifdef EDM_ML_DEBUG
0552         edm::LogVerbatim("HcalCalib") << "track: ieta " << ietahit << " iphi: " << iphihit << " depth: " << depthhit
0553                                       << " energydepos: " << enehit;
0554 #endif
0555         for (HBHERecHitCollection::const_iterator hhit2 = Hithbhe.begin(); hhit2 != Hithbhe.end(); hhit2++) {
0556           recal = RecalibFactor(hhit2->detid());
0557           int iphihit2 = (hhit2->id()).iphi();
0558           int ietahit2 = (hhit2->id()).ieta();
0559           int depthhit2 = (hhit2->id()).depth();
0560           float enehit2 = hhit2->energy() * recal;
0561 
0562           if (iphihitNoise == iphihit2 && ietahitNoise == ietahit2 && depthhitNoise == depthhit2 && enehit2 > 0.) {
0563             eHcalConeNoise += hhit2->energy() * recal;
0564             UsedCellsNoise++;
0565 #ifdef EDM_ML_DEBUG
0566             edm::LogVerbatim("HcalCalib") << "Noise: ieta " << ietahit2 << " iphi: " << iphihit2
0567                                           << " depth: " << depthhit2 << " energydepos: " << enehit2;
0568 #endif
0569           }
0570         }
0571       }
0572     }  //end of all HBHE  hits loop
0573 
0574     for (HFRecHitCollection::const_iterator hhit = Hithf.begin(); hhit != Hithf.end(); hhit++) {
0575       recal = RecalibFactor(hhit->detid());
0576 
0577       GlobalPoint pos = gHcal->getPosition(hhit->detid());
0578       //float phihit = pos.phi();
0579       //float etahit = pos.eta();
0580 
0581       int iphihit = (hhit->id()).iphi();
0582       int ietahit = (hhit->id()).ieta();
0583       int depthhit = (hhit->id()).depth();
0584       float enehit = hhit->energy() * recal;
0585 
0586       //Set noise RecHit opposite to track hits
0587       int iphihitNoise = iphihit > 36 ? iphihit - 36 : iphihit + 36;
0588       int ietahitNoise = -ietahit;
0589       int depthhitNoise = depthhit;
0590 
0591       double distAtHcal = getDistInPlaneSimple(gPointHcal, pos);
0592 
0593       if (distAtHcal < clusterConeSize_ && MaxHit.hitenergy > 0.)
0594       //if(dr<radius_ && enehit>0.)
0595       {
0596         eHcalCone += enehit;
0597         UsedCells++;
0598 
0599         int DIETA = 100;
0600         if (MaxHit.ietahitm * (hhit->id()).ieta() > 0) {
0601           DIETA = MaxHit.ietahitm - (hhit->id()).ieta();
0602         }
0603         if (MaxHit.ietahitm * (hhit->id()).ieta() < 0) {
0604           DIETA = MaxHit.ietahitm - (hhit->id()).ieta();
0605           DIETA = DIETA > 0 ? DIETA - 1 : DIETA + 1;
0606         }
0607         int DIPHI = abs(MaxHit.iphihitm - (hhit->id()).iphi());
0608         DIPHI = DIPHI > 36 ? 72 - DIPHI : DIPHI;
0609 
0610         if (abs(DIETA) <= 2 && (abs(DIPHI) <= 2 || ((abs(MaxHit.ietahitm) > 20 && abs(DIPHI) <= 4) &&
0611                                                     !((abs(MaxHit.ietahitm) == 21 || abs(MaxHit.ietahitm) == 22) &&
0612                                                       abs((hhit->id()).ieta()) <= 20 && abs(DIPHI) > 2)))) {
0613           e5x5 += hhit->energy();
0614         }
0615         if (abs(DIETA) <= 1 &&
0616             (abs(DIPHI) <= 1 || ((abs(MaxHit.ietahitm) > 20 && abs(DIPHI) <= 2) &&
0617                                  !(abs(MaxHit.ietahitm) == 21 && abs((hhit->id()).ieta()) <= 20 && abs(DIPHI) > 1)))) {
0618           e3x3 += hhit->energy();
0619         }
0620 
0621         for (HFRecHitCollection::const_iterator hhit2 = Hithf.begin(); hhit2 != Hithf.end(); hhit2++) {
0622           recal = RecalibFactor(hhit2->detid());
0623 
0624           int iphihit2 = (hhit2->id()).iphi();
0625           int ietahit2 = (hhit2->id()).ieta();
0626           int depthhit2 = (hhit2->id()).depth();
0627           float enehit2 = hhit2->energy() * recal;
0628 
0629           if (iphihitNoise == iphihit2 && ietahitNoise == ietahit2 && depthhitNoise == depthhit2 && enehit2 > 0.01) {
0630             eHcalConeNoise += hhit2->energy() * recal;
0631             UsedCellsNoise++;
0632           }
0633         }
0634       }
0635     }  //end of all HF hits loop
0636 
0637     /*  ----------------      --------------------         ----------------------------------------------       -------- */
0638 
0639     /* ------------- -   Track-MC matching  (if any tracks are in event)    ------------    - */
0640 
0641     nTracks = 0;
0642 
0643     delRmc[0] = 5;
0644 
0645     float delR_track_particle = 100;
0646 
0647     for (reco::TrackCollection::const_iterator track1 = generalTracks->begin(); track1 != generalTracks->end();
0648          track1++) {
0649       delR_track_particle = deltaR(etaParticle, phiParticle, track1->eta(), track1->phi());
0650 
0651       trackEta[nTracks] = track1->eta();
0652       trackPhi[nTracks] = track1->phi();
0653       trackP[nTracks] = sqrt(track1->px() * track1->px() + track1->py() * track1->py() + track1->pz() * track1->pz());
0654 
0655       delRmc[nTracks] = delR_track_particle;
0656       numValidTrkHits[nTracks] = track1->hitPattern().numberOfValidHits();
0657       numValidTrkStrips[nTracks] = track1->hitPattern().numberOfValidStripTECHits();
0658       numLayers[nTracks] = track1->hitPattern().trackerLayersWithMeasurement();  //layers crossed
0659       trkQual[nTracks] = track1->quality(reco::TrackBase::highPurity);
0660 
0661       nTracks++;
0662     }
0663 
0664     /*        ------------------          ------------------------------ ------- */
0665 
0666     int dieta_M_P = 100;
0667     int diphi_M_P = 100;
0668     if (MaxHit.ietahitm * ietatrue > 0) {
0669       dieta_M_P = abs(MaxHit.ietahitm - ietatrue);
0670     }
0671     if (MaxHit.ietahitm * ietatrue < 0) {
0672       dieta_M_P = abs(MaxHit.ietahitm - ietatrue) - 1;
0673     }
0674     diphi_M_P = abs(MaxHit.iphihitm - iphitrue);
0675 
0676     diphi_M_P = diphi_M_P > 36 ? 72 - diphi_M_P : diphi_M_P;
0677     iDr = sqrt(diphi_M_P * diphi_M_P + dieta_M_P * dieta_M_P);
0678 
0679 #ifdef EDM_ML_DEBUG
0680     if (iDr > 15) {
0681       edm::LogVerbatim("HcalCalib") << "diphi: " << diphi_M_P << "  dieta: " << dieta_M_P << "   iDr: " << iDr
0682                                     << " ietatrue:" << ietatrue << "  iphitrue:" << iphitrue;
0683       edm::LogVerbatim("HcalCalib") << "M ieta: " << MaxHit.ietahitm << "  M iphi: " << MaxHit.iphihitm;
0684     }
0685 #endif
0686 
0687     Bool_t passCuts = kFALSE;
0688     passCuts = kTRUE;
0689     //if(eEcalCone < energyECALmip && iDr<2.) passCuts = kTRUE;
0690     //if(MaxHit.hitenergy>0.) passCuts = kTRUE;
0691 
0692     if (passCuts) {
0693       /*      
0694       enHcal -> Fill(ietatrue,  eHcalCone);
0695       nCells -> Fill(ietatrue,  UsedCells);
0696       enHcalNoise -> Fill(ietatrue,  eHcalConeNoise);
0697       nCellsNoise -> Fill(ietatrue,  UsedCellsNoise); 
0698       */
0699 
0700       //e3x3=0; e5x5=0;
0701 
0702       iEta = ietatrue;
0703       iPhi = iphitrue;
0704 
0705       //iEta = MaxHit.ietahitm;
0706       //iPhi = MaxHit.iphihitm;
0707       delR = MaxHit.dr;
0708       eCentHit = MaxHit.hitenergy;
0709 
0710       eParticle = mom_MC;
0711       //eTrack = mom_MC;
0712       //phiTrack = phiParticle;
0713       //etaTrack = etaParticle;
0714 
0715       pfTree->Fill();
0716     }
0717 
0718   }  //Hep:MC
0719 }
0720 
0721 void HcalCorrPFCalculation::beginJob() {
0722   pfTree = fs->make<TTree>("pfTree", "Tree for pf info");
0723 
0724   pfTree->Branch("nTracks", &nTracks, "nTracks/I");
0725   pfTree->Branch("trackEta", trackEta, "trackEta[nTracks]/F");
0726   pfTree->Branch("trackPhi", trackPhi, "trackPhi[nTracks]/F");
0727   pfTree->Branch("trackP", trackP, "trackP[nTracks]/F");
0728 
0729   pfTree->Branch("delRmc", delRmc, "delRmc[nTracks]/F");
0730   pfTree->Branch("numValidTrkHits", numValidTrkHits, "numValidTrkHits[nTracks]/I");
0731   pfTree->Branch("numValidTrkStrips", numValidTrkStrips, "numValidTrkStrips[nTracks]/I");
0732   pfTree->Branch("numLayers", numLayers, "numLayers[nTracks]/I");
0733   pfTree->Branch("trkQual", trkQual, "trkQual[nTracks]/O");
0734 
0735   pfTree->Branch("eEcalCone", &eEcalCone, "eEcalCone/F");
0736   pfTree->Branch("eHcalCone", &eHcalCone, "eHcalCone/F");
0737   pfTree->Branch("eHcalConeNoise", &eHcalConeNoise, "eHcalConeNoise/F");
0738 
0739   pfTree->Branch("UsedCellsNoise", &UsedCellsNoise, "UsedCellsNoise/I");
0740   pfTree->Branch("UsedCells", &UsedCells, "UsedCells/I");
0741 
0742   pfTree->Branch("eCentHit", &eCentHit, "eCentHit/F");
0743 
0744   pfTree->Branch("eParticle", &eParticle, "eParticle/F");
0745   pfTree->Branch("etaParticle", &etaParticle, "etaParticle/F");
0746   pfTree->Branch("phiParticle", &phiParticle, "phiParticle/F");
0747 
0748   pfTree->Branch("etaGPoint", &etaGPoint, "etaGPoint/F");
0749   pfTree->Branch("phiGPoint", &phiGPoint, "phiGPoint/F");
0750 
0751   pfTree->Branch("xAtHcal", &xAtHcal, "xAtHcal/F");
0752   pfTree->Branch("yAtHcal", &yAtHcal, "yAtHcal/F");
0753   pfTree->Branch("zAtHcal", &zAtHcal, "zAtHcal/F");
0754 
0755   pfTree->Branch("eECAL09cm", &eECAL09cm, "eECAL09cm/F");
0756   pfTree->Branch("eECAL40cm", &eECAL40cm, "eECAL40cm/F");
0757   pfTree->Branch("eECAL", &eECAL, "eECAL/F");
0758 
0759   pfTree->Branch("e3x3 ", &e3x3, "e3x3/F");
0760   pfTree->Branch("e5x5", &e5x5, "e5x5/F");
0761 
0762   pfTree->Branch("iDr", &iDr, "iDr/F");
0763   pfTree->Branch("delR", &delR, "delR/F");
0764 
0765   pfTree->Branch("iEta", &iEta, "iEta/I");
0766   pfTree->Branch("iPhi", &iPhi, "iPhi/I");
0767 
0768   //  pfTree->Branch("numValidTrkHits", &numValidTrkHits, "numValidTrkHits/I");
0769   // pfTree->Branch("numValidTrkStrips", &numValidTrkStrips, "numValidTrkStrips/I");
0770   // pfTree->Branch("trkQual", &trkQual, "trkQual/");
0771   // pfTree->Branch("numLayers", &numLayers, "numLayers/I");
0772 }
0773 
0774 #include "FWCore/Framework/interface/MakerMacros.h"
0775 
0776 //define this as a plug-in
0777 DEFINE_FWK_MODULE(HcalCorrPFCalculation);
[ 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!