ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/UHHAnalysis/SFrameTools/src/EventCalc.cxx
Revision: 1.8
Committed: Tue Aug 21 15:49:48 2012 UTC (12 years, 8 months ago) by peiffer
Content type: text/plain
Branch: MAIN
Changes since 1.7: +1 -0 lines
Log Message:
small fixes

File Contents

# User Rev Content
1 rkogler 1.1 // Dear emacs, this is -*- c++ -*-
2    
3     #include "include/EventCalc.h"
4    
5     #include <iostream>
6    
7     EventCalc* EventCalc::m_instance = NULL;
8    
9     EventCalc* EventCalc::Instance()
10     {
11     // Get a pointer to the object handler.
12     // This is the only way to access this class,
13     // since it's a singleton. This method is accessible
14     // from everywhere.
15    
16     if (m_instance == NULL){
17     m_instance = new EventCalc();
18     }
19     return m_instance;
20     }
21    
22     EventCalc::EventCalc() : m_logger( "EventCalc" )
23     {
24     // constructor: initialise all variables
25     m_logger << DEBUG << "Constructor called." << SLogger::endmsg;
26     m_bcc = NULL;
27     m_lumi = NULL;
28 peiffer 1.5 m_primlep = NULL;
29 peiffer 1.6 m_ttgen = NULL;
30 rkogler 1.1 }
31    
32     EventCalc::~EventCalc()
33     {
34     // default destructor
35     }
36    
37     void EventCalc::Reset()
38     {
39     // reset: set all booleans to false
40     // this has to be done at the beginning of each event
41     // after a call to reset all quantities will be re-calculated
42     // when they are accessed
43    
44     // (re-)set the pointers using the ObjectHandler
45     ObjectHandler* objs = ObjectHandler::Instance();
46     m_bcc = objs->GetBaseCycleContainer();
47     m_lumi = objs->GetLumiHandler();
48    
49     // reset booleans
50     b_HT = false;
51 rkogler 1.2 b_HTlep = false;
52 peiffer 1.6 b_Reconstruction = false;
53 rkogler 1.2
54 mmeyer 1.7 m_TotalWeight = 1.;
55    
56 peiffer 1.5 m_primlep = NULL;
57 peiffer 1.6 m_ttgen = NULL;
58 mmeyer 1.7
59 rkogler 1.1 }
60    
61     BaseCycleContainer* EventCalc::GetBaseCycleContainer()
62     {
63     // return the pointer to the container with all objects
64     if (!m_bcc){
65     m_logger << WARNING << "Pointer to BaseCycleContainer is NULL." << SLogger::endmsg;
66     }
67     return m_bcc;
68     }
69    
70     LuminosityHandler* EventCalc::GetLumiHandler()
71     {
72     // return the pointer to the container with all objects
73     if (!m_lumi){
74     m_logger << WARNING << "Pointer to LumiHandler is NULL." << SLogger::endmsg;
75     }
76     return m_lumi;
77     }
78    
79     double EventCalc::GetHT()
80     {
81     // calculate HT, which is defined as the scalar sum of all
82     // jets, leptons and missing transverse momentum in the event
83     if (!b_HT){
84    
85     b_HT = true;
86     m_HT = 0;
87    
88     // add lepton pt and MET
89     m_HT += GetHTlep();
90    
91     // sum over pt of all jets
92     if(m_bcc->jets){
93     for(unsigned int i=0; i<m_bcc->jets->size(); ++i){
94     m_HT += m_bcc->jets->at(i).pt();
95     }
96     }
97    
98     }
99     return m_HT;
100     }
101    
102     double EventCalc::GetHTlep()
103     {
104     // calculate HT_lep, which is defined as the scalar sum of all
105     // leptons and missing transverse momentum in the event
106     if (!b_HTlep){
107    
108     b_HTlep = true;
109     m_HTlep=0;
110    
111     // sum over pt of all electrons
112     if(m_bcc->electrons){
113     for(unsigned int i=0; i<m_bcc->electrons->size(); ++i){
114     m_HTlep += m_bcc->electrons->at(i).pt();
115     }
116     }
117    
118     // sum over pt of all muons
119     if(m_bcc->muons){
120     for(unsigned int i=0; i<m_bcc->muons->size(); ++i){
121     m_HTlep += m_bcc->muons->at(i).pt();
122     }
123     }
124    
125     // sum over pt of all taus
126     if(m_bcc->taus){
127     for(unsigned int i=0; i<m_bcc->taus->size(); ++i){
128     m_HTlep += m_bcc->taus->at(i).pt();
129     }
130     }
131    
132     // add MET
133     if(m_bcc->met) m_HTlep += m_bcc->met->pt();
134    
135     }
136     return m_HTlep;
137     }
138    
139 peiffer 1.5 Particle* EventCalc::GetPrimaryLepton(){
140    
141     if(!m_primlep){
142     double ptmax = -999;
143     for(unsigned int i=0; i<m_bcc->electrons->size(); ++i){
144     if(m_bcc->electrons->at(i).pt()>ptmax){
145     ptmax = m_bcc->electrons->at(i).pt();
146     m_primlep = &m_bcc->electrons->at(i);
147     }
148     }
149     for(unsigned int i=0; i<m_bcc->muons->size(); ++i){
150     if(m_bcc->muons->at(i).pt()>ptmax){
151     ptmax = m_bcc->muons->at(i).pt();
152     m_primlep = &m_bcc->muons->at(i);
153     }
154     }
155     }
156    
157     return m_primlep;
158     }
159    
160 peiffer 1.6 TTbarGen* EventCalc::GetTTbarGen(){
161    
162     if(!m_ttgen){
163     m_ttgen = new TTbarGen();
164     }
165    
166     return m_ttgen;
167     }
168    
169 peiffer 1.5
170     std::vector<LorentzVector> EventCalc::NeutrinoReconstruction(const LorentzVector lepton, const LorentzVector met){
171    
172    
173     TVector3 lepton_pT = toVector(lepton);
174     lepton_pT.SetZ(0);
175    
176     TVector3 neutrino_pT = toVector(met);
177     neutrino_pT.SetZ(0);
178    
179     const float mass_w = 80.399;
180     float mu = mass_w * mass_w / 2 + lepton_pT * neutrino_pT;
181    
182     float A = - (lepton_pT * lepton_pT);
183     float B = mu * lepton.pz();
184     float C = mu * mu - lepton.e() * lepton.e() * (neutrino_pT * neutrino_pT);
185    
186     float discriminant = B * B - A * C;
187    
188     std::vector<LorentzVector> solutions;
189    
190     if (0 >= discriminant)
191     {
192     // Take only real part of the solution
193     //
194     LorentzVectorXYZE solution (0,0,0,0);
195     solution.SetPx(met.Px());
196     solution.SetPy(met.Py());
197     solution.SetPz(-B / A);
198     solution.SetE(toVector(solution).Mag());
199    
200     solutions.push_back(toPtEtaPhi(solution));
201    
202     //_solutions = 0 > discriminant ? 0 : 1;
203     }
204     else
205     {
206     discriminant = sqrt(discriminant);
207    
208     LorentzVectorXYZE solution (0,0,0,0);
209     solution.SetPx(met.Px());
210     solution.SetPy(met.Py());
211     solution.SetPz((-B - discriminant) / A);
212     solution.SetE(toVector(solution).Mag());
213    
214     solutions.push_back(toPtEtaPhi(solution));
215    
216     LorentzVectorXYZE solution2 (0,0,0,0);
217     solution2.SetPx(met.Px());
218     solution2.SetPy(met.Py());
219     solution2.SetPz((-B + discriminant) / A);
220     solution2.SetE(toVector(solution2).Mag());
221    
222     solutions.push_back(toPtEtaPhi(solution2));
223    
224     //_solutions = 2;
225     }
226    
227     return solutions;
228     }
229    
230     void EventCalc::FillHighMassTTbarHypotheses(){
231    
232 peiffer 1.6 if(b_Reconstruction) return;
233     b_Reconstruction=true;
234    
235 peiffer 1.5 //clear hypothesis list
236     m_bcc->recoHyps->clear();
237    
238     //find primary charged lepton
239     Particle* lepton = GetPrimaryLepton();
240    
241     //reconstruct neutrino
242     std::vector<LorentzVector> neutrinos = NeutrinoReconstruction( lepton->v4(), m_bcc->met->v4());
243    
244     ReconstructionHypothesis hyp;
245    
246     hyp.set_lepton(*lepton);
247    
248     //loop over neutrino solutions and jet assignments to fill hyotheses
249     for(unsigned int i=0; i< neutrinos.size();++i){
250    
251     hyp.set_neutrino_v4(neutrinos[i]);
252     LorentzVector wlep_v4 = lepton->v4()+neutrinos[i];
253    
254     unsigned int n_jets = m_bcc->jets->size();
255     unsigned int max_j = myPow(3, n_jets);
256 peiffer 1.8 if(max_j>10) max_j=10; //avoid crashes in events with many jets
257 peiffer 1.5 for (unsigned int j=0; j < max_j; j++) {
258     LorentzVector tophad_v4(0,0,0,0);
259     LorentzVector toplep_v4 = wlep_v4;
260     int hadjets=0;
261     int lepjets=0;
262     int num = j;
263     hyp.clear_jetindices();
264     for (unsigned int k=0; k<n_jets; k++) {
265     // num is the k-th digit of j if you
266     // write j in a base-3 system. According
267     // to the value of this digit (which takes
268     // values from 0 to 2,
269     // in all possible combinations with the other digits),
270     // decide how to treat the jet.
271    
272     if(num%3==0) {
273     tophad_v4 = tophad_v4 + m_bcc->jets->at(k).v4();
274     hyp.add_tophad_jet_index(k);
275     hadjets++;
276     }
277    
278     if(num%3==1) {
279     toplep_v4 = toplep_v4 + m_bcc->jets->at(k).v4();
280     hyp.add_toplep_jet_index(k);
281     lepjets++;
282     }
283     //if(num%3==2); //do not take this jet
284    
285     //shift the trigits of num to the right:
286     num /= 3;
287     }
288    
289 peiffer 1.6 //search jet with highest pt assigned to leptonic top
290     float maxpt=-999;
291     int maxind=-1;
292     for(unsigned int i=0; i<hyp.toplep_jets_indices().size(); ++i){
293     float pt = m_bcc->jets->at(hyp.toplep_jets_indices().at(i)).pt();
294     if(pt>maxpt){
295     maxpt=pt;
296     maxind=hyp.toplep_jets_indices().at(i);
297     }
298     }
299     hyp.set_blep_index(maxind);
300    
301 peiffer 1.5
302 peiffer 1.6 //fill only hypotheses with at least one jet assigned to each top quark
303 peiffer 1.5 if(hadjets>0 && lepjets>0){
304     hyp.set_tophad_v4(tophad_v4);
305     hyp.set_toplep_v4(toplep_v4);
306     m_bcc->recoHyps->push_back(hyp);
307     }
308     }
309     }
310    
311     }
312    
313    
314 peiffer 1.3 void EventCalc::PrintEventContent(){
315 rkogler 1.1
316 peiffer 1.3 m_logger << INFO << "----------------- event content -----------------" << SLogger::endmsg;
317     m_logger << INFO << "run: " << m_bcc->run << " lumi block:" << m_bcc->luminosityBlock << " event: " << m_bcc->event << SLogger::endmsg;
318     m_logger << INFO << "MET = " << m_bcc->met->pt() << " METphi = " << m_bcc->met->phi() << " HTlep = " << GetHTlep() << SLogger::endmsg;
319     if(m_bcc->electrons){m_logger << INFO << "Electrons:" << SLogger::endmsg;
320     for(unsigned int i=0; i<m_bcc->electrons->size(); ++i){
321 peiffer 1.5 m_logger << INFO << " " << i+1 << " pt = " << m_bcc->electrons->at(i).pt() <<" eta = " << m_bcc->electrons->at(i).eta()
322     << " supercluster eta = " <<m_bcc->electrons->at(i).supercluster_eta() << " IP wrt bsp = " << fabs(m_bcc->electrons->at(i).gsfTrack_dxy_vertex(m_bcc->pvs->at(0).x(), m_bcc->pvs->at(0).y()))
323     << " pass conv veto = " << m_bcc->electrons->at(i).passconversionveto() << " mvaTrigV0 = " << m_bcc->electrons->at(i).mvaTrigV0()
324     << " dEtaIn = " << m_bcc->electrons->at(i).dEtaIn() << " sigmaIEtaIEta = " << m_bcc->electrons->at(i).sigmaIEtaIEta()
325     << " HoverE = " << m_bcc->electrons->at(i).HoverE() << " EcalEnergy = " << m_bcc->electrons->at(i).EcalEnergy()
326     << " EoverPIn = " << m_bcc->electrons->at(i).EoverPIn() << " trackMomentumAtVtx = " << m_bcc->electrons->at(i).EcalEnergy()/m_bcc->electrons->at(i).EoverPIn()
327     << SLogger::endmsg;
328 peiffer 1.3 }
329     }
330     if(m_bcc->muons){m_logger << INFO << "Muons:" << SLogger::endmsg;
331     for(unsigned int i=0; i<m_bcc->muons->size(); ++i){
332     m_logger << INFO << " " << i+1 << " pt = " << m_bcc->muons->at(i).pt() <<" eta = " << m_bcc->muons->at(i).eta() << SLogger::endmsg;
333     }
334     }
335     if(m_bcc->taus){m_logger << INFO << "Taus:" << SLogger::endmsg;
336     for(unsigned int i=0; i<m_bcc->taus->size(); ++i){
337     m_logger << INFO << " " << i+1 << " pt = " << m_bcc->taus->at(i).pt() <<" eta = " << m_bcc->taus->at(i).eta() << SLogger::endmsg;
338     }
339     }
340     if(m_bcc->jets){m_logger << INFO << "Jets:" << SLogger::endmsg;
341     for(unsigned int i=0; i<m_bcc->jets->size(); ++i){
342 peiffer 1.4 m_logger << INFO << " " << i+1 << " pt = " << m_bcc->jets->at(i).pt() <<" eta = " << m_bcc->jets->at(i).eta() << " genjet_pt = " << m_bcc->jets->at(i).genjet_pt() << " genjet_eta = " << m_bcc->jets->at(i).genjet_eta() <<SLogger::endmsg;
343 peiffer 1.3 }
344     }
345     if(m_bcc->topjets){m_logger << INFO << "TopJets:" << SLogger::endmsg;
346     for(unsigned int i=0; i<m_bcc->topjets->size(); ++i){
347     m_logger << INFO << " " << i+1 << " pt = " << m_bcc->topjets->at(i).pt() <<" eta = " << m_bcc->topjets->at(i).eta() << SLogger::endmsg;
348     }
349     }
350     if(m_bcc->photons){m_logger << INFO << "Photons:" << SLogger::endmsg;
351     for(unsigned int i=0; i<m_bcc->photons->size(); ++i){
352     m_logger << INFO << " " << i+1 << " pt = " << m_bcc->photons->at(i).pt() <<" eta = " << m_bcc->photons->at(i).eta() << SLogger::endmsg;
353     }
354     }
355     }
356 mmeyer 1.7
357    
358     void EventCalc::ProduceWeight(double weight)
359     {
360     m_TotalWeight = m_TotalWeight * weight;
361     }
362    
363     double EventCalc::GetWeight()
364     {
365    
366     return m_TotalWeight;
367     }
368