00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023 #include <math.h>
00024 #include "GaudiKernel/MsgStream.h"
00025 #include "GaudiKernel/AlgFactory.h"
00026 #include "MdcxReco/MdcxAddHits.h"
00027 #include "MdcxReco/MdcxHel.h"
00028 #include "MdcxReco/MdcxHit.h"
00029 #include "MdcxReco/MdcxFittedHel.h"
00030 #include "MdcxReco/MdcxParameters.h"
00031 #include "CLHEP/Alist/AIterator.h"
00032 using std::cout;
00033 using std::endl;
00034 #include "AIDA/IHistogram1D.h"
00035 #include "AIDA/IHistogram2D.h"
00036
00037 extern AIDA::IHistogram1D* g_addHitCut;
00038 extern AIDA::IHistogram2D* g_addHitCut2d;
00039
00040 MdcxAddHits::MdcxAddHits(HepAList<MdcxFittedHel> &trkl,
00041 const HepAList<MdcxHit> &hitl, float addit):
00042 ncalc(0), nadded(0), sumpull(0.0), thot(0), tuot(0)
00043 {
00044 ntracks = trkl.length();
00045 nhits = hitl.length();
00046 addcut = addit;
00047
00048
00049 if ( (ntracks < 1) || (nhits < 1) ) return;
00050 kkk = 0;
00051 kkl = 0;
00052 while(hitl[kkl]) {
00053 hitl[kkl]->SetUsedOnHel(0);
00054 hhhh.append(hitl[kkl++]);
00055 }
00056 while(trkl[kkk]) {
00057 MdcxHel* thel = trkl[kkk];
00058 tttt.append(thel);
00059 const HepAList<MdcxHit>& dchitlist = trkl[kkk++]->XHitList();
00060 thot += dchitlist.length();
00061 kkl = 0;
00062 while (dchitlist[kkl]) dchitlist[kkl++]->SetUsedOnHel(1);
00063 }
00064 }
00065
00066 MdcxAddHits::MdcxAddHits(HepAList<MdcxHel> &trkl, const HepAList<MdcxHit> &hitl, float addit):
00067 ncalc(0), nadded(0), sumpull(0.0)
00068 {
00069
00070
00071
00072
00073 ntracks = trkl.length();
00074 nhits = hitl.length();
00075 addcut = addit;
00076
00077
00078 if ( (ntracks < 1) || (nhits < 1) ) return;
00079 hhhh = hitl;
00080 tttt = trkl;
00081 }
00082
00083 MdcxAddHits::~MdcxAddHits() {
00084 }
00085
00086 const HepAList<MdcxHit>& MdcxAddHits::GetAssociates(int whichtrack) {
00087 int debug = MdcxParameters::debug;
00088
00089 listoasses.removeAll();
00090 if ((whichtrack >= ntracks) || (whichtrack < 0)) {
00091 cout << "asked for associates of track " << whichtrack
00092 << " while there are only " << ntracks << " tracks in list " << endl;
00093 return listoasses;
00094 }
00095
00096 double m_2pi = 2.0*M_PI;
00097 MdcxHel* temphel = tttt[whichtrack];
00098
00099 if(5 == debug) temphel->print();
00100
00101 double lmax = temphel->Lmax();
00102
00103
00104 double gvin = MdcxParameters::firstMdcAxialRadius;
00105 double gvout = MdcxParameters::maxMdcRadius;
00106 double phiex = 0.38;
00107 double phim = -M_PI;
00108 double phip = M_PI;
00109 double rmin = fabs(temphel->D0());
00110 double rmax, pc, rc, rhel;
00111 if (temphel->Ominfl()) {
00112 rmax = fabs(temphel->D0() + 2.0/temphel->Omega());
00113 double xc = temphel->Xc();
00114 double yc = temphel->Yc();
00115 pc = atan2(yc,xc);
00116 rc = fabs(temphel->D0() + 1.0/temphel->Omega());
00117 rhel = fabs(1.0/temphel->Omega());
00118 } else {
00119 rmax = 1000.;
00120 pc = temphel->Phi0();
00121 rc = 1000.;
00122 rhel = 1000.;
00123 }
00124 if(5 == debug) {
00125 std::cout<< "==Test add hit phi: rmin >?"<<rmin << " gvin "<<gvin
00126 << " rmax >?"<<rmax << " gvout "<<gvout << std::endl;
00127 }
00128
00129 if ((rmin<gvin) && (rmax>gvout)) {
00130
00131 if(5 == debug) std::cout<<" case A (exiter) "<<std::endl;
00132
00133 double len = 0;
00134 double lstep = m_2pi*rhel/16.;
00135 if (lstep>10.) lstep = 10.;
00136 double xl = temphel->Xh(len);
00137 double yl = temphel->Yh(len);
00138 double phil = atan2(yl, xl);
00139 double rl = sqrt(xl*xl + yl*yl);
00140 int nstep = 0;
00141 double philast = phil;
00142 int isin = 0;
00143 int notout = 1;
00144 while ((notout) && (nstep<20)) {
00145 len += lstep;
00146 nstep++;
00147 xl = temphel->Xh(len);
00148 yl = temphel->Yh(len);
00149 phil = atan2(yl, xl);
00150 rl = sqrt(xl*xl + yl*yl);
00151 if ((rl<gvin) && (!isin)) {
00152 philast = phil;
00153 } else if ((rl>gvin) && (!isin)) {
00154 isin = 1; phim = philast; phip = philast;
00155 }
00156 if (isin) {
00157 double deltap = phil - philast;
00158 if (deltap > M_PI) phil -= m_2pi;
00159 if (deltap < -M_PI) phil += m_2pi;
00160 philast = phil;
00161 if (rl > gvout) notout = 0;
00162 if (phil < phim) phim = phil;
00163 if (phil > phip) phip = phil;
00164 }
00165 }
00166 } else if ((rmin>gvin) && (rmax>gvout)) {
00167 if(5 == debug) std::cout<<" case B (albedo) "<<std::endl;
00168
00169 double len=0; double lstep=m_2pi*rhel/16.; if (lstep>10.)lstep=10.;
00170 double xl=temphel->Xh(len); double yl=temphel->Yh(len);
00171 double phil=atan2(yl,xl); double rl=sqrt(xl*xl+yl*yl);
00172 phim=phil; phip=phil; double phis=phil; double deltap,dp1,dp2;
00173 int nstep=0; double philast=phil;
00174 while ((rl<gvout)&&(nstep<20)){
00175 len+=lstep; nstep++; xl=temphel->Xh(len); yl=temphel->Yh(len);
00176 phil=atan2(yl,xl); rl=sqrt(xl*xl+yl*yl);
00177 deltap=phil-philast; if (deltap>M_PI)phil-=m_2pi;
00178 if (deltap<-M_PI)phil+=m_2pi; philast=phil;
00179 if (phil<phim)phim=phil; if (phil>phip)phip=phil;
00180 }
00181 dp1=fabs(phis-phim); dp2=fabs(phip-phis); deltap=dp1;
00182 if(dp2>dp1)deltap=dp2; phim=phis-deltap; phip=phis+deltap;
00183 } else if ((rmin<gvin) && (rmax<gvout)) {
00184
00185 if(5 == debug) std::cout<<" case C (looper) "<<std::endl;
00186 if (rc>rhel){
00187 double alp=asin(rhel/rc); phim=pc-alp; phip=pc+alp;
00188 }else{
00189
00190 double len=0;
00191 double lstep=m_2pi*rhel/16.;
00192 if (lstep>10.)lstep=10.;
00193 if(5 == debug) std::cout<< "ini step "<<m_2pi*rhel/16.<<" lstep " << lstep <<"cm"<<std::endl;
00194 double xl=temphel->Xh(len); double yl=temphel->Yh(len);
00195 double phil=atan2(yl,xl); double rl=sqrt(xl*xl+yl*yl);
00196 int nstep=0; double philast=phil; int isin=0; int notout=1;
00197 while ((notout)&&(nstep<33)){
00198 len+=lstep; nstep++; xl=temphel->Xh(len); yl=temphel->Yh(len);
00199 phil=atan2(yl,xl); rl=sqrt(xl*xl+yl*yl);
00200 if(5 == debug) {
00201 std::cout<< nstep<<" rl "<<rl<< " gvin " <<gvin<< " isin "<<isin;
00202 }
00203 if ((rl<gvin)&&(!isin)){
00204 philast=phil;
00205 }else if((rl>gvin)&&(!isin)){
00206 isin=1; phim=philast; phip=philast;
00207 }
00208 if (isin){
00209 double deltap=phil-philast; if (deltap>M_PI)phil-=m_2pi;
00210 if (deltap<-M_PI)phil+=m_2pi; philast=phil; if (rl<gvin)notout=0;
00211 if (phil<phim)phim=phil; if (phil>phip)phip=phil;
00212 }
00213
00214 if(len > M_PI*rhel*0.75) {
00215 if(5 == debug) {
00216 std::cout<< " len "<<len<<" >pi*R_helix "<<M_PI*rhel<< " rhel "<<rhel<< " break"<<std::endl;
00217 }
00218 break;
00219 }
00220 if(5 == debug) {
00221 std::cout<< " philast "<<philast<<" phim "<<phim << " phip "<<phip <<" len "<<len<<std::endl;
00222 }
00223
00224 }
00225 }
00226 } else if ((rmin>gvin) && (rmax<gvout)) {
00227
00228 if (rc > rhel) {
00229 double alp = asin(rhel/rc);
00230 phim = pc - alp;
00231 phip = pc + alp;
00232 }
00233
00234 }
00235 phim -= phiex;
00236 phip += phiex;
00237
00238 if(5 == debug) {
00239 std::cout<<"add hits phim "<<phim <<" phip "<<phip<< std::endl;
00240 }
00241
00242
00243 kkl = 0;
00244 if(5 == debug) std::cout<<"===== try to add hits:"<< std::endl;
00245 while(hhhh[kkl]) {
00246 MdcxHit* temphit = hhhh[kkl++];
00247
00248 if(temphel->Doca_Len()<0) {
00249 if(5 == debug) std::cout<< " len<0 " << temphel->Doca_Len()<< " continue"<<std::endl;
00250 continue;
00251 }
00252
00253 double factor=1.;
00254
00255 if(5 == debug) {
00256 temphit->print(std::cout);
00257 std::cout<< " pull "<<temphit->pull(*temphel)
00258 << " nsig "<<factor * MdcxParameters::nSigAddHitTrk<< " len "<<temphel->Doca_Len() <<std::endl;
00259 }
00260
00261 if((!temphit->GetUsedOnHel())&&(fabs(temphit->pull(*temphel))< factor * MdcxParameters::nSigAddHitTrk) ) {
00262
00263
00264 double pw = temphit->pw();
00265 if((phip-pw) > m_2pi) pw += m_2pi;
00266 if((phim-pw) < -m_2pi) pw -= m_2pi;
00267 if(5 == debug) {
00268 std::cout<< "phi wire "<<pw << " phim "<<phim<< " phip "<<phip<<" len "<<temphel->Doca_Len();
00269 }
00270 if ( (pw>phim) && (pw<phip) ) {
00271 if (5 == debug){
00272 std::cout<< " used "<<temphit->GetUsedOnHel()
00273 <<" pull " << fabs(temphit->pull(*temphel))
00274 <<" <? nSig " << MdcxParameters::nSigAddHitTrk
00275 << std::endl;
00276 }
00277 temphit->SetConstErr(0);
00278 double pull = temphit->pull( *temphel );
00279 ncalc++;
00280 sumpull += fabs(pull);
00281
00282 if(g_addHitCut) g_addHitCut->fill(fabs(pull));
00283 int layer = temphit->Layer();
00284 if(g_addHitCut2d) g_addHitCut2d->fill(layer,fabs(pull));
00285 if(5 == debug) {
00286 std::cout<< " pull "<<pull
00287 << " addcut "<< MdcxParameters::nSigAddHitTrk
00288 << "* factor:"<< factor <<"="<<factor * MdcxParameters::nSigAddHitTrk
00289 << " len "<< temphel->Doca_Len()
00290 << " >? lmax "<< lmax
00291 << " >? maxTkLen "<< MdcxParameters::maxTrkLength
00292 << std::endl;
00293 }
00294
00295 if(fabs(pull) < factor * MdcxParameters::nSigAddHitTrk) {
00296
00297 double len = temphel->Doca_Len();
00298
00299
00300
00301
00302 if (((len>0.0)&&(len<lmax))||((len<0.0)&&(len>-MdcxParameters::maxTrkLength))) {
00303
00304 listoasses.append(temphit);
00305 temphit->SetUsedOnHel(1);
00306 nadded++;
00307 if(debug>2){
00308 temphit->print(std::cout);
00309 std::cout<<"Added "<< std::endl;
00310 }
00311 }
00312 }
00313 temphit->SetConstErr(1);
00314 }
00315 }
00316 }
00317
00318 return listoasses;
00319 }