1 | /*
|
---|
2 | Cut Functions: Return 1 when cut is passed. Return 0 when cut is failed.
|
---|
3 | */
|
---|
4 |
|
---|
5 | #ifdef __CLING__
|
---|
6 | #include <map>
|
---|
7 | #include <vector>
|
---|
8 | R__LOAD_LIBRARY(libDelphes)
|
---|
9 | #include "classes/DelphesClasses.h"
|
---|
10 | #include "external/ExRootAnalysis/ExRootTreeReader.h"
|
---|
11 | #endif
|
---|
12 |
|
---|
13 | #include "external/fastjet/PseudoJet.hh"
|
---|
14 | #include "external/fastjet/JetDefinition.hh"
|
---|
15 | #include "external/fastjet/ClusterSequence.hh"
|
---|
16 |
|
---|
17 | namespace using fastjet;
|
---|
18 |
|
---|
19 |
|
---|
20 |
|
---|
21 | //------------------------------------------------------------------------------
|
---|
22 |
|
---|
23 |
|
---|
24 |
|
---|
25 | void jetMass_s(void){
|
---|
26 |
|
---|
27 | gSystem->Load("libDelphes");
|
---|
28 |
|
---|
29 | TChain chain("Delphes");
|
---|
30 |
|
---|
31 | Processes smbkg;
|
---|
32 |
|
---|
33 | const int file_error_tolerence = 2;
|
---|
34 | char inputFile[150];
|
---|
35 | char outputFile[150];
|
---|
36 |
|
---|
37 | const double luminosity = 2000.;
|
---|
38 | double weight = 0.;
|
---|
39 |
|
---|
40 | int mgluino = 2000;
|
---|
41 | char inputFile[150];
|
---|
42 | char outputFile[150];
|
---|
43 |
|
---|
44 |
|
---|
45 | for(int msquark = 3000; msquark < 9000; msquark = msquark + 2500){
|
---|
46 | sprintf(inputFile, "/media/john/EC7A174B7A1711C6/Linux/Stops100TeV/%d_%d/001.root", msquark, mgluino);
|
---|
47 | sprintf(outputFile, "OutputGridFiles/JetMass/point_%d_mstop_%d_mgluino_analysis.root", msquark, mgluino);
|
---|
48 |
|
---|
49 | TFile* outFile = new TFile(outputFile, "RECREATE");
|
---|
50 |
|
---|
51 | TH1F *scalarHT = new TH1F("scalarHT","Scalar Sum HT", 100, .0, 50000.);
|
---|
52 | TH1F *missingET = new TH1F("missingET","MET", 100, .0, 10000.);
|
---|
53 | TH1F *numberOfJet = new TH1F("numberOfJet","Number Of Jet", 20, .0, 20.);
|
---|
54 | TH1F *massjet[6];
|
---|
55 | for(int i = 0; i < 6; i++){
|
---|
56 | char name[20];
|
---|
57 | char title[100];
|
---|
58 | sprintf(name, "massjet%i", i+1);
|
---|
59 | sprintf(title, "Jet %i Mass", i+1);
|
---|
60 | massjet[i] = new TH1F(name,title, 100, .0, 2000.0);
|
---|
61 | }
|
---|
62 | /*
|
---|
63 | double Rparam = 0.5;
|
---|
64 | fastjet::Strategy strategy = fastjet::Best;
|
---|
65 | fastjet::RecombinationScheme recombScheme = fastjet::E_scheme;
|
---|
66 | fastjet::JetDefinition *jetDef = NULL;
|
---|
67 | jetDef = new fastjet::JetDefinition(fastjet::antikt_algorithm, Rparam,
|
---|
68 | recombScheme, strategy);
|
---|
69 |
|
---|
70 | // Fastjet input
|
---|
71 | std::vector <fastjet::PseudoJet> fjInputs;
|
---|
72 |
|
---|
73 |
|
---|
74 | fjInputs.resize(0);
|
---|
75 | */
|
---|
76 |
|
---|
77 | for(auto iterator = smbkg.Stops.begin(); iterator != smbkg.Stops.end(); iterator++) {
|
---|
78 | int msquark_s = atoi(iterator->first);
|
---|
79 | if(msquark_s == msquark){
|
---|
80 | weight = 1000.*luminosity*iterator->second;
|
---|
81 | break;
|
---|
82 | }
|
---|
83 | }
|
---|
84 | chain.Add(inputFile);
|
---|
85 |
|
---|
86 | ExRootTreeReader *treeReader = new ExRootTreeReader(&chain);
|
---|
87 | Long64_t numberOfEntries = treeReader->GetEntries();
|
---|
88 |
|
---|
89 | branchEvent = treeReader->UseBranch("Event");
|
---|
90 | branchJet = treeReader->UseBranch("Jet");
|
---|
91 | branchEFlowMuon = treeReader->UseBranch("EFlowMuon");
|
---|
92 | branchMuon = treeReader->UseBranch("Muon");
|
---|
93 | branchElectron = treeReader->UseBranch("Electron");
|
---|
94 | branchParticle = treeReader->UseBranch("Particle");
|
---|
95 | branchMissingET = treeReader->UseBranch("MissingET");
|
---|
96 | branchScalarHT = treeReader->UseBranch("ScalarHT");
|
---|
97 |
|
---|
98 | double n_entries = numberOfEntries;
|
---|
99 | if(n_entries < 1) {
|
---|
100 | cout << inputFile << " contains insufficient entries!" << endl;
|
---|
101 | continue;
|
---|
102 | }
|
---|
103 | weight = weight/n_entries;
|
---|
104 | cout << "Reading: " << inputFile << " with bin weight " << weight << endl;
|
---|
105 |
|
---|
106 | const double inclusive_jet_pt = 0.0;
|
---|
107 | const double PI = 3.14159265359;
|
---|
108 | const double jetR = 0.5;
|
---|
109 |
|
---|
110 | const double scalarHTcut = 3500;
|
---|
111 | const double missingETcut = 500;
|
---|
112 |
|
---|
113 |
|
---|
114 | // Loop over all events
|
---|
115 | for(Int_t entry = 0; entry < numberOfEntries; ++entry){
|
---|
116 | treeReader->ReadEntry(entry);
|
---|
117 |
|
---|
118 | if(CutScalarHT(branchScalarHT, scalarHTcut) < 1) continue;
|
---|
119 | if(CohenMET(branchMissingET, missingETcut) < 1) continue;
|
---|
120 | if(branchJet->GetEntriesFast() < 6) continue;
|
---|
121 |
|
---|
122 | MissingET *met = (MissingET*) branchMissingET->At(0);
|
---|
123 | ScalarHT *sht = (ScalarHT*) branchScalarHT->At(0);
|
---|
124 |
|
---|
125 | numberOfJet->Fill(branchJet->GetEntriesFast(),weight);
|
---|
126 | scalarHT->Fill(sht->HT,weight);
|
---|
127 |
|
---|
128 | missingET->Fill(met->MET,weight);
|
---|
129 |
|
---|
130 |
|
---|
131 | for(int njet = 0; njet < 6; ++njet){
|
---|
132 | Jet* jet = branchJet->At(njet);
|
---|
133 |
|
---|
134 | massjet[njet]->Fill(jet->Mass,weight);
|
---|
135 | }
|
---|
136 | }
|
---|
137 |
|
---|
138 |
|
---|
139 | outFile->Write();
|
---|
140 | outFile->Close();
|
---|
141 |
|
---|
142 | cout << endl << "Output File - " << outputFile << endl;
|
---|
143 |
|
---|
144 | delete outFile;
|
---|
145 | }
|
---|
146 |
|
---|
147 | }
|
---|
148 |
|
---|
149 | //-------------------------------------------------------------------------------------------------------------------
|
---|
150 | int CohenJetPT(TClonesArray *br_Jet, const int njetreq = 2, const double ptJetCut = 1000., const double etaJetCut = 2.5){
|
---|
151 |
|
---|
152 | int nJetsEvent = br_Jet->GetEntriesFast();
|
---|
153 | int jetsPassed = 0;
|
---|
154 |
|
---|
155 | if(nJetsEvent < njetreq) return 0;
|
---|
156 |
|
---|
157 | Jet* jet;
|
---|
158 | for(int i=0; i < nJetsEvent; i++){
|
---|
159 | jet = (Jet*) br_Jet->At(i);
|
---|
160 | if(fabs(jet->Eta) > etaJetCut) continue;
|
---|
161 | if(jet->PT > ptJetCut) jetsPassed += 1;
|
---|
162 | if(jetsPassed == njetreq) return 1;
|
---|
163 | }
|
---|
164 | return 0;
|
---|
165 | }
|
---|
166 |
|
---|
167 | //-------------------------------------------------------------------------------------------------------------------
|
---|
168 | int CohenMuonInJet(TClonesArray *br_Jet, TClonesArray *br_EFlowMuon, const double ptMuInJetCut = 1000., const int nJetMuons = 1, const int nLeadJets = 2, const double jetR = 0.5, const double etaCut = 2.5){
|
---|
169 |
|
---|
170 | const double PI = 3.14159265359;
|
---|
171 | int muonInJet = 0;
|
---|
172 | int nMuonEntries = br_EFlowMuon->GetEntriesFast();
|
---|
173 | int jetsLoopSize = nLeadJets;
|
---|
174 |
|
---|
175 | Muon *muon;
|
---|
176 | Jet *jet;
|
---|
177 |
|
---|
178 | if(br_Jet->GetEntriesFast() < nLeadJets) jetsLoopSize = br_Jet->GetEntriesFast();
|
---|
179 |
|
---|
180 | for(int i=0; i < jetsLoopSize; i++) {
|
---|
181 | jet = (Jet*) br_Jet->At(i);
|
---|
182 | if (fabs(jet->Eta) > etaCut) continue;
|
---|
183 |
|
---|
184 | for (int j = 0; j < nMuonEntries; ++j) {
|
---|
185 | muon = (Muon *) br_EFlowMuon->At(j);
|
---|
186 | if(muon->PT > ptMuInJetCut && fabs(muon->Eta) < etaCut) {
|
---|
187 | double dphi = fabs(muon->Phi - jet->Phi);
|
---|
188 | if (dphi > PI) dphi = 2*PI - dphi;
|
---|
189 | double deta = fabs(muon->Eta - jet->Eta);
|
---|
190 | if (sqrt(dphi*dphi + deta*deta) < jetR) muonInJet += 1;
|
---|
191 | if (muonInJet == nJetMuons) return 1;
|
---|
192 | }
|
---|
193 | }
|
---|
194 | }
|
---|
195 | return 0;
|
---|
196 | }
|
---|
197 |
|
---|
198 | //-------------------------------------------------------------------------------------------------------------------
|
---|
199 | int CohenNoIsoLepton(TClonesArray *br_Electron, TClonesArray *br_Muon, const double ptLepton = 35., const double ptIsolationRatio = 0.1, const int isoLTolerance = 0, const double etaCut = 2.5){
|
---|
200 |
|
---|
201 | int nElectronEntries = br_Electron->GetEntriesFast();
|
---|
202 | int nMuonEntries = br_Muon->GetEntriesFast();
|
---|
203 | int nIsoLeptons = 0;
|
---|
204 | int counter = 0;
|
---|
205 | Electron *electron;
|
---|
206 | Muon *muon;
|
---|
207 |
|
---|
208 | for (int i = 0; i < nElectronEntries; ++i) {
|
---|
209 | electron = (Electron *) br_Electron->At(i);
|
---|
210 | if (electron->PT > ptLepton && fabs(electron->Eta) < etaCut) nIsoLeptons += 1 ;
|
---|
211 | ++counter;
|
---|
212 | }
|
---|
213 | for (int i = 0; i < nMuonEntries; ++i) {
|
---|
214 | muon = (Muon *) br_Muon->At(i);
|
---|
215 | if (muon->PT > ptLepton && fabs(muon->Eta) < etaCut) nIsoLeptons += 1 ;
|
---|
216 | ++counter;
|
---|
217 | }
|
---|
218 |
|
---|
219 | if (nIsoLeptons > isoLTolerance) return 0;
|
---|
220 | if (counter == nMuonEntries + nElectronEntries) return 1;
|
---|
221 | }
|
---|
222 |
|
---|
223 | //-------------------------------------------------------------------------------------------------------------------
|
---|
224 | int CohenPhiIsoMET(TClonesArray *br_MissingET, TClonesArray *br_Jet, const double phiCut = 1.0, const double JetpTCut = 200., const double etaCut = 2.5){
|
---|
225 |
|
---|
226 | const double PI = 3.14159265359;
|
---|
227 | double jetdPhi = 0.;
|
---|
228 | int nJetsEvent = br_Jet->GetEntriesFast();
|
---|
229 |
|
---|
230 | MissingET *met = (MissingET*) br_MissingET->At(0);
|
---|
231 | Jet *jet;
|
---|
232 |
|
---|
233 | for (int i = 0; i < nJetsEvent; ++i){
|
---|
234 | jet = (Jet*) br_Jet->At(i);
|
---|
235 | if (jet->PT < JetpTCut || fabs(jet->Eta) > etaCut) continue;
|
---|
236 | jetdPhi = fabs(jet->Phi - met->Phi);
|
---|
237 | if (jetdPhi > PI) jetdPhi = 2*PI - jetdPhi;
|
---|
238 | if (jetdPhi < phiCut) return 0;
|
---|
239 | if (i == nJetsEvent) return 1;
|
---|
240 | }
|
---|
241 | }
|
---|
242 |
|
---|
243 | //-------------------------------------------------------------------------------------------------------------------
|
---|
244 |
|
---|
245 | int CohenMET(TClonesArray *br_MissingET, const double eMetCut = 3000.0){
|
---|
246 |
|
---|
247 | MissingET *met = (MissingET*) br_MissingET->At(0);
|
---|
248 | if (met->MET > eMetCut) return 1;
|
---|
249 | return 0;
|
---|
250 | }
|
---|
251 | //-------------------------------------------------------------------------------------------------------------------
|
---|
252 |
|
---|
253 | int CutScalarHT(TClonesArray *br_ScalarHT, const double scalarHTCut = 2000.0){
|
---|
254 |
|
---|
255 | ScalarHT *sht = (ScalarHT*) br_ScalarHT->At(0);
|
---|
256 | if (sht->HT > scalarHTCut) return 1;
|
---|
257 | return 0;
|
---|
258 | }
|
---|
259 |
|
---|
260 | //-------------------------------------------------------------------------------------------------------------------
|
---|
261 |
|
---|
262 |
|
---|
263 | int ZJetDistribution(TClonesArray *br_Jet, const double ZVarCut = 0.8, const double jetPtCut = 200., const int minJetsEvent = 4, const int leadingJets = 2) {
|
---|
264 |
|
---|
265 | int jetsLoopSize = minJetsEvent;
|
---|
266 | double z = -1.;
|
---|
267 | double x = 0. ,y = 0.;
|
---|
268 |
|
---|
269 | if(br_Jet->GetEntriesFast() < minJetsEvent) jetsLoopSize = br_Jet->GetEntriesFast();
|
---|
270 | Jet *lastjet = (Jet*) br_Jet->At(jetsLoopSize-1);
|
---|
271 | if(lastjet->PT < jetPtCut) return 0;
|
---|
272 |
|
---|
273 | for(int i=0; i < jetsLoopSize; i++){
|
---|
274 | Jet *jet = (Jet*) br_Jet->At(i);
|
---|
275 | if (jet->PT > jetPtCut && i < leadingJets) y += jet->PT;
|
---|
276 | else if(jet->PT > jetPtCut) x += jet->PT;
|
---|
277 | }
|
---|
278 | z = 2.*x/y;
|
---|
279 | if(z > ZVarCut) return 1;
|
---|
280 | return 0;
|
---|
281 | }
|
---|
282 |
|
---|
283 | //------------------------------------------------------------------------------------------------------------------------
|
---|
284 | class Processes
|
---|
285 | {
|
---|
286 | public:
|
---|
287 | Processes();
|
---|
288 | ~Processes();
|
---|
289 | XsecOf(const char *, const char *);
|
---|
290 | Processes::BinsOf(const char, std::map <char *, double> );
|
---|
291 |
|
---|
292 | std::vector <char *> process;
|
---|
293 | std::vector <char *> process_sm;
|
---|
294 | std::vector <char *> process_susy;
|
---|
295 | std::map <char *, double> ttB;
|
---|
296 | std::map <char *, double> tt;
|
---|
297 | std::map <char *, double> tB;
|
---|
298 | std::map <char *, double> Stops;
|
---|
299 | };
|
---|
300 |
|
---|
301 | void Processes::Processes()
|
---|
302 | {
|
---|
303 | process_sm.push_back("ttB");
|
---|
304 | process_sm.push_back("tt");
|
---|
305 | process_sm.push_back("tB");
|
---|
306 | process_susy.push_back("Stops");
|
---|
307 |
|
---|
308 | process.insert(process.end(), process_sm.begin(), process_sm.end());
|
---|
309 | process.insert(process.end(), process_susy.begin(), process_susy.end());
|
---|
310 |
|
---|
311 | ttB["0-1500"] = 206.01;
|
---|
312 | ttB["1500-3000"] = 12.58;
|
---|
313 | ttB["3000-5500"] = 1.18;
|
---|
314 | ttB["5500-9000"] = 0.092;
|
---|
315 | ttB["9000-100000"] = 0.009;
|
---|
316 |
|
---|
317 | tt["0-1000"] = 29141.30;
|
---|
318 | tt["1000-2000"] = 1777.28;
|
---|
319 | tt["2000-3500"] = 185.22;
|
---|
320 | tt["3500-5500"] = 18.92;
|
---|
321 | tt["5500-8500"] = 2.39;
|
---|
322 | tt["8500-100000"] = 0.277;
|
---|
323 |
|
---|
324 | tB["0-1000"] = 3399.65;
|
---|
325 | tB["1000-2000"] = 165.25;
|
---|
326 | tB["2000-3500"] = 15.57;
|
---|
327 | tB["3500-6000"] = 1.59;
|
---|
328 | tB["6000-9000"] = 0.11;
|
---|
329 | tB["9000-100000"] = 0.013;
|
---|
330 |
|
---|
331 | Stops["1500"] = 0.535;
|
---|
332 | Stops["2000"] = 0.123;
|
---|
333 | Stops["2500"] = 0.038;
|
---|
334 | Stops["3000"] = 0.0142;
|
---|
335 | Stops["3500"] = 0.00593;
|
---|
336 | Stops["4000"] = 0.00276;
|
---|
337 | Stops["4500"] = 0.00136;
|
---|
338 | Stops["5000"] = 0.00071;
|
---|
339 | Stops["5500"] = 0.00038;
|
---|
340 | Stops["6000"] = 0.00022;
|
---|
341 | Stops["6500"] = 0.00013;
|
---|
342 | Stops["7000"] = 0.000077;
|
---|
343 | Stops["7500"] = 0.000047;
|
---|
344 | Stops["8000"] = 0.000029;
|
---|
345 | Stops["8500"] = 0.000019;
|
---|
346 | Stops["9000"] = 0.000012;
|
---|
347 |
|
---|
348 | }
|
---|
349 |
|
---|
350 | void Processes::~Processes() {
|
---|
351 | }
|
---|
352 |
|
---|
353 | double Processes::XsecOf(const char * proc, const char * bin){
|
---|
354 | if(proc == "ttB") return Processes::ttB[bin];
|
---|
355 | if(proc == "tt") return Processes::tt[bin];
|
---|
356 | if(proc == "tB") return Processes::tB[bin];
|
---|
357 | if(proc == "Stops") return Processes::Stops[bin];
|
---|
358 | return -1.;
|
---|
359 | }
|
---|
360 |
|
---|
361 | int Processes::BinsOf(const char * proc, std::map <char *, double>* bins){
|
---|
362 | if(proc == "ttB") &bins = Processes::ttB; return 0;
|
---|
363 | if(proc == "tt") &bins = Processes::tt; return 0;
|
---|
364 | if(proc == "tB") &bins = Processes::tB; return 0;
|
---|
365 | if(proc == "Stops") &bins = Processes::Stops; return 0;
|
---|
366 | return -1.;
|
---|
367 | }
|
---|
368 |
|
---|
369 |
|
---|
370 | //-------------------------------------------------------------------------------------------------------------------------------------------
|
---|
371 |
|
---|
372 | bool check_File (const char * file) {
|
---|
373 | ifstream f(file);
|
---|
374 | if (f.good()) {
|
---|
375 | f.close();
|
---|
376 | return true;
|
---|
377 | } else {
|
---|
378 | f.close();
|
---|
379 | return false;
|
---|
380 | }
|
---|
381 | }
|
---|
382 |
|
---|
383 |
|
---|
384 |
|
---|
385 |
|
---|
386 |
|
---|