// $Id$ // // File: JEventProcessor_trkassassin.cc // Created: Fri May 5 06:46:06 EDT 2017 // Creator: davidl (on Linux gluon119.jlab.org 2.6.32-642.3.1.el6.x86_64 x86_64) // #include using namespace std; #include "JEventProcessor_trkassassin.h" using namespace jana; #include #include #include #include #include #include #include // Routine used to create our JEventProcessor #include #include extern "C"{ void InitPlugin(JApplication *app){ InitJANAPlugin(app); app->AddProcessor(new JEventProcessor_trkassassin()); } } // "C" //------------------ // JEventProcessor_trkassassin (Constructor) //------------------ JEventProcessor_trkassassin::JEventProcessor_trkassassin() { } //------------------ // ~JEventProcessor_trkassassin (Destructor) //------------------ JEventProcessor_trkassassin::~JEventProcessor_trkassassin() { } //------------------ // init //------------------ jerror_t JEventProcessor_trkassassin::init(void) { tTrkParms = new TTree("trkparms", "Tracking Parameters"); tTrkParms->Branch("trkparms", "trkparms", &tparms); tTrkParms->Branch("Ncdcpulls", &Ncdcpulls, "Ncdcpulls/I"); tTrkParms->Branch("Nfdcpulls", &Nfdcpulls, "Nfdcpulls/I"); tTrkParms->Branch("cdcpulls", cdcpulls, "cdcpulls[Ncdcpulls]/F"); tTrkParms->Branch("fdcpulls", fdcpulls, "fdcpulls[Nfdcpulls]/F"); tHitParms = new TTree("hitparms", "FDC Hit Parameters"); tHitParms->Branch("hitparms", "hitparms", &hparms); tTrkHypoth = new TTree("trkhypoth", "Track Hypotheses"); tTrkHypoth->Branch("trkhypoth", "trkhy", &trkhy); Nwirebased_per_candidate = new TH1D("Nwirebased_per_candidate", "Num. Wire-based Tracks per Candidate", 20, 0.0, 20.0); Ntimebased_per_candidate = new TH1D("Ntimebased_per_candidate", "Num. Time-based Tracks per Candidate", 20, 0.0, 20.0); Ntimebased_per_wirebased = new TH1D("Ntimebased_per_wirebased", "Num. Time-based Tracks per Wire-based track", 20, 0.0, 20.0); return NOERROR; } //------------------ // brun //------------------ jerror_t JEventProcessor_trkassassin::brun(JEventLoop *eventLoop, int32_t runnumber) { // This is called whenever the run number changes return NOERROR; } //------------------ // evnt //------------------ jerror_t JEventProcessor_trkassassin::evnt(JEventLoop *loop, uint64_t eventnumber) { vector tbts; vector wbts; vector candidates; vector cdctrackhits; vector fdcpseudos; loop->Get(tbts); loop->Get(wbts); loop->Get(candidates); loop->Get(cdctrackhits); loop->Get(fdcpseudos); const DTrackHitSelector *hitselector = NULL; loop->GetSingle(hitselector); map > can_map; // value is pid for successfully fit tbt // Count successfully fit tracks for each candidate map wbt_count; map tbt_count; map tbt_wbt_count; for(auto wbt : wbts){ const DTrackCandidate *can = NULL; wbt->GetSingle(can); wbt_count[can]++; } vector all_hitparms; for(auto tbt : tbts){ // Get pointers to candidate and wire-based // track from which this came vector cans; tbt->Get(cans); const DTrackCandidate *can = cans.empty() ? NULL:cans[0]; tbt_count[can]++; const DTrackWireBased *wbt = NULL; tbt->GetSingle(wbt); tbt_wbt_count[wbt]++; can_map[can].insert(tbt->PID()); // Get list of DFDCPseudo hits from each stage set can_hits; set wbt_hits; set tbt_hits; // Candidates and Wire-based tracks do not fill // the pulls vector (this is only done when the // Kalman smoother is called). Use associated // objects to get which hits were used there. vector vcan_hits; vector vwbt_hits; vector vtbt_hits; can->Get(vcan_hits); wbt->Get(vwbt_hits); tbt->Get(vtbt_hits); for(auto &h : vcan_hits) can_hits.insert(h); for(auto &h : vwbt_hits) wbt_hits.insert(h); for(auto &h : vtbt_hits) tbt_hits.insert(h); // Use entries in pulls vector for Time-based tracks // This seems to be a subset of the hits in the // asoociated objects and Simon isn't sure why. // (They should be the same). The assumption here // is this more restrictive set represents the // purest track hits. // for(auto &p : tbt->pulls){ // if(p.fdc_hit!=NULL) tbt_hits.insert(p.fdc_hit); // } // Swim reference trajctory for candidate parameters // using Time-based (and presumably wire-based) mass. DReferenceTrajectory rt(*(tbt->rt)); rt.SetMass(tbt->mass()); rt.FastSwim(can->position(),can->momentum(),can->charge(),2000.0,0.,370.); // vector vwbt_input_hits; // hitselector->GetFDCHits(DTrackHitSelector::kWireBased, &rt, fdcpseudos, vwbt_input_hits, can->Ndof+3); // // // Get list of output hits from timebased tracking // vector vtbt_hits; // tbt->Get(vtbt_hits); // // // Convert to STL sets so we can easily check if hit appears in both // set wbt_input_hits; // set tbt_hits; // for(auto h : vwbt_input_hits) wbt_input_hits.insert(h); // for(auto h : vtbt_hits ) tbt_hits.insert(h); // // for(auto &p : can->pulls){ // if(p.fdc_hit!=NULL) can_hits.insert(p.fdc_hit); // } // for(auto &p : wbt->pulls){ // if(p.fdc_hit!=NULL) wbt_hits.insert(p.fdc_hit); // } // for(auto &p : tbt->pulls){ // if(p.fdc_hit!=NULL) tbt_hits.insert(p.fdc_hit); // } // _DBG_<<"can_hits.size()=" << can_hits.size() << endl; // _DBG_<<"wbt_hits.size()=" << wbt_hits.size() << endl; // _DBG_<<"tbt_hits.size()=" << tbt_hits.size() << endl; // Calculate average values from candidate hits Float_t x_avg = 0.0; Float_t y_avg = 0.0; Float_t u_avg = 0.0; Float_t v_avg = 0.0; Float_t t_u_avg = 0.0; Float_t t_v_avg = 0.0; for(auto h : can_hits){ x_avg += h->xy.X(); y_avg += h->xy.Y(); u_avg += h->u; v_avg += h->v; t_u_avg += h->t_u; t_v_avg += h->t_v; } if(!can_hits.empty()){ Float_t N = (Float_t)can_hits.size(); x_avg /= N; y_avg /= N; u_avg /= N; v_avg /= N; t_u_avg /= N; t_v_avg /= N; } // Loop over all FDC hits and create hitparms for each hitparms hp; hp.candidate_p = can->momentum().Mag(); hp.candidate_theta = can->momentum().Theta(); hp.candidate_phi = can->momentum().Phi(); for(auto h : fdcpseudos){ hp.on_candidate = can_hits.count(h); hp.on_wirebased = wbt_hits.count(h); hp.on_timebased = tbt_hits.count(h); Float_t xrel = h->xy.X() - x_avg; Float_t yrel = h->xy.Y() - y_avg; Float_t urel = h->u - u_avg; Float_t vrel = h->v - v_avg; Float_t turel = h->t_u - h->time; Float_t tvrel = h->t_v - h->time; hp.xy_dist_from_avg = sqrt( xrel*xrel + yrel*yrel ); hp.uv_dist_from_avg = sqrt( urel*urel + vrel*vrel ); hp.u_tdiff_from_avg = h->t_u - t_u_avg; hp.v_tdiff_from_avg = h->t_v - t_v_avg; hp.tu_minus_tv = h->t_u - h->t_v; hp.phiu_minus_phiv = h->phi_u - h->phi_v; hp.cos_phiu_minus_phiv = cos(hp.phiu_minus_phiv); hp.w = h->w; hp.dw = h->dw; hp.w_c = h->w_c; hp.s = h->s; hp.ds = h->ds; hp.centroid_u_pos = h->cluster_u.pos; hp.centroid_u_q = h->cluster_u.q; hp.centroid_u_q_from_pulse_height = h->cluster_u.q_from_pulse_height; hp.centroid_u_numstrips = h->cluster_u.numstrips; hp.centroid_v_pos = h->cluster_v.pos; hp.centroid_v_q = h->cluster_v.q; hp.centroid_v_q_from_pulse_height = h->cluster_v.q_from_pulse_height; hp.centroid_v_numstrips = h->cluster_v.numstrips; hp.t = h->time; hp.t_dist = sqrt( turel*turel + tvrel*tvrel ); hp.dE = h->dE; hp.q = h->q; hp.covxx = h->covxx; hp.covxy = h->covxy; hp.covyy = h->covyy; hp.status = h->status; hp.layer = h->wire->layer; hp.wire = h->wire->wire; // Additional parameters based on DTrackHitSelectorALT2 are filled // by other routine below. bool ok = FillAdditionalSelectorParms(h, hp, &rt, (double)(can->Ndof+3)); if(ok) all_hitparms.push_back(hp); } } // Mass hypotheses to check for vector hypotheses; hypotheses.push_back(Positron); hypotheses.push_back(PiPlus); hypotheses.push_back(KPlus); hypotheses.push_back(Proton); hypotheses.push_back(Electron); hypotheses.push_back(PiMinus); hypotheses.push_back(KMinus); hypotheses.push_back(AntiProton); // Lock ROOT mutex and fill ROOT objects japp->RootFillLock(this); for(auto c : candidates){ Nwirebased_per_candidate->Fill(wbt_count[c]); Ntimebased_per_candidate->Fill(tbt_count[c]); } for(auto w : wbts){ Ntimebased_per_wirebased->Fill(tbt_wbt_count[w]); } for(auto tbt : tbts){ tparms.chisq = tbt->chisq; tparms.Ndof = tbt->Ndof; // pulls Ncdcpulls = 0; Nfdcpulls = 0; for(auto &p : tbt->pulls){ if(p.err == 0.0) continue; if(p.cdc_hit){ if(Ncdcpulls < kMaxPulls) cdcpulls[Ncdcpulls++] = p.resi/p.err; } if(p.fdc_hit){ if(Nfdcpulls < kMaxPulls) fdcpulls[Nfdcpulls++] = p.resi/p.err; } } // Number of hits at wire-based and candidate stages const DTrackWireBased *wbt = NULL; const DTrackCandidate *can = NULL; tbt->GetSingle(wbt); wbt->GetSingle(can); vector cdchits; vector fdchits; tbt->Get(cdchits); tbt->Get(fdchits); tparms.Ncdchits_timebased = cdchits.size(); tparms.Nfdchits_timebased = fdchits.size(); wbt->Get(cdchits); wbt->Get(fdchits); tparms.Ncdchits_wirebased = cdchits.size(); tparms.Nfdchits_wirebased = fdchits.size(); can->Get(cdchits); can->Get(fdchits); tparms.Ncdchits_candidate = cdchits.size(); tparms.Nfdchits_candidate = fdchits.size(); tTrkParms->Fill(); } for(auto &hp : all_hitparms) { hparms = hp; tHitParms->Fill(); } for(auto p : can_map){ auto c = p.first; auto &s = p.second; vector cdchits; vector fdchits; c->Get(cdchits); c->Get(fdchits); trkhy.Ncdchits_candidate = cdchits.size(); trkhy.Nfdchits_candidate = fdchits.size(); trkhy.chisq = c->chisq; trkhy.Ndof = c->Ndof; trkhy.p = c->momentum().Mag(); trkhy.theta = c->momentum().Theta(); for(auto h : hypotheses){ trkhy.has_tbt = s.count(h); trkhy.pid = h; tTrkHypoth->Fill(); } } japp->RootFillUnLock(this); return NOERROR; } //------------------ // FillAdditionalSelectorParms //------------------ bool JEventProcessor_trkassassin::FillAdditionalSelectorParms(const DFDCPseudo *hit, hitparms &hp, DReferenceTrajectory *rt, double N) { // The variance on the residual due to measurement error. double var_anode = 0.25/12.0; // scale factor reflects field-sense wire separation const double VAR_CATHODE_STRIPS=0.000225; double var_cathode = VAR_CATHODE_STRIPS; // To estimate the impact of errors in the track momentum on the variance of the residual, // use a helical approximation. const DReferenceTrajectory::swim_step_t *my_step=&rt->swim_steps[0]; double Bz0=my_step->B.z(); double z0=my_step->origin.z(); double a=-0.003*Bz0*rt->q; double p=my_step->mom.Mag(); //double p_sq=p*p; double p_over_a=p/a; double a_over_p=1./p_over_a; double lambda=M_PI_2-my_step->mom.Theta(); double cosl=cos(lambda); //double cosl2=cosl*cosl; double sinl=sin(lambda); double sinl2=sinl*sinl; double tanl=tan(lambda); double tanl2=tanl*tanl; double pt_over_a=cosl*p_over_a; double phi=my_step->mom.Phi(); double cosphi=cos(phi); double sinphi=sin(phi); double var_lambda=0.,var_phi=0.,var_lambda_res=0.; //double mass=rt->GetMass(); double var_x0=0.01,var_y0=0.01; double var_pt_over_pt_sq=0.; // // Loop over hits // bool most_downstream_hit=true; // Find the DOCA to this wire double s; double doca = rt->DistToRT(hit->wire, &s); double MAX_DOCA=2.5; if(!isfinite(doca)) return false; if(!isfinite(s))return false; if (s<=0.) return false; if (doca>MAX_DOCA)return false; const DReferenceTrajectory::swim_step_t *last_step = rt->GetLastSwimStep(); // Position along trajectory DVector3 pos=rt->GetLastDOCAPoint(); double dz=pos.z()-z0; // Direction variables for wire double cosa=hit->wire->udir.y(); double sina=hit->wire->udir.x(); // Cathode Residual double u=rt->GetLastDistAlongWire(); double u_cathodes = hit->s; double resic = u - u_cathodes; // Get "measured" distance to wire. // For matching purposes this is assumed to be half a cell size double dist=0.25; // Take into account non-normal incidence to FDC plane double pz=last_step->mom.z(); double tx=last_step->mom.x()/pz; double ty=last_step->mom.y()/pz; double tu=tx*cosa-ty*sina; double alpha=atan(tu); double cosalpha=cos(alpha); // Compensate for the fact that the field at the "origin" of the // track does not correspond to the average Bz used to compute pt double Bz=last_step->B.z(); double Bratio=Bz/Bz0; double invBratio=1./Bratio; pt_over_a*=invBratio; p_over_a*=invBratio; a_over_p*=Bratio; // Anode Residual double resi = dist - doca/cosalpha; // Initialize some probability-related variables needed later double probability=0.,chisq=0.; //if (fit_type!=kHelical){ if (false){ // Correct for deflection of avalanche position due to Lorentz force double sign=(pos.x()*cosa-pos.y()*sina-hit->w)<0?1:-1.; double ucor=0.1458*Bz*(1.-0.048*last_step->B.Perp())*sign*doca; resic-=ucor; } else{ // Cathode variance due to Lorentz deflection double max_deflection=0.1458*Bz*(1.-0.048*last_step->B.Perp())*0.5; var_cathode=VAR_CATHODE_STRIPS+max_deflection*max_deflection/3.; } double var_tot=var_anode+var_cathode; // Variance in angles due to multiple scattering var_lambda = last_step->itheta02; var_phi=var_lambda*(1.+tanl2); // if (most_downstream_hit){ // // Fractional variance in the curvature k due to resolution and multiple scattering // double s_sq=s*s; // double var_k_over_k_sq_res=var_tot*p_over_a*p_over_a // *0.0720/double(N+4)/(s_sq*s_sq)/cosl2; // // // double one_over_beta=sqrt(1.+mass*mass/p_sq); // double var_pt_factor=0.016*one_over_beta/(cosl*0.003*last_step->B.z()); // double var_k_over_k_sq_ms=var_pt_factor*var_pt_factor*last_step->invX0/s; // // Fractional variance in pt // var_pt_over_pt_sq=var_k_over_k_sq_ms+var_k_over_k_sq_res; // // // Variance in dip angle due to measurement error // var_lambda_res=12.0*var_tot*double(N-1)/double(N*(N+1)) // *sinl2*sinl2/s_sq; // // most_downstream_hit=false; // } // Include error in lambda due to measurements var_lambda+=var_lambda_res; // Variance in position due to multiple scattering double var_pos_ms=last_step->itheta02s2/48.; // Variances in x and y due to uncertainty in track parameters double as_over_p=s*a_over_p; double sin_as_over_p=sin(as_over_p); double cos_as_over_p=cos(as_over_p); double one_minus_cos_as_over_p=1-cos_as_over_p; double diff1=sin_as_over_p-as_over_p*cos_as_over_p; double diff2=one_minus_cos_as_over_p-as_over_p*sin_as_over_p; double pdx_dp=pt_over_a*(cosphi*diff1-sinphi*diff2); double dx_ds=cosl*(cosphi*cos_as_over_p-sinphi*sin_as_over_p); double ds_dcosl=dz*cosl/(sinl*sinl2); double dx_dcosl =p_over_a*(cosphi*sin_as_over_p-sinphi*one_minus_cos_as_over_p) +dx_ds*ds_dcosl; double dx_dphi=-pt_over_a*(sinphi*sin_as_over_p+cosphi*one_minus_cos_as_over_p); double var_z0=2.*tanl2*(var_tot)*double(2*N-1)/double(N*(N+1)); double var_x=var_x0+pdx_dp*pdx_dp*var_pt_over_pt_sq+var_pos_ms +dx_dcosl*dx_dcosl*sinl2*var_lambda+dx_dphi*dx_dphi*var_phi +dx_ds*dx_ds*var_z0/sinl2; double pdy_dp=pt_over_a*(sinphi*diff1+cosphi*diff2); double dy_ds=cosl*(sinphi*cos_as_over_p+cosphi*sin_as_over_p); double dy_dcosl =p_over_a*(sinphi*sin_as_over_p+cosphi*one_minus_cos_as_over_p) +dy_ds*ds_dcosl; double dy_dphi=pt_over_a*(cosphi*sin_as_over_p-sinphi*one_minus_cos_as_over_p); double var_y=var_y0+pdy_dp*pdy_dp*var_pt_over_pt_sq+var_pos_ms +dy_dcosl*dy_dcosl*sinl2*var_lambda+dy_dphi*dy_dphi*var_phi +dy_ds*dy_ds*var_z0/sinl2; // Rotate from global coordinate system into FDC local system double cos2a=cosa*cosa; double sin2a=sina*sina; double var_d=(cos2a*var_x+sin2a*var_y)/(cosalpha*cosalpha); double var_u=cos2a*var_y+sin2a*var_x; //if (fit_type!=kHelical){ if (false){ // Factors take into account improved resolution after wire-based fit var_d*=0.1; var_u*=0.1; } // Calculate chisq chisq = resi*resi/(var_d+var_anode)+resic*resic/(var_u+var_cathode); // Probability of this hit being on the track probability = TMath::Prob(chisq,2); hp.doca = doca; hp.z = pos.z(); hp.dz = dz; hp.u = u; hp.u_cathodes = u_cathodes; hp.resic = resic; hp.pz = pz; hp.tx = tx; hp.ty = ty; hp.tu = tu; hp.alpha = alpha; hp.cosalpha = cosalpha; hp.Bratio = Bratio; hp.resi = resi; hp.chisq = chisq; hp.prob = probability; return true; } //------------------ // erun //------------------ jerror_t JEventProcessor_trkassassin::erun(void) { // This is called whenever the run number changes, before it is // changed to give you a chance to clean up before processing // events from the next run number. return NOERROR; } //------------------ // fini //------------------ jerror_t JEventProcessor_trkassassin::fini(void) { // Called before program exit after event processing is finished. return NOERROR; }