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