BOSS 7.1.1
BESIII Offline Software System
Loading...
Searching...
No Matches
MdcTrackList Class Reference

#include <MdcTrackList.h>

+ Inheritance diagram for MdcTrackList:

Public Member Functions

 MdcTrackList (const MdcTrackParams &tkPar)
 
 ~MdcTrackList ()
 
int createFromSegs (MdcSegList *segs, const MdcHitMap *, const MdcDetector *, TrkContext &, double bunchTime)
 
int finishCircle (MdcTrack &track, const MdcHitMap *, const MdcDetector *)
 
int finishHelix (MdcTrack &track, const MdcHitMap *, const MdcDetector *)
 
int pickHits (MdcTrack *, const MdcHitMap *, const MdcDetector *, bool pickAmb=true)
 
void dumpSeed (const MdcSeg *seed, bool goodOnly)
 
void dumpCircle (const MdcTrack *)
 
void dumpAxFill (const MdcTrack *)
 
void dumpAxCombine (const MdcTrack *)
 
void dumpD0 (const TrkExchangePar &)
 
void dumpStFill ()
 
void dumpStCombine (const MdcTrack *)
 
void dumpHelix (const MdcTrack *)
 
void dropMultiHotInLayer (const MdcTrack *tk)
 
- Public Member Functions inherited from MdcTrackListBase
 MdcTrackListBase (const MdcTrackParams &tkPar)
 
virtual ~MdcTrackListBase ()
 
int nTrack () const
 
void setPlot (int plotFlag)
 
void newParams (const MdcTrackParams &tkPar)
 
void plot () const
 
void store (RecMdcTrackCol *, RecMdcHitCol *)
 
int arbitrateHits ()
 
void dropMultiHotInLayer (const MdcTrack *tk)
 
void remove (MdcTrack *atrack)
 
void setD0Cut (double d0Cut)
 
void setZ0Cut (double z0Cut)
 
void setPtCut (double ptCut)
 
void setNoInner (bool noInnerFlag)
 

Additional Inherited Members

- Static Public Attributes inherited from MdcTrackListBase
static double m_d0Cut = -999.
 
static double m_z0Cut = -999.
 
static double m_ptCut = -999.
 
- Protected Attributes inherited from MdcTrackListBase
MdcTrackParams tkParam
 
bool m_noInner
 

Detailed Description

Definition at line 31 of file MdcTrackList.h.

Constructor & Destructor Documentation

◆ MdcTrackList()

MdcTrackList::MdcTrackList ( const MdcTrackParams & tkPar)

Definition at line 91 of file MdcTrackList.cxx.

91 :
92 MdcTrackListBase(tkPar) {
93 // *************************************************************************
94 return;
95 }
MdcTrackListBase(const MdcTrackParams &tkPar)

◆ ~MdcTrackList()

MdcTrackList::~MdcTrackList ( )

Definition at line 98 of file MdcTrackList.cxx.

98{}

Member Function Documentation

◆ createFromSegs()

int MdcTrackList::createFromSegs ( MdcSegList * segs,
const MdcHitMap * map,
const MdcDetector * gm,
TrkContext & context,
double bunchTime )
virtual

Implements MdcTrackListBase.

Definition at line 103 of file MdcTrackList.cxx.

105 {
106 // *************************************************************************
107 //Forms tracks out of list of superlayer segments
108 int nTracks = 0;
109
110 m_debug = tkParam.lPrint;//yzhang debug
111
112 // Create empty list of stereo segs (to save time)
113 MdcSegGrouperSt stereoSegs(gm, tkParam.lPrint);
114 stereoSegs.setNoInner(m_noInner);//yzhang add 2016-10-12
115
116 // *** Create r-phi track
117
118 // Create list of axial segments, treated as on tracks from origin
119#ifdef MDCPATREC_TIMETEST
120 TAU_PROFILE_TIMER(timer8,"fill ax seg", "int ()", TAU_DEFAULT);
121 TAU_PROFILE_START(timer8);
122#endif
123 MdcSegGrouperAx origSegs(gm, tkParam.lPrint);
124 origSegs.fillWithSegs(segs);
125 //std::cout << "Plot segs after ax fillWithSegs " << std::endl;//yzhang debug
126 //segs->plot(0);//yzhang debug
127#ifdef MDCPATREC_TIMETEST
128 TAU_PROFILE_STOP(timer8);
129#endif
130 MdcSeg *seed;
131 bool goodOnly = true;
132 bool useBad = (segs->count() < 400); // if true, use non-good seeds
133 //bool useBad = false;
134 segs->resetSeed(gm);
135 MdcTrack *trialTrack = 0;
136
137 while (1) {
138 seed = segs->getSeed(0,goodOnly);
139 if (seed == 0 && goodOnly && useBad) {
140 segs->resetSeed(gm);
141 goodOnly = false;
142 continue;
143 }
144 else if (seed == 0 && (!goodOnly || !useBad)) {
145 break;
146 }
147
148 if (3 == m_debug) dumpSeed(seed, goodOnly);//yzhang debug
149
150 // Pass through first superlayer of MDC required, Field layer No5 = 12.135 cm
151 if ( fabs( ((MdcSegInfoAxialO *) seed->info())->curv()) > 0.0417) continue;//curv cut yzhang
152 delete trialTrack;
153 trialTrack = 0;
154
155 //---------Combine Ax segs--------
156#ifdef MDCPATREC_TIMETEST
157 TAU_PROFILE_TIMER(timer3,"combine ax seg", "int ()", TAU_DEFAULT);
158 TAU_PROFILE_START(timer3);
159#endif
160 int success = origSegs.combineSegs(trialTrack, seed, context, bunchTime,
162#ifdef MDCPATREC_TIMETEST
163 TAU_PROFILE_STOP(timer3);
164#endif
165
166
167 if (3 <= m_debug){
168 cout<<" ***** Ax.combineSegs success? " << success<<"\n";
169 dumpAxCombine(trialTrack);//yzhang debug
170 }
171 if (!success) continue;
172
173
174 //--------Finish circle--------
175#ifdef MDCPATREC_TIMETEST
176 TAU_PROFILE_TIMER(timer4,"circle fitting", "int ()", TAU_DEFAULT);
177 TAU_PROFILE_START(timer4);
178#endif
179 success = finishCircle(*trialTrack, map, gm);
180#ifdef MDCPATREC_TIMETEST
181 TAU_PROFILE_STOP(timer4);
182#endif
183
184 if (3 <= m_debug){
185 cout<<"finishCircle success? " << success<<"\n";
186 dumpCircle(trialTrack);//yzhang debug
187 }
188
189 if (!success) continue;
190 //--------Make sure it really did come from origin--------
191 const TrkFit* tkFit = trialTrack->track().fitResult();
192 assert (tkFit != 0);
193 TrkExchangePar par = tkFit->helix(0.0);
194 //dumpD0(par);
195
196 //------------------d0 cut-------------------
197 // 2010-03-31 , m_d0Cut from 0 to epsilon
198
199 //std::cout<< __FILE__ <<" d0 cut------------"<< fabs(par.d0())<<" d0cut "<< m_d0Cut << std::endl;
200 if ( (m_d0Cut > Constants::epsilon) && (fabs(par.d0()) > m_d0Cut) ){
201 if (tkParam.lPrint>1)cout<<__FILE__<<" Killing track by d0:" <<par.d0()<<">"<<m_d0Cut << endl;
202 continue;
203 }
204
205 //------------------pt cut-------------------
206 if ( (fabs(m_ptCut)>0.) && (fabs(tkFit->pt(0.)) < m_ptCut) ){
207 if (tkParam.lPrint>1)cout<<__FILE__<<" Killing track by pt:" <<tkFit->pt(0.)<<"<"<<m_ptCut << endl;
208 continue;
209 }
210
211 // *** End r-phi track-finding
212 if (3 <= m_debug) std::cout << " *** End r-phi track-finding "<<std::endl;
213
214 //--------Add stereo to taste--------
215#ifdef MDCPATREC_TIMETEST
216 TAU_PROFILE_TIMER(timer5,"fill st seg", "int ()", TAU_DEFAULT);
217 TAU_PROFILE_START(timer5);
218#endif
219 stereoSegs.fillWithSegs(segs, trialTrack);
220#ifdef MDCPATREC_TIMETEST
221 TAU_PROFILE_STOP(timer5);
222#endif
223 if (3 <= m_debug){
224 //dumpStFill();//yzhang debug
225 std::cout<<std::endl<<"----------------------------------------"<< std::endl;
226 std::cout<<" Segment list after St.fillWithSegs "<< std::endl;
227 stereoSegs.dumpSegList();
228 }
229
230#ifdef MDCPATREC_TIMETEST
231 TAU_PROFILE_TIMER(timer6,"combine st seg", "int ()", TAU_DEFAULT);
232 TAU_PROFILE_START(timer6);
233#endif
234 success = stereoSegs.combineSegs(trialTrack, 0, context, bunchTime,
236#ifdef MDCPATREC_TIMETEST
237 TAU_PROFILE_STOP(timer6);
238#endif
239
240 if (3 <= m_debug){
241 cout<<" ***** St.combineSegs success? " << success<<"\n";
242 dumpStCombine(trialTrack);
243 }
244
245
246 //--------Finish 3-d track--------
247 if (success) {
248#ifdef MDCPATREC_TIMETEST
249 TAU_PROFILE_TIMER(timer7,"helix fitting", "int ()", TAU_DEFAULT);
250 TAU_PROFILE_START(timer7);
251#endif
252 success = finishHelix(*trialTrack, map, gm);
253#ifdef MDCPATREC_TIMETEST
254 TAU_PROFILE_STOP(timer7);
255#endif
256 } // end if (success -- st)
257
258 if (3 == m_debug){
259 dumpHelix(trialTrack);
260 }
261 if (!success) continue;
262
263 //------------------d0 cut after helix fitting-------------------
264 //yzhang add 2011-08-01
265 double d0par = trialTrack->track().fitResult()->helix(0.).d0();
266 if ( (m_d0Cut > Constants::epsilon) && (fabs(d0par) > m_d0Cut) ){
267 if (tkParam.lPrint>1) {cout<<__FILE__
268 <<" Killing track by d0 after 3d fit:" <<d0par<<">"<<m_d0Cut << endl;}
269 continue;
270 }
271
272 double z0par = trialTrack->track().fitResult()->helix(0.).z0();
273 if ( (m_z0Cut > Constants::epsilon) && (fabs(z0par) > m_z0Cut) ){
274 if (tkParam.lPrint>1) {cout<<__FILE__
275 <<" Killing track by z0:" <<z0par<<">"<<m_z0Cut << endl;}
276 continue;
277 }
278
279
280 nTracks++;
281 append(trialTrack); // Add to list of Tracks
282
283 trialTrack = 0;
284
285 if (3 == m_debug) std::cout << " ***** End one track-finding *****"<<std::endl;
286 } // end while(1)
287
288 delete trialTrack;
289 return nTracks;
290
291}
static const double epsilon
Definition Constants.h:46
int count() const
MdcSeg * getSeed(int iview, int goodOnly)
void resetSeed(const MdcDetector *gm)
MdcSegInfo * info() const
Definition MdcSeg.h:52
static double m_d0Cut
static double m_ptCut
static double m_z0Cut
MdcTrackParams tkParam
void dumpHelix(const MdcTrack *)
void dumpCircle(const MdcTrack *)
int finishHelix(MdcTrack &track, const MdcHitMap *, const MdcDetector *)
void dumpAxCombine(const MdcTrack *)
void dumpSeed(const MdcSeg *seed, bool goodOnly)
int finishCircle(MdcTrack &track, const MdcHitMap *, const MdcDetector *)
void dumpStCombine(const MdcTrack *)
TrkRecoTrk & track()
Definition MdcTrack.h:27
virtual double pt(double fltL=0.) const =0
double z0() const
double d0() const
virtual TrkExchangePar helix(double fltL) const =0
const TrkFit * fitResult() const
char * append(const char *str, EvtCyclic3::Index i)

◆ dropMultiHotInLayer()

void MdcTrackList::dropMultiHotInLayer ( const MdcTrack * tk)

Definition at line 1152 of file MdcTrackList.cxx.

1152 {
1153 double tdr[43];
1154 double tdr_wire[43];
1155 for(int i=0; i<43; i++){tdr[i]=9999.;}
1156
1157 // make flag
1158 TrkHotList::hot_iterator hotIter= tk->track().hits()->hotList().begin();
1159 while (hotIter!=tk->track().hits()->hotList().end()) {
1160 MdcHitOnTrack* hot = const_cast<MdcHitOnTrack*> (&(*hotIter->mdcHitOnTrack()));
1161
1162 //driftTime(tof,z)
1163 double dt = hot->mdcHit()->driftTime(0.,0.);
1164 int layer = hot->mdcHit()->layernumber();
1165 int wire = hot->mdcHit()->wirenumber();
1166 if (dt < tdr[layer]) {
1167 tdr[layer] = dt;
1168 tdr_wire[layer] = wire;
1169 }
1170 hotIter++;
1171 }
1172
1173 std::cout<<" tdr wire ";
1174 for(int i=0;i<43;i++){
1175 std::cout<<i<<" "<<tdr[i]<<" "<<tdr_wire<<" ";
1176 }
1177 std::cout<<" "<< std::endl;
1178 // inactive multi hit
1179 hotIter= tk->track().hits()->hotList().begin();
1180 while (hotIter!=tk->track().hits()->hotList().end()) {
1181 int layer = hotIter->mdcHitOnTrack()->mdcHit()->layernumber();
1182 int wire = hotIter->mdcHitOnTrack()->mdcHit()->wirenumber();
1183 double dt = hotIter->mdcHitOnTrack()->mdcHit()->driftTime(0.,0.);
1184
1185 if ((tdr[layer] <9998.) && (tdr_wire[layer]!=wire)){
1186 MdcHitOnTrack* hot = const_cast<MdcHitOnTrack*> (&(*hotIter->mdcHitOnTrack()));
1187 hot->setActivity(false);
1188
1189 std::cout<<__FILE__<<" inactive "<< layer<<" "<<wire<<" dt "<<dt << std::endl;
1190 }
1191 hotIter++;
1192 }
1193}
TGraph2DErrors * dt
Definition McCor.cxx:45
virtual const MdcHit * mdcHit() const
unsigned layernumber() const
Definition MdcHit.h:61
unsigned wirenumber() const
Definition MdcHit.h:62
double driftTime(double tof, double z) const
Definition MdcHit.cxx:142
const TrkHotList & hotList() const
void setActivity(bool turnOn)
hot_iterator end() const
Definition TrkHotList.h:45
hot_iterator begin() const
Definition TrkHotList.h:44
TrkHitList * hits()
Definition TrkRecoTrk.h:107

◆ dumpAxCombine()

void MdcTrackList::dumpAxCombine ( const MdcTrack * trialTrack)

Definition at line 1058 of file MdcTrackList.cxx.

1058 {
1059 if (NULL == trialTrack) return;
1060 std::cout<<std::endl<< "-------------------------------------"<<std::endl;
1061 std::cout<<"Track and hitList after AxCombine "<<std::endl;
1062 trialTrack->track().printAll(cout);//yzhang debug
1063 TrkHotList::hot_iterator hotIter= trialTrack->track().hits()->hotList().begin();
1064 while (hotIter!=trialTrack->track().hits()->hotList().end()) {
1065 cout <<"(" <<((MdcHit*)(hotIter->hit()))->layernumber()
1066 <<","<<((MdcHit*)(hotIter->hit()))->wirenumber()
1067 <<":"<<hotIter->isActive()<<") ";
1068 hotIter++;
1069 }
1070 std::cout << std::endl;
1071 std::cout<< "-------------------------------------"<<std::endl;
1072}
#define NULL
virtual void printAll(std::ostream &) const

Referenced by createFromSegs().

◆ dumpAxFill()

void MdcTrackList::dumpAxFill ( const MdcTrack * trialTrack)

Definition at line 1045 of file MdcTrackList.cxx.

1045 {
1046 std::cout << "ax fill: "<<std::endl;
1047 if(!trialTrack){
1048 trialTrack->track().printAll(cout);//yzhang debug
1049 }
1050}

◆ dumpCircle()

void MdcTrackList::dumpCircle ( const MdcTrack * trialTrack)

Definition at line 1074 of file MdcTrackList.cxx.

1074 {
1075 if(NULL == trialTrack) return;
1076 if (!trialTrack->track().fitResult()) return;
1077 /*
1078 double omega,r,pt;
1079 omega =trialTrack->track().fitResult()->helix(0.0).omega();
1080 if (omega == 0) pt = 0;
1081 else pt = (-1) * 1/(omega * 333.576 );
1082 std::cout<<"in MdcTrackList Circle Pt = "<< pt <<std::endl;//yzhang deubg
1083
1084 if (omega == 0) r=0;
1085 else r = 1/ omega;
1086 std::cout<<"in MdcTrackList Circle R = "<< r <<std::endl;//yzhang deubg
1087 */
1088 std::cout<<std::endl<< "-------------------------------------"<<std::endl;
1089 std::cout << "Track and hitList after finishCircle" << std::endl;
1090 trialTrack->track().printAll(cout);
1091 TrkHotList::hot_iterator hotIter= trialTrack->track().hits()->hotList().begin();
1092 while (hotIter!=trialTrack->track().hits()->hotList().end()) {
1093 cout <<"(" <<((MdcHit*)(hotIter->hit()))->layernumber()
1094 <<","<<((MdcHit*)(hotIter->hit()))->wirenumber()
1095 <<":"<<hotIter->isActive()<<") ";
1096 hotIter++;
1097 }
1098 cout <<endl;
1099 std::cout<< "-------------------------------------"<<std::endl;
1100}

Referenced by createFromSegs().

◆ dumpD0()

void MdcTrackList::dumpD0 ( const TrkExchangePar & par)

Definition at line 1102 of file MdcTrackList.cxx.

1102 {
1103 //yzhang hist
1104 //m_hd0->fill(par.d0());
1105 //m_d0 = par.d0();
1106 // m_tuple1->write();//yzhang hist
1107 std::cout <<std::endl<< " Dump d0() " << par.d0()<<"\n";//yzhang debug
1108
1109 //zhangy hist
1110}

◆ dumpHelix()

void MdcTrackList::dumpHelix ( const MdcTrack * trialTrack)

Definition at line 1134 of file MdcTrackList.cxx.

1134 {
1135 std::cout<< std::endl<<"-------------------------------------"<<std::endl;
1136 std::cout<< "Track and hitList after finishHelix " << std::endl;
1137 trialTrack->track().printAll(cout);
1138 TrkHotList::hot_iterator hotIter= trialTrack->track().hits()->hotList().begin();
1139 int tmplay = -1;
1140 while (hotIter!=trialTrack->track().hits()->hotList().end()) {
1141 int layer = ((MdcHit*)(hotIter->hit()))->layernumber();
1142 if( (layer%4) ==0 ) { if( tmplay != layer ) cout<<endl; }
1143 cout <<"(" <<layer <<","<<((MdcHit*)(hotIter->hit()))->wirenumber()
1144 <<":"<<hotIter->isActive() <<") ";
1145 hotIter++;
1146 tmplay = layer;
1147 }
1148 cout <<endl;
1149 std::cout<< "-------------------------------------"<<std::endl;
1150}

Referenced by createFromSegs().

◆ dumpSeed()

void MdcTrackList::dumpSeed ( const MdcSeg * seed,
bool goodOnly )

Definition at line 1052 of file MdcTrackList.cxx.

1052 {
1053 std::cout << std::endl<<"Dump seed segment goodOnly="<<goodOnly<<" ";
1054 seed->plotSegAll();
1055 std::cout<< std::endl;
1056}
void plotSegAll() const
Definition MdcSeg.cxx:159

Referenced by createFromSegs().

◆ dumpStCombine()

void MdcTrackList::dumpStCombine ( const MdcTrack * trialTrack)

Definition at line 1117 of file MdcTrackList.cxx.

1117 {
1118 std::cout<<std::endl<< "-------------------------------------"<<std::endl;
1119 std::cout<<"Track and hitList after StCombine "<<std::endl;
1120 if(trialTrack)trialTrack->track().printAll(cout);
1121 TrkHotList::hot_iterator hotIter= trialTrack->track().hits()->hotList().begin();
1122 int tmplay = -1;
1123 while (hotIter!=trialTrack->track().hits()->hotList().end()) {
1124 int layer = ((MdcHit*)(hotIter->hit()))->layernumber();
1125 if( (layer%4) ==0 ) { if( tmplay != layer ) cout<<endl; }
1126 cout <<"(" <<layer <<","<<((MdcHit*)(hotIter->hit()))->wirenumber()
1127 <<" act:"<<hotIter->isActive() <<" lr:"<<hotIter->ambig() <<") ";
1128 hotIter++;
1129 tmplay=layer;
1130 }
1131 cout <<endl;
1132 std::cout<< "-------------------------------------"<<std::endl;
1133}

Referenced by createFromSegs().

◆ dumpStFill()

void MdcTrackList::dumpStFill ( )

Definition at line 1111 of file MdcTrackList.cxx.

1111 {
1112 std::cout << "Plot segs after st fillWithSegs " << std::endl;
1113 cout <<endl;
1114
1115}

◆ finishCircle()

int MdcTrackList::finishCircle ( MdcTrack & track,
const MdcHitMap * map,
const MdcDetector * gm )

Definition at line 952 of file MdcTrackList.cxx.

953 {
954 //************************************************************************
955 TrkRecoTrk& trk = mdcTrk.track();
956 if(9==tkParam.lPrint){
957 std::cout << " finishCircle "<< std::endl;
958 trk.print(std::cout);
959 }
960
961 const TrkFit* tkFit = trk.fitResult();
962 if(9==tkParam.lPrint){ cout << "Before circle fit, nactive " << tkFit->nActive() << endl;}
963 assert (tkFit != 0);
964 TrkHitList* hitList = trk.hits();
965 assert (hitList != 0);
966 TrkCircleMaker circMaker;
967 circMaker.setFlipAndDrop(trk, false, false); // reset as a precaution
968 //circMaker.setFactor(trk, 1.);//nSigma cut factor//yzhang FIXME 2010-09-13
969
970 TrkErrCode fitResult = hitList->fit();
971 if (fitResult.failure()) {
972 trk.status()->addHistory(TrkErrCode(TrkErrCode::fail,15,"finishCircle"),"MdcTrkRecon");
973 if (tkParam.lPrint > 1) {
974 cout << "First circle fit failed: " << endl;
975 fitResult.print(std::cout);
976 }
977 return 0;
978 }
979 trk.status()->addHistory(TrkErrCode(TrkErrCode::succeed,15,"finishCircle"),"MdcTrkRecon");
980
981 if(9==tkParam.lPrint){ cout << "After circle fit, nactive " << tkFit->nActive() << endl;}
982 double chisqperDOF;
983 int nDOF = tkFit->nActive() - 3;
984 if (nDOF > 3){
985 chisqperDOF = tkFit->chisq() / nDOF;
986 }else{
987 chisqperDOF = tkFit->chisq();
988 }
989
990 if(g_cirTkChi2) { g_cirTkChi2->fill( chisqperDOF ); }//yzhang hist cut
991 int success = (chisqperDOF <= tkParam.maxChisq && tkFit->nActive() >= 3);
992
993 if(9==tkParam.lPrint){
994 std::cout<<__FILE__<<" "<<__LINE__
995 << " success "<<success
996 << " chisqperDOF "<< chisqperDOF<<"<? maxChisq "<< tkParam.maxChisq
997 << " nAct "<<tkFit->nActive()<<">=3 "
998 << std::endl;
999 mdcTrk.track().hots()->printAll(std::cout);
1000 }
1001 bool lcurler(fabs(tkFit->helix(0).omega()) > 3.4);//yzhang FIXME, ??
1002 pickHits(&mdcTrk, map, gm, lcurler);
1003
1004 if(9==tkParam.lPrint){
1005 std::cout<< __FILE__ << " " << __LINE__ << " nHit after pickHit "<<hitList->nHit() <<std::endl;
1006 }
1007 //if(hitList->nHit()<=0) return 0;
1008
1009 circMaker.setFlipAndDrop(trk, true, true);
1010 fitResult = hitList->fit();
1011 if (fitResult.failure()) {
1012 if(9==tkParam.lPrint){
1013 cout << "Second circle fit failed: " << endl;
1014 fitResult.print(std::cout);
1015 }
1016 return 0;
1017 }
1018 if(9==tkParam.lPrint){
1019 cout << "Final fit: " << endl << trk << endl;
1020 }
1021 circMaker.setFlipAndDrop(trk, false, false);
1022
1023 nDOF = tkFit->nActive() - 3;
1024 if (nDOF > 3) {
1025 chisqperDOF = tkFit->chisq() / nDOF;
1026 }
1027 else {
1028 chisqperDOF = tkFit->chisq();
1029 }
1030 if(g_cirTkChi2) { g_cirTkChi2->fill( chisqperDOF ); }//yzhang hist cut
1031 success = (chisqperDOF <= tkParam.maxChisq) && (tkFit->nActive() >= 3);
1032
1033 if(9==tkParam.lPrint){
1034 cout << "chisqperDOF "<<chisqperDOF<<"=?" << tkParam.maxChisq<< endl;
1035 cout << "nActive "<<tkFit->nActive()<<">= 3"<< endl;
1036 }
1037
1038 if(9==tkParam.lPrint){
1039 trk.printAll(cout);
1040 }
1041
1042 return success;
1043}
AIDA::IHistogram1D * g_cirTkChi2
Definition MdcHistItem.h:28
int pickHits(MdcTrack *, const MdcHitMap *, const MdcDetector *, bool pickAmb=true)
virtual double chisq() const =0
void print(std::ostream &ostr) const
int failure() const
Definition TrkErrCode.h:61
double omega() const
virtual void addHistory(const TrkErrCode &status, const char *modulename)
virtual int nActive() const =0
TrkErrCode fit()
unsigned nHit() const
Definition TrkHitList.h:44
virtual void print(std::ostream &) const
const TrkFitStatus * status() const
bool setFlipAndDrop(TrkRecoTrk &, bool allowFlips, bool allowDrops) const

Referenced by createFromSegs().

◆ finishHelix()

int MdcTrackList::finishHelix ( MdcTrack & track,
const MdcHitMap * map,
const MdcDetector * gm )

Definition at line 864 of file MdcTrackList.cxx.

865 {
866 //***********************************************************************
867 int success = 1;
868
869 TrkRecoTrk& trk = mdcTrk.track();
870 TrkErrCode fitResult;
871 TrkHelixMaker helMaker;
872 const TrkFit* tkFit = trk.fitResult();
873 assert (tkFit != 0);
874 TrkHitList* hitList = trk.hits();
875 assert (hitList != 0);
876 TrkExchangePar par = tkFit->helix(0.0);
877
878
879 helMaker.setFlipAndDrop(trk, true, true);
880 fitResult = hitList->fit();
881
882 if (!fitResult.success() && (3 == tkParam.lPrint)) {
883 cout << "Helix fit failure: " << endl;
884 fitResult.print(cout);
885 }
886 helMaker.setFlipAndDrop(trk, false, false);
887
888 if (!fitResult.success()) return 0;
889
890 bool lcurler(fabs(tkFit->helix(0).omega()) > 3.4);//yzhang FIXME, ??
891 pickHits(&mdcTrk, map, gm, lcurler);//yzhang add 2010-05-10
892
893 if(3==tkParam.lPrint) std::cout<< __FILE__ << " " << __LINE__ << " nHit after pickHit "<<hitList->nHit() <<std::endl;
894 if(3==tkParam.lPrint) hitList->hotList().printAll(std::cout);
895 //if(hitList->nHit()<=0) return 0;
896
897 helMaker.setFlipAndDrop(trk, true, true);
898 fitResult = hitList->fit();
899 if (fitResult.failure()) {
900 if(3==tkParam.lPrint){
901 cout << "Second helix fit failed: " << endl;
902 fitResult.print(std::cout);
903 }
904 return 0;
905 }
906 if(3==tkParam.lPrint){ cout << "Final fit: " << endl << trk << endl; }
907 helMaker.setFlipAndDrop(trk, false, false);
908
909 double chisqperDOF = 0.;
910 int nACT = tkFit->nActive();
911 int nDOF = nACT - 5;
912 if (nDOF > 5) {
913 chisqperDOF = tkFit->chisq() / nDOF;
914 } else {
915 chisqperDOF = tkFit->chisq();
916 }
917 if(g_3dTkChi2) { g_3dTkChi2->fill( chisqperDOF ); }//yzhang hist cut
918 if(g_fitNAct) { g_fitNAct->fill(nACT ); }//yzhang hist cut
919
920 int goodMatch = 1;
921 if (fitResult.failure() || chisqperDOF > tkParam.maxChisq || nACT < tkParam.minHits ) {
922 goodMatch = 0;
923 if (3 == tkParam.lPrint) {
924 std::cout<<" goodMatch=0; chi2/dof "<<chisqperDOF <<" >?maxChisq"<<tkParam.maxChisq
925 <<" nAct:"<<nACT <<"<?minHits"<<tkParam.minHits << std::endl;
926 }
927 }
928 //yzhang add
929 ////yzhang FIXME
930 //if (tkParam.lUseQualCuts) {
931 //double tem2 = (float) trk.hits()->nHit() - nACT;
932 //if (tem2 >= tkParam.maxNmissTrack) goodMatch = 0;
933 //if (tem2 / float(trk.hits()->nHit()) > tkParam.maxNmissNorm)
934 //goodMatch = 0;
935 //}
936 //zhangy add
937
938 if (goodMatch) {
939 if(3 <= m_debug){std::cout<<" ***** finishHelix success!"<< std::endl;}
940 trk.fitResult()->positionErr(0.0);
941 } else { // Not a good match
942 success = 0;
943 if(3 <= m_debug){std::cout<<" ***** finishHelix failure!"<< std::endl;}
944 } // end if goodmatch
945
946 return success;
947}
AIDA::IHistogram1D * g_fitNAct
Definition MdcHistItem.h:38
AIDA::IHistogram1D * g_3dTkChi2
Definition MdcHistItem.h:39
virtual BesPointErr positionErr(double fltL) const =0
int success() const
Definition TrkErrCode.h:62
void printAll(std::ostream &o) const

Referenced by createFromSegs().

◆ pickHits()

int MdcTrackList::pickHits ( MdcTrack * trk,
const MdcHitMap * map,
const MdcDetector * gm,
bool pickAmb = true )

!! change to using omega itself

Definition at line 296 of file MdcTrackList.cxx.

297 {
298
299 /***************************************************************************/
300
301 // Selects candidate hits along track;
302 // sorts them into "active" (small residual) and "inactive" (large resid);
303 // throws away hits separated from main group by large gaps. //?? FIXME
304 // pickAm => pick the ambiguity for hits already on track
305 // Return # of active hits found
306
307 bool is2d = trk->track().status()->is2d();
308 if(6==tkParam.lPrint){
309 cout << "Before pickHits";
310 if (is2d) cout<<" 2d:";
311 else cout<<" 3d:";
312 cout<< endl;
313 }
314
315 int nFound = 0;
316 int nCand = 0; // cells tried
317 double rMin, rMax; // min & max search radii for this track
318 double rEnter, rExit; // radii at which track enters, exits layer
319 BesAngle phiEnter, phiExit;//yzhang change
320 int wireLow, wireHigh;
321 double phiLow, phiHigh;
322 int nCell;
323 MdcHit *thisHit;
324 HepAList<MdcHitOnTrack> localHistory; //temporary list of added hits
325 //until bogus hits are chucked
326 double cellWidth;
327 double goodDriftCut; // missing hits with predDrift > goodDriftCut don't count in figuring gaps
328 double aresid = 0.; // = abs(resid)
329 int firstInputHit = -1; //first hit in firstInputLayer/lastInputLayer region
330 int lCurl = 0; // reached curl point
331
332 //****************************************************/
333 const MdcLayer *firstInputLayer = trk->firstLayer();
334 const MdcLayer *lastInputLayer = trk->lastLayer();
335
336 double bunchTime = trk->track().trackT0();
337 const TrkFit* tkFit = trk->track().fitResult();
338 assert (tkFit != 0);
339 const TrkFitStatus* tkStatus = trk->track().status();
340 assert (tkStatus != 0);
341 TrkHitList* hitList = trk->track().hits();
342 assert (hitList != 0);
343 TrkExchangePar par = tkFit->helix(0.0);
344 double d0 = par.d0();
345 double curv = 0.5 * par.omega(); //!!! change to using omega itself
346 double sinPhi0 = sin(par.phi0());
347 double cosPhi0 = cos(par.phi0());
348
349 // *** Set min and max radius for track
350 rMin = gm->firstLayer()->rIn();
351 double absd0 = fabs(d0);
352 if (absd0 > rMin) rMin = absd0 + Constants::epsilon;
353
354 bool willCurl = false;
355 double rCurl = fabs(d0 + 1./curv);
356 rMax = gm->lastLayer()->rOut();
357 //std::cout<< __FILE__ <<" "<< __LINE__<<" rCurl "<< rCurl <<" rMax "<< rMax << std::endl;
358 if (rCurl < rMax) {
359 willCurl = true;
360 rMax = rCurl - Constants::epsilon;
361 }
362 //std::cout<<" willCurl "<< willCurl << std::endl;
363 bool reachedLastInput = false;
364 bool hasCurled = false; // hit found past curl point
365
366 //yzhang add skip layer with hot, 2011-05-03
367 bool isHotOnLayer[43];
369 for(int ii=0; ii<43; ii++) isHotOnLayer[ii]=false;
370 for (TrkHitList::hot_iterator ihit(hitList->begin()); ihit != hitList->end(); ++ihit){
371 isHotOnLayer[ihit->layerNumber()]=true;
372 }
373 }
374
375 // *** Loop through layers in view, looking for the hits
376 // assumes axial inner
377 for (const MdcLayer *layer = gm->firstLayer(); layer != 0;
378 layer = ((!lCurl) ? ( (hasCurled) ? gm->prevLayer(layer) :
379 gm->nextLayer(layer)) : layer) ) {
380
381
382 if (lCurl) {
383 lCurl = 0;
384 hasCurled = true;
385 }
386 if (tkStatus->is2d() && layer->view() != 0) continue;
387
388 if(tkParam.pickSkipExistLayer && isHotOnLayer[layer->layNum()]) continue;//yzhang add 2011-05-03
389
390 //std::cout<< __FILE__ <<" "<< __LINE__ << " lCurl "<< lCurl
391 //<<" hasCurled "<< hasCurled <<" layer "<< layer->layNum() << std::endl;
392 bool lContinue = true;
393 if(tkParam.pickUitlLastLayer) {//yzhang change 2010-09-10
394 if (layer == lastInputLayer) reachedLastInput = true;
395 }
396
397 // *** Find enter and exit points
398 if (hasCurled) {
399 rEnter = layer->rOut();
400 if (rEnter < rMin) {
401 //std::cout<< __FILE__ <<" "<< __LINE__
402 //<<" rEnter<rMin "<<rEnter<<" "<<rMin<< std::endl;
403 break;
404 }
405 rExit = layer->rIn();
406 if (rExit < rMin) rExit = rMin;
407 if (rEnter > rMax) rEnter = rMax;
408
409 //std::cout<< __FILE__ <<" hasCurled: rEnter "<< rEnter
410 //<<" rExit " << rExit << " rMin "<<rMin<<" rMax "<<rMax<< std::endl;
411 } else {
412 rEnter = layer->rIn();
413 rExit = layer->rOut();
414 //std::cout<< __FILE__ <<" NOT hasCurled: rEnter "<< rEnter
415 //<<" rExit " << rExit << " rMin "<<rMin<<" rMax "<<rMax<< std::endl;
416 if (rExit < rMin) continue;
417
418 if (willCurl) {
419 if (rEnter > rMax) {
420 //std::cout<< __FILE__<<" Reached curl point before entering layer"<< std::endl;
421 // Reached curl point before entering layer
422 hasCurled = 1;
423 continue;
424 }
425 if (rExit > rMax) {
426 lCurl = 1;
427 rExit = rMax;
428 }
429 } else { // not a potential curler
430 //std::cout<< __FILE__ <<" "<< __LINE__ <<" not a potential curler"<< std::endl;
431 if (rEnter > rMax) {
432 //std::cout<< __FILE__ <<" rEnter> rMax "<< rEnter << std::endl;
433 break;
434 }
435 if (rExit > rMax) rExit = rMax;
436 }
437 } // end if curled already
438
439 nCell = layer->nWires();
440 cellWidth = Constants::twoPi * layer->rMid() / nCell;//??
441 // don't count outer xmm of cell
442 goodDriftCut = 0.5 * cellWidth * M_SQRT2 + tkParam.pickHitMargin;
443 //goodDriftCut = 0.5 * cellWidth - tkParam.pickHitMargin;//yzhang change 2012-08-17
444 double deltaPhiCellWidth = 0.5 * (cellWidth * M_SQRT2)/layer->rMid();
445
446 //**** Find entrance and exit phi's
447 BesAngle phiTemp(0.0);
448 int ierror = trk->projectToR(rEnter, phiTemp, hasCurled);
449 phiEnter = phiTemp.posRad();
450 if (ierror != 0) {
451 if(6==tkParam.lPrint) std::cout<< " ErrMsg(warning) "
452 << "Error in MdcTrackList::pickHits projection, ierror " <<ierror<< std::endl;
453 continue;//yzhang 2011-04-14
454 }
455 ierror = trk->projectToR(rExit, phiTemp, hasCurled);
456 phiExit = phiTemp.posRad();
457 if (ierror != 0) {
458 if(6==tkParam.lPrint) std::cout<< " ErrMsg(warning) "
459 << "Error in MdcTrackList::pickHits projection, ierror "<<ierror << std::endl;
460 continue;//yzhang 2011-04-14
461 }
462
463
464 if(6==tkParam.lPrint){
465 std::cout<< endl<<"--pickHit of layer "<<layer->layNum()<<"--"<<std::endl;
466 std::cout<< " track phiEnter "<< phiEnter.deg()<<" phiExit "<<phiExit.deg()<<" degree"<< std::endl;
467 std::cout<< " cell width "<< 360./layer->nWires()<<" dPhiz "<<layer->dPhiz()*Constants::radToDegrees <<" deltaPhiCellWidth "<<0.5 * (cellWidth * M_SQRT2)/layer->rMid() * Constants::radToDegrees<<std::endl;
468 std::cout<< " maxInactiveResid "<< tkParam.maxInactiveResid <<" pickHitPhiFactor "<<tkParam.pickHitPhiFactor<< std::endl;
469 std::cout<< " goodDriftCut "<< goodDriftCut <<"=sqrt(2)*0.5*"<<cellWidth<<"+"<<tkParam.pickHitMargin<< std::endl;
470 }
471
472 double Bz = trk->track().bField().bFieldZ();
473 //std::cout<< " facotr "<<tkParam.pickHitPhiFactor<< " dPhiz "<<layer->dPhiz()
474 //<< " factor*dPhiz "<<layer->dPhiz()*tkParam.pickHitPhiFactor<< std::endl;
475 double t_phiHit = -999.;
476 if (curv*Bz <= 0.0) {
477 //positive track in minus Bz
478 phiLow = phiEnter;
479 phiHigh = phiExit;
480 // Allow some slop in phi
481 phiLow -= tkParam.maxInactiveResid / rEnter;
482 phiHigh += tkParam.maxInactiveResid / rExit;
483 } else {
484 phiLow = phiExit;
485 phiHigh = phiEnter;
486 phiLow -= tkParam.maxInactiveResid / rExit;
487 phiHigh += tkParam.maxInactiveResid / rEnter;
488 }
489 //yzhang fix cross x axis bug 2011-04-10
490 if((phiHigh>0 && phiLow<0)){
491 phiLow += Constants::twoPi;
492 }else if((phiHigh<0 && phiLow>0)){
493 phiHigh += Constants::twoPi;
494 }
495
496 double lowFloat = nCell /Constants::twoPi * (phiLow - layer->phiOffset()) + 0.5;
497 double highFloat = nCell /Constants::twoPi * (phiHigh - layer->phiOffset()) + 0.5;
498 wireLow = (lowFloat >= 0.0) ? int(lowFloat) : int(floor(lowFloat));
499 wireHigh = (highFloat >= 0.0) ? int(highFloat) : int(floor(highFloat));
500
501 if(6==tkParam.lPrint){
502 std::cout << " wireLow "<<wireLow << " wireHigh "<<wireHigh <<" phiLow "<<phiLow*180./Constants::pi << " phiHigh "<<phiHigh*180./Constants::pi << std::endl;
503 }
504 // *** Loop through the predicted hit wires
505 int tempDiff = 0;
506 if(g_pickHitWire) {
507 int tempN = Constants::maxCell[layer->layNum()];
508 if(wireLow>tempN/2. && wireHigh<tempN/2.){
509 g_pickHitWire->fill(wireHigh+tempN - wireLow);
510 tempDiff = wireHigh+tempN - wireLow +1;
511 }else{
512 g_pickHitWire->fill(wireHigh - wireLow);
513 tempDiff = wireHigh - wireLow +1;
514 }
515 }//yzhang hist cut
516
517 if(wireLow>wireHigh) wireLow -= nCell;//yzhang 2011-05-16
518 long t_iHit = 0;
519 for (int thisWire = wireLow; thisWire <= wireHigh; thisWire++) {
520 int corrWire = mdcWrapWire(thisWire, nCell);
521 thisHit = map->hitWire(layer->layNum(), corrWire);
522
523 if(6==tkParam.lPrint){
524 if(thisHit != 0) {cout<<endl<<"test hit "; thisHit->print(std::cout);}
525 else std::cout << endl<<"test hit ("<<layer->layNum()<<","<<corrWire<<")"<< std::endl;
526 }
527
528 double tof = 0.;
529 double z = 0.;
530 double driftDist = 0.;
531 // Calculate expected drift distance
532 double delx, dely;
533 double resid = 0., predDrift = 0.;
534 int ambig = 0;
535 const MdcHitOnTrack *alink = 0;
536 double mcTkId = -9999; //yzhang for tuple 2011-06-28
537 if (thisHit == 0 ) {
538 delx = -d0 * sinPhi0 - layer->xWire(corrWire);
539 dely = d0 * cosPhi0 - layer->yWire(corrWire);
540 predDrift = cosPhi0 * dely - sinPhi0 * delx +
541 curv * (delx * delx + dely * dely);
542 if(6==tkParam.lPrint) cout<<"No hit. predDrift="<<predDrift<<endl;
543 continue;
544 } else {
545 if(m_tuplePick) mcTkId = thisHit->digi()->getTrackIndex();
546 // look for existing hit
547 if(thisHit->getHitOnTrack(&(trk->track())) ==0){
548 alink = 0;//yzhang temp
549 }else{
550 alink = thisHit->getHitOnTrack(&(trk->track()))->mdcHitOnTrack();//yzhang temp
551 if(6==tkParam.lPrint) { cout << " existing hot; " ;}
552 }
553
554 if (alink == 0 || pickAm) {
555 if ((!tkStatus->is2d()) && layer->view() != 0){
556 double rc = 9999.;
557 if(fabs(par.omega())>Constants::epsilon) rc = 1./fabs(par.omega());
558 double rw = layer->rMid();
559 double alpha = acos(1 - rw*rw/(2*rc*rc));
560 double fltLen = rw;
561 if(fabs(1 - rw*rw/(2*rc*rc))<1) fltLen = rc * alpha;
562 z = par.z0() + fltLen* par.tanDip();
563 tof = fltLen / Constants::c;
564 double x = layer->getWire(corrWire)->xWireDC(z);
565 double y = layer->getWire(corrWire)->yWireDC(z);
566 delx = -d0 * sinPhi0 - x;
567 dely = d0 * cosPhi0 - y;
568 if(m_tuplePick) t_phiHit = thisHit->phi(z);
569 }else{
570 delx = -d0 * sinPhi0 - thisHit->x();
571 dely = d0 * cosPhi0 - thisHit->y();
572 if(m_tuplePick) t_phiHit = thisHit->phi();
573 }
574 predDrift = cosPhi0 * dely - sinPhi0 * delx +
575 curv * (delx * delx + dely * dely);
576
577 // predDrift = predDrift * (1. - curv() * predDrift);
578 ambig = (predDrift >= 0) ? 1 : -1;
579 if (hasCurled) ambig = -ambig;
580 double entranceAngle=0.;
581 driftDist = thisHit->driftDist(tof+bunchTime,ambig,entranceAngle,0.,z);
582 if (alink == 0) {
583 // FIXME: is this ambig here VVVVV OK for incoming tracks?
584 resid = ambig * fabs(driftDist) - predDrift;
585 aresid = fabs(resid);
586 //if (aresid > tkParam.maxInactiveResid ) ambig = 0;//yzhang delete 2012-10-09
587 }
588 } else {
589 ambig = alink->ambig();
590 if(6==tkParam.lPrint) cout << " pick ambig for hot"<<endl;
591 if(m_tuplePick) t_phiHit = par.phi0()+par.omega()*alink->fltLen();
592 }
593 }
594
595 //yzhang 2012-08-20 , guess phi of this hit on z
596 double zGuess = par.z0() + layer->rMid() * par.tanDip();
597 double phiDCz = layer->getWire(corrWire)->phiDC(zGuess);
598 BesAngle phiDCzMax(phiDCz + deltaPhiCellWidth);
599 BesAngle phiDCzMin(phiDCz - deltaPhiCellWidth);
600
601 if(m_tuplePick && alink==0) {
602 double sigma = 999.;
603 if (thisHit != 0 &&alink==0) {
604 double entranceAngle = 0.;
605 sigma = thisHit->sigma(driftDist, ambig, entranceAngle, atan(par.tanDip()), z);
606 }
607 m_pickPredDrift[t_iHit] = predDrift;
608 m_pickDrift[t_iHit] = driftDist;
609 m_pickDriftTruth[t_iHit] = haveDigiDrift[thisHit->layernumber()][thisHit->wirenumber()];
610 if(curv*Bz<=0.){
611 //positive track under minus Bz
612 if((phiDCzMin-phiExit>0) || (phiDCzMax-phiEnter<0)) m_pickPhizOk[t_iHit] = 0;
613 else m_pickPhizOk[t_iHit] = 1;
614 }else{
615 if((phiDCzMin-phiEnter>0) || (phiDCzMax-phiExit<0)) m_pickPhizOk[t_iHit] = 0;
616 else m_pickPhizOk[t_iHit] = 1;
617 }
618 m_pickZ[t_iHit] = z;
619 m_pickWire[t_iHit]=thisHit->wirenumber();
620 m_pickResid[t_iHit] = aresid;
621 if(abs(sigma)>0.000001) m_pickSigma[t_iHit] = sigma;
622 else m_pickSigma[t_iHit] = 999.;
623 double t_phiLowCut=-999.;
624 double t_phiHighCut= -999.;
625 if(t_phiHit > -998.){
626 t_phiLowCut = (phiEnter-t_phiHit)*rEnter;
627 t_phiHighCut = (phiExit-t_phiHit)*rExit;
628 }
629 m_pickPhiLowCut[t_iHit] = t_phiLowCut;
630 m_pickPhiHighCut[t_iHit] = t_phiHighCut;
631 m_pickDriftCut[t_iHit] = goodDriftCut;
632 m_pickMcTk[t_iHit] = mcTkId;
633 m_pickPt[t_iHit] = tkFit->pt();
634 m_pickCurv[t_iHit] = curv;
635 t_iHit++;
636 }
637
638 if(curv*Bz<=0.){
639 //positive track
640 if((phiDCzMin-phiExit>0) || (phiDCzMax-phiEnter<0)) {
641 if(6==tkParam.lPrint){ std::cout<<" CUT by phiDCz not in phi En Ex range, curv>=0"<<std::endl; }
642 continue;
643 }
644 }else{
645 //negtive track
646 if((phiDCzMin-phiEnter>0) || (phiDCzMax-phiExit<0)) {
647 if(6==tkParam.lPrint){ std::cout<<" CUT by phiDCz not in phi En Ex range, curv<0"<<std::endl; }
648 continue;
649 }
650 }
651
652 // Cases
653 // (0) pred drift dist > goodDriftCut, drop Hit
654 if (ambig != 0 && fabs(predDrift) > goodDriftCut){//yzhang add 2012-08-17
655 if(6==tkParam.lPrint){cout<<" predDrift "<<predDrift<<">goodDriftCut "<<goodDriftCut<<endl;}
656 continue;
657 }
658
659 // (1) No good hit, but track near cell-edge => do nothing and continue
660 if (ambig == 0 && fabs(predDrift) > goodDriftCut){//yzhang 2009-11-03 add 3factor, 2011-05-30 from 3 to 2,2012-08-17 from 2 to 1
661 if(6==tkParam.lPrint){
662 cout<<" ambig==0 && |predDirft| "<<fabs(predDrift) <<"> goodDriftCut "<< goodDriftCut<<endl;
663 cout<<" No good hit, but track near cell-edge " << endl;
664 }
665 continue;
666 }
667
668
669 // (2) Hit found:
670 if (ambig != 0) {
671 //yzhang changed 2009-10-19
672 //if resid> maxActiveSimga * sigma do not add hits to track
673 double entranceAngle = 0.;
674 double sigma = thisHit->sigma(driftDist, ambig, entranceAngle, atan(par.tanDip()), z);
675 double factor = tkParam.pickHitFactor;
676 //yzhang delete 2012-10-09
677 //if(!is2d fabs(par.tanDip())<2) factor = (2-par.tanDip())/2 * factor;
678 double residCut = tkParam.maxActiveSigma * factor * sigma;//yzhang 2012-08-21
679 //if(6==tkParam.lPrint){
680 //std::cout<< "aresid "<<aresid<<" residCut "<<residCut<<" sigma "<<sigma <<" tanDip "<< par.tanDip() <<" factor "<<factor <<" tkParam.maxActiveSigma "<<tkParam.maxActiveSigma<<std::endl;
681 //}
682
683 if (alink == 0 && (aresid <= residCut) ) {
684 if(6==tkParam.lPrint){
685 cout << " (2) New hit found " << endl;//yzhang debug
686 }
687 //yzhang 2012-08-17 delete
688 int isActive = 1;
689 //if (aresid > residCut) isActive = 0;
690 //if(6==tkParam.lPrint) {
691 // if (aresid > residCut)
692 // std::cout << "notACT, resid "<<aresid<<" >residCut " << residCut<< std::endl;
693 //}
694 MdcRecoHitOnTrack tempHot(*thisHit, ambig, bunchTime);
695 tempHot.setActivity(isActive);
696 // Don't add active hits if they're in use.
697 if (thisHit->usedHit()){
698 tempHot.setUsability(false);
699 if(6==tkParam.lPrint) std::cout<< " thisHit used, setUsability false " << std::endl;
700 }
701 // very crude starting point for flight length !!!! improve
702 double flt = layer->rMid();
703 if (hasCurled) flt = Constants::twoPi * rCurl - layer->rMid();
704 flt += 0.000001 * (thisHit->x() + thisHit->y());
705 tempHot.setFltLen(flt);
706 if(6==tkParam.lPrint) {
707 std::cout<< " aresid "<<aresid<<"<=residCut " <<residCut<<" nSig "<<aresid/sigma<< " nSigCut "<<(tkParam.maxActiveSigma*factor)<<" factor "<<factor<<" maxActiveSigma "<<tkParam.maxActiveSigma<<" tanDip "<<par.tanDip()<<std::endl;
708 std::cout<< " Append Hot " << std::endl;
709 }
710 alink = (MdcHitOnTrack*) hitList->appendHot(&tempHot);
711 }else{
712 if(6==tkParam.lPrint){
713 if(alink!=0){
714 std::cout<< "Exist hot found"<<std::endl;
715 }else if(aresid > residCut){
716 thisHit->print(std::cout);
717 std::cout<< " Drop hit. aresid "<<aresid<<">residCut " <<residCut<<" nSig "<<aresid/sigma<< " nSigCut "<<(tkParam.maxActiveSigma*factor)<<" factor "<<factor<<" maxActiveSigma "<<tkParam.maxActiveSigma<<" tanDip "<<par.tanDip()<<std::endl;
718 }
719 }
720 }
721 if (!localHistory.member(const_cast<MdcHitOnTrack*>(alink))) {
722 localHistory.append(const_cast<MdcHitOnTrack*>(alink));
723 if (hasCurled) trk->setHasCurled();
724 nFound++;
725 if(6==tkParam.lPrint) std::cout<<" nFound="<<nFound<<" nCand "<<nCand<<std::endl;
726 if (layer == firstInputLayer && firstInputHit < 0) {
727 firstInputHit = nCand;
728 }
729 } else {
730 if(6==tkParam.lPrint) std::cout << "ErrMsg(warning) "<< "would have inserted identical HOT "
731 "twice in history buffer" << std::endl;
732 }
733 }
734
735 // (3) No hit found; if beyond known good region, should we continue?
736 else if (ambig == 0 && reachedLastInput) {
737 if(6==tkParam.lPrint) std::cout << "ambig==0 "<<std::endl;
738 lContinue = false;
739 int nSuccess = 0;
740 int last = 8;
741 if (nCand < 8) last = nCand;
742 for (int i = 0; i < last; i++) {
743 int j = nCand - 1 - i;
744 if (localHistory[j] != 0) {
745 //std::cout<< __FILE__ <<" localHistory["<<j<<"]!=0 nSuccess++ "<< std::endl;
746 nSuccess++;
747 }
748 if (i == 2 && nSuccess >= 2) lContinue = true;
749 if (i == 4 && nSuccess >= 3) lContinue = true;
750 if (i == 7 && nSuccess >= 5) lContinue = true;
751 if(6==tkParam.lPrint) cout <<i<< " (3) No hit found; if beyond known good region " << endl;//yzhang debug
752 if (lContinue) {
753 if(6==tkParam.lPrint) std::cout<<" pickHits break by lContinue. i="<<i<<" nSuccess="<<nSuccess<< std::endl;
754 break;
755 }
756 }
757
758 if(6==tkParam.lPrint) cout << " (3) No hit found; if beyond known good region " << endl;//yzhang debug
759 // if lContinue == false => there's been a gap, so quit
760 if (!lContinue) {
761 //std::cout<< __FILE__ <<" "<< __LINE__ <<" break by !lContinue "<< std::endl;
762 break;
763 }
764 localHistory.append(0);
765 }
766
767 nCand++;
768 // Update last layer:
769 if (ambig != 0 && reachedLastInput) {
770 if (trk->hasCurled() == 0) {
771 if (thisHit->layer()->rEnd() > trk->lastLayer()->rEnd()) {
772 trk->setLastLayer(thisHit->layer());
773 }
774 }
775 else {
776 if (thisHit->layer()->rEnd() < trk->lastLayer()->rEnd()) {
777 trk->setLastLayer(thisHit->layer());
778 }
779 }
780 }
781
782 } // end loop over hit wires
783 if(t_iHit>0 && m_tuplePick) {
784 m_pickNCell = tempDiff;
785 m_pickNCellWithDigi = t_iHit;
786 m_pickLayer = layer->layNum();
787 m_pickIs2D = is2d;
788 m_tuplePick->write();
789 }
790 if (!lContinue && reachedLastInput) {
791 //std::cout<< __FILE__ <<" break by !lContinue && reachedLastInput "<< std::endl;
792 break;
793 }
794 } // end loop over layers
795
796 // Now look back at hits in early layer and see if there are any to be thrown away there.
797 bool lContinue = true;
798 for (int ihit = firstInputHit; ihit >= 0; ihit--) {
799 if (localHistory[ihit] != 0) {
800 if (lContinue) {
801 // Update firstLayer
802 const MdcHitOnTrack *mdcHit = localHistory[ihit]->mdcHitOnTrack();
803 if (mdcHit != 0) {
804 if (mdcHit->layer()->rEnd() < trk->firstLayer()->rEnd()) {
805 trk->setFirstLayer(mdcHit->layer());
806 }
807 }
808 continue; // no gap yet
809 } else {
810 // gap found; delete hits
811 if(6==tkParam.lPrint){
812 std::cout << " gap found; delete hits. ";
813 }
814 if (!localHistory[ihit]->isUsable()) {
815 if(6==tkParam.lPrint){
816 cout << "about to delete hit for unusable HOT:" << endl;
817 localHistory[ihit]->print(std::cout);
818 }
819 hitList->removeHit(localHistory[ihit]->hit());
820 }
821 if(6==tkParam.lPrint){
822 cout << " current contents of localHistory: "
823 <<localHistory.length()<<"hot" << endl;
824 //for (int i=0; i<localHistory.length();++i) {
825 //cout << " hit @ " << localHistory[i]->hit() << "; hot @ " << localHistory[i] << endl;
826 //}
827 }
828 nFound--;
829 }
830 }
831 else if (localHistory[ihit] == 0) {
832 if(6==tkParam.lPrint){ cout << " localHistory= 0 " << endl; }
833 int nSuccess = 0;
834 lContinue = false;
835 int last = 8;
836 if (nCand < 8) last = nCand;
837 for (int i = 0; i < last; i++) {
838 int j = ihit + 1 + i;
839 if (localHistory[j] != 0) nSuccess++;
840 if (i == 2 && nSuccess >= 2) lContinue = true;
841 if (i == 4 && nSuccess >= 3) lContinue = true;
842 // if (i == 7 && nSuccess >= 5) lContinue = true;
843 if (lContinue) break;
844 }
845 }
846 }
847 if(6==tkParam.lPrint){
848 std::cout<< "nFound="<<nFound <<" < "<< tkParam.pickHitFract*nCand<<" pickHitFract? "<< tkParam.pickHitFract<<"*"<<nCand << std::endl;
849 }
850 if (nFound < tkParam.pickHitFract * nCand) nFound = 0;//yzhang delete 2010-05-10 use pickHitFract=0.
851 if(6==tkParam.lPrint){ cout << " localHistory= 0 " << endl; }
852
853 if(6==tkParam.lPrint || 3==tkParam.lPrint){
854 cout << "After pickHits found " << nFound <<" hit(s)"<< endl;
855 hitList->hotList().print(std::cout);
856 std::cout<< std::endl;
857 }
858
859 return nFound;
860}
double sin(const BesAngle a)
Definition BesAngle.h:210
double cos(const BesAngle a)
Definition BesAngle.h:213
TTree * sigma
Double_t x[10]
const double alpha
NTuple::Array< long > m_pickIs2D
NTuple::Array< double > m_pickPt
NTuple::Item< long > m_pickLayer
NTuple::Array< long > m_pickMcTk
AIDA::IHistogram1D * g_pickHitWire
Definition MdcHistItem.h:40
NTuple::Item< long > m_pickNCellWithDigi
NTuple::Array< double > m_pickPhiLowCut
double haveDigiDrift[43][288]
NTuple::Array< double > m_pickDriftCut
NTuple::Array< double > m_pickCurv
NTuple::Array< double > m_pickDrift
NTuple::Array< double > m_pickResid
NTuple::Array< int > m_pickPhizOk
NTuple::Array< double > m_pickSigma
NTuple::Array< long > m_pickWire
NTuple::Tuple * m_tuplePick
NTuple::Array< double > m_pickDriftTruth
NTuple::Array< double > m_pickZ
NTuple::Array< double > m_pickPredDrift
NTuple::Item< long > m_pickNCell
NTuple::Array< double > m_pickPhiHighCut
double posRad() const
Definition BesAngle.h:124
double deg() const
Definition BesAngle.h:121
static const double pi
Definition Constants.h:38
static const double twoPi
Definition Constants.h:39
static const double c
Definition Constants.h:43
static const int maxCell[43]
Definition Constants.h:40
static const double radToDegrees
Definition Constants.h:41
const MdcLayer * prevLayer(int lay) const
Definition MdcDetector.h:39
const MdcLayer * lastLayer() const
Definition MdcDetector.h:37
const MdcLayer * firstLayer() const
Definition MdcDetector.h:36
const MdcLayer * nextLayer(int lay) const
Definition MdcDetector.h:38
int ambig() const
const MdcLayer * layer() const
const MdcDigi * digi() const
Definition MdcHit.h:55
void print(std::ostream &o) const
Definition MdcHit.cxx:121
double x() const
Definition MdcHit.h:76
double phi() const
Definition MdcHit.h:75
double sigma(double, int, double, double, double) const
Definition MdcHit.cxx:184
double y() const
Definition MdcHit.h:77
double driftDist(double, int, double, double, double) const
Definition MdcHit.cxx:156
const MdcLayer * layer() const
Definition MdcHit.h:56
double rEnd(void) const
Definition MdcLayer.h:37
double rOut(void) const
Definition MdcLayer.h:39
double rIn(void) const
Definition MdcLayer.h:38
double pickHitPhiFactor
double maxInactiveResid
void setFirstLayer(const MdcLayer *l)
Definition MdcTrack.h:35
int projectToR(double radius, BesAngle &phiIntersect, int lCurl=0) const
Definition MdcTrack.cxx:63
int hasCurled() const
Definition MdcTrack.h:30
void setHasCurled(bool c=true)
Definition MdcTrack.h:34
const MdcLayer * firstLayer() const
Definition MdcTrack.h:31
const MdcLayer * lastLayer() const
Definition MdcTrack.h:32
void setLastLayer(const MdcLayer *l)
Definition MdcTrack.h:36
int getTrackIndex() const
Definition RawData.cxx:50
double phi0() const
double tanDip() const
bool is2d() const
const TrkHitOnTrk * getHitOnTrack(const TrkRecoTrk *trk) const
bool usedHit(void) const
Definition TrkFundHit.h:57
hot_iterator begin() const
Definition TrkHitList.h:45
TrkHitOnTrk * appendHot(const TrkHitOnTrk *theHot)
hot_iterator end() const
Definition TrkHitList.h:46
bool removeHit(const TrkFundHit *theHit)
double fltLen() const
Definition TrkHitOnTrk.h:91
virtual const MdcHitOnTrack * mdcHitOnTrack() const
void print(std::ostream &o) const
double trackT0() const
const BField & bField() const
Definition TrkRecoTrk.h:82
double y[1000]
int mdcWrapWire(int wireIn, int nCell)
Definition mdcWrapWire.h:6

Referenced by finishCircle(), and finishHelix().


The documentation for this class was generated from the following files: