ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/grimes/L1Menu/src/MenuSample.cpp
Revision: 1.9
Committed: Thu Jul 4 13:02:22 2013 UTC (11 years, 10 months ago) by grimes
Branch: MAIN
CVS Tags: HEAD
Changes since 1.8: +0 -0 lines
State: FILE REMOVED
Log Message:
Renamed ReducedMenuSample to ReducedSample, and MenuSample to FullSample

File Contents

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