BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
TBuilderCurl.cxx
Go to the documentation of this file.
1//-----------------------------------------------------------------------------
2// $Id: TBuilderCurl.cxx,v 1.24 2022/01/28 22:04:42 maqm Exp $
3//-----------------------------------------------------------------------------
4// Filename : TBuilderCurl.cc
5// Section : Tracking
6// Owner : Yoshi Iwasaki
7// Email : yoshihito.iwasaki@kek.jp
8//-----------------------------------------------------------------------------
9// Description : A class to build a curl track.
10// See http://bsunsrv1.kek.jp/~yiwasaki/tracking/
11//-----------------------------------------------------------------------------
12
13// #ifdef TRKRECO_DEBUG_DETAIL
14// #define TRKRECO_DEBUG
15// #endif
16#include "TrkReco/TBuilderCurl.h"
17#include "TrackUtil/Lpav.h"
18#include "TrkReco/TMLink.h"
19#include "TrkReco/TTrack.h"
20// #include "TrkReco/Lpav.h"
21// #include "TrkReco/TSvdFinder.h"
22// #include "TrkReco/TSvdHit.h"
23
24#include "TrkReco/TMLine.h"
25#include "TrkReco/TRobustLineFitter.h"
26
27// Following 3 parameters : (0,0,0) is best!
28// Other cases are for the development.
29#define DEBUG_CURL_DUMP 0
30#define DEBUG_CURL_GNUPLOT 0
31#define DEBUG_CURL_MC 0
32
33#if ( DEBUG_CURL_GNUPLOT + DEBUG_CURL_MC )
34# include "TrkReco/TMDCWireHitMC.h"
35# include "TrkReco/TTrackHEP.h"
36// #include MDC_H
37# include "MdcTables/MdcTables.h"
38# include <algo.h>
39#endif
40
41#include "GaudiKernel/Bootstrap.h"
42#include "GaudiKernel/IDataProviderSvc.h"
43#include "GaudiKernel/IInterface.h"
44#include "GaudiKernel/IMessageSvc.h"
45#include "GaudiKernel/ISvcLocator.h"
46#include "GaudiKernel/Kernel.h"
47#include "GaudiKernel/MsgStream.h"
48#include "GaudiKernel/Service.h"
49#include "GaudiKernel/SmartDataPtr.h"
50#include "GaudiKernel/StatusCode.h"
51
52#if ( DEBUG_CURL_DUMP + DEBUG_CURL_GNUPLOT + DEBUG_CURL_MC )
53// #include "tables/belletdf.h"
54// #include "tables/bestdf.h"
55static int debugMcFlag = 1;
56#endif
57
58TBuilderCurl::TBuilderCurl( const std::string& name )
59 : TBuilder0( name ), _fitter( "TBuilderCurl Fitter" ), m_pmgnIMF( nullptr ) {
60
61#if 0
62// if(m_param.ON_CORRECTION){
63// _fitter.sag(true);
64// _fitter.propagation(true);
65// _fitter.tof(true);
66// }
67 if (m_param.ON_CORRECTION & 1) _fitter.sag(true);
68 if (m_param.ON_CORRECTION & 2) _fitter.propagation(true);
69 if (m_param.ON_CORRECTION & 4) _fitter.tof(true);
70#endif
71}
72
73TBuilderCurl::~TBuilderCurl() {
74 // if(m_param.SVD_RECONSTRUCTION){
75 // delete m_svdFinder;
76 // delete m_svdAssociator;
77 // }
78}
79
80void TBuilderCurl::setParam( const TCurlFinderParameter& p ) {
81 m_param.Z_CUT = p.Z_CUT;
82 m_param.Z_DIFF_FOR_LAST_ATTEND = p.Z_DIFF_FOR_LAST_ATTEND;
83 m_param.SELECTOR_MAX_SIGMA = p.SELECTOR_MAX_SIGMA;
84 m_param.SELECTOR_MAX_IMPACT = p.SELECTOR_MAX_IMPACT;
85 m_param.SELECTOR_STRANGE_PZ = p.SELECTOR_STRANGE_PZ;
86 m_param.SELECTOR_REPLACE_DZ = p.SELECTOR_REPLACE_DZ;
87 m_param.SVD_RECONSTRUCTION = p.SVD_RECONSTRUCTION;
88 m_param.MIN_SVD_ELECTRONS = p.MIN_SVD_ELECTRONS;
89 m_param.ON_CORRECTION = p.ON_CORRECTION;
90 m_param.CURL_VERSION = p.CURL_VERSION;
91#if 1
92 // if(m_param.ON_CORRECTION){
93 // _fitter.sag(true);
94 // _fitter.propagation(true);
95 // _fitter.tof(true);
96 // }
97 if ( m_param.ON_CORRECTION & 1 ) _fitter.sag( true );
98 if ( m_param.ON_CORRECTION & 2 ) _fitter.propagation( true );
99 if ( m_param.ON_CORRECTION & 4 ) _fitter.tof( true );
100#endif
101 // if(m_param.SVD_RECONSTRUCTION){
102 ////m_svdFinder = new TSvdFinder(-1.*(m_param.MIN_SVD_ELECTRONS),
103 //// m_param.MIN_SVD_ELECTRONS);
104 // m_svdAssociator = new TSvdAssociator(-1.*(m_param.MIN_SVD_ELECTRONS),
105 // m_param.MIN_SVD_ELECTRONS);
106 // }
107}
108
109// Stereo Finder For Curl Tracks : by jtanaka == from here ==
110
111void TBuilderCurl::resetHelixFit( THelixFitter* fit ) const {
112 fit->fit2D( false );
113 fit->sag( false );
114 fit->propagation( false );
115 fit->tof( false );
116 fit->freeT0( false );
117}
118
119//.............................................
120//...............Main Part.....................
121//.............................................
122
123/*
124bool TBuilderCurl::buildStereo( TTrack& track, double& dZ, double& tanL ) const {
125 if(!(m_param.SVD_RECONSTRUCTION))return false;
126
127 #if 0
128 double q = 1.;
129 if(track.helix().kappa()<0.)q = -1.;
130 double h[5], chisq;
131 if(m_svdFinder->recTrk(track.helix().center().x(),
132 track.helix().center().y(),
133 fabs(track.helix().radius()),
134 q,
135 0.2,
136 h[0], h[1], h[2], h[3], h[4], chisq)){
137 #if 0
138 std::cout << "SVD = "
139 << h[0] << " "
140 << h[1] << " "
141 << h[2] << " "
142 << h[3] << " "
143 << h[4];
144 std::cout << ", chi2 = " << chisq << " pt=" << 1./h[2] << std::endl;
145 #endif
146 dZ = h[3];
147 tanL = h[4];
148 return true;
149 }else{
150 //std::cout << "SVD Fail....." << std::endl;
151 return false;
152 }
153 #else
154 Helix th(track.helix());
155 th.pivot(HepPoint3D(0.,0.,0.));
156 AList<TSvdHit> cand;
157 if(m_svdAssociator->recTrk(th,dZ,tanL,0.5,-1.0,cand,0.5)){
158 //std::cout << "SVD " << dZ << " " << tanL << std::endl;
159 return true;
160 }else{
161 //std::cout << "SVD..." << std::endl;
162 return false;
163 }
164 #endif
165}
166*/
167
168TTrack* TBuilderCurl::buildStereo( TTrack& track, const AList<TMLink>& stereoList ) const {
169 return NULL;
170}
171
172TTrack* TBuilderCurl::buildStereo( TTrack& track, const AList<TMLink>& stereoList,
173 const AList<TMLink>& allStereoList ) const {
174#if ( DEBUG_CURL_DUMP + DEBUG_CURL_GNUPLOT + DEBUG_CURL_MC )
175 Belle_event_Manager& evtMgr = Belle_event_Manager::get_manager();
176 debugMcFlag = 1;
177 if ( evtMgr.count() != 0 && evtMgr[0].ExpMC() != 2 ) debugMcFlag = 0; // not MC
178#endif
179
180 AList<TMLink> list = stereoList;
181 unsigned nList = list.length();
182
183 //...gets stereo wires from track.links
184 for ( unsigned i = 0, size = track.links().length(); i < size; ++i )
185 {
186 unsigned superID = ( track.links() )[i]->wire()->superLayerId();
187 if ( superID == 0 || superID == 1 || superID == 5 || superID == 6 || superID == 7 ||
188 superID == 8 )
189 {
190 int ok = 1;
191 for ( unsigned j = 0; j < nList; ++j )
192 {
193 TMLink* l = list[j];
194 if ( l->hit()->wire()->id() == ( track.links() )[i]->hit()->wire()->id() )
195 {
196 ok = 0;
197 break;
198 }
199 }
200 if ( ok == 1 ) list.append( ( ( track.links() )[i] ) );
201 }
202 // set LR 2
203 ( track.links() )[i]->leftRight( 2 ); // returns left-right. 0:left, 1:right, 2:wire
204 }
205 for ( unsigned i = 0, size = list.length(); i < size; ++i )
206 {
207 track._links.remove( list[i] );
208 // set LR 2
209 list[i]->leftRight( 2 );
210 }
211
212 if ( list.length() < 2 || list.length() + track.nLinks() < 5 ) return NULL;
213#if DEBUG_CURL_DUMP
214 unsigned debug_stereo_counter1 = 0;
215 for ( unsigned i = 0; i < track.links().length(); ++i )
216 {
217 unsigned superID = ( track.links() )[i]->hit()->wire()->superLayerId();
218 if ( superID == 0 || superID == 1 || superID == 5 || superID == 6 || superID == 7 ||
219 superID == 8 )
220 ++debug_stereo_counter1;
221 }
222 std::cout << "(TBuilderCurl)Fitted Track:";
223 std::cout << ", A# = " << track.links().length() - debug_stereo_counter1;
224 std::cout << ", S# = " << debug_stereo_counter1 << "(==0)";
225 std::cout << ", A+S# = " << track.links().length();
226 std::cout << ", Cand Stereo Wires # = " << list.length() << std::endl;
227 double debugChg = -1.;
228 if ( track.charge() > 0. ) debugChg = 1.;
229 if ( debugChg > 0. ) std::cout << "(TBuilderCurl) Positive Track" << std::endl;
230 else std::cout << "(TBuilderCurl) Negative Track" << std::endl;
231#endif
232
233 //...calculates a circle.
234 double xc2D;
235 double yc2D;
236 double r2D;
237 AList<TMLink> tmpAxialList = track.links();
238 bool err2D = fitWDD( xc2D, yc2D, r2D, tmpAxialList );
239
240 //...using SVD //....Liuqg
241 // double dZSVD, tanLSVD;
242 // bool initWithSVD = buildStereo(track, dZSVD, tanLSVD);
243
244 //...sets arc and z pairs of each stereo hit wire.
245 setArcZ( track, list );
246 AList<TMLink> removeList;
247 for ( unsigned i = 0, size = list.length(); i < size; ++i )
248 {
249 if ( list[i]->arcZ( 0 ).x() == -999. && list[i]->arcZ( 1 ).x() == -999. )
250 removeList.append( list[i] );
251 }
252 list.remove( removeList );
253 if ( list.length() < 2 || list.length() + track.nLinks() < 5 )
254 { // stereo >=2 && axial >= 3
255#if DEBUG_CURL_DUMP
256 std::cout << "(TBuilderCurl) Fail...few wires which can be set Arc-Z # = "
257 << list.length() << std::endl;
258#endif
259 return NULL;
260 }
261#if DEBUG_CURL_DUMP
262 std::cout << "(TBuilderCurl) Cand Stereo Wires which can be set Arc-Z # = " << list.length()
263 << std::endl;
264 plotArcZ( list, 0., 0., 0 );
265#endif
266
267#if 0
268 //...separates
269 AList<TMLink> layer[18];
270 for(unsigned i = 0, size = list.length(); i < size; ++i){
271 layer[list[i]->wire()->layer()->axialStereoLayerId()].append(*list[i]);
272 }
273
274# if DEBUG_CURL_DUMP
275 for(int i=0;i<18;++i){
276 if(layer[i].length())
277 std::cout << "(TBuilderCurl) Cand Stereo Wires which can be set Arc-Z # = "
278 << layer[i].length()
279 << " on " << i << " Layer." << std::endl;
280 }
281# endif
282#endif
283
284 // makes a line.
285 AList<TMLink> goodL;
286 AList<HepPoint3D> goodP;
287 double minChi2 = 9999.;
288 double goodA = 9999.;
289 double goodB = 9999.;
290 makeLine( track, list, allStereoList, goodL, minChi2, goodA, goodB, goodP );
291 HepAListDeleteAll( goodP );
292
293#if DEBUG_CURL_DUMP
294 std::cout << "(TBuilderCurl) a = " << goodA << ", b = " << goodB
295 << " (after makeLine-function)" << std::endl;
296#endif
297
298 // refits
299 static TRobustLineFitter rfitter( "Can you work well?" );
300 TMLine newLine0( goodL );
301 if ( rfitter.fit( newLine0 ) != 0 )
302 {
303#if DEBUG_CURL_DUMP
304 std::cout << "(TBuilderCurl) Fail...linefitting...fail." << std::endl;
305#endif
306 return NULL;
307 }
308
309#if DEBUG_CURL_DUMP
310 std::cout << "(TBuilderCurl) a = " << newLine0.a() << ", b = " << newLine0.b()
311 << " (after robustline-fit)" << std::endl;
312#endif
313
314 //...appends at last chance
315 unsigned nGoodL = goodL.length();
316 list.remove( goodL );
317 AList<TMLink> goodL2;
318 if ( m_param.CURL_VERSION == 0 )
319 { // default
320 if ( fabs( newLine0.b() ) < 10. )
321 {
322 appendPoints( list, goodL2, newLine0.a(), newLine0.b(), track,
323 m_param.Z_DIFF_FOR_LAST_ATTEND );
324 }
325 else
326 { // in case of bad result of robustLineFitter. (2001/04/04)
327 appendPoints( list, goodL2, goodA, goodB, track, m_param.Z_DIFF_FOR_LAST_ATTEND );
328 }
329 }
330 else
331 {
332 // same with b20010409_2122 at least
333 appendPoints( list, goodL2, newLine0.a(), newLine0.b(), track,
334 m_param.Z_DIFF_FOR_LAST_ATTEND );
335 }
336 goodL.append( goodL2 );
337 TMLine newLine( goodL );
338
339 //...refits
340 if ( rfitter.fit( newLine ) != 0 )
341 {
342#if DEBUG_CURL_DUMP
343 std::cout << "(TBuilderCurl) Fail...appending and re-fitting...fail." << std::endl;
344#endif
345 return NULL;
346 }
347
348#if DEBUG_CURL_DUMP
349 std::cout << "(TBuilderCurl) a = " << newLine.a() << ", b = " << newLine.b()
350 << " (after last-append + re-robustline-fit)" << std::endl;
351#endif
352 //...makes 3D tracks
353 const AList<TMLink>& good = newLine.links();
354 track.TTrackBase::append( good );
355 if ( !check( track ) )
356 {
357#if DEBUG_CURL_DUMP
358 std::cout << "(TBuilderCurl) Fail...checking wire numbers...fail." << std::endl;
359#endif
360 return NULL;
361 }
362
363 //...sets initial values
364 Vector a = track.helix().a();
365 if ( err2D )
366 {
367 double tmpA[3];
368 double tmpQ = track.charge() > 0. ? 1. : -1.;
369 tmpA[1] = fmod( atan2( tmpQ * yc2D, tmpQ * xc2D ) + 4. * M_PI, 2. * M_PI );
370 // tmpA[2] = tmpQ*Helix::ConstantAlpha/r2D;
371 // tmpA[2] = tmpQ*333.564095/r2D;
372 tmpA[2] = tmpQ * ( 333.564095 / ( -1000 * ( getPmgnIMF()->getReferField() ) ) ) / r2D;
373 // cout<<"TBuilderCurl::tmpA[2]:
374 // "<<(333.564095/(-1000*(getPmgnIMF()->getReferField())))<<endl;
375 tmpA[0] = xc2D / cos( tmpA[1] ) - tmpQ * r2D;
376 // std::cout << yc2D/sin(tmpA[1])-tmpQ*r2D << std::endl;
377 // std::cout << a[0] << " -0- " << tmpA[0] << std::endl;
378 // std::cout << a[1] << " -1- " << tmpA[1] << std::endl;
379 // std::cout << a[2] << " -2- " << tmpA[2] << std::endl;
380 a[0] = tmpA[0];
381 a[1] = tmpA[1];
382 a[2] = tmpA[2];
383 }
384 a[3] = newLine.b();
385 a[4] = newLine.a();
386 // Liuqg
387 /* if(initWithSVD){ // use SVD initial values if possible.
388 a[3] = dZSVD;
389 a[4] = tanLSVD;
390 }
391 */
392 if ( m_param.CURL_VERSION == 0 )
393 {
394 if ( fabs( a[3] ) > 10. && fabs( goodB ) < 10. )
395 { // 50cm --> 10cm @ 2001/04/04
396 // use initial values when results of RobustFit is bad.
397 a[3] = goodB;
398 a[4] = goodA;
399 }
400 }
401 else
402 {
403 if ( fabs( a[3] ) > 50. && fabs( goodB ) < 50. )
404 {
405 a[3] = goodB;
406 a[4] = goodA;
407 track.TTrackBase::remove( goodL2 );
408 }
409 }
410
411 track._helix->a( a );
412
413#if DEBUG_CURL_DUMP
414 std::cout << "(TBuilderCurl) Created Line: y = " << newLine.a() << " * x + "
415 << newLine.b() << ", size = " << good.length() << std::endl;
416 plotArcZ( const_cast<AList<TMLink>&>( good ), newLine.a(), newLine.b(), 0 );
417#endif
418
419 /* std::cout << "1st : "
420 << track.helix().a()[0] << " " << track.helix().a()[1] << " "
421 << track.helix().a()[2] << " " << track.helix().a()[3] << " "
422 << track.helix().a()[4] << std::endl; */
423
424#if 1
425 // if(m_param.ON_CORRECTION){
426 // _fitter.sag();
427 // _fitter.propagation();
428 // _fitter.tof();
429 // }
430 if ( m_param.ON_CORRECTION & 1 ) _fitter.sag( true );
431 if ( m_param.ON_CORRECTION & 2 ) _fitter.propagation( true );
432 if ( m_param.ON_CORRECTION & 4 ) _fitter.tof( true );
433#endif
434
435 //...fits
436 AList<TMLink> bad;
437 int err = _fitter.fit( track );
438 if ( err < 0 )
439 {
440#if DEBUG_CURL_DUMP
441 std::cout << "(TBuilderCurl) Fail fitting(0)...error code = " << err << std::endl;
442#endif
443 return NULL;
444 }
445 else if ( fabs( track.helix().a()[3] ) > fabs( a[3] ) &&
446 ( fabs( track.helix().a()[3] ) > 50. || // 50cm
447 fabs( track.helix().a()[2] ) > 100. || // 0.01GeV
448 fabs( track.helix().a()[2] ) < 0.1 ) )
449 { // 10GeV
450 // in strange case, set "correction" of wires OFF
451 // and then fit with the initial values.
452 if ( m_param.ON_CORRECTION )
453 {
454 if ( m_param.ON_CORRECTION & 1 ) _fitter.sag( true );
455 if ( m_param.ON_CORRECTION & 2 ) _fitter.propagation( true );
456 if ( m_param.ON_CORRECTION & 4 ) _fitter.tof( true );
457 //_fitter.sag();
458 //_fitter.propagation();
459 //_fitter.tof();
460 track._helix->a( a );
461 // std::cout << "ON --> OFF" << std::endl;
462 }
463 else { return NULL; }
464 }
465
466 /* std::cout << "2nd : "
467 << track.helix().a()[0] << " " << track.helix().a()[1] << " "
468 << track.helix().a()[2] << " " << track.helix().a()[3] << " "
469 << track.helix().a()[4] << std::endl; */
470
471 track.refine( bad, m_param.SELECTOR_MAX_SIGMA * 100. );
472 if ( !check( track ) )
473 {
474#if DEBUG_CURL_DUMP
475 std::cout << "(TBuilderCurl) Fail checking(1)..." << std::endl;
476#endif
477 return NULL;
478 }
479 err = _fitter.fit( track );
480 if ( err < 0 )
481 {
482#if DEBUG_CURL_DUMP
483 std::cout << "(TBuilderCurl) Fail fitting(1)...error code = " << err << std::endl;
484#endif
485 return NULL;
486 }
487 track.refine( bad, m_param.SELECTOR_MAX_SIGMA * 10. );
488 if ( !check( track ) )
489 {
490#if DEBUG_CURL_DUMP
491 std::cout << "(TBuilderCurl) Fail checking(2)..." << std::endl;
492#endif
493 return NULL;
494 }
495 err = _fitter.fit( track );
496 if ( err < 0 )
497 {
498#if DEBUG_CURL_DUMP
499 std::cout << "(TBuilderCurl) Fail fitting(2)...error code = " << err << std::endl;
500#endif
501 return NULL;
502 }
503 track.refine( bad, m_param.SELECTOR_MAX_SIGMA );
504 if ( !check( track ) )
505 {
506#if DEBUG_CURL_DUMP
507 std::cout << "(TBuilderCurl) Fail checking(3)..." << std::endl;
508#endif
509 return NULL;
510 }
511 err = _fitter.fit( track );
512 if ( err < 0 )
513 {
514#if DEBUG_CURL_DUMP
515 std::cout << "(TBuilderCurl) Fail fitting(3)...error code = " << err << std::endl;
516#endif
517 return NULL;
518 }
519
520#if 0
521 for(int i=0;i<track.links().length();++i){
522 if(track.links()[i]->pull() > 36.){
523 std::cout << track.links()[i]->wire()->id()
524 << " :+: " << track.links()[i]->pull() << std::endl;
525 }
526 if(track.links()[i]->hit()->state() & WireHitInvalidForFit)
527 std::cout << "Not Valid For Fit!" << std::endl;
528 //if(!(track.links()[i]->hit()->state() & WireHitFittingValid))
529 //std::cout << "No-Valid For Fit!" << std::endl;
530 }
531#endif
532 track.refine( bad, m_param.SELECTOR_MAX_SIGMA );
533 if ( !check( track ) )
534 {
535#if DEBUG_CURL_DUMP
536 std::cout << "(TBuilderCurl) Fail checking(4)..." << std::endl;
537#endif
538 return NULL;
539 }
540 err = _fitter.fit( track );
541 if ( err < 0 )
542 {
543#if DEBUG_CURL_DUMP
544 std::cout << "(TBuilderCurl) Fail fitting(4)...error code = " << err << std::endl;
545#endif
546 return NULL;
547 }
548
549#if 0
550 for(int i=0;i<track.links().length();++i){
551 if(track.links()[i]->pull() > 36.){
552 std::cout << track.links()[i]->wire()->id()
553 << " :*: " << track.links()[i]->pull() << std::endl;
554 }
555 if(track.links()[i]->hit()->state() & WireHitInvalidForFit)
556 std::cout << "Not Valid For Fit!" << std::endl;
557 //if(!(track.links()[i]->hit()->state() & WireHitFittingValid))
558 //std::cout << "No-Valid For Fit!" << std::endl;
559 }
560#endif
561
562 /* std::cout << "3rd : "
563 << track.helix().a()[0] << " " << track.helix().a()[1] << " "
564 << track.helix().a()[2] << " " << track.helix().a()[3] << " "
565 << track.helix().a()[4] << std::endl; */
566
567 //...tests
568 if ( track.nLinks() < 5 )
569 { // axial + stereo >= 5
570#if DEBUG_CURL_DUMP
571 std::cout << "(TBuilderCurl) Success fitting, but pre-selection...fail." << std::endl;
572#endif
573 return NULL;
574 }
575 if ( fabs( track.impact() ) > m_param.SELECTOR_MAX_IMPACT ||
576 track.pt() < 0.005 || // Pt >= 5MeV
577 fabs( track.pz() ) > m_param.SELECTOR_STRANGE_PZ )
578 {
579#if DEBUG_CURL_DUMP
580 std::cout << "(TBuilderCurl) Success fitting, but selection...fail." << std::endl;
581 std::cout << "(TBuilderCurl) impact = " << track.impact()
582 << std::endl;
583 std::cout << "(TBuilderCurl) pt = " << track.pt()
584 << std::endl;
585 std::cout << "(TBuilderCurl) pz = " << track.pz()
586 << std::endl;
587 std::cout << "(TBuilderCurl) dz = " << track.helix().a()[3]
588 << std::endl;
589#endif
590 return NULL;
591 }
592
593 //...replaces init helix if dz is bad.
594 if ( fabs( track.helix().a()[3] ) > m_param.SELECTOR_REPLACE_DZ &&
595 fabs( a[3] ) < fabs( track.helix().a()[3] ) )
596 { track._helix->a( a ); }
597#if DEBUG_CURL_DUMP
598 std::cout << "(TBuilderCurl) Success Build Stereo!!!" << std::endl;
599#endif
600 return &track;
601}
602
603//.............................................
604//.............Set Arc and Z Part..............
605//.............................................
606
607void TBuilderCurl::setArcZ( TTrack& track, AList<TMLink>& list ) const {
608 if ( track.nLinks() < 4 )
609 {
610 track.stereoHitForCurl( list );
611 return;
612 }
613 // Liuqg
614 AList<TMLink> alayer[5];
615 AList<TMLink> slayer[6];
616 for ( unsigned i = 0, size = track.nLinks(); i < size; ++i )
617 {
618 unsigned id = ( track.links() )[i]->wire()->superLayerId();
619 if ( id == 2 ) alayer[0].append( ( track.links() )[i] );
620 else if ( id == 3 ) alayer[1].append( ( track.links() )[i] );
621 else if ( id == 4 ) alayer[2].append( ( track.links() )[i] );
622 else if ( id == 9 ) alayer[3].append( ( track.links() )[i] );
623 else if ( id == 10 ) alayer[4].append( ( track.links() )[i] );
624 }
625
626 for ( unsigned i = 0, size = list.length(); i < size; ++i )
627 {
628 unsigned id = list[i]->wire()->superLayerId();
629 if ( id == 0 ) slayer[0].append( list[i] );
630 else if ( id == 1 ) slayer[1].append( list[i] );
631 else if ( id == 5 ) slayer[2].append( list[i] );
632 else if ( id == 6 ) slayer[3].append( list[i] );
633 else if ( id == 7 ) slayer[4].append( list[i] );
634 else if ( id == 8 ) slayer[5].append( list[i] );
635 }
636#if 0
637 std::cout << "Stereo = "
638 << slayer[0].length() << " "
639 << slayer[1].length() << " "
640 << slayer[2].length() << " "
641 << slayer[3].length() << " "
642 << slayer[4].length() << std::endl;
643 std::cout << "Axial = "
644 << alayer[0].length() << " "
645 << alayer[1].length() << " "
646 << alayer[2].length() << " "
647 << alayer[3].length() << " "
648 << alayer[4].length() << " "
649 << alayer[5].length() << std::endl;
650#endif
651 unsigned ip = 0;
652 if ( slayer[0].length() >= 1 )
653 {
654 if ( alayer[0].length() + alayer[1].length() >= 4 )
655 { setArcZ( track, slayer[0], alayer[0], alayer[1], ip ); }
656 else if ( alayer[0].length() + alayer[1].length() + alayer[2].length() >= 4 )
657 { setArcZ( track, slayer[0], alayer[0], alayer[1], alayer[2], ip ); }
658 else if ( alayer[0].length() + alayer[1].length() + alayer[2].length() +
659 alayer[3].length() >=
660 4 )
661 { setArcZ( track, slayer[0], alayer[0], alayer[1], alayer[2], alayer[3], ip ); }
662 else if ( alayer[0].length() + alayer[1].length() + alayer[2].length() +
663 alayer[3].length() + alayer[4].length() >=
664 4 )
665 { setArcZ( track, slayer[0], alayer[0], alayer[1], alayer[2], alayer[3], alayer[4], ip ); }
666 }
667 if ( slayer[1].length() >= 1 )
668 {
669 if ( alayer[0].length() + alayer[1].length() >= 4 )
670 { setArcZ( track, slayer[1], alayer[0], alayer[1], ip ); }
671 else if ( alayer[0].length() + alayer[1].length() + alayer[2].length() >= 4 )
672 { setArcZ( track, slayer[1], alayer[0], alayer[1], alayer[2], ip ); }
673 else if ( alayer[0].length() + alayer[1].length() + alayer[2].length() +
674 alayer[3].length() >=
675 4 )
676 { setArcZ( track, slayer[1], alayer[0], alayer[1], alayer[2], alayer[3], ip ); }
677 else if ( alayer[0].length() + alayer[1].length() + alayer[2].length() +
678 alayer[3].length() + alayer[4].length() >=
679 4 )
680 { setArcZ( track, slayer[1], alayer[0], alayer[1], alayer[2], alayer[3], alayer[4], ip ); }
681 }
682 if ( slayer[2].length() >= 1 )
683 {
684 if ( alayer[1].length() + alayer[2].length() >= 4 )
685 { setArcZ( track, slayer[2], alayer[1], alayer[2], ip ); }
686 else if ( alayer[0].length() + alayer[1].length() + alayer[2].length() >= 4 )
687 { setArcZ( track, slayer[2], alayer[0], alayer[1], alayer[2], ip ); }
688 else if ( alayer[0].length() + alayer[1].length() + alayer[2].length() +
689 alayer[3].length() >=
690 4 )
691 { setArcZ( track, slayer[2], alayer[0], alayer[1], alayer[2], alayer[3], ip ); }
692 else if ( alayer[0].length() + alayer[1].length() + alayer[2].length() +
693 alayer[3].length() + alayer[4].length() >=
694 4 )
695 { setArcZ( track, slayer[2], alayer[0], alayer[1], alayer[2], alayer[3], alayer[4], ip ); }
696 }
697 if ( slayer[3].length() >= 1 )
698 {
699 if ( alayer[1].length() + alayer[2].length() >= 4 )
700 { setArcZ( track, slayer[3], alayer[1], alayer[2], ip ); }
701 else if ( alayer[0].length() + alayer[1].length() + alayer[2].length() >= 4 )
702 { setArcZ( track, slayer[3], alayer[0], alayer[1], alayer[2], ip ); }
703 else if ( alayer[0].length() + alayer[1].length() + alayer[2].length() +
704 alayer[3].length() >=
705 4 )
706 { setArcZ( track, slayer[3], alayer[0], alayer[1], alayer[2], alayer[3], ip ); }
707 else if ( alayer[0].length() + alayer[1].length() + alayer[2].length() +
708 alayer[3].length() + alayer[4].length() >=
709 4 )
710 { setArcZ( track, slayer[3], alayer[0], alayer[1], alayer[2], alayer[3], alayer[4], ip ); }
711 }
712 if ( slayer[4].length() >= 1 )
713 {
714 if ( alayer[1].length() + alayer[2].length() >= 4 )
715 { setArcZ( track, slayer[4], alayer[1], alayer[2], ip ); }
716 else if ( alayer[0].length() + alayer[1].length() + alayer[2].length() >= 4 )
717 { setArcZ( track, slayer[4], alayer[0], alayer[1], alayer[2], ip ); }
718 else if ( alayer[0].length() + alayer[1].length() + alayer[2].length() +
719 alayer[3].length() >=
720 4 )
721 { setArcZ( track, slayer[4], alayer[0], alayer[1], alayer[2], alayer[3], ip ); }
722 else if ( alayer[0].length() + alayer[1].length() + alayer[2].length() +
723 alayer[3].length() + alayer[4].length() >=
724 4 )
725 { setArcZ( track, slayer[4], alayer[0], alayer[1], alayer[2], alayer[3], alayer[4], ip ); }
726 }
727 if ( slayer[5].length() >= 1 )
728 {
729 if ( alayer[1].length() + alayer[2].length() >= 4 )
730 { setArcZ( track, slayer[5], alayer[1], alayer[2], ip ); }
731 else if ( alayer[1].length() + alayer[2].length() + alayer[3].length() >= 4 )
732 { setArcZ( track, slayer[5], alayer[1], alayer[2], alayer[3], ip ); }
733 else if ( alayer[0].length() + alayer[1].length() + alayer[2].length() +
734 alayer[3].length() >=
735 4 )
736 { setArcZ( track, slayer[5], alayer[0], alayer[1], alayer[2], alayer[3], ip ); }
737 else if ( alayer[0].length() + alayer[1].length() + alayer[2].length() +
738 alayer[3].length() + alayer[4].length() >=
739 4 )
740 { setArcZ( track, slayer[5], alayer[0], alayer[1], alayer[2], alayer[3], alayer[4], ip ); }
741 }
742}
743//.............................................
744//.............Append Points(last).............
745//.............................................
746
747unsigned TBuilderCurl::appendPoints( AList<TMLink>& list, AList<TMLink>& line, double a,
748 double b, TTrack& track, double z_cut ) const {
749 unsigned size = list.length();
750 if ( size == 0 ) return 0;
751 unsigned counter( 0 );
752 for ( unsigned i = 0; i < size; ++i )
753 {
754 for ( unsigned j = 0; j < 4; ++j )
755 {
756 if ( j <= 1 && list[i]->arcZ( j ).x() == -999. ) continue;
757 else if ( j > 1 && list[i]->arcZ( j ).x() == -999. ) break;
758 double y = a * list[i]->arcZ( j ).x() + b;
759 if ( fabs( y - list[i]->arcZ( j ).y() ) < z_cut )
760 {
761 list[i]->position( list[i]->arcZ( j ) );
762 line.append( list[i] );
763 ++counter;
764 break;
765 }
766 }
767 }
768 return counter;
769}
770
771//.............................................
772//..............Line Fit Part..................
773//.............................................
774
775int doLineFit( AList<TMLink>& points, double& m_a, double& m_b, double& chi2, double& nhits,
776 int ipC = 1 ) // Liuqg, IP Constraint
777{
778 m_a = m_b = nhits = 0.;
779 chi2 = 1.e+10;
780 unsigned n = points.length();
781 double sum = double( n );
782 double sumX = 0., sumY = 0., sumX2 = 0., sumXY = 0., sumY2 = 0.;
783 for ( unsigned i = 0; i < n; i++ )
784 {
785 TMLink& l = *points[i];
786
787 double x = l.position().x();
788 double y = l.position().y();
789 sumX += x;
790 sumY += y;
791 sumX2 += x * x;
792 sumXY += x * y;
793 sumY2 += y * y;
794 }
795 if ( ipC != 0 ) sum += 1.0; // IP Constraint
796 nhits = sum;
797
798 double m_det = sum * sumX2 - sumX * sumX;
799 if ( m_det == 0. && sum != 2. ) { return -1; }
800 else if ( m_det == 0. && sum == 2. )
801 {
802 double x0 = points[0]->position().x();
803 double y0 = points[0]->position().y();
804 double x1 = points[1]->position().x();
805 double y1 = points[1]->position().y();
806 if ( x0 == x1 ) return -1;
807 m_a = ( y0 - y1 ) / ( x0 - x1 );
808 m_b = -m_a * x1 + y1;
809 chi2 = 0.;
810 return 0;
811 }
812 chi2 = 0.;
813 m_a = ( sumXY * sum - sumX * sumY ) / m_det;
814 m_b = ( sumX2 * sumY - sumX * sumXY ) / m_det;
815
816 for ( unsigned i = 0; i < n; i++ )
817 {
818 TMLink& l = *points[i];
819
820 double x = l.position().x();
821 double y = l.position().y();
822 double d = y - ( x * m_a + m_b );
823 chi2 += d * d;
824 }
825
826 return 0;
827}
828
829void TBuilderCurl::fitLine( AList<TMLink>& tmpLine, double& min_chi2, double& good_a,
830 double& good_b, AList<TMLink>& goodLine,
831 AList<HepPoint3D>& goodPosition, int& overCounter ) const {
832#if 0
833 // OLD
834 unsigned size = tmpLine.length();
835 TMLine line(tmpLine);
836 if(!(line.fit())){
837 double chi2 = line.chi2()/(static_cast<double>(size-2));
838 if(chi2 < min_chi2 && fabs(line.b()) < m_param.Z_CUT){
839 min_chi2 = chi2;
840 good_a = line.a();
841 good_b = line.b();
842 goodLine = tmpLine;
843 HepAListDeleteAll(goodPosition);
844 for(unsigned i=0;i<size;++i)
845 goodPosition.append(new HepPoint3D(const_cast<HepPoint3D&>(tmpLine[i]->position())));
846 }
847 }
848 ++overCounter;
849#else
850 unsigned size = tmpLine.length();
851 double ta, tb, tc, tn;
852 if ( !( doLineFit( tmpLine, ta, tb, tc, tn, 1 ) ) && tn > 2. )
853 { // with IP
854 double chi2 = tc / ( tn - 2. );
855 if ( chi2 < min_chi2 && fabs( tb ) < m_param.Z_CUT )
856 {
857 min_chi2 = chi2;
858 good_a = ta;
859 good_b = tb;
860 goodLine = tmpLine;
861 HepAListDeleteAll( goodPosition );
862 for ( unsigned i = 0; i < size; ++i )
863 goodPosition.append(
864 new HepPoint3D( const_cast<HepPoint3D&>( tmpLine[i]->position() ) ) );
865 }
866 }
867 ++overCounter;
868#endif
869}
870
871//.............................................
872//.................Check Part..................
873//.............................................
874
875unsigned TBuilderCurl::check( const TTrack& track ) const {
876 unsigned nAhits( 0 ), nShits( 0 );
877 for ( unsigned i = 0, size = track.nLinks(); i < size; ++i )
878 {
879 if ( !( track.links()[i]->hit()->state() & WireHitFittingValid ) ) continue;
880 if ( track.links()[i]->wire()->stereo() ) ++nShits;
881 else ++nAhits;
882 }
883 if ( nAhits >= 3 && nShits >= 2 ) return 1; // hard coding
884 return 0;
885}
886
887//.............................................
888//.............Checker(Debuger)................
889//.............................................
890
891int checkBorder( AList<TMLink>& layer0, AList<TMLink>& layer1, AList<TMLink>& layer2 ) {
892 AList<TMLink> list = layer0;
893 list.append( layer1 );
894 list.append( layer2 );
895 int size = list.length();
896 if ( size <= 1 ) return 0;
897 int layerId = list[0]->hit()->wire()->layerId();
898 int maxLocalId = 79;
899 if ( layerId >= 15 ) maxLocalId = 127;
900 if ( layerId >= 23 ) maxLocalId = 159;
901 if ( layerId >= 32 ) maxLocalId = 207;
902 if ( layerId >= 41 ) maxLocalId = 255;
903 int HId = (int)( 0.8 * ( maxLocalId + 1 ) );
904 int LId = (int)( 0.2 * ( maxLocalId + 1 ) );
905 int low = 0;
906 int high = 0;
907 for ( int i = 0; i > size; ++i )
908 {
909 if ( list[i]->hit()->wire()->localId() < LId ) low = 1;
910 else if ( list[i]->hit()->wire()->localId() > HId ) high = 1;
911 if ( low == 1 && high == 1 ) return 1;
912 }
913 return 0;
914}
915
916int checkBorder( AList<TMLink>& layer0, AList<TMLink>& layer1, AList<TMLink>& layer2,
917 AList<TMLink>& layer3 ) {
918 AList<TMLink> list = layer0;
919 list.append( layer1 );
920 list.append( layer2 );
921 list.append( layer3 );
922 int size = list.length();
923 if ( size <= 1 ) return 0;
924 int layerId = list[0]->hit()->wire()->layerId();
925 int maxLocalId = 79;
926 if ( layerId >= 15 ) maxLocalId = 127;
927 if ( layerId >= 23 ) maxLocalId = 159;
928 if ( layerId >= 32 ) maxLocalId = 207;
929 if ( layerId >= 41 ) maxLocalId = 255;
930 int HId = (int)( 0.8 * ( maxLocalId + 1 ) );
931 int LId = (int)( 0.2 * ( maxLocalId + 1 ) );
932 int low = 0;
933 int high = 0;
934 for ( int i = 0; i > size; ++i )
935 {
936 if ( list[i]->hit()->wire()->localId() < LId ) low = 1;
937 else if ( list[i]->hit()->wire()->localId() > HId ) high = 1;
938 if ( low == 1 && high == 1 ) return 1;
939 }
940 return 0;
941}
942
943int offsetBorder( TMLink* l ) {
944 int layerId = l->hit()->wire()->layerId();
945 int maxLocalId = 79;
946 if ( layerId >= 15 ) maxLocalId = 127;
947 if ( layerId >= 23 ) maxLocalId = 159;
948 if ( layerId >= 32 ) maxLocalId = 207;
949 if ( layerId >= 41 ) maxLocalId = 255;
950 return maxLocalId + 1;
951}
952
953void makeList( AList<TMLink>& layer, AList<TMLink>& list, double q, int border, int checkB,
954 TMLink* layer0 ) {
955 int n = layer.length();
956 if ( checkB == 0 )
957 {
958 for ( int i = 0; i < n; ++i )
959 {
960 if ( q < 0. )
961 {
962 if ( layer[i]->hit()->wire()->localId() >= border ) list.append( layer[i] );
963 }
964 else
965 {
966 if ( layer[i]->hit()->wire()->localId() <= border ) list.append( layer[i] );
967 }
968 }
969 }
970 else
971 {
972 // difficult!! --> puls offset
973 int offset = offsetBorder( layer0 );
974 if ( border * 2 < offset ) border += offset;
975 for ( int i = 0; i < n; ++i )
976 {
977 int lId = layer[i]->hit()->wire()->localId();
978 if ( lId * 2 < offset ) lId += offset;
979 if ( q < 0. )
980 {
981 if ( lId >= border ) list.append( layer[i] );
982 }
983 else
984 {
985 if ( lId <= border ) list.append( layer[i] );
986 }
987 }
988 }
989}
990
991int TBuilderCurl::sortByLocalId( AList<TMLink>& list ) const {
992 int size = list.length();
993 if ( size <= 1 ) return 0;
994 int layerId = list[0]->hit()->wire()->layerId();
995 int maxLocalId = 0; // add initialize 25-05-15
996 /* if(layerId < 15)maxLocalId = 79;
997 else if(layerId >= 15)maxLocalId = 127;
998 else if(layerId >= 23)maxLocalId = 159;
999 else if(layerId >= 32)maxLocalId = 207;
1000 else if(layerId >= 41)maxLocalId = 255;
1001 */
1002 // Liuqg in this class, the parameter of superlayer changed to layerid.
1003 if ( layerId == 0 ) maxLocalId = 39;
1004 else if ( layerId == 1 ) maxLocalId = 43;
1005 else if ( layerId == 2 ) maxLocalId = 47;
1006 else if ( layerId == 3 ) maxLocalId = 55;
1007 else if ( layerId == 4 ) maxLocalId = 63;
1008 else if ( layerId == 5 ) maxLocalId = 71;
1009 else if ( layerId < 8 ) maxLocalId = 79;
1010 else if ( layerId < 24 ) maxLocalId = 159;
1011 else if ( layerId < 28 ) maxLocalId = 175;
1012 else if ( layerId < 32 ) maxLocalId = 207;
1013 else if ( layerId < 43 ) maxLocalId = 239;
1014 int flag = 0;
1015 for ( int i = 0; i < size; ++i )
1016 {
1017 if ( list[i]->hit()->wire()->localId() == 0 || list[i]->hit()->wire()->localId() == 1 ||
1018 list[i]->hit()->wire()->localId() == maxLocalId - 1 ||
1019 list[i]->hit()->wire()->localId() == maxLocalId )
1020 {
1021 flag = 1;
1022 break;
1023 }
1024 }
1025 if ( flag == 0 ) return 0;
1026 int maxDif = (int)( 0.5 * ( maxLocalId + 1 ) );
1027 AList<TMLink> fList; // former
1028 AList<TMLink> lList; // later
1029 for ( int i = 0; i < size; ++i )
1030 {
1031 if ( list[i]->hit()->wire()->localId() < maxDif ) lList.append( list[i] );
1032 else fList.append( list[i] );
1033 }
1034 list.removeAll();
1035 list.append( fList );
1036 list.append( lList );
1037 if ( fList.length() >= 1 && lList.length() >= 1 ) return 1;
1038 else return 0;
1039}
1040
1041//..................................................
1042//..................................................
1043//....................MC............................
1044//..................................................
1045//..................................................
1046
1047TTrack* TBuilderCurl::buildStereoMC( TTrack& track, const AList<TMLink>& stereoList ) const {
1048#if DEBUG_CURL_MC
1049 AList<TMLink> list = stereoList;
1050
1051 HepPoint3D center = track.helix().center();
1052 double r = fabs( track.helix().curv() );
1053 for ( unsigned i = 0, size = list.length(); i < size; ++i )
1054 {
1055 HepPoint3D point(
1056 ( list[i]->hit()->mc()->datcdc()->m_xin + list[i]->hit()->mc()->datcdc()->m_xout ) *
1057 0.5,
1058 ( list[i]->hit()->mc()->datcdc()->m_yin + list[i]->hit()->mc()->datcdc()->m_yout ) *
1059 0.5,
1060 0. );
1061 double cosdPhi = -center.dot( ( point - center ).unit() ) / center.mag();
1062 double dPhi;
1063 if ( fabs( cosdPhi ) < 1.0 ) { dPhi = acos( cosdPhi ); }
1064 else if ( cosdPhi >= 1.0 ) { dPhi = 0.; }
1065 else { dPhi = M_PI; }
1066 list[i]->position( HepPoint3D(
1067 r * dPhi,
1068 ( list[i]->hit()->mc()->datcdc()->m_zin + list[i]->hit()->mc()->datcdc()->m_zout ) *
1069 0.5,
1070 0. ) );
1071 /* std::cout << list[i]->wire()->id() << ": "
1072 << point.x() << " "
1073 << point.y() << " "
1074 << (list[i]->hit()->mc()->datcdc()->m_zin+
1075 list[i]->hit()->mc()->datcdc()->m_zout)*0.5 << std::endl; */
1076 }
1077 // std::cout << "A# = " << track.links().length() << ", S# = " << list.length() << std::endl;
1078 double xc2D;
1079 double yc2D;
1080 double r2D;
1081 bool err2D = fitWDD( xc2D, yc2D, r2D, track.links() ); // axial only
1082
1083 TLine0 newLine( list );
1084 if ( newLine.fit() != 0 ) return NULL;
1085 const AList<TMLink>& good = newLine.links();
1086 track.append( good );
1087 Vector a = track.helix().a();
1088 a[3] = newLine.b();
1089 a[4] = newLine.a();
1090 if ( err2D )
1091 {
1092 double tmpA[3];
1093 double tmpQ = track.charge() > 0. ? 1. : -1.;
1094 tmpA[1] = fmod( atan2( tmpQ * yc2D, tmpQ * xc2D ) + 4. * M_PI, 2. * M_PI );
1095 // tmpA[2] = tmpQ*Helix::ConstantAlpha/r2D;
1096 // tmpA[2] = tmpQ*333.564095/r2D;
1097 tmpA[2] = tmpQ * ( 333.564095 / ( -1000 * ( getPmgnIMF()->getReferField() ) ) ) / r2D;
1098 tmpA[0] = xc2D / cos( tmpA[1] ) - tmpQ * r2D;
1099 a[0] = tmpA[0];
1100 a[1] = tmpA[1];
1101 a[2] = tmpA[2];
1102 }
1103 track._helix->a( a );
1104# if 0
1105 std::cout << track.helix().a()[0] << " " << track.helix().a()[1] << " "
1106 << track.helix().a()[2] << " " << track.helix().a()[3] << " "
1107 << track.helix().a()[4] << std::endl;
1108# endif
1109 //...fits
1110 AList<TMLink> bad;
1111 int err = _fitter.fit( track );
1112 if ( err < m_param.ERROR_FOR_HELIX_FIT ) return NULL;
1113 track.refine( bad, m_param.SELECTOR_MAX_SIGMA * 100. );
1114 err = _fitter.fit( track );
1115 if ( err < m_param.ERROR_FOR_HELIX_FIT ) return NULL;
1116 track.refine( bad, m_param.SELECTOR_MAX_SIGMA * 10. );
1117 err = _fitter.fit( track );
1118 if ( err < m_param.ERROR_FOR_HELIX_FIT ) return NULL;
1119 track.refine( bad, m_param.SELECTOR_MAX_SIGMA );
1120 err = _fitter.fit( track );
1121 if ( err < m_param.ERROR_FOR_HELIX_FIT ) return NULL;
1122# if 0
1123 std::cout << track.helix().a()[0] << " " << track.helix().a()[1] << " "
1124 << track.helix().a()[2] << " " << track.helix().a()[3] << " "
1125 << track.helix().a()[4] << std::endl;
1126# endif
1127 // track._helix->a(a); // re-write
1128
1129 //...checks
1130 if ( track.nLinks() < m_param.SELECTOR_NLINKS ) return NULL;
1131 if ( fabs( track.impact() ) > m_param.SELECTOR_MAX_IMPACT ||
1132 track.pt() < m_param.SELECTOR_MIN_PT )
1133 return NULL;
1134 return &track;
1135#else
1136 return NULL;
1137#endif
1138}
1139
1140//..................................................
1141//....................Plot..........................
1142//..................................................
1143
1144void TBuilderCurl::plotArcZ( AList<TMLink>& tmpLine, double a, double b,
1145 const int flag ) const {
1146#if DEBUG_CURL_GNUPLOT
1147 // #if 1
1148 if ( a == 9999. || b == 9999. )
1149 {
1150 a = 0.;
1151 b = 0.;
1152 }
1153 int nCounter = 0;
1154 double gmaxX = 0., gminX = 0.;
1155 double gmaxY = 0., gminY = 0.;
1156 FILE * gnuplot, *data;
1157 if ( ( data = fopen( "you_can_remove_this.dat", "w" ) ) != NULL )
1158 {
1159 for ( int ii = 0; ii < tmpLine.length(); ii++ )
1160 {
1161 ++nCounter;
1162 fprintf( data, "%lf, %lf\n", tmpLine[ii]->position().x(), tmpLine[ii]->position().y() );
1163 if ( flag )
1164 std::cout << " Wire ID = " << tmpLine[ii]->hit()->wire()->id()
1165 << " Arc = " << tmpLine[ii]->position().x()
1166 << " Z = " << tmpLine[ii]->position().y();
1167 // if(flag && !debugMcFlag)std::cout << std::endl;
1168 // if(flag && debugMcFlag){
1169 // std::cout << " Z(true) = " << (tmpLine[ii]->hit()->mc()->datcdc()->m_zin+
1170 // tmpLine[ii]->hit()->mc()->datcdc()->m_zout)*0.5;
1171 // std::cout << " HepTrackID = " << tmpLine[ii]->hit()->mc()->hep()->id() << std::endl;
1172 // }
1173 if ( gmaxX < tmpLine[ii]->position().x() ) gmaxX = tmpLine[ii]->position().x();
1174 if ( gminX > tmpLine[ii]->position().x() ) gminX = tmpLine[ii]->position().x();
1175 if ( gmaxY < tmpLine[ii]->position().y() ) gmaxY = tmpLine[ii]->position().y();
1176 if ( gminY > tmpLine[ii]->position().y() ) gminY = tmpLine[ii]->position().y();
1177 }
1178 fclose( data );
1179 }
1180 if ( ( data = fopen( "you_can_remove_this.dat2", "w" ) ) != NULL )
1181 {
1182# if 0
1183 if(debugMcFlag){
1184 for(int ii=0;ii<tmpLine.length();ii++){
1185 double z = (tmpLine[ii]->hit()->mc()->datcdc()->m_zin+
1186 tmpLine[ii]->hit()->mc()->datcdc()->m_zout)*0.5;
1187 if(tmpLine[ii]->arcZ(0).x() != -999.){
1188 ++nCounter;
1189 fprintf(data,"%lf, %lf\n",tmpLine[ii]->arcZ(0).x(),z);
1190 if(gmaxX < tmpLine[ii]->arcZ(0).x())
1191 gmaxX = tmpLine[ii]->arcZ(0).x();
1192 if(gminX > tmpLine[ii]->arcZ(0).x())
1193 gminX = tmpLine[ii]->arcZ(0).x();
1194 if(gmaxY < z)
1195 gmaxY = z;
1196 if(gminY > z)
1197 gminY = z;
1198 }
1199 }
1200 }
1201# endif
1202 fclose( data );
1203 }
1204 if ( ( data = fopen( "you_can_remove_this.dat3", "w" ) ) != NULL )
1205 {
1206 for ( int ii = 0; ii < tmpLine.length(); ii++ )
1207 {
1208 for ( int jj = 0; jj < 4; ++jj )
1209 {
1210 if ( tmpLine[ii]->arcZ( jj ).x() != -999. )
1211 {
1212 ++nCounter;
1213 fprintf( data, "%lf, %lf\n", tmpLine[ii]->arcZ( jj ).x(),
1214 tmpLine[ii]->arcZ( jj ).y() );
1215 if ( gmaxX < tmpLine[ii]->arcZ( jj ).x() ) gmaxX = tmpLine[ii]->arcZ( jj ).x();
1216 if ( gminX > tmpLine[ii]->arcZ( jj ).x() ) gminX = tmpLine[ii]->arcZ( jj ).x();
1217 if ( gmaxY < tmpLine[ii]->arcZ( jj ).y() ) gmaxY = tmpLine[ii]->arcZ( jj ).y();
1218 if ( gminY > tmpLine[ii]->arcZ( jj ).y() ) gminY = tmpLine[ii]->arcZ( jj ).y();
1219 }
1220 else { break; }
1221 }
1222 }
1223 fclose( data );
1224 }
1225 if ( gmaxX < 0. ) gmaxX = 0.;
1226 if ( gminX > 0. ) gminX = 0.;
1227 if ( gmaxY < 0. ) gmaxY = 0.;
1228 if ( gminY > 0. ) gminY = 0.;
1229 if ( nCounter && ( gnuplot = popen( "gnuplot", "w" ) ) != NULL )
1230 {
1231 fprintf( gnuplot, "set nokey \n" );
1232 fprintf( gnuplot, "set size 0.721,1.0 \n" );
1233 fprintf( gnuplot, "set xrange [%f:%f] \n", gminX, gmaxX );
1234 fprintf( gnuplot, "set yrange [%f:%f] \n", gminY, gmaxY );
1235 if ( a == 0. && b == 0. )
1236 {
1237 fprintf( gnuplot, "plot \"you_can_remove_this.dat3\", \"you_can_remove_this.dat\", "
1238 "\"you_can_remove_this.dat2\" \n" );
1239 }
1240 else
1241 {
1242 fprintf( gnuplot,
1243 "plot \"you_can_remove_this.dat3\", \"you_can_remove_this.dat\", "
1244 "\"you_can_remove_this.dat2\", %lf*x+%lf \n",
1245 a, b );
1246 }
1247 fflush( gnuplot );
1248 char tmp[8];
1249 gets( tmp );
1250 pclose( gnuplot );
1251 }
1252#endif
1253 return;
1254}
1255
1256//
1257//
1258//
1259
1260unsigned findMaxLocalId( unsigned layerId ) {
1261 /* unsigned maxLocalId = 79;
1262 if(superLayerId == 3)maxLocalId = 127;
1263 else if(superLayerId == 5)maxLocalId = 159;
1264 else if(superLayerId == 7)maxLocalId = 207;
1265 else if(superLayerId == 9)maxLocalId = 255;
1266 return maxLocalId;
1267 */
1268 // changed to BESiii, Liuqg
1269 unsigned maxLocalId = 39;
1270 if ( layerId == 1 ) maxLocalId = 43;
1271 else if ( layerId == 2 ) maxLocalId = 47;
1272 else if ( layerId == 3 ) maxLocalId = 55;
1273 else if ( layerId == 4 ) maxLocalId = 63;
1274 else if ( layerId == 5 ) maxLocalId = 71;
1275 else if ( layerId == 6 || layerId == 7 ) maxLocalId = 79;
1276 else if ( layerId == 20 || layerId == 21 || layerId == 22 || layerId == 23 )
1277 maxLocalId = 159;
1278 else if ( layerId == 24 || layerId == 25 || layerId == 26 || layerId == 27 )
1279 maxLocalId = 175;
1280 else if ( layerId == 28 || layerId == 29 || layerId == 30 || layerId == 31 )
1281 maxLocalId = 207;
1282 else if ( layerId == 32 || layerId == 33 || layerId == 34 || layerId == 35 )
1283 maxLocalId = 239;
1284
1285 return maxLocalId;
1286}
1287
1288unsigned isIsolation( unsigned localId, unsigned maxLocalId, unsigned layerId, int lr,
1289 const AList<TMLink>& allStereoList ) {
1290 unsigned findId;
1291 if ( lr == 1 )
1292 { // R : ox, Dose a wire exist at "x"?
1293 findId = maxLocalId;
1294 if ( localId != 0 ) findId = localId - 1;
1295 }
1296 else if ( lr == -1 )
1297 { // L : xo, Dose a wire exist at "x"?
1298 findId = 0;
1299 if ( localId != maxLocalId ) findId = localId + 1;
1300 }
1301 else { return 1; }
1302 unsigned isolate = 1;
1303 for ( int i = 0; i < allStereoList.length(); ++i )
1304 {
1305 if ( allStereoList[i]->wire()->layerId() == layerId &&
1306 allStereoList[i]->wire()->localId() == findId )
1307 {
1308 isolate = 0;
1309 break;
1310 }
1311 }
1312 return isolate;
1313}
1314
1315void findTwoHits( AList<TMLink>& twoOnLayer, const AList<TMLink>& hitsOnLayer,
1316 const AList<TMLink>& allStereoList ) {
1317 //...finds "two" seq. and isolated hits.
1318 //...and then sets LR for selected hits.
1319 if ( hitsOnLayer.length() == 0 || hitsOnLayer.length() > 3 ) return;
1320 twoOnLayer.removeAll();
1321 if ( hitsOnLayer.length() == 1 )
1322 {
1323 if ( hitsOnLayer[0]->wire()->superLayerId() == 0 ) return; // origin is == 1, Liuqg 060921
1324 unsigned maxLocalId = findMaxLocalId( hitsOnLayer[0]->wire()->layerId() );
1325 unsigned R = isIsolation( hitsOnLayer[0]->wire()->localId(), maxLocalId,
1326 hitsOnLayer[0]->wire()->layerId(), 1, allStereoList );
1327 unsigned L = isIsolation( hitsOnLayer[0]->wire()->localId(), maxLocalId,
1328 hitsOnLayer[0]->wire()->layerId(), -1, allStereoList );
1329 if ( R == 1 && L == 0 )
1330 {
1331 unsigned nextLocalId = hitsOnLayer[0]->wire()->localIdForPlus() + 1;
1332 L = isIsolation( nextLocalId, maxLocalId, hitsOnLayer[0]->wire()->layerId(), -1,
1333 allStereoList );
1334 if ( L == 1 )
1335 { // xuox
1336 hitsOnLayer[0]->leftRight( 1 ); // R
1337 hitsOnLayer[0]->position( hitsOnLayer[0]->arcZ( 1 ) );
1338 twoOnLayer.append( hitsOnLayer[0] );
1339 }
1340 }
1341 else if ( R == 0 && L == 1 )
1342 {
1343 unsigned nextLocalId = hitsOnLayer[0]->wire()->localIdForMinus() + 1;
1344 R = isIsolation( nextLocalId, maxLocalId, hitsOnLayer[0]->wire()->layerId(), 1,
1345 allStereoList );
1346 if ( R == 1 )
1347 { // xoux
1348 hitsOnLayer[0]->leftRight( 0 ); // L
1349 hitsOnLayer[0]->position( hitsOnLayer[0]->arcZ( 0 ) );
1350 twoOnLayer.append( hitsOnLayer[0] );
1351 }
1352 }
1353 }
1354 if ( hitsOnLayer.length() == 2 )
1355 {
1356 if ( hitsOnLayer[0]->wire()->localIdForPlus() + 1 == hitsOnLayer[1]->wire()->localId() )
1357 {
1358 unsigned maxLocalId = findMaxLocalId( hitsOnLayer[0]->wire()->layerId() );
1359 unsigned R = isIsolation( hitsOnLayer[0]->wire()->localId(), maxLocalId,
1360 hitsOnLayer[0]->wire()->layerId(), 1, allStereoList );
1361 unsigned L = isIsolation( hitsOnLayer[1]->wire()->localId(), maxLocalId,
1362 hitsOnLayer[1]->wire()->layerId(), -1, allStereoList );
1363 if ( R == 1 && L == 1 )
1364 { // xoox
1365 hitsOnLayer[0]->leftRight( 1 ); // R
1366 hitsOnLayer[0]->position( hitsOnLayer[0]->arcZ( 1 ) );
1367 hitsOnLayer[1]->leftRight( 0 ); // L
1368 hitsOnLayer[1]->position( hitsOnLayer[1]->arcZ( 0 ) );
1369 twoOnLayer.append( hitsOnLayer[0] );
1370 twoOnLayer.append( hitsOnLayer[1] );
1371 }
1372 }
1373 }
1374 if ( hitsOnLayer.length() == 3 )
1375 {
1376 if ( hitsOnLayer[0]->wire()->localIdForPlus() + 1 == hitsOnLayer[1]->wire()->localId() &&
1377 hitsOnLayer[1]->wire()->localIdForPlus() + 1 != hitsOnLayer[2]->wire()->localId() )
1378 {
1379 unsigned maxLocalId = findMaxLocalId( hitsOnLayer[0]->wire()->layerId() );
1380 unsigned R = isIsolation( hitsOnLayer[0]->wire()->localId(), maxLocalId,
1381 hitsOnLayer[0]->wire()->layerId(), 1, allStereoList );
1382 unsigned L = isIsolation( hitsOnLayer[1]->wire()->localId(), maxLocalId,
1383 hitsOnLayer[1]->wire()->layerId(), -1, allStereoList );
1384 if ( R == 1 && L == 1 )
1385 { // oxoox
1386 hitsOnLayer[0]->leftRight( 1 ); // R
1387 hitsOnLayer[0]->position( hitsOnLayer[0]->arcZ( 1 ) );
1388 hitsOnLayer[1]->leftRight( 0 ); // L
1389 hitsOnLayer[1]->position( hitsOnLayer[1]->arcZ( 0 ) );
1390 twoOnLayer.append( hitsOnLayer[0] );
1391 twoOnLayer.append( hitsOnLayer[1] );
1392 }
1393 }
1394 else if ( hitsOnLayer[0]->wire()->localIdForPlus() + 1 !=
1395 hitsOnLayer[1]->wire()->localId() &&
1396 hitsOnLayer[1]->wire()->localIdForPlus() + 1 ==
1397 hitsOnLayer[2]->wire()->localId() )
1398 {
1399 unsigned maxLocalId = findMaxLocalId( hitsOnLayer[1]->wire()->layerId() );
1400 unsigned R = isIsolation( hitsOnLayer[1]->wire()->localId(), maxLocalId,
1401 hitsOnLayer[1]->wire()->layerId(), 1, allStereoList );
1402 unsigned L = isIsolation( hitsOnLayer[2]->wire()->localId(), maxLocalId,
1403 hitsOnLayer[2]->wire()->layerId(), -1, allStereoList );
1404 if ( R == 1 && L == 1 )
1405 { // xooxo
1406 hitsOnLayer[1]->leftRight( 1 ); // R
1407 hitsOnLayer[1]->position( hitsOnLayer[1]->arcZ( 1 ) );
1408 hitsOnLayer[2]->leftRight( 0 ); // L
1409 hitsOnLayer[2]->position( hitsOnLayer[2]->arcZ( 0 ) );
1410 twoOnLayer.append( hitsOnLayer[1] );
1411 twoOnLayer.append( hitsOnLayer[2] );
1412 }
1413 }
1414 }
1415 /* if(twoOnLayer.length() != 0){
1416 std::cout << "TWO " << twoOnLayer.length() << std::endl;
1417 } */
1418}
1419
1420void setLR( AList<TMLink>& hitsOnLayer, unsigned LR = 0 ) {
1421 // LR = 0 : L
1422 // = 1 : R
1423 for ( unsigned i = 0; i < hitsOnLayer.length(); ++i )
1424 {
1425 if ( LR == 0 )
1426 {
1427 hitsOnLayer[i]->leftRight( 0 ); // L
1428 hitsOnLayer[i]->position( hitsOnLayer[i]->arcZ( 0 ) );
1429 }
1430 else
1431 {
1432 hitsOnLayer[i]->leftRight( 1 ); // R
1433 hitsOnLayer[i]->position( hitsOnLayer[i]->arcZ( 1 ) );
1434 }
1435 }
1436}
1437
1438bool moveLR( AList<TMLink>& hitsOnLayer ) {
1439 unsigned nHits = hitsOnLayer.length();
1440 if ( nHits == 0 ) return false;
1441 // ex) LLLL --> LLLR --> LLRR --> LRRR --> RRRR
1442 if ( hitsOnLayer[nHits - 1]->leftRight() == 1 ) return false; // All R
1443 for ( unsigned i = 0; i < nHits; ++i )
1444 {
1445 if ( hitsOnLayer[i]->leftRight() == 0 )
1446 { // L
1447 hitsOnLayer[i]->leftRight( 1 ); // R
1448 hitsOnLayer[i]->position( hitsOnLayer[i]->arcZ( 1 ) );
1449 return true;
1450 }
1451 }
1452 return false;
1453}
1454
1455void selectGoodWires( const AList<TMLink>& allWires, AList<TMLink>& goodWires ) {
1456 goodWires.removeAll();
1457 for ( int i = 0; i < allWires.length(); ++i )
1458 {
1459 if ( allWires[i]->position().x() != -999. ) { goodWires.append( allWires[i] ); }
1460 }
1461}
1462
1463void TBuilderCurl::fitLine2( const AList<TMLink>& tmpLine, double& min_chi2, double& good_a,
1464 double& good_b, AList<TMLink>& goodLine,
1465 AList<HepPoint3D>& goodPosition, int& overCounter ) const {
1466 AList<TMLink> goodWires;
1467 selectGoodWires( tmpLine, goodWires );
1468 if ( goodWires.length() >= 3 )
1469 fitLine( goodWires, min_chi2, good_a, good_b, goodLine, goodPosition, overCounter );
1470}
1471
1472void calVirtualCircle( const TMLink& hit, const TTrack& track, const int LR,
1473 HepPoint3D& center, double& radius ) {
1474 if ( abs( LR ) != 1 ) return;
1475 double Q = track.charge();
1476 int isOuter = 1;
1477 if ( Q > 0. && LR == 1 ) isOuter = -1; // Inner
1478 else if ( Q < 0. && LR == -1 ) isOuter = -1; // Inner
1479 radius = fabs( track.radius() );
1480 center = track.helix().center();
1481 HepPoint3D wire = hit.wire()->xyPosition();
1482 center.setZ( 0. );
1483 wire.setZ( 0. );
1484 // new center(virtual)
1485 center = wire + ( radius + ( (double)isOuter ) * ( hit.hit()->drift() ) ) *
1486 ( center - wire ).unit();
1487}
1488
1489void moveLR( AList<TMLink>& hits, const AList<TMLink>& hitsOnLayerOrg, const TTrack& track ) {
1490 AList<TMLink> hitsOnLayer = hitsOnLayerOrg;
1491 hitsOnLayer.remove( hits );
1492 if ( hitsOnLayer.length() == 0 ) return;
1493
1494 unsigned nHits = hits.length();
1495 if ( nHits == 0 ) return;
1496
1497 //...finds "ref" from hits.
1498 // ex) LLLL|, LLL|R, LL|RR, L|RRR, |RRRR
1499 int LR = -1; // L
1500 TMLink ref;
1501 if ( hits[nHits - 1]->leftRight() == 1 )
1502 { // All R
1503 LR = 1; // R
1504 ref = *hits[nHits - 1];
1505 }
1506 for ( unsigned i = 0; i < nHits; ++i )
1507 {
1508 if ( hits[i]->leftRight() == 0 )
1509 { // L
1510 ref = *hits[i];
1511 break;
1512 }
1513 }
1514
1515 HepPoint3D center;
1516 double radius;
1517 calVirtualCircle( ref, track, LR, center, radius );
1518
1519 double Q = track.charge();
1520 for ( int i = 0; i < hitsOnLayer.length(); ++i )
1521 {
1522 int isOuter = 1;
1523 if ( ( hitsOnLayer[i]->wire()->xyPosition() - center ).mag() - radius < 0. ) isOuter = -1;
1524 if ( Q > 0. && isOuter == 1 )
1525 {
1526 hitsOnLayer[i]->position( hitsOnLayer[i]->arcZ( 0 ) ); // L
1527 hitsOnLayer[i]->leftRight( 0 ); // L
1528 }
1529 else if ( Q > 0. && isOuter == -1 )
1530 {
1531 hitsOnLayer[i]->position( hitsOnLayer[i]->arcZ( 1 ) ); // R
1532 hitsOnLayer[i]->leftRight( 1 ); // R
1533 }
1534 else if ( Q < 0. && isOuter == 1 )
1535 {
1536 hitsOnLayer[i]->position( hitsOnLayer[i]->arcZ( 1 ) ); // R
1537 hitsOnLayer[i]->leftRight( 1 ); // R
1538 }
1539 else if ( Q < 0. && isOuter == -1 )
1540 {
1541 hitsOnLayer[i]->position( hitsOnLayer[i]->arcZ( 0 ) ); // L
1542 hitsOnLayer[i]->leftRight( 0 ); // L
1543 }
1544 }
1545}
1546
1547void TBuilderCurl::makeLine( TTrack& track, AList<TMLink>& list,
1548 const AList<TMLink>& allStereoList, AList<TMLink>& goodLine,
1549 double& min_chi2, double& good_a, double& good_b,
1550 AList<HepPoint3D>& goodPosition ) const {
1551 if ( list.length() == 0 ) return;
1552
1553 AList<TMLink> layer[24]; // Liuqg, origin is 18.
1554 AList<TMLink> layerOrg[24]; // Liuqg, origin is 18.
1555 for ( unsigned i = 0, size = list.length(); i < size; ++i )
1556 {
1557 layer[list[i]->wire()->layer()->axialStereoLayerId()].append( *list[i] );
1558 layerOrg[list[i]->wire()->layer()->axialStereoLayerId()].append( *list[i] );
1559 }
1560
1561 AList<TMLink> fixedWires[6]; // each superlayer origin is 5,Liuqg 060920
1562 AList<TMLink> nonFixedWires[6]; // each superlayer origin is 5,Liuqg 060920
1563 AList<TMLink> allFixedWires;
1564 for ( unsigned i = 0; i < 24; ++i )
1565 {
1566 if ( layer[i].length() )
1567 {
1568 layer[i].sort( SortByWireId );
1569 sortByLocalId( layer[i] ); // chiisai-jun but kyoukai fukin ha sukoshi kufuu shite iru.
1570 AList<TMLink> tmp;
1571 findTwoHits( tmp, layer[i], allStereoList );
1572 if ( tmp.length() )
1573 {
1574 layer[i].removeAll();
1575 allFixedWires.append( tmp );
1576 // the following parameters have been changed, liuqg 060919
1577 if ( i < 4 ) fixedWires[0].append( tmp );
1578 else if ( i < 8 ) fixedWires[1].append( tmp );
1579 else if ( i < 12 ) fixedWires[2].append( tmp );
1580 else if ( i < 16 ) fixedWires[3].append( tmp );
1581 else if ( i < 20 ) fixedWires[4].append( tmp );
1582 else fixedWires[5].append( tmp );
1583 }
1584 else
1585 {
1586 if ( i < 4 ) nonFixedWires[0].append( layer[i] );
1587 else if ( i < 8 ) nonFixedWires[1].append( layer[i] );
1588 else if ( i < 12 ) nonFixedWires[2].append( layer[i] );
1589 else if ( i < 16 ) nonFixedWires[3].append( layer[i] );
1590 else if ( i < 20 ) nonFixedWires[4].append( layer[i] );
1591 else nonFixedWires[5].append( layer[i] );
1592 }
1593 }
1594 }
1595
1596#if DEBUG_CURL_DUMP
1597 std::cout << "(TBuilderCurl) 1st fixed & non-fixed wires selection:" << std::endl;
1598 std::cout << "(TBuilderCurl) all fixed wires # = " << allFixedWires.length() << std::endl;
1599 std::cout << "(TBuilderCurl) fixed wires # = (" << fixedWires[0].length() << ", "
1600 << fixedWires[1].length() << ", " << fixedWires[2].length() << ", "
1601 << fixedWires[3].length() << ", " << fixedWires[4].length() << ", "
1602 << fixedWires[5].length() << ")" << std::endl;
1603 std::cout << "(TBuilderCurl) non fixed wires # = (" << nonFixedWires[0].length() << ", "
1604 << nonFixedWires[1].length() << ", " << nonFixedWires[2].length() << ", "
1605 << nonFixedWires[3].length() << ", " << nonFixedWires[4].length() << ", "
1606 << nonFixedWires[5].length() << ")" << std::endl;
1607# if 0 /* in detail */
1608 for(unsigned i=0;i<allFixedWires.length();++i)
1609 std::cout << i << ": LocalID/LayerID = " << allFixedWires[i]->wire()->localId()
1610 << "/" << allFixedWires[i]->wire()->layerId()
1611 << ", LR = " << allFixedWires[i]->leftRight() << std::endl;
1612# endif /* in detail */
1613#endif
1614
1615 int createdLine = 0;
1616 int overCounter = 0;
1617 AList<TMLink> goodWires;
1618#if 1 /* fastest finder */
1619 if ( allFixedWires.length() >= 5 )
1620 {
1621 selectGoodWires( allFixedWires, goodWires );
1622 if ( goodWires.length() >= 5 )
1623 {
1624 fitLine( goodWires, min_chi2, good_a, good_b, goodLine, goodPosition, overCounter );
1625 if ( fabs( good_b ) < 10. ) createdLine = 1;
1626 }
1627 }
1628#endif /* fastest finder */
1629#if 1 /* faster finder */
1630 if ( createdLine == 0 )
1631 {
1632 // Q > 0 : Outer = L, Inner = R
1633 // Q < 0 : Outer = R, Inner = L
1634 double Q = track.charge();
1635 unsigned isIncreased = 0;
1636 for ( int sl = 0; sl < 6; ++sl )
1637 { // origin is 5, Liuqg 060919
1638 if ( fixedWires[sl].length() >= 1 && nonFixedWires[sl].length() >= 1 )
1639 {
1640 isIncreased = 1;
1641 unsigned bestNCorrectLR = 0;
1642 double bestR = 1.e9; // add initialize 25-05-15
1643 HepPoint3D bestC;
1644 for ( int i = 0; i < fixedWires[sl].length(); ++i )
1645 {
1646 int LR = fixedWires[sl][i]->leftRight() == 0 ? -1 : 1;
1647 HepPoint3D center;
1648 double radius;
1649 calVirtualCircle( *fixedWires[sl][i], track, LR, center, radius );
1650 unsigned nCorrectLR = 0;
1651 for ( int j = 0; j < fixedWires[sl].length(); ++j )
1652 {
1653 if ( i != j )
1654 {
1655 int tmpIsOuter = 1;
1656 if ( ( fixedWires[sl][j]->wire()->xyPosition() - center ).mag() - radius < 0. )
1657 tmpIsOuter = -1;
1658 if ( Q > 0. && tmpIsOuter == 1 && fixedWires[sl][j]->leftRight() == 0 )
1659 ++nCorrectLR;
1660 else if ( Q > 0. && tmpIsOuter == -1 && fixedWires[sl][j]->leftRight() == 1 )
1661 ++nCorrectLR;
1662 else if ( Q < 0. && tmpIsOuter == 1 && fixedWires[sl][j]->leftRight() == 1 )
1663 ++nCorrectLR;
1664 else if ( Q < 0. && tmpIsOuter == -1 && fixedWires[sl][j]->leftRight() == 0 )
1665 ++nCorrectLR;
1666 }
1667 }
1668 if ( i == 0 || nCorrectLR > bestNCorrectLR )
1669 {
1670 bestNCorrectLR = nCorrectLR;
1671 bestR = radius;
1672 bestC = center;
1673 }
1674 if ( bestNCorrectLR == fixedWires[sl].length() - 1 ) break;
1675 }
1676 for ( int i = 0; i < nonFixedWires[sl].length(); ++i )
1677 {
1678 int isOuter = 1;
1679 if ( ( nonFixedWires[sl][i]->wire()->xyPosition() - bestC ).mag() - bestR < 0. )
1680 isOuter = -1;
1681 if ( Q > 0. && isOuter == 1 )
1682 {
1683 nonFixedWires[sl][i]->position( nonFixedWires[sl][i]->arcZ( 0 ) ); // L
1684 nonFixedWires[sl][i]->leftRight( 0 ); // L
1685 }
1686 else if ( Q > 0. && isOuter == -1 )
1687 {
1688 nonFixedWires[sl][i]->position( nonFixedWires[sl][i]->arcZ( 1 ) ); // R
1689 nonFixedWires[sl][i]->leftRight( 1 ); // R
1690 }
1691 else if ( Q < 0. && isOuter == 1 )
1692 {
1693 nonFixedWires[sl][i]->position( nonFixedWires[sl][i]->arcZ( 1 ) ); // R
1694 nonFixedWires[sl][i]->leftRight( 1 ); // R
1695 }
1696 else if ( Q < 0. && isOuter == -1 )
1697 {
1698 nonFixedWires[sl][i]->position( nonFixedWires[sl][i]->arcZ( 0 ) ); // L
1699 nonFixedWires[sl][i]->leftRight( 0 ); // L
1700 }
1701 }
1702 allFixedWires.append( nonFixedWires[sl] );
1703 fixedWires[sl].append( nonFixedWires[sl] );
1704 nonFixedWires[sl].removeAll();
1705 }
1706 }
1707
1708# if DEBUG_CURL_DUMP
1709 std::cout << "(TBuilderCurl) 2nd fixed & non-fixed wires selection:" << std::endl;
1710 std::cout << "(TBuilderCurl) all fixed wires # = " << allFixedWires.length()
1711 << std::endl;
1712 std::cout << "(TBuilderCurl) fixed wires # = (" << fixedWires[0].length() << ", "
1713 << fixedWires[1].length() << ", " << fixedWires[2].length() << ", "
1714 << fixedWires[3].length() << ", " << fixedWires[4].length() << ", "
1715 << fixedWires[5].length() << ")" << std::endl;
1716 std::cout << "(TBuilderCurl) non fixed wires # = (" << nonFixedWires[0].length() << ", "
1717 << nonFixedWires[1].length() << ", " << nonFixedWires[2].length() << ", "
1718 << nonFixedWires[3].length() << ", " << nonFixedWires[4].length() << ", "
1719 << nonFixedWires[5].length() << ")" << std::endl;
1720# if 0 /* in detail */
1721 for(unsigned i=0;i<allFixedWires.length();++i)
1722 std::cout << i << ": LocalID/LayerID = " << allFixedWires[i]->wire()->localId()
1723 << "/" << allFixedWires[i]->wire()->layerId()
1724 << ", LR = " << allFixedWires[i]->leftRight() << std::endl;
1725# endif /* in detail */
1726# endif
1727
1728 if ( isIncreased == 1 && allFixedWires.length() >= 5 )
1729 {
1730 selectGoodWires( allFixedWires, goodWires );
1731 if ( goodWires.length() >= 5 )
1732 {
1733 fitLine( goodWires, min_chi2, good_a, good_b, goodLine, goodPosition, overCounter );
1734 if ( fabs( good_b ) < 10. ) createdLine = 1;
1735 }
1736 }
1737 }
1738#endif /* faster finder */
1739#if 1 /* slow finder */
1740 // 2000/1/27...very slow but unlike an infinity loop
1741 if ( createdLine == 0 )
1742 {
1743 // nonFixed Wires
1744 int maxNonFixedLayerIndex[6] = { -1, -1, -1, -1, -1, -1 }; // origin is 5, Liuqg 060919
1745 int maxLength[6] = { 0, 0, 0, 0, 0, 0 }; // origin is 5, Liuqg 060919
1746 for ( int l = 0; l < 24; ++l )
1747 {
1748 unsigned sl = 5; // superlayer id, changed by Liuqg
1749 if ( l < 4 ) sl = 0;
1750 else if ( l < 8 ) sl = 1;
1751 else if ( l < 12 ) sl = 2;
1752 else if ( l < 16 ) sl = 3;
1753 else if ( l < 20 ) sl = 4;
1754 layer[l].remove( fixedWires[sl] );
1755 setLR( layer[l] ); // set All L
1756 if ( layer[l].length() && layer[l].length() > maxLength[sl] )
1757 {
1758 maxLength[sl] = layer[l].length();
1759 maxNonFixedLayerIndex[sl] = l;
1760 }
1761 }
1762
1763 unsigned index = 0;
1764 unsigned nonFixedSuperLayerIndex[6] = { 0, 0, 0, 0, 0, 0 }; // origin is 5, Liuqg 060919
1765 unsigned isIncreased = 0;
1766 for ( unsigned i = 0; i < 6; ++i )
1767 { // origin is 5, Liuqg 060919
1768 allFixedWires.append( nonFixedWires[i] );
1769 if ( nonFixedWires[i].length() )
1770 {
1771 isIncreased = 1;
1772 nonFixedSuperLayerIndex[index] = i;
1773 ++index;
1774 }
1775 }
1776
1777 if ( isIncreased )
1778 {
1779 do {
1780 moveLR( layer[maxNonFixedLayerIndex[nonFixedSuperLayerIndex[0]]],
1781 nonFixedWires[nonFixedSuperLayerIndex[0]], track );
1782 if ( index > 1 )
1783 {
1784 setLR( nonFixedWires[nonFixedSuperLayerIndex[1]] );
1785 do {
1786 moveLR( layer[maxNonFixedLayerIndex[nonFixedSuperLayerIndex[1]]],
1787 nonFixedWires[nonFixedSuperLayerIndex[1]], track );
1788 if ( index > 2 )
1789 {
1790 setLR( nonFixedWires[nonFixedSuperLayerIndex[2]] );
1791 do {
1792 moveLR( layer[maxNonFixedLayerIndex[nonFixedSuperLayerIndex[2]]],
1793 nonFixedWires[nonFixedSuperLayerIndex[2]], track );
1794 if ( index > 3 )
1795 {
1796 setLR( nonFixedWires[nonFixedSuperLayerIndex[3]] );
1797 do {
1798 moveLR( layer[maxNonFixedLayerIndex[nonFixedSuperLayerIndex[3]]],
1799 nonFixedWires[nonFixedSuperLayerIndex[3]], track );
1800 if ( index > 4 )
1801 {
1802 setLR( nonFixedWires[nonFixedSuperLayerIndex[4]] );
1803 do {
1804 moveLR( layer[maxNonFixedLayerIndex[nonFixedSuperLayerIndex[4]]],
1805 nonFixedWires[nonFixedSuperLayerIndex[4]], track );
1806 if ( index > 5 )
1807 {
1808 setLR( nonFixedWires[nonFixedSuperLayerIndex[5]] );
1809 do {
1810 moveLR( layer[maxNonFixedLayerIndex[nonFixedSuperLayerIndex[5]]],
1811 nonFixedWires[nonFixedSuperLayerIndex[5]], track );
1812 fitLine2( allFixedWires, min_chi2, good_a, good_b, goodLine,
1813 goodPosition, overCounter );
1814 if ( overCounter > 10000 ) goto kokohe;
1815 } while ( moveLR(
1816 layer[maxNonFixedLayerIndex[nonFixedSuperLayerIndex[5]]] ) );
1817 }
1818 else
1819 {
1820 fitLine2( allFixedWires, min_chi2, good_a, good_b, goodLine,
1821 goodPosition, overCounter );
1822 if ( overCounter > 10000 ) goto kokohe;
1823 }
1824 } while (
1825 moveLR( layer[maxNonFixedLayerIndex[nonFixedSuperLayerIndex[4]]] ) );
1826 }
1827 else
1828 {
1829 fitLine2( allFixedWires, min_chi2, good_a, good_b, goodLine,
1830 goodPosition, overCounter );
1831 if ( overCounter > 10000 ) goto kokohe;
1832 }
1833 } while (
1834 moveLR( layer[maxNonFixedLayerIndex[nonFixedSuperLayerIndex[3]]] ) );
1835 }
1836 else
1837 {
1838 fitLine2( allFixedWires, min_chi2, good_a, good_b, goodLine, goodPosition,
1839 overCounter );
1840 if ( overCounter > 10000 ) goto kokohe;
1841 }
1842 } while ( moveLR( layer[maxNonFixedLayerIndex[nonFixedSuperLayerIndex[2]]] ) );
1843 }
1844 else
1845 {
1846 fitLine2( allFixedWires, min_chi2, good_a, good_b, goodLine, goodPosition,
1847 overCounter );
1848 if ( overCounter > 10000 ) goto kokohe;
1849 }
1850 } while ( moveLR( layer[maxNonFixedLayerIndex[nonFixedSuperLayerIndex[1]]] ) );
1851 }
1852 else
1853 {
1854 fitLine2( allFixedWires, min_chi2, good_a, good_b, goodLine, goodPosition,
1855 overCounter );
1856 if ( overCounter > 10000 ) goto kokohe;
1857 }
1858 } while ( moveLR( layer[maxNonFixedLayerIndex[nonFixedSuperLayerIndex[0]]] ) );
1859 kokohe:;
1860 }
1861 else if ( allFixedWires.length() >= 3 )
1862 {
1863 fitLine2( allFixedWires, min_chi2, good_a, good_b, goodLine, goodPosition, overCounter );
1864 }
1865 }
1866#endif /* slow finder */
1867 for ( unsigned i = 0, size = goodLine.length(); i < size; ++i )
1868 { goodLine[i]->position( *( goodPosition[i] ) ); }
1869#if DEBUG_CURL_DUMP
1870 std::cout << "(TBuilderCurl) make a line from All SuperLayers." << std::endl;
1871 plotArcZ( goodLine, good_a, good_b, 0 );
1872#endif
1873}
1874
1875bool TBuilderCurl::fitWDD( double& xc, double& yc, double& r, AList<TMLink>& list ) const {
1876 if ( list.length() <= 3 ) return false;
1877 Lpav circle;
1878 // MDC
1879 for ( int i = 0; i < list.length(); ++i )
1880 {
1881 circle.add_point( list[i]->wire()->xyPosition().x(), list[i]->wire()->xyPosition().y(),
1882 1.0 );
1883 }
1884 circle.add_point( 0., 0., 1.0 ); // IP Constraint
1885 if ( circle.fit() < 0.0 || circle.kappa() == 0.0 ) return false;
1886 xc = circle.center()[0];
1887 yc = circle.center()[1];
1888 r = circle.radius();
1889 const int maxIte = 2;
1890 for ( int ite = 0; ite < maxIte; ++ite )
1891 {
1892 Lpav circle2;
1893 circle2.clear();
1894 // MDC
1895 for ( int i = 0; i < list.length(); ++i )
1896 {
1897 double R = sqrt( ( list[i]->wire()->xyPosition().x() - xc ) *
1898 ( list[i]->wire()->xyPosition().x() - xc ) +
1899 ( list[i]->wire()->xyPosition().y() - yc ) *
1900 ( list[i]->wire()->xyPosition().y() - yc ) );
1901 if ( R == 0. ) continue;
1902 double U = 1. / R;
1903 double dir = R > r ? -1. : 1.;
1904 double X = xc + ( list[i]->wire()->xyPosition().x() - xc ) * U *
1905 ( R + dir * list[i]->hit()->drift() );
1906 double Y = yc + ( list[i]->wire()->xyPosition().y() - yc ) * U *
1907 ( R + dir * list[i]->hit()->drift() );
1908 circle2.add_point( X, Y, 1.0 );
1909 }
1910 circle2.add_point( 0., 0., 1.0 ); // IP Constraint
1911 if ( circle2.fit() < 0.0 || circle2.kappa() == 0.0 ) return false;
1912 xc = circle2.center()[0];
1913 yc = circle2.center()[1];
1914 r = circle2.radius();
1915 // std::cout << xc << ", " << yc << " : " << r << std::endl;
1916 }
1917 return true;
1918}
1919
1920int TBuilderCurl::stereoHit( double& xc, double& yc, double& r, double& q,
1921 AList<TMLink>& list ) const {
1922 if ( list.length() == 0 ) return -1;
1923
1924 HepPoint3D center( xc, yc, 0. );
1925 HepPoint3D tmp( -999., -999., 0. );
1926 for ( unsigned i = 0, size = list.length(); i < size; ++i )
1927 {
1928 TMDCWireHit& h = *const_cast<TMDCWireHit*>( list[i]->hit() );
1929 HepVector3D X = 0.5 * ( h.wire()->forwardPosition() + h.wire()->backwardPosition() );
1930 HepVector3D x = HepVector3D( X.x(), X.y(), 0. );
1931 HepVector3D w = x - center;
1932 HepVector3D V = h.wire()->direction();
1933 HepVector3D v = HepVector3D( V.x(), V.y(), 0. );
1934 double vmag2 = v.mag2();
1935 double vmag = sqrt( vmag2 );
1936 //...temporary
1937 for ( unsigned j = 0; j < 4; ++j ) list[i]->arcZ( tmp, j );
1938
1939 //...stereo?
1940 if ( vmag == 0. ) continue;
1941
1942 double drift = h.drift( WireHitLeft );
1943 double R[2] = { r + drift, r - drift };
1944 double wv = w.dot( v );
1945 double d2[2];
1946 d2[0] = vmag2 * R[0] * R[0] +
1947 ( wv * wv - vmag2 * w.mag2() ); //...= v^2*(r^2 - w^2*sin()^2)...outer
1948 d2[1] = vmag2 * R[1] * R[1] +
1949 ( wv * wv - vmag2 * w.mag2() ); //...= v^2*(r^2 - w^2*sin()^2)...inner
1950 //...No crossing in R/Phi plane...
1951 if ( d2[0] < 0. && d2[1] < 0. ) continue;
1952
1953 bool ok_inner( true );
1954 bool ok_outer( true );
1955 double d[2] = { -1., -1. };
1956 //...outer
1957 if ( d2[0] >= 0. ) { d[0] = sqrt( d2[0] ); }
1958 else { ok_outer = false; }
1959 if ( d2[1] >= 0. ) { d[1] = sqrt( d2[1] ); }
1960 else { ok_inner = false; }
1961
1962 //...Cal. length and z to crossing points...
1963 double l[2][2];
1964 double z[2][2];
1965 //...outer
1966 if ( ok_outer )
1967 {
1968 l[0][0] = ( -wv + d[0] ) /
1969 vmag2; //...= (-wvcos()+d)/v/v = (-wcos() + (r^2 - w^2*sin()^2)^0.5)/v
1970 l[1][0] = ( -wv - d[0] ) /
1971 vmag2; //...= (-wvcos()+d)/v/v = (-wcos() - (r^2 - w^2*sin()^2)^0.5)/v
1972 z[0][0] = X.z() + l[0][0] * V.z();
1973 z[1][0] = X.z() + l[1][0] * V.z();
1974 }
1975 //...inner
1976 if ( ok_inner )
1977 {
1978 l[0][1] = ( -wv + d[1] ) /
1979 vmag2; //...= (-wvcos()+d)/v/v = (-wcos() + (r^2 - w^2*sin()^2)^0.5)/v
1980 l[1][1] = ( -wv - d[1] ) /
1981 vmag2; //...= (-wvcos()+d)/v/v = (-wcos() - (r^2 - w^2*sin()^2)^0.5)/v
1982 z[0][1] = X.z() + l[0][1] * V.z();
1983 z[1][1] = X.z() + l[1][1] * V.z();
1984 }
1985
1986 //...Cal. xy position of crossing points...
1987 HepVector3D p[2][2];
1988 if ( ok_outer )
1989 {
1990 p[0][0] = x + l[0][0] * v;
1991 p[1][0] = x + l[1][0] * v;
1992 HepVector3D tmp_pc = p[0][0] - center;
1993 HepVector3D pc0 = HepVector3D( tmp_pc.x(), tmp_pc.y(), 0. );
1994 p[0][0] -= drift / pc0.mag() * pc0;
1995 tmp_pc = p[1][0] - center;
1996 HepVector3D pc1 = HepVector3D( tmp_pc.x(), tmp_pc.y(), 0. );
1997 p[1][0] -= drift / pc1.mag() * pc1;
1998 }
1999 if ( ok_inner )
2000 {
2001 p[0][1] = x + l[0][1] * v;
2002 p[1][1] = x + l[1][1] * v;
2003 HepVector3D tmp_pc = p[0][1] - center;
2004 HepVector3D pc0 = HepVector3D( tmp_pc.x(), tmp_pc.y(), 0. );
2005 p[0][1] += drift / pc0.mag() * pc0;
2006 tmp_pc = p[1][1] - center;
2007 HepVector3D pc1 = HepVector3D( tmp_pc.x(), tmp_pc.y(), 0. );
2008 p[1][1] += drift / pc1.mag() * pc1;
2009 }
2010
2011 //...Check r-phi...
2012 bool ok_xy[2][2];
2013 if ( ok_outer )
2014 {
2015 ok_xy[0][0] = true;
2016 ok_xy[1][0] = true;
2017 }
2018 else
2019 {
2020 ok_xy[0][0] = false;
2021 ok_xy[1][0] = false;
2022 }
2023 if ( ok_inner )
2024 {
2025 ok_xy[0][1] = true;
2026 ok_xy[1][1] = true;
2027 }
2028 else
2029 {
2030 ok_xy[0][1] = false;
2031 ok_xy[1][1] = false;
2032 }
2033 if ( ok_outer )
2034 {
2035 if ( q * ( center.x() * p[0][0].y() - center.y() * p[0][0].x() ) < 0. )
2036 ok_xy[0][0] = false;
2037 if ( q * ( center.x() * p[1][0].y() - center.y() * p[1][0].x() ) < 0. )
2038 ok_xy[1][0] = false;
2039 }
2040 if ( ok_inner )
2041 {
2042 if ( q * ( center.x() * p[0][1].y() - center.y() * p[0][1].x() ) < 0. )
2043 ok_xy[0][1] = false;
2044 if ( q * ( center.x() * p[1][1].y() - center.y() * p[1][1].x() ) < 0. )
2045 ok_xy[1][1] = false;
2046 }
2047 if ( !ok_inner && ok_outer && ( !ok_xy[0][0] ) && ( !ok_xy[1][0] ) ) { continue; }
2048 if ( ok_inner && !ok_outer && ( !ok_xy[0][1] ) && ( !ok_xy[1][1] ) ) { continue; }
2049
2050 //...Check z position...
2051 // Liuqg060925, change temporary! these should be changed to bes3, reference is in
2052 // TTrack::szPosition!!!
2053 if ( ok_xy[0][0] )
2054 {
2055 if ( z[0][0] < h.wire()->backwardPosition().z() ||
2056 z[0][0] > h.wire()->forwardPosition().z() )
2057 ok_xy[0][0] = false;
2058 }
2059 if ( ok_xy[1][0] )
2060 {
2061 if ( z[1][0] < h.wire()->backwardPosition().z() ||
2062 z[1][0] > h.wire()->forwardPosition().z() )
2063 ok_xy[1][0] = false;
2064 }
2065 if ( ok_xy[0][1] )
2066 {
2067 if ( z[0][1] < h.wire()->backwardPosition().z() ||
2068 z[0][1] > h.wire()->forwardPosition().z() )
2069 ok_xy[0][1] = false;
2070 }
2071 if ( ok_xy[1][1] )
2072 {
2073 if ( z[1][1] < h.wire()->backwardPosition().z() ||
2074 z[1][1] > h.wire()->forwardPosition().z() )
2075 ok_xy[1][1] = false;
2076 }
2077 if ( ( !ok_xy[0][0] ) && ( !ok_xy[1][0] ) && ( !ok_xy[0][1] ) && ( !ok_xy[1][1] ) )
2078 { continue; }
2079
2080 double cosdPhi, dPhi;
2081 // unsigned index = 0;
2082 unsigned indexL = 0;
2083 unsigned indexR = 0;
2084 // std::cout << "Stereo " << std::endl;
2085 // Q > 0 : Outer = L, Inner = R
2086 // Q < 0 : Outer = R, Inner = L
2087 if ( ok_xy[0][0] )
2088 {
2089 //...cal. arc length...
2090 cosdPhi = -center.dot( ( p[0][0] - center ).unit() ) / center.mag();
2091 if ( fabs( cosdPhi ) < 1.0 ) { dPhi = acos( cosdPhi ); }
2092 else if ( cosdPhi >= 1.0 ) { dPhi = 0.; }
2093 else { dPhi = M_PI; }
2094 // list[i]->arcZ(HepPoint3D(r*dPhi, z[0][0], 0.), index);
2095 // std::cout << r*dPhi << ", " << z[0][0] << std::endl;
2096 //++index;
2097 if ( q > 0 )
2098 {
2099 // std::cout << "Outer L" << std::endl;
2100 if ( indexL == 0 ) list[i]->arcZ( HepPoint3D( r * dPhi, z[0][0], 0. ), 0 );
2101 else list[i]->arcZ( HepPoint3D( r * dPhi, z[0][0], 0. ), 2 );
2102 ++indexL;
2103 }
2104 else
2105 {
2106 // std::cout << "Outer R" << std::endl;
2107 if ( indexR == 0 ) list[i]->arcZ( HepPoint3D( r * dPhi, z[0][0], 0. ), 1 );
2108 else list[i]->arcZ( HepPoint3D( r * dPhi, z[0][0], 0. ), 3 );
2109 ++indexR;
2110 }
2111 }
2112 if ( ok_xy[1][0] )
2113 {
2114 //...cal. arc length...
2115 cosdPhi = -center.dot( ( p[1][0] - center ).unit() ) / center.mag();
2116 if ( fabs( cosdPhi ) < 1.0 ) { dPhi = acos( cosdPhi ); }
2117 else if ( cosdPhi >= 1.0 ) { dPhi = 0.; }
2118 else { dPhi = M_PI; }
2119 // list[i]->arcZ(HepPoint3D(r*dPhi, z[1][0], 0.), index);
2120 // std::cout << r*dPhi << ", " << z[1][0] << std::endl;
2121 //++index;
2122 if ( q > 0 )
2123 {
2124 // std::cout << "Outer L" << std::endl;
2125 if ( indexL == 0 ) list[i]->arcZ( HepPoint3D( r * dPhi, z[1][0], 0. ), 0 );
2126 else list[i]->arcZ( HepPoint3D( r * dPhi, z[1][0], 0. ), 2 );
2127 ++indexL;
2128 }
2129 else
2130 {
2131 // std::cout << "Outer R" << std::endl;
2132 if ( indexR == 0 ) list[i]->arcZ( HepPoint3D( r * dPhi, z[1][0], 0. ), 1 );
2133 else list[i]->arcZ( HepPoint3D( r * dPhi, z[1][0], 0. ), 3 );
2134 ++indexR;
2135 }
2136 }
2137 if ( ok_xy[0][1] )
2138 {
2139 //...cal. arc length...
2140 cosdPhi = -center.dot( ( p[0][1] - center ).unit() ) / center.mag();
2141 if ( fabs( cosdPhi ) < 1.0 ) { dPhi = acos( cosdPhi ); }
2142 else if ( cosdPhi >= 1.0 ) { dPhi = 0.; }
2143 else { dPhi = M_PI; }
2144 // list[i]->arcZ(HepPoint3D(r*dPhi, z[0][1], 0.), index);
2145 // std::cout << r*dPhi << ", " << z[0][1] << std::endl;
2146 //++index;
2147 if ( q > 0 )
2148 {
2149 // std::cout << "Inner R" << std::endl;
2150 if ( indexR == 0 ) list[i]->arcZ( HepPoint3D( r * dPhi, z[0][1], 0. ), 1 );
2151 else list[i]->arcZ( HepPoint3D( r * dPhi, z[0][1], 0. ), 3 );
2152 ++indexR;
2153 }
2154 else
2155 {
2156 // std::cout << "Inner L" << std::endl;
2157 if ( indexL == 0 ) list[i]->arcZ( HepPoint3D( r * dPhi, z[0][1], 0. ), 0 );
2158 else list[i]->arcZ( HepPoint3D( r * dPhi, z[0][1], 0. ), 2 );
2159 ++indexL;
2160 }
2161 }
2162 if ( ok_xy[1][1] )
2163 {
2164 //...cal. arc length...
2165 cosdPhi = -center.dot( ( p[1][1] - center ).unit() ) / center.mag();
2166 if ( fabs( cosdPhi ) < 1.0 ) { dPhi = acos( cosdPhi ); }
2167 else if ( cosdPhi >= 1.0 ) { dPhi = 0.; }
2168 else { dPhi = M_PI; }
2169 // list[i]->arcZ(HepPoint3D(r*dPhi, z[1][1], 0.), index);
2170 // std::cout << r*dPhi << ", " << z[1][1] << std::endl;
2171 //++index;
2172 if ( q > 0 )
2173 {
2174 // std::cout << "Inner R" << std::endl;
2175 if ( indexR == 0 ) list[i]->arcZ( HepPoint3D( r * dPhi, z[1][1], 0. ), 1 );
2176 else list[i]->arcZ( HepPoint3D( r * dPhi, z[1][1], 0. ), 3 );
2177 ++indexR;
2178 }
2179 else
2180 {
2181 // std::cout << "Inner L" << std::endl;
2182 if ( indexL == 0 ) list[i]->arcZ( HepPoint3D( r * dPhi, z[1][1], 0. ), 0 );
2183 else list[i]->arcZ( HepPoint3D( r * dPhi, z[1][1], 0. ), 2 );
2184 ++indexL;
2185 }
2186 }
2187 }
2188 return 0;
2189}
2190
2191void TBuilderCurl::setArcZ( TTrack& track, AList<TMLink>& slayer, AList<TMLink>& alayer0,
2192 AList<TMLink>& alayer1, unsigned ip ) const {
2193 AList<TMLink> tmp = alayer0;
2194 tmp.append( alayer1 );
2195 double xc, yc, r;
2196 if ( fitWDD( xc, yc, r, tmp ) )
2197 {
2198 double q = track.charge();
2199 stereoHit( xc, yc, r, q, slayer );
2200 }
2201}
2202
2203void TBuilderCurl::setArcZ( TTrack& track, AList<TMLink>& slayer, AList<TMLink>& alayer0,
2204 AList<TMLink>& alayer1, AList<TMLink>& alayer2,
2205 unsigned ip ) const {
2206 AList<TMLink> tmp = alayer0;
2207 tmp.append( alayer1 );
2208 tmp.append( alayer2 );
2209 double xc, yc, r;
2210 if ( fitWDD( xc, yc, r, tmp ) )
2211 {
2212 double q = track.charge();
2213 stereoHit( xc, yc, r, q, slayer );
2214 }
2215}
2216
2217void TBuilderCurl::setArcZ( TTrack& track, AList<TMLink>& slayer, AList<TMLink>& alayer0,
2218 AList<TMLink>& alayer1, AList<TMLink>& alayer2,
2219 AList<TMLink>& alayer3, unsigned ip ) const {
2220 AList<TMLink> tmp = alayer0;
2221 tmp.append( alayer1 );
2222 tmp.append( alayer2 );
2223 tmp.append( alayer3 );
2224 double xc, yc, r;
2225 if ( fitWDD( xc, yc, r, tmp ) )
2226 {
2227 double q = track.charge();
2228 stereoHit( xc, yc, r, q, slayer );
2229 }
2230}
2231
2232void TBuilderCurl::setArcZ( TTrack& track, AList<TMLink>& slayer, AList<TMLink>& alayer0,
2233 AList<TMLink>& alayer1, AList<TMLink>& alayer2,
2234 AList<TMLink>& alayer3, AList<TMLink>& alayer4,
2235 unsigned ip ) const {
2236 AList<TMLink> tmp = alayer0;
2237 tmp.append( alayer1 );
2238 tmp.append( alayer2 );
2239 tmp.append( alayer3 );
2240 tmp.append( alayer4 );
2241 double xc, yc, r;
2242 if ( fitWDD( xc, yc, r, tmp ) )
2243 {
2244 double q = track.charge();
2245 stereoHit( xc, yc, r, q, slayer );
2246 }
2247}
2248
2249/*void
2250TBuilderCurl::setArcZ(TTrack &track, AList<TMLink> &slayer,
2251 AList<TMLink> &alayer0,AList<TMLink> &alayer1,
2252 AList<TMLink> &alayer2,AList<TMLink> &alayer3,
2253 AList<TMLink> &alayer4,AList<TMLink> &alayer5,
2254 unsigned ip) const {
2255 AList<TMLink> tmp = alayer0;
2256 tmp.append(alayer1);
2257 tmp.append(alayer2);
2258 tmp.append(alayer3);
2259 tmp.append(alayer4);
2260 tmp.append(alayer5);
2261 double xc,yc,r;
2262 if(fitWDD(xc,yc,r,tmp)){
2263 double q = track.charge();
2264 stereoHit(xc,yc,r,q,slayer);
2265 }
2266}
2267*/
2268
2269IBesMagFieldSvc* TBuilderCurl::getPmgnIMF() const {
2270 if ( nullptr == m_pmgnIMF )
2271 {
2272 StatusCode sc = Gaudi::svcLocator()->service( "MagneticFieldSvc", m_pmgnIMF );
2273 if ( sc.isFailure() )
2274 { std::cout << "Unable to open Magnetic field service" << std::endl; }
2275 }
2276 return m_pmgnIMF;
2277}
2278
2279#undef DEBUG_CURL_DUMP
2280#undef DEBUG_CURL_GNUPLOT
2281#undef DEBUG_CURL_MC
2282
2283// End === Stereo Finder For Curl Tracks : by jtanaka ===
HepVector3D
HepGeom::Vector3D< double > HepVector3D
Definition Analysis/VertexFit/src/BField.cxx:13
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition CosmicGenerator.h:36
n
const Int_t n
Definition DataBase/tau_mode.c:57
data
TTree * data
Definition DataBase/tau_mode.c:5
x
Double_t x[10]
Definition DataBase/tau_mode.c:49
w
double w
Definition EvtPhokharaDef.hh:28
unit
*******INTEGER m_nBinMax INTEGER m_NdiMax !No of bins in histogram for cell exploration division $ !Last vertex $ !Last active cell $ !Last cell in buffer $ !No of sampling when dividing cell $ !No of function total $ !Flag for random ceel for $ !Flag for type of for WtMax $ !Flag which decides whether vertices are included in the sampling $ entire domain is hyp !Maximum effective eevents per saves r n generator level $ !Flag for chat level in !Latex Output unit
Definition FoamA.h:90
m_b
*********DOUBLE PRECISION m_pi INTEGER m_lenwt !max no of aux weights INTEGER m_phmax !maximum photon multiplicity ISR FSR *DOUBLE COMPLEX m_Pauli4 DOUBLE COMPLEX m_AmpBorn DOUBLE COMPLEX m_AmpBoxy DOUBLE COMPLEX m_AmpBorn1 DOUBLE COMPLEX m_AmpBorn2 DOUBLE COMPLEX m_AmpExpo2p DOUBLE COMPLEX m_Rmat DOUBLE COMPLEX m_BoxGZut !DOUBLE COMPLEX m_F1finPair2 !DOUBLE PRECISION m_Vcut DOUBLE PRECISION m_Alfinv DOUBLE PRECISION m_Lorin1 DOUBLE PRECISION m_Lorin2 DOUBLE PRECISION m_b
Definition GPS.h:30
cos
double cos(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:187
Vector
CLHEP::HepVector Vector
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrackUtil/Helix.h:38
WireHitFittingValid
#define WireHitFittingValid
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMDCWireHit.h:29
WireHitLeft
#define WireHitLeft
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMDCWireHit.h:21
WireHitInvalidForFit
#define WireHitInvalidForFit
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMDCWireHit.h:55
SortByWireId
int SortByWireId(const void *a, const void *b)
Sorter.
Definition TMLink.cxx:379
q
****INTEGER imax DOUBLE PRECISION m_pi *DOUBLE PRECISION m_amfin DOUBLE PRECISION m_Chfin DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_sinw2 DOUBLE PRECISION m_GFermi DOUBLE PRECISION m_MfinMin DOUBLE PRECISION m_ta2 INTEGER m_out INTEGER m_KeyFSR INTEGER m_KeyQCD *COMMON c_Semalib $ !copy of input $ !CMS energy $ !beam mass $ !final mass $ !beam charge $ !final charge $ !smallest final mass $ !Z mass $ !Z width $ !EW mixing angle $ !Gmu Fermi $ alphaQED at q
Definition KKsem.h:33
v
**********Class see also m_nmax DOUBLE PRECISION m_amel DOUBLE PRECISION m_x2 DOUBLE PRECISION m_alfinv DOUBLE PRECISION m_Xenph INTEGER m_KeyWtm INTEGER m_idyfs DOUBLE PRECISION m_zini DOUBLE PRECISION m_q2 DOUBLE PRECISION m_Wt_KF DOUBLE PRECISION m_WtCut INTEGER m_KFfin *COMMON c_KarLud $ !Input CMS energy[GeV] $ !CMS energy after beam spread beam strahlung[GeV] $ !Beam energy spread[GeV] $ !z boost due to beam spread $ !electron beam mass *ff pair spectrum $ !minimum v
Definition KarLud.h:35
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition Reconstruction/TrkReco/include/TrkReco/TBuilderCurl.h:36
checkBorder
int checkBorder(AList< TMLink > &layer0, AList< TMLink > &layer1, AList< TMLink > &layer2)
Definition TBuilderCurl.cxx:891
selectGoodWires
void selectGoodWires(const AList< TMLink > &allWires, AList< TMLink > &goodWires)
Definition TBuilderCurl.cxx:1455
makeList
void makeList(AList< TMLink > &layer, AList< TMLink > &list, double q, int border, int checkB, TMLink *layer0)
Definition TBuilderCurl.cxx:953
moveLR
bool moveLR(AList< TMLink > &hitsOnLayer)
Definition TBuilderCurl.cxx:1438
calVirtualCircle
void calVirtualCircle(const TMLink &hit, const TTrack &track, const int LR, HepPoint3D &center, double &radius)
Definition TBuilderCurl.cxx:1472
findMaxLocalId
unsigned findMaxLocalId(unsigned layerId)
Definition TBuilderCurl.cxx:1260
doLineFit
int doLineFit(AList< TMLink > &points, double &m_a, double &m_b, double &chi2, double &nhits, int ipC=1)
Definition TBuilderCurl.cxx:775
isIsolation
unsigned isIsolation(unsigned localId, unsigned maxLocalId, unsigned layerId, int lr, const AList< TMLink > &allStereoList)
Definition TBuilderCurl.cxx:1288
setLR
void setLR(AList< TMLink > &hitsOnLayer, unsigned LR=0)
Definition TBuilderCurl.cxx:1420
offsetBorder
int offsetBorder(TMLink *l)
Definition TBuilderCurl.cxx:943
findTwoHits
void findTwoHits(AList< TMLink > &twoOnLayer, const AList< TMLink > &hitsOnLayer, const AList< TMLink > &allStereoList)
Definition TBuilderCurl.cxx:1315
M_PI
#define M_PI
Definition TConstant.h:4
Helix::center
const HepPoint3D & center(void) const
returns position of helix center(z = 0.);
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrackUtil/Helix.h:197
Helix::a
const HepVector & a(void) const
returns helix parameters.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrackUtil/Helix.h:217
Lpar::center
HepVector center() const
Definition TrackUtil/src/Lpar.cxx:241
Lpav::add_point
void add_point(double x, double y, double w=1)
TBuilder0::name
const std::string & name(void) const
returns name.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TBuilder0.h:126
TBuilder0::TBuilder0
TBuilder0(const std::string &name)
Constructor.
Definition TBuilder0.cxx:31
TBuilder0::fit
virtual int fit(TTrackBase &) const
fits a track using a private fitter.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TBuilder0.h:128
TBuilderCurl::~TBuilderCurl
virtual ~TBuilderCurl()
Destructor.
Definition TBuilderCurl.cxx:73
TBuilderCurl::TBuilderCurl
TBuilderCurl(const std::string &name)
Constructor.
Definition TBuilderCurl.cxx:58
TBuilderCurl::buildStereo
TTrack * buildStereo(TTrack &track, const AList< TMLink > &) const
appends stereo hits to a track.
Definition TBuilderCurl.cxx:168
TBuilderCurl::buildStereoMC
TTrack * buildStereoMC(TTrack &track, const AList< TMLink > &) const
Definition TBuilderCurl.cxx:1047
TBuilderCurl::setParam
void setParam(const TCurlFinderParameter &)
Definition TBuilderCurl.cxx:80
THelixFitter
A class to fit a TTrackBase object to a helix.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/THelixFitter.h:47
TLine0
A class to represent a track in tracking.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TLine0.h:30
TLine0::b
double b(void) const
returns coefficient b.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TLine0.h:135
TLine0::a
double a(void) const
returns coefficient a.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TLine0.h:128
TMDCWireHit::drift
float drift(unsigned) const
returns drift distance.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMDCWireHit.h:218
TMDCWireHit::wire
const TMDCWire *const wire(void) const
returns a pointer to a TMDCWire.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMDCWireHit.h:212
TMDCWire::layerId
unsigned layerId(void) const
returns layer id.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMDCWire.h:207
TMDCWire::direction
const HepVector3D & direction(void) const
returns direction vector of the wire.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMDCWire.h:264
TMDCWire::xyPosition
const HepPoint3D & xyPosition(void) const
returns middle position of a wire. z componet is 0.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMDCWire.h:255
TMDCWire::forwardPosition
const HepPoint3D & forwardPosition(void) const
returns position in forward endplate.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMDCWire.h:242
TMDCWire::backwardPosition
const HepPoint3D & backwardPosition(void) const
returns position in backward endplate.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMDCWire.h:244
TMLine
A class to represent a track in tracking.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMLine.h:40
TMLine::a
double a(void) const
returns coefficient a.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMLine.h:146
TMLine::b
double b(void) const
returns coefficient b.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TMLine.h:153
TRobustLineFitter::fit
virtual int fit(TTrackBase &) const
Definition TRobustLineFitter.cxx:42
TTrackBase::fit
virtual int fit(void)
fits itself by a default fitter. Error was happened if return value is not zero.
Definition TTrackBase.cxx:344
TTrackBase::refine
virtual void refine(AList< TMLink > &list, double maxSigma)
Definition TTrackBase.cxx:164
TTrackBase::append
void append(TMLink &)
appends a TMLink.
Definition TTrackBase.cxx:346
TTrackBase::links
const AList< TMLink > & links(unsigned mask=0) const
Definition TTrackBase.cxx:289
TTrackBase::nLinks
unsigned nLinks(unsigned mask=0) const
returns # of masked TMLinks assigned to this track object.
Definition TTrackBase.cxx:296
TTrack
A class to represent a track in tracking.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TTrack.h:128
TTrack::helix
const Helix & helix(void) const
returns helix parameter.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TTrack.h:434
TTrack::stereoHitForCurl
int stereoHitForCurl(TMLink &link, AList< HepPoint3D > &arcZList) const
TTrack::radius
double radius(void) const
returns signed radius.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TTrack.h:487
TTrack::impact
double impact(void) const
returns signed impact parameter to the origin.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TTrack.h:483
TTrack::charge
double charge(void) const
returns charge.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TTrack.h:455
rb::R
complex_t R(double Q2, double M2, double G, double Mp2, double Mm2)
Definition TUtil.h:22