applyDataLumiZeroBias.cc File Reference

#include <cstdlib>
#include <cassert>
#include <iostream>
#include <vector>
#include <map>
#include <filesystem>
#include "Core/Objects/interface/Jet.h"
#include "Core/Objects/interface/Event.h"
#include "Core/CommonTools/interface/terminal.h"
#include "Core/CommonTools/interface/Looper.h"
#include "Core/CommonTools/interface/MetaInfo.h"
#include "Core/CommonTools/interface/toolbox.h"
#include "Core/CommonTools/interface/ControlPlots.h"
#include <TROOT.h>
#include <TFile.h>
#include <TChain.h>
#include <TH2.h>
#include "TriggerEff.h"

Include dependency graph for applyDataLumiZeroBias.cc:

Functions
void	applyDataPrescalesZeroBias (const fs::path &input, const fs::path &output, auto total_lumi, auto maxpt, int nSplit=1, int nNow=0)
int	main (int argc, char *argv[])

Function Documentation

applyDataPrescalesZeroBias()

void applyDataPrescalesZeroBias	(	const fs::path &	input,
		const fs::path &	output,
		auto	total_lumi,
		auto	maxpt,
		int	nSplit = 1,
		int	nNow = 0
	)

Normalise with prescales.

Parameters

input	name of input root file
output	name of output root file
nSplit	number of jobs/tasks/cores
nNow	index of job/task/core

{
    assert(fs::exists(input));
    TChain * oldchain = new TChain("events"); 
    oldchain->Add(input.c_str());
 
    vector<RecJet> * recJets = nullptr;
    Event * event = nullptr;
    Trigger * trigger = nullptr;
    oldchain->SetBranchAddress("recJets", &recJets);
    oldchain->SetBranchAddress("event", &event);
    oldchain->SetBranchAddress("trigger", &trigger);
 
    if (fs::exists(output)) 
        cerr << red << output << " will be overwritten\n" << normal;
    auto newfile = TFile::Open(output.c_str(), "RECREATE");
    auto newtree = oldchain->CloneTree(0);
 
    MetaInfo metainfo(newtree);
    if (nNow == 0) metainfo.Print();
    metainfo.AddCorrection("normalised");
 
    int year = metainfo.year();
    assert(year > 2015 && year < 2019);
 
    assert(total_lumi > 0);
    auto inv_total_lumi = 1./total_lumi;
 
    newfile->cd();
 
    // a few control plots
    ControlPlots::isMC = false;
    vector<ControlPlots> corr ;
    for (int i=0; i<metainfo.GetNRecEvWgts(); i++)
        corr.push_back(ControlPlots(metainfo.GetRecEvWgt(i)));
 
    cout << "looping over events:" << endl;
 
    Looper looper(__func__, oldchain, nSplit, nNow);
    while (looper.Next()) {
 
        // we find the leading jet in tracker acceptance
        auto leadingInTk = recJets->begin();
        while (leadingInTk != recJets->end()
                && abs(leadingInTk->Rapidity()) >= 2.5)
            ++leadingInTk;
 
        // removed events alreadu covered by jet triggers
        if (leadingInTk != recJets->end() && leadingInTk->CorrPt() >= maxpt) continue;
 
        // get prescale
        auto preHLT   = trigger->PreHLT.front();
        auto preL1min = trigger->PreL1min.front();
        auto preL1max = trigger->PreL1max.front();
        assert(preL1min == preL1max);
        if (preL1min != preL1max)
            cerr << "\x1B[31m\e[1m" << preL1min << ' ' <<  preL1max << "\x1B[30m\e[0m\n";
 
        auto prescale = preHLT * preL1min;
 
        // sanity checks:
 
        // 1) we assume that the same prescale is indeed constant for a given LS
        {
            static vector<map<pair<int,int>,int>> prescales(1); // this could be written in a simpler way...
            pair<int,int> runlum = {event->runNo, event->lumi};
            if (prescales.front().count(runlum)) 
                assert(prescales.front().at(runlum) == prescale);
            else 
                prescales.front()[runlum] = prescale;
        }
 
        // 2) check that the data have not yet been rescaled
        assert(event->genWgts.size() == 0 && event->recWgts.size() == 1 && event->recWgts.front() == 1);
 
        // setting the inverse of the eff lumi to the event including correction of for the trigger efficiency
        for (size_t i = 0; i < corr.size(); ++i)
            event->recWgts.at(i) *= prescale * inv_total_lumi;
 
        newtree->Fill();
        for(size_t i=0; i<corr.size(); i++){
            if(recJets->size()==0) continue;
            corr.at(i)(*recJets, event->recWgts.at(i));
        }
    }
 
    cout << "saving" << endl;
    newtree->AutoSave();
    for (size_t i = 0; i < corr.size(); ++i)
        corr.at(i).Write(newfile);
    auto controlplots = newfile->mkdir("controlplots");
    controlplots->cd();
    newfile->Close();
    delete oldchain;
}

main()

int main	(	int	argc,
		char *	argv[]
	)

{
    gROOT->SetBatch();
 
    if (argc < 4) {
        cout << argv[0] << " input output lumi_file turnon_file trigger_curves [nSplit [nNow]]\n"
             << "where\tinput = data n-tuple\n"
             << "     \toutput = new n-tuple\n"
             << "     \ttotal_lumi = lumi in /pb (can be found in `$DAS_WORKAREA/tables/lumi`)\n"
             << "     \tmaxpt = max pt to be covered by ZeroBias trigger (suggestion: lowest pt threshold in `$DAS_WORKAREA/tables/trigger`)\n"
             << flush;
        return EXIT_SUCCESS;
    }
 
    fs::path input = argv[1],
             output = argv[2];
    auto total_lumi = atof(argv[3]),
         maxpt = atof(argv[4]);
 
    int nNow = 0, nSplit = 1;
    if (argc > 5) nSplit = atoi(argv[5]);
    if (argc > 6) nNow = atoi(argv[6]);
 
    applyDataPrescalesZeroBias(input, output, total_lumi, maxpt, nSplit, nNow);
    return EXIT_SUCCESS;
}