| File: | /home/sdobbs/work/clang/halld_recon/src/libraries/TRACKING/DTrackHitSelectorALT2.cc |
| Warning: | line 815, column 5 Value stored to 'var_lambda' is never read |
Press '?' to see keyboard shortcuts
Keyboard shortcuts:
| 1 | // $Id$ |
| 2 | // |
| 3 | // File: DTrackHitSelectorALT2.cc |
| 4 | // Created: Fri Feb 6 08:22:58 EST 2009 |
| 5 | // Creator: davidl (on Darwin harriet.jlab.org 9.6.0 i386) |
| 6 | // |
| 7 | |
| 8 | #include <cmath> |
| 9 | using namespace std; |
| 10 | |
| 11 | #include <TROOT.h> |
| 12 | |
| 13 | #include <TRACKING/DReferenceTrajectory.h> |
| 14 | |
| 15 | #include "DTrackHitSelectorALT2.h" |
| 16 | |
| 17 | #ifndef ansi_escape((char)0x1b) |
| 18 | #define ansi_escape((char)0x1b) ((char)0x1b) |
| 19 | #define ansi_bold((char)0x1b)<<"[1m" ansi_escape((char)0x1b)<<"[1m" |
| 20 | #define ansi_normal((char)0x1b)<<"[0m" ansi_escape((char)0x1b)<<"[0m" |
| 21 | #define ansi_red((char)0x1b)<<"[31m" ansi_escape((char)0x1b)<<"[31m" |
| 22 | #define ansi_green((char)0x1b)<<"[32m" ansi_escape((char)0x1b)<<"[32m" |
| 23 | #define ansi_blue((char)0x1b)<<"[34m" ansi_escape((char)0x1b)<<"[34m" |
| 24 | #endif // ansi_escape |
| 25 | |
| 26 | #define ONE_OVER_SQRT120.288675 0.288675 |
| 27 | #define ONE_OVER_120.08333333333333 0.08333333333333 |
| 28 | #define EPS1e-6 1e-6 |
| 29 | |
| 30 | bool static DTrackHitSelector_cdchit_cmp(pair<double,const DCDCTrackHit *>a, |
| 31 | pair<double,const DCDCTrackHit *>b){ |
| 32 | if (a.second->wire->ring!=b.second->wire->ring) |
| 33 | return (a.second->wire->ring>b.second->wire->ring); |
| 34 | return (a.first>b.first); |
| 35 | } |
| 36 | bool static DTrackHitSelector_fdchit_cmp(pair<double,const DFDCPseudo *>a, |
| 37 | pair<double,const DFDCPseudo *>b){ |
| 38 | if (a.second->wire->layer!=b.second->wire->layer) |
| 39 | return (a.second->wire->layer>b.second->wire->layer); |
| 40 | return (a.first>b.first); |
| 41 | } |
| 42 | bool static DTrackHitSelector_gemhit_cmp(pair<double,const DGEMPoint *>a, |
| 43 | pair<double,const DGEMPoint *>b){ |
| 44 | if (a.second->detector!=b.second->detector) |
| 45 | return (a.second->detector>b.second->detector); |
| 46 | return (a.first>b.first); |
| 47 | } |
| 48 | bool static DTrackHitSelector_trdhit_cmp(pair<double,const DTRDPoint *>a, |
| 49 | pair<double,const DTRDPoint *>b){ |
| 50 | if (a.second->detector!=b.second->detector) |
| 51 | return (a.second->detector>b.second->detector); |
| 52 | return (a.first>b.first); |
| 53 | } |
| 54 | |
| 55 | bool static DTrackHitSelector_cdchit_in_cmp(const DCDCTrackHit *a, const DCDCTrackHit *b){ |
| 56 | if (a->wire->ring != b->wire->ring) return (a->wire->ring < b->wire->ring); |
| 57 | return (a->wire->straw < b->wire->straw); |
| 58 | } |
| 59 | |
| 60 | bool static DTrackHitSelector_fdchit_in_cmp(const DFDCPseudo *a, const DFDCPseudo *b){ |
| 61 | if (a->wire->layer != b->wire->layer) return (a->wire->layer < b->wire->layer); |
| 62 | return (a->wire->wire < b->wire->wire); |
| 63 | } |
| 64 | |
| 65 | |
| 66 | //--------------------------------- |
| 67 | // DTrackHitSelectorALT2 (Constructor) |
| 68 | //--------------------------------- |
| 69 | DTrackHitSelectorALT2::DTrackHitSelectorALT2(jana::JEventLoop *loop, int32_t runnumber):DTrackHitSelector(loop) |
| 70 | { |
| 71 | HS_DEBUG_LEVEL = 0; |
| 72 | MAKE_DEBUG_TREES = false; |
| 73 | MIN_HIT_PROB_CDC = 0.01; |
| 74 | MIN_HIT_PROB_FDC = 0.01; |
| 75 | MIN_HIT_PROB_TRD = 0.01; |
| 76 | MIN_HIT_PROB_GEM = 0.01; |
| 77 | MIN_FDC_SIGMA_ANODE_CANDIDATE = 0.1000; |
| 78 | MIN_FDC_SIGMA_CATHODE_CANDIDATE = 0.1000; |
| 79 | MIN_FDC_SIGMA_ANODE_WIREBASED = 0.0100; |
| 80 | MIN_FDC_SIGMA_CATHODE_WIREBASED = 0.0100; |
| 81 | MAX_DOCA=2.5; |
| 82 | |
| 83 | gPARMS->SetDefaultParameter("TRKFIT:MAX_DOCA",MAX_DOCA,"Maximum doca for associating hit with track"); |
| 84 | |
| 85 | gPARMS->SetDefaultParameter("TRKFIT:HS_DEBUG_LEVEL", HS_DEBUG_LEVEL, "Debug verbosity level for hit selector used in track fitting (0=no debug messages)"); |
| 86 | gPARMS->SetDefaultParameter("TRKFIT:MAKE_DEBUG_TREES", MAKE_DEBUG_TREES, "Create a TTree with debugging info on hit selection for the FDC and CDC"); |
| 87 | gPARMS->SetDefaultParameter("TRKFIT:MIN_HIT_PROB_CDC", MIN_HIT_PROB_CDC, "Minimum probability a CDC hit may have to be associated with a track to be included in list passed to fitter"); |
| 88 | gPARMS->SetDefaultParameter("TRKFIT:MIN_HIT_PROB_FDC", MIN_HIT_PROB_FDC, "Minimum probability a FDC hit may have to be associated with a track to be included in list passed to fitter"); |
| 89 | gPARMS->SetDefaultParameter("TRKFIT:MIN_FDC_SIGMA_ANODE_CANDIDATE", MIN_FDC_SIGMA_ANODE_CANDIDATE, "Minimum sigma used for FDC anode hits on track candidates"); |
| 90 | gPARMS->SetDefaultParameter("TRKFIT:MIN_FDC_SIGMA_CATHODE_CANDIDATE", MIN_FDC_SIGMA_CATHODE_CANDIDATE, "Minimum sigma used for FDC cathode hits on track candidates"); |
| 91 | gPARMS->SetDefaultParameter("TRKFIT:MIN_FDC_SIGMA_ANODE_WIREBASED", MIN_FDC_SIGMA_ANODE_WIREBASED, "Minimum sigma used for FDC anode hits on wire-based tracks"); |
| 92 | gPARMS->SetDefaultParameter("TRKFIT:MIN_FDC_SIGMA_CATHODE_WIREBASED", MIN_FDC_SIGMA_CATHODE_WIREBASED, "Minimum sigma used for FDC cathode hits on wire-based tracks"); |
| 93 | |
| 94 | cdchitsel = NULL__null; |
| 95 | fdchitsel = NULL__null; |
| 96 | if(MAKE_DEBUG_TREES){ |
| 97 | loop->GetJApplication()->Lock(); |
| 98 | |
| 99 | cdchitsel= (TTree*)gROOT(ROOT::GetROOT())->FindObject("cdchitsel"); |
| 100 | if(!cdchitsel){ |
| 101 | cdchitsel = new TTree("cdchitsel", "CDC Hit Selector"); |
| 102 | cdchitsel->Branch("H", &cdchitdbg, "fit_type/I:p/F:theta:mass:sigma:x:y:z:s:itheta02:itheta02s:itheta02s2:dist:doca:resi:chisq:prob:sig_phi:sig_lambda:sig_pt"); |
| 103 | }else{ |
| 104 | _DBG__std::cerr<<"libraries/TRACKING/DTrackHitSelectorALT2.cc" <<":"<<104<<std::endl; |
| 105 | jerr<<" !!! WARNING !!!"<<endl; |
| 106 | jerr<<"It appears that the cdchitsel TTree is already defined."<<endl; |
| 107 | jerr<<"This is probably means you are running with multiple threads."<<endl; |
| 108 | jerr<<"To avoid complication, filling of the hit selector debug"<<endl; |
| 109 | jerr<<"trees will be disabled for this thread."<<endl; |
| 110 | _DBG__std::cerr<<"libraries/TRACKING/DTrackHitSelectorALT2.cc" <<":"<<110<<std::endl; |
| 111 | cdchitsel = NULL__null; |
| 112 | } |
| 113 | |
| 114 | fdchitsel= (TTree*)gROOT(ROOT::GetROOT())->FindObject("fdchitsel"); |
| 115 | if(!fdchitsel){ |
| 116 | fdchitsel = new TTree("fdchitsel", "FDC Hit Selector"); |
| 117 | fdchitsel->Branch("H", &fdchitdbg, "fit_type/I:p/F:theta:mass:sigma_anode:sigma_cathode:x:y:z:s:itheta02:itheta02s:itheta02s2:dist:doca:resi:u:u_cathodes:resic:chisq:prob:sig_phi:sig_lambda:sig_pt"); |
| 118 | }else{ |
| 119 | _DBG__std::cerr<<"libraries/TRACKING/DTrackHitSelectorALT2.cc" <<":"<<119<<std::endl; |
| 120 | jerr<<" !!! WARNING !!!"<<endl; |
| 121 | jerr<<"It appears that the fdchitsel TTree is already defined."<<endl; |
| 122 | jerr<<"This is probably means you are running with multiple threads."<<endl; |
| 123 | jerr<<"To avoid complication, filling of the hit selector debug"<<endl; |
| 124 | jerr<<"trees will be disabled for this thread."<<endl; |
| 125 | _DBG__std::cerr<<"libraries/TRACKING/DTrackHitSelectorALT2.cc" <<":"<<125<<std::endl; |
| 126 | fdchitsel = NULL__null; |
| 127 | } |
| 128 | |
| 129 | loop->GetJApplication()->Unlock(); |
| 130 | } |
| 131 | |
| 132 | DApplication* dapp = dynamic_cast<DApplication*>(loop->GetJApplication()); |
| 133 | bfield = dapp->GetBfield(runnumber); // this should be run number based! |
| 134 | |
| 135 | } |
| 136 | |
| 137 | //--------------------------------- |
| 138 | // ~DTrackHitSelectorALT2 (Destructor) |
| 139 | //--------------------------------- |
| 140 | DTrackHitSelectorALT2::~DTrackHitSelectorALT2() |
| 141 | { |
| 142 | |
| 143 | } |
| 144 | |
| 145 | |
| 146 | //--------------------------------- |
| 147 | // GetTRDHits |
| 148 | //--------------------------------- |
| 149 | void DTrackHitSelectorALT2::GetTRDHits(const vector<DTrackFitter::Extrapolation_t> &extrapolations, const vector<const DTRDPoint*> &trdhits_in, vector<const DTRDPoint*> &trdhits_out) const { |
| 150 | // Vector of pairs storing the hit with the probability it is on the track |
| 151 | vector<pair<double,const DTRDPoint*> >trdhits_tmp; |
| 152 | |
| 153 | for (unsigned int k=0;k<extrapolations.size();k++){ |
| 154 | DVector3 pos=extrapolations[k].position; |
| 155 | for (unsigned int j=0;j<trdhits_in.size();j++){ |
| 156 | const DTRDPoint *hit=trdhits_in[j]; |
| 157 | if (int(k)==hit->detector){ |
| 158 | double dx=hit->x-pos.X(); |
| 159 | double dy=hit->y-pos.Y(); |
| 160 | double varx=1.,vary=1.; |
| 161 | double chisq=dx*dx/varx+dy*dy/vary; |
| 162 | // Use chi-sq probability function with Ndof=2 to calculate probability |
| 163 | double probability = TMath::Prob(chisq, 2); |
| 164 | if(probability>=MIN_HIT_PROB_GEM){ |
| 165 | pair<double,const DTRDPoint*>myhit; |
| 166 | myhit.first=probability; |
| 167 | myhit.second=hit; |
| 168 | trdhits_tmp.push_back(myhit); |
| 169 | } |
| 170 | } |
| 171 | } |
| 172 | } |
| 173 | // Order according to layer number and probability,then put the hits in the |
| 174 | // output list with the following algorithm: put hits with the highest |
| 175 | // probability in a given layer in the output list. |
| 176 | sort(trdhits_tmp.begin(),trdhits_tmp.end(),DTrackHitSelector_trdhit_cmp); |
| 177 | int old_layer=1000; |
| 178 | for (unsigned int i=0;i<trdhits_tmp.size();i++){ |
| 179 | if (trdhits_tmp[i].second->detector!=old_layer){ |
| 180 | trdhits_out.push_back(trdhits_tmp[i].second); |
| 181 | } |
| 182 | old_layer=trdhits_tmp[i].second->detector; |
| 183 | } |
| 184 | } |
| 185 | |
| 186 | //--------------------------------- |
| 187 | // GetGEMHits |
| 188 | //--------------------------------- |
| 189 | void DTrackHitSelectorALT2::GetGEMHits(const vector<DTrackFitter::Extrapolation_t> &extrapolations, const vector<const DGEMPoint*> &gemhits_in, vector<const DGEMPoint*> &gemhits_out) const { |
| 190 | // Vector of pairs storing the hit with the probability it is on the track |
| 191 | vector<pair<double,const DGEMPoint*> >gemhits_tmp; |
| 192 | |
| 193 | for (unsigned int k=0;k<extrapolations.size();k++){ |
| 194 | DVector3 pos=extrapolations[k].position; |
| 195 | for (unsigned int j=0;j<gemhits_in.size();j++){ |
| 196 | const DGEMPoint *hit=gemhits_in[j]; |
| 197 | if (int(k)==hit->detector){ |
| 198 | double dx=hit->x-pos.X(); |
| 199 | double dy=hit->y-pos.Y(); |
| 200 | double varx=1.,vary=1.; |
| 201 | double chisq=dx*dx/varx+dy*dy/vary; |
| 202 | // Use chi-sq probability function with Ndof=2 to calculate probability |
| 203 | double probability = TMath::Prob(chisq, 2); |
| 204 | if(probability>=MIN_HIT_PROB_GEM){ |
| 205 | pair<double,const DGEMPoint*>myhit; |
| 206 | myhit.first=probability; |
| 207 | myhit.second=hit; |
| 208 | gemhits_tmp.push_back(myhit); |
| 209 | } |
| 210 | } |
| 211 | } |
| 212 | |
| 213 | } |
| 214 | // Order according to layer number and probability,then put the hits in the |
| 215 | // output list with the following algorithm: put hits with the highest |
| 216 | // probability in a given layer in the output list. |
| 217 | sort(gemhits_tmp.begin(),gemhits_tmp.end(),DTrackHitSelector_gemhit_cmp); |
| 218 | int old_layer=1000; |
| 219 | for (unsigned int i=0;i<gemhits_tmp.size();i++){ |
| 220 | if (gemhits_tmp[i].second->detector!=old_layer){ |
| 221 | gemhits_out.push_back(gemhits_tmp[i].second); |
| 222 | } |
| 223 | old_layer=gemhits_tmp[i].second->detector; |
| 224 | } |
| 225 | |
| 226 | } |
| 227 | |
| 228 | //--------------------------------- |
| 229 | // GetCDCHits |
| 230 | //--------------------------------- |
| 231 | void DTrackHitSelectorALT2::GetCDCHits(double Bz,double q, |
| 232 | const vector<DTrackFitter::Extrapolation_t> &extrapolations, const vector<const DCDCTrackHit*> &cdchits_in, vector<const DCDCTrackHit*> &cdchits_out, int N) const |
| 233 | { |
| 234 | // Vector of pairs storing the hit with the probability it is on the track |
| 235 | vector<pair<double,const DCDCTrackHit*> >cdchits_tmp; |
| 236 | |
| 237 | // Sort so innermost ring is first and outermost is last |
| 238 | vector<const DCDCTrackHit*> cdchits_in_sorted = cdchits_in; |
| 239 | sort(cdchits_in_sorted.begin(),cdchits_in_sorted.end(),DTrackHitSelector_cdchit_in_cmp); |
| 240 | |
| 241 | // The variance on the residual due to measurement error. |
| 242 | double var=0.25*2.4336*ONE_OVER_120.08333333333333; |
| 243 | |
| 244 | // To estimate the impact of errors in the track momentum on the variance |
| 245 | // of the residual, use a helical approximation. |
| 246 | double a=-0.003*Bz*q; |
| 247 | DVector3 mom=extrapolations[extrapolations.size()/2].momentum; |
| 248 | double p=mom.Mag(); |
| 249 | double p_over_a=p/a; |
| 250 | double a_over_p=1./p_over_a; |
| 251 | double lambda=M_PI_21.57079632679489661923-mom.Theta(); |
| 252 | double tanl=tan(lambda); |
| 253 | double tanl2=tanl*tanl; |
| 254 | double sinl=sin(lambda); |
| 255 | double cosl=cos(lambda); |
| 256 | double cosl2=cosl*cosl; |
| 257 | double sinl2=sinl*sinl; |
| 258 | double pt_over_a=cosl*p_over_a; |
| 259 | |
| 260 | // variances |
| 261 | double var_lambda_res=0.; |
| 262 | double var_phi_radial=0.; |
| 263 | double var_x0=0.0,var_y0=0.0; |
| 264 | double var_k=0.; |
| 265 | |
| 266 | // Loop over all the CDC hits looking for matches with the track |
| 267 | bool outermost_hit=true; |
| 268 | vector<const DCDCTrackHit*>::const_reverse_iterator iter; |
| 269 | for(iter=cdchits_in_sorted.rbegin(); iter!=cdchits_in_sorted.rend(); iter++){ |
| 270 | const DCDCTrackHit *hit = *iter; |
| 271 | DVector3 origin=hit->wire->origin; |
| 272 | DVector3 dir=hit->wire->udir; |
| 273 | double ux=dir.x(); |
| 274 | double uy=dir.y(); |
| 275 | double uz=dir.z(); |
| 276 | double z0=origin.z(); |
| 277 | double d2=0.,d2_old=1.e6; |
| 278 | double s=0.; |
| 279 | DVector3 old_trackpos; |
| 280 | for (unsigned int i=0;i<extrapolations.size();i++){ |
| 281 | DVector3 trackpos=extrapolations[i].position; |
| 282 | double dz=trackpos.z()-z0; |
| 283 | DVector3 wirepos=origin+dz/uz*dir; |
| 284 | d2=(trackpos-wirepos).Perp2(); |
| 285 | if (d2>d2_old){ |
| 286 | if (d2<4){ // has to be reasonably close to the wire in question |
| 287 | double phi=extrapolations[i].momentum.Phi(); |
| 288 | double sinphi=sin(phi); |
| 289 | double cosphi=cos(phi); |
| 290 | // Variables relating wire direction and track direction |
| 291 | double my_ux=ux*sinl-cosl*cosphi; |
| 292 | double my_uy=uy*sinl-cosl*sinphi; |
| 293 | double denom=my_ux*my_ux+my_uy*my_uy; |
| 294 | // For simplicity make a linear approximation for the path increment |
| 295 | double ds=((old_trackpos.x()-origin.x()-ux*dz)*my_ux |
| 296 | +(old_trackpos.y()-origin.y()-uy*dz)*my_uy)/denom; |
| 297 | wirepos+=ds*dir; |
| 298 | double dx=old_trackpos.x()+mom.X()/mom.Mag()*ds-wirepos.x(); |
| 299 | double dy=old_trackpos.y()+mom.Y()/mom.Mag()*ds-wirepos.y(); |
| 300 | d2_old=d2; |
| 301 | d2=dx*dx+dy*dy; |
| 302 | if (d2>d2_old){ |
| 303 | // linear approximation did not work... |
| 304 | d2=d2_old; |
| 305 | } |
| 306 | // Variance in dip angle due to multiple scattering |
| 307 | double var_lambda = extrapolations[i].theta2ms_sum/3.; |
| 308 | // variance in phi due to multiple scattering |
| 309 | double var_phi=var_lambda*(1.+tanl2); |
| 310 | |
| 311 | if (outermost_hit){ |
| 312 | // variance in the curvature k due to resolution and multiple scattering |
| 313 | double s_sq=s*s; |
| 314 | double s_4=s_sq*s_sq; |
| 315 | double sum_s_theta_ms=extrapolations[i].s_theta_ms_sum; |
| 316 | double var_k_ms=(4./3.)*sum_s_theta_ms*sum_s_theta_ms/s_4; |
| 317 | double var_k_res=var*0.0720/double(N+4)/s_4/(cosl2*cosl2); |
| 318 | var_k=var_k_ms+var_k_res; |
| 319 | |
| 320 | // Variance in dip angle due to measurement error |
| 321 | var_lambda_res=12.0*var*double(N-1)/double(N*(N+1)) |
| 322 | *sinl2*sinl2/s_sq; |
| 323 | |
| 324 | // Variance in phi, crude approximation |
| 325 | double tan_s_2R=tan(s*cosl/(2.*pt_over_a)); |
| 326 | var_phi_radial=tan_s_2R*tan_s_2R*pt_over_a*pt_over_a*var_k; |
| 327 | |
| 328 | outermost_hit=false; |
| 329 | } |
| 330 | // Include error in lambda due to measurements |
| 331 | var_lambda+=var_lambda_res; |
| 332 | |
| 333 | // Include uncertainty in phi due to uncertainty in the center of |
| 334 | // the circle |
| 335 | var_phi+=var_phi_radial; |
| 336 | |
| 337 | // Variance in position due to multiple scattering |
| 338 | double var_pos_ms=extrapolations[i].s_theta_ms_sum/3.; |
| 339 | |
| 340 | // Hit doca and residual |
| 341 | double doca=sqrt(d2); |
| 342 | double dist=0.39; |
| 343 | double resi=dist-doca; |
| 344 | |
| 345 | // Variances in x and y due to uncertainty in track parameters |
| 346 | double as_over_p=s*a_over_p; |
| 347 | double sin_as_over_p=sin(as_over_p); |
| 348 | double cos_as_over_p=cos(as_over_p); |
| 349 | double one_minus_cos_as_over_p=1-cos_as_over_p; |
| 350 | double diff1=sin_as_over_p-as_over_p*cos_as_over_p; |
| 351 | double diff2=one_minus_cos_as_over_p-as_over_p*sin_as_over_p; |
| 352 | double dx_dk=2.*pt_over_a*pt_over_a*(sinphi*diff2-cosphi*diff1); |
| 353 | double dx_dcosl=p_over_a*(cosphi*sin_as_over_p-sinphi*one_minus_cos_as_over_p); |
| 354 | double dx_dphi=-pt_over_a*(sinphi*sin_as_over_p+cosphi*one_minus_cos_as_over_p); |
| 355 | double var_x=var_x0+dx_dk*dx_dk*var_k+var_pos_ms |
| 356 | +dx_dcosl*dx_dcosl*sinl*sinl*var_lambda+dx_dphi*dx_dphi*var_phi; |
| 357 | |
| 358 | double dy_dk=-2.*pt_over_a*pt_over_a*(sinphi*diff1+cosphi*diff2); |
| 359 | double dy_dcosl=p_over_a*(sinphi*sin_as_over_p+cosphi*one_minus_cos_as_over_p); |
| 360 | double dy_dphi=pt_over_a*(cosphi*sin_as_over_p-sinphi*one_minus_cos_as_over_p); |
| 361 | double var_y=var_y0+dy_dk*dy_dk*var_k+var_pos_ms |
| 362 | +dy_dcosl*dy_dcosl*sinl*sinl*var_lambda+dy_dphi*dy_dphi*var_phi; |
| 363 | |
| 364 | double dd_dx=dx/doca; |
| 365 | double dd_dy=dy/doca; |
| 366 | double var_d=dd_dx*dd_dx*var_x+dd_dy*dd_dy*var_y; |
| 367 | double chisq=resi*resi/(var+var_d); |
| 368 | |
| 369 | // Use chi-sq probability function with Ndof=1 to calculate probability |
| 370 | double probability = TMath::Prob(chisq, 1); |
| 371 | if(probability>=MIN_HIT_PROB_CDC){ |
| 372 | pair<double,const DCDCTrackHit*>myhit; |
| 373 | myhit.first=probability; |
| 374 | myhit.second=hit; |
| 375 | cdchits_tmp.push_back(myhit); |
| 376 | } |
| 377 | |
| 378 | // Optionally fill debug tree |
| 379 | if(cdchitsel){ |
| 380 | cdchitdbg.fit_type = kWireBased; |
| 381 | cdchitdbg.p = p; |
| 382 | cdchitdbg.theta = extrapolations[i].momentum.Theta(); |
| 383 | //cdchitdbg.mass = mass; |
| 384 | cdchitdbg.sigma = sqrt(var+var_d); |
| 385 | cdchitdbg.x = old_trackpos.X(); |
| 386 | cdchitdbg.y = old_trackpos.Y(); |
| 387 | cdchitdbg.z = old_trackpos.Z(); |
| 388 | cdchitdbg.s = s; |
| 389 | cdchitdbg.itheta02 = extrapolations[i].theta2ms_sum; |
| 390 | //cdchitdbg.itheta02s = last_step->itheta02s; |
| 391 | //cdchitdbg.itheta02s2 = last_step->itheta02s2; |
| 392 | cdchitdbg.dist = dist; |
| 393 | cdchitdbg.doca = doca; |
| 394 | cdchitdbg.resi = resi; |
| 395 | cdchitdbg.chisq = chisq; |
| 396 | cdchitdbg.prob = probability; |
| 397 | cdchitdbg.sig_phi=sqrt(var_phi); |
| 398 | cdchitdbg.sig_lambda=sqrt(var_lambda); |
| 399 | cdchitdbg.sig_pt=fabs(2.*pt_over_a)*sqrt(var_k); |
| 400 | |
| 401 | cdchitsel->Fill(); |
| 402 | } |
| 403 | |
| 404 | if(HS_DEBUG_LEVEL>10){ |
| 405 | _DBG_std::cerr<<"libraries/TRACKING/DTrackHitSelectorALT2.cc" <<":"<<405<<" "; |
| 406 | if (probability>=MIN_HIT_PROB_CDC) jerr<<"---> "; |
| 407 | jerr<<"s="<<s<<" t=" <<hit->tdrift << " doca="<<doca<<" dist="<<dist<<" resi="<<resi<<" sigma="/*<<sigma_total*/<<" prob="<<probability<<endl; |
| 408 | } |
| 409 | } |
| 410 | break; |
| 411 | } |
| 412 | d2_old=d2; |
| 413 | old_trackpos=trackpos; |
| 414 | s=extrapolations[i].s; |
| 415 | } |
| 416 | |
| 417 | } |
| 418 | |
| 419 | // Order according to ring number and probability, then put the hits in the |
| 420 | // output list with the following algorithm: hits with the highest |
| 421 | // probability in a given ring are automatically put in the output list, |
| 422 | // but if there is more than one hit in a given ring, only those hits |
| 423 | // that are within +/-1 of the straw # of the most probable hit are added |
| 424 | // to the list. |
| 425 | sort(cdchits_tmp.begin(),cdchits_tmp.end(),DTrackHitSelector_cdchit_cmp); |
| 426 | int old_straw=1000,old_ring=1000; |
| 427 | for (unsigned int i=0;i<cdchits_tmp.size();i++){ |
| 428 | if (cdchits_tmp[i].second->wire->ring!=old_ring || |
| 429 | abs(cdchits_tmp[i].second->wire->straw-old_straw)==1){ |
| 430 | cdchits_out.push_back(cdchits_tmp[i].second); |
| 431 | } |
| 432 | old_straw=cdchits_tmp[i].second->wire->straw; |
| 433 | old_ring=cdchits_tmp[i].second->wire->ring; |
| 434 | } |
| 435 | } |
| 436 | |
| 437 | //--------------------------------- |
| 438 | // GetCDCHits |
| 439 | //--------------------------------- |
| 440 | void DTrackHitSelectorALT2::GetCDCHits(fit_type_t fit_type, const DReferenceTrajectory *rt, const vector<const DCDCTrackHit*> &cdchits_in, vector<const DCDCTrackHit*> &cdchits_out, int N) const |
| 441 | { |
| 442 | // Vector of pairs storing the hit with the probability it is on the track |
| 443 | vector<pair<double,const DCDCTrackHit*> >cdchits_tmp; |
| 444 | |
| 445 | /// Determine the probability that for each CDC hit that it came from the |
| 446 | /// track with the given trajectory. |
| 447 | /// |
| 448 | /// This will calculate a probability for each CDC hit that |
| 449 | /// it came from the track represented by the given |
| 450 | /// DReference trajectory. The probability is based on |
| 451 | /// the residual between the distance of closest approach |
| 452 | /// of the trajectory to the wire and the drift time for |
| 453 | /// time-based tracks and the distance to the wire for |
| 454 | /// wire-based tracks. |
| 455 | |
| 456 | // Sort so innermost ring is first and outermost is last |
| 457 | vector<const DCDCTrackHit*> cdchits_in_sorted = cdchits_in; |
| 458 | sort(cdchits_in_sorted.begin(),cdchits_in_sorted.end(),DTrackHitSelector_cdchit_in_cmp); |
| 459 | |
| 460 | // Calculate beta of particle. |
| 461 | //double my_mass=rt->GetMass(); |
| 462 | // double one_over_beta =sqrt(1.0+my_mass*my_mass/rt->swim_steps[0].mom.Mag2()); |
| 463 | |
| 464 | // The variance on the residual due to measurement error. |
| 465 | double var=0.25*2.4336*ONE_OVER_120.08333333333333; |
| 466 | |
| 467 | // To estimate the impact of errors in the track momentum on the variance of the residual, |
| 468 | // use a helical approximation. |
| 469 | const DReferenceTrajectory::swim_step_t *my_step=&rt->swim_steps[0]; |
| 470 | double Bz0=my_step->B.z(); |
| 471 | double a=-0.003*Bz0*rt->q; |
| 472 | double p=my_step->mom.Mag(); |
| 473 | double p_over_a=p/a; |
| 474 | double a_over_p=1./p_over_a; |
| 475 | double lambda=M_PI_21.57079632679489661923-my_step->mom.Theta(); |
| 476 | double cosl=cos(lambda); |
| 477 | double sinl=sin(lambda); |
| 478 | double sinl2=sinl*sinl; |
| 479 | double cosl2=cosl*cosl; |
| 480 | double tanl=tan(lambda); |
| 481 | double tanl2=tanl*tanl; |
| 482 | double pt_over_a=cosl*p_over_a; |
| 483 | double phi=my_step->mom.Phi(); |
| 484 | double cosphi=cos(phi); |
| 485 | double sinphi=sin(phi); |
| 486 | |
| 487 | double var_lambda_res=0.; |
| 488 | double var_phi_radial=0.; |
| 489 | double var_x0=0.0,var_y0=0.0; |
| 490 | double mass=rt->GetMass(); |
| 491 | double one_over_beta=sqrt(1.+mass*mass/(p*p)); |
| 492 | double var_pt_factor=0.016*one_over_beta/(cosl*a); |
| 493 | double var_pt_over_pt_sq=0.; |
| 494 | |
| 495 | // Keep track of straws and rings |
| 496 | int old_straw=1000,old_ring=1000; |
| 497 | |
| 498 | // Loop over hits |
| 499 | bool outermost_hit=true; |
| 500 | vector<const DCDCTrackHit*>::const_reverse_iterator iter; |
| 501 | for(iter=cdchits_in_sorted.rbegin(); iter!=cdchits_in_sorted.rend(); iter++){ |
| 502 | const DCDCTrackHit *hit = *iter; |
| 503 | |
| 504 | // Skip hit if it is on the same wire as the previous hit |
| 505 | if (hit->wire->ring == old_ring && hit->wire->straw==old_straw){ |
| 506 | //_DBG_ << "ring " << hit->wire->ring << " straw " << hit->wire->straw << endl; |
| 507 | continue; |
| 508 | } |
| 509 | old_ring=hit->wire->ring; |
| 510 | old_straw=hit->wire->straw; |
| 511 | |
| 512 | // Find the DOCA to this wire |
| 513 | double s=0.; |
| 514 | double doca = rt->DistToRT(hit->wire, &s); |
| 515 | |
| 516 | if(!isfinite(doca)) continue; |
| 517 | if(!isfinite(s))continue; |
| 518 | if (s<1.) continue; |
| 519 | if (doca>MAX_DOCA) continue; |
| 520 | |
| 521 | const DReferenceTrajectory::swim_step_t *last_step = rt->GetLastSwimStep(); |
| 522 | |
| 523 | // Position along trajectory |
| 524 | DVector3 pos=rt->GetLastDOCAPoint(); |
| 525 | |
| 526 | // Compensate for the fact that the field at the "origin" of the |
| 527 | // track does not correspond to the average Bz used to compute pt |
| 528 | double Bratio=last_step->B.z()/Bz0; |
| 529 | double invBratio=1./Bratio; |
| 530 | pt_over_a*=invBratio; |
| 531 | p_over_a*=invBratio; |
| 532 | a_over_p*=Bratio; |
| 533 | |
| 534 | // Variances in angles due to multiple scattering |
| 535 | double var_lambda = last_step->itheta02/3.; |
| 536 | double var_phi=var_lambda*(1.+tanl2); |
| 537 | |
| 538 | if (outermost_hit){ |
| 539 | // Fractional variance in the curvature k due to resolution and multiple scattering |
| 540 | double s_sq=s*s; |
| 541 | double var_k_over_k_sq_res=var*p_over_a*p_over_a*0.0720/double(N+4)/(s_sq*s_sq*sinl2)/cosl2; |
| 542 | double var_k_over_k_sq_ms=var_pt_factor*var_pt_factor*last_step->invX0/s; |
| 543 | // Fractional variance in pt |
| 544 | var_pt_over_pt_sq=var_k_over_k_sq_ms+var_k_over_k_sq_res; |
| 545 | |
| 546 | // Variance in dip angle due to measurement error |
| 547 | var_lambda_res=12.0*var*double(N-1)/double(N*(N+1))*cosl2*cosl2/s_sq; |
| 548 | |
| 549 | // Variance in phi, crude approximation |
| 550 | double tan_s_2R=tan(s*cosl/(2.*pt_over_a)); |
| 551 | var_phi_radial=tan_s_2R*tan_s_2R*var_pt_over_pt_sq; |
| 552 | |
| 553 | outermost_hit=false; |
| 554 | } |
| 555 | // Include error in lambda due to measurements |
| 556 | var_lambda+=var_lambda_res; |
| 557 | |
| 558 | // Include uncertainty in phi due to uncertainty in the radius of |
| 559 | // the circle |
| 560 | var_phi+=var_phi_radial; |
| 561 | |
| 562 | // Get "measured" distance to wire. |
| 563 | double dist=0.39; |
| 564 | |
| 565 | // Residual |
| 566 | double resi = dist - doca; |
| 567 | |
| 568 | // Variance in position due to multiple scattering |
| 569 | double var_pos_ms=last_step->itheta02s2/3.; |
| 570 | |
| 571 | // Variances in x and y due to uncertainty in track parameters |
| 572 | double as_over_p=s*a_over_p; |
| 573 | double sin_as_over_p=sin(as_over_p); |
| 574 | double cos_as_over_p=cos(as_over_p); |
| 575 | double one_minus_cos_as_over_p=1-cos_as_over_p; |
| 576 | double diff1=sin_as_over_p-as_over_p*cos_as_over_p; |
| 577 | double diff2=one_minus_cos_as_over_p-as_over_p*sin_as_over_p; |
| 578 | double pdx_dp=pt_over_a*(cosphi*diff1-sinphi*diff2); |
| 579 | double dx_dcosl=p_over_a*(cosphi*sin_as_over_p-sinphi*one_minus_cos_as_over_p); |
| 580 | double dx_dphi=-pt_over_a*(sinphi*sin_as_over_p+cosphi*one_minus_cos_as_over_p); |
| 581 | double var_x=var_x0+pdx_dp*pdx_dp*var_pt_over_pt_sq+var_pos_ms |
| 582 | +dx_dcosl*dx_dcosl*sinl*sinl*var_lambda+dx_dphi*dx_dphi*var_phi; |
| 583 | |
| 584 | double pdy_dp=pt_over_a*(sinphi*diff1+cosphi*diff2); |
| 585 | double dy_dcosl=p_over_a*(sinphi*sin_as_over_p+cosphi*one_minus_cos_as_over_p); |
| 586 | double dy_dphi=pt_over_a*(cosphi*sin_as_over_p-sinphi*one_minus_cos_as_over_p); |
| 587 | double var_y=var_y0+pdy_dp*pdy_dp*var_pt_over_pt_sq+var_pos_ms |
| 588 | +dy_dcosl*dy_dcosl*sinl*sinl*var_lambda+dy_dphi*dy_dphi*var_phi; |
| 589 | |
| 590 | |
| 591 | DVector3 origin=hit->wire->origin; |
| 592 | DVector3 dir=hit->wire->udir; |
| 593 | double uz=dir.z(); |
| 594 | double z0=origin.z(); |
| 595 | DVector3 wirepos=origin+(pos.z()-z0)/uz*dir; |
| 596 | double dd_dx=(pos.x()-wirepos.x())/doca; |
| 597 | double dd_dy=(pos.y()-wirepos.y())/doca; |
| 598 | double var_d=dd_dx*dd_dx*var_x+dd_dy*dd_dy*var_y; |
| 599 | |
| 600 | double chisq=resi*resi/(var+var_d); |
| 601 | |
| 602 | // Use chi-sq probability function with Ndof=1 to calculate probability |
| 603 | double probability = TMath::Prob(chisq, 1); |
| 604 | if(probability>=MIN_HIT_PROB_CDC){ |
| 605 | pair<double,const DCDCTrackHit*>myhit; |
| 606 | myhit.first=probability; |
| 607 | myhit.second=hit; |
| 608 | cdchits_tmp.push_back(myhit); |
| 609 | } |
| 610 | |
| 611 | // Optionally fill debug tree |
| 612 | if(cdchitsel){ |
| 613 | cdchitdbg.fit_type = fit_type; |
| 614 | cdchitdbg.p = p; |
| 615 | cdchitdbg.theta = rt->swim_steps[0].mom.Theta(); |
| 616 | cdchitdbg.mass = mass; |
| 617 | cdchitdbg.sigma = sqrt(var+var_d); |
| 618 | cdchitdbg.x = pos.X(); |
| 619 | cdchitdbg.y = pos.Y(); |
| 620 | cdchitdbg.z = pos.Z(); |
| 621 | cdchitdbg.s = s; |
| 622 | cdchitdbg.itheta02 = last_step->itheta02; |
| 623 | cdchitdbg.itheta02s = last_step->itheta02s; |
| 624 | cdchitdbg.itheta02s2 = last_step->itheta02s2; |
| 625 | cdchitdbg.dist = dist; |
| 626 | cdchitdbg.doca = doca; |
| 627 | cdchitdbg.resi = resi; |
| 628 | cdchitdbg.chisq = chisq; |
| 629 | cdchitdbg.prob = probability; |
| 630 | cdchitdbg.sig_phi=sqrt(var_phi); |
| 631 | cdchitdbg.sig_lambda=sqrt(var_lambda); |
| 632 | cdchitdbg.sig_pt=sqrt(var_pt_over_pt_sq); |
| 633 | |
| 634 | cdchitsel->Fill(); |
| 635 | } |
| 636 | |
| 637 | if(HS_DEBUG_LEVEL>10){ |
| 638 | _DBG_std::cerr<<"libraries/TRACKING/DTrackHitSelectorALT2.cc" <<":"<<638<<" "; |
| 639 | if (probability>=MIN_HIT_PROB_CDC) jerr<<"---> "; |
| 640 | jerr<<"s="<<s<<" t=" <<hit->tdrift << " doca="<<doca<<" dist="<<dist<<" resi="<<resi<<" sigma="/*<<sigma_total*/<<" prob="<<probability<<endl; |
| 641 | } |
| 642 | } |
| 643 | |
| 644 | // Order according to ring number and probability, then put the hits in the |
| 645 | // output list with the following algorithm: hits with the highest |
| 646 | // probability in a given ring are automatically put in the output list, |
| 647 | // but if there is more than one hit in a given ring, only those hits |
| 648 | // that are within +/-1 of the straw # of the most probable hit are added |
| 649 | // to the list. |
| 650 | sort(cdchits_tmp.begin(),cdchits_tmp.end(),DTrackHitSelector_cdchit_cmp); |
| 651 | old_straw=1000,old_ring=1000; |
| 652 | for (unsigned int i=0;i<cdchits_tmp.size();i++){ |
| 653 | if (cdchits_tmp[i].second->wire->ring!=old_ring || |
| 654 | abs(cdchits_tmp[i].second->wire->straw-old_straw)==1){ |
| 655 | cdchits_out.push_back(cdchits_tmp[i].second); |
| 656 | } |
| 657 | old_straw=cdchits_tmp[i].second->wire->straw; |
| 658 | old_ring=cdchits_tmp[i].second->wire->ring; |
| 659 | } |
| 660 | } |
| 661 | |
| 662 | //--------------------------------- |
| 663 | // GetFDCHits |
| 664 | //--------------------------------- |
| 665 | void DTrackHitSelectorALT2::GetFDCHits(fit_type_t fit_type, const DReferenceTrajectory *rt, const vector<const DFDCPseudo*> &fdchits_in, vector<const DFDCPseudo*> &fdchits_out,int N) const |
| 666 | { |
| 667 | // Vector of pairs storing the hit with the probability it is on the track |
| 668 | vector<pair<double,const DFDCPseudo*> >fdchits_tmp; |
| 669 | |
| 670 | /// Determine the probability that for each FDC hit that it came from the |
| 671 | /// track with the given trajectory. |
| 672 | /// |
| 673 | /// This will calculate a probability for each FDC hit that |
| 674 | /// it came from the track represented by the given |
| 675 | /// DReference trajectory. The probability is based on |
| 676 | /// the residual between the distance of closest approach |
| 677 | /// of the trajectory to the wire and the drift distance |
| 678 | /// and the distance along the wire. |
| 679 | |
| 680 | // Sort so innermost ring is first and outermost is last |
| 681 | vector<const DFDCPseudo*> fdchits_in_sorted = fdchits_in; |
| 682 | sort(fdchits_in_sorted.begin(),fdchits_in_sorted.end(),DTrackHitSelector_fdchit_in_cmp); |
| 683 | |
| 684 | // The variance on the residual due to measurement error. |
| 685 | double var_anode = 0.25*ONE_OVER_120.08333333333333; |
| 686 | const double VAR_CATHODE_STRIPS=0.000225; |
| 687 | double var_cathode = VAR_CATHODE_STRIPS; |
| 688 | |
| 689 | // To estimate the impact of errors in the track momentum on the variance of the residual, |
| 690 | // use a helical approximation. |
| 691 | const DReferenceTrajectory::swim_step_t *my_step=&rt->swim_steps[0]; |
| 692 | double Bz0=my_step->B.z(); |
| 693 | double z0=my_step->origin.z(); |
| 694 | double a=-0.003*Bz0*rt->q; |
| 695 | double p=my_step->mom.Mag(); |
| 696 | double p_sq=p*p; |
| 697 | double p_over_a=p/a; |
| 698 | double a_over_p=1./p_over_a; |
| 699 | double lambda=M_PI_21.57079632679489661923-my_step->mom.Theta(); |
| 700 | double cosl=cos(lambda); |
| 701 | double cosl2=cosl*cosl; |
| 702 | double sinl=sin(lambda); |
| 703 | double sinl2=sinl*sinl; |
| 704 | double tanl=tan(lambda); |
| 705 | double tanl2=tanl*tanl; |
| 706 | double pt_over_a=cosl*p_over_a; |
| 707 | double phi=my_step->mom.Phi(); |
| 708 | double cosphi=cos(phi); |
| 709 | double sinphi=sin(phi); |
| 710 | double var_lambda=0.,var_phi=0.,var_lambda_res=0.; |
| 711 | double var_phi_radial=0.; |
| 712 | double mass=rt->GetMass(); |
| 713 | |
| 714 | double var_x0=0.0,var_y0=0.0; |
| 715 | double var_pt_over_pt_sq=0.; |
| 716 | |
| 717 | // Loop over hits |
| 718 | bool most_downstream_hit=true; |
| 719 | vector<const DFDCPseudo*>::const_reverse_iterator iter; |
| 720 | for(iter=fdchits_in_sorted.rbegin(); iter!=fdchits_in_sorted.rend(); iter++){ |
| 721 | const DFDCPseudo *hit = *iter; |
| 722 | |
| 723 | // Find the DOCA to this wire |
| 724 | double s=0.; |
| 725 | double doca = rt->DistToRT(hit->wire, &s); |
| 726 | |
| 727 | if(!isfinite(doca)) continue; |
| 728 | if(!isfinite(s))continue; |
| 729 | if (s<1.) continue; |
| 730 | if (doca>MAX_DOCA)continue; |
| 731 | |
| 732 | const DReferenceTrajectory::swim_step_t *last_step = rt->GetLastSwimStep(); |
| 733 | |
| 734 | // Position along trajectory |
| 735 | DVector3 pos=rt->GetLastDOCAPoint(); |
| 736 | double dz=pos.z()-z0; |
| 737 | |
| 738 | // Direction variables for wire |
| 739 | double cosa=hit->wire->udir.y(); |
| 740 | double sina=hit->wire->udir.x(); |
| 741 | |
| 742 | // Cathode Residual |
| 743 | double u=rt->GetLastDistAlongWire(); |
| 744 | double u_cathodes = hit->s; |
| 745 | double resic = u - u_cathodes; |
| 746 | |
| 747 | // Get "measured" distance to wire. |
| 748 | double dist=0.25; |
| 749 | |
| 750 | // Take into account non-normal incidence to FDC plane |
| 751 | double pz=last_step->mom.z(); |
| 752 | double tx=last_step->mom.x()/pz; |
| 753 | double ty=last_step->mom.y()/pz; |
| 754 | double tu=tx*cosa-ty*sina; |
| 755 | double alpha=atan(tu); |
| 756 | double cosalpha=cos(alpha); |
| 757 | |
| 758 | // Compensate for the fact that the field at the "origin" of the |
| 759 | // track does not correspond to the average Bz used to compute pt |
| 760 | double Bz=last_step->B.z(); |
| 761 | double Bratio=Bz/Bz0; |
| 762 | double invBratio=1./Bratio; |
| 763 | pt_over_a*=invBratio; |
| 764 | p_over_a*=invBratio; |
| 765 | a_over_p*=Bratio; |
| 766 | |
| 767 | // Anode Residual |
| 768 | double resi = dist - doca/cosalpha; |
| 769 | |
| 770 | // Initialize some probability-related variables needed later |
| 771 | double probability=0.,chisq=0.; |
| 772 | |
| 773 | if (fit_type!=kHelical){ |
| 774 | // Correct for deflection of avalanche position due to Lorentz force |
| 775 | double sign=(pos.x()*cosa-pos.y()*sina-hit->w)<0?1:-1.; |
| 776 | double ucor=0.1458*Bz*(1.-0.048*last_step->B.Perp())*sign*doca; |
| 777 | resic-=ucor; |
| 778 | } |
| 779 | else{ |
| 780 | // Cathode variance due to Lorentz deflection |
| 781 | double max_deflection=0.1458*Bz*(1.-0.048*last_step->B.Perp())*0.5; |
| 782 | var_cathode=VAR_CATHODE_STRIPS+max_deflection*max_deflection/12.; |
| 783 | } |
| 784 | double var_tot=var_anode+var_cathode; |
| 785 | |
| 786 | // Variance in angles due to multiple scattering |
| 787 | var_lambda = last_step->itheta02/3.; |
| 788 | var_phi=var_lambda*(1.+tanl2); |
| 789 | |
| 790 | if (most_downstream_hit){ |
| 791 | // Fractional variance in the curvature k due to resolution and multiple scattering |
| 792 | double s_sq=s*s; |
| 793 | double var_k_over_k_sq_res=var_tot*p_over_a*p_over_a |
| 794 | *0.0720/double(N+4)/(s_sq*s_sq)/cosl2; |
| 795 | |
| 796 | |
| 797 | double one_over_beta=sqrt(1.+mass*mass/p_sq); |
| 798 | double var_pt_factor=0.016*one_over_beta/(cosl*0.003*last_step->B.z()); |
| 799 | double var_k_over_k_sq_ms=var_pt_factor*var_pt_factor*last_step->invX0/s; |
| 800 | // Fractional variance in pt |
| 801 | var_pt_over_pt_sq=var_k_over_k_sq_ms+var_k_over_k_sq_res; |
| 802 | |
| 803 | // Variance in dip angle due to measurement error |
| 804 | var_lambda_res=12.0*var_tot*double(N-1)/double(N*(N+1)) |
| 805 | *sinl2*sinl2/s_sq; |
| 806 | |
| 807 | // Variance in phi, crude approximation |
| 808 | double tan_s_2R=tan(s*cosl/(2.*pt_over_a)); |
| 809 | var_phi_radial=tan_s_2R*tan_s_2R*var_pt_over_pt_sq; |
| 810 | |
| 811 | most_downstream_hit=false; |
| 812 | } |
| 813 | |
| 814 | // Include error in lambda due to measurements |
| 815 | var_lambda+=var_lambda_res; |
Value stored to 'var_lambda' is never read | |
| 816 | |
| 817 | // Include uncertainty in phi due to uncertainty in the center of |
| 818 | // the circle |
| 819 | var_phi+=var_phi_radial; |
| 820 | |
| 821 | var_phi=0.; |
| 822 | var_lambda=0.; |
| 823 | |
| 824 | // Variance in position due to multiple scattering |
| 825 | double var_pos_ms=last_step->itheta02s2/3.; |
| 826 | |
| 827 | // Variances in x and y due to uncertainty in track parameters |
| 828 | double as_over_p=s*a_over_p; |
| 829 | double sin_as_over_p=sin(as_over_p); |
| 830 | double cos_as_over_p=cos(as_over_p); |
| 831 | double one_minus_cos_as_over_p=1-cos_as_over_p; |
| 832 | double diff1=sin_as_over_p-as_over_p*cos_as_over_p; |
| 833 | double diff2=one_minus_cos_as_over_p-as_over_p*sin_as_over_p; |
| 834 | double pdx_dp=pt_over_a*(cosphi*diff1-sinphi*diff2); |
| 835 | double dx_ds=cosl*(cosphi*cos_as_over_p-sinphi*sin_as_over_p); |
| 836 | double ds_dcosl=dz*cosl/(sinl*sinl2); |
| 837 | double dx_dcosl |
| 838 | =p_over_a*(cosphi*sin_as_over_p-sinphi*one_minus_cos_as_over_p) |
| 839 | +dx_ds*ds_dcosl; |
| 840 | double dx_dphi=-pt_over_a*(sinphi*sin_as_over_p+cosphi*one_minus_cos_as_over_p); |
| 841 | double var_z0=2.*tanl2*(var_tot)*double(2*N-1)/double(N*(N+1)); |
| 842 | |
| 843 | double var_x=var_x0+pdx_dp*pdx_dp*var_pt_over_pt_sq+var_pos_ms |
| 844 | +dx_dcosl*dx_dcosl*sinl2*var_lambda+dx_dphi*dx_dphi*var_phi |
| 845 | +dx_ds*dx_ds*var_z0/sinl2; |
| 846 | |
| 847 | double pdy_dp=pt_over_a*(sinphi*diff1+cosphi*diff2); |
| 848 | double dy_ds=cosl*(sinphi*cos_as_over_p+cosphi*sin_as_over_p); |
| 849 | double dy_dcosl |
| 850 | =p_over_a*(sinphi*sin_as_over_p+cosphi*one_minus_cos_as_over_p) |
| 851 | +dy_ds*ds_dcosl; |
| 852 | double dy_dphi=pt_over_a*(cosphi*sin_as_over_p-sinphi*one_minus_cos_as_over_p); |
| 853 | double var_y=var_y0+pdy_dp*pdy_dp*var_pt_over_pt_sq+var_pos_ms |
| 854 | +dy_dcosl*dy_dcosl*sinl2*var_lambda+dy_dphi*dy_dphi*var_phi |
| 855 | +dy_ds*dy_ds*var_z0/sinl2; |
| 856 | |
| 857 | // Rotate from global coordinate system into FDC local system |
| 858 | double cos2a=cosa*cosa; |
| 859 | double sin2a=sina*sina; |
| 860 | double var_d=(cos2a*var_x+sin2a*var_y)/(cosalpha*cosalpha); |
| 861 | double var_u=cos2a*var_y+sin2a*var_x; |
| 862 | |
| 863 | if (fit_type!=kHelical){ |
| 864 | // Factors take into account improved resolution after wire-based fit |
| 865 | var_d*=0.1; |
| 866 | var_u*=0.1; |
| 867 | } |
| 868 | |
| 869 | // Calculate chisq |
| 870 | chisq = resi*resi/(var_d+var_anode)+resic*resic/(var_u+var_cathode); |
| 871 | |
| 872 | // Probability of this hit being on the track |
| 873 | probability = TMath::Prob(chisq,2); |
| 874 | |
| 875 | if(probability>=MIN_HIT_PROB_FDC){ |
| 876 | pair<double,const DFDCPseudo*>myhit; |
| 877 | myhit.first=probability; |
| 878 | myhit.second=hit; |
| 879 | fdchits_tmp.push_back(myhit); |
| 880 | } |
| 881 | |
| 882 | // Optionally fill debug tree |
| 883 | if(fdchitsel){ |
| 884 | fdchitdbg.fit_type = fit_type; |
| 885 | fdchitdbg.p = p; |
| 886 | fdchitdbg.theta = rt->swim_steps[0].mom.Theta(); |
| 887 | fdchitdbg.mass = mass; |
| 888 | fdchitdbg.sigma_anode = sqrt(var_d+var_anode); |
| 889 | fdchitdbg.sigma_cathode = sqrt(var_u+var_cathode); |
| 890 | fdchitdbg.x = pos.X(); |
| 891 | fdchitdbg.y = pos.Y(); |
| 892 | fdchitdbg.z = pos.Z(); |
| 893 | fdchitdbg.s = s; |
| 894 | fdchitdbg.itheta02 = last_step->itheta02; |
| 895 | fdchitdbg.itheta02s = last_step->itheta02s; |
| 896 | fdchitdbg.itheta02s2 = last_step->itheta02s2; |
| 897 | fdchitdbg.dist = dist; |
| 898 | fdchitdbg.doca = doca; |
| 899 | fdchitdbg.resi = resi; |
| 900 | fdchitdbg.u = u; |
| 901 | fdchitdbg.u_cathodes = u_cathodes; |
| 902 | fdchitdbg.resic = resic; |
| 903 | fdchitdbg.chisq = chisq; |
| 904 | fdchitdbg.prob = probability; |
| 905 | fdchitdbg.sig_phi=sqrt(var_phi); |
| 906 | fdchitdbg.sig_lambda=sqrt(var_lambda); |
| 907 | fdchitdbg.sig_pt=sqrt(var_pt_over_pt_sq); |
| 908 | |
| 909 | fdchitsel->Fill(); |
| 910 | } |
| 911 | |
| 912 | if(HS_DEBUG_LEVEL>10){ |
| 913 | _DBG_std::cerr<<"libraries/TRACKING/DTrackHitSelectorALT2.cc" <<":"<<913<<" "; |
| 914 | if(probability>=MIN_HIT_PROB_FDC)jerr<<"----> "; |
| 915 | jerr<<"s="<<s<<" doca="<<doca<<" dist="<<dist<<" resi="<<resi<<" resic="<<resic<<" chisq="<<chisq<<" prob="<<probability<<endl; |
| 916 | } |
| 917 | } |
| 918 | // Order according to layer number and probability,then put the hits in the |
| 919 | // output list with the following algorithm: hits with the highest |
| 920 | // probability in a given layer are automatically put in the output list, |
| 921 | // but if there is more than one hit in a given layer, only those hits |
| 922 | // that are within +/-1 of the wire # of the most probable hit are added |
| 923 | // to the list. |
| 924 | sort(fdchits_tmp.begin(),fdchits_tmp.end(),DTrackHitSelector_fdchit_cmp); |
| 925 | int old_layer=1000,old_wire=1000; |
| 926 | for (unsigned int i=0;i<fdchits_tmp.size();i++){ |
| 927 | if (fdchits_tmp[i].second->wire->layer!=old_layer || |
| 928 | abs(fdchits_tmp[i].second->wire->wire-old_wire)==1){ |
| 929 | fdchits_out.push_back(fdchits_tmp[i].second); |
| 930 | } |
| 931 | old_wire=fdchits_tmp[i].second->wire->wire; |
| 932 | old_layer=fdchits_tmp[i].second->wire->layer; |
| 933 | } |
| 934 | } |
| 935 | void DTrackHitSelectorALT2::GetFDCHits(double Bz,double q, |
| 936 | const vector<DTrackFitter::Extrapolation_t> &extrapolations, const vector<const DFDCPseudo*> &fdchits_in, vector<const DFDCPseudo*> &fdchits_out,int N) const |
| 937 | { |
| 938 | // Vector of pairs storing the hit with the probability it is on the track |
| 939 | vector<pair<double,const DFDCPseudo*> >fdchits_tmp; |
| 940 | |
| 941 | /// Determine the probability that for each FDC hit that it came from the |
| 942 | /// track with the given trajectory. The probability is based on |
| 943 | /// the residual between the distance of closest approach |
| 944 | /// of the trajectory to the wire and the drift distance |
| 945 | /// and the distance along the wire. |
| 946 | |
| 947 | // Sort so innermost plane is first and outermost is last |
| 948 | vector<const DFDCPseudo*> fdchits_in_sorted = fdchits_in; |
| 949 | sort(fdchits_in_sorted.begin(),fdchits_in_sorted.end(),DTrackHitSelector_fdchit_in_cmp); |
| 950 | |
| 951 | // The variance on the residual due to measurement error. |
| 952 | double var_anode = 0.25*ONE_OVER_120.08333333333333; |
| 953 | const double VAR_CATHODE_STRIPS=0.000225; |
| 954 | double var_cathode = VAR_CATHODE_STRIPS+0.15*0.15*ONE_OVER_120.08333333333333; |
| 955 | double var_tot=var_anode+var_cathode; |
| 956 | |
| 957 | // To estimate the impact of errors in the track momentum on the variance of the residual, |
| 958 | // use a helical approximation. |
| 959 | double a=-0.003*Bz*q; |
| 960 | DVector3 mom=extrapolations[extrapolations.size()/2].momentum; |
| 961 | double p=mom.Mag(); |
| 962 | double p_over_a=p/a; |
| 963 | double a_over_p=1./p_over_a; |
| 964 | double lambda=M_PI_21.57079632679489661923-mom.Theta(); |
| 965 | double tanl=tan(lambda); |
| 966 | double tanl2=tanl*tanl; |
| 967 | double sinl=sin(lambda); |
| 968 | double cosl=cos(lambda); |
| 969 | double cosl2=cosl*cosl; |
| 970 | double sinl2=sinl*sinl; |
| 971 | double pt_over_a=cosl*p_over_a; |
| 972 | double z0=extrapolations[0].position.z(); |
| 973 | // variances |
| 974 | double var_lambda=0.,var_phi=0.,var_lambda_res=0.,var_k=0.; |
| 975 | double var_phi_radial=0.; |
| 976 | double var_x0=0.0,var_y0=0.0; |
| 977 | double var_z0=2.*tanl2*(var_tot)*double(2*N-1)/double(N*(N+1)); |
| 978 | |
| 979 | // Loop over hits |
| 980 | bool most_downstream_hit=true; |
| 981 | int last_extrapolation_index=extrapolations.size()-1; |
| 982 | vector<const DFDCPseudo*>::const_reverse_iterator iter; |
| 983 | for(iter=fdchits_in_sorted.rbegin(); iter!=fdchits_in_sorted.rend(); iter++){ |
| 984 | const DFDCPseudo *hit = *iter; |
| 985 | for (int k=last_extrapolation_index;k>=0;k--){ |
| 986 | // Position along trajectory |
| 987 | DVector3 pos=extrapolations[k].position; |
| 988 | double dz=pos.z()-z0; |
| 989 | double s=extrapolations[k].s; |
| 990 | if (fabs(pos.z()-hit->wire->origin.z())<0.5){ |
| 991 | // Variance in dip angle due to multiple scattering |
| 992 | var_lambda = extrapolations[k].theta2ms_sum/3.; |
| 993 | // the above expression seems to lead to overestimation of the uncertainty in the dip angle after the wire-based pass.. |
| 994 | var_lambda*=0.01; |
| 995 | // Variance in phi due to multiple scattering |
| 996 | var_phi=var_lambda*(1.+tanl2); |
| 997 | |
| 998 | // azimuthal angle |
| 999 | double phi=extrapolations[k].momentum.Phi(); |
| 1000 | double sinphi=sin(phi); |
| 1001 | double cosphi=cos(phi); |
| 1002 | |
| 1003 | if (most_downstream_hit){ |
| 1004 | // variance in the curvature k due to resolution and multiple scattering |
| 1005 | double s_sq=s*s; |
| 1006 | double s_4=s_sq*s_sq; |
| 1007 | double sum_s_theta_ms=extrapolations[k].s_theta_ms_sum; |
| 1008 | double var_k_ms=(4./3.)*sum_s_theta_ms*sum_s_theta_ms/s_4; |
| 1009 | double var_k_res=var_tot*0.0720/double(N+4)/s_4/(cosl2*cosl2); |
| 1010 | var_k=var_k_ms+var_k_res; |
| 1011 | |
| 1012 | // Variance in dip angle due to measurement error |
| 1013 | var_lambda_res=12.0*var_tot*double(N-1)/double(N*(N+1)) |
| 1014 | *sinl2*sinl2/s_sq; |
| 1015 | |
| 1016 | // Variance in phi, crude approximation |
| 1017 | double tan_s_2R=tan(s*cosl/(2.*pt_over_a)); |
| 1018 | var_phi_radial=tan_s_2R*tan_s_2R*pt_over_a*pt_over_a*var_k; |
| 1019 | |
| 1020 | most_downstream_hit=false; |
| 1021 | } |
| 1022 | // Include error in lambda due to measurements |
| 1023 | var_lambda+=var_lambda_res; |
| 1024 | |
| 1025 | // Include uncertainty in phi due to uncertainty in the center of |
| 1026 | // the circle |
| 1027 | var_phi+=var_phi_radial; |
| 1028 | |
| 1029 | // Variance in position due to multiple scattering |
| 1030 | double var_pos_ms=extrapolations[k].s_theta_ms_sum/3.; |
| 1031 | |
| 1032 | // doca to wire |
| 1033 | double x=pos.x(); |
| 1034 | double y=pos.y(); |
| 1035 | double dx=hit->xy.X()-x; |
| 1036 | double dy=hit->xy.Y()-y; |
| 1037 | double doca=sqrt(dx*dx+dy*dy); |
| 1038 | |
| 1039 | // Direction variables for wire |
| 1040 | double cosa=hit->wire->udir.y(); |
| 1041 | double sina=hit->wire->udir.x(); |
| 1042 | |
| 1043 | // Cathode Residual |
| 1044 | double u=x*sina+y*cosa; |
| 1045 | double u_cathodes = hit->s; |
| 1046 | double resic = u - u_cathodes; |
| 1047 | |
| 1048 | // Get "measured" distance to wire. |
| 1049 | double dist=0.25; |
| 1050 | |
| 1051 | // Take into account non-normal incidence to FDC plane |
| 1052 | double pz=extrapolations[k].momentum.z(); |
| 1053 | double tx=extrapolations[k].momentum.x()/pz; |
| 1054 | double ty=extrapolations[k].momentum.y()/pz; |
| 1055 | double tu=tx*cosa-ty*sina; |
| 1056 | double alpha=atan(tu); |
| 1057 | double cosalpha=cos(alpha); |
| 1058 | |
| 1059 | // Anode Residual |
| 1060 | double resi = dist - doca/cosalpha; |
| 1061 | |
| 1062 | // Initialize some probability-related variables needed later |
| 1063 | double probability=0.,chisq=0.; |
| 1064 | |
| 1065 | // Variances in x and y due to uncertainty in track parameters |
| 1066 | double as_over_p=s*a_over_p; |
| 1067 | double sin_as_over_p=sin(as_over_p); |
| 1068 | double cos_as_over_p=cos(as_over_p); |
| 1069 | double one_minus_cos_as_over_p=1-cos_as_over_p; |
| 1070 | double diff1=sin_as_over_p-as_over_p*cos_as_over_p; |
| 1071 | double diff2=one_minus_cos_as_over_p-as_over_p*sin_as_over_p; |
| 1072 | double dx_ds=cosl*(cosphi*cos_as_over_p-sinphi*sin_as_over_p); |
| 1073 | double ds_dcosl=dz*cosl/(sinl*sinl2); |
| 1074 | double dx_dcosl |
| 1075 | =p_over_a*(cosphi*sin_as_over_p-sinphi*one_minus_cos_as_over_p) |
| 1076 | +dx_ds*ds_dcosl; |
| 1077 | double dx_dphi=-pt_over_a*(sinphi*sin_as_over_p+cosphi*one_minus_cos_as_over_p); |
| 1078 | double dx_dk=2.*pt_over_a*pt_over_a*(sinphi*diff2-cosphi*diff1); |
| 1079 | |
| 1080 | double var_x=var_x0+dx_dk*dx_dk*var_k+var_pos_ms |
| 1081 | +dx_dcosl*dx_dcosl*sinl2*var_lambda+dx_dphi*dx_dphi*var_phi |
| 1082 | +dx_ds*dx_ds*var_z0/sinl2; |
| 1083 | |
| 1084 | double dy_dk=-2.*pt_over_a*pt_over_a*(sinphi*diff1+cosphi*diff2); |
| 1085 | double dy_ds=cosl*(sinphi*cos_as_over_p+cosphi*sin_as_over_p); |
| 1086 | double dy_dcosl |
| 1087 | =p_over_a*(sinphi*sin_as_over_p+cosphi*one_minus_cos_as_over_p) |
| 1088 | +dy_ds*ds_dcosl; |
| 1089 | double dy_dphi=pt_over_a*(cosphi*sin_as_over_p-sinphi*one_minus_cos_as_over_p); |
| 1090 | double var_y=var_y0+dy_dk*dy_dk*var_k+var_pos_ms |
| 1091 | +dy_dcosl*dy_dcosl*sinl2*var_lambda+dy_dphi*dy_dphi*var_phi |
| 1092 | +dy_ds*dy_ds*var_z0/sinl2; |
| 1093 | |
| 1094 | // Rotate from global coordinate system into FDC local system |
| 1095 | double cos2a=cosa*cosa; |
| 1096 | double sin2a=sina*sina; |
| 1097 | double var_d=(cos2a*var_x+sin2a*var_y)/(cosalpha*cosalpha); |
| 1098 | double var_u=cos2a*var_y+sin2a*var_x; |
| 1099 | |
| 1100 | // Factors take into account improved resolution after wire-based fit |
| 1101 | //var_d*=0.1; |
| 1102 | //var_u*=0.1; |
| 1103 | |
| 1104 | // Calculate chisq |
| 1105 | chisq = resi*resi/(var_d+var_anode)+resic*resic/(var_u+var_cathode); |
| 1106 | |
| 1107 | // Probability of this hit being on the track |
| 1108 | probability = TMath::Prob(chisq,2); |
| 1109 | |
| 1110 | if(probability>=MIN_HIT_PROB_FDC){ |
| 1111 | pair<double,const DFDCPseudo*>myhit; |
| 1112 | myhit.first=probability; |
| 1113 | myhit.second=hit; |
| 1114 | fdchits_tmp.push_back(myhit); |
| 1115 | } |
| 1116 | // Optionally fill debug tree |
| 1117 | if(fdchitsel){ |
| 1118 | fdchitdbg.fit_type = kWireBased; |
| 1119 | fdchitdbg.p = p; |
| 1120 | fdchitdbg.theta = extrapolations[k].momentum.Theta(); |
| 1121 | //fdchitdbg.mass = mass; |
| 1122 | fdchitdbg.sigma_anode = sqrt(var_d+var_anode); |
| 1123 | fdchitdbg.sigma_cathode = sqrt(var_u+var_cathode); |
| 1124 | fdchitdbg.x = pos.X(); |
| 1125 | fdchitdbg.y = pos.Y(); |
| 1126 | fdchitdbg.z = pos.Z(); |
| 1127 | fdchitdbg.s = s; |
| 1128 | fdchitdbg.itheta02 = extrapolations[k].theta2ms_sum; |
| 1129 | //fdchitdbg.itheta02s = last_step->itheta02s; |
| 1130 | //fdchitdbg.itheta02s2 = last_step->itheta02s2; |
| 1131 | fdchitdbg.dist = dist; |
| 1132 | fdchitdbg.doca = doca; |
| 1133 | fdchitdbg.resi = resi; |
| 1134 | fdchitdbg.u = u; |
| 1135 | fdchitdbg.u_cathodes = u_cathodes; |
| 1136 | fdchitdbg.resic = resic; |
| 1137 | fdchitdbg.chisq = chisq; |
| 1138 | fdchitdbg.prob = probability; |
| 1139 | fdchitdbg.sig_phi=sqrt(var_phi); |
| 1140 | fdchitdbg.sig_lambda=sqrt(var_lambda); |
| 1141 | fdchitdbg.sig_pt=fabs(2.*pt_over_a)*sqrt(var_k); |
| 1142 | |
| 1143 | fdchitsel->Fill(); |
| 1144 | } |
| 1145 | |
| 1146 | if(HS_DEBUG_LEVEL>10){ |
| 1147 | _DBG_std::cerr<<"libraries/TRACKING/DTrackHitSelectorALT2.cc" <<":"<<1147<<" "; |
| 1148 | if(probability>=MIN_HIT_PROB_FDC)jerr<<"----> "; |
| 1149 | jerr<<"s="<<s<<" doca="<<doca<<" dist="<<dist<<" resi="<<resi<<" resic="<<resic<<" chisq="<<chisq<<" prob="<<probability<<endl; |
| 1150 | } |
| 1151 | |
| 1152 | } |
| 1153 | } |
| 1154 | |
| 1155 | } |
| 1156 | // Order according to layer number and probability,then put the hits in the |
| 1157 | // output list with the following algorithm: hits with the highest |
| 1158 | // probability in a given layer are automatically put in the output list, |
| 1159 | // but if there is more than one hit in a given layer, only those hits |
| 1160 | // that are within +/-1 of the wire # of the most probable hit are added |
| 1161 | // to the list. |
| 1162 | sort(fdchits_tmp.begin(),fdchits_tmp.end(),DTrackHitSelector_fdchit_cmp); |
| 1163 | int old_layer=1000,old_wire=1000; |
| 1164 | for (unsigned int i=0;i<fdchits_tmp.size();i++){ |
| 1165 | if (fdchits_tmp[i].second->wire->layer!=old_layer || |
| 1166 | abs(fdchits_tmp[i].second->wire->wire-old_wire)==1){ |
| 1167 | fdchits_out.push_back(fdchits_tmp[i].second); |
| 1168 | } |
| 1169 | old_wire=fdchits_tmp[i].second->wire->wire; |
| 1170 | old_layer=fdchits_tmp[i].second->wire->layer; |
| 1171 | } |
| 1172 | } |