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