ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/grimes/L1Menu/src/FullSample.cpp
Revision: 1.4
Committed: Fri Jul 5 14:21:56 2013 UTC (11 years, 10 months ago) by grimes
Branch: MAIN
Changes since 1.3: +1 -0 lines
Log Message:
Set the event weight in FullSamples properly

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 "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 FullSamplePrivateMembers
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 FullSamplePrivateMembers( FullSample* 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::FullSamplePrivateMembers::PHIBINS=18;
66 const double l1menu::FullSamplePrivateMembers::PHIBIN[]={10,30,50,70,90,110,130,150,170,190,210,230,250,270,290,310,330,350};
67 const size_t l1menu::FullSamplePrivateMembers::ETABINS=23;
68 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.};
69
70
71 double l1menu::FullSamplePrivateMembers::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::FullSamplePrivateMembers::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::FullSamplePrivateMembers::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::FullSamplePrivateMembers::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 && !event.Taujet.at(i) )
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::FullSamplePrivateMembers::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 && !event.Taujet.at(i) )
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::FullSamplePrivateMembers::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 currentEvent.setWeight( inputNtuple.event_->puWeight );
157 analysisDataFormat.Run=inputNtuple.event_->run;
158 analysisDataFormat.LS=inputNtuple.event_->lumi;
159 analysisDataFormat.Event=inputNtuple.event_->event;
160
161 /* =======================================================================================================
162 / Select the input source information
163 / ---------------------------------------------------------------------------------------------------------
164 / case 0: Use Original L1ExtraTree that came with the event
165 / case 1: Use reEmulated L1ExtraTree that was produced (note this may not be present or may be identical to the Original tree)
166 / 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)
167 /
168 / case 10: Use Original L1Tree (GMT/GT) that came with the event
169 / case 11: Use reEmulated L1Tree (GMT/GT) that was produced (note this may not be present or may be identical to the Original tree)
170 / case 12: Use Original L1Tree (GMT/GCT) that was produced
171 / case 13: Use reEmulated L1Tree (GMT/GCT) that was produced (note this may not be present or may be identical to the Original tree)
172 /
173 / case 21: Use the L1ExtraUpgradeTree (Assuming Stage 1 Quantities filled in tree)
174 / case 22: Use the L1ExtraUpgradeTree (Assuming Stage 2 Quantities filled in tree)
175 / ======================================================================================================= */
176 switch( selectDataInput )
177 {
178 case 22: //Select from L1ExtraUpgradeTree (Stage 2)
179
180 // NOTES: Stage 1 has EG Relaxed and EG Isolated. The isolated EG are a subset of the Relaxed.
181 // so sort through the relaxed list and flag those that also appear in the isolated list.
182 for( unsigned int i=0; i<inputNtuple.l1upgrade_->nEG; i++ )
183 {
184
185 analysisDataFormat.Bxel.push_back( inputNtuple.l1upgrade_->egBx.at( i ) );
186 analysisDataFormat.Etel.push_back( inputNtuple.l1upgrade_->egEt.at( i ) );
187 analysisDataFormat.Phiel.push_back( phiINjetCoord( inputNtuple.l1upgrade_->egPhi.at( i ) ) ); //PROBLEM: real value, trigger wants bin convert with phiINjetCoord
188 analysisDataFormat.Etael.push_back( etaINjetCoord( inputNtuple.l1upgrade_->egEta.at( i ) ) ); //PROBLEM: real value, trigger wants bin convert with etaINjetCoord
189
190 // Check whether this EG is located in the isolation list
191 bool isolated=false;
192 bool fnd=false;
193 unsigned int isoEG=0;
194 while( !fnd && isoEG < inputNtuple.l1upgrade_->nIsoEG )
195 {
196 if( inputNtuple.l1upgrade_->isoEGPhi.at( isoEG )==inputNtuple.l1upgrade_->egPhi.at( i )
197 && inputNtuple.l1upgrade_->isoEGEta.at( isoEG )==inputNtuple.l1upgrade_->egEta.at( i ) )
198 {
199 isolated=true;
200 fnd=true;
201 }
202 isoEG++;
203 }
204 analysisDataFormat.Isoel.push_back( isolated );
205 analysisDataFormat.Nele++;
206 }
207
208 // Note: Taus are in the jet list. Decide what to do with them. For now
209 // leave them the there as jets (not even flagged..)
210 for( unsigned int i=0; i<inputNtuple.l1upgrade_->nJets; i++ )
211 {
212
213 // For each jet look for a possible duplicate if so remove it.
214 bool duplicate=false;
215 for( unsigned int j=0; j<i; j++ )
216 {
217 if( inputNtuple.l1upgrade_->jetBx.at( i )==inputNtuple.l1upgrade_->jetBx.at( j )
218 && inputNtuple.l1upgrade_->jetEt.at( i )==inputNtuple.l1upgrade_->jetEt.at( j )
219 && inputNtuple.l1upgrade_->jetEta.at( i )==inputNtuple.l1upgrade_->jetEta.at( j )
220 && inputNtuple.l1upgrade_->jetPhi.at( i )==inputNtuple.l1upgrade_->jetPhi.at( j ) )
221 {
222 duplicate=true;
223 //printf("Duplicate jet found and removed \n");
224 }
225 }
226
227 if( !duplicate )
228 {
229 analysisDataFormat.Bxjet.push_back( inputNtuple.l1upgrade_->jetBx.at( i ) );
230 analysisDataFormat.Etjet.push_back( inputNtuple.l1upgrade_->jetEt.at( i ) );
231 analysisDataFormat.Phijet.push_back( phiINjetCoord( inputNtuple.l1upgrade_->jetPhi.at( i ) ) ); //PROBLEM: real value, trigger wants bin convert with phiINjetCoord
232 analysisDataFormat.Etajet.push_back( etaINjetCoord( inputNtuple.l1upgrade_->jetEta.at( i ) ) ); //PROBLEM: real value, trigger wants bin convert with etaINjetCoord
233 analysisDataFormat.Taujet.push_back( false );
234 analysisDataFormat.isoTaujet.push_back( false );
235 //analysisDataFormat.Fwdjet.push_back(false); //COMMENT OUT IF JET ETA FIX
236
237 //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));
238 // Eta Jet Fix. Some Jets with eta>3 has appeared in central jet list. Move them by hand
239 // This is a problem in Stage 2 Jet code.
240 (fabs( inputNtuple.l1upgrade_->jetEta.at( i ) )>=3.0) ? analysisDataFormat.Fwdjet.push_back( true ) : analysisDataFormat.Fwdjet.push_back( false );
241
242 analysisDataFormat.Njet++;
243 }
244 }
245
246 for( unsigned int i=0; i<inputNtuple.l1upgrade_->nFwdJets; i++ )
247 {
248
249 analysisDataFormat.Bxjet.push_back( inputNtuple.l1upgrade_->fwdJetBx.at( i ) );
250 analysisDataFormat.Etjet.push_back( inputNtuple.l1upgrade_->fwdJetEt.at( i ) );
251 analysisDataFormat.Phijet.push_back( phiINjetCoord( inputNtuple.l1upgrade_->fwdJetPhi.at( i ) ) ); //PROBLEM: real value, trigger wants bin convert with phiINjetCoord
252 analysisDataFormat.Etajet.push_back( etaINjetCoord( inputNtuple.l1upgrade_->fwdJetEta.at( i ) ) ); //PROBLEM: real value, trigger wants bin convert with etaINjetCoord
253 analysisDataFormat.Taujet.push_back( false );
254 analysisDataFormat.isoTaujet.push_back( false );
255 analysisDataFormat.Fwdjet.push_back( true );
256
257 analysisDataFormat.Njet++;
258 }
259
260 // NOTES: Stage 1 has Tau Relaxed and TauIsolated. The isolated Tau are a subset of the Relaxed.
261 // so sort through the relaxed list and flag those that also appear in the isolated list.
262
263 for( unsigned int i=0; i<inputNtuple.l1upgrade_->nTau; i++ )
264 {
265
266 // remove duplicates
267 bool duplicate=false;
268 for( unsigned int j=0; j<i; j++ )
269 {
270 if( inputNtuple.l1upgrade_->tauBx.at( i )==inputNtuple.l1upgrade_->tauBx.at( j )
271 && inputNtuple.l1upgrade_->tauEt.at( i )==inputNtuple.l1upgrade_->tauEt.at( j )
272 && inputNtuple.l1upgrade_->tauEta.at( i )==inputNtuple.l1upgrade_->tauEta.at( j )
273 && inputNtuple.l1upgrade_->tauPhi.at( i )==inputNtuple.l1upgrade_->tauPhi.at( j ) )
274 {
275 duplicate=true;
276 //printf("Duplicate jet found and removed \n");
277 }
278 }
279
280 if( !duplicate )
281 {
282 analysisDataFormat.Bxjet.push_back( inputNtuple.l1upgrade_->tauBx.at( i ) );
283 analysisDataFormat.Etjet.push_back( inputNtuple.l1upgrade_->tauEt.at( i ) );
284 analysisDataFormat.Phijet.push_back( phiINjetCoord( inputNtuple.l1upgrade_->tauPhi.at( i ) ) ); //PROBLEM: real value, trigger wants bin convert with phiINjetCoord
285 analysisDataFormat.Etajet.push_back( etaINjetCoord( inputNtuple.l1upgrade_->tauEta.at( i ) ) ); //PROBLEM: real value, trigger wants bin convert with etaINjetCoord
286 analysisDataFormat.Taujet.push_back( true );
287 analysisDataFormat.Fwdjet.push_back( false );
288
289 bool isolated=false;
290 bool fnd=false;
291 unsigned int isoTau=0;
292 while( !fnd && isoTau < inputNtuple.l1upgrade_->nIsoTau )
293 {
294 if( inputNtuple.l1upgrade_->isoTauPhi.at( isoTau )==inputNtuple.l1upgrade_->tauPhi.at( i )
295 && inputNtuple.l1upgrade_->isoTauEta.at( isoTau )==inputNtuple.l1upgrade_->tauEta.at( i ) )
296 {
297 isolated=true;
298 fnd=true;
299 }
300 isoTau++;
301 }
302 analysisDataFormat.isoTaujet.push_back( isolated );
303
304 analysisDataFormat.Njet++;
305 } // duplicate check
306 }
307
308 // Fill energy sums (Are overflow flags accessible in l1extra?)
309 for( unsigned int i=0; i<inputNtuple.l1upgrade_->nMet; i++ )
310 {
311 //if(inputNtuple.l1upgrade_->metBx.at(i)==0) {
312 analysisDataFormat.ETT=inputNtuple.l1upgrade_->et.at( i );
313 analysisDataFormat.ETM=inputNtuple.l1upgrade_->met.at( i );
314 analysisDataFormat.PhiETM=inputNtuple.l1upgrade_->metPhi.at( i );
315 }
316 analysisDataFormat.OvETT=0; //not available in l1extra
317 analysisDataFormat.OvETM=0; //not available in l1extra
318
319 for( unsigned int i=0; i<inputNtuple.l1upgrade_->nMht; i++ )
320 {
321 if( inputNtuple.l1upgrade_->mhtBx.at( i )==0 )
322 {
323 analysisDataFormat.HTT=calculateHTT( analysisDataFormat ); //inputNtuple.l1upgrade_->ht.at(i) ;
324 analysisDataFormat.HTM=calculateHTM( analysisDataFormat ); //inputNtuple.l1upgrade_->mht.at(i) ;
325 analysisDataFormat.PhiHTM=0.; //inputNtuple.l1upgrade_->mhtPhi.at(i) ;
326 }
327 }
328 analysisDataFormat.OvHTM=0; //not available in l1extra
329 analysisDataFormat.OvHTT=0; //not available in l1extra
330
331 // Get the muon information from reEmul GMT
332 for( int i=0; i<inputNtuple.gmtEmu_->N; i++ )
333 {
334
335 analysisDataFormat.Bxmu.push_back( inputNtuple.gmtEmu_->CandBx[i] );
336 analysisDataFormat.Ptmu.push_back( inputNtuple.gmtEmu_->Pt[i] );
337 analysisDataFormat.Phimu.push_back( inputNtuple.gmtEmu_->Phi[i] );
338 analysisDataFormat.Etamu.push_back( inputNtuple.gmtEmu_->Eta[i] );
339 analysisDataFormat.Qualmu.push_back( inputNtuple.gmtEmu_->Qual[i] );
340 analysisDataFormat.Isomu.push_back( false );
341 analysisDataFormat.Nmu++;
342 }
343
344 break;
345
346 default:
347 throw std::runtime_error( "---Not a valid input source FULL STOP! " );
348
349 break;
350 }
351
352 return;
353 }
354
355 void l1menu::FullSamplePrivateMembers::fillL1Bits()
356 {
357 bool* PhysicsBits=currentEvent.physicsBits();
358
359 // I really don't think this if statement is correct. Surely it
360 // should be "if( gt_ )"? - M. Grimes.
361 if( !inputNtuple.gt_ )
362 {
363 for( Int_t ibit=0; ibit<128; ibit++ )
364 {
365 PhysicsBits[ibit]=0;
366 if( ibit<64 )
367 {
368 PhysicsBits[ibit]=(inputNtuple.gt_->tw1[2]>>ibit)&1;
369 }
370 else
371 {
372 PhysicsBits[ibit]=(inputNtuple.gt_->tw2[2]>>(ibit-64))&1;
373 }
374 }
375 }
376 else
377 {
378 PhysicsBits[0]=1; //set zero bias on if no gt information
379 }
380 }
381
382
383 l1menu::FullSample::FullSample()
384 : pImple_( new FullSamplePrivateMembers( this ) )
385 {
386 // No operation besides the initialiser list
387 }
388
389 l1menu::FullSample::~FullSample()
390 {
391 delete pImple_;
392 }
393
394 l1menu::FullSample::FullSample( const l1menu::FullSample& otherFullSample )
395 : pImple_( new FullSamplePrivateMembers(*otherFullSample.pImple_) )
396 {
397 // No operation besides the initialiser list
398 }
399
400 l1menu::FullSample::FullSample( l1menu::FullSample&& otherFullSample ) noexcept
401 : pImple_( otherFullSample.pImple_ )
402 {
403 otherFullSample.pImple_=NULL;
404 }
405
406 l1menu::FullSample& l1menu::FullSample::operator=( const l1menu::FullSample& otherFullSample )
407 {
408 *pImple_=*otherFullSample.pImple_;
409 return *this;
410 }
411
412 l1menu::FullSample& l1menu::FullSample::operator=( l1menu::FullSample&& otherFullSample ) noexcept
413 {
414 pImple_=otherFullSample.pImple_;
415 otherFullSample.pImple_=NULL;
416 return *this;
417 }
418
419 void l1menu::FullSample::loadFile( const std::string& filename )
420 {
421 pImple_->sumOfWeights=-1;
422 pImple_->inputNtuple.Open( filename );
423 }
424
425 void l1menu::FullSample::loadFilesFromList( const std::string& filenameOfList )
426 {
427 pImple_->sumOfWeights=-1;
428 pImple_->inputNtuple.OpenWithList( filenameOfList );
429 }
430
431 const l1menu::L1TriggerDPGEvent& l1menu::FullSample::getFullEvent( size_t eventNumber ) const
432 {
433 // Make sure the event number requested is valid. Use static_cast to get rid
434 // of the "comparison between signed and unsigned" compiler warning.
435 if( eventNumber>static_cast<size_t>(pImple_->inputNtuple.GetEntries()) ) throw std::runtime_error( "Requested event number is out of range" );
436
437 pImple_->inputNtuple.LoadTree(eventNumber);
438 pImple_->inputNtuple.GetEntry(eventNumber);
439 // This next call fills pImple_->currentEvent with the information in pImple_->inputNtuple
440 pImple_->fillDataStructure( 22 );
441 pImple_->fillL1Bits();
442
443 return pImple_->currentEvent;
444 }
445
446 size_t l1menu::FullSample::numberOfEvents() const
447 {
448 return static_cast<size_t>( pImple_->inputNtuple.GetEntries() );
449 }
450
451 const l1menu::IEvent& l1menu::FullSample::getEvent( size_t eventNumber ) const
452 {
453 // This returns a derived class so just delegate to that
454 return getFullEvent( eventNumber );
455 }
456
457 std::unique_ptr<l1menu::ICachedTrigger> l1menu::FullSample::createCachedTrigger( const l1menu::ITrigger& trigger ) const
458 {
459 return std::unique_ptr<l1menu::ICachedTrigger>( new CachedTriggerImplementation(trigger) );
460 }
461
462 float l1menu::FullSample::eventRate() const
463 {
464 return pImple_->eventRate;
465 }
466
467 void l1menu::FullSample::setEventRate( float rate )
468 {
469 pImple_->eventRate=rate;
470 }
471
472 float l1menu::FullSample::sumOfWeights() const
473 {
474 if( pImple_->sumOfWeights==-1 )
475 {
476 pImple_->sumOfWeights=0;
477 for( int eventNumber=0; eventNumber<pImple_->inputNtuple.GetEntries(); ++eventNumber )
478 {
479 pImple_->inputNtuple.LoadTree(eventNumber);
480 pImple_->inputNtuple.GetEntry(eventNumber);
481 pImple_->sumOfWeights+=pImple_->inputNtuple.event_->puWeight;
482 }
483 }
484
485 return pImple_->sumOfWeights;
486 }
487
488 std::unique_ptr<const l1menu::IMenuRate> l1menu::FullSample::rate( const l1menu::TriggerMenu& menu ) const
489 {
490 return std::unique_ptr<const l1menu::IMenuRate>();
491 }