ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/grimes/L1Menu/src/MenuSample.cpp
Revision: 1.1
Committed: Sun May 5 20:37:32 2013 UTC (12 years ago) by grimes
Branch: MAIN
Log Message:
First draft of the L1 Trigger rate estimation code

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