CMS 3D CMS Logo

/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/RecoTracker/SpecialSeedGenerators/src/CtfSpecialSeedGenerator.cc

Go to the documentation of this file.
00001 #include "RecoTracker/SpecialSeedGenerators/interface/CtfSpecialSeedGenerator.h"
00002 //#include "RecoTracker/SpecialSeedGenerators/interface/CosmicLayerTriplets.h"
00003 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
00004 #include "Geometry/Records/interface/IdealGeometryRecord.h"
00005 #include "DataFormats/GeometrySurface/interface/RectangularPlaneBounds.h"
00006 #include "TrackingTools/GeomPropagators/interface/StraightLinePlaneCrossing.h"
00007 
00008 #include "RecoTracker/TkTrackingRegions/interface/TrackingRegionProducerFactory.h"
00009 #include "RecoTracker/TkTrackingRegions/interface/TrackingRegion.h"
00010 #include "RecoTracker/TkTrackingRegions/interface/OrderedHitsGeneratorFactory.h"
00011 #include "RecoTracker/TkSeedingLayers/interface/OrderedSeedingHits.h"
00012 
00013 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00014 
00015 #include "RecoTracker/SpecialSeedGenerators/interface/ClusterChecker.h"
00016 
00017 using namespace ctfseeding;
00018 
00019 CtfSpecialSeedGenerator::CtfSpecialSeedGenerator(const edm::ParameterSet& conf): 
00020   conf_(conf),
00021   requireBOFF(conf.getParameter<bool>("requireBOFF")),
00022   theMaxSeeds(conf.getParameter<int32_t>("maxSeeds"))
00023 {
00024         useScintillatorsConstraint = conf_.getParameter<bool>("UseScintillatorsConstraint");
00025         edm::LogVerbatim("CtfSpecialSeedGenerator") << "Constructing CtfSpecialSeedGenerator";
00026         produces<TrajectorySeedCollection>();
00027         theSeedBuilder =0; 
00028         theRegionProducer =0;
00029 }
00030 
00031 CtfSpecialSeedGenerator::~CtfSpecialSeedGenerator(){
00032 }
00033 
00034 void CtfSpecialSeedGenerator::endRun(edm::Run const&, edm::EventSetup const&){
00035     if (theSeedBuilder)    { delete theSeedBuilder;    theSeedBuilder = 0; }
00036     if (theRegionProducer) { delete theRegionProducer; theRegionProducer = 0; }
00037     std::vector<OrderedHitsGenerator*>::iterator iGen;  
00038     for (iGen = theGenerators.begin(); iGen != theGenerators.end(); iGen++){
00039         delete (*iGen);
00040     }
00041     theGenerators.clear();
00042 }
00043 
00044 void CtfSpecialSeedGenerator::beginRun(edm::Run const&, const edm::EventSetup& iSetup){
00045         std::string builderName = conf_.getParameter<std::string>("TTRHBuilder");
00046         iSetup.get<TransientRecHitRecord>().get(builderName,theBuilder);
00047 
00048         iSetup.get<IdealMagneticFieldRecord>().get(theMagfield);
00049         iSetup.get<TrackerDigiGeometryRecord>().get(theTracker);
00050 
00051         edm::LogVerbatim("CtfSpecialSeedGenerator") << "Initializing...";
00052         if (useScintillatorsConstraint){
00053                 edm::ParameterSet upperScintPar = conf_.getParameter<edm::ParameterSet>("UpperScintillatorParameters");
00054                 edm::ParameterSet lowerScintPar = conf_.getParameter<edm::ParameterSet>("LowerScintillatorParameters");
00055                 RectangularPlaneBounds upperBounds(upperScintPar.getParameter<double>("WidthInX"),
00056                                              upperScintPar.getParameter<double>("LenghtInZ"),
00057                                              1);
00058                 GlobalPoint upperPosition(upperScintPar.getParameter<double>("GlobalX"),
00059                                         upperScintPar.getParameter<double>("GlobalY"),
00060                                         upperScintPar.getParameter<double>("GlobalZ"));
00061                 edm::LogVerbatim("CtfSpecialSeedGenerator")
00062                         << "Upper Scintillator position x, y, z " << upperPosition.x()
00063                         << ", " << upperPosition.y() << ", " << upperPosition.z();
00064                 RectangularPlaneBounds lowerBounds(lowerScintPar.getParameter<double>("WidthInX"),
00065                                              lowerScintPar.getParameter<double>("LenghtInZ"),
00066                                              1);
00067                 GlobalPoint lowerPosition(lowerScintPar.getParameter<double>("GlobalX"),
00068                                      lowerScintPar.getParameter<double>("GlobalY"),
00069                                      lowerScintPar.getParameter<double>("GlobalZ"));
00070                 edm::LogVerbatim("CtfSpecialSeedGenerator")
00071                         << "Lower Scintillator position x, y, z " << lowerPosition.x()
00072                         << ", " << lowerPosition.y() << ", " << lowerPosition.z() ;
00073                 TkRotation<float> rot(1,0,0,0,0,1,0,1,0);
00074                 upperScintillator = BoundPlane::build(upperPosition, rot, &upperBounds);
00075                 lowerScintillator = BoundPlane::build(lowerPosition, rot, &lowerBounds);
00076         } 
00077         edm::ParameterSet regfactoryPSet = conf_.getParameter<edm::ParameterSet>("RegionFactoryPSet");
00078         std::string regfactoryName = regfactoryPSet.getParameter<std::string>("ComponentName");
00079         theRegionProducer = TrackingRegionProducerFactory::get()->create(regfactoryName,regfactoryPSet);
00080         
00081         edm::ESHandle<Propagator>  propagatorAlongHandle;
00082         iSetup.get<TrackingComponentsRecord>().get("PropagatorWithMaterial",propagatorAlongHandle);
00083         edm::ESHandle<Propagator>  propagatorOppositeHandle;
00084         iSetup.get<TrackingComponentsRecord>().get("PropagatorWithMaterialOpposite",propagatorOppositeHandle);
00085 
00086 /*      edm::ParameterSet hitsfactoryOutInPSet = conf_.getParameter<edm::ParameterSet>("OrderedHitsFactoryOutInPSet");
00087         std::string hitsfactoryOutInName = hitsfactoryOutInPSet.getParameter<std::string>("ComponentName");
00088         hitsGeneratorOutIn = OrderedHitsGeneratorFactory::get()->create( hitsfactoryOutInName, hitsfactoryOutInPSet);
00089         std::string propagationDirection = hitsfactoryOutInPSet.getUntrackedParameter<std::string>("PropagationDirection", 
00090                                                                                                     "alongMomentum");
00091         if (propagationDirection == "alongMomentum") outInPropagationDirection = alongMomentum;
00092         else outInPropagationDirection = oppositeToMomentum;
00093         edm::LogVerbatim("CtfSpecialSeedGenerator") << "hitsGeneratorOutIn done";
00094 
00095         edm::ParameterSet hitsfactoryInOutPSet = conf_.getParameter<edm::ParameterSet>("OrderedHitsFactoryInOutPSet");
00096         std::string hitsfactoryInOutName = hitsfactoryInOutPSet.getParameter<std::string>("ComponentName");
00097         hitsGeneratorInOut = OrderedHitsGeneratorFactory::get()->create( hitsfactoryInOutName, hitsfactoryInOutPSet);
00098 
00099         propagationDirection = hitsfactoryInOutPSet.getUntrackedParameter<std::string>("PropagationDirection",
00100                                                                                         "alongMomentum");
00101         if (propagationDirection == "alongMomentum") inOutPropagationDirection = alongMomentum;
00102         else inOutPropagationDirection = oppositeToMomentum;
00103         edm::LogVerbatim("CtfSpecialSeedGenerator") << "hitsGeneratorInOut done";
00104         if (!hitsGeneratorOutIn || !hitsGeneratorInOut) 
00105                 throw cms::Exception("CtfSpecialSeedGenerator") << "Only corcrete implementation GenericPairOrTripletGenerator of OrderedHitsGenerator is allowed ";
00106 */
00107         std::vector<edm::ParameterSet> pSets = conf_.getParameter<std::vector<edm::ParameterSet> >("OrderedHitsFactoryPSets");
00108         std::vector<edm::ParameterSet>::const_iterator iPSet;
00109         for (iPSet = pSets.begin(); iPSet != pSets.end(); iPSet++){
00110                 std::string hitsfactoryName = iPSet->getParameter<std::string>("ComponentName");
00111                 theGenerators.push_back(OrderedHitsGeneratorFactory::get()->create( hitsfactoryName, *iPSet));
00112                 std::string propagationDirection = iPSet->getParameter<std::string>("PropagationDirection");
00113                 if (propagationDirection == "alongMomentum") thePropDirs.push_back(alongMomentum);
00114                 else thePropDirs.push_back(oppositeToMomentum);
00115                 std::string navigationDirection = iPSet->getParameter<std::string>("NavigationDirection");              
00116                 if (navigationDirection == "insideOut") theNavDirs.push_back(insideOut);
00117                 else theNavDirs.push_back(outsideIn);
00118                 edm::LogVerbatim("CtfSpecialSeedGenerator") << "hitsGenerator done";
00119         } 
00120         bool setMomentum = conf_.getParameter<bool>("SetMomentum");
00121         std::vector<int> charges;
00122         if (setMomentum){
00123                 charges = conf_.getParameter<std::vector<int> >("Charges");
00124         }
00125         theSeedBuilder = new SeedFromGenericPairOrTriplet(theMagfield.product(), 
00126                                                           theTracker.product(), 
00127                                                           theBuilder.product(),
00128                                                           propagatorAlongHandle.product(),
00129                                                           propagatorOppositeHandle.product(),
00130                                                           charges,
00131                                                           setMomentum,
00132                                                           conf_.getParameter<double>("ErrorRescaling"));
00133         double p = 1;
00134         if (setMomentum) {
00135                 p = conf_.getParameter<double>("SeedMomentum");
00136                 theSeedBuilder->setMomentumTo(p);
00137         }
00138 
00139 }
00140 
00141 void CtfSpecialSeedGenerator::produce(edm::Event& e, const edm::EventSetup& iSetup)
00142 {
00143   // get Inputs
00144   std::auto_ptr<TrajectorySeedCollection> output(new TrajectorySeedCollection);
00145   
00146   //check on the number of clusters
00147   ClusterChecker check(conf_);
00148   if ( !requireBOFF || (theMagfield->inTesla(GlobalPoint(0,0,0)).mag() == 0.00) ) {
00149       size_t clustsOrZero = check.tooManyClusters(e);
00150       if (!clustsOrZero){
00151           bool ok = run(iSetup, e, *output);
00152           if (!ok) { ; } // nothing to do
00153       } else edm::LogError("TooManyClusters") << "Found too many clusters (" << clustsOrZero << "), bailing out.\n";
00154   }
00155   
00156   
00157   edm::LogVerbatim("CtfSpecialSeedGenerator") << " number of seeds = "<< output->size();
00158   e.put(output);
00159 }
00160 
00161 bool CtfSpecialSeedGenerator::run(const edm::EventSetup& iSetup,
00162                                                const edm::Event& e,
00163                                                TrajectorySeedCollection& output){
00164         std::vector<TrackingRegion*> regions = theRegionProducer->regions(e, iSetup);
00165         std::vector<TrackingRegion*>::const_iterator iReg;
00166         bool ok = true;
00167         for (iReg = regions.begin(); iReg != regions.end(); iReg++){
00168                 if(!theSeedBuilder->momentumFromPSet()) theSeedBuilder->setMomentumTo((*iReg)->ptMin());
00169                 std::vector<OrderedHitsGenerator*>::const_iterator iGen;
00170                 int i = 0;
00171                 for (iGen = theGenerators.begin(); iGen != theGenerators.end(); iGen++){ 
00172                   ok = buildSeeds(iSetup, 
00173                              e, 
00174                              (*iGen)->run(**iReg, e, iSetup),
00175                              theNavDirs[i], 
00176                              thePropDirs[i], 
00177                              output);
00178                   i++;
00179                   if (!ok) break;
00180                 }
00181                 if (!ok) break;
00182         }
00183         //clear memory
00184         for (std::vector<TrackingRegion*>::iterator iReg = regions.begin(); iReg != regions.end(); iReg++){
00185                 delete *iReg;
00186         }
00187         return ok;
00188 }
00189 
00190 bool CtfSpecialSeedGenerator::buildSeeds(const edm::EventSetup& iSetup,
00191                                          const edm::Event& e,
00192                                          const OrderedSeedingHits& osh,
00193                                          const NavigationDirection& navdir,
00194                                          const PropagationDirection& dir,
00195                                          TrajectorySeedCollection& output){ 
00196   //SeedFromGenericPairOrTriplet seedBuilder(conf_, magfield.product(), tracker.product(), theBuilder.product());
00197   edm::LogInfo("CtfSpecialSeedGenerator")<<"osh.size() " << osh.size();
00198   for (unsigned int i = 0; i < osh.size(); i++){
00199         SeedingHitSet shs = osh[i];
00200         if (preliminaryCheck(shs,iSetup)){
00201                 std::vector<TrajectorySeed*> seeds = theSeedBuilder->seed(shs, 
00202                                                                         dir,
00203                                                                         navdir, 
00204                                                                         iSetup);
00205                 for (std::vector<TrajectorySeed*>::const_iterator iSeed = seeds.begin(); iSeed != seeds.end(); iSeed++){
00206                   if (!*iSeed) {edm::LogError("CtfSpecialSeedGenerator")<<"a seed pointer is null. skipping.";continue;}
00207                         if (postCheck(**iSeed)){
00208                                 output.push_back(**iSeed);
00209                         }
00210                         delete *iSeed;
00211                         edm::LogVerbatim("CtfSpecialSeedGenerator") << "Seed built";
00212                 }
00213         }
00214   }      
00215   if ((theMaxSeeds > 0) && (output.size() > size_t(theMaxSeeds))) {
00216     edm::LogWarning("TooManySeeds") << "Too many seeds ("<< output.size() <<"), bailing out.\n";
00217     output.clear(); 
00218     return false;
00219   }
00220   return true;
00221 }
00222 //checks the hits are on diffrent layers
00223 bool CtfSpecialSeedGenerator::preliminaryCheck(const SeedingHitSet& shs, const edm::EventSetup &es ){
00224 
00225         edm::ESHandle<TrackerTopology> tTopo;
00226         es.get<IdealGeometryRecord>().get(tTopo);
00227 
00228         std::vector<std::pair<unsigned int, unsigned int> > vSubdetLayer;
00229         //std::vector<std::string> vSeedLayerNames;
00230         bool checkHitsAtPositiveY       = conf_.getParameter<bool>("SeedsFromPositiveY");
00231         //***top-bottom
00232         bool checkHitsAtNegativeY       = conf_.getParameter<bool>("SeedsFromNegativeY");
00233         //***
00234         bool checkHitsOnDifferentLayers = conf_.getParameter<bool>("CheckHitsAreOnDifferentLayers");
00235       unsigned int nHits = shs.size();
00236       for (unsigned int iHit=0; iHit < nHits; ++iHit) {
00237                 //hits for the seeds must be at positive y
00238             const TrackingRecHit * trh = shs[iHit]->hit();
00239                 TransientTrackingRecHit::RecHitPointer recHit = theBuilder->build(trh);
00240                 GlobalPoint hitPos = recHit->globalPosition();
00241                 //GlobalPoint point = 
00242                 //  theTracker->idToDet(iHits->geographicalId() )->surface().toGlobal(iHits->localPosition());
00243                 if (checkHitsAtPositiveY){ if (hitPos.y() < 0) return false;}
00244                 //***top-bottom
00245                 if (checkHitsAtNegativeY){ if (hitPos.y() > 0) return false;}
00246                 //***
00247                 //std::string name = iHits->seedinglayer().name(); 
00248                 //hits for the seeds must be in different layers
00249                 unsigned int subid=(*trh).geographicalId().subdetId();
00250                 unsigned int layer = tTopo->layer( (*trh).geographicalId());
00251                 std::vector<std::pair<unsigned int, unsigned int> >::const_iterator iter;
00252                 //std::vector<std::string>::const_iterator iNames;
00253                 if (checkHitsOnDifferentLayers){
00254                         
00255                         for (iter = vSubdetLayer.begin(); iter != vSubdetLayer.end(); iter++){
00256                                 if (iter->first == subid && iter->second == layer) return false;
00257                         }
00258                         /*
00259                         for (iNames = vSeedLayerNames.begin(); iNames != vSeedLayerNames.end(); iNames++){
00260                                 if (*iNames == name) return false;
00261                         }
00262                         */
00263                 }
00264                 //vSeedLayerNames.push_back(iHits->seedinglayer().name());
00265                 vSubdetLayer.push_back(std::make_pair(subid, layer));   
00266         }
00267         return true;
00268 }
00269 
00270 
00271 bool CtfSpecialSeedGenerator::postCheck(const TrajectorySeed& seed){
00272         if (!useScintillatorsConstraint) return true; 
00273         
00274         PTrajectoryStateOnDet pstate = seed.startingState();
00275         TrajectoryStateOnSurface theTSOS = trajectoryStateTransform::transientState(pstate,
00276                                                                       &(theTracker->idToDet(DetId(pstate.detId()))->surface()),
00277                                                                       &(*theMagfield)); 
00278         FreeTrajectoryState* state = theTSOS.freeState();       
00279         StraightLinePlaneCrossing planeCrossingLower( Basic3DVector<float>(state->position()), 
00280                                                       Basic3DVector<float>(state->momentum()),
00281                                                       alongMomentum);
00282         StraightLinePlaneCrossing planeCrossingUpper( Basic3DVector<float>(state->position()),
00283                                                       Basic3DVector<float>(state->momentum()),
00284                                                       oppositeToMomentum);
00285         std::pair<bool,StraightLinePlaneCrossing::PositionType> positionUpper = 
00286           planeCrossingUpper.position(*upperScintillator);
00287         std::pair<bool,StraightLinePlaneCrossing::PositionType> positionLower = 
00288           planeCrossingLower.position(*lowerScintillator);
00289         if (!(positionUpper.first && positionLower.first)) {
00290                  edm::LogVerbatim("CtfSpecialSeedGenerator::checkDirection") 
00291                         << "Scintillator plane not crossed";
00292                  return false;
00293         }
00294         LocalPoint positionUpperLocal = upperScintillator->toLocal((GlobalPoint)(positionUpper.second));
00295         LocalPoint positionLowerLocal = lowerScintillator->toLocal((GlobalPoint)(positionLower.second));
00296         if (upperScintillator->bounds().inside(positionUpperLocal) && 
00297             lowerScintillator->bounds().inside(positionLowerLocal)) {
00298                  edm::LogVerbatim("CtfSpecialSeedGenerator::checkDirection") 
00299                                 << "position on Upper scintillator " 
00300                                 << positionUpper.second;
00301                  edm::LogVerbatim("CtfSpecialSeedGenerator::checkDirection") 
00302                                 << "position on Lower scintillator " 
00303                                 << positionLower.second;
00304                  
00305                  return true;
00306         }
00307         edm::LogVerbatim("CtfSpecialSeedGenerator::checkDirection") 
00308                 << "scintillator not crossed in bounds: position on Upper scintillator " 
00309                 << positionUpper.second << " position on Lower scintillator " << positionLower.second;
00310         return false;
00311 }
00312 
00313