ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/grimes/L1Menu/src/MenuSample.cpp
Revision: 1.2
Committed: Mon May 6 00:05:38 2013 UTC (12 years ago) by grimes
Branch: MAIN
Changes since 1.1: +2 -1 lines
Log Message:
Fixed a bug where I hadn't fully changed a namespace name

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