XANATracks.cpp Source File

00001 
00002 #include <Utilities/Configuration/interface/Architecture.h>
00003 
00004 #include <XANADOO/XANATracks/interface/XANATrack.h>
00005 #include <XANADOO/XANATracks/interface/XANATrackHit.h>
00006 #include <XANADOO/XANAUtilities/interface/XANAToolkit.h>
00007 
00008 #include <CommonReco/PatternTools/interface/TTrack.h>
00009 #include <CommonReco/PatternTools/interface/TrackFinder.h>
00010 #include <TrackerReco/GtfPattern/interface/ModularKfReconstructor.h>
00011 #include <TrackerReco/TkSeedGenerator/interface/CombinatorialSeedGeneratorFromPixel.h>
00012 #include <TrackerReco/GtfPattern/interface/CombinatorialTrajectoryBuilder.h>
00013 #include <CommonReco/TrackFitters/interface/KFFittingSmoother.h>
00014 #include <CommonReco/PatternTools/interface/TrajectoryCleanerBySharedHits.h>
00015 #include <CommonReco/PatternTools/interface/RefittedRecTrack.h>
00016 #include "CommonDet/PatternPrimitives/interface/FreeTrajectoryState.h"
00017 #include "CommonReco/PatternTools/interface/ImpactPointExtrapolator.h"
00018 #include <CommonReco/MaterialEffects/interface/PropagatorWithMaterial.h>
00019 #include <CommonReco/MaterialEffects/interface/MultipleScatteringUpdator.h>
00020 #include <CommonReco/GeomPropagators/interface/AnalyticalPropagator.h>
00021 #include <CommonDet/PatternPrimitives/interface/FTSFromVertexToPointFactory.h>
00022 #include "CommonDet/BasicDet/interface/DetUnit.h"
00023 #include "CommonDet/BasicDet/interface/DetType.h"
00024 
00025 #include <CARF/G3Event/interface/G3EventProxy.h>
00026 #include <CARF/G3Event/interface/G3SetUp.h>
00027 #include <CARF/Reco/interface/RecCollection.h>
00028 #include <CARF/Reco/interface/RecQuery.h>
00029 
00030 #include <Utilities/Notification/interface/Observer.h>
00031 #include <Utilities/Notification/interface/PackageInitializer.h>
00032 #include <Utilities/UI/interface/SimpleConfigurable.h> 
00033  
00034 #include <CLHEP/Geometry/Point3D.h>
00035 #include <vector>
00036 
00037 struct TrackObserver : public Observer<G3EventProxy *> {
00038 
00039   int nbTracks_ ;
00040   double momentumSum_ ;
00041   SimpleConfigurable<bool> produceReference_;
00042   SimpleConfigurable<string> trackReconstructorName_;
00043   SimpleConfigurable<bool> writeTrackHits_;
00044  
00045   TrackObserver()
00046    : nbTracks_(0), momentumSum_(0)
00047    { 
00048     // initialize observer
00049     Observer<G3EventProxy *>::init() ;      
00050    
00051     // steering
00052     produceReference_=SimpleConfigurable<bool>(false,"XANATracks:ProduceReference");
00053     writeTrackHits_ = SimpleConfigurable<bool>(false,"XANATracks:WriteTrackHits");
00054     trackReconstructorName_=
00055      SimpleConfigurable<string>("CombinatorialTrackFinder","XANATracks:TrackReconstructor");
00056     // initialisation of tracking algo must be done in observing G3SetUp       
00057    }
00058    
00059   ~TrackObserver()
00060    {
00061     xana::out(8)<<"[OVAL] number of Tracks : "<<nbTracks_<<endl ;
00062     xana::out(8)<<"[OVAL] mean energy : "<<momentumSum_/nbTracks_<<endl ;
00063    }
00064     
00065   void upDate( G3EventProxy * ev )
00066    { 
00067     
00068     if (!ev) return ;    
00069     
00070     RecCollection<TTrack> tracks( RecQuery(trackReconstructorName_.value()) );
00071     RecCollection<TTrack>::const_iterator track;
00072     
00073     // to produce reference
00074     double momentum;
00075     if (produceReference_) {
00076        for (track=tracks.begin(); track!=tracks.end(); track++) {
00077         TrajectoryStateOnSurface ts = (*track).impactPointState() ;
00078         if (ts.isValid()) {
00079           momentum = ts.globalMomentum().mag() ;
00080           nbTracks_++ ;
00081           momentumSum_ += momentum ;
00082           xana::out(8)<<"[OVAL] new Track with momentum "<<momentum<<flush ;
00083           xana::out(8)<<" and impact parameter "<<(*track).impactParameter3D().value()<<endl ;
00084          }  
00085       }
00086       return ;
00087     } 
00088     
00089     // construct XANATracks from TTracks
00090     vector<XANATrack> xTracks;
00091     XANATrack *xTrack=0;
00092     for (track=tracks.begin(); track!=tracks.end(); track++) {  
00093       TrajectoryStateOnSurface ts = (*track).impactPointState() ;
00094       if (ts.isValid()) {
00095         momentum = ts.globalMomentum().mag() ;
00096         nbTracks_++ ;
00097         momentumSum_ += momentum ;
00098       }        
00099       TrajectoryStateOnSurface tsi = (*track).innermostState() ;
00100       TrajectoryStateOnSurface tso = (*track).outermostState() ;
00101       HepVector3D momentumAtVertex((*track).momentumAtVertex().x(), 
00102        (*track).momentumAtVertex().y(),(*track).momentumAtVertex().z());
00103       HepPoint3D firstHit(tsi.globalPosition().x(),tsi.globalPosition().y(),
00104         tsi.globalPosition().z()); 
00105       HepVector3D momentumAtFirst(tsi.globalMomentum().x(),tsi.globalMomentum().y(),
00106         tsi.globalMomentum().z());
00107       HepPoint3D lastHit(tso.globalPosition().x(),tso.globalPosition().y(),
00108         tso.globalPosition().z()); 
00109       HepVector3D momentumAtLast(tso.globalMomentum().x(),tso.globalMomentum().y(),
00110         tso.globalMomentum().z());
00111       xTrack = new XANATrack;
00112       xTrack->setCharge((*track).charge()); 
00113       xTrack->setChi2OverDof((*track).normalisedChiSquared());
00114       xTrack->setNumberOfHits((*track).foundHits());
00115       xTrack->setNumberOfLostHits((*track).lostHits());
00116       xTrack->setImpactParameter((*track).impactParameter3D().value());
00117       xTrack->setLongImpactParameter((*track).zImpactParameter().value());
00118       xTrack->setTransImpactParameter((*track).transverseImpactParameter().value());
00119       xTrack->setMomentumAtVertex(momentumAtVertex);
00120       xTrack->setPositionAtFirstPoint(firstHit);
00121       xTrack->setMomentumAtFirstPoint(momentumAtFirst);
00122       xTrack->setPositionAtLastPoint(lastHit);
00123       xTrack->setMomentumAtLastPoint(momentumAtLast);
00124       xTrack->setAlgoName(trackReconstructorName_.value());
00125       
00126       // hits (if requested)
00127       XANATrackHit *thit=0;
00128       if (writeTrackHits_) {
00129 
00130          std::vector<TrajectoryMeasurement> MyTrMeas = (*track).measurements();
00131          std::vector<TrajectoryMeasurement>::const_iterator IterTrMeas;
00132 
00133          xana::out(2) << " standard track " <<std::endl; 
00134          xana::out(2) << " number of measurements " << MyTrMeas.size()<<std::endl;
00135          xana::out(2) << " number of valid hits "<< (*track).foundHits()<<std::endl;
00136 
00137          for (IterTrMeas=MyTrMeas.begin(); IterTrMeas!=MyTrMeas.end(); IterTrMeas++)
00138            {
00139             // taking the rec hit 
00140             RecHit IterRecHit = IterTrMeas->recHit();                
00141             // checking if it is valid or not  
00142             if( !(IterRecHit.isValid()) ){xana::out(2) << "non valid hit" << endl;}
00143 
00144             // taking the updated state 
00145             TrajectoryStateOnSurface IterTSN = IterTrMeas->updatedState();                
00146             // checking if it is valid or not  
00147             if( !(IterTSN.isValid()) ){xana::out(2) << "non valid updated state" << endl;}
00148 
00149 
00150             // analysis only if both are valid 
00151             if (!(IterRecHit.isValid()) && (IterTSN.isValid())){xana::out(2)<<"non valid hit but valid updated state"<< endl;}
00152             if ((IterRecHit.isValid()) && !(IterTSN.isValid())){xana::out(2)<<"non valid updated state but valid hit"<< endl;}
00153             if(  IterRecHit.isValid() && IterTSN.isValid() )
00154             {
00155               HepPoint3D hitposition( (IterRecHit).globalPosition().x(), 
00156                                       (IterRecHit).globalPosition().y(), 
00157                                       (IterRecHit).globalPosition().z());
00158 
00159               HepPoint3D updposition( (IterTSN).globalPosition().x(),
00160                                       (IterTSN).globalPosition().y(),
00161                                       (IterTSN).globalPosition().z());
00162 
00163               HepVector3D updmomentum( (IterTSN).globalMomentum().x(),
00164                                        (IterTSN).globalMomentum().y(),
00165                                        (IterTSN).globalMomentum().z());
00166 
00167               thit = new XANATrackHit( hitposition, updposition, updmomentum ); 
00168               
00169               thit->setIsStereo((IterRecHit).det().detUnits().front()->type().isStereo());
00170 
00171 //                check
00172                   xana::out(1) << "NEW STEP - standard" << endl; 
00173                   xana::out(1) << "new hit" << endl;    
00174                   xana::out(1) << (IterRecHit).globalPosition().x() << endl;     
00175                   xana::out(1) << (IterRecHit).globalPosition().y() << endl;          
00176                   xana::out(1) << (IterRecHit).globalPosition().z() << endl;          
00177                   xana::out(1) << "new updated state: position" << endl;        
00178                   xana::out(1) << (IterTSN).globalPosition().x() << endl;     
00179                   xana::out(1) << (IterTSN).globalPosition().y() << endl;     
00180                   xana::out(1) << (IterTSN).globalPosition().z() << endl;     
00181                   xana::out(1) << "new updated state: momentum" << endl;        
00182                   xana::out(1) << (IterTSN).globalMomentum().x() << endl;     
00183                   xana::out(1) << (IterTSN).globalMomentum().y() << endl;     
00184                   xana::out(1) << (IterTSN).globalMomentum().z() << endl;   
00185 
00186               xTrack->addHit(thit);           
00187               delete thit;
00188 
00189            } // valid hit and upd state
00190            } // loop on measurements
00191       } // work on hits or not
00192       xTracks.push_back(*xTrack);
00193     }
00194      
00195     // read XANATracks 
00196     vector<XANATrack>::iterator it;
00197     for (it=xTracks.begin(); it!=xTracks.end(); it++)  
00198      {
00199       xana::out(8)<<"[OVAL] new Track with momentum "<<it->getMomentumAtVertex().mag()<<flush ;
00200       xana::out(8)<<" and impact parameter "<<it->getImpactParameter()<<endl ;
00201      }
00202    }
00203       
00204 } ;
00205  
00206 static PKBuilder<TrackObserver> trackObserver("TrackObserver") ;
00207 
00208