ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/grimes/L1Menu/src/MenuSample.cpp
Revision: 1.5
Committed: Wed Jun 5 11:02:36 2013 UTC (11 years, 11 months ago) by grimes
Branch: MAIN
Changes since 1.4: +2 -2 lines
Log Message:
Fixes to the number of muons (so muon triggers now work); error when loading ReducedMenuSamples from disk; made L1_SingleIsoMu an alias for L1_SingleMu since there's no isolation on muons and they're the same

File Contents

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