CGEM BOSS 6.6.5.i
BESIII Offline Software System
Loading...
Searching...
No Matches
Public Member Functions | Protected Attributes | List of all members
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)
 

Protected Attributes

float addcut
 
int ncalc
 
int nadded
 
double sumpull
 
int kkk
 
int kkl
 
int thot
 
int tuot
 
int ntracks
 
int nhits
 
HepAList< MdcxHel > tttt
 
HepAList< MdcxHit > hhhh
 
HepAList< MdcxHit > listoasses
 

Detailed Description

Definition at line 33 of file MdcxAddHits.h.

Constructor & Destructor Documentation

◆ MdcxAddHits() [1/2]

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

Definition at line 40 of file MdcxAddHits.cxx.

41 :
42 ncalc(0), nadded(0), sumpull(0.0), thot(0), tuot(0)
43{
44 ntracks = trkl.length();
45 nhits = hitl.length();
46 addcut = addit;
47 // cout << "MdcxAddHits called with " << ntracks << " tracks, "
48 // << nhits << " hits, addit = " << addit << endl;
49 if ( (ntracks < 1) || (nhits < 1) ) return;
50 kkk = 0;
51 kkl = 0;
52 while(hitl[kkl]) {
53 hitl[kkl]->SetUsedOnHel(0);
54 hhhh.append(hitl[kkl++]);
55 }
56 while(trkl[kkk]) {
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::sumpull
double sumpull
Definition MdcxAddHits.h:53
MdcxAddHits::hhhh
HepAList< MdcxHit > hhhh
Definition MdcxAddHits.h:57
MdcxAddHits::tttt
HepAList< MdcxHel > tttt
Definition MdcxAddHits.h:56

◆ MdcxAddHits() [2/2]

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

Definition at line 66 of file MdcxAddHits.cxx.

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

◆ ~MdcxAddHits()

MdcxAddHits::~MdcxAddHits ( )
virtual

Definition at line 83 of file MdcxAddHits.cxx.

83 {
84}

Member Function Documentation

◆ GetAssociates()

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

zoujh

FIXME epsilon ? yzhang

FIXME

Definition at line 86 of file MdcxAddHits.cxx.

86 {
87 int debug = MdcxParameters::debug;
88
89 listoasses.removeAll(); /// zoujh
90 if ((whichtrack >= ntracks) || (whichtrack < 0)) {
91 cout << "asked for associates of track " << whichtrack
92 << " while there are only " << ntracks << " tracks in list " << endl;
93 return listoasses;
94 }
95
96 double m_2pi = 2.0*M_PI;
97 MdcxHel* temphel = tttt[whichtrack];
98
99 if(5 == debug) temphel->print();
100
101 double lmax = temphel->Lmax();
102
103 // calc. phim, phip
104 double gvin = MdcxParameters::firstMdcAxialRadius;//yzhang 2010-05-05
105 double gvout = MdcxParameters::maxMdcRadius;
106 double phiex = 0.38; ///FIXME epsilon ? yzhang
107 double phim = -M_PI;
108 double phip = M_PI;
109 double rmin = fabs(temphel->D0());
110 double rmax, pc, rc, rhel;
111 if (temphel->Ominfl()) {
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 // case A (exiter)
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 }//end while
166 } else if ((rmin>gvin) && (rmax>gvout)) {
167 if(5 == debug) std::cout<<" case B (albedo) "<<std::endl;
168 // case B (albedo)
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 // case C (looper)
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 // going to have to step to find min, max
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 //yzhang add 2011-10-10
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 //zhangy add
224 }//end while
225 }//end rc<>rhel
226 } else if ((rmin>gvin) && (rmax<gvout)) {
227 // case D (contained)
228 if (rc > rhel) {
229 double alp = asin(rhel/rc);
230 phim = pc - alp;
231 phip = pc + alp;
232 }
233 // if rc<rhel and case D, it's an orbiter. Use max phim, phip, so do nothing
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 // now try to add hits
243 kkl = 0;
244 if(5 == debug) std::cout<<"===== try to add hits:"<< std::endl;
245 while(hhhh[kkl]) {
246 MdcxHit* temphit = hhhh[kkl++];
247 //yzhang add 2011-10-11
248 if(temphel->Doca_Len()<0) {
249 if(5 == debug) std::cout<< " len<0 " << temphel->Doca_Len()<< " continue"<<std::endl;
250 continue;
251 }
252 //zhangy
253 double factor=1.;
254 //if ((0 != temphit->type()) && (temphit->Layer()>7)) factor = MdcxParameters::addHitFactor;//yzhang FIXME 2009-11-03
255 if(5 == debug) {
256 temphit->print(std::cout);
257 std::cout<< " pull "<<temphit->pull(*temphel)
258 << " nsig "<<factor * MdcxParameters::nSigAddHitTrk<< " len "<<temphel->Doca_Len() <<std::endl;
259 }
260 //yzhang 2009-10-21 16:55:25 ///FIXME
261 if((!temphit->GetUsedOnHel())&&(fabs(temphit->pull(*temphel))< factor * MdcxParameters::nSigAddHitTrk) ) {
262
263 //if( (!temphit->GetUsedOnHel()) && (fabs(temphit->d())<1.2) ) //delete
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){
272 std::cout<< " used "<<temphit->GetUsedOnHel()
273 <<" pull " << fabs(temphit->pull(*temphel))
274 <<" <? nSig " << MdcxParameters::nSigAddHitTrk
275 << std::endl;
276 }
277 temphit->SetConstErr(0);
278 double pull = temphit->pull( *temphel );
279 ncalc++;
280 sumpull += fabs(pull);
281 //cout << "MdcxAddHits: hit " << kkl-1 << " trk " << whichtrack << " pull " << pull << endl;
282 if(g_addHitCut) g_addHitCut->fill(fabs(pull));
283 int layer = temphit->Layer();
284 if(g_addHitCut2d) g_addHitCut2d->fill(layer,fabs(pull));
285 if(5 == debug) {
286 std::cout<< " pull "<<pull
287 << " addcut "<< MdcxParameters::nSigAddHitTrk
288 << "* factor:"<< factor <<"="<<factor * MdcxParameters::nSigAddHitTrk
289 << " len "<< temphel->Doca_Len()
290 << " >? lmax "<< lmax
291 << " >? maxTkLen "<< MdcxParameters::maxTrkLength
292 << std::endl;
293 }
294 // yzhang 2009-10-21 22:55:26
295 if(fabs(pull) < factor * MdcxParameters::nSigAddHitTrk) {
296 //if(fabs(pull) < addcut)
297 double len = temphel->Doca_Len();
298 //{ hhhh[kkl-1]->print(std::cout);
299 // cout << " trk " << whichtrack << " pull " << pull
300 // << " d() " << temphit->d() << " len " << len << endl; }
301 //int theflag = temphel->GetTurnFlag(); //zoujh ???
302 if (((len>0.0)&&(len<lmax))||((len<0.0)&&(len>-MdcxParameters::maxTrkLength))) { ///FIXME
303 //templist.append(temphit);
304 listoasses.append(temphit);
305 temphit->SetUsedOnHel(1);
306 nadded++;
307 if(debug>2){
308 temphit->print(std::cout);
309 std::cout<<"Added "<< std::endl;
310 }
311 }
312 }
313 temphit->SetConstErr(1);
314 }//end phim,phip cuts
315 }
316 }
317 //listoasses = templist;
318 return listoasses;
319}
g_addHitCut
AIDA::IHistogram1D * g_addHitCut
Definition MdcxHistItem.h:58
g_addHitCut2d
AIDA::IHistogram2D * g_addHitCut2d
Definition MdcxHistItem.h:59
M_PI
#define M_PI
Definition TConstant.h:4
MdcxAddHits::listoasses
HepAList< MdcxHit > listoasses
Definition MdcxAddHits.h:58
MdcxHel::Phi0
double Phi0() const
Definition MdcxHel.h:54
MdcxHel::D0
double D0() const
Definition MdcxHel.h:53
MdcxHel::Lmax
double Lmax() const
Definition MdcxHel.cxx:142
MdcxHel::Yc
double Yc() const
Definition MdcxHel.cxx:68
MdcxHel::Xc
double Xc() const
Definition MdcxHel.cxx:59
MdcxHel::Omega
double Omega() const
Definition MdcxHel.h:55
MdcxHel::Ominfl
int Ominfl() const
Definition MdcxHel.h:76
MdcxHel::print
void print() const
Definition MdcxHel.cxx:357
MdcxHel::Yh
double Yh(double l) const
Definition MdcxHel.cxx:94
MdcxHel::Xh
double Xh(double l) const
Definition MdcxHel.cxx:85
MdcxHel::Doca_Len
double Doca_Len() const
Definition MdcxHel.h:65
MdcxHit::pull
float pull(MdcxHel &hel) const
Definition MdcxHit.cxx:174
MdcxHit::Layer
int Layer() const
Definition MdcxHit.h:67
MdcxHit::pw
float pw() const
Definition MdcxHit.h:78
MdcxHit::SetConstErr
void SetConstErr(int i)
Definition MdcxHit.h:85
MdcxHit::print
void print(std::ostream &o, int i=0) const
Definition MdcxHit.cxx:205
MdcxHit::GetUsedOnHel
int GetUsedOnHel() const
Definition MdcxHit.h:103
MdcxHit::SetUsedOnHel
void SetUsedOnHel(int i)
Definition MdcxHit.h:102
MdcxParameters::firstMdcAxialRadius
static const double firstMdcAxialRadius
Definition MdcxParameters.h:34
MdcxParameters::maxTrkLength
static const double maxTrkLength
Definition MdcxParameters.h:13
MdcxParameters::maxMdcRadius
static const double maxMdcRadius
MDC Geometry.
Definition MdcxParameters.h:33
MdcxParameters::debug
static int debug
Definition MdcxParameters.h:44
MdcxParameters::nSigAddHitTrk
static double nSigAddHitTrk
Definition MdcxParameters.h:22

Member Data Documentation

◆ addcut

float MdcxAddHits::addcut
protected

Definition at line 51 of file MdcxAddHits.h.

Referenced by MdcxAddHits(), and MdcxAddHits().

◆ hhhh

HepAList<MdcxHit> MdcxAddHits::hhhh
protected

Definition at line 57 of file MdcxAddHits.h.

Referenced by GetAssociates(), MdcxAddHits(), and MdcxAddHits().

◆ kkk

int MdcxAddHits::kkk
protected

Definition at line 54 of file MdcxAddHits.h.

Referenced by MdcxAddHits().

◆ kkl

int MdcxAddHits::kkl
protected

Definition at line 54 of file MdcxAddHits.h.

Referenced by GetAssociates(), and MdcxAddHits().

◆ listoasses

HepAList<MdcxHit> MdcxAddHits::listoasses
protected

Definition at line 58 of file MdcxAddHits.h.

Referenced by GetAssociates().

◆ nadded

int MdcxAddHits::nadded
protected

Definition at line 52 of file MdcxAddHits.h.

Referenced by GetAssociates().

◆ ncalc

int MdcxAddHits::ncalc
protected

Definition at line 52 of file MdcxAddHits.h.

Referenced by GetAssociates().

◆ nhits

int MdcxAddHits::nhits
protected

Definition at line 55 of file MdcxAddHits.h.

Referenced by MdcxAddHits(), and MdcxAddHits().

◆ ntracks

int MdcxAddHits::ntracks
protected

Definition at line 55 of file MdcxAddHits.h.

Referenced by GetAssociates(), MdcxAddHits(), and MdcxAddHits().

◆ sumpull

double MdcxAddHits::sumpull
protected

Definition at line 53 of file MdcxAddHits.h.

Referenced by GetAssociates().

◆ thot

int MdcxAddHits::thot
protected

Definition at line 54 of file MdcxAddHits.h.

Referenced by MdcxAddHits().

◆ tttt

HepAList<MdcxHel> MdcxAddHits::tttt
protected

Definition at line 56 of file MdcxAddHits.h.

Referenced by GetAssociates(), MdcxAddHits(), and MdcxAddHits().

◆ tuot

int MdcxAddHits::tuot
protected

Definition at line 54 of file MdcxAddHits.h.

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