BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
MdcxAddHits Class Reference

#include <MdcxAddHits.h>

Public Member Functions

 MdcxAddHits (HepAList< MdcxFittedHel > &f, const HepAList< MdcxHit > &h, float addit=1.5)
 MdcxAddHits (HepAList< MdcxHel > &f, const HepAList< MdcxHit > &h, float addit=1.5)
virtual ~MdcxAddHits ()
const HepAList< MdcxHit > & GetAssociates (int i=0)
 MdcxAddHits (HepAList< MdcxFittedHel > &f, const HepAList< MdcxHit > &h, float addit=1.5)
 MdcxAddHits (HepAList< MdcxHel > &f, const HepAList< MdcxHit > &h, float addit=1.5)
virtual ~MdcxAddHits ()
const HepAList< MdcxHit > & GetAssociates (int i=0)
 MdcxAddHits (HepAList< MdcxFittedHel > &f, const HepAList< MdcxHit > &h, float addit=1.5)
 MdcxAddHits (HepAList< MdcxHel > &f, const HepAList< MdcxHit > &h, float addit=1.5)
virtual ~MdcxAddHits ()
const HepAList< MdcxHit > & GetAssociates (int i=0)

Protected Attributes

float addcut
int ncalc
int nadded
double sumpull
int kkk
int kkl
int thot
int tuot
int ntracks
int nhits
HepAList< MdcxHeltttt
HepAList< MdcxHithhhh
HepAList< MdcxHitlistoasses

Detailed Description

Constructor & Destructor Documentation

◆ MdcxAddHits() [1/6]

MdcxAddHits::MdcxAddHits ( HepAList< MdcxFittedHel > & f,
const HepAList< MdcxHit > & h,
float addit = 1.5 )

Definition at line 39 of file MdcxAddHits.cxx.

41 : ncalc( 0 ), nadded( 0 ), sumpull( 0.0 ), thot( 0 ), tuot( 0 ) {
42 ntracks = trkl.length();
43 nhits = hitl.length();
44 addcut = addit;
45 // cout << "MdcxAddHits called with " << ntracks << " tracks, "
46 // << nhits << " hits, addit = " << addit << endl;
47 if ( ( ntracks < 1 ) || ( nhits < 1 ) ) return;
48 kkk = 0;
49 kkl = 0;
50 while ( hitl[kkl] )
51 {
52 hitl[kkl]->SetUsedOnHel( 0 );
53 hhhh.append( hitl[kkl++] );
54 }
55 while ( trkl[kkk] )
56 {
57 MdcxHel* thel = trkl[kkk];
58 tttt.append( thel );
59 const HepAList<MdcxHit>& dchitlist = trkl[kkk++]->XHitList();
60 thot += dchitlist.length();
61 kkl = 0;
62 while ( dchitlist[kkl] ) dchitlist[kkl++]->SetUsedOnHel( 1 );
63 }
64}

◆ MdcxAddHits() [2/6]

MdcxAddHits::MdcxAddHits ( HepAList< MdcxHel > & f,
const HepAList< MdcxHit > & h,
float addit = 1.5 )

Definition at line 66 of file MdcxAddHits.cxx.

67 : ncalc( 0 ), nadded( 0 ), sumpull( 0.0 ) {
68 //
69 // c-tor designed to work with MdcxHitAdder; all MdcxHit's in hitl are
70 // assumed to be fresh and fair game.
71 //
72 ntracks = trkl.length();
73 nhits = hitl.length();
74 addcut = addit;
75 // cout << "MdcxAddHits called with " << ntracks << " tracks, "
76 // << nhits << " hits, addit = " << addit << endl;
77 if ( ( ntracks < 1 ) || ( nhits < 1 ) ) return;
78 hhhh = hitl;
79 tttt = trkl;
80}

◆ ~MdcxAddHits() [1/3]

MdcxAddHits::~MdcxAddHits ( )
virtual

Definition at line 82 of file MdcxAddHits.cxx.

82{}

◆ MdcxAddHits() [3/6]

MdcxAddHits::MdcxAddHits ( HepAList< MdcxFittedHel > & f,
const HepAList< MdcxHit > & h,
float addit = 1.5 )

◆ MdcxAddHits() [4/6]

MdcxAddHits::MdcxAddHits ( HepAList< MdcxHel > & f,
const HepAList< MdcxHit > & h,
float addit = 1.5 )

◆ ~MdcxAddHits() [2/3]

virtual MdcxAddHits::~MdcxAddHits ( )
virtual

◆ MdcxAddHits() [5/6]

MdcxAddHits::MdcxAddHits ( HepAList< MdcxFittedHel > & f,
const HepAList< MdcxHit > & h,
float addit = 1.5 )

◆ MdcxAddHits() [6/6]

MdcxAddHits::MdcxAddHits ( HepAList< MdcxHel > & f,
const HepAList< MdcxHit > & h,
float addit = 1.5 )

◆ ~MdcxAddHits() [3/3]

virtual MdcxAddHits::~MdcxAddHits ( )
virtual

Member Function Documentation

◆ GetAssociates() [1/3]

const HepAList< MdcxHit > & MdcxAddHits::GetAssociates ( int i = 0)

zoujh

FIXME epsilon ? yzhang

FIXME

Definition at line 84 of file MdcxAddHits.cxx.

84 {
85 int debug = MdcxParameters::debug;
86
87 listoasses.removeAll(); /// zoujh
88 if ( ( whichtrack >= ntracks ) || ( whichtrack < 0 ) )
89 {
90 cout << "asked for associates of track " << whichtrack << " while there are only "
91 << ntracks << " tracks in list " << endl;
92 return listoasses;
93 }
94
95 double m_2pi = 2.0 * M_PI;
96 MdcxHel* temphel = tttt[whichtrack];
97
98 if ( 5 == debug ) temphel->print();
99
100 double lmax = temphel->Lmax();
101
102 // calc. phim, phip
103 double gvin = MdcxParameters::firstMdcAxialRadius; // yzhang 2010-05-05
104 double gvout = MdcxParameters::maxMdcRadius;
105 double phiex = 0.38; /// FIXME epsilon ? yzhang
106 double phim = -M_PI;
107 double phip = M_PI;
108 double rmin = fabs( temphel->D0() );
109 double rmax, pc, rc, rhel;
110 if ( temphel->Ominfl() )
111 {
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 }
119 else
120 {
121 rmax = 1000.;
122 pc = temphel->Phi0();
123 rc = 1000.;
124 rhel = 1000.;
125 }
126 if ( 5 == debug )
127 {
128 std::cout << "==Test add hit phi: rmin >?" << rmin << " gvin " << gvin << " rmax >?"
129 << rmax << " gvout " << gvout << std::endl;
130 }
131
132 if ( ( rmin < gvin ) && ( rmax > gvout ) )
133 {
134 // case A (exiter)
135 if ( 5 == debug ) std::cout << " case A (exiter) " << std::endl;
136
137 double len = 0;
138 double lstep = m_2pi * rhel / 16.;
139 if ( lstep > 10. ) lstep = 10.;
140 double xl = temphel->Xh( len );
141 double yl = temphel->Yh( len );
142 double phil = atan2( yl, xl );
143 double rl = sqrt( xl * xl + yl * yl );
144 int nstep = 0;
145 double philast = phil;
146 int isin = 0;
147 int notout = 1;
148 while ( ( notout ) && ( nstep < 20 ) )
149 {
150 len += lstep;
151 nstep++;
152 xl = temphel->Xh( len );
153 yl = temphel->Yh( len );
154 phil = atan2( yl, xl );
155 rl = sqrt( xl * xl + yl * yl );
156 if ( ( rl < gvin ) && ( !isin ) ) { philast = phil; }
157 else if ( ( rl > gvin ) && ( !isin ) )
158 {
159 isin = 1;
160 phim = philast;
161 phip = philast;
162 }
163 if ( isin )
164 {
165 double deltap = phil - philast;
166 if ( deltap > M_PI ) phil -= m_2pi;
167 if ( deltap < -M_PI ) phil += m_2pi;
168 philast = phil;
169 if ( rl > gvout ) notout = 0;
170 if ( phil < phim ) phim = phil;
171 if ( phil > phip ) phip = phil;
172 }
173 } // end while
174 }
175 else if ( ( rmin > gvin ) && ( rmax > gvout ) )
176 {
177 if ( 5 == debug ) std::cout << " case B (albedo) " << std::endl;
178 // case B (albedo)
179 double len = 0;
180 double lstep = m_2pi * rhel / 16.;
181 if ( lstep > 10. ) lstep = 10.;
182 double xl = temphel->Xh( len );
183 double yl = temphel->Yh( len );
184 double phil = atan2( yl, xl );
185 double rl = sqrt( xl * xl + yl * yl );
186 phim = phil;
187 phip = phil;
188 double phis = phil;
189 double deltap, dp1, dp2;
190 int nstep = 0;
191 double philast = phil;
192 while ( ( rl < gvout ) && ( nstep < 20 ) )
193 {
194 len += lstep;
195 nstep++;
196 xl = temphel->Xh( len );
197 yl = temphel->Yh( len );
198 phil = atan2( yl, xl );
199 rl = sqrt( xl * xl + yl * yl );
200 deltap = phil - philast;
201 if ( deltap > M_PI ) phil -= m_2pi;
202 if ( deltap < -M_PI ) phil += m_2pi;
203 philast = phil;
204 if ( phil < phim ) phim = phil;
205 if ( phil > phip ) phip = phil;
206 }
207 dp1 = fabs( phis - phim );
208 dp2 = fabs( phip - phis );
209 deltap = dp1;
210 if ( dp2 > dp1 ) deltap = dp2;
211 phim = phis - deltap;
212 phip = phis + deltap;
213 }
214 else if ( ( rmin < gvin ) && ( rmax < gvout ) )
215 {
216 // case C (looper)
217 if ( 5 == debug ) std::cout << " case C (looper) " << std::endl;
218 if ( rc > rhel )
219 {
220 double alp = asin( rhel / rc );
221 phim = pc - alp;
222 phip = pc + alp;
223 }
224 else
225 {
226 // going to have to step to find min, max
227 double len = 0;
228 double lstep = m_2pi * rhel / 16.;
229 if ( lstep > 10. ) lstep = 10.;
230 if ( 5 == debug )
231 std::cout << "ini step " << m_2pi * rhel / 16. << " lstep " << lstep << "cm"
232 << std::endl;
233 double xl = temphel->Xh( len );
234 double yl = temphel->Yh( len );
235 double phil = atan2( yl, xl );
236 double rl = sqrt( xl * xl + yl * yl );
237 int nstep = 0;
238 double philast = phil;
239 int isin = 0;
240 int notout = 1;
241 while ( ( notout ) && ( nstep < 33 ) )
242 {
243 len += lstep;
244 nstep++;
245 xl = temphel->Xh( len );
246 yl = temphel->Yh( len );
247 phil = atan2( yl, xl );
248 rl = sqrt( xl * xl + yl * yl );
249 if ( 5 == debug )
250 { std::cout << nstep << " rl " << rl << " gvin " << gvin << " isin " << isin; }
251 if ( ( rl < gvin ) && ( !isin ) ) { philast = phil; }
252 else if ( ( rl > gvin ) && ( !isin ) )
253 {
254 isin = 1;
255 phim = philast;
256 phip = philast;
257 }
258 if ( isin )
259 {
260 double deltap = phil - philast;
261 if ( deltap > M_PI ) phil -= m_2pi;
262 if ( deltap < -M_PI ) phil += m_2pi;
263 philast = phil;
264 if ( rl < gvin ) notout = 0;
265 if ( phil < phim ) phim = phil;
266 if ( phil > phip ) phip = phil;
267 }
268 // yzhang add 2011-10-10
269 if ( len > M_PI * rhel * 0.75 )
270 {
271 if ( 5 == debug )
272 {
273 std::cout << " len " << len << " >pi*R_helix " << M_PI * rhel << " rhel " << rhel
274 << " break" << std::endl;
275 }
276 break;
277 }
278 if ( 5 == debug )
279 {
280 std::cout << " philast " << philast << " phim " << phim << " phip " << phip
281 << " len " << len << std::endl;
282 }
283 // zhangy add
284 } // end while
285 } // end rc<>rhel
286 }
287 else if ( ( rmin > gvin ) && ( rmax < gvout ) )
288 {
289 // case D (contained)
290 if ( rc > rhel )
291 {
292 double alp = asin( rhel / rc );
293 phim = pc - alp;
294 phip = pc + alp;
295 }
296 // if rc<rhel and case D, it's an orbiter. Use max phim, phip, so do nothing
297 }
298 phim -= phiex;
299 phip += phiex;
300
301 if ( 5 == debug ) { std::cout << "add hits phim " << phim << " phip " << phip << std::endl; }
302
303 // now try to add hits
304 kkl = 0;
305 if ( 5 == debug ) std::cout << "===== try to add hits:" << std::endl;
306 while ( hhhh[kkl] )
307 {
308 MdcxHit* temphit = hhhh[kkl++];
309 // yzhang add 2011-10-11
310 if ( temphel->Doca_Len() < 0 )
311 {
312 if ( 5 == debug )
313 std::cout << " len<0 " << temphel->Doca_Len() << " continue" << std::endl;
314 continue;
315 }
316 // zhangy
317 double factor = 1.;
318 // if ((0 != temphit->type()) && (temphit->Layer()>7)) factor =
319 // MdcxParameters::addHitFactor;//yzhang FIXME 2009-11-03
320 if ( 5 == debug )
321 {
322 temphit->print( std::cout );
323 std::cout << " pull " << temphit->pull( *temphel ) << " nsig "
324 << factor * MdcxParameters::nSigAddHitTrk << " len " << temphel->Doca_Len()
325 << std::endl;
326 }
327 // yzhang 2009-10-21 16:55:25 ///FIXME
328 if ( ( !temphit->GetUsedOnHel() ) &&
329 ( fabs( temphit->pull( *temphel ) ) < factor * MdcxParameters::nSigAddHitTrk ) )
330 {
331
332 // if( (!temphit->GetUsedOnHel()) && (fabs(temphit->d())<1.2) ) //delete
333 double pw = temphit->pw();
334 if ( ( phip - pw ) > m_2pi ) pw += m_2pi;
335 if ( ( phim - pw ) < -m_2pi ) pw -= m_2pi;
336 if ( 5 == debug )
337 {
338 std::cout << "phi wire " << pw << " phim " << phim << " phip " << phip << " len "
339 << temphel->Doca_Len();
340 }
341 if ( ( pw > phim ) && ( pw < phip ) )
342 {
343 if ( 5 == debug )
344 {
345 std::cout << " used " << temphit->GetUsedOnHel() << " pull "
346 << fabs( temphit->pull( *temphel ) ) << " <? nSig "
347 << MdcxParameters::nSigAddHitTrk << std::endl;
348 }
349 temphit->SetConstErr( 0 );
350 double pull = temphit->pull( *temphel );
351 ncalc++;
352 sumpull += fabs( pull );
353 // cout << "MdcxAddHits: hit " << kkl-1 << " trk " << whichtrack << " pull " << pull <<
354 // endl;
355 if ( g_addHitCut ) g_addHitCut->fill( fabs( pull ) );
356 int layer = temphit->Layer();
357 if ( g_addHitCut2d ) g_addHitCut2d->fill( layer, fabs( pull ) );
358 if ( 5 == debug )
359 {
360 std::cout << " pull " << pull << " addcut " << MdcxParameters::nSigAddHitTrk
361 << "* factor:" << factor << "=" << factor * MdcxParameters::nSigAddHitTrk
362 << " len " << temphel->Doca_Len() << " >? lmax " << lmax << " >? maxTkLen "
363 << MdcxParameters::maxTrkLength << std::endl;
364 }
365 // yzhang 2009-10-21 22:55:26
366 if ( fabs( pull ) < factor * MdcxParameters::nSigAddHitTrk )
367 {
368 // if(fabs(pull) < addcut)
369 double len = temphel->Doca_Len();
370 //{ hhhh[kkl-1]->print(std::cout);
371 // cout << " trk " << whichtrack << " pull " << pull
372 // << " d() " << temphit->d() << " len " << len << endl; }
373 // int theflag = temphel->GetTurnFlag(); //zoujh ???
374 if ( ( ( len > 0.0 ) && ( len < lmax ) ) ||
375 ( ( len < 0.0 ) && ( len > -MdcxParameters::maxTrkLength ) ) )
376 { /// FIXME
377 // templist.append(temphit);
378 listoasses.append( temphit );
379 temphit->SetUsedOnHel( 1 );
380 nadded++;
381 if ( debug > 2 )
382 {
383 temphit->print( std::cout );
384 std::cout << "Added " << std::endl;
385 }
386 }
387 }
388 temphit->SetConstErr( 1 );
389 } // end phim,phip cuts
390 }
391 }
392 // listoasses = templist;
393 return listoasses;
394}
character *LEPTONflag integer iresonances real zeta5 real * alp
AIDA::IHistogram1D * g_addHitCut
AIDA::IHistogram2D * g_addHitCut2d
#define M_PI
Definition TConstant.h:4
double Lmax() const
Definition MdcxHel.cxx:137
double Yc() const
Definition MdcxHel.cxx:76
double Xc() const
Definition MdcxHel.cxx:67
void print() const
Definition MdcxHel.cxx:403
double Yh(double l) const
Definition MdcxHel.cxx:98
double Xh(double l) const
Definition MdcxHel.cxx:89
float pull(MdcxHel &hel) const
Definition MdcxHit.cxx:168
void print(std::ostream &o, int i=0) const

◆ GetAssociates() [2/3]

const HepAList< MdcxHit > & MdcxAddHits::GetAssociates ( int i = 0)

◆ GetAssociates() [3/3]

const HepAList< MdcxHit > & MdcxAddHits::GetAssociates ( int i = 0)

Member Data Documentation

◆ addcut

float MdcxAddHits::addcut
protected

◆ hhhh

HepAList< MdcxHit > MdcxAddHits::hhhh
protected

◆ kkk

int MdcxAddHits::kkk
protected

◆ kkl

int MdcxAddHits::kkl
protected

◆ listoasses

HepAList< MdcxHit > MdcxAddHits::listoasses
protected

◆ nadded

int MdcxAddHits::nadded
protected

◆ ncalc

int MdcxAddHits::ncalc
protected

◆ nhits

int MdcxAddHits::nhits
protected

◆ ntracks

int MdcxAddHits::ntracks
protected

◆ sumpull

double MdcxAddHits::sumpull
protected

◆ thot

int MdcxAddHits::thot
protected

◆ tttt

HepAList< MdcxHel > MdcxAddHits::tttt
protected

◆ tuot

int MdcxAddHits::tuot
protected

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