ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/grimes/L1Menu/src/MenuSample.cpp
Revision: 1.3
Committed: Wed May 15 10:41:36 2013 UTC (11 years, 11 months ago) by grimes
Branch: MAIN
Changes since 1.2: +2 -0 lines
Log Message:
Got everything pretty much working for a scram build

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