CGEM BOSS 6.6.5.f
BESIII Offline Software System
Loading...
Searching...
No Matches
FTTrack.cxx
Go to the documentation of this file.
1// -*- C++ -*-
2//
3// Package: MdcFastTrkAlg
4// Module: FTTrack
5//
6// Description: track class for MdcFastTrkAlg
7//
8//
9
10#include <stdlib.h>
11#include "MdcFastTrkAlg/FTTrack.h"
12#include "MdcFastTrkAlg/FTSegment.h"
13#include "MdcFastTrkAlg/FTFinder.h"
14#include "MdcFastTrkAlg/FTWire.h"
15#include "MdcFastTrkAlg/MdcParameter.h"
16
17#ifndef OnlineMode
18#include "AIDA/IAxis.h"
19#include "AIDA/IHistogram1D.h"
20#include "AIDA/IHistogram2D.h"
21#include "AIDA/IHistogram3D.h"
22#include "AIDA/IHistogramFactory.h"
23
24extern IHistogram1D* g_sigmaxy;
25extern IHistogram1D* g_sigmaz;
26extern IHistogram1D* g_chi2xy;
27extern IHistogram1D* g_chi2sz;
28
29#endif
30
31//const HepPoint3D
32//FTTrack::m_ftFinder->pivot(0,0,0);
33
35
36int
38{
39 IMessageSvc *msgSvc;
40 Gaudi::svcLocator()->service("MessageSvc", msgSvc);
41
42 MsgStream log(msgSvc, "FTFinder");
43
44 // static const float alpha(333.564095);
45 if (fabs(_kappa) > 1.2/param->_minPt) return 0;
46
47 _la = new Lpav;
48 int n = _axial_segments.length();
49 for (int i = 0; i^n; i++){
50 FTList<FTWire *> & hits = _axial_segments[i]->wireHits();
51 int m = hits.length();
52 for (int j = 0; j^m; j++){
53 FTWire & h = * hits[j];
55 hits.remove(j);
56 continue;
57 }
58 // log << MSG::DEBUG << "layer: "<< h.layer().layerId() << " phi: " << h.phi() << " df_distance: "<< h.distance() << endreq;
59 double par = h.distance()/(0.25*h.layer().csize());
60 _la->add_point((double)h.x(),(double)h.y(),exp(-par*par));
61 }
62 }
63
64 double chi2 = _la->fit();
65 HepVector a = _la->Hpar(m_ftFinder->pivot);
66 log << MSG::DEBUG << " chi2/_la cut(1): " << chi2 << " / " << _la->nc()
67 <<" a1(" << param->_minDr << "): " << a(1)
68 <<" a3(" << 1.05/param->_minPt << "):" << a(3) <<endreq;
69 if (chi2/_la->nc() > 1.) return 0;
70 if (fabs(a(3))>(1.05/param->_minPt)) return 0;
71 if (fabs(a(1))>param->_minDr) return 0;
72 _la->clear();
73 log << MSG::DEBUG << " passed chi2/_la, a(3) and a(1) cut" <<endreq;
74 int layer0 = 0;
75 for (int i = 0; i^n; i++){
76 FTList<FTWire *> & hits = _axial_segments[i]->wireHits();
77 if (!_axial_segments[i]->superLayer().superLayerId()) layer0 = 1;
78 int m = hits.length();
79 for (int j = 0; j^m; j++){
80 FTWire & h = * hits[j];
81 const float x = h.x();
82 const float y = h.y();
83 double d0 = _la->d((double)x,(double)y);
84 float delta = h.distance()/h.layer().r();
85 if (fabs(d0) > 0.7*h.layer().csize()) continue; // remove bad hits
86 if (d0>0){ // left or right
87 d0 = -d0;
88 }else{
89 delta = -delta;
90 }
91 _la->add_point(x - delta*y, y + delta*x, 1);
92 }
93 }
94 if (!layer0){ // salvage hits from complecated segments in layer0
95 FTList<FTSegment *> & salvage =
96 m_ftFinder->superLayer(0)->complecated_segments();
97 HepVector center = _la->center();
98 const float xc = center(1);
99 const float yc = center(2);
100 const float rc = a(1)+(-1. / 2.99792458 /m_pmgnIMF->getReferField())/a(3); // rho+drho(signed)
101 int nn = salvage.length();
102 for (int i = 0; i^nn; i++){
103 int salvaged = 0;
104 FTList<FTWire *> & hits = salvage[i]->wireHits();
105 int m = hits.length();
106 for (int j = 0; j^m; j++){
107 FTWire & h = * hits[j];
108 float x = h.x();
109 float y = h.y();
110 float r = h.layer().r();
111 if ((y*xc-x*yc)/(r*rc)<0.707) break;
112 double d0 = _la->d((double)x,(double)y);
113 float delta = h.distance()/r;
114 if (fabs(d0) > 0.7*h.layer().csize()) continue; // remove bad hits
115 if (d0>0){ // left or right
117 d0 = -d0;
118 }else{
120 delta = -delta;
121 }
122 _la->add_point(x - delta*y, y + delta*x, 1);
123 salvaged = 1;
124 }
125 if (salvaged){
126 _axial_segments.append(salvage[i]);
127 break;
128 }
129 }
130 }
131 chi2 = _la->fit(); // refit using drift distance
135 _za = new zav;
136 _SigmaS = 0.;
137 _SigmaZ = 0.;
138 _SigmaSS = 0.;
139 _SigmaSZ = 0.;
140 return 1;
141}
142
143int
144FTTrack::r_phiReFit(float vx, float vy, int vtx_flag)
145{
146 if (vtx_flag) _la->fit(vx, vy, 20*_la->chisq()/_la->nc());
147 int n = _axial_segments.length();
148 _la->clear();
149 for (int i = 0; i^n; i++){
150 FTList<FTWire *> & hits = _axial_segments[i]->wireHits();
151 int m = hits.length();
152 for (int j = 0; j^m; j++){
153 FTWire & h = * hits[j];
154 const float x = h.x();
155 const float y = h.y();
156 double d0 = _la->d((double)x,(double)y);
157 float cellsize = h.layer().csize();
158 if (m_ftFinder->evtTiming){
159 float time = h.time() + m_ftFinder->evtTiming;
160 if(time<-100) continue; // discard the bad TDC time
161 h.distance(m_ftFinder->t2x(h.layer(),time));
162 //h.distance(time*40*0.0001);
163 }
164 float delta = h.distance()/h.layer().r();
165 if (fabs(d0) > 0.5*cellsize) continue; // remove bad hits
166 if (d0>0){ // left or right
168 d0 = -d0;
169 }else{
171 delta = -delta;
172 }
173 h.setChi2(h.distance()+d0);
174#ifndef OnlineMode
175 //fill histogram
176 g_sigmaxy->fill(h.distance()+d0, 1.0);
177#endif
178 _la->add_point(x - delta*y, y + delta*x, 1);
179 }
180 }
181 HepVector b = _la->Hpar(m_ftFinder->pivot);
182 double chi2 = _la->fit(); // refit using drift distance
183#ifndef OnlineMode
184 g_chi2xy->fill(chi2/_la->nc(), 1.0);
185#endif
186
187 return 1;
188}
189
190int
191FTTrack::r_phi2Fit(float vx, float vy, int vtx_flag)
192{
193 if (vtx_flag) _la->fit(vx, vy, 20*_la->chisq()/_la->nc());
194 int n = _axial_segments.length();
195 _la->clear();
196 int k=0,l=0;
197 float temp=0;
198 for (int i = 0; i^n; i++){
199 FTList<FTWire *> & hits = _axial_segments[i]->wireHits();
200 int m = hits.length();
201 for (int j = 0; j^m; j++){
202 k++;
203 FTWire & h = * hits[j];
204 const float x = h.x();
205 const float y = h.y();
206 double d0 = _la->d((double)x,(double)y);
207 float cellsize = h.layer().csize();
208 if (m_ftFinder->evtTiming){
209 float time = h.time() + m_ftFinder->evtTiming ;
210 if(time<-100) continue;
211 h.distance(m_ftFinder->t2x(h.layer(),time));
212 // h.distance(time*40*0.0001);
213 }
214 float delta = h.distance()/h.layer().r();
215 if (fabs(d0) > 0.5*cellsize) continue; // remove bad hits
216 if (h.distance()<0.1 || h.distance()>0.6) continue;
217 if (d0>0){ // left or right
219 d0 = -d0;
220 }else{
222 delta = -delta;
223 }
224#ifndef OnlineMode
225 //fill histogram
226 //g_sigmaxy->fill(h.distance()+d0, 1.0);
227#endif
228 _la->add_point(x - delta*y, y + delta*x, 1);
229 if(temp<fabs(h.distance()+d0)){
230 temp=fabs(h.distance()+d0);
231 l=k;
232 }
233 }
234 }
235 HepVector b = _la->Hpar(m_ftFinder->pivot);
236 double chi2 = _la->fit(); // refit using drift distance
237#ifndef OnlineMode
238 // g_chi2xy->fill(chi2/_la->nc(), 1.0);
239#endif
240 if(k>5){
241 return l;
242 }
243 else{
244 return -99;
245 }
246}
247
248int
249FTTrack::r_phi3Fit(int l, float vx, float vy, int vtx_flag)
250{
251 if (vtx_flag) _la->fit(vx, vy, 20*_la->chisq()/_la->nc());
252 int n = _axial_segments.length();
253 _la->clear();
254 int k=0;
255 for (int i = 0; i^n; i++){
256 FTList<FTWire *> & hits = _axial_segments[i]->wireHits();
257 int m = hits.length();
258 for (int j = 0; j^m; j++){
259 k++;
260 FTWire & h = * hits[j];
261 const float x = h.x();
262 const float y = h.y();
263 double d0 = _la->d((double)x,(double)y);
264 float cellsize = h.layer().csize();
265 if (m_ftFinder->evtTiming){
266 float time = h.time() + m_ftFinder->evtTiming ;
267 if(time<-100) continue;
268 h.distance(m_ftFinder->t2x(h.layer(),time));
269 //h.distance(time*40*0.0001);
270 }
271 float delta = h.distance()/h.layer().r();
272 if (fabs(d0) > 0.5*cellsize) continue; // remove bad hits
273 if (d0>0){ // left or right
275 d0 = -d0;
276 }else{
278 delta = -delta;
279 }
280 if (k==l){
281 //delete hits[j];
282 m=hits.remove(j);
283 continue;
284 }
285 _la->add_point(x - delta*y, y + delta*x, 1);
286 // h.setChi2(h.distance()+d0);
287 }
288 }
289 HepVector b = _la->Hpar(m_ftFinder->pivot);
290 double chi2 = _la->fit(); // refit using drift distance
291 return 1;
292}
293
294int
295FTTrack::r_phi4Fit(float vx, float vy, int vtx_flag)
296{
297 if (vtx_flag) _la->fit(vx, vy, 20*_la->chisq()/_la->nc());
298 int n = _axial_segments.length();
299 _la->clear();
300 for (int i = 0; i^n; i++){
301 FTList<FTWire *> & hits = _axial_segments[i]->wireHits();
302 int m = hits.length();
303 for (int j = 0; j^m; j++){
304 FTWire & h = * hits[j];
305 const float x = h.x();
306 const float y = h.y();
307 double d0 = _la->d((double)x,(double)y);
308 float cellsize = h.layer().csize();
309 if (m_ftFinder->evtTiming){
310 float time = h.time() + m_ftFinder->evtTiming ;
311 if(time<-100) continue;
312 h.distance(m_ftFinder->t2x(h.layer(),time));
313 //h.distance(time*40*0.0001);
314 }
315 float delta = h.distance()/h.layer().r();
316 if (fabs(d0) > 0.5*cellsize) continue; // remove bad hits
317 if (d0>0){ // left or right
319 d0 = -d0;
320 }else{
322 delta = -delta;
323 }
324
325 _la->add_point(x - delta*y, y + delta*x, 1);
326 h.setChi2(h.distance()+d0);
327#ifndef OnlineMode
328 //fill histogram
329 // g_sigmaxy->fill(h.distance()+d0, 1.0);
330#endif
331 }
332 }
333 HepVector b = _la->Hpar(m_ftFinder->pivot);
334 double chi2 = _la->fit(); // refit using drift distance
335#ifndef OnlineMode
336 // g_chi2xy->fill(chi2/_la->nc(), 1.0);
337#endif
338 return 1;
339}
340
341int
345 for (int i = 0; i^n; i++){
346 FTList<FTSegment *> & segments = *(*_stereo_segments_by_superLayer)[i];
347 int m = segments.length();
348 float min_D_z = 9998.;
349 float S = 0.;
350 float Z = 0.;
351 FTSegment * selected = NULL;
352 for (int j = 0; j^m; j++){
353 FTSegment * s = segments[j];
354 float s_tmp = s->s();
355 float z_tmp = s->z();
356 double D_z = fabs(d_z(s_tmp,z_tmp));
357 if (D_z < min_D_z){
358 selected = s;
359 min_D_z = D_z;
360 S = s_tmp;
361 Z = z_tmp;
362 }
363 }
364 if (selected){
366 _SigmaS += S;
367 _SigmaZ += Z;
368 _SigmaSZ += S*Z;
369 _SigmaSS += S*S;
370 }
371 }
374 return 0;
375}
376
377int
379 IMessageSvc *msgSvc;
380 Gaudi::svcLocator()->service("MessageSvc", msgSvc);
381
382 MsgStream log(msgSvc, "FTFinder");
383
384 HepVector a = _la->Hpar(m_ftFinder->pivot);
385 int n = _stereo_segments->length();
386 log<<MSG::DEBUG<<"number of stereo segments: "<< n << endreq;
387 if (n < (param->_nseg)){
388 a(4) = -9999.;
389 a(5) = -9999.;
390 _helix = new Helix(m_ftFinder->pivot,a);
391 log<<MSG::DEBUG<<"cut by _nseg" << param->_nseg << endreq;
392 return 0;
393 }
394 FTList<double> zList(50);
395 FTList<double> sList(50);
396 FTList<FTWire *> hList(50);
397 for (int i = 0; i^n; i++){
398 FTList<FTWire *> & hits = (*_stereo_segments)[i]->wireHits();
399 int m = hits.length();
400 for (int j = 0; j^m; j++){
401 FTWire & h = * hits[j];
402 double z;
403 if (!(h.z(*_la,z))) continue;
404 double s = _la->s(h.layer().r());
405 float cellsize = h.layer().csize();
406 log<<MSG::DEBUG<<"cellsize: " << cellsize << endreq;
407
408 if (m_ftFinder->evtTiming){
409 float time = h.time() + m_ftFinder->evtTiming;
410 if(time<-100) continue; //discard the bad TDC time
411 double distance = m_ftFinder->t2x(h.layer(),time);
412 h.distance(distance);
413 log<<MSG::DEBUG<<"m_ftFinder->evtTiming: "<< m_ftFinder->evtTiming << " TDC time: " << h.time() << " distance: "<< distance << endreq;
414 }
415
416 double par = h.distance()/(0.25*cellsize);
417 _za->add(s,z,exp(-par*par));
418 sList.append(s);
419 zList.append(z);
420 hList.append(&h);
421 }
422 }
423 if (hList.length() < (param->_nlength)){
424 a(4) = -9999.;
425 a(5) = -9999.;
426 _helix = new Helix(m_ftFinder->pivot,a);
427 log << MSG::DEBUG << "cut by _nlength: " << param->_nlength << endreq;
428 return 0;
429 }
430 double chi2 = _za->calculate();
431 _za->clear();
432 n = hList.length();
433 for (int i = 0; i^n; i++){
434 double d = _za->d(sList[i],zList[i]);
435 float z_distance = hList[i]->distance_z();
436 if (fabs(fabs(d)-z_distance) > (param->_z_cut1)){
437 log<<MSG::DEBUG<<"cut by _z_cut1: "<< param->_z_cut1 << endreq;
438 zList.remove2(i);
439 sList.remove2(i);
440 n = hList.remove(i);
441 continue;
442 }
443 _za->add(sList[i], (d>0) ? zList[i]-z_distance : zList[i]+z_distance, 1.);
444 }
445
446 if (_za->nc() < (param->_nc)){
447 a(4) = -9999.;
448 a(5) = -9999.;
449 _helix = new Helix(m_ftFinder->pivot,a);
450 log<<MSG::DEBUG<<"cut by _nc: " << param->_nc << endreq;
451 return 0;
452 }
453 chi2 = _za->calculate();
454 _za->clear();
455 for (int i = 0; i^n; i++){
456 double d = _za->d(sList[i],zList[i]);
457 float z_distance = hList[i]->distance_z();
458 hList[i]->setChi2(fabs(d)-z_distance);
459 #ifndef OnlineMode
460 //fill ntuple
461 g_sigmaz->fill(fabs(d)-z_distance, 1.0);
462#endif
463 if (fabs(fabs(d)-z_distance) > (param->_z_cut2)) continue;
464 _za->add(sList[i], (d>0) ? zList[i]-z_distance : zList[i]+z_distance, 1.);
465 }
466
467 if (_za->nc() < (param->_nc)){
468 a(4) = -9999.;
469 a(5) = -9999.;
470 _helix = new Helix(m_ftFinder->pivot,a);
471 log<<MSG::DEBUG<<"cut by _nc" << param->_nc << endreq;
472 return 0;
473 }
474 chi2 = _za->calculate();
475#ifndef OnlineMode
476 g_chi2sz->fill(chi2/_za->nc(), 1.0);
477#endif
478 a(4) = _za->b(); // dz
479 a(5) = _za->a(); // tanLambda
480 _helix = new Helix(m_ftFinder->pivot,a);
481 return 1;
482}
483
485{
486 int nhits = 0;
487 int n = _axial_segments.length();
488 for (int i = 0; i^n; i++) {
489 nhits += _axial_segments[i]->wireHits().length();
490 }
492 for (int j = 0; j^n; j++) {
493 nhits += (*_stereo_segments)[j]->wireHits().length();
494 }
495 return nhits;
496}
497
498#ifndef OnlineMode
500 int n = axial_segments().length();
501 for (int i = 0; i^n; i++) axial_segments()[i]->printout();
503 for (int i = 0; i^n; i++) stereo_segments()[i]->printout();
504}
505#endif
int selected
const Int_t n
Double_t x[10]
Double_t time
EvtComplex exp(const EvtComplex &c)
Definition: EvtComplex.hh:252
XmlRpcServer s
Definition: HelloServer.cpp:11
FTSuperLayer * superLayer(int id) const
returns superlayer
float t2x(const FTLayer &l, const float t) const
convert t to x
const float r(void) const
returns r form origin
double csize(void) const
returns cell size
void deleteAll(void)
delete all object and clear(allocated memory remains same)
void remove2(int)
remove objects by index
int remove(int &)
remove objects by index and returns decremented index and length
int length(void) const
returns the length of the list
int append(const T &x)
append an object into the end of the list
FTList< FTSegment * > & complecated_segments(void)
returns complecated segments
int r_phiReFit(float vx, float xy, int vtx_flag)
do r-phi refit
Definition: FTTrack.cxx:144
int get_nhits(void)
calculate the wire hits number
Definition: FTTrack.cxx:484
FTList< FTSegment * > & axial_segments(void) const
returns axial segments
int r_phi4Fit(float vx, float xy, int vtx_flag)
Definition: FTTrack.cxx:295
int r_phiFit(void)
do r-phi circle fit
Definition: FTTrack.cxx:37
void printout()
Definition: FTTrack.cxx:499
FTList< FTList< FTSegment * > * > * _stereo_segments_by_superLayer
FTList< FTSegment * > & stereo_segments(void) const
returns stereo_segments
int r_phi3Fit(int l, float vx, float xy, int vtx_flag)
Definition: FTTrack.cxx:249
int r_phi2Fit(float vx, float xy, int vtx_flag)
Definition: FTTrack.cxx:191
int s_zFit(void)
do s-z linear fit
Definition: FTTrack.cxx:378
int linkStereoSegments(void)
link stereo segments by tanLambda
Definition: FTTrack.cxx:342
const float y(void) const
returns position y
int z(const Lpav &la, double &z) const
returns z for track la
float distance(void) const
returns drift distance
unsigned stateAND(const unsigned mask) const
returns state bit
const FTLayer & layer(void) const
returns layer
const float x(void) const
returns position x
void stateOR(const unsigned mask)
set state bit
float time(void) const
rerurns TDC time(after t0 subtraction)
void setChi2(float chi2)
set residual fit chi2
virtual double getReferField()=0
double d(double x, double y) const
HepVector Hpar(const HepPoint3D &pivot) const
double s(double x, double y) const
void add_point(double x, double y, double w=1)
static MdcParameter * instance()
Definition: MdcParameter.cxx:6
void add(double, double, double)
double d(double s, double z) const