ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/grimes/L1Menu/src/MenuSample.cpp
Revision: 1.7
Committed: Fri Jun 28 14:30:08 2013 UTC (11 years, 10 months ago) by grimes
Branch: MAIN
Changes since 1.6: +3 -29 lines
Log Message:
Added ICachedTrigger to speed up processing

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