19 #include <JANA/JApplication.h>
20 #include <JANA/JEventLoop.h>
46 pthread_mutex_init(&mutex, NULL);
63 TDirectory *
dir = (TDirectory*)gROOT->FindObject(
"TRACKING");
64 if(!dir)dir =
new TDirectoryFile(
"TRACKING",
"TRACKING");
67 ttrkres =
new TTree(
"track",
"Track Res.");
68 ttrkres->Branch(
"E",
"trackres",&trkres_ptr);
81 double locSIGMA_CDC = 0.0150;
82 double locSIGMA_FDC_ANODE = 0.0200;
83 double locSIGMA_FDC_CATHODE = 0.0200;
85 gPARMS->SetDefaultParameter(
"TRKFIT:SIGMA_CDC", locSIGMA_CDC);
86 gPARMS->SetDefaultParameter(
"TRKFIT:SIGMA_FDC_ANODE", locSIGMA_FDC_ANODE);
87 gPARMS->SetDefaultParameter(
"TRKFIT:SIGMA_FDC_CATHODE", locSIGMA_FDC_CATHODE);
91 pthread_mutex_lock(&mutex);
95 SIGMA_CDC = locSIGMA_CDC;
96 SIGMA_FDC_ANODE = locSIGMA_FDC_ANODE;
97 SIGMA_FDC_CATHODE = locSIGMA_FDC_CATHODE;
99 pthread_mutex_unlock(&mutex);
109 vector<const DMCTrajectoryPoint*> trajpoints;
110 vector<const DCDCTrackHit*> cdchits;
111 vector<const DFDCPseudo*> fdchits;
114 loop->Get(trajpoints);
117 loop->GetSingle(mcthrown);
125 for(
unsigned int i=0; i<cdchits.size(); i++){
133 m.
traj = FindTrajectoryPoint(wire, m.
radlen, m.
s, trajpoints);
136 m.
B =
sqrt(Bx*Bx + By*By + Bz*Bz);
142 for(
unsigned int i=0; i<fdchits.size(); i++){
146 m.
err = SIGMA_FDC_ANODE;
150 m.
traj = FindTrajectoryPoint(wire, m.
radlen, m.
s, trajpoints);
153 m.
B =
sqrt(Bx*Bx + By*By + Bz*Bz);
158 if(meas.size()<5)
return NOERROR;
160 double deltak, pt_res;
161 GetPtRes(meas, deltak, pt_res);
164 GetThetaRes(meas, theta_res);
166 double pt =
sqrt(pow((
double)meas[0].traj->px,2.0) + pow((
double)meas[0].traj->py,2.0));
167 double p =
sqrt(pow((
double)meas[0].traj->pz,2.0) + pow(pt,2.0));
168 double theta = asin(pt/p);
169 if(theta<0.0)theta+=2.0*M_PI;
174 japp->RootWriteLock();
176 trkres.event = eventnumber;
177 trkres.recon.SetXYZ(meas[0].traj->px, meas[0].traj->py, meas[0].traj->pz);
178 trkres.thrown.SetXYZ(dthrown.X(), dthrown.Y(), dthrown.Z());
179 trkres.deltak = deltak;
180 trkres.pt_res = pt_res;
181 trkres.p_res = (pt_res*pt + fabs(pt/tan(theta)*theta_res))/
sin(theta)/p;
182 trkres.theta_res = theta_res;
199 double min_dist = 1.0E6;
202 double best_radlen=0.0;
204 for(
unsigned int i=0; i<trajpoints.size(); i++){
207 double d_dot_udir = d.Dot(wire->
udir);
208 double dist =
sqrt(d.Mag2() - d_dot_udir*d_dot_udir);
209 double dist_past_end = fabs(d_dot_udir) - wire->
L/2.0;
210 if(dist_past_end>0.0) dist =
sqrt(dist*dist + dist_past_end*dist_past_end);
213 radlen += (double)traj->
step/(
double)traj->
radlen;
219 best_radlen = radlen;
225 radlen = best_radlen;
242 for(
unsigned int i=0; i<meas.size(); i++){
243 double s = meas[i].s - meas[0].s;
244 double e = meas[i].err;
250 Vs4 += s*s*s*s/(e*
e);
253 double Vss = (Vs2 - Vs1*Vs1)/w;
254 double Vs2s2 = (Vs4 - Vs2*Vs2)/w;
255 double Vss2 = (Vs3 - Vs2*Vs1)/w;
257 double deltak_res =
sqrt(4.0/w*Vss/(Vss*Vs2s2 - Vss2*Vss2));
260 double radlen = meas[meas.size()-1].radlen;
261 double L = meas[meas.size()-1].s - meas[0].s;
262 double pt =
sqrt(pow((
double)meas[0].traj->px,2.0) + pow((
double)meas[0].traj->py,2.0));
263 double p =
sqrt(pow((
double)meas[0].traj->pz,2.0) + pow(pt,2.0));
264 double beta = p/
sqrt(p*p + pow(0.13957,2.0));
265 double lambda = M_PI/2.0 - asin(pt/p);
266 double deltak_ms = 0.016*
sqrt(radlen)/(L*p*beta*pow(cos(lambda), 2.0));
268 deltak =
sqrt(deltak_res*deltak_res + deltak_ms*deltak_ms);
272 Bavg/=(double)meas.size();
273 double k = 0.003*Bavg/pt;
291 double err2_avg = 0.0;
292 for(
unsigned int i=0; i<meas.size(); i++){
293 double s = meas[i].s - meas[0].s;
294 double e = meas[i].err;
299 s_avg/=(double)meas.size();
300 err2_avg/=(double)(meas.size()-2);
303 for(
unsigned int i=0; i<meas.size(); i++){
304 double s = meas[i].s - meas[0].s;
306 sum += pow(s - s_avg, 2.0);
319 double dtheta_res =
sqrt(err2_avg/sum);
322 double radlen = meas[meas.size()-1].radlen;
323 double pt =
sqrt(pow((
double)meas[0].traj->px,2.0) + pow((
double)meas[0].traj->py,2.0));
324 double p =
sqrt(pow((
double)meas[0].traj->pz,2.0) + pow(pt,2.0));
325 double beta = p/
sqrt(p*p + pow(0.13957,2.0));
326 double dtheta_ms = 0.0136*
sqrt(radlen)/(beta*p)*(1.0+0.038*log(radlen))/
sqrt(3.0);
328 theta_res =
sqrt(dtheta_res*dtheta_res + dtheta_ms*dtheta_ms);
jerror_t evnt(JEventLoop *eventLoop, uint64_t eventnumber)
Called every event.
void GetPtRes(vector< meas_t > &meas, double &deltak, double &pt_res)
jerror_t fini(void)
Called after last event of last event source has been processed.
DMagneticFieldMap * GetBfield(unsigned int run_number=1)
jerror_t erun(void)
Called everytime run number changes, provided brun has been called.
jerror_t init(void)
Called once at program start.
const DMCTrajectoryPoint * traj
const DMCTrajectoryPoint * FindTrajectoryPoint(const DCoordinateSystem *wire, double &radlen, double &s, vector< const DMCTrajectoryPoint * > trajpoints)
void GetThetaRes(vector< meas_t > &meas, double &theta_res)
const DVector3 & momentum(void) const
jerror_t brun(JEventLoop *eventLoop, int32_t runnumber)
Called everytime a new run number is detected.
~DEventProcessor_trkres_tree()
DEventProcessor_trkres_tree()