18 #include <JANA/JApplication.h>
19 #include <JANA/JEventLoop.h>
26 #include <PID/DParticle.h>
55 target.origin.SetXYZ(0.0, 0.0, 65.0);
56 target.sdir.SetXYZ(1.0, 0.0, 0.0);
57 target.tdir.SetXYZ(0.0, 1.0, 0.0);
58 target.udir.SetXYZ(0.0, 0.0, 1.0);
61 pthread_mutex_init(&mutex, NULL);
85 TDirectory *
dir = (TDirectory*)gROOT->FindObject(
"TRACKING");
86 if(!dir)dir =
new TDirectoryFile(
"TRACKING",
"TRACKING");
89 cdchits =
new TTree(
"cdchit",
"CDC hits");
90 cdchits->Branch(
"cdchit",
"dchit",&cdchit_ptr);
92 fdchits =
new TTree(
"fdchit",
"FDC hits");
93 fdchits->Branch(
"fdchit",
"dchit",&fdchit_ptr);
95 ttrack =
new TTree(
"track",
"Track");
96 ttrack->Branch(
"track",
"trackpar",&trk_ptr);
110 _DBG_<<
"Cannot get DApplication from JEventLoop! (are you using a JApplication based program perhaps?)"<<endl;
111 return RESOURCE_UNAVAILABLE;
133 sprintf(str,
"%3.4f%%", 100.0*(
double)NLRbad/(
double)(NLRbad+NLRgood));
135 cout<<endl<<setprecision(4);
136 cout<<
" NLRgood: "<<NLRgood<<endl;
137 cout<<
" NLRbad: "<<NLRbad<<endl;
138 cout<<
" NLRfit==0: "<<NLRfit_unknown<<endl;
139 cout<<
"Percentage bad: "<<str<<endl;
140 cout<<
" Nevents: "<<
Nevents<<endl;
151 vector<const DChargedTrack*> chargedtracks;
152 vector<const DTrackCandidate*> candidates;
153 vector<const DMCThrown*> mcthrowns;
154 vector<const DMCTrackHit*> mctrackhits;
156 loop->Get(chargedtracks);
157 loop->Get(candidates);
158 loop->Get(mcthrowns);
159 loop->Get(mctrackhits);
164 if(chargedtracks.size() <1 || mcthrowns.size() !=1 || candidates.size()<1 || chargedtracks[0]->hypotheses.size()<1){
165 _DBG_<<
"chargedtracks.size()="<<chargedtracks.size()<<
" mcthrowns.size()="<<mcthrowns.size()<<
" candidates.size()="<<candidates.size()<<endl;
175 for(
unsigned int i=0; i<chargedtracks.size(); i++){
176 if(chargedtracks[i]->hypotheses.size()<1)
continue;
177 DVector3 mom = chargedtracks[i]->hypotheses[0]->momentum();
179 if(mom.Mag()<1.0E-9)
continue;
183 double cos_phi = mom.Dot(mc_mom)/mom.Mag()/mc_mom.Mag();
184 if(cos_phi>1.)cos_phi=1.0;
185 if(cos_phi<-1.)cos_phi=-1.0;
187 double delta_pt = (mom.Pt()-mc_mom.Pt())/mc_mom.Pt();
188 double delta_theta = (mom.Theta() - mc_mom.Theta())*1000.0;
189 double delta_phi = acos(cos_phi)*1000.0;
190 double chisq = pow(delta_pt/0.04, 2.0) + pow(delta_theta/20.0, 2.0) + pow(delta_phi/20.0, 2.0);
194 recon = chargedtracks[i]->hypotheses[0];
199 vector<const DCDCTrackHit *> cdctrackhits;
200 vector<const DFDCPseudo *> fdcpseudohits;
201 recon->Get(cdctrackhits);
202 recon->Get(fdcpseudohits);
210 japp->RootWriteLock();
213 int NLRcorrect_this_track = 0;
214 int NLRincorrect_this_track = 0;
215 for(
unsigned int i=0; i<cdctrackhits.size(); i++){
223 hit_info.
FindLR(mctrackhits);
224 cdchit.eventnumber = eventnumber;
226 cdchit.layer = cdctrackhit->
wire->
ring;
227 cdchit.t = cdctrackhit->
tdrift;
228 cdchit.tof = hit_info.
tof;
229 cdchit.doca = hit_info.
doca;
230 cdchit.resi = hit_info.
dist - hit_info.
doca;
232 cdchit.u_pseudo = 0.0;
233 cdchit.u_lorentz = 0.0;
235 cdchit.trk_chisq = recon->
chisq;
236 cdchit.trk_Ndof = recon->
Ndof;
238 cdchit.LRfit = hit_info.
LRfit;
239 cdchit.pos_doca = hit_info.
pos_doca;
240 cdchit.pos_wire = hit_info.
pos_wire;
245 if(cdchit.LRis_correct){
246 NLRcorrect_this_track++;
249 NLRincorrect_this_track++;
260 for(
unsigned int i=0; i<fdcpseudohits.size(); i++){
261 const DFDCPseudo *fdcpseudohit = fdcpseudohits[i];
273 hit_info.
FindLR(mctrackhits, lorentz_def);
275 fdchit.eventnumber = eventnumber;
276 fdchit.wire = fdcpseudohit->
wire->
wire;
277 fdchit.layer = fdcpseudohit->
wire->
layer;
278 fdchit.t = fdcpseudohit->
time;
279 fdchit.tof = hit_info.
tof;
280 fdchit.doca = hit_info.
doca;
281 fdchit.resi = hit_info.
dist - hit_info.
doca;
282 fdchit.u = hit_info.
u;
283 fdchit.u_pseudo = fdcpseudohit->
s;
285 fdchit.resic = hit_info.
u - (fdcpseudohit->
s + hit_info.
u_lorentz) ;
286 fdchit.trk_chisq = recon->
chisq;
287 fdchit.trk_Ndof = recon->
Ndof;
289 fdchit.LRfit = hit_info.
LRfit;
290 fdchit.pos_doca = hit_info.
pos_doca;
291 fdchit.pos_wire = hit_info.
pos_wire;
296 if(fdchit.LRis_correct){
297 NLRcorrect_this_track++;
300 NLRincorrect_this_track++;
310 double doca_tgt = rt->
DistToRT(&target);
315 TVector3 thrown_mom(tmp.X(), tmp.Y(), tmp.Z());
317 TVector3 recon_mom(tmp.X(), tmp.Y(), tmp.Z());
319 TVector3 thrown_pos(tmp.X(), tmp.Y(), tmp.Z());
321 TVector3 can_pos(tmp.X(), tmp.Y(), tmp.Z());
324 trk.eventnumber = eventnumber;
325 trk.pthrown = thrown_mom;
326 trk.pfit = recon_mom;
327 trk.z_thrown = thrown_pos.Z();
328 trk.z_fit = tgt_doca.Z();
329 trk.z_can = can_pos.Z();
330 trk.r_fit = doca_tgt;
331 trk.NLRcorrect = NLRcorrect_this_track;
332 trk.NLRincorrect = NLRincorrect_this_track;
359 DVector3 pos_doca_dvec3(pos_doca.X(), pos_doca.Y(), pos_doca.Z());
360 DVector3 mom_doca_dvec3(mom_doca.X(), mom_doca.Y(), mom_doca.Z());
363 doca = rt->DistToRT(wire, &s);
364 rt->GetLastDOCAPoint(pos_doca_dvec3, mom_doca_dvec3);
365 DVector3 shift = wire->udir.Cross(mom_doca_dvec3);
366 u = rt->GetLastDistAlongWire();
367 DVector3 pos_wire_devc3 = wire->origin +
u*wire->udir;
368 pos_wire.SetXYZ(pos_wire_devc3.X(), pos_wire_devc3.Y(), pos_wire_devc3.Z());
371 double mass = 0.13957;
372 double beta = 1.0/
sqrt(1.0 + pow(mass/mom_doca.Mag(), 2.0))*2.998E10;
374 dist = (tdrift - tof)*55E-4;
378 DVector3 shift = wire->udir.Cross(mom_doca_dvec3);
380 double LRsign = shift.Dot(pos_doca_dvec3-pos_wire_devc3)<0.0 ? +1.0:-1.0;
381 double alpha = mom_doca.Angle(TVector3(0,0,1));
392 LRis_correct =
false;
394 TVector3 pos_truth(mctrackhit->
r*cos(mctrackhit->
phi), mctrackhit->
r*
sin(mctrackhit->
phi), mctrackhit->
z);
395 TVector3 pos_diff_truth = pos_truth-pos_wire;
396 TVector3 pos_diff_traj = pos_doca-pos_wire;
398 TVector3 shift_tvec3(shift.X(), shift.Y(), shift.Z());
399 LRfit = shift_tvec3.Dot(pos_diff_traj)<0.0 ? -1:1;
400 LRis_correct = pos_diff_truth.Dot(pos_diff_traj)>0.0;
float chisq
Chi-squared for the track (not chisq/dof!)
sprintf(text,"Post KinFit Cut")
DVector3 GetLastDOCAPoint(void) const
const DVector3 & position(void) const
const DFDCWire * wire
DFDCWire for this wire.
DLorentzDeflections * GetLorentzDeflections(unsigned int run_number=1)
class DFDCPseudo: definition for a reconstructed point in the FDC
static const DMCTrackHit * GetMCTrackHit(const DCoordinateSystem *wire, double rdrift, vector< const DMCTrackHit * > &mctrackhits, int trackno_filter=-1)
~DEventProcessor_track_hists()
jerror_t brun(JEventLoop *loop, int32_t runnumber)
double DistToRT(double x, double y, double z) const
const DCoordinateSystem * wire
float z
coordinates of hit in cm and rad
jerror_t evnt(JEventLoop *loop, uint64_t eventnumber)
Invoked via DEventProcessor virtual method.
double s
< wire position computed from cathode data, assuming the avalanche occurs at the wire ...
void FindLR(vector< const DMCTrackHit * > &mctrackhits, const DLorentzDeflections *lorentz_def=NULL)
int Ndof
Number of degrees of freedom in the fit.
<A href="index.html#legend"> <IMG src="CORE.png" width="100"> </A>
double time
time corresponding to this pseudopoint.
jerror_t fini(void)
Invoked via DEventProcessor virtual method.
double GetLorentzCorrection(double x, double y, double z, double alpha, double dx) const
DReferenceTrajectory * rt
const DVector3 & momentum(void) const
jerror_t erun(void)
Invoked via DEventProcessor virtual method.
jerror_t init(void)
Invoked via DEventProcessor virtual method.
DEventProcessor_track_hists()