Inheritance diagram for LHEAnalyzer:

Public Member Functions
	LHEAnalyzer (const edm::ParameterSet &params)
virtual	~LHEAnalyzer ()
Protected Member Functions
virtual void	analyze (const edm::Event &event, const edm::EventSetup &es)
Private Member Functions
void	fillDJRSched (unsigned int min, unsigned int max)
Private Attributes
unsigned int	binsDelta
unsigned int	binsPt
double	defaultDeltaCut
double	defaultPtCut
boost::ptr_vector< JetClustering >	deltaClustering
std::vector< unsigned int >	djrSched
std::vector< TH1 * >	histoDelta
std::vector< TH1 * >	histoPt
JetInput	jetInput
double	maxDelta
unsigned int	maxDJR
double	maxEta
double	maxPt
double	minDelta
unsigned int	minDJR
double	minPt
std::auto_ptr< JetClustering >	ptClustering
double	ptFraction
edm::InputTag	sourceLabel
bool	useEt

Detailed Description

Definition at line 38 of file LHEAnalyzer.cc.

Constructor & Destructor Documentation

LHEAnalyzer::LHEAnalyzer ( const edm::ParameterSet & params ) [explicit]

Definition at line 77 of file LHEAnalyzer.cc.

References edm::ParameterSet::addParameter(), binsDelta, binsPt, defaultDeltaCut, defaultPtCut, deltaClustering, fillDJRSched(), edm::ParameterSet::getParameter(), h, histoDelta, histoPt, i, maxDelta, maxDJR, maxPt, minDelta, minDJR, minPt, ptClustering, ptFraction, tmp, and useEt.

                                                      :
        sourceLabel(params.getParameter<edm::InputTag>("src")),
        jetInput(params.getParameter<edm::ParameterSet>("jetInput")),
        defaultDeltaCut(params.getParameter<double>("defaultDeltaCut")),
        defaultPtCut(params.getParameter<double>("defaultPtCut")),
        maxEta(params.getParameter<double>("maxEta")),
        useEt(params.getParameter<bool>("useEt")),
        ptFraction(params.getUntrackedParameter<double>("ptFraction", 0.75)),
        binsDelta(params.getParameter<unsigned int>("binsDelta")),
        minDelta(params.getParameter<double>("minDelta")),
        maxDelta(params.getParameter<double>("maxDelta")),
        binsPt(params.getParameter<unsigned int>("binsPt")),
        minPt(params.getParameter<double>("minPt")),
        maxPt(params.getParameter<double>("maxPt")),
        minDJR(params.getParameter<unsigned int>("minDJR")),
        maxDJR(params.getParameter<unsigned int>("maxDJR"))
{
        edm::ParameterSet jetClusPSet =
                params.getParameter<edm::ParameterSet>("jetClustering");

        for(unsigned int i = 0; i <= binsDelta; i++) {
                double deltaCut =
                        minDelta + (maxDelta - minDelta) * i / binsDelta;
                jetClusPSet.addParameter("coneRadius", deltaCut);
                edm::ParameterSet tmp;
                tmp.addParameter("algorithm", jetClusPSet);
                deltaClustering.push_back(
                        new JetClustering(tmp, defaultPtCut * ptFraction));
        }

        jetClusPSet.addParameter("coneRadius", defaultDeltaCut);
        edm::ParameterSet tmp;
        tmp.addParameter("algorithm", jetClusPSet);
        ptClustering.reset(new JetClustering(tmp, minPt * ptFraction));

        fillDJRSched(minDJR <= 0 ? 1 : minDJR, maxDJR - 1);

        edm::Service<TFileService> fs;
        for(unsigned int i = minDJR; i < maxDJR; i++) {
                std::ostringstream ss, ss2;
                ss << (i + 1) << "#rightarrow" << i << " jets";
                ss2 << i;
                TH1 *h = fs->make<TH1D>(("delta" + ss2.str()).c_str(),
                                        ("DJR " + ss.str()).c_str(),
                                        binsDelta, minDelta, maxDelta);
                h->SetXTitle("p_{T} [GeV/c^2]");
                h->SetYTitle("#delta#sigma [mb]");

                if (i == 0) {
                        h->Delete();
                        h = 0;
                }
                histoDelta.push_back(h);

                std::string what = useEt ? "E" : "p";
                h = fs->make<TH1D>(("pt" + ss2.str()).c_str(),
                                   ("DJR " + ss.str()).c_str(), binsPt,
                                   std::log10(minPt), std::log10(maxPt));
                h->SetXTitle(("log_{10}(" + what +
                                        "_{T} [GeV/c^{2}])").c_str());
                h->SetYTitle("#delta#sigma [mb]");

                histoPt.push_back(h);
        }
}

LHEAnalyzer::~LHEAnalyzer ( ) [virtual]

Definition at line 143 of file LHEAnalyzer.cc.

{
}

Member Function Documentation

void LHEAnalyzer::analyze	(	const edm::Event &	event,
		const edm::EventSetup &	es
	)		`[protected, virtual]`

Implements edm::EDAnalyzer.

Definition at line 159 of file LHEAnalyzer.cc.

References binsDelta, defaultPtCut, deltaClustering, djrSched, first, MCTruth::genEvent, h, histoDelta, histoPt, i, j, jetInput, analyzePatCleaning_cfg::jets, maxDJR, minDJR, query::result, edm::second(), python::multivaluedict::sort(), sourceLabel, pileupCalc::upper, and useEt.

{
        using boost::bind;
        typedef JetClustering::Jet Jet;

        edm::Handle<edm::HepMCProduct> hepmc;
        event.getByLabel(sourceLabel, hepmc);

        std::auto_ptr<HepMC::GenEvent> clonedEvent;
        const HepMC::GenEvent *genEvent = hepmc->GetEvent();
        if (!genEvent->signal_process_vertex()) {
                clonedEvent.reset(new HepMC::GenEvent(*genEvent));
                const HepMC::GenVertex *signalVertex =
                        LHEEvent::findSignalVertex(clonedEvent.get());
                clonedEvent->set_signal_process_vertex(
                        const_cast<HepMC::GenVertex*>(signalVertex));
                genEvent = clonedEvent.get();
        }

        JetInput::ParticleVector particles = jetInput(genEvent);

        std::vector<Jet> ptJets = (*ptClustering)(particles);
        std::sort(ptJets.begin(), ptJets.end(),
                  bind(std::greater<double>(),
                       bind(useEt ? &Jet::et : &Jet::pt, _1),
                       bind(useEt ? &Jet::et : &Jet::pt, _2)));

        typedef std::pair<int, int> Pair;
        std::vector<Pair> deltaJets(maxDJR - minDJR + 1,
                                    Pair(-1, binsDelta + 1));

        for(std::vector<unsigned int>::const_iterator djr = djrSched.begin();
            djr != djrSched.end(); ++djr) {
//std::cout << "DJR schedule " << (*djr + 1) << " -> " << *djr << std::endl;
                int result = -1;
                for(;;) {
//for(int i = minDJR; i <= maxDJR; i++)
//std::cout << "+++ " << i << ": (" << deltaJets[i - minDJR].first << ", " << deltaJets[i - minDJR].second << ")" << std::endl;
                        int upper = binsDelta + 1;
                        for(int i = *djr; i >= (int)minDJR; i--) {
                                if (deltaJets[i - minDJR].second <=
                                                        (int)binsDelta) {
                                        upper = deltaJets[i - minDJR].second;
                                        break;
                                }
                        }
                        int lower = -1;
                        for(int i = *djr + 1; i <= (int)maxDJR; i++) {
                                if (deltaJets[i - minDJR].first >= 0) {
                                        lower = deltaJets[i - minDJR].first;
                                        break;
                                }
                        }
//std::cout << "\t" << lower << " < " << upper << std::endl;

                        result = (lower + upper + 2) / 2 - 1;
                        if (result == lower)
                                break;
                        else if (result < lower) {
                                result = -1;
                                break;
                        }

                        std::vector<Jet> jets =
                                deltaClustering[result](particles);
                        unsigned int nJets = 0;
                        for(std::vector<Jet>::const_iterator iter =
                                jets.begin(); iter != jets.end(); ++iter)
                                if ((useEt ? iter->et() : iter->pt())
                                                        > defaultPtCut)
                                        nJets++;

//std::cout << "\t---(" << *djr << ")--> bin " << result << ": " << nJets << " jets" << std::endl;

                        if (nJets < minDJR)
                                nJets = minDJR;
                        else if (nJets > maxDJR)
                                nJets = maxDJR;

                        for(int j = nJets; j >= (int)minDJR; j--) {
                                if (deltaJets[j - minDJR].first < 0 ||
                                    result > deltaJets[j - minDJR].first)
                                        deltaJets[j - minDJR].first = result;
                        }
                        for(int j = nJets; j <= (int)maxDJR; j++) {
                                if (deltaJets[j - minDJR].second <
                                                        (int)binsDelta ||
                                    result < deltaJets[j - minDJR].second)
                                        deltaJets[j - minDJR].second = result;
                        }
                }

//std::cout << "final " << *djr << ": " << result << std::endl;
                TH1 *h = histoDelta[*djr - minDJR];
                h->Fill(h->GetBinCenter(result + 1));

                h = histoPt[*djr - minDJR];
                if (*djr >= ptJets.size())
                        h->Fill(-999.0);
                else if (useEt)
                        h->Fill(std::log10(ptJets[*djr].et()));
                else
                        h->Fill(std::log10(ptJets[*djr].pt()));
        }

        if (minDJR <= 0) {
                TH1 *h = histoPt[0];
                if (minDJR >= ptJets.size())
                        h->Fill(-999.0);
                else if (useEt)
                        h->Fill(std::log10(ptJets[minDJR].et()));
                else
                        h->Fill(std::log10(ptJets[minDJR].pt()));
        }
}

void LHEAnalyzer::fillDJRSched	(	unsigned int	min,
		unsigned int	max
	)		`[private]`

Definition at line 147 of file LHEAnalyzer.cc.

References djrSched, and max().

Referenced by LHEAnalyzer().

{
        unsigned int middle = (min + max) / 2;

        djrSched.push_back(middle);

        if (min < middle)
                fillDJRSched(min, middle - 1);
        if (middle < max)
                fillDJRSched(middle + 1, max);
}