00001
00002
00003 #include "MdcTrkRecon/MdcSeg.h"
00004 #include <stdlib.h>
00005 #include <math.h>
00006 #include "MdcGeom/BesAngle.h"
00007 #include "MdcTrkRecon/mdcWrapAng.h"
00008 #include "MdcTrkRecon/mdcWrapWire.h"
00009 #include "MdcTrkRecon/MdcLine.h"
00010 #include "MdcTrkRecon/MdcSegParams.h"
00011 #include "MdcData/MdcHit.h"
00012 #include "MdcGeom/MdcSuperLayer.h"
00013 #include "MdcGeom/MdcLayer.h"
00014 #include "MdcTrkRecon/MdcSegInfoSterO.h"
00015 #include "MdcTrkRecon/MdcSegUsage.h"
00016 #include "MdcTrkRecon/MdcMap.h"
00017 #include "MdcData/MdcHitUse.h"
00018 #include "MdcData/MdcHitMap.h"
00019
00020
00021 #include "AIDA/IHistogram1D.h"
00022 extern AIDA::IHistogram1D* g_nSigAdd;
00023
00024
00025
00026 extern int haveDigiTk[43][288];
00027 extern double haveDigiPt[43][288];
00028 extern double haveDigiTheta[43][288];
00029 extern double haveDigiPhi[43][288];
00030 extern int haveDigiAmbig[43][288];
00031
00032 MdcSegParams* MdcSeg::segParam = 0;
00033 const double twoPi = Constants::twoPi;
00034
00035 MdcSeg::MdcSeg(double bt) {
00036
00037 _info = 0;
00038 _bunchTime = bt;
00039 }
00040
00041
00042 MdcSeg::~MdcSeg() {
00043
00044 if (_info != 0) delete _info;
00045 reset();
00046 }
00047
00048
00049 MdcSeg::MdcSeg(const MdcSeg& other):
00050 GmsListLink(), _slayer(other._slayer), _phi(other._phi), _slope(other._slope), _chisq(other._chisq), _qual(other._qual), _pattern(other._pattern), _info(other._info), _bunchTime(other._bunchTime)
00051
00052 {
00053 HepAListDeleteAll(_theList);
00054 for(int i=0; i<other.nHit(); i++){
00055 _theList.append(other.hit(i));
00056 }
00057 for(int j=0; j<3; j++){
00058 _errmat[0] = other._errmat[0];
00059 _errmat[1] = other._errmat[1];
00060 _errmat[2] = other._errmat[2];
00061 }
00062 segParam = other.segParam;
00063 }
00064
00065
00066 MdcSeg&
00067 MdcSeg::operator = (const MdcSeg& other) {
00068
00069 if(&other != this){
00070
00071 HepAListDeleteAll(_theList);
00072 for(int i=0; i<other.nHit(); i++){
00073 _theList.append(other.hit(i));
00074 }
00075 _slayer = other._slayer;
00076 _phi = other._phi;
00077 _slope= other._slope;
00078 _errmat[0] = other._errmat[0];
00079 _errmat[1] = other._errmat[1];
00080 _errmat[2] = other._errmat[2];
00081 _chisq = other._chisq;
00082 _qual = other._qual;
00083 _pattern = other._pattern;
00084 _info = other._info;
00085 _bunchTime = other._bunchTime;
00086 segParam = other.segParam;
00087 }
00088
00089 return *this;
00090 }
00091
00092
00093
00094
00095 void
00096 MdcSeg::setInfo(MdcSegInfo *newInfo) {
00097
00098 delete _info;
00099 _info = newInfo;
00100 }
00101
00102
00103 void
00104 MdcSeg::setValues(int nInPatt, int nhit, MdcHit *hits[],
00105 MdcLine *span, const MdcSuperLayer *slay, int ambig[]) {
00106
00107 _qual = 0;
00108 if (nInPatt == 4) _qual |= segFullFlag;
00109 _phi = BesAngle(span->intercept);
00110 _slope = span->slope;
00111 _chisq = span->chisq;
00112 _errmat[0] = span->errmat[0];
00113 _errmat[1] = span->errmat[1];
00114 _errmat[2] = span->errmat[2];
00115 reset();
00116 _slayer = slay;
00117 for (int i = 0; i < nhit; i++) {
00118 MdcHitUse *alink = new MdcHitUse(*(hits[i]), superlayer()->rad0(),
00119 ambig[i]);
00120 append(alink);
00121 }
00122
00123 return;
00124 }
00125
00126
00127 void
00128 MdcSeg::setValues(int nInPatt, double inPhi, double inSlope,
00129 double chi2, double inError[3], const MdcSuperLayer *slay) {
00130
00131
00132 _qual = 0;
00133 if (nInPatt == 4) _qual |= segFullFlag;
00134 _phi = inPhi;
00135 _slope = inSlope;
00136 _chisq = chi2;
00137 _errmat[0] = inError[0];
00138 _errmat[1] = inError[1];
00139 _errmat[2] = inError[2];
00140 _slayer = slay;
00141 reset();
00142
00143 return;
00144 }
00145
00146
00147 void
00148 MdcSeg::markHits(const MdcMap<const MdcHit*,MdcSegUsage*>& usedHits) const {
00149
00150 for (int i = 0; i < nHit(); i++) {
00151 MdcHitUse *alink = _theList[i];
00152 MdcSegUsage *x ;
00153 if ( usedHits.get( alink->mdcHit() , x) ) x->setUsedAmbig( alink->ambig() );
00154 }
00155 }
00156
00157
00158 void
00159 MdcSeg::plotSegAll() const {
00160
00161
00162
00163 for (int ihit=0 ; ihit< nHit() ; ihit++){
00164 hit(ihit)->mdcHit()->print(std::cout);
00165 std::cout << setw(2)<<hit(ihit)->ambig()<<" ";
00166 }
00167
00168 cout<<setiosflags(ios::fixed);
00169 if (info()!=NULL){
00170 std::cout<< " phi " << setprecision(2) << phi()
00171 << " slope " <<std::setw(2)<< setprecision(2) <<slope()<<" ";
00172 if(superlayer()->whichView()==0){
00173 std::cout <<setprecision (2) <<"phi0="<<info()->par(0);
00174 cout<<setprecision(5)<<" cpa="<<info()->par(1);
00175 }else{
00176 std::cout <<setprecision(2)<<"z0="<<info()->par(0)
00177 <<setprecision(2)<<" ct="<<info()->par(1);
00178 }
00179 if(fabs(info()->arc())>0.0001){
00180 std::cout<<setprecision(2)<<" arc="<<info()->arc();
00181 }
00182 std::cout<<setprecision(3)<<" chi2="<<_chisq;
00183 }else{
00184 std::cout<< " phi " << setprecision(2) << phi()
00185 << " slope " <<std::setw(2)<< setprecision(2) <<slope()
00186 << " chi2 "<<setprecision(3) <<chisq();
00187 }
00188
00189 cout<<setprecision(6);
00190 cout<<setiosflags(ios::scientific);
00191 }
00192
00193 void
00194 MdcSeg::plotSeg() const {
00195
00196 std::cout<<superlayer()->slayNum()<<" pat"<<segPattern()<<" ";
00197 for (int ihit=0 ; ihit< nHit() ; ihit++){
00198 hit(ihit)->mdcHit()->print(std::cout);
00199 std::cout <<hit(ihit)->ambig()<<" ";
00200 }
00201 if (info()!=NULL){
00202 cout<< " . ";
00203 }
00204
00205 }
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229 int
00230 MdcSeg::addHits(MdcLine *span, MdcHit** , const MdcHitMap* map,
00231 double corr) {
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242 int cell[2], ambig[2];
00243 double phiwire[2];
00244 int cellused[4] = {0};
00245 int lAdded = 0;
00246 int nhits = nHit();
00247 int m_debug = 0;
00248
00249
00250 int firstnum = superlayer()->firstLayer()->layNum();
00251 for (int i = 0; i < nHit(); i++) {
00252 const MdcHit* dHit = hit(i)->mdcHit();
00253 int laynum = dHit->layernumber();
00254 cellused[laynum - firstnum ] = dHit->wirenumber();
00255 }
00256
00257
00258
00259 for (int layIndex = 0; layIndex < superlayer()->nLayers(); layIndex++) {
00260 const MdcLayer *layer = superlayer()->layer(layIndex);
00261 int laynum = layer->layNum();
00262
00263
00264 double rinv = 1. / layer->rMid();
00265 double ncellinv = 1. / (double) layer->nWires();
00266 double phiproj = phi() + (layer->rMid() - superlayer()->rad0()) * slope();
00267
00268 BesAngle tmp(phiproj - layer->phiOffset());
00269 cell[0] = (int) floor(layer->nWires() *
00270 tmp.rad() / twoPi + 0.5);
00271 cell[0] = mdcWrapWire( cell[0], layer->nWires() );
00272 phiwire[0] = mdcWrapAng( phi(), cell[0] * twoPi * ncellinv +
00273 layer->phiOffset() );
00274
00275 ambig[0] = (phiwire[0] < phiproj) ? 1 : -1;
00276
00277
00278 ambig[1] = -ambig[0];
00279 cell[1] = mdcWrapWire( cell[0] + ambig[0], layer->nWires() );
00280 phiwire[1] = mdcWrapAng( phi(), cell[1] * twoPi * ncellinv +
00281 layer->phiOffset() );
00282
00283 if(m_debug) std::cout<< " loop over the two possible wires " << std::endl;
00284
00285 for (int iroad = 0; iroad < 2; iroad++) {
00286 if (cellused[laynum - firstnum] == cell[iroad]) continue;
00287 if(m_debug) std::cout<< "possible wires "<<laynum<<" "<<cell[iroad]<<std::endl;
00288 if (map->hitWire(laynum, cell[iroad]) != 0) {
00289 MdcHit *ahit = map->hitWire(laynum, cell[iroad]);
00290
00291 if (ahit->usedHit()) {
00292 if(m_debug) std::cout<< "hit used continue " <<std::endl;
00293 continue;
00294 }
00295
00296
00297 double phidrift = ahit->driftDist(bunchTime(), ambig[iroad]) * corr * rinv;
00298 double phihit = mdcWrapAng( phi(), phidrift * ambig[iroad] + ahit->phi());
00299
00300
00301 double sigphi = corr * ahit->sigma(bunchTime(), 0) * rinv;
00302
00303 if ( g_nSigAdd && fabs(sigphi)>0.0001 ) { g_nSigAdd->fill(fabs(phihit - phiproj) / sigphi); }
00304 if ( fabs(phihit - phiproj) > sigphi * segParam->nsigAddHit ) {
00305 if(m_debug) std::cout<< fabs(phihit-phiproj) <<"> add hit sigma "
00306 << sigphi<< "*"<< segParam->nsigAddHit <<"="<<sigphi*segParam->nsigAddHit<<std::endl;
00307 continue;
00308 }
00309
00310
00311 lAdded = 1;
00312 span->sigma[nhits] = sigphi;
00313 span->x[nhits] = layer->rMid() - superlayer()->rad0();
00314 span->y[nhits] = mdcWrapAng(span->y[0], phihit);
00315
00316
00317 MdcHitUse *alink = new MdcHitUse(*ahit, superlayer()->rad0(),
00318 ambig[iroad]);
00319 append(alink);
00320
00321
00322 nhits++;
00323 }
00324
00325 }
00326
00327 }
00328
00329 if (lAdded) {
00330 span->fit(nhits);
00331 BesAngle tmpAngle(span->intercept);
00332 span->intercept = tmpAngle;
00333 _phi = span->intercept;
00334 _slope = span->slope;
00335 _chisq = span->chisq;
00336 _errmat[0] = span->errmat[0];
00337 _errmat[1] = span->errmat[1];
00338 _errmat[2] = span->errmat[2];
00339 }
00340
00341 return nhits;
00342 }
00343
00344
00345 void
00346 MdcSeg::reset() {
00347
00348 HepAListDeleteAll( _theList );
00349 }
00350
00351
00352 void
00353 MdcSeg::append(MdcHitUse* theHitUse) {
00354
00355 _theList.append(theHitUse);
00356 }
00357
00358
00359 void
00360 MdcSeg::remove(MdcHitUse* theHitUse) {
00361
00362 _theList.remove(theHitUse);
00363 delete theHitUse;
00364 }
00365
00366
00367 int
00368 MdcSeg::nHit() const {
00369
00370 return _theList.length();
00371 }
00372
00373
00374
00375 double
00376 MdcSeg::testCombSeg(const MdcSeg* testSeg)const{
00377
00378 int tkId= -1;
00379 for (int i=0; i<nHit(); i++){
00380 const MdcHit* h = hit(i)->mdcHit();
00381 unsigned int l = h->layernumber();
00382 unsigned int w = h->wirenumber();
00383
00384 if ( haveDigiTk[l][w]<0 || haveDigiTk[l][w]>100 ) continue;
00385 if (tkId<0){
00386 tkId = haveDigiTk[l][w];
00387 }else if (haveDigiTk[l][w] != tkId){
00388 return -1;
00389 }
00390 }
00391
00392 double nSame = 0.;
00393 for (int i=0; i<testSeg->nHit(); i++){
00394 const MdcHit* h = testSeg->hit(i)->mdcHit();
00395 unsigned int l = h->layernumber();
00396 unsigned int w = h->wirenumber();
00397 if (haveDigiTk[l][w] == tkId){
00398 ++nSame;
00399 }
00400 }
00401
00402 return nSame/testSeg->nHit();
00403 }
00404
00405
00406 double
00407 MdcSeg::testCombSegPt()const{
00408
00409 double truthPt = -1.;
00410 for (int i=0; i<nHit(); i++){
00411 const MdcHit* h = hit(i)->mdcHit();
00412 unsigned int l = h->layernumber();
00413 unsigned int w = h->wirenumber();
00414 if ( haveDigiPt[l][w]<0 ) continue;
00415
00416 if (truthPt<0){ truthPt = haveDigiPt[l][w]; }
00417 }
00418
00419 return truthPt;
00420 }
00421
00422
00423 double
00424 MdcSeg::testCombSegTheta()const{
00425
00426 double truthTheta = -999.;
00427 for (int i=0; i<nHit(); i++){
00428 const MdcHit* h = hit(i)->mdcHit();
00429 unsigned int l = h->layernumber();
00430 unsigned int w = h->wirenumber();
00431 if ( haveDigiTheta[l][w]<-998. ) continue;
00432
00433 if (truthTheta<-998.){ truthTheta = haveDigiTheta[l][w]; }
00434 }
00435
00436 return truthTheta;
00437 }
00438
00439
00440 double
00441 MdcSeg::testCombSegPhi()const{
00442
00443 double truthPhi = -999.;
00444 for (int i=0; i<nHit(); i++){
00445 const MdcHit* h = hit(i)->mdcHit();
00446 unsigned int l = h->layernumber();
00447 unsigned int w = h->wirenumber();
00448 if ( haveDigiPhi[l][w]<-998. ) continue;
00449
00450 if (truthPhi<-998.){ truthPhi = haveDigiPhi[l][w]; }
00451 }
00452
00453 return truthPhi;
00454 }
00455
00456
00457 double
00458 MdcSeg::testCombSegAmbig()const{
00459
00460 double ambigOk = 0.;
00461 for (int i=0; i<nHit(); i++){
00462 const MdcHit* h = hit(i)->mdcHit();
00463 unsigned int l = h->layernumber();
00464 unsigned int w = h->wirenumber();
00465 if( hit(i)->ambig() == haveDigiAmbig[l][w] ) ambigOk++;
00466 }
00467
00468 return ambigOk/nHit();
00469 }
00470