Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-02-28 01:32: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   int npart = 0;
0254 
0255   GlobalPoint initpos(0, 0, 0);
0256 
0257   HepMC::GenEvent* myGenEvent = new HepMC::GenEvent(*(evtMC->GetEvent()));
0258   for (HepMC::GenEvent::particle_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end(); ++p) {
0259     phiParticle = (*p)->momentum().phi();
0260     etaParticle = (*p)->momentum().eta();
0261     double pt = (*p)->momentum().perp();
0262     mom_MC = (*p)->momentum().rho();
0263     if (pt > maxPt) {
0264       npart++;
0265       maxPt = pt; /*phi_MC = phiParticle; eta_MC = etaParticle;*/
0266     }
0267     GlobalVector mom((*p)->momentum().x(), (*p)->momentum().y(), (*p)->momentum().z());
0268     int charge = -1;
0269 
0270     if (abs((*p)->pdg_id()) == 211)
0271       charge = (*p)->pdg_id() / abs((*p)->pdg_id());  // pions only !!!
0272     else
0273       continue;
0274 
0275     /*  Propogate the partoicle to Hcal */
0276 
0277     const FreeTrajectoryState* freetrajectorystate_ = new FreeTrajectoryState(initpos, mom, charge, &(*theMagField));
0278 
0279     TrackDetMatchInfo info = trackAssociator_.associate(ev, c, *freetrajectorystate_, parameters_);
0280 
0281     GlobalPoint barrelMC(0, 0, 0), endcapMC(0, 0, 0), forwardMC1(0, 0, 0), forwardMC2(0, 0, 0);
0282 
0283     GlobalPoint gPointHcal(0., 0., 0.);
0284 
0285     /*   
0286       xTrkHcal=info.trkGlobPosAtHcal.x();
0287       yTrkHcal=info.trkGlobPosAtHcal.y();
0288       zTrkHcal=info.trkGlobPosAtHcal.z();
0289       //GlobalPoint gPointHcal(xTrkHcal,yTrkHcal,zTrkHcal);
0290       
0291       GlobalVector trackMomAtHcal = info.trkMomAtHcal;
0292  
0293       if (xTrkHcal==0 && yTrkHcal==0 && zTrkHcal==0) continue;
0294       */
0295 
0296     if (fabs(etaParticle) < 1.392) {
0297       Cylinder* cylinder = new Cylinder(181.1, Surface::PositionType(0, 0, 0), Surface::RotationType());
0298 
0299       TrajectoryStateOnSurface steppingHelixstateinfo_ = stepPropF->propagate(*freetrajectorystate_, (*cylinder));
0300       if (steppingHelixstateinfo_.isValid()) {
0301         barrelMC = steppingHelixstateinfo_.freeState()->position();
0302       }
0303     }
0304 
0305     doHF = kFALSE;
0306     if (fabs(etaParticle) >= 1.392 && fabs(etaParticle) < 5.191) {
0307       Surface::RotationType rot(GlobalVector(1, 0, 0), GlobalVector(0, 1, 0));
0308 
0309       Surface::PositionType pos(0., 0., 400.5 * fabs(etaParticle) / etaParticle);
0310       PlaneBuilder::ReturnType aPlane = PlaneBuilder().plane(pos, rot);
0311 
0312       TrajectoryStateOnSurface steppingHelixstateinfo_ = stepPropF->propagate(*freetrajectorystate_, (*aPlane));
0313       if (steppingHelixstateinfo_.isValid() && fabs(steppingHelixstateinfo_.freeState()->position().eta()) < 3.0)
0314         endcapMC = steppingHelixstateinfo_.freeState()->position();
0315       if (steppingHelixstateinfo_.isValid() && fabs(steppingHelixstateinfo_.freeState()->position().eta()) > 3.0)
0316         doHF = kTRUE;
0317     }
0318     if (doHF) {
0319       if (abs(etaParticle) > 5.191)
0320         continue;
0321 
0322       if (abs(etaParticle) > 2.99) {
0323         Surface::RotationType rot(GlobalVector(1, 0, 0), GlobalVector(0, 1, 0));
0324 
0325         Surface::PositionType pos1(0., 0., 1125 * fabs(etaParticle) / etaParticle);
0326         //  Surface::PositionType pos1(0., 0.,1115*fabs(etaParticle)/etaParticle);
0327         Surface::PositionType pos2(0., 0., 1137 * fabs(etaParticle) / etaParticle);
0328         PlaneBuilder::ReturnType aPlane1 = PlaneBuilder().plane(pos1, rot);
0329         PlaneBuilder::ReturnType aPlane2 = PlaneBuilder().plane(pos2, rot);
0330 
0331         TrajectoryStateOnSurface steppingHelixstateinfo_ = stepPropF->propagate(*freetrajectorystate_, (*aPlane1));
0332         if (steppingHelixstateinfo_.isValid())
0333           forwardMC1 = steppingHelixstateinfo_.freeState()->position();
0334 
0335         steppingHelixstateinfo_ = stepPropF->propagate(*freetrajectorystate_, (*aPlane2));
0336         if (steppingHelixstateinfo_.isValid())
0337           forwardMC2 = steppingHelixstateinfo_.freeState()->position();
0338       }
0339     }
0340     /*   ------------       ------------    -----------------------------  */
0341 
0342     /*  Finding the closest cell at Hcal  */
0343     Int_t iphitrue = -10;
0344     Int_t ietatrue = 100;
0345     HcalDetId tempId;
0346 
0347     if (abs(etaParticle) < 1.392) {
0348       gPointHcal = barrelMC;
0349       tempId = gHcal->getClosestCell(gPointHcal);
0350     }
0351     if (abs(etaParticle) >= 1.392 && abs(etaParticle) < 3.0) {
0352       gPointHcal = endcapMC;
0353       tempId = gHcal->getClosestCell(gPointHcal);
0354     }
0355     if (abs(etaParticle) >= 3.0 && abs(etaParticle) < 5.191) {
0356       gPointHcal = forwardMC1;
0357       tempId = gHcal->getClosestCell(gPointHcal);
0358       //tempId = gHcal->CaloSubdetectorGeometry::getClosestCell(gPointHcal);
0359     }
0360 
0361     tempId = gHcal->getClosestCell(gPointHcal);
0362 
0363     ietatrue = tempId.ieta();
0364     iphitrue = tempId.iphi();
0365 
0366     etaGPoint = gPointHcal.eta();
0367     phiGPoint = gPointHcal.phi();
0368 
0369     //xAtHcal = gPointHcal.x();
0370     //yAtHcal = gPointHcal.y();
0371     //zAtHcal = gPointHcal.z();
0372     /*       -----------------   ------------------------      */
0373 
0374     if (gPointHcal.x() == 0 && gPointHcal.y() == 0 && gPointHcal.z() == 0) {
0375 #ifdef EDM_ML_DEBUG
0376       edm::LogVerbatim("HcalCalib") << "gPointHcal is Zero!";
0377 #endif
0378       continue;
0379     }
0380 
0381     float etahcal = gPointHcal.eta();
0382     if (abs(etahcal) > 5.192)
0383       continue;
0384 #ifdef EDM_ML_DEBUG
0385     if (std::abs(etahcal) > 3.0 && std::abs(etahcal) < 5.191) {
0386       edm::LogVerbatim("HcalCalib") << gPointHcal.x() << "   " << gPointHcal.y() << "   " << gPointHcal.z() << "    "
0387                                     << gPointHcal.eta() << "  " << gPointHcal.phi() << "   " << ietatrue << "   "
0388                                     << iphitrue;
0389       if (ietatrue == 100 || iphitrue == -10) {
0390         edm::LogVerbatim("HcalCalib") << "ietatrue: " << ietatrue << "   iphitrue: " << iphitrue
0391                                       << "  etahcal: " << etahcal << "  phihcal: " << gPointHcal.phi();
0392       }
0393     }
0394 #endif
0395     /*   -------------   Calculate Ecal Energy using TrackAssociator  ---------------------- */
0396 
0397     //float etaecal=info.trkGlobPosAtEcal.eta();
0398 
0399     xTrkEcal = info.trkGlobPosAtEcal.x();
0400     yTrkEcal = info.trkGlobPosAtEcal.y();
0401     zTrkEcal = info.trkGlobPosAtEcal.z();
0402 
0403     GlobalPoint gPointEcal(xTrkEcal, yTrkEcal, zTrkEcal);
0404 
0405     eECAL = ecalEnergyInCone(gPointEcal, ecalCone_, Hitecal, geo);
0406     eECAL09cm = ecalEnergyInCone(gPointEcal, 9, Hitecal, geo);
0407     eECAL40cm = ecalEnergyInCone(gPointEcal, 40, Hitecal, geo);
0408 
0409     eEcalCone = eECAL;
0410     //if(abs(etaecal)<1.5) enEcalB -> Fill(eEcalCone);
0411     //if(abs(etaecal)>1.5 && abs(etaecal)<3.1) enEcalE -> Fill(eEcalCone);
0412 
0413     /*    ------------------------------              --------------------------       ----------------- */
0414 
0415     /*    ----------------- Find the Hottest Hcal RecHit   --------------------------- */
0416     MaxHit_struct MaxHit;
0417 
0418     MaxHit.hitenergy = -100.;
0419     Float_t recal = 1.0;
0420 
0421     //Hcal:
0422     eHcalCone = 0.;
0423     eHcalConeNoise = 0.;
0424     UsedCells = 0;
0425     UsedCellsNoise = 0;
0426     e3x3 = 0.;
0427     e5x5 = 0.;
0428 
0429     for (HBHERecHitCollection::const_iterator hhit = Hithbhe.begin(); hhit != Hithbhe.end(); hhit++)
0430     //for (HcalRecHitCollection::const_iterator hhit=Hithcal.begin(); hhit!=Hithcal.end(); hhit++)
0431     {
0432       recal = RecalibFactor(hhit->detid());
0433       GlobalPoint pos = gHcal->getPosition(hhit->detid());
0434 
0435       int iphihit = (hhit->id()).iphi();
0436       int ietahit = (hhit->id()).ieta();
0437       int depthhit = (hhit->id()).depth();
0438       float enehit = hhit->energy() * recal;
0439 
0440       if (depthhit != 1)
0441         continue;
0442 
0443       double distAtHcal = getDistInPlaneSimple(gPointHcal, pos);
0444 
0445       if (distAtHcal < clusterConeSize_) {
0446         for (HBHERecHitCollection::const_iterator hhit2 = Hithbhe.begin(); hhit2 != Hithbhe.end(); hhit2++)
0447         //for (HcalRecHitCollection::const_iterator hhit2=Hithcal.begin(); hhit2!=Hithcal.end(); hhit2++)
0448         {
0449           int iphihit2 = (hhit2->id()).iphi();
0450           int ietahit2 = (hhit2->id()).ieta();
0451           int depthhit2 = (hhit2->id()).depth();
0452           float enehit2 = hhit2->energy() * recal;
0453 
0454           if (iphihit == iphihit2 && ietahit == ietahit2 && depthhit != depthhit2)
0455             enehit = enehit + enehit2;
0456         }
0457 
0458         if (enehit > MaxHit.hitenergy) {
0459           MaxHit.hitenergy = enehit;
0460           MaxHit.ietahitm = (hhit->id()).ieta();
0461           MaxHit.iphihitm = (hhit->id()).iphi();
0462           MaxHit.dr = distAtHcal;
0463           //MaxHit.depthhit  = (hhit->id()).depth();
0464           MaxHit.depthhit = 1;
0465         }
0466       }
0467     }
0468 
0469     for (HFRecHitCollection::const_iterator hhit = Hithf.begin(); hhit != Hithf.end(); hhit++) {
0470       recal = RecalibFactor(hhit->detid());
0471 
0472       GlobalPoint pos = gHcal->getPosition(hhit->detid());
0473 
0474       int iphihit = (hhit->id()).iphi();
0475       int ietahit = (hhit->id()).ieta();
0476       int depthhit = (hhit->id()).depth();
0477       float enehit = hhit->energy() * recal;
0478 
0479       double distAtHcal = getDistInPlaneSimple(gPointHcal, pos);
0480 
0481       if (distAtHcal < associationConeSize_) {
0482         for (HFRecHitCollection::const_iterator hhit2 = Hithf.begin(); hhit2 != Hithf.end(); hhit2++) {
0483           int iphihit2 = (hhit2->id()).iphi();
0484           int ietahit2 = (hhit2->id()).ieta();
0485           int depthhit2 = (hhit2->id()).depth();
0486           float enehit2 = hhit2->energy() * recal;
0487 
0488           if (iphihit == iphihit2 && ietahit == ietahit2 && depthhit != depthhit2)
0489             enehit = enehit + enehit2;
0490         }
0491 
0492         if (enehit > MaxHit.hitenergy) {
0493           MaxHit.hitenergy = enehit;
0494           MaxHit.ietahitm = (hhit->id()).ieta();
0495           MaxHit.iphihitm = (hhit->id()).iphi();
0496           MaxHit.dr = distAtHcal;
0497           MaxHit.depthhit = 1;
0498         }
0499       }
0500     }
0501 
0502     /*    ----------------------       ----------------              --------------------------------------------   */
0503 
0504     /*    ----------- Collect Hcal Energy in a Cone (and also 3x3 and 5x5 around the hottest cell)*/
0505 
0506     for (HBHERecHitCollection::const_iterator hhit = Hithbhe.begin(); hhit != Hithbhe.end(); hhit++)
0507     //    for (HcalRecHitCollection::const_iterator hhit=Hithcal.begin(); hhit!=Hithcal.end(); hhit++)
0508     {
0509       recal = RecalibFactor(hhit->detid());
0510 #ifdef EDM_ML_DEBUG
0511       edm::LogVerbatim("HcalCalib") << "recal: " << recal;
0512 #endif
0513       GlobalPoint pos = gHcal->getPosition(hhit->detid());
0514 
0515       int iphihit = (hhit->id()).iphi();
0516       int ietahit = (hhit->id()).ieta();
0517       int depthhit = (hhit->id()).depth();
0518       float enehit = hhit->energy() * recal;
0519 
0520       //if (depthhit!=1) continue;
0521 
0522       //Set noise RecHit opposite to track hits
0523       int iphihitNoise = iphihit > 36 ? iphihit - 36 : iphihit + 36;
0524       int ietahitNoise = -ietahit;
0525       int depthhitNoise = depthhit;
0526 
0527       double distAtHcal = getDistInPlaneSimple(gPointHcal, pos);
0528       if (distAtHcal < clusterConeSize_ && MaxHit.hitenergy > 0.) {
0529         eHcalCone += enehit;
0530         UsedCells++;
0531 
0532         int DIETA = 100;
0533         if (MaxHit.ietahitm * (hhit->id()).ieta() > 0) {
0534           DIETA = MaxHit.ietahitm - (hhit->id()).ieta();
0535         }
0536         if (MaxHit.ietahitm * (hhit->id()).ieta() < 0) {
0537           DIETA = MaxHit.ietahitm - (hhit->id()).ieta();
0538           DIETA = DIETA > 0 ? DIETA - 1 : DIETA + 1;
0539         }
0540         int DIPHI = abs(MaxHit.iphihitm - (hhit->id()).iphi());
0541         DIPHI = DIPHI > 36 ? 72 - DIPHI : DIPHI;
0542 
0543         if (abs(DIETA) <= 2 && (abs(DIPHI) <= 2 || ((abs(MaxHit.ietahitm) > 20 && abs(DIPHI) <= 4) &&
0544                                                     !((abs(MaxHit.ietahitm) == 21 || abs(MaxHit.ietahitm) == 22) &&
0545                                                       abs((hhit->id()).ieta()) <= 20 && abs(DIPHI) > 2)))) {
0546           e5x5 += hhit->energy();
0547         }
0548         if (abs(DIETA) <= 1 &&
0549             (abs(DIPHI) <= 1 || ((abs(MaxHit.ietahitm) > 20 && abs(DIPHI) <= 2) &&
0550                                  !(abs(MaxHit.ietahitm) == 21 && abs((hhit->id()).ieta()) <= 20 && abs(DIPHI) > 1)))) {
0551           e3x3 += hhit->energy();
0552         }
0553 #ifdef EDM_ML_DEBUG
0554         edm::LogVerbatim("HcalCalib") << "track: ieta " << ietahit << " iphi: " << iphihit << " depth: " << depthhit
0555                                       << " energydepos: " << enehit;
0556 #endif
0557         for (HBHERecHitCollection::const_iterator hhit2 = Hithbhe.begin(); hhit2 != Hithbhe.end(); hhit2++) {
0558           recal = RecalibFactor(hhit2->detid());
0559           int iphihit2 = (hhit2->id()).iphi();
0560           int ietahit2 = (hhit2->id()).ieta();
0561           int depthhit2 = (hhit2->id()).depth();
0562           float enehit2 = hhit2->energy() * recal;
0563 
0564           if (iphihitNoise == iphihit2 && ietahitNoise == ietahit2 && depthhitNoise == depthhit2 && enehit2 > 0.) {
0565             eHcalConeNoise += hhit2->energy() * recal;
0566             UsedCellsNoise++;
0567 #ifdef EDM_ML_DEBUG
0568             edm::LogVerbatim("HcalCalib") << "Noise: ieta " << ietahit2 << " iphi: " << iphihit2
0569                                           << " depth: " << depthhit2 << " energydepos: " << enehit2;
0570 #endif
0571           }
0572         }
0573       }
0574     }  //end of all HBHE  hits loop
0575 
0576     for (HFRecHitCollection::const_iterator hhit = Hithf.begin(); hhit != Hithf.end(); hhit++) {
0577       recal = RecalibFactor(hhit->detid());
0578 
0579       GlobalPoint pos = gHcal->getPosition(hhit->detid());
0580       //float phihit = pos.phi();
0581       //float etahit = pos.eta();
0582 
0583       int iphihit = (hhit->id()).iphi();
0584       int ietahit = (hhit->id()).ieta();
0585       int depthhit = (hhit->id()).depth();
0586       float enehit = hhit->energy() * recal;
0587 
0588       //Set noise RecHit opposite to track hits
0589       int iphihitNoise = iphihit > 36 ? iphihit - 36 : iphihit + 36;
0590       int ietahitNoise = -ietahit;
0591       int depthhitNoise = depthhit;
0592 
0593       double distAtHcal = getDistInPlaneSimple(gPointHcal, pos);
0594 
0595       if (distAtHcal < clusterConeSize_ && MaxHit.hitenergy > 0.)
0596       //if(dr<radius_ && enehit>0.)
0597       {
0598         eHcalCone += enehit;
0599         UsedCells++;
0600 
0601         int DIETA = 100;
0602         if (MaxHit.ietahitm * (hhit->id()).ieta() > 0) {
0603           DIETA = MaxHit.ietahitm - (hhit->id()).ieta();
0604         }
0605         if (MaxHit.ietahitm * (hhit->id()).ieta() < 0) {
0606           DIETA = MaxHit.ietahitm - (hhit->id()).ieta();
0607           DIETA = DIETA > 0 ? DIETA - 1 : DIETA + 1;
0608         }
0609         int DIPHI = abs(MaxHit.iphihitm - (hhit->id()).iphi());
0610         DIPHI = DIPHI > 36 ? 72 - DIPHI : DIPHI;
0611 
0612         if (abs(DIETA) <= 2 && (abs(DIPHI) <= 2 || ((abs(MaxHit.ietahitm) > 20 && abs(DIPHI) <= 4) &&
0613                                                     !((abs(MaxHit.ietahitm) == 21 || abs(MaxHit.ietahitm) == 22) &&
0614                                                       abs((hhit->id()).ieta()) <= 20 && abs(DIPHI) > 2)))) {
0615           e5x5 += hhit->energy();
0616         }
0617         if (abs(DIETA) <= 1 &&
0618             (abs(DIPHI) <= 1 || ((abs(MaxHit.ietahitm) > 20 && abs(DIPHI) <= 2) &&
0619                                  !(abs(MaxHit.ietahitm) == 21 && abs((hhit->id()).ieta()) <= 20 && abs(DIPHI) > 1)))) {
0620           e3x3 += hhit->energy();
0621         }
0622 
0623         for (HFRecHitCollection::const_iterator hhit2 = Hithf.begin(); hhit2 != Hithf.end(); hhit2++) {
0624           recal = RecalibFactor(hhit2->detid());
0625 
0626           int iphihit2 = (hhit2->id()).iphi();
0627           int ietahit2 = (hhit2->id()).ieta();
0628           int depthhit2 = (hhit2->id()).depth();
0629           float enehit2 = hhit2->energy() * recal;
0630 
0631           if (iphihitNoise == iphihit2 && ietahitNoise == ietahit2 && depthhitNoise == depthhit2 && enehit2 > 0.01) {
0632             eHcalConeNoise += hhit2->energy() * recal;
0633             UsedCellsNoise++;
0634           }
0635         }
0636       }
0637     }  //end of all HF hits loop
0638 
0639     /*  ----------------      --------------------         ----------------------------------------------       -------- */
0640 
0641     /* ------------- -   Track-MC matching  (if any tracks are in event)    ------------    - */
0642 
0643     nTracks = 0;
0644 
0645     delRmc[0] = 5;
0646 
0647     float delR_track_particle = 100;
0648 
0649     for (reco::TrackCollection::const_iterator track1 = generalTracks->begin(); track1 != generalTracks->end();
0650          track1++) {
0651       delR_track_particle = deltaR(etaParticle, phiParticle, track1->eta(), track1->phi());
0652 
0653       trackEta[nTracks] = track1->eta();
0654       trackPhi[nTracks] = track1->phi();
0655       trackP[nTracks] = sqrt(track1->px() * track1->px() + track1->py() * track1->py() + track1->pz() * track1->pz());
0656 
0657       delRmc[nTracks] = delR_track_particle;
0658       numValidTrkHits[nTracks] = track1->hitPattern().numberOfValidHits();
0659       numValidTrkStrips[nTracks] = track1->hitPattern().numberOfValidStripTECHits();
0660       numLayers[nTracks] = track1->hitPattern().trackerLayersWithMeasurement();  //layers crossed
0661       trkQual[nTracks] = track1->quality(reco::TrackBase::highPurity);
0662 
0663       nTracks++;
0664     }
0665 
0666     /*        ------------------          ------------------------------ ------- */
0667 
0668     int dieta_M_P = 100;
0669     int diphi_M_P = 100;
0670     if (MaxHit.ietahitm * ietatrue > 0) {
0671       dieta_M_P = abs(MaxHit.ietahitm - ietatrue);
0672     }
0673     if (MaxHit.ietahitm * ietatrue < 0) {
0674       dieta_M_P = abs(MaxHit.ietahitm - ietatrue) - 1;
0675     }
0676     diphi_M_P = abs(MaxHit.iphihitm - iphitrue);
0677 
0678     diphi_M_P = diphi_M_P > 36 ? 72 - diphi_M_P : diphi_M_P;
0679     iDr = sqrt(diphi_M_P * diphi_M_P + dieta_M_P * dieta_M_P);
0680 
0681 #ifdef EDM_ML_DEBUG
0682     if (iDr > 15) {
0683       edm::LogVerbatim("HcalCalib") << "diphi: " << diphi_M_P << "  dieta: " << dieta_M_P << "   iDr: " << iDr
0684                                     << " ietatrue:" << ietatrue << "  iphitrue:" << iphitrue;
0685       edm::LogVerbatim("HcalCalib") << "M ieta: " << MaxHit.ietahitm << "  M iphi: " << MaxHit.iphihitm;
0686     }
0687 #endif
0688 
0689     Bool_t passCuts = kFALSE;
0690     passCuts = kTRUE;
0691     //if(eEcalCone < energyECALmip && iDr<2.) passCuts = kTRUE;
0692     //if(MaxHit.hitenergy>0.) passCuts = kTRUE;
0693 
0694     if (passCuts) {
0695       /*      
0696       enHcal -> Fill(ietatrue,  eHcalCone);
0697       nCells -> Fill(ietatrue,  UsedCells);
0698       enHcalNoise -> Fill(ietatrue,  eHcalConeNoise);
0699       nCellsNoise -> Fill(ietatrue,  UsedCellsNoise); 
0700       */
0701 
0702       //e3x3=0; e5x5=0;
0703 
0704       iEta = ietatrue;
0705       iPhi = iphitrue;
0706 
0707       //iEta = MaxHit.ietahitm;
0708       //iPhi = MaxHit.iphihitm;
0709       delR = MaxHit.dr;
0710       eCentHit = MaxHit.hitenergy;
0711 
0712       eParticle = mom_MC;
0713       //eTrack = mom_MC;
0714       //phiTrack = phiParticle;
0715       //etaTrack = etaParticle;
0716 
0717       pfTree->Fill();
0718     }
0719 
0720   }  //Hep:MC
0721 }
0722 
0723 void HcalCorrPFCalculation::beginJob() {
0724   pfTree = fs->make<TTree>("pfTree", "Tree for pf info");
0725 
0726   pfTree->Branch("nTracks", &nTracks, "nTracks/I");
0727   pfTree->Branch("trackEta", trackEta, "trackEta[nTracks]/F");
0728   pfTree->Branch("trackPhi", trackPhi, "trackPhi[nTracks]/F");
0729   pfTree->Branch("trackP", trackP, "trackP[nTracks]/F");
0730 
0731   pfTree->Branch("delRmc", delRmc, "delRmc[nTracks]/F");
0732   pfTree->Branch("numValidTrkHits", numValidTrkHits, "numValidTrkHits[nTracks]/I");
0733   pfTree->Branch("numValidTrkStrips", numValidTrkStrips, "numValidTrkStrips[nTracks]/I");
0734   pfTree->Branch("numLayers", numLayers, "numLayers[nTracks]/I");
0735   pfTree->Branch("trkQual", trkQual, "trkQual[nTracks]/O");
0736 
0737   pfTree->Branch("eEcalCone", &eEcalCone, "eEcalCone/F");
0738   pfTree->Branch("eHcalCone", &eHcalCone, "eHcalCone/F");
0739   pfTree->Branch("eHcalConeNoise", &eHcalConeNoise, "eHcalConeNoise/F");
0740 
0741   pfTree->Branch("UsedCellsNoise", &UsedCellsNoise, "UsedCellsNoise/I");
0742   pfTree->Branch("UsedCells", &UsedCells, "UsedCells/I");
0743 
0744   pfTree->Branch("eCentHit", &eCentHit, "eCentHit/F");
0745 
0746   pfTree->Branch("eParticle", &eParticle, "eParticle/F");
0747   pfTree->Branch("etaParticle", &etaParticle, "etaParticle/F");
0748   pfTree->Branch("phiParticle", &phiParticle, "phiParticle/F");
0749 
0750   pfTree->Branch("etaGPoint", &etaGPoint, "etaGPoint/F");
0751   pfTree->Branch("phiGPoint", &phiGPoint, "phiGPoint/F");
0752 
0753   pfTree->Branch("xAtHcal", &xAtHcal, "xAtHcal/F");
0754   pfTree->Branch("yAtHcal", &yAtHcal, "yAtHcal/F");
0755   pfTree->Branch("zAtHcal", &zAtHcal, "zAtHcal/F");
0756 
0757   pfTree->Branch("eECAL09cm", &eECAL09cm, "eECAL09cm/F");
0758   pfTree->Branch("eECAL40cm", &eECAL40cm, "eECAL40cm/F");
0759   pfTree->Branch("eECAL", &eECAL, "eECAL/F");
0760 
0761   pfTree->Branch("e3x3 ", &e3x3, "e3x3/F");
0762   pfTree->Branch("e5x5", &e5x5, "e5x5/F");
0763 
0764   pfTree->Branch("iDr", &iDr, "iDr/F");
0765   pfTree->Branch("delR", &delR, "delR/F");
0766 
0767   pfTree->Branch("iEta", &iEta, "iEta/I");
0768   pfTree->Branch("iPhi", &iPhi, "iPhi/I");
0769 
0770   //  pfTree->Branch("numValidTrkHits", &numValidTrkHits, "numValidTrkHits/I");
0771   // pfTree->Branch("numValidTrkStrips", &numValidTrkStrips, "numValidTrkStrips/I");
0772   // pfTree->Branch("trkQual", &trkQual, "trkQual/");
0773   // pfTree->Branch("numLayers", &numLayers, "numLayers/I");
0774 }
0775 
0776 #include "FWCore/Framework/interface/MakerMacros.h"
0777 
0778 //define this as a plug-in
0779 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!