Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-02-28 01:32:04 

0001 //
0002 // Package:    ValidIsoTrkCalib
0003 // Class:      ValidIsoTrkCalib
0004 //
0005 
0006 /*
0007  Description: Validation for Isolated tracks Calibration
0008  
0009  Implementation:
0010 See the twiki page for details:
0011 https://twiki.cern.ch/twiki/bin/view/CMS/ValidIsoTrkCalib
0012 
0013 */
0014 
0015 //
0016 // Original Author:  Andrey Pozdnyakov
0017 //         Created:  Tue Nov  4 01:16:05 CET 2008
0018 //
0019 
0020 // system include files
0021 
0022 #include <memory>
0023 
0024 // user include files
0025 #include "FWCore/Framework/interface/Frameworkfwd.h"
0026 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0027 #include "FWCore/Framework/interface/Event.h"
0028 #include "FWCore/Framework/interface/MakerMacros.h"
0029 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0030 #include "FWCore/ServiceRegistry/interface/Service.h"
0031 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0032 
0033 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0034 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0035 #include "TrackingTools/TrackAssociator/interface/TrackDetectorAssociator.h"
0036 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0037 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0038 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0039 #include "DataFormats/HcalIsolatedTrack/interface/IsolatedPixelTrackCandidate.h"
0040 #include "DataFormats/HcalIsolatedTrack/interface/IsolatedPixelTrackCandidateFwd.h"
0041 
0042 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0043 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0044 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0045 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0046 #include "Geometry/Records/interface/CaloTopologyRecord.h"
0047 #include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
0048 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0049 #include "CondFormats/HcalObjects/interface/HcalRespCorrs.h"
0050 #include "CondFormats/DataRecord/interface/HcalRespCorrsRcd.h"
0051 
0052 #include "Calibration/HcalCalibAlgos/interface/CommonUsefulStuff.h"
0053 
0054 #include "TROOT.h"
0055 #include "TFile.h"
0056 #include "TTree.h"
0057 #include "TH1F.h"
0058 
0059 #include <fstream>
0060 #include <map>
0061 
0062 //#define EDM_ML_DEBUG
0063 
0064 class ValidIsoTrkCalib : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0065 public:
0066   explicit ValidIsoTrkCalib(const edm::ParameterSet&);
0067 
0068   //  double getDistInPlaneSimple(const GlobalPoint caloPoint, const GlobalPoint rechitPoint);
0069 
0070 private:
0071   void beginJob() override;
0072   void analyze(const edm::Event&, const edm::EventSetup&) override;
0073 
0074   // ----------member data ---------------------------
0075 
0076   //Variables from HcalIsotrackAnalyzer
0077 
0078   TrackDetectorAssociator trackAssociator_;
0079   TrackAssociatorParameters parameters_;
0080   double taECALCone_;
0081   double taHCALCone_;
0082 
0083   const CaloGeometry* geo;
0084   // nothing is done with these tags, so I leave it - cowden
0085   const bool takeGenTracks_;
0086   const edm::InputTag genhbheLabel_;
0087   const double associationConeSize_;
0088   const bool allowMissingInputs_;
0089   const std::string AxB_;
0090   const double calibrationConeSize_;
0091   const int MinNTrackHitsBarrel;
0092   const int MinNTECHitsEndcap;
0093   const double energyECALmip;
0094   const double energyMinIso;
0095   const double energyMaxIso;
0096   const double maxPNear;
0097 
0098   const edm::EDGetTokenT<reco::TrackCollection> tok_genTrack_;
0099   const edm::EDGetTokenT<HBHERecHitCollection> tok_hbhe_;
0100   const edm::EDGetTokenT<HORecHitCollection> tok_ho_;
0101   const edm::EDGetTokenT<reco::IsolatedPixelTrackCandidateCollection> tok_track_;
0102   const edm::EDGetTokenT<reco::TrackCollection> tok_track1_;
0103 
0104   const edm::ESGetToken<HcalRespCorrs, HcalRespCorrsRcd> tok_recalibCorrs_;
0105   const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tok_geom_;
0106 
0107   edm::InputTag genhoLabel_;
0108   std::vector<edm::InputTag> genecalLabel_;
0109   std::string outputFileName_;
0110 
0111   int gen, iso, pix;
0112   float genPt[500], genPhi[500], genEta[500];
0113   float isoPt[500], isoPhi[500], isoEta[500];
0114   float pixPt[500], pixPhi[500], pixEta[500];
0115 
0116   int NisoTrk;
0117   float trackPt, trackE, trackEta, trackPhi;
0118   float ptNear;
0119   float ptrack, rvert;
0120   //float eecal, ehcal;
0121 
0122   Float_t emEnergy;
0123   Float_t targetE;
0124 
0125   TTree *tTree, *fTree;
0126 
0127   Float_t xTrkEcal;
0128   Float_t yTrkEcal;
0129   Float_t zTrkEcal;
0130 
0131   Float_t xTrkHcal;
0132   Float_t yTrkHcal;
0133   Float_t zTrkHcal;
0134 
0135   int Nhits;
0136   float eClustBefore;  //Calo energy before calibration
0137   float eClustAfter;   //After calibration
0138   float eTrack;        //Track energy
0139   float etaTrack;
0140   float phiTrack;
0141   float eECAL;  // Energy deposited in ECAL
0142   int numHits;  //number of rechits
0143 
0144   float eBeforeDepth1;
0145   float eAfterDepth1;
0146   float eBeforeDepth2;
0147   float eAfterDepth2;
0148   float eCentHitBefore;
0149   float eCentHitAfter;
0150   float CentHitFactor;
0151   int iEta;
0152   int iPhi;
0153   int iEtaTr;
0154   int iPhiTr;
0155   float iDr, delR;
0156   int dietatr;
0157   int diphitr;
0158 
0159   float iTime;
0160   float HTime[100];
0161   float e3x3Before;
0162   float e3x3After;
0163   float e5x5Before;
0164   float e5x5After;
0165   int eventNumber;
0166   int runNumber;
0167   float PtNearBy;
0168   float numVH, numVS, numValidTrkHits, numValidTrkStrips;
0169 
0170   const HcalRespCorrs* respRecalib;
0171 
0172   TH1F* nTracks;
0173 };
0174 
0175 ValidIsoTrkCalib::ValidIsoTrkCalib(const edm::ParameterSet& iConfig)
0176     : takeGenTracks_(iConfig.getUntrackedParameter<bool>("takeGenTracks")),
0177       genhbheLabel_(iConfig.getParameter<edm::InputTag>("genHBHE")),
0178       associationConeSize_(iConfig.getParameter<double>("associationConeSize")),
0179       allowMissingInputs_(iConfig.getUntrackedParameter<bool>("allowMissingInputs", true)),
0180       AxB_(iConfig.getParameter<std::string>("AxB")),
0181       calibrationConeSize_(iConfig.getParameter<double>("calibrationConeSize")),
0182       MinNTrackHitsBarrel(iConfig.getParameter<int>("MinNTrackHitsBarrel")),
0183       MinNTECHitsEndcap(iConfig.getParameter<int>("MinNTECHitsEndcap")),
0184       energyECALmip(iConfig.getParameter<double>("energyECALmip")),
0185       energyMinIso(iConfig.getParameter<double>("energyMinIso")),
0186       energyMaxIso(iConfig.getParameter<double>("energyMaxIso")),
0187       maxPNear(iConfig.getParameter<double>("maxPNear")),
0188       tok_genTrack_(consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("genTracksLabel"))),
0189       tok_hbhe_(consumes<HBHERecHitCollection>(iConfig.getParameter<edm::InputTag>("hbheInput"))),
0190       tok_ho_(consumes<HORecHitCollection>(iConfig.getParameter<edm::InputTag>("hoInput"))),
0191       tok_track_(consumes<reco::IsolatedPixelTrackCandidateCollection>(
0192           iConfig.getParameter<edm::InputTag>("HcalIsolTrackInput"))),
0193       tok_track1_(consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("trackInput"))),
0194       tok_recalibCorrs_(esConsumes(edm::ESInputTag("", "recalibrate"))),
0195       tok_geom_(esConsumes()) {
0196   usesResource(TFileService::kSharedResource);
0197 
0198   edm::ParameterSet parameters = iConfig.getParameter<edm::ParameterSet>("TrackAssociatorParameters");
0199   edm::ConsumesCollector iC = consumesCollector();
0200   parameters_.loadParameters(parameters, iC);
0201   trackAssociator_.useDefaultPropagator();
0202 }
0203 
0204 // ------------ method called to for each event  ------------
0205 void ValidIsoTrkCalib::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0206   try {
0207     respRecalib = &iSetup.getData(tok_recalibCorrs_);
0208 
0209     edm::LogVerbatim("ValidIsoTrkCalib") << "  Loaded:  OK ";
0210 
0211   } catch (const cms::Exception& e) {
0212     edm::LogWarning("ValidIsoTrkCalib") << "   Not Found!! ";
0213   }
0214 
0215   const edm::Handle<reco::TrackCollection>& generalTracks = iEvent.getHandle(tok_genTrack_);
0216 
0217   const edm::Handle<reco::TrackCollection>& isoProdTracks = iEvent.getHandle(tok_track1_);
0218 
0219   const edm::Handle<reco::IsolatedPixelTrackCandidateCollection>& isoPixelTracks = iEvent.getHandle(tok_track_);
0220   //edm::Handle<reco::TrackCollection> isoPixelTracks;
0221 
0222   /*
0223   edm::Handle<EcalRecHitCollection> ecal;
0224   iEvent.getByLabel(eLabel_,ecal);
0225   const EcalRecHitCollection Hitecal = *(ecal.product());
0226   */
0227 
0228   const edm::Handle<HBHERecHitCollection>& hbhe = iEvent.getHandle(tok_hbhe_);
0229   const HBHERecHitCollection Hithbhe = *(hbhe.product());
0230 
0231   geo = &iSetup.getData(tok_geom_);
0232 
0233   const HcalGeometry* gHcal = static_cast<const HcalGeometry*>(geo->getSubdetectorGeometry(DetId::Hcal, HcalBarrel));
0234   //Note: even though it says HcalBarrel, we actually get the whole Hcal detector geometry!
0235 
0236   // Lumi_n=iEvent.luminosityBlock();
0237   parameters_.useEcal = true;
0238   parameters_.useHcal = true;
0239   parameters_.useCalo = false;
0240   parameters_.useMuon = false;
0241   //parameters_.dREcal = taECALCone_;
0242   //parameters_.dRHcal = taHCALCone_;
0243 #ifdef EDM_ML_DEBUG
0244   edm::LogVerbatim("ValidIsoTrkCalib") << "Hello World. TrackCollectionSize: " << isoPixelTracks->size();
0245 #endif
0246   if (isoPixelTracks->empty())
0247     return;
0248 
0249   for (reco::TrackCollection::const_iterator trit = isoProdTracks->begin(); trit != isoProdTracks->end(); trit++) {
0250     reco::IsolatedPixelTrackCandidateCollection::const_iterator isoMatched = isoPixelTracks->begin();
0251     //reco::TrackCollection::const_iterator isoMatched=isoPixelTracks->begin();
0252     bool matched = false;
0253 
0254     //for (reco::IsolatedPixelTrackCandidateCollection::const_iterator trit = isoPixelTracks->begin(); trit!=isoPixelTracks->end(); trit++)
0255     for (reco::IsolatedPixelTrackCandidateCollection::const_iterator it = isoPixelTracks->begin();
0256          it != isoPixelTracks->end();
0257          it++)
0258     //for (reco::TrackCollection::const_iterator it = isoPixelTracks->begin(); it!=isoPixelTracks->end(); it++)
0259     {
0260       if (abs((trit->pt() - it->pt()) / it->pt()) < 0.005 && abs(trit->eta() - it->eta()) < 0.01) {
0261         isoMatched = it;
0262         matched = true;
0263         break;
0264       }
0265     }
0266     // CUT
0267 
0268     if (!matched)
0269       continue;
0270     if (isoMatched->maxPtPxl() > maxPNear)
0271       continue;
0272 
0273     ptNear = isoMatched->maxPtPxl();
0274 #ifdef EDM_ML_DEBUG
0275     edm::LogVerbatim("ValidIsoTrkCalib") << "Point 0.1  isoMatch. ptnear: " << ptNear;
0276 #endif
0277     // CUT
0278     if (trit->hitPattern().numberOfValidHits() < MinNTrackHitsBarrel)
0279       continue;
0280     if (fabs(trit->eta()) > 1.47 && trit->hitPattern().numberOfValidStripTECHits() < MinNTECHitsEndcap)
0281       continue;
0282 
0283 #ifdef EDM_ML_DEBUG
0284     edm::LogVerbatim("ValidIsoTrkCalib") << "Point 0.2.1 after numofvalidhits HB: "
0285                                          << trit->hitPattern().numberOfValidHits();
0286     edm::LogVerbatim("ValidIsoTrkCalib") << "Point 0.2.2 after numofvalidstrips HE: "
0287                                          << trit->hitPattern().numberOfValidStripTECHits();
0288 #endif
0289     numVH = trit->hitPattern().numberOfValidHits();
0290     numVS = trit->hitPattern().numberOfValidStripTECHits();
0291 
0292     trackE = sqrt(trit->px() * trit->px() + trit->py() * trit->py() + trit->pz() * trit->pz() + 0.14 * 0.14);
0293     trackPt = trit->pt();
0294     trackEta = trit->eta();
0295     trackPhi = trit->phi();
0296 
0297     emEnergy = isoMatched->energyIn();
0298 
0299 #ifdef EDM_ML_DEBUG
0300     edm::LogVerbatim("ValidIsoTrkCalib") << "Point 0.3.  Matched :: pt: " << trit->pt() << " wholeEnergy: " << trackE
0301                                          << "  emEnergy: " << emEnergy << "  eta: " << trackEta << " phi: " << trackPhi;
0302     edm::LogVerbatim("ValidIsoTrkCalib") << "Point 0.4.  EM energy in cone: " << emEnergy;
0303 #endif
0304     TrackDetMatchInfo info = trackAssociator_.associate(
0305         iEvent,
0306         iSetup,
0307         trackAssociator_.getFreeTrajectoryState(&iSetup.getData(parameters_.bFieldToken), *trit),
0308         parameters_);
0309 
0310     //float etaecal=info.trkGlobPosAtEcal.eta();
0311     //float phiecal=info.trkGlobPosAtEcal.phi();
0312     //  float etahcal=info.trkGlobPosAtHcal.eta();
0313     // float phihcal=info.trkGlobPosAtHcal.phi();
0314 
0315     xTrkEcal = info.trkGlobPosAtEcal.x();
0316     yTrkEcal = info.trkGlobPosAtEcal.y();
0317     zTrkEcal = info.trkGlobPosAtEcal.z();
0318 
0319     xTrkHcal = info.trkGlobPosAtHcal.x();
0320     yTrkHcal = info.trkGlobPosAtHcal.y();
0321     zTrkHcal = info.trkGlobPosAtHcal.z();
0322 
0323     if (xTrkEcal == 0 && yTrkEcal == 0 && zTrkEcal == 0) {
0324       edm::LogVerbatim("ValidIsoTrkCalib") << "zero point at Ecal";
0325       continue;
0326     }
0327     if (xTrkHcal == 0 && yTrkHcal == 0 && zTrkHcal == 0) {
0328       edm::LogVerbatim("ValidIsoTrkCalib") << "zero point at Hcal";
0329       continue;
0330     }
0331 
0332     /*GlobalVector trackMomAtEcal = info.trkMomAtEcal;
0333       GlobalVector trackMomAtHcal = info.trkMomAtHcal;
0334     
0335       PxTrkHcal = trackMomAtHcal.x();
0336       PyTrkHcal = trackMomAtHcal.y();
0337       PzTrkHcal = trackMomAtHcal.z();
0338       */
0339 
0340     GlobalPoint gPointEcal(xTrkEcal, yTrkEcal, zTrkEcal);
0341     GlobalPoint gPointHcal(xTrkHcal, yTrkHcal, zTrkHcal);
0342 
0343     int iphitrue = -10;
0344     int ietatrue = 100;
0345     const HcalDetId tempId = gHcal->getClosestCell(gPointHcal);
0346     ietatrue = tempId.ieta();
0347     iphitrue = tempId.iphi();
0348 
0349     MaxHit_struct MaxHit;
0350 
0351     MaxHit.hitenergy = -100.;
0352 
0353     //container for used recHits
0354     std::vector<DetId> usedHits;
0355     //
0356     usedHits.clear();
0357 #ifdef EDM_ML_DEBUG
0358     edm::LogVerbatim("ValidIsoTrkCalib") << "Point 1. Entrance to HBHECollection";
0359 #endif
0360     //float dddeta = 1000.;
0361     //float dddphi = 1000.;
0362     //int iphitrue = 1234;
0363     //int ietatrue = 1234;
0364 
0365     GlobalPoint gPhot;
0366 
0367     for (HBHERecHitCollection::const_iterator hhit = Hithbhe.begin(); hhit != Hithbhe.end(); hhit++) {
0368       //check that this hit was not considered before and push it into usedHits
0369       bool hitIsUsed = false;
0370       for (uint32_t i = 0; i < usedHits.size(); i++) {
0371         if (usedHits[i] == hhit->id())
0372           hitIsUsed = true;
0373       }
0374       if (hitIsUsed)
0375         continue;
0376       usedHits.push_back(hhit->id());
0377       //
0378 
0379       // rof 16.05.2008 start: include the possibility for recalibration
0380       float recal = 1;
0381       // rof end
0382 
0383       GlobalPoint pos = geo->getPosition(hhit->detid());
0384       //float phihit = pos.phi();
0385       //float etahit = pos.eta();
0386 
0387       int iphihitm = (hhit->id()).iphi();
0388       int ietahitm = (hhit->id()).ieta();
0389       int depthhit = (hhit->id()).depth();
0390       float enehit = hhit->energy() * recal;
0391 
0392       if (depthhit != 1)
0393         continue;
0394 
0395 #ifdef EDM_ML_DEBUG
0396       float dphi = fabs(info.trkGlobPosAtHcal.phi() - pos.phi());
0397       if (dphi > 4. * atan(1.))
0398         dphi = 8. * atan(1.) - dphi;
0399       float deta = fabs(info.trkGlobPosAtHcal.eta() - pos.eta());
0400       float dr = sqrt(dphi * dphi + deta * deta);
0401 #endif
0402 
0403       //double distAtHcal =  getDistInPlaneTrackDir(gPointHcal, trackMomAtHcal, pos);
0404       double distAtHcal = getDistInPlaneSimple(gPointHcal, pos);
0405 
0406       if (distAtHcal < associationConeSize_) {
0407         for (HBHERecHitCollection::const_iterator hhit2 = Hithbhe.begin(); hhit2 != Hithbhe.end(); hhit2++) {
0408           int iphihitm2 = (hhit2->id()).iphi();
0409           int ietahitm2 = (hhit2->id()).ieta();
0410           int depthhit2 = (hhit2->id()).depth();
0411           float enehit2 = hhit2->energy() * recal;
0412 
0413           if (iphihitm == iphihitm2 && ietahitm == ietahitm2 && depthhit != depthhit2)
0414             enehit = enehit + enehit2;
0415         }
0416 
0417 #ifdef EDM_ML_DEBUG
0418         edm::LogVerbatim("ValidIsoTrkCalib") << "IN CONE ieta: " << ietahitm << "  iphi: " << iphihitm
0419                                              << " depthhit: " << depthhit << "  dr: " << dr << " energy: " << enehit;
0420 #endif
0421         //Find a Hit with Maximum Energy
0422 
0423         if (enehit > MaxHit.hitenergy) {
0424           MaxHit.hitenergy = enehit;
0425           MaxHit.ietahitm = (hhit->id()).ieta();
0426           MaxHit.iphihitm = (hhit->id()).iphi();
0427           MaxHit.dr = distAtHcal;
0428           //MaxHit.depthhit  = (hhit->id()).depth();
0429           MaxHit.depthhit = 1;
0430 
0431           //gPhot = geo->getPosition(hhit->detid());
0432         }
0433       }
0434     }  //end of all HBHE hits cycle
0435 
0436     usedHits.clear();
0437 
0438 #ifdef EDM_ML_DEBUG
0439     edm::LogVerbatim("ValidIsoTrkCalib") << "Hottest ieta: " << MaxHit.ietahitm << "  iphi: " << MaxHit.iphihitm
0440                                          << "  dr: " << MaxHit.dr;
0441     edm::LogVerbatim("ValidIsoTrkCalib") << "Track   ieta: " << ietatrue << "  iphi: " << iphitrue;
0442     edm::LogVerbatim("ValidIsoTrkCalib") << "Point 3.  MaxHit :::En " << MaxHit.hitenergy
0443                                          << "  ieta: " << MaxHit.ietahitm << "  iphi: " << MaxHit.iphihitm;
0444 #endif
0445 
0446     Bool_t passCuts = kFALSE;
0447     if (trackE > energyMinIso && trackE < energyMaxIso && emEnergy < energyECALmip && MaxHit.hitenergy > 0. &&
0448         abs(MaxHit.ietahitm) < 29)
0449       passCuts = kTRUE;
0450 
0451 #ifdef EDM_ML_DEBUG
0452     edm::LogVerbatim("ValidIsoTrkCalib") << "Pont 0.1.1.  trackE:" << trackE << "  emEn: " << emEnergy;
0453 #endif
0454     numHits = 0;
0455 
0456     eClustBefore = 0.0;
0457     eClustAfter = 0.0;
0458     eBeforeDepth1 = 0.0;
0459     eAfterDepth1 = 0.0;
0460     eBeforeDepth2 = 0.0;
0461     eAfterDepth2 = 0.0;
0462     CentHitFactor = 0.0;
0463     e3x3After = 0.0;
0464     e3x3Before = 0.0;
0465     e5x5After = 0.0;
0466     e5x5Before = 0.0;
0467 
0468     for (HBHERecHitCollection::const_iterator hhit = Hithbhe.begin(); hhit != Hithbhe.end(); hhit++) {
0469       //check that this hit was not considered before and push it into usedHits
0470       bool hitIsUsed = false;
0471       for (uint32_t i = 0; i < usedHits.size(); i++) {
0472         if (usedHits[i] == hhit->id())
0473           hitIsUsed = true;
0474       }
0475       if (hitIsUsed)
0476         continue;
0477       usedHits.push_back(hhit->id());
0478 
0479       int DIETA = 100;
0480       if (MaxHit.ietahitm * (hhit->id()).ieta() > 0) {
0481         DIETA = MaxHit.ietahitm - (hhit->id()).ieta();
0482       }
0483       if (MaxHit.ietahitm * (hhit->id()).ieta() < 0) {
0484         DIETA = MaxHit.ietahitm - (hhit->id()).ieta();
0485         DIETA = DIETA > 0 ? DIETA - 1 : DIETA + 1;
0486       }
0487 
0488       int DIPHI = abs(MaxHit.iphihitm - (hhit->id()).iphi());
0489       DIPHI = DIPHI > 36 ? 72 - DIPHI : DIPHI;
0490 
0491       int numbercell = 100;  //always collect Wide clastor!
0492 
0493       //if(AxB_=="5x5" || AxB_=="3x3" || AxB_=="7x7"|| AxB_=="Cone")
0494 
0495       //if(AxB_=="3x3") numbercell = 1;
0496       //if(AxB_=="5x5") numbercell = 2;
0497       //if(AxB_=="Cone") numbercell = 1000;
0498 
0499       if (abs(DIETA) <= numbercell &&
0500           (abs(DIPHI) <= numbercell || (abs(MaxHit.ietahitm) >= 20 && abs(DIPHI) <= numbercell + 1))) {
0501         const GlobalPoint pos2 = geo->getPosition(hhit->detid());
0502 
0503         if (passCuts && hhit->energy() > 0) {
0504           float factor = 0.0;
0505           // factor = CalibFactors[hhit->id()];
0506           factor = respRecalib->getValues(hhit->id())->getValue();
0507 
0508 #ifdef EDM_ML_DEBUG
0509           edm::LogVerbatim("ValidIsoTrkCalib") << " calib factors: " << factor;
0510 #endif
0511 
0512           if (hhit->id().ieta() == MaxHit.ietahitm && hhit->id().iphi() == MaxHit.iphihitm)
0513             CentHitFactor = factor;
0514 
0515           if (hhit->id().ieta() == MaxHit.ietahitm && hhit->id().iphi() == MaxHit.iphihitm)
0516             iTime = hhit->time();
0517 
0518           if (AxB_ != "3x3" && AxB_ != "5x5" && AxB_ != "Cone")
0519             edm::LogWarning(" AxB ") << "   Not supported: " << AxB_;
0520 
0521           if (abs(DIETA) <= 2 && (abs(DIPHI) <= 2 || ((abs(MaxHit.ietahitm) > 20 && abs(DIPHI) <= 4) &&
0522                                                       !((abs(MaxHit.ietahitm) == 21 || abs(MaxHit.ietahitm) == 22) &&
0523                                                         abs((hhit->id()).ieta()) <= 20 && abs(DIPHI) > 2)))) {
0524             e5x5Before += hhit->energy();
0525             e5x5After += hhit->energy() * factor;
0526           }
0527 
0528           if (abs(DIETA) <= 1 && (abs(DIPHI) <= 1 || ((abs(MaxHit.ietahitm) > 20 && abs(DIPHI) <= 2) &&
0529                                                       !(abs(MaxHit.ietahitm) == 21 && abs((hhit->id()).ieta()) <= 20 &&
0530                                                         abs(DIPHI) > 1)))) {
0531             e3x3Before += hhit->energy();
0532             e3x3After += hhit->energy() * factor;
0533           }
0534 
0535           if (AxB_ == "5x5") {
0536             if (abs(DIETA) <= 2 && (abs(DIPHI) <= 2 || (abs(MaxHit.ietahitm) > 20 && abs(DIPHI) <= 4))) {
0537               if (abs(MaxHit.ietahitm) == 21 && abs((hhit->id()).ieta()) <= 20 && abs(DIPHI) > 3)
0538                 continue;
0539 
0540               HTime[numHits] = hhit->time();
0541               numHits++;
0542 
0543               eClustBefore += hhit->energy();
0544               eClustAfter += hhit->energy() * factor;
0545 
0546               if ((hhit->id().depth() == 1) && (abs(hhit->id().ieta()) > 17) && (abs(hhit->id().ieta()) < 29)) {
0547                 eBeforeDepth1 += hhit->energy();
0548                 eAfterDepth1 += hhit->energy() * factor;
0549               } else if ((hhit->id().depth() == 2) && (abs(hhit->id().ieta()) > 17) && (abs(hhit->id().ieta()) < 29)) {
0550                 eBeforeDepth2 += hhit->energy();
0551                 eAfterDepth2 += hhit->energy() * factor;
0552               }
0553             }
0554           }  //end of 5x5
0555 
0556           if (AxB_ == "3x3") {
0557             if (abs(DIETA) <= 1 && (abs(DIPHI) <= 1 || (abs(MaxHit.ietahitm) > 20 && abs(DIPHI) <= 2))) {
0558               if (abs(MaxHit.ietahitm) == 21 && abs((hhit->id()).ieta()) <= 20 && abs(DIPHI) > 2)
0559                 continue;
0560 
0561               HTime[numHits] = hhit->time();
0562               numHits++;
0563 
0564               eClustBefore += hhit->energy();
0565               eClustAfter += hhit->energy() * factor;
0566 
0567               if ((hhit->id().depth() == 1) && (abs(hhit->id().ieta()) > 17) && (abs(hhit->id().ieta()) < 29)) {
0568                 eBeforeDepth1 += hhit->energy();
0569                 eAfterDepth1 += hhit->energy() * factor;
0570               } else if ((hhit->id().depth() == 2) && (abs(hhit->id().ieta()) > 17) && (abs(hhit->id().ieta()) < 29)) {
0571                 eBeforeDepth2 += hhit->energy();
0572                 eAfterDepth2 += hhit->energy() * factor;
0573               }
0574             }
0575           }  //end of 3x3
0576 
0577           if (AxB_ == "Cone" && getDistInPlaneSimple(gPointHcal, pos2) < calibrationConeSize_) {
0578             HTime[numHits] = hhit->time();
0579             numHits++;
0580 
0581             eClustBefore += hhit->energy();
0582             eClustAfter += hhit->energy() * factor;
0583 
0584             if ((hhit->id().depth() == 1) && (abs(hhit->id().ieta()) > 17) && (abs(hhit->id().ieta()) < 29)) {
0585               eBeforeDepth1 += hhit->energy();
0586               eAfterDepth1 += hhit->energy() * factor;
0587             } else if ((hhit->id().depth() == 2) && (abs(hhit->id().ieta()) > 17) && (abs(hhit->id().ieta()) < 29)) {
0588               eBeforeDepth2 += hhit->energy();
0589               eAfterDepth2 += hhit->energy() * factor;
0590             }
0591 
0592           }  //end of Cone
0593 
0594         }  //end of passCuts
0595 
0596       }  //end of DIETA DIPHI
0597 
0598     }  //end of associatedcone HBHE hits cycle
0599 
0600     int dieta_M_P = 100;
0601     int diphi_M_P = 100;
0602     if (MaxHit.ietahitm * ietatrue > 0) {
0603       dieta_M_P = abs(MaxHit.ietahitm - ietatrue);
0604     }
0605     if (MaxHit.ietahitm * ietatrue < 0) {
0606       dieta_M_P = abs(MaxHit.ietahitm - ietatrue) - 1;
0607     }
0608     diphi_M_P = abs(MaxHit.iphihitm - iphitrue);
0609     diphi_M_P = diphi_M_P > 36 ? 72 - diphi_M_P : diphi_M_P;
0610 
0611     if (passCuts)
0612 
0613     {
0614       eventNumber = iEvent.id().event();
0615       runNumber = iEvent.id().run();
0616 
0617       eCentHitBefore = MaxHit.hitenergy;
0618       eCentHitAfter = MaxHit.hitenergy * CentHitFactor;
0619       eECAL = emEnergy;
0620       numValidTrkHits = numVH;
0621       numValidTrkStrips = numVS;
0622       PtNearBy = ptNear;
0623 
0624       eTrack = trackE;
0625       phiTrack = trackPhi;
0626       etaTrack = trackEta;
0627 
0628       iEta = MaxHit.ietahitm;
0629       iPhi = MaxHit.iphihitm;
0630 
0631       iEtaTr = ietatrue;
0632       iPhiTr = iphitrue;
0633       iDr = sqrt(diphi_M_P * diphi_M_P + dieta_M_P * dieta_M_P);
0634       delR = MaxHit.dr;
0635       dietatr = dieta_M_P;
0636       diphitr = diphi_M_P;
0637 
0638       fTree->Fill();
0639     }
0640 
0641   }  //end of isoProdTracks cycle
0642 
0643   /* ------------------   Some stuff for general tracks  ----------   ----*/
0644 #ifdef EDM_ML_DEBUG
0645   edm::LogVerbatim("ValidIsoTrkCalib") << " generalTracks Size: " << generalTracks->size();
0646 #endif
0647   int n = generalTracks->size();
0648   nTracks->Fill(n);
0649 
0650   if (takeGenTracks_ && iEvent.id().event() % 10 == 1) {
0651     gen = generalTracks->size();
0652     iso = isoProdTracks->size();
0653     pix = isoPixelTracks->size();
0654 
0655     genPt[0] = -33;
0656     genPhi[0] = -33;
0657     genEta[0] = -33;
0658 
0659     isoPt[0] = -33;
0660     isoPhi[0] = -33;
0661     isoEta[0] = -33;
0662 
0663     pixPt[0] = -33;
0664     pixPhi[0] = -33;
0665     pixEta[0] = -33;
0666 
0667     Int_t gencount = 0, isocount = 0, pixcount = 0;
0668     for (reco::TrackCollection::const_iterator gentr = generalTracks->begin(); gentr != generalTracks->end(); gentr++) {
0669       genPt[gencount] = gentr->pt();
0670       genPhi[gencount] = gentr->phi();
0671       genEta[gencount] = gentr->eta();
0672       gencount++;
0673     }
0674 
0675     for (reco::TrackCollection::const_iterator isotr = isoProdTracks->begin(); isotr != isoProdTracks->end(); isotr++) {
0676       isoPt[isocount] = isotr->pt();
0677       isoPhi[isocount] = isotr->phi();
0678       isoEta[isocount] = isotr->eta();
0679       isocount++;
0680     }
0681 
0682     for (reco::IsolatedPixelTrackCandidateCollection::const_iterator pixtr = isoPixelTracks->begin();
0683          pixtr != isoPixelTracks->end();
0684          pixtr++) {
0685       pixPt[pixcount] = pixtr->pt();
0686       pixPhi[pixcount] = pixtr->phi();
0687       pixEta[pixcount] = pixtr->eta();
0688       pixcount++;
0689     }
0690   }
0691 
0692   tTree->Fill();
0693 }
0694 
0695 // ------------ method called once each job just before starting event loop  ------------
0696 void ValidIsoTrkCalib::beginJob() {
0697   edm::Service<TFileService> fs;
0698 
0699   nTracks = fs->make<TH1F>("nTracks", "general;number of general tracks", 11, -0.5, 10.5);
0700 
0701   tTree = fs->make<TTree>("tTree", "Tree for gen info");
0702 
0703   fTree = fs->make<TTree>("fTree", "Tree for IsoTrack Calibration");
0704 
0705   fTree->Branch("eventNumber", &eventNumber, "eventNumber/I");
0706   fTree->Branch("runNumber", &runNumber, "runNumber/I");
0707 
0708   fTree->Branch("eClustBefore", &eClustBefore, "eClustBefore/F");
0709   fTree->Branch("eClustAfter", &eClustAfter, "eClustAfter/F");
0710   fTree->Branch("eTrack", &eTrack, "eTrack/F");
0711   fTree->Branch("etaTrack", &etaTrack, "etaTrack/F");
0712   fTree->Branch("phiTrack", &phiTrack, "phiTrack/F");
0713 
0714   fTree->Branch("numHits", &numHits, "numHits/I");
0715   fTree->Branch("eECAL", &eECAL, "eECAL/F");
0716   fTree->Branch("PtNearBy", &PtNearBy, "PtNearBy/F");
0717   fTree->Branch("numValidTrkHits", &numValidTrkHits, "numValidTrkHits/F");
0718   fTree->Branch("numValidTrkStrips", &numValidTrkStrips, "numValidTrkStrips/F");
0719 
0720   fTree->Branch("eBeforeDepth1", &eBeforeDepth1, "eBeforeDepth1/F");
0721   fTree->Branch("eBeforeDepth2", &eBeforeDepth2, "eBeforeDepth2/F");
0722   fTree->Branch("eAfterDepth1", &eAfterDepth1, "eAfterDepth1/F");
0723   fTree->Branch("eAfterDepth2", &eAfterDepth2, "eAfterDepth2/F");
0724 
0725   fTree->Branch("e3x3Before", &e3x3Before, "e3x3Before/F");
0726   fTree->Branch("e3x3After", &e3x3After, "e3x3After/F");
0727   fTree->Branch("e5x5Before", &e5x5Before, "e5x5Before/F");
0728   fTree->Branch("e5x5After", &e5x5After, "e5x5After/F");
0729 
0730   fTree->Branch("eCentHitAfter", &eCentHitAfter, "eCentHitAfter/F");
0731   fTree->Branch("eCentHitBefore", &eCentHitBefore, "eCentHitBefore/F");
0732   fTree->Branch("iEta", &iEta, "iEta/I");
0733   fTree->Branch("iPhi", &iPhi, "iPhi/I");
0734 
0735   fTree->Branch("iEtaTr", &iEtaTr, "iEtaTr/I");
0736   fTree->Branch("iPhiTr", &iPhiTr, "iPhiTr/I");
0737   fTree->Branch("dietatr", &dietatr, "dietatr/I");
0738   fTree->Branch("diphitr", &diphitr, "diphitr/I");
0739   fTree->Branch("iDr", &iDr, "iDr/F");
0740   fTree->Branch("delR", &delR, "delR/F");
0741 
0742   fTree->Branch("iTime", &iTime, "iTime/F");
0743   fTree->Branch("HTime", HTime, "HTime[numHits]/F");
0744 
0745   fTree->Branch("xTrkEcal", &xTrkEcal, "xTrkEcal/F");
0746   fTree->Branch("yTrkEcal", &yTrkEcal, "yTrkEcal/F");
0747   fTree->Branch("zTrkEcal", &zTrkEcal, "zTrkEcal/F");
0748   fTree->Branch("xTrkHcal", &xTrkHcal, "xTrkHcal/F");
0749   fTree->Branch("yTrkHcal", &yTrkHcal, "yTrkHcal/F");
0750   fTree->Branch("zTrkHcal", &zTrkHcal, "zTrkHcal/F");
0751 
0752   if (takeGenTracks_) {
0753     tTree->Branch("gen", &gen, "gen/I");
0754     tTree->Branch("iso", &iso, "iso/I");
0755     tTree->Branch("pix", &pix, "pix/I");
0756     tTree->Branch("genPt", genPt, "genPt[gen]/F");
0757     tTree->Branch("genPhi", genPhi, "genPhi[gen]/F");
0758     tTree->Branch("genEta", genEta, "genEta[gen]/F");
0759 
0760     tTree->Branch("isoPt", isoPt, "isoPt[iso]/F");
0761     tTree->Branch("isoPhi", isoPhi, "isoPhi[iso]/F");
0762     tTree->Branch("isoEta", isoEta, "isoEta[iso]/F");
0763 
0764     tTree->Branch("pixPt", pixPt, "pixPt[pix]/F");
0765     tTree->Branch("pixPhi", pixPhi, "pixPhi[pix]/F");
0766     tTree->Branch("pixEta", pixEta, "pixEta[pix]/F");
0767   }
0768 }
0769 
0770 //define this as a plug-in
0771 DEFINE_FWK_MODULE(ValidIsoTrkCalib);
[ 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!