WPrimeAnalysis_obs/MatchAll.cc at master · sifuluo/WPrimeAnalysis_obs · 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
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
#include "Utilities/Analyzer.cc"
#include "Utilities/JESTools.cc"
#include "Utilities/GenTools.cc"
#include "Utilities/ROOTMini.cc"

#include <TROOT.h>
#include <TClonesArray.h>
#include <TLorentzVector.h>

#include <TH1.h>
#include <TH2.h>
#include <TProfile.h>
#include <TProfile2D.h>
#include <TEfficiency.h>
#include <TString.h>
#include <TTree.h>

#include <vector>
#include <utility>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <cmath>
#include <string>
#include <map>
#include <algorithm>

using namespace std;

void MatchAll(int SampleType = 0, int irun = 1, int OptionCode = 0, int debug = -1) {
  cout << "start" << endl;
  Analyzer *a = new Analyzer(SampleType, irun, 30);
  TString savepath = "results/";
  TString savename = "MatchAll";
  double dRToMatch = 0.4;
  cout << "Running " << savename << endl;

  a->SetupROOTMini();
  a->SetOutput(savepath, savename);
  a->DebugMode(debug);
  a->CDOut();
  // Histograms

  vector<int>* Conflict = new vector<int>;
  vector<double>* DeltaR = new vector<double>;
  vector<double>* PtRatio = new vector<double>;
  cout << "before init" <<endl;
  a->Tree_Init(); // Init tree and initialize gen and reco branch
  cout << "done init" <<endl;
  a->t->Branch("Conflict",&Conflict);
  a->t->Branch("DeltaR",&DeltaR);
  a->t->Branch("PtRatio",&PtRatio);

  for (Int_t entry = a->GetStartEntry(); entry < a->GetEndEntry(); ++entry) {
    a->ReadEvent(entry);
    if (a->AssignGenParticles() == -1) continue;
    if (a->RecoPass == -1) continue;
    a->Tree_Reco();
    Conflict->clear();
    DeltaR->clear();
    PtRatio->clear();
    a->CountEvent("PassPreSelection");
    vector<TLorentzVector> Jets = a->LVJets;
    vector<TLorentzVector> GenJets = a->LVGenOutSort; // LF0, LF1, HadB, LepB, W' B
    vector<int> tmpvec;
    for (unsigned igen = 0; igen < GenJets.size(); ++igen) {
      map<double,int> matchcontainer; // Sort by deltaR before put in vector
      TLorentzVector lvgen = GenJets.at(igen);
      if (lvgen.Pt() == 0) {
        tmpvec.push_back(-2);
        continue;
      }
      for (unsigned ij = 0; ij < Jets.size(); ++ij) {
        TLorentzVector lvreco = Jets.at(ij);
        double dr = lvreco.DeltaR(lvgen);
        if (dr < dRToMatch) {
          matchcontainer[dr] = ij;
        }
      }
      for (const auto &it: matchcontainer) {
        tmpvec.push_back(it.second);
      }
      tmpvec.push_back(-2);
    }

    int indexgen = 0;
    int indexconf = 1;
    Conflict->resize(tmpvec.size(),0);
    for (unsigned iv = 0; iv < tmpvec.size(); ++iv) {
      bool has_conf = false;
      int indexreco = tmpvec[iv];
      if (indexreco == -2) {
        Conflict->at(iv) = -2;
        DeltaR->push_back(-2);
        PtRatio->push_back(-2);
        indexgen++;
        continue;
      }
      for (unsigned iv2 = iv + 1; iv2 < tmpvec.size(); ++iv2) {
        int indexreco2 = tmpvec[iv2];
        if (indexreco2 == indexreco) {
          has_conf = true;
          Conflict->at(iv2) = indexconf;
        }
      }
      if (has_conf) {
        Conflict->at(iv) = indexconf;
        indexconf++;
      }
      TLorentzVector jet = Jets.at(indexreco);
      TLorentzVector gen = GenJets.at(indexgen);
      DeltaR->push_back(jet.DeltaR(gen));
      PtRatio->push_back(jet.Pt() / gen.Pt());
    }
    // cout << Form("vector sizes: tmpvec = %d, Conflict = %d, DeltaR = %d",tmpvec.size(),Conflict->size(), DeltaR->size())<<endl;
    a->Tree_Fill();
  }

  a->Tree_Save();
  a->SaveOutput();


}