20 #ifdef PROFILE_TRK_TIMES
22 static map<string, prof_time::time_diffs> prof_times;
35 fit_status = kFitNotDone;
36 unsigned int run_number = (loop->GetJEvent()).GetRunNumber();
39 CORRECT_FOR_ELOSS=
true;
40 gPARMS->SetDefaultParameter(
"TRKFIT:CORRECT_FOR_ELOSS",CORRECT_FOR_ELOSS);
45 _DBG_<<
"Cannot get DApplication from JEventLoop! (are you using a JApplication based program?)"<<endl;
52 MATERIAL_MAP_MODEL =
"DGeometry";
53 gPARMS->SetDefaultParameter(
"TRKFIT:MATERIAL_MAP_MODEL",MATERIAL_MAP_MODEL);
54 if(MATERIAL_MAP_MODEL==
"DRootGeom"){
58 vector<Extrapolation_t>myvector;
59 extrapolations.emplace(
SYS_BCAL,myvector);
60 extrapolations.emplace(
SYS_TOF,myvector);
61 extrapolations.emplace(
SYS_FCAL,myvector);
62 extrapolations.emplace(
SYS_FDC,myvector);
63 extrapolations.emplace(
SYS_CDC,myvector);
64 extrapolations.emplace(
SYS_START,myvector);
65 extrapolations.emplace(
SYS_DIRC,myvector);
67 extrapolations[
SYS_TOF].reserve(1);
68 extrapolations[
SYS_BCAL].reserve(300);
70 extrapolations[
SYS_FDC].reserve(24);
71 extrapolations[
SYS_CDC].reserve(200);
77 #ifdef PROFILE_TRK_TIMES
80 prof_times[
"Ntracks"] = tdiff_zero;
96 #ifdef PROFILE_TRK_TIMES
102 fit_type = kWireBased;
105 cdchits_used_in_fit.clear();
106 fdchits_used_in_fit.clear();
107 ClearExtrapolations();
110 fit_status = kFitNotDone;
112 #ifdef PROFILE_TRK_TIMES
122 cdchits.push_back(cdchit);
123 fit_status = kFitNotDone;
131 for(
unsigned int i=0; i<cdchits.size(); i++)this->cdchits.push_back(cdchits[i]);
132 fit_status = kFitNotDone;
140 fdchits.push_back(fdchit);
141 fit_status = kFitNotDone;
149 for(
unsigned int i=0; i<fdchits.size(); i++)this->fdchits.push_back(fdchits[i]);
150 fit_status = kFitNotDone;
158 #ifdef PROFILE_TRK_TIMES
161 input_params.setPosition(pos);
162 input_params.setMomentum(mom);
163 input_params.setPID(
IDTrack(q, mass));
164 input_params.setTime(t0);
165 input_params.setT0(t0,0.,t0_det);
169 #ifdef PROFILE_TRK_TIMES
181 #ifdef PROFILE_TRK_TIMES
185 SetInputParameters(starting_params);
188 #ifdef PROFILE_TRK_TIMES
199 const map<
DetectorSystem_t,vector<DTrackFitter::Extrapolation_t> >&extrapolations,
201 double mass,
int N,
double t0,
206 fit_type = save_type;
209 if(mass<0.0)mass = starting_params.
mass();
211 double q=starting_params.
charge();
214 vector<const DTrackHitSelector *> hitselectors;
215 loop->Get(hitselectors);
216 if(hitselectors.size()<1){
217 _DBG_<<
"Unable to get a DTrackHitSelector object! NO Charged track fitting will be done!"<<endl;
218 return fit_status = kFitNotDone;
223 vector<const DCDCTrackHit*> cdctrackhits;
224 vector<const DFDCPseudo*> fdcpseudos;
225 loop->Get(cdctrackhits);
226 loop->Get(fdcpseudos);
231 if (extrapolations.at(
SYS_CDC).size()>0){
232 vector<Extrapolation_t>extraps=extrapolations.at(
SYS_CDC);
233 DVector3 mypos=extraps[extraps.size()/2].position;
234 double Bz=GetDMagneticFieldMap()->GetBz(mypos.x(),mypos.y(),mypos.z());
235 hitselector->
GetCDCHits(Bz,q,extraps,cdctrackhits,
this,N);
238 if (extrapolations.at(
SYS_FDC).size()>0){
239 vector<Extrapolation_t>extraps=extrapolations.at(
SYS_FDC);
240 DVector3 mypos=extraps[extraps.size()/2].position;
241 double Bz=GetDMagneticFieldMap()->GetBz(mypos.x(),mypos.y(),mypos.z());
242 hitselector->
GetFDCHits(Bz,q,extraps,fdcpseudos,
this,N);
245 if (got_hits==
false){
246 return fit_status = kFitNotDone;
250 fit_params.setPID(
IDTrack(q, mass));
252 #ifdef PROFILE_TRK_TIMES
253 start_time.TimeDiffNow(prof_times,
"Find Hits");
259 fit_status = FitTrack(pos, mom,q, mass,t0,t0_det);
261 #ifdef PROFILE_TRK_TIMES
262 start_time.TimeDiffNow(prof_times,
"Find Hits and Fit Track");
273 double mass,
int N,
double t0,
286 #ifdef PROFILE_TRK_TIMES
287 prof_times[
"Ntracks"].real += 1.0;
295 fit_type = save_type;
298 if(mass<0.0)mass = starting_params.
mass();
314 double q=starting_params.
charge();
321 vector<const DTrackHitSelector *> hitselectors;
322 loop->Get(hitselectors);
323 if(hitselectors.size()<1){
324 _DBG_<<
"Unable to get a DTrackHitSelector object! NO Charged track fitting will be done!"<<endl;
325 return fit_status = kFitNotDone;
330 vector<const DCDCTrackHit*> cdctrackhits;
331 vector<const DFDCPseudo*> fdcpseudos;
332 loop->Get(cdctrackhits);
333 loop->Get(fdcpseudos);
335 hitselector->
GetAllHits(input_type, rt, cdctrackhits, fdcpseudos,
this,N);
340 if (fdchits.size()+cdchits.size()==0)
return fit_status=kFitNotDone;
343 fit_params.setPID(
IDTrack(q, mass));
345 #ifdef PROFILE_TRK_TIMES
350 fit_status = FitTrack(pos, mom,q, mass,t0,t0_det);
352 #ifdef PROFILE_TRK_TIMES
353 start_time.
TimeDiffNow(prof_times,
"Find Hits and Fit Track");
366 vector<const DCDCTrackHit*> cdchits;
367 starting_params.Get(cdchits);
368 if(cdchits.size()>0){
370 first_wire = cdchits[cdchits.size()/2]->wire;
372 vector<const DFDCPseudo*> fdchits;
373 starting_params.Get(fdchits);
374 if(fdchits.size()!=0){
376 first_wire = fdchits[fdchits.size()/2]->wire;
381 return RESOURCE_UNAVAILABLE;
388 pos,mom,starting_params.
charge(), 1000.0, first_wire);
392 target.
origin.SetXYZ(0.0, 0.0, 65.0);
393 target.
sdir.SetXYZ(1.0, 0.0, 0.0);
394 target.
tdir.SetXYZ(0.0, 1.0, 0.0);
395 target.
udir.SetXYZ(0.0, 0.0, 1.0);
413 double Z_over_A,
double I){
414 double rho_Z_over_A=density*Z_over_A;
416 return CalcDensityEffect(p,mass,rho_Z_over_A,LnI);
423 double rho_Z_over_A,
double LnI){
425 double betagamma=p/mass;
426 return CalcDensityEffect(betagamma,rho_Z_over_A,LnI);
431 double rho_Z_over_A,
double LnI){
432 double X=log10(betagamma);
434 double Cbar=2.*(LnI-log(28.816
e-9*
sqrt(rho_Z_over_A)))+1.;
435 if (rho_Z_over_A>0.01){
437 if (Cbar<=3.681) X0=0.2;
438 else X0=0.326*Cbar-1.;
442 if (Cbar<=5.215) X0=0.2;
443 else X0=0.326*Cbar-1.5;
449 if (Cbar<=9.5) X0=1.6;
450 else if (Cbar>9.5 && Cbar<=10.) X0=1.7;
451 else if (Cbar>10 && Cbar<=10.5) X0=1.8;
452 else if (Cbar>10.5 && Cbar<=11.) X0=1.9;
453 else if (Cbar>11.0 && Cbar<=12.25) X0=2.;
454 else if (Cbar>12.25 && Cbar<=13.804){
465 delta=4.606*X-Cbar+(Cbar-4.606*
X0)*pow((X1-X)/(X1-
X0),3.);
474 const vector<Extrapolation_t>&extraps,
477 if (extraps.size()<2)
return false;
479 for (
unsigned int j=1;j<extraps.size();j++){
480 if (extraps[j].position.Perp()>R){
498 DVector2 x2(prevpos.X(),prevpos.Y());
500 double A = dx.Mod2();
501 double B = 2.0*(x1.X()*dx.X() + x1.Y()*dx.Y());
502 double C = x1.Mod2() - R*R;
504 double sqrt_D=
sqrt(B*B-4.0*A*C);
505 double one_over_denom=0.5/A;
506 double alpha1 = (-B + sqrt_D)*one_over_denom;
507 double alpha2 = (-B - sqrt_D)*one_over_denom;
508 double alpha = alpha1;
509 if(alpha1<0.0 || alpha1>1.0)alpha=alpha2;
510 if (isfinite(alpha)){
514 double ds=(1.-
alpha)*delta.Mag();
515 double v=(extrap.
s-prev.
s)/(extrap.
t-prev.
t);
528 const vector<Extrapolation_t>&extraps,
532 return ExtrapolateToRadius(R,extraps,pos,mom,s,t);
538 const vector<Extrapolation_t>&extrapolations,
540 DVector3 *position_along_wire)
const{
541 if (extrapolations.size()<3)
return 1000.;
544 double z0w=wire->
origin.z();
545 double ux=wire->
udir.x();
546 double uy=wire->
udir.y();
547 double uz=wire->
udir.z();
549 double doca_old=1000.,doca=1000.;
550 for (
unsigned int i=1;i<extrapolations.size();i++){
551 DVector3 trackpos=extrapolations[i].position;
552 double z=trackpos.z();
555 doca=(wirepos-trackpos).Perp();
557 DVector3 trackdir=extrapolations[i-1].momentum;
560 double lambda=M_PI_2-trackdir.Theta();
561 double sinl=
sin(lambda);
562 double cosl=cos(lambda);
563 double phi=trackdir.Phi();
564 double sinphi=
sin(phi);
565 double cosphi=cos(phi);
566 double my_ux=ux*sinl-cosl*cosphi;
567 double my_uy=uy*sinl-cosl*sinphi;
568 double denom=my_ux*my_ux+my_uy*my_uy;
569 double ds=((trackpos.X()-wire->
origin.X()-ux*dzw)*my_ux
570 +(trackpos.Y()-wire->
origin.Y()-uy*dzw)*my_uy)/denom;
571 if (fabs(ds)<2.*fabs(extrapolations[i].s-extrapolations[i-1].s)){
572 trackpos+=ds*trackdir;
573 wirepos=wire->
origin+((trackpos.z()-z0w)/uz)*wire->
udir;
575 doca=(wirepos-trackpos).Perp()*cosstereo;
577 if (pos!=NULL) *pos=trackpos;
578 if (mom!=NULL) *mom=extrapolations[i-1].momentum;
579 if (position_along_wire!=NULL) *position_along_wire=wirepos;
595 #ifdef PROFILE_TRK_TIMES
596 void DTrackFitter::GetProfilingTimes(map<string, prof_time::time_diffs> &my_prof_times)
const
598 my_prof_times = prof_times;
jerror_t CorrectForELoss(const DKinematicData &starting_params, DReferenceTrajectory *rt, DVector3 &pos, DVector3 &mom, double mass)
double CalcDensityEffect(double p, double mass, double density, double Z_over_A, double I)
bool ExtrapolateToRadius(double R, const vector< Extrapolation_t > &extraps, DVector3 &pos, DVector3 &mom, double &t, double &s) const
void TimeDiffNow(std::map< std::string, time_diffs > &prof_times, std::string what)
void FastSwim(const DVector3 &pos, const DVector3 &mom, double q, double smax=2000.0, double zmin=-100., double zmax=1000.0)
void AddHits(vector< const DCDCTrackHit * > cdchits)
DMagneticFieldMap * GetBfield(unsigned int run_number=1)
void SetMass(double mass)
const DVector3 & position(void) const
virtual fit_status_t FitTrack(void)=0
class DFDCPseudo: definition for a reconstructed point in the FDC
The DTrackHitSelector class is a base class for algorithms that will select hits from the drift chamb...
void SetDGeometry(const DGeometry *geom)
fit_status_t FindHitsAndFitTrack(const DKinematicData &starting_params, const DReferenceTrajectory *rt, JEventLoop *loop, double mass=-1.0, int N=0, double t0=QuietNaN, DetectorSystem_t t0_det=SYS_NULL)
mass<0 means get it from starting_params
virtual void GetCDCHits(fit_type_t fit_type, const DReferenceTrajectory *rt, const vector< const DCDCTrackHit * > &cdchits_in, vector< const DCDCTrackHit * > &cdchits_out, int N=20) const =0
DRootGeom * GetRootGeom(unsigned int run_number)
DGeometry * GetDGeometry(unsigned int run_number)
bool FDCSortByZincreasing(DTrackCandidate_factory_FDC::DFDCTrkHit *const &hit1, DTrackCandidate_factory_FDC::DFDCTrkHit *const &hit2)
double charge(void) const
void AddHit(const DCDCTrackHit *cdchit)
const DVector3 & momentum(void) const
This utility class is to help in doing detailed profiling of code. It is uses the unix itimer system ...
void SetPLossDirection(direction_t direction)
virtual void GetFDCHits(fit_type_t fit_type, const DReferenceTrajectory *rt, const vector< const DFDCPseudo * > &fdchits_in, vector< const DFDCPseudo * > &fdchits_out, int N=20) const =0
static Particle_t IDTrack(float locCharge, float locMass)
bool CDCSortByRincreasing(const DCDCTrackHit *const &hit1, const DCDCTrackHit *const &hit2)
double DistToWire(const DCoordinateSystem *wire, const vector< Extrapolation_t > &extrapolations, DVector3 *pos=NULL, DVector3 *mom=NULL, DVector3 *position_along_wire=NULL) const
void GetAllHits(fit_type_t fit_type, const DReferenceTrajectory *rt, const vector< const DCDCTrackHit * > &cdchits_in, const vector< const DFDCPseudo * > &fdchits_in, DTrackFitter *fitter, int N=20) const