MA5SandBox: cms_top_18_003.cpp

File cms_top_18_003.cpp, 12.4 KB (added by Benjamin Fuks, 6 years ago)
Line 
1#include "SampleAnalyzer/User/Analyzer/cms_top_18_003.h"
2using namespace MA5;
3
4// -----------------------------------------------------------------------------
5// Initialize
6// function called one time at the beginning of the analysis
7// -----------------------------------------------------------------------------
8bool cms_top_18_003::Initialize(const MA5::Configuration& cfg, const std::map<std::string,std::string>& parameters)
9{
10 // Information on the analysis, authors, ...
11 INFO << " <><><><><><><><><><><><><><><><><><><><><><><><><>" << endmsg;
12 INFO << " <> Analysis: CMS-TOP-18-003, arXiv:xxxxxxxx <>" << endmsg;
13 INFO << " <> (four top - SS2L & multilepton) <>" << endmsg;
14 INFO << " <> Recaster: B. Fuks, L. Darme <>" << endmsg;
15 INFO << " <> Contact: fuks@lpthe.jussieu.fr <>" << endmsg;
16 INFO << " <> Based on MadAnalysis 5 v1.6 <>" << endmsg;
17 INFO << " <> DOI: XX.YYYY/ZZZ <>" << endmsg;
18 INFO << " <><><><><><><><><><><><><><><><><><><><><><><><><>" << endmsg;
19
20 // Declaration of the signal regions
21 Manager()->AddRegionSelection("SR1");
22 Manager()->AddRegionSelection("SR2");
23 Manager()->AddRegionSelection("SR3");
24 Manager()->AddRegionSelection("SR4");
25 Manager()->AddRegionSelection("SR5");
26 Manager()->AddRegionSelection("SR6");
27 Manager()->AddRegionSelection("SR7");
28 Manager()->AddRegionSelection("SR8");
29 Manager()->AddRegionSelection("SR9");
30 Manager()->AddRegionSelection("SR10");
31 Manager()->AddRegionSelection("SR11");
32 Manager()->AddRegionSelection("SR12");
33 Manager()->AddRegionSelection("SR13");
34 Manager()->AddRegionSelection("SR14");
35
36 // Declaration of the preselection cuts
37 Manager()->AddCut("baseline");
38 Manager()->AddCut("mee > 12 GeV");
39 Manager()->AddCut("veto_3l");
40 Manager()->AddCut("InSR");
41
42 // Histograms
43 Manager()->AddHisto("njets", 11,0.5,11.5);
44 Manager()->AddHisto("nb", 6,0.5, 6.5);
45 Manager()->AddHisto("HT", 16,0.,1600.);
46 Manager()->AddHisto("MET", 15,0., 600.);
47
48 // Declaration of the selection cuts
49 std::string TwoLeps[] = { "SR1", "SR2", "SR3", "SR4", "SR5", "SR6","SR7","SR8"};
50 std::string MoreLeps[] = { "SR9", "SR10", "SR11", "SR12", "SR13", "SR14"};
51 Manager()->AddCut("2 leptons", TwoLeps);
52 Manager()->AddCut("more leptons", MoreLeps);
53
54 std::string TwoBs[] = { "SR1", "SR2", "SR3", "SR9", "SR10", "SR11"};
55 std::string ThreeBs[] = { "SR4", "SR5", "SR6", "SR7"};
56 std::string MoreThreeBs[] = {"SR12", "SR13", "SR14"};
57 Manager()->AddCut("2 bjets", TwoBs);
58 Manager()->AddCut("3 bjets", ThreeBs);
59 Manager()->AddCut("more4 bjets", "SR8");
60 Manager()->AddCut("more3 bjets", MoreThreeBs);
61
62 std::string More8Js[] = { "SR3", "SR7"};
63 std::string FiveJs[] = { "SR4", "SR9", "SR13"};
64 std::string SixJs[] = { "SR1", "SR5", "SR10"};
65 std::string SevenJs[] = { "SR2", "SR6"};
66 Manager()->AddCut("4 jets","SR12");
67 Manager()->AddCut("5 jets",FiveJs);
68 Manager()->AddCut("6 jets",SixJs);
69 Manager()->AddCut("7 jets", SevenJs);
70 Manager()->AddCut("more5 jets", "SR8");
71 Manager()->AddCut("more6 jets", "SR14");
72 Manager()->AddCut("more7 jets", "SR11");
73 Manager()->AddCut("more8 jets", More8Js);
74
75 return true;
76}
77
78// -----------------------------------------------------------------------------
79// Finalize
80// function called one time at the end of the analysis
81// -----------------------------------------------------------------------------
82void cms_top_18_003::Finalize(const SampleFormat& summary, const std::vector<SampleFormat>& files)
83{
84}
85
86// -----------------------------------------------------------------------------
87// Execute
88// function called each time one event is read
89// -----------------------------------------------------------------------------
90bool cms_top_18_003::Execute(SampleFormat& sample, const EventFormat& event)
91{
92 // Event weight
93 double myWeight;
94 if(Configuration().IsNoEventWeight()) myWeight=1.;
95 else if(event.mc()->weight()!=0.) myWeight=event.mc()->weight();
96 else
97 {
98 WARNING << "Found one event with a zero weight. Skipping...\n";
99 return false;
100 }
101 Manager()->InitializeForNewEvent(myWeight);
102
103 // Security for empty events
104 if (event.rec()==0) return true;
105
106 // Electrons
107 std::vector<const RecLeptonFormat*> LooseElectrons;
108 for(unsigned int ii=0; ii<event.rec()->electrons().size(); ii++)
109 {
110 const RecLeptonFormat *myElec = &(event.rec()->electrons()[ii]);
111 if(IsLooseLepton(myElec, event, 5., 2.5))
112 LooseElectrons.push_back(myElec);
113 }
114
115 // Muons
116 std::vector<const RecLeptonFormat*> LooseMuons;
117 for(unsigned int ii=0; ii<event.rec()->muons().size(); ii++)
118 {
119 const RecLeptonFormat *myMuon = &(event.rec()->muons()[ii]);
120 if(IsLooseLepton(myMuon, event, 7., 2.4))
121 LooseMuons.push_back(myMuon);
122 }
123
124 // Jets
125 std::vector <const RecJetFormat*> IsoJets, SignalJets;
126 unsigned int nj=0, nb=0;
127 double HT=0.;
128 for(unsigned int ii=0; ii<event.rec()->jets().size(); ii++)
129 {
130 const RecJetFormat * myJet = &(event.rec()->jets()[ii]);
131 double eta = std::fabs(myJet->eta());
132 double pt = myJet->pt();
133 double iso = true, isL = false;
134 for (unsigned int jj=0; jj<LooseElectrons.size(); jj++)
135 {
136 if(myJet->dr(LooseElectrons[jj])<0.1) { isL = true; break; }
137 if(myJet->dr(LooseElectrons[jj])<0.4) { iso = false; break; }
138 }
139 if(iso && !isL)
140 {
141 for (unsigned int jj=0; jj<LooseMuons.size(); jj++)
142 {
143 if(myJet->dr(LooseMuons[jj])<0.1) { isL = true; break; }
144 if(myJet->dr(LooseMuons[jj])<0.4) { iso = false; break; }
145 }
146 }
147 if(!iso && !isL /*&& pt>5.*/) IsoJets.push_back(myJet);
148 else if(!isL && eta<2.4 && iso)
149 {
150 if(pt>25. && myJet->btag()) nb++;
151 if(pt>40.) { nj++; HT+=myJet->pt(); SignalJets.push_back(myJet);}
152 }
153 }
154
155 // Lepton isolation
156 std::vector<const RecLeptonFormat*> IsoElectrons, IsoMuons;
157 std::vector<const RecLeptonFormat*> SignalLeptons, SignalElectrons, SignalMuons;
158 for (unsigned int jj=0; jj<LooseElectrons.size(); jj++)
159 if(IsIsolated(LooseElectrons[jj],IsoJets,event,0.12,0.80,7.2))
160 {
161 if(LooseElectrons[jj]->pt()>20.)
162 {
163 SignalLeptons.push_back(LooseElectrons[jj]);
164 SignalElectrons.push_back(LooseElectrons[jj]);
165 }
166 IsoElectrons.push_back(LooseElectrons[jj]);
167 }
168 for (unsigned int jj=0; jj<LooseMuons.size(); jj++)
169 if(IsIsolated(LooseMuons[jj],IsoJets,event,0.16,0.76,7.2))
170 {
171 if(LooseMuons[jj]->pt()>20.)
172 {
173 SignalLeptons.push_back(LooseMuons[jj]);
174 SignalMuons.push_back(LooseMuons[jj]);
175 }
176 IsoMuons.push_back(LooseMuons[jj]);
177 }
178 SORTER->sort(SignalLeptons, PTordering);
179 unsigned int nl = SignalLeptons.size();
180
181 // MET
182 double MET = event.rec()->MET().pt();
183
184 // Cut-0: baseline
185 bool base_jets = (nj>=2) && (nb>=2) && (HT>300.);
186 bool base_leps = (nl>=2) && (SignalLeptons[0]->pt()>25.) && (SignalLeptons[1]->pt()>20.); // LD -- Added Trail lepton
187 bool SS_leps = false;
188 unsigned int iSS=99;
189 for(unsigned int ii=1; ii<SignalLeptons.size(); ii++)
190 if(SignalLeptons[0]->charge()==SignalLeptons[ii]->charge())
191 {
192 if(iSS==99) {SS_leps=true; iSS=ii; }
193 else {SS_leps=false; break; }
194 }
195 if(!Manager()->ApplyCut( base_jets && base_leps && SS_leps && MET > 50.,"baseline")) return true;
196
197 // cut-1: Mee > 12 GeV
198 bool mee_12=false;
199 for(unsigned int ii=0; ii<SignalElectrons.size(); ii++)
200 for(unsigned int jj=ii+1; jj<SignalElectrons.size(); jj++)
201 {
202 if(SignalElectrons[ii]->charge()==SignalElectrons[jj]->charge())
203 {
204 double mee = (SignalElectrons[ii]->momentum()+SignalElectrons[jj]->momentum()).M();
205 if(mee<12.) {mee_12=true; break;}
206 }
207 }
208 if(!Manager()->ApplyCut(!mee_12,"mee > 12 GeV")) return true;
209
210 // cut-3: 3rd lepton
211 bool veto_3l=false;
212 for(unsigned int jj=0; jj<IsoMuons.size();jj++)
213 {
214 if(IsoMuons[jj]->dr(SignalLeptons[0])<0.1) continue;
215 if(IsoMuons[jj]->dr(SignalLeptons[iSS])<0.1) continue;
216 for(unsigned int ii=0;ii<SignalMuons.size();ii++)
217 {
218 if(SignalMuons[ii]->dr(SignalLeptons[0])>0.1 &&
219 SignalMuons[ii]->dr(SignalLeptons[iSS])>0.1) continue;
220 if(IsoMuons[jj]->charge()==SignalMuons[ii]->charge()) continue;
221 double mll = (SignalMuons[ii]->momentum()+IsoMuons[jj]->momentum()).M();
222 double ThirdmupT= SignalMuons[ii]->pt();
223 if((ThirdmupT >5) && (mll<12. || (mll>76. && mll<106.)) ) { veto_3l=true; break; } // LD -- Added pt limit
224 }
225 if (veto_3l) break;
226 }
227 if(!veto_3l)
228 {
229 for(unsigned int jj=0; jj<IsoElectrons.size();jj++)
230 {
231 if(IsoElectrons[jj]->dr(SignalLeptons[0])<0.1) continue;
232 if(IsoElectrons[jj]->dr(SignalLeptons[iSS])<0.1) continue;
233 for(unsigned int ii=0;ii<SignalElectrons.size();ii++)
234 {
235 if(SignalElectrons[ii]->dr(SignalLeptons[0])>0.1 &&
236 SignalElectrons[ii]->dr(SignalLeptons[iSS])>0.1) continue;
237 if(IsoElectrons[jj]->charge()==SignalElectrons[ii]->charge()) continue;
238 double mll = (SignalElectrons[ii]->momentum()+IsoElectrons[jj]->momentum()).M();
239 double ThirdepT= SignalElectrons[ii]->pt();
240 if((ThirdepT >7) && (mll<12. || (mll>76. && mll<106.)) ) { veto_3l=true; break; }// LD -- Added pt limit
241 }
242 if (veto_3l) break;
243 }
244 }
245 if(!Manager()->ApplyCut(!veto_3l,"veto_3l")) return true;
246 if(!Manager()->ApplyCut( (nl==2 && nb<4 && (nb+nj)>7) || (nl==2 && nb>3 && nj>4) ||
247 (nl>2 && nb==2 && nj>4) || (nl>2 && nb>2 && nj>3) ,"InSR")) return true;
248
249
250 // Crazy observable
251 // SORTER->sort(SignalJets, PTordering);
252 // MALorentzVector HTvec(0.,0.,0.,0.);
253 // for(unsigned ii=1;ii<SignalJets.size();ii++)
254 // HTvec=HTvec+SignalJets[ii]->momentum();
255 // double crazy = std::fabs(SignalJets[0]->phi() - HTvec.Phi());
256 // if(crazy>M_PI) crazy=2.*M_PI-crazy;
257
258 // Histograms
259 Manager()->FillHisto("njets",nj);
260 Manager()->FillHisto("nb" ,nb);
261 Manager()->FillHisto("HT" ,HT);
262 Manager()->FillHisto("MET" ,MET);
263 // Manager()->FillHisto("test1" ,crazy);
264
265
266
267 // SR definitions : leptons
268 if(!Manager()->ApplyCut(nl==2,"2 leptons")) return true;
269 if(!Manager()->ApplyCut(nl> 2,"more leptons")) return true;
270
271 // SR definitions : b-jets
272 if(!Manager()->ApplyCut(nb==2,"2 bjets")) return true;
273 if(!Manager()->ApplyCut(nb==3,"3 bjets")) return true;
274 if(!Manager()->ApplyCut(nb> 3,"more4 bjets")) return true;
275 if(!Manager()->ApplyCut(nb> 2,"more3 bjets")) return true;
276
277 // SR definitions : jets
278 if(!Manager()->ApplyCut(nj==4,"4 jets")) return true;
279 if(!Manager()->ApplyCut(nj==5,"5 jets")) return true;
280 if(!Manager()->ApplyCut(nj==6,"6 jets")) return true;
281 if(!Manager()->ApplyCut(nj==7,"7 jets")) return true;
282 if(!Manager()->ApplyCut(nj> 4,"more5 jets")) return true;
283 if(!Manager()->ApplyCut(nj> 5,"more6 jets")) return true;
284 if(!Manager()->ApplyCut(nj> 6,"more7 jets")) return true;
285 if(!Manager()->ApplyCut(nj> 7,"more8 jets")) return true;
286
287return true;
288}
289
290
291// -----------------------------------------------------------------------------
292// User-defined methods
293// -----------------------------------------------------------------------------
294bool cms_top_18_003::IsLooseLepton(const RecLeptonFormat* Lep, const EventFormat& event,
295 const double &ptTH, const double &etaTH)
296{
297 double eta = std::fabs(Lep->eta());
298 double pt = Lep->pt();
299
300 double DR = 10./std::min(std::max(pt,50.),200.);
301 double imini = (PHYSICS->Isol->eflow->sumIsolation(Lep,
302 event.rec(),DR,0.,IsolationEFlow::ALL_COMPONENTS))/pt;
303
304 if( eta<etaTH && pt>ptTH && imini<0.4) return true;
305 return false;
306}
307
308
309bool cms_top_18_003::IsIsolated(const RecLeptonFormat* Lep,
310 std::vector<const RecJetFormat*> IsoJets, const EventFormat& event,
311 const double &I1, const double &I2, const double &I3)
312{
313 double pt = Lep->pt();
314 // I1
315 double DR = 10./std::min(std::max(pt,50.),200.);
316 double imini = (PHYSICS->Isol->eflow->sumIsolation(Lep,
317 event.rec(),DR,0.,IsolationEFlow::ALL_COMPONENTS))/pt;
318 // Closest jet
319 unsigned int ijet=9999;
320 double drmin=1e+09;
321 for(unsigned int ii=0; ii<IsoJets.size(); ii++)
322 if(Lep->dr(IsoJets[ii])<drmin) { drmin=Lep->dr(IsoJets[ii]); ijet=ii; }
323 // I2
324 double ptratio=1.;
325 if (ijet!=9999) ptratio = pt/IsoJets[ijet]->pt();
326 // I3
327 double ptrel = 0.;
328 if (ijet!=9999)
329 {
330 MAVector3 axis(IsoJets[ijet]->momentum().Vect() - Lep->momentum().Vect());
331 ptrel = (axis.Cross(Lep->momentum().Vect())).Mag()/axis.Mag();
332 }
333 // output
334 if (imini<I1 && (ptratio>I2 || ptrel>I3)) return true;
335 return false;
336}