CGEM BOSS 6.6.5.f
BESIII Offline Software System
Loading...
Searching...
No Matches
MdcHoughFinder/MdcHoughFinder-00-00-12/src/HoughHit.cxx
Go to the documentation of this file.
1#include "MdcHoughFinder/HoughHit.h"
2#include "MdcHoughFinder/HoughGlobal.h"
3#include "RawEvent/RawDataUtil.h"
4#include "Identifier/MdcID.h"
5#include "MdcGeomSvc/IMdcGeomSvc.h"
6#include "MdcGeomSvc/MdcGeoWire.h"
7#include "MdcGeomSvc/MdcGeoLayer.h"
8#include "GaudiKernel/SvcFactory.h"
9#include "GaudiKernel/ISvcLocator.h"
10#include "GaudiKernel/StatusCode.h"
11#include <assert.h>
12
13const MdcCalibFunSvc* HoughHit::_calibPtr=NULL;
14MdcGeomSvc* HoughHit::_mdcGeomSvc=NULL;
15CgemGeomSvc* HoughHit::_cgemGeomSvc=NULL;
16double HoughHit::_bunchTime;
17int HoughHit::_npart=360;
18//MdcDetector* HoughHit::_det=NULL;
19
20HoughHit::HoughHit(){
21 _det=Global::m_gm;
22 _digiPtr=NULL;
23 _layerPtr=NULL;
24 _wirePtr=NULL;
25 _clusterPtr=NULL;
26 _hitMap = NULL;
27}
28HoughHit::HoughHit(const MdcDigi* const aDigi){
29 _det= Global::m_gm;
30 _digiPtr=aDigi;
31 _id = aDigi->identify();
32 _layer = MdcID::layer(_id);
33 _wire = MdcID::wire(_id);
34 _layerPtr = _det->Layer(_id);
35 _wirePtr = _det->Wire(_id);
36 _rawTime = RawDataUtil::MdcTime(_digiPtr->getTimeChannel());
37 _charge = _digiPtr->getChargeChannel();
38 _driftTime = driftTime();
39 _driftDist = calDriftDist(_bunchTime,0);
40 _slayerType=slayerType(_layer);
41 _lr = 0;
42 _rmid = _wirePtr->rMid();
43
44 assert(_mdcGeomSvc!=NULL);
45 const MdcGeoWire* wire = _mdcGeomSvc->Wire(_layer,_wire);
46 assert(wire!=NULL);
47 HepPoint3D eastP = wire->Backward()/10.;
48 HepPoint3D westP = wire->Forward()/10.;
49 _eastPoint = eastP;
50 _westPoint = westP;
51 _midPoint = (eastP+westP)/2.;
52
53 _type = MIDPOINT;
54 conformalTrans( _midPoint.x(), _midPoint.y(),_driftDist);
55
56 _detectorType = MDC;
57 _clusterPtr = NULL;
58 _trkid = aDigi->getTrackIndex();
59 _clusterid = -999;
60 _hitMap = NULL;
61 _used = 0;
62}
63
64HoughHit::HoughHit(const RecCgemCluster* const cluster){
65 _det= Global::m_gm;
66 _digiPtr=NULL;
67 //_id = aDigi->identify(); //FIXME
68 _layer = cluster->getlayerid();
69 _wire = cluster->getsheetid(); //FIXME
70 //_wire = cluster->getclusterid(); //FIXME
71 _layerPtr = NULL;
72 _wirePtr = NULL;
73 _rawTime = 0; //FIXME
74 _charge = cluster->getenergydeposit(); //FIXME
75 _driftTime = 0; //FIXME
76 _driftDist = 0.00001; //FIXME
77 _slayerType= cluster->getflag(); //FIXME
78 _lr = 0;
79 _rmid = _cgemGeomSvc->getCgemLayer(_layer)->getMiddleROfGapD()/10.;//mm 2 cm; //FIXME
80
81 assert(_cgemGeomSvc!=NULL);
82 //const MdcGeoWire* wire = _mdcGeomSvc->Wire(_layer,_wire);
83 //assert(wire!=NULL);
84 //HepPoint3D point(0,0,0); //FIXME
85 HepPoint3D point((_cgemGeomSvc->getCgemLayer(_layer)->getMiddleROfGapD()/10.)*cos(cluster->getrecphi()),(_cgemGeomSvc->getCgemLayer(_layer)->getMiddleROfGapD()/10.)*sin(cluster->getrecphi()),cluster->getRecZ()/10.); //FIXME
86 //if(_slayerType==2)cout<<"cluster "<<cluster->getclusterid()<<" layer "<<_layer<<" phi,z,V"<<cluster->getrecphi()<<" ,"<<cluster->getRecZ()<<" ,"<<cluster->getrecv()<<endl;
87 HepPoint3D eastP = point;
88 HepPoint3D westP = point;
89 _eastPoint = eastP;
90 _westPoint = westP;
91 _midPoint = (eastP+westP)/2.;
92
93 _type = MIDPOINT;
94 conformalTrans( _midPoint.x(), _midPoint.y(),_driftDist);
95 //conformalTrans( (_cgemGeomSvc->getCgemLayer(_layer)->getMiddleROfGapD()/10.)*cos(cluster->getrecphi()), (_cgemGeomSvc->getCgemLayer(_layer)->getMiddleROfGapD()/10.)*sin(cluster->getrecphi()),_driftDist);
96
97 _detectorType = CGEM;
98 _clusterPtr = cluster;
99 _trkid = 999;
100 _clusterid = cluster->getclusterid();
101 _hitMap = NULL;
102 _used = 0;
103}
104
105HoughHit::HoughHit(const HoughHit& other) :
106 _digiPtr(other._digiPtr),
107 _det(other._det),
108 _layerPtr(other._layerPtr), _wirePtr(other._wirePtr), _id(other._id),
109 _layer(other._layer), _wire(other._wire),
110 _rawTime(other._rawTime), _charge(other._charge),
111 _eastPoint(other._eastPoint),_westPoint(other._westPoint),_midPoint(other._midPoint),
112 _leftPoint(other._leftPoint),_rightPoint(other._rightPoint),
113 _u(other._u),_v(other._v),_lr(other._lr),
114 _type(other._type),_detectorType(other._detectorType),
115 _rmid(other._rmid), _slayerType(other._slayerType),
116 _cirlist(other._cirlist),_style(other._style),
117 _driftTime(other._driftTime), _driftDist(other._driftDist),
118 CF_drift(other.CF_drift),
119 vec_cfcir(other.vec_cfcir),
120 //truth
121 _truthU(other._truthU),
122 _truthV(other._truthV),
123 _truthR(other._truthR),
124 _truthlr(other._truthlr),
125 _truthId(other._truthId),
126 _truthPoint(other._truthPoint),
127 _truthDrift(other._truthDrift),
128 _deltad(other._deltad),
129 _flightLength(other._flightLength),
130 _clusterPtr(other._clusterPtr),
131 _trkid(other._trkid),
132 _hitid(other._hitid),
133 _clusterid(other._clusterid),
134 _hitMap(other._hitMap),
135 _used(other._used)
136{ }
137
138HoughHit& HoughHit::operator=(const HoughHit& other){
139 _digiPtr=(other._digiPtr);
140 _det=(other._det);
141 _layerPtr=(other._layerPtr); _wirePtr=(other._wirePtr); _id=(other._id);
142 _layer=(other._layer); _wire=(other._wire);
143 _rawTime=(other._rawTime); _charge=(other._charge);
144 _eastPoint=(other._eastPoint);_westPoint=(other._westPoint);_midPoint=(other._midPoint);
145 _leftPoint=(other._leftPoint);_rightPoint=(other._rightPoint);
146 _u=(other._u);_v=(other._v);_lr=(other._lr);
147 _type=(other._type);_detectorType=(other._detectorType);
148 _rmid=(other._rmid); _slayerType=(other._slayerType);
149 _cirlist=(other._cirlist);_style=(other._style);
150 _driftTime=(other._driftTime); _driftDist=(other._driftDist);
151 CF_drift = (other.CF_drift);
152 vec_cfcir = (other.vec_cfcir);
153 //truth
154 _truthU= (other._truthU);
155 _truthV= (other._truthV);
156 _truthR= (other._truthR);
157 _truthlr= (other._truthlr);
158 _truthId= (other._truthId);
159 _truthDrift= (other._truthDrift);
160 _truthPoint= (other._truthPoint);
161 _deltad= (other._deltad);
162 _flightLength= (other._flightLength);
163 _clusterPtr= (other._clusterPtr);
164 _trkid = (other._trkid);
165 _hitid = (other._hitid);
166 _clusterid = (other._clusterid);
167 _hitMap = (other._hitMap);
168 _used = (other._used);
169}
170
171void HoughHit::setTruthInfo(const MdcMcHit*& aMcHit){
172 if(!aMcHit) return;
173 _truthDrift = aMcHit->getDriftDistance()/10.;//mm to cm
174 _truthPoint.setX(aMcHit->getPositionX()/10.);//mm to cm
175 _truthPoint.setY(aMcHit->getPositionY()/10.);
176 _truthPoint.setZ(aMcHit->getPositionZ()/10.);
177
178 int mcLR = aMcHit->getPositionFlag();
179 // if (mcLR == 0) mcLR = -1;//FIXME
180 // for truth reserve from digi
181 if (mcLR == 1) mcLR = -1;//
182 if (mcLR == 0) mcLR = 1;//
183 _truthId = aMcHit->getTrackIndex();
184 _truthlr = mcLR;
185
186 //fix
187 _truthU = getConformal_u( _truthPoint.x(), _truthPoint.y(), _truthDrift );
188 _truthV = getConformal_v( _truthPoint.y(), _truthPoint.x(), _truthDrift );
189 _truthR = getConformal_r( _truthPoint.x(), _truthPoint.y(), _truthDrift );
190
191 //std::cout<<__FILE__<<" "<<_layer<<","<<_wire<<" "<<_truthId<<" truth "<<_truthPoint<<" truthU "<<_truthU<<std::endl;
192
193 _type = DIGIWITHTRUTH;
194}
195
196void HoughHit::setTruthInfo(const CgemMcHit* aMcHit){
197 if(!aMcHit) return;
198 _truthDrift = 0.00001;//mm to cm
199 _truthPoint.setX((aMcHit->GetPositionXOfPrePoint()+aMcHit->GetPositionXOfPostPoint())/2/10.);//mm to cm
200 _truthPoint.setY((aMcHit->GetPositionYOfPrePoint()+aMcHit->GetPositionYOfPostPoint())/2/10.);
201 _truthPoint.setZ((aMcHit->GetPositionZOfPrePoint()+aMcHit->GetPositionZOfPostPoint())/2/10.);
202
203 //int mcLR = aMcHit->getPositionFlag();
204 // if (mcLR == 0) mcLR = -1;//FIXME
205 // for truth reserve from digi
206 //if (mcLR == 1) mcLR = -1;//
207 //if (mcLR == 0) mcLR = 1;//
208 _truthId = aMcHit->GetTrackID();
209 _truthlr = 0;
210
211 //fix
212 _truthU = getConformal_u( _truthPoint.x(), _truthPoint.y(), _truthDrift );
213 _truthV = getConformal_v( _truthPoint.y(), _truthPoint.x(), _truthDrift );
214 _truthR = getConformal_r( _truthPoint.x(), _truthPoint.y(), _truthDrift );
215
216 //std::cout<<__FILE__<<" "<<_layer<<","<<_wire<<" "<<_truthId<<" truth "<<_truthPoint<<" truthU "<<_truthU<<std::endl;
217
218 _type = CLUSTERWITHTRUTH;
219}
220
221void HoughHit::conformalTrans(double x,double y,double r){
222 _u = 2*x/(x*x+y*y-r*r);
223 _v = 2*y/(x*x+y*y-r*r);
224 CF_drift = 2*r/(x*x+y*y-r*r);
225}
226
227double HoughHit::getConformal_u(double x,double y,double r) {
228 return 2*x/(x*x+y*y-r*r);
229}
230double HoughHit::getConformal_v(double x,double y,double r) {
231 return 2*x/(x*x+y*y-r*r);
232}
233double HoughHit::getConformal_r(double x,double y,double r) {
234 return 2*r/(x*x+y*y-r*r);
235}
236
237double HoughHit::driftTime() const {
238 double tprop = _calibPtr->getTprop(_layer,0);
239 double T0Walk = _calibPtr->getT0(_layer,_wire) + _calibPtr->getTimeWalk(_layer, _charge);
240 //tof in ns, driftTime in ns, _T0Walk in ns
241 double driftT = _rawTime - T0Walk - 1.e9*_bunchTime- tprop;
242 return driftT;
243}
244double HoughHit::driftTime(double tof, double z) const {
245 double tprop = _calibPtr->getTprop(_layer,z*10.);
246 double T0Walk = _calibPtr->getT0(_layer,_wire) + _calibPtr->getTimeWalk(_layer, _charge);
247 //tof in ns, driftTime in ns, _T0Walk in ns
248 double driftT = _rawTime - T0Walk - 1.e9*tof - tprop;
249
250 return driftT;
251}
252
253double HoughHit::calDriftDist(double bunchTime, int ambig) const {
254
255 // double crudeTof = 0; //FIXME
256 return calDriftDist(bunchTime+crudeTof(), ambig, 0., 0., 0. );
257}
258
259double HoughHit::calDriftDist(double tof, int ambig, double entranceAngle, double /*dipAngle*/, double z) const {
260
261 double driftD;
262 //drift time ns, layer id begin with 0, entrance angle rads,
263 //lr ambig: wire ambig 1,-1,0 -> Calib 0,1,2
264 int lrCalib=2;
265 if (ambig==1) lrCalib = 0;
266 else if (ambig==-1) lrCalib = 1;
267
268 // tof in s, driftDist in cm, dirftTime in ns
269 if (fabs(z)>150. || fabs(driftTime(tof,z))>1500.){
270 return 9999.;
271 }
272 driftD = 0.1 * _calibPtr->driftTimeToDist(driftTime(tof,z),_layer,_wire,lrCalib,entranceAngle);//to cm
273 //std::cout<<"driftDist "<<"("<<_layer <<","<<_wire <<") dd "<<driftD<<" dt "<<driftTime(tof,z) <<" lr "<<lrCalib <<" eAng "<<entranceAngle <<" tof "<<tof*1.e9<<" z "<<z <<" t0walk "<<_T0Walk<<" rawT "<<_rawTime <<" tprop "<< _rawTime - driftTime(tof,z)- _T0Walk-1.e9*tof<<std::endl;
274
275 if (abs(driftD)<0.00001) driftD = 0.00001;
276 return driftD;
277}
278
279void HoughHit::print() const{
280 std::cout<<"("<<_layer<<","<<_wire<<") "<<std::endl;
281}
282
283void HoughHit::printAll() const{
284 //std::cout<<"("<<_layer<<","<<_wire<<") id "<<this->digi()->getTrackIndex()<<" xyz "<<_midPoint<<endl;
285 if(_detectorType==MDC){
286 std::cout<<"("<<_layer<<","<<_wire<<") id "<<this->digi()->getTrackIndex()<<" xyz "<<_midPoint<<endl;
287 }else{
288 cout <<"cgem "
289 <<_clusterPtr->getclusterid()<<" layer "<<_clusterPtr->getlayerid()
290 <<" uv "<<_u<<" "<<_v
291 <<" flag "<<(_clusterPtr)->getflag()
292 <<" recphi "<<(_clusterPtr)->getrecphi()<<" recv "<<(_clusterPtr)->getrecv()<<" recz "<<(_clusterPtr)->getRecZ()
293 <<" sheet "<<(_clusterPtr)->getsheetid()
294 <<" energy "<<(_clusterPtr)->getenergydeposit()
295 <<endl;
296 }
297}
298void HoughHit::printTruth() const{
299 std::cout<<"("<<_layer<<","<<_wire<<") truth "<<_truthId<<" xyz "<<_truthPoint<<" uv "<<_truthU<<","<<_truthV<<") "<<" cir list "<<_cirlist<<" ambig "<<_truthlr<<std::endl;
300}
301
302int HoughHit::slayerType(int layer){
303 //get nominal shift cell of this layer
304 double nomShift= _mdcGeomSvc->Layer(layer)->nomShift();
305
306 if(nomShift>0) return 1;
307 else if(nomShift<0) return -1;
308 else return 0;
309
310 return -999;
311}
312
313//calcu drift cir
314void HoughHit::makeCir(int n,double phi_begin,double phi_last,double r){
315 vec_cfcir.clear();
316 vector<CFCir>().swap(vec_cfcir);
317 for(int i =0; i<n; i++){
318 double phi_slice = (phi_last-phi_begin)/n;
319 double phi = phi_begin + (1/2.+i)*phi_slice;
320 double x = _u + r*cos(phi);
321 double y = _v + r*sin(phi);
322 vec_cfcir.push_back( CFCir(x,y,phi,n,_u,_v,r) );
323 }
324}
325
326void HoughHit::buildMap(int x_bin, double x_min, double x_max, int y_bin, double y_min, double y_max, int nPoint, int charge)
327{
328 if(_hitMap!=NULL)delete _hitMap;
329 stringstream ssname;
330 ssname<<"id"<<_hitid<<",layer"<<_layer<<",wire(sheet)"<<_wire;
331 string sname = ssname.str();
332 const char * name = sname.c_str() ;
333 //cout<<name<<endl;
334 _hitMap = new TH2D(name,name, x_bin, x_min, x_max, y_bin, y_min, y_max);
335 int N = x_bin*nPoint;
336 double delta_alpha = (x_max - x_min )/N;
337 double alpha = x_min-delta_alpha/2;
338 for(int i=0;i<N;i++){
339 alpha += delta_alpha;
340 if(alpha>M_PI)alpha -= M_PI;
341 double rho1 = _u*cos(alpha) + _v*sin(alpha) + CF_drift;
342 double rho2 = _u*cos(alpha) + _v*sin(alpha) - CF_drift;
343 //cout<<alpha<<" "<<rho1<<" "<<rho2<<endl;
344 double slantOfLine = _v*cos(alpha) - _u*sin(alpha);
345 int chargeOfHitOnCir1 = (slantOfLine/fabs(slantOfLine))*(rho1/fabs(rho1));
346 int chargeOfHitOnCir2 = (slantOfLine/fabs(slantOfLine))*(rho2/fabs(rho2));
347 if(chargeOfHitOnCir1 != charge && y_min<(rho1) && (rho1)<y_max)
348 //if( y_min<(rho1) && (rho1)<y_max)
349 _hitMap->Fill(alpha,rho1);
350 if(chargeOfHitOnCir2 != charge && y_min<(rho2) && (rho2)<y_max)
351 //if( y_min<(rho2) && (rho2)<y_max)
352 _hitMap->Fill(alpha,rho2);
353 }
354}
355
356void HoughHit::clearMap()
357{
358 if(_hitMap!=NULL)delete _hitMap;
359 _hitMap = NULL;
360}
n
const Int_t n
Definition: DataBase/tau_mode.c:65
x
Double_t x[10]
Definition: DataBase/tau_mode.c:57
abs
double abs(const EvtComplex &c)
Definition: EvtComplex.hh:212
alpha
const double alpha
Definition: FastVertexFit.cxx:4
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
HoughHit::conformalTrans
void conformalTrans(double x, double y, double r)
Definition: MdcHoughFinder/MdcHoughFinder-00-00-12/src/HoughHit.cxx:221
HoughHit::makeCir
void makeCir(int n, double phi_begin, double phi_last, double r)
Definition: MdcHoughFinder/MdcHoughFinder-00-00-12/src/HoughHit.cxx:314
HoughHit::setTruthInfo
void setTruthInfo(const MdcMcHit *&mcHit)
HoughHit::getConformal_v
double getConformal_v(double, double, double)
Definition: MdcHoughFinder/MdcHoughFinder-00-00-12/src/HoughHit.cxx:230
HoughHit::calDriftDist
double calDriftDist(double, int, double, double, double) const
Definition: MdcHoughFinder/MdcHoughFinder-00-00-12/src/HoughHit.cxx:259
HoughHit::buildMap
void buildMap(int x_bin, double x_min, double x_max, int y_bin, double y_min, double y_max, int nPoint, int charge)
Definition: MdcHoughFinder/MdcHoughFinder-00-00-12/src/HoughHit.cxx:326
HoughHit::getConformal_r
double getConformal_r(double, double, double)
Definition: MdcHoughFinder/MdcHoughFinder-00-00-12/src/HoughHit.cxx:233
HoughHit::operator=
HoughHit & operator=(const HoughHit &other)
Definition: HoughTransAlg/HoughTransAlg-00-00-14/src/HoughHit.cxx:209
HoughHit::getConformal_u
double getConformal_u(double, double, double)
Definition: MdcHoughFinder/MdcHoughFinder-00-00-12/src/HoughHit.cxx:227
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
MdcDetector::Layer
const MdcLayer * Layer(unsigned id) const
Definition: InstallArea/include/MdcGeom/MdcGeom/MdcDetector.h:33
MdcDetector::Wire
const MdcSWire * Wire(unsigned id) const
Definition: InstallArea/include/MdcGeom/MdcGeom/MdcDetector.h:28
MdcGeomSvc::Wire
const MdcGeoWire *const Wire(unsigned id)
Definition: MdcGeomSvc.cxx:770
MdcGeomSvc::Layer
const MdcGeoLayer *const Layer(unsigned id)
Definition: MdcGeomSvc.cxx:786
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::getTrackIndex
int getTrackIndex() const
Definition: RawData.cxx:50
RawData::getTimeChannel
unsigned int getTimeChannel() const
Definition: RawData.cxx:40
EvtCyclic3::other
Index other(Index i, Index j)
Definition: EvtCyclic3.cc:118