wprime/CreateAuxHists.cc at master · sifuluo/wprime · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
#include <TChain.h>
#include <TString.h>
#include <vector>
#include <TLorentzVector.h>
#include <iostream>
#include <cmath>
#include <array>

#include <TApplication.h>

#include "Utilities/NanoAODReader.cc"
#include "Utilities/JetScale.cc"
#include "Utilities/Permutations.cc"
#include "Utilities/Hypothesis.cc"

void CreateAuxHists(int sampleyear = 3, int sampletype = 22, int ifile = -1, string infile = "All") {
  dlib.AppendAndrewDatasets();
  // infile = "/eos/user/p/pflanaga/andrewsdata/skimmed_samples/wprime_500/2017/003C756A-BC7B-5A48-8F2E-A61D3CDC7C32.root";
  // sampletype = dlib.AddDataset_NGTCXS("WP500");
  if (sampletype < 0) {
    dlib.ListDatasets();
    cout << "Please enter the index for a dataset:" <<endl;
    cin >> sampletype;
  }
  Configs* conf = new Configs(sampleyear, sampletype, ifile);
  conf->DASInput = false;
  conf->bTagEffHistCreation = false;
  conf->AuxHistCreation = true;
  conf->UseMergedAuxHist = true;
  // conf->ProgressInterval = 1;
  conf->InputFile = infile;
  // conf->EntryMax = 1000;

  NanoAODReader* r = new NanoAODReader(conf);
  bTagEff* bTE = new bTagEff(conf);
  r->SetbTag(bTE);
  JetScale *JS = new JetScale(conf);
  GenHypothesis *gh = new GenHypothesis(conf->WPType);
  Permutations *PM = new Permutations(conf);

  Long64_t nentries = r->GetEntries();
  if (conf->EntryMax > 0 && conf->EntryMax < nentries) nentries = conf->EntryMax;
  Progress* progress = new Progress(nentries,conf->ProgressInterval);
  for (unsigned ievt = 0; ievt < nentries; ++ievt) {
    progress->Print(ievt);
    if (r->ReadEvent(ievt) < 0) continue;
    r->BranchSizeCheck();
    if (conf->EntryMax > 0 && ievt == conf->EntryMax) break;
    //check jet multiplicity
    int nj = r->RegionAssociations.GetNJets(0); // Nominal region
    int nb = r->RegionAssociations.GetbTags(0); // Nominal region
    for (unsigned ij = 0; ij < r->Jets.size(); ++ij) {
      if (r->Jets[ij].genJetIdx != -1 && ((nj >= 5 && nb >= 3) || (nj >= 6 && nb >= 3))) { // Fitter takes SR inputs
        JS->FillJet(r->Jets[ij],r->GenJets[r->Jets[ij].genJetIdx], r->EventWeights[0].first);
      }
    }
    if (conf->Type != 2) continue;

    if (nj != 5 && nj != 6) continue;
    gh->SetGenParts(r->GenParts);
    gh->MatchToJets(r->GenJets, r->Jets);
    if (!gh->ValidReco || !gh->ValidGen) continue;
    gh->MatchToLeptons(r->Leptons);
    gh->CreateHypothesisSet(r->TheLepton, r->Met);
    Hypothesis Tar = gh->TarRecoHypo;
    vector<TLorentzVector> Neus;
    if (JS->SolveNeutrinos(Tar.Lep, Tar.Neu, Neus) == 1) {
      if (fabs(172.186 - (Tar.Lep + Neus[0] + Tar.Jets[3]).M()) < fabs(172.186 - (Tar.Lep + Neus[1] + Tar.Jets[3]).M())) Tar.Neu = Neus[0];
      else Tar.Neu = Neus[1];
    }
    else Tar.Neu = TLorentzVector();
    PM->FillHypo(Tar, gh->OutJets, r->EventWeights[0].first);
    JS->FillHypo(Tar, r->EventWeights[0].first);
  }
  JS->PostProcess();
  JS->Clear();
  PM->PostProcess();
  PM->Clear();
  gh->Summarize();

  progress->JobEnd();

  gApplication->Terminate();

}