ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/grimes/L1Menu/src/FullSample.cpp
Revision: 1.4
Committed: Fri Jul 5 14:21:56 2013 UTC (11 years, 10 months ago) by grimes
Branch: MAIN
Changes since 1.3: +1 -0 lines
Log Message:
Set the event weight in FullSamples properly

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