00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 using namespace std;
00013
00014 #include <iostream>
00015 #include <vector>
00016 #include <CLHEP/Vector/ThreeVector.h>
00017 #include <CLHEP/Geometry/Point3D.h>
00018
00019 #include "Identifier/MucID.h"
00020 #include "MucGeomSvc/MucGeomSvc.h"
00021 #include "MucGeomSvc/MucConstant.h"
00022 #include "MucRecEvent/MucRec2DRoad.h"
00023 #include "MucRecEvent/MucRecLineFit.h"
00024 #include "MucRecEvent/MucRecQuadFit.h"
00025
00026 const double muckInfinity = 1e30;
00027
00028 MucRec2DRoad::MucRec2DRoad(const int& part,
00029 const int& seg,
00030 const int& orient,
00031 const float& xVertex,
00032 const float& yVertex,
00033 const float& zVertex,
00034 const int& fittingMethod)
00035 : m_VertexPos(xVertex, yVertex, zVertex),
00036 m_VertexSigma(0.0, 0.0, 0.0),
00037 m_Part(part), m_Seg(seg), m_Orient(orient),
00038 m_Chi2(0.0), m_DOF(0),
00039 m_MaxHitsPerGap(0), m_LastGap(0), m_TotalHits(0),
00040 m_FitOK(false), m_QuadFitOK(false),
00041 m_HitDistance(5, float(muckInfinity)),
00042 m_pHits(0), m_fittingMethod(fittingMethod)
00043 { }
00044
00045
00046 MucRec2DRoad&
00047 MucRec2DRoad::operator=(const MucRec2DRoad& orig)
00048 {
00049
00050 if ( this != &orig ) {
00051 m_VertexPos = orig.m_VertexPos;
00052 m_VertexSigma = orig.m_VertexSigma;
00053 m_Index = orig.m_Index;
00054 m_Part = orig.m_Part;
00055 m_Seg = orig.m_Seg;
00056 m_Orient = orig.m_Orient;
00057 m_Chi2 = orig.m_Chi2;
00058 m_DOF = orig.m_DOF;
00059 m_FitOK = orig.m_FitOK;
00060 m_MaxHitsPerGap = orig.m_MaxHitsPerGap;
00061 m_LastGap = orig.m_LastGap;
00062 m_TotalHits = orig.m_TotalHits;
00063 m_HitDistance = orig.m_HitDistance;
00064 m_pHits = orig.m_pHits;
00065 m_fittingMethod = orig.m_fittingMethod;
00066 }
00067
00068 return *this;
00069 }
00070
00071
00072 MucRec2DRoad::MucRec2DRoad(const MucRec2DRoad& source)
00073 : m_VertexPos(source.m_VertexPos),
00074 m_VertexSigma(source.m_VertexSigma),
00075 m_Index(source.m_Index),
00076 m_Part(source.m_Part), m_Seg(source.m_Seg), m_Orient(source.m_Orient),
00077 m_Chi2(source.m_Chi2), m_DOF(source.m_DOF),
00078 m_MaxHitsPerGap(source.m_MaxHitsPerGap),
00079 m_LastGap(source.m_LastGap), m_TotalHits(source.m_TotalHits),
00080 m_FitOK(source.m_FitOK),
00081 m_HitDistance(source.m_HitDistance),
00082 m_pHits(source.m_pHits),
00083 m_fittingMethod(source.m_fittingMethod)
00084 { }
00085
00086
00087 MucRec2DRoad::~MucRec2DRoad()
00088 { }
00089
00090
00091 void
00092 MucRec2DRoad::SetIndex(const int& index)
00093 {
00094 if (index >= 0) m_Index = index;
00095 }
00096
00097
00098
00099 void
00100 MucRec2DRoad::AttachHit(MucRecHit* hit)
00101 {
00102
00103 if (!hit) {
00104 cout << "MucRec2DRoad::AttachHit-E1 null hit pointer!" << endl;
00105 return ;
00106 }
00107
00108 int gap = hit->Gap();
00109 if ( (gap < 0) || (gap >= (int)MucID::getGapMax()) ) {
00110
00111 cout << "MucRec2DRoad::AttachHit-W3 bad gap number = " << gap
00112 << endl;
00113 return;
00114 }
00115
00116
00117
00118 for (int i = 0; i < (int)m_pHits.size(); i++) {
00119 if (m_pHits[i] == hit) return;
00120 }
00121 m_pHits.push_back(hit);
00122
00123
00124
00125
00126
00127 CountHits();
00128
00129
00130
00131 float a, sa, b, sb, chi;
00132 int n;
00133 int status = SimpleFit(a, b, sa, sb, chi, n);
00134 if (status < 0) {
00135
00136 }
00137
00138
00139
00140 }
00141
00142
00143
00144 void
00145 MucRec2DRoad::AttachHitNoFit(MucRecHit* hit)
00146 {
00147
00148 if (!hit) {
00149 cout << "MucRec2DRoad::AttachHit-E1 null hit pointer!" << endl;
00150 return ;
00151 }
00152
00153 int gap = hit->Gap();
00154 if ( (gap < 0) || (gap >= (int)MucID::getGapMax()) ) {
00155
00156 cout << "MucRec2DRoad::AttachHit-W3 bad gap number = " << gap
00157 << endl;
00158 return;
00159 }
00160
00161
00162
00163 for (int i = 0; i < (int)m_pHits.size(); i++) {
00164 if (m_pHits[i] == hit) return;
00165 }
00166 m_pHits.push_back(hit);
00167
00168
00169
00170
00171
00172 CountHits();
00173
00174
00175 }
00176
00177
00178
00179 void
00180 MucRec2DRoad::SetMaxNSkippedGaps(const int& numGaps)
00181 {
00182 m_MaxNSkippedGaps = numGaps;
00183 CountHits();
00184 }
00185
00186
00187 int
00188 MucRec2DRoad::SimpleFit(float& slope,
00189 float& intercept,
00190 float& sigmaSlope,
00191 float& sigmaIntercept,
00192 float& chisq,
00193 int& ndof)
00194 {
00195
00196 if(m_pHits.size()>100){
00197 cout<<"MucRec2DRoad: too many hits in this track!"<<endl;
00198 return -1;
00199 }
00200
00201 float px[100];
00202 float py[100];
00203 float pw[100];
00204 int npoints = 0;
00205
00206 vector<MucRecHit*>::const_iterator iHit;
00207
00208 float weight[100];
00209
00210
00211
00212
00213 for (int i = 0; i < m_pHits.size(); i++) {
00214 weight[i] = 1;
00215
00216 for(int j = 0; j < m_pHits.size(); j++){
00217
00218 if(j == i) continue;
00219 if(m_pHits[i]->Part() == m_pHits[j]->Part() &&
00220 m_pHits[i]->Seg() == m_pHits[j]->Seg() &&
00221 m_pHits[i]->Gap() == m_pHits[j]->Gap() )
00222 {
00223 int deltaStrip = fabs(m_pHits[i]->Strip()- m_pHits[j]->Strip());
00224
00225
00226 if(deltaStrip == 0){
00227 cout<<"deltaStrip == 0 ? check it"<<endl;
00228 }
00229 else{
00230 weight[i] *= (deltaStrip+1) * (deltaStrip+1);
00231 }
00232
00233 }
00234 }
00235 }
00236
00237 for (iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++) {
00238 if (*iHit) {
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257 Hep3Vector center = (*iHit)->GetCenterPos();
00258 Hep3Vector sigma = (*iHit)->GetCenterSigma();
00259
00260
00261 if (m_Part == 1) {
00262 if ( m_Orient == 0) {
00263 px[npoints] = center.z();
00264 py[npoints] = sqrt(center.x()*center.x()
00265 + center.y()*center.y());
00266 if(m_Seg==2) py[npoints] = center.y();
00267
00268 pw[npoints] = 4.0 * 1.0/(sigma.y()*sigma.y()) / weight[npoints];
00269 }
00270 else {
00271 px[npoints] = center.x();
00272 py[npoints] = center.y();
00273 pw[npoints] = 4.0 * 1.0/(sigma.x()*sigma.x()) / weight[npoints];
00274 }
00275 }
00276 else {
00277 if ( m_Orient == 0) {
00278 px[npoints] = center.z();
00279 py[npoints] = center.y();
00280 pw[npoints] = 4.0 * 1.0/(sigma.y()*sigma.y()) / weight[npoints];
00281 }
00282 else {
00283 px[npoints] = center.z();
00284 py[npoints] = center.x();
00285 pw[npoints] = 4.0 * 1.0/(sigma.x()*sigma.x()) / weight[npoints];
00286 }
00287 }
00288
00289
00290 npoints++;
00291
00292 if(npoints > 99)
00293 { cout<<"MucRec2DRoad: to many hits in this track, ignore it!"<<endl;
00294 return -1;
00295 }
00296 }
00297 }
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326 ndof = npoints - 2;
00327 if (ndof < 0) return -1;
00328
00329 if (npoints == 1) {
00330 if (m_Part == 1) {
00331 if ( m_Orient == 0) {
00332 px[npoints] = m_VertexPos.z();
00333 py[npoints] = sqrt(m_VertexPos.x()*m_VertexPos.x()
00334 + m_VertexPos.y()*m_VertexPos.y());
00335 if(m_Seg==2) py[npoints] = m_VertexPos.y();
00336 }
00337 else {
00338 px[npoints] = m_VertexPos.x();
00339 py[npoints] = m_VertexPos.y();
00340 }
00341 }
00342 else {
00343 if ( m_Orient == 0) {
00344 px[npoints] = m_VertexPos.z();
00345 py[npoints] = m_VertexPos.y();
00346 }
00347 else {
00348 px[npoints] = m_VertexPos.z();
00349 py[npoints] = m_VertexPos.x();
00350 }
00351 }
00352 pw[npoints] = 1.0;
00353 npoints++;
00354 }
00355 else {
00356 if (npoints == 0 ) {
00357 return -1;
00358 }
00359 }
00360
00361
00362 MucRecLineFit fit;
00363 int status = fit.LineFit(px, py, pw, npoints,
00364 &slope, &intercept, &chisq,
00365 &sigmaSlope, &sigmaIntercept);
00366
00367
00368 float tempslope, tempintercept,tempchisq, tempsigmaslope, sigmaPos;
00369 int status4 = fit.LineFit(px, py, pw,m_Part,m_Seg,m_Orient, npoints,
00370 &tempslope, &tempintercept, &tempchisq,
00371 &tempsigmaslope, &sigmaPos);
00372
00373 MucRecQuadFit quadfit;
00374 float quad_a, quad_b, quad_c, sigmaquad_a, sigmaquad_b, sigmaquad_c, chisq_quad;
00375 int whichhalf, status2;
00376
00377 if(m_fittingMethod == 2){
00378 status2 = quadfit.QuadFit(px, py, pw, npoints,
00379 &quad_a, &quad_b, &quad_c, &whichhalf, &chisq_quad,
00380 &sigmaquad_a, &sigmaquad_b, &sigmaquad_c);
00381
00382 }
00383
00384
00385 if (slope > 1.0e10 || slope < -1.0e10) {
00386 if (m_Seg > 4) slope *= -1.0;
00387 }
00388
00389 m_SimpleSlope = slope;
00390 m_SimpleSlopeSigma = sigmaSlope;
00391 m_SimpleIntercept = intercept;
00392 m_SimpleInterceptSigma = sigmaIntercept;
00393 m_SimplePosSigma = sigmaPos;
00394 m_SimpleChi2 = chisq;
00395 m_SimpleDOF = ndof;
00396 m_SimpleFitOK = (status == 0) && (chisq < 1000.0);
00397 m_QuadFitOK = (status2 == 1);
00398
00399 m_SimpleQuad_a = quad_a;
00400 m_SimpleQuad_b = quad_b;
00401 m_SimpleQuad_c = quad_c;
00402 m_SimpleQuad_whichhalf = whichhalf;
00403 m_SimpleQuad_aSigma = sigmaquad_a;
00404 m_SimpleQuad_bSigma = sigmaquad_b;
00405 m_SimpleQuad_cSigma = sigmaquad_c;
00406
00407 return status;
00408 }
00409
00410
00412
00413 int
00414 MucRec2DRoad::SimpleFitNoCurrentGap(
00415 int currentgap,
00416 float& slope,
00417 float& intercept,
00418 float& sigmaSlope,
00419 float& sigmaIntercept,
00420 float& chisq,
00421 int& ndof)
00422 {
00423
00424 float px[100];
00425 float py[100];
00426 float pw[100];
00427 int npoints = 0;
00428
00429 int notused = 0;
00430
00431 vector<MucRecHit*>::const_iterator iHit;
00432
00433 float pw2[100];
00434 for(int i = 0; i< 100; i++) {
00435 pw2[i] = 1;
00436
00437 }
00438
00439 for (iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++) {
00440 if (*iHit) {
00441 if( (*iHit)->Gap()==currentgap ) continue;
00442
00443 Hep3Vector center = (*iHit)->GetCenterPos();
00444 Hep3Vector sigma = (*iHit)->GetCenterSigma();
00445
00446
00447 if (m_Part == 1) {
00448 if ( m_Orient == 0) {
00449 px[npoints] = center.z();
00450 py[npoints] = sqrt(center.x()*center.x()
00451 + center.y()*center.y());
00452 if(m_Seg==2) py[npoints] = center.y();
00453
00454 pw[npoints] = 1.0/(sigma.y()*sigma.y())/pw2[npoints]/pw2[npoints];
00455 }
00456 else {
00457 px[npoints] = center.x();
00458 py[npoints] = center.y();
00459 pw[npoints] = 1.0/(sigma.x()*sigma.x())/pw2[npoints]/pw2[npoints];
00460 }
00461 }
00462 else {
00463 if ( m_Orient == 0) {
00464 px[npoints] = center.z();
00465 py[npoints] = center.y();
00466 pw[npoints] = 1.0/(sigma.y()*sigma.y())/pw2[npoints]/pw2[npoints];
00467 }
00468 else {
00469 px[npoints] = center.z();
00470 py[npoints] = center.x();
00471 pw[npoints] = 1.0/(sigma.x()*sigma.x())/pw2[npoints]/pw2[npoints];
00472 }
00473 }
00474
00475
00476
00477 npoints++;
00478
00479 if(npoints > 99)
00480 { cout<<"MucRec2DRoad: to many hits in this track, ignore it!"<<endl;
00481 return -1;
00482 }
00483 }
00484 }
00485
00486 ndof = npoints - 2;
00487 if (ndof < 0) return -1;
00488
00489 if (npoints == 1) {
00490 if (m_Part == 1) {
00491 if ( m_Orient == 0) {
00492 px[npoints] = m_VertexPos.z();
00493 py[npoints] = sqrt(m_VertexPos.x()*m_VertexPos.x()
00494 + m_VertexPos.y()*m_VertexPos.y());
00495 if(m_Seg==2) py[npoints] = m_VertexPos.y();
00496 }
00497 else {
00498 px[npoints] = m_VertexPos.x();
00499 py[npoints] = m_VertexPos.y();
00500 }
00501 }
00502 else {
00503 if ( m_Orient == 0) {
00504 px[npoints] = m_VertexPos.z();
00505 py[npoints] = m_VertexPos.y();
00506 }
00507 else {
00508 px[npoints] = m_VertexPos.z();
00509 py[npoints] = m_VertexPos.x();
00510 }
00511 }
00512 pw[npoints] = 1.0;
00513 npoints++;
00514 }
00515 else {
00516 if (npoints == 0 ) {
00517 return -1;
00518 }
00519
00520 }
00521
00522
00523 MucRecLineFit fit;
00524 int status = fit.LineFit(px, py, pw, npoints,
00525 &slope, &intercept, &chisq,
00526 &sigmaSlope, &sigmaIntercept);
00527
00528 MucRecQuadFit quadfit;
00529 float quad_a, quad_b, quad_c, sigmaquad_a, sigmaquad_b, sigmaquad_c, chisq_quad;
00530 int whichhalf, status2;
00531
00532 if(m_fittingMethod == 2){
00533 status2 = quadfit.QuadFit(px, py, pw, npoints,
00534 &quad_a, &quad_b, &quad_c, &whichhalf, &chisq_quad,
00535 &sigmaquad_a, &sigmaquad_b, &sigmaquad_c);
00536
00537 }
00538
00539
00540 if (slope > 1.0e10 || slope < -1.0e10) {
00541 if (m_Seg > 4) slope *= -1.0;
00542 }
00543
00545 m_SimpleSlope = slope;
00546 m_SimpleSlopeSigma = sigmaSlope;
00547 m_SimpleIntercept = intercept;
00548 m_SimpleInterceptSigma = sigmaIntercept;
00549
00550 m_SimpleChi2 = chisq;
00551 m_SimpleDOF = ndof;
00552 m_SimpleFitOK = (status == 0) && (chisq < 1000.0);
00553 m_QuadFitOK = (status2 == 1);
00554
00555 m_SimpleQuad_a = quad_a;
00556 m_SimpleQuad_b = quad_b;
00557 m_SimpleQuad_c = quad_c;
00558 m_SimpleQuad_whichhalf = whichhalf;
00559 m_SimpleQuad_aSigma = sigmaquad_a;
00560 m_SimpleQuad_bSigma = sigmaquad_b;
00561 m_SimpleQuad_cSigma = sigmaquad_c;
00562
00563
00564 return status;
00565 }
00566
00567
00568
00569
00570
00571
00572
00573
00574 int
00575 MucRec2DRoad::Fit(const float& x,
00576 const float& y,
00577 const float& z,
00578 float& slope, float& intercept,
00579 float& sigmaSlope, float& sigmaIntercept,
00580 float& chisq, int& ndof)
00581 {
00582 int status;
00583
00584
00585 if (!m_SimpleFitOK) {
00586
00587 float a, sa, b, sb, chi;
00588 int n;
00589 status = SimpleFit(a, b, sa, sb, chi, n);
00590 if (status < 0) {
00591
00592
00593 return status;
00594 }
00595 }
00596
00597
00598 float px[100];
00599 float py[100];
00600 float pw[100];
00601 int npoints = 0;
00602 float dx, dy, dr;
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614
00615
00616
00617
00618
00619 vector<MucRecHit*>::const_iterator iHit;
00620
00621
00622
00623
00624
00625
00626
00627 for (iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++) {
00628 if (*iHit) {
00629
00630 Hep3Vector center = (*iHit)->GetCenterPos();
00631 Hep3Vector sigma = (*iHit)->GetCenterSigma();
00632
00633 if (m_Part == 1) {
00634 if ( m_Orient == 0) {
00635 px[npoints] = center.z();
00636 dx = px[npoints] - z;
00637 py[npoints] = sqrt(center.x()*center.x()
00638 + center.y()*center.y());
00639 if(m_Seg==2) py[npoints] = center.y();
00640 dy = py[npoints] - sqrt(x*x + y*y);
00641 }
00642 else {
00643 px[npoints] = center.x();
00644 dx = px[npoints] - x;
00645 py[npoints] = center.y();
00646 dy = py[npoints] - y;
00647 }
00648 }
00649 else {
00650 if ( m_Orient == 0) {
00651 px[npoints] = center.z();
00652 dx = px[npoints] - z;
00653 py[npoints] = center.y();
00654 dy = py[npoints] - y;
00655 }
00656 else {
00657 px[npoints] = center.z();
00658 dx = px[npoints] - z;
00659 py[npoints] = center.x();
00660 dy = py[npoints] - x;
00661 }
00662 }
00663
00664 dr = sqrt(dx*dx + dy*dy);
00665 pw[npoints] = WeightFunc(m_SimpleChi2, dr);
00666
00667
00668
00669
00670
00671
00672
00673 npoints++;
00674
00675 if(npoints > 99)
00676 { cout<<"MucRec2DRoad: to many hits in this track, ignore it!"<<endl;
00677 return -1;
00678 }
00679
00680 }
00681 }
00682
00683
00684 ndof = npoints - 2;
00685 if (npoints == 1) {
00686 if (m_Part == 1) {
00687 if ( m_Orient == 0) {
00688 px[npoints] = m_VertexPos.z();
00689 py[npoints] = sqrt(m_VertexPos.x()*m_VertexPos.x()
00690 + m_VertexPos.y()*m_VertexPos.y());
00691 if(m_Seg==2) py[npoints] = m_VertexPos.y();
00692 }
00693 else {
00694 px[npoints] = m_VertexPos.x();
00695 py[npoints] = m_VertexPos.y();
00696 }
00697 }
00698 else {
00699 if ( m_Orient == 0) {
00700 px[npoints] = m_VertexPos.z();
00701 py[npoints] = m_VertexPos.y();
00702 }
00703 else {
00704 px[npoints] = m_VertexPos.z();
00705 py[npoints] = m_VertexPos.x();
00706 }
00707 }
00708 pw[npoints] = 1.0;
00709 npoints++;
00710 }
00711 else {
00712 if (npoints == 0) {
00713 return -1;
00714 }
00715 }
00716
00717 MucRecLineFit fit;
00718
00719 status = fit.LineFit(px, py, pw, npoints,
00720 &slope, &intercept, &chisq,
00721 &sigmaSlope, &sigmaIntercept);
00722 m_DOF = ndof;
00723 m_Chi2 = chisq;
00724
00725 if (status < 0) {
00726
00727 }
00728
00729 return status;
00730 }
00731
00732
00733 void
00734 MucRec2DRoad::Project(const int& gap,
00735 float& x, float& y, float& z, float& x2, float& y2, float& z2)
00736 {
00737 float sigmaX, sigmaY, sigmaZ;
00738
00739 x = 0.0; sigmaX = 0.0; x2 = 0.0;
00740 y = 0.0; sigmaY = 0.0; y2 = 0.0;
00741 z = 0.0; sigmaZ = 0.0; z2 = 0.0;
00742
00743 if ( (gap < 0) || (gap >= (int)MucID::getGapMax()) ) {
00744 cout << "MucRec2DRoad::Project-E1 invalid gap = " << gap << endl;
00745 return;
00746 }
00747
00748
00749 float x0, y0, z0, sigmaX0, sigmaY0, sigmaZ0;
00750 float phi = m_Seg*0.25*kPi;
00751 MucGeoGeneral *geom = MucGeoGeneral::Instance();
00752
00753 if (m_Part == 1) {
00754 if (m_Orient == 0) {
00755 geom->FindIntersection(m_Part, m_Seg, gap,
00756 m_SimpleSlope*cos(phi), m_SimpleSlope*sin(phi), 1.0,
00757 m_SimpleIntercept*cos(phi),m_SimpleIntercept*sin(phi), 0.0,
00758 m_SimpleSlopeSigma, 0.0, 0.0,
00759 m_SimpleInterceptSigma, 0.0, 0.0,
00760 x0, y0, z0,
00761 sigmaX0, sigmaY0, sigmaZ0);
00762
00763 if(m_SimpleSlope>10000){
00764 geom->FindIntersection(m_Part, m_Seg, gap,
00765 m_SimpleSlope*cos(phi), m_SimpleSlope*sin(phi), 1.0,
00766 0.0, 0.0, m_SimpleIntercept,
00767 m_SimpleSlopeSigma, 0.0, 0.0,
00768 m_SimpleInterceptSigma, 0.0, 0.0,
00769 x0, y0, z0,
00770 sigmaX0, sigmaY0, sigmaZ0);
00771
00772 }
00773
00774 }
00775 else {
00776 geom->FindIntersection(m_Part, m_Seg, gap,
00777 1.0, m_SimpleSlope, 0.0,
00778 0.0, m_SimpleIntercept, 0.0,
00779 0.0, m_SimpleSlopeSigma, 0.0,
00780 0.0, m_SimpleInterceptSigma, 0.0,
00781 x0, y0, z0,
00782 sigmaX0, sigmaY0, sigmaZ0);
00783
00784
00785 }
00786 }
00787 else {
00788 if (m_Orient == 0) {
00789 geom->FindIntersection(m_Part, m_Seg, gap,
00790 m_SimpleSlope, m_SimpleSlope, 1.0,
00791 0.0, m_SimpleIntercept, 0.0,
00792 0.0, m_SimpleSlopeSigma, 0.0,
00793 0.0, m_SimpleInterceptSigma, 0.0,
00794 x0, y0, z0,
00795 sigmaX0, sigmaY0, sigmaZ0);
00796 }
00797 else {
00798 geom->FindIntersection(m_Part, m_Seg, gap,
00799 m_SimpleSlope, m_SimpleSlope, 1.0,
00800 m_SimpleIntercept, 0.0, 0.0,
00801 m_SimpleSlopeSigma, 0.0, 0.0,
00802 m_SimpleInterceptSigma, 0.0, 0.0,
00803 x0, y0, z0,
00804 sigmaX0, sigmaY0, sigmaZ0);
00805 }
00806
00807
00808
00809
00810
00811
00812 }
00813
00814
00815
00816 float a,b,sa,sb,chi2;
00817 int ndof;
00818
00819 int status = Fit(x0, y0, z0, a, b, sa, sb, chi2, ndof);
00820
00821
00822
00823
00824
00825
00826
00827
00828
00829
00830
00831
00832 if (m_Part == 1) {
00833 if (m_Orient == 0) {
00834 geom->FindIntersection(m_Part, m_Seg, gap,
00835 a*cos(phi), a*sin(phi), 1.0,
00836
00837 b*cos(phi), b*sin(phi), 0.0,
00838 sa, 0.0, 0.0,
00839 sb, 0.0, 0.0,
00840 x, y, z,
00841 sigmaX, sigmaY, sigmaZ);
00842
00843 if(fabs(a)>10000){
00844 geom->FindIntersection(m_Part, m_Seg, gap,
00845 a*cos(phi), a*sin(phi), 1.0,
00846 0.0, 0.0, b,
00847
00848
00849 sa, 0.0, 0.0,
00850 sb, 0.0, 0.0,
00851 x, y, z,
00852 sigmaX, sigmaY, sigmaZ);
00853
00854 }
00855 }
00856 else {
00857 geom->FindIntersection(m_Part, m_Seg, gap,
00858 1.0, a, 0.0,
00859 0.0, b, 0.0,
00860 0.0, sa, 0.0,
00861 0.0, sb, 0.0,
00862 x, y, z,
00863 sigmaX, sigmaY, sigmaZ);
00864
00865 if(m_fittingMethod == 2 && m_QuadFitOK ){
00866 float a, b, c;
00867 float sa, sb, sc, chi2x; int ydof; int whichhalf;
00868
00869 GetSimpleFitParams(a, b, c, whichhalf, sa, sb, sc, chi2x, ydof);
00870 geom->FindIntersection(m_Part, m_Seg, gap,
00871 10E30, 0.0, m_SimpleIntercept, 0.0,
00872 a, b, c,
00873 whichhalf,
00874 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
00875 x, y, z, x2, y2, z2,
00876 sigmaX, sigmaY, sigmaZ);
00877
00878
00879 }
00880
00881 }
00882 }
00883 else {
00884 if (m_Orient == 0) {
00885 geom->FindIntersection(m_Part, m_Seg, gap,
00886 a, a, 1.0,
00887 0.0, b, 0.0,
00888 0.0, sa, 0.0,
00889 0.0, sb, 0.0,
00890 x, y, z,
00891 sigmaX, sigmaY, sigmaZ);
00892 }
00893 else {
00894 geom->FindIntersection(m_Part, m_Seg, gap,
00895 a, a, 1.0,
00896 b, 0.0, 0.0,
00897 sa, 0.0, 0.0,
00898 sb, 0.0, 0.0,
00899 x, y, z,
00900 sigmaX, sigmaY, sigmaZ);
00901 }
00902 }
00903
00904 return;
00905 }
00906
00907
00908 int
00909 MucRec2DRoad::GetIndex() const
00910 {
00911 return m_Index;
00912 }
00913
00914
00915 int
00916 MucRec2DRoad::GetPart() const { return m_Part; }
00917
00918
00919 int
00920 MucRec2DRoad::GetSeg() const { return m_Seg; }
00921
00922
00923 int
00924 MucRec2DRoad::GetOrient() const { return m_Orient; }
00925
00926
00927 void
00928 MucRec2DRoad::GetVertexPos(float& x, float& y, float& z) const
00929 {
00930 x = m_VertexPos.x();
00931 y = m_VertexPos.y();
00932 z = m_VertexPos.z();
00933
00934 return;
00935 }
00936
00937
00938 int
00939 MucRec2DRoad::GetLastGap() const { return m_LastGap; }
00940
00941
00942 int
00943 MucRec2DRoad::GetMaxHitsPerGap() const { return m_MaxHitsPerGap; }
00944
00945
00946 int
00947 MucRec2DRoad::GetTotalHits() const { return m_TotalHits; }
00948
00949
00950 int
00951 MucRec2DRoad::GetNGapsWithHits() const
00952 {
00953 const int ngap = (int)MucID::getGapMax();
00954 int gap, count = 0;
00955 vector<int> firedGap;
00956 for ( gap = 0; gap < ngap; gap++) {
00957 firedGap.push_back(0);
00958 }
00959
00960 vector<MucRecHit*>::const_iterator iHit;
00961 for (iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++) {
00962 if (*iHit) {
00963 gap = (*iHit)->Gap();
00964 firedGap[gap] = 1;
00965 }
00966 }
00967
00968 for ( gap = 0; gap < ngap; gap++) {
00969 count += firedGap[gap];
00970 }
00971
00972 return count;
00973 }
00974
00975
00976 int
00977 MucRec2DRoad::GetHitsPerGap(const int& gap) const
00978 {
00979 if ( (gap < 0) || (gap >= (int)MucID::getGapMax()) ) {
00980 cout << "MucRec2DRoad::GetHitsPerGap-E1 invalid gap = " << gap << endl;
00981 return 0;
00982 }
00983
00984 vector<MucRecHit*>::const_iterator iHit;
00985 int hitsInGap = 0;
00986
00987 for (iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++) {
00988
00989 if ( !(*iHit) ) {
00990 cout << "MucRec2DRoad::GetHitsPerGap()-E2 null hit pointer !" << endl;
00991 return 0;
00992 }
00993 else {
00994 if( gap == (*iHit)->Gap() ) hitsInGap += 1;
00995 }
00996 }
00997
00998 return hitsInGap;
00999 }
01000
01001
01002 bool
01003 MucRec2DRoad::HasHitInGap(const int& gap) const
01004 {
01005 if ( (gap < 0) || (gap >= (int)MucID::getGapMax()) ) {
01006 cout << "MucRec2DRoad::HasHitInGap-E2 invalid gap = " << gap << endl;
01007 return false;
01008 }
01009
01010 bool found = false;
01011 vector<MucRecHit*>::const_iterator iHit;
01012
01013 for (iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++) {
01014 if (*iHit) {
01015 if ( (*iHit)->Gap() == gap ) {
01016 found = true;
01017 }
01018 }
01019 }
01020
01021 return found;
01022 }
01023
01024
01025 int
01026 MucRec2DRoad::GetNSharedHits(const MucRec2DRoad* road2) const
01027 {
01028 if (!road2) {
01029 return 0;
01030 }
01031
01032 int count = 0;
01033 vector<MucRecHit*>::const_iterator iHit1;
01034 vector<MucRecHit*>::const_iterator iHit2;
01035 MucRecHit *hit1, *hit2;
01036
01037 for( iHit1 = m_pHits.begin(); iHit1 != m_pHits.end(); iHit1++){
01038 for( iHit2 = road2->m_pHits.begin();
01039 iHit2 != road2->m_pHits.end(); iHit2++){
01040 hit1 = (*iHit1);
01041 hit2 = (*iHit2);
01042
01043 if ( (hit1 != 0) && (hit2 != 0) ) {
01044 if (hit1->GetID() == hit2->GetID()) {
01045 count++;
01046 }
01047 }
01048 }
01049 }
01050
01051 return count;
01052 }
01053
01054
01055 float
01056 MucRec2DRoad::GetSlope() const { return m_SimpleSlope; }
01057
01058
01059 float
01060 MucRec2DRoad::GetIntercept() const { return m_SimpleIntercept; }
01061
01062
01063 int
01064 MucRec2DRoad::GetDegreesOfFreedom() const { return m_SimpleDOF; }
01065
01066
01067 float
01068 MucRec2DRoad::GetReducedChiSquare() const
01069 {
01070 if ( (!m_SimpleFitOK) || (m_SimpleDOF < 0) ) {
01071 return -1.0;
01072 }
01073 else if (m_SimpleDOF == 0) {
01074 return 0.0;
01075 }
01076 else {
01077 return (m_SimpleChi2 / m_SimpleDOF);
01078 }
01079 }
01080
01081
01082 void
01083 MucRec2DRoad::GetSimpleFitParams(float& slope, float& intercept,
01084 float& sigmaSlope, float& sigmaIntercept,
01085 float& chisq, int& ndof) const
01086 {
01087 slope = m_SimpleSlope;
01088 intercept = m_SimpleIntercept;
01089 sigmaSlope = m_SimpleSlopeSigma;
01090 sigmaIntercept = m_SimpleInterceptSigma;
01091 chisq = m_SimpleChi2;
01092 ndof = m_SimpleDOF;
01093
01094 return;
01095 }
01096
01097
01098 void
01099 MucRec2DRoad::GetPosSigma(float &possigma)const
01100 {
01101 possigma = m_SimplePosSigma;
01102
01103 }
01104
01105
01106
01107 void
01108 MucRec2DRoad::GetSimpleFitParams(float& a, float& b, float& c, int& whichhalf,
01109 float& sigmaa, float& sigmab, float& sigmac,
01110 float& chisq, int& ndof) const
01111 {
01112 a = m_SimpleQuad_a;
01113 b = m_SimpleQuad_b;
01114 c = m_SimpleQuad_c;
01115 whichhalf = m_SimpleQuad_whichhalf;
01116
01117 sigmaa = m_SimpleQuad_aSigma;
01118 sigmab = m_SimpleQuad_bSigma;
01119 sigmac = m_SimpleQuad_cSigma;
01120
01121 chisq = m_SimpleChi2;
01122 ndof = m_SimpleDOF;
01123
01124 return;
01125 }
01126
01127
01128
01129
01130 MucRecHit*
01131 MucRec2DRoad::GetHit(const int& gap) const
01132 {
01133 if ( (gap < 0) || (gap >= (int)MucID::getGapMax()) ) {
01134 cout << "MucRec2DRoad::Hit-E1 invalid gap = " << gap << endl;
01135 return 0;
01136 }
01137
01138 vector<MucRecHit*>::const_iterator iHit;
01139
01140 for (iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++) {
01141 if (*iHit) {
01142 if ( (*iHit)->Gap() == gap) {
01143 return (*iHit);
01144 }
01145 }
01146 }
01147
01148 return 0L;
01149 }
01150
01151
01152
01153 float
01154 MucRec2DRoad::GetHitDistance(const MucRecHit* hit)
01155 {
01156
01157 if (!hit) {
01158 cout << "MucRec2DRoad::GetHitDistance-E1 null hit pointer!" << endl;
01159 return muckInfinity;
01160 }
01161
01162 int gap = hit->Gap();
01163 if ( (gap < 0) || (gap >= (int)MucID::getGapMax()) ) {
01164 cout << "MucRec2DRoad::HitDistance-W2 bad gap number = " << gap << endl;
01165 return muckInfinity;
01166 }
01167
01168 if ( hit->GetGap()->Orient() != m_Orient) {
01169
01170 cout << "MucRec2DRoad::GetHitDistance-W3 "
01171 << " hit has wrong orientation = "
01172 << hit->GetGap()->Orient()
01173 << endl;
01174 return muckInfinity;
01175 }
01176
01177 HepPoint3D r = hit->GetCenterPos();
01178 HepPoint3D rl = hit->GetGap()->TransformToGap(r);
01179
01180
01181
01182
01183
01184 float delta, delta1,delta2;
01185 float x0, y0, z0;
01186 float x2, y2, z2;
01187
01188
01189
01190
01191 Project(gap, x0, y0, z0, x2, y2, z2);
01192
01193
01194
01195
01196 if(x0==y0&&x0==z0&&x0==0) {x0+=-9999;y0+=-9999;z0+=-9999;}
01197 if(x2==y2&&x2==z2&&x2==0) {x2+=-9999;y2+=-9999;z2+=-9999;}
01198 HepPoint3D s = HepPoint3D(x0, y0, z0);
01199 HepPoint3D s_2 = HepPoint3D(x2, y2, z2);
01200 HepPoint3D sl = hit->GetGap()->TransformToGap(s);
01201 HepPoint3D s_2l = hit->GetGap()->TransformToGap(s_2);
01202
01203
01204
01205
01206
01207 if ( m_Part == 0 ) {
01208 if ( m_Orient == 0 ) {
01209 delta1 = fabs( sl.y() - rl.y() );
01210 delta2= fabs( s_2l.y() - rl.y() );
01211 }
01212 else {
01213 delta1 = fabs( sl.x() - rl.x() );
01214 delta2= fabs( s_2l.x() - rl.x() );
01215 }
01216 }
01217 else {
01218 if ( m_Orient == 0 ) {
01219 delta1 = fabs( sl.y() - rl.y() );
01220 delta2= fabs( s_2l.y() - rl.y() );
01221 }
01222 else {
01223 delta1 = fabs( sl.x() - rl.x() );
01224 delta2= fabs( s_2l.x() - rl.x() );
01225 }
01226 }
01227
01228
01229
01230
01231
01232
01233 if(delta1 < delta2) return delta1;
01234 else return delta2;
01235 }
01236
01237
01238 float
01239 MucRec2DRoad::GetSearchWindowSize( const int& gap) const
01240 {
01241 if ( (gap < 0) || (gap >= (int)MucID::getGapMax()) ) {
01242 cout << "MucRec2DRoad::GetSearchWindowSize-E1 invalid gap = " << gap << endl;
01243 return 0.0;
01244 }
01245
01246
01247
01248
01249
01250
01251
01252 MucGeoGeneral *geom = MucGeoGeneral::Instance();
01253 float x1, y1, z1, x2, y2, z2, dx, dy, dr;
01254 float sigmaX, sigmaY, sigmaZ;
01255
01256 if (m_Part == 0) {
01257 if (m_Orient == 0) {
01258 geom->FindIntersection(m_Part, 0, m_LastGap,
01259 1.0, m_SimpleSlope, 0.0,
01260 0.0, m_SimpleIntercept, 0.0,
01261 0.0, m_SimpleSlopeSigma, 0.0,
01262 0.0, m_SimpleInterceptSigma, 0.0,
01263 x1, y1, z1,
01264 sigmaX, sigmaY, sigmaZ);
01265 geom->FindIntersection(m_Part, 0, gap,
01266 1.0, m_SimpleSlope, 0.0,
01267 0.0, m_SimpleIntercept, 0.0,
01268 0.0, m_SimpleSlopeSigma, 0.0,
01269 0.0, m_SimpleInterceptSigma, 0.0,
01270 x2, y2, z2,
01271 sigmaX, sigmaY, sigmaZ);
01272 dx = z2 - z1;
01273 dy = sqrt(x2*x2 + y2*y2) - sqrt(x1*x1 + y1*y1);
01274 }
01275 else {
01276 geom->FindIntersection(m_Part, m_Seg, m_LastGap,
01277 m_SimpleSlope, 0.0, 1.0,
01278 m_SimpleIntercept, 0.0, 0.0,
01279 m_SimpleSlopeSigma, 0.0, 0.0,
01280 m_SimpleInterceptSigma, 0.0, 0.0,
01281 x1, y1, z1,
01282 sigmaX, sigmaY, sigmaZ);
01283 geom->FindIntersection(m_Part, m_Seg, gap,
01284 m_SimpleSlope, 0.0, 1.0,
01285 m_SimpleIntercept, 0.0, 0.0,
01286 m_SimpleSlopeSigma, 0.0, 0.0,
01287 m_SimpleInterceptSigma, 0.0, 0.0,
01288 x2, y2, z2,
01289 sigmaX, sigmaY, sigmaZ);
01290 dx = x2 - x1;
01291 dy = y2 - y1;
01292 }
01293 }
01294 else {
01295 if (m_Orient == 0) {
01296 geom->FindIntersection(m_Part, m_Seg, m_LastGap,
01297 0.0, m_SimpleSlope, 1.0,
01298 0.0, m_SimpleIntercept, 0.0,
01299 0.0, m_SimpleSlopeSigma, 0.0,
01300 0.0, m_SimpleInterceptSigma, 0.0,
01301 x1, y1, z1,
01302 sigmaX, sigmaY, sigmaZ);
01303 geom->FindIntersection(m_Part, m_Seg, gap,
01304 0.0, m_SimpleSlope, 1.0,
01305 0.0, m_SimpleIntercept, 0.0,
01306 0.0, m_SimpleSlopeSigma, 0.0,
01307 0.0, m_SimpleInterceptSigma, 0.0,
01308 x2, y2, z2,
01309 sigmaX, sigmaY, sigmaZ);
01310 dx = z2 - z1;
01311 dy = y2 - y1;
01312 }
01313 else {
01314 geom->FindIntersection(m_Part, m_Seg, m_LastGap,
01315 m_SimpleSlope, 0.0, 1.0,
01316 m_SimpleIntercept, 0.0, 0.0,
01317 m_SimpleSlopeSigma, 0.0, 0.0,
01318 m_SimpleInterceptSigma, 0.0, 0.0,
01319 x1, y1, z1,
01320 sigmaX, sigmaY, sigmaZ);
01321 geom->FindIntersection(m_Part, m_Seg, gap,
01322 m_SimpleSlope, 0.0, 1.0,
01323 m_SimpleIntercept, 0.0, 0.0,
01324 m_SimpleSlopeSigma, 0.0, 0.0,
01325 m_SimpleInterceptSigma, 0.0, 0.0,
01326 x2, y2, z2,
01327 sigmaX, sigmaY, sigmaZ);
01328 dx = z2 - z1;
01329 dy = x2 - x1;
01330 }
01331 }
01332
01333 dr = sqrt(dx*dx + dy*dy);
01334
01335 return WindowFunc(m_SimpleChi2,dr);
01336 }
01337
01338
01339 void
01340 MucRec2DRoad::CountHits()
01341 {
01342 m_MaxHitsPerGap = 0;
01343 m_TotalHits = 0;
01344 m_LastGap = 0;
01345
01346 vector<MucRecHit*>::const_iterator iHit;
01347 const int ngap = (int)MucID::getGapNum(m_Part);
01348 vector<int> HitsPerGap;
01349 for (int gap = 0; gap < ngap; gap++) {
01350 HitsPerGap.push_back(0);
01351 }
01352
01353
01354 for (iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++) {
01355 if ( !(*iHit) ) {
01356 cout << "MucRec2DRoad::CountHits()-E1 null hit pointer !" << endl;
01357 return;
01358 }
01359 else {
01360 int gap = (*iHit)->Gap();
01361 HitsPerGap[gap]++;
01362
01363 }
01364 }
01365
01366
01367
01368
01369
01370 int StopGap = ngap;
01371
01372
01373
01374
01375
01376
01377
01378
01379
01380
01381
01382
01383
01384
01385
01386 for (int gap = 0; gap < StopGap; gap++) {
01387 m_TotalHits += HitsPerGap[gap];
01388 if(m_MaxHitsPerGap < HitsPerGap[gap]) m_MaxHitsPerGap = HitsPerGap[gap];
01389 if(HitsPerGap[gap] > 0) m_LastGap = gap;
01390 }
01391 }
01392
01393
01394 bool
01395 MucRec2DRoad::HasHit(MucRecHit* hit) const
01396 {
01397
01398 vector<MucRecHit*>::const_iterator iHit;
01399 bool HitExist = false;
01400
01401
01402
01403
01404 Identifier id = hit->GetID();
01405
01406 for ( iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++ ) {
01407 if ( *iHit ) {
01408 if ( (*iHit)->GetID() == id ) HitExist = true;
01409 }
01410 if (HitExist) break;
01411 }
01412
01413 return HitExist;
01414 }
01415
01416
01417 vector<Identifier>
01418 MucRec2DRoad::GetHitsID() const
01419 {
01420 vector<Identifier> idCon;
01421
01422 vector<MucRecHit*>::const_iterator iHit;
01423 for ( iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++) {
01424 if (*iHit) {
01425 Identifier id = (*iHit)->GetID();
01426 idCon.push_back(id);
01427
01428
01429
01430
01431
01432
01433 }
01434 }
01435 return idCon;
01436 }
01437
01438 vector<MucRecHit*>
01439 MucRec2DRoad::GetHits() const
01440 {
01441 return m_pHits;
01442
01443 }
01444
01445
01446
01447
01448 void
01449 MucRec2DRoad::PrintHitsInfo() const
01450 {
01451 vector<MucRecHit*>::const_iterator iHit;
01452 for ( iHit = m_pHits.begin(); iHit != m_pHits.end(); iHit++) {
01453 if (*iHit) {
01454 float xl, yl, zl;
01455 (*iHit)->GetStrip()->GetCenterPos(xl, yl, zl);
01456 HepPoint3D vl(xl, yl, zl);
01457 HepPoint3D vg = (*iHit)->GetGap()->TransformToGlobal(vl);
01458
01459 cout << " orient " << m_Orient
01460 << " part " << (*iHit)->Part()
01461 << " seg " << (*iHit)->Seg()
01462 << " gap " << (*iHit)->Gap()
01463 << " strip " << (*iHit)->Strip()
01464 << " pos (" << vg.x() << ", " << vg.y() << ", " << vg.z() << ")"
01465 << endl;
01466 }
01467 }
01468
01469 }
01470
01471
01472 float
01473 MucRec2DRoad::WeightFunc(const float& chisq, const float& distance) const
01474 {
01475 return 1.0;
01476 }
01477
01478
01479 float
01480 MucRec2DRoad::WindowFunc(const float& chisq, const float& distance) const
01481 {
01482 return 1.0;
01483 }
01484