ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/grimes/L1Menu/src/MenuSample.cpp
Revision: 1.6
Committed: Thu Jun 6 15:05:04 2013 UTC (11 years, 11 months ago) by grimes
Branch: MAIN
Changes since 1.5: +8 -2 lines
Log Message:
Improved memory performance. Had to copy L1UpgradeNtuple.C into this package to fix some leaks.

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.6
8     #include "l1menu/IEvent.h"
9     #include "L1UpgradeNtuple.h"
10 grimes 1.2 #include "UserCode/L1TriggerUpgrade/interface/L1AnalysisDataFormat.h"
11 grimes 1.6 #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 grimes 1.4
17     // I'll use the unnamed namespace for things that are only used in this file
18     namespace
19     {
20     /** Implementation of the l1menu::IEvent interface.
21     *
22     * @author Mark Grimes (mark.grimes@bristol.ac.uk)
23     * @date 21/May/2013
24     */
25     class EventImplementation : public l1menu::IEvent
26     {
27     public:
28     //
29     // These are the methods required by the IEvent interface
30     //
31     virtual L1Analysis::L1AnalysisDataFormat& rawEvent() { return rawEvent_; }
32     virtual const L1Analysis::L1AnalysisDataFormat& rawEvent() const { return rawEvent_; }
33     virtual bool* physicsBits() { return physicsBits_; }
34     virtual const bool* physicsBits() const { return physicsBits_; }
35     virtual float weight() const { return 1; }
36     protected:
37     L1Analysis::L1AnalysisDataFormat rawEvent_;
38     bool physicsBits_[128];
39     }; // end of the EventImplementation class
40    
41     } // end of the unnamed namespace
42 grimes 1.1
43     namespace l1menu
44     {
45     class MenuSamplePrivateMembers
46     {
47     private:
48     static const size_t PHIBINS;
49     static const double PHIBIN[];
50     static const size_t ETABINS;
51     static const double ETABIN[];
52    
53     double degree( double radian );
54     int phiINjetCoord( double phi );
55     int etaINjetCoord( double eta );
56     double calculateHTT( const L1Analysis::L1AnalysisDataFormat& event );
57     double calculateHTM( const L1Analysis::L1AnalysisDataFormat& event );
58     public:
59     void fillDataStructure( int selectDataInput );
60 grimes 1.4 void fillL1Bits();
61 grimes 1.1 L1UpgradeNtuple inputNtuple;
62 grimes 1.4 ::EventImplementation currentEvent;
63 grimes 1.1 };
64     }
65    
66     const size_t l1menu::MenuSamplePrivateMembers::PHIBINS=18;
67     const double l1menu::MenuSamplePrivateMembers::PHIBIN[]={10,30,50,70,90,110,130,150,170,190,210,230,250,270,290,310,330,350};
68     const size_t l1menu::MenuSamplePrivateMembers::ETABINS=23;
69     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.};
70    
71    
72     double l1menu::MenuSamplePrivateMembers::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::MenuSamplePrivateMembers::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::MenuSamplePrivateMembers::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::MenuSamplePrivateMembers::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 )
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::MenuSamplePrivateMembers::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 )
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::MenuSamplePrivateMembers::fillDataStructure( int selectDataInput )
150     {
151 grimes 1.4 // Use a reference for ease of use
152     L1Analysis::L1AnalysisDataFormat& analysisDataFormat=currentEvent.rawEvent();
153    
154     analysisDataFormat.Reset();
155 grimes 1.1
156     // Grab standard event information
157 grimes 1.4 analysisDataFormat.Run=inputNtuple.event_->run;
158     analysisDataFormat.LS=inputNtuple.event_->lumi;
159     analysisDataFormat.Event=inputNtuple.event_->event;
160 grimes 1.1
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 grimes 1.4 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 grimes 1.1
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 grimes 1.4 analysisDataFormat.Isoel.push_back( isolated );
205     analysisDataFormat.Nele++;
206 grimes 1.1 }
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 grimes 1.4 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 grimes 1.1
237 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));
238 grimes 1.1 // 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 grimes 1.4 (fabs( inputNtuple.l1upgrade_->jetEta.at( i ) )>=3.0) ? analysisDataFormat.Fwdjet.push_back( true ) : analysisDataFormat.Fwdjet.push_back( false );
241 grimes 1.1
242 grimes 1.4 analysisDataFormat.Njet++;
243 grimes 1.1 }
244     }
245    
246     for( unsigned int i=0; i<inputNtuple.l1upgrade_->nFwdJets; i++ )
247     {
248    
249 grimes 1.4 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 grimes 1.1
257 grimes 1.4 analysisDataFormat.Njet++;
258 grimes 1.1 }
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 grimes 1.4 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 grimes 1.1
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 grimes 1.4 analysisDataFormat.isoTaujet.push_back( isolated );
303 grimes 1.1
304 grimes 1.4 analysisDataFormat.Njet++;
305 grimes 1.1 } // 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 grimes 1.4 analysisDataFormat.ETT=inputNtuple.l1upgrade_->et.at( i );
313     analysisDataFormat.ETM=inputNtuple.l1upgrade_->met.at( i );
314     analysisDataFormat.PhiETM=inputNtuple.l1upgrade_->metPhi.at( i );
315 grimes 1.1 }
316 grimes 1.4 analysisDataFormat.OvETT=0; //not available in l1extra
317     analysisDataFormat.OvETM=0; //not available in l1extra
318 grimes 1.1
319     for( unsigned int i=0; i<inputNtuple.l1upgrade_->nMht; i++ )
320     {
321     if( inputNtuple.l1upgrade_->mhtBx.at( i )==0 )
322     {
323 grimes 1.4 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 grimes 1.1 }
327     }
328 grimes 1.4 analysisDataFormat.OvHTM=0; //not available in l1extra
329     analysisDataFormat.OvHTT=0; //not available in l1extra
330 grimes 1.1
331     // Get the muon information from reEmul GMT
332     for( int i=0; i<inputNtuple.gmtEmu_->N; i++ )
333     {
334    
335 grimes 1.4 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 grimes 1.5 analysisDataFormat.Nmu++;
342 grimes 1.1 }
343    
344     break;
345    
346     default:
347 grimes 1.5 throw std::runtime_error( "---Not a valid input source FULL STOP! " );
348 grimes 1.1
349     break;
350     }
351    
352     return;
353     }
354    
355 grimes 1.4 void l1menu::MenuSamplePrivateMembers::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 grimes 1.1 l1menu::MenuSample::MenuSample()
384     : pImple_( new MenuSamplePrivateMembers )
385     {
386     // No operation besides the initialiser list
387     }
388    
389     l1menu::MenuSample::~MenuSample()
390     {
391     delete pImple_;
392     }
393    
394     l1menu::MenuSample::MenuSample( const l1menu::MenuSample& otherMenuSample )
395     : pImple_( new MenuSamplePrivateMembers(*otherMenuSample.pImple_) )
396     {
397     // No operation besides the initialiser list
398 grimes 1.4 }
399    
400     l1menu::MenuSample::MenuSample( l1menu::MenuSample&& otherMenuSample ) noexcept
401     : pImple_( otherMenuSample.pImple_ )
402     {
403     otherMenuSample.pImple_=NULL;
404 grimes 1.1 }
405    
406     l1menu::MenuSample& l1menu::MenuSample::operator=( const l1menu::MenuSample& otherMenuSample )
407     {
408     *pImple_=*otherMenuSample.pImple_;
409     return *this;
410     }
411    
412 grimes 1.4 l1menu::MenuSample& l1menu::MenuSample::operator=( l1menu::MenuSample&& otherMenuSample ) noexcept
413     {
414     pImple_=otherMenuSample.pImple_;
415     otherMenuSample.pImple_=NULL;
416     return *this;
417     }
418    
419 grimes 1.1 void l1menu::MenuSample::loadFile( const std::string& filename )
420     {
421     pImple_->inputNtuple.Open( filename );
422     }
423    
424     size_t l1menu::MenuSample::numberOfEvents() const
425     {
426     return static_cast<size_t>( pImple_->inputNtuple.GetEntries() );
427     }
428    
429 grimes 1.4 const l1menu::IEvent& l1menu::MenuSample::getEvent( size_t eventNumber ) const
430 grimes 1.1 {
431     // Make sure the event number requested is valid. Use static_cast to get rid
432     // of the "comparison between signed and unsigned" compiler warning.
433     if( eventNumber>static_cast<size_t>(pImple_->inputNtuple.GetEntries()) ) throw std::runtime_error( "Requested event number is out of range" );
434    
435     pImple_->inputNtuple.LoadTree(eventNumber);
436     pImple_->inputNtuple.GetEntry(eventNumber);
437     // This next call fills pImple_->currentEvent with the information in pImple_->inputNtuple
438     pImple_->fillDataStructure( 22 );
439 grimes 1.4 pImple_->fillL1Bits();
440 grimes 1.1
441     return pImple_->currentEvent;
442     }