MA5SandBox: atlas_susy_2018_06.cpp

File atlas_susy_2018_06.cpp, 18.5 KB (added by Benjamin Fuks, 5 years ago)
Line 
1#include "SampleAnalyzer/User/Analyzer/atlas_susy_2018_06.h"
2using namespace MA5;
3using namespace std;
4
5namespace atlas_susy_2018_06_func
6{
7template<typename T1,typename T2> std::vector<const T1*>
8Removal(std::vector<const T1*> &v1, std::vector<const T2*> &v2, const MAdouble64 &drmin);
9
10template<typename T1,typename T2> std::vector<const T1*>
11RemovalLeptons(std::vector<const T1*> &v1, std::vector<const T2*> &v2);
12
13vector<int> LeptonPick(vector<const RecLeptonFormat*> leptons);
14MAdouble64 GetHBoost(vector<const RecLeptonFormat*> leptons,const RecParticleFormat* missing);
15}
16
17
18
19// -----------------------------------------------------------------------------
20// Initialize
21// function called one time at the beginning of the analysis
22// -----------------------------------------------------------------------------
23bool atlas_susy_2018_06::Initialize(const MA5::Configuration& cfg, const std::map<std::string,std::string>& parameters)
24{
25 INFO << " <><><><><><><><><><><><><><><><><><><><><><><><><><><><>" << endmsg;
26 INFO << " <> Analysis: ATLAS-SUSY-2018-06 <>" << endmsg;
27 INFO << " <> arXiv:1912.08479 <>" << endmsg;
28 INFO << " <> (chargino-neutralino) <>" << endmsg;
29 INFO << " <> Recasters: J.W.Kim & J.H.Kim & H.Jang & T.G.Lee <>" << endmsg;
30 INFO << " <> Contact: adg9334@naver.com <>" << endmsg;
31 INFO << " <> Based on MadAnalysis 5 v1.7 <>" << endmsg;
32 INFO << " <> For more information, see <>" << endmsg;
33 INFO << " <> http://madanalysis.irmp.ucl.ac.be/wiki/ <>" << endmsg;
34 INFO << " <> /PhysicsAnalysisDatabase <>" << endmsg;
35 INFO << " <><><><><><><><><><><><><><><><><><><><><><><><><><><><>" << endmsg;
36
37
38 //Add Signal Region
39
40
41 //Signal Region for Cutflows
42 Manager()->AddRegionSelection("SR_low");
43 Manager()->AddRegionSelection("SR_ISR");
44
45 Manager()->AddCut("3Leptons"); // Exact three leptons
46 Manager()->AddCut("SFOS"); // Same flavor, opposite-charge pair
47 //Manager()->AddCut("Dilepton Trigger");
48 Manager()->AddCut("B-veto"); // b-jet veto
49 Manager()->AddCut("Three Lepton Mass"); // Invariant mass of three leptons system
50 Manager()->AddCut("low-LeptonPT > 60, 40, 30","SR_low");
51 Manager()->AddCut("ISR-LeptonPT > 25, 25, 20","SR_ISR");
52 Manager()->AddCut("Dilepton InvMass"); // Invariant mass of SFOS pair
53
54 Manager()->AddCut("low-Jet-veto","SR_low");
55 Manager()->AddCut("low-HBoost","SR_low");
56 Manager()->AddCut("low-PTsoft/(PTsoft+Meff)","SR_low");
57 Manager()->AddCut("low-Meff/HBoost","SR_low");
58
59 Manager()->AddCut("ISR-Njet > 0","SR_ISR");
60 Manager()->AddCut("ISR-Njet < 4","SR_ISR");
61 Manager()->AddCut("ISR-DeltaPhi(MET,Jets) > 2.0","SR_ISR");
62 Manager()->AddCut("ISR-R(MET,Jets)","SR_ISR");
63 Manager()->AddCut("ISR-JetsPT > 100","SR_ISR");
64 Manager()->AddCut("ISR-MET > 80","SR_ISR");
65 Manager()->AddCut("Transverse Mass");
66 Manager()->AddCut("ISR-PTsoft < 25","SR_ISR");
67
68//For diagrams. If you want to draw diagrams, uncomment the paragraph below and the last part of Execute and comment the paragraph above.
69/*
70 //Signal Region for diagrams
71 Manager()->AddRegionSelection("SR_low");
72 Manager()->AddRegionSelection("SR_ISR");
73 Manager()->AddRegionSelection("SR-low-MT");
74 Manager()->AddRegionSelection("SR-low-HBoost");
75 Manager()->AddRegionSelection("SR-low-R(Meff,HBoost)");
76 Manager()->AddRegionSelection("SR-low-R(PTsoft,(PTsoft+Meff))");
77
78 Manager()->AddRegionSelection("SR-ISR-MT");
79 Manager()->AddRegionSelection("SR-ISR-R(MET,Jets)");
80 Manager()->AddRegionSelection("SR-ISR-PTsoft");
81 Manager()->AddRegionSelection("SR-ISR-PTjets");
82
83 std::string SR_low[] = {"SR_low","SR-low-MT","SR-low-HBoost","SR-low-R(Meff,HBoost)","SR-low-R(PTsoft,(PTsoft+Meff))"};
84 std::string SR_low_HBoost[] = {"SR_low","SR-low-MT","SR-low-R(Meff,HBoost)","SR-low-R(PTsoft,(PTsoft+Meff))"};
85 std::string SR_low_R_Meff_HBoost[] = {"SR_low","SR-low-MT","SR-low-HBoost","SR-low-R(PTsoft,(PTsoft+Meff))"};
86 std::string SR_low_R_PTsoft_Meff[] = {"SR_low","SR-low-MT","SR-low-HBoost","SR-low-R(Meff,HBoost)"};
87
88 std::string SR_ISR[] = {"SR_ISR","SR-ISR-MT","SR-ISR-R(MET,Jets)","SR-ISR-PTsoft","SR-ISR-PTjets"};
89 std::string SR_ISR_R_MET_Jets[] = {"SR_ISR","SR-ISR-MT","SR-ISR-PTsoft","SR-ISR-PTjets"};
90 std::string SR_ISR_PTsoft[] = {"SR_ISR","SR-ISR-MT","SR-ISR-R(MET,Jets)","SR-ISR-PTjets"};
91 std::string SR_ISR_PTjets[] = {"SR_ISR","SR-ISR-MT","SR-ISR-R(MET,Jets)","SR-ISR-PTsoft"};
92
93 std::string SR_MT[] = {"SR_low","SR-low-HBoost","SR-low-R(Meff,HBoost)","SR-low-R(PTsoft,(PTsoft+Meff))"
94 ,"SR_ISR","SR-ISR-R(MET,Jets)","SR-ISR-PTsoft","SR-ISR-PTjets"};
95
96
97 Manager()->AddCut("3Leptons"); // Exact three leptons
98 Manager()->AddCut("SFOS"); // Same flavor, opposite-charge pair
99 //Manager()->AddCut("Dilepton Trigger");
100 Manager()->AddCut("B-veto"); // b-jet veto
101 Manager()->AddCut("Three Lepton Mass"); // Invariant mass of three leptons system
102 Manager()->AddCut("low-LeptonPT > 60, 40, 30",SR_low);
103 Manager()->AddCut("ISR-LeptonPT > 25, 25, 20",SR_ISR);
104 Manager()->AddCut("Dilepton InvMass"); // Invariant mass of SFOS pair
105
106 Manager()->AddCut("low-Jet-veto",SR_low);
107 Manager()->AddCut("low-HBoost",SR_low_HBoost);
108 Manager()->AddCut("low-PTsoft/(PTsoft+Meff)",SR_low_R_PTsoft_Meff);
109 Manager()->AddCut("low-Meff/HBoost",SR_low_R_Meff_HBoost);
110
111 Manager()->AddCut("ISR-Njet > 0",SR_ISR);
112 Manager()->AddCut("ISR-Njet < 4",SR_ISR);
113 Manager()->AddCut("ISR-DeltaPhi(MET,Jets) > 2.0",SR_ISR);
114 Manager()->AddCut("ISR-R(MET,Jets)",SR_ISR_R_MET_Jets);
115 Manager()->AddCut("ISR-JetsPT > 100",SR_ISR_PTjets);
116 Manager()->AddCut("ISR-MET > 80",SR_ISR);
117 Manager()->AddCut("Transverse Mass",SR_MT);
118 Manager()->AddCut("ISR-PTsoft < 25",SR_ISR_PTsoft);
119
120
121 Manager()->AddHisto("SR-low-MT",20,0,500,"SR-low-MT");
122 Manager()->AddHisto("SR-low-HBoost",20,200,700,"SR-low-HBoost");
123 Manager()->AddHisto("SR-low-R(Meff,HBoost)",15,0.3,1.05,"SR-low-R(Meff,HBoost)");
124 Manager()->AddHisto("SR-low-R(PTsoft,(PTsoft+Meff))",15,0,0.15,"SR-low-R(PTsoft,(PTsoft+Meff))");
125
126 Manager()->AddHisto("SR-ISR-MT",50,0,500,"SR-ISR-MT");
127 Manager()->AddHisto("SR-ISR-R(MET,Jets)",18,0.1,1,"SR-ISR-R(MET,Jets)");
128 Manager()->AddHisto("SR-ISR-PTsoft",20,0,100,"SR-ISR-PTsoft");
129 Manager()->AddHisto("SR-ISR-PTjets",18,0,500,"SR-ISR-PTjets");
130
131 //Auxiliery figure
132 Manager()->AddHisto("SR-low-PTlep1",17,60,400,"SR_low");
133 Manager()->AddHisto("SR-low-PTlep2",13,40,300,"SR_low");
134 Manager()->AddHisto("SR-low-PTlep3",12,30,150,"SR_low");
135 Manager()->AddHisto("SR-low-MET",20,0,500,"SR_low");
136 Manager()->AddHisto("SR-ISR-PTlep1",10,25,225,"SR_ISR");
137 Manager()->AddHisto("SR-ISR-PTlep2",15,25,175,"SR_ISR");
138 Manager()->AddHisto("SR-ISR-PTlep3",8,20,100,"SR_ISR");
139*/
140 return true;
141}
142
143// -----------------------------------------------------------------------------
144// Finalize
145// function called one time at the end of the analysis
146// -----------------------------------------------------------------------------
147void atlas_susy_2018_06::Finalize(const SampleFormat& summary, const std::vector<SampleFormat>& files)
148{
149}
150
151// -----------------------------------------------------------------------------
152// Execute
153// function called each time one event is read
154// -----------------------------------------------------------------------------
155bool atlas_susy_2018_06::Execute(SampleFormat& sample, const EventFormat& event)
156{
157 //Check Weight
158 double myWeight;
159 if(Configuration().IsNoEventWeight()) myWeight=1. ;
160 else if(event.mc()->weight()!=0.) myWeight=event.mc()->weight();
161 else { return false; }
162 Manager()->InitializeForNewEvent(myWeight);
163
164 if(event.rec()==0) return true;
165 //Objects Isolation
166 //Signal Electrons
167 vector<const RecLeptonFormat*> SignalElectrons;
168 for(unsigned int i=0; i<event.rec()->electrons().size();i++)
169 {
170 const RecLeptonFormat *Lep = &(event.rec()->electrons()[i]);
171
172 //Kinamatics
173 double eta = Lep->abseta();
174 double pt = Lep->pt();
175
176 double iso_dR = 0.2;
177 double iso_tracks = PHYSICS->Isol->eflow->relIsolation(Lep,event.rec(),iso_dR,0,IsolationEFlow::TRACK_COMPONENT);
178 double iso_all = PHYSICS->Isol->eflow->relIsolation(Lep,event.rec(),iso_dR,0,IsolationEFlow::ALL_COMPONENTS);
179 bool iso = (iso_tracks<0.06 && iso_all<0.06);
180 if(eta<2.47 && pt>10. && iso) { SignalElectrons.push_back(Lep);}
181 }
182 //Signal Muon
183 vector<const RecLeptonFormat*> SignalMuons;
184 for(unsigned int i=0; i<event.rec()->muons().size();i++)
185 {
186 const RecLeptonFormat *Lep = &(event.rec()->muons()[i]);
187
188 //Kinamatics
189 double eta = Lep->abseta();
190 double pt = Lep->pt();
191
192 double iso_dR = 0.2;
193 if(pt<=33) {iso_dR=0.3;}
194 else if(pt<50) {iso_dR=-0.0059*pt + 0.4941;}
195
196 double iso_tracks = PHYSICS->Isol->eflow->relIsolation(Lep,event.rec(),iso_dR,0,IsolationEFlow::TRACK_COMPONENT);
197 double iso_all = PHYSICS->Isol->eflow->relIsolation(Lep,event.rec(),0.2,0,IsolationEFlow::ALL_COMPONENTS);
198 bool iso = (iso_tracks<0.04 && iso_all<0.15);
199
200 if(eta<2.4 && pt>10. && iso){SignalMuons.push_back(Lep);}
201 }
202 //Signal Jet
203 vector<const RecJetFormat*> Jets;
204 for(unsigned int i=0; i<event.rec()->jets().size();i++)
205 {
206 const RecJetFormat *Jet = &(event.rec()->jets()[i]);
207 double eta = Jet->abseta();
208 double pt = Jet->pt();
209 if(pt>20 && eta<2.4) {Jets.push_back(Jet);}
210 }
211
212 // Overlap removal
213 SignalElectrons = atlas_susy_2018_06_func::Removal(SignalElectrons, SignalMuons,0.2);
214 Jets = PHYSICS->Isol->JetCleaning(Jets,SignalElectrons,0.2);
215 Jets = PHYSICS->Isol->JetCleaning(Jets,SignalMuons,0.2);
216 SignalElectrons = atlas_susy_2018_06_func::RemovalLeptons(SignalElectrons,Jets);
217 SignalMuons = atlas_susy_2018_06_func::RemovalLeptons(SignalMuons,Jets);
218
219
220
221 vector<const RecLeptonFormat*> SignalLeptons;
222 for(unsigned int i=0; i<SignalMuons.size();i++)
223 {SignalLeptons.push_back(SignalMuons[i]);}
224 for(unsigned int i=0; i<SignalElectrons.size();i++)
225 {SignalLeptons.push_back(SignalElectrons[i]);}
226 SORTER->sort(SignalLeptons,PTordering);
227
228 //Missing Vector
229 const RecParticleFormat *MET = &(event.rec()->MET());
230
231 //////////////////////////////////////////////////////////////////////
232 //Cut
233 //////////////////////////////////////////////////////////////////////
234
235 // Nleptons = 3
236 bool Nlepton = (SignalLeptons.size()==3);
237 if(!Manager()->ApplyCut(Nlepton,"3Leptons")) return true;
238
239 //SFOS
240 vector<int> PairLeptons = atlas_susy_2018_06_func::LeptonPick(SignalLeptons);
241 bool SFOS = (SignalLeptons[PairLeptons[0]]->isMuon()==SignalLeptons[PairLeptons[1]]->isMuon());
242 SFOS = SFOS && (SignalLeptons[PairLeptons[0]]->charge()!=SignalLeptons[PairLeptons[1]]->charge());
243 if(!Manager()->ApplyCut(SFOS,"SFOS")) return true;
244
245 //if(!Manager()->ApplyCut("Dilepton Trigger")); // Not Applied
246
247 //B-veto
248 bool Bveto = true;
249 for(unsigned int i=0;i<Jets.size();i++)
250 if(Jets[i]->btag()) {Bveto = false; break;}
251 if(!Manager()->ApplyCut(Bveto,"B-veto")) return true;
252
253 // Mlll
254 MALorentzVector ThreeLeptonsSystem;
255 ThreeLeptonsSystem = SignalLeptons[0]->momentum()+SignalLeptons[1]->momentum()+SignalLeptons[2]->momentum();
256 MAdouble64 mlll = ThreeLeptonsSystem.M();
257 if(!Manager()->ApplyCut(mlll>105.,"Three Lepton Mass")) return true;
258
259 //Leptons PT
260 bool LowLeptonsPT = SignalLeptons[0]->pt()>60. && SignalLeptons[1]->pt()>40. && SignalLeptons[2]->pt()>30.;
261 bool ISRLeptonsPT = SignalLeptons[0]->pt()>25. && SignalLeptons[1]->pt()>25. && SignalLeptons[2]->pt()>20.;
262 if(!Manager()->ApplyCut(LowLeptonsPT,"low-LeptonPT > 60, 40, 30")) return true;
263 if(!Manager()->ApplyCut(ISRLeptonsPT,"ISR-LeptonPT > 25, 25, 20")) return true;
264
265 //Mll
266 MAdouble64 mll = (SignalLeptons[PairLeptons[0]]->momentum()+SignalLeptons[PairLeptons[1]]->momentum()).M();
267 if(!Manager()->ApplyCut(75. <= mll && mll <= 105.,"Dilepton InvMass")) return true;
268
269
270 /////// LOW-region ///////
271
272 //Number of jets
273 int Jet_size = Jets.size();
274 if(!Manager()->ApplyCut(Jet_size==0,"low-Jet-veto")) return true;
275
276 //low HBoost
277 MAdouble64 HBoost = atlas_susy_2018_06_func::GetHBoost(SignalLeptons,MET);
278
279 if(!Manager()->ApplyCut(HBoost>250.,"low-HBoost")) return true;
280
281 //low - PT soft
282 MAdouble64 low_PTsoft = (ThreeLeptonsSystem+MET->momentum()).Pt();
283 MAdouble64 meff = SignalLeptons[0]->pt()+SignalLeptons[1]->pt()+SignalLeptons[2]->pt()+MET->pt();
284 MAdouble64 R_PTsoft = low_PTsoft/(low_PTsoft+meff);
285 if(!Manager()->ApplyCut(R_PTsoft < 0.05,"low-PTsoft/(PTsoft+Meff)")) return true;
286 if(!Manager()->ApplyCut(meff/HBoost > 0.9,"low-Meff/HBoost")) return true;
287
288
289 /////// ISR-region ///////
290
291 //Number of jets
292 if(!Manager()->ApplyCut(Jet_size>0,"ISR-Njet > 0")) return true;
293 if(!Manager()->ApplyCut(Jet_size<4,"ISR-Njet < 4")) return true;
294
295 //ISR Delta Phi between MET and Jets
296 MALorentzVector Jets_momentum;
297 for(int i=0;i<Jet_size;i++)
298 {Jets_momentum += Jets[i]->momentum();}
299 double DeltaPhi_MET_Jets = abs((MET->momentum()).DeltaPhi(Jets_momentum));
300 if(!Manager()->ApplyCut(DeltaPhi_MET_Jets >= 2.0,"ISR-DeltaPhi(MET,Jets) > 2.0")) return true;
301
302 //ISR Ratio PT betwwen MET and Jets
303 double R_MET_Jets;
304 double R_MET_Jets_x = MET->momentum().Px()*Jets_momentum.Px();
305 double R_MET_Jets_y = MET->momentum().Py()*Jets_momentum.Py();
306 R_MET_Jets = abs(R_MET_Jets_x+R_MET_Jets_y)/pow(Jets_momentum.Pt(),2);
307
308 if(!Manager()->ApplyCut(0.55<=R_MET_Jets && R_MET_Jets <= 1.,"ISR-R(MET,Jets)")) return true;
309
310 //PT jet
311 bool ISR_JetPT = (Jets_momentum.Pt() > 100.);
312 if(!Manager()->ApplyCut(ISR_JetPT,"ISR-JetsPT > 100")) return true;
313
314 //MET PT
315 if(!Manager()->ApplyCut(MET->pt()>=80.,"ISR-MET > 80")) return true;
316
317
318 // MT
319 const RecLeptonFormat* thirdlepton = SignalLeptons[PairLeptons[2]];
320 double MT = sqrt(2*thirdlepton->pt()*MET->pt()*(1-cos(thirdlepton->dphi_0_2pi(MET->momentum()))));
321
322 if(!Manager()->ApplyCut(MT>100.,"Transverse Mass")) return true;
323
324 //ISR PT soft
325 MAdouble64 ISR_PTsoft = (ThreeLeptonsSystem+Jets_momentum+MET->momentum()).Pt();
326
327 if(!Manager()->ApplyCut(ISR_PTsoft<25.,"ISR-PTsoft < 25")) return true;
328
329 //For diagrams, If you want to draw, uncomment the paragraph below
330
331 //Fill Histogram
332/*
333 Manager()->FillHisto("SR-low-MT",MT);
334 Manager()->FillHisto("SR-low-HBoost",HBoost);
335 Manager()->FillHisto("SR-low-R(Meff,HBoost)",(meff/HBoost));
336 Manager()->FillHisto("SR-low-R(PTsoft,(PTsoft+Meff))",R_PTsoft);
337
338 Manager()->FillHisto("SR-ISR-MT",MT);
339 Manager()->FillHisto("SR-ISR-R(MET,Jets)",R_MET_Jets);
340 Manager()->FillHisto("SR-ISR-PTsoft",ISR_PTsoft);
341 Manager()->FillHisto("SR-ISR-PTjets",Jets_momentum.Pt());
342
343
344 Manager()->FillHisto("SR-low-PTlep1",SignalLeptons[0]->pt());
345 Manager()->FillHisto("SR-low-PTlep2",SignalLeptons[1]->pt());
346 Manager()->FillHisto("SR-low-PTlep3",SignalLeptons[2]->pt());
347 Manager()->FillHisto("SR-low-MET",MET->pt());
348 Manager()->FillHisto("SR-ISR-PTlep1",SignalLeptons[0]->pt());
349 Manager()->FillHisto("SR-ISR-PTlep2",SignalLeptons[1]->pt());
350 Manager()->FillHisto("SR-ISR-PTlep3",SignalLeptons[2]->pt());
351*/
352
353 return true; // end of anlaysis
354}
355
356
357////////////////////////////////////////////////////////
358//
359// Function Define
360//
361////////////////////////////////////////////////////////
362template<typename T1,typename T2> std::vector<const T1*>
363atlas_susy_2018_06_func::Removal(std::vector<const T1*> &v1, std::vector<const T2*> &v2, const MAdouble64 &drmin)
364{
365 std::vector<bool> mask(v1.size(),false);
366 for(MAuint32 j=0;j<v1.size();j++)
367 {for(MAuint32 i=0;i<v2.size();i++)
368 if(v2[i]->dr(v1[j]) < drmin)
369 {
370 mask[j]=true;
371 break;
372 }
373 }
374 std::vector<const T1*> cleaned_v1;
375 for(MAuint32 i=0;i<v1.size();i++)
376 {if(!mask[i]) cleaned_v1.push_back(v1[i]);}
377
378 return cleaned_v1;
379}
380
381
382template<typename T1,typename T2> std::vector<const T1*>
383atlas_susy_2018_06_func::RemovalLeptons(std::vector<const T1*> &v1, std::vector<const T2*> &v2)
384{
385 std::vector<bool> mask(v1.size(),false);
386 for(MAuint32 j=0;j<v1.size();j++)
387 {
388 double pt = v1[j]->pt();
389 double drmin = 0.2;
390 if(pt <= 25) drmin = 0.4;
391 else if(pt <= 50) drmin = -0.008*pt+0.6 ;
392 else continue;
393
394 for(MAuint32 i=0;i<v2.size();i++)
395 if(v2[i]->dr(v1[j]) < drmin)
396 {
397 mask[j]=true;
398 break;
399 }
400 }
401
402 std::vector<const T1*> cleaned_v1;
403 for(MAuint32 i=0;i<v1.size();i++)
404 if(!mask[i]) cleaned_v1.push_back(v1[i]);
405
406 return cleaned_v1;
407}
408
409vector<int> atlas_susy_2018_06_func::LeptonPick(vector<const RecLeptonFormat*> leptons)
410{
411 int size = leptons.size();
412 vector<int> pair(2);
413 double Zmass=91.1876;
414 double Invmass;
415 double deltaMass = 1000000;
416 for(int i=0;i<size-1;i++)
417 for(int j=i+1;j<size;j++)
418 {
419 if(leptons[i]->isMuon()==leptons[j]->isMuon())
420 if(leptons[i]->charge()!=leptons[j]->charge())
421 {
422 double temp_invmass = (leptons[i]->momentum()+leptons[j]->momentum()).M();
423 if(abs(Zmass-temp_invmass)<deltaMass)
424 {
425 Invmass = temp_invmass;
426 deltaMass = abs(Zmass-Invmass);
427 pair[0]=i;
428 pair[1]=j;
429 }
430
431 }
432 }
433
434 for(int i=0;i<size;i++)
435 if(i!=pair[0] && i!=pair[1]) pair.push_back(i);
436
437 return pair;
438}
439
440MAdouble64 atlas_susy_2018_06_func::GetHBoost(vector<const RecLeptonFormat*> leptons,const RecParticleFormat* missing)
441{
442 MAdouble64 HBoost=0;
443
444 MALorentzVector LeptonsMomentum = leptons[0]->momentum()+leptons[1]->momentum()+leptons[2]->momentum();
445 MAdouble64 Missing_Pz = (LeptonsMomentum.Pz()*missing->pt())/sqrt(pow(LeptonsMomentum.Pt(),2)+pow(LeptonsMomentum.M(),2));
446
447 MAdouble64 Missing_E = sqrt(pow(missing->pt(),2)+pow(Missing_Pz,2));
448 MALorentzVector MissingMomentum(missing->px(),missing->py(),Missing_Pz,Missing_E);
449 //cout<<MissingMomentum.M()<<endl;
450 MAdouble64 ETotal = LeptonsMomentum.E()+MissingMomentum.P();
451
452 MAdouble64 Boost_x = (LeptonsMomentum.Px()+MissingMomentum.Px())/ETotal;
453 MAdouble64 Boost_y = (LeptonsMomentum.Py()+MissingMomentum.Py())/ETotal;
454 MAdouble64 Boost_z = (LeptonsMomentum.Pz()+MissingMomentum.Pz())/ETotal;
455
456 MABoost Boost(-Boost_x,-Boost_y,-Boost_z);
457
458 for(unsigned int i=0;i<leptons.size();i++)
459 {
460 MALorentzVector temp_mom = leptons[i]->momentum();
461 Boost.boost(temp_mom);
462 //HBoost +=(Boost*leptons[i]->momentum()).P();
463 HBoost += temp_mom.P();
464 }
465 HBoost += (Boost*MissingMomentum).P();
466
467 return HBoost;
468}