CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
MuonCSCSeedFromRecHits.cc
Go to the documentation of this file.
12 #include <iomanip>
13 
14 
17 {
18 }
19 
20 
22 {
24  if(theRhits.size() == 1)
25  {
26  //return createSeed(100., 100., theRhits[0]);
27  makeDefaultSeed(result);
28  return result;
29  }
30  //@@ doesn't handle overlap between ME11 and ME12 correctly
31  // sort by station
32  MuonRecHitContainer station1Hits, station2Hits, station3Hits, station4Hits;
33  for ( MuonRecHitContainer::const_iterator iter = theRhits.begin(), end = theRhits.end();
34  iter != end; ++iter)
35  {
36  int station = CSCDetId((*iter)->geographicalId().rawId()).station();
37  if(station == 1)
38  {
39  station1Hits.push_back(*iter);
40  }
41  else if(station == 2)
42  {
43  station2Hits.push_back(*iter);
44  }
45  else if(station == 3)
46  {
47  station3Hits.push_back(*iter);
48  }
49  else if(station == 4)
50  {
51  station4Hits.push_back(*iter);
52  }
53  }
54 
55  //std::cout << "Station hits " << station1Hits.size() << " "
56  // << station2Hits.size() << " "
57  // << station3Hits.size() << std::endl;
58 
59  // see whether station 2 or station 3 is better
60  MuonRecHitContainer * betterSecondHits = &station2Hits;
61  MuonRecHitContainer * notAsGoodSecondHits = &station3Hits;
62  if(!station2Hits.empty() && !station3Hits.empty())
63  {
64  // swap if station 3 has better quailty
65  if(segmentQuality(station3Hits[0]) < segmentQuality(station2Hits[0]))
66  {
67  betterSecondHits = &station3Hits;
68  notAsGoodSecondHits = &station2Hits;
69  }
70  }
71 
72  // now try to make pairs
73  if(makeSeed(station1Hits, *betterSecondHits, result))
74  {
75  return result;
76  }
77  if(makeSeed(station1Hits, *notAsGoodSecondHits, result))
78  {
79  return result;
80  }
81  if(makeSeed(station2Hits, station3Hits, result))
82  {
83  return result;
84  }
85  if(makeSeed(station1Hits, station4Hits, result))
86  {
87  return result;
88  }
89 
90 
91 
92  // no luck
93  makeDefaultSeed(result);
94  return result;
95 }
96 
97 
99  TrajectorySeed & seed) const
100 {
101  for ( MuonRecHitContainer::const_iterator itr1 = hits1.begin(), end1 = hits1.end();
102  itr1 != end1; ++itr1)
103  {
104  CSCDetId cscId1((*itr1)->geographicalId().rawId());
105  //int type1 = CSCChamberSpecs::whatChamberType(cscId1.station(), cscId1.ring());
106 
107  for ( MuonRecHitContainer::const_iterator itr2 = hits2.begin(), end2 = hits2.end();
108  itr2 != end2; ++itr2)
109  {
110 
111  CSCDetId cscId2((*itr2)->geographicalId().rawId());
112  //int type2 = CSCChamberSpecs::whatChamberType(cscId2.station(), cscId2.ring());
113 
114  // take the first pair that comes along. Probably want to rank them later
115  std::vector<double> pts = thePtExtractor->pT_extract(*itr1, *itr2);
116 
117  double pt = pts[0];
118  double sigmapt = pts[1];
119  double minpt = 3.;
120 
121  // if too small, probably an error. Keep trying.
122  if(fabs(pt) > minpt)
123  {
124  double maxpt = 2000.;
125  if(pt > maxpt) {
126  pt = maxpt;
127  sigmapt = maxpt;
128  }
129  if(pt < -maxpt) {
130  pt = -maxpt;
131  sigmapt = maxpt;
132  }
133 
134  // get the position and direction from the higher-quality segment
136  seed = createSeed(pt, sigmapt, bestSeg);
137 
138  //std::cout << "FITTED TIMESPT " << pt << " dphi " << dphi << " eta " << eta << std::endl;
139  return true;
140  }
141 
142  }
143  }
144 
145  // guess it didn't find one
146 //std::cout << "NOTHING FOUND ! " << std::endl;
147  return false;
148 }
149 
150 
151 //typedef MuonTransientTrackingRecHit::MuonRecHitContainer MuonRecHitContainer;
153 {
154  int Nchi2 = 0;
155  int quality = 0;
156  int nhits = segment->recHits().size();
157  if ( segment->chi2()/(nhits*2.-4.) > 3. ) Nchi2 = 1;
158  if ( segment->chi2()/(nhits*2.-4.) > 9. ) Nchi2 = 2;
159 
160  if ( nhits > 4 ) quality = 1 + Nchi2;
161  if ( nhits == 4 ) quality = 3 + Nchi2;
162  if ( nhits == 3 ) quality = 5 + Nchi2;
163 
164  float dPhiGloDir = fabs ( deltaPhi(segment->globalPosition().phi(), segment->globalDirection().phi()) );
165 
166  if ( dPhiGloDir > .2 ) ++quality;
167  // add a penalty for being ME1A
168  if(segment->isCSC() && CSCDetId(segment->geographicalId()).ring() == 4) ++quality;
169  return quality;
170 }
171 
172 
173 
176 {
177  MuonRecHitPointer me1=0, meit=0;
178  float dPhiGloDir = .0; // +v
179  float bestdPhiGloDir = M_PI; // +v
180  int quality1 = 0, quality = 0; // +v I= 5,6-p. / II= 4p. / III= 3p.
181 
182  for ( MuonRecHitContainer::const_iterator iter = endcapHits.begin(); iter!= endcapHits.end(); iter++ ){
183  if ( !(*iter)->isCSC() ) continue;
184 
185  // tmp compar. Glob-Dir for the same tr-segm:
186 
187  meit = *iter;
188 
189  quality = segmentQuality(meit);
190 
191  dPhiGloDir = fabs ( deltaPhi(meit->globalPosition().phi(), meit->globalDirection().phi()) );
192 
193  if(!me1){
194  me1 = meit;
195  quality1 = quality;
196  bestdPhiGloDir = dPhiGloDir;
197  }
198 
199  if(me1) {
200 
201  if ( !me1->isValid() ) {
202  me1 = meit;
203  quality1 = quality;
204  bestdPhiGloDir = dPhiGloDir;
205  }
206 
207  if ( me1->isValid() && quality < quality1 ) {
208  me1 = meit;
209  quality1 = quality;
210  bestdPhiGloDir = dPhiGloDir;
211  }
212 
213  if ( me1->isValid() && bestdPhiGloDir > .03 ) {
214  if ( dPhiGloDir < bestdPhiGloDir - .01 && quality == quality1 ) {
215  me1 = meit;
216  quality1 = quality;
217  bestdPhiGloDir = dPhiGloDir;
218  }
219  }
220  }
221 
222  } // iter
223 
224  return me1;
225 }
226 
227 
229 {
230  //Search ME1 ...
232  bool good=false;
233 
234  if(me1 && me1->isValid() )
235  {
236  good = createDefaultEndcapSeed(me1, seed);
237  }
238 }
239 
240 
241 
242 
243 
244 bool
246  TrajectorySeed & seed) const {
247  //float momentum = computeDefaultPt(last);
248  std::vector<double> momentum = thePtExtractor->pT_extract(last, last);
249  seed = createSeed(momentum[0],momentum[1],last);
250  return true;
251 }
252 
253 
254 
256 {
257  for ( MuonRecHitContainer::const_iterator iter = theRhits.begin(); iter!= theRhits.end(); iter++ ){
258  if ( !(*iter)->isCSC() ) continue;
259  for ( MuonRecHitContainer::const_iterator iter2 = iter + 1; iter2 != theRhits.end(); ++iter2)
260  {
261  if( !(*iter2)->isCSC() ) continue;
262 
263  CSCDetId cscId1((*iter)->geographicalId().rawId());
264  CSCDetId cscId2((*iter2)->geographicalId().rawId());
265  double dphi = deltaPhi((**iter).globalPosition().phi(), (**iter2).globalPosition().phi());
266 
267  int type1 = cscId1.iChamberType();
268  int type2 = cscId2.iChamberType();
269 
270  // say the lower station first
271  if(type1 < type2)
272  {
273  std::cout << "HITPAIRA," << type1 << type2 << "," <<
274  dphi << "," << (**iter).globalPosition().eta() << std::endl;
275  }
276  if(type2 < type1)
277  {
278  std::cout << "HITPAIRB," << type2 << type1 << "," <<
279  -dphi << "," << (**iter2).globalPosition().eta() << std::endl;
280  }
281  if(type1 == type2)
282  {
283  std::cout << "HITPAIRSAMESTATION," << type1 << cscId1.ring() << cscId2.ring()
284  << "," << dphi << "," << (**iter).globalPosition().eta() << std::endl;
285  }
286 
287  }
288  }
289 }
virtual TrajectorySeed seed() const
double deltaPhi(float phi1, float phi2)
Definition: VectorUtil.h:30
MuonTransientTrackingRecHit::MuonRecHitContainer theRhits
MuonTransientTrackingRecHit::MuonRecHitContainer MuonRecHitContainer
const MuonSeedPtExtractor * thePtExtractor
ConstMuonRecHitPointer bestEndcapHit(const MuonRecHitContainer &endcapHits) const
bool createDefaultEndcapSeed(ConstMuonRecHitPointer last, TrajectorySeed &seed) const
tuple result
Definition: query.py:137
int segmentQuality(ConstMuonRecHitPointer muon) const
#define end
Definition: vmac.h:38
void makeDefaultSeed(TrajectorySeed &seed) const
#define M_PI
Definition: BFit3D.cc:3
TrajectorySeed createSeed(float ptmean, float sptmean, MuonTransientTrackingRecHit::ConstMuonRecHitPointer last) const
tuple cout
Definition: gather_cfg.py:41
virtual std::vector< double > pT_extract(MuonTransientTrackingRecHit::ConstMuonRecHitPointer firstHit, MuonTransientTrackingRecHit::ConstMuonRecHitPointer secondHit) const
bool makeSeed(const MuonRecHitContainer &hits1, const MuonRecHitContainer &hits2, TrajectorySeed &seed) const