ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/grimes/L1Menu/src/FullSample.cpp
Revision: 1.5
Committed: Wed Jul 24 11:48:55 2013 UTC (11 years, 9 months ago) by grimes
Branch: MAIN
CVS Tags: HEAD
Changes since 1.4: +2 -1 lines
Log Message:
Big commit of lots of changes before migrating to Git.

File Contents

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