DTSegment/src/DTCombinatorialPatternReco4D.cc

0001 /** \file
0002  *
0003  * \author Stefano Lacaprara - INFN Legnaro <stefano.lacaprara@pd.infn.it>
0004  * \author Riccardo Bellan - INFN TO <riccardo.bellan@cern.ch>
0005  */
0006 #include "RecoLocalMuon/DTSegment/src/DTCombinatorialPatternReco4D.h"
0007
0008 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0009 #include "FWCore/Framework/interface/EventSetup.h"
0010
0011 #include "RecoLocalMuon/DTSegment/src/DTSegmentUpdator.h"
0012 // For the 2D reco I use thei reconstructor!
0013 #include "RecoLocalMuon/DTSegment/src/DTCombinatorialPatternReco.h"
0014 #include "RecoLocalMuon/DTSegment/src/DTSegmentCand.h"
0015
0016 #include "DataFormats/Common/interface/OwnVector.h"
0017 #include "DataFormats/DTRecHit/interface/DTRangeMapAccessor.h"
0018 #include "DataFormats/MuonDetId/interface/DTWireId.h"
0019 #include "DataFormats/DTRecHit/interface/DTRecHit1DPair.h"
0020 #include "DataFormats/DTRecHit/interface/DTSLRecSegment2D.h"
0021 #include "DataFormats/DTRecHit/interface/DTChamberRecSegment2D.h"
0022
0023 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0024
0025 using namespace std;
0026 using namespace edm;
0027
0028 // TODO
0029 // Throw an exception if a theta segment container is requested and in the event
0030 // there isn't it. (Or launch a "lazy" reco on demand)
0031
0032 DTCombinatorialPatternReco4D::DTCombinatorialPatternReco4D(const ParameterSet& pset, ConsumesCollector cc)
0033     : DTRecSegment4DBaseAlgo(pset), theAlgoName("DTCombinatorialPatternReco4D"), theDTGeometryToken(cc.esConsumes()) {
0034   // debug parameter
0035   debug = pset.getUntrackedParameter<bool>("debug");
0036
0037   //do you want the T0 correction?
0038   applyT0corr = pset.getParameter<bool>("performT0SegCorrection");
0039
0040   computeT0corr = pset.existsAs<bool>("computeT0Seg") ? pset.getParameter<bool>("computeT0Seg") : true;
0041
0042   // the updator
0043   theUpdator = new DTSegmentUpdator(pset, cc);
0044
0045   // the input type.
0046   // If true the instructions in setDTRecSegment2DContainer will be schipped and the
0047   // theta segment will be recomputed from the 1D rechits
0048   // If false the theta segment will be taken from the Event. Caveat: in this case the
0049   // event must contain the 2D segments!
0050   allDTRecHits = pset.getParameter<bool>("AllDTRecHits");
0051
0052   // Get the concrete 2D-segments reconstruction algo from the factory
0053   // For the 2D reco I use this reconstructor!
0054   the2DAlgo = new DTCombinatorialPatternReco(pset.getParameter<ParameterSet>("Reco2DAlgoConfig"), cc);
0055 }
0056
0057 DTCombinatorialPatternReco4D::~DTCombinatorialPatternReco4D() {
0058   delete the2DAlgo;
0059   delete theUpdator;
0060 }
0061
0062 void DTCombinatorialPatternReco4D::setES(const EventSetup& setup) {
0063   theDTGeometry = setup.getHandle(theDTGeometryToken);
0064   the2DAlgo->setES(setup);
0065   theUpdator->setES(setup);
0066 }
0067
0068 void DTCombinatorialPatternReco4D::setChamber(const DTChamberId& chId) {
0069   // Set the chamber
0070   theChamber = theDTGeometry->chamber(chId);
0071 }
0072
0073 void DTCombinatorialPatternReco4D::setDTRecHit1DContainer(Handle<DTRecHitCollection> all1DHits) {
0074   theHitsFromPhi1.clear();
0075   theHitsFromPhi2.clear();
0076   theHitsFromTheta.clear();
0077
0078   DTRecHitCollection::range rangeHitsFromPhi1 =
0079       all1DHits->get(DTRangeMapAccessor::layersBySuperLayer(DTSuperLayerId(theChamber->id(), 1)));
0080   DTRecHitCollection::range rangeHitsFromPhi2 =
0081       all1DHits->get(DTRangeMapAccessor::layersBySuperLayer(DTSuperLayerId(theChamber->id(), 3)));
0082
0083   vector<DTRecHit1DPair> hitsFromPhi1(rangeHitsFromPhi1.first, rangeHitsFromPhi1.second);
0084   vector<DTRecHit1DPair> hitsFromPhi2(rangeHitsFromPhi2.first, rangeHitsFromPhi2.second);
0085   if (debug)
0086     cout << "Number of DTRecHit1DPair in the SL 1 (Phi 1): " << hitsFromPhi1.size() << endl
0087          << "Number of DTRecHit1DPair in the SL 3 (Phi 2): " << hitsFromPhi2.size() << endl;
0088
0089   theHitsFromPhi1 = hitsFromPhi1;
0090   theHitsFromPhi2 = hitsFromPhi2;
0091
0092   if (allDTRecHits) {
0093     DTRecHitCollection::range rangeHitsFromTheta =
0094         all1DHits->get(DTRangeMapAccessor::layersBySuperLayer(DTSuperLayerId(theChamber->id(), 2)));
0095
0096     vector<DTRecHit1DPair> hitsFromTheta(rangeHitsFromTheta.first, rangeHitsFromTheta.second);
0097     if (debug)
0098       cout << "Number of DTRecHit1DPair in the SL 2 (Theta): " << hitsFromTheta.size() << endl;
0099     theHitsFromTheta = hitsFromTheta;
0100   }
0101 }
0102
0103 void DTCombinatorialPatternReco4D::setDTRecSegment2DContainer(Handle<DTRecSegment2DCollection> all2DSegments) {
0104   theSegments2DTheta.clear();
0105
0106   if (!allDTRecHits) {
0107     //Extract the DTRecSegment2DCollection range for the theta SL
0108     DTRecSegment2DCollection::range rangeTheta = all2DSegments->get(DTSuperLayerId(theChamber->id(), 2));
0109
0110     // Fill the DTRecSegment2D container for the theta SL
0111     vector<DTSLRecSegment2D> segments2DTheta(rangeTheta.first, rangeTheta.second);
0112
0113     if (debug)
0114       cout << "Number of 2D-segments in the second SL (Theta): " << segments2DTheta.size() << endl;
0115     theSegments2DTheta = segments2DTheta;
0116   }
0117 }
0118
0119 OwnVector<DTRecSegment4D> DTCombinatorialPatternReco4D::reconstruct() {
0120   OwnVector<DTRecSegment4D> result;
0121
0122   if (debug) {
0123     cout << "Segments in " << theChamber->id() << endl;
0124     cout << "Reconstructing of the Phi segments" << endl;
0125   }
0126
0127   vector<std::shared_ptr<DTHitPairForFit>> pairPhiOwned;
0128   vector<DTSegmentCand*> resultPhi = buildPhiSuperSegmentsCandidates(pairPhiOwned);
0129
0130   if (debug)
0131     cout << "There are " << resultPhi.size() << " Phi cand" << endl;
0132
0133   if (allDTRecHits) {
0134     // take the theta SL of this chamber
0135     const DTSuperLayer* sl = theChamber->superLayer(2);
0136     // sl points to 0 if the theta SL was not found
0137     if (sl) {
0138       // reconstruct the theta segments
0139       if (debug)
0140         cout << "Reconstructing of the Theta segments" << endl;
0141       OwnVector<DTSLRecSegment2D> thetaSegs = the2DAlgo->reconstruct(sl, theHitsFromTheta);
0142       vector<DTSLRecSegment2D> segments2DTheta(thetaSegs.begin(), thetaSegs.end());
0143       theSegments2DTheta = segments2DTheta;
0144     }
0145   }
0146
0147   bool hasZed = false;
0148
0149   // has this chamber the Z-superlayer?
0150   if (!theSegments2DTheta.empty()) {
0151     hasZed = !theSegments2DTheta.empty();
0152     if (debug)
0153       cout << "There are " << theSegments2DTheta.size() << " Theta cand" << endl;
0154   } else {
0155     if (debug)
0156       cout << "No Theta SL" << endl;
0157   }
0158
0159   // Now I want to build the concrete DTRecSegment4D.
0160   if (debug)
0161     cout << "Building of the concrete DTRecSegment4D" << endl;
0162   if (!resultPhi.empty()) {
0163     for (vector<DTSegmentCand*>::const_iterator phi = resultPhi.begin(); phi != resultPhi.end(); ++phi) {
0164       std::unique_ptr<DTChamberRecSegment2D> superPhi(**phi);
0165
0166       theUpdator->update(superPhi.get(), false);
0167       if (debug)
0168         cout << "superPhi: " << *superPhi << endl;
0169
0170       if (hasZed) {
0171         // Create all the 4D-segment combining the Z view with the Phi one
0172         // loop over the Z segments
0173         for (vector<DTSLRecSegment2D>::const_iterator zed = theSegments2DTheta.begin(); zed != theSegments2DTheta.end();
0174              ++zed) {
0175           if (debug)
0176             cout << "Theta: " << *zed << endl;
0177           // Important!!
0178           DTSuperLayerId ZedSegSLId(zed->geographicalId().rawId());
0179
0180           // Put the theta segment poistion in its 3D place.
0181           // note: (superPhi is in the CHAMBER local frame)
0182           const DTSuperLayer* zSL = theChamber->superLayer(ZedSegSLId);
0183
0184           // FIXME: should rather extrapolate for Y!
0185           LocalPoint zPos(
0186               zed->localPosition().x(), (zSL->toLocal(theChamber->toGlobal(superPhi->localPosition()))).y(), 0.);
0187
0188           const LocalPoint posZInCh = theChamber->toLocal(zSL->toGlobal(zPos));
0189           // FIXME: zed->localDirection() is in 2D. Should add the phi direction in the orthogonal plane as well!!
0190           const LocalVector dirZInCh = theChamber->toLocal(zSL->toGlobal(zed->localDirection()));
0191
0192           DTRecSegment4D* newSeg = new DTRecSegment4D(*superPhi, *zed, posZInCh, dirZInCh);
0193
0194           if (debug)
0195             cout << "Created a 4D seg " << *newSeg << endl;
0196
0197           /// 4d segment: I have the pos along the wire => further update!
0198           theUpdator->update(newSeg, false, false);
0199           if (debug)
0200             cout << "     seg updated " << *newSeg << endl;
0201
0202           if (!applyT0corr && computeT0corr)
0203             theUpdator->calculateT0corr(newSeg);
0204           if (applyT0corr)
0205             theUpdator->update(newSeg, true, false);
0206
0207           result.push_back(newSeg);
0208         }
0209       } else {
0210         // Only phi
0211
0212         DTRecSegment4D* newSeg = new DTRecSegment4D(*superPhi);
0213
0214         if (debug)
0215           cout << "Created a 4D segment using only the 2D Phi segment " << *newSeg << endl;
0216
0217         //update the segment with the t0 and possibly vdrift correction
0218         if (!applyT0corr && computeT0corr)
0219           theUpdator->calculateT0corr(newSeg);
0220         if (applyT0corr)
0221           theUpdator->update(newSeg, true, false);
0222
0223         result.push_back(newSeg);
0224       }
0225     }
0226   } else {
0227     // DTRecSegment4D from zed projection only (unlikely, not so useful, but...)
0228     if (hasZed) {
0229       for (vector<DTSLRecSegment2D>::const_iterator zed = theSegments2DTheta.begin(); zed != theSegments2DTheta.end();
0230            ++zed) {
0231         if (debug)
0232           cout << "Theta: " << *zed << endl;
0233
0234         // Important!!
0235         DTSuperLayerId ZedSegSLId(zed->geographicalId().rawId());
0236
0237         const LocalPoint posZInCh =
0238             theChamber->toLocal(theChamber->superLayer(ZedSegSLId)->toGlobal(zed->localPosition()));
0239         const LocalVector dirZInCh =
0240             theChamber->toLocal(theChamber->superLayer(ZedSegSLId)->toGlobal(zed->localDirection()));
0241
0242         DTRecSegment4D* newSeg = new DTRecSegment4D(*zed, posZInCh, dirZInCh);
0243
0244         if (debug)
0245           cout << "Created a 4D segment using only the 2D Theta segment " << *newSeg << endl;
0246
0247         if (!applyT0corr && computeT0corr)
0248           theUpdator->calculateT0corr(newSeg);
0249         if (applyT0corr)
0250           theUpdator->update(newSeg, true, false);
0251
0252         result.push_back(newSeg);
0253       }
0254     }
0255   }
0256   // finally delete the candidates!
0257   for (vector<DTSegmentCand*>::iterator phi = resultPhi.begin(); phi != resultPhi.end(); ++phi)
0258     delete *phi;
0259
0260   return result;
0261 }
0262
0263 vector<DTSegmentCand*> DTCombinatorialPatternReco4D::buildPhiSuperSegmentsCandidates(
0264     vector<std::shared_ptr<DTHitPairForFit>>& pairPhiOwned) {
0265   DTSuperLayerId slId;
0266
0267   if (!theHitsFromPhi1.empty())
0268     slId = theHitsFromPhi1.front().wireId().superlayerId();
0269   else if (!theHitsFromPhi2.empty())
0270     slId = theHitsFromPhi2.front().wireId().superlayerId();
0271   else {
0272     if (debug)
0273       cout << "DTCombinatorialPatternReco4D::buildPhiSuperSegmentsCandidates: "
0274            << "No Hits in the two Phi SL" << endl;
0275     return vector<DTSegmentCand*>();
0276   }
0277
0278   const DTSuperLayer* sl = theDTGeometry->superLayer(slId);
0279
0280   vector<std::shared_ptr<DTHitPairForFit>> pairPhi1 = the2DAlgo->initHits(sl, theHitsFromPhi1);
0281   // same sl!! Since the fit will be in the sl phi 1!
0282   vector<std::shared_ptr<DTHitPairForFit>> pairPhi2 = the2DAlgo->initHits(sl, theHitsFromPhi2);
0283   // copy the pairPhi2 in the pairPhi1 vector
0284   copy(pairPhi2.begin(), pairPhi2.end(), back_inserter(pairPhi1));
0285
0286   pairPhiOwned.swap(pairPhi1);
0287   // Build the segment candidate
0288   return the2DAlgo->buildSegments(sl, pairPhiOwned);
0289 }
0290
0291 /// Build a 4d segment with a zed component made by just one hits
0292 DTRecSegment4D* DTCombinatorialPatternReco4D::segmentSpecialZed(const DTRecSegment4D* seg) {
0293   // Get the zed projection
0294   //if (!seg->hasZed()) return seg;
0295   const DTSLRecSegment2D* zedSeg = seg->zSegment();
0296   std::vector<DTRecHit1D> hits = zedSeg->specificRecHits();
0297
0298   // pick up a hit "in the middle", where the single hit will be put.
0299   int nHits = hits.size();
0300   DTRecHit1D middle = hits[static_cast<int>(nHits / 2.)];
0301
0302   // Need to extrapolate pos to the middle layer z
0303   LocalPoint posInSL = zedSeg->localPosition();
0304   LocalVector dirInSL = zedSeg->localDirection();
0305   LocalPoint posInMiddleLayer = posInSL + dirInSL * (-posInSL.z()) / cos(dirInSL.theta());
0306
0307   // create a hit with position and error as the Zed projection one's
0308   auto hit = std::make_unique<DTRecHit1D>(
0309       middle.wireId(), middle.lrSide(), middle.digiTime(), posInMiddleLayer, zedSeg->localPositionError());
0310
0311   std::vector<DTRecHit1D> newHits(1, *hit);
0312
0313   // create a new zed segment with that single hits, but same dir and pos
0314   LocalPoint pos(zedSeg->localPosition());
0315   LocalVector dir(zedSeg->localDirection());
0316   AlgebraicSymMatrix cov(zedSeg->covMatrix());
0317   double chi2(zedSeg->chi2());
0318   //cout << "zed " << *zedSeg << endl;
0319   auto newZed = std::make_unique<DTSLRecSegment2D>(zedSeg->superLayerId(), pos, dir, cov, chi2, newHits);
0320   //cout << "newZed " << *newZed << endl;
0321
0322   // create a 4d segment with the special zed
0323   DTRecSegment4D* result = new DTRecSegment4D(*seg->phiSegment(), *newZed, seg->localPosition(), seg->localDirection());
0324   // delete the input segment
0325   delete seg;
0326
0327   // return it
0328   return result;
0329 }