CMSSW Reference Manual

00001 #include "RecoTracker/SpecialSeedGenerators/interface/SeedGeneratorForCosmics.h"
00002 #include "FWCore/Framework/interface/EventSetup.h"
00003 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00004 #include "RecoTracker/TkHitPairs/interface/CosmicLayerPairs.h"
00005 #include "RecoPixelVertexing/PixelTriplets/interface/CosmicLayerTriplets.h"
00006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00007 #include "DataFormats/GeometryCommonDetAlgo/interface/GlobalError.h"
00008 #include "RecoTracker/TkSeedGenerator/interface/SeedFromConsecutiveHits.h"
00009 #include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
00010 #include "RecoTracker/TransientTrackingRecHit/interface/TkTransientTrackingRecHitBuilder.h"
00011 #include "TrackingTools/Records/interface/TransientRecHitRecord.h" 
00012 #include "RecoTracker/TkSeedGenerator/interface/FastHelix.h"
00013 void 
00014 SeedGeneratorForCosmics::init(const SiStripRecHit2DCollection &collstereo,
00015                               const SiStripRecHit2DCollection &collrphi ,
00016                               const SiStripMatchedRecHit2DCollection &collmatched,
00017                               const edm::EventSetup& iSetup)
00018 {
00019   iSetup.get<IdealMagneticFieldRecord>().get(magfield);
00020   iSetup.get<TrackerDigiGeometryRecord>().get(tracker);
00021   thePropagatorAl=    new PropagatorWithMaterial(alongMomentum,0.1057,&(*magfield) );
00022   thePropagatorOp=    new PropagatorWithMaterial(oppositeToMomentum,0.1057,&(*magfield) );
00023   theUpdator=       new KFUpdator();
00024   
00025   // get the transient builder
00026   //
00027 
00028   edm::ESHandle<TransientTrackingRecHitBuilder> theBuilder;
00029 
00030   iSetup.get<TransientRecHitRecord>().get(builderName,theBuilder);
00031   TTTRHBuilder = theBuilder.product();
00032   LogDebug("CosmicSeedFinder")<<" Hits built with  "<<hitsforseeds<<" hits";
00033  
00034     CosmicLayerPairs cosmiclayers(geometry);
00035 
00036     cosmiclayers.init(collstereo,collrphi,collmatched,iSetup);
00037     thePairGenerator=new CosmicHitPairGenerator(cosmiclayers,iSetup);
00038     HitPairs.clear();
00039     if ((hitsforseeds=="pairs")||(hitsforseeds=="pairsandtriplets")){
00040       thePairGenerator->hitPairs(region,HitPairs,iSetup);
00041   }
00042 
00043     CosmicLayerTriplets cosmiclayers2;
00044     cosmiclayers2.init(collstereo,collrphi,collmatched,geometry,iSetup);
00045     theTripletGenerator=new CosmicHitTripletGenerator(cosmiclayers2,iSetup);
00046     HitTriplets.clear();
00047     if ((hitsforseeds=="triplets")||(hitsforseeds=="pairsandtriplets")){
00048       theTripletGenerator->hitTriplets(region,HitTriplets,iSetup);
00049     }
00050 }
00051 
00052 SeedGeneratorForCosmics::SeedGeneratorForCosmics(edm::ParameterSet const& conf):
00053   conf_(conf),
00054   maxSeeds_(conf.getParameter<int32_t>("maxSeeds"))
00055 {  
00056 
00057   float ptmin=conf_.getParameter<double>("ptMin");
00058   float originradius=conf_.getParameter<double>("originRadius");
00059   float halflength=conf_.getParameter<double>("originHalfLength");
00060   float originz=conf_.getParameter<double>("originZPosition");
00061   seedpt = conf_.getParameter<double>("SeedPt");
00062 
00063   builderName = conf_.getParameter<std::string>("TTRHBuilder");   
00064   geometry=conf_.getUntrackedParameter<std::string>("GeometricStructure","STANDARD");
00065   region=GlobalTrackingRegion(ptmin,originradius,
00066                               halflength,originz);
00067   hitsforseeds=conf_.getUntrackedParameter<std::string>("HitsForSeeds","pairs");
00068   edm::LogInfo("SeedGeneratorForCosmics")<<" PtMin of track is "<<ptmin<< 
00069     " The Radius of the cylinder for seeds is "<<originradius <<"cm"  << " The set Seed Momentum" <<  seedpt;
00070 
00071   //***top-bottom
00072   positiveYOnly=conf_.getParameter<bool>("PositiveYOnly");
00073   negativeYOnly=conf_.getParameter<bool>("NegativeYOnly");
00074   //***
00075 
00076 
00077 }
00078 
00079 void SeedGeneratorForCosmics::run(TrajectorySeedCollection &output,const edm::EventSetup& iSetup){
00080   seeds(output,iSetup,region);
00081   delete thePairGenerator;
00082   delete theTripletGenerator;
00083   delete thePropagatorAl;
00084   delete thePropagatorOp;
00085   delete theUpdator;
00086 }
00087 bool SeedGeneratorForCosmics::seeds(TrajectorySeedCollection &output,
00088                                     const edm::EventSetup& iSetup,
00089                                     const TrackingRegion& region){
00090   LogDebug("CosmicSeedFinder")<<"Number of triplets "<<HitTriplets.size();
00091   LogDebug("CosmicSeedFinder")<<"Number of pairs "<<HitPairs.size();
00092 
00093   for (uint it=0;it<HitTriplets.size();it++){
00094     
00095     //const TrackingRecHit *hit = &(
00096     //    const TrackingRecHit* hit = it->hits();
00097 
00098   //   GlobalPoint inner = tracker->idToDet(HitTriplets[it].inner().RecHit()->
00099 //                                       geographicalId())->surface().
00100 //       toGlobal(HitTriplets[it].inner().RecHit()->localPosition());
00101     //const TrackingRecHit  *innerhit =  &(*HitTriplets[it].inner());
00102     //const TrackingRecHit  *middlehit =  &(*HitTriplets[it].middle());
00103 
00104     GlobalPoint inner = tracker->idToDet((*(HitTriplets[it].inner())).geographicalId())->surface().
00105       toGlobal((*(HitTriplets[it].inner())).localPosition());
00106 
00107     GlobalPoint middle = tracker->idToDet((*(HitTriplets[it].middle())).geographicalId())->surface().
00108       toGlobal((*(HitTriplets[it].middle())).localPosition());
00109 
00110     GlobalPoint outer = tracker->idToDet((*(HitTriplets[it].outer())).geographicalId())->surface().
00111       toGlobal((*(HitTriplets[it].outer())).localPosition());   
00112 
00113     // TransientTrackingRecHit::ConstRecHitPointer outrhit=TTTRHBuilder->build(HitPairs[is].outer())
00114 
00115     TransientTrackingRecHit::ConstRecHitPointer outrhit= TTTRHBuilder->build(HitTriplets[it].outer());
00116     //***top-bottom
00117     TransientTrackingRecHit::ConstRecHitPointer innrhit = TTTRHBuilder->build(HitTriplets[it].inner());
00118     if (positiveYOnly && (outrhit->globalPosition().y()<0 || innrhit->globalPosition().y()<0
00119                           || outrhit->globalPosition().y() < innrhit->globalPosition().y()
00120                           ) ) continue;
00121     if (negativeYOnly && (outrhit->globalPosition().y()>0 || innrhit->globalPosition().y()>0
00122                           || outrhit->globalPosition().y() > innrhit->globalPosition().y()
00123                           ) ) continue;
00124     //***
00125 
00126     edm::OwnVector<TrackingRecHit> hits;
00127     hits.push_back((*(HitTriplets[it].outer())).clone());
00128     FastHelix helix(inner, middle, outer,iSetup);
00129     GlobalVector gv=helix.stateAtVertex().parameters().momentum();
00130     float ch=helix.stateAtVertex().parameters().charge();
00131     float Mom = sqrt( gv.x()*gv.x() + gv.y()*gv.y() + gv.z()*gv.z() ); 
00132     if(Mom > 1000 || isnan(Mom))  continue;   // ChangedByDaniele 
00133 
00134     if (gv.y()>0){
00135       gv=-1.*gv;
00136       ch=-1.*ch;
00137     }
00138 
00139     GlobalTrajectoryParameters Gtp(outer,
00140                                    gv,int(ch), 
00141                                    &(*magfield));
00142     FreeTrajectoryState CosmicSeed(Gtp,
00143                                    CurvilinearTrajectoryError(AlgebraicSymMatrix55(AlgebraicMatrixID())));  
00144     if((outer.y()-inner.y())>0){
00145       const TSOS outerState =
00146         thePropagatorAl->propagate(CosmicSeed,
00147                                    tracker->idToDet((*(HitTriplets[it].outer())).geographicalId())->surface());
00148       if ( outerState.isValid()) {
00149         LogDebug("CosmicSeedFinder") <<"outerState "<<outerState;
00150         const TSOS outerUpdated= theUpdator->update( outerState,*outrhit);
00151         if ( outerUpdated.isValid()) {
00152           LogDebug("CosmicSeedFinder") <<"outerUpdated "<<outerUpdated;
00153           
00154           std::auto_ptr<PTrajectoryStateOnDet> PTraj(
00155             transformer.persistentState(outerUpdated,(*(HitTriplets[it].outer())).geographicalId().rawId())
00156           );
00157           output.push_back(TrajectorySeed(*PTraj,hits,alongMomentum));
00158 
00159           if ((maxSeeds_ > 0) && (output.size() > size_t(maxSeeds_))) {
00160             edm::LogWarning("SeedGeneratorForCosmics") << "Found too many seeds, bailing out.\n";
00161             output.clear(); 
00162             return false;
00163           }
00164         }
00165       }
00166     } else {
00167       const TSOS outerState =
00168         thePropagatorOp->propagate(CosmicSeed,
00169                                    tracker->idToDet((*(HitTriplets[it].outer())).geographicalId())->surface());
00170       if ( outerState.isValid()) {
00171         LogDebug("CosmicSeedFinder") <<"outerState "<<outerState;
00172         const TSOS outerUpdated= theUpdator->update( outerState,*outrhit);
00173         if ( outerUpdated.isValid()) {
00174           LogDebug("CosmicSeedFinder") <<"outerUpdated "<<outerUpdated;
00175           
00176           std::auto_ptr<PTrajectoryStateOnDet> PTraj(
00177                 transformer.persistentState(outerUpdated, (*(HitTriplets[it].outer())).geographicalId().rawId())
00178           );
00179           output.push_back(TrajectorySeed(*PTraj,hits,oppositeToMomentum));
00180 
00181           if ((maxSeeds_ > 0) && (output.size() > size_t(maxSeeds_))) {
00182             edm::LogWarning("SeedGeneratorForCosmics") << "Found too many seeds, bailing out.\n";
00183             output.clear(); 
00184             return false;
00185           }
00186         }
00187       }
00188     }
00189   }
00190   
00191 
00192   for(uint is=0;is<HitPairs.size();is++){
00193 
00194     
00195 
00196     GlobalPoint inner = tracker->idToDet((*(HitPairs[is].inner())).geographicalId())->surface().toGlobal((*(HitPairs[is].inner())).localPosition());
00197     GlobalPoint outer = tracker->idToDet((*(HitPairs[is].outer())).geographicalId())->surface().toGlobal((*(HitPairs[is].outer())).localPosition());
00198     
00199     LogDebug("CosmicSeedFinder") <<"inner point of the seed "<<inner <<" outer point of the seed "<<outer; 
00200     //RC const TransientTrackingRecHit* outrhit=TTTRHBuilder->build(HitPairs[is].outer().RecHit());  
00201     TransientTrackingRecHit::ConstRecHitPointer outrhit = TTTRHBuilder->build((HitPairs[is].outer()));
00202     //***top-bottom
00203     TransientTrackingRecHit::ConstRecHitPointer innrhit = TTTRHBuilder->build((HitPairs[is].inner()));
00204     if (positiveYOnly && (outrhit->globalPosition().y()<0 || innrhit->globalPosition().y()<0
00205                           || outrhit->globalPosition().y() < innrhit->globalPosition().y()
00206                           ) ) continue;
00207     if (negativeYOnly && (outrhit->globalPosition().y()>0 || innrhit->globalPosition().y()>0
00208                           || outrhit->globalPosition().y() > innrhit->globalPosition().y()
00209                           ) ) continue;
00210     //***
00211 
00212     edm::OwnVector<TrackingRecHit> hits;
00213     hits.push_back((*(HitPairs[is].outer())).clone());
00214     //    hits.push_back(HitPairs[is].inner()->clone());
00215 
00216     for (int i=0;i<2;i++){
00217       //FIRST STATE IS CALCULATED CONSIDERING THAT THE CHARGE CAN BE POSITIVE OR NEGATIVE
00218       int predsign=(2*i)-1;
00219       if((outer.y()-inner.y())>0){
00220         GlobalTrajectoryParameters Gtp(outer,
00221                                        (inner-outer)*(seedpt/(inner-outer).mag()),
00222                                        predsign, 
00223                                        &(*magfield));
00224         
00225         FreeTrajectoryState CosmicSeed(Gtp,
00226                                        CurvilinearTrajectoryError(AlgebraicSymMatrix55(AlgebraicMatrixID())));
00227         
00228         
00229         LogDebug("CosmicSeedFinder") << " FirstTSOS "<<CosmicSeed;
00230         //First propagation
00231         const TSOS outerState =
00232           thePropagatorAl->propagate(CosmicSeed,
00233                                      tracker->idToDet((*(HitPairs[is].outer())).geographicalId())->surface());
00234         if ( outerState.isValid()) {
00235           LogDebug("CosmicSeedFinder") <<"outerState "<<outerState;
00236           const TSOS outerUpdated= theUpdator->update( outerState,*outrhit);
00237           if ( outerUpdated.isValid()) {
00238             LogDebug("CosmicSeedFinder") <<"outerUpdated "<<outerUpdated;
00239             
00240             std::auto_ptr<PTrajectoryStateOnDet> PTraj(
00241               transformer.persistentState(outerUpdated, (*(HitPairs[is].outer())).geographicalId().rawId())
00242             );
00243             
00244             output.push_back( TrajectorySeed(*PTraj,hits,alongMomentum));
00245 
00246             if ((maxSeeds_ > 0) && (output.size() > size_t(maxSeeds_))) {
00247               edm::LogWarning("SeedGeneratorForCosmics") << "Found too many seeds, bailing out.\n";
00248               output.clear(); 
00249               return false;
00250             }
00251             
00252           }else      edm::LogWarning("CosmicSeedFinder") << " SeedForCosmics first update failed ";
00253         }else      edm::LogWarning("CosmicSeedFinder") << " SeedForCosmics first propagation failed ";
00254       
00255       
00256       }
00257       else{
00258         GlobalTrajectoryParameters Gtp(outer,
00259                                        (outer-inner)*(seedpt/(outer-inner).mag()),
00260                                        predsign, 
00261                                        &(*magfield));
00262         FreeTrajectoryState CosmicSeed(Gtp,
00263                                        CurvilinearTrajectoryError(AlgebraicSymMatrix55(AlgebraicMatrixID())));
00264         LogDebug("CosmicSeedFinder") << " FirstTSOS "<<CosmicSeed;
00265         //First propagation
00266         const TSOS outerState =
00267           thePropagatorAl->propagate(CosmicSeed,
00268                                      tracker->idToDet((*(HitPairs[is].outer())).geographicalId())->surface());
00269         if ( outerState.isValid()) {
00270           
00271           LogDebug("CosmicSeedFinder") <<"outerState "<<outerState;
00272           const TSOS outerUpdated= theUpdator->update( outerState,*outrhit);
00273           if ( outerUpdated.isValid()) {
00274           LogDebug("CosmicSeedFinder") <<"outerUpdated "<<outerUpdated;
00275 
00276           std::auto_ptr<PTrajectoryStateOnDet> PTraj(
00277             transformer.persistentState(outerUpdated,(*(HitPairs[is].outer())).geographicalId().rawId())
00278           );
00279           output.push_back(TrajectorySeed(*PTraj,hits,oppositeToMomentum));
00280 
00281           if ((maxSeeds_ > 0) && (output.size() > size_t(maxSeeds_))) {
00282             edm::LogWarning("SeedGeneratorForCosmics") << "Found too many seeds, bailing out.\n";
00283             output.clear(); 
00284             return false;
00285           }
00286 
00287           }else      edm::LogWarning("CosmicSeedFinder") << " SeedForCosmics first update failed ";
00288         }else      edm::LogWarning("CosmicSeedFinder") << " SeedForCosmics first propagation failed ";
00289       }
00290       
00291     }
00292   }
00293   return true;
00294 }
SeedGeneratorForCosmics.cc
