CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_4_1_8_patch12/src/RecoTracker/SpecialSeedGenerators/src/SimpleCosmicBONSeeder.cc

Go to the documentation of this file.
00001 //
00002 // Package:         RecoTracker/TkSeedGenerator
00003 // Class:           GlobalPixelLessSeedGenerator
00004 // 
00005 // From CosmicSeedGenerator+SimpleCosmicBONSeeder, with changes by Giovanni
00006 // to seed Cosmics with B != 0
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   //***top-bottom
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 // Functions that gets called by framework every event
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     // FIXME: these should come from ES too!!
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         //FIXME: inner gives a SEGV
00140 #if 0
00141         //const TransientTrackingRecHit::ConstRecHitPointer &ihit1 = trip1.inner();
00142         //const TransientTrackingRecHit::ConstRecHitPointer &ihit2 = trip2.inner();
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) { // 1 Downgoing
00157             if (oy2 - iy2 > 0) { // 2 Downgoing
00158                 // sort by inner, or by outer if inners are the same
00159                 return (iy1 != iy2 ? (iy1 > iy2) : (oy1 > oy2));
00160             } else return true; // else prefer downgoing
00161         } else if (oy2 - iy2 > 0) {
00162             return false; // prefer downgoing
00163         } else { // both upgoing
00164             // sort by inner, or by outer
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(); // this caches by itself
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(); // this caches by itself
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(); // this caches by itself
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                     //***top-bottom
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();                      // clear
00276                             //OrderedHitTriplets().swap(hitTriplets); // really clear   
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)) { // debug the momentum here too
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; // should not happen
00310     int subdet = detid.subdetId(); 
00311     if (subdet < 3) { // pixel
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(); // status on inner hit
00394         std::cout << "FastHelix P = " << gv   << "\n";
00395         std::cout << "FastHelix Q = " << helix.stateAtVertex().parameters().charge() << "\n";
00396     }
00397 
00398     // My attempt (with different approx from FastHelix)
00399     // 1) fit the circle
00400     FastCircle theCircle(inner,middle,outer);
00401     double rho = theCircle.rho();
00402 
00403     // 2) Get the PT
00404     GlobalVector tesla = magfield->inTesla(middle);
00405     double pt = 0.01 * rho * (0.3*tesla.z());
00406 
00407     // 3) Get the PX,PY at OUTER hit (VERTEX)
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     // 4) Get the PZ through pz = pT*(dz/d(R*phi)))
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 //            std::swap(const_cast<ctfseeding::SeedingHit &>(trip.inner()), 
00462 //                      const_cast<ctfseeding::SeedingHit &>(trip.outer()) );
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         // First use FastHelix out of the box
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 //      OrderedHitTriplet::Hits seedHits;
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