CMS 3D CMS Logo

/data/doxygen/doxygen-1.7.3/gen/CMSSW_4_2_8/src/RecoTracker/SpecialSeedGenerators/src/SeedGeneratorForCosmics.cc

Go to the documentation of this file.
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   // get the transient builder
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   //***top-bottom
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     //const TrackingRecHit *hit = &(
00090     //    const TrackingRecHit* hit = it->hits();
00091 
00092   //   GlobalPoint inner = tracker->idToDet(HitTriplets[it].inner().RecHit()->
00093 //                                       geographicalId())->surface().
00094 //       toGlobal(HitTriplets[it].inner().RecHit()->localPosition());
00095     //const TrackingRecHit  *innerhit =  &(*HitTriplets[it].inner());
00096     //const TrackingRecHit  *middlehit =  &(*HitTriplets[it].middle());
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     // TransientTrackingRecHit::ConstRecHitPointer outrhit=TTTRHBuilder->build(HitPairs[is].outer())
00108 
00109     TransientTrackingRecHit::ConstRecHitPointer outrhit= TTTRHBuilder->build(HitTriplets[it].outer()->hit());
00110     //***top-bottom
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;   // ChangedByDaniele 
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     //RC const TransientTrackingRecHit* outrhit=TTTRHBuilder->build(HitPairs[is].outer().RecHit());  
00195     TransientTrackingRecHit::ConstRecHitPointer outrhit = TTTRHBuilder->build(HitPairs[is].outer()->hit());
00196     //***top-bottom
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     //    hits.push_back(HitPairs[is].inner()->clone());
00209 
00210     for (int i=0;i<2;i++){
00211       //FIRST STATE IS CALCULATED CONSIDERING THAT THE CHARGE CAN BE POSITIVE OR NEGATIVE
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         //First propagation
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         //First propagation
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 }