---

XANAMuons.cpp

Go to the documentation of this file.

#include <Utilities/Configuration/interface/Architecture.h>

#include <XANADOO/XANAMuons/interface/XANAMuonCandidate.h>
#include <XANADOO/XANAMuons/interface/XANAMuonTrack.h>
#include <XANADOO/XANATracks/interface/XANATrack.h>
#include <XANADOO/XANAUtilities/interface/XANAToolkit.h>

#include "Trigger/L1GlobalMuonTrigger/interface/L1MuGMTSetup.h"
#include "Trigger/L1GlobalMuonTrigger/interface/L1MuGlobalMuonTrigger.h"
#include "Trigger/L1GlobalMuonTrigger/interface/L1MuGMTCand.h"

#include <CommonReco/PersistentTrack/interface/TTrack.h>
#include <CommonReco/PersistentTrack/interface/TrackReconstructor.h>
#include <CommonReco/PatternTools/interface/TrackFinder.h>
#include <TrackerReco/GtfPattern/interface/ModularKfReconstructor.h>
#include <TrackerReco/TkSeedGenerator/interface/CombinatorialSeedGeneratorFromPixel.h>
#include <TrackerReco/GtfPattern/interface/CombinatorialTrajectoryBuilder.h>
#include <CommonReco/TrackFitters/interface/KFFittingSmoother.h>
#include <CommonReco/PatternTools/interface/TrajectoryCleanerBySharedHits.h>

#include <CARF/G3Event/interface/G3EventProxy.h>
#include <CARF/G3Event/interface/G3SetUp.h>
#include <CARF/Reco/interface/AutoRecCollection.h>
#include <CARF/Reco/interface/RecItr.h>

#include <Utilities/Notification/interface/Observer.h>
#include <Utilities/Notification/interface/PackageInitializer.h>
#include <Utilities/UI/interface/SimpleConfigurable.h> 
#include "Utilities/Notification/interface/Singleton.h"
 
#include <CLHEP/Geometry/Point3D.h>
#include <vector>

struct MuonObserver : public Observer<G3EventProxy *>, public Observer<G3SetUp *> {

  int nbMuons_ ;
  double momentumSum_ ;
  SimpleConfigurable<bool> produceReference_;
  SimpleConfigurable<string> muonReconstructorName_;
  TrackReconstructor *muonReconstructor_;
 
  MuonObserver()
   : nbMuons_(0), momentumSum_(0), muonReconstructor_(0)
   { 
    // initialize observer
    Observer<G3EventProxy *>::init() ;      
    Observer<G3SetUp *>::init() ;
    
    // steering
    produceReference_=SimpleConfigurable<bool>(false,"XANAMuons:ProduceReference");
    muonReconstructorName_=
     SimpleConfigurable<string>("GlobalMuonReconstructor","XANAMuons:MuonReconstructor");
    // initialisation of muoning algo must be done in observing G3SetUp       
   }
   
  ~MuonObserver()
   {
    delete muonReconstructor_;
    xana::out(1)<<"[OVAL] number of muons : "<<nbMuons_<<endl ;
    xana::out(1)<<"[OVAL] mean energy : "<<momentumSum_/nbMuons_<<endl ;
   }

  // observing G3SetUp needed for traker initialisation
  // algos can also be initialised in an update(G3SetUp*) method if the o
  // observer is also observing G3SetUp. Then initialisation of RecUnits is
  // done a bit later in the execution, but god knows exactly when! 
  void upDate( G3SetUp *su )
   {    
    // defines muon reconstruction algo  
    // muon reconstructor initialisation can be done in the event loop as well
    L1MuGMTSetup* setup = Singleton<L1MuGMTSetup>::instance();

   }
    
  void upDate( G3EventProxy * ev )
   { 
    
    if (!ev) return ;    
    
    RecItr<TTrack> muon(ev->recEvent(),muonReconstructorName_.value().c_str());
    
    // to produce reference
    double momentum=0.;
    if (produceReference_) {
      while (muon.next()) {
        TrajectoryStateOnSurface ts = (*muon).impactPointState() ;
        if (ts.isValid()) {
          momentum = ts.globalMomentum().mag() ;
          nbMuons_++ ;
          momentumSum_ += momentum ;
          cout<<"[OVAL] new muon with momentum "<<momentum<<flush ;
          cout<<" and impact parameter "<<muon->impactParameter3D().value()<<endl ;
         }  
      }
      return ;
    } 
    
    // construct XANAmuons from Tmuons
    vector<XANAMuonCandidate> xmuons;
    XANAMuonCandidate *xmuon=0;
    while (muon.next()) { 
      TrajectoryStateOnSurface ts = (*muon).impactPointState() ;
      if (ts.isValid()) {
        momentum = ts.globalMomentum().mag() ;
        nbMuons_++ ;
        momentumSum_ += momentum ;
      }        
      cout << "Found new muon with momentum : "<<momentum<<endl ; 
      TrajectoryStateOnSurface tsi = (*muon).innermostState() ;
      TrajectoryStateOnSurface tso = (*muon).outermostState() ;
      HepVector3D momentumAtVertex(muon->momentumAtVertex().x(), 
       muon->momentumAtVertex().y(),muon->momentumAtVertex().z());
      HepPoint3D firstHit(tsi.globalPosition().x(),tsi.globalPosition().y(),
        tsi.globalPosition().z()); 
      HepVector3D momentumAtFirst(tsi.globalMomentum().x(),tsi.globalMomentum().y(),
        tsi.globalMomentum().z());
      HepPoint3D lastHit(tso.globalPosition().x(),tso.globalPosition().y(),
        tso.globalPosition().z()); 
      HepVector3D momentumAtLast(tso.globalMomentum().x(),tso.globalMomentum().y(),
        tso.globalMomentum().z());
      xmuon = new XANAMuonCandidate();
      HepLorentzVector momentum(muon->momentumAtVertex().x(), 
       muon->momentumAtVertex().y(),muon->momentumAtVertex().z(),
       muon->momentumAtVertex().mag());
      xmuon->setMomentum(momentum);
      xmuon->setPt(momentum.perp());
      xmuon->setCharge(muon->charge());
      XANAMuonTrack *xmuontrack = new XANAMuonTrack();
      xmuontrack->setCharge(muon->charge()); 
      xmuontrack->setChi2OverDof(muon->normalisedChiSquared());
      xmuontrack->setNumberOfHits(muon->foundHits());
      xmuontrack->setNumberOfLostHits(muon->lostHits());
      xmuontrack->setImpactParameter(muon->impactParameter3D().value());
      xmuontrack->setLongImpactParameter(muon->zImpactParameter().value());
      xmuontrack->setTransImpactParameter(muon->transverseImpactParameter().value());
      xmuontrack->setMomentumAtVertex(momentumAtVertex);
      xmuontrack->setPositionAtFirstPoint(firstHit);
      xmuontrack->setMomentumAtFirstPoint(momentumAtFirst);
      xmuontrack->setPositionAtLastPoint(lastHit);
      xmuontrack->setMomentumAtLastPoint(momentumAtLast);
      xmuontrack->setAlgoName(muonReconstructorName_.value());
      xmuon->setMuonTrack(xmuontrack);
      xmuons.push_back(*xmuon);
    }
     
    // read XANAmuons 
    vector<XANAMuonCandidate>::iterator it;
    for (it=xmuons.begin(); it!=xmuons.end(); it++)  
     {
      xana::out(1)<<"[OVAL] new muon with transverse momentum "<<it->getMuonTrack()->getMomentumAtVertex().perp()<<flush ;
      xana::out(1)<<" and impact parameter "<<it->getMuonTrack()->getImpactParameter()<<endl ;
     }
   }
      
} ;
 
static PKBuilder<MuonObserver> muonObserver("MuonObserver") ;

---