Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Make running doAll_Zp more flexible and setup for condor jobs #57

Merged
merged 18 commits into from
Sep 19, 2022
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
32 changes: 26 additions & 6 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -27,23 +27,43 @@ For more dedtails, refer to `cpp/README.md`.

### Example instructions

Edit `cpp/doAll_Zp.C` with an appropriate file (or hopefully the default one still exists).
Edit `cpp/doAll_Zp.C` and/or `cpp/ScanChain_Zp.C` with an appropriate file (or hopefully the default one still exists).

**To run locally:**

```bash
pushd cpp/
root -b -q -l -n doAll_Zp.C
popd
sh cpp/runOutput_Zp.sh DIRECTORY/FOR/OUTPUT/FILES YEAR RUNDATA RUNBKG RUNSIGNAL RUNBFF SAMPLE ADDITIONALBOOLEANFLAGS
```

For example, to run all data, bkg and signal but not the BFF samples, for all years and with all additional flags to their default values, the command to save the files in a folder called "temp_data" would be:

```bash
sh cpp/runOutput_Zp.sh temp_data all 1 1 1 0 all
```

One should check `cpp/doAll_Zp.C` for details on the (additional) arguments and their meaning.

This loops and creates a number of output files of the form `output_"process"_"year".root` containing histograms.
Optionally, the file also contains a `RooDataSet` to be used as input for fitting (see below).

To produce plots:
**To run on condor:**

```bash
voms-proxy-init -voms cms
export X509_USER_PROXY=/ABSOLUTE/PATH/TO/PROXY/LOCATION
sh utils/condor/runOutput_Zp_onCondor.sh FOLDER/FOR/OUTPUT/FILES
```

This script will package the current state of the repository and send it to condor jobs running the `cpp/runOutput_Zp.sh` with the arguments included in the different lines of `utils/condor/runOutput_Zp_onCondor.sub` (please edit this file to control what condor jobs you send). More details on this are included in the description of [PR#57](https://github.com/cmstas/ZPrimeSnT/pull/57).

The output of your jobs will be found under `/ceph/cms/store/user/$USER/ZPrimeSnTOutput/FOLDER/FOR/OUTPUT/FILES` and the plotting logs under `utils/condor/plotting_logs`.

**To produce plots:**
```bash
python python/stack_plots.py
```

To produce cutflow table:
**To produce cutflow table:**
```bash
python python/make_cutflow_table.py
```
Expand Down
125 changes: 110 additions & 15 deletions cpp/ScanChain_Zp.C
Original file line number Diff line number Diff line change
Expand Up @@ -107,7 +107,7 @@ using namespace duplicate_removal;
using namespace RooFit;


int ScanChain(TChain *ch, double genEventSumw, TString year, TString process, int prefireWeight=1, int topPtWeight=1, int PUWeight=1, int muonRecoSF=1, int muonIdSF=1, int muonIsoSF=1, int triggerSF=1, int bTagSF=1, int JECUnc=0, int JERUnc=0, int runBFFSync=0, const char* outdir="temp_data") {
int ScanChain(TChain *ch, double genEventSumw, TString year, TString process, int prefireWeight=1, int topPtWeight=1, int PUWeight=1, int muonRecoSF=1, int muonIdSF=1, int muonIsoSF=1, int muonResUnc=0, int triggerSF=1, int bTagSF=1, int JECUnc=0, int JERUnc=0, int UnclEnUnc=0, int runBFFSync=0, const char* outdir="temp_data") {
// Event weights / scale factors:
// 0: Do not apply
// 1: Apply central value
Expand All @@ -125,10 +125,12 @@ int ScanChain(TChain *ch, double genEventSumw, TString year, TString process, in
muonRecoSF = 0;
muonIdSF = 0;
muonIsoSF = 0;
muonResUnc=0;
triggerSF = 0;
bTagSF = 0;
JECUnc = 0;
JERUnc = 0;
UnclEnUnc = 0;
}

int mdir = mkdir(outdir,0755);
Expand Down Expand Up @@ -336,7 +338,7 @@ int ScanChain(TChain *ch, double genEventSumw, TString year, TString process, in

// Define RooDataSet's for fit
RooRealVar mfit("mfit", "mfit", 175.0, 6500.0);
RooRealVar roow("roow", "roow", 0.0, 100.0);
RooRealVar roow("roow", "roow", -10000.0, 10000.0);
map<TString, RooDataSet> roods;
for ( unsigned int imll=0; imll < mllbin.size(); imll++ ) {
for ( unsigned int inb=0; inb < nbtag.size(); inb++ ) {
Expand Down Expand Up @@ -407,6 +409,10 @@ int ScanChain(TChain *ch, double genEventSumw, TString year, TString process, in
plot_names.push_back("mu2_dxy");
plot_names.push_back("mu1_dz");
plot_names.push_back("mu2_dz");
plot_names.push_back("mu1_relPtErr");
plot_names.push_back("mu2_relPtErr");
plot_names.push_back("mu1_tunepRelPt");
plot_names.push_back("mu2_tunepRelPt");
plot_names.push_back("mu1_trkRelIso");
plot_names.push_back("mu1_trkAbsIso");
plot_names.push_back("mu2_trkRelIso");
Expand All @@ -427,6 +433,10 @@ int ScanChain(TChain *ch, double genEventSumw, TString year, TString process, in
plot_names.push_back("mu2_dxy");
plot_names.push_back("mu1_dz");
plot_names.push_back("mu2_dz");
plot_names.push_back("mu1_relPtErr");
plot_names.push_back("mu2_relPtErr");
plot_names.push_back("mu1_tunepRelPt");
plot_names.push_back("mu2_tunepRelPt");
plot_names.push_back("mu1_trkRelIso");
plot_names.push_back("mu1_trkAbsIso");
plot_names.push_back("mu2_trkRelIso");
Expand Down Expand Up @@ -572,6 +582,10 @@ int ScanChain(TChain *ch, double genEventSumw, TString year, TString process, in
plot_names.push_back("mu2_dxy");
plot_names.push_back("mu1_dz");
plot_names.push_back("mu2_dz");
plot_names.push_back("mu1_relPtErr");
plot_names.push_back("mu2_relPtErr");
plot_names.push_back("mu1_tunepRelPt");
plot_names.push_back("mu2_tunepRelPt");
plot_names.push_back("mu1_trkRelIso");
plot_names.push_back("mu1_trkAbsIso");
plot_names.push_back("mu2_trkRelIso");
Expand Down Expand Up @@ -651,6 +665,9 @@ int ScanChain(TChain *ch, double genEventSumw, TString year, TString process, in
if ( triggerSF!=0 ) set_triggerSF();
if ( bTagSF!=0 ) set_allbTagEff();

// Setting up muon momentum resolution
TRandom3 rnd_muonMomRes(12345);

// Setting up btagging scale factors
BTagCalibration_v2* btagCalib;
BTagCalibrationReader_v2* btagReaderTight = new BTagCalibrationReader_v2(BTagEntry_v2::OP_TIGHT, "central", {"up", "down"});
Expand Down Expand Up @@ -826,7 +843,7 @@ int ScanChain(TChain *ch, double genEventSumw, TString year, TString process, in


if ( isMC ) {
if(removeSpikes && weight*factor>1e2) continue;
if(removeSpikes && fabs(weight*factor)>1e2) continue;
if (process=="DYbb") {
int nGluons = 0;
for ( unsigned int p=0; p<nt.nLHEPart(); p++ ) {
Expand Down Expand Up @@ -943,6 +960,16 @@ int ScanChain(TChain *ch, double genEventSumw, TString year, TString process, in
continue;
if ( isinf(nt.PuppiMET_phiJERDown()) || isnan(nt.PuppiMET_phiJERDown()) )
continue;
// Unclustered up
if ( isinf(nt.PuppiMET_ptUnclusteredUp()) || isnan(nt.PuppiMET_ptUnclusteredUp()) )
continue;
if ( isinf(nt.PuppiMET_phiUnclusteredUp()) || isnan(nt.PuppiMET_phiUnclusteredUp()) )
continue;
// Unclustered down
if ( isinf(nt.PuppiMET_ptUnclusteredDown()) || isnan(nt.PuppiMET_ptUnclusteredDown()) )
continue;
if ( isinf(nt.PuppiMET_phiUnclusteredDown()) || isnan(nt.PuppiMET_phiUnclusteredDown()) )
continue;
}
}

Expand All @@ -969,6 +996,14 @@ int ScanChain(TChain *ch, double genEventSumw, TString year, TString process, in
puppimet_pt = nt.PuppiMET_ptJERDown();
puppimet_phi = nt.PuppiMET_phiJERDown();
}
if ( UnclEnUnc==2 ) {
puppimet_pt = nt.PuppiMET_ptUnclusteredUp();
puppimet_phi = nt.PuppiMET_phiUnclusteredUp();
}
if ( UnclEnUnc==-2 ) {
puppimet_pt = nt.PuppiMET_ptUnclusteredDown();
puppimet_phi = nt.PuppiMET_phiUnclusteredDown();
}
}

if ( !removeCorrections ) {
Expand Down Expand Up @@ -999,9 +1034,29 @@ int ScanChain(TChain *ch, double genEventSumw, TString year, TString process, in
float triggerWeight = -1.0;
//
for ( unsigned int mu = 0; mu < nt.nMuon(); mu++ ) {
Muon_pt.push_back(nt.Muon_pt().at(mu));
Muon_p4.push_back(LorentzVector(nt.Muon_pt().at(mu),nt.Muon_eta().at(mu),nt.Muon_phi().at(mu),nt.Muon_mass().at(mu)));
Muon_tkRelIso.push_back(nt.Muon_tkRelIso().at(mu));
if ( isMC && muonResUnc != 0 ) { // 2 means that variation are to be applied
// https://twiki.cern.ch/twiki/bin/view/CMS/MuonUL2018#Momentum_Resolution
// https://twiki.cern.ch/twiki/bin/view/CMS/MuonUL2017#Momentum_Resolution
// https://twiki.cern.ch/twiki/bin/view/CMS/MuonUL2016#Momentum_Resolution
float a, b, c;
if ( fabs(Muon_p4[mu].Eta()) < 1.2 ) {
a = ( year == "2018" ? 0.0141926 : ( year == "2017" ? 0.0145351 : 0.0146577 ) );
b = ( year == "2018" ? 4.23456e-05 : ( year == "2017" ? 4.24022e-05 : 4.48517e-05 ) );
c = ( year == "2018" ? -9.91644e-09 : ( year == "2017" ? -1.01461e-08 : -9.87206e-09 ) );
}
else {
a = ( year == "2018" ? 0.016883 : ( year == "2017" ? 0.0168087 : 0.0155794 ) );
b = ( year == "2018" ? 6.21438e-05 : ( year == "2017" ? 5.57382e-05 : 7.22117e-05 ) );
c = ( year == "2018" ? -9.79418e-09 : ( year == "2017" ? -8.89713e-09 : -1.27255e-08 ) );
}
float muonP = Muon_p4[mu].P();
float muonResSmearParam = a + b * muonP + c * muonP * muonP;
double muonResSmearFactor = rnd_muonMomRes.Gaus(0,muonResSmearParam*0.46);
Muon_p4[mu].SetPt( Muon_p4[mu].Pt()*(1 + 0.5*muonResUnc*muonResSmearFactor) ); // 0.5*muonResUnc = sign of variation
}
Muon_pt.push_back(Muon_p4[mu].Pt());
Muon_tkRelIso.push_back(nt.Muon_tkRelIso().at(mu)); // We don't change the relative isolation based on the muon momentum resolution, as this seems to be out of the scope of this uncertainty
//
bool mu_trk_and_global = ( nt.Muon_isGlobal().at(mu) && nt.Muon_isTracker().at(mu) );
bool mu_id = ( nt.Muon_highPtId().at(mu) >= 2 );
Expand Down Expand Up @@ -1163,6 +1218,18 @@ int ScanChain(TChain *ch, double genEventSumw, TString year, TString process, in
plot_names.push_back("mu2_dz");
variable.insert({"mu2_dz", fabs(nt.Muon_dz().at(1))});

plot_names.push_back("mu1_relPtErr");
variable.insert({"mu1_relPtErr", nt.Muon_ptErr().at(0) / Muon_pt.at(0)});

plot_names.push_back("mu2_relPtErr");
variable.insert({"mu2_relPtErr", nt.Muon_ptErr().at(1) / Muon_pt.at(1)});

plot_names.push_back("mu1_tunepRelPt");
variable.insert({"mu1_tunepRelPt", nt.Muon_tunepRelPt().at(0)});

plot_names.push_back("mu2_tunepRelPt");
variable.insert({"mu2_tunepRelPt", nt.Muon_tunepRelPt().at(1)});

plot_names.push_back("mu1_trkRelIso");
variable.insert({"mu1_trkRelIso", Muon_tkRelIso.at(0)});

Expand Down Expand Up @@ -1452,6 +1519,10 @@ int ScanChain(TChain *ch, double genEventSumw, TString year, TString process, in
variable["mu2_dxy"] = fabs(nt.Muon_dxy().at(cand_muons_id[1]));
variable["mu1_dz"] = fabs(nt.Muon_dz().at(cand_muons_id[0]));
variable["mu2_dz"] = fabs(nt.Muon_dz().at(cand_muons_id[1]));
variable["mu1_relPtErr"] = nt.Muon_ptErr().at(cand_muons_id[0]) / Muon_pt.at(cand_muons_id[0]);
variable["mu2_relPtErr"] = nt.Muon_ptErr().at(cand_muons_id[1]) / Muon_pt.at(cand_muons_id[1]);
variable["mu1_tunepRelPt"] = nt.Muon_tunepRelPt().at(cand_muons_id[0]);
variable["mu2_tunepRelPt"] = nt.Muon_tunepRelPt().at(cand_muons_id[1]);
variable["mu1_trkRelIso"] = Muon_tkRelIso.at(cand_muons_id[0]);
variable["mu1_trkAbsIso"] = Muon_tkRelIso.at(cand_muons_id[0]) * Muon_pt.at(cand_muons_id[0]);
variable["mu2_trkRelIso"] = Muon_tkRelIso.at(cand_muons_id[1]);
Expand Down Expand Up @@ -1500,6 +1571,10 @@ int ScanChain(TChain *ch, double genEventSumw, TString year, TString process, in
variable["mu2_dxy"] = fabs(nt.Muon_dxy().at(cand_muons_id_pteta[1]));
variable["mu1_dz"] = fabs(nt.Muon_dz().at(cand_muons_id_pteta[0]));
variable["mu2_dz"] = fabs(nt.Muon_dz().at(cand_muons_id_pteta[1]));
variable["mu1_relPtErr"] = nt.Muon_ptErr().at(cand_muons_id_pteta[0]) / Muon_pt.at(cand_muons_id_pteta[0]);
variable["mu2_relPtErr"] = nt.Muon_ptErr().at(cand_muons_id_pteta[1]) / Muon_pt.at(cand_muons_id_pteta[1]);
variable["mu1_tunepRelPt"] = nt.Muon_tunepRelPt().at(cand_muons_id_pteta[0]);
variable["mu2_tunepRelPt"] = nt.Muon_tunepRelPt().at(cand_muons_id_pteta[1]);
variable["mu1_trkRelIso"] = Muon_tkRelIso.at(cand_muons_id_pteta[0]);
variable["mu1_trkAbsIso"] = Muon_tkRelIso.at(cand_muons_id_pteta[0]) * Muon_pt.at(cand_muons_id_pteta[0]);
variable["mu2_trkRelIso"] = Muon_tkRelIso.at(cand_muons_id_pteta[1]);
Expand Down Expand Up @@ -1549,6 +1624,10 @@ int ScanChain(TChain *ch, double genEventSumw, TString year, TString process, in
variable["mu2_dxy"] = fabs(nt.Muon_dxy().at(cand_muons[1]));
variable["mu1_dz"] = fabs(nt.Muon_dz().at(cand_muons[0]));
variable["mu2_dz"] = fabs(nt.Muon_dz().at(cand_muons[1]));
variable["mu1_relPtErr"] = nt.Muon_ptErr().at(cand_muons[0]) / Muon_pt.at(cand_muons[0]);
variable["mu2_relPtErr"] = nt.Muon_ptErr().at(cand_muons[1]) / Muon_pt.at(cand_muons[1]);
variable["mu1_tunepRelPt"] = nt.Muon_tunepRelPt().at(cand_muons[0]);
variable["mu2_tunepRelPt"] = nt.Muon_tunepRelPt().at(cand_muons[1]);
variable["mu1_trkRelIso"] = Muon_tkRelIso.at(cand_muons[0]);
variable["mu1_trkAbsIso"] = Muon_tkRelIso.at(cand_muons[0]) * Muon_pt.at(cand_muons[0]);
variable["mu2_trkRelIso"] = Muon_tkRelIso.at(cand_muons[1]);
Expand Down Expand Up @@ -1732,6 +1811,10 @@ int ScanChain(TChain *ch, double genEventSumw, TString year, TString process, in
variable["mu2_dxy"] = fabs(nt.Muon_dxy().at(subleadingMu_idx));
variable["mu1_dz"] = fabs(nt.Muon_dz().at(leadingMu_idx));
variable["mu2_dz"] = fabs(nt.Muon_dz().at(subleadingMu_idx));
variable["mu1_relPtErr"] = nt.Muon_ptErr().at(leadingMu_idx) / Muon_pt.at(leadingMu_idx);
variable["mu2_relPtErr"] = nt.Muon_ptErr().at(subleadingMu_idx) / Muon_pt.at(subleadingMu_idx);
variable["mu1_tunepRelPt"] = nt.Muon_tunepRelPt().at(leadingMu_idx);
variable["mu2_tunepRelPt"] = nt.Muon_tunepRelPt().at(subleadingMu_idx);
variable["mu1_trkRelIso"] = Muon_tkRelIso.at(leadingMu_idx);
variable["mu1_trkAbsIso"] = Muon_tkRelIso.at(leadingMu_idx) * Muon_pt.at(leadingMu_idx);
variable["mu2_trkRelIso"] = Muon_tkRelIso.at(subleadingMu_idx);
Expand Down Expand Up @@ -1769,6 +1852,18 @@ int ScanChain(TChain *ch, double genEventSumw, TString year, TString process, in
plot_names.push_back("mu2_dz");
variable.insert({"mu2_dz", fabs(nt.Muon_dz().at(subleadingMu_idx))});

plot_names.push_back("mu1_relPtErr");
variable.insert({"mu1_relPtErr", nt.Muon_ptErr().at(leadingMu_idx) / Muon_pt.at(leadingMu_idx)});

plot_names.push_back("mu2_relPtErr");
variable.insert({"mu2_relPtErr", nt.Muon_ptErr().at(subleadingMu_idx) / Muon_pt.at(subleadingMu_idx)});

plot_names.push_back("mu1_tunepRelPt");
variable.insert({"mu1_tunepRelPt", nt.Muon_tunepRelPt().at(leadingMu_idx)});

plot_names.push_back("mu2_tunepRelPt");
variable.insert({"mu2_tunepRelPt", nt.Muon_tunepRelPt().at(subleadingMu_idx)});

plot_names.push_back("mu1_trkRelIso");
variable.insert({"mu1_trkRelIso", Muon_tkRelIso.at(leadingMu_idx)});

Expand Down Expand Up @@ -2545,16 +2640,16 @@ int ScanChain(TChain *ch, double genEventSumw, TString year, TString process, in
if ( removeDataDuplicates && !(isMC) )
cout << "Number of duplicates found: " << nDuplicates << endl;

// Avoid histograms with unphysical negative bin content (due to negative GEN weights)
map<TString, TH1F*>::iterator it;
for ( it = histos.begin(); it != histos.end(); it++ ) {
for ( unsigned int b=1; b<(it->second)->GetNbinsX()+1; b++ ) {
if ( (it->second)->GetBinContent(b)<0.0) {
(it->second)->SetBinContent(b,0.0);
(it->second)->SetBinError(b,0.0);
}
}
}
//// Avoid histograms with unphysical negative bin content (due to negative GEN weights)
//map<TString, TH1F*>::iterator it;
//for ( it = histos.begin(); it != histos.end(); it++ ) {
// for ( unsigned int b=1; b<(it->second)->GetNbinsX()+1; b++ ) {
// if ( (it->second)->GetBinContent(b)<0.0) {
// (it->second)->SetBinContent(b,0.0);
// (it->second)->SetBinError(b,0.0);
// }
// }
//}

if ( writeOutYields_BeforeSel ) {
TString model = ((TObjString *)process.Tokenize("_")->At(0))->String();
Expand Down
38 changes: 29 additions & 9 deletions cpp/configuration_Zp.h
Original file line number Diff line number Diff line change
Expand Up @@ -119,38 +119,58 @@ void histoDefinition(map<TString, int> &nbins, map<TString, float> &low, map<TSt
nbins.insert({"mu1_dxy", 200});
low.insert({"mu1_dxy", 0.001});
high.insert({"mu1_dxy", 0.021});
title.insert({"mu1_dxy", "d_{xy} (leading #mu)"});
title.insert({"mu1_dxy", "|d_{xy}| (leading #mu) [cm]"});

nbins.insert({"mu2_dxy", 200});
low.insert({"mu2_dxy", 0.001});
high.insert({"mu2_dxy", 0.021});
title.insert({"mu2_dxy", "d_{xy} (subleading #mu)"});
title.insert({"mu2_dxy", "|d_{xy}| (subleading #mu) [cm]"});

nbins.insert({"mu1_dz", 100});
low.insert({"mu1_dz", 0.001});
high.insert({"mu1_dz", 0.101});
title.insert({"mu1_dz", "d_{z} (leading #mu)"});
title.insert({"mu1_dz", "|d_{z}| (leading #mu) [cm]"});

nbins.insert({"mu2_dz", 100});
low.insert({"mu2_dz", 0.001});
high.insert({"mu2_dz", 0.101});
title.insert({"mu2_dz", "d_{z} (subleading #mu)"});
title.insert({"mu2_dz", "|d_{z}| (subleading #mu) [cm]"});

nbins.insert({"mu1_relPtErr", 100});
low.insert({"mu1_relPtErr", 0.001});
high.insert({"mu1_relPtErr", 1.001});
title.insert({"mu1_relPtErr", "Relative p_{T} error (leading #mu)"});

nbins.insert({"mu2_relPtErr", 100});
low.insert({"mu2_relPtErr", 0.001});
high.insert({"mu2_relPtErr", 1.001});
title.insert({"mu2_relPtErr", "Relative p_{T} error (subleading #mu)"});

nbins.insert({"mu1_tunepRelPt", 100});
low.insert({"mu1_tunepRelPt", 0.001});
high.insert({"mu1_tunepRelPt", 2.001});
title.insert({"mu1_tunepRelPt", "TuneP relative p_{T} (leading #mu)"});

nbins.insert({"mu2_tunepRelPt", 100});
low.insert({"mu2_tunepRelPt", 0.001});
high.insert({"mu2_tunepRelPt", 2.001});
title.insert({"mu2_tunepRelPt", "TuneP relative p_{T} (subleading #mu)"});

nbins.insert({"mu1_trkRelIso", 100});
low.insert({"mu1_trkRelIso", 0.001});
high.insert({"mu1_trkRelIso", 0.051});
title.insert({"mu1_trkRelIso", "Track iso./p_{T} (leading #mu)"});

nbins.insert({"mu1_trkAbsIso", 100});
low.insert({"mu1_trkAbsIso", 0.1});
high.insert({"mu1_trkAbsIso", 5.1});
title.insert({"mu1_trkAbsIso", "Track iso. (leading #mu)"});

nbins.insert({"mu2_trkRelIso", 100});
low.insert({"mu2_trkRelIso", 0.001});
high.insert({"mu2_trkRelIso", 0.051});
title.insert({"mu2_trkRelIso", "Track iso./p_{T} (subleading #mu)"});

nbins.insert({"mu1_trkAbsIso", 100});
low.insert({"mu1_trkAbsIso", 0.1});
high.insert({"mu1_trkAbsIso", 5.1});
title.insert({"mu1_trkAbsIso", "Track iso. (leading #mu)"});

nbins.insert({"mu2_trkAbsIso", 100});
low.insert({"mu2_trkAbsIso", 0.1});
high.insert({"mu2_trkAbsIso", 5.1});
Expand Down
Loading