ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/grimes/L1Menu/src/MenuSample.cpp
Revision: 1.7
Committed: Fri Jun 28 14:30:08 2013 UTC (11 years, 10 months ago) by grimes
Branch: MAIN
Changes since 1.6: +3 -29 lines
Log Message:
Added ICachedTrigger to speed up processing

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