ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/grimes/L1Menu/src/MenuSample.cpp
Revision: 1.4
Committed: Tue May 28 23:14:03 2013 UTC (11 years, 11 months ago) by grimes
Branch: MAIN
Changes since 1.3: +127 -55 lines
Log Message:
Numerous changes

File Contents

# Content
1 #include "l1menu/MenuSample.h"
2
3 #include <iostream>
4 #include <stdexcept>
5 #include <cmath>
6
7 #include <TSystem.h>
8 #include "UserCode/L1TriggerUpgrade/macros/L1UpgradeNtuple.h"
9 #include "UserCode/L1TriggerUpgrade/interface/L1AnalysisDataFormat.h"
10 #include "l1menu/IEvent.h"
11
12 // I'll use the unnamed namespace for things that are only used in this file
13 namespace
14 {
15 /** Implementation of the l1menu::IEvent interface.
16 *
17 * @author Mark Grimes (mark.grimes@bristol.ac.uk)
18 * @date 21/May/2013
19 */
20 class EventImplementation : public l1menu::IEvent
21 {
22 public:
23 //
24 // These are the methods required by the IEvent interface
25 //
26 virtual L1Analysis::L1AnalysisDataFormat& rawEvent() { return rawEvent_; }
27 virtual const L1Analysis::L1AnalysisDataFormat& rawEvent() const { return rawEvent_; }
28 virtual bool* physicsBits() { return physicsBits_; }
29 virtual const bool* physicsBits() const { return physicsBits_; }
30 virtual float weight() const { return 1; }
31 protected:
32 L1Analysis::L1AnalysisDataFormat rawEvent_;
33 bool physicsBits_[128];
34 }; // end of the EventImplementation class
35
36 } // end of the unnamed namespace
37
38 namespace l1menu
39 {
40 class MenuSamplePrivateMembers
41 {
42 private:
43 static const size_t PHIBINS;
44 static const double PHIBIN[];
45 static const size_t ETABINS;
46 static const double ETABIN[];
47
48 double degree( double radian );
49 int phiINjetCoord( double phi );
50 int etaINjetCoord( double eta );
51 double calculateHTT( const L1Analysis::L1AnalysisDataFormat& event );
52 double calculateHTM( const L1Analysis::L1AnalysisDataFormat& event );
53 public:
54 void fillDataStructure( int selectDataInput );
55 void fillL1Bits();
56 L1UpgradeNtuple inputNtuple;
57 ::EventImplementation currentEvent;
58 };
59 }
60
61 const size_t l1menu::MenuSamplePrivateMembers::PHIBINS=18;
62 const double l1menu::MenuSamplePrivateMembers::PHIBIN[]={10,30,50,70,90,110,130,150,170,190,210,230,250,270,290,310,330,350};
63 const size_t l1menu::MenuSamplePrivateMembers::ETABINS=23;
64 const double l1menu::MenuSamplePrivateMembers::ETABIN[]={-5.,-4.5,-4.,-3.5,-3.,-2.172,-1.74,-1.392,-1.044,-0.696,-0.348,0,0.348,0.696,1.044,1.392,1.74,2.172,3.,3.5,4.,4.5,5.};
65
66
67 double l1menu::MenuSamplePrivateMembers::degree( double radian )
68 {
69 if( radian<0 ) return 360.+(radian/M_PI*180.);
70 else return radian/M_PI*180.;
71 }
72
73 int l1menu::MenuSamplePrivateMembers::phiINjetCoord( double phi )
74 {
75 size_t phiIdx=0;
76 double phidegree=degree( phi );
77 for( size_t idx=0; idx<PHIBINS; idx++ )
78 {
79 if( phidegree>=PHIBIN[idx] and phidegree<PHIBIN[idx+1] ) phiIdx=idx;
80 else if( phidegree>=PHIBIN[PHIBINS-1] || phidegree<=PHIBIN[0] ) phiIdx=idx;
81 }
82 phiIdx=phiIdx+1;
83 if( phiIdx==18 ) phiIdx=0;
84 return int( phiIdx );
85 }
86
87 int l1menu::MenuSamplePrivateMembers::etaINjetCoord( double eta )
88 {
89 size_t etaIdx=0.;
90 for( size_t idx=0; idx<ETABINS; idx++ )
91 {
92 if( eta>=ETABIN[idx] and eta<ETABIN[idx+1] ) etaIdx=idx;
93 }
94 return int( etaIdx );
95 }
96
97 double l1menu::MenuSamplePrivateMembers::calculateHTT( const L1Analysis::L1AnalysisDataFormat& event )
98 {
99 double httValue=0.;
100
101 // Calculate our own HT and HTM from the jets that survive the double jet removal.
102 for( int i=0; i<event.Njet; i++ )
103 {
104 if( event.Bxjet.at( i )==0 )
105 {
106 if( event.Etajet.at( i )>4 and event.Etajet.at( i )<17 )
107 {
108 httValue+=event.Etjet.at( i );
109 } //in proper eta range
110 } //correct beam crossing
111 } //loop over cleaned jets
112
113 return httValue;
114 }
115
116 double l1menu::MenuSamplePrivateMembers::calculateHTM( const L1Analysis::L1AnalysisDataFormat& event )
117 {
118 double htmValue=0.;
119 double htmValueX=0.;
120 double htmValueY=0.;
121
122 // Calculate our own HT and HTM from the jets that survive the double jet removal.
123 for( int i=0; i<event.Njet; i++ )
124 {
125 if( event.Bxjet.at( i )==0 )
126 {
127 if( event.Etajet.at( i )>4 and event.Etajet.at( i )<17 )
128 {
129
130 // Get the phi angle towers are 0-17 (this is probably not real mapping but OK for just magnitude of HTM
131 float phi=2*M_PI*(event.Phijet.at( i )/18.);
132 htmValueX+=cos( phi )*event.Etjet.at( i );
133 htmValueY+=sin( phi )*event.Etjet.at( i );
134
135 } //in proper eta range
136 } //correct beam crossing
137 } //loop over cleaned jets
138
139 htmValue=sqrt( htmValueX*htmValueX+htmValueY*htmValueY );
140
141 return htmValue;
142 }
143
144 void l1menu::MenuSamplePrivateMembers::fillDataStructure( int selectDataInput )
145 {
146 // Use a reference for ease of use
147 L1Analysis::L1AnalysisDataFormat& analysisDataFormat=currentEvent.rawEvent();
148
149 analysisDataFormat.Reset();
150
151 // Grab standard event information
152 analysisDataFormat.Run=inputNtuple.event_->run;
153 analysisDataFormat.LS=inputNtuple.event_->lumi;
154 analysisDataFormat.Event=inputNtuple.event_->event;
155
156 /* =======================================================================================================
157 / Select the input source information
158 / ---------------------------------------------------------------------------------------------------------
159 / case 0: Use Original L1ExtraTree that came with the event
160 / case 1: Use reEmulated L1ExtraTree that was produced (note this may not be present or may be identical to the Original tree)
161 / case 2: Use Original L1ExtraTree that came with the event except for muons which get from GMT. (For old Ntuple with no quality flag in L1Extra)
162 /
163 / case 10: Use Original L1Tree (GMT/GT) that came with the event
164 / case 11: Use reEmulated L1Tree (GMT/GT) that was produced (note this may not be present or may be identical to the Original tree)
165 / case 12: Use Original L1Tree (GMT/GCT) that was produced
166 / case 13: Use reEmulated L1Tree (GMT/GCT) that was produced (note this may not be present or may be identical to the Original tree)
167 /
168 / case 21: Use the L1ExtraUpgradeTree (Assuming Stage 1 Quantities filled in tree)
169 / case 22: Use the L1ExtraUpgradeTree (Assuming Stage 2 Quantities filled in tree)
170 / ======================================================================================================= */
171 switch( selectDataInput )
172 {
173 case 22: //Select from L1ExtraUpgradeTree (Stage 2)
174
175 // NOTES: Stage 1 has EG Relaxed and EG Isolated. The isolated EG are a subset of the Relaxed.
176 // so sort through the relaxed list and flag those that also appear in the isolated list.
177 for( unsigned int i=0; i<inputNtuple.l1upgrade_->nEG; i++ )
178 {
179
180 analysisDataFormat.Bxel.push_back( inputNtuple.l1upgrade_->egBx.at( i ) );
181 analysisDataFormat.Etel.push_back( inputNtuple.l1upgrade_->egEt.at( i ) );
182 analysisDataFormat.Phiel.push_back( phiINjetCoord( inputNtuple.l1upgrade_->egPhi.at( i ) ) ); //PROBLEM: real value, trigger wants bin convert with phiINjetCoord
183 analysisDataFormat.Etael.push_back( etaINjetCoord( inputNtuple.l1upgrade_->egEta.at( i ) ) ); //PROBLEM: real value, trigger wants bin convert with etaINjetCoord
184
185 // Check whether this EG is located in the isolation list
186 bool isolated=false;
187 bool fnd=false;
188 unsigned int isoEG=0;
189 while( !fnd && isoEG < inputNtuple.l1upgrade_->nIsoEG )
190 {
191 if( inputNtuple.l1upgrade_->isoEGPhi.at( isoEG )==inputNtuple.l1upgrade_->egPhi.at( i )
192 && inputNtuple.l1upgrade_->isoEGEta.at( isoEG )==inputNtuple.l1upgrade_->egEta.at( i ) )
193 {
194 isolated=true;
195 fnd=true;
196 }
197 isoEG++;
198 }
199 analysisDataFormat.Isoel.push_back( isolated );
200 analysisDataFormat.Nele++;
201 }
202
203 // Note: Taus are in the jet list. Decide what to do with them. For now
204 // leave them the there as jets (not even flagged..)
205 for( unsigned int i=0; i<inputNtuple.l1upgrade_->nJets; i++ )
206 {
207
208 // For each jet look for a possible duplicate if so remove it.
209 bool duplicate=false;
210 for( unsigned int j=0; j<i; j++ )
211 {
212 if( inputNtuple.l1upgrade_->jetBx.at( i )==inputNtuple.l1upgrade_->jetBx.at( j )
213 && inputNtuple.l1upgrade_->jetEt.at( i )==inputNtuple.l1upgrade_->jetEt.at( j )
214 && inputNtuple.l1upgrade_->jetEta.at( i )==inputNtuple.l1upgrade_->jetEta.at( j )
215 && inputNtuple.l1upgrade_->jetPhi.at( i )==inputNtuple.l1upgrade_->jetPhi.at( j ) )
216 {
217 duplicate=true;
218 //printf("Duplicate jet found and removed \n");
219 }
220 }
221
222 if( !duplicate )
223 {
224 analysisDataFormat.Bxjet.push_back( inputNtuple.l1upgrade_->jetBx.at( i ) );
225 analysisDataFormat.Etjet.push_back( inputNtuple.l1upgrade_->jetEt.at( i ) );
226 analysisDataFormat.Phijet.push_back( phiINjetCoord( inputNtuple.l1upgrade_->jetPhi.at( i ) ) ); //PROBLEM: real value, trigger wants bin convert with phiINjetCoord
227 analysisDataFormat.Etajet.push_back( etaINjetCoord( inputNtuple.l1upgrade_->jetEta.at( i ) ) ); //PROBLEM: real value, trigger wants bin convert with etaINjetCoord
228 analysisDataFormat.Taujet.push_back( false );
229 analysisDataFormat.isoTaujet.push_back( false );
230 //analysisDataFormat.Fwdjet.push_back(false); //COMMENT OUT IF JET ETA FIX
231
232 //if(fabs(inputNtuple.l1upgrade_->jetEta.at(i))>=3.0) printf("Et %f Eta %f iEta %f Phi %f iPhi %f \n",analysisDataFormat.Etjet.at(analysisDataFormat.Njet),inputNtuple.l1upgrade_->jetEta.at(i),analysisDataFormat.Etajet.at(analysisDataFormat.Njet),inputNtuple.l1upgrade_->jetPhi.at(i),analysisDataFormat.Phijet.at(analysisDataFormat.Njet));
233 // Eta Jet Fix. Some Jets with eta>3 has appeared in central jet list. Move them by hand
234 // This is a problem in Stage 2 Jet code.
235 (fabs( inputNtuple.l1upgrade_->jetEta.at( i ) )>=3.0) ? analysisDataFormat.Fwdjet.push_back( true ) : analysisDataFormat.Fwdjet.push_back( false );
236
237 analysisDataFormat.Njet++;
238 }
239 }
240
241 for( unsigned int i=0; i<inputNtuple.l1upgrade_->nFwdJets; i++ )
242 {
243
244 analysisDataFormat.Bxjet.push_back( inputNtuple.l1upgrade_->fwdJetBx.at( i ) );
245 analysisDataFormat.Etjet.push_back( inputNtuple.l1upgrade_->fwdJetEt.at( i ) );
246 analysisDataFormat.Phijet.push_back( phiINjetCoord( inputNtuple.l1upgrade_->fwdJetPhi.at( i ) ) ); //PROBLEM: real value, trigger wants bin convert with phiINjetCoord
247 analysisDataFormat.Etajet.push_back( etaINjetCoord( inputNtuple.l1upgrade_->fwdJetEta.at( i ) ) ); //PROBLEM: real value, trigger wants bin convert with etaINjetCoord
248 analysisDataFormat.Taujet.push_back( false );
249 analysisDataFormat.isoTaujet.push_back( false );
250 analysisDataFormat.Fwdjet.push_back( true );
251
252 analysisDataFormat.Njet++;
253 }
254
255 // NOTES: Stage 1 has Tau Relaxed and TauIsolated. The isolated Tau are a subset of the Relaxed.
256 // so sort through the relaxed list and flag those that also appear in the isolated list.
257
258 for( unsigned int i=0; i<inputNtuple.l1upgrade_->nTau; i++ )
259 {
260
261 // remove duplicates
262 bool duplicate=false;
263 for( unsigned int j=0; j<i; j++ )
264 {
265 if( inputNtuple.l1upgrade_->tauBx.at( i )==inputNtuple.l1upgrade_->tauBx.at( j )
266 && inputNtuple.l1upgrade_->tauEt.at( i )==inputNtuple.l1upgrade_->tauEt.at( j )
267 && inputNtuple.l1upgrade_->tauEta.at( i )==inputNtuple.l1upgrade_->tauEta.at( j )
268 && inputNtuple.l1upgrade_->tauPhi.at( i )==inputNtuple.l1upgrade_->tauPhi.at( j ) )
269 {
270 duplicate=true;
271 //printf("Duplicate jet found and removed \n");
272 }
273 }
274
275 if( !duplicate )
276 {
277 analysisDataFormat.Bxjet.push_back( inputNtuple.l1upgrade_->tauBx.at( i ) );
278 analysisDataFormat.Etjet.push_back( inputNtuple.l1upgrade_->tauEt.at( i ) );
279 analysisDataFormat.Phijet.push_back( phiINjetCoord( inputNtuple.l1upgrade_->tauPhi.at( i ) ) ); //PROBLEM: real value, trigger wants bin convert with phiINjetCoord
280 analysisDataFormat.Etajet.push_back( etaINjetCoord( inputNtuple.l1upgrade_->tauEta.at( i ) ) ); //PROBLEM: real value, trigger wants bin convert with etaINjetCoord
281 analysisDataFormat.Taujet.push_back( true );
282 analysisDataFormat.Fwdjet.push_back( false );
283
284 bool isolated=false;
285 bool fnd=false;
286 unsigned int isoTau=0;
287 while( !fnd && isoTau < inputNtuple.l1upgrade_->nIsoTau )
288 {
289 if( inputNtuple.l1upgrade_->isoTauPhi.at( isoTau )==inputNtuple.l1upgrade_->tauPhi.at( i )
290 && inputNtuple.l1upgrade_->isoTauEta.at( isoTau )==inputNtuple.l1upgrade_->tauEta.at( i ) )
291 {
292 isolated=true;
293 fnd=true;
294 }
295 isoTau++;
296 }
297 analysisDataFormat.isoTaujet.push_back( isolated );
298
299 analysisDataFormat.Njet++;
300 } // duplicate check
301 }
302
303 // Fill energy sums (Are overflow flags accessible in l1extra?)
304 for( unsigned int i=0; i<inputNtuple.l1upgrade_->nMet; i++ )
305 {
306 //if(inputNtuple.l1upgrade_->metBx.at(i)==0) {
307 analysisDataFormat.ETT=inputNtuple.l1upgrade_->et.at( i );
308 analysisDataFormat.ETM=inputNtuple.l1upgrade_->met.at( i );
309 analysisDataFormat.PhiETM=inputNtuple.l1upgrade_->metPhi.at( i );
310 }
311 analysisDataFormat.OvETT=0; //not available in l1extra
312 analysisDataFormat.OvETM=0; //not available in l1extra
313
314 for( unsigned int i=0; i<inputNtuple.l1upgrade_->nMht; i++ )
315 {
316 if( inputNtuple.l1upgrade_->mhtBx.at( i )==0 )
317 {
318 analysisDataFormat.HTT=calculateHTT( analysisDataFormat ); //inputNtuple.l1upgrade_->ht.at(i) ;
319 analysisDataFormat.HTM=calculateHTM( analysisDataFormat ); //inputNtuple.l1upgrade_->mht.at(i) ;
320 analysisDataFormat.PhiHTM=0.; //inputNtuple.l1upgrade_->mhtPhi.at(i) ;
321 }
322 }
323 analysisDataFormat.OvHTM=0; //not available in l1extra
324 analysisDataFormat.OvHTT=0; //not available in l1extra
325
326 // Get the muon information from reEmul GMT
327 for( int i=0; i<inputNtuple.gmtEmu_->N; i++ )
328 {
329
330 analysisDataFormat.Bxmu.push_back( inputNtuple.gmtEmu_->CandBx[i] );
331 analysisDataFormat.Ptmu.push_back( inputNtuple.gmtEmu_->Pt[i] );
332 analysisDataFormat.Phimu.push_back( inputNtuple.gmtEmu_->Phi[i] );
333 analysisDataFormat.Etamu.push_back( inputNtuple.gmtEmu_->Eta[i] );
334 analysisDataFormat.Qualmu.push_back( inputNtuple.gmtEmu_->Qual[i] );
335 analysisDataFormat.Isomu.push_back( false );
336 }
337
338 break;
339
340 default:
341 std::cout << "---Not a valid input source FULL STOP! " << std::endl;
342
343 break;
344 }
345
346 return;
347 }
348
349 void l1menu::MenuSamplePrivateMembers::fillL1Bits()
350 {
351 bool* PhysicsBits=currentEvent.physicsBits();
352
353 // I really don't think this if statement is correct. Surely it
354 // should be "if( gt_ )"? - M. Grimes.
355 if( !inputNtuple.gt_ )
356 {
357 for( Int_t ibit=0; ibit<128; ibit++ )
358 {
359 PhysicsBits[ibit]=0;
360 if( ibit<64 )
361 {
362 PhysicsBits[ibit]=(inputNtuple.gt_->tw1[2]>>ibit)&1;
363 }
364 else
365 {
366 PhysicsBits[ibit]=(inputNtuple.gt_->tw2[2]>>(ibit-64))&1;
367 }
368 }
369 }
370 else
371 {
372 PhysicsBits[0]=1; //set zero bias on if no gt information
373 }
374 }
375
376
377 l1menu::MenuSample::MenuSample()
378 : pImple_( new MenuSamplePrivateMembers )
379 {
380 // No operation besides the initialiser list
381 }
382
383 l1menu::MenuSample::~MenuSample()
384 {
385 delete pImple_;
386 }
387
388 l1menu::MenuSample::MenuSample( const l1menu::MenuSample& otherMenuSample )
389 : pImple_( new MenuSamplePrivateMembers(*otherMenuSample.pImple_) )
390 {
391 // No operation besides the initialiser list
392 }
393
394 l1menu::MenuSample::MenuSample( l1menu::MenuSample&& otherMenuSample ) noexcept
395 : pImple_( otherMenuSample.pImple_ )
396 {
397 otherMenuSample.pImple_=NULL;
398 }
399
400 l1menu::MenuSample& l1menu::MenuSample::operator=( const l1menu::MenuSample& otherMenuSample )
401 {
402 *pImple_=*otherMenuSample.pImple_;
403 return *this;
404 }
405
406 l1menu::MenuSample& l1menu::MenuSample::operator=( l1menu::MenuSample&& otherMenuSample ) noexcept
407 {
408 pImple_=otherMenuSample.pImple_;
409 otherMenuSample.pImple_=NULL;
410 return *this;
411 }
412
413 void l1menu::MenuSample::loadFile( const std::string& filename )
414 {
415 pImple_->inputNtuple.Open( filename );
416 }
417
418 size_t l1menu::MenuSample::numberOfEvents() const
419 {
420 return static_cast<size_t>( pImple_->inputNtuple.GetEntries() );
421 }
422
423 const l1menu::IEvent& l1menu::MenuSample::getEvent( size_t eventNumber ) const
424 {
425 // Make sure the event number requested is valid. Use static_cast to get rid
426 // of the "comparison between signed and unsigned" compiler warning.
427 if( eventNumber>static_cast<size_t>(pImple_->inputNtuple.GetEntries()) ) throw std::runtime_error( "Requested event number is out of range" );
428
429 pImple_->inputNtuple.LoadTree(eventNumber);
430 pImple_->inputNtuple.GetEntry(eventNumber);
431 // This next call fills pImple_->currentEvent with the information in pImple_->inputNtuple
432 pImple_->fillDataStructure( 22 );
433 pImple_->fillL1Bits();
434
435 return pImple_->currentEvent;
436 }