12 #include <JANA/JApplication.h>
13 #include <JANA/JFactory.h>
46 double z_lower_limit = 38.5;
47 double z_upper_limit = 98.5;
49 gPARMS->SetParameter(
"TRKFIT:USE_SC_TIME",
false);
51 if(gPARMS->Exists(
"TRKFIT:USE_SC_TIME"))
52 gPARMS->GetParameter(
"TRKFIT:USE_SC_TIME", USE_SC_TIME);
58 cout <<
"=========================================================================="<< endl;
59 cout <<
"TRKFIT: USE_SC_TIME = 0; WARNING SC TIME WILL NOT BE USED IN TRACK FITTING"<< endl;
60 cout <<
"This is required in ST_ZEff plugin which calculate SC efficiency "<< endl;
61 cout <<
"=========================================================================="<< endl;
64 TDirectory *
main = gDirectory;
65 gDirectory->mkdir(
"st_efficiency")->cd();
67 h_fom =
new TH1D(
"h_fom",
" FOM ; FOM ; Counts", 1000,0.0, 1.0);
68 h_N_Hit_in_track =
new TH1D(
"h_N_Hit_in_track",
" Ndof + 5; Number of Hits per track ; Counts", 70,0.0, 70.0);
69 h1_qVectorSize=
new TH1D(
" h1_qVectorSize",
" number of q tracks; number of q tracks; Counts", 50,0.0, 50.0);
70 h1_qVectorSize_ACuts=
new TH1D(
" h1_qVectorSize_ACuts",
" number of q tracks; number of q tracks; Counts", 50,0.0, 50.0);
71 h1_tDiff=
new TH1D(
"h1_tDiff",
" SC_time - RF; t(ns) ; Counts", 600,-30.0, 30.0);
72 h1_RFtime=
new TH1D(
"h1_RFtime",
" t0; t0(ns) ; Counts", 600,-30.0, 30.0);
73 h1_Centered_RFtime=
new TH1D(
"h1_Centered_RFtime",
" RF_C; RF_C(ns) ; Counts", 600,-30.0, 30.0);
74 h1_SC_ShiftedTime=
new TH1D(
"h1_SC_ShiftedTime",
" SC_shifted Time; Time(ns) ; Counts", 600,-30.0, 30.0);
76 h_z_v =
new TH1D(
"h_z_v",
" z_v; Z (cm); Counts", 1000,0.0, 100.0);
77 h1_st_pred_id =
new TH1D(
"h1_st_prd_id",
"h1_st_prd_id; Sector; Predicted Hit Counts", 31, -0.5, 30.5);
79 h2_z_vs_r =
new TH2I(
"h2_z_vs_r",
"Z vs R; Z (cm); R (cm)", NoBins_z, z_lower_limit, z_upper_limit, 100, 0.0, 10.0);
80 h2_x_vs_y =
new TH2I (
"h2_x_vs_y",
"X vs Y; X (cm); Y (cm)", 1000, -20., 20., 500, -10., 10.);
85 h_N_trck_prd_z_ss[i] =
new TH1D(Form(
"h_N_trck_prd_z_ss_%i",i+1), Form(
"z #in interval %i; Sector; Predicted Hit Counts", i+1), 31, -0.5, 30.5);
86 h_N_recd_hit_z_ss[i] =
new TH1D(Form(
"h_N_recd_hit_z_ss_%i",i+1), Form(
"z #in interval %i; Sector; Recorded Hit Counts", i+1), 31, -0.5, 30.5);
87 h_N_recd_hit_z_ss_ACC[i] =
new TH1D(Form(
"h_N_recd_hit_z_ss_ACC_%i",i+1), Form(
"z #in interval %i; Sector; ACC Hit Counts", i+1), 31, -0.5, 30.5);
88 h_N_trck_prd_z_ss_eff[i] =
new TH1D(Form(
"h_N_trck_prd_z_ss_eff_%i",i+1), Form(
"z #in interval %i; Sector; SS Effeiciency ", i+1), 31, -0.5, 30.5);
92 h_N_trck_prd_z_bs[i] =
new TH1D(Form(
"h_N_trck_prd_z_bs_%i",i+1), Form(
"z #in interval %i; Sector; Predicted Hit Counts", i+1), 31, -0.5, 30.5);
93 h_N_recd_hit_z_bs[i] =
new TH1D(Form(
"h_N_recd_hit_z_bs_%i",i+1), Form(
"z #in interval %i; Sector; Recorded Hit Counts", i+1), 31, -0.5, 30.5);
94 h_N_recd_hit_z_bs_ACC[i] =
new TH1D(Form(
"h_N_recd_hit_z_bs_ACC_%i",i+1), Form(
"z #in interval %i; Sector; ACC Hit Counts", i+1), 31, -0.5, 30.5);
95 h_N_trck_prd_z_bs_eff[i] =
new TH1D(Form(
"h_N_trck_prd_z_bs_eff_%i",i+1), Form(
"z #in interval %i; Sector; BS Effeiciency ", i+1), 31, -0.5, 30.5);
99 h_N_trck_prd_z_ns[i] =
new TH1D(Form(
"h_N_trck_prd_z_ns_%i",i+1), Form(
"z #in interval %i; Sector; Predicted Hit Counts", i+1), 31, -0.5, 30.5);
100 h_N_recd_hit_z_ns[i] =
new TH1D(Form(
"h_N_recd_hit_z_ns_%i",i+1), Form(
"z #in interval %i; Sector; Recorded Hit Counts", i+1), 31, -0.5, 30.5);
101 h_N_recd_hit_z_ns_ACC[i] =
new TH1D(Form(
"h_N_recd_hit_z_ns_ACC_%i",i+1), Form(
"z #in interval %i; Sector; ACC Hit Counts", i+1), 31, -0.5, 30.5);
102 h_N_trck_prd_z_ns_eff[i] =
new TH1D(Form(
"h_N_trck_prd_z_ns_eff_%i",i+1), Form(
"z #in interval %i; Sector; NS Effeiciency ", i+1), 31, -0.5, 30.5);
106 gDirectory->cd(
"../");
118 map<string,double> target_params;
119 if (eventLoop->GetCalib(
"/TARGET/target_parms", target_params))
120 jout <<
"Error loading /TARGET/target_parms/ !" << endl;
121 if (target_params.find(
"TARGET_Z_POSITION") != target_params.end())
122 z_target_center = target_params[
"TARGET_Z_POSITION"];
124 jerr <<
"Unable to get TARGET_Z_POSITION from /TARGET/target_parms !" << endl;
128 _DBG_<<
"Cannot get DApplication from JEventLoop! (are you using a JApplication based program?)"<<endl;
133 double theta = sc_norm[0][sc_norm[0].size()-2].Theta();
134 sc_angle_corr = 1./cos(M_PI_2 - theta);
138 if(eventLoop->GetCalib(
"START_COUNTER/propagation_time_corr", propagation_time_corr))
139 jout <<
"Error loading /START_COUNTER/propagation_time_corr !" << endl;
141 if(eventLoop->GetCalib(
"START_COUNTER/PTC_Boundary", PTC_Boundary))
142 jout <<
"Error loading /START_COUNTER/PTC_Boundary !" << endl;
162 vector<const DSCHit*> st_hits;
163 vector<const DChargedTrack*> chargedTrackVector;
165 eventLoop->Get(st_hits);
166 eventLoop->Get(chargedTrackVector);
167 vector<DVector3> sc_track_position;
170 eventLoop->GetSingle(locTrigger);
173 double speed_light = 29.9792458;
175 eventLoop->GetSingle(locEventRFBunch);
180 vector<const DParticleID* > locParticleID_algos;
181 eventLoop->Get(locParticleID_algos);
182 if(locParticleID_algos.size() < 1)
184 _DBG_<<
"Unable to get a DParticleID object! NO PID will be done!"<<endl;
185 return RESOURCE_UNAVAILABLE;
187 auto locParticleID = locParticleID_algos[0];
194 eventLoop->GetSingle(locDetectorMatches);
196 japp->RootFillLock(
this);
197 sc_track_position.clear();
198 for (uint32_t i = 0; i < chargedTrackVector.size(); i++)
201 const DChargedTrack *thisChargedTrack = chargedTrackVector[i];
207 float trackingFOMCut = 2.699796E-3;
208 float NHitsTrack = 14.0;
214 if(tb_track->
Ndof + 5 < NHitsTrack)
continue;
215 if(tb_track->
FOM < trackingFOMCut)
continue;
217 shared_ptr<const DSCHitMatchParams> locBestSCHitMatchParams;
218 bool foundSC = locParticleID->Get_BestSCMatchParams(tb_track, locDetectorMatches, locBestSCHitMatchParams);
219 if (!foundSC)
continue;
227 double z_v = vertex.z();
228 double r_v = vertex.Perp();
229 bool z_vertex_cut = fabs(z_target_center - z_v) <= 15.0;
230 bool r_vertex_cut = r_v < 1.0;
233 if (!z_vertex_cut)
continue;
234 if (!r_vertex_cut)
continue;
238 shared_ptr<const DTOFHitMatchParams> locTOFHitMatchParams;
239 bool foundTOF = locParticleID->Get_BestTOFMatchParams(tb_track, locDetectorMatches, locTOFHitMatchParams);
242 shared_ptr<const DBCALShowerMatchParams> locBCALHitMatchParams;
243 bool foundBCAL = locParticleID->Get_BestBCALMatchParams(tb_track, locDetectorMatches, locBCALHitMatchParams);
245 shared_ptr<const DFCALShowerMatchParams> locFCALHitMatchParams;
246 bool foundFCAL = locParticleID->Get_BestFCALMatchParams(tb_track, locDetectorMatches, locFCALHitMatchParams);
248 if (!(foundBCAL || (foundFCAL && foundTOF)))
continue;
252 vector<shared_ptr<const DSCHitMatchParams>> st_params;
255 if (!st_match)
continue;
257 DVector3 IntersectionPoint, IntersectionMomentum,locProjMom, locPaddleNorm;
258 shared_ptr<DSCHitMatchParams> locSCHitMatchParams;
260 bool st_match_pid = locParticleID->Cut_MatchDistance(extrapolations, st_params[0]->dSCHit, st_params[0]->dSCHit->t, locSCHitMatchParams,
true, &IntersectionPoint, &IntersectionMomentum);
262 if(!st_match_pid)
continue;
264 int st_pred_id=locBestSCHitMatchParams->dSCHit->sector;
265 double t0 = locEventRFBunch->
dTime;
269 int st_pred_id_index = st_pred_id - 1;
272 sc_track_position.push_back(IntersectionPoint);
276 double locCenteredRFTime = t0;
279 double locCenterToVertexRFTime = (tb_track->
z() - z_target_center)*(1.0/speed_light);
280 double locVertexRFTime = locCenteredRFTime + locCenterToVertexRFTime;
282 int sc_index = st_pred_id_index;
284 double sc_pos_soss, sc_pos_eoss, sc_pos_eobs, sc_pos_eons;
286 sc_pos_soss = sc_pos[sc_index][0].z();
287 sc_pos_eoss = sc_pos[sc_index][1].z();
288 sc_pos_eobs = sc_pos[sc_index][11].z();
289 sc_pos_eons = sc_pos[sc_index][12].z();
292 incpt_ss = propagation_time_corr[sc_index][0];
293 slope_ss = propagation_time_corr[sc_index][1];
294 incpt_bs = propagation_time_corr[sc_index][2];
295 slope_bs = propagation_time_corr[sc_index][3];
296 incpt_ns = propagation_time_corr[sc_index][4];
297 slope_ns = propagation_time_corr[sc_index][5];
299 Bound1 = PTC_Boundary[0][0];
300 Bound2 = PTC_Boundary[1][0];
301 for (uint32_t i = 0; i < sc_track_position.size(); i++)
304 h_z_v->Fill(sc_track_position[i].z());
305 h2_z_vs_r->Fill(sc_track_position[i].z(),sc_track_position[i].Perp());
306 h2_x_vs_y->Fill(sc_track_position[i].
x(), sc_track_position[i].
y());
307 double locSCzIntersection = sc_track_position[i].z();
316 z_ss[k] = sc_pos_soss + ss_interval * k ;
318 if ((z_ss[k] <= locSCzIntersection) && (locSCzIntersection < (z_ss[k]+ss_interval)))
323 for (uint32_t j = 0; j < st_hits.size(); j++)
325 int hit_sector = st_hits[j]->sector;
326 if ((st_pred_id == hit_sector)||(st_pred_id == hit_sector-1) || (st_pred_id == hit_sector+1) ||(st_pred_id == hit_sector+29)||(st_pred_id == hit_sector-29))
329 double FlightTime = locBestSCHitMatchParams->dFlightTime;
331 double st_corr_FlightTime = st_hits[j]->t - FlightTime;
333 double pathlength = locSCzIntersection - sc_pos_soss;
334 double corr_time = st_corr_FlightTime - (incpt_ss + (slope_ss * pathlength));
337 double t_diff =corr_time - t0;
341 if (!((-10 < t_diff) && (t_diff <10)) && (-20 < t_diff) && (t_diff < 20))
346 if ((-10 < t_diff) && (t_diff <10))
358 z_bs[p] = sc_pos_eoss + bs_interval * p ;
360 if ((z_bs[p] <= locSCzIntersection) && (locSCzIntersection < (z_bs[p]+bs_interval)))
364 for (uint32_t j = 0; j < st_hits.size(); j++)
366 int hit_sector = st_hits[j]->sector;
367 if ((st_pred_id == hit_sector)||(st_pred_id == hit_sector-1) || (st_pred_id == hit_sector+1) ||(st_pred_id == hit_sector+29)||(st_pred_id == hit_sector-29))
371 double FlightTime = locBestSCHitMatchParams->dFlightTime;
373 double st_corr_FlightTime = st_hits[j]->t - FlightTime;
374 double SS_Length = sc_pos_eoss - sc_pos_soss;
376 double path_bs = SS_Length + (locSCzIntersection - sc_pos_eoss)*sc_angle_corr;
377 double corr_time_bs= st_corr_FlightTime - (incpt_bs + (slope_bs * path_bs));
380 double t_diff = corr_time_bs - t0;
387 if (!((-10 < t_diff) && (t_diff <10)) && (-20 < t_diff) && (t_diff < 20))
392 if ((-10 < t_diff) && (t_diff <10))
403 z_ns[k] = sc_pos_eobs + ns_interval * k ;
405 if ((z_ns[k] <= locSCzIntersection) && (locSCzIntersection <= (z_ns[k]+ns_interval)))
408 for (uint32_t j = 0; j < st_hits.size(); j++)
410 int hit_sector = st_hits[j]->sector;
411 if ((st_pred_id == hit_sector)||(st_pred_id == hit_sector-1) || (st_pred_id == hit_sector+1) ||(st_pred_id == hit_sector+29)||(st_pred_id == hit_sector-29))
414 double FlightTime = locBestSCHitMatchParams->dFlightTime;
415 double SS_Length = sc_pos_eoss - sc_pos_soss;
417 double st_corr_FlightTime = st_hits[j]->t - FlightTime;
419 double path_ns = SS_Length + (locSCzIntersection - sc_pos_eoss)*sc_angle_corr;
420 if (path_ns <= Bound1)
422 Corr_Time_ns= st_corr_FlightTime - (incpt_ss + (slope_ss * path_ns));
425 else if ((Bound1 < path_ns)&&(path_ns <= Bound2))
427 Corr_Time_ns = st_corr_FlightTime - (incpt_bs + (slope_bs * path_ns));
432 Corr_Time_ns = st_corr_FlightTime - (incpt_ns + (slope_ns * path_ns));
435 double t_diff =Corr_Time_ns - t0;
437 if (!((-10 < t_diff) && (t_diff <10)) && (-20 < t_diff) && (t_diff < 20))
442 if ((-10 < t_diff) && (t_diff <10))
454 japp->RootFillUnLock(
this);
static TH1D * h_N_trck_prd_z_ns_eff[Nof_ns_intervals]
static TH1D * h_N_recd_hit_z_bs[Nof_bs_intervals]
jerror_t brun(jana::JEventLoop *eventLoop, int32_t runnumber)
Called everytime a new run number is detected.
DRFTime_factory * locRFTimeFactory
static TH1D * h_N_recd_hit_z_ss[Nof_ss_intervals]
static TH1D * h_N_trck_prd_z_bs_eff[Nof_bs_intervals]
bool Get_SCMatchParams(const DTrackingData *locTrack, vector< shared_ptr< const DSCHitMatchParams > > &locMatchParams) const
const uint32_t Nof_ss_intervals
jerror_t fini(void)
Called after last event of last event source has been processed.
const DChargedTrackHypothesis * Get_BestTrackingFOM(void) const
static TH1D * h_N_Hit_in_track
const DVector3 & position(void) const
const DTrackTimeBased * Get_TrackTimeBased(void) const
static TH1D * h1_qVectorSize
uint32_t Get_L1FrontPanelTriggerBits(void) const
static TH1D * h_N_trck_prd_z_ss[Nof_ss_intervals]
static TH1D * h1_qVectorSize_ACuts
static TH1D * h_N_recd_hit_z_ss_ACC[Nof_ss_intervals]
jerror_t evnt(jana::JEventLoop *eventLoop, uint64_t eventnumber)
Called every event.
static TH1D * h_N_trck_prd_z_ns[Nof_ns_intervals]
JEventProcessor_ST_ZEff()
static TH1D * h1_Centered_RFtime
DGeometry * GetDGeometry(unsigned int run_number)
static TH1D * h_N_trck_prd_z_bs[Nof_bs_intervals]
int Ndof
Number of degrees of freedom in the fit.
static TH1D * h_N_recd_hit_z_ns[Nof_ns_intervals]
const uint32_t Nof_ns_intervals
static TH1D * h_N_recd_hit_z_bs_ACC[Nof_bs_intervals]
const DVector3 & momentum(void) const
jerror_t init(void)
Called once at program start.
map< DetectorSystem_t, vector< DTrackFitter::Extrapolation_t > > extrapolations
unsigned int dNumParticleVotes
double Step_TimeToNearInputTime(double locTimeToStep, double locTimeToStepTo) const
jerror_t erun(void)
Called everytime run number changes, provided brun has been called.
static TH1D * h_N_recd_hit_z_ns_ACC[Nof_ns_intervals]
~JEventProcessor_ST_ZEff()
static TH1D * h1_SC_ShiftedTime
static TH1D * h1_st_pred_id
int main(int argc, char *argv[])
const uint32_t Nof_bs_intervals
bool GetStartCounterGeom(vector< vector< DVector3 > > &pos, vector< vector< DVector3 > > &norm) const
static TH1D * h_N_trck_prd_z_ss_eff[Nof_ss_intervals]