00001
00002
00003
00004
00005
00006
00007
00008 #include "RecoTracker/SpecialSeedGenerators/interface/SimpleCosmicBONSeeder.h"
00009 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
00010 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
00011 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h"
00012 #include "RecoTracker/TkSeedGenerator/interface/FastLine.h"
00013 typedef TransientTrackingRecHit::ConstRecHitPointer SeedingHit;
00014
00015 #include <numeric>
00016
00017 using namespace std;
00018 SimpleCosmicBONSeeder::SimpleCosmicBONSeeder(edm::ParameterSet const& conf) :
00019 conf_(conf),
00020 theLsb(conf.getParameter<edm::ParameterSet>("TripletsPSet")),
00021 writeTriplets_(conf.getParameter<bool>("writeTriplets")),
00022 seedOnMiddle_(conf.existsAs<bool>("seedOnMiddle") ? conf.getParameter<bool>("seedOnMiddle") : false),
00023 rescaleError_(conf.existsAs<double>("rescaleError") ? conf.getParameter<double>("rescaleError") : 1.0),
00024 tripletsVerbosity_(conf.getParameter<edm::ParameterSet>("TripletsPSet").getUntrackedParameter<uint32_t>("debugLevel",0)),
00025 seedVerbosity_(conf.getUntrackedParameter<uint32_t>("seedDebugLevel",0)),
00026 helixVerbosity_(conf.getUntrackedParameter<uint32_t>("helixDebugLevel",0)),
00027 check_(conf.getParameter<edm::ParameterSet>("ClusterCheckPSet")),
00028 maxTriplets_(conf.getParameter<int32_t>("maxTriplets")),
00029 maxSeeds_(conf.getParameter<int32_t>("maxSeeds"))
00030 {
00031 edm::ParameterSet regionConf = conf_.getParameter<edm::ParameterSet>("RegionPSet");
00032 float ptmin = regionConf.getParameter<double>("ptMin");
00033 float originradius = regionConf.getParameter<double>("originRadius");
00034 float halflength = regionConf.getParameter<double>("originHalfLength");
00035 float originz = regionConf.getParameter<double>("originZPosition");
00036 region_ = GlobalTrackingRegion(ptmin, originradius, halflength, originz);
00037 pMin_ = regionConf.getParameter<double>("pMin");
00038
00039 builderName = conf_.getParameter<std::string>("TTRHBuilder");
00040
00041
00042 positiveYOnly=conf_.getParameter<bool>("PositiveYOnly");
00043 negativeYOnly=conf_.getParameter<bool>("NegativeYOnly");
00044
00045
00046 produces<TrajectorySeedCollection>();
00047 if (writeTriplets_) produces<edm::OwnVector<TrackingRecHit> >("cosmicTriplets");
00048
00049 layerTripletNames_ = conf_.getParameter<edm::ParameterSet>("TripletsPSet").getParameter<std::vector<std::string> >("layerList");
00050
00051 if (conf.existsAs<edm::ParameterSet>("ClusterChargeCheck")) {
00052 edm::ParameterSet cccc = conf.getParameter<edm::ParameterSet>("ClusterChargeCheck");
00053 checkCharge_ = cccc.getParameter<bool>("checkCharge");
00054 matchedRecHitUsesAnd_ = cccc.getParameter<bool>("matchedRecHitsUseAnd");
00055 chargeThresholds_.resize(7,0);
00056 edm::ParameterSet ccct = cccc.getParameter<edm::ParameterSet>("Thresholds");
00057 chargeThresholds_[StripSubdetector::TIB] = ccct.getParameter<int32_t>("TIB");
00058 chargeThresholds_[StripSubdetector::TID] = ccct.getParameter<int32_t>("TID");
00059 chargeThresholds_[StripSubdetector::TOB] = ccct.getParameter<int32_t>("TOB");
00060 chargeThresholds_[StripSubdetector::TEC] = ccct.getParameter<int32_t>("TEC");
00061 } else {
00062 checkCharge_ = false;
00063 }
00064 if (conf.existsAs<edm::ParameterSet>("HitsPerModuleCheck")) {
00065 edm::ParameterSet hpmcc = conf.getParameter<edm::ParameterSet>("HitsPerModuleCheck");
00066 checkMaxHitsPerModule_ = hpmcc.getParameter<bool>("checkHitsPerModule");
00067 maxHitsPerModule_.resize(7,std::numeric_limits<int32_t>::max());
00068 edm::ParameterSet hpmct = hpmcc.getParameter<edm::ParameterSet>("Thresholds");
00069 maxHitsPerModule_[StripSubdetector::TIB] = hpmct.getParameter<int32_t>("TIB");
00070 maxHitsPerModule_[StripSubdetector::TID] = hpmct.getParameter<int32_t>("TID");
00071 maxHitsPerModule_[StripSubdetector::TOB] = hpmct.getParameter<int32_t>("TOB");
00072 maxHitsPerModule_[StripSubdetector::TEC] = hpmct.getParameter<int32_t>("TEC");
00073 } else {
00074 checkMaxHitsPerModule_ = false;
00075 }
00076 if (checkCharge_ || checkMaxHitsPerModule_) {
00077 goodHitsPerSeed_ = conf.getParameter<int32_t>("minimumGoodHitsInSeed");
00078 } else {
00079 goodHitsPerSeed_ = 0;
00080 }
00081
00082 }
00083
00084
00085 void SimpleCosmicBONSeeder::produce(edm::Event& ev, const edm::EventSetup& es)
00086 {
00087 std::auto_ptr<TrajectorySeedCollection> output(new TrajectorySeedCollection());
00088 std::auto_ptr<edm::OwnVector<TrackingRecHit> > outtriplets(new edm::OwnVector<TrackingRecHit>());
00089
00090 es.get<IdealMagneticFieldRecord>().get(magfield);
00091 if (magfield->inTesla(GlobalPoint(0,0,0)).mag() > 0.01) {
00092 size_t clustsOrZero = check_.tooManyClusters(ev);
00093 if (clustsOrZero) {
00094 edm::LogError("TooManyClusters") << "Found too many clusters (" << clustsOrZero << "), bailing out.\n";
00095 } else {
00096 init(es);
00097 bool tripletsOk = triplets(ev,es);
00098 if (tripletsOk) {
00099
00100 bool seedsOk = seeds(*output,es);
00101 if (!seedsOk) { }
00102
00103 if (writeTriplets_) {
00104 for (OrderedHitTriplets::const_iterator it = hitTriplets.begin(); it != hitTriplets.end(); ++it) {
00105 const TrackingRecHit * hit1 = it->inner()->hit();
00106 const TrackingRecHit * hit2 = it->middle()->hit();
00107 const TrackingRecHit * hit3 = it->outer()->hit();
00108 outtriplets->push_back(hit1->clone());
00109 outtriplets->push_back(hit2->clone());
00110 outtriplets->push_back(hit3->clone());
00111 }
00112 }
00113 }
00114 done();
00115 }
00116 }
00117
00118 if (writeTriplets_) {
00119 ev.put(outtriplets, "cosmicTriplets");
00120 }
00121 ev.put(output);
00122 }
00123
00124 void
00125 SimpleCosmicBONSeeder::init(const edm::EventSetup& iSetup)
00126 {
00127 iSetup.get<TrackerDigiGeometryRecord>().get(tracker);
00128 iSetup.get<TransientRecHitRecord>().get(builderName,TTTRHBuilder);
00129
00130
00131 thePropagatorAl = new PropagatorWithMaterial(alongMomentum,0.1057,&(*magfield) );
00132 thePropagatorOp = new PropagatorWithMaterial(oppositeToMomentum,0.1057,&(*magfield) );
00133 theUpdator = new KFUpdator();
00134
00135 }
00136
00137 struct HigherInnerHit {
00138 bool operator()(const OrderedHitTriplet &trip1, const OrderedHitTriplet &trip2) const {
00139
00140 #if 0
00141
00142
00143 const TransientTrackingRecHit::ConstRecHitPointer &ihit1 = trip1.middle();
00144 const TransientTrackingRecHit::ConstRecHitPointer &ihit2 = trip2.middle();
00145 const TransientTrackingRecHit::ConstRecHitPointer &ohit1 = trip1.outer();
00146 const TransientTrackingRecHit::ConstRecHitPointer &ohit2 = trip2.outer();
00147 #endif
00148 TransientTrackingRecHit::ConstRecHitPointer ihit1 = trip1.inner();
00149 TransientTrackingRecHit::ConstRecHitPointer ihit2 = trip2.inner();
00150 TransientTrackingRecHit::ConstRecHitPointer ohit1 = trip1.outer();
00151 TransientTrackingRecHit::ConstRecHitPointer ohit2 = trip2.outer();
00152 float iy1 = ihit1->globalPosition().y();
00153 float oy1 = ohit1->globalPosition().y();
00154 float iy2 = ihit2->globalPosition().y();
00155 float oy2 = ohit2->globalPosition().y();
00156 if (oy1 - iy1 > 0) {
00157 if (oy2 - iy2 > 0) {
00158
00159 return (iy1 != iy2 ? (iy1 > iy2) : (oy1 > oy2));
00160 } else return true;
00161 } else if (oy2 - iy2 > 0) {
00162 return false;
00163 } else {
00164
00165 return (iy1 != iy2 ? (iy1 < iy2) : (oy1 < oy2));
00166 }
00167 }
00168 };
00169
00170 bool SimpleCosmicBONSeeder::triplets(const edm::Event& e, const edm::EventSetup& es) {
00171 using namespace ctfseeding;
00172
00173 hitTriplets.clear();
00174 hitTriplets.reserve(0);
00175 SeedingLayerSets lss = theLsb.layers(es);
00176 SeedingLayerSets::const_iterator iLss;
00177
00178 double minRho = region_.ptMin() / ( 0.003 * magfield->inTesla(GlobalPoint(0,0,0)).z() );
00179
00180 for (iLss = lss.begin(); iLss != lss.end(); iLss++){
00181 SeedingLayers ls = *iLss;
00182 if (ls.size() != 3){
00183 throw cms::Exception("CtfSpecialSeedGenerator") << "You are using " << ls.size() <<" layers in set instead of 3 ";
00184 }
00185
00187 std::vector<SeedingHit> innerHits = region_.hits(e, es, &ls[0]);
00188 std::vector<SeedingHit> middleHits = region_.hits(e, es, &ls[1]);
00189 std::vector<SeedingHit> outerHits = region_.hits(e, es, &ls[2]);
00190 std::vector<SeedingHit>::const_iterator iOuterHit,iMiddleHit,iInnerHit;
00191
00192 if (tripletsVerbosity_ > 0) {
00193 std::cout << "GenericTripletGenerator iLss = " << layerTripletNames_[iLss - lss.begin()]
00194 << " (" << (iLss - lss.begin()) << "): # = "
00195 << innerHits.size() << "/" << middleHits.size() << "/" << outerHits.size() << std::endl;
00196 }
00197
00199 typedef TransientTrackingRecHit::RecHitPointer TTRH;
00200 std::vector<TTRH> innerTTRHs, middleTTRHs, outerTTRHs;
00201
00203 std::vector<bool> innerOk( innerHits.size(), true);
00204 std::vector<bool> middleOk(middleHits.size(), true);
00205 std::vector<bool> outerOk( outerHits.size(), true);
00206
00207 size_t sizBefore = hitTriplets.size();
00209 int idx = 0;
00210 for (iOuterHit = outerHits.begin(), idx = 0; iOuterHit != outerHits.end(); ++idx, ++iOuterHit){
00211 outerTTRHs.push_back(ls[2].hitBuilder()->build((**iOuterHit).hit()));
00212 if (checkCharge_ && !checkCharge(outerTTRHs.back()->hit())) outerOk[idx] = false;
00213 }
00214 for (iMiddleHit = middleHits.begin(), idx = 0; iMiddleHit != middleHits.end(); ++idx, ++iMiddleHit){
00215 middleTTRHs.push_back(ls[1].hitBuilder()->build((**iMiddleHit).hit()));
00216 if (checkCharge_ && !checkCharge(middleTTRHs.back()->hit())) middleOk[idx] = false;
00217 }
00218 for (iInnerHit = innerHits.begin(), idx = 0; iInnerHit != innerHits.end(); ++idx, ++iInnerHit){
00219 innerTTRHs.push_back(ls[0].hitBuilder()->build((**iInnerHit).hit()));
00220 if (checkCharge_ && !checkCharge(innerTTRHs.back()->hit())) innerOk[idx] = false;
00221 }
00222 if (checkMaxHitsPerModule_) {
00223 checkNoisyModules(innerTTRHs, innerOk);
00224 checkNoisyModules(middleTTRHs, middleOk);
00225 checkNoisyModules(outerTTRHs, outerOk);
00226 }
00227
00228 for (iOuterHit = outerHits.begin(); iOuterHit != outerHits.end(); iOuterHit++){
00229 idx = iOuterHit - outerHits.begin();
00230 TTRH & outerTTRH = outerTTRHs[idx];
00231 GlobalPoint outerpos = outerTTRH->globalPosition();
00232 bool outerok = outerOk[idx];
00233 if (outerok < goodHitsPerSeed_ - 2) {
00234 if (tripletsVerbosity_ > 2)
00235 std::cout << "Skipping at first hit: " << (outerok) << " < " << (goodHitsPerSeed_ - 2) << std::endl;
00236 continue;
00237 }
00238
00239 for (iMiddleHit = middleHits.begin(); iMiddleHit != middleHits.end(); iMiddleHit++){
00240 idx = iMiddleHit - middleHits.begin();
00241 TTRH & middleTTRH = middleTTRHs[idx];
00242 GlobalPoint middlepos = middleTTRH->globalPosition();
00243 bool middleok = middleOk[idx];
00244 if (outerok+middleok < goodHitsPerSeed_ - 1) {
00245 if (tripletsVerbosity_ > 2)
00246 std::cout << "Skipping at second hit: " << (outerok+middleok) << " < " << (goodHitsPerSeed_ - 1) << std::endl;
00247 continue;
00248 }
00249
00250 for (iInnerHit = innerHits.begin(); iInnerHit != innerHits.end(); iInnerHit++){
00251 idx = iInnerHit - innerHits.begin();
00252 TTRH & innerTTRH = innerTTRHs[idx];
00253 GlobalPoint innerpos = innerTTRH->globalPosition();
00254 bool innerok = innerOk[idx];
00255 if (outerok+middleok+innerok < goodHitsPerSeed_) {
00256 if (tripletsVerbosity_ > 2)
00257 std::cout << "Skipping at third hit: " << (outerok+middleok+innerok) << " < " << (goodHitsPerSeed_) << std::endl;
00258 continue;
00259 }
00260
00261
00262 if (positiveYOnly && (innerpos.y()<0 || middlepos.y()<0 || outerpos.y()<0
00263 || outerpos.y() < innerpos.y()
00264 ) ) continue;
00265 if (negativeYOnly && (innerpos.y()>0 || middlepos.y()>0 || outerpos.y()>0
00266 || outerpos.y() > innerpos.y()
00267 ) ) continue;
00268
00269
00270 if (tripletsVerbosity_ > 2) std::cout << "Trying seed with: " << innerpos << " + " << middlepos << " + " << outerpos << std::endl;
00271 if (goodTriplet(innerpos,middlepos,outerpos,minRho)) {
00272 OrderedHitTriplet oht(*iInnerHit,*iMiddleHit,*iOuterHit);
00273 hitTriplets.push_back(oht);
00274 if ((maxTriplets_ > 0) && (hitTriplets.size() > size_t(maxTriplets_))) {
00275 hitTriplets.clear();
00276
00277 edm::LogError("TooManyTriplets") << "Found too many triplets, bailing out.\n";
00278 return false;
00279 }
00280 if (tripletsVerbosity_ > 3) {
00281 std::cout << " accepted seed #" << (hitTriplets.size()-1) << " w/: "
00282 << innerpos << " + " << middlepos << " + " << outerpos << std::endl;
00283 }
00284 if (tripletsVerbosity_ == 2) {
00285 std::cout << " good seed #" << (hitTriplets.size()-1) << " w/: "
00286 << innerpos << " + " << middlepos << " + " << outerpos << std::endl;
00287 }
00288 if (tripletsVerbosity_ > 3 && (helixVerbosity_ > 0)) {
00289 pqFromHelixFit(innerpos,middlepos,outerpos,es);
00290 }
00291 }
00292 }
00293 }
00294 }
00295 if ((tripletsVerbosity_ > 0) && (hitTriplets.size() > sizBefore)) {
00296 std::cout << " iLss = " << layerTripletNames_[iLss - lss.begin()]
00297 << " (" << (iLss - lss.begin()) << "): # = "
00298 << innerHits.size() << "/" << middleHits.size() << "/" << outerHits.size()
00299 << ": Found " << (hitTriplets.size() - sizBefore) << " seeds [running total: " << hitTriplets.size() << "]"
00300 << std::endl ;
00301 }
00302
00303 }
00304 std::sort(hitTriplets.begin(),hitTriplets.end(),HigherInnerHit());
00305 return true;
00306 }
00307 bool SimpleCosmicBONSeeder::checkCharge(const TrackingRecHit *hit) const {
00308 DetId detid(hit->geographicalId());
00309 if (detid.det() != DetId::Tracker) return false;
00310 int subdet = detid.subdetId();
00311 if (subdet < 3) {
00312 return true;
00313 } else {
00314 if (typeid(*hit) == typeid(SiStripMatchedRecHit2D)) {
00315 const SiStripMatchedRecHit2D *mhit = static_cast<const SiStripMatchedRecHit2D *>(hit);
00316 if (matchedRecHitUsesAnd_) {
00317 return checkCharge(*mhit->monoHit(), subdet) && checkCharge(*mhit->stereoHit(), subdet);
00318 } else {
00319 return checkCharge(*mhit->monoHit(), subdet) || checkCharge(*mhit->stereoHit(), subdet);
00320 }
00321 } else if (typeid(*hit) == typeid(SiStripRecHit2D)) {
00322 return checkCharge(static_cast<const SiStripRecHit2D &>(*hit), subdet);
00323 } else {
00324 return true;
00325 }
00326 }
00327 }
00328 bool SimpleCosmicBONSeeder::checkCharge(const SiStripRecHit2D &hit, int subdetid) const {
00329 const SiStripCluster *clust = (hit.cluster().isNonnull() ? hit.cluster().get() : hit.cluster_regional().get());
00330 int charge = std::accumulate(clust->amplitudes().begin(), clust->amplitudes().end(), int(0));
00331 if (tripletsVerbosity_ > 1) {
00332 std::cerr << "Hit on " << subdetid << ", charge = " << charge << ", threshold = " << chargeThresholds_[subdetid]
00333 << ", detid = " << hit.geographicalId().rawId() << ", firstStrip = " << clust->firstStrip() << std::endl;
00334 } else if ((tripletsVerbosity_ == 1) && (charge < chargeThresholds_[subdetid])) {
00335 std::cerr << "Hit on " << subdetid << ", charge = " << charge << ", threshold = " << chargeThresholds_[subdetid]
00336 << ", detid = " << hit.geographicalId().rawId() << ", firstStrip = " << clust->firstStrip() << std::endl;
00337 }
00338 return charge > chargeThresholds_[subdetid];
00339 }
00340
00341 void SimpleCosmicBONSeeder::checkNoisyModules(const std::vector<TransientTrackingRecHit::RecHitPointer> &hits, std::vector<bool> &oks) const {
00342 typedef TransientTrackingRecHit::RecHitPointer TTRH;
00343 std::vector<TTRH>::const_iterator it = hits.begin(), start = it, end = hits.end();
00344 std::vector<bool>::iterator ok = oks.begin(), okStart = ok, okEnd = oks.end();
00345 while (start < end) {
00346 DetId lastid = (*start)->geographicalId();
00347 for (it = start + 1; (it < end) && ((*it)->geographicalId() == lastid); ++it) {
00348 ++ok;
00349 }
00350 if ( (it - start) > maxHitsPerModule_[lastid.subdetId()] ) {
00351 if (tripletsVerbosity_ > 0) {
00352 std::cerr << "SimpleCosmicBONSeeder: Marking noisy module " << lastid.rawId() << ", it has " << (it-start) << " rechits"
00353 << " (threshold is " << maxHitsPerModule_[lastid.subdetId()] << ")" << std::endl;
00354 }
00355 std::fill(okStart,ok,false);
00356 } else if (tripletsVerbosity_ > 0) {
00357 if ( (it - start) > std::min(4,maxHitsPerModule_[lastid.subdetId()]/4) ) {
00358 std::cerr << "SimpleCosmicBONSeeder: Not marking noisy module " << lastid.rawId() << ", it has " << (it-start) << " rechits"
00359 << " (threshold is " << maxHitsPerModule_[lastid.subdetId()] << ")" << std::endl;
00360 }
00361 }
00362 start = it; okStart = ok;
00363 }
00364 }
00365
00366 bool SimpleCosmicBONSeeder::goodTriplet(const GlobalPoint &inner, const GlobalPoint & middle, const GlobalPoint & outer, const double &minRho) const {
00367 float dyOM = outer.y() - middle.y(), dyIM = inner.y() - middle.y();
00368 if ((dyOM * dyIM > 0) && (fabs(dyOM)>10) && (fabs(dyIM)>10)) {
00369 if (tripletsVerbosity_ > 2) std::cout << " fail for non coherent dy" << std::endl;
00370 return false;
00371 }
00372 float dzOM = outer.z() - middle.z(), dzIM = inner.z() - middle.z();
00373 if ((dzOM * dzIM > 0) && (fabs(dzOM)>50) && (fabs(dzIM)>50)) {
00374 if (tripletsVerbosity_ > 2) std::cout << " fail for non coherent dz" << std::endl;
00375 return false;
00376 }
00377 if (minRho > 0) {
00378 FastCircle theCircle(inner,middle,outer);
00379 if (theCircle.rho() < minRho) {
00380 if (tripletsVerbosity_ > 2) std::cout << " fail for pt cut" << std::endl;
00381 return false;
00382 }
00383 }
00384 return true;
00385 }
00386
00387 std::pair<GlobalVector,int>
00388 SimpleCosmicBONSeeder::pqFromHelixFit(const GlobalPoint &inner, const GlobalPoint & middle, const GlobalPoint & outer,
00389 const edm::EventSetup& iSetup) const {
00390 if (helixVerbosity_ > 0) {
00391 std::cout << "DEBUG PZ =====" << std::endl;
00392 FastHelix helix(inner,middle,outer,iSetup);
00393 GlobalVector gv=helix.stateAtVertex().parameters().momentum();
00394 std::cout << "FastHelix P = " << gv << "\n";
00395 std::cout << "FastHelix Q = " << helix.stateAtVertex().parameters().charge() << "\n";
00396 }
00397
00398
00399
00400 FastCircle theCircle(inner,middle,outer);
00401 double rho = theCircle.rho();
00402
00403
00404 GlobalVector tesla = magfield->inTesla(middle);
00405 double pt = 0.01 * rho * (0.3*tesla.z());
00406
00407
00408 double dx1 = outer.x()-theCircle.x0();
00409 double dy1 = outer.y()-theCircle.y0();
00410 double py = pt*dx1/rho, px = -pt*dy1/rho;
00411 if(px*(middle.x() - outer.x()) + py*(middle.y() - outer.y()) < 0.) {
00412 px *= -1.; py *= -1.;
00413 }
00414
00415
00416 double dz = inner.z() - outer.z();
00417 double sinphi = ( dx1*(inner.y()-theCircle.y0()) - dy1*(inner.x()-theCircle.x0())) / (rho * rho);
00418 double dphi = std::abs(std::asin(sinphi));
00419 double pz = pt * dz / (dphi * rho);
00420
00421 int myq = ((theCircle.x0()*py - theCircle.y0()*px) / tesla.z()) > 0. ? +1 : -1;
00422
00423 std::pair<GlobalVector,int> mypq(GlobalVector(px,py,pz),myq);
00424
00425 if (helixVerbosity_ > 1) {
00426 std::cout << "Gio: pt = " << pt << std::endl;
00427 std::cout << "Gio: dz = " << dz << ", sinphi = " << sinphi << ", dphi = " << dphi << ", dz/drphi = " << (dz/dphi/rho) << std::endl;
00428 }
00429 if (helixVerbosity_ > 0) {
00430 std::cout << "Gio's fit P = " << mypq.first << "\n";
00431 std::cout << "Gio's fit Q = " << myq << "\n";
00432 }
00433
00434 return mypq;
00435 }
00436
00437 bool SimpleCosmicBONSeeder::seeds(TrajectorySeedCollection &output, const edm::EventSetup& iSetup)
00438 {
00439 typedef TrajectoryStateOnSurface TSOS;
00440
00441 for (size_t it=0;it<hitTriplets.size();it++){
00442 const OrderedHitTriplet &trip = hitTriplets[it];
00443
00444 GlobalPoint inner = tracker->idToDet((*(trip.inner())).geographicalId())->surface().
00445 toGlobal((*(trip.inner())).localPosition());
00446
00447 GlobalPoint middle = tracker->idToDet((*(trip.middle())).geographicalId())->surface().
00448 toGlobal((*(trip.middle())).localPosition());
00449
00450 GlobalPoint outer = tracker->idToDet((*(trip.outer())).geographicalId())->surface().
00451 toGlobal((*(trip.outer())).localPosition());
00452
00453 if (seedVerbosity_ > 1)
00454 std::cout << "Processing triplet " << it << ": " << inner << " + " << middle << " + " << outer << std::endl;
00455
00456 if ( (outer.y()-inner.y())*outer.y() < 0 ) {
00457 std::swap(inner,outer);
00458 std::swap(const_cast<TransientTrackingRecHit::ConstRecHitPointer &>(trip.inner()),
00459 const_cast<TransientTrackingRecHit::ConstRecHitPointer &>(trip.outer()) );
00460
00461
00462
00463 if (seedVerbosity_ > 1) {
00464 std::cout << "The seed was going away from CMS! swapped in <-> out" << std::endl;
00465 std::cout << "Processing swapped triplet " << it << ": " << inner << " + " << middle << " + " << outer << std::endl;
00466 }
00467 }
00468
00469
00470 std::pair<GlobalVector,int> pq = pqFromHelixFit(inner,middle,outer,iSetup);
00471 GlobalVector gv = pq.first;
00472 float ch = pq.second;
00473 float Mom = sqrt( gv.x()*gv.x() + gv.y()*gv.y() + gv.z()*gv.z() );
00474
00475 if(Mom > 10000 || isnan(Mom)) {
00476 if (seedVerbosity_ > 1)
00477 std::cout << "Processing triplet " << it << ": fail for momentum." << std::endl;
00478 continue;
00479 }
00480
00481 if (gv.perp() < region_.ptMin()) {
00482 if (seedVerbosity_ > 1)
00483 std::cout << "Processing triplet " << it << ": fail for pt = " << gv.perp() << " < ptMin = " << region_.ptMin() << std::endl;
00484 continue;
00485 }
00486
00487 const Propagator * propagator = 0;
00488 if((outer.y()-inner.y())>0){
00489 if (seedVerbosity_ > 1)
00490 std::cout << "Processing triplet " << it << ": downgoing." << std::endl;
00491 propagator = thePropagatorAl;
00492 } else {
00493 gv = -1*gv; ch = -1.*ch;
00494 propagator = thePropagatorOp;
00495 if (seedVerbosity_ > 1)
00496 std::cout << "Processing triplet " << it << ": upgoing." << std::endl;
00497 }
00498
00499 if (seedVerbosity_ > 1) {
00500 if (( gv.z() * (outer.z()-inner.z()) > 0 ) && ( fabs(outer.z()-inner.z()) > 5) && (fabs(gv.z()) > .01)) {
00501 std::cout << "ORRORE: outer.z()-inner.z() = " << (outer.z()-inner.z()) << ", gv.z() = " << gv.z() << std::endl;
00502 }
00503 }
00504
00505 GlobalTrajectoryParameters Gtp(outer,
00506 gv,int(ch),
00507 &(*magfield));
00508 FreeTrajectoryState CosmicSeed(Gtp,
00509 CurvilinearTrajectoryError(AlgebraicSymMatrix55(AlgebraicMatrixID())));
00510 CosmicSeed.rescaleError(100);
00511 if (seedVerbosity_ > 2) {
00512 std::cout << "Processing triplet " << it << ". start from " << std::endl;
00513 std::cout << " X = " << outer << ", P = " << gv << std::endl;
00514 std::cout << " Cartesian error (X,P) = \n" << CosmicSeed.cartesianError().matrix() << std::endl;
00515 }
00516
00517 edm::OwnVector<TrackingRecHit> hits;
00518
00519 OrderedHitTriplet::RecHits seedHits;
00520 seedHits.push_back(trip.outer());
00521 seedHits.push_back(trip.middle());
00522 seedHits.push_back(trip.inner());
00523 TSOS propagated, updated;
00524 bool fail = false;
00525 for (size_t ih = 0; ih < 3; ++ih) {
00526 if ((ih == 2) && seedOnMiddle_) {
00527 if (seedVerbosity_ > 2)
00528 std::cout << "Stopping at middle hit, as requested." << std::endl;
00529 break;
00530 }
00531 if (seedVerbosity_ > 2)
00532 std::cout << "Processing triplet " << it << ", hit " << ih << "." << std::endl;
00533 if (ih == 0) {
00534 propagated = propagator->propagate(CosmicSeed, tracker->idToDet((*seedHits[ih]).geographicalId())->surface());
00535 } else {
00536 propagated = propagator->propagate(updated, tracker->idToDet((*seedHits[ih]).geographicalId())->surface());
00537 }
00538 if (!propagated.isValid()) {
00539 if (seedVerbosity_ > 1)
00540 std::cout << "Processing triplet " << it << ", hit " << ih << ": failed propagation." << std::endl;
00541 fail = true; break;
00542 } else {
00543 if (seedVerbosity_ > 2)
00544 std::cout << "Processing triplet " << it << ", hit " << ih << ": propagated state = " << propagated;
00545 }
00546 const TransientTrackingRecHit::ConstRecHitPointer & tthp = seedHits[ih];
00547 TransientTrackingRecHit::RecHitPointer newtth = tthp->clone(propagated);
00548 hits.push_back(newtth->hit()->clone());
00549 updated = theUpdator->update(propagated, *newtth);
00550 if (!updated.isValid()) {
00551 if (seedVerbosity_ > 1)
00552 std::cout << "Processing triplet " << it << ", hit " << ih << ": failed update." << std::endl;
00553 fail = true; break;
00554 } else {
00555 if (seedVerbosity_ > 2)
00556 std::cout << "Processing triplet " << it << ", hit " << ih << ": updated state = " << updated;
00557 }
00558 }
00559 if (!fail && updated.isValid() && (updated.globalMomentum().perp() < region_.ptMin())) {
00560 if (seedVerbosity_ > 1)
00561 std::cout << "Processing triplet " << it <<
00562 ": failed for final pt " << updated.globalMomentum().perp() << " < " << region_.ptMin() << std::endl;
00563 fail = true;
00564 }
00565 if (!fail && updated.isValid() && (updated.globalMomentum().mag() < pMin_)) {
00566 if (seedVerbosity_ > 1)
00567 std::cout << "Processing triplet " << it <<
00568 ": failed for final p " << updated.globalMomentum().perp() << " < " << pMin_ << std::endl;
00569 fail = true;
00570 }
00571 if (!fail) {
00572 if (rescaleError_ != 1.0) {
00573 if (seedVerbosity_ > 2) {
00574 std::cout << "Processing triplet " << it << ", rescale error by " << rescaleError_ << ": state BEFORE rescaling " << updated;
00575 std::cout << " Cartesian error (X,P) before rescaling= \n" << updated.cartesianError().matrix() << std::endl;
00576 }
00577 updated.rescaleError(rescaleError_);
00578 }
00579 if (seedVerbosity_ > 0) {
00580 std::cout << "Processed triplet " << it << ": success (saved as #"<<output.size()<<") : "
00581 << inner << " + " << middle << " + " << outer << std::endl;
00582 std::cout << " pt = " << updated.globalMomentum().perp() <<
00583 " eta = " << updated.globalMomentum().eta() <<
00584 " phi = " << updated.globalMomentum().phi() <<
00585 " ch = " << updated.charge() << std::endl;
00586 if (seedVerbosity_ > 1) {
00587 std::cout << " State:" << updated;
00588 } else {
00589 std::cout << " X = " << updated.globalPosition() << ", P = " << updated.globalMomentum() << std::endl;
00590 }
00591 std::cout << " Cartesian error (X,P) = \n" << updated.cartesianError().matrix() << std::endl;
00592 }
00593
00594 std::auto_ptr<PTrajectoryStateOnDet> PTraj(transformer.persistentState(updated,
00595 (*(seedOnMiddle_ ? trip.middle() : trip.inner())).geographicalId().rawId()));
00596 output.push_back(TrajectorySeed(*PTraj,hits,
00597 ( (outer.y()-inner.y()>0) ? alongMomentum : oppositeToMomentum) ));
00598 if ((maxSeeds_ > 0) && (output.size() > size_t(maxSeeds_))) {
00599 output.clear();
00600 edm::LogError("TooManySeeds") << "Found too many seeds, bailing out.\n";
00601 return false;
00602 }
00603 }
00604 }
00605 return true;
00606 }
00607
00608 void SimpleCosmicBONSeeder::done(){
00609 delete thePropagatorAl;
00610 delete thePropagatorOp;
00611 delete theUpdator;
00612 }
00613
00614