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 std::auto_ptr<PTrajectoryStateOnDet> PTraj(
00149 transformer.persistentState(outerUpdated,(*(HitTriplets[it].outer())).geographicalId().rawId())
00150 );
00151 output.push_back(TrajectorySeed(*PTraj,hits,alongMomentum));
00152
00153 if ((maxSeeds_ > 0) && (output.size() > size_t(maxSeeds_))) {
00154 edm::LogError("TooManySeeds") << "Found too many seeds, bailing out.\n";
00155 output.clear();
00156 return false;
00157 }
00158 }
00159 }
00160 } else {
00161 const TSOS outerState =
00162 thePropagatorOp->propagate(CosmicSeed,
00163 tracker->idToDet((*(HitTriplets[it].outer())).geographicalId())->surface());
00164 if ( outerState.isValid()) {
00165 LogDebug("CosmicSeedFinder") <<"outerState "<<outerState;
00166 const TSOS outerUpdated= theUpdator->update( outerState,*outrhit);
00167 if ( outerUpdated.isValid()) {
00168 LogDebug("CosmicSeedFinder") <<"outerUpdated "<<outerUpdated;
00169
00170 std::auto_ptr<PTrajectoryStateOnDet> PTraj(
00171 transformer.persistentState(outerUpdated, (*(HitTriplets[it].outer())).geographicalId().rawId())
00172 );
00173 output.push_back(TrajectorySeed(*PTraj,hits,oppositeToMomentum));
00174
00175 if ((maxSeeds_ > 0) && (output.size() > size_t(maxSeeds_))) {
00176 edm::LogError("TooManySeeds") << "Found too many seeds, bailing out.\n";
00177 output.clear();
00178 return false;
00179 }
00180 }
00181 }
00182 }
00183 }
00184
00185
00186 for(unsigned int is=0;is<HitPairs.size();is++){
00187
00188
00189
00190 GlobalPoint inner = tracker->idToDet((*(HitPairs[is].inner())).geographicalId())->surface().toGlobal((*(HitPairs[is].inner())).localPosition());
00191 GlobalPoint outer = tracker->idToDet((*(HitPairs[is].outer())).geographicalId())->surface().toGlobal((*(HitPairs[is].outer())).localPosition());
00192
00193 LogDebug("CosmicSeedFinder") <<"inner point of the seed "<<inner <<" outer point of the seed "<<outer;
00194
00195 TransientTrackingRecHit::ConstRecHitPointer outrhit = TTTRHBuilder->build(HitPairs[is].outer()->hit());
00196
00197 TransientTrackingRecHit::ConstRecHitPointer innrhit = TTTRHBuilder->build(HitPairs[is].inner()->hit());
00198 if (positiveYOnly && (outrhit->globalPosition().y()<0 || innrhit->globalPosition().y()<0
00199 || outrhit->globalPosition().y() < innrhit->globalPosition().y()
00200 ) ) continue;
00201 if (negativeYOnly && (outrhit->globalPosition().y()>0 || innrhit->globalPosition().y()>0
00202 || outrhit->globalPosition().y() > innrhit->globalPosition().y()
00203 ) ) continue;
00204
00205
00206 edm::OwnVector<TrackingRecHit> hits;
00207 hits.push_back(HitPairs[is].outer()->hit()->clone());
00208
00209
00210 for (int i=0;i<2;i++){
00211
00212 int predsign=(2*i)-1;
00213 if((outer.y()-inner.y())>0){
00214 GlobalTrajectoryParameters Gtp(outer,
00215 (inner-outer)*(seedpt/(inner-outer).mag()),
00216 predsign,
00217 &(*magfield));
00218
00219 FreeTrajectoryState CosmicSeed(Gtp,
00220 CurvilinearTrajectoryError(AlgebraicSymMatrix55(AlgebraicMatrixID())));
00221
00222
00223 LogDebug("CosmicSeedFinder") << " FirstTSOS "<<CosmicSeed;
00224
00225 const TSOS outerState =
00226 thePropagatorAl->propagate(CosmicSeed,
00227 tracker->idToDet((*(HitPairs[is].outer())).geographicalId())->surface());
00228 if ( outerState.isValid()) {
00229 LogDebug("CosmicSeedFinder") <<"outerState "<<outerState;
00230 const TSOS outerUpdated= theUpdator->update( outerState,*outrhit);
00231 if ( outerUpdated.isValid()) {
00232 LogDebug("CosmicSeedFinder") <<"outerUpdated "<<outerUpdated;
00233
00234 std::auto_ptr<PTrajectoryStateOnDet> PTraj(
00235 transformer.persistentState(outerUpdated, (*(HitPairs[is].outer())).geographicalId().rawId())
00236 );
00237
00238 output.push_back( TrajectorySeed(*PTraj,hits,alongMomentum));
00239
00240 if ((maxSeeds_ > 0) && (output.size() > size_t(maxSeeds_))) {
00241 edm::LogError("TooManySeeds") << "Found too many seeds, bailing out.\n";
00242 output.clear();
00243 return false;
00244 }
00245
00246 }else edm::LogWarning("CosmicSeedFinder") << " SeedForCosmics first update failed ";
00247 }else edm::LogWarning("CosmicSeedFinder") << " SeedForCosmics first propagation failed ";
00248
00249
00250 }
00251 else{
00252 GlobalTrajectoryParameters Gtp(outer,
00253 (outer-inner)*(seedpt/(outer-inner).mag()),
00254 predsign,
00255 &(*magfield));
00256 FreeTrajectoryState CosmicSeed(Gtp,
00257 CurvilinearTrajectoryError(AlgebraicSymMatrix55(AlgebraicMatrixID())));
00258 LogDebug("CosmicSeedFinder") << " FirstTSOS "<<CosmicSeed;
00259
00260 const TSOS outerState =
00261 thePropagatorAl->propagate(CosmicSeed,
00262 tracker->idToDet((*(HitPairs[is].outer())).geographicalId())->surface());
00263 if ( outerState.isValid()) {
00264
00265 LogDebug("CosmicSeedFinder") <<"outerState "<<outerState;
00266 const TSOS outerUpdated= theUpdator->update( outerState,*outrhit);
00267 if ( outerUpdated.isValid()) {
00268 LogDebug("CosmicSeedFinder") <<"outerUpdated "<<outerUpdated;
00269
00270 std::auto_ptr<PTrajectoryStateOnDet> PTraj(
00271 transformer.persistentState(outerUpdated,(*(HitPairs[is].outer())).geographicalId().rawId())
00272 );
00273 output.push_back(TrajectorySeed(*PTraj,hits,oppositeToMomentum));
00274
00275 if ((maxSeeds_ > 0) && (output.size() > size_t(maxSeeds_))) {
00276 edm::LogError("TooManySeeds") << "Found too many seeds, bailing out.\n";
00277 output.clear();
00278 return false;
00279 }
00280
00281 }else edm::LogWarning("CosmicSeedFinder") << " SeedForCosmics first update failed ";
00282 }else edm::LogWarning("CosmicSeedFinder") << " SeedForCosmics first propagation failed ";
00283 }
00284
00285 }
00286 }
00287 return true;
00288 }