86 {
88
90 if ((whichtrack >=
ntracks) || (whichtrack < 0)) {
91 cout << "asked for associates of track " << whichtrack
92 <<
" while there are only " <<
ntracks <<
" tracks in list " << endl;
94 }
95
96 double m_2pi = 2.0*
M_PI;
98
99 if(5 == debug) temphel->
print();
100
101 double lmax = temphel->
Lmax();
102
103
106 double phiex = 0.38;
109 double rmin = fabs(temphel->
D0());
110 double rmax, pc, rc, rhel;
112 rmax = fabs(temphel->
D0() + 2.0/temphel->
Omega());
113 double xc = temphel->
Xc();
114 double yc = temphel->
Yc();
115 pc = atan2(yc,xc);
116 rc = fabs(temphel->
D0() + 1.0/temphel->
Omega());
117 rhel = fabs(1.0/temphel->
Omega());
118 } else {
119 rmax = 1000.;
120 pc = temphel->
Phi0();
121 rc = 1000.;
122 rhel = 1000.;
123 }
124 if(5 == debug) {
125 std::cout<< "==Test add hit phi: rmin >?"<<rmin << " gvin "<<gvin
126 << " rmax >?"<<rmax << " gvout "<<gvout << std::endl;
127 }
128
129 if ((rmin<gvin) && (rmax>gvout)) {
130
131 if(5 == debug) std::cout<<" case A (exiter) "<<std::endl;
132
133 double len = 0;
134 double lstep = m_2pi*rhel/16.;
135 if (lstep>10.) lstep = 10.;
136 double xl = temphel->
Xh(len);
137 double yl = temphel->
Yh(len);
138 double phil = atan2(yl, xl);
139 double rl = sqrt(xl*xl + yl*yl);
140 int nstep = 0;
141 double philast = phil;
142 int isin = 0;
143 int notout = 1;
144 while ((notout) && (nstep<20)) {
145 len += lstep;
146 nstep++;
147 xl = temphel->
Xh(len);
148 yl = temphel->
Yh(len);
149 phil = atan2(yl, xl);
150 rl = sqrt(xl*xl + yl*yl);
151 if ((rl<gvin) && (!isin)) {
152 philast = phil;
153 } else if ((rl>gvin) && (!isin)) {
154 isin = 1; phim = philast; phip = philast;
155 }
156 if (isin) {
157 double deltap = phil - philast;
158 if (deltap >
M_PI) phil -= m_2pi;
159 if (deltap < -
M_PI) phil += m_2pi;
160 philast = phil;
161 if (rl > gvout) notout = 0;
162 if (phil < phim) phim = phil;
163 if (phil > phip) phip = phil;
164 }
165 }
166 } else if ((rmin>gvin) && (rmax>gvout)) {
167 if(5 == debug) std::cout<<" case B (albedo) "<<std::endl;
168
169 double len=0; double lstep=m_2pi*rhel/16.; if (lstep>10.)lstep=10.;
170 double xl=temphel->
Xh(len);
double yl=temphel->
Yh(len);
171 double phil=atan2(yl,xl); double rl=sqrt(xl*xl+yl*yl);
172 phim=phil; phip=phil; double phis=phil; double deltap,dp1,dp2;
173 int nstep=0; double philast=phil;
174 while ((rl<gvout)&&(nstep<20)){
175 len+=lstep; nstep++; xl=temphel->
Xh(len); yl=temphel->
Yh(len);
176 phil=atan2(yl,xl); rl=sqrt(xl*xl+yl*yl);
177 deltap=phil-philast;
if (deltap>
M_PI)phil-=m_2pi;
178 if (deltap<-
M_PI)phil+=m_2pi; philast=phil;
179 if (phil<phim)phim=phil; if (phil>phip)phip=phil;
180 }
181 dp1=fabs(phis-phim); dp2=fabs(phip-phis); deltap=dp1;
182 if(dp2>dp1)deltap=dp2; phim=phis-deltap; phip=phis+deltap;
183 } else if ((rmin<gvin) && (rmax<gvout)) {
184
185 if(5 == debug) std::cout<<" case C (looper) "<<std::endl;
186 if (rc>rhel){
187 double alp=asin(rhel/rc); phim=pc-
alp; phip=pc+
alp;
188 }else{
189
190 double len=0;
191 double lstep=m_2pi*rhel/16.;
192 if (lstep>10.)lstep=10.;
193 if(5 == debug) std::cout<< "ini step "<<m_2pi*rhel/16.<<" lstep " << lstep <<"cm"<<std::endl;
194 double xl=temphel->
Xh(len);
double yl=temphel->
Yh(len);
195 double phil=atan2(yl,xl); double rl=sqrt(xl*xl+yl*yl);
196 int nstep=0; double philast=phil; int isin=0; int notout=1;
197 while ((notout)&&(nstep<33)){
198 len+=lstep; nstep++; xl=temphel->
Xh(len); yl=temphel->
Yh(len);
199 phil=atan2(yl,xl); rl=sqrt(xl*xl+yl*yl);
200 if(5 == debug) {
201 std::cout<< nstep<<" rl "<<rl<< " gvin " <<gvin<< " isin "<<isin;
202 }
203 if ((rl<gvin)&&(!isin)){
204 philast=phil;
205 }else if((rl>gvin)&&(!isin)){
206 isin=1; phim=philast; phip=philast;
207 }
208 if (isin){
209 double deltap=phil-philast;
if (deltap>
M_PI)phil-=m_2pi;
210 if (deltap<-
M_PI)phil+=m_2pi; philast=phil;
if (rl<gvin)notout=0;
211 if (phil<phim)phim=phil; if (phil>phip)phip=phil;
212 }
213
214 if(len >
M_PI*rhel*0.75) {
215 if(5 == debug) {
216 std::cout<<
" len "<<len<<
" >pi*R_helix "<<
M_PI*rhel<<
" rhel "<<rhel<<
" break"<<std::endl;
217 }
218 break;
219 }
220 if(5 == debug) {
221 std::cout<< " philast "<<philast<<" phim "<<phim << " phip "<<phip <<" len "<<len<<std::endl;
222 }
223
224 }
225 }
226 } else if ((rmin>gvin) && (rmax<gvout)) {
227
228 if (rc > rhel) {
229 double alp = asin(rhel/rc);
232 }
233
234 }
235 phim -= phiex;
236 phip += phiex;
237
238 if(5 == debug) {
239 std::cout<<"add hits phim "<<phim <<" phip "<<phip<< std::endl;
240 }
241
242
244 if(5 == debug) std::cout<<"===== try to add hits:"<< std::endl;
247
249 if(5 == debug) std::cout<<
" len<0 " << temphel->
Doca_Len()<<
" continue"<<std::endl;
250 continue;
251 }
252
253 double factor=1.;
254
255 if(5 == debug) {
256 temphit->
print(std::cout);
257 std::cout<<
" pull "<<temphit->
pull(*temphel)
259 }
260
262
263
264 double pw = temphit->
pw();
265 if((phip-pw) > m_2pi) pw += m_2pi;
266 if((phim-pw) < -m_2pi) pw -= m_2pi;
267 if(5 == debug) {
268 std::cout<<
"phi wire "<<pw <<
" phim "<<phim<<
" phip "<<phip<<
" len "<<temphel->
Doca_Len();
269 }
270 if ( (pw>phim) && (pw<phip) ) {
271 if (5 == debug){
273 <<
" pull " << fabs(temphit->
pull(*temphel))
275 << std::endl;
276 }
278 double pull = temphit->
pull( *temphel );
281
283 int layer = temphit->
Layer();
285 if(5 == debug) {
286 std::cout<< " pull "<<pull
290 << " >? lmax "<< lmax
292 << std::endl;
293 }
294
296
298
299
300
301
303
307 if(debug>2){
308 temphit->
print(std::cout);
309 std::cout<<"Added "<< std::endl;
310 }
311 }
312 }
314 }
315 }
316 }
317
319}
character *LEPTONflag integer iresonances real zeta5 real * alp
AIDA::IHistogram1D * g_addHitCut
AIDA::IHistogram2D * g_addHitCut2d
HepAList< MdcxHit > listoasses
double Yh(double l) const
double Xh(double l) const
float pull(MdcxHel &hel) const
void print(std::ostream &o, int i=0) const
static const double firstMdcAxialRadius
static const double maxTrkLength
static const double maxMdcRadius
MDC Geometry.
static double nSigAddHitTrk