Skip to content

Commit

Permalink
modernize EopTreeWriter
Browse files Browse the repository at this point in the history
  • Loading branch information
mmusich committed Dec 9, 2021
1 parent a6366f4 commit d060680
Showing 1 changed file with 111 additions and 105 deletions.
216 changes: 111 additions & 105 deletions Alignment/OfflineValidation/plugins/EopTreeWriter.cc
Original file line number Diff line number Diff line change
Expand Up @@ -25,40 +25,35 @@
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "FWCore/ServiceRegistry/interface/Service.h"
#include "FWCore/Utilities/interface/EDGetToken.h"

// user include files
#include <DataFormats/TrackReco/interface/Track.h>
#include "DataFormats/TrackReco/interface/TrackFwd.h"
#include <TMath.h>
#include <TH1.h>
#include "TTree.h"

#include "Alignment/OfflineValidation/interface/EopVariables.h"
#include "CommonTools/UtilAlgos/interface/TFileService.h"
#include "CommonTools/Utils/interface/TFileDirectory.h"

#include "TrackingTools/TrackAssociator/interface/TrackDetectorAssociator.h"
#include "TrackingTools/TrackAssociator/interface/TrackAssociatorParameters.h"
#include "TrackingTools/GeomPropagators/interface/PropagationExceptions.h"
#include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h"

#include "DataFormats/CaloTowers/interface/CaloTowerDefs.h"
#include "DataFormats/Common/interface/Handle.h"
#include "DataFormats/DTRecHit/interface/DTRecHitCollection.h"
#include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
#include "DataFormats/HcalIsolatedTrack/interface/IsolatedPixelTrackCandidate.h"
#include "DataFormats/HcalIsolatedTrack/interface/IsolatedPixelTrackCandidateFwd.h"
#include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
#include "DataFormats/JetReco/interface/CaloJet.h"
#include "DataFormats/CaloTowers/interface/CaloTowerDefs.h"
#include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"

#include "DataFormats/RecoCandidate/interface/RecoCaloTowerCandidate.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/TrackReco/interface/TrackFwd.h"
#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
#include "Geometry/Records/interface/CaloGeometryRecord.h"
#include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
#include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h"
#include "TrackingTools/GeomPropagators/interface/PropagationExceptions.h"
#include "TrackingTools/TrackAssociator/interface/TrackAssociatorParameters.h"
#include "TrackingTools/TrackAssociator/interface/TrackDetectorAssociator.h"

#include "DataFormats/HcalIsolatedTrack/interface/IsolatedPixelTrackCandidateFwd.h"
#include "DataFormats/HcalIsolatedTrack/interface/IsolatedPixelTrackCandidate.h"
#include "DataFormats/RecoCandidate/interface/RecoCaloTowerCandidate.h"
#include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
#include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
#include "DataFormats/DTRecHit/interface/DTRecHitCollection.h"

#include "Alignment/OfflineValidation/interface/EopVariables.h"
// ROOT includes
#include "TH1.h"
#include "TTree.h"

//
// class decleration
Expand All @@ -67,13 +62,13 @@
class EopTreeWriter : public edm::one::EDAnalyzer<edm::one::SharedResources> {
public:
explicit EopTreeWriter(const edm::ParameterSet&);
~EopTreeWriter() override;
~EopTreeWriter() override = default;

static void fillDescriptions(edm::ConfigurationDescriptions&);

private:
void beginJob() override;
void analyze(const edm::Event&, const edm::EventSetup&) override;
void endJob() override;

double getDistInCM(double eta1, double phi1, double eta2, double phi2);

// ----------member data ---------------------------
Expand All @@ -85,15 +80,13 @@ class EopTreeWriter : public edm::one::EDAnalyzer<edm::one::SharedResources> {
EopVariables* treeMemPtr_;
TrackDetectorAssociator trackAssociator_;
TrackAssociatorParameters parameters_;
};

//
// constants, enums and typedefs
//

//
// static data member definitions
//
edm::EDGetTokenT<reco::TrackCollection> theTrackCollectionToken_;
edm::EDGetTokenT<reco::IsolatedPixelTrackCandidateCollection> isoPixelTkToken_;
edm::EDGetTokenT<EcalRecHitCollection> ecalRecHitTokenAlCaToken_;
edm::EDGetTokenT<EcalRecHitCollection> ecalRecHitTokenEBRecoToken_;
edm::EDGetTokenT<EcalRecHitCollection> ecalRecHitTokenEERecoToken_;
};

//
// constructors and destructor
Expand All @@ -111,11 +104,15 @@ EopTreeWriter::EopTreeWriter(const edm::ParameterSet& iConfig)
tree_ = fs_->make<TTree>("EopTree", "EopTree");
treeMemPtr_ = new EopVariables;
tree_->Branch("EopVariables", &treeMemPtr_); // address of pointer!
}

EopTreeWriter::~EopTreeWriter() {
// do anything here that needs to be done at destruction time
// (e.g. close files, deallocate resources etc.)
// do the consumes
theTrackCollectionToken_ = consumes<reco::TrackCollection>(src_);
isoPixelTkToken_ =
consumes<reco::IsolatedPixelTrackCandidateCollection>(edm::InputTag("IsoProd", "HcalIsolatedTrackCollection"));

ecalRecHitTokenAlCaToken_ = consumes<EcalRecHitCollection>(edm::InputTag("IsoProd", "IsoTrackEcalRecHitCollection"));
ecalRecHitTokenEBRecoToken_ = consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit", "EcalRecHitsEB"));
ecalRecHitTokenEERecoToken_ = consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit", "EcalRecHitsEE"));
}

//
Expand All @@ -128,71 +125,63 @@ void EopTreeWriter::analyze(const edm::Event& iEvent, const edm::EventSetup& iSe

// get geometry
const CaloGeometry* geo = &iSetup.getData(geometryToken_);
// const CaloSubdetectorGeometry* towerGeometry =
// geo->getSubdetectorGeometry(DetId::Calo, CaloTowerDetId::SubdetId);

// temporary collection of EB+EE recHits
std::unique_ptr<EcalRecHitCollection> tmpEcalRecHitCollection(new EcalRecHitCollection);
std::vector<edm::InputTag> ecalLabels_;
bool ecalInAlca = iEvent.getHandle(ecalRecHitTokenAlCaToken_).isValid();
bool ecalInReco =
iEvent.getHandle(ecalRecHitTokenEBRecoToken_) && iEvent.getHandle(ecalRecHitTokenEERecoToken_).isValid();

edm::Handle<EcalRecHitCollection> tmpEc;
bool ecalInAlca = iEvent.getByLabel(edm::InputTag("IsoProd", "IsoTrackEcalRecHitCollection"), tmpEc);
bool ecalInReco = iEvent.getByLabel(edm::InputTag("ecalRecHit", "EcalRecHitsEB"), tmpEc) &&
iEvent.getByLabel(edm::InputTag("ecalRecHit", "EcalRecHitsEE"), tmpEc);
std::vector<edm::EDGetTokenT<EcalRecHitCollection> > ecalTokens_;
if (ecalInAlca)
ecalLabels_.push_back(edm::InputTag("IsoProd", "IsoTrackEcalRecHitCollection"));
ecalTokens_.push_back(ecalRecHitTokenAlCaToken_);
else if (ecalInReco) {
ecalLabels_.push_back(edm::InputTag("ecalRecHit", "EcalRecHitsEB"));
ecalLabels_.push_back(edm::InputTag("ecalRecHit", "EcalRecHitsEE"));
} else
ecalTokens_.push_back(ecalRecHitTokenEBRecoToken_);
ecalTokens_.push_back(ecalRecHitTokenEERecoToken_);
} else {
throw cms::Exception("MissingProduct", "can not find EcalRecHits");
}

std::vector<edm::InputTag>::const_iterator i;
for (i = ecalLabels_.begin(); i != ecalLabels_.end(); i++) {
edm::Handle<EcalRecHitCollection> ec;
iEvent.getByLabel(*i, ec);
for (const auto& i : ecalTokens_) {
edm::Handle<EcalRecHitCollection> ec = iEvent.getHandle(i);
for (EcalRecHitCollection::const_iterator recHit = (*ec).begin(); recHit != (*ec).end(); ++recHit) {
tmpEcalRecHitCollection->push_back(*recHit);
}
}

edm::Handle<reco::TrackCollection> tracks;
iEvent.getByLabel(src_, tracks);
const auto& tracks = iEvent.get(theTrackCollectionToken_);

edm::Handle<reco::IsolatedPixelTrackCandidateCollection> isoPixelTracks;
edm::Handle<reco::IsolatedPixelTrackCandidateCollection> tmpPix;
bool pixelInAlca = iEvent.getByLabel(edm::InputTag("IsoProd", "HcalIsolatedTrackCollection"), tmpPix);
if (pixelInAlca)
iEvent.getByLabel(edm::InputTag("IsoProd", "HcalIsolatedTrackCollection"), isoPixelTracks);

Double_t trackemc1;
Double_t trackemc3;
Double_t trackemc5;
Double_t trackhac1;
Double_t trackhac3;
Double_t trackhac5;
Double_t maxPNearby;
Double_t dist;
Double_t EnergyIn;
Double_t EnergyOut;
edm::Handle<reco::IsolatedPixelTrackCandidateCollection> isoPixelTracks = iEvent.getHandle(isoPixelTkToken_);
bool pixelInAlca = isoPixelTracks.isValid();

double trackemc1;
double trackemc3;
double trackemc5;
double trackhac1;
double trackhac3;
double trackhac5;
double maxPNearby;
double dist;
double EnergyIn;
double EnergyOut;

parameters_.useMuon = false;

if (pixelInAlca)
if (isoPixelTracks->empty())
return;

for (reco::TrackCollection::const_iterator track = tracks->begin(); track != tracks->end(); ++track) {
for (const auto& track : tracks) {
bool noChargedTracks = true;

if (track->p() < 9.)
if (track.p() < 9.)
continue;

trackAssociator_.useDefaultPropagator();
TrackDetMatchInfo info = trackAssociator_.associate(
iEvent,
iSetup,
trackAssociator_.getFreeTrajectoryState(&iSetup.getData(parameters_.bFieldToken), *track),
trackAssociator_.getFreeTrajectoryState(&iSetup.getData(parameters_.bFieldToken), track),
parameters_);

trackemc1 = info.nXnEnergy(TrackDetMatchInfo::EcalRecHits, 0);
Expand All @@ -210,10 +199,11 @@ void EopTreeWriter::analyze(const edm::Event& iEvent, const edm::EventSetup& iSe

maxPNearby = -10;
dist = 50;
for (reco::TrackCollection::const_iterator track1 = tracks->begin(); track1 != tracks->end(); track1++) {
if (track == track1)
for (const auto& track1 : tracks) {
if (&track == &track1) {
continue;
TrackDetMatchInfo info1 = trackAssociator_.associate(iEvent, iSetup, *track1, parameters_);
}
TrackDetMatchInfo info1 = trackAssociator_.associate(iEvent, iSetup, track1, parameters_);
double etaecal1 = info1.trkGlobPosAtEcal.eta();
double phiecal1 = info1.trkGlobPosAtEcal.phi();

Expand All @@ -224,13 +214,13 @@ void EopTreeWriter::analyze(const edm::Event& iEvent, const edm::EventSetup& iSe

if (ecDist < 40.) {
//calculate maximum P and sum P near seed track
if (track1->p() > maxPNearby) {
maxPNearby = track1->p();
if (track1.p() > maxPNearby) {
maxPNearby = track1.p();
dist = ecDist;
}

//apply loose isolation criteria
if (track1->p() > 5.) {
if (track1.p() > 5.) {
noChargedTracks = false;
break;
}
Expand Down Expand Up @@ -258,19 +248,19 @@ void EopTreeWriter::analyze(const edm::Event& iEvent, const edm::EventSetup& iSe
}
}

treeMemPtr_->fillVariables(track->charge(),
track->innerOk(),
track->outerRadius(),
track->numberOfValidHits(),
track->numberOfLostHits(),
track->chi2(),
track->normalizedChi2(),
track->p(),
track->pt(),
track->ptError(),
track->theta(),
track->eta(),
track->phi(),
treeMemPtr_->fillVariables(track.charge(),
track.innerOk(),
track.outerRadius(),
track.numberOfValidHits(),
track.numberOfLostHits(),
track.chi2(),
track.normalizedChi2(),
track.p(),
track.pt(),
track.ptError(),
track.theta(),
track.eta(),
track.phi(),
trackemc1,
trackemc3,
trackemc5,
Expand All @@ -287,21 +277,25 @@ void EopTreeWriter::analyze(const edm::Event& iEvent, const edm::EventSetup& iSe
}
}

// ------------ method called once each job just before starting event loop ------------
void EopTreeWriter::beginJob() {}

// ------------ method called once each job just after ending the event loop ------------
void EopTreeWriter::endJob() {
delete treeMemPtr_;
treeMemPtr_ = nullptr;
}

//*************************************************************
double EopTreeWriter::getDistInCM(double eta1, double phi1, double eta2, double phi2) {
//*************************************************************

static constexpr float EEBoundary = 1.479; // psedo-rapidity
static constexpr float EBRadius = 129; // in cm
static constexpr float EEIPdis = 317; // in cm

double deltaPhi = phi1 - phi2;
while (deltaPhi > TMath::Pi())
deltaPhi -= 2 * TMath::Pi();
while (deltaPhi <= -TMath::Pi())
deltaPhi += 2 * TMath::Pi();
while (deltaPhi > M_PI)
deltaPhi -= 2 * M_PI;
while (deltaPhi <= -M_PI)
deltaPhi += 2 * M_PI;
double dR;
// double Rec;
double theta1 = 2 * atan(exp(-eta1));
Expand All @@ -320,13 +314,25 @@ double EopTreeWriter::getDistInCM(double eta1, double phi1, double eta2, double
// else Rec=317; //distance from IP to ECAL endcap
//|vect| times tg of acos(scalar product)
// dR=fabs((Rec/sin(theta1))*tan(acos(sin(theta1)*sin(theta2)*(sin(phi1)*sin(phi2)+cos(phi1)*cos(phi2))+cos(theta1)*cos(theta2))));
if (fabs(eta1) < 1.479)
dR = 129 * sqrt((cotantheta1 - cotantheta2) * (cotantheta1 - cotantheta2) + deltaPhi * deltaPhi);
else
dR = 317 *

if (fabs(eta1) < EEBoundary) {
dR = EBRadius * sqrt((cotantheta1 - cotantheta2) * (cotantheta1 - cotantheta2) + deltaPhi * deltaPhi);
} else {
dR = EEIPdis *
sqrt(tan(theta1) * tan(theta1) + tan(theta2) * tan(theta2) - 2 * tan(theta1) * tan(theta2) * cos(deltaPhi));
}
return dR;
}

//*************************************************************
void EopTreeWriter::fillDescriptions(edm::ConfigurationDescriptions& descriptions)
//*************************************************************
{
edm::ParameterSetDescription desc;
desc.setComment("Generate tree for Tracker Alignment E/p validation");
desc.add<edm::InputTag>("src", edm::InputTag("generalTracks"));
descriptions.addWithDefaultLabel(desc);
}

//define this as a plug-in
DEFINE_FWK_MODULE(EopTreeWriter);

0 comments on commit d060680

Please sign in to comment.