EDAnalyzer to exercise and demonstrate usage of Calibratable tree. More...

#include <CalibratableTest.h>

Inheritance diagram for CalibratableTest:

Public Member Functions
	CalibratableTest (const edm::ParameterSet &)
template<class T >
void	getCollection (edm::Handle< T > &c, const edm::InputTag &tag, const edm::Event &event) const
	~CalibratableTest ()
Private Member Functions
virtual void	analyze (const edm::Event &, const edm::EventSetup &)
virtual void	beginJob ()
double	deltaR (const double &eta1, const double &eta2, const double &phi1, const double &phi2)
virtual void	endJob ()
void	extractCandidate (const reco::PFCandidate &cand)
virtual void	fillTreeAndReset ()
std::vector< unsigned >	findCandidatesInDeltaR (const reco::PFSimParticle &pft, const std::vector< reco::PFCandidate > &cands, const double &deltaR)
std::vector< unsigned >	findPrimarySimParticles (const std::vector< reco::PFSimParticle > &sims)
Private Attributes
pftools::Calibratable *	calib_
edm::Handle < reco::PFClusterCollection > *	clustersEcal_
edm::Handle < reco::PFClusterCollection > *	clustersHcal_
int	debug_
double	deltaRCandToSim_
edm::Service< TFileService >	fileservice_
edm::InputTag	inputTagCandidates_
edm::InputTag	inputTagClustersEcal_
edm::InputTag	inputTagClustersHcal_
edm::InputTag	inputTagSimParticles_
unsigned	nEventFails_
unsigned	nEventWrites_
unsigned	nParticleFails_
unsigned	nParticleWrites_
edm::Handle < reco::PFCandidateCollection > *	pfCandidates_
edm::Handle < reco::PFSimParticleCollection > *	simParticles_
bool	thisEventPasses_
bool	thisParticlePasses_
TTree *	tree_

Detailed Description

EDAnalyzer to exercise and demonstrate usage of Calibratable tree.

Author:: Jamie Ballin, Imperial College London

Date:: November 2008

OBJECTIVE: this analyzer will create a tree of PFClusterTools/Calibratable instances. Each entry of the tree represents particle flow information relating to one pion in a multi-pion event. The use of monte carlo is assumed as sim particle information is extracted too.

USAGE: This is an analyzer, not a producer: the tree it produces is stored in a seperate ROOT file using the CMSSW TFileService.

Consult test/CalibratableTest_FastSim.py for a sample CMSSW configuration. Run the example from that directory with, cmsRun CalibratableTest_FastSim.py This will produce a file with a tree inside. Browse the contents in bare root with a TBrowser.

NOTE: This analyzer does not exercise the complete functionality of the Calibratable class. For instance, tracks and rechits are not stored. Contact me (the author) for more details, and for an analyzer that we use to extract particle information to calibrate the particle flow algorithm. A MUCH more involved example from my private analysis is to be found in, UserCode/JamieBallin/interface/DipionDelegate.h

General details about the usage of Calibratable may be found at, https://twiki.cern.ch/twiki/bin/view/CMS/PFClusterToolsPackage

Definition at line 57 of file CalibratableTest.h.

Constructor & Destructor Documentation

CalibratableTest::CalibratableTest ( const edm::ParameterSet & parameters ) [explicit]

Definition at line 28 of file CalibratableTest.cc.

References calib_, gather_cfg::cout, debug_, deltaRCandToSim_, fileservice_, edm::ParameterSet::getParameter(), inputTagCandidates_, inputTagClustersEcal_, inputTagClustersHcal_, inputTagSimParticles_, and tree_.

                                                                    :
        debug_(0), nParticleWrites_(0), nParticleFails_(0), nEventWrites_(0), nEventFails_(0),
                         deltaRCandToSim_(0.4) {

        std::cout << __PRETTY_FUNCTION__ << std::endl;
        
        /* This procedure is GENERIC to storing any dictionary enable class in a ROOT tree. */
        tree_ = fileservice_->make<TTree>("CalibratableTest", "");
        calib_ = new Calibratable();
        tree_->Branch("Calibratable", "pftools::Calibratable", &calib_, 32000, 2);
        
        inputTagCandidates_= parameters.getParameter<InputTag>("PFCandidates");
        inputTagSimParticles_= parameters.getParameter<InputTag>("PFSimParticles");
        inputTagClustersEcal_= parameters.getParameter<InputTag>("PFClustersEcal");
        inputTagClustersHcal_= parameters.getParameter<InputTag>("PFClustersHcal");
        deltaRCandToSim_ = parameters.getParameter<double>("deltaRCandToSim");
        debug_= parameters.getParameter<int>("debug");
}

CalibratableTest::~CalibratableTest ( )

Definition at line 47 of file CalibratableTest.cc.

                                    {

}

Member Function Documentation

void CalibratableTest::analyze	(	const edm::Event &	event,
		const edm::EventSetup &	iSetup
	)		`[private, virtual]`

Implements edm::EDAnalyzer.

Definition at line 57 of file CalibratableTest.cc.

                                             {
        if (debug_ > 1)
                std::cout << __PRETTY_FUNCTION__ << "\n";

        //Extract new collection references
        pfCandidates_ = new Handle<PFCandidateCollection>;
        simParticles_ = new Handle<PFSimParticleCollection>;
        clustersEcal_ = new Handle<PFClusterCollection>;
        clustersHcal_ = new Handle<PFClusterCollection>;

        getCollection(*pfCandidates_, inputTagCandidates_, event);
        getCollection(*simParticles_, inputTagSimParticles_, event);
        getCollection(*clustersEcal_, inputTagClustersEcal_, event);
        getCollection(*clustersHcal_, inputTagClustersHcal_, event);

        //Reset calibratable branch
        thisEventPasses_ = true;
        thisParticlePasses_ = true;
        calib_->reset();

        if (debug_ > 1)
                std::cout << "\tStarting event..."<< std::endl;

        PFSimParticleCollection sims = **simParticles_;
        PFCandidateCollection candidates = **pfCandidates_;
        PFClusterCollection clustersEcal = **clustersEcal_;
        PFClusterCollection clustersHcal = **clustersHcal_;

        if (sims.size() == 0) {
                std::cout << "\tNo sim particles found!" << std::endl;
                thisEventPasses_ = false;
        }

        //Find primary pions in the event
        std::vector<unsigned> primarySims = findPrimarySimParticles(sims);
        if (debug_) {
                std::cout << "\tFound "<< primarySims.size()
                                << " primary sim particles, "<< (**pfCandidates_).size() << " pfCandidates\n";
        }
        for (std::vector<unsigned>::const_iterator cit = primarySims.begin(); cit
                        != primarySims.end(); ++cit) {
                //There will be one write to the tree for each pion found.
                if (debug_ > 1)
                        std::cout << "\t**Starting particle...**\n";
                const PFSimParticle& sim = sims[*cit];
                //One sim per calib =>
                calib_->sim_numEvent_ = 1;
                calib_->sim_isMC_ = true;
                calib_->sim_energyEvent_ = sim.trajectoryPoint(PFTrajectoryPoint::ClosestApproach).momentum().E();
                calib_->sim_phi_ = sim.trajectoryPoint(PFTrajectoryPoint::ClosestApproach).momentum().Phi();
                calib_->sim_eta_ = sim.trajectoryPoint(PFTrajectoryPoint::ClosestApproach).momentum().Eta();

                if (sim.nTrajectoryPoints() > PFTrajectoryPoint::ECALEntrance) {
                        calib_->sim_etaEcal_ = sim.trajectoryPoint(PFTrajectoryPoint::ECALEntrance).positionREP().Eta();
                        calib_->sim_phiEcal_ = sim.trajectoryPoint(PFTrajectoryPoint::ECALEntrance).positionREP().Phi();
                }
                if (sim.nTrajectoryPoints() > PFTrajectoryPoint::HCALEntrance) {
                        calib_->sim_etaHcal_ = sim.trajectoryPoint(PFTrajectoryPoint::HCALEntrance).positionREP().Eta();
                        calib_->sim_phiHcal_ = sim.trajectoryPoint(PFTrajectoryPoint::HCALEntrance).positionREP().Phi();
                }

                // Find candidates near this sim particle
                std::vector<unsigned> matchingCands = findCandidatesInDeltaR(sim,
                                candidates, deltaRCandToSim_);
                if (debug_ > 3)
                        std::cout << "\t\tFound candidates near sim, found "
                                        << matchingCands.size()<< " of them.\n";
                if (matchingCands.size() == 0)
                        thisParticlePasses_ = false;
                for (std::vector<unsigned>::const_iterator mcIt = matchingCands.begin(); mcIt
                                != matchingCands.end(); ++mcIt) {
                        const PFCandidate& theCand = candidates[*mcIt];
                        extractCandidate(theCand);
                }
                //Finally,
                fillTreeAndReset();
                
        }

        delete pfCandidates_;
        delete simParticles_;
        delete clustersEcal_;
        delete clustersHcal_;

        if (thisEventPasses_)
                ++nEventWrites_;
        else
                ++nEventFails_;

}

void CalibratableTest::beginJob ( void ) [private, virtual]

Reimplemented from edm::EDAnalyzer.

Definition at line 51 of file CalibratableTest.cc.

References gather_cfg::cout, and debug_.

                                {
        if (debug_ > 1)
                std::cout << __PRETTY_FUNCTION__ << "\n";

}

double CalibratableTest::deltaR	(	const double &	eta1,
		const double &	eta2,
		const double &	phi1,
		const double &	phi2
	)		`[private]`

Definition at line 305 of file CalibratableTest.cc.

References Geom::deltaPhi(), M_PI, funct::pow(), and mathSSE::sqrt().

Referenced by findCandidatesInDeltaR().

                                                        {
        double deltaEta = fabs(eta1 - eta2);
        double deltaPhi = fabs(phi1 - phi2);
        if (deltaPhi > M_PI) {
                deltaPhi = 2 * M_PI- deltaPhi;
        }
        return sqrt(pow(deltaEta, 2) + pow(deltaPhi, 2));
}

void CalibratableTest::endJob ( void ) [private, virtual]

Reimplemented from edm::EDAnalyzer.

Definition at line 291 of file CalibratableTest.cc.

References gather_cfg::cout, debug_, nEventFails_, nEventWrites_, nParticleFails_, and nParticleWrites_.

                              {

        if (debug_> 0) {
                std::cout << __PRETTY_FUNCTION__ << std::endl;

                std::cout << "\tnParticleWrites: "<< nParticleWrites_
                                << ", nParticleFails: "<< nParticleFails_ << "\n";
                std::cout << "\tnEventWrites: "<< nEventWrites_ << ", nEventFails: "
                                << nEventFails_ << "\n";
                std::cout << "Leaving "<< __PRETTY_FUNCTION__ << "\n";
        }

}

void CalibratableTest::extractCandidate ( const reco::PFCandidate & cand ) [private]

Definition at line 171 of file CalibratableTest.cc.

References Association::block, calib_, pftools::Calibratable::cands_, pftools::Calibratable::cluster_ecal_, pftools::Calibratable::cluster_hcal_, gather_cfg::cout, debug_, ECAL, reco::PFBlock::elements(), asciidump::elements, reco::PFCluster::energy(), pftools::CandidateWrapper::energy_, pftools::CandidateWrapper::energyEcal_, pftools::CandidateWrapper::energyHcal_, pftools::CandidateWrapper::eta_, HCAL, reco::PFCluster::layer(), pftools::CandidateWrapper::phi_, reco::PFCluster::positionREP(), and pftools::CandidateWrapper::type_.

Referenced by analyze().

                                                               {
        if (debug_ > 3)
                std::cout << "\tCandidate: "<< cand << "\n";

        //There may be several PFCandiates per sim particle. So we create a mini-class
        //to represent each one. CandidateWrapper is defined in Calibratable.
        //It's very easy to use, as we shall see...
        CandidateWrapper cw;
        cw.energy_ = cand.energy();
        cw.eta_ = cand.eta();
        cw.phi_ = cand.phi();
        cw.type_ = cand.particleId();
        cw.energyEcal_ = cand.ecalEnergy();
        cw.energyHcal_ = cand.hcalEnergy();
        if (debug_ > 4)
                std::cout << "\t\tECAL energy = " << cand.ecalEnergy()
                                << ", HCAL energy = " << cand.hcalEnergy() << "\n";

        //Now, extract block elements from the pfCandidate:
        PFCandidate::ElementsInBlocks eleInBlocks = cand.elementsInBlocks();
        if (debug_ > 3)
                std::cout << "\tLooping over elements in blocks, "
                                << eleInBlocks.size() << " of them."<< std::endl;
        for (PFCandidate::ElementsInBlocks::iterator bit = eleInBlocks.begin(); bit
                        != eleInBlocks.end(); ++bit) {

                /* 
                 * Find PFClusters associated with this candidate.
                 */
                
                //Extract block reference
                PFBlockRef blockRef((*bit).first);
                //Extract index
                unsigned indexInBlock((*bit).second);
                //Dereference the block (what a palava!)
                const PFBlock& block = *blockRef;
                //And finally get a handle on the elements
                const edm::OwnVector<reco::PFBlockElement> & elements = block.elements();
                //get references to the candidate's track, ecal clusters and hcal clusters
                switch (elements[indexInBlock].type()) {
                case PFBlockElement::ECAL: {
                        reco::PFClusterRef clusterRef = elements[indexInBlock].clusterRef();
                        const PFCluster theRealCluster = *clusterRef;
                        CalibratableElement d(theRealCluster.energy(),
                                        theRealCluster.positionREP().eta(), theRealCluster.positionREP().phi(), theRealCluster.layer() );
                        calib_->cluster_ecal_.push_back(d);
                        if (debug_ > 4)
                                std::cout << "\t\tECAL cluster: "<< theRealCluster << "\n";

                        break;
                }

                case PFBlockElement::HCAL: {
                        reco::PFClusterRef clusterRef = elements[indexInBlock].clusterRef();
                        const PFCluster theRealCluster = *clusterRef;
                        CalibratableElement d(theRealCluster.energy(),
                                        theRealCluster.positionREP().eta(), theRealCluster.positionREP().phi(), theRealCluster.layer() );
                        calib_->cluster_hcal_.push_back(d);
                        if (debug_ > 4)
                                std::cout << "\t\tHCAL cluster: "<< theRealCluster << "\n";

                        break;
                }

                default:
                        if (debug_ > 3)
                                std::cout << "\t\tOther block type: "
                                                << elements[indexInBlock].type() << "\n";
                        break;
                }

        }
        //For each candidate found,
        calib_->cands_.push_back(cw);
}

void CalibratableTest::fillTreeAndReset ( ) [private, virtual]

Definition at line 274 of file CalibratableTest.cc.

References calib_, gather_cfg::cout, debug_, nParticleFails_, nParticleWrites_, pftools::Calibratable::recompute(), pftools::Calibratable::reset(), thisEventPasses_, thisParticlePasses_, and tree_.

Referenced by analyze().

                                        {
        if (thisEventPasses_ && thisParticlePasses_) {
                ++nParticleWrites_;
                calib_->recompute();
                if (debug_> 4) {
                        //print a summary
                        std::cout << *calib_;
                }
                tree_->Fill();
        } else {
                ++nParticleFails_;
        }
        if (debug_ > 1)
                std::cout << "\t**Finished particle.**\n";
        calib_->reset();
}

std::vector< unsigned > CalibratableTest::findCandidatesInDeltaR	(	const reco::PFSimParticle &	pft,
		const std::vector< reco::PFCandidate > &	cands,
		const double &	deltaR
	)		`[private]`

Definition at line 247 of file CalibratableTest.cc.

References deltaR(), reco::LeafCandidate::eta(), getHLTprescales::index, reco::LeafCandidate::phi(), reco::PFTrajectoryPoint::positionREP(), and reco::PFTrack::trajectoryPoint().

Referenced by analyze().

                                         {

        unsigned index(0);
        std::vector<unsigned> answers;

        double trEta = pft.trajectoryPoint(PFTrajectoryPoint::ECALEntrance).positionREP().Eta();
        double trPhi = pft.trajectoryPoint(PFTrajectoryPoint::ECALEntrance).positionREP().Phi();

        for (std::vector<PFCandidate>::const_iterator cit = cands.begin(); cit
                        != cands.end(); ++cit) {

                PFCandidate cand = *cit;
                double cEta = cand.eta();
                double cPhi = cand.phi();

                if (deltaR(cEta, trEta, cPhi, trPhi) < deltaRCut) {
                        //accept
                        answers.push_back(index);
                }

                ++index;
        }
        return answers;
}

std::vector< unsigned > CalibratableTest::findPrimarySimParticles ( const std::vector< reco::PFSimParticle > & sims ) [private]

Definition at line 149 of file CalibratableTest.cc.

References abs, reco::PFSimParticle::daughterIds(), getHLTprescales::index, reco::PFSimParticle::motherId(), and reco::PFSimParticle::pdgCode().

Referenced by analyze().

                                                      {
        std::vector<unsigned> answers;
        unsigned index(0);
        for (std::vector<PFSimParticle>::const_iterator cit = sims.begin(); cit
                        != sims.end(); ++cit) {
                PFSimParticle theSim = *cit;
                //TODO: what about rejected events?
                if (theSim.motherId() >= 0)
                        continue;
                int particleId = abs(theSim.pdgCode());
                if (particleId != 211)
                        continue;
                //TODO: ...particularly interacting pions?
                if (theSim.daughterIds().size() > 0)
                        continue;
                answers.push_back(index);
                ++index;
        }
        return answers;
}

template<class T >

void CalibratableTest::getCollection	(	edm::Handle< T > &	c,
		const edm::InputTag &	tag,
		const edm::Event &	event
	)		const

Definition at line 147 of file CalibratableTest.h.

References gather_cfg::cout, and edm::HandleBase::isValid().

Referenced by analyze().

                                                                     {

        try {
                event.getByLabel(tag, c);
                if(!c.isValid()) {
                        std::cout << "Warning! Collection for label " << tag << " is not valid!" << std::endl;
                }
        }
        catch (cms::Exception& err) {
                std::cout << "Couldn't get collection\n";
                std::ostringstream err;
                //LogError("Error getting collection!") << err;
        }
}