Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-10-25 09:34:41 

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