ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/grimes/L1Menu/src/FullSample.cpp
Revision: 1.5
Committed: Wed Jul 24 11:48:55 2013 UTC (11 years, 9 months ago) by grimes
Branch: MAIN
CVS Tags: HEAD
Changes since 1.4: +2 -1 lines
Log Message:
Big commit of lots of changes before migrating to Git.

File Contents

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