00001
00013 #include "RecoLocalMuon/DTSegment/src/DTRefitAndCombineReco4D.h"
00014
00015 #include "RecoLocalMuon/DTSegment/src/DTSegmentUpdator.h"
00016 #include "RecoLocalMuon/DTSegment/src/DTSegmentCand.h"
00017
00018 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00019 #include "FWCore/Framework/interface/EventSetup.h"
00020 #include "FWCore/Utilities/interface/Exception.h"
00021
00022 #include "DataFormats/Common/interface/OwnVector.h"
00023 #include "DataFormats/DTRecHit/interface/DTChamberRecSegment2D.h"
00024 #include "DataFormats/DTRecHit/interface/DTRecSegment2DCollection.h"
00025
00026 #include "Geometry/Records/interface/MuonGeometryRecord.h"
00027
00028 using namespace std;
00029 using namespace edm;
00030
00031 DTRefitAndCombineReco4D::DTRefitAndCombineReco4D(const ParameterSet& pset):
00032 DTRecSegment4DBaseAlgo(pset), theAlgoName("DTRefitAndCombineReco4D"){
00033
00034
00035 debug = pset.getUntrackedParameter<bool>("debug");
00036
00037
00038 theUpdator = new DTSegmentUpdator(pset);
00039
00040
00041 theMaxChi2forPhi = pset.getParameter<double>("MaxChi2forPhi");
00042 }
00043
00044 void DTRefitAndCombineReco4D::setES(const EventSetup& setup){
00045 setup.get<MuonGeometryRecord>().get(theDTGeometry);
00046 theUpdator->setES(setup);
00047
00048 }
00049
00050 void DTRefitAndCombineReco4D::setChamber(const DTChamberId &chId){
00051
00052 theChamber = theDTGeometry->chamber(chId);
00053 }
00054
00055 void
00056 DTRefitAndCombineReco4D::setDTRecSegment2DContainer(Handle<DTRecSegment2DCollection> allHits) {
00057 theSegments2DPhi1.clear();
00058 theSegments2DTheta.clear();
00059 theSegments2DPhi2.clear();
00060
00061
00062
00063
00064 const DTChamberId chId = theChamber->id();
00065
00066
00067 DTRecSegment2DCollection::range rangePhi1 = allHits->get(DTSuperLayerId(chId,1));
00068 DTRecSegment2DCollection::range rangeTheta = allHits->get(DTSuperLayerId(chId,2));
00069 DTRecSegment2DCollection::range rangePhi2 = allHits->get(DTSuperLayerId(chId,3));
00070
00071
00072 vector<DTSLRecSegment2D> segments2DPhi1(rangePhi1.first,rangePhi1.second);
00073 vector<DTSLRecSegment2D> segments2DTheta(rangeTheta.first,rangeTheta.second);
00074 vector<DTSLRecSegment2D> segments2DPhi2(rangePhi2.first,rangePhi2.second);
00075
00076 if(debug)
00077 cout << "Number of 2D-segments in the first SL (Phi)" << segments2DPhi1.size() << endl
00078 << "Number of 2D-segments in the second SL (Theta)" << segments2DTheta.size() << endl
00079 << "Number of 2D-segments in the third SL (Phi)" << segments2DPhi2.size() << endl;
00080
00081 theSegments2DPhi1 = segments2DPhi1;
00082 theSegments2DTheta = segments2DTheta;
00083 theSegments2DPhi2 = segments2DPhi2;
00084 }
00085
00086
00087
00088 OwnVector<DTRecSegment4D>
00089 DTRefitAndCombineReco4D::reconstruct(){
00090 OwnVector<DTRecSegment4D> result;
00091
00092 if (debug) cout << "Segments in " << theChamber->id() << endl;
00093
00094 vector<DTChamberRecSegment2D> resultPhi = refitSuperSegments();
00095
00096 if (debug) cout << "There are " << resultPhi.size() << " Phi cand" << endl;
00097
00098 bool hasZed=false;
00099
00100
00101 if (theSegments2DTheta.size()){
00102 hasZed = theSegments2DTheta.size()>0;
00103 if (debug) cout << "There are " << theSegments2DTheta.size() << " Theta cand" << endl;
00104 } else {
00105 if (debug) cout << "No Theta SL" << endl;
00106 }
00107
00108
00109 if (resultPhi.size()) {
00110 for (vector<DTChamberRecSegment2D>::const_iterator phi=resultPhi.begin();
00111 phi!=resultPhi.end(); ++phi) {
00112
00113 if (hasZed) {
00114
00115
00116
00117 for(vector<DTSLRecSegment2D>::const_iterator zed = theSegments2DTheta.begin();
00118 zed != theSegments2DTheta.end(); ++zed){
00119
00120
00121 DTSuperLayerId ZedSegSLId(zed->geographicalId().rawId());
00122
00123 const LocalPoint posZInCh = theChamber->toLocal( theChamber->superLayer(ZedSegSLId)->toGlobal(zed->localPosition() )) ;
00124 const LocalVector dirZInCh = theChamber->toLocal( theChamber->superLayer(ZedSegSLId)->toGlobal(zed->localDirection() )) ;
00125
00126 DTRecSegment4D* newSeg = new DTRecSegment4D(*phi,*zed,posZInCh,dirZInCh);
00127
00128
00130 theUpdator->update(newSeg);
00131 if (debug) cout << "Created a 4D seg " << endl;
00132 result.push_back(newSeg);
00133 }
00134 } else {
00135
00136 DTRecSegment4D* newSeg = new DTRecSegment4D(*phi);
00137 if (debug) cout << "Created a 4D segment using only the 2D Phi segment" << endl;
00138 result.push_back(newSeg);
00139 }
00140 }
00141 } else {
00142
00143 if (hasZed) {
00144 for(vector<DTSLRecSegment2D>::const_iterator zed = theSegments2DTheta.begin();
00145 zed != theSegments2DTheta.end(); ++zed){
00146
00147
00148 DTSuperLayerId ZedSegSLId(zed->geographicalId().rawId());
00149
00150 const LocalPoint posZInCh = theChamber->toLocal( theChamber->superLayer(ZedSegSLId)->toGlobal(zed->localPosition() )) ;
00151 const LocalVector dirZInCh = theChamber->toLocal( theChamber->superLayer(ZedSegSLId)->toGlobal(zed->localDirection() )) ;
00152
00153 DTRecSegment4D* newSeg = new DTRecSegment4D( *zed,posZInCh,dirZInCh);
00154
00155
00156 if (debug) cout << "Created a 4D segment using only the 2D Theta segment" << endl;
00157 result.push_back(newSeg);
00158 }
00159 }
00160 }
00161
00162 return result;
00163 }
00164
00165 vector<DTChamberRecSegment2D> DTRefitAndCombineReco4D::refitSuperSegments(){
00166 vector<DTChamberRecSegment2D> result;
00167
00168
00169 for(vector<DTSLRecSegment2D>::const_iterator segment2DPhi1 = theSegments2DPhi1.begin();
00170 segment2DPhi1 != theSegments2DPhi1.end(); ++segment2DPhi1){
00171 for(vector<DTSLRecSegment2D>::const_iterator segment2DPhi2 = theSegments2DPhi2.begin();
00172 segment2DPhi2 != theSegments2DPhi2.end(); ++segment2DPhi2){
00173
00174
00175 if(segment2DPhi1->chamberId() != segment2DPhi2->chamberId())
00176 throw cms::Exception("refitSuperSegments")
00177 <<"he phi segments have different chamber id"<<std::endl;
00178
00179
00180 vector<DTRecHit1D> recHitsSeg2DPhi1 = segment2DPhi1->specificRecHits();
00181 vector<DTRecHit1D> recHitsSeg2DPhi2 = segment2DPhi2->specificRecHits();
00182
00183 copy(recHitsSeg2DPhi2.begin(),recHitsSeg2DPhi2.end(),back_inserter(recHitsSeg2DPhi1));
00184
00185 const DTChamberId chId = segment2DPhi1->chamberId();
00186
00187
00188 DTChamberRecSegment2D superPhi(chId,recHitsSeg2DPhi1);
00189
00190
00191 theUpdator->fit(&superPhi);
00192
00193
00194 if (superPhi.chi2() > theMaxChi2forPhi)
00195 result.push_back(superPhi);
00196 }
00197 }
00198
00199
00200
00201 return result;
00202 }