OfflineValidation/plugins/EopTreeWriter.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    EopTreeWriter
0004 // Class:      EopTreeWriter
0005 //
0006 /**\class EopTreeWriter EopTreeWriter.cc Alignment/OfflineValidation/plugins/EopTreeWriter.cc
0007
0008  Description: <one line class summary>
0009
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Holger Enderle
0015 //         Created:  Thu Dec  4 11:22:48 CET 2008
0016 // $Id: EopTreeWriter.cc,v 1.2 2011/11/30 07:45:28 mussgill Exp $
0017 //
0018 //
0019
0020 // framework include files
0021 #include "FWCore/Framework/interface/Frameworkfwd.h"
0022 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0023 #include "FWCore/Framework/interface/EventSetup.h"
0024 #include "FWCore/Framework/interface/Event.h"
0025 #include "FWCore/Framework/interface/MakerMacros.h"
0026 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0027 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0028 #include "FWCore/ServiceRegistry/interface/Service.h"
0029 #include "FWCore/Utilities/interface/EDGetToken.h"
0030
0031 // user include files
0032 #include "Alignment/OfflineValidation/interface/EopVariables.h"
0033 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0034 #include "CommonTools/Utils/interface/TFileDirectory.h"
0035 #include "DataFormats/CaloTowers/interface/CaloTowerDefs.h"
0036 #include "DataFormats/Common/interface/Handle.h"
0037 #include "DataFormats/DTRecHit/interface/DTRecHitCollection.h"
0038 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0039 #include "DataFormats/HcalIsolatedTrack/interface/IsolatedPixelTrackCandidate.h"
0040 #include "DataFormats/HcalIsolatedTrack/interface/IsolatedPixelTrackCandidateFwd.h"
0041 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0042 #include "DataFormats/JetReco/interface/CaloJet.h"
0043 #include "DataFormats/RecoCandidate/interface/RecoCaloTowerCandidate.h"
0044 #include "DataFormats/TrackReco/interface/Track.h"
0045 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0046 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0047 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0048 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0049 #include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h"
0050 #include "TrackingTools/GeomPropagators/interface/PropagationExceptions.h"
0051 #include "TrackingTools/TrackAssociator/interface/TrackAssociatorParameters.h"
0052 #include "TrackingTools/TrackAssociator/interface/TrackDetectorAssociator.h"
0053
0054 // ROOT includes
0055 #include "TH1.h"
0056 #include "TTree.h"
0057
0058 //
0059 // class decleration
0060 //
0061
0062 class EopTreeWriter : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0063 public:
0064   explicit EopTreeWriter(const edm::ParameterSet&);
0065   ~EopTreeWriter() override = default;
0066
0067   static void fillDescriptions(edm::ConfigurationDescriptions&);
0068
0069 private:
0070   void analyze(const edm::Event&, const edm::EventSetup&) override;
0071   void endJob() override;
0072   double getDistInCM(double eta1, double phi1, double eta2, double phi2);
0073
0074   // ----------member data ---------------------------
0075   edm::InputTag src_;
0076   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> geometryToken_;
0077
0078   edm::Service<TFileService> fs_;
0079   TTree* tree_;
0080   EopVariables* treeMemPtr_;
0081   TrackDetectorAssociator trackAssociator_;
0082   TrackAssociatorParameters parameters_;
0083
0084   edm::EDGetTokenT<reco::TrackCollection> theTrackCollectionToken_;
0085   edm::EDGetTokenT<reco::IsolatedPixelTrackCandidateCollection> isoPixelTkToken_;
0086   edm::EDGetTokenT<EcalRecHitCollection> ecalRecHitTokenAlCaToken_;
0087   edm::EDGetTokenT<EcalRecHitCollection> ecalRecHitTokenEBRecoToken_;
0088   edm::EDGetTokenT<EcalRecHitCollection> ecalRecHitTokenEERecoToken_;
0089 };
0090
0091 //
0092 // constructors and destructor
0093 //
0094 EopTreeWriter::EopTreeWriter(const edm::ParameterSet& iConfig)
0095     : src_(iConfig.getParameter<edm::InputTag>("src")), geometryToken_(esConsumes()) {
0096   usesResource(TFileService::kSharedResource);
0097   //now do what ever initialization is needed
0098
0099   // TrackAssociator parameters
0100   edm::ParameterSet parameters = iConfig.getParameter<edm::ParameterSet>("TrackAssociatorParameters");
0101   edm::ConsumesCollector iC = consumesCollector();
0102   parameters_.loadParameters(parameters, iC);
0103
0104   tree_ = fs_->make<TTree>("EopTree", "EopTree");
0105   treeMemPtr_ = new EopVariables;
0106   tree_->Branch("EopVariables", &treeMemPtr_);  // address of pointer!
0107
0108   // do the consumes
0109   theTrackCollectionToken_ = consumes<reco::TrackCollection>(src_);
0110   isoPixelTkToken_ =
0111       consumes<reco::IsolatedPixelTrackCandidateCollection>(edm::InputTag("IsoProd", "HcalIsolatedTrackCollection"));
0112
0113   ecalRecHitTokenAlCaToken_ = consumes<EcalRecHitCollection>(edm::InputTag("IsoProd", "IsoTrackEcalRecHitCollection"));
0114   ecalRecHitTokenEBRecoToken_ = consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit", "EcalRecHitsEB"));
0115   ecalRecHitTokenEERecoToken_ = consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit", "EcalRecHitsEE"));
0116 }
0117
0118 //
0119 // member functions
0120 //
0121
0122 // ------------ method called to for each event  ------------
0123 void EopTreeWriter::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0124   using namespace edm;
0125
0126   // get geometry
0127   const CaloGeometry* geo = &iSetup.getData(geometryToken_);
0128
0129   // temporary collection of EB+EE recHits
0130   std::unique_ptr<EcalRecHitCollection> tmpEcalRecHitCollection(new EcalRecHitCollection);
0131   bool ecalInAlca = iEvent.getHandle(ecalRecHitTokenAlCaToken_).isValid();
0132   bool ecalInReco =
0133       iEvent.getHandle(ecalRecHitTokenEBRecoToken_) && iEvent.getHandle(ecalRecHitTokenEERecoToken_).isValid();
0134
0135   std::vector<edm::EDGetTokenT<EcalRecHitCollection> > ecalTokens_;
0136   if (ecalInAlca)
0137     ecalTokens_.push_back(ecalRecHitTokenAlCaToken_);
0138   else if (ecalInReco) {
0139     ecalTokens_.push_back(ecalRecHitTokenEBRecoToken_);
0140     ecalTokens_.push_back(ecalRecHitTokenEERecoToken_);
0141   } else {
0142     throw cms::Exception("MissingProduct", "can not find EcalRecHits");
0143   }
0144
0145   for (const auto& i : ecalTokens_) {
0146     edm::Handle<EcalRecHitCollection> ec = iEvent.getHandle(i);
0147     for (EcalRecHitCollection::const_iterator recHit = (*ec).begin(); recHit != (*ec).end(); ++recHit) {
0148       tmpEcalRecHitCollection->push_back(*recHit);
0149     }
0150   }
0151
0152   const auto& tracks = iEvent.get(theTrackCollectionToken_);
0153
0154   edm::Handle<reco::IsolatedPixelTrackCandidateCollection> isoPixelTracks = iEvent.getHandle(isoPixelTkToken_);
0155   bool pixelInAlca = isoPixelTracks.isValid();
0156
0157   double trackemc1;
0158   double trackemc3;
0159   double trackemc5;
0160   double trackhac1;
0161   double trackhac3;
0162   double trackhac5;
0163   double maxPNearby;
0164   double dist;
0165   double EnergyIn;
0166   double EnergyOut;
0167
0168   parameters_.useMuon = false;
0169
0170   if (pixelInAlca)
0171     if (isoPixelTracks->empty())
0172       return;
0173
0174   for (const auto& track : tracks) {
0175     bool noChargedTracks = true;
0176
0177     if (track.p() < 9.)
0178       continue;
0179
0180     trackAssociator_.useDefaultPropagator();
0181     TrackDetMatchInfo info = trackAssociator_.associate(
0182         iEvent,
0183         iSetup,
0184         trackAssociator_.getFreeTrajectoryState(&iSetup.getData(parameters_.bFieldToken), track),
0185         parameters_);
0186
0187     trackemc1 = info.nXnEnergy(TrackDetMatchInfo::EcalRecHits, 0);
0188     trackemc3 = info.nXnEnergy(TrackDetMatchInfo::EcalRecHits, 1);
0189     trackemc5 = info.nXnEnergy(TrackDetMatchInfo::EcalRecHits, 2);
0190     trackhac1 = info.nXnEnergy(TrackDetMatchInfo::HcalRecHits, 0);
0191     trackhac3 = info.nXnEnergy(TrackDetMatchInfo::HcalRecHits, 1);
0192     trackhac5 = info.nXnEnergy(TrackDetMatchInfo::HcalRecHits, 2);
0193
0194     if (trackhac3 < 5.)
0195       continue;
0196
0197     double etaecal = info.trkGlobPosAtEcal.eta();
0198     double phiecal = info.trkGlobPosAtEcal.phi();
0199
0200     maxPNearby = -10;
0201     dist = 50;
0202     for (const auto& track1 : tracks) {
0203       if (&track == &track1) {
0204         continue;
0205       }
0206       TrackDetMatchInfo info1 = trackAssociator_.associate(iEvent, iSetup, track1, parameters_);
0207       double etaecal1 = info1.trkGlobPosAtEcal.eta();
0208       double phiecal1 = info1.trkGlobPosAtEcal.phi();
0209
0210       if (etaecal1 == 0 && phiecal1 == 0)
0211         continue;
0212
0213       double ecDist = getDistInCM(etaecal, phiecal, etaecal1, phiecal1);
0214
0215       if (ecDist < 40.) {
0216         //calculate maximum P and sum P near seed track
0217         if (track1.p() > maxPNearby) {
0218           maxPNearby = track1.p();
0219           dist = ecDist;
0220         }
0221
0222         //apply loose isolation criteria
0223         if (track1.p() > 5.) {
0224           noChargedTracks = false;
0225           break;
0226         }
0227       }
0228     }
0229     EnergyIn = 0;
0230     EnergyOut = 0;
0231     if (noChargedTracks) {
0232       for (std::vector<EcalRecHit>::const_iterator ehit = tmpEcalRecHitCollection->begin();
0233            ehit != tmpEcalRecHitCollection->end();
0234            ehit++) {
0235         ////////////////////// FIND ECAL CLUSTER ENERGY
0236         // R-scheme of ECAL CLUSTERIZATION
0237         const GlobalPoint& posH = geo->getPosition((*ehit).detid());
0238         double phihit = posH.phi();
0239         double etahit = posH.eta();
0240
0241         double dHitCM = getDistInCM(etaecal, phiecal, etahit, phihit);
0242
0243         if (dHitCM < 9.0) {
0244           EnergyIn += ehit->energy();
0245         }
0246         if (dHitCM > 15.0 && dHitCM < 35.0) {
0247           EnergyOut += ehit->energy();
0248         }
0249       }
0250
0251       treeMemPtr_->fillVariables(track.charge(),
0252                                  track.innerOk(),
0253                                  track.outerRadius(),
0254                                  track.numberOfValidHits(),
0255                                  track.numberOfLostHits(),
0256                                  track.chi2(),
0257                                  track.normalizedChi2(),
0258                                  track.p(),
0259                                  track.pt(),
0260                                  track.ptError(),
0261                                  track.theta(),
0262                                  track.eta(),
0263                                  track.phi(),
0264                                  trackemc1,
0265                                  trackemc3,
0266                                  trackemc5,
0267                                  trackhac1,
0268                                  trackhac3,
0269                                  trackhac5,
0270                                  maxPNearby,
0271                                  dist,
0272                                  EnergyIn,
0273                                  EnergyOut);
0274
0275       tree_->Fill();
0276     }
0277   }
0278 }
0279
0280 // ------------ method called once each job just after ending the event loop  ------------
0281 void EopTreeWriter::endJob() {
0282   delete treeMemPtr_;
0283   treeMemPtr_ = nullptr;
0284 }
0285
0286 //*************************************************************
0287 double EopTreeWriter::getDistInCM(double eta1, double phi1, double eta2, double phi2) {
0288   //*************************************************************
0289
0290   static constexpr float EEBoundary = 1.479;  // psedo-rapidity
0291   static constexpr float EBRadius = 129;      // in cm
0292   static constexpr float EEIPdis = 317;       // in cm
0293
0294   double deltaPhi = phi1 - phi2;
0295   while (deltaPhi > M_PI)
0296     deltaPhi -= 2 * M_PI;
0297   while (deltaPhi <= -M_PI)
0298     deltaPhi += 2 * M_PI;
0299   double dR;
0300   // double Rec;
0301   double theta1 = 2 * atan(exp(-eta1));
0302   double theta2 = 2 * atan(exp(-eta2));
0303   double cotantheta1;
0304   if (cos(theta1) == 0)
0305     cotantheta1 = 0;
0306   else
0307     cotantheta1 = 1 / tan(theta1);
0308   double cotantheta2;
0309   if (cos(theta2) == 0)
0310     cotantheta2 = 0;
0311   else
0312     cotantheta2 = 1 / tan(theta2);
0313   // if (fabs(eta1)<1.479) Rec=129; //radius of ECAL barrel
0314   // else Rec=317; //distance from IP to ECAL endcap
0315   //|vect| times tg of acos(scalar product)
0316   // dR=fabs((Rec/sin(theta1))*tan(acos(sin(theta1)*sin(theta2)*(sin(phi1)*sin(phi2)+cos(phi1)*cos(phi2))+cos(theta1)*cos(theta2))));
0317
0318   if (fabs(eta1) < EEBoundary) {
0319     dR = EBRadius * sqrt((cotantheta1 - cotantheta2) * (cotantheta1 - cotantheta2) + deltaPhi * deltaPhi);
0320   } else {
0321     dR = EEIPdis *
0322          sqrt(tan(theta1) * tan(theta1) + tan(theta2) * tan(theta2) - 2 * tan(theta1) * tan(theta2) * cos(deltaPhi));
0323   }
0324   return dR;
0325 }
0326
0327 //*************************************************************
0328 void EopTreeWriter::fillDescriptions(edm::ConfigurationDescriptions& descriptions)
0329 //*************************************************************
0330 {
0331   edm::ParameterSetDescription desc;
0332   desc.setComment("Generate tree for Tracker Alignment E/p validation");
0333   desc.add<edm::InputTag>("src", edm::InputTag("generalTracks"));
0334
0335   // track association (take defaults)
0336   edm::ParameterSetDescription psd0;
0337   TrackAssociatorParameters::fillPSetDescription(psd0);
0338   desc.add<edm::ParameterSetDescription>("TrackAssociatorParameters", psd0);
0339
0340   descriptions.addWithDefaultLabel(desc);
0341 }
0342
0343 //define this as a plug-in
0344 DEFINE_FWK_MODULE(EopTreeWriter);