ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/grimes/L1Menu/src/MenuSample.cpp
Revision: 1.4
Committed: Tue May 28 23:14:03 2013 UTC (11 years, 11 months ago) by grimes
Branch: MAIN
Changes since 1.3: +127 -55 lines
Log Message:
Numerous changes

File Contents

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