00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #include "TrkReco/TCircle.h"
00014 #include "TrkReco/TMDCWire.h"
00015 #include "TrkReco/TMDCWireHit.h"
00016 #include "TrkReco/TMDCWireHitMC.h"
00017 #include "TrkReco/TMLink.h"
00018
00019 const TCircleFitter
00020 TCircle::_fitter = TCircleFitter("TCircle Default Circle Fitter");
00021
00022 TCircle::TCircle(const AList<TMLink> & a)
00023 : TTrackBase(a), _radius(0.), _charge(0.) {
00024
00025
00026 fitter(& TCircle::_fitter);
00027 }
00028
00029 TCircle::~TCircle() {
00030 }
00031
00032 void
00033 TCircle::dump(const std::string & msg, const std::string & pre) const {
00034 bool def = false;
00035 if (msg == "") def = true;
00036
00037 if (def || msg.find("circle") != std::string::npos || msg.find("detail") != std::string::npos) {
00038 std::cout << pre;
00039 std::cout << "#links=" << _links.length();
00040 if (_fitted) {
00041 std::cout << ",charge=" << _charge;
00042 std::cout << ",center=" << _center;
00043 std::cout << ",radius=" << _radius;
00044 std::cout << std::endl << pre;
00045 std::cout << "pt=" << pt();
00046 std::cout << ",impact=" << impact();
00047 std::cout << std::endl;
00048 }
00049 else {
00050 std::cout << ", not fitted yet" << std::endl;
00051 }
00052 }
00053 if (! def) TTrackBase::dump(msg, pre);
00054 }
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110 double
00111 TCircle::weight(const TMLink & l) const {
00112
00113
00114 int maxLink = 0;
00115 int localID[6];
00116 int layerID[6];
00117 int LayerID = l.hit()->wire()->layerId();
00118 for(int i=0;i<6;i++){
00119 if(l.neighbor(i)){
00120 maxLink = i;
00121 localID[i] = l.neighbor(i)->hit()->wire()->localId();
00122 layerID[i] = l.neighbor(i)->hit()->wire()->layerId();
00123 }else{
00124 break;
00125 }
00126 }
00127 if(maxLink != 1)return 1.0;
00128 if(layerID[0] == LayerID &&
00129 layerID[1] == LayerID)return 0.5;
00130 if(layerID[0]+layerID[1] != LayerID*2)return 1.0;
00131
00132 return 1.0;
00133 }
00134
00135 int
00136 TCircle::fitForCurl(int ipConst) {
00137 unsigned n = _links.length();
00138 if (n < 3) return -1;
00139
00140
00141 unsigned flagIP = 1;
00142 unsigned layerID = _links[0]->hit()->wire()->layerId();
00143 for(unsigned i = 0; i < n; i++){
00144 if(layerID != _links[i]->hit()->wire()->layerId()){
00145 flagIP = 0;
00146 break;
00147 }
00148 }
00149 if(ipConst != 0)flagIP = 1;
00150 if(flagIP == 1){
00151 _circle.add_point(0., 0., 0.5);
00152
00153 }
00154 for (unsigned i = 0; i < n; i++) {
00155 TMLink * l = _links[i];
00156 if (l == 0) continue;
00157
00158 const TMDCWireHit * h = l->hit();
00159 if (h == 0) continue;
00160
00161
00162 HepPoint3D point;
00163 point = h->xyPosition();
00164
00165 double weight = 1.0;
00166 weight = this->weight(* l);
00167
00168 _circle.add_point(point.x(), point.y(), weight);
00169
00170 }
00171
00172
00173 if(_circle.fit()<0.0 || _circle.kappa()==0.0) return -1;
00174 HepVector v(_circle.center());
00175 _center.setX(v(1));
00176 _center.setY(v(2));
00177 _radius = _circle.radius();
00178
00179
00180 int qSum = 0;
00181 for (unsigned i = 0; i < n; i++) {
00182 TMLink * l = _links[i];
00183 if (l == 0) continue;
00184
00185 const TMDCWireHit * h = l->hit();
00186 if (h == 0) continue;
00187
00188 float q = (_center.cross(h->xyPosition())).z();
00189 if (q > 0.) qSum += 1;
00190 else qSum -= 1;
00191 }
00192 if (qSum >= 0) _charge = +1.;
00193 else _charge = -1.;
00194 _radius *= _charge;
00195
00196 _fitted = true;
00197 return 0;
00198 }