CGEM BOSS 6.6.5.f
BESIII Offline Software System
Loading...
Searching...
No Matches
HoughTransAlg/HoughTransAlg-00-00-14/src/HoughHit.cxx
Go to the documentation of this file.
1#include "HoughTransAlg/HoughHit.h"
2#include "Identifier/MdcID.h"
3#include "RawEvent/RawDataUtil.h"
4#include "MdcGeom/Constants.h"
5#include "TrkBase/TrkPoca.h"
6#include "TrkBase/HelixTraj.h"
7#include "TrkFitter/TrkCircleTraj.h"
8#include "TrkBase/TrkExchangePar.h"
9
10#include <iomanip>
11#include <string>
12
13using namespace std;
14
15MdcGeomSvc* HoughHit::m_mdcGeomSvc = NULL;
16MdcCalibFunSvc* HoughHit::m_mdcCalibFunSvc = NULL;
17CgemGeomSvc* HoughHit::m_cgemGeomSvc = NULL;
18CgemCalibFunSvc* HoughHit::m_cgemCalibFunSvc = NULL;
19const MdcDetector* HoughHit::m_mdcDetector = NULL;
20HoughHit::HoughHit(){
21 //m_det=Global::m_gm;
22 m_cgemCluster=NULL;
23 m_mdcDigi=NULL;
24 m_cgemMcHit=NULL;
25 m_mdcMcHit=NULL;
26 //m_hitMap = NULL;
27}
28
29HoughHit::HoughHit(const MdcDigi* mdcDigi, double bunchTime, int hitID)
30{
31 m_hitID = hitID;
32 m_hitType = MDCHIT;
33 m_cgemCluster = NULL;
34 m_mdcDigi = mdcDigi;
35 m_cgemMcHit=NULL;
36 m_mdcMcHit=NULL;
37 m_layer = MdcID::layer(mdcDigi->identify());
38 m_wire = MdcID::wire(mdcDigi->identify());
39 m_flag = flag();
40 m_use = 0;
41 m_bunchTime = bunchTime;
42 m_rawTime = RawDataUtil::MdcTime(mdcDigi->getTimeChannel());
43 m_depositEnergy = mdcDigi->getChargeChannel();
44 const MdcGeoWire* wire = m_mdcGeomSvc->Wire(m_layer,m_wire);
45 m_hitPosition = (wire->Forward() + wire->Backward())/10/2;
46 m_westPoint = wire->Forward();
47 m_eastPoint = wire->Backward();
48 m_driftDist = driftDistance();
49 m_residual = 9999;
50 //m_hitMap = NULL;
51 //m_sLeft = 0;
52 //m_zLeft = 0;
53 //m_sRight = 0;
54 //m_zRight = 0;
55 m_trkID.clear();
56 m_sz.clear();
57 m_pairHit = NULL;
58 m_halfCircle = 0;
59 m_position.clear();
60 m_mdcHit = NULL;
61 //MdcHit mdcHit = *(new MdcHit(m_mdcDigi,m_mdcDetector));
62 //mdcHit.setCalibSvc(m_mdcCalibFunSvc);
63 //mdcHit.setCountPropTime(true);
64}
65
66HoughHit::HoughHit(const RecCgemCluster* cgemCluster, double bunchTime, int hitID)
67{
68 m_hitID = hitID;
69 m_hitType = CGEMHIT;
70 m_cgemCluster = cgemCluster;
71 m_mdcDigi=NULL;
72 m_cgemMcHit=NULL;
73 m_mdcMcHit=NULL;
74 m_layer = cgemCluster->getlayerid();
75 m_wire = cgemCluster->getsheetid();
76 //m_wire = cgemCluster->getclusterid();
77 //m_wire = -1;
78 m_flag = cgemCluster->getflag();
79 m_use = 0;
80 m_bunchTime = bunchTime;
81 m_rawTime = 0; //FIXME
82 m_depositEnergy = cgemCluster->getenergydeposit(); //FIXME
83 m_hitPosition.setX((m_cgemGeomSvc->getCgemLayer(m_layer)->getMiddleROfGapD()/10.)*cos(cgemCluster->getrecphi()));
84 m_hitPosition.setY((m_cgemGeomSvc->getCgemLayer(m_layer)->getMiddleROfGapD()/10.)*sin(cgemCluster->getrecphi()));
85 m_hitPosition.setZ(cgemCluster->getRecZ()/10.);
86 m_westPoint = HepPoint3D(0,0,0);
87 m_eastPoint = HepPoint3D(0,0,0);
88 m_driftDist = 0;
89 m_residual = 9999;
90 //m_hitMap = NULL;
91 //m_sLeft = 0;
92 //m_zLeft = 0;
93 //m_sRight = 0;
94 //m_zRight = 0;
95 m_trkID.clear();
96 m_sz.clear();
97 m_pairHit = NULL;
98 m_halfCircle = 0;
99 m_position.clear();
100 m_mdcHit = NULL;
101}
102
103HoughHit::HoughHit(const MdcMcHit* mdcMcHit, double bunchTime, int hitID)
104{
105 m_hitID = hitID;
106 m_hitType = MDCMCHIT;
107 m_cgemCluster = NULL;
108 m_mdcDigi = NULL;
109 m_cgemMcHit = NULL;
110 m_mdcMcHit = mdcMcHit;
111 m_layer = MdcID::layer(mdcMcHit->identify());
112 m_wire = MdcID::wire(mdcMcHit->identify());
113 m_flag = flag();
114 m_use = mdcMcHit->getPositionFlag();
115 m_bunchTime = bunchTime;
116 m_rawTime = 0;
117 m_depositEnergy = mdcMcHit->getDepositEnergy();
118 m_hitPosition.setX(mdcMcHit->getPositionX()/10.);
119 m_hitPosition.setY(mdcMcHit->getPositionY()/10.);
120 m_hitPosition.setZ(mdcMcHit->getPositionZ()/10.);
121 m_westPoint = HepPoint3D(0,0,0);
122 m_eastPoint = HepPoint3D(0,0,0);
123 m_driftDist = mdcMcHit->getDriftDistance()/10.;
124 m_residual = 9999;
125 //m_hitMap = NULL;
126 //m_sLeft = 0;
127 //m_zLeft = 0;
128 //m_sRight = 0;
129 //m_zRight = 0;
130 m_trkID.push_back(mdcMcHit->getTrackIndex());
131 //cout<<m_layer<<" "<<mdcMcHit->getTrackIndex()<<endl;
132 m_sz.clear();
133 m_pairHit = NULL;
134 m_halfCircle = 0;
135 m_position.clear();
136 m_mdcHit = NULL;
137}
138
139HoughHit::HoughHit(const CgemMcHit* cgemMcHit, double bunchTime, int hitID)
140{
141 m_hitID = hitID;
142 m_hitType = CGEMMCHIT;
143 m_cgemCluster = NULL;
144 m_mdcDigi = NULL;
145 m_cgemMcHit = cgemMcHit;
146 m_mdcMcHit = NULL;
147 m_layer = cgemMcHit->GetLayerID();
148 //m_wire = cgemMcHit->GetParentID();
149 //m_flag = cgemMcHit->GetPDGCode();
150 m_wire = -1;
151 m_flag = 2;
152 m_use = 0;
153 m_bunchTime = bunchTime;
154 m_rawTime = 0;
155 m_depositEnergy = 0;
156 m_hitPosition.setX((cgemMcHit->GetPositionXOfPrePoint()+cgemMcHit->GetPositionXOfPostPoint())/2/10.);
157 m_hitPosition.setY((cgemMcHit->GetPositionYOfPrePoint()+cgemMcHit->GetPositionYOfPostPoint())/2/10.);
158 m_hitPosition.setZ((cgemMcHit->GetPositionZOfPrePoint()+cgemMcHit->GetPositionZOfPostPoint())/2/10.);
159 m_westPoint = HepPoint3D(0,0,0);
160 m_eastPoint = HepPoint3D(0,0,0);
161 m_driftDist = 0;
162 m_residual = 9999;
163 //m_hitMap = NULL;
164 //m_sLeft = 0;
165 //m_zLeft = 0;
166 //m_sRight = 0;
167 //m_zRight = 0;
168 m_trkID.push_back(cgemMcHit->GetTrackID());
169 //cout<<m_layer<<" "<<cgemMcHit->GetTrackID()<<endl;
170 m_sz.clear();
171 m_pairHit = NULL;
172 m_halfCircle = 0;
173 m_position.clear();
174 m_mdcHit = NULL;
175}
176
177HoughHit::HoughHit(const HoughHit& other):
178 m_hitID(other.m_hitID),
179 m_hitType(other.m_hitType),
180 m_cgemCluster(other.m_cgemCluster),
181 m_mdcDigi(other.m_mdcDigi),
182 m_cgemMcHit(other.m_cgemMcHit),
183 m_mdcMcHit(other.m_mdcMcHit),
184 m_layer(other.m_layer),
185 m_wire(other.m_wire),
186 m_flag(other.m_flag),
187 m_use(other.m_use),
188 m_bunchTime(other.m_bunchTime),
189 m_rawTime(other.m_rawTime),
190 m_depositEnergy(other.m_depositEnergy),
191 m_hitPosition(other.m_hitPosition),
192 m_westPoint(other.m_westPoint),
193 m_eastPoint(other.m_eastPoint),
194 m_driftDist(other.m_driftDist),
195 m_residual(other.m_residual),
196 //m_hitMap(other.m_hitMap),
197 //m_sLeft(other.m_sLeft),
198 //m_zLeft(other.m_zLeft),
199 //m_sRight(other.m_sRight),
200 //m_zRight(other.m_zRight),
201 m_trkID(other.m_trkID),
202 m_sz(other.m_sz),
203 m_pairHit(other.m_pairHit),
204 m_halfCircle(other.m_halfCircle),
205 m_position(other.m_position),
206 m_mdcHit(other.m_mdcHit)
207{}
208
209HoughHit& HoughHit::operator=(const HoughHit& other)
210{
211 m_hitID = other.m_hitID;
212 m_hitType = other.m_hitType;
213 m_cgemCluster = other.m_cgemCluster;
214 m_mdcDigi = other.m_mdcDigi;
215 m_cgemMcHit = other.m_cgemMcHit;
216 m_mdcMcHit = other.m_mdcMcHit;
217 m_layer = other.m_layer;
218 m_wire = other.m_wire;
219 m_flag = other.m_flag;
220 m_use = other.m_use;
221 m_bunchTime = other.m_bunchTime;
222 m_rawTime = other.m_rawTime;
223 m_depositEnergy = other.m_depositEnergy;
224 m_hitPosition = other.m_hitPosition;
225 m_westPoint = other.m_westPoint;
226 m_eastPoint = other.m_eastPoint;
227 m_driftDist = other.m_driftDist;
228 m_residual = other.m_residual;
229 //m_hitMap = other.m_hitMap;
230 //m_sLeft = other.m_sLeft;
231 //m_zLeft = other.m_zLeft;
232 //m_sRight = other.m_sRight;
233 //m_zRight = other.m_zRight;
234 m_trkID = other.m_trkID;
235 m_sz = other.m_sz;
236 m_pairHit = other.m_pairHit;
237 m_halfCircle = other.m_halfCircle;
238 m_position = other.m_position;
239 m_mdcHit =other.m_mdcHit;
240}
241
242/*
243void HoughHit::buildMap(int x_bin, double x_min, double x_max, int y_bin, double y_min, double y_max, int nPoint, int charge)
244{
245 if(m_hitMap!=NULL)delete m_hitMap;
246 stringstream ssname;
247 ssname<<"id"<<m_hitID<<",layer"<<m_layer<<",wire(sheet)"<<m_wire;
248 string sname = ssname.str();
249 const char * name = sname.c_str();
250 m_hitMap = new TH2D(name,name, x_bin, x_min, x_max, y_bin, y_min, y_max);
251 int N = x_bin*nPoint;
252 double delta_alpha = (x_max - x_min )/N;
253 double alpha = x_min-delta_alpha/2;
254 double u = 2*m_hitPosition.x()/(m_hitPosition.x()*m_hitPosition.x()+m_hitPosition.y()*m_hitPosition.y()-m_driftDist*m_driftDist);
255 double v = 2*m_hitPosition.y()/(m_hitPosition.x()*m_hitPosition.x()+m_hitPosition.y()*m_hitPosition.y()-m_driftDist*m_driftDist);
256 double w = 2*m_driftDist/(m_hitPosition.x()*m_hitPosition.x()+m_hitPosition.y()*m_hitPosition.y()-m_driftDist*m_driftDist);
257 for(int i=0;i<N;i++){
258 alpha += delta_alpha;
259 if(alpha>M_PI)alpha -= M_PI;
260 double rho1 = u*cos(alpha) + v*sin(alpha) + w;
261 double rho2 = u*cos(alpha) + v*sin(alpha) - w;
262 //cout<<alpha<<" "<<rho1<<" "<<rho2<<endl;
263 double slantOfLine = v*cos(alpha) - u*sin(alpha);
264 int chargeOfHitOnCir1 = (slantOfLine/fabs(slantOfLine))*(rho1/fabs(rho1));
265 if(chargeOfHitOnCir1 != charge && y_min<(rho1) && (rho1)<y_max){
266 //if( y_min<(rho1) && (rho1)<y_max){
267 m_hitMap->Fill(alpha,rho1);
268 }
269 int chargeOfHitOnCir2 = (slantOfLine/fabs(slantOfLine))*(rho2/fabs(rho2));
270 if(chargeOfHitOnCir2 != charge && y_min<(rho2) && (rho2)<y_max){
271 //if( y_min<(rho2) && (rho2)<y_max){
272 if(m_hitType == CGEMHIT)continue;
273 //int bin = m_hitMap->FindBin(alpha,rho2);
274 //int bin = m_hitMap->FindFixBin(alpha,rho2);
275 //if(m_hitMap->GetBinContent(bin)>0)continue;
276 m_hitMap->Fill(alpha,rho2);
277 }
278 }
279}
280
281void HoughHit::clearMap()
282{
283 if(m_hitMap!=NULL)delete m_hitMap;
284 m_hitMap = NULL;
285}
286*/
287
288int HoughHit::flag()
289{
290 double nomShift = m_mdcGeomSvc->Layer(m_layer)->nomShift();
291 if(nomShift>0) return 1;
292 else if(nomShift<0) return -1;
293 else return 0;
294}
295
296double HoughHit::driftTime()
297{
298 double ToF = 1.e9*sqrt(m_hitPosition.x()*m_hitPosition.x()+m_hitPosition.y()*m_hitPosition.y())/Constants::c;
299 double Tprop = m_mdcCalibFunSvc->getTprop(m_layer,m_hitPosition.z()*10.);
300 double Twalk = m_mdcCalibFunSvc->getTimeWalk(m_layer, m_depositEnergy);
301 double T0 = m_mdcCalibFunSvc->getT0(m_layer,m_wire);
302 double driftTime = m_rawTime - 1.e9*m_bunchTime - ToF - Tprop - Twalk -T0;
303 //cout<<m_rawTime - Twalk - T0 - ToF - Tprop<<endl;
304 //cout<<"("<<m_layer<<", "<<m_wire<<"): ";
305 //cout<<"bunchT = "<<1.e9*m_bunchTime<<" , rawT "<<m_rawTime<<" , Tof "<<ToF<<" , Tprop "<<Tprop<<" , T0walk "<<Twalk+T0<<" , driftT "<<driftTime<<" , rmid "<<sqrt(m_hitPosition.x()*m_hitPosition.x()+m_hitPosition.y()*m_hitPosition.y())<<" , z "<<m_hitPosition.z()*10.0<<endl;
306 //cout<<endl;
307 return driftTime;
308}
309
310double HoughHit::driftDistance()
311{
312 //if(fabs(m_hitPosition.z())>150. || driftTime>1500.)return -1;
313 double driftT = driftTime();
314 if(fabs(driftT)>1500.||fabs(m_hitPosition.z()>150.)){
315 m_use = -1;
316 return 9999.;
317 }
318 double driftDist(-1);
319 int lrCalib=2;
320 //if (m_flag==1) lrCalib = 0;
321 //else if (m_flag==-1) lrCalib = 1;
322 double entranceAngle = 0;
323 driftDist = 0.1 * m_mdcCalibFunSvc->driftTimeToDist(driftT,m_layer,m_wire,lrCalib,entranceAngle);
324 return driftDist;
325}
326
327void HoughHit::print()
328{
329 string hitType;
330 switch(m_hitType){
331 case CGEMHIT : hitType = "CGEMHIT"; break;
332 case MDCHIT : hitType = "MDCHIT"; break;
333 case CGEMMCHIT : hitType = "CGEMMCHIT"; break;
334 case MDCMCHIT : hitType = "MDCMCHIT"; break;
335 default : hitType = "";
336 }
337 cout<<"ID:" <<setw(6)<<m_hitID;
338 cout<<"(" <<setw(2)<<m_layer;
339 cout<<"," <<setw(3)<<m_wire<<") ";
340 cout<<" flag:" <<setw(3)<<m_flag;
341 if(m_hitType==CGEMHIT||m_hitType==MDCHIT){
342 if(m_pairHit!=NULL){
343 cout<<" half:"<<setw(2)<<m_pairHit->getHalfCircle();
344 cout<<" trkID:"<<setw(3)<<(m_pairHit->getTrkID())[0];
345 cout<<" McHit:("<<setw(10)<<m_pairHit->getHitPosition().x()<<", "<<setw(10)<<m_pairHit->getHitPosition().y()<<", "<<setw(10)<<m_pairHit->getHitPosition().z()<<") ";
346 if(m_hitType==MDCHIT){
347 cout<<" McLR:"<<setw(3)<<m_pairHit->getMdcMcHit()->getPositionFlag();
348 cout<<" drift:"<<setw(10)<<m_pairHit->getMdcMcHit()->getDriftDistance();
349 }
350 //if(m_hitType==MDCHIT)
351 //cout<<" trkID:"<<setw(3)<<m_pairHit->getMdcMcHit()->getTrackIndex();
352 }
353 }
354 if(m_hitType==CGEMMCHIT||m_hitType==MDCMCHIT){
355 cout<<" half:"<<setw(2)<<getHalfCircle();
356 cout<<" trkID:"<<setw(3)<<m_trkID[0];
357 cout<<" McHit:("<<setw(10)<<m_hitPosition.x()<<", "<<setw(10)<<m_hitPosition.y()<<", "<<setw(10)<<m_hitPosition.z()<<") ";
358 //if(m_hitType==MDCMCHIT)
359 //cout<<" trkID:"<<setw(3)<<m_mdcMcHit->getTrackIndex();
360 if(m_hitType==MDCMCHIT){
361 cout<<" McLR:"<<setw(3)<<m_mdcMcHit->getPositionFlag();
362 cout<<" drift:"<<setw(10)<<m_mdcMcHit->getDriftDistance();
363 }
364 }
365
366 //cout<<setw(7)<<"hitID:" <<setw(6)<<m_hitID
367 //<<setw(10)<<"hitType:" <<setw(10)<<hitType
368 //<<setw(8)<<"layer:" <<setw(4)<<m_layer
369 //<<setw(6)<<"wire:" <<setw(6)<<m_wire
370 //<<setw(10)<<"flag:" <<setw(10)<<m_flag
371 //<<setw(5)<<"use:" <<setw(3)<<m_use
372 //<<endl
373 //<<setw(12)<<"bunchTime:"<<setw(12)<<m_bunchTime
374 //<<setw(12)<<"rawTime:" <<setw(12)<<m_rawTime
375 //<<setw(12)<<"charge:" <<setw(12)<<m_depositEnergy
376 //<<setw(12)<<"driftTime:"<<setw(12)<<driftTime()
377 //<<setw(11)<<"driftDist:"<<setw(12)<<m_driftDist
378 //<<setw(12)<<"residual:" <<setw(12)<<m_residual
379 //<<endl
380 //<<setw(10)<<"sLeft:" <<setw(10)<<m_sLeft
381 //<<setw(10)<<"zLeft:" <<setw(10)<<m_zLeft
382 //<<setw(10)<<"sRight:" <<setw(10)<<m_sRight
383 //<<setw(10)<<"zRight:" <<setw(10)<<m_zRight
384 //<<endl
385 //<<endl
386 //cout<<m_trkID.size()<<" trkID:";
387 //for(int i=0;i<m_trkID.size();i++)cout<<setw(4)<<m_trkID[i]<<endl;
388 //;
389 //if(m_layer>3)cout<<"TrkID:"<<setw(4)<<getDigi()->getTrackIndex();
390 //cout<<endl;
391 //;
392 cout<<endl;
393}
394
395
396// erase trk id and residual from the Hough hit
397void HoughHit::dropTrkID(int trkID)
398{
399 vector<int>::iterator it0 = m_trkID.begin();
400 vector<int>::iterator result = find(m_trkID.begin(),m_trkID.end(),trkID);
401 if(result!=m_trkID.end()) {
402 m_trkID.erase(result);
403 vector<double>::iterator itRes = m_vecResid.begin()+(result-it0);
404 if(itRes!=m_vecResid.end()) m_vecResid.erase(itRes);
405 }
406}
407
408vector<HepPoint3D>& HoughHit::VHitPosition(HoughTrack* track)
409{
410 if(getFlag()==0)return m_position;
411 m_position.clear();
412 //if(getPairHit()==NULL)return m_position;
413 //if(getPairHit()->getHalfCircle()!=1)return m_position;
414 double z(0);
415 double xc = track->center().x();
416 double yc = track->center().y();
417 double rTrack = track->radius();// signed FIXME
418 if(getHitType()==HoughHit::MDCHIT){
419 resetSZ();
420 double drift = getDriftDist();
421 const MdcGeoWire* wire = getMdcGeomSvc()->Wire(getLayer(),getWire());
422 HepPoint3D westPoint = wire->Forward();
423 HepPoint3D eastPoint = wire->Backward();
424 double xEast = eastPoint.x()/10.0;
425 double xWest = westPoint.x()/10.0;
426 double yEast = eastPoint.y()/10.0;
427 double yWest = westPoint.y()/10.0;
428 double zEast = eastPoint.z()/10.0;
429 double zWest = westPoint.z()/10.0;
430 //cout<<"wast:x,y,z: "<<xWest<<", "<<yWest<<", "<<zWest<<endl;
431 //cout<<"east:x,y,z: "<<xEast<<", "<<yEast<<", "<<zEast<<endl;
432 //cout<<"Xc, Yc, Rc: "<<xc<<", "<<yc<<", "<<rTrack<<endl;
433 double west2east = sqrt((xEast-xWest)*(xEast-xWest)+(yEast-yWest)*(yEast-yWest));
434
435 double slope = (yEast-yWest)/(xEast-xWest);
436 double intercept = (yWest-slope*xWest+yEast-slope*xEast)/2;
437 double a = 1+slope*slope;
438 double b = -2*(xc+slope*yc-slope*intercept);
439 double c1 = xc*xc+(yc-intercept)*(yc-intercept)-(rTrack+drift)*(rTrack+drift);
440 double c2 = xc*xc+(yc-intercept)*(yc-intercept)-(rTrack-drift)*(rTrack-drift);
441 //double c1 = intercept*(2*yc-intercept)+drift*(drift+rTrack);
442 //double c2 = intercept*(2*yc-intercept)+drift*(drift-rTrack);
443 double delta1 = (b*b-4*a*c1);
444 double delta2 = (b*b-4*a*c2);
445 //cout<<"a,b,c1,c2: "<<a<<", "<<b<<", "<<c1<<", "<<c2<<endl;
446 //cout<<"delta: "<<delta1<<", "<<delta2<<endl;
447 if(delta1>=0){
448 double x1 = (-b+sqrt(delta1))/(2*a);
449 double x2 = (-b-sqrt(delta1))/(2*a);
450 double y1 = slope*x1+intercept;
451 double y2 = slope*x2+intercept;
452 if((x1-xWest)*(x1-xEast)<0){
453 double l = sqrt((x1-xWest)*(x1-xWest)+(y1-yWest)*(y1-yWest));
454 z = zWest + l/west2east*fabs((zEast-zWest));
455 HepPoint3D position(x1,y1,z);
456 m_position.push_back(position);
457 }
458 if((x2-xWest)*(x2-xEast)<0){
459 double l = sqrt((x2-xWest)*(x2-xWest)+(y2-yWest)*(y2-yWest));
460 z = zWest + l/west2east*fabs((zEast-zWest));
461 HepPoint3D position(x2,y2,z);
462 m_position.push_back(position);
463 }
464 }
465
466 if( delta2>=0){
467 double x1 = (-b+sqrt(delta2))/(2*a);
468 double x2 = (-b-sqrt(delta2))/(2*a);
469 double y1 = slope*x1+intercept;
470 double y2 = slope*x2+intercept;
471 if((x1-xWest)*(x1-xEast)<0){
472 double l = sqrt((x1-xWest)*(x1-xWest)+(y1-yWest)*(y1-yWest));
473 z = zWest + l/west2east*fabs((zEast-zWest));
474 HepPoint3D position(x1,y1,z);
475 m_position.push_back(position);
476 }
477 if((x2-xWest)*(x2-xEast)<0){
478 //double l = zWest + sqrt((x2-xWest)*(x2-xWest)+(y2-yWest)*(y2-yWest));
479 //z = l/west2east*fabs((zEast-zWest));
480 double l = sqrt((x2-xWest)*(x2-xWest)+(y2-yWest)*(y2-yWest));
481 z = zWest + l/west2east*fabs((zEast-zWest));
482 HepPoint3D position(x2,y2,z);
483 m_position.push_back(position);
484 }
485 }
486 }
487 return m_position;
488}
489///*
490double HoughHit::residual(HoughTrack* track)
491{
492 double residual(9999.);
493 if(m_flag==0){
494 HepPoint3D hitPoint = getHitPosition();
495 if(getHitType() == HoughHit::MDCHIT || getHitType() == HoughHit::MDCMCHIT){
496 //HepPoint3D circleCenter = m_helix.center();
497 HepPoint3D circleCenter = track->center();
498 //double distance = circleCenter.perp(hitPoint);
499 double distance = (circleCenter-hitPoint).perp();
500 //double Rc = m_helix.radius();
501 double Rc = fabs(track->radius());
502 double driftDist(0);
503 //if(getHitType() == HoughHit::MDCMCHIT)driftDist = 0;
504 if(getHitType() == HoughHit::MDCHIT)driftDist = getDriftDist();
505 //residual = fabs(distance - Rc) - driftDist;
506 residual = driftDist-fabs(distance - Rc);
507 //cout<<"driftDist, distance, Rc = "<<driftDist<<", "<<distance<<", "<<Rc<<endl;
508 //double distance = (m_helix.center()).perp(hit.getHitPosition());
509 //res = fabs(distance - m_helix.radius()) - hit.getDriftDist();
510 }else if(getHitType() == HoughHit::CGEMHIT || getHitType() == HoughHit::CGEMMCHIT){
511 double rCgem = hitPoint.perp();
512 //double phiTrkFlt = m_helix.IntersectCylinder(rCgem);
513 double phiTrkFlt = track->IntersectCylinder(rCgem);
514 phiTrkFlt = track->judgeHalf(this)*phiTrkFlt;
515 //HepPoint3D crossPoint = m_helix.x(phiTrkFlt);
516 HepPoint3D crossPoint = track->x(phiTrkFlt);
517 double phiCrossPoint = crossPoint.phi();
518 double phiMeasure = hitPoint.phi();
519 residual = phiMeasure - phiCrossPoint;
520 while(residual<-M_PI)residual += 2*M_PI;
521 while(residual> M_PI)residual -= 2*M_PI;
522 residual = rCgem*residual;
523 }
524 }else{
525 double zTruth = getPairHit()->getHitPosition().z();
526 double minRes(9999);
527 if(getHitType() == HoughHit::MDCHIT || getHitType() == HoughHit::MDCMCHIT){
528 vector<HepPoint3D> position = VHitPosition(track);
529 for(vector<HepPoint3D>::iterator posIt = position.begin();posIt!=position.end();posIt++){
530 HepPoint3D point = *posIt;
531 double s = track->flightArc(point);
532 double zTrk = track->dz()+s*track->tanl();
533 double res = point.z() - zTrk;
534 //res = point.z() - zTruth;
535 if(fabs(res)<minRes){
536 minRes=fabs(res);
537 residual = res;
538 }
539 }
540 }else if(getHitType() == HoughHit::CGEMHIT || getHitType() == HoughHit::CGEMMCHIT){
541 double s = track->flightArc(m_hitPosition);
542 double zTrk = track->dz()+s*track->tanl();
543 residual = m_hitPosition.z() - zTrk;
544 //residual = m_hitPosition.z() - zTruth;
545 }
546 }
547 return residual;
548}
549
550void HoughHit::updateVHit(HoughTrack* track)
551{
552 if(m_flag!=0&&m_hitType==HoughHit::MDCHIT){
553 double minRes(9999);
554 HepPoint3D hitPoint(m_hitPosition);
555 vector<HepPoint3D> position = VHitPosition(track);
556 for(vector<HepPoint3D>::iterator posIt = position.begin();posIt!=position.end();posIt++){
557 HepPoint3D point = *posIt;
558 double s = track->flightArc(point);
559 double zTrk = track->dz()+s*track->tanl();
560 double res = point.z() - zTrk;
561 if(fabs(res)<minRes){
562 minRes=fabs(res);
563 hitPoint = point;
564 }
565 }
566 m_hitPosition = hitPoint;
567 m_driftDist = driftDistance();
568 }
569}
570//*/
571double HoughHit::residual(HoughTrack* track, HepPoint3D& positionOntrack, HepPoint3D& positionOnDetector)
572{
573 double residual(9999.);
574 HepPoint3D pivot(0,0,0);
575 track->pivot(pivot);
576 if(m_hitType==MDCHIT){
577 const MdcDetector* mdcDetector(m_mdcDetector);
578 MdcHit mdcHit(m_mdcDigi,mdcDetector);
579 const Trajectory* wire(mdcHit.hitTraj());
580
581 double d0 = -track->dr();
582 double phi0 = track->phi0()+M_PI/2;
583 double omega = track->kappa()/fabs(track->alpha());
584 double z0 = track->dz();
585 double tanl = track->tanl();
586
587 //HepVector a(5,0);
588 //a(1) = d0;
589 //a(2) = phi0;
590 //a(3) = omega;
591 //a(4) = z0;
592 //a(5) = tanl;
593 //HepSymMatrix Ea(5,0);
594 //const HepVector pvec(a);
595 //const HepSymMatrix pcov(Ea);
596 //const HepPoint3D refpoint(track->pivot());
597 //const HelixTraj trk(pvec, pcov, -99999., 99999., refpoint);
598
599 double hitWireLength = m_hitPosition.z()-m_westPoint.z();
600 const TrkExchangePar trkExchangePar(d0,phi0,omega,z0,tanl);
601 double trkFlightLength = track->flightLength(m_hitPosition);
602 const HelixTraj helixTraj(trkExchangePar);
603 TrkPoca trkPoca(helixTraj,trkFlightLength,*wire,hitWireLength);
604 trkFlightLength = trkPoca.flt1();
605 hitWireLength = trkPoca.flt2();
606 double doca = trkPoca.doca();
607 //double trkFlightArc = track->flightArc(m_hitPosition);
608 //const TrkCircleTraj circleTraj(trkExchangePar);
609 //TrkPoca trkPoca(circleTraj,trkFlightArc,*wire,hitWireLength);
610 //trkFlightArc = trkPoca.flt1();
611 //hitWireLength = trkPoca.flt2();
612 //double doca = trkPoca.doca();
613
614 Hep3Vector wireDirection, trackDirection;
615 //HepPoint3D hitOnWire, hitOntrack;
616 wire->getInfo(hitWireLength, positionOnDetector, wireDirection);
617 helixTraj.getInfo(trkFlightLength, positionOntrack, trackDirection);
618 residual = m_driftDist - fabs(doca);
619 //cout<<"("<<m_layer<<","<<m_wire<<") ";
620 //cout<<setw(12)<<m_driftDist
621 //<<setw(12)<<doca
622 //<<setw(12)<<residual
623 //<<endl;
624
625 //cout<<"("<<m_layer<<","<<m_wire<<") "
626 //<<setw(12)<<getPairHit()->getHitPosition().z()
627 //<<setw(12)<<positionOntrack.z()
628 //<<setw(12)<<positionOnDetector.z()
629 //<<setw(12)<<tanl/sqrt(1+tanl*tanl)*trkFlightLength+z0
630 //<<setw(12)<<tanl*(track->flightArc(m_hitPosition))+z0
631 //<<endl;
632 }
633 else if(m_hitType==CGEMHIT){
634 double rCgem = m_hitPosition.perp();
635 //double phiTrkFlt = m_helix.IntersectCylinder(rCgem);
636 double phiTrkFlt = track->IntersectCylinder(rCgem);
637 phiTrkFlt = track->judgeHalf(this)*phiTrkFlt;
638 //HepPoint3D crossPoint = m_helix.x(phiTrkFlt);
639 positionOntrack = track->x(phiTrkFlt);
640 positionOnDetector = track->x(phiTrkFlt);
641 if(m_flag==0){
642 double phiCrossPoint = positionOntrack.phi();
643 double phiMeasure = m_hitPosition.phi();
644 double dphi = phiMeasure - phiCrossPoint;
645 while(dphi<-M_PI)dphi += 2*M_PI;
646 while(dphi> M_PI)dphi -= 2*M_PI;
647 residual = rCgem*dphi;
648 }
649 else{
650 residual = (m_cgemCluster->getrecv() - m_cgemGeomSvc->getReadoutPlane(m_layer,m_wire)->getVFromPhiZ(positionOntrack.phi(),positionOntrack.z()*10.))/10.;
651 }
652 }
653 return residual;
654}
s
XmlRpcServer s
Definition: HelloServer.cpp:11
sin
double sin(const BesAngle a)
Definition: InstallArea/include/MdcGeom/MdcGeom/BesAngle.h:210
cos
double cos(const BesAngle a)
Definition: InstallArea/include/MdcGeom/MdcGeom/BesAngle.h:213
M_PI
#define M_PI
Definition: TConstant.h:4
CgemGeomSvc::getReadoutPlane
CgemGeoReadoutPlane * getReadoutPlane(int iLayer, int iSheet) const
Definition: Cgem/CgemGeomSvc/CgemGeomSvc-00-00-31-p01/CgemGeomSvc/CgemGeomSvc.h:220
HelixTraj::getInfo
virtual void getInfo(double fltLen, HepPoint3D &pos, Hep3Vector &dir) const
Definition: HelixTraj.cxx:181
Helix::flightLength
double flightLength(HepPoint3D &hit) const
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-12/src/Helix.cxx:812
Helix::center
const HepPoint3D & center(void) const
returns position of helix center(z = 0.);
Definition: InstallArea/include/TrackUtil/TrackUtil/Helix.h:212
Helix::x
HepPoint3D x(double dPhi=0.) const
returns position after rotating angle dPhi in phi direction.
Definition: Reconstruction/TrackUtil/bak_TrackUtil-00-00-11/src/Helix.cxx:132
Helix::radius
double radius(void) const
returns radious of helix.
Definition: InstallArea/include/TrackUtil/TrackUtil/Helix.h:224
Helix::dr
double dr(void) const
returns an element of parameters.
Definition: InstallArea/include/TrackUtil/TrackUtil/Helix.h:236
Helix::pivot
const HepPoint3D & pivot(void) const
returns pivot position.
Definition: InstallArea/include/TrackUtil/TrackUtil/Helix.h:218
HoughHit::VHitPosition
vector< HepPoint3D > & VHitPosition(HoughTrack *track)
Definition: HoughTransAlg/HoughTransAlg-00-00-14/src/HoughHit.cxx:408
HoughHit::operator=
HoughHit & operator=(const HoughHit &other)
Definition: HoughTransAlg/HoughTransAlg-00-00-14/src/HoughHit.cxx:209
MdcCalibFunSvc::getT0
double getT0(int layid, int cellid) const
Definition: MdcCalibFunSvc.cxx:507
MdcCalibFunSvc::driftTimeToDist
double driftTimeToDist(double drifttime, int layid, int cellid, int lr, double entrance=0.0) const
Definition: MdcCalibFunSvc.cxx:149
MdcCalibFunSvc::getTprop
double getTprop(int lay, double z) const
Definition: MdcCalibFunSvc.cxx:143
MdcCalibFunSvc::getTimeWalk
double getTimeWalk(int layid, double Q) const
Definition: MdcCalibFunSvc.cxx:514
MdcGeomSvc::Wire
const MdcGeoWire *const Wire(unsigned id)
Definition: MdcGeomSvc.cxx:770
MdcGeomSvc::Layer
const MdcGeoLayer *const Layer(unsigned id)
Definition: MdcGeomSvc.cxx:786
MdcHit::hitTraj
const Trajectory * hitTraj() const
Definition: MdcHit.cxx:231
MdcID::layer
static int layer(const Identifier &id)
Values of different levels (failure returns 0)
Definition: bak_Identifier-01-02-14/src/MdcID.cxx:49
MdcID::wire
static int wire(const Identifier &id)
Definition: bak_Identifier-01-02-14/src/MdcID.cxx:54
RawData::identify
virtual Identifier identify() const
Definition: RawData.cxx:15
RawData::getChargeChannel
unsigned int getChargeChannel() const
Definition: RawData.cxx:45
RawData::getTimeChannel
unsigned int getTimeChannel() const
Definition: RawData.cxx:40
c1
TCanvas * c1
Definition: tau_mode.c:75