Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-12-13 13:13:41

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
0336   edm::ParameterSetDescription psd0;
0337   psd0.add<bool>("accountForTrajectoryChangeCalo", false);
0338   psd0.add<bool>("propagateAllDirections", true);
0339   psd0.add<bool>("truthMatch", false);
0340   psd0.add<bool>("useCalo", false);
0341   psd0.add<bool>("useEcal", true);
0342   psd0.add<bool>("useGEM", false);
0343   psd0.add<bool>("useHO", true);
0344   psd0.add<bool>("useHcal", true);
0345   psd0.add<bool>("useME0", false);
0346   psd0.add<bool>("useMuon", true);
0347   psd0.add<bool>("usePreshower", false);
0348   psd0.add<double>("dREcal", 9999.0);
0349   psd0.add<double>("dREcalPreselection", 0.05);
0350   psd0.add<double>("dRHcal", 9999.0);
0351   psd0.add<double>("dRHcalPreselection", 0.2);
0352   psd0.add<double>("dRMuon", 9999.0);
0353   psd0.add<double>("dRMuonPreselection", 0.2);
0354   psd0.add<double>("dRPreshowerPreselection", 0.2);
0355   psd0.add<double>("muonMaxDistanceSigmaX", 0.0);
0356   psd0.add<double>("muonMaxDistanceSigmaY", 0.0);
0357   psd0.add<double>("muonMaxDistanceX", 5.0);
0358   psd0.add<double>("muonMaxDistanceY", 5.0);
0359   psd0.add<double>("trajectoryUncertaintyTolerance", -1.0);
0360   psd0.add<edm::InputTag>("CSCSegmentCollectionLabel", edm::InputTag("cscSegments"));
0361   psd0.add<edm::InputTag>("CaloTowerCollectionLabel", edm::InputTag("towerMaker"));
0362   psd0.add<edm::InputTag>("DTRecSegment4DCollectionLabel", edm::InputTag("dt4DSegments"));
0363   psd0.add<edm::InputTag>("EBRecHitCollectionLabel", edm::InputTag("IsoProd", "IsoTrackEcalRecHitCollection"));
0364   psd0.add<edm::InputTag>("EERecHitCollectionLabel", edm::InputTag("IsoProd", "IsoTrackEcalRecHitCollection"));
0365   psd0.add<edm::InputTag>("GEMSegmentCollectionLabel", edm::InputTag("gemSegments"));
0366   psd0.add<edm::InputTag>("HBHERecHitCollectionLabel", edm::InputTag("IsoProd", "IsoTrackHBHERecHitCollection"));
0367   psd0.add<edm::InputTag>("HORecHitCollectionLabel", edm::InputTag("IsoProd", "IsoTrackHORecHitCollection"));
0368   psd0.add<edm::InputTag>("ME0SegmentCollectionLabel", edm::InputTag("me0Segments"));
0369   desc.add<edm::ParameterSetDescription>("TrackAssociatorParameters", psd0);
0370 
0371   descriptions.addWithDefaultLabel(desc);
0372 }
0373 
0374 //define this as a plug-in
0375 DEFINE_FWK_MODULE(EopTreeWriter);