Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:28:45 

0001 //
0002 // Package:         RecoTracker/TkSeedGenerator
0003 // Class:           GlobalPixelLessSeedGenerator
0004 //
0005 // From CosmicSeedGenerator+SimpleCosmicBONSeeder, with changes by Giovanni
0006 // to seed Cosmics with B != 0
0007 
0008 #include "RecoTracker/SpecialSeedGenerators/interface/SimpleCosmicBONSeeder.h"
0009 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
0010 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h"
0011 #include "FWCore/Utilities/interface/isFinite.h"
0012 #include "FWCore/Framework/interface/ConsumesCollector.h"
0013 #include "TrackingTools/TransientTrackingRecHit/interface/SeedingLayerSetsHits.h"
0014 typedef SeedingHitSet::ConstRecHitPointer SeedingHit;
0015 
0016 #include <numeric>
0017 
0018 namespace {
0019   std::string seedingLayersToString(const SeedingLayerSetsHits::SeedingLayerSet &layer) {
0020     return layer[0].name() + "+" + layer[1].name() + "+" + layer[2].name();
0021   }
0022 }  // namespace
0023 
0024 using namespace std;
0025 SimpleCosmicBONSeeder::SimpleCosmicBONSeeder(edm::ParameterSet const &conf)
0026     : seedingLayerToken_(consumes<SeedingLayerSetsHits>(conf.getParameter<edm::InputTag>("TripletsSrc"))),
0027       magfieldToken_(esConsumes()),
0028       trackerToken_(esConsumes()),
0029       ttrhBuilderToken_(esConsumes(edm::ESInputTag("", conf.getParameter<std::string>("TTRHBuilder")))),
0030       writeTriplets_(conf.getParameter<bool>("writeTriplets")),
0031       seedOnMiddle_(conf.existsAs<bool>("seedOnMiddle") ? conf.getParameter<bool>("seedOnMiddle") : false),
0032       rescaleError_(conf.existsAs<double>("rescaleError") ? conf.getParameter<double>("rescaleError") : 1.0),
0033       tripletsVerbosity_(conf.getUntrackedParameter<uint32_t>("TripletsDebugLevel", 0)),
0034       seedVerbosity_(conf.getUntrackedParameter<uint32_t>("seedDebugLevel", 0)),
0035       helixVerbosity_(conf.getUntrackedParameter<uint32_t>("helixDebugLevel", 0)),
0036       check_(conf.getParameter<edm::ParameterSet>("ClusterCheckPSet"), consumesCollector()),
0037       maxTriplets_(conf.getParameter<int32_t>("maxTriplets")),
0038       maxSeeds_(conf.getParameter<int32_t>("maxSeeds")) {
0039   edm::ParameterSet regionConf = conf.getParameter<edm::ParameterSet>("RegionPSet");
0040   float ptmin = regionConf.getParameter<double>("ptMin");
0041   float originradius = regionConf.getParameter<double>("originRadius");
0042   float halflength = regionConf.getParameter<double>("originHalfLength");
0043   float originz = regionConf.getParameter<double>("originZPosition");
0044   region_ = GlobalTrackingRegion(ptmin, originradius, halflength, originz);
0045   pMin_ = regionConf.getParameter<double>("pMin");
0046 
0047   //***top-bottom
0048   positiveYOnly = conf.getParameter<bool>("PositiveYOnly");
0049   negativeYOnly = conf.getParameter<bool>("NegativeYOnly");
0050   //***
0051 
0052   produces<TrajectorySeedCollection>();
0053   if (writeTriplets_)
0054     produces<edm::OwnVector<TrackingRecHit>>("cosmicTriplets");
0055 
0056   if (conf.existsAs<edm::ParameterSet>("ClusterChargeCheck")) {
0057     edm::ParameterSet cccc = conf.getParameter<edm::ParameterSet>("ClusterChargeCheck");
0058     checkCharge_ = cccc.getParameter<bool>("checkCharge");
0059     matchedRecHitUsesAnd_ = cccc.getParameter<bool>("matchedRecHitsUseAnd");
0060     chargeThresholds_.resize(7, 0);
0061     edm::ParameterSet ccct = cccc.getParameter<edm::ParameterSet>("Thresholds");
0062     chargeThresholds_[StripSubdetector::TIB] = ccct.getParameter<int32_t>("TIB");
0063     chargeThresholds_[StripSubdetector::TID] = ccct.getParameter<int32_t>("TID");
0064     chargeThresholds_[StripSubdetector::TOB] = ccct.getParameter<int32_t>("TOB");
0065     chargeThresholds_[StripSubdetector::TEC] = ccct.getParameter<int32_t>("TEC");
0066   } else {
0067     checkCharge_ = false;
0068   }
0069   if (conf.existsAs<edm::ParameterSet>("HitsPerModuleCheck")) {
0070     edm::ParameterSet hpmcc = conf.getParameter<edm::ParameterSet>("HitsPerModuleCheck");
0071     checkMaxHitsPerModule_ = hpmcc.getParameter<bool>("checkHitsPerModule");
0072     maxHitsPerModule_.resize(7, std::numeric_limits<int32_t>::max());
0073     edm::ParameterSet hpmct = hpmcc.getParameter<edm::ParameterSet>("Thresholds");
0074     maxHitsPerModule_[StripSubdetector::TIB] = hpmct.getParameter<int32_t>("TIB");
0075     maxHitsPerModule_[StripSubdetector::TID] = hpmct.getParameter<int32_t>("TID");
0076     maxHitsPerModule_[StripSubdetector::TOB] = hpmct.getParameter<int32_t>("TOB");
0077     maxHitsPerModule_[StripSubdetector::TEC] = hpmct.getParameter<int32_t>("TEC");
0078   } else {
0079     checkMaxHitsPerModule_ = false;
0080   }
0081   if (checkCharge_ || checkMaxHitsPerModule_) {
0082     goodHitsPerSeed_ = conf.getParameter<int32_t>("minimumGoodHitsInSeed");
0083   } else {
0084     goodHitsPerSeed_ = 0;
0085   }
0086 }
0087 
0088 // Functions that gets called by framework every event
0089 void SimpleCosmicBONSeeder::produce(edm::Event &ev, const edm::EventSetup &es) {
0090   auto output = std::make_unique<TrajectorySeedCollection>();
0091   auto outtriplets = std::make_unique<edm::OwnVector<TrackingRecHit>>();
0092 
0093   magfield = &es.getData(magfieldToken_);
0094   if (magfield->inTesla(GlobalPoint(0, 0, 0)).mag() > 0.01) {
0095     size_t clustsOrZero = check_.tooManyClusters(ev);
0096     if (clustsOrZero) {
0097       edm::LogError("TooManyClusters") << "Found too many clusters (" << clustsOrZero << "), bailing out.\n";
0098     } else {
0099       init(es);
0100       bool tripletsOk = triplets(ev);
0101       if (tripletsOk) {
0102         bool seedsOk = seeds(*output);
0103         if (!seedsOk) {
0104         }
0105 
0106         if (writeTriplets_) {
0107           for (OrderedHitTriplets::const_iterator it = hitTriplets.begin(); it != hitTriplets.end(); ++it) {
0108             const TrackingRecHit *hit1 = it->inner()->hit();
0109             const TrackingRecHit *hit2 = it->middle()->hit();
0110             const TrackingRecHit *hit3 = it->outer()->hit();
0111             outtriplets->push_back(hit1->clone());
0112             outtriplets->push_back(hit2->clone());
0113             outtriplets->push_back(hit3->clone());
0114           }
0115         }
0116       }
0117       done();
0118     }
0119   }
0120 
0121   if (writeTriplets_) {
0122     ev.put(std::move(outtriplets), "cosmicTriplets");
0123   }
0124   ev.put(std::move(output));
0125 }
0126 
0127 void SimpleCosmicBONSeeder::init(const edm::EventSetup &iSetup) {
0128   tracker = &iSetup.getData(trackerToken_);
0129   cloner = dynamic_cast<TkTransientTrackingRecHitBuilder const &>(iSetup.getData(ttrhBuilderToken_)).cloner();
0130   // FIXME: these should come from ES too!!
0131   thePropagatorAl = new PropagatorWithMaterial(alongMomentum, 0.1057, magfield);
0132   thePropagatorOp = new PropagatorWithMaterial(oppositeToMomentum, 0.1057, magfield);
0133   theUpdator = new KFUpdator();
0134 }
0135 
0136 struct HigherInnerHit {
0137   bool operator()(const OrderedHitTriplet &trip1, const OrderedHitTriplet &trip2) const {
0138     //FIXME: inner gives a SEGV
0139 #if 0
0140         //const SeedingHitSet::ConstRecHitPointer &ihit1 = trip1.inner();
0141         //const SeedingHitSet::ConstRecHitPointer &ihit2 = trip2.inner();
0142         const SeedingHitSet::ConstRecHitPointer &ihit1 = trip1.middle();
0143         const SeedingHitSet::ConstRecHitPointer &ihit2 = trip2.middle();
0144         const SeedingHitSet::ConstRecHitPointer &ohit1 = trip1.outer();
0145         const SeedingHitSet::ConstRecHitPointer &ohit2 = trip2.outer();
0146 #endif
0147     SeedingHitSet::ConstRecHitPointer ihit1 = trip1.inner();
0148     SeedingHitSet::ConstRecHitPointer ihit2 = trip2.inner();
0149     SeedingHitSet::ConstRecHitPointer ohit1 = trip1.outer();
0150     SeedingHitSet::ConstRecHitPointer ohit2 = trip2.outer();
0151     float iy1 = ihit1->globalPosition().y();
0152     float oy1 = ohit1->globalPosition().y();
0153     float iy2 = ihit2->globalPosition().y();
0154     float oy2 = ohit2->globalPosition().y();
0155     if (oy1 - iy1 > 0) {    // 1 Downgoing
0156       if (oy2 - iy2 > 0) {  // 2 Downgoing
0157         // sort by inner, or by outer if inners are the same
0158         return (iy1 != iy2 ? (iy1 > iy2) : (oy1 > oy2));
0159       } else
0160         return true;  // else prefer downgoing
0161     } else if (oy2 - iy2 > 0) {
0162       return false;  // prefer downgoing
0163     } else {         // both upgoing
0164       // sort by inner, or by outer
0165       return (iy1 != iy2 ? (iy1 < iy2) : (oy1 < oy2));
0166     }
0167   }
0168 };
0169 
0170 bool SimpleCosmicBONSeeder::triplets(const edm::Event &e) {
0171   hitTriplets.clear();
0172   hitTriplets.reserve(0);
0173   edm::Handle<SeedingLayerSetsHits> hlayers;
0174   e.getByToken(seedingLayerToken_, hlayers);
0175   const SeedingLayerSetsHits &layers = *hlayers;
0176   if (layers.numberOfLayersInSet() != 3)
0177     throw cms::Exception("CtfSpecialSeedGenerator")
0178         << "You are using " << layers.numberOfLayersInSet() << " layers in set instead of 3 ";
0179 
0180   double minRho = region_.ptMin() / (0.003 * magfield->inTesla(GlobalPoint(0, 0, 0)).z());
0181 
0182   for (SeedingLayerSetsHits::LayerSetIndex layerIndex = 0; layerIndex < layers.size(); ++layerIndex) {
0183     SeedingLayerSetsHits::SeedingLayerSet ls = layers[layerIndex];
0184     /// ctfseeding SeedinHits and their iterators
0185     auto innerHits = region_.hits(ls[0]);
0186     auto middleHits = region_.hits(ls[1]);
0187     auto outerHits = region_.hits(ls[2]);
0188 
0189     if (tripletsVerbosity_ > 0) {
0190       std::cout << "GenericTripletGenerator iLss = " << seedingLayersToString(ls) << " (" << layerIndex
0191                 << "): # = " << innerHits.size() << "/" << middleHits.size() << "/" << outerHits.size() << std::endl;
0192     }
0193 
0194     /// Transient Tracking RecHits (not anymore....)
0195     typedef SeedingHitSet::ConstRecHitPointer TTRH;
0196     std::vector<TTRH> innerTTRHs, middleTTRHs, outerTTRHs;
0197 
0198     /// Checks on the cluster charge and on noisy modules
0199     std::vector<bool> innerOk(innerHits.size(), true);
0200     std::vector<bool> middleOk(middleHits.size(), true);
0201     std::vector<bool> outerOk(outerHits.size(), true);
0202 
0203     size_t sizBefore = hitTriplets.size();
0204     /// Now actually filling in the charges for all the clusters
0205     int idx = 0;
0206     for (auto iOuterHit = outerHits.begin(); iOuterHit != outerHits.end(); ++idx, ++iOuterHit) {
0207       outerTTRHs.push_back(&(**iOuterHit));
0208       if (checkCharge_ && !checkCharge(outerTTRHs.back()->hit()))
0209         outerOk[idx] = false;
0210     }
0211     idx = 0;
0212     for (auto iMiddleHit = middleHits.begin(); iMiddleHit != middleHits.end(); ++idx, ++iMiddleHit) {
0213       middleTTRHs.push_back(&(**iMiddleHit));
0214       if (checkCharge_ && !checkCharge(middleTTRHs.back()->hit()))
0215         middleOk[idx] = false;
0216     }
0217     idx = 0;
0218     for (auto iInnerHit = innerHits.begin(); iInnerHit != innerHits.end(); ++idx, ++iInnerHit) {
0219       innerTTRHs.push_back(&(**iInnerHit));
0220       if (checkCharge_ && !checkCharge(innerTTRHs.back()->hit()))
0221         innerOk[idx] = false;
0222     }
0223     if (checkMaxHitsPerModule_) {
0224       checkNoisyModules(innerTTRHs, innerOk);
0225       checkNoisyModules(middleTTRHs, middleOk);
0226       checkNoisyModules(outerTTRHs, outerOk);
0227     }
0228 
0229     for (auto iOuterHit = outerHits.begin(); iOuterHit != outerHits.end(); iOuterHit++) {
0230       idx = iOuterHit - outerHits.begin();
0231       TTRH &outerTTRH = outerTTRHs[idx];
0232       GlobalPoint outerpos = outerTTRH->globalPosition();  // this caches by itself
0233       bool outerok = outerOk[idx];
0234       if (outerok < goodHitsPerSeed_ - 2) {
0235         if (tripletsVerbosity_ > 2)
0236           std::cout << "Skipping at first hit: " << (outerok) << " < " << (goodHitsPerSeed_ - 2) << std::endl;
0237         continue;
0238       }
0239 
0240       for (auto iMiddleHit = middleHits.begin(); iMiddleHit != middleHits.end(); iMiddleHit++) {
0241         idx = iMiddleHit - middleHits.begin();
0242         TTRH &middleTTRH = middleTTRHs[idx];
0243         GlobalPoint middlepos = middleTTRH->globalPosition();  // this caches by itself
0244         bool middleok = middleOk[idx];
0245         if (outerok + middleok < goodHitsPerSeed_ - 1) {
0246           if (tripletsVerbosity_ > 2)
0247             std::cout << "Skipping at second hit: " << (outerok + middleok) << " < " << (goodHitsPerSeed_ - 1)
0248                       << std::endl;
0249           continue;
0250         }
0251 
0252         for (auto iInnerHit = innerHits.begin(); iInnerHit != innerHits.end(); iInnerHit++) {
0253           idx = iInnerHit - innerHits.begin();
0254           TTRH &innerTTRH = innerTTRHs[idx];
0255           GlobalPoint innerpos = innerTTRH->globalPosition();  // this caches by itself
0256           bool innerok = innerOk[idx];
0257           if (outerok + middleok + innerok < goodHitsPerSeed_) {
0258             if (tripletsVerbosity_ > 2)
0259               std::cout << "Skipping at third hit: " << (outerok + middleok + innerok) << " < " << (goodHitsPerSeed_)
0260                         << std::endl;
0261             continue;
0262           }
0263 
0264           //***top-bottom
0265           if (positiveYOnly &&
0266               (innerpos.y() < 0 || middlepos.y() < 0 || outerpos.y() < 0 || outerpos.y() < innerpos.y()))
0267             continue;
0268           if (negativeYOnly &&
0269               (innerpos.y() > 0 || middlepos.y() > 0 || outerpos.y() > 0 || outerpos.y() > innerpos.y()))
0270             continue;
0271           //***
0272 
0273           if (tripletsVerbosity_ > 2)
0274             std::cout << "Trying seed with: " << innerpos << " + " << middlepos << " + " << outerpos << std::endl;
0275           if (goodTriplet(innerpos, middlepos, outerpos, minRho)) {
0276             OrderedHitTriplet oht(&(**iInnerHit), &(**iMiddleHit), &(**iOuterHit));
0277             hitTriplets.push_back(oht);
0278             if ((maxTriplets_ > 0) && (hitTriplets.size() > size_t(maxTriplets_))) {
0279               hitTriplets.clear();  // clear
0280               //OrderedHitTriplets().swap(hitTriplets); // really clear
0281               edm::LogError("TooManyTriplets") << "Found too many triplets, bailing out.\n";
0282               return false;
0283             }
0284             if (tripletsVerbosity_ > 3) {
0285               std::cout << " accepted seed #" << (hitTriplets.size() - 1) << " w/: " << innerpos << " + " << middlepos
0286                         << " + " << outerpos << std::endl;
0287             }
0288             if (tripletsVerbosity_ == 2) {
0289               std::cout << " good seed #" << (hitTriplets.size() - 1) << " w/: " << innerpos << " + " << middlepos
0290                         << " + " << outerpos << std::endl;
0291             }
0292             if (tripletsVerbosity_ > 3 && (helixVerbosity_ > 0)) {  // debug the momentum here too
0293               pqFromHelixFit(innerpos, middlepos, outerpos);
0294             }
0295           }
0296         }
0297       }
0298     }
0299     if ((tripletsVerbosity_ > 0) && (hitTriplets.size() > sizBefore)) {
0300       std::cout << "                        iLss = " << seedingLayersToString(ls) << " (" << layerIndex
0301                 << "): # = " << innerHits.size() << "/" << middleHits.size() << "/" << outerHits.size() << ": Found "
0302                 << (hitTriplets.size() - sizBefore) << " seeds [running total: " << hitTriplets.size() << "]"
0303                 << std::endl;
0304     }
0305   }
0306   std::sort(hitTriplets.begin(), hitTriplets.end(), HigherInnerHit());
0307   return true;
0308 }
0309 bool SimpleCosmicBONSeeder::checkCharge(const TrackingRecHit *hit) const {
0310   DetId detid(hit->geographicalId());
0311   if (detid.det() != DetId::Tracker)
0312     return false;  // should not happen
0313   int subdet = detid.subdetId();
0314   if (subdet < 3) {  // pixel
0315     return true;
0316   } else {
0317     if (typeid(*hit) == typeid(SiStripMatchedRecHit2D)) {
0318       const SiStripMatchedRecHit2D *mhit = static_cast<const SiStripMatchedRecHit2D *>(hit);
0319       if (matchedRecHitUsesAnd_) {
0320         return checkCharge(mhit->monoHit(), subdet) && checkCharge(mhit->stereoHit(), subdet);
0321       } else {
0322         return checkCharge(mhit->monoHit(), subdet) || checkCharge(mhit->stereoHit(), subdet);
0323       }
0324     } else if (typeid(*hit) == typeid(SiStripRecHit2D)) {
0325       return checkCharge(static_cast<const SiStripRecHit2D &>(*hit), subdet);
0326     } else {
0327       return true;
0328     }
0329   }
0330 }
0331 
0332 // to be fixed to use OmniCluster
0333 bool SimpleCosmicBONSeeder::checkCharge(const SiStripRecHit2D &hit, int subdetid) const {
0334   const SiStripCluster *clust = hit.cluster().get();
0335   int charge = std::accumulate(clust->amplitudes().begin(), clust->amplitudes().end(), int(0));
0336   if (tripletsVerbosity_ > 1) {
0337     std::cerr << "Hit on " << subdetid << ", charge = " << charge << ", threshold = " << chargeThresholds_[subdetid]
0338               << ", detid = " << hit.geographicalId().rawId() << ", firstStrip = " << clust->firstStrip() << std::endl;
0339   } else if ((tripletsVerbosity_ == 1) && (charge < chargeThresholds_[subdetid])) {
0340     std::cerr << "Hit on " << subdetid << ", charge = " << charge << ", threshold = " << chargeThresholds_[subdetid]
0341               << ", detid = " << hit.geographicalId().rawId() << ", firstStrip = " << clust->firstStrip() << std::endl;
0342   }
0343   return charge > chargeThresholds_[subdetid];
0344 }
0345 
0346 void SimpleCosmicBONSeeder::checkNoisyModules(const std::vector<SeedingHitSet::ConstRecHitPointer> &hits,
0347                                               std::vector<bool> &oks) const {
0348   typedef SeedingHitSet::ConstRecHitPointer TTRH;
0349   std::vector<TTRH>::const_iterator it = hits.begin(), start = it, end = hits.end();
0350   std::vector<bool>::iterator ok = oks.begin(), okStart = ok;
0351   while (start < end) {
0352     DetId lastid = (*start)->geographicalId();
0353     for (it = start + 1; (it < end) && ((*it)->geographicalId() == lastid); ++it) {
0354       ++ok;
0355     }
0356     if ((it - start) > maxHitsPerModule_[lastid.subdetId()]) {
0357       if (tripletsVerbosity_ > 0) {
0358         std::cerr << "SimpleCosmicBONSeeder: Marking noisy module " << lastid.rawId() << ", it has " << (it - start)
0359                   << " rechits"
0360                   << " (threshold is " << maxHitsPerModule_[lastid.subdetId()] << ")" << std::endl;
0361       }
0362       std::fill(okStart, ok, false);
0363     } else if (tripletsVerbosity_ > 0) {
0364       if ((it - start) > std::min(4, maxHitsPerModule_[lastid.subdetId()] / 4)) {
0365         std::cerr << "SimpleCosmicBONSeeder: Not marking noisy module " << lastid.rawId() << ", it has " << (it - start)
0366                   << " rechits"
0367                   << " (threshold is " << maxHitsPerModule_[lastid.subdetId()] << ")" << std::endl;
0368       }
0369     }
0370     start = it;
0371     okStart = ok;
0372   }
0373 }
0374 
0375 bool SimpleCosmicBONSeeder::goodTriplet(const GlobalPoint &inner,
0376                                         const GlobalPoint &middle,
0377                                         const GlobalPoint &outer,
0378                                         const double &minRho) const {
0379   float dyOM = outer.y() - middle.y(), dyIM = inner.y() - middle.y();
0380   if ((dyOM * dyIM > 0) && (fabs(dyOM) > 10) && (fabs(dyIM) > 10)) {
0381     if (tripletsVerbosity_ > 2)
0382       std::cout << "  fail for non coherent dy" << std::endl;
0383     return false;
0384   }
0385   float dzOM = outer.z() - middle.z(), dzIM = inner.z() - middle.z();
0386   if ((dzOM * dzIM > 0) && (fabs(dzOM) > 50) && (fabs(dzIM) > 50)) {
0387     if (tripletsVerbosity_ > 2)
0388       std::cout << "  fail for non coherent dz" << std::endl;
0389     return false;
0390   }
0391   if (minRho > 0) {
0392     FastCircle theCircle(inner, middle, outer);
0393     if (theCircle.rho() < minRho) {
0394       if (tripletsVerbosity_ > 2)
0395         std::cout << "  fail for pt cut" << std::endl;
0396       return false;
0397     }
0398   }
0399   return true;
0400 }
0401 
0402 std::pair<GlobalVector, int> SimpleCosmicBONSeeder::pqFromHelixFit(const GlobalPoint &inner,
0403                                                                    const GlobalPoint &middle,
0404                                                                    const GlobalPoint &outer) const {
0405   if (helixVerbosity_ > 0) {
0406     std::cout << "DEBUG PZ =====" << std::endl;
0407     FastHelix helix(inner, middle, outer, magfield->nominalValue(), &*magfield);
0408     GlobalVector gv = helix.stateAtVertex().momentum();  // status on inner hit
0409     std::cout << "FastHelix P = " << gv << "\n";
0410     std::cout << "FastHelix Q = " << helix.stateAtVertex().charge() << "\n";
0411   }
0412 
0413   // My attempt (with different approx from FastHelix)
0414   // 1) fit the circle
0415   FastCircle theCircle(inner, middle, outer);
0416   double rho = theCircle.rho();
0417 
0418   // 2) Get the PT
0419   GlobalVector tesla = magfield->inTesla(middle);
0420   double pt = 0.01 * rho * (0.3 * tesla.z());
0421 
0422   // 3) Get the PX,PY at OUTER hit (VERTEX)
0423   double dx1 = outer.x() - theCircle.x0();
0424   double dy1 = outer.y() - theCircle.y0();
0425   double py = pt * dx1 / rho, px = -pt * dy1 / rho;
0426   if (px * (middle.x() - outer.x()) + py * (middle.y() - outer.y()) < 0.) {
0427     px *= -1.;
0428     py *= -1.;
0429   }
0430 
0431   // 4) Get the PZ through pz = pT*(dz/d(R*phi)))
0432   double dz = inner.z() - outer.z();
0433   double sinphi = (dx1 * (inner.y() - theCircle.y0()) - dy1 * (inner.x() - theCircle.x0())) / (rho * rho);
0434   double dphi = std::abs(std::asin(sinphi));
0435   double pz = pt * dz / (dphi * rho);
0436 
0437   int myq = ((theCircle.x0() * py - theCircle.y0() * px) / tesla.z()) > 0. ? +1 : -1;
0438 
0439   std::pair<GlobalVector, int> mypq(GlobalVector(px, py, pz), myq);
0440 
0441   if (helixVerbosity_ > 1) {
0442     std::cout << "Gio: pt = " << pt << std::endl;
0443     std::cout << "Gio: dz = " << dz << ", sinphi = " << sinphi << ", dphi = " << dphi
0444               << ", dz/drphi = " << (dz / dphi / rho) << std::endl;
0445   }
0446   if (helixVerbosity_ > 0) {
0447     std::cout << "Gio's fit P = " << mypq.first << "\n";
0448     std::cout << "Gio's fit Q = " << myq << "\n";
0449   }
0450 
0451   return mypq;
0452 }
0453 
0454 bool SimpleCosmicBONSeeder::seeds(TrajectorySeedCollection &output) {
0455   typedef TrajectoryStateOnSurface TSOS;
0456 
0457   for (size_t it = 0; it < hitTriplets.size(); it++) {
0458     OrderedHitTriplet &trip = hitTriplets.at(it);
0459 
0460     GlobalPoint inner =
0461         tracker->idToDet((*(trip.inner())).geographicalId())->surface().toGlobal((*(trip.inner())).localPosition());
0462 
0463     GlobalPoint middle =
0464         tracker->idToDet((*(trip.middle())).geographicalId())->surface().toGlobal((*(trip.middle())).localPosition());
0465 
0466     GlobalPoint outer =
0467         tracker->idToDet((*(trip.outer())).geographicalId())->surface().toGlobal((*(trip.outer())).localPosition());
0468 
0469     if (seedVerbosity_ > 1)
0470       std::cout << "Processing triplet " << it << ": " << inner << " + " << middle << " + " << outer << std::endl;
0471 
0472     if ((outer.y() - inner.y()) * outer.y() < 0) {
0473       std::swap(inner, outer);
0474       trip = OrderedHitTriplet(trip.outer(), trip.middle(), trip.inner());
0475 
0476       if (seedVerbosity_ > 1) {
0477         std::cout << "The seed was going away from CMS! swapped in <-> out" << std::endl;
0478         std::cout << "Processing swapped triplet " << it << ": " << inner << " + " << middle << " + " << outer
0479                   << std::endl;
0480       }
0481     }
0482 
0483     // First use FastHelix out of the box
0484     std::pair<GlobalVector, int> pq = pqFromHelixFit(inner, middle, outer);
0485     GlobalVector gv = pq.first;
0486     float ch = pq.second;
0487     float Mom = sqrt(gv.x() * gv.x() + gv.y() * gv.y() + gv.z() * gv.z());
0488 
0489     if (Mom > 10000 || edm::isNotFinite(Mom)) {
0490       if (seedVerbosity_ > 1)
0491         std::cout << "Processing triplet " << it << ": fail for momentum." << std::endl;
0492       continue;
0493     }
0494 
0495     if (gv.perp() < region_.ptMin()) {
0496       if (seedVerbosity_ > 1)
0497         std::cout << "Processing triplet " << it << ": fail for pt = " << gv.perp() << " < ptMin = " << region_.ptMin()
0498                   << std::endl;
0499       continue;
0500     }
0501 
0502     const Propagator *propagator = nullptr;
0503     if ((outer.y() - inner.y()) > 0) {
0504       if (seedVerbosity_ > 1)
0505         std::cout << "Processing triplet " << it << ":  downgoing." << std::endl;
0506       propagator = thePropagatorAl;
0507     } else {
0508       gv = -1 * gv;
0509       ch = -1. * ch;
0510       propagator = thePropagatorOp;
0511       if (seedVerbosity_ > 1)
0512         std::cout << "Processing triplet " << it << ":  upgoing." << std::endl;
0513     }
0514 
0515     if (seedVerbosity_ > 1) {
0516       if ((gv.z() * (outer.z() - inner.z()) > 0) && (fabs(outer.z() - inner.z()) > 5) && (fabs(gv.z()) > .01)) {
0517         std::cout << "ORRORE: outer.z()-inner.z() = " << (outer.z() - inner.z()) << ", gv.z() = " << gv.z()
0518                   << std::endl;
0519       }
0520     }
0521 
0522     GlobalTrajectoryParameters Gtp(outer, gv, int(ch), &(*magfield));
0523     FreeTrajectoryState CosmicSeed(Gtp, CurvilinearTrajectoryError(AlgebraicSymMatrix55(AlgebraicMatrixID())));
0524     CosmicSeed.rescaleError(100);
0525     if (seedVerbosity_ > 2) {
0526       std::cout << "Processing triplet " << it << ". start from " << std::endl;
0527       std::cout << "    X  = " << outer << ", P = " << gv << std::endl;
0528       std::cout << "    Cartesian error (X,P) = \n" << CosmicSeed.cartesianError().matrix() << std::endl;
0529     }
0530 
0531     edm::OwnVector<TrackingRecHit> hits;
0532     OrderedHitTriplet seedHits(trip.outer(), trip.middle(), trip.inner());
0533     TSOS propagated, updated;
0534     bool fail = false;
0535     for (size_t ih = 0; ih < 3; ++ih) {
0536       if ((ih == 2) && seedOnMiddle_) {
0537         if (seedVerbosity_ > 2)
0538           std::cout << "Stopping at middle hit, as requested." << std::endl;
0539         break;
0540       }
0541       if (seedVerbosity_ > 2)
0542         std::cout << "Processing triplet " << it << ", hit " << ih << "." << std::endl;
0543       if (ih == 0) {
0544         propagated = propagator->propagate(CosmicSeed, tracker->idToDet((*seedHits[ih]).geographicalId())->surface());
0545       } else {
0546         propagated = propagator->propagate(updated, tracker->idToDet((*seedHits[ih]).geographicalId())->surface());
0547       }
0548       if (!propagated.isValid()) {
0549         if (seedVerbosity_ > 1)
0550           std::cout << "Processing triplet " << it << ", hit " << ih << ": failed propagation." << std::endl;
0551         fail = true;
0552         break;
0553       } else {
0554         if (seedVerbosity_ > 2)
0555           std::cout << "Processing triplet " << it << ", hit " << ih << ": propagated state = " << propagated;
0556       }
0557       SeedingHitSet::ConstRecHitPointer tthp = seedHits[ih];
0558       auto newtth = static_cast<SeedingHitSet::RecHitPointer>(cloner(*tthp, propagated));
0559       updated = theUpdator->update(propagated, *newtth);
0560       hits.push_back(newtth);
0561       if (!updated.isValid()) {
0562         if (seedVerbosity_ > 1)
0563           std::cout << "Processing triplet " << it << ", hit " << ih << ": failed update." << std::endl;
0564         fail = true;
0565         break;
0566       } else {
0567         if (seedVerbosity_ > 2)
0568           std::cout << "Processing triplet " << it << ", hit " << ih << ": updated state = " << updated;
0569       }
0570     }
0571     if (!fail && updated.isValid() && (updated.globalMomentum().perp() < region_.ptMin())) {
0572       if (seedVerbosity_ > 1)
0573         std::cout << "Processing triplet " << it << ": failed for final pt " << updated.globalMomentum().perp() << " < "
0574                   << region_.ptMin() << std::endl;
0575       fail = true;
0576     }
0577     if (!fail && updated.isValid() && (updated.globalMomentum().mag() < pMin_)) {
0578       if (seedVerbosity_ > 1)
0579         std::cout << "Processing triplet " << it << ": failed for final p " << updated.globalMomentum().perp() << " < "
0580                   << pMin_ << std::endl;
0581       fail = true;
0582     }
0583     if (!fail) {
0584       if (rescaleError_ != 1.0) {
0585         if (seedVerbosity_ > 2) {
0586           std::cout << "Processing triplet " << it << ", rescale error by " << rescaleError_
0587                     << ": state BEFORE rescaling " << updated;
0588           std::cout << "    Cartesian error (X,P) before rescaling= \n"
0589                     << updated.cartesianError().matrix() << std::endl;
0590         }
0591         updated.rescaleError(rescaleError_);
0592       }
0593       if (seedVerbosity_ > 0) {
0594         std::cout << "Processed  triplet " << it << ": success (saved as #" << output.size() << ") : " << inner << " + "
0595                   << middle << " + " << outer << std::endl;
0596         std::cout << "    pt = " << updated.globalMomentum().perp() << "    eta = " << updated.globalMomentum().eta()
0597                   << "    phi = " << updated.globalMomentum().phi() << "    ch = " << updated.charge() << std::endl;
0598         if (seedVerbosity_ > 1) {
0599           std::cout << "    State:" << updated;
0600         } else {
0601           std::cout << "    X  = " << updated.globalPosition() << ", P = " << updated.globalMomentum() << std::endl;
0602         }
0603         std::cout << "    Cartesian error (X,P) = \n" << updated.cartesianError().matrix() << std::endl;
0604       }
0605 
0606       PTrajectoryStateOnDet const &PTraj = trajectoryStateTransform::persistentState(
0607           updated, (*(seedOnMiddle_ ? trip.middle() : trip.inner())).geographicalId().rawId());
0608       output.push_back(TrajectorySeed(PTraj, hits, ((outer.y() - inner.y() > 0) ? alongMomentum : oppositeToMomentum)));
0609       if ((maxSeeds_ > 0) && (output.size() > size_t(maxSeeds_))) {
0610         output.clear();
0611         edm::LogError("TooManySeeds") << "Found too many seeds, bailing out.\n";
0612         return false;
0613       }
0614     }
0615   }
0616   return true;
0617 }
0618 
0619 void SimpleCosmicBONSeeder::done() {
0620   delete thePropagatorAl;
0621   delete thePropagatorOp;
0622   delete theUpdator;
0623 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!