Report problems to ATLAS LXR Team (with time and IP address indicated)

The LXR Cross Referencer
[ source navigation ]
[ diff markup ]
[ identifier search ]
[ general search ]
 
 
Architecture: [ linux ]
Version: [ head ] [ nightly ] [ GaudiDev ]
  Links to LXR source navigation pages for stable releases [ 12.*.* ]   [ 13.*.* ]   [ 14.*.* ]   [ 15.*.* ]  
001 // ================================================
002 // EventHeaderMaker class Implementation
003 // ================================================
004 //
005 // THIS TEXT TO BE REPLACED BY ATLAS STANDARD FORMAT
006 //
007 // Namespace Atlfast
008 //
009 
010 #include "AtlfastAlgs/EventHeaderMaker.h"
011 #include "AtlfastAlgs/MissingMomentum.h"
012 #include "AtlfastAlgs/GlobalEventData.h"
013 
014 #include "AtlfastEvent/EventHeader.h"
015 #include "AtlfastEvent/IKinematic.h"
016 #include "AtlfastEvent/MsgStreamDefs.h"
017 #include "AtlfastEvent/ReconstructedParticle.h"
018 #include "AtlfastEvent/Cell.h"
019 #include "AtlfastEvent/Cluster.h"
020 #include "AtlfastEvent/KtCluster.h"
021 #include "AtlfastEvent/ParticleCodes.h"
022 #include "AtlfastEvent/CollectionDefs.h"
023 #include "AtlfastEvent/MCparticleCollection.h"
024 #include "AtlfastEvent/ReconstructedParticle.h"
025 #include "AtlfastEvent/Jet.h"
026 
027 #include "AtlfastUtils/HeaderPrinter.h"
028 #include "AtlfastUtils/TesIO.h"
029 
030 // standard library includes
031 #include <cmath> 
032 #include <string>
033 #include <vector>
034 #include <algorithm>
035 #include <assert.h>
036 
037 // Gaudi includes
038 #include "GaudiKernel/DataSvc.h"
039 #include "GaudiKernel/ISvcLocator.h"
040 #include "GaudiKernel/MsgStream.h"
041 
042 
043 //Generator includes
044 #include "GeneratorObjects/McEventCollection.h"
045 
046 
047 #include "CLHEP/Units/SystemOfUnits.h"
048 
049 /** @brief Makes the Atlfast Event Header and computes missing momentum.
050    *
051    *  This must be run as the final algorithm in the
052    *  Atlfast suite as it requires the information from
053    *  all the earlier makers.
054    *
055    * @image html EventHeaderMaker.jpg "EventHeaderMaker sequence"
056    */
057 
058 namespace Atlfast {
059 
060   //--------------------------------
061   // Constructors and destructors
062   //-------------------------------
063   EventHeaderMaker::EventHeaderMaker 
064   ( const std::string& name, ISvcLocator* pSvcLocator ) 
065     : Algorithm( name, pSvcLocator ),
066       m_tesIO(0),
067       m_visibleToAtlas(0)
068   {
069     
070     // Default paths for entities in the TES
071     m_electronLocation         = "/Event/AtlfastIsolatedElectrons";
072     m_photonLocation           = "/Event/AtlfastIsolatedPhotons";
073     m_isolatedMuonLocation     = "/Event/AtlfastIsolatedMuons";
074     m_nonIsolatedMuonLocation  = "/Event/AtlfastNonIsolatedMuons";
075     m_jetLocation              = "/Event/AtlfastJets";
076     m_cellLocation             = "/Event/AtlfastCells";
077     m_clusterLocation          = "/Event/AtlfastClusters";
078     m_outputLocation           = "/Event/AtlfastEventHeader";
079     m_missingMomentumLocation  = "/Event/AtlfastMissingMomentum";
080     m_beamEnergy               = 7000.0*GeV;
081     m_testMode                 = false;
082     //m_addHalos                 = false;
083     m_adjustMissETForIsolation = true;
084     
085     declareProperty( "MissingMomentumLocation", m_missingMomentumLocation ) ;
086     declareProperty( "ElectronLocation",        m_electronLocation ) ;
087     declareProperty( "PhotonLocation",          m_photonLocation );
088     declareProperty( "IsolatedMuonLocation",    m_isolatedMuonLocation );
089     declareProperty( "NonIsolatedMuonLocation", m_nonIsolatedMuonLocation );
090     declareProperty( "JetLocation",             m_jetLocation );
091     declareProperty( "OutputLocation",          m_outputLocation );
092     declareProperty( "BeamEnergy",              m_beamEnergy );
093     //declareProperty( "AddHalos",                m_addHalos );
094 
095     declareProperty( "TestMode",                m_testMode ) ;
096     
097     
098     
099   }
100   
101   
102   //-------------------------------------
103   // Destructor 
104   //-------------------------------------
105   EventHeaderMaker::~EventHeaderMaker(){
106     if (m_tesIO) delete m_tesIO;
107   }
108 
109   //---------------------------------
110   // initialise() 
111   //---------------------------------
112   StatusCode EventHeaderMaker::initialize()
113   {
114     // .. put any initialisation in here...
115     MsgStream log(messageService(), name());
116     log << MSG::DEBUG << "in initialize()" << endreq;
117     
118     //set up the list of particles which are considered as "escaping"
119     m_escapingParticles.push_back(ParticleCodes::NU_E);
120     m_escapingParticles.push_back(ParticleCodes::NU_MU);
121     m_escapingParticles.push_back(ParticleCodes::NU_TAU);
122 
123 
124     //get the Global Event Data from singleton pattern
125     GlobalEventData* ged = GlobalEventData::Instance();
126     m_visibleToAtlas           = ged ->visibleToAtlas();
127     m_lumi                     = ged -> lumi();
128     m_barrelForwardEta         = ged -> barrelForwardEta();
129     m_mcLocation               = ged -> mcLocation();
130     m_adjustMissETForIsolation = ged->adjustMissETForIsolation();
131 
132     m_tesIO = new TesIO(m_mcLocation, ged->justHardScatter());
133 
134     HeaderPrinter hp("Atlfast EventHeader Maker:", log);
135     hp.add("TestMode                    ", m_testMode);
136     hp.add("EndCap Begins at (eta)      ", m_barrelForwardEta);
137     hp.add("Luminosity                  ", m_lumi);
138     hp.add("Beam Energy                 ", m_beamEnergy);
139     hp.add("MC location                 ", m_mcLocation);
140     hp.add("Electron Location           ", m_electronLocation);
141     hp.add("Photon   Location           ", m_photonLocation);
142     hp.add("Isolated Muon Location      ", m_isolatedMuonLocation);
143     hp.add("Non-Isolated Muon Location  ", m_nonIsolatedMuonLocation);
144     hp.add("Jet Location                ", m_jetLocation);
145     hp.add("Cell Location               ", m_cellLocation);
146     hp.add("Cluster Location            ", m_clusterLocation);
147     hp.add("Unused cell+cluster sum     ", m_missingMomentumLocation);
148     hp.add("Output Location             ", m_outputLocation);    
149     hp.print();
150     
151 
152     
153     
154     return StatusCode::SUCCESS ;
155   }
156   
157   //---------------------------------
158   // finalise() 
159   //---------------------------------
160   StatusCode EventHeaderMaker::finalize()
161   {
162     // .. put any finalisation in here...
163     MsgStream log( messageService(), name() ) ;
164     log << MSG::INFO << "finalizing " << endreq;
165     
166     return StatusCode::SUCCESS ;
167   }
168   
169   //----------------------------------------------
170   // execute() method called once per event
171   //----------------------------------------------
172   StatusCode EventHeaderMaker::execute( )
173   {
174     
175     //................................
176     // make a message logging stream
177     
178     MsgStream log( messageService(), name() ) ;
179     log << MSG::DEBUG << "Executing " << endreq;
180     std::string mess;
181     
182       
183 
184     //
185     // Check each particle list in turn and get the number of entries
186     //
187 
188    /** numberInList - finds out how many entries are in a TES list */
189 
190     int nElectrons         = 
191       numberInList<ReconstructedParticleCollection>(m_electronLocation);
192     int nPhotons           = 
193       numberInList<ReconstructedParticleCollection>(m_photonLocation);
194     int nIsolatedMuons     = 
195       numberInList<ReconstructedParticleCollection>(m_isolatedMuonLocation);
196     int nNonIsolatedMuons  = 
197       numberInList<ReconstructedParticleCollection>(m_nonIsolatedMuonLocation);
198     int nJets              = numberInList<JetCollection>(m_jetLocation);
199     /** numberJetFlavor - finds out how many numbers of jets by flavor */
200     int nBJets             = numberJetFlavor(m_jetLocation,
201                                              ParticleCodes::BQUARK);
202     int nCJets             = numberJetFlavor(m_jetLocation,
203                                              ParticleCodes::CQUARK);
204     int nTauJets           = numberJetFlavor(m_jetLocation,
205                                              ParticleCodes::TAU);
206     nTauJets              += numberJetFlavor(m_jetLocation,
207                                              -ParticleCodes::TAU);
208     
209     log << MSG::DEBUG << "Electrons: " << nElectrons << endreq;
210     log << MSG::DEBUG << "Photons: " << nPhotons << endreq;
211     log << MSG::DEBUG << "IsolatedMuons: " << nIsolatedMuons << endreq;
212     log << MSG::DEBUG << "NonIsolatedMuons: " << nNonIsolatedMuons << endreq;
213     log << MSG::DEBUG << "Jets: " << nJets << endreq;
214     
215     
216     // Event Shape variable calculations
217     
218     double jetCircularity = 0.0; //FIXME not written yet
219     double eventCircularity = 0.0; //FIXME not written yet
220     double thrust = 0.0; //FIXME not written yet
221     double oblateness = 0.0; //FIXME not written yet
222     
223     // Cumulative ET sum
224     double sumET = 0.0;
225 
226     //missing momentum and sumET
227 
228     /** calculate the reconstructed missing momentum */
229     HepLorentzVector pMiss = missingMomentum(log,sumET);
230     
231     log << MSG::DEBUG << "Missing momentum " <<pMiss<<endreq; 
232     
233     /** calculated the momentum of escaping particles from MC truth */
234     HepLorentzVector pEscaped = escapedMomentum(log);
235     
236     log << MSG::DEBUG << "Escaped momentum " << pEscaped<<endreq;
237 
238     //get a container of weight containers for this event
239     MCweightContainerCollection weightCollections;
240     TesIoStat stat = m_tesIO->getMCweightContainers(weightCollections);
241     if(!stat){
242       log << MSG::ERROR <<"Problem getting b-quarks" << endreq; 
243       return stat;
244     }
245 
246     EventHeader* header = 
247       new EventHeader(
248                       nElectrons, 
249                       nIsolatedMuons, 
250                       nNonIsolatedMuons, 
251                       nPhotons, 
252                       nJets,
253                       nBJets, 
254                       nCJets, 
255                       nTauJets, 
256                       jetCircularity, 
257                       eventCircularity,
258                       thrust, 
259                       oblateness, 
260                       pMiss,
261                       sumET,
262                       pEscaped,
263                       weightCollections
264                       );
265     
266     stat = m_tesIO->store(header,m_outputLocation) ;
267     std::string message = stat? "Stored Evt Header":"Failed to Evt Header";
268     log<<MSG::DEBUG<<message<<endreq;
269     return stat.status();
270     
271   }
272   
273   //-------------------------------------------------------
274   // missingMomentum
275   //-------------------------------------------------------
276   HepLorentzVector EventHeaderMaker::missingMomentum(MsgStream& log, double &sumET) {
277 
278     HepLorentzVector unused(0., 0., 0., 0.);
279     const MissingMomentum* mm(0);
280     TesIoStat stat = m_tesIO->getDH(mm,m_missingMomentumLocation );
281     if(!stat){
282       log << MSG::WARNING<<"Could Not Find Missing momentum in TES"<<endreq;
283       return unused;
284     }
285     log << MSG::DEBUG<<"Found Missing Momentum in TES"<<endreq;
286 
287     unused = mm->momentum();
288 
289     HepLorentzVector temp(0.0,0.0,0.0,0.0);
290 
291     temp = unused
292       +  momentumSum<ReconstructedParticleCollection> (m_electronLocation)
293       +  momentumSum<ReconstructedParticleCollection> (m_photonLocation)
294       +  momentumSum<ReconstructedParticleCollection> (m_isolatedMuonLocation)
295       +  momentumSum<JetCollection>                   (m_jetLocation)
296       +  momentumSum<ReconstructedParticleCollection> (m_nonIsolatedMuonLocation)
297       -  assocMomentumSum<JetCollection, ReconstructedParticle> (m_jetLocation);
298 
299     
300     double unused_sumET = mm->sumET();
301     sumET = unused_sumET +
302       +  transverseEnergySum<ReconstructedParticleCollection> (m_electronLocation)
303       +  transverseEnergySum<ReconstructedParticleCollection> (m_photonLocation)
304       +  transverseEnergySum<ReconstructedParticleCollection> (m_isolatedMuonLocation)
305       +  transverseEnergySum<JetCollection>                   (m_jetLocation)
306       +  transverseEnergySum<ReconstructedParticleCollection> (m_nonIsolatedMuonLocation)
307       -  assocTransverseEnergySum<JetCollection, ReconstructedParticle> (m_jetLocation);
308 
309     
310     // Add halos if requested and if the missing ET has not already been adjusted
311     // to treat isolation problem
312     /*if (m_addHalos && !m_adjustMissETForIsolation){
313       log << MSG::DEBUG <<"Adding halo values to missing ET, sumET" << endreq; 
314       temp += haloVectorSum(m_electronLocation) + haloVectorSum(m_photonLocation);
315       sumET += haloSum(m_electronLocation) + haloSum(m_photonLocation);
316     }*/ 
317     
318     log << MSG::DEBUG<<"Total Imbalanced momentum = " << temp 
319         << ", total sum of scalar ETs = " << sumET << endreq;
320     
321     if(m_testMode){
322       /** returns the summed momentum of a list of IKinematic objects
323         *  in the TES */
324       HepLorentzVector totCluster = 
325         momentumSum<IClusterCollection> (m_clusterLocation);
326       /** assocMomentumSum - returns the summed momentum of a list of associated
327        * IKinematic objects in the TES */
328 
329       HepLorentzVector clusterAssE = 
330         assocMomentumSum<ReconstructedParticleCollection, Cluster>
331           (m_electronLocation);
332         HepLorentzVector clusterAssE2 = 
333           assocMomentumSum<ReconstructedParticleCollection, KtCluster>
334           (m_electronLocation);
335 
336         HepLorentzVector clusterAssP = 
337           assocMomentumSum<ReconstructedParticleCollection, Cluster>
338           (m_photonLocation);
339         HepLorentzVector clusterAssP2 = 
340           assocMomentumSum<ReconstructedParticleCollection, KtCluster>
341           (m_photonLocation);
342 
343         HepLorentzVector clusterAssM = 
344           assocMomentumSum<ReconstructedParticleCollection, Cluster>
345           (m_isolatedMuonLocation);
346         HepLorentzVector clusterAssM2 = 
347           assocMomentumSum<ReconstructedParticleCollection, KtCluster>
348           (m_isolatedMuonLocation);
349 
350         HepLorentzVector clusterAssJ = 
351           assocMomentumSum<JetCollection, Cluster>(m_jetLocation);
352         HepLorentzVector clusterAssJ2 = 
353           assocMomentumSum<JetCollection, KtCluster>(m_jetLocation);
354 
355       HepLorentzVector clusterDiff = totCluster-
356         clusterAssJ-
357         clusterAssE-
358         clusterAssM-
359         clusterAssP-
360         clusterAssJ2-
361         clusterAssE2-
362         clusterAssM2-
363         clusterAssP2;
364 
365       HepLorentzVector totCell = 
366         momentumSum<ITwoCptCellCollection>(m_cellLocation);
367 
368       HepLorentzVector cellAssCluster = 
369         assocMomentumSum<IClusterCollection, TwoCptCell>(m_clusterLocation);
370 
371       HepLorentzVector cellDiff = totCell-cellAssCluster;
372 
373       HepLorentzVector diffS  = clusterDiff+cellDiff-unused;
374       
375 
376       
377       log<<MSG::DEBUG<<"Total cluster    "<<totCluster<<endreq;
378       log<<MSG::DEBUG<<"Clus Ass Jet     "<<clusterAssJ+clusterAssJ2<<endreq;
379       log<<MSG::DEBUG<<"Clus Ass El      "<<clusterAssE+clusterAssE2<<endreq;
380       log<<MSG::DEBUG<<"Clus Ass Mu      "<<clusterAssM+clusterAssM2<<endreq;
381       log<<MSG::DEBUG<<"Clus Ass Ph      "<<clusterAssP+clusterAssP2<<endreq;
382       log<<MSG::DEBUG<<"Clus diff        "<<clusterDiff<<endreq;
383       log<<MSG::DEBUG<<"Total Cell       "<<totCell<<endreq;
384       log<<MSG::DEBUG<<"Cell Ass Clus    "<<cellAssCluster<<endreq;
385       log<<MSG::DEBUG<<"Clus diff        "<<cellDiff<<endreq;
386       log<<MSG::DEBUG<<"Cell+Clus diff   "<<cellDiff<<endreq;
387       log<<MSG::DEBUG<<"Unused           "<<unused<<endreq;
388       log<<MSG::DEBUG<<"diff             "<<diffS<<endreq;
389 
390       log<<MSG::DEBUG<<endreq;
391       
392       assert( diffS.rho() <0.001);
393 
394 
395 
396 
397 
398       HepLorentzVector temp2 =  
399         momentumSum<ReconstructedParticleCollection>    (m_electronLocation)
400         +  momentumSum<ReconstructedParticleCollection> (m_photonLocation)
401         +  momentumSum<ReconstructedParticleCollection> (m_isolatedMuonLocation)
402         +  momentumSum<JetCollection>  (m_jetLocation)
403         +  momentumSum<IClusterCollection> (m_clusterLocation)
404         +  momentumSum<ITwoCptCellCollection> (m_cellLocation)
405         +  momentumSum<ReconstructedParticleCollection> (m_nonIsolatedMuonLocation)
406         
407         - assocMomentumSum<ReconstructedParticleCollection, Cluster>
408         (m_electronLocation)   
409         - assocMomentumSum<ReconstructedParticleCollection, KtCluster>
410         (m_electronLocation)   
411         
412         - assocMomentumSum<ReconstructedParticleCollection, Cluster>
413         (m_photonLocation)
414         - assocMomentumSum<ReconstructedParticleCollection, KtCluster>
415         (m_photonLocation)
416         
417         - assocMomentumSum<ReconstructedParticleCollection, Cluster>
418         (m_nonIsolatedMuonLocation)
419         - assocMomentumSum<ReconstructedParticleCollection, KtCluster>
420         (m_nonIsolatedMuonLocation)
421         
422         -  assocMomentumSum<JetCollection, Cluster>
423         (m_jetLocation)
424         -  assocMomentumSum<JetCollection, KtCluster>
425         (m_jetLocation)
426         
427         -  assocMomentumSum<JetCollection, ReconstructedParticle>
428         (m_jetLocation)
429 
430         - assocMomentumSum<IClusterCollection, Cell>
431         (m_clusterLocation);
432 
433       
434       assert( (temp-temp2).rho() <0.001);
435     }
436     
437     // calculate missing momentum from the reconstructed sum of momenta:
438     //
439     temp.setPx( -temp.px() );
440     temp.setPy( -temp.py() );
441     temp.setPz( -temp.pz() );
442     temp.setE(  2.0*m_beamEnergy - temp.e() );
443     
444     return temp;
445   }
446   
447   
448   
449   //-------------------------------------------------------
450   // escapedMomentum
451   //-------------------------------------------------------
452   HepLorentzVector EventHeaderMaker::escapedMomentum(MsgStream& log) {
453 
454     HepLorentzVector temp(0.0,0.0,0.0,0.0);
455 
456     MCparticleCollection  mcParticles;
457     TesIoStat stat = m_tesIO->getMC( mcParticles, m_visibleToAtlas ) ;
458     std::string mess = 
459       stat? "MC particles retrieved":"MC particles retrieve failed";
460     log<<MSG::DEBUG<<mess<<endreq;
461     
462 
463     MCparticleCollection::const_iterator part = mcParticles.begin();
464     for (;part != mcParticles.end(); ++part) {
465       temp.setPx(temp.px()+(*part)->momentum().px());
466       temp.setPy(temp.py()+(*part)->momentum().py());
467       temp.setPz(temp.pz()+(*part)->momentum().pz());
468       temp.setE(temp.e()+(*part)->momentum().e());
469     }
470 
471     return temp;
472   }
473 
474 } // end of namespace bracket
475 
476 
477 
478 
479 
480 
481 
482 
483 
484 
485 
486
[ source navigation ] [ diff markup ] [ identifier search ] [ general search ]
Due to the LXR bug, the updates fail sometimes to remove references to deleted files. The Saturday's full rebuilds fix these problems
This page was automatically generated by the LXR engine. Valid HTML 4.01!