// $Id$ // // File: DTrackHitSelectorALT1.cc // Created: Fri Feb 6 08:22:58 EST 2009 // Creator: davidl (on Darwin harriet.jlab.org 9.6.0 i386) // #include #include #include "DTrackHitSelectorALT1.h" #ifndef ansi_escape #define ansi_escape ((char)0x1b) #define ansi_bold ansi_escape<<"[1m" #define ansi_normal ansi_escape<<"[0m" #define ansi_red ansi_escape<<"[31m" #define ansi_green ansi_escape<<"[32m" #define ansi_blue ansi_escape<<"[34m" #endif // ansi_escape #define ONE_OVER_SQRT12 0.288675 bool static DTrackHitSelector_cdchit_cmp(paira, pairb){ if (a.second->wire->ring!=b.second->wire->ring) return (a.second->wire->ring>b.second->wire->ring); return (a.first>b.first); } bool static DTrackHitSelector_fdchit_cmp(paira, pairb){ if (a.second->wire->layer!=b.second->wire->layer) return (a.second->wire->layer>b.second->wire->layer); return (a.first>b.first); } //--------------------------------- // DTrackHitSelectorALT1 (Constructor) //--------------------------------- DTrackHitSelectorALT1::DTrackHitSelectorALT1(jana::JEventLoop *loop):DTrackHitSelector(loop) { HS_DEBUG_LEVEL = 0; MAKE_DEBUG_TREES = false; gPARMS->SetDefaultParameter("TRKFIT:HS_DEBUG_LEVEL", HS_DEBUG_LEVEL, "Debug verbosity level for hit selector used in track fitting (0=no debug messages)"); gPARMS->SetDefaultParameter("TRKFIT:MAKE_DEBUG_TREES", MAKE_DEBUG_TREES, "Create a TTree with debugging info on hit selection for the FDC and CDC"); cdchitsel = NULL; fdchitsel = NULL; if(MAKE_DEBUG_TREES){ loop->GetJApplication()->Lock(); cdchitsel= (TTree*)gROOT->FindObject("cdchitsel"); if(!cdchitsel){ cdchitsel = new TTree("cdchitsel", "CDC Hit Selector"); cdchitsel->Branch("H", &cdchitdbg, "fit_type/I:p/F:theta:mass:sigma:mom_factor:x:y:z:s:s_factor:itheta02:itheta02s:itheta02s2:dist:doca:resi:sigma_total:chisq:prob"); }else{ _DBG__; jerr<<" !!! WARNING !!!"<FindObject("fdchitsel"); if(!fdchitsel){ fdchitsel = new TTree("fdchitsel", "FDC Hit Selector"); fdchitsel->Branch("H", &fdchitdbg, "fit_type/I:p/F:theta:mass:sigma_anode:sigma_cathode:mom_factor_anode:mom_factor_cathode:x:y:z:s:s_factor_anode:s_factor_cathode:itheta02:itheta02s:itheta02s2:dist:doca:resi:u:u_cathodes:resic:sigma_anode_total:sigma_cathode_total:chisq:prob:prob_anode:prob_cathode:pull_anode:pull_cathode"); }else{ _DBG__; jerr<<" !!! WARNING !!!"<GetJApplication()->Unlock(); } } //--------------------------------- // ~DTrackHitSelectorALT1 (Destructor) //--------------------------------- DTrackHitSelectorALT1::~DTrackHitSelectorALT1() { } //--------------------------------- // GetCDCHits //--------------------------------- void DTrackHitSelectorALT1::GetCDCHits(fit_type_t fit_type, DReferenceTrajectory *rt, const vector &cdchits_in, vector &cdchits_out) const { // Vector of pairs storing the hit with the probability it is on the track vector >cdchits_tmp; /// Determine the probability that for each CDC hit that it came from the /// track with the given trajectory. /// /// This will calculate a probability for each CDC hit that /// it came from the track represented by the given /// DReference trajectory. The probability is based on /// the residual between the distance of closest approach /// of the trajectory to the wire and the drift time for /// time-based tracks and the distance to the wire for /// wire-based tracks. // Calculate beta of particle. double p = rt->swim_steps[0].mom.Mag(); double my_mass=rt->GetMass(); double one_over_beta =sqrt(1.0+my_mass*my_mass/(p*p)); // The error on the residual. This will be different based on the // quality of the track and whether MULS is on or not etc. // In principle, this could also depend on the momentum parameters // of the track. double sigma; switch(fit_type){ case kTimeBased: sigma = 0.8*ONE_OVER_SQRT12; break; case kWireBased: sigma = 1.6*ONE_OVER_SQRT12; break; case kHelical: default: sigma = 8.0*ONE_OVER_SQRT12; } // Low-momentum tracks are more poorly defined than high-momentum tracks. // We account for that here by increasing the error as a function of momentum double g = 0.350/sqrt(log(2.0)); // total guess double p_over_g=p/g; //sigma *= 1.0 + exp(-pow(rt->swim_steps[0].mom.Mag()/g,2.0)); double mom_factor = 1.0+exp(-p_over_g*p_over_g); sigma *= mom_factor; int old_straw=1000,old_ring=1000; // Minimum probability of hit belonging to wire and still be accepted double MIN_HIT_PROB = 0.05; vector::const_iterator iter; for(iter=cdchits_in.begin(); iter!=cdchits_in.end(); iter++){ const DCDCTrackHit *hit = *iter; // Skip hit if it is on the same wire as the previous hit if (fit_type!=kTimeBased){ if (hit->wire->ring == old_ring && hit->wire->straw==old_straw) continue; old_ring=hit->wire->ring; old_straw=hit->wire->straw; } // Find the DOCA to this wire double s; double doca = rt->DistToRT(hit->wire, &s); // Get "measured" distance to wire. For time-based tracks // this is calculated from the drift time. For all other // tracks, this is assumed to be half a cell size double dist; if(fit_type == kTimeBased){ // Distance using drift time // NOTE: Right now we assume pions for the TOF // and a constant drift velocity of 55um/ns double tof = s*one_over_beta/29.98; dist = (hit->tdrift - tof)*55E-4; }else{ dist = 0.4; // =0.8/2.0; half cell-size } // For time-based and wire-based tracks, the fit was // weighted for multiple scattering by material times // angle giving preference to the begining of the // track. Take this into account here by enhancing the // error for hits further from the vertex double sigma_total = sigma; double s_factor = 1.0; if(fit_type == kTimeBased || fit_type == kWireBased){ s_factor = 1.0 + s/50.0; // double error at 50cm out (guess for now) sigma_total *= s_factor; } // Residual double resi = dist - doca; //double chisq = pow(resi/sigma_total, 2.0); double chisq=resi*resi/(sigma_total*sigma_total); // Use chi-sq probability function with Ndof=1 to calculate probability double probability = TMath::Prob(chisq, 1); if(probability>=MIN_HIT_PROB){ pairmyhit; myhit.first=probability; myhit.second=hit; cdchits_tmp.push_back(myhit); } // Optionally fill debug tree if(cdchitsel){ DVector3 pos = rt->GetLastDOCAPoint(); const DReferenceTrajectory::swim_step_t *last_step = rt->GetLastSwimStep(); cdchitdbg.fit_type = fit_type; cdchitdbg.p = p; cdchitdbg.theta = rt->swim_steps[0].mom.Theta(); cdchitdbg.mass = my_mass; cdchitdbg.sigma = sigma; cdchitdbg.mom_factor = mom_factor; cdchitdbg.x = pos.X(); cdchitdbg.y = pos.Y(); cdchitdbg.z = pos.Z(); cdchitdbg.s = s; cdchitdbg.s_factor = s_factor; cdchitdbg.itheta02 = last_step->itheta02; cdchitdbg.itheta02s = last_step->itheta02s; cdchitdbg.itheta02s2 = last_step->itheta02s2; cdchitdbg.dist = dist; cdchitdbg.doca = doca; cdchitdbg.resi = resi; cdchitdbg.sigma_total = sigma_total; cdchitdbg.chisq = chisq; cdchitdbg.prob = probability; cdchitsel->Fill(); static bool printed_first = false; if(!printed_first){ _DBG_<<"=== Printing first entry for CDC hit selector debug tree ==="<time - tof)*55E-4; }else{ dist = 0.25; //= 0.5/2.0; half cell-size } // Anode Residual double resi = dist - doca; // Cathode Residual double u=rt->GetLastDistAlongWire(); double u_cathodes = hit->s; double resic = u - u_cathodes; // Probability of this hit being on the track //double chisq = pow(resi/sigma_anode, 2.0) + pow(resic/sigma_cathode, 2.0); double sigma_anode_total = sigma_anode*mom_factor*mass_factor; double sigma_cathode_total = sigma_cathode*mom_factor*mass_factor; double pull_anode = resi/sigma_anode_total; double pull_cathode = resic/sigma_cathode_total; double chisq = pull_anode*pull_anode + pull_cathode*pull_cathode; double probability = TMath::Prob(chisq, 2); if(probability>=MIN_HIT_PROB){ pairmyhit; myhit.first=probability; myhit.second=hit; fdchits_tmp.push_back(myhit); } // Optionally fill debug tree if(fdchitsel){ DVector3 pos = rt->GetLastDOCAPoint(); fdchitdbg.fit_type = fit_type; fdchitdbg.p = p; fdchitdbg.theta = rt->swim_steps[0].mom.Theta(); fdchitdbg.mass = my_mass; fdchitdbg.sigma_anode = sigma_anode; fdchitdbg.sigma_cathode = sigma_cathode; fdchitdbg.mom_factor_anode = mom_factor; fdchitdbg.mom_factor_cathode = mom_factor; fdchitdbg.x = pos.X(); fdchitdbg.y = pos.Y(); fdchitdbg.z = pos.Z(); fdchitdbg.s = s; fdchitdbg.s_factor_anode = 1.0; fdchitdbg.s_factor_cathode = 1.0; fdchitdbg.itheta02 = last_step->itheta02; fdchitdbg.itheta02s = last_step->itheta02s; fdchitdbg.itheta02s2 = last_step->itheta02s2; fdchitdbg.dist = dist; fdchitdbg.doca = doca; fdchitdbg.resi = resi; fdchitdbg.u = u; fdchitdbg.u_cathodes = u_cathodes; fdchitdbg.resic = resic; fdchitdbg.sigma_anode_total = sigma_anode_total; fdchitdbg.sigma_cathode_total = sigma_cathode_total; fdchitdbg.chisq = chisq; fdchitdbg.prob = probability; fdchitdbg.prob_anode = TMath::Prob(pull_anode*pull_anode, 1); fdchitdbg.prob_cathode = TMath::Prob(pull_cathode*pull_cathode, 1); fdchitdbg.pull_anode = pull_anode; fdchitdbg.pull_cathode = pull_cathode; //_DBG_<<"fdchitdbg.pull_anode="<