00001
00002
00003
00004 #include <iostream>
00005 #include <fstream>
00006 #include <math.h>
00007 #include <stdlib.h>
00008
00009 #include "CLHEP/String/Strings.h"
00010 #include "TrkReco/TCurlFinder.h"
00011 #include "TrkReco/TMDCWire.h"
00012 #include "TrkReco/TMDCWireHit.h"
00013 #include "TrkReco/TMLink.h"
00014 #include "TrkReco/TCircle.h"
00015 #include "TrkReco/TTrack.h"
00016 #include "TrkReco/TSegmentCurl.h"
00017
00018
00019
00020
00021
00022 #include "MdcTables/MdcTables.h"
00023
00024
00025
00026 #define DEBUG_CURL_DUMP 0
00027 #define DEBUG_CURL_SEGMENT 0
00028 #define DEBUG_CURL_GNUPLOT 0
00029 #define DEBUG_CURL_MC 0
00030
00031 #if (DEBUG_CURL_DUMP+DEBUG_CURL_GNUPLOT+DEBUG_CURL_MC)
00032 #include "TrkReco/TMDCWireHitMC.h"
00033 #include "TrkReco/TTrackHEP.h"
00034 #include <set>
00035 #include <algo.h>
00036
00037 static int debugMcFlag = 1;
00038 #endif
00039
00040
00041
00042
00043 TCurlFinder::TCurlFinder(void)
00044 : m_builder("CurlBuilder"),
00045 m_fitter("TCurlFinder Fitter"),
00046 m_debugCdcFrame(false),
00047 m_debugPlotFlag(0),
00048 m_debugFileNumber(0) {
00049
00050
00051
00052 }
00053
00054 TCurlFinder::TCurlFinder(const unsigned min_segment,
00055 const unsigned min_salvage,
00056 const double bad_distance_for_salvage,
00057 const double good_distance_for_salvage,
00058 const unsigned min_sequence,
00059 const unsigned min_fullwire,
00060 const double range_for_axial_search,
00061 const double range_for_stereo_search,
00062 const unsigned superlayer_for_stereo_search,
00063 const double range_for_axial_last2d_search,
00064 const double range_for_stereo_last2d_search,
00065 const double trace2d_distance,
00066 const double trace2d_first_distance,
00067 const double trace3d_distance,
00068 const unsigned determine_one_track,
00069
00070 const double selector_max_impact,
00071 const double selector_max_sigma,
00072 const double selector_strange_pz,
00073 const double selector_replace_dz,
00074
00075 const unsigned stereo_2dfind,
00076 const unsigned merge_exe,
00077 const double merge_ratio,
00078 const double merge_z_diff,
00079 const double mask_distance,
00080 const double ratio_used_wire,
00081 const double range_for_stereo1,
00082 const double range_for_stereo2,
00083 const double range_for_stereo3,
00084 const double range_for_stereo4,
00085 const double range_for_stereo5,
00086 const double range_for_stereo6,
00087
00088 const double z_cut,
00089 const double z_diff_for_last_attend,
00090 const unsigned svd_reconstruction,
00091 const double min_svd_electrons,
00092 const unsigned on_correction,
00093 const unsigned output_2dtracks,
00094 const unsigned curl_version,
00095
00096
00097 const double minimum_seedLength,
00098 const double minimum_2DTrackLength,
00099 const double minimum_3DTrackLength,
00100 const double minimum_closeHitsLength,
00101 const double MIN_RADIUS_OF_STRANGE_TRACK,
00102 const double ULTIMATE_MIN_RADIUS_OF_STRANGE_TRACK)
00103 : m_builder("CurlBuilder"),
00104 m_fitter("TCurlFinder Fitter"),
00105 m_debugCdcFrame(false),
00106 m_debugPlotFlag(0),
00107 m_debugFileNumber(0) {
00108 StatusCode scmgn = Gaudi::svcLocator()->service("MagneticFieldSvc",m_pmgnIMF);
00109 if(scmgn!=StatusCode::SUCCESS) {
00110 std::cout<< __FILE__<<" Unable to open Magnetic field service"<<std::endl;
00111 }
00112
00113 m_param.MIN_SEGMENT = min_segment;
00114 m_param.MIN_SALVAGE = min_salvage;
00115 m_param.BAD_DISTANCE_FOR_SALVAGE = bad_distance_for_salvage;
00116 m_param.GOOD_DISTANCE_FOR_SALVAGE = good_distance_for_salvage;
00117 m_param.MIN_SEQUENCE = min_sequence;
00118 m_param.MAX_FULLWIRE = min_fullwire;
00119 m_param.RANGE_FOR_AXIAL_SEARCH = range_for_axial_search;
00120 m_param.RANGE_FOR_STEREO_SEARCH = range_for_stereo_search;
00121 m_param.SUPERLAYER_FOR_STEREO_SEARCH = superlayer_for_stereo_search;
00122 m_param.RANGE_FOR_AXIAL_LAST2D_SEARCH = range_for_axial_last2d_search;
00123 m_param.RANGE_FOR_STEREO_LAST2D_SEARCH = range_for_stereo_last2d_search;
00124 m_param.TRACE2D_DISTANCE = trace2d_distance;
00125 m_param.TRACE2D_FIRST_SUPERLAYER = trace2d_first_distance;
00126 m_param.TRACE3D_DISTANCE = trace3d_distance;
00127 m_param.DETERMINE_ONE_TRACK = determine_one_track;
00128
00129 m_param.SELECTOR_MAX_IMPACT = selector_max_impact;
00130 m_param.SELECTOR_MAX_SIGMA = selector_max_sigma;
00131 m_param.SELECTOR_STRANGE_PZ = selector_strange_pz;
00132 m_param.SELECTOR_REPLACE_DZ = selector_replace_dz;
00133
00134 m_param.STEREO_2DFIND = stereo_2dfind;
00135 m_param.MERGE_EXE = merge_exe;
00136 m_param.MERGE_RATIO = merge_ratio;
00137 m_param.MERGE_Z_DIFF = merge_z_diff;
00138 m_param.MASK_DISTANCE = mask_distance;
00139 m_param.RATIO_USED_WIRE = ratio_used_wire;
00140 m_param.RANGE_FOR_STEREO_FIRST = range_for_stereo1;
00141 m_param.RANGE_FOR_STEREO_SECOND = range_for_stereo2;
00142 m_param.RANGE_FOR_STEREO_THIRD = range_for_stereo3;
00143 m_param.RANGE_FOR_STEREO_FORTH = range_for_stereo4;
00144 m_param.RANGE_FOR_STEREO_FIFTH = range_for_stereo5;
00145 m_param.RANGE_FOR_STEREO_SIXTH = range_for_stereo6;
00146
00147 m_param.Z_CUT = z_cut;
00148 m_param.Z_DIFF_FOR_LAST_ATTEND = z_diff_for_last_attend;
00149 m_param.SVD_RECONSTRUCTION = svd_reconstruction;
00150 m_param.MIN_SVD_ELECTRONS = min_svd_electrons;
00151 m_param.ON_CORRECTION = on_correction;
00152 m_param.OUTPUT_2DTRACKS = output_2dtracks;
00153 m_param.CURL_VERSION = curl_version;
00154 m_param.now();
00155
00156 m_param.minimum_seedLength = minimum_seedLength;
00157 m_param.minimum_2DTrackLength = minimum_2DTrackLength;
00158 m_param.minimum_3DTrackLength = minimum_3DTrackLength;
00159 m_param.minimum_closeHitsLength = minimum_closeHitsLength;
00160 m_param.MIN_RADIUS_OF_STRANGE_TRACK = MIN_RADIUS_OF_STRANGE_TRACK;
00161 m_param.ULTIMATE_MIN_RADIUS_OF_STRANGE_TRACK = ULTIMATE_MIN_RADIUS_OF_STRANGE_TRACK;
00162
00163
00164
00165 m_builder.setParam(m_param);
00166 }
00167
00168 TCurlFinder::~TCurlFinder(void) {
00169 }
00170
00171
00172
00173
00174 std::string
00175 TCurlFinder::name(void) const {
00176 return std::string("Curling Track Finder");
00177 }
00178
00179 std::string
00180 TCurlFinder::version(void) const {
00181 return std::string("3.00");
00182 }
00183
00184
00185
00186
00187 #if defined(__GNUG__)
00188 int
00189 sortBySequentialLength( const TSegmentCurl **a, const TSegmentCurl **b ) {
00190 if( (*a)->maxSeq() < (*b)->maxSeq() ){
00191 return 1;
00192 }else if( (*a)->maxSeq() == (*b)->maxSeq() ){
00193 return 0;
00194 }else{
00195 return -1;
00196 }
00197 }
00198 #else
00199 extern "C" int
00200 sortBySequentialLength( const void *av, const void *bv ) {
00201 const TSegmentCurl **a((const TSegmentCurl**)av);
00202 const TSegmentCurl **b((const TSegmentCurl**)bv);
00203 if( (*a)->maxSeq() < (*b)->maxSeq() ){
00204 return 1;
00205 }else if( (*a)->maxSeq() == (*b)->maxSeq() ){
00206 return 0;
00207 }else{
00208 return -1;
00209 }
00210 }
00211 #endif
00212
00213 #if defined(__GNUG__)
00214 int
00215 sortByArcLength( const TMLink **a, const TMLink **b ) {
00216 if( (*a)->position().x() > (*b)->position().x() ){
00217 return 1;
00218 }else if( (*a)->position().x() == (*b)->position().x() ){
00219 return 0;
00220 }else{
00221 return -1;
00222 }
00223 }
00224 #else
00225 int
00226 sortByArcLength( const void*av, const void*bv ) {
00227 const TMLink **a((const TMLink**)av);
00228 const TMLink **b((const TMLink**)bv);
00229 if( (*a)->position().x() > (*b)->position().x() ){
00230 return 1;
00231 }else if( (*a)->position().x() == (*b)->position().x() ){
00232 return 0;
00233 }else{
00234 return -1;
00235 }
00236 }
00237 #endif
00238
00239
00240
00241 double
00242 TCurlFinder::distance(const double x, const double y) const {
00243 return sqrt(x*x+y*y);
00244 }
00245
00246 unsigned
00247 TCurlFinder::offset(const unsigned layerId) const {
00248
00249
00250 if( layerId == 0 || layerId == 2 || layerId == 4 ||
00251 layerId == 6 || layerId == 8 || layerId == 9 ||
00252 layerId == 11 || layerId == 13 || layerId == 15 ||
00253 layerId == 17 || layerId == 18 || layerId == 20 ||
00254 layerId == 22 || layerId == 24 || layerId == 26 ||
00255 layerId == 27 || layerId == 29 || layerId == 31 ||
00256 layerId == 33 || layerId == 35 || layerId == 36 ||
00257 layerId == 38 || layerId == 40 || layerId == 42 ||
00258 layerId == 44 || layerId == 45 || layerId == 47 ||
00259 layerId == 49 ) return 1;
00260 return 0;
00261 }
00262
00263 unsigned
00264 TCurlFinder::layerId(const double &R) const {
00265
00266 double r = R*10.;
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287 if(r < 73.0 || r > 772.2)return 43;
00288 if(r <= 85.0)return 0; if(r <= 97.0)return 1; if(r <= 109.0)return 2;
00289 if(r <= 121.0)return 3; if(r <= 133.0)return 4; if(r <= 145.0)return 5;
00290 if(r <= 157.0)return 6; if(r <= 179.0)return 7; if(r <= 205.2)return 8;
00291 if(r <= 221.4)return 9; if(r <= 237.6)return 10; if(r <= 253.8)return 11;
00292 if(r <= 270.0)return 12; if(r <= 286.2)return 13; if(r <= 302.4)return 14;
00293 if(r <= 318.6)return 15; if(r <= 334.8)return 16; if(r <= 351.0)return 17;
00294 if(r <= 367.2)return 18; if(r <= 387.45)return 19; if(r <= 407.7)return 20;
00295 if(r <= 423.9)return 21; if(r <= 440.1)return 22; if(r <= 456.3)return 23;
00296 if(r <= 472.5)return 24; if(r <= 488.7)return 25; if(r <= 504.9)return 26;
00297 if(r <= 521.1)return 27; if(r <= 537.3)return 28; if(r <= 554.5)return 29;
00298 if(r <= 569.7)return 30; if(r <= 585.9)return 31; if(r <= 602.1)return 32;
00299 if(r <= 618.3)return 33; if(r <= 634.5)return 34; if(r <= 654.75)return 35;
00300 if(r <= 675.0)return 36; if(r <= 691.2)return 37; if(r <= 707.4)return 38;
00301 if(r <= 723.6)return 39; if(r <= 739.8)return 40; if(r <= 756.0)return 41;
00302 if(r <= 772.2)return 42;
00303 }
00304
00305 unsigned
00306 TCurlFinder::maxLocalLayerId(const unsigned superLayerId) const {
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316 if(superLayerId == 10)return 2;
00317 if(superLayerId >= 0 && superLayerId < 10)return 3;
00318
00319 std::cerr << "Error in the CurlFinder(maxLocalLayerId). superLayerId = "
00320 << superLayerId << std::endl;
00321 return 0;
00322 }
00323
00324 int
00325 TCurlFinder::nextSuperAxialLayerId(const unsigned superLayerId, const int in) const {
00326
00327
00328
00329
00330
00331
00332
00333
00334 if(superLayerId > 10 || superLayerId < 0){
00335 std::cerr << "Error in the CurlFinder(nextSuperAxialLayerId)." << std::endl;
00336 return -1;
00337 }
00338
00339 if(in == 0){
00340 if(superLayerId == 2 || superLayerId == 3 ||
00341 superLayerId == 4 || superLayerId == 9 ||
00342 superLayerId == 10){
00343 return superLayerId;
00344 }else{
00345
00346 return -1;
00347 }
00348 }
00349
00350 if((superLayerId == 3 ) && in == 1)return 2;
00351 if(superLayerId == 4){
00352 if(in == 1)return 3; if(in == 2)return 2;
00353 }
00354 if(superLayerId == 5 || superLayerId == 6 ||
00355 superLayerId == 7 || superLayerId == 8 ||
00356 superLayerId == 9){
00357 if(in == 1)return 4; if(in == 2)return 3;
00358 if(in == 3)return 2;
00359 }
00360 if(superLayerId == 10){
00361 if(in == 1)return 9; if(in == 2)return 4;
00362 if(in == 3)return 3; if(in == 4)return 2;
00363 }
00364
00365 if(superLayerId == 0 || superLayerId == 1){
00366 if(in == -1)return 2; if(in == -2)return 3;
00367 if(in == -3)return 4; if(in == -4)return 9;
00368 if(in == -5)return 10;
00369 }
00370 if( superLayerId == 2){
00371 if(in == -1)return 3; if(in == -2)return 4;
00372 if(in == -3)return 9; if(in == -4)return 10;
00373 }
00374 if(superLayerId == 3){
00375 if(in == -1)return 4; if(in == -2)return 9;
00376 if(in == -3)return 10;
00377 }
00378 if(superLayerId == 4){
00379 if(in == -1)return 9; if(in == -2)return 10;
00380 }
00381 if(superLayerId == 5 || superLayerId == 6 ||
00382 superLayerId == 7 || superLayerId == 8 ||
00383 superLayerId == 9){
00384 if(in == -1)return 10;
00385 }
00386
00387 return -1;
00388 }
00389
00390 int
00391 TCurlFinder::nextSuperStereoLayerId(const unsigned superLayerId, const int in) const {
00392
00393
00394
00395
00396
00397
00398
00399
00400 if(superLayerId > 10 || superLayerId < 0){
00401 std::cerr << "Error in the CurlFinder(nextSuperStereoLayerId)." << std::endl;
00402 return -1;
00403 }
00404
00405 if(in == 0){
00406 if(superLayerId == 0 || superLayerId == 1 ||
00407 superLayerId == 5 || superLayerId == 6 ||
00408 superLayerId == 7 || superLayerId == 8){
00409 return superLayerId;
00410 }else{
00411 return -1;
00412 }
00413 }
00414
00415 if((superLayerId == 1 || superLayerId == 2 ||
00416 superLayerId == 3 || superLayerId == 4) && in == 1)return 0;
00417 if(superLayerId == 5){
00418 if(in == 1)return 1; if(in == 2)return 0;
00419 }
00420 if(superLayerId == 6){
00421 if(in == 1)return 5; if(in == 2)return 1;
00422 if(in == 3)return 0;
00423 }
00424 if(superLayerId == 7){
00425 if(in == 1)return 6; if(in == 2)return 5;
00426 if(in == 3)return 1; if(in == 4)return 0;
00427 }
00428 if(superLayerId == 8){
00429 if(in == 1)return 7; if(in == 2)return 6;
00430 if(in == 3)return 5; if(in == 4)return 1;
00431 if(in == 5)return 0;
00432 }
00433 if(superLayerId == 9 || superLayerId == 10){
00434 if(in == 1)return 8; if(in == 2)return 7;
00435 if(in == 3)return 6; if(in == 4)return 5;
00436 if(in == 5)return 1; if(in == 6)return 0;
00437 }
00438
00439 if(superLayerId == 0){
00440 if(in == -1)return 1; if(in == -2)return 5;
00441 if(in == -3)return 6; if(in == -4)return 7;
00442 if(in == -5)return 8;
00443 }
00444 if(superLayerId == 1 || superLayerId == 2 ||
00445 superLayerId == 3 || superLayerId == 4){
00446 if(in == -1)return 5; if(in == -2)return 6;
00447 if(in == -3)return 7; if(in == -4)return 8;
00448 }
00449 if(superLayerId == 5){
00450 if(in == -1)return 6; if(in == -2)return 7;
00451 if(in == -3)return 8;
00452 }
00453 if(superLayerId == 6){
00454 if(in == -1)return 7; if(in == -2)return 8;
00455 }
00456 if(superLayerId == 7){
00457 if(in == -1)return 8;
00458 }
00459
00460 return -1;
00461 }
00462
00463 unsigned
00464 TCurlFinder::nAxialHits(const double &r) const {
00465
00466 const double eps = 0.2;
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502 if(r < 20.52-eps) return 0;
00503 else if(r < 22.14-eps) return 1;
00504 else if(r < 23.76-eps) return 2;
00505 else if(r < 25.38-eps) return 3;
00506 else if(r < 27.0-eps) return 4;
00507 else if(r < 28.62-eps) return 5;
00508 else if(r < 30.24-eps) return 6;
00509 else if(r < 31.86-eps) return 7;
00510 else if(r < 33.48-eps) return 8;
00511 else if(r < 35.1-eps) return 9;
00512 else if(r < 36.72-eps) return 10;
00513 else if(r < 38.34-eps) return 11;
00514 else if(r < 67.5-eps) return 12;
00515 else if(r < 69.12-eps) return 13;
00516 else if(r < 70.74-eps) return 14;
00517 else if(r < 72.36-eps) return 15;
00518 else if(r < 73.98-eps) return 16;
00519 else if(r < 75.6-eps) return 17;
00520 else if(r < 77.22-eps) return 18;
00521 }
00522
00523
00524
00525
00526 void
00527 TCurlFinder::makeList(AList<TMLink> &madeList,
00528 const AList<TSegmentCurl> &originalList,
00529 const AList<TMLink> &removeList) {
00530
00531
00532
00533 madeList.removeAll();
00534 for(unsigned i = 0, size = originalList.length(); i < size; ++i)
00535 madeList.append(originalList[i]->list());
00536 madeList.remove(removeList);
00537 }
00538
00539 void
00540 TCurlFinder::makeList(AList<TMLink> &madeList,
00541 const AList<TMLink> &originalList,
00542 const AList<TMLink> &removeList) {
00543
00544
00545
00546 madeList.removeAll();
00547 madeList.append(originalList);
00548 madeList.remove(removeList);
00549 }
00550
00551 void
00552 TCurlFinder::clear(void) {
00553
00554
00555 HepAListDeleteAll(m_allAxialHitsOriginal);
00556 HepAListDeleteAll(m_allStereoHitsOriginal);
00557 HepAListDeleteAll(m_segmentList);
00558 HepAListDeleteAll(m_allCircles);
00559 HepAListDeleteAll(m_allTracks);
00560
00561 m_unusedAxialHitsOriginal.removeAll();
00562 m_unusedStereoHitsOriginal.removeAll();
00563 m_unusedAxialHits.removeAll();
00564 m_unusedStereoHits.removeAll();
00565 m_removedHits.removeAll();
00566 m_circles.removeAll();
00567 m_tracks.removeAll();
00568 m_2dTracks.removeAll();
00569
00570 for(int i=0;i<19;++i)
00571 m_unusedAxialHitsOnEachLayer[i].removeAll();
00572 for(int i=0;i<24;++i)
00573 m_unusedStereoHitsOnEachLayer[i].removeAll();
00574 for(int i=0;i<5;++i)
00575 m_unusedAxialHitsOnEachSuperLayer[i].removeAll();
00576 for(int i=0;i<6;++i)
00577 m_unusedStereoHitsOnEachSuperLayer[i].removeAll();
00578 m_hitsOnInnerSuperLayer.removeAll();
00579 }
00580
00581
00582
00583
00584 int
00585 TCurlFinder::doit(const AList<TMDCWireHit> & axialHits,
00586 const AList<TMDCWireHit> & stereoHits,
00587 AList<TTrack> & tracks,
00588 AList<TTrack> & tracks2D) {
00589 #if (DEBUG_CURL_DUMP+DEBUG_CURL_GNUPLOT+DEBUG_CURL_MC)
00590 Belle_event_Manager &evtMgr = Belle_event_Manager::get_manager();
00591 debugMcFlag = 1;
00592 if(evtMgr.count() != 0 &&
00593 evtMgr[0].ExpMC() != 2)debugMcFlag = 0;
00594 m_debugCdcFrame = false;
00595 #endif
00596 #if !(DEBUG_CURL_MC)
00597 #if DEBUG_CURL_DUMP
00598 std::cout << "(TCurlFinder)Plot Menu : All Off = 0, Interactive = 1, All On = 2" << std::endl;
00599 cin >> m_debugPlotFlag;
00600 #endif
00601
00602
00603 makeWireHitsListsSegments(axialHits, stereoHits);
00604 #if DEBUG_CURL_SEGMENT
00605 std::cout << "(TCurlFinder)# of segment = " << m_segmentList.length() << std::endl;
00606 debugCheckSegments0();
00607 debugCheckSegments1();
00608 debugCheckSegments2();
00609 #endif
00610
00611 if(checkSortSegments() == 0)return 0;
00612 #if DEBUG_CURL_DUMP
00613 if(m_debugPlotFlag){
00614 int noPlot = 1;
00615 if(m_debugPlotFlag == 1){
00616 std::cout << "(TCurlFinder) Do you want to see Segment Plot? : yes = 1, no = other #" << std::endl;
00617 cin >> noPlot;
00618 }
00619 if(noPlot == 1){
00620 for(int i=0;i<m_segmentList.length();++i)
00621 plotSegment(m_segmentList[i]->list(),0);
00622 }
00623 }
00624 #endif
00625
00626 makeCurlTracks(tracks,tracks2D);
00627 #else
00628 makeWithMC(axialHits, stereoHits, tracks);
00629 #endif
00630
00631
00632 unsigned n = tracks2D.length();
00633 for (unsigned i = 0; i < n; i++)
00634 tracks2D[i]->quality(TrackQuality2D);
00635
00636 return 0;
00637 }
00638
00639
00640
00641
00642
00643 void
00644 TCurlFinder::makeWireHitsListsSegments(const AList<TMDCWireHit> &axialList,
00645 const AList<TMDCWireHit> &stereoList) {
00646
00647
00648
00649
00650
00651
00652
00653 unsigned size = axialList.length();
00654 for(unsigned i=0;i<size;++i){
00655 if(axialList[i]->reccdc()->tdc > 500) continue;
00656 if(axialList[i]->state() & WireHitFindingValid){
00657 m_allAxialHitsOriginal.append(new TMLink(0,axialList[i]));
00658 if(axialList[i]->wire()->superLayerId() <= 3)
00659 m_hitsOnInnerSuperLayer.append(axialList[i]);
00660
00661 }
00662 }
00663 size = m_allAxialHitsOriginal.length();
00664 for(unsigned i=0;i<size;++i){
00665 if(!(m_allAxialHitsOriginal[i]->hit()->state() & WireHitUsed)){
00666 if(m_allAxialHitsOriginal[i]->hit()->state() & WireHitInvalidForFit){
00667 unsigned newState = m_allAxialHitsOriginal[i]->hit()->state()&(~WireHitInvalidForFit);
00668 m_allAxialHitsOriginal[i]->hit()->state(newState);
00669 }
00670 m_unusedAxialHitsOriginal.append(m_allAxialHitsOriginal[i]);
00671 }
00672 }
00673 m_unusedAxialHits = m_unusedAxialHitsOriginal;
00674
00675 size = stereoList.length();
00676 for(unsigned i=0;i<size;++i){
00677 if (stereoList[i]->reccdc()->tdc > 500) continue;
00678 if(stereoList[i]->state() & WireHitFindingValid){
00679 m_allStereoHitsOriginal.append(new TMLink(0,stereoList[i]));
00680 if(stereoList[i]->wire()->superLayerId() <= 3)
00681 m_hitsOnInnerSuperLayer.append(stereoList[i]);
00682 }
00683 }
00684 size = m_allStereoHitsOriginal.length();
00685 for(unsigned i=0;i<size;++i){
00686 if(!(m_allStereoHitsOriginal[i]->hit()->state() & WireHitUsed)){
00687 if(m_allStereoHitsOriginal[i]->hit()->state() & WireHitInvalidForFit){
00688 unsigned newState = m_allStereoHitsOriginal[i]->hit()->state()&(~WireHitInvalidForFit);
00689 m_allStereoHitsOriginal[i]->hit()->state(newState);
00690 }
00691 m_unusedStereoHitsOriginal.append(m_allStereoHitsOriginal[i]);
00692 }
00693 }
00694 m_unusedStereoHits = m_unusedStereoHitsOriginal;
00695
00696
00697 size = m_unusedAxialHitsOriginal.length();
00698 for(unsigned i=0;i<size;++i){
00699 m_unusedAxialHitsOnEachLayer[m_unusedAxialHitsOriginal[i]->hit()->wire()->
00700 axialStereoLayerId()].append(m_unusedAxialHitsOriginal[i]);
00701 }
00702 size = m_unusedStereoHitsOriginal.length();
00703 for(unsigned i=0;i<size;++i){
00704 m_unusedStereoHitsOnEachLayer[m_unusedStereoHitsOriginal[i]->hit()->wire()->
00705 axialStereoLayerId()].append(m_unusedStereoHitsOriginal[i]);
00706 }
00707
00708
00709 linkNeighboringWires(m_unusedAxialHitsOnEachLayer,19);
00710 linkNeighboringWires(m_unusedStereoHitsOnEachLayer,24);
00711
00712
00713 createSuperLayer();
00714
00715 m_segmentList.removeAll();
00716 for(unsigned i=0;i<5;++i)
00717 if(m_unusedAxialHitsOnEachSuperLayer[i].length() > 0)
00718 createSegments(m_unusedAxialHitsOnEachSuperLayer[i]);
00719 for(unsigned i=0;i<6;++i)
00720 if(m_unusedStereoHitsOnEachSuperLayer[i].length() > 0)
00721 createSegments(m_unusedStereoHitsOnEachSuperLayer[i]);
00722
00723 }
00724
00725 void
00726 TCurlFinder::linkNeighboringWires(AList<TMLink> *list, const unsigned num) {
00727
00728
00729
00730
00731 for(int i=0;i<num;++i){
00732 if(list[i].length() == 0)continue;
00733 for(int j=0;j<list[i].length();++j){
00734
00735
00736 if(num == 19){
00737 if( i == 0 || i == 4 || i == 8 ||
00738 i == 12 || i == 16)goto outer;
00739 }else if(num == 24){
00740 if( i == 0 || i == 4 ||
00741 i == 8 || i == 12 || i == 16 || i == 20 )goto outer;
00742 }
00743
00744 for(int k=0;k<list[i-1].length();++k){
00745 if(list[i-1][k]->hit()->wire()->localId() ==
00746 list[i][j]->wire()->neighbor(1)->localId())
00747 setNeighboringWires(list[i][j], list[i-1][k]);
00748 else if(list[i-1][k]->hit()->wire()->localId() ==
00749 list[i][j]->wire()->neighbor(0)->localId())
00750 setNeighboringWires(list[i][j], list[i-1][k]);
00751 }
00752 outer:
00753
00754 if(num == 19){
00755 if( i == 3 || i == 7 || i == 11 ||
00756 i == 15 || i == 18)goto same;
00757 }else if(num == 24){
00758 if( i == 3 || i == 7 || i == 11 ||
00759 i == 15 || i == 19 || i == 23 )goto same;
00760 }
00761 for(int k=0;k<list[i+1].length();++k){
00762 if(list[i+1][k]->hit()->wire()->localId() ==
00763 list[i][j]->hit()->wire()->neighbor(4)->localId())
00764 setNeighboringWires(list[i][j], list[i+1][k]);
00765 else if(list[i+1][k]->wire()->localId() ==
00766 list[i][j]->wire()->neighbor(5)->localId())
00767 setNeighboringWires(list[i][j], list[i+1][k]);
00768 }
00769 same:
00770
00771 for(int k=0;k<list[i].length();++k){
00772 if(list[i][k]->hit()->wire()->localId() ==
00773 list[i][j]->hit()->wire()->localIdForPlus()+1){
00774 setNeighboringWires(list[i][j], list[i][k]);
00775 }else if(list[i][k]->hit()->wire()->localId() ==
00776 list[i][j]->hit()->wire()->localIdForMinus()-1){
00777 setNeighboringWires(list[i][j], list[i][k]);
00778 }
00779 }
00780 }
00781 }
00782 }
00783
00784 void
00785 TCurlFinder::setNeighboringWires(TMLink *list, const TMLink *next) {
00786
00787
00788 for(int i=0;i<6;++i){
00789 if(!list->neighbor(i)){
00790 list->neighbor(i,const_cast<TMLink*>(next));
00791 break;
00792 }
00793 }
00794 }
00795
00796 void
00797 TCurlFinder::createSuperLayer(void) {
00798
00799 for(int i=0;i<4;++i){
00800 if(m_unusedAxialHitsOnEachLayer[i].length() > 0)
00801 m_unusedAxialHitsOnEachSuperLayer[0].append(m_unusedAxialHitsOnEachLayer[i]);
00802 if(m_unusedAxialHitsOnEachLayer[i+4].length() > 0)
00803 m_unusedAxialHitsOnEachSuperLayer[1].append(m_unusedAxialHitsOnEachLayer[i+4]);
00804 if(m_unusedAxialHitsOnEachLayer[i+8].length() > 0)
00805 m_unusedAxialHitsOnEachSuperLayer[2].append(m_unusedAxialHitsOnEachLayer[i+8]);
00806 if(m_unusedAxialHitsOnEachLayer[i+12].length() > 0)
00807 m_unusedAxialHitsOnEachSuperLayer[3].append(m_unusedAxialHitsOnEachLayer[i+12]);
00808 if(m_unusedAxialHitsOnEachLayer[i+16].length() > 0 && i+16 < 19)
00809 m_unusedAxialHitsOnEachSuperLayer[4].append(m_unusedAxialHitsOnEachLayer[i+16]);
00810 }
00811 for(int i=0;i<4;++i){
00812 if(m_unusedStereoHitsOnEachLayer[i].length() > 0)
00813 m_unusedStereoHitsOnEachSuperLayer[0].append(m_unusedStereoHitsOnEachLayer[i]);
00814 if(m_unusedStereoHitsOnEachLayer[i+4].length() > 0)
00815 m_unusedStereoHitsOnEachSuperLayer[1].append(m_unusedStereoHitsOnEachLayer[i+4]);
00816 if(m_unusedStereoHitsOnEachLayer[i+8].length() > 0)
00817 m_unusedStereoHitsOnEachSuperLayer[2].append(m_unusedStereoHitsOnEachLayer[i+8]);
00818 if(m_unusedStereoHitsOnEachLayer[i+12].length() > 0)
00819 m_unusedStereoHitsOnEachSuperLayer[3].append(m_unusedStereoHitsOnEachLayer[i+12]);
00820 if(m_unusedStereoHitsOnEachLayer[i+16].length() > 0)
00821 m_unusedStereoHitsOnEachSuperLayer[4].append(m_unusedStereoHitsOnEachLayer[i+16]);
00822 if(m_unusedStereoHitsOnEachLayer[i+20].length() > 0)
00823 m_unusedStereoHitsOnEachSuperLayer[5].append(m_unusedStereoHitsOnEachLayer[i+20]);
00824 }
00825 }
00826
00827 void
00828 TCurlFinder::createSegments(AList<TMLink> &list) {
00829
00830
00831 AList<TMLink> seedStock;
00832 do{
00833 TSegmentCurl *segment = new TSegmentCurl(list[0]->hit()->wire()->superLayerId(),
00834 maxLocalLayerId(list[0]->hit()->wire()
00835 ->superLayerId()));
00836
00837 segment->append(list[0]);
00838 TMLink *seed = list[0];
00839 list.remove(seed);
00840 next:
00841 searchSegment(seed, list, seedStock, segment);
00842 if(seedStock.length() > 0){
00843 seed = seedStock[0];
00844 seedStock.remove(seed);
00845 goto next;
00846 }else if(segment->size() >= m_param.MIN_SEGMENT){
00847 segment->update();
00848 m_segmentList.append(segment);
00849 #if DEBUG_CURL_DUMP
00850 std::cout << "Segment # = " << m_segmentList.length() << std::endl;
00851 segment->dump();
00852 #endif
00853 }else{
00854 delete segment;
00855 }
00856 }while(list.length() > 0);
00857 }
00858
00859 void
00860 TCurlFinder::searchSegment(TMLink *seed, AList<TMLink> &list,
00861 AList<TMLink> &seedStock, TSegmentCurl *segment) {
00862 for(int i=0;i<6;++i){
00863 if(seed->neighbor(i)){
00864 if(!findLink(seed->neighbor(i),list))continue;
00865 segment->append(seed->neighbor(i));
00866 seedStock.append(seed->neighbor(i));
00867 list.remove(seed->neighbor(i));
00868 }else{
00869 break;
00870 }
00871 }
00872 }
00873
00874 TMLink *
00875 TCurlFinder::findLink(const TMLink *seed, const AList<TMLink> &list) {
00876
00877
00878
00879 unsigned size = list.length();
00880 if(size == 0)return NULL;
00881 for(unsigned i=0;i<size;++i){
00882 if(seed == list[i])return list[i];
00883 }
00884 return NULL;
00885 }
00886
00887
00888
00889
00890
00891 int
00892 TCurlFinder::checkSortSegments(void) {
00893
00894 #if DEBUG_CURL_DUMP
00895 std::cout << "(TCurlFinder)checking and sorting segments..." << std::endl;
00896 #endif
00897 unsigned length = m_segmentList.length();
00898 if(length == 0)return 0;
00899 checkExceptionalSegmentsType03();
00900
00901 checkExceptionalSegmentsType01();
00902 m_segmentList.sort(sortBySequentialLength);
00903 #if DEBUG_CURL_DUMP
00904 std::cout << "(TCurlFinder)...done check and sort of segments." << std::endl;
00905 #endif
00906 return 1;
00907 }
00908
00909 void
00910 TCurlFinder::checkExceptionalSegmentsType03(void) {
00911 int max = m_param.MAX_FULLWIRE;
00912 int nMinWires;
00913 if(max == 7)nMinWires = 21;
00914 else if(max == 6)nMinWires = 19;
00915 else if(max == 5)nMinWires = 18;
00916 else if(max == 4)nMinWires = 16;
00917 else if(max == 3)nMinWires = 14;
00918 else if(max == 2)nMinWires = 12;
00919 else if(max == 1)nMinWires = 10;
00920 else if(max == 0)nMinWires = 7;
00921
00922 AList<TSegmentCurl> removeList;
00923 for(unsigned i = 0, length = m_segmentList.length(); i < length; ++i){
00924 if(m_segmentList[i]->size() >= nMinWires){
00925 unsigned nWires = m_segmentList[i]->size();
00926 unsigned n6Wires = 0;
00927 for(unsigned j=0;j<nWires;++j){
00928 if(((m_segmentList[i]->list())[j])->neighbor(5))++n6Wires;
00929 if(n6Wires > max)break;
00930 }
00931 if(n6Wires <= max)continue;
00932 removeList.append(m_segmentList[i]);
00933 #if DEBUG_CURL_SEGMENT
00934 writeSegment(m_segmentList[i]->list(),3);
00935 #endif
00936 }
00937 }
00938 if(removeList.length() >= 1){
00939 #if DEBUG_CURL_DUMP
00940 std::cout << "(TCurlFinder)removing large segments: # = " << removeList.length() << std::endl;
00941 #endif
00942 m_segmentList.remove(removeList);
00943 HepAListDeleteAll(removeList);
00944 }
00945 }
00946
00947 void
00948 TCurlFinder::checkExceptionalSegmentsType02(void) {
00949 int max = 10;
00950 int hmax = 5;
00951 AList<TSegmentCurl> removeList;
00952 for(unsigned i = 0, length = m_segmentList.length(); i < length; ++i){
00953 int lSize = max*3+hmax*3;
00954 int lNum = 3;
00955 if(m_segmentList[i]->superLayerId() == 1 ||
00956 m_segmentList[i]->superLayerId() == 3){
00957 lSize = max*2+hmax;
00958 lNum = 2;
00959 }
00960 if(m_segmentList[i]->superLayerId() == 5 ||
00961 m_segmentList[i]->superLayerId() == 7 ||
00962 m_segmentList[i]->superLayerId() == 9){
00963 lSize = max*2+hmax*2;
00964 lNum = 2;
00965 }
00966 if(m_segmentList[i]->superLayerId() == 4 ||
00967 m_segmentList[i]->superLayerId() == 6 ||
00968 m_segmentList[i]->superLayerId() == 8 ||
00969 m_segmentList[i]->superLayerId() == 10)lSize = max*3+hmax*2;
00970 if(m_segmentList[i]->size() < lSize)continue;
00971 int nL = 0;
00972 for(unsigned j=0,size=m_segmentList[i]->maxLocalLayerId();j<size;++j){
00973 if(m_segmentList[i]->sizeOfLayer(j) >= max)++nL;
00974 }
00975 if(nL < lNum)continue;
00976 removeList.append(m_segmentList[i]);
00977
00978 #if DEBUG_CURL_SEGMENT
00979
00980 #endif
00981 }
00982 if(removeList.length() >= 1){
00983 #if DEBUG_CURL_DUMP
00984 std::cout << "(TCurlFinder)removing large segments: # = " << removeList.length() << std::endl;
00985 #endif
00986 m_segmentList.remove(removeList);
00987 HepAListDeleteAll(removeList);
00988 }
00989 }
00990
00991 void
00992 TCurlFinder::checkExceptionalSegmentsType01(void) {
00993 for(unsigned i = 0, length = m_segmentList.length(); i < length; ++i){
00994 if(m_segmentList[i]->maxLocalLayerId() != m_segmentList[i]->layerIdOfMaxSeq() &&
00995 m_segmentList[i]->maxSeq() >= m_param.MIN_SEQUENCE){
00996 unsigned innerHits = 0;
00997 if(m_segmentList[i]->layerIdOfMaxSeq() == 0)continue;
00998 TSegmentCurl *outer = new TSegmentCurl(m_segmentList[i]->superLayerId(),
00999 m_segmentList[i]->maxLocalLayerId());
01000 for(unsigned j = 0, size = m_segmentList[i]->size(); j < size; ++j){
01001 if(m_segmentList[i]->layerIdOfMaxSeq()+1 <=
01002 (m_segmentList[i]->list())[j]->hit()->wire()->localLayerId() &&
01003 (m_segmentList[i]->list())[j]->hit()->wire()->localLayerId() <=
01004 m_segmentList[i]->maxLocalLayerId()){
01005 outer->append((m_segmentList[i]->list())[j]);
01006 }else if(m_segmentList[i]->layerIdOfMaxSeq()-1 >=
01007 (m_segmentList[i]->list())[j]->hit()->wire()->localLayerId()){
01008 ++innerHits;
01009 }
01010 }
01011 if(innerHits != 0 && outer->size() != 0){
01012 #if DEBUG_CURL_DUMP
01013 std::cout << "(TCurlFinder)removing some wires in the segment." << std::endl;
01014 #endif
01015 #if DEBUG_CURL_SEGMENT
01016
01017 #endif
01018 m_segmentList[i]->remove(const_cast< AList<TMLink>& >(outer->list()));
01019 outer->removeAll();
01020 delete outer;
01021 }else{
01022 outer->removeAll();
01023 delete outer;
01024 }
01025 }
01026 }
01027 }
01028
01029
01030
01031
01032
01033
01034
01035
01036
01037
01038 void
01039 TCurlFinder::makeCurlTracks(AList<TTrack> &tracks,
01040 AList<TTrack> &tracks2D) {
01041 AList<TSegmentCurl> segmentList = m_segmentList;
01042
01043
01044 for(unsigned i = 0, size = m_segmentList.length(); i < size; ++i){
01045 TCircle *circle = make2DTrack(segmentList[i]->list(), segmentList, 1);
01046 if(circle){
01047 AList<TMLink> tmp(circle->links());
01048 #if DEBUG_CURL_DUMP
01049 std::cout << "(TCurlFinder) 2D:Created Circle!!!" << std::endl;
01050 if(m_debugPlotFlag){
01051 int noPlot = 1;
01052 if(m_debugPlotFlag == 1){
01053 std::cout << "(TCurlFinder) Do you want to see Circle Plot(2D)? : yes = 1, no = other #" << std::endl;
01054 cin >> noPlot;
01055 }
01056 if(noPlot == 1)plotCircle(*circle,0);
01057 }
01058 #endif
01059
01060
01061
01062
01063
01064
01065
01066
01067
01068
01069 if(TCircle *dividedCircle = dividing2DTrack(circle)){
01070 #if DEBUG_CURL_DUMP
01071 std::cout << "(TCurlFinder) 2D:dividing...good...2 Circles!!!" << std::endl;
01072 #endif
01073 TTrack *track1(NULL), *track2(NULL);
01074 int ok2d[2] = { 0, 0 };
01075 int ok3d[2] = { 0, 0 };
01076
01077
01078 if(trace2DTrack(circle) &&
01079 check2DCircle(circle) &&
01080 circle->fitForCurl(1) != -1){
01081 ok2d[0] = 1;
01082 #if DEBUG_CURL_DUMP
01083 std::cout << "(TCurlFinder) 2D:Success Circle Fit!!!" << std::endl;
01084 #endif
01085 track1 = make3DTrack(circle, segmentList);
01086 }
01087
01088
01089 if(trace2DTrack(dividedCircle) &&
01090 check2DCircle(dividedCircle) &&
01091 dividedCircle->fitForCurl(1) != -1){
01092 ok2d[1] = 1;
01093 #if DEBUG_CURL_DUMP
01094 std::cout << "(TCurlFinder) 2D:Success Circle Fit!!!" << std::endl;
01095 #endif
01096 track2 = make3DTrack(dividedCircle, segmentList);
01097 }
01098 if(track1 && track2){
01099 #if DEBUG_CURL_DUMP
01100 std::cout << "(TCurlFinder) 3D:Create Track!!! in track1 && track2" << std::endl;
01101 #endif
01102 salvage3DTrack(track1,true);
01103 salvage3DTrack(track2,true);
01104 #if DEBUG_CURL_DUMP
01105 std::cout << "(TCurlFinder) 1:dz = " << track1->helix().dz()
01106 << ", 2:dz = " << track2->helix().dz() << std::endl;
01107 #endif
01108 if(m_param.DETERMINE_ONE_TRACK){
01109 if(fabs(track1->helix().dz()) < fabs(track2->helix().dz())){
01110 m_tracks.remove(track2);
01111 if(merge3DTrack(track1, tracks))
01112 if(check3DTrack(track1) &&
01113 trace3DTrack(track1)){
01114 ok3d[0] = 1;
01115 ok3d[1] = 1;
01116 mask3DTrack(track1,tmp);
01117 tmp.append(track1->links());
01118 }
01119 }else{
01120 m_tracks.remove(track1);
01121 if(merge3DTrack(track2, tracks))
01122 if(check3DTrack(track2) &&
01123 trace3DTrack(track2)){
01124 ok3d[0] = 1;
01125 ok3d[1] = 1;
01126 mask3DTrack(track2,tmp);
01127 tmp.append(track2->links());
01128 }
01129 }
01130 }else{
01131 int isSaved[2] = { 0, 0 };
01132 if(merge3DTrack(track1, tracks)){
01133 if(check3DTrack(track1) &&
01134 trace3DTrack(track1)){
01135 ok3d[0] = 1;
01136 mask3DTrack(track1,tmp);
01137 tmp.append(track1->links());
01138 isSaved[0] = 1;
01139 }
01140 }
01141 if(merge3DTrack(track2, tracks)){
01142 if(check3DTrack(track2) &&
01143 trace3DTrack(track2)){
01144 ok3d[1] = 1;
01145 mask3DTrack(track2,tmp);
01146 tmp.append(track2->links());
01147 isSaved[1] = 1;
01148 }
01149 }
01150 if(isSaved[0] == 1 && isSaved[1] == 1){
01151 track1->daughter(track2);
01152 track2->daughter(track1);
01153 }
01154 }
01155 }else if(track1){
01156 #if DEBUG_CURL_DUMP
01157 std::cout << "(TCurlFinder) 3D:Create Track!!! in track1" << std::endl;
01158 #endif
01159 salvage3DTrack(track1,true);
01160 if(merge3DTrack(track1, tracks))
01161 if(check3DTrack(track1) &&
01162 trace3DTrack(track1)){
01163 ok3d[0] = 1;
01164 mask3DTrack(track1,tmp);
01165 tmp.append(track1->links());
01166 }
01167 }else if(track2){
01168 #if DEBUG_CURL_DUMP
01169 std::cout << "(TCurlFinder) 3D:Create Track!!! in track2" << std::endl;
01170 #endif
01171 salvage3DTrack(track2,true);
01172 if(merge3DTrack(track2, tracks))
01173 if(check3DTrack(track2) &&
01174 trace3DTrack(track2)){
01175 ok3d[1] = 1;
01176 mask3DTrack(track2,tmp);
01177 tmp.append(track2->links());
01178 }
01179 }
01180
01181 if(m_param.OUTPUT_2DTRACKS){
01182
01183 if(ok2d[0] == 1 && ok3d[0] == 0){
01184 removeStereo(*circle);
01185 double chi2_2d;
01186 int ndf_2d;
01187 if(fitWDD(*circle,chi2_2d,ndf_2d)){
01188 TTrack *trk2d = new TTrack(*circle);
01189 trk2d->_ndf = ndf_2d;
01190 trk2d->_chi2 = chi2_2d;
01191 m_2dTracks.append(trk2d);
01192 m_allTracks.append(trk2d);
01193 }
01194 #if DEBUG_CURL_DUMP
01195 else{
01196 std::cout << "(TCurlFinder) 2D:fit with drift information!!!" << std::endl;
01197 }
01198 #endif
01199 }
01200 if(ok2d[1] == 1 && ok3d[1] == 0){
01201 removeStereo(*dividedCircle);
01202 double chi2_2d;
01203 int ndf_2d;
01204 if(fitWDD(*dividedCircle,chi2_2d,ndf_2d)){
01205 TTrack *trk2d = new TTrack(*dividedCircle);
01206 trk2d->_ndf = ndf_2d;
01207 trk2d->_chi2 = chi2_2d;
01208 m_2dTracks.append(trk2d);
01209 m_allTracks.append(trk2d);
01210 }
01211 #if DEBUG_CURL_DUMP
01212 else{
01213 std::cout << "(TCurlFinder) 2D:fit with drift information!!!" << std::endl;
01214 }
01215 #endif
01216 }
01217 }
01218
01219 }else{
01220 #if DEBUG_CURL_DUMP
01221 std::cout << "(TCurlFinder) 2D:dividing...no good...1 Circles!!!" << std::endl;
01222 #endif
01223 int ok2d = 0;
01224 int ok3d = 0;
01225
01226
01227 if(trace2DTrack(circle) &&
01228 check2DCircle(circle) &&
01229 circle->fitForCurl(1) != -1){
01230 #if DEBUG_CURL_DUMP
01231 std::cout << "(TCurlFinder) 2D:Success Circle Fit!!!" << std::endl;
01232 #endif
01233 ok2d = 1;
01234 TTrack *track3 = make3DTrack(circle, segmentList);
01235 if(track3){
01236 #if DEBUG_CURL_DUMP
01237 std::cout << "(TCurlFinder) 3D:Create Track!!! in track3" << std::endl;
01238 #endif
01239 salvage3DTrack(track3,true);
01240 if(merge3DTrack(track3, tracks))
01241 if(check3DTrack(track3) &&
01242 trace3DTrack(track3)){
01243 ok3d = 1;
01244 mask3DTrack(track3,tmp);
01245 tmp.append(track3->links());
01246 }
01247 }
01248
01249 if(m_param.OUTPUT_2DTRACKS){
01250
01251 if(ok2d == 1 && ok3d == 0){
01252 removeStereo(*circle);
01253 double chi2_2d;
01254 int ndf_2d;
01255 if(fitWDD(*circle,chi2_2d,ndf_2d)){
01256 TTrack *trk2d = new TTrack(*circle);
01257 trk2d->_ndf = ndf_2d;
01258 trk2d->_chi2 = chi2_2d;
01259 m_2dTracks.append(trk2d);
01260 m_allTracks.append(trk2d);
01261 }
01262 #if DEBUG_CURL_DUMP
01263 else{
01264 std::cout << "(TCurlFinder) 2D:fit with drift information!!!" << std::endl;
01265 }
01266 #endif
01267 }
01268 }
01269 }
01270 }
01271 m_unusedAxialHits.remove(tmp);
01272 m_unusedStereoHits.remove(tmp);
01273 for(unsigned ii=0, nsize=m_segmentList.length();ii<nsize;++ii){
01274 m_segmentList[ii]->remove(tmp);
01275 if(m_segmentList[ii]->list().length() < m_param.MIN_SEGMENT)
01276 m_segmentList[ii]->removeAll();
01277 }
01278 }
01279 segmentList[i]->removeAll();
01280 }
01281
01282
01283
01284
01285
01286 assignTracks();
01287
01288 #if DEBUG_CURL_DUMP
01289 std::cout << "(TCurlFinder)MDC Rec Track # 3D = " << m_tracks.length()
01290 << ", 2D = " << m_2dTracks.length() << std::endl;
01291 std::cout << "3D Track List" << std::endl;
01292 for(int j=0;j<m_tracks.length();++j){
01293 m_tracks[j]->setFinderType(3);
01294 unsigned nA = 0, nS = 0;
01295 unsigned nAOK = 0, nSOK = 0;
01296 for(unsigned i=0,size=m_tracks[j]->nLinks();i<size;++i){
01297 if(m_tracks[j]->links()[i]->wire()->stereo())++nS;
01298 else ++nA;
01299 if(
01300 (m_tracks[j]->links()[i]->hit()->state() & WireHitFittingValid)){
01301 if(m_tracks[j]->links()[i]->wire()->stereo())++nSOK;
01302 else ++nAOK;
01303 }
01304 }
01305 std::cout << "(TCurlFinder) #" << j << ": wire info...A+S: " << m_tracks[j]->nLinks()
01306 << ", A: " << nAOK << "/" << nA
01307 << ", S: " << nSOK << "/" << nS << std::endl;
01308 if(m_tracks[j]->daughter())
01309 std::cout << "(TCurlFinder) Relation = EXIST" << std::endl;
01310 else
01311 std::cout << "(TCurlFinder) Relation = NO EXIST" << std::endl;
01312 }
01313 std::cout << "2D Track List" << std::endl;
01314 for(int j=0;j<m_2dTracks.length();++j){
01315 unsigned nA = 0, nS = 0;
01316 unsigned nAOK = 0, nSOK = 0;
01317 for(unsigned i=0,size=m_2dTracks[j]->nLinks();i<size;++i){
01318 if(m_2dTracks[j]->links()[i]->wire()->stereo())++nS;
01319 else ++nA;
01320 if(
01321 (m_2dTracks[j]->links()[i]->hit()->state() & WireHitFittingValid)){
01322 if(m_2dTracks[j]->links()[i]->wire()->stereo())++nSOK;
01323 else ++nAOK;
01324 }
01325 }
01326 std::cout << "(TCurlFinder) #" << j << ": wire info...A+S: " << m_2dTracks[j]->nLinks()
01327 << ", A: " << nAOK << "/" << nA
01328 << ", S: " << nSOK << "/" << nS
01329 << ", Chi2: " << m_2dTracks[j]->chi2()
01330 << ", Ndf: " << m_2dTracks[j]->ndf() << std::endl;
01331 if(m_2dTracks[j]->daughter())
01332 std::cout << "(TCurlFinder) Relation = EXIST" << std::endl;
01333 else
01334 std::cout << "(TCurlFinder) Relation = NO EXIST" << std::endl;
01335 }
01336 #endif
01337
01338
01339 m_allTracks.remove(m_tracks);
01340 checkRelation(m_tracks);
01341 tracks.append(m_tracks);
01342 if(m_param.OUTPUT_2DTRACKS){
01343
01344 m_allTracks.remove(m_2dTracks);
01345 tracks2D.append(m_2dTracks);
01346 }
01347 }
01348
01349 void
01350 TCurlFinder::check2DTracks(void)
01351 {
01352 if(m_2dTracks.length() == 0)return;
01353 AList<TMLink> allWires_3Dtrks;
01354 for(int i=0;i<m_tracks.length();++i){
01355 allWires_3Dtrks.append(m_tracks[i]->links());
01356 }
01357
01358
01359 for(int i=0;i<m_2dTracks.length();++i){
01360 AList<TMLink> usedWires;
01361 for(int j=0;j<m_2dTracks[i]->nLinks();++j){
01362 int ok = 1;
01363 for(int k=0;k<allWires_3Dtrks.length();++k){
01364 if(m_2dTracks[i]->links()[j]->wire()->id() ==
01365 allWires_3Dtrks[k]->wire()->id()){
01366 ok = 0;
01367 break;
01368 }
01369 }
01370 if(ok == 0){
01371 usedWires.append(m_2dTracks[i]->links()[j]);
01372 }
01373 }
01374
01375
01376 m_2dTracks[i]->remove(usedWires);
01377 }
01378 }
01379
01380 void
01381 TCurlFinder::checkRelation(AList<TTrack> &list) {
01382 unsigned nT = list.length();
01383 if(nT <= 1)return;
01384 for(unsigned i=0;i<nT;++i){
01385 if(list[i]->daughter()){
01386 int isHere = 0;
01387 for(unsigned j=0;j<nT;++j){
01388 if(i != j &&
01389 list[i]->daughter() == list[j]){
01390 isHere = 1;
01391 break;
01392 }
01393 }
01394 if(isHere == 0){
01395 list[i]->daughter(NULL);
01396 }
01397 }
01398 }
01399 }
01400
01401 TCircle *
01402 TCurlFinder::dividing2DTrack(TCircle *circle) {
01403 AList<TMLink> positive, negative;
01404 for(unsigned i = 0, size = circle->nLinks(); i < size; ++i){
01405 if(circle->center().x()*circle->links()[i]->hit()->wire()->xyPosition().y() -
01406 circle->center().y()*circle->links()[i]->hit()->wire()->xyPosition().x() > 0.){
01407 positive.append(circle->links()[i]);
01408 }else{
01409 negative.append(circle->links()[i]);
01410 }
01411 }
01412 if(positive.length() > negative.length()){
01413 circle->remove(negative);
01414 circle->property(1.,fabs(circle->radius()),circle->center());
01415 if(negative.length() >= 3){
01416 TCircle *new_circle = new TCircle(negative);
01417 m_allCircles.append(new_circle);
01418 new_circle->property(-1.,-1.*fabs(circle->radius()),circle->center());
01419 return new_circle;
01420 }else{
01421 return NULL;
01422 }
01423 }else{
01424 circle->remove(positive);
01425 circle->property(-1.,-1.*fabs(circle->radius()),circle->center());
01426 if(positive.length() >= 3){
01427 TCircle *new_circle = new TCircle(positive);
01428 m_allCircles.append(new_circle);
01429 new_circle->property(1.,fabs(circle->radius()),circle->center());
01430 return new_circle;
01431 }else{
01432 return NULL;
01433 }
01434 }
01435 }
01436
01437 void
01438 TCurlFinder::assignTracks(void) {
01439
01440 for(int i=0,size=m_tracks.length();i<size;++i) {
01441 m_tracks[i]->assign(WireHitCurlFinder);
01442 m_tracks[i]->finder(TrackCurlFinder);
01443
01444 }
01445
01446
01447 for(int i=0,size=m_2dTracks.length();i<size;++i) {
01448 m_2dTracks[i]->assign(WireHitCurlFinder);
01449 m_2dTracks[i]->finder(TrackCurlFinder);
01450 }
01451 }
01452
01453 extern "C" int
01454 TCurlFinder_doubleCompare(const void *i, const void *j){
01455 if(*(static_cast<const double*>(i)) > *(static_cast<const double*>(j)))return 1;
01456 if(*(static_cast<const double*>(i)) < *(static_cast<const double*>(j)))return -1;
01457 return 0;
01458 }
01459
01460 int
01461 TCurlFinder::trace2DTrack(TCircle *circle) {
01462
01463 unsigned nSize = circle->links().length();
01464 if(nSize == 0)return 0;
01465 double r = fabs(circle->radius());
01466 if(r < 0.01)return 0;
01467 double cx = circle->center().x();
01468 double cy = circle->center().y();
01469 double th = atan2(-cy,-cx);
01470 if(th < 0.)th += 2.*M_PI;
01471
01472 unsigned innerOK = 0;
01473 double *angle = new double [circle->links().length()];
01474 for(unsigned i=0,size=nSize;i<size;++i){
01475 double th_r = atan2(circle->links()[i]->wire()->xyPosition().y()-cy,
01476 circle->links()[i]->wire()->xyPosition().x()-cx);
01477 if(th_r < 0.)th_r += 2.*M_PI;
01478 double diff = th_r-th+2.*M_PI;
01479 if(th_r > th)diff = th_r-th;
01480 if(circle->links()[i]->wire()->superLayerId() <= 2)innerOK = 1;
01481 angle[i] = diff;
01482 }
01483 qsort(angle,nSize,sizeof(double),TCurlFinder_doubleCompare);
01484 double maxDiffAngle = 0.;
01485 unsigned maxIndex = 0;
01486 for(unsigned i=0,size=nSize;i<size-1;++i){
01487 if(angle[i+1]-angle[i] > maxDiffAngle){
01488 maxDiffAngle = angle[i+1]-angle[i];
01489 maxIndex = i;
01490 }
01491 }
01492 delete [] angle;
01493
01494
01495
01496
01497 if(innerOK == 1)return 1;
01498
01499 double q = circle->radius() > 0. ? 1. : -1;
01500 for(unsigned i=0,size=m_hitsOnInnerSuperLayer.length();i<size;++i){
01501 double mag = distance(m_hitsOnInnerSuperLayer[i]->wire()->xyPosition().x()-cx,
01502 m_hitsOnInnerSuperLayer[i]->wire()->xyPosition().y()-cy);
01503 if(fabs(mag-r) < m_param.TRACE2D_FIRST_SUPERLAYER &&
01504 q*(cx*m_hitsOnInnerSuperLayer[i]->wire()->xyPosition().y()-
01505 cy*m_hitsOnInnerSuperLayer[i]->wire()->xyPosition().x()) > 0.){
01506 return 1;
01507 }
01508 }
01509 return 0;
01510 }
01511
01512 bool
01513 TCurlFinder::check2DCircle(TCircle *circle) {
01514 unsigned nA(nAxialHits(fabs(circle->radius())*2.0));
01515
01516 unsigned nMA = static_cast<unsigned>(floor(m_param.RATIO_USED_WIRE*static_cast<double>(nA)));
01517 if(nMA < 3)nMA = 3;
01518
01519 unsigned nAhits(0), nShits(0);
01520 for(unsigned i=0,size=circle->nLinks();i<size;++i){
01521 if((circle->links()[i])->wire()->axial())++nAhits;
01522 else ++nShits;
01523 }
01524 #if DEBUG_CURL_DUMP
01525 if(nAhits < nMA){
01526 std::cout << "(TCurlFinder) 2D:Fail...checking axial wires # = "
01527 << nAhits << " < " << nMA << std::endl;
01528 }
01529 #endif
01530 if(nAhits >= nMA)return true;
01531 return false;
01532 }
01533
01534 bool
01535 TCurlFinder::check3DTrack(TTrack *track) {
01536 trace3DTrack(track);
01537 unsigned nA = 0, nS = 0;
01538 for(unsigned i=0,size=track->nLinks();i<size;++i){
01539 if(!(track->links()[i]->hit()->state() & WireHitFittingValid))continue;
01540 if(track->links()[i]->wire()->stereo())++nS;
01541 else ++nA;
01542 if(nA >= 3 && nS >= 2)return true;
01543 }
01544 m_tracks.remove(track);
01545 #if DEBUG_CURL_DUMP
01546 std::cout << "(TCurlFinder) 3D:Checked...Fail...removing this track. Valid Axial # = "
01547 << nA << ", Stereo # = " << nS << std::endl;
01548 #endif
01549 return false;
01550 }
01551
01552 int
01553 TCurlFinder::trace3DTrack(TTrack *track) {
01554
01555 unsigned nSize = track->links().length();
01556 if(nSize == 0){
01557 m_tracks.remove(track);
01558 return 0;
01559 }
01560 double r = fabs(track->helix().radius());
01561 if(r < 0.01){
01562 m_tracks.remove(track);
01563 return 0;
01564 }
01565 double cx = track->helix().center().x();
01566 double cy = track->helix().center().y();
01567 double th = atan2(-cy,-cx);
01568 if(th < 0.)th += 2.*M_PI;
01569
01570 double *angle = new double [track->links().length()];
01571 for(unsigned i=0,size=nSize;i<size;++i){
01572 double th_r = atan2(track->links()[i]->positionOnTrack().y()-cy,
01573 track->links()[i]->positionOnTrack().x()-cx);
01574 if(th_r < 0.)th_r += 2.*M_PI;
01575 double diff = th_r-th+2.*M_PI;
01576 if(th_r > th)diff = th_r-th;
01577 angle[i] = diff;
01578 }
01579 qsort(angle,nSize,sizeof(double),TCurlFinder_doubleCompare);
01580 double maxDiffAngle = 0.;
01581 unsigned maxIndex = 0;
01582 for(unsigned i=0,size=nSize;i<size-1;++i){
01583 if(angle[i+1]-angle[i] > maxDiffAngle){
01584 maxDiffAngle = angle[i+1]-angle[i];
01585 maxIndex = i;
01586 }
01587 }
01588 delete [] angle;
01589
01590
01591 if(r*maxDiffAngle > m_param.TRACE3D_DISTANCE){
01592 m_tracks.remove(track);
01593 return 0;
01594 }else{
01595 return 1;
01596 }
01597 }
01598
01599 void
01600 TCurlFinder::mask3DTrack(TTrack *track,
01601 AList<TMLink> &maskList) {
01602 double r(fabs(track->helix().radius()));
01603 double cx(track->helix().center().x());
01604 double cy(track->helix().center().y());
01605
01606 AList<TMLink> list(m_unusedAxialHits);
01607 list.append(m_unusedStereoHits);
01608 list.remove(track->links());
01609 list.sort(SortByWireId);
01610
01611 AList<TMLink> removeList;
01612 for(unsigned i=0,size=list.length();i<size;++i){
01613 double d = distance(*track, *(list[i]));
01614 if(d < m_param.MASK_DISTANCE){
01615 HepPoint3D tmp(d, 0., 0.);
01616 list[i]->position(tmp);
01617 removeList.append(list[i]);
01618 }
01619 }
01620
01621 int pLayerId1 = static_cast<int>(layerId(2.*r));
01622 if(pLayerId1 != 50)pLayerId1 -= 1;
01623 int pLayerId2 = pLayerId1+2;
01624
01625 AList<TMLink> preCand, cand;
01626 while(removeList.length()){
01627 preCand.removeAll();
01628 preCand.append(removeList[0]);
01629 if(removeList.length() >= 2){
01630 for(unsigned j=1,size=removeList.length();j<size;++j){
01631 if(removeList[0]->wire()->layerId() == removeList[j]->wire()->layerId()){
01632 for(unsigned k=0,num=preCand.length();k<num;++k){
01633 if(preCand[k]->wire()->localIdForPlus()+1 == removeList[j]->wire()->localId()){
01634 preCand.append(removeList[j]);
01635 break;
01636 }
01637 }
01638 }
01639 }
01640 #if 1
01641
01642 if(preCand[0]->wire()->layerId() >= pLayerId1 &&
01643 preCand[0]->wire()->layerId() <= pLayerId2){
01644 cand.append(preCand);
01645 }else if(preCand.length() == 2){
01646 cand.append(preCand);
01647 }else if(preCand.length() == 1){
01648 cand.append(preCand[0]);
01649 }
01650 #else
01651 if(preCand.length() == 1){
01652 if(preCand[0]->position().x() < MASK_DISTANCE)cand.append(preCand[0]);
01653 }else{
01654 if(preCand[0]->wire()->layerId() >= pLayerId1 &&
01655 preCand[0]->wire()->layerId() <= pLayerId2){
01656 cand.append(preCand);
01657 }else if(preCand.length() == 2){
01658 cand.append(preCand);
01659 }
01660 }
01661 #endif
01662 }else{
01663 cand.append(removeList[0]);
01664 }
01665 removeList.remove(removeList[0]);
01666 removeList.remove(cand);
01667 }
01668 maskList.append(cand);
01669 }
01670
01671 TTrack*
01672 TCurlFinder::merge3DTrack(TTrack *track, AList<TTrack> &confTracks) {
01673 if(!m_param.MERGE_EXE)return track;
01674
01675 AList<TTrack> tracks(confTracks);
01676 tracks.append(m_tracks);
01677 tracks.remove(track);
01678 if(tracks.length() == 0)return track;
01679
01680 double r = track->helix().radius();
01681 double cx = track->helix().center().x();
01682 double cy = track->helix().center().y();
01683
01684 unsigned bestIndex = 0;
01685 double bestDiff = 1.0e+20;
01686 double R, cX, cY;
01687 for(unsigned i=0,size=tracks.length();i<size;++i){
01688 R = fabs(tracks[i]->helix().radius());
01689 cX = tracks[i]->helix().center().x();
01690 cY = tracks[i]->helix().center().y();
01691 if(fabs(r)*(1.-m_param.MERGE_RATIO) <= R && R <= fabs(r)*(1.+m_param.MERGE_RATIO)){
01692 if(fabs(cx-cX) <= fabs(r)*m_param.MERGE_RATIO && fabs(cy-cY) <= fabs(r)*m_param.MERGE_RATIO){
01693 double diff = fabs((fabs(r)-fabs(R))*(cx-cX)*(cy-cY));
01694 if(diff < bestDiff){
01695 bestDiff = diff;
01696 bestIndex = i;
01697 }
01698 }
01699 }
01700 }
01701 if(bestDiff == 1.0e20)return track;
01702 R = tracks[bestIndex]->helix().radius();
01703 cX = tracks[bestIndex]->helix().center().x();
01704 cY = tracks[bestIndex]->helix().center().y();
01705 if(r*R >= 0.){
01706 if(fabs(track->helix().dz()-tracks[bestIndex]->helix().dz()) < m_param.MERGE_Z_DIFF){
01707 if(track->nLinks() > tracks[bestIndex]->nLinks()){
01708 m_tracks.remove(tracks[bestIndex]);
01709 return track;
01710 }else{
01711 m_tracks.remove(track);
01712 #if DEBUG_CURL_DUMP
01713 std::cout << "(TCurlFinder) 3D:Merged...removing this track.(type1)" << std::endl;
01714 #endif
01715 return NULL;
01716 }
01717 }
01718 }else{
01719 bool newTrack(false), oldTrack(false);
01720 unsigned newCounter(0), oldCounter(0);
01721 for(unsigned i=0, size=m_hitsOnInnerSuperLayer.length();
01722 i<size;++i){
01723 if(!oldTrack){
01724 if((R > 0. &&
01725 cX*(m_hitsOnInnerSuperLayer[i]->wire()->xyPosition().y()-cY)-
01726 cY*(m_hitsOnInnerSuperLayer[i]->wire()->xyPosition().x()-cX) > 0.) ||
01727 (R < 0. &&
01728 cX*(m_hitsOnInnerSuperLayer[i]->wire()->xyPosition().y()-cY)-
01729 cY*(m_hitsOnInnerSuperLayer[i]->wire()->xyPosition().x()-cX) < 0.)){
01730 double dist = distance(m_hitsOnInnerSuperLayer[i]->wire()->xyPosition().x()-cX,
01731 m_hitsOnInnerSuperLayer[i]->wire()->xyPosition().y()-cY);
01732 if(dist < fabs(R)){
01733 if(fabs(fabs(R)-dist-m_hitsOnInnerSuperLayer[i]->drift()) < 0.5){
01734 ++oldCounter;
01735 if(oldCounter >= 3)oldTrack = true;
01736 }
01737 }else{
01738 if(fabs(dist-fabs(R)-m_hitsOnInnerSuperLayer[i]->drift()) < 0.5){
01739 ++oldCounter;
01740 if(oldCounter >= 3)oldTrack = true;
01741 }
01742 }
01743 }
01744 }
01745 if(!newTrack){
01746 if((r > 0. &&
01747 cx*(m_hitsOnInnerSuperLayer[i]->wire()->xyPosition().y()-cy)-
01748 cy*(m_hitsOnInnerSuperLayer[i]->wire()->xyPosition().x()-cx) > 0.) ||
01749 (r < 0. &&
01750 cx*(m_hitsOnInnerSuperLayer[i]->wire()->xyPosition().y()-cy)-
01751 cy*(m_hitsOnInnerSuperLayer[i]->wire()->xyPosition().x()-cx) < 0.)){
01752 double dist = distance(m_hitsOnInnerSuperLayer[i]->wire()->xyPosition().x()-cx,
01753 m_hitsOnInnerSuperLayer[i]->wire()->xyPosition().y()-cy);
01754 if(dist < fabs(r)){
01755 if(fabs(fabs(r)-dist-m_hitsOnInnerSuperLayer[i]->drift()) < 0.5){
01756 ++newCounter;
01757 if(newCounter >= 3)newTrack = true;
01758 }
01759 }else{
01760 if(fabs(dist-fabs(r)-m_hitsOnInnerSuperLayer[i]->drift()) < 0.5){
01761 ++newCounter;
01762 if(newCounter >= 3)newTrack = true;
01763 }
01764 }
01765 }
01766 }
01767 if(oldTrack && newTrack)break;
01768 }
01769 if(oldTrack && !newTrack){
01770 m_tracks.remove(track);
01771 #if DEBUG_CURL_DUMP
01772 std::cout << "(TCurlFinder) 3D:Merged...removing this track.(type2)" << std::endl;
01773 #endif
01774 return NULL;
01775 }else if(!oldTrack && newTrack){
01776 m_tracks.remove(tracks[bestIndex]);
01777 return track;
01778 }else if(!oldTrack && !newTrack){
01779 m_tracks.remove(track);
01780 m_tracks.remove(tracks[bestIndex]);
01781 #if DEBUG_CURL_DUMP
01782 std::cout << "(TCurlFinder) 3D:Merged...removing this track.(type3)" << std::endl;
01783 #endif
01784 return NULL;
01785 }else if(oldTrack && newTrack){
01786 if(fabs(track->helix().dz()) > fabs(tracks[bestIndex]->helix().dz()) &&
01787 fabs(track->helix().dz()) > fabs(tracks[bestIndex]->helix().dz())+m_param.MERGE_Z_DIFF){
01788 m_tracks.remove(track);
01789 #if DEBUG_CURL_DUMP
01790 std::cout << "(TCurlFinder) 3D:Merged...removing this track.(type4)" << std::endl;
01791 #endif
01792 return NULL;
01793 }else if(fabs(tracks[bestIndex]->helix().dz()) > fabs(track->helix().dz()) &&
01794 fabs(tracks[bestIndex]->helix().dz()) > fabs(track->helix().dz())+m_param.MERGE_Z_DIFF){
01795 m_tracks.remove(tracks[bestIndex]);
01796 return track;
01797 }
01798 }
01799 }
01800 return track;
01801 }
01802
01803 void
01804 TCurlFinder::salvage3DTrack(TTrack *track, bool half) {
01805 if(track->nLinks() >= m_param.MIN_SALVAGE){
01806 AList<TMLink> list = m_unusedAxialHits;
01807 list.append(m_unusedStereoHits);
01808 list.remove(track->links());
01809
01810 if(half){
01811 double q = track->charge();
01812 double x = track->helix().center().x();
01813 double y = track->helix().center().y();
01814 AList<TMLink> removeList;
01815 const double Bz = -1000*m_pmgnIMF->getReferField();
01816 for(unsigned i = 0, size = list.length(); i < size; ++i){
01817 if(q*Bz*(x*list[i]->wire()->xyPosition().y()-y*list[i]->wire()->xyPosition().x())<0.){
01818
01819 #ifdef TRKRECO_DEBUG_DETAIL
01820 std::cout<<" "<<__FILE__<<" "<<__LINE__<<" removeList append "<<i<<std::endl;
01821 #endif
01822 removeList.append(list[i]);
01823 }
01824 }
01825 list.remove(removeList);
01826 }
01827
01828 AList<TMLink> badCand, goodCand;
01829 double dist;
01830 for(unsigned j=0, nLinks=track->nLinks();j<nLinks;++j){
01831 dist = distance(*track, *(track->links()[j]));
01832 #ifdef TRKRECO_DEBUG_DETAIL
01833 std::cout<<" "<<__FILE__<<" "<<__LINE__<<" "<<j<<" distance "<<dist<<std::endl;
01834 #endif
01835 if(dist > m_param.BAD_DISTANCE_FOR_SALVAGE)badCand.append(track->links()[j]);
01836 }
01837 for(unsigned j=0, nList=list.length();j<nList;++j){
01838 dist = distance(*track, *(list[j]));
01839 #ifdef TRKRECO_DEBUG_DETAIL
01840 std::cout<<" "<<__FILE__<<" "<<__LINE__<<" "<<j<<" distance "<<dist<<std::endl;
01841 #endif
01842 if(dist < m_param.GOOD_DISTANCE_FOR_SALVAGE)goodCand.append(list[j]);
01843 }
01844 track->TTrackBase::remove(badCand);
01845 track->TTrackBase::append(goodCand);
01846 if(m_fitter.fit(*track) < 0){
01847 track->TTrackBase::remove(goodCand);
01848 track->TTrackBase::append(badCand);
01849 m_fitter.fit(*track);
01850 }
01851 }
01852 return;
01853 }
01854
01855 double
01856 TCurlFinder::distance(const TTrack &track, const TMLink &link) const {
01857
01858 if((link.wire()->superLayerId() > 1 && link.wire()->superLayerId() <5) || link.wire()->superLayerId() >8 ){
01859
01860 double d = distance(track.helix().center().x()-link.xyPosition().x(),
01861 track.helix().center().y()-link.xyPosition().y());
01862 double diff = fabs(d - fabs(track.helix().radius()));
01863 return fabs(link.hit()->drift()-diff);
01864 }
01865
01866 HepPoint3D xc(track.helix().center());
01867 HepPoint3D xw(link.xyPosition());
01868 HepPoint3D xt(track.helix().x());
01869 HepVector3D v0(xt-xc), v1(xw-xc);
01870 double vCrs(v0.x() * v1.y() - v0.y() * v1.x());
01871 double vDot(v0.x() * v1.x() + v0.y() * v1.y());
01872 double dPhi = atan2(vCrs, vDot);
01873 Vector a(track.helix().a());
01874 double kappa = a[2];
01875 double phi0 = a[1];
01876
01877 const double Bz = -1000*m_pmgnIMF->getReferField();
01878
01879 double rho = m_param.ALPHA_SAME_WITH_HELIX/kappa/Bz;
01880 double tanLambda = a[4];
01881 HepVector3D v = link.wire()->direction();
01882 Vector c(3);
01883 c = HepPoint3D(link.wire()->backwardPosition()-(v*link.wire()->backwardPosition())*v);
01884
01885 HepDiagMatrix e(3,1);
01886 Matrix t(3, 3);
01887 t[0][0] = v.x() * v.x();
01888 t[0][1] = v.x() * v.y();
01889 t[0][2] = v.x() * v.z();
01890 t[1][0] = t[0][1];
01891 t[1][1] = v.y() * v.y();
01892 t[1][2] = v.y() * v.z();
01893 t[2][0] = t[0][2];
01894 t[2][1] = t[1][2];
01895 t[2][2] = v.z() * v.z();
01896 t -= e;
01897
01898 double factor = 1.;
01899 unsigned nTrial = 0;
01900
01901
01902 Vector x(3);
01903 Vector dXdPhi(3);
01904 Vector d2Xd2Phi(3);
01905 double fOld = 0.;
01906 const double convergence = 1.0e-5;
01907 while(nTrial < 100){
01908 x = track.helix().x(dPhi);
01909 double cosPhi = cos(phi0+dPhi);
01910 double sinPhi = sin(phi0+dPhi);
01911 dXdPhi[0] = rho*sinPhi;
01912 dXdPhi[1] = - rho*cosPhi;
01913 dXdPhi[2] = - rho*tanLambda;
01914
01915
01916 double f = dot(c,dXdPhi)+dot(x,(t*dXdPhi));
01917
01918 if(fabs(f) < convergence)break;
01919 if(nTrial > 0){
01920 double eval = (1.-0.25*factor)*fabs(fOld)-fabs(f);
01921 if(eval <= 0.)factor *= 0.5;
01922 }
01923
01924 d2Xd2Phi[0] = rho*cosPhi;
01925 d2Xd2Phi[1] = rho*sinPhi;
01926 d2Xd2Phi[2] = 0.;
01927 double df = dot(c, d2Xd2Phi)+
01928 dot(dXdPhi, (t * dXdPhi))+
01929 dot(x, (t * d2Xd2Phi));
01930 dPhi -= factor * f / df;
01931
01932 fOld = f;
01933 ++nTrial;
01934 }
01935
01936 double beta = v*(track.helix().x(dPhi)-link.wire()->backwardPosition());
01937 return fabs((link.wire()->backwardPosition()+beta*v-track.helix().x(dPhi)).mag()-
01938 link.hit()->drift());
01939 }
01940
01941
01942
01943
01944
01945 TCircle *
01946 TCurlFinder::make2DTrack(const AList<TMLink> &seed,
01947 const AList<TSegmentCurl> &segmentList,
01948 const unsigned ip) {
01949
01950 if(seed.length() < m_param.minimum_seedLength)return NULL;
01951 TCircle *circle = new TCircle(seed);
01952 m_allCircles.append(circle);
01953 int errorFlag = circle->fitForCurl(ip);
01954
01955
01956 if( (fabs(circle->radius()) > m_param.MIN_RADIUS_OF_STRANGE_TRACK) || (fabs(circle->radius()) < m_param.ULTIMATE_MIN_RADIUS_OF_STRANGE_TRACK) )return NULL;
01957 int searchDirection = 1;
01958 int searchPath(searchDirection);
01959 bool searchZero(false);
01960 bool changeDirection(false);
01961 unsigned superLayerId = seed[0]->hit()->wire()->superLayerId();
01962
01963 AList<TMLink> cand, tmpList;
01964 AList<TMLink> preAxialCand, preStereoCand;
01965 makeList(tmpList, segmentList, seed);
01966 for(unsigned i=0, size=tmpList.length();i<size;++i){
01967 if(tmpList[i]->wire()){
01968 if(tmpList[i]->wire()->axial()) preAxialCand.append(tmpList[i]);
01969 else preStereoCand.append(tmpList[i]);
01970 }
01971 }
01972 #if DEBUG_CURL_DUMP
01973 std::cout << "(TCurlFinder) 2D: Superlayer of seed = " << superLayerId << std::endl;
01974 #endif
01975
01976 bool appendFlag = false;
01977 nextStep:
01978 #if DEBUG_CURL_DUMP
01979 std::cout << "(TCurlFinder) 2D: SearchPath = " << searchPath
01980 << " Search SelfSuperlayer = " << (int)(searchZero)
01981 << " Change Direction of Search = " << (int)(changeDirection) << std::endl;
01982 #endif
01983 if(preAxialCand.length() == 0 && preStereoCand.length() == 0){
01984 if(circle->links().length() >= 3){
01985 if(m_unusedAxialHits.length() == 0)
01986 if(errorFlag == -1)return NULL;
01987 else return circle;
01988 else goto salvage;
01989 }else{
01990 if(m_unusedAxialHits.length() == 0)return NULL;
01991 else goto salvage;
01992 }
01993 }
01994 searchAxialCand(cand, preAxialCand, circle,
01995 searchPath, superLayerId, m_param.RANGE_FOR_AXIAL_SEARCH);
01996 if(cand.length() > 0){
01997 appendFlag = true;
01998 for(unsigned i = 0, size = cand.length(); i < size; ++i)circle->append(*cand[i]);
01999 errorFlag = circle->fitForCurl(ip);
02000 preAxialCand.remove(circle->links());
02001 if(m_param.STEREO_2DFIND){
02002 searchStereoCand(cand, preStereoCand, circle,
02003 searchPath, superLayerId, m_param.RANGE_FOR_STEREO_SEARCH);
02004 if(cand.length() > 0){
02005 appendFlag = true;
02006 for(unsigned i = 0, size = cand.length(); i < size; ++i)circle->append(*cand[i]);
02007 errorFlag = circle->fitForCurl(ip);
02008 preStereoCand.remove(circle->links());
02009 }
02010 }
02011 if(searchDirection == 1)++searchPath;
02012 else --searchPath;
02013 goto nextStep;
02014 }else{
02015 if(m_param.STEREO_2DFIND){
02016 searchStereoCand(cand, preStereoCand, circle,
02017 searchPath, superLayerId, m_param.RANGE_FOR_STEREO_SEARCH);
02018 if(cand.length() > 0){
02019 appendFlag = true;
02020 for(unsigned i = 0, size = cand.length(); i < size; ++i)circle->append(*cand[i]);
02021 errorFlag = circle->fitForCurl(ip);
02022 preStereoCand.remove(circle->links());
02023 if(searchDirection == 1)++searchPath;
02024 else --searchPath;
02025 goto nextStep;
02026 }else if((searchPath == 1 || searchPath == -1) && !searchZero){
02027 searchPath = 0;
02028 searchZero = true;
02029 goto nextStep;
02030 }else if((searchPath == 1 || searchPath == -1) && searchZero && !changeDirection){
02031 searchPath *= -1;
02032 searchDirection *= -1;
02033 changeDirection = true;
02034 goto nextStep;
02035 }else{
02036 if(circle->links().length() >= 3){
02037 if(m_unusedAxialHits.length() == 0)
02038 if(errorFlag == -1)return NULL;
02039 else return circle;
02040 else goto salvage;
02041 }else{
02042 if(m_unusedAxialHits.length() == 0)return NULL;
02043 else goto salvage;
02044 }
02045 }
02046 }else{
02047 if((searchPath == 1 || searchPath == -1) && !searchZero){
02048 searchPath = 0;
02049 searchZero = true;
02050 goto nextStep;
02051 }else if((searchPath == 1 || searchPath == -1) && searchZero && !changeDirection){
02052 searchPath *= -1;
02053 searchDirection *= -1;
02054 changeDirection = true;
02055 goto nextStep;
02056 }else{
02057 if(circle->links().length() >= 3){
02058 if(m_unusedAxialHits.length() == 0)
02059 if(errorFlag == -1)return NULL;
02060 else return circle;
02061 else goto salvage;
02062 }else{
02063 if(m_unusedAxialHits.length() == 0)return NULL;
02064 else goto salvage;
02065 }
02066 }
02067 }
02068 }
02069
02070 salvage:
02071 cand.removeAll();
02072 searchHits(cand, m_unusedAxialHits, circle, m_param.RANGE_FOR_AXIAL_LAST2D_SEARCH);
02073 if(m_param.STEREO_2DFIND){
02074 searchHits(cand, m_unusedStereoHits, circle, m_param.RANGE_FOR_STEREO_LAST2D_SEARCH);
02075 }
02076 if(checkAppendHits(circle->links(), cand)){
02077 circle->append(cand);
02078 if(circle->nLinks() >= 3)
02079 if(circle->fitForCurl(ip) == -1)return NULL;
02080 else return circle;
02081 else return NULL;
02082 }else if(circle->nLinks() >= 3){
02083 return circle;
02084 }else return NULL;
02085 }
02086
02087 TMLink *
02088 findIsolatedCloseHits(TMLink *link);
02089
02090 void
02091 TCurlFinder::searchAxialCand(AList<TMLink> &cand,
02092 const AList<TMLink> &preCand,
02093 const TCircle *circle,
02094 const int depth,
02095 const unsigned superLayerID,
02096 const double searchError) {
02097 cand.removeAll();
02098 int innerSuperLayerId = nextSuperAxialLayerId(superLayerID, depth);
02099 if(innerSuperLayerId < 0)return;
02100 for(unsigned i = 0, size = preCand.length(); i < size; ++i){
02101 if(preCand[i]->hit()->wire()->superLayerId() ==
02102 (static_cast<unsigned>(innerSuperLayerId))){
02103 #if 0
02104 if(searchHits(preCand[i], circle, searchError))cand.append(preCand[i]);
02105 #else
02106 if(searchHits(preCand[i], circle, searchError)){
02107 cand.remove(preCand[i]);
02108 cand.append(preCand[i]);
02109 TMLink * cand2 = findIsolatedCloseHits(preCand[i]);
02110 if(cand2){
02111 for(unsigned j = 0; j < size; ++j){
02112 if(preCand[j]->wire()->id() == cand2->wire()->id()){
02113 cand.remove(cand2);
02114 cand.append(cand2);
02115 #if 0
02116 std::cout << "Axial Appending....";
02117 std::cout << " layerID = " << cand2->wire()->layerId();
02118 std::cout << " localID = " << cand2->wire()->localId() << std::endl;
02119 if(searchHits(cand2, circle, searchError)){
02120 std::cout << " But this can be added by default!" << std::endl;
02121 }else{
02122 std::cout << " Good!! this cannot be added by default!" << std::endl;
02123 }
02124 #endif
02125 break;
02126 }
02127 }
02128 }
02129 }
02130 #endif
02131 }
02132 }
02133 }
02134
02135 void
02136 TCurlFinder::searchStereoCand(AList<TMLink> &cand,
02137 const AList<TMLink> &preCand,
02138 const TCircle *circle,
02139 const int depth,
02140 const unsigned superLayerID,
02141 const double searchError) {
02142 cand.removeAll();
02143 int innerSuperLayerId = nextSuperStereoLayerId(superLayerID, depth);
02144 if(innerSuperLayerId < 0 || innerSuperLayerId > m_param.SUPERLAYER_FOR_STEREO_SEARCH)return;
02145 for(unsigned i = 0, size = preCand.length(); i < size; ++i){
02146 if(preCand[i]->hit()->wire()->superLayerId() ==
02147 (static_cast<unsigned>(innerSuperLayerId))){
02148 if(searchHits(preCand[i], circle, searchError))cand.append(preCand[i]);
02149 }
02150 }
02151 }
02152
02153 unsigned
02154 TCurlFinder::searchHits(const TMLink *link, const TCircle *circle,
02155 const double searchError) const {
02156
02157
02158
02159
02160 double dist = distance(link->hit()->wire()->xyPosition().x() - circle->center().x(),
02161 link->hit()->wire()->xyPosition().y() - circle->center().y());
02162 double radius = fabs(circle->radius());
02163
02164
02165 if(radius - searchError < dist &&
02166 radius + searchError > dist)return 1;
02167 return 0;
02168 }
02169
02170 unsigned
02171 TCurlFinder::searchHits(AList<TMLink> &cand,
02172 const AList<TMLink> &preCand,
02173 const TCircle *circle,
02174 const double searchError) const {
02175 unsigned numBefore = cand.length();
02176 for(unsigned i = 0, size = preCand.length(); i < size; ++i){
02177 if(searchHits(preCand[i], circle, searchError)){
02178 #if 0
02179 cand.append(preCand[i]);
02180 #else
02181 cand.remove(preCand[i]);
02182 cand.append(preCand[i]);
02183 TMLink * cand2 = findIsolatedCloseHits(preCand[i]);
02184 if(cand2){
02185 for(unsigned j = 0; j < size; ++j){
02186 if(preCand[j]->wire()->id() == cand2->wire()->id()){
02187 cand.remove(cand2);
02188 cand.append(cand2);
02189 break;
02190 }
02191 }
02192 }
02193 #endif
02194 }
02195 }
02196 if(numBefore == cand.length())return 0;
02197 return 1;
02198 }
02199
02200 unsigned
02201 TCurlFinder::checkAppendHits(const AList<TMLink> &link,
02202 AList<TMLink> &cand) const {
02203 if(cand.length() == 0)return 0;
02204 AList<TMLink> tmp;
02205 for(unsigned i = 0, size1 = cand.length(), size2 = link.length(); i < size1; ++i){
02206 for(unsigned j = 0; j < size2; ++j){
02207 if((cand[i])->wire()->id() == (link[j])->wire()->id()){
02208 tmp.append(cand[i]);
02209 break;
02210 }
02211 }
02212 }
02213 cand.remove(tmp);
02214 if(cand.length() > 0)return 1;
02215 return 0;
02216 }
02217
02218 void
02219 TCurlFinder::removeStereo(TCircle &c) const
02220 {
02221 AList<TMLink> stereoList;
02222 for(int i=0;i<c.links().length();++i){
02223 if(c.links()[i]->wire()->stereo()){
02224 stereoList.append(c.links()[i]);
02225 }
02226 }
02227 if(stereoList.length()>0)c.remove(stereoList);
02228 }
02229
02230 bool
02231 TCurlFinder::fitWDD(TCircle &c,
02232 double &chi2,
02233 int &ndf) const
02234 {
02235 if(c.links().length() <= 3)return false;
02236 Lpav circle;
02237
02238 for(int i=0;i<c.links().length();++i){
02239 circle.add_point((c.links()[i])->wire()->xyPosition().x(),
02240 (c.links()[i])->wire()->xyPosition().y(),1.0);
02241 }
02242 circle.add_point(0.,0.,1.0);
02243 if (circle.fit() < 0.0 || circle.kappa() == 0.0) return false;
02244 double xc = circle.center()[0];
02245 double yc = circle.center()[1];
02246 double r = circle.radius();
02247 const int maxIte = 2;
02248 for(int ite=0;ite<maxIte;++ite){
02249 Lpav circle2;
02250 circle2.clear();
02251
02252 for(int i=0;i<c.links().length();++i){
02253 if(!((c.links()[i])->hit()->state() & WireHitFittingValid))continue;
02254 double R = sqrt(((c.links()[i])->wire()->xyPosition().x()-xc)*((c.links()[i])->wire()->xyPosition().x()-xc)+
02255 ((c.links()[i])->wire()->xyPosition().y()-yc)*((c.links()[i])->wire()->xyPosition().y()-yc));
02256 if(R == 0.)continue;
02257 double U = 1./R;
02258 double dir = R > r ? -1. : 1.;
02259 double X = xc+((c.links()[i])->wire()->xyPosition().x()-xc)*U*(R+dir*(c.links()[i])->hit()->drift());
02260 double Y = yc+((c.links()[i])->wire()->xyPosition().y()-yc)*U*(R+dir*(c.links()[i])->hit()->drift());
02261 circle2.add_point(X,Y,1.0);
02262 }
02263 circle2.add_point(0.,0.,1.0);
02264 if (circle2.fit() < 0.0 || circle2.kappa() == 0.0) return false;
02265 xc = circle2.center()[0];
02266 yc = circle2.center()[1];
02267 r = circle2.radius();
02268
02269 }
02270
02271
02272 double totalChi2 = 0.;
02273 int totalNHit = 0;
02274 for(int i=0;i<c.links().length();++i){
02275 if(!((c.links()[i])->hit()->state() & WireHitFittingValid))continue;
02276 double xw = (c.links()[i])->wire()->xyPosition().x();
02277 double yw = (c.links()[i])->wire()->xyPosition().y();
02278 double R = sqrt((xw-xc)*(xw-xc)+(yw-yc)*(yw-yc));
02279 if(R == 0.)continue;
02280 double U = 1./R;
02281 double X = xc+(xw-xc)*U*r;
02282 double Y = yc+(yw-yc)*U*r;
02283 double zlr = xw*Y-yw*X;
02284 unsigned leftRight = zlr > 0. ? WireHitRight : WireHitLeft;
02285 double pChi2 = sqrt((X-xw)*(X-xw)+(Y-yw)*(Y-yw))-(c.links()[i])->hit()->drift();
02286
02287
02288 if((c.links()[i])->hit()->dDrift() != 0.){
02289 pChi2 *= pChi2/((c.links()[i])->hit()->dDrift()*(c.links()[i])->hit()->dDrift());
02290 totalChi2 += pChi2;
02291
02292 ++totalNHit;
02293 }else pChi2 = 1.0e+10;
02294 (c.links()[i])->update(HepPoint3D(X,Y,0.),HepPoint3D(xw,yw,0.),leftRight,pChi2);
02295
02296 }
02297 chi2 = totalChi2;
02298 if(totalNHit <= 3)return false;
02299 ndf = totalNHit-3;
02300
02301 HepPoint3D center(xc,yc,0.);
02302 double charge = 0.;
02303
02304 int qSum = 0;
02305 for(int i=0;i<c.links().length();++i){
02306 TMLink * l = c.links()[i];
02307 if(l == 0)continue;
02308 const TMDCWireHit * h = l->hit();
02309 if(h == 0)continue;
02310 double q = (center.cross(h->xyPosition())).z();
02311 if(q > 0.)qSum += 1;
02312 else qSum -= 1;
02313 }
02314 if(qSum >= 0)charge = +1.;
02315 else charge = -1.;
02316 r *= charge;
02317
02318 c.property(charge,r,center);
02319
02320 return true;
02321 }
02322
02323
02324
02325
02326
02327 TTrack *
02328 TCurlFinder::make3DTrack(const TCircle *circle, AList<TSegmentCurl> &segmentList) {
02329 unsigned size = segmentList.length();
02330 if(TTrack *track = make3DTrack(circle)){
02331 m_tracks.append(track);
02332 #if 0
02333 std::cout << "MDC Helix+Pt: " << track->helix().dr() << ", "
02334 << track->helix().phi0() << ", "
02335 << track->helix().kappa() << ", "
02336 << track->helix().dz() << ", "
02337 << track->helix().tanl()
02338 << ": " << 10000./2.9979258/15./track->helix().kappa() << std::endl;
02339 #endif
02340 return track;
02341 }
02342 return NULL;
02343 }
02344
02345 TTrack *
02346 TCurlFinder::make3DTrack(const TCircle *circle) {
02347 TTrack *track = new TTrack(*circle);
02348 m_allTracks.append(track);
02349
02350
02351 if(track->links().length() < m_param.minimum_2DTrackLength){
02352 #if DEBUG_CURL_DUMP
02353 std::cout << "(TCurlFinder) 3D:Fail...inital hit wire # < 3." << std::endl;
02354 #endif
02355 return NULL;
02356 }
02357 AList<TMLink> allStereoHits(m_unusedStereoHits);
02358 allStereoHits.remove(track->links());
02359 AList<TMLink> closeHits;
02360 findCloseHits(allStereoHits, *track, closeHits);
02361
02362
02363 if ( closeHits.length() < m_param.minimum_closeHitsLength ){
02364 #if DEBUG_CURL_DUMP
02365 std::cout << "(TCurlFinder) 3D:Fail...stereohit wire # < 5." << std::endl;
02366 #endif
02367 return NULL;
02368 }
02369
02370 if(!m_builder.buildStereo(*track, closeHits, m_allStereoHitsOriginal)){
02371 #if DEBUG_CURL_DUMP
02372 std::cout << "(TCurlFinder) 3D:Fail...can not build stereo." << std::endl;
02373 #endif
02374 return NULL;
02375 }
02376
02377
02378 if( (fabs(circle->radius()) > m_param.MIN_RADIUS_OF_STRANGE_TRACK) || (fabs(circle->radius()) < m_param.ULTIMATE_MIN_RADIUS_OF_STRANGE_TRACK) ){
02379 #if DEBUG_CURL_DUMP
02380 std::cout << "(TCurlFinder) 3D:Fail...success 3D fit, but large radius > "
02381 << m_param.MIN_RADIUS_OF_STRANGE_TRACK << "." << std::endl;
02382 #endif
02383 return NULL;
02384 }
02385
02386 if(track->links().length() >= m_param.minimum_3DTrackLength){
02387 #if DEBUG_CURL_DUMP
02388 std::cout << "(TCurlFinder) 3D:Success...can build stereo!!!" << std::endl;
02389 #endif
02390 return track;
02391 }else{
02392 #if DEBUG_CURL_DUMP
02393 std::cout << "(TCurlFinder) 3D:Fail...success 3D fit, but final hit wire # < 3." << std::endl;
02394 #endif
02395 return NULL;
02396 }
02397 }
02398
02399 TMLink *
02400 findIsolatedCloseHits(TMLink *link)
02401 {
02402 int nNeighbor = 0;
02403 int nIsolated = 0;
02404 TMLink * isolatedLink[2] = {NULL,NULL};
02405 unsigned layerID = link->wire()->layerId();
02406 unsigned localID = link->wire()->localId();
02407 for(int i=0;i<6;++i){
02408 if(link->neighbor(i)){
02409 if(link->neighbor(i)->wire()->layerId() == layerID){
02410 ++nNeighbor;
02411 int isolated = 1;
02412 int testEach = 0;
02413 for(int j=0;j<6;++j){
02414 if(link->neighbor(i)->neighbor(j)){
02415 if(link->neighbor(i)->neighbor(j)->wire()->layerId() == layerID &&
02416 link->neighbor(i)->neighbor(j)->wire()->localId() != localID){
02417 isolated = 0;
02418 }
02419 #if 1
02420 else if(link->neighbor(i)->neighbor(j)->wire()->layerId() == layerID &&
02421 link->neighbor(i)->neighbor(j)->wire()->localId() == localID){
02422 testEach = 1;
02423 }
02424 #endif
02425 }else break;
02426 }
02427 if(isolated == 1){
02428 if(nIsolated < 2)
02429 isolatedLink[nIsolated] = link->neighbor(i);
02430 ++nIsolated;
02431 }
02432 #if 1
02433 if(testEach == 0){
02434 std::cout << "Why?? Neighborhood info. dose not exist!!" << std::endl;
02435 }
02436 #endif
02437 }
02438 }else break;
02439 }
02440 #if 0
02441 std::cout << "isolated/neighbor # = " << nIsolated << "/" << nNeighbor << std::endl;
02442 std::cout << "layer ID = " << layerID << " ";
02443 std::cout << "local ID = " << localID << " --> ";
02444 if(isolatedLink[0])std::cout << isolatedLink[0]->wire()->localId() << " ";
02445 if(isolatedLink[1])std::cout << isolatedLink[1]->wire()->localId() << " ";
02446 std::cout << std::endl;
02447 #endif
02448 if(nIsolated == 1 &&
02449 nNeighbor == 1 && isolatedLink[0])return isolatedLink[0];
02450 else return NULL;
02451 }
02452
02453 void
02454 TCurlFinder::findCloseHits(AList<TMLink> &links,
02455 TTrack &track, AList<TMLink> &list) {
02456
02457
02458
02459 double dRcut[11] = {m_param.RANGE_FOR_STEREO_FIRST, m_param.RANGE_FOR_STEREO_SECOND,
02460 0., 0.,
02461 0., m_param.RANGE_FOR_STEREO_THIRD,
02462 m_param.RANGE_FOR_STEREO_FORTH, m_param.RANGE_FOR_STEREO_FIFTH,
02463 m_param.RANGE_FOR_STEREO_SIXTH, 0.,
02464 0.};
02465 double r = fabs(track.helix().curv());
02466 double q = track.charge();
02467 double x = track.helix().center().x();
02468 double y = track.helix().center().y();
02469 for(unsigned i = 0, size = links.length(); i < size; ++i){
02470 if(fabs((links[i]->wire()->xyPosition() - track.helix().center()).mag() - r) <
02471 dRcut[links[i]->wire()->superLayerId()]){
02472 if(q*(x*links[i]->wire()->xyPosition().y()-y*links[i]->wire()->xyPosition().x()) > 0.){
02473 list.remove(links[i]);
02474 list.append(links[i]);
02475 TMLink *cand = findIsolatedCloseHits(links[i]);
02476 if(cand){
02477 list.remove(cand);
02478 list.append(cand);
02479 }
02480 }
02481 }
02482 }
02483 return;
02484 }
02485
02486
02487
02488
02489
02490
02491
02492 int
02493 TCurlFinder::makeWithMC(const AList<TMDCWireHit> & axialHits,
02494 const AList<TMDCWireHit> & stereoHits,
02495 AList<TTrack> & tracks) {
02496 #if DEBUG_CURL_MC
02497 #define MAX_INDEX_MAKEMC 100
02498 #if DEBUG_CURL_DUMP
02499 std::cout << "(TCurlFinder)Now making tracks using MC info..." << std::endl;
02500 #endif
02501 int index[MAX_INDEX_MAKEMC];
02502 for(unsigned i = 0; i < MAX_INDEX_MAKEMC; ++i)index[i] = 9999;
02503
02504 int counter(0);
02505 bool first(true);
02506
02507 for(unsigned i = 0, size = axialHits.length(); i < size; ++i){
02508 if(axialHits[i]->mc() &&
02509 axialHits[i]->mc()->hep() &&
02510 !(axialHits[i]->state() & WireHitUsed)){
02511 int flag(1);
02512 for(unsigned j = 0; j < MAX_INDEX_MAKEMC; ++j){
02513 if(index[j] != 9999 && index[j] == axialHits[i]->mc()->hep()->id()){
02514 flag = 0;
02515 break;
02516 }
02517 }
02518 if(flag){
02519 index[counter] = axialHits[i]->mc()->hep()->id();
02520 ++counter;
02521 }
02522 }
02523 }
02524 #if DEBUG_CURL_DUMP
02525 std::cout << "(TCurlFinder)Found " << counter
02526 << " tracks with MC information." << std::endl;
02527 #endif
02528 for(unsigned j = 0; j < counter; ++j){
02529 AList<TMLink> axialList;
02530 AList<TMLink> stereoList;
02531 int axialCounter(0);
02532 int stereoCounter(0);
02533
02534 for(unsigned i = 0, size = axialHits.length(); i < size; ++i){
02535 if(index[j] == axialHits[i]->mc()->hep()->id() &&
02536 !(axialHits[i]->state() & WireHitUsed)){
02537 axialList.append(new TMLink(0, axialHits[i]));
02538 ++axialCounter;
02539 }
02540 }
02541 if(axialCounter < 3){
02542 HepAListDeleteAll(axialList);
02543 continue;
02544 }
02545
02546 for(unsigned i = 0, size = stereoHits.length(); i < size; ++i){
02547 if(index[j] == stereoHits[i]->mc()->hep()->id() &&
02548 !(stereoHits[i]->state() & WireHitUsed)){
02549 stereoList.append(new TMLink(0, stereoHits[i]));
02550 ++stereoCounter;
02551 }
02552 }
02553 if(stereoCounter < 2){
02554 HepAListDeleteAll(axialList);
02555 HepAListDeleteAll(stereoList);
02556 continue;
02557 }
02558 #if DEBUG_CURL_DUMP
02559 std::cout << "(TCurlFinder)#" << j << " : Use "
02560 << axialCounter << " axial hit wires and "
02561 << stereoCounter << " stereo hit wires" << std::endl;
02562 std::cout << "(TCurlFinder)Particle Type(LUND) = "
02563 << axialList[0]->hit()->mc()->hep()->pType() << std::endl;
02564 #endif
02565
02566 m_unusedAxialHitsOriginal.append(axialList);
02567 m_unusedAxialHitsOriginal.append(stereoList);
02568 TCircle *circle = new TCircle(axialList);
02569 m_allCircles.append(circle);
02570 circle->fitForCurl();
02571 double charge = 1.;
02572 if(axialList[0]->hit()->mc()->hep()->pType() < 0)charge = -1.;
02573 if(fabs(axialList[0]->hit()->mc()->hep()->pType()) == 11 ||
02574 fabs(axialList[0]->hit()->mc()->hep()->pType()) == 13 ||
02575 fabs(axialList[0]->hit()->mc()->hep()->pType()) == 15)charge *= -1.;
02576 circle->property(charge, charge*fabs(circle->radius()), circle->center());
02577
02578 AList<TMLink> removeList;
02579 double x = circle->center().x();
02580 double y = circle->center().y();
02581
02582 for(unsigned i = 0, size = axialList.length();
02583 i < size; ++i){
02584 if(charge*(x*axialList[i]->xyPosition().y()-
02585 y*axialList[i]->xyPosition().x())< 0.){
02586 removeList.append(axialList[i]);
02587
02588 }
02589 }
02590 circle->remove(removeList);
02591 if(circle->nLinks() < 3)continue;
02592
02593 circle->fitForCurl(1);
02594 x = circle->center().x();
02595 y = circle->center().y();
02596 removeList.removeAll();
02598 for(unsigned i = 0, size = stereoList.length();
02599 i < size; ++i){
02600 if(charge*(x*stereoList[i]->xyPosition().y()-
02601 y*stereoList[i]->xyPosition().x())< 0.){
02602 removeList.append(stereoList[i]);
02603 }
02604 }
02605 stereoList.remove(removeList);
02606 if(stereoList.length() < 2)continue;
02607
02608 TTrack *track = new TTrack(*circle);
02609 m_allTracks.append(track);
02610 if(m_builder.buildStereoMC(*track, stereoList)){
02611 if(track->links().length() >= 5){
02612 track->assign(WireHitCurlFinder, TrackCurlFinder | TrackValid | Track3D);
02613 tracks.append(track);
02614 m_allTracks.remove(track);
02615 }else{
02616 std::cout << "Can not reconstruct with MC information!" << std::endl;
02617 }
02618 }else{
02619 std::cout << "Can not reconstruct with MC information!" << std::endl;
02620 }
02621 }
02622 #endif
02623 return 0;
02624 }
02625
02626 void
02627 TCurlFinder::makeCdcFrame(void) {
02628 #if DEBUG_CURL_GNUPLOT+DEBUG_CURL_SEGMENT
02629
02630 double X = 0.;
02631 double Y = 0.;
02632 double R[12] = {8.3, 16.9, 21.7, 31.3, 36.1, 44.1,
02633 50.5, 58.5, 64.9, 72.9, 79.3, 87.4};
02634 double step = 300.;
02635 double dStep = 2.*M_PI/step;
02636 FILE *data;
02637 std::string nameHead = "tmp.cdc_";
02638 for(int j=0;j<12;++j){
02639 std::string nameFile = nameHead+"0"+itostring(j);
02640 if(j>=10)nameFile = nameHead+itostring(j);
02641 if((data = fopen(nameFile,"w")) != NULL){
02642 for(int i=0;i<step;++i){
02643 double x = X + R[j] * cos(dStep*static_cast<double>(i));
02644 double y = Y + R[j] * sin(dStep*static_cast<double>(i));
02645 fprintf(data,"%lf, %lf\n",x,y);
02646 }
02647 fclose(data);
02648 }
02649 }
02650
02651 if((data = fopen("tmp_wires.dat","w")) != NULL){
02652 AList<TMLink> list = m_unusedAxialHitsOriginal;
02653 list.append(m_unusedStereoHitsOriginal);
02654 for(int i=0;i<list.length();i++){
02655 double x = list[i]->hit()->wire()->xyPosition().x();
02656 double y = list[i]->hit()->wire()->xyPosition().y();
02657 fprintf(data,"%lf, %lf\n",x,y);
02658 }
02659 fclose(data);
02660 }
02661 #endif
02662 return;
02663 }
02664
02665 void
02666 TCurlFinder::plotSegment(const AList<TMLink>& list, const int flag) {
02667 #if DEBUG_CURL_GNUPLOT
02668 if(!m_debugCdcFrame){
02669 makeCdcFrame();
02670 m_debugCdcFrame = true;
02671 }
02672 double gmaxX = 90. ,gminX = -90.;
02673 double gmaxY = 90. ,gminY = -90.;
02674 FILE *gnuplot, *data;
02675 if((data = fopen("tmp.dat","w")) != NULL){
02676 if(flag)std::cout << "Wire ID = ";
02677 for(int i=0;i<list.length();i++){
02678 double x = list[i]->hit()->wire()->xyPosition().x();
02679 double y = list[i]->hit()->wire()->xyPosition().y();
02680 fprintf(data,"%lf, %lf\n",x,y);
02681 if(flag)std::cout << list[i]->hit()->wire()->id() << ", ";
02682 }
02683 if(flag)std::cout << std::endl;
02684 fclose(data);
02685 }
02686 if((gnuplot = popen("gnuplot","w")) != NULL){
02687 fprintf(gnuplot,"set nokey \n");
02688 fprintf(gnuplot,"set size 0.721,1.0 \n");
02689 fprintf(gnuplot,"set xrange [%f:%f] \n",gminX,gmaxX);
02690 fprintf(gnuplot,"set yrange [%f:%f] \n",gminY,gmaxY);
02691 std::string longName = "plot \"tmp_wires.dat\", \"tmp.dat\"";
02692 std::string nameHead = ",\"tmp.cdc_";
02693 for(int j=0;j<12;++j){
02694 std::string nameFile = nameHead+"0"+itostring(j)+"\"w l 0";
02695 if(j>=10)nameFile = nameHead+itostring(j)+"\"w l 0";
02696 longName += nameFile;
02697 }
02698 longName += " \n";
02699 fprintf(gnuplot,longName);
02700 fflush(gnuplot);
02701 char tmp[8];
02702 gets(tmp);
02703 pclose(gnuplot);
02704 }
02705 #endif
02706 return;
02707 }
02708
02709 void
02710 TCurlFinder::plotCircle(const TCircle& circle, const int flag) {
02711 #if DEBUG_CURL_GNUPLOT
02712
02713 if(!m_debugCdcFrame){
02714 makeCdcFrame();
02715 m_debugCdcFrame = true;
02716 }
02717 double gmaxX = 90. ,gminX = -90.;
02718 double gmaxY = 90. ,gminY = -90.;
02719 FILE *gnuplot, *data;
02720 if((data = fopen("tmp.dat1","w")) != NULL){
02721 if(flag)std::cout << "Axial Wire ID ==> " << std::endl;
02722 for(int i=0;i<circle.nLinks();++i){
02723 if(circle.links()[i]->hit()->wire()->axial()){
02724 double x = circle.links()[i]->hit()->wire()->xyPosition().x();
02725 double y = circle.links()[i]->hit()->wire()->xyPosition().y();
02726 fprintf(data,"%lf, %lf\n",x,y);
02727 if(flag){
02728
02729
02730
02731
02732
02733 }
02734 }
02735 }
02736 if(flag)std::cout << std::endl;
02737 fclose(data);
02738 }
02739 if((data = fopen("tmp.dat2","w")) != NULL){
02740 if(flag)std::cout << "Stereo Wire ID ==> " << std::endl;
02741 for(int i=0;i<circle.nLinks();++i){
02742 if(circle.links()[i]->hit()->wire()->stereo()){
02743 double x = circle.links()[i]->hit()->wire()->xyPosition().x();
02744 double y = circle.links()[i]->hit()->wire()->xyPosition().y();
02745 fprintf(data,"%lf, %lf\n",x,y);
02746 if(flag){
02747
02748
02749
02750
02751
02752 }
02753 }
02754 }
02755 if(flag)std::cout << std::endl;
02756 fclose(data);
02757 }
02758 double X = circle.center().x();
02759 double Y = circle.center().y();
02760 double R = fabs(circle.radius());
02761 double step = 300.;
02762 double dStep = 2.*M_PI/step;
02763 if((data = fopen("tmp.dat3","w")) != NULL){
02764 for(int i=0;i<step;++i){
02765 double x = X + R * cos(dStep*static_cast<double>(i));
02766 double y = Y + R * sin(dStep*static_cast<double>(i));
02767 fprintf(data,"%lf, %lf\n",x,y);
02768 }
02769 fclose(data);
02770 }
02771 if((gnuplot = popen("gnuplot","w")) != NULL){
02772 fprintf(gnuplot,"set nokey \n");
02773 fprintf(gnuplot,"set size 0.721,1.0 \n");
02774 fprintf(gnuplot,"set xrange [%f:%f] \n",gminX,gmaxX);
02775 fprintf(gnuplot,"set yrange [%f:%f] \n",gminY,gmaxY);
02776 std::string longName = "plot \"tmp_wires.dat\", \"tmp.dat1\", \"tmp.dat3\" w l, \"tmp.dat2\"";
02777 std::string nameHead = ",\"tmp.cdc_";
02778 for(int j=0;j<12;++j){
02779 std::string nameFile = nameHead+"0"+std::string(j)+"\"w l 0";
02780 if(j>=10)nameFile = nameHead+std::string(j)+"\"w l 0";
02781 longName += nameFile;
02782 }
02783 longName += " \n";
02784 fprintf(gnuplot,longName);
02785 fflush(gnuplot);
02786 char tmp[8];
02787 gets(tmp);
02788 pclose(gnuplot);
02789 }
02790 #endif
02791 return;
02792 }
02793
02794 void
02795 TCurlFinder::plotTrack(const TTrack& track, const int flag) {
02796 #if DEBUG_CURL_GNUPLOT
02797 if(!m_debugCdcFrame){
02798 makeCdcFrame();
02799 m_debugCdcFrame = true;
02800 }
02801 double gmaxX = 90. ,gminX = -90.;
02802 double gmaxY = 90. ,gminY = -90.;
02803 FILE *gnuplot, *data;
02804 if((data = fopen("tmp.dat1","w")) != NULL){
02805 if(flag)std::cout << "Axial Wire ID ==> " << std::endl;
02806 for(int i=0;i<track.nLinks();++i){
02807 if(track.links()[i]->hit()->wire()->axial()){
02808 double x = track.links()[i]->hit()->wire()->xyPosition().x();
02809 double y = track.links()[i]->hit()->wire()->xyPosition().y();
02810 fprintf(data,"%lf, %lf\n",x,y);
02811 if(flag){
02812 if(debugMcFlag){
02813 std::cout << " A:" << track.links()[i]->hit()->wire()->id() << ", ";
02814 std::cout << ", HepTrackID = " << track.links()[i]->hit()->mc()->hep()->id();
02815 std::cout << ", HepLundID = " << track.links()[i]->hit()->mc()->hep()->pType() << std::endl;
02816 }else std::cout << " A:" << track.links()[i]->hit()->wire()->id() << std::endl;
02817 }
02818 }
02819 }
02820 if(flag)std::cout << std::endl;
02821 fclose(data);
02822 }
02823 if((data = fopen("tmp.dat2","w")) != NULL){
02824 if(flag)std::cout << "Stereo Wire ID ==> " << std::endl;
02825 for(int i=0;i<track.nLinks();++i){
02826 if(track.links()[i]->hit()->wire()->stereo()){
02827 double x = track.links()[i]->hit()->wire()->xyPosition().x();
02828 double y = track.links()[i]->hit()->wire()->xyPosition().y();
02829 fprintf(data,"%lf, %lf\n",x,y);
02830 if(flag){
02831 if(debugMcFlag){
02832 std::cout << " S:" << track.links()[i]->hit()->wire()->id() << ", ";
02833 std::cout << ", HepTrackID = " << track.links()[i]->hit()->mc()->hep()->id();
02834 std::cout << ", HepLundID = " << track.links()[i]->hit()->mc()->hep()->pType() << std::endl;
02835 }else std::cout << " S:" << track.links()[i]->hit()->wire()->id() << std::endl;
02836 }
02837 }
02838 }
02839 if(flag)std::cout << std::endl;
02840 fclose(data);
02841 }
02842 double X = track.helix().center().x();
02843 double Y = track.helix().center().y();
02844 double R = fabs(track.helix().radius());
02845 double step = 300.;
02846 double dStep = 2.*M_PI/step;
02847 if((data = fopen("tmp.dat3","w")) != NULL){
02848 for(int i=0;i<step;++i){
02849 double x = X + R * cos(dStep*static_cast<double>(i));
02850 double y = Y + R * sin(dStep*static_cast<double>(i));
02851 fprintf(data,"%lf, %lf\n",x,y);
02852 }
02853 fclose(data);
02854 }
02855 if((gnuplot = popen("gnuplot","w")) != NULL){
02856 fprintf(gnuplot,"set nokey \n");
02857 fprintf(gnuplot,"set size 0.721,1.0 \n");
02858 fprintf(gnuplot,"set xrange [%f:%f] \n",gminX,gmaxX);
02859 fprintf(gnuplot,"set yrange [%f:%f] \n",gminY,gmaxY);
02860 std::string longName = "plot \"tmp_wires.dat\", \"tmp.dat1\", \"tmp.dat3\" w l, \"tmp.dat2\"";
02861 std::string nameHead = ",\"tmp.cdc_";
02862 for(int j=0;j<12;++j){
02863 std::string nameFile = nameHead+"0"+itostring(j)+"\"w l 0";
02864 if(j>=10)nameFile = nameHead+itostring(j)+"\"w l 0";
02865 longName += nameFile;
02866 }
02867 longName += " \n";
02868 fprintf(gnuplot,longName);
02869 fflush(gnuplot);
02870 char tmp[8];
02871 gets(tmp);
02872 pclose(gnuplot);
02873 }
02874 #endif
02875 return;
02876 }
02877
02878 void
02879 TCurlFinder::writeSegment(const AList<TMLink>& list, const int type) {
02880 #if DEBUG_CURL_SEGMENT
02881 if(!m_debugCdcFrame){
02882 makeCdcFrame();
02883 m_debugCdcFrame = true;
02884 }
02885 double gmaxX = 90. ,gminX = -90.;
02886 double gmaxY = 90. ,gminY = -90.;
02887
02888 FILE *data;
02889 std::string nameHead = "tmp.segment_";
02890 std::string nameFile = nameHead+itostring(type)+"_"+itostring(m_debugFileNumber);
02891 ++m_debugFileNumber;
02892 if((data = fopen(nameFile,"w")) != NULL){
02893 for(int i=0;i<list.length();i++){
02894 double x = list[i]->hit()->wire()->xyPosition().x();
02895 double y = list[i]->hit()->wire()->xyPosition().y();
02896 fprintf(data,"%lf, %lf\n",x,y);
02897 }
02898 fclose(data);
02899 }
02900 #endif
02901 return;
02902 }
02903
02904 void
02905 TCurlFinder::dumpType1(TTrack *track) {
02906 #if DEBUG_CURL_DUMP
02907 for(int j=0;j<track->nLinks();++j){
02908 std::cout << "Used Wire Info...";
02909 if(track->links()[j]->hit()->wire()->axial()){
02910 std::cout << "A:" << track->links()[j]->hit()->wire()->id() << ", ";
02911 }else{
02912 std::cout << "S:" << track->links()[j]->hit()->wire()->id() << ", ";
02913 }
02914 if(debugMcFlag){
02915 std::cout << ", HepTrackID = " << track->links()[j]->hit()->mc()->hep()->id();
02916 std::cout << ", HepLundID = " << track->links()[j]->hit()->mc()->hep()->pType();
02917 }
02918 double dist = distance(*track, *(track->links()[j]));
02919 if(dist > 2.)std::cout << ": Large Distance( >2cm ) = " << dist;
02920 std::cout << std::endl;
02921 }
02922 AList<TMLink> list=m_unusedAxialHits;
02923 list.append(m_unusedStereoHits);
02924 for(unsigned j=0, nList=list.length();j<nList;++j){
02925 double dist = distance(*track, *(list[j]));
02926 std::cout << "Close Wire Info in ALL( <0.5cm )...";
02927 if(dist < 0.5){
02928 if(list[j]->hit()->wire()->axial())
02929 std::cout << "CA:" << list[j]->hit()->wire()->id() << ", ";
02930 else
02931 std::cout << "CS:" << list[j]->hit()->wire()->id() << ", ";
02932 if(debugMcFlag){
02933 std::cout << ", HepTrackID = " << list[j]->hit()->mc()->hep()->id();
02934 std::cout << ", HepLundID = " << list[j]->hit()->mc()->hep()->pType();
02935 }
02936 std::cout << ", Distance = " << dist << std::endl;
02937 }
02938 }
02939 #endif
02940 return;
02941 }
02942
02943 void
02944 TCurlFinder::dumpType2(TTrack *track) {
02945 #if DEBUG_CURL_DUMP
02946 unsigned size = track->nLinks();
02947 if(size == 0)return;
02948
02949 set< int, less<int> > uniqueHepID;
02950 vector<int> hepID;
02951 vector<double> ratio;
02952 for(int i=0;i<size;++i){
02953 uniqueHepID.insert(track->links()[i]->hit()->mc()->hep()->id());
02954 hepID.push_back(track->links()[i]->hit()->mc()->hep()->id());
02955
02956 }
02957
02958 set< int, less<int> >::iterator u = uniqueHepID.begin();
02959 vector<int>::size_type sizeInt;
02960 for(unsigned i=0;i<uniqueHepID.size();++i){
02961 sizeInt = 0;
02962 count(hepID.begin(), hepID.end(), *u, sizeInt);
02963 ratio.push_back((static_cast<double>(sizeInt)/static_cast<double>(size)));
02964
02965 ++u;
02966 }
02967
02968 vector<double>::iterator m = max_element(ratio.begin(), ratio.end());
02969 int maxIndex = 0;
02970 ::distance(ratio.begin(), m, maxIndex);
02971 u = uniqueHepID.begin();
02972 advance(u,maxIndex);
02973
02974 std::cout << "Ratio " << ratio[maxIndex] << std::endl;
02975 for(int i=0;i<size;++i){
02976 if(track->links()[i]->hit()->wire()->axial())std::cout << "A ";
02977 else std::cout << "S ";
02978
02979 double dist = distance(*track, *(track->links()[i]));
02980 if(*u != track->links()[i]->hit()->mc()->hep()->id()){
02981 std::cout << "Bad " << dist << std::endl;
02982 }else{
02983 std::cout << "Good " << dist << std::endl;
02984 }
02985 }
02986 #endif
02987 return;
02988 }
02989
02990
02991 #if 0
02992
02993
02994
02995
02996
02997
02998
02999
03000
03001
03002
03003
03004
03005 Gen_hepevt *
03006 cluster2hep(Recsvd_cluster *clus) {
03007 AList<Datsvd_hit> m_datsvd_hit;
03008 Datsvd_hit_Manager& svdHitMgr = Datsvd_hit_Manager::get_manager();
03009 for(Datsvd_hit_Manager::iterator
03010 it = svdHitMgr.begin(),
03011 end = svdHitMgr.end();
03012 it != end; ++it){
03013 m_datsvd_hit.append(it);
03014 }
03015
03016
03017 unsigned size = m_datsvd_hit.length();
03018 int startPosition = -1;
03019 int direction = 1;
03020 #if 0
03021 for(unsigned i=0;i<size;++i){
03022
03023 std::cout << i << ": " << clus->hit().get_ID() << " <--> " << m_datsvd_hit[i]->get_ID()
03024 << ", RLA = " << m_datsvd_hit[i]->rla() << ", LSA = " << m_datsvd_hit[i]->amp()
03025 << ", Width = " << clus->width()
03026 << ", Cluster LSA = " << clus->lsa() << std::endl;
03027 }
03028 #endif
03029 for(unsigned i=0;i<size;++i){
03030 if(m_datsvd_hit[i]->amp() == 0)std::cout << "DatSVD_Hit:amp == 0" << std::endl;
03031
03032 if(m_datsvd_hit[i]->rla() == 0)std::cout << "DatSVD_Hit:rla == 0" << std::endl;
03033 if(m_datsvd_hit[i]->rla() == 81920)std::cout << "DatSVD_Hit:rla == 81920" << std::endl;
03034 if(clus->hit().get_ID() == m_datsvd_hit[i]->get_ID()){
03035 startPosition = i;
03036 if(static_cast<double>(m_datsvd_hit[i]->amp()-1) > clus->lsa())direction = -1;
03037 break;
03038 }
03039 }
03040 if(startPosition == -1)return NULL;
03041 int width = clus->width();
03042 #if 0
03043 std::cout << "Start # = " << startPosition
03044 << ", Width = " << width
03045 << ", Direction = " << direction <<std::endl;
03046 #endif
03047
03048 int* hepID = new int[width];
03049 set< int, less<int> > uniqueHepID;
03050
03051 for(int i=startPosition;i<startPosition+width;++i)hepID[i-startPosition] = 0;
03052
03053 Datsvd_mcdata_Manager& svdMcHitMgr = Datsvd_mcdata_Manager::get_manager();
03054
03055 if(direction == 1){
03056 for(int i=startPosition;i<startPosition+width;++i){
03057 for(Datsvd_mcdata_Manager::iterator
03058 it = svdMcHitMgr.begin(),
03059 end = svdMcHitMgr.end();
03060 it != end; ++it){
03061 if(it->Hep()){
03062 if(m_datsvd_hit[i]->rla() == it->rla() &&
03063 it->Hep().get_ID() != 0){
03064 hepID[i-startPosition] = it->Hep().get_ID();
03065 uniqueHepID.insert(it->Hep().get_ID());
03066 break;
03067 }
03068 }
03069 }
03070 }
03071 }else{
03072 int reverse = 0;
03073 for(int i=startPosition;i<startPosition+width;++i){
03074 ++reverse;
03075
03076
03077 for(Datsvd_mcdata_Manager::iterator
03078 it = svdMcHitMgr.begin(),
03079 end = svdMcHitMgr.end();
03080 it != end; ++it){
03081 if(it->Hep()){
03082 if(m_datsvd_hit[startPosition+1-reverse]->rla() == it->rla() &&
03083 it->Hep().get_ID() != 0){
03084 hepID[i-startPosition] = it->Hep().get_ID();
03085 uniqueHepID.insert(it->Hep().get_ID());
03086 break;
03087 }
03088 }
03089 }
03090 }
03091 }
03092
03093 unsigned num = uniqueHepID.size();
03094 int* counter = new int[num];
03095 set< int, less<int> >::iterator u = uniqueHepID.begin();
03096 for(int i=0;i<num;++i) counter[i] = 0;
03097
03098 for(int i=0;i<num;++i){
03099 for(int j=0;j<width;++j){
03100 if(*u == hepID[j]){
03101 counter[i] += 1;
03102 }
03103 }
03104 ++u;
03105 }
03106
03107 Gen_hepevt_Manager& genMgr = Gen_hepevt_Manager::get_manager();
03108 #if 0
03109 u = uniqueHepID.begin();
03110 for(int i=0;i<num;++i){
03111 std::cout << i << ": TrackID = "<< *u - 1
03112 << ", Count = " << counter[i]
03113 << ", LundID = " << genMgr[*u-1].idhep() << std::endl;
03114 ++u;
03115 }
03116 #endif
03117
03118 delete [] hepID;
03119 delete [] counter;
03120 if(num == 1){
03121
03122 return &(genMgr[*(uniqueHepID.begin())-1]);
03123 }else if(num >= 2){
03124
03125 return NULL;
03126 }else{
03127 return NULL;
03128 }
03129 }
03130 #endif
03131
03132 void
03133 TCurlFinder::debugCheckSegments1(void) {
03134 #if DEBUG_CURL_SEGMENT
03135
03136
03137 std::cout << "(TCurlFinder)checking consistency of segement..." << std::endl;
03138 unsigned nA = m_unusedAxialHitsOriginal.length();
03139 if(nA >= 2){
03140 for(unsigned i=0;i<nA-1;++i){
03141 int superLayerId = (int)(m_unusedAxialHitsOriginal[i]->wire()->superLayerId());
03142 int layerId = (int)(m_unusedAxialHitsOriginal[i]->wire()->layerId());
03143 int localId = (int)(m_unusedAxialHitsOriginal[i]->wire()->localId());
03144 int localIdP = (int)(m_unusedAxialHitsOriginal[i]->wire()->localIdForPlus());
03145 int localIdM = (int)(m_unusedAxialHitsOriginal[i]->wire()->localIdForMinus());
03146 for(unsigned j=i+1;j<nA;++j){
03147 int superLayerId2 = (int)(m_unusedAxialHitsOriginal[j]->wire()->superLayerId());
03148 int layerId2 = (int)(m_unusedAxialHitsOriginal[j]->wire()->layerId());
03149 int localId2 = (int)(m_unusedAxialHitsOriginal[j]->wire()->localId());
03150 if(superLayerId == superLayerId2){
03151 if(layerId2 == layerId){
03152 if(localIdP+1 == localId2 || localIdM-1 == localId2)
03153 debugCheckSegments(localId, layerId,
03154 localId2,layerId2);
03155 }else if(layerId2 == layerId-1 || layerId2 == layerId+1){
03156 if(offset(layerId) == offset(layerId2)){
03157 std::cout << "(TCurlFinder: Waring) Offset is same at the same superlayer!!" << std::endl;
03158 }else if(offset(layerId) > offset(layerId2)){
03159 if(localId == localId2 || localIdP+1 == localId2)
03160 debugCheckSegments(localId, layerId,
03161 localId2,layerId2);
03162 }else{
03163 if(localId == localId2 || localIdM-1 == localId2)
03164 debugCheckSegments(localId, layerId,
03165 localId2,layerId2);
03166 }
03167 }
03168 }
03169 }
03170 }
03171 }
03172 unsigned nS = m_unusedStereoHitsOriginal.length();
03173 if(nS >= 2){
03174 for(unsigned i=0;i<nS-1;++i){
03175 int superLayerId = (int)(m_unusedStereoHitsOriginal[i]->wire()->superLayerId());
03176 int layerId = (int)(m_unusedStereoHitsOriginal[i]->wire()->layerId());
03177 int localId = (int)(m_unusedStereoHitsOriginal[i]->wire()->localId());
03178 int localIdP = (int)(m_unusedStereoHitsOriginal[i]->wire()->localIdForPlus());
03179 int localIdM = (int)(m_unusedStereoHitsOriginal[i]->wire()->localIdForMinus());
03180 for(unsigned j=i+1;j<nS;++j){
03181 int superLayerId2 = (int)(m_unusedStereoHitsOriginal[j]->wire()->superLayerId());
03182 int layerId2 = (int)(m_unusedStereoHitsOriginal[j]->wire()->layerId());
03183 int localId2 = (int)(m_unusedStereoHitsOriginal[j]->wire()->localId());
03184 if(superLayerId == superLayerId2){
03185 if(layerId2 == layerId){
03186 if(localIdP+1 == localId2 || localIdM-1 == localId2)
03187 debugCheckSegments(localId, layerId,
03188 localId2,layerId2);
03189 }else if(layerId2 == layerId-1 || layerId2 == layerId+1){
03190 if(offset(layerId) == offset(layerId2)){
03191 std::cout << "(TCurlFinder: Waring) Offset is same at the same superlayer!!" << std::endl;
03192 }else if(offset(layerId) > offset(layerId2)){
03193 if(localId == localId2 || localIdP+1 == localId2)
03194 debugCheckSegments(localId, layerId,
03195 localId2,layerId2);
03196 }else{
03197 if(localId == localId2 || localIdM-1 == localId2)
03198 debugCheckSegments(localId, layerId,
03199 localId2,layerId2);
03200 }
03201 }
03202 }
03203 }
03204 }
03205 }
03206 std::cout << "(TCurlFinder)...done check of segement!" << std::endl;
03207 std::cout << "(TCurlFinder)...If no warning message exists, check of segement is complete!" << std::endl;
03208 std::cout << "(TCurlFinder)...Note: a segment size should be 1 or 2 to use this debugger." << std::endl;
03209 #endif
03210 return;
03211 }
03212
03213 void
03214 TCurlFinder::debugCheckSegments(const double localId, const double layerId,
03215 const double localId2,const double layerId2) {
03216 #if DEBUG_CURL_DUMP
03217 unsigned nSeg = m_segmentList.length();
03218 unsigned nFound = 0;
03219 for(unsigned i=0;i<nSeg;++i){
03220 unsigned nWire = m_segmentList[i]->list().length();
03221 unsigned mFound = 0;
03222 for(unsigned j=0;j<nWire;++j){
03223 if(((m_segmentList[i]->list())[j])->wire()->layerId() == layerId &&
03224 ((m_segmentList[i]->list())[j])->wire()->localId() == localId)++mFound;
03225 if(((m_segmentList[i]->list())[j])->wire()->layerId() == layerId2 &&
03226 ((m_segmentList[i]->list())[j])->wire()->localId() == localId2)++mFound;
03227 }
03228 if(mFound != 0 && mFound != 2){
03229 std::cout << "(TCurlFinder: Warning) Segment is inconsistency(0)!! mFound = " << mFound << std::endl;
03230 }
03231 if(mFound == 2)++nFound;
03232 }
03233 if(nFound != 1)
03234 std::cout << "(TCurlFinder: Warning) Segment is inconsistency(1)!! nFound = " << nFound << std::endl;
03235 #endif
03236 return;
03237 }
03238
03239 void
03240 TCurlFinder::debugCheckSegments0(void) {
03241 #if DEBUG_CURL_SEGMENT
03242 unsigned nSeg = m_segmentList.length();
03243 unsigned nWire = 0;
03244 for(unsigned i=0;i<nSeg;++i)nWire += m_segmentList[i]->list().length();
03245
03246 unsigned nWireOriginal = m_unusedAxialHitsOriginal.length()+
03247 m_unusedStereoHitsOriginal.length();
03248
03249 std::cout << "(TCurlFinder: SelfChecker) Segment Parts" << std::endl;
03250 std::cout << " MIN_SEGMENT = " << m_param.MIN_SEGMENT << std::endl;
03251 std::cout << " Wire # of Orinal List = " << nWireOriginal << std::endl;
03252 std::cout << " Wire # of Segments = " << nWire << std::endl;
03253 std::cout << " If MIN_SEGMENT <= 1, above numbers should be same." << std::endl;
03254 std::cout << " If MIN_SEGMENT > 1, former >= latter." << std::endl;
03255 #endif
03256 return;
03257 }
03258
03259 void
03260 TCurlFinder::debugCheckSegments2(void) {
03261 #if DEBUG_CURL_SEGMENT
03262
03263 #define DEBUG_TMP_N_CURL 50
03264
03265 unsigned nSeg = m_segmentList.length();
03266 unsigned nWire[DEBUG_TMP_N_CURL] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
03267 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
03268 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
03269 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
03270 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
03271 for(unsigned i=0;i<nSeg;++i){
03272 if(m_segmentList[i]->list().length() < DEBUG_TMP_N_CURL)
03273 ++(nWire[m_segmentList[i]->list().length()]);
03274 else
03275 ++(nWire[DEBUG_TMP_N_CURL-1]);
03276 }
03277 std::ifstream fin("tmp.wire.data");
03278 unsigned nTotalWire[DEBUG_TMP_N_CURL] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
03279 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
03280 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
03281 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
03282 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
03283 if(fin){
03284 for(int i=0;i<DEBUG_TMP_N_CURL;++i)fin >> nTotalWire[i];
03285 }else{
03286 std::cout << "(TCurlFinder) tmp.wire.data does not exist!" << std::endl;
03287 }
03288 for(int i=0;i<DEBUG_TMP_N_CURL;++i)nTotalWire[i] += nWire[i];
03289 std::ofstream fout("tmp.wire.data");
03290 if(fout){
03291 fout << nTotalWire[0];
03292 for(int i=1;i<DEBUG_TMP_N_CURL;++i)fout << " " << nTotalWire[i];
03293 }else{
03294 std::cout << "(TCurlFinder) tmp.wire.data can not be made!" << std::endl;
03295 }
03296 #endif
03297 return;
03298 }
03299
03300 #undef DEBUG_CURL_DUMP
03301 #undef DEBUG_CURL_SEGMENT
03302 #undef DEBUG_CURL_GNUPLOT
03303 #undef DEBUG_CURL_MC