---

XANAVertex.cpp

Go to the documentation of this file.

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

#include <XANADOO/XANAVertex/interface/XANAVertex.h>
#include <XANADOO/XANATracks/interface/XANATrack.h>
#include <XANADOO/XANAUtilities/interface/XANAToolkit.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 "CommonReco/PatternTools/interface/RefittedRecTrack.h"

#include "Vertex/CARFVertex/interface/RecoVertex.h"

#include "CARF/Reco/interface/AutoRecCollection.h"
#include "CARF/G3Event/interface/EventObservers.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 <CLHEP/Geometry/Point3D.h>
#include <vector>

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

  int nbVertices_ ;
  SimpleConfigurable<string> vertexReconstructorName_;
  SimpleConfigurable<string> trackReconstructorName_;
  TrackReconstructor *trackReconstructor_;
 
  VertexObserver()
   : nbVertices_(0), trackReconstructor_(0)
   { 
    // initialize observer
    Observer<G3EventProxy *>::init() ;      
    Observer<G3SetUp *>::init() ;      

    // steering
    vertexReconstructorName_=
     SimpleConfigurable<string>("PVRVertex","XANAVertex:VertexReconstructor");
    trackReconstructorName_=
//     SimpleConfigurable<string>("FKfTracks","XANAVertex:TrackReconstructor");
     SimpleConfigurable<string>("FKfTracks","XANAVertex:TrackReconstructor");
    // initialisation of tracking algo must be done in observing G3SetUp       
   }
   
  ~VertexObserver()
   {
    delete trackReconstructor_;
    xana::out(1)<<"[OVAL] number of vertices : "<<nbVertices_<<endl ;
   }

  void upDate( G3SetUp *su )
   {    
    // defines Tracking algo  
    // beware! not to add to reconstructor! this is not a RecUnit 
    // this will use default propagator, trajectory builder, etc..
    TrackFinder *myTrackFinder = new ModularKfReconstructor(
     new CombinatorialSeedGeneratorFromPixel(0.8,0.1,15.),
     new CombinatorialTrajectoryBuilder,
     new KFFittingSmoother,
     new TrajectoryCleanerBySharedHits);
    trackReconstructor_ = 
     new TrackReconstructor(myTrackFinder,trackReconstructorName_.value().c_str()) ;
   }
    
  void upDate( G3EventProxy * ev )
   { 
    
    if (!ev) return ;    
  
    RecItr<TTrack> track(ev->recEvent(),trackReconstructorName_.value().c_str());
    while (track.next()) { 
      xana::out(1) << "new reconstructed track " << std::endl; 
      TrajectoryStateOnSurface ts = (*track).impactPointState() ;
      xana::out(1) << "(*track).impactPointState() " <<(*track).impactPointState() << std::endl;      
      TrajectoryStateOnSurface tsi = (*track).innermostState() ;
      xana::out(1) << "(*track).innermostState() " <<(*track).innermostState() << std::endl;      
      TrajectoryStateOnSurface tso = (*track).outermostState() ;
      xana::out(1) << "(*track).outermostState() " <<(*track).outermostState() << std::endl;      
      HepVector3D momentumAtVertex(track->momentumAtVertex().x(), 
       track->momentumAtVertex().y(),track->momentumAtVertex().z());
      xana::out(1) << "(*track).momentumAtVertex() " <<(*track).momentumAtVertex() << std::endl;      
      HepVector3D momentumAtFirst(tsi.globalMomentum().x(),tsi.globalMomentum().y(),
        tsi.globalMomentum().z());
      HepVector3D momentumAtLast(tso.globalMomentum().x(),tso.globalMomentum().y(),
        tso.globalMomentum().z());
   }    
    AutoRecCollection<RecoVertex, G3EventProxy *> vertices(vertexReconstructorName_.value().c_str());
        
    // construct XANAVertex from RecoVertex
    vector<XANAVertex> xVertices;
    XANAVertex *xVertex=0;
    XANATrack *xTrack=0;
    for(AutoRecCollection<RecoVertex, G3EventProxy *>::const_iterator iv=vertices.begin(); 
     iv!=vertices.end(); iv++) {
      HepPoint3D position((*iv)->position().x(),(*iv)->position().y(), 
       (*iv)->position().z()); 
      xVertex = new XANAVertex(position, (*iv)->totalChiSquared(), 
       (*iv)->degreesOfFreedom(), vertexReconstructorName_.value().c_str());
      // then add constituent tracks 
      xana::out(1) << "vertex found, loop on constituent tracks " << std::endl;
      vector <RefittedRecTrack> refTracks = (*iv)->refittedTracks();
      for  (vector<RefittedRecTrack>::const_iterator it=refTracks.begin();
       it!=refTracks.end(); it++) {
        xana::out(1) << "new refitted track " << std::endl;
        TrajectoryStateOnSurface tsi = it->innermostState() ;
        xana::out(1) << "(*track).innermostState() " <<(*track).innermostState() << std::endl;      
        TrajectoryStateOnSurface tso = it->outermostState() ;
        xana::out(1) << "(*track).outermostState() " <<(*track).outermostState() << std::endl;      
        // acessing track params if vertex exists is buggy but tsi, tso OK
        GlobalVector momentum = it->momentumAtVertex();
        xana::out(1) << "tsi.momentum " << tsi.globalMomentum() << std::endl;
        HepVector3D momentumAtVertex;
         momentumAtVertex=HepVector3D(momentum.x(),
         momentum.y(), momentum.z());
        xana::out(1) << "momentum at vertex " << momentumAtVertex << std::endl;        
      } 
      // this is needed otherwise access to originaltracks params is buggy
      vector<RecTrack> recTracks = (*iv)->originalTracks();
      for  (vector<RecTrack>::const_iterator it=recTracks.begin();
       it!=recTracks.end(); it++) {
        xana::out(1) << "new original track " << std::endl;
        TrajectoryStateOnSurface tsi = it->innermostState() ;
        xana::out(1) << "(*track).innermostState() " <<(*track).innermostState() << std::endl;      
        TrajectoryStateOnSurface tso = it->outermostState() ;
        xana::out(1) << "(*track).outermostState() " <<(*track).outermostState() << std::endl;      
        GlobalVector momentum = it->momentumAtVertex();
        HepVector3D momentumAtVertex;
         momentumAtVertex=HepVector3D(momentum.x(),
         momentum.y(), momentum.z());
        xana::out(1) << "momentum at vertex " << momentum << std::endl;
        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());
        xTrack = new XANATrack; 
        xana::out(1) << "track created " << std::endl;
        xTrack->setCharge(it->charge()); 
        xTrack->setChi2OverDof(it->normalisedChiSquared());
        xTrack->setNumberOfHits(it->foundHits());
        xTrack->setNumberOfLostHits(it->lostHits());
        xTrack->setImpactParameter(it->impactParameter3D().value());
        xTrack->setLongImpactParameter(it->zImpactParameter().value());
        xTrack->setTransImpactParameter(it->transverseImpactParameter().value());
        xTrack->setMomentumAtVertex(momentumAtVertex);
        xTrack->setPositionAtFirstPoint(firstHit);
        xTrack->setMomentumAtFirstPoint(momentumAtFirst);
        xTrack->setPositionAtLastPoint(lastHit);
        xTrack->setMomentumAtLastPoint(momentumAtLast);
        xTrack->setAlgoName(trackReconstructorName_.value());
        xana::out(1) << "adding track " << std::endl;
        xVertex->addTrack(xTrack);
      } 
    }
     
    // read XANAVertices 
    vector<XANAVertex>::iterator iv;
    for (iv=xVertices.begin(); iv!=xVertices.end(); iv++)  
     {
      xana::out(1)<<"[OVAL] new vertex with position "<<iv->getPosition()<<endl ;
     }
   }
      
} ;
 
static PKBuilder<VertexObserver> VertexObserver("VertexObserver") ;

---