Fork me on GitHub

Ticket #923: jetMass_s.C

File jetMass_s.C, 11.5 KB (added by john, 8 years ago)
Line 
1/*
2Cut Functions: Return 1 when cut is passed. Return 0 when cut is failed.
3*/
4
5#ifdef __CLING__
6#include <map>
7#include <vector>
8R__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
17namespace using fastjet;
18
19
20
21//------------------------------------------------------------------------------
22
23
24
25void 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//-------------------------------------------------------------------------------------------------------------------
150int 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//-------------------------------------------------------------------------------------------------------------------
168int 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//-------------------------------------------------------------------------------------------------------------------
199int 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//-------------------------------------------------------------------------------------------------------------------
224int 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
245int 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
253int 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
263int 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//------------------------------------------------------------------------------------------------------------------------
284class Processes
285{
286public:
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
301void 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
350void Processes::~Processes() {
351}
352
353double 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
361int 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
372bool 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