DTSegment/src/DTRefitAndCombineReco4D.cc

0001 /** \class DTRefitAndCombineReco4D
0002  *
0003  * Algo for reconstructing 4d segment in DT refitting the 2D phi SL hits and combining
0004  * the results with the theta view.
0005  *
0006  * \author Stefano Lacaprara - INFN Legnaro <stefano.lacaprara@pd.infn.it>
0007  * \author Riccardo Bellan - INFN TO <riccardo.bellan@cern.ch>
0008  *
0009  */
0010
0011 #include "RecoLocalMuon/DTSegment/src/DTRefitAndCombineReco4D.h"
0012
0013 #include "RecoLocalMuon/DTSegment/src/DTSegmentUpdator.h"
0014 #include "RecoLocalMuon/DTSegment/src/DTSegmentCand.h"
0015
0016 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0017 #include "FWCore/Framework/interface/EventSetup.h"
0018 #include "FWCore/Framework/interface/ConsumesCollector.h"
0019 #include "FWCore/Utilities/interface/Exception.h"
0020
0021 #include "DataFormats/Common/interface/OwnVector.h"
0022 #include "DataFormats/DTRecHit/interface/DTChamberRecSegment2D.h"
0023 #include "DataFormats/DTRecHit/interface/DTRecSegment2DCollection.h"
0024
0025 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0026
0027 using namespace std;
0028 using namespace edm;
0029
0030 DTRefitAndCombineReco4D::DTRefitAndCombineReco4D(const ParameterSet& pset, ConsumesCollector cc)
0031     : DTRecSegment4DBaseAlgo(pset), theAlgoName("DTRefitAndCombineReco4D"), theDTGeometryToken(cc.esConsumes()) {
0032   // debug parameter
0033   debug = pset.getUntrackedParameter<bool>("debug");
0034
0035   // the updator
0036   theUpdator = new DTSegmentUpdator(pset, cc);
0037
0038   // the max allowd chi^2 for the fit of th combination of two phi segments
0039   theMaxChi2forPhi = pset.getParameter<double>("MaxChi2forPhi");
0040 }
0041
0042 void DTRefitAndCombineReco4D::setES(const EventSetup& setup) {
0043   theDTGeometry = setup.getHandle(theDTGeometryToken);
0044   theUpdator->setES(setup);
0045   //  the2DAlgo->setES(setup);
0046 }
0047
0048 void DTRefitAndCombineReco4D::setChamber(const DTChamberId& chId) {
0049   // Set the chamber
0050   theChamber = theDTGeometry->chamber(chId);
0051 }
0052
0053 void DTRefitAndCombineReco4D::setDTRecSegment2DContainer(Handle<DTRecSegment2DCollection> allHits) {
0054   theSegments2DPhi1.clear();
0055   theSegments2DTheta.clear();
0056   theSegments2DPhi2.clear();
0057
0058   // Get the chamber
0059   //    const DTChamber *chamber = theDTGeometry->chamber(chId);
0060
0061   const DTChamberId chId = theChamber->id();
0062
0063   //Extract the DTRecSegment2DCollection ranges for the three different SL
0064   DTRecSegment2DCollection::range rangePhi1 = allHits->get(DTSuperLayerId(chId, 1));
0065   DTRecSegment2DCollection::range rangeTheta = allHits->get(DTSuperLayerId(chId, 2));
0066   DTRecSegment2DCollection::range rangePhi2 = allHits->get(DTSuperLayerId(chId, 3));
0067
0068   // Fill the DTSLRecSegment2D containers for the three different SL
0069   vector<DTSLRecSegment2D> segments2DPhi1(rangePhi1.first, rangePhi1.second);
0070   vector<DTSLRecSegment2D> segments2DTheta(rangeTheta.first, rangeTheta.second);
0071   vector<DTSLRecSegment2D> segments2DPhi2(rangePhi2.first, rangePhi2.second);
0072
0073   if (debug)
0074     cout << "Number of 2D-segments in the first  SL (Phi)" << segments2DPhi1.size() << endl
0075          << "Number of 2D-segments in the second SL (Theta)" << segments2DTheta.size() << endl
0076          << "Number of 2D-segments in the third  SL (Phi)" << segments2DPhi2.size() << endl;
0077
0078   theSegments2DPhi1 = segments2DPhi1;
0079   theSegments2DTheta = segments2DTheta;
0080   theSegments2DPhi2 = segments2DPhi2;
0081 }
0082
0083 OwnVector<DTRecSegment4D> DTRefitAndCombineReco4D::reconstruct() {
0084   OwnVector<DTRecSegment4D> result;
0085
0086   if (debug)
0087     cout << "Segments in " << theChamber->id() << endl;
0088
0089   vector<DTChamberRecSegment2D> resultPhi = refitSuperSegments();
0090
0091   if (debug)
0092     cout << "There are " << resultPhi.size() << " Phi cand" << endl;
0093
0094   bool hasZed = false;
0095
0096   // has this chamber the Z-superlayer?
0097   if (!theSegments2DTheta.empty()) {
0098     hasZed = !theSegments2DTheta.empty();
0099     if (debug)
0100       cout << "There are " << theSegments2DTheta.size() << " Theta cand" << endl;
0101   } else {
0102     if (debug)
0103       cout << "No Theta SL" << endl;
0104   }
0105
0106   // Now I want to build the concrete DTRecSegment4D.
0107   if (!resultPhi.empty()) {
0108     for (vector<DTChamberRecSegment2D>::const_iterator phi = resultPhi.begin(); phi != resultPhi.end(); ++phi) {
0109       if (hasZed) {
0110         // Create all the 4D-segment combining the Z view with the Phi one
0111         // loop over the Z segments
0112         for (vector<DTSLRecSegment2D>::const_iterator zed = theSegments2DTheta.begin(); zed != theSegments2DTheta.end();
0113              ++zed) {
0114           //>> Important!!
0115           DTSuperLayerId ZedSegSLId(zed->geographicalId().rawId());
0116
0117           const LocalPoint posZInCh =
0118               theChamber->toLocal(theChamber->superLayer(ZedSegSLId)->toGlobal(zed->localPosition()));
0119           const LocalVector dirZInCh =
0120               theChamber->toLocal(theChamber->superLayer(ZedSegSLId)->toGlobal(zed->localDirection()));
0121
0122           DTRecSegment4D* newSeg = new DTRecSegment4D(*phi, *zed, posZInCh, dirZInCh);
0123           //<<
0124
0125           /// 4d segment: I have the pos along the wire => further update!
0126           theUpdator->update(newSeg, false, true);
0127           if (debug)
0128             cout << "Created a 4D seg " << endl;
0129           result.push_back(newSeg);
0130         }
0131       } else {
0132         // Only phi
0133         DTRecSegment4D* newSeg = new DTRecSegment4D(*phi);
0134         if (debug)
0135           cout << "Created a 4D segment using only the 2D Phi segment" << endl;
0136         result.push_back(newSeg);
0137       }
0138     }
0139   } else {
0140     // DTRecSegment4D from zed projection only (unlikely, not so useful, but...)
0141     if (hasZed) {
0142       for (vector<DTSLRecSegment2D>::const_iterator zed = theSegments2DTheta.begin(); zed != theSegments2DTheta.end();
0143            ++zed) {
0144         // Important!!
0145         DTSuperLayerId ZedSegSLId(zed->geographicalId().rawId());
0146
0147         const LocalPoint posZInCh =
0148             theChamber->toLocal(theChamber->superLayer(ZedSegSLId)->toGlobal(zed->localPosition()));
0149         const LocalVector dirZInCh =
0150             theChamber->toLocal(theChamber->superLayer(ZedSegSLId)->toGlobal(zed->localDirection()));
0151
0152         DTRecSegment4D* newSeg = new DTRecSegment4D(*zed, posZInCh, dirZInCh);
0153         //<<
0154
0155         if (debug)
0156           cout << "Created a 4D segment using only the 2D Theta segment" << endl;
0157         result.push_back(newSeg);
0158       }
0159     }
0160   }
0161
0162   return result;
0163 }
0164
0165 vector<DTChamberRecSegment2D> DTRefitAndCombineReco4D::refitSuperSegments() {
0166   vector<DTChamberRecSegment2D> result;
0167
0168   //double-loop over all the DTSLRecSegment2D in order to make all the possible pairs
0169   for (vector<DTSLRecSegment2D>::const_iterator segment2DPhi1 = theSegments2DPhi1.begin();
0170        segment2DPhi1 != theSegments2DPhi1.end();
0171        ++segment2DPhi1) {
0172     for (vector<DTSLRecSegment2D>::const_iterator segment2DPhi2 = theSegments2DPhi2.begin();
0173          segment2DPhi2 != theSegments2DPhi2.end();
0174          ++segment2DPhi2) {
0175       // check the id
0176       if (segment2DPhi1->chamberId() != segment2DPhi2->chamberId())
0177         throw cms::Exception("refitSuperSegments") << "he phi segments have different chamber id" << std::endl;
0178
0179       // create a super phi starting from 2 phi
0180       vector<DTRecHit1D> recHitsSeg2DPhi1 = segment2DPhi1->specificRecHits();
0181       vector<DTRecHit1D> recHitsSeg2DPhi2 = segment2DPhi2->specificRecHits();
0182       // copy the recHitsSeg2DPhi2 in the recHitsSeg2DPhi1 container
0183       copy(recHitsSeg2DPhi2.begin(), recHitsSeg2DPhi2.end(), back_inserter(recHitsSeg2DPhi1));
0184
0185       const DTChamberId chId = segment2DPhi1->chamberId();
0186
0187       // create the super phi
0188       DTChamberRecSegment2D superPhi(chId, recHitsSeg2DPhi1);
0189
0190       // refit it!
0191       theUpdator->fit(&superPhi, false, false);
0192
0193       // cut on the chi^2
0194       if (superPhi.chi2() > theMaxChi2forPhi)
0195         result.push_back(superPhi);
0196     }
0197   }
0198   // TODO clean the container!!!
0199   // there are some possible repetition!
0200   // maybe using the cleaner, previous a conversion from DTChamberRecSegment2D to DTSegmentCandidate
0201   return result;
0202 }