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