CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_0/src/RecoHI/HiMuonAlgos/src/HITrackVertexMaker.cc

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 //
00003 // Package:    TestMuL1L2
00004 // Class:      TestMuL1L2
00005 //
00006 // \class TestMuL1L2 TestMuL1L2.cc TestMuonL1L2/TestMuL1L2/src/TestMuL1L2.cc
00007 //
00008 // Original Author:  Dong Ho Moon
00009 //         Created:  Wed May  9 06:22:36 CEST 2007
00010 // $Id: HITrackVertexMaker.cc,v 1.7 2011/01/10 00:19:47 dlange Exp $
00011 //
00012 // added CAIR error cut
00013  
00014 #include "RecoHI/HiMuonAlgos/interface/HITrackVertexMaker.h" 
00015 
00016 // C++ Headers
00017 
00018 #include <memory>
00019 #include<iostream>
00020 #include<iomanip>
00021 #include<vector>
00022 #include<cmath>
00023 
00024 // ES include files 
00025 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
00026 #include "RecoTracker/Record/interface/CkfComponentsRecord.h"
00027 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h"
00028 
00029 // Muon trigger includes 
00030 #include "L1Trigger/GlobalMuonTrigger/interface/L1MuGlobalMuonTrigger.h"
00031 #include "DataFormats/L1GlobalMuonTrigger/interface/L1MuGMTReadoutRecord.h"
00032 #include "DataFormats/L1GlobalMuonTrigger/interface/L1MuGMTReadoutCollection.h"
00033 #include "DataFormats/L1GlobalMuonTrigger/interface/L1MuRegionalCand.h"
00034 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidate.h"
00035 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidateFwd.h"
00036 
00037 // Tracking includes
00038 #include "DataFormats/TrackReco/interface/Track.h"
00039 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
00040 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidate.h"
00041 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidateFwd.h"
00042 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
00043 #include "TrackingTools/MaterialEffects/interface/PropagatorWithMaterial.h"
00044 #include "TrackingTools/GeomPropagators/interface/StateOnTrackerBound.h"
00045 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
00046 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
00047 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
00048 #include "TrackingTools/TrajectoryFiltering/interface/RegionalTrajectoryFilter.h"
00049 #include "TrackingTools/PatternTools/interface/TrajectoryMeasurement.h"
00050 #include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
00051 #include "RecoTracker/TkNavigation/interface/NavigationSchoolFactory.h"
00052 #include "DataFormats/TrackReco/interface/TrackBase.h"
00053 #include "RecoVertex/VertexPrimitives/interface/TransientVertex.h"
00054 #include "DataFormats/VertexReco/interface/VertexFwd.h"
00055 #include "DataFormats/VertexReco/interface/Vertex.h"
00056 #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexFitter.h"
00057 //#include "TrackingTools/PatternTools/interface/TrajectoryStateClosestToBeamLineBuilder.h"
00058 #include "TrackingTools/PatternTools/interface/TSCBLBuilderNoMaterial.h"
00059 #include "RecoTracker/TkNavigation/interface/SimpleNavigationSchool.h"
00060 #include "TrackingTools/DetLayers/interface/NavigationSetter.h"
00061 #include "TrackingTools/DetLayers/interface/NavigationSchool.h"
00062 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
00063 #include "DataFormats/TrackReco/interface/TrackBase.h"
00064 #include "TrackingTools/PatternTools/interface/TrajectorySmoother.h"
00065 #include "TrackingTools/TrackFitters/interface/TrajectoryFitter.h"
00066 
00067 // Geometry includes
00068 
00069 #include "Geometry/CommonDetUnit/interface/GeomDetType.h"
00070 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
00071 #include "Geometry/CommonDetUnit/interface/GeomDetUnit.h"
00072 #include "Geometry/TrackerGeometryBuilder/interface/GluedGeomDet.h"
00073 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
00074 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
00075 #include "Geometry/TrackerNumberingBuilder/interface/GeometricDet.h"
00076 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
00077 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
00078 
00079 // RecoHIMuon includes
00080 
00081 #include "RecoHI/HiMuonAlgos/interface/HICFTSfromL1orL2.h"
00082 #include "RecoHI/HiMuonAlgos/interface/HICConst.h"
00083 #include "RecoHI/HiMuonAlgos/interface/FastMuPropagator.h"
00084 #include "RecoHI/HiMuonAlgos/interface/FmpConst.h"
00085 #include "RecoHI/HiMuonAlgos/interface/HICTkOuterStartingLayerFinder.h"
00086 #include "RecoHI/HiMuonAlgos/interface/DiMuonSeedGeneratorHIC.h"
00087 #include "RecoHI/HiMuonAlgos/interface/HICTrajectoryBuilder.h"
00088 #include "RecoHI/HiMuonAlgos/interface/HICSimpleNavigationSchool.h"
00089 #include "RecoHI/HiMuonAlgos/interface/HICMuonUpdator.h"
00090 #include "RecoHI/HiMuonAlgos/interface/HICConst.h"
00091 
00092 // Global  points
00093 
00094 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
00095 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
00096 
00097 // Vertex reco modification
00098 #include "TrackingTools/PatternTools/interface/TwoTrackMinimumDistance.h"
00099 #include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h"
00100 
00101 
00102 //Constructor
00103 
00104 using namespace reco;
00105 using namespace std;
00106 
00107 //#define DEBUG
00108 
00109 namespace cms{
00110 HITrackVertexMaker::HITrackVertexMaker(const edm::ParameterSet& ps1, const edm::EventSetup& es1)
00111 {
00112    
00113    L2candTag_          = ps1.getParameter< edm::InputTag > ("L2CandTag");
00114    rphirecHitsTag      = ps1.getParameter< edm::InputTag > ("rphiRecHits");
00115    builderName         = ps1.getParameter< std::string >   ("TTRHBuilder");
00116    primaryVertexTag    = ps1.getParameter< edm::InputTag > ("PrimaryVertexTag");
00117 #ifdef DEBUG
00118    std::cout<<" Start HI TrackVertexMaker constructor "<<std::endl;
00119 #endif
00120    pset_ = ps1;
00121    std::cout<<" No initialization "<<std::endl;
00122    eventCount = 0;
00123 #ifdef DEBUG
00124     std::cout<<" HICTrajectoryBuilder constructed "<<std::endl;
00125 #endif
00126 }
00127 
00128 
00129 
00130 //Destructor
00131 
00132 HITrackVertexMaker::~HITrackVertexMaker()
00133 {
00134 //  std::cout<<" HITrackVertexMaker::destructor "<<std::endl;
00135 } 
00136 
00137 bool HITrackVertexMaker::produceTracks(const edm::Event& e1, const edm::EventSetup& es1, HICConst* theHICConst, FmpConst* theFmpConst)
00138 {
00139    bool dimuon = false;
00140    
00141    edm::Handle<RecoChargedCandidateCollection> L2mucands;
00142 //   edm::Handle<TrackCollection> L2mucands;
00143    e1.getByLabel (L2candTag_,L2mucands);
00144       
00145 #ifdef DEBUG
00146    cout<<" Number of muon candidates "<<L2mucands->size()<<endl;
00147    if( L2mucands->size() < 2 ) cout<<"L2 muon failed"<<endl;
00148 #endif
00149 
00150    if( L2mucands->size() < 2 ) return dimuon;
00151    
00152 #ifdef DEBUG
00153    cout<<"L2 muon accepted"<<endl;
00154 #endif
00155 //   std::cout<<" Just do nothing for L3 but initiate ESHandles "<<std::endl;
00156    edm::Handle<reco::VertexCollection> vertexcands;
00157    e1.getByLabel (primaryVertexTag,vertexcands);
00158 
00159 #ifdef DEBUG   
00160    cout<<" Number of vertices primary  "<<vertexcands->size()<<endl;
00161    if(vertexcands->size()<1) cout<<" Primary vertex failed "<<endl;
00162 #endif
00163 
00164    if(vertexcands->size()<1) return dimuon;
00165 
00166 #ifdef DEBUG_COUNT
00167    cout<<" Accepted for L3 propagation  "<<endl;  
00168 #endif
00169    
00170 
00171    int iv = 0;
00172    for (reco::VertexCollection::const_iterator ipvertex=vertexcands->begin();ipvertex!=vertexcands->end();ipvertex++)
00173    {
00174 //     cout<<" Vertex position from pixels "<<(*ipvertex).position().z()<<endl;
00175      if (iv == 0) {theHICConst->setVertex((*ipvertex).position().z()); theFmpConst->setVertex((*ipvertex).position().z());} 
00176      iv++;
00177    } 
00178 
00179 //   cout << " Vertex is set to (found by pixel finder)"<<theHICConst->zvert<<endl;
00180    
00181    eventCount++;
00182 // ============================ Event accepted for L3   
00183 // Initialization from Records
00184 // 
00185    std::string updatorName = "KFUpdator";
00186    std::string propagatorAlongName    = "PropagatorWithMaterial";
00187    std::string propagatorOppositeName = "PropagatorWithMaterialOpposite";
00188    double theChiSquareCut = 500.;
00189    double nsig = 3.;
00190    double ptmin=1.;
00191 
00192   
00193    edm::ESHandle<MagneticField>                  magfield;
00194    edm::ESHandle<TransientTrackingRecHitBuilder> recHitBuilderHandle;
00195    edm::ESHandle<MeasurementTracker>             measurementTrackerHandle;
00196    edm::ESHandle<GeometricSearchTracker>         tracker; 
00197    edm::ESHandle<Propagator>                     propagatorAlongHandle;
00198    edm::ESHandle<Propagator>                     propagatorOppositeHandle;
00199    edm::ESHandle<TrajectoryStateUpdator>         updatorHandle;
00200 
00201    es1.get<TrackerRecoGeometryRecord>().get( tracker );
00202    es1.get<IdealMagneticFieldRecord>().get(magfield);
00203    es1.get<CkfComponentsRecord>().get("",measurementTrackerHandle);
00204    es1.get<TransientRecHitRecord>().get(builderName,recHitBuilderHandle); 
00205    es1.get<TrackingComponentsRecord>().get(propagatorAlongName,propagatorAlongHandle);
00206    es1.get<TrackingComponentsRecord>().get(propagatorOppositeName,propagatorOppositeHandle);
00207    es1.get<TrackingComponentsRecord>().get(updatorName,updatorHandle);
00208 
00209 // Initialization of navigation school
00210 
00211    if(eventCount == 1) {
00212 
00213    int lost=-1;
00214    int lostf=-1;
00215    theNavigationSchoolV.push_back(new HICSimpleNavigationSchool(&(*tracker), &(*magfield), lost, lostf)); 
00216    for(int i=11;i>-1;i--) {
00217      theNavigationSchoolV.push_back(new HICSimpleNavigationSchool(&(*tracker), &(*magfield), i, lostf));
00218    } 
00219    lost=-1;    
00220    for(int i=12;i>-1;i--) {
00221      theNavigationSchoolV.push_back(new HICSimpleNavigationSchool(&(*tracker), &(*magfield), lost, i));
00222    }
00223 
00224    } // initialization at the first event
00225        
00226     
00227 
00228    MinPtTrajectoryFilter theMinPtFilter(ptmin);
00229    HICMeasurementEstimator theEstimator(&(*tracker), &(*magfield), theChiSquareCut, nsig);
00230    
00231    
00232    HICTrajectoryBuilder theTrajectoryBuilder(        pset_, es1,
00233                                                      updatorHandle.product(),
00234                                                      propagatorAlongHandle.product(),
00235                                                      propagatorOppositeHandle.product(),
00236                                                      &theEstimator,
00237                                                      recHitBuilderHandle.product(),
00238                                                      measurementTrackerHandle.product(),
00239                                                      &theMinPtFilter);
00240 
00241    
00242 
00243 
00244  
00245   measurementTrackerHandle->update(e1);
00246   
00247 #ifdef DEBUG
00248   std::cout<<" After first tracker update "<<std::endl;
00249 #endif
00250 
00251 
00252 // For trajectory builder   
00253    int  theLowMult = 1;
00254    theEstimator.setHICConst(theHICConst);
00255    theEstimator.setMult(theLowMult);
00256 
00257 
00258    theTrajectoryBuilder.settracker(measurementTrackerHandle.product());
00259    
00260 //============================   
00261 //    FastMuPropagator* theFmp = new FastMuPropagator(&(*magfield),theFmpConst); 
00262 //    StateOnTrackerBound state(theFmp);  
00263       
00264     FastMuPropagator theFmp(&(*magfield),theFmpConst);
00265     StateOnTrackerBound state(&theFmp); 
00266 
00267     TrajectoryStateOnSurface tsos;
00268     
00269   
00270     HICFTSfromL1orL2 vFts(&(*magfield));
00271     
00272     
00273     int NumOfSigma=4;
00274     HICTkOuterStartingLayerFinder TkOSLF(NumOfSigma, &(*magfield), &(*tracker), theHICConst);
00275    
00276     int mult = 1;
00277     DiMuonSeedGeneratorHIC Seed(rphirecHitsTag,&(*magfield),&(*tracker), theHICConst, builderName, mult);
00278 
00279 //    vector<FreeTrajectoryState> theFts = vFts.createFTSfromStandAlone((*mucands));
00280 
00281    vector<FreeTrajectoryState> theFts = vFts.createFTSfromL2((*L2mucands)); 
00282 
00283 #ifdef DEBUG
00284     cout<<" Size of the freeTS "<<theFts.size()<<endl;
00285 #endif 
00286    
00287    
00288     
00289    DiMuonSeedGeneratorHIC::SeedContainer myseeds;  
00290    map<DetLayer*, DiMuonSeedGeneratorHIC::SeedContainer> seedmap;
00291    vector<Trajectory> theTmpTrajectories0;
00292  
00293 //   map< FreeTrajectoryState*, Trajectory>  
00294 
00295    vector<FreeTrajectoryState*> theFoundFts;
00296    map<FreeTrajectoryState*, vector<Trajectory> >  theMapFtsTraj;
00297 
00298 
00299    for(vector<FreeTrajectoryState>::iterator ifts=theFts.begin(); ifts!=theFts.end(); ifts++)
00300    {
00301       theTmpTrajectories0.clear();
00302 #ifdef DEBUG
00303      cout<<" cycle on Muon Trajectory State " <<(*ifts).parameters().position().perp()<<
00304                                           " " <<(*ifts).parameters().position().z()   <<endl;
00305 #endif
00306 
00307      tsos=state((*ifts));          
00308      if(tsos.isValid())
00309      {
00310 
00311 //        vector<Trajectory> theTmpTrajectories0;
00312 #ifdef DEBUG
00313         cout<<" Position "<<tsos.globalPosition().perp()<<" "<<tsos.globalPosition().phi()<<
00314               " "<<tsos.globalPosition().z()<<" "<<tsos.globalMomentum().perp()<<endl;
00315 #endif
00316 
00317 // Start to find starting layers
00318         FreeTrajectoryState* ftsnew=tsos.freeTrajectoryState();
00319         vector<DetLayer*> seedlayers = TkOSLF.startingLayers((*ftsnew));
00320         
00321 #ifdef DEBUG
00322         std::cout<<" The size of the starting layers "<<seedlayers.size()<<std::endl;
00323 #endif
00324 
00325         if( seedlayers.size() == 0 ) {
00326 #ifdef DEBUG
00327           cout<<" Starting layers failed for muon candidate "<<endl;
00328 #endif
00329           continue;
00330         }
00331         map<DetLayer*, DiMuonSeedGeneratorHIC::SeedContainer> seedmap = Seed.produce(e1 ,es1, 
00332                                                           (*ftsnew), tsos, (*ifts), 
00333                                                                 recHitBuilderHandle.product(),
00334                                                             measurementTrackerHandle.product(), 
00335                                                             &seedlayers);
00336 
00337  
00338 
00339        for( vector<DetLayer*>::iterator it = seedlayers.begin(); it != seedlayers.end(); it++)
00340        {
00341 
00342        DiMuonSeedGeneratorHIC::SeedContainer seeds = (*seedmap.find(*it)).second;
00343 
00344 #ifdef DEBUG
00345        std::cout<<" Layer::Position z "<<(**it).surface().position().z()<<
00346                   " Number of seeds in layer "<<seeds.size()<<std::endl;            
00347 #endif
00348 
00349   if(seeds.size() == 0) {
00350 #ifdef DEBUG
00351     std::cout<<" No seeds are found: do not continue with this Detlayer "<<std::endl;   
00352 #endif
00353     continue;
00354   }
00355     // set the first navigation (all layers without gap)
00356        
00357        NavigationSetter setter( *(theNavigationSchoolV[0]));
00358     
00359        for(DiMuonSeedGeneratorHIC::SeedContainer::iterator iseed=seeds.begin();
00360                                                            iseed!=seeds.end();iseed++)
00361        {
00362          std::vector<TrajectoryMeasurement> theV = (*iseed).measurements();
00363 
00364 #ifdef DEBUG
00365          std::cout<< "RecHIT Layer position r "<<
00366          theV[0].recHit()->globalPosition().perp()<<" phi "<<
00367          theV[0].recHit()->globalPosition().phi()<<" z "<<
00368          theV[0].recHit()->globalPosition().z()<<" momentum "<<
00369          theV[0].updatedState().freeTrajectoryState()->parameters().momentum().perp()<<" "<<
00370          theV[0].updatedState().freeTrajectoryState()->parameters().momentum().z()<<std::endl;
00371 #endif
00372 
00373          vector<Trajectory> theTmpTrajectories = theTrajectoryBuilder.trajectories(*iseed); 
00374     
00375 #ifdef DEBUG
00376         cout<<" Number of found trajectories "<<theTmpTrajectories.size()<<endl;         
00377 #endif
00378         if(theTmpTrajectories.size()>0) {
00379 
00380         theTmpTrajectories0.insert(theTmpTrajectories0.end(),
00381                                    theTmpTrajectories.begin(),
00382                                    theTmpTrajectories.end());
00383 #ifdef DEBUG
00384         std::cout<<"We found trajectories for at least one seed "<<theTmpTrajectories0.size()<<std::endl; 
00385 #endif
00386         break;
00387         }  // There are trajectories
00388        } // seeds           
00389 
00390     if( theTmpTrajectories0.size() > 0 ) {
00391 #ifdef DEBUG
00392         std::cout<<"We found trajectories for at least one seed "<<theTmpTrajectories0.size()<<"break layer cycle "<<std::endl; 
00393 #endif
00394      break;  
00395     }
00396    } // seedlayer
00397 
00398     if( theTmpTrajectories0.size() > 0 ) {
00399 #ifdef DEBUG
00400         std::cout<<"We found trajectories for at least one seed "<<theTmpTrajectories0.size()
00401                  <<"continue"<<std::endl; 
00402 #endif
00403      theFoundFts.push_back(&(*ifts));
00404      theMapFtsTraj[&(*ifts)] = theTmpTrajectories0;
00405      continue;
00406     } 
00407 
00408 #ifdef DEBUG
00409         std::cout<<"No trajectories for this FTS "<<theTmpTrajectories0.size()<<
00410                    "try second path"<<std::endl; 
00411 #endif
00412 
00413    
00414 // No trajectories found for this Muon FTS although seeds exist. 
00415 // Try to allow another trajectory map with lost 1 layer (only for barrel layers)
00416 
00417        for( vector<DetLayer*>::iterator it = seedlayers.begin(); it != seedlayers.end(); it++)
00418        {
00419 
00420           DiMuonSeedGeneratorHIC::SeedContainer seeds = (*seedmap.find(*it)).second;
00421 
00422           if(seeds.size() == 0) {
00423 #ifdef DEBUG
00424     std::cout<<" Second path: No seeds are found: do not continue with this Detlayer "<<std::endl;   
00425 #endif
00426              continue;
00427            }
00428 
00429        if((**it).location() == GeomDetEnumerators::endcap) {
00430        for(unsigned int j=13; j<theNavigationSchoolV.size();j++){
00431        NavigationSetter setter( *(theNavigationSchoolV[j]));
00432        for(DiMuonSeedGeneratorHIC::SeedContainer::iterator iseed=seeds.begin();
00433                                                            iseed!=seeds.end();iseed++)
00434        {
00435          std::vector<TrajectoryMeasurement> theV = (*iseed).measurements();
00436          vector<Trajectory> theTmpTrajectories = theTrajectoryBuilder.trajectories(*iseed);
00437 
00438          if(theTmpTrajectories.size()>0) {
00439 
00440          theTmpTrajectories0.insert(theTmpTrajectories0.end(),
00441                                     theTmpTrajectories.begin(),
00442                                     theTmpTrajectories.end());
00443 #ifdef DEBUG
00444          std::cout<<"Second path: We found trajectories for at least one seed in barrel layer "<<
00445                                     theTmpTrajectories0.size()<<std::endl; 
00446 #endif
00447           break;
00448          }  // There are trajectories
00449         } // seeds
00450 
00451          if( theTmpTrajectories0.size() > 0 ) {
00452 #ifdef DEBUG
00453          std::cout<<"Second path: no trajectories: we try next barrel layer "<<
00454                                   theTmpTrajectories0.size()<<std::endl; 
00455 #endif
00456                     break;
00457          }   
00458        } // navigation maps
00459        } else {
00460        for(int j=1; j<13;j++){
00461        NavigationSetter setter(*(theNavigationSchoolV[j]));
00462        for(DiMuonSeedGeneratorHIC::SeedContainer::iterator iseed=seeds.begin();
00463                                                            iseed!=seeds.end();iseed++)
00464        {
00465          std::vector<TrajectoryMeasurement> theV = (*iseed).measurements();
00466          vector<Trajectory> theTmpTrajectories = theTrajectoryBuilder.trajectories(*iseed);
00467 
00468          if(theTmpTrajectories.size()>0) {
00469 
00470          theTmpTrajectories0.insert(theTmpTrajectories0.end(),
00471                                     theTmpTrajectories.begin(),
00472                                     theTmpTrajectories.end());
00473 #ifdef DEBUG
00474          std::cout<<"Second path: We found trajectories for at least one seed in barrel layer "<<
00475                                   theTmpTrajectories0.size()<<std::endl; 
00476 #endif
00477           break;
00478          }  // There are trajectories
00479         } // seeds
00480 
00481          if( theTmpTrajectories0.size() > 0 ) {
00482 #ifdef DEBUG
00483          std::cout<<
00484    "Second path:  We found trajectories for at least one seed in barrel/endcap layer"<<
00485                                   theTmpTrajectories0.size()<<std::endl; 
00486 #endif
00487                     break;
00488          }   
00489      } // navigation maps
00490     } // barrel or endcap seed
00491          if( theTmpTrajectories0.size() > 0 ) {
00492 #ifdef DEBUG
00493          std::cout<<
00494               "Second path: We found trajectories for at least one seed in barrel/endcap layer "
00495                                <<
00496                                   theTmpTrajectories0.size()<<std::endl; 
00497 #endif     
00498               theFoundFts.push_back(&(*ifts));     
00499               theMapFtsTraj[&(*ifts)] = theTmpTrajectories0;
00500               break;
00501          }  
00502     } // seedlayer
00503    } // tsos. isvalid
00504  } // Muon Free trajectory state
00505         
00506 //
00507 // start fitting procedure
00508 //
00509 #ifdef DEBUG
00510     if(theFoundFts.size()>0 ) {
00511      std::cout<<" Event reconstruction finished with "<<theFoundFts.size()  
00512               <<std::endl;} 
00513            else {std::cout<<" Event reconstruction::no tracks found"<<std::endl;}
00514 #endif
00515     if(theFoundFts.size()<2)  return dimuon;
00516 
00517 // Look for vertex constraints
00518 
00519     edm::ESHandle<GlobalTrackingGeometry> globTkGeomHandle;
00520     es1.get<GlobalTrackingGeometryRecord>().get(globTkGeomHandle);
00521 
00522     reco::BeamSpot::CovarianceMatrix matrix;
00523     matrix(2,2) = 0.001;
00524     matrix(3,3) = 0.001;
00525 
00526     reco::BeamSpot bs( reco::BeamSpot::Point(0., 0., theHICConst->zvert),
00527                                                      0.1,
00528                                                      0.,
00529                                                      0.,
00530                                                      0.,
00531                                                     matrix
00532                      );
00533 
00534 
00535     reco::TrackBase::TrackAlgorithm Algo = reco::TrackBase::undefAlgorithm;
00536 
00537 // For trajectory refitting
00538         vector<reco::Track> firstTrack;
00539         vector<reco::TransientTrack> firstTransTracks;
00540         vector<reco::TrackRef> firstTrackRefs;
00541         vector<reco::Track> secondTrack;
00542         vector<reco::TransientTrack> secondTransTracks;
00543         vector<reco::TrackRef> secondTrackRefs;
00544         
00545     for(vector<FreeTrajectoryState*>::iterator it = theFoundFts.begin(); it!= theFoundFts.end(); it++)
00546     {
00547         vector<Trajectory> first = (*theMapFtsTraj.find(*it)).second;
00548 
00549         for(vector<Trajectory>::iterator im=first.begin();im!=first.end(); im++) {
00550 
00551           TrajectoryStateOnSurface innertsos;
00552           if (im->direction() == alongMomentum) {
00553            innertsos = im->firstMeasurement().updatedState();
00554           } else {
00555             innertsos = im->lastMeasurement().updatedState();
00556           }
00557 
00558 
00559          // CMSSW31X
00560           TSCBLBuilderNoMaterial tscblBuilder;
00561          // CMSSW22X
00562          //TrajectoryStateClosestToBeamLineBuilder tscblBuilder;
00563          TrajectoryStateClosestToBeamLine tscbl = tscblBuilder(*(innertsos.freeState()),bs);
00564 
00565          if (tscbl.isValid()==false) {
00566             //cout<<" false track "<<endl; 
00567          continue;
00568          }
00569 
00570          GlobalPoint v = tscbl.trackStateAtPCA().position();
00571          math::XYZPoint  pos( v.x(), v.y(), v.z() );
00572     
00573 //    cout<<" Position of track close to vertex "<<v.perp()<<" "<<v.z()<<" Primary vertex "<<theHICConst->zvert<<
00574 //                                     " charge "<<tscbl.trackStateAtPCA().charge()<<endl;
00575    
00576          if(v.perp() > 0.1 ) continue;
00577          if(fabs(v.z() - theHICConst->zvert ) > 0.4 ) continue;
00578 
00579           GlobalVector p = tscbl.trackStateAtPCA().momentum();
00580           math::XYZVector mom( p.x(), p.y(), p.z() );
00581 
00582           Track theTrack(im->chiSquared(),
00583                          im->ndof(), 
00584                          pos, mom, tscbl.trackStateAtPCA().charge(), 
00585                          tscbl.trackStateAtPCA().curvilinearError(),Algo);
00586            TransientTrack tmpTk( theTrack, &(*magfield), globTkGeomHandle );
00587 
00588          
00589            firstTrack.push_back( theTrack );
00590            firstTransTracks.push_back( tmpTk );
00591         }
00592 
00593     if( firstTrack.size() == 0 ) continue;
00594 
00595     for(vector<FreeTrajectoryState*>::iterator jt = it+1; jt!= theFoundFts.end(); jt++)
00596     {
00597         vector<Trajectory> second = (*theMapFtsTraj.find(*jt)).second;
00598        // cout<<" Number of trajectories first "<<first.size()<<" second "<<second.size()<<endl;
00599 
00600         for(vector<Trajectory>::iterator im=second.begin();im!=second.end(); im++) {
00601 
00602           TrajectoryStateOnSurface innertsos;
00603 
00604           if (im->direction() == alongMomentum) {
00605            innertsos = im->firstMeasurement().updatedState();
00606           } else {
00607             innertsos = im->lastMeasurement().updatedState();
00608           }
00609 
00610          TSCBLBuilderNoMaterial  tscblBuilder;
00611          //TrajectoryStateClosestToBeamLineBuilder tscblBuilder;
00612          TrajectoryStateClosestToBeamLine tscbl = tscblBuilder(*(innertsos.freeState()),bs);
00613 
00614          if (tscbl.isValid()==false) {
00615             //cout<<" false track "<<endl; 
00616          continue;
00617          }
00618 
00619          GlobalPoint v = tscbl.trackStateAtPCA().position();
00620          math::XYZPoint  pos( v.x(), v.y(), v.z() );
00621     
00622 //    cout<<" Position of track close to vertex "<<v.perp()<<" "<<v.z()<<" Primary vertex "<<theHICConst->zvert<<
00623 //                                     " charge "<<tscbl.trackStateAtPCA().charge()<<endl;
00624    
00625          if(v.perp() > 0.1 ) continue;
00626          if(fabs(v.z() - theHICConst->zvert ) > 0.4 ) continue;
00627 
00628           GlobalVector p = tscbl.trackStateAtPCA().momentum();
00629           math::XYZVector mom( p.x(), p.y(), p.z() );
00630 
00631           Track theTrack(im->chiSquared(),
00632                          im->ndof(), 
00633                          pos, mom, tscbl.trackStateAtPCA().charge(), 
00634                              tscbl.trackStateAtPCA().curvilinearError(),Algo);
00635            TransientTrack tmpTk( theTrack, &(*magfield), globTkGeomHandle );
00636 
00637          
00638            secondTrack.push_back( theTrack );
00639            secondTransTracks.push_back( tmpTk );
00640         }
00641         if( secondTrack.size() == 0 ) continue;
00642 
00643     } // FTS
00644   } // FTS 
00645 
00646 
00647    if( firstTrack.size() < 1 || secondTrack.size() < 1 ){ 
00648 #ifdef DEBUG
00649      cout<<" No enough tracks to get vertex "<<endl; 
00650 #endif
00651      return dimuon; 
00652    }
00653 
00654 
00655    bool useRefTrax=true;
00656    KalmanVertexFitter theFitter(useRefTrax);
00657    TransientVertex theRecoVertex;
00658 //
00659 // Create possible two particles vertices
00660 //
00661    vector<reco::TransientTrack> theTwoTransTracks;
00662 
00663    vector<TransientVertex> theVertexContainer;
00664 
00665    for(vector<reco::TransientTrack>::iterator iplus = firstTransTracks.begin(); 
00666                                               iplus != firstTransTracks.end(); iplus++)
00667    {
00668        for(vector<reco::TransientTrack>::iterator iminus = secondTransTracks.begin(); 
00669                iminus != secondTransTracks.end(); iminus++)
00670        {
00671            // To chech CAIR error before the vertex fitting
00672            TwoTrackMinimumDistance ttmd;
00673            bool CAIR_ST = false;
00674            GlobalTrajectoryParameters sta = (*iminus).impactPointState().globalParameters();
00675            GlobalTrajectoryParameters stb = (*iplus).impactPointState().globalParameters();
00676            ClosestApproachInRPhi theIniAlgo;
00677            
00678            CAIR_ST = theIniAlgo.calculate( sta, stb );
00679            //cout<<"%%%%% CAIR_ST : "<<CAIR_ST<<" %%%%%"<<endl;
00680            if(CAIR_ST == 0) continue;
00681            
00682            theTwoTransTracks.clear();
00683            theTwoTransTracks.push_back(*iplus);
00684            theTwoTransTracks.push_back(*iminus);
00685            theRecoVertex = theFitter.vertex(theTwoTransTracks);
00686            if( !theRecoVertex.isValid() ) {
00687                continue;
00688            } 
00689 
00690      //   cout<<" Vertex is found "<<endl;
00691      //   cout<<" Chi2 = "<<theRecoVertex.totalChiSquared()<<
00692      //           " r= "<<theRecoVertex.position().perp()<<
00693      //           " z= "<<theRecoVertex.position().z()<<endl;
00694 
00695 // Additional cuts       
00696            if ( theRecoVertex.totalChiSquared() > 0.0002 ) {
00697            //    cout<<" Vertex is failed with Chi2 : "<<theRecoVertex.totalChiSquared()<<endl; 
00698                continue;
00699            }
00700            if ( theRecoVertex.position().perp() > 0.08 ) {
00701            //    cout<<" Vertex is failed with r position : "<<theRecoVertex.position().perp()<<endl; 
00702                continue;
00703            }    
00704            if ( fabs(theRecoVertex.position().z()-theHICConst->zvert) > 0.06 ) {
00705            //    cout<<" Vertex is failed with z position : "<<theRecoVertex.position().z()<<endl; 
00706                continue;
00707            }
00708            double quality = theRecoVertex.normalisedChiSquared();
00709            std::vector<reco::TransientTrack> tracks = theRecoVertex.originalTracks();
00710 
00711            for (std::vector<reco::TransientTrack>::iterator ivb = tracks.begin(); ivb != tracks.end(); ivb++)
00712            {
00713                quality = quality * (*ivb).chi2() /(*ivb).ndof();
00714            }
00715            if( quality > 70. ) {
00716            //    cout<<" Vertex failed quality cut "<<quality<<endl; 
00717                continue;
00718            }
00719            theVertexContainer.push_back(theRecoVertex);
00720 
00721            dimuon = true;
00722            break;
00723        } // iminus
00724       if(dimuon) break; 
00725   } // iplus
00726     return dimuon;
00727 
00728 } 
00729 }
00730 
00731