18 #include <JANA/JApplication.h>
19 #include <JANA/JEventLoop.h>
20 #include <JANA/JCalibration.h>
36 #define ADJACENT_MATCH_RADIUS 1.0
37 #define MATCH_RADIUS 2.0
65 pthread_mutex_init(&mutex, NULL);
80 cout <<
"initializing" <<endl;
87 alignments.resize(24);
89 for (
unsigned int i=0;i<24;i++){
90 alignments[i].A(kDx)=0.;
91 alignments[i].A(kDy)=0.;
92 alignments[i].A(kDPhi)=0.;
93 alignments[i].E(kDx,kDx)=1.0;
94 alignments[i].E(kDy,kDy)=1.0;
95 alignments[i].E(kDPhi,kDPhi)=1.0;
99 TDirectory *
dir =
new TDirectoryFile(
"FDC",
"FDC");
103 fdctree =
new TTree(
"fdc",
"FDC algnments");
104 fdcbranch = fdctree->Branch(
"T",
"FDC_branch",&fdc_ptr);
122 double endplate_dz,endplate_rmin,endplate_rmax;
123 dgeom->
GetCDCEndplate(endplate_z,endplate_dz,endplate_rmin,endplate_rmax);
124 endplate_z+=0.5*endplate_dz;
126 japp->RootWriteLock();
128 Hqratio_vs_wire= (TH2F *)gROOT->FindObject(
"Hqratio_vs_wire");
129 if (!Hqratio_vs_wire)
130 Hqratio_vs_wire=
new TH2F(
"Hqratio_vs_wire",
"Charge ratio vs wire number",
131 2304,0.5,2304.5,100,-0.5,0.5);
133 Hdelta_z_vs_wire= (TH2F*)gROOT->FindObject(
"Hdelta_z_vs_wire");
134 if (!Hdelta_z_vs_wire)
135 Hdelta_z_vs_wire=
new TH2F(
"Hdelta_z_vs_wire",
"wire z offset vs wire number",
136 2304,0.5,2304.5,100,-0.1,0.1);
138 Hprob = (TH1F*)gROOT->FindObject(
"Hprob");
140 Hprob=
new TH1F(
"Hprob",
"Confidence level for wire-based fit",
143 Hreduced_chi2 = (TH1F*)gROOT->FindObject(
"Hreduced_chi2");
145 Hreduced_chi2=
new TH1F(
"Hreduced_chi2",
"chi^2/ndf",
147 Htime_prob = (TH1F*)gROOT->FindObject(
"Htime_prob");
149 Htime_prob=
new TH1F(
"Htime_prob",
"Confidence level for time-based fit",
152 Hures_vs_layer=(TH2F*)gROOT->FindObject(
"Hures_vs_layer");
153 if (!Hures_vs_layer){
154 Hures_vs_layer=
new TH2F(
"Hures_vs_layer",
"wire-based residuals",
155 24,0.5,24.5,200,-0.5,0.5);
157 Hvres_vs_layer=(TH2F*)gROOT->FindObject(
"Hvres_vs_layer");
158 if (!Hvres_vs_layer){
159 Hvres_vs_layer=
new TH2F(
"Hvres_vs_layer",
"residual for position along wire",
160 24,0.5,24.5,200,-0.5,0.5);
163 Hcand_ty_vs_tx=(TH2F*)gROOT->FindObject(
"Hcand_ty_vs_tx");
164 if (!Hcand_ty_vs_tx){
165 Hcand_ty_vs_tx=
new TH2F(
"Hcand_ty_vs_tx",
"candidate ty vs tx",100,-1,1,
168 Hty_vs_tx=(TH2F*)gROOT->FindObject(
"Hty_vs_tx");
170 Hty_vs_tx=
new TH2F(
"Hty_vs_tx",
"wire-based ty vs tx",100,-1,1,
173 Htime_ty_vs_tx=(TH2F*)gROOT->FindObject(
"Htime_ty_vs_tx");
174 if (!Htime_ty_vs_tx){
175 Htime_ty_vs_tx=
new TH2F(
"Htime_ty_vs_tx",
"time-based ty vs tx",100,-1,1,
178 Htime_y_vs_x=(TH3F*)gROOT->FindObject(
"Htime_y_vs_x");
180 Htime_y_vs_x=
new TH3F(
"Htime_y_vs_x",
"time-based y vs x",24,0.5,24.5,100,-48.5,48.5,
184 Hdrift_time=(TH2F*)gROOT->FindObject(
"Hdrift_time");
186 Hdrift_time=
new TH2F(
"Hdrift_time",
187 "doca vs drift time",201,-21,381,100,0,1);
189 Hdrift_integral=(TH1F*)gROOT->FindObject(
"Hdrift_integral");
190 if (!Hdrift_integral){
191 Hdrift_integral=
new TH1F(
"Hdrift_integral",
192 "drift time integral",140,-20,260);
194 Hz_target=(TH1F*)gROOT->FindObject(
"Hz_target");
196 Hz_target=
new TH1F(
"Hz_target",
197 "Target position",260,0,130);
200 Hres_vs_drift_time=(TH2F*)gROOT->FindObject(
"Hres_vs_drift_time");
201 if (!Hres_vs_drift_time){
202 Hres_vs_drift_time=
new TH2F(
"Hres_vs_drift_time",
"Residual vs drift time",320,-20,300,1000,-1,1);
204 Hdv_vs_dE=(TH2F*)gROOT->FindObject(
"Hdv_vs_dE");
206 Hdv_vs_dE=
new TH2F(
"Hdv_vs_dE",
"dv vs energy dep",100,0,20
e-6,200,-1,1);
209 Hxcand_prob = (TH1F *)gROOT->FindObject(
"Hscand_prob");
211 Hxcand_prob=
new TH1F(
"Hxcand_prob",
"x candidate prob",100,0,1);
213 Hycand_prob = (TH1F *)gROOT->FindObject(
"Hycand_prob");
215 Hycand_prob=
new TH1F(
"Hycand_prob",
"y candidate prob",100,0,1);
218 Hfcal_match=(TH1F*)gROOT->FindObject(
"Hfcal_match");
220 Hfcal_match=
new TH1F(
"Hfcal_match",
221 "match dr for FCAL",200,0,10);
224 Hbcal_match=(TH1F*)gROOT->FindObject(
"Hbcal_match");
226 Hbcal_match=
new TH1F(
"Hbcal_match",
227 "match dr for BCAL",200,0,10);
231 Htheta=(TH1F*)gROOT->FindObject(
"Htheta");
233 Htheta=
new TH1F(
"Htheta",
236 HdEdx=(TH1F*)gROOT->FindObject(
"HdEdx");
238 HdEdx=
new TH1F(
"HdEdx",
"dEdx",100,0,10
e-6);
243 JCalibration *jcalib = dapp->GetJCalibration(0);
244 vector< map<string, float> > tvals;
245 if (jcalib->Get(
"FDC/fdc_drift_Bzero", tvals)==
false){
249 japp->RootFillLock(
this);
251 for(
unsigned int i=0; i<tvals.size(); i++){
252 map<string, float> &row = tvals[i];
253 map<string,float>::iterator iter = row.begin();
254 fdc_drift_table[i] = iter->second;
255 Hdrift_integral->Fill(2.*i-20,fdc_drift_table[i]/0.5);
258 japp->RootFillUnLock(
this);
262 jerr <<
" FDC time-to-distance table not available... bailing..." << endl;
297 vector<const DFCALShower*>fcalshowers;
298 loop->Get(fcalshowers);
299 vector<const DBCALShower*>bcalshowers;
301 vector<const DFDCPseudo*>pseudos;
304 for (
unsigned int i=0;i<pseudos.size();i++){
305 vector<const DFDCCathodeCluster*>cathode_clusters;
306 pseudos[i]->GetT(cathode_clusters);
308 unsigned int wire_number
309 =96*(pseudos[i]->wire->layer-1)+pseudos[i]->wire->wire;
311 for (
unsigned int j=0;j<cathode_clusters[0]->members.size();j++){
312 double q=cathode_clusters[0]->members[j]->q;
315 for (
unsigned int j=0;j<cathode_clusters[1]->members.size();j++){
316 double q=cathode_clusters[1]->members[j]->q;
324 japp->RootFillLock(
this);
326 Hqratio_vs_wire->Fill(wire_number,ratio-1.);
327 Hdelta_z_vs_wire->Fill(wire_number,0.5336*(1.-ratio)/(1.+ratio));
328 japp->RootFillUnLock(
this);
331 if (pseudos.size()>2 && (bcalshowers.size()>0 || fcalshowers.size()>0)){
333 vector<const DFDCPseudo*>packages[4];
334 for (
unsigned int i=0;i<pseudos.size();i++){
335 packages[(pseudos[i]->wire->layer-1)/6].push_back(pseudos[i]);
339 vector<segment_t>segments[4];
340 for (
unsigned int i=0;i<4;i++){
341 FindSegments(packages[i],segments[i]);
345 vector<vector<const DFDCPseudo *> >LinkedSegments;
346 LinkSegments(segments,LinkedSegments);
349 for (
unsigned int k=0;k<LinkedSegments.size();k++){
350 vector<const DFDCPseudo *>hits=LinkedSegments[k];
353 double var_x,var_tx,cov_x_tx,var_y,var_ty,cov_y_ty,chi2x,chi2y;
354 int myndf=hits.size()-2;
355 DMatrix4x1 S=FitLine(hits,var_x,cov_x_tx,var_tx,chi2x,var_y,cov_y_ty,
357 double probx=TMath::Prob(chi2x,myndf);
358 double proby=TMath::Prob(chi2y,myndf);
362 japp->RootFillLock(
this);
364 Hxcand_prob->Fill(probx);
365 Hycand_prob->Fill(proby);
367 japp->RootFillUnLock(
this);
370 bool got_match=
false;
372 double dz=0.,outer_time=0.;
374 for (
unsigned int i=0;i<fcalshowers.size();i++){
375 double fcal_z=fcalshowers[i]->getPosition().z();
376 double x=
S(state_x)+fcal_z*
S(state_tx);
377 double y=
S(state_y)+fcal_z*
S(state_ty);
378 double dx=fcalshowers[i]->getPosition().x()-
x;
379 double dy=fcalshowers[i]->getPosition().y()-
y;
380 double dr=
sqrt(dx*dx+dy*dy);
384 dz=fcal_z-endplate_z;
385 outer_time=fcalshowers[i]->getTime();
392 japp->RootFillLock(
this);
394 Hfcal_match->Fill(drmin);
396 japp->RootFillUnLock(
this);
405 for (
unsigned int i=0;i<bcalshowers.size();i++){
406 double bcal_z=bcalshowers[i]->z;
407 double R2=bcalshowers[i]->x*bcalshowers[i]->x
408 +bcalshowers[i]->y*bcalshowers[i]->y;
409 double x0=
S(state_x);
410 double y0=
S(state_y);
411 double tx=
S(state_tx);
412 double ty=
S(state_ty);
413 double myz=(-x0*tx-y0*ty
414 +
sqrt(R2*(tx*tx+ty*ty)-(x0*ty-y0*tx)*(x0*ty-y0*tx)))
418 double dx=bcalshowers[i]->x-
x;
419 double dy=bcalshowers[i]->y-
y;
420 double dr=
sqrt(dx*dx+dy*dy);
424 dz=bcal_z-endplate_z;
425 outer_time=bcalshowers[i]->t;
431 japp->RootFillLock(
this);
433 Hbcal_match->Fill(drmin);
435 japp->RootFillUnLock(
this);
444 japp->RootFillLock(
this);
446 Hcand_ty_vs_tx->Fill(
S(state_tx),
S(state_ty));
449 double tanl=1./
sqrt(
S(state_tx)*
S(state_tx)+
S(state_ty)*
S(state_ty));
450 double theta=90-180/M_PI*atan(tanl);
451 double sinl=
sin(atan(tanl));
457 vector<double>dElist;
458 for (
unsigned int i=0;i<hits.size();i++){
459 dElist.push_back(hits[i]->dE);
461 sort(dElist.begin(),dElist.end(),
dEdxSort);
462 unsigned int N=dElist.size()/2;
463 for (
unsigned int i=0;i<N;i++){
471 double one_over_tx2=1./(
S(state_tx)*
S(state_tx));
473 =one_over_tx2*(var_x+
S(state_x)*
S(state_x)*one_over_tx2*var_tx
474 -2.*
S(state_x)/
S(state_tx)*cov_x_tx);
475 double one_over_ty2=1./(
S(state_ty)*
S(state_ty));
477 =one_over_ty2*(var_y+
S(state_y)*
S(state_y)*one_over_ty2*var_ty
478 -2.*
S(state_x)/
S(state_ty)*cov_y_ty);
481 =(-
S(state_x)/(
S(state_tx)*varz_x)-
S(state_y)/(
S(state_ty)*varz_y))/
482 (1./varz_x+1/varz_y);
484 Hz_target->Fill(z_target);
486 japp->RootFillUnLock(
this);
490 mT0=outer_time-dz/(29.98*sinl);
504 vector<const DFDCPseudo *>&hits){
505 unsigned int num_hits=hits.size();
508 vector<update_t>updates(num_hits);
509 vector<update_t>smoothed_updates(num_hits);
512 double anneal_factor=1.;
514 anneal_factor=pow(1e6,(
double(NEVENTS-myevt))/(NEVENTS-1.));
515 if (myevt>NEVENTS) anneal_factor=1.;
525 deque<trajectory_t>trajectory;
527 S(state_x)+=endplate_z*
S(state_tx);
528 S(state_y)+=endplate_z*
S(state_ty);
529 SetReferenceTrajectory(endplate_z,S,trajectory,hits);
533 C0(state_x,state_x)=
C0(state_y,state_y)=1.;
534 C0(state_tx,state_tx)=
C0(state_ty,state_ty)=0.01;
537 double chi2=1e16,chi2_old=1e16;
538 unsigned int ndof=0,ndof_old=0;
545 if (KalmanFilter(anneal_factor,S,C,hits,trajectory,updates,chi2,ndof)
549 if (chi2>chi2_old || fabs(chi2_old-chi2)<0.1 || iter==
ITER_MAX)
break;
556 Smooth(S,C,trajectory,hits,updates,smoothed_updates);
561 japp->RootFillLock(
this);
563 double reduced_chi2=chi2_old/double(ndof_old);
564 double prob=TMath::Prob(chi2_old,ndof_old);
566 Hreduced_chi2->Fill(reduced_chi2);
569 Hty_vs_tx->Fill(Sbest(state_tx),Sbest(state_ty));
571 for (
unsigned int i=0;i<smoothed_updates.size();i++){
572 unsigned int layer=hits[i]->wire->layer;
573 Hures_vs_layer->Fill(layer,smoothed_updates[i].ures);
574 Hvres_vs_layer->Fill(layer,smoothed_updates[i].vres);
576 Hres_vs_drift_time->Fill(smoothed_updates[i].drift_time,
577 smoothed_updates[i].ures);
578 Hdrift_time->Fill(smoothed_updates[i].drift_time,
579 smoothed_updates[i].doca);
581 Hdv_vs_dE->Fill(hits[i]->dE,smoothed_updates[i].vres);
586 japp->RootFillUnLock(
this);
589 FindOffsets(hits,smoothed_updates);
592 japp->RootWriteLock();
596 double dxr=alignments[
layer].A(kDx);
597 double dyr=alignments[
layer].A(kDy);
598 fdc.dPhi=alignments[
layer].A(kDPhi);
599 double cosdphi=cos(fdc.dPhi);
600 double sindphi=
sin(fdc.dPhi);
601 double dx=dxr*cosdphi-dyr*sindphi;
602 double dy=dxr*sindphi+dyr*cosdphi;
629 return VALUE_OUT_OF_RANGE;
636 vector<vector<const DFDCPseudo *> >&LinkedSegments){
637 vector<const DFDCPseudo *>myhits;
638 for (
unsigned int i=0;i<4;i++){
639 for (
unsigned int j=0;j<segments[i].size();j++){
640 if (segments[i][j].matched==
false){
641 myhits.assign(segments[i][j].hits.begin(),segments[i][j].hits.end());
643 unsigned i_plus_1=i+1;
645 double tx=segments[i][j].S(state_tx);
646 double ty=segments[i][j].S(state_ty);
647 double x0=segments[i][j].S(state_x);
648 double y0=segments[i][j].S(state_y);
650 for (
unsigned int k=0;k<segments[i_plus_1].size();k++){
651 if (segments[i_plus_1][k].matched==
false){
652 double z=segments[i_plus_1][k].hits[0]->wire->origin.z();
655 if ((proj-segments[i_plus_1][k].hits[0]->xy).Mod()<
MATCH_RADIUS){
656 segments[i_plus_1][k].matched=
true;
657 myhits.insert(myhits.end(),segments[i_plus_1][k].hits.begin(),
658 segments[i_plus_1][k].hits.end());
660 unsigned int i_plus_2=i_plus_1+1;
662 tx=segments[i_plus_1][k].S(state_tx);
663 ty=segments[i_plus_1][k].S(state_ty);
664 x0=segments[i_plus_1][k].S(state_x);
665 y0=segments[i_plus_1][k].S(state_y);
667 for (
unsigned int m=0;m<segments[i_plus_2].size();m++){
668 if (segments[i_plus_2][m].matched==
false){
669 z=segments[i_plus_2][m].hits[0]->wire->origin.z();
670 proj.Set(x0+tx*z,y0+ty*z);
672 if ((proj-segments[i_plus_2][m].hits[0]->xy).Mod()<
MATCH_RADIUS){
673 segments[i_plus_2][m].matched=
true;
674 myhits.insert(myhits.end(),segments[i_plus_2][m].hits.begin(),
675 segments[i_plus_2][m].hits.end());
677 unsigned int i_plus_3=i_plus_2+1;
679 tx=segments[i_plus_2][m].S(state_tx);
680 ty=segments[i_plus_2][m].S(state_ty);
681 x0=segments[i_plus_2][m].S(state_x);
682 y0=segments[i_plus_2][m].S(state_y);
684 for (
unsigned int n=0;n<segments[i_plus_3].size();n++){
685 if (segments[i_plus_3][n].matched==
false){
686 z=segments[i_plus_3][n].hits[0]->wire->origin.z();
687 proj.Set(x0+tx*z,y0+ty*z);
689 if ((proj-segments[i_plus_3][n].hits[0]->xy).Mod()<
MATCH_RADIUS){
690 segments[i_plus_3][n].matched=
true;
691 myhits.insert(myhits.end(),segments[i_plus_3][n].hits.begin(),
692 segments[i_plus_3][n].hits.end());
709 LinkedSegments.push_back(myhits);
720 vector<segment_t>&segments){
721 if (points.size()==0)
return RESOURCE_UNAVAILABLE;
722 vector<int>used(points.size());
726 double old_z=points[0]->wire->origin.z();
727 vector<unsigned int>x_list;
729 for (
unsigned int i=0;i<points.size();i++){
730 used.push_back(
false);
731 if (points[i]->wire->origin.z()!=old_z){
734 old_z=points[i]->wire->origin.z();
736 x_list.push_back(points.size());
738 unsigned int start=0;
740 while (start<x_list.size()-1){
742 for (
unsigned int i=x_list[start];i<x_list[start+1];i++){
750 vector<const DFDCPseudo*>neighbors;
751 neighbors.push_back(points[i]);
752 unsigned int match=0;
753 double delta,delta_min=1000.;
754 for (
unsigned int k=0;k<x_list.size()-1;k++){
757 for (
unsigned int m=x_list[k];m<x_list[k+1];m++){
758 delta=(XY-points[m]->xy).Mod();
765 && used[match]==
false
767 XY=points[match]->xy;
769 neighbors.push_back(points[match]);
772 unsigned int num_neighbors=neighbors.size();
776 for (
unsigned int k=0;k<points.size();k++){
778 for (
unsigned int j=0;j<num_neighbors;j++){
779 delta=(points[k]->xy-neighbors[j]->xy).Mod();
782 abs(neighbors[j]->wire->wire-points[k]->wire->wire)<=1
783 && neighbors[j]->wire->origin.z()==points[k]->wire->origin.z()){
785 neighbors.push_back(points[k]);
792 if (neighbors.size()>4){
795 mysegment.
S=FitLine(neighbors);
796 mysegment.
hits=neighbors;
797 segments.push_back(mysegment);
803 while (start<x_list.size()-1){
804 if (used[x_list[start]]==
false)
break;
816 double var_x,var_tx,cov_x_tx,chi2x;
817 double var_y,var_ty,cov_y_ty,chi2y;
819 return FitLine(fdchits,var_x,cov_x_tx,var_tx,chi2x,var_y,cov_y_ty,var_ty,
828 double &var_x,
double &cov_x_tx,
829 double &var_tx,
double &chi2x,
830 double &var_y,
double &cov_y_ty,
831 double &var_ty,
double &chi2y){
845 for (
unsigned int i=0;i<fdchits.size();i++){
846 double cosa=fdchits[i]->wire->udir.y();
847 double sina=fdchits[i]->wire->udir.x();
848 double x=fdchits[i]->xy.X();
849 double y=fdchits[i]->xy.Y();
850 double z=fdchits[i]->wire->origin.z();
851 double sig2x=cosa*cosa/12+sina*sina*sig2v;
852 double sig2y=sina*sina/12+cosa*cosa*sig2v;
853 double one_over_var1=1/sig2y;
854 double one_over_var2=1/sig2x;
857 S1z+=z*one_over_var1;
858 S1y+=y*one_over_var1;
859 S1zz+=z*z*one_over_var1;
860 S1zy+=z*y*one_over_var1;
863 S2z+=z*one_over_var2;
864 S2x+=x*one_over_var2;
865 S2zz+=z*z*one_over_var2;
866 S2zx+=z*x*one_over_var2;
868 double D1=S1*S1zz-S1z*S1z;
869 double y_intercept=(S1zz*S1y-S1z*S1zy)/D1;
870 double y_slope=(S1*S1zy-S1z*S1y)/D1;
871 double D2=S2*S2zz-S2z*S2z;
872 double x_intercept=(S2zz*S2x-S2z*S2zx)/D2;
873 double x_slope=(S2*S2zx-S2z*S2x)/D2;
888 for (
unsigned int i=0;i<fdchits.size();i++){
889 double cosa=fdchits[i]->wire->udir.y();
890 double sina=fdchits[i]->wire->udir.x();
891 double sig2x=cosa*cosa/12+sina*sina*sig2v;
892 double sig2y=sina*sina/12+cosa*cosa*sig2v;
893 double one_over_var1=1/sig2y;
894 double one_over_var2=1/sig2x;
896 double z=fdchits[i]->wire->origin.z();
897 double dx=fdchits[i]->xy.X()-(x_intercept+x_slope*z);
898 double dy=fdchits[i]->xy.Y()-(y_intercept+y_slope*z);
900 chi2x+=dx*dx*one_over_var2;
901 chi2y+=dy*dy*one_over_var1;
903 return DMatrix4x1(x_intercept,y_intercept,x_slope,y_slope);
910 deque<trajectory_t>&trajectory,
911 vector<const DFDCPseudo *>&hits,
912 vector<strip_update_t>updates,
913 vector<strip_update_t>&smoothed_updates
919 unsigned int max=trajectory.size()-1;
920 S=(trajectory[
max].Skk);
921 C=(trajectory[
max].Ckk);
922 JT=(trajectory[
max].J.Transpose());
925 for (
unsigned int m=max-1;m>0;m--){
926 if (trajectory[m].h_id==0){
927 A=trajectory[m].Ckk*JT*C.
Invert();
928 Ss=trajectory[m].Skk+A*(Ss-
S);
929 Cs=trajectory[m].Ckk+A*(Cs-C)*A.
Transpose();
931 else if (trajectory[m].h_id>0){
932 unsigned int id=trajectory[m].h_id-1;
933 A=updates[id].C*JT*C.
Invert();
935 Ss=updates[id].S+A*(Ss-
S);
944 double cosa=hits[id]->wire->udir.y();
945 double sina=hits[id]->wire->udir.x();
948 double x=Ss(state_x);
949 double y=Ss(state_y);
950 double tx=Ss(state_tx);
951 double ty=Ss(state_ty);
954 unsigned int layer=hits[id]->wire->layer-1;
959 double sindphi=
sin(A(kDPhi));
960 double cosdphi=cos(A(kDPhi));
963 double cospsi=cosa*cosdphi+sina*sindphi;
964 double sinpsi=sina*cosdphi-cosa*sindphi;
971 double upred=x*cospsi-y*sinpsi-dx*cosa+dy*sina;
972 double vpred=x*sinpsi+y*cospsi-dx*sina-dy*cosa;
973 double tu=tx*cospsi-ty*sinpsi;
974 double tv=tx*sinpsi-ty*cospsi;
978 double alpha=atan(tu);
979 double cosalpha=cos(alpha);
980 double sinalpha=
sin(alpha);
983 double uwire=hits[id]->w;
984 double v=hits[id]->s;
985 double d=(upred-uwire)*cosalpha;
986 smoothed_updates[id].vres=v-vpred+tv*d*sinalpha;
987 smoothed_updates[id].ures=(d>0?1.:-1.)*updates[id].drift-d;
989 smoothed_updates[id].id=trajectory[m].h_id;
990 smoothed_updates[id].drift=updates[id].drift;
991 smoothed_updates[id].drift_time=updates[id].drift_time;
992 smoothed_updates[id].S=Ss;
993 smoothed_updates[id].C=Cs;
994 smoothed_updates[id].R=updates[id].R-updates[id].H*dC*updates[id].H_T;
998 JT=trajectory[m].J.Transpose();
1008 A=trajectory[0].Ckk*JT*C.
Invert();
1009 Ss=trajectory[0].Skk+A*(Ss-
S);
1010 Cs=trajectory[0].Ckk+A*(Cs-C)*A.
Transpose();
1017 deque<trajectory_t>&trajectory,
1018 vector<const DFDCPseudo *>&hits,
1019 vector<update_t>updates,
1020 vector<update_t>&smoothed_updates
1026 unsigned int max=trajectory.size()-1;
1027 S=(trajectory[
max].Skk);
1028 C=(trajectory[
max].Ckk);
1029 JT=(trajectory[
max].J.Transpose());
1032 for (
unsigned int m=max-1;m>0;m--){
1033 if (trajectory[m].h_id==0){
1034 A=trajectory[m].Ckk*JT*C.
Invert();
1035 Ss=trajectory[m].Skk+A*(Ss-
S);
1036 Cs=trajectory[m].Ckk+A*(Cs-C)*A.
Transpose();
1038 else if (trajectory[m].h_id>0){
1039 unsigned int id=trajectory[m].h_id-1;
1040 A=updates[id].C*JT*C.
Invert();
1042 Ss=updates[id].S+A*(Ss-
S);
1043 Cs=updates[id].C+dC;
1051 double cosa=hits[id]->wire->udir.y();
1052 double sina=hits[id]->wire->udir.x();
1055 double x=Ss(state_x);
1056 double y=Ss(state_y);
1057 double tx=Ss(state_tx);
1058 double ty=Ss(state_ty);
1061 unsigned int layer=hits[id]->wire->layer-1;
1066 double sindphi=
sin(A(kDPhi));
1067 double cosdphi=cos(A(kDPhi));
1070 double cospsi=cosa*cosdphi+sina*sindphi;
1071 double sinpsi=sina*cosdphi-cosa*sindphi;
1078 double upred=x*cospsi-y*sinpsi-dx*cosa+dy*sina;
1079 double vpred=x*sinpsi+y*cospsi-dx*sina-dy*cosa;
1080 double tu=tx*cospsi-ty*sinpsi;
1081 double tv=tx*sinpsi-ty*cospsi;
1085 double alpha=atan(tu);
1086 double cosalpha=cos(alpha);
1087 double sinalpha=
sin(alpha);
1090 double uwire=hits[id]->w;
1091 double v=hits[id]->s;
1092 double d=(upred-uwire)*cosalpha;
1093 smoothed_updates[id].vres=v-vpred+tv*d*sinalpha;
1094 smoothed_updates[id].ures=(d>0?1.:-1.)*updates[id].drift-d;
1095 smoothed_updates[id].doca=fabs(d);
1097 smoothed_updates[id].id=trajectory[m].h_id;
1098 smoothed_updates[id].drift=updates[id].drift;
1099 smoothed_updates[id].drift_time=updates[id].drift_time;
1100 smoothed_updates[id].S=Ss;
1101 smoothed_updates[id].C=Cs;
1102 smoothed_updates[id].R=updates[id].R-updates[id].H*dC*updates[id].H_T;
1104 S=trajectory[m].Skk;
1105 C=trajectory[m].Ckk;
1106 JT=trajectory[m].J.Transpose();
1116 A=trajectory[0].Ckk*JT*C.
Invert();
1117 Ss=trajectory[0].Skk+A*(Ss-
S);
1118 Cs=trajectory[0].Ckk+A*(Cs-C)*A.
Transpose();
1127 vector<const DFDCPseudo *>&hits,
1128 deque<trajectory_t>&trajectory,
1129 vector<strip_update_t>&updates,
1130 double &chi2,
unsigned int &ndof){
1147 trajectory[0].Skk=
S;
1148 trajectory[0].Ckk=C;
1149 for (
unsigned int k=1;k<trajectory.size();k++){
1150 if (C(0,0)<=0. || C(1,1)<=0. || C(2,2)<=0. || C(3,3)<=0.)
1151 return UNRECOVERABLE_ERROR;
1154 S=trajectory[k].S+J*(S-S0);
1158 trajectory[k].Skk=
S;
1159 trajectory[k].Ckk=C;
1166 if (trajectory[k].h_id>0){
1167 unsigned int id=trajectory[k].h_id-1;
1169 double cosa=hits[id]->wire->udir.y();
1170 double sina=hits[id]->wire->udir.x();
1173 double x=
S(state_x);
1174 double y=
S(state_y);
1175 double tx=
S(state_tx);
1176 double ty=
S(state_ty);
1177 if (isnan(x) || isnan(y))
return UNRECOVERABLE_ERROR;
1180 unsigned int layer=hits[id]->wire->layer-1;
1185 double sindphi=
sin(A(kDPhi));
1186 double cosdphi=cos(A(kDPhi));
1189 double cospsi=cosa*cosdphi+sina*sindphi;
1190 double sinpsi=sina*cosdphi-cosa*sindphi;
1197 double upred=x*cospsi-y*sinpsi-dx*cosa+dy*sina;
1198 double vpred=x*sinpsi+y*cospsi-dx*sina-dy*cosa;
1199 double tu=tx*cospsi-ty*sinpsi;
1200 double tv=tx*sinpsi-ty*cospsi;
1204 double alpha=atan(tu);
1205 double cosalpha=cos(alpha);
1206 double sinalpha=
sin(alpha);
1209 for (
int m=trajectory[k].num_hits-1;m>=0;m--){
1210 unsigned int my_id=
id+m;
1211 double uwire=hits[my_id]->w;
1212 double v=hits[my_id]->s;
1217 double drift_time=hits[my_id]->time-trajectory[k].t-mT0;
1218 double drift=GetDriftDistance(drift_time);
1219 updates[my_id].drift=drift;
1220 updates[my_id].drift_time=drift_time;
1222 double du=upred-uwire;
1223 double d=du*cosalpha;
1225 Mdiff=v-vpred+tv*d*sinalpha;
1231 double sigma=3.7e-8/hits[my_id]->dE+0.0062;
1232 V=anneal_factor*sigma*
sigma;
1235 double sinalpha_cosalpha=sinalpha*cosalpha;
1236 H_T(state_x)=sinpsi-tv*cospsi*sinalpha_cosalpha;
1237 H_T(state_y)=cospsi+tv*sinpsi*sinalpha_cosalpha;
1239 double temp=tv*cosalpha*(cosalpha*cosalpha-sinalpha*sinalpha);
1240 H_T(state_tx)=-d*(sinpsi*sinalpha+cospsi*
temp);
1241 H_T(state_ty)=-d*(cospsi*sinalpha-sinpsi*
temp);
1244 H(state_x)=H_T(state_x);
1245 H(state_y)=H_T(state_y);
1246 H(state_tx)=H_T(state_tx);
1247 H(state_ty)=H_T(state_ty);
1250 updates[my_id].H_T=H_T;
1252 double Vtemp=V+H*C*H_T;
1259 G_T(kDx)=-sina+tv*cosa*sinalpha_cosalpha;
1260 G_T(kDy)=-cosa-tv*sina*sinalpha_cosalpha;
1262 G_T(kDPhi)=-x*cospsi+y*sinpsi
1263 +tu*d*sinalpha-tv*(x*sinpsi+y*cospsi)*sinalpha_cosalpha
1264 -tu*tv*d*cosalpha*(cosalpha*cosalpha-sinalpha*sinalpha)
1270 G(kDPhi)=G_T(kDPhi);
1272 Vtemp=Vtemp+G*E*G_T;
1276 K=(1./Vtemp)*(C*H_T);
1288 updates[my_id].R=scale*V;
1292 chi2+=scale*Mdiff*Mdiff/V;
1299 chi2*=anneal_factor;
1313 vector<const DFDCPseudo *>&hits,
1314 deque<trajectory_t>&trajectory,
1315 vector<update_t>&updates,
1316 double &chi2,
unsigned int &ndof){
1335 trajectory[0].Skk=
S;
1336 trajectory[0].Ckk=C;
1337 for (
unsigned int k=1;k<trajectory.size();k++){
1338 if (C(0,0)<=0. || C(1,1)<=0. || C(2,2)<=0. || C(3,3)<=0.)
1339 return UNRECOVERABLE_ERROR;
1342 S=trajectory[k].S+J*(S-S0);
1346 trajectory[k].Skk=
S;
1347 trajectory[k].Ckk=C;
1354 if (trajectory[k].h_id>0){
1355 unsigned int id=trajectory[k].h_id-1;
1357 double cosa=hits[id]->wire->udir.y();
1358 double sina=hits[id]->wire->udir.x();
1361 double x=
S(state_x);
1362 double y=
S(state_y);
1363 double tx=
S(state_tx);
1364 double ty=
S(state_ty);
1365 if (isnan(x) || isnan(y))
return UNRECOVERABLE_ERROR;
1368 unsigned int layer=hits[id]->wire->layer-1;
1373 double sindphi=
sin(A(kDPhi));
1374 double cosdphi=cos(A(kDPhi));
1377 double cospsi=cosa*cosdphi+sina*sindphi;
1378 double sinpsi=sina*cosdphi-cosa*sindphi;
1385 double upred=x*cospsi-y*sinpsi-dx*cosa+dy*sina;
1386 double vpred=x*sinpsi+y*cospsi-dx*sina-dy*cosa;
1387 double tu=tx*cospsi-ty*sinpsi;
1388 double tv=tx*sinpsi-ty*cospsi;
1392 double alpha=atan(tu);
1393 double cosalpha=cos(alpha);
1394 double sinalpha=
sin(alpha);
1397 for (
int m=trajectory[k].num_hits-1;m>=0;m--){
1398 unsigned int my_id=
id+m;
1399 double uwire=hits[my_id]->w;
1400 double v=hits[my_id]->s;
1407 double drift_time=hits[my_id]->time-trajectory[k].t-mT0;
1408 double drift=GetDriftDistance(drift_time);
1409 updates[my_id].drift=drift;
1410 updates[my_id].drift_time=drift_time;
1412 double du=upred-uwire;
1413 double d=du*cosalpha;
1415 double sign=(du>0)?1.:-1.;
1419 Mdiff(0)=sign*drift-d;
1420 Mdiff(1)=v-vpred+tv*d*sinalpha;
1426 V(0,0)=anneal_factor*GetDriftVariance(drift_time);
1428 double sigma=3.7e-8/hits[my_id]->dE+0.0062;
1429 V(1,1)=anneal_factor*sigma*
sigma;
1432 double sinalpha_cosalpha=sinalpha*cosalpha;
1433 H_T(state_x,0)=cospsi*cosalpha;
1434 H_T(state_y,0)=-sinpsi*cosalpha;
1436 H_T(state_x,1)=sinpsi-tv*cospsi*sinalpha_cosalpha;
1437 H_T(state_y,1)=cospsi+tv*sinpsi*sinalpha_cosalpha;
1439 double temp=d*sinalpha_cosalpha;
1440 H_T(state_tx,0)=-temp*cospsi;
1441 H_T(state_ty,0)=+temp*sinpsi;
1443 temp=tv*cosalpha*(cosalpha*cosalpha-sinalpha*sinalpha);
1444 H_T(state_tx,1)=-d*(sinpsi*sinalpha+cospsi*
temp);
1445 H_T(state_ty,1)=-d*(cospsi*sinalpha-sinpsi*
temp);
1448 H(0,state_x)=H_T(state_x,0);
1449 H(0,state_y)=H_T(state_y,0);
1451 H(1,state_x)=H_T(state_x,1);
1452 H(1,state_y)=H_T(state_y,1);
1454 H(0,state_tx)=H_T(state_tx,0);
1455 H(0,state_ty)=H_T(state_ty,0);
1458 updates[my_id].H_T=H_T;
1467 G_T(kDx,0)=-cosa*cosalpha;
1468 G_T(kDy,0)=+sina*cosalpha;
1470 G_T(kDx,1)=-sina+tv*cosa*sinalpha_cosalpha;
1471 G_T(kDy,1)=-cosa-tv*sina*sinalpha_cosalpha;
1473 G_T(kDPhi,0)=(x*sinpsi+y*cospsi)*cosalpha;
1474 G_T(kDPhi,1)=-x*cospsi+y*sinpsi
1475 +tu*d*sinalpha-tv*(x*sinpsi+y*cospsi)*sinalpha_cosalpha
1476 -tu*tv*d*cosalpha*(cosalpha*cosalpha-sinalpha*sinalpha)
1480 G(0,kDx)=G_T(kDx,0);
1481 G(0,kDy)=G_T(kDy,0);
1482 G(1,kDx)=G_T(kDx,1);
1483 G(1,kDy)=G_T(kDy,1);
1484 G(0,kDPhi)=G_T(kDPhi,0);
1485 G(1,kDPhi)=G_T(kDPhi,1);
1487 Vtemp=Vtemp+G*E*G_T;
1507 updates[my_id].R=RC;
1518 chi2*=anneal_factor;
1527 vector<const DFDCPseudo *>&pseudos){
1529 DMatrix4x4 J(1.,0.,1.,0., 0.,1.,0.,1., 0.,0.,1.,0., 0.,0.,0.,1.);
1542 trajectory.push_front(temp);
1545 unsigned int itrajectory=0;
1546 for (
unsigned int i=0;i<pseudos.size();i++){
1547 zhit=pseudos[i]->wire->origin.z();
1549 if (fabs(zhit-old_zhit)<
EPS){
1552 temp.
S=trajectory[0].S;
1553 temp.
t=trajectory[0].t;
1555 temp.
z=trajectory[0].z;
1558 trajectory.push_front(temp);
1567 temp.
J(state_x,state_tx)=dz;
1568 temp.
J(state_y,state_ty)=dz;
1570 temp.
t=(t+=dz*
sqrt(1+
S(state_tx)*
S(state_tx)+
S(state_ty)*
S(state_ty))
1573 temp.
S(state_x)=
S(state_x)+
S(state_tx)*dz;
1574 temp.
S(state_y)=
S(state_y)+
S(state_ty)*dz;
1575 temp.
S(state_tx)=
S(state_tx);
1576 temp.
S(state_ty)=
S(state_ty);
1588 trajectory.push_front(temp);
1600 temp.
J(state_x,state_tx)=-dz;
1601 temp.
J(state_y,state_ty)=-dz;
1603 temp.
t=(t+=dz*
sqrt(1+
S(state_tx)*
S(state_tx)+
S(state_ty)*
S(state_ty))
1606 temp.
S(state_x)=
S(state_x)+
S(state_tx)*dz;
1607 temp.
S(state_y)=
S(state_y)+
S(state_ty)*dz;
1608 temp.
S(state_tx)=
S(state_tx);
1609 temp.
S(state_ty)=
S(state_ty);
1611 trajectory.push_front(temp);
1615 for (
unsigned int i=0;i<trajectory.size();i++){
1616 printf(
" x %f y %f z %f hit %d\n",trajectory[i].
S(state_x),
1617 trajectory[i].
S(state_y),trajectory[i].z,trajectory[i].h_id);
1628 double tpar[5]={0.0162,-0.001064,3.59e-5,-5.04e-7,2.46e-9};
1629 double sigma=tpar[0]+0.008;
1630 for (
int i=1;i<5;i++) sigma+=tpar[i]*pow(t,i);
1637 #define FDC_T0_OFFSET 20.
1643 double d=fdc_drift_table[id];
1647 double dd=fdc_drift_table[
id+1]-fdc_drift_table[id];
1655 vector<update_t>smoothed_updates){
1659 unsigned int num_hits=hits.size();
1661 for (
unsigned int i=0;i<num_hits;i++){
1662 if (smoothed_updates[i].
id>0){
1663 unsigned int id=smoothed_updates[i].id-1;
1664 double x=smoothed_updates[i].S(state_x);
1665 double y=smoothed_updates[i].S(state_y);
1666 double tx=smoothed_updates[i].S(state_tx);
1667 double ty=smoothed_updates[i].S(state_ty);
1669 double cosa=hits[id]->wire->udir.y();
1670 double sina=hits[id]->wire->udir.x();
1671 double uwire=hits[id]->w;
1672 double v=hits[id]->s;
1675 unsigned int layer=hits[id]->wire->layer-1;
1680 double sindphi=
sin(A(kDPhi));
1681 double cosdphi=cos(A(kDPhi));
1684 double cospsi=cosa*cosdphi+sina*sindphi;
1685 double sinpsi=sina*cosdphi-cosa*sindphi;
1692 double upred=x*cospsi-y*sinpsi-dx*cosa+dy*sina;
1693 double vpred=x*sinpsi+y*cospsi-dx*sina-dy*cosa;
1694 double tu=tx*cospsi-ty*sinpsi;
1695 double tv=tx*sinpsi-ty*cospsi;
1696 double du=upred-uwire;
1697 double sign=(du>0)?1.:-1.;
1701 double alpha=atan(tu);
1702 double cosalpha=cos(alpha);
1703 double sinalpha=
sin(alpha);
1706 G_T(kDx,0)=-cosa*cosalpha;
1707 G_T(kDy,0)=+sina*cosalpha;
1709 double sinalpha_cosalpha=sinalpha*cosalpha;
1710 G_T(kDx,1)=-sina+tv*cosa*sinalpha_cosalpha;
1711 G_T(kDy,1)=-cosa-tv*sina*sinalpha_cosalpha;
1713 double d=du*cosalpha;
1714 G_T(kDPhi,0)=(x*sinpsi+y*cospsi)*cosalpha;
1715 G_T(kDPhi,1)=-x*cospsi+y*sinpsi
1716 +tu*d*sinalpha-tv*(x*sinpsi+y*cospsi)*sinalpha_cosalpha
1717 -tu*tv*d*cosalpha*(cosalpha*cosalpha-sinalpha*sinalpha)
1721 G(0,kDx)=G_T(kDx,0);
1722 G(0,kDy)=G_T(kDy,0);
1723 G(1,kDx)=G_T(kDx,1);
1724 G(1,kDy)=G_T(kDy,1);
1725 G(0,kDPhi)=G_T(kDPhi,0);
1726 G(1,kDPhi)=G_T(kDPhi,1);
1729 DMatrix2x2 InvV=(smoothed_updates[i].R+G*E*G_T).Invert();
1733 Mdiff(0)=sign*smoothed_updates[i].drift-du*cosalpha;
1734 Mdiff(1)=v-vpred+tv*du*cosalpha*sinalpha;
1740 if (Etemp(0,0)>0 && Etemp(1,1)>0 && Etemp(2,2)>0){
1741 alignments[
layer].E=Etemp;
1742 alignments[
layer].A=A+Ka*Mdiff;
1752 vector<strip_update_t>smoothed_updates){
1756 unsigned int num_hits=hits.size();
1758 for (
unsigned int i=0;i<num_hits;i++){
1759 if (smoothed_updates[i].
id>0){
1760 unsigned int id=smoothed_updates[i].id-1;
1761 double x=smoothed_updates[i].S(state_x);
1762 double y=smoothed_updates[i].S(state_y);
1763 double tx=smoothed_updates[i].S(state_tx);
1764 double ty=smoothed_updates[i].S(state_ty);
1766 double cosa=hits[id]->wire->udir.y();
1767 double sina=hits[id]->wire->udir.x();
1768 double uwire=hits[id]->w;
1769 double v=hits[id]->s;
1772 unsigned int layer=hits[id]->wire->layer-1;
1777 double sindphi=
sin(A(kDPhi));
1778 double cosdphi=cos(A(kDPhi));
1781 double cospsi=cosa*cosdphi+sina*sindphi;
1782 double sinpsi=sina*cosdphi-cosa*sindphi;
1789 double upred=x*cospsi-y*sinpsi-dx*cosa+dy*sina;
1790 double vpred=x*sinpsi+y*cospsi-dx*sina-dy*cosa;
1791 double tu=tx*cospsi-ty*sinpsi;
1792 double tv=tx*sinpsi-ty*cospsi;
1793 double du=upred-uwire;
1797 double alpha=atan(tu);
1798 double cosalpha=cos(alpha);
1799 double sinalpha=
sin(alpha);
1802 double sinalpha_cosalpha=sinalpha*cosalpha;
1803 G_T(kDx)=-sina+tv*cosa*sinalpha_cosalpha;
1804 G_T(kDy)=-cosa-tv*sina*sinalpha_cosalpha;
1806 double d=du*cosalpha;
1807 G_T(kDPhi)=-x*cospsi+y*sinpsi
1808 +tu*d*sinalpha-tv*(x*sinpsi+y*cospsi)*sinalpha_cosalpha
1809 -tu*tv*d*cosalpha*(cosalpha*cosalpha-sinalpha*sinalpha)
1815 G(kDPhi)=G_T(kDPhi);
1818 double InvV=1./(smoothed_updates[i].R+G*E*G_T);
1821 double Mdiff=v-vpred+tv*d*sinalpha;
1827 if (Etemp(0,0)>0 && Etemp(1,1)>0 && Etemp(2,2)>0){
1828 alignments[
layer].E=Etemp;
1829 alignments[
layer].A=A+dA;
jerror_t erun(void)
Invoked via DEventProcessor virtual method.
#define ADJACENT_MATCH_RADIUS
static vector< vector< DFDCWire * > > fdcwires
jerror_t KalmanFilter(double anneal_factor, DMatrix4x1 &S, DMatrix4x4 &C, vector< const DFDCPseudo * > &hits, deque< trajectory_t > &trajectory, vector< strip_update_t > &updates, double &chi2, unsigned int &ndof)
DEventProcessor_fdc_hists()
jerror_t evnt(JEventLoop *loop, uint64_t eventnumber)
Invoked via DEventProcessor virtual method.
jerror_t FindSegments(vector< const DFDCPseudo * > &pseudos, vector< segment_t > &segments)
~DEventProcessor_fdc_hists()
double Chi2(const DMatrix2x1 &R) const
vector< const DFDCPseudo * > hits
DGeometry * GetDGeometry(unsigned int run_number)
double GetDriftVariance(double t)
Double_t sigma[NCHANNELS]
jerror_t fini(void)
Invoked via DEventProcessor virtual method.
jerror_t LinkSegments(vector< segment_t >segments[4], vector< vector< const DFDCPseudo * > > &LinkedSegments)
jerror_t Smooth(DMatrix4x1 &Ss, DMatrix4x4 &Cs, deque< trajectory_t > &trajectory, vector< const DFDCPseudo * > &hits, vector< strip_update_t >updates, vector< strip_update_t > &smoothed_updates)
jerror_t brun(JEventLoop *loop, int32_t runnumber)
bool GetFDCWires(vector< vector< DFDCWire * > > &fdcwires) const
jerror_t DoFilter(DMatrix4x1 &S, vector< const DFDCPseudo * > &fdchits)
bool GetCDCEndplate(double &z, double &dz, double &rmin, double &rmax) const
double GetDriftDistance(double t)
printf("string=%s", string)
DMatrix4x1 FitLine(vector< const DFDCPseudo * > &fdchits)
jerror_t SetReferenceTrajectory(double z, DMatrix4x1 &S, deque< trajectory_t > &trajectory, vector< const DFDCPseudo * > &wires)
for(auto locVertexInfo:dStepVertexInfos)
bool dEdxSort(double a, double b)
jerror_t init(void)
Invoked via DEventProcessor virtual method.
jerror_t FindOffsets(vector< const DFDCPseudo * > &hits, vector< update_t >smoothed_updates)