Hall-D Software  alpha
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DEventProcessor_fdc_hists.cc
Go to the documentation of this file.
1 // $Id: DEventProcessor_fdc_hists.cc 2774 2007-07-19 15:59:02Z davidl $
2 //
3 // File: DEventProcessor_fdc_hists.cc
4 // Created: Sun Apr 24 06:45:21 EDT 2005
5 // Creator: davidl (on Darwin Harriet.local 7.8.0 powerpc)
6 //
7 
8 #include <iostream>
9 #include <cmath>
10 using namespace std;
11 
12 #include <TThread.h>
13 #include <TMath.h>
14 #include <TROOT.h>
15 
17 
18 #include <JANA/JApplication.h>
19 #include <JANA/JEventLoop.h>
20 #include <JANA/JCalibration.h>
21 
22 #include <DANA/DApplication.h>
23 #include <TRACKING/DMCThrown.h>
24 #include <TRACKING/DMCTrackHit.h>
25 #include <FDC/DFDCGeometry.h>
26 #include <FDC/DFDCHit.h>
27 #include <FDC/DFDCPseudo.h>
28 #include <FDC/DFDCSegment.h>
29 #include <FDC/DFDCCathodeCluster.h>
30 
31 #include <BCAL/DBCALShower.h>
32 #include <DVector2.h>
33 
34 #define EPS 1e-3
35 #define ITER_MAX 10
36 #define ADJACENT_MATCH_RADIUS 1.0
37 #define MATCH_RADIUS 2.0
38 
39 // The executable should define the ROOTfile global variable. It will
40 // be automatically linked when dlopen is called.
41 //extern TFile *ROOTfile;
42 
43 // Routine used to create our DEventProcessor
44 extern "C"{
45 void InitPlugin(JApplication *app){
46  InitJANAPlugin(app);
47  app->AddProcessor(new DEventProcessor_fdc_hists());
48 }
49 } // "C"
50 
51 
52 bool dEdxSort(double a,double b){
53  return (a<b);
54 }
55 
56 
57 //------------------
58 // DEventProcessor_fdc_hists
59 //------------------
61 {
62  fdc_ptr = &fdc;
63  fdchit_ptr = &fdchit;
64 
65  pthread_mutex_init(&mutex, NULL);
66 }
67 
68 //------------------
69 // ~DEventProcessor_fdc_hists
70 //------------------
72 {
73 }
74 
75 //------------------
76 // init
77 //------------------
79 {
80  cout << "initializing" <<endl;
81 
82  mT0=0.;
83  myevt=0;
84 
85  DoAlign=false;
86  //DoAlign=true;
87  alignments.resize(24);
88  if (DoAlign){
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;
96  }
97  }
98 
99  TDirectory *dir = new TDirectoryFile("FDC","FDC");
100  dir->cd();
101 
102  // Create Tree
103  fdctree = new TTree("fdc","FDC algnments");
104  fdcbranch = fdctree->Branch("T","FDC_branch",&fdc_ptr);
105 
106  // Go back up to the parent directory
107  dir->cd("../");
108 
109  return NOERROR;
110 }
111 
112 //------------------
113 // brun
114 //------------------
115 jerror_t DEventProcessor_fdc_hists::brun(JEventLoop *loop, int32_t runnumber)
116 {
117  DApplication* dapp=dynamic_cast<DApplication*>(loop->GetJApplication());
118  const DGeometry *dgeom = dapp->GetDGeometry(runnumber);
119  dgeom->GetFDCWires(fdcwires);
120 
121  // Get the position of the CDC downstream endplate from DGeometry
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;
125 
126  japp->RootWriteLock(); //ACQUIRE ROOT LOCK
127 
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);
132 
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);
137 
138  Hprob = (TH1F*)gROOT->FindObject("Hprob");
139  if (!Hprob)
140  Hprob=new TH1F("Hprob","Confidence level for wire-based fit",
141  100,0.0,1.);
142 
143  Hreduced_chi2 = (TH1F*)gROOT->FindObject("Hreduced_chi2");
144  if (!Hreduced_chi2)
145  Hreduced_chi2=new TH1F("Hreduced_chi2","chi^2/ndf",
146  1000,0,100);
147  Htime_prob = (TH1F*)gROOT->FindObject("Htime_prob");
148  if (!Htime_prob)
149  Htime_prob=new TH1F("Htime_prob","Confidence level for time-based fit",
150  100,0,1);
151 
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);
156  }
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);
161  }
162 
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,
166  100,-1,1);
167  }
168  Hty_vs_tx=(TH2F*)gROOT->FindObject("Hty_vs_tx");
169  if (!Hty_vs_tx){
170  Hty_vs_tx=new TH2F("Hty_vs_tx","wire-based ty vs tx",100,-1,1,
171  100,-1,1);
172  }
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,
176  100,-1,1);
177  }
178  Htime_y_vs_x=(TH3F*)gROOT->FindObject("Htime_y_vs_x");
179  if (!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,
181  100,-48.5,48.5);
182  }
183 
184  Hdrift_time=(TH2F*)gROOT->FindObject("Hdrift_time");
185  if (!Hdrift_time){
186  Hdrift_time=new TH2F("Hdrift_time",
187  "doca vs drift time",201,-21,381,100,0,1);
188  }
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);
193  }
194  Hz_target=(TH1F*)gROOT->FindObject("Hz_target");
195  if (!Hz_target){
196  Hz_target=new TH1F("Hz_target",
197  "Target position",260,0,130);
198  }
199 
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);
203  }
204  Hdv_vs_dE=(TH2F*)gROOT->FindObject("Hdv_vs_dE");
205  if (!Hdv_vs_dE){
206  Hdv_vs_dE=new TH2F("Hdv_vs_dE","dv vs energy dep",100,0,20e-6,200,-1,1);
207  }
208 
209  Hxcand_prob = (TH1F *)gROOT->FindObject("Hscand_prob");
210  if (!Hxcand_prob){
211  Hxcand_prob=new TH1F("Hxcand_prob","x candidate prob",100,0,1);
212  }
213  Hycand_prob = (TH1F *)gROOT->FindObject("Hycand_prob");
214  if (!Hycand_prob){
215  Hycand_prob=new TH1F("Hycand_prob","y candidate prob",100,0,1);
216  }
217 
218  Hfcal_match=(TH1F*)gROOT->FindObject("Hfcal_match");
219  if (!Hfcal_match){
220  Hfcal_match=new TH1F("Hfcal_match",
221  "match dr for FCAL",200,0,10);
222  }
223 
224  Hbcal_match=(TH1F*)gROOT->FindObject("Hbcal_match");
225  if (!Hbcal_match){
226  Hbcal_match=new TH1F("Hbcal_match",
227  "match dr for BCAL",200,0,10);
228  }
229 
230 
231  Htheta=(TH1F*)gROOT->FindObject("Htheta");
232  if (!Htheta){
233  Htheta=new TH1F("Htheta",
234  "#theta",100,0,25);
235  }
236  HdEdx=(TH1F*)gROOT->FindObject("HdEdx");
237  if (!HdEdx){
238  HdEdx=new TH1F("HdEdx","dEdx",100,0,10e-6);
239  }
240 
241  japp->RootUnLock(); //RELEASE ROOT LOCK
242 
243  JCalibration *jcalib = dapp->GetJCalibration(0); // need run number here
244  vector< map<string, float> > tvals;
245  if (jcalib->Get("FDC/fdc_drift_Bzero", tvals)==false){
246 
247  // FILL HISTOGRAMS
248  // Since we are filling histograms local to this plugin, it will not interfere with other ROOT operations: can use plugin-wide ROOT fill lock
249  japp->RootFillLock(this); //ACQUIRE ROOT FILL LOCK
250 
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);
256  }
257 
258  japp->RootFillUnLock(this); //RELEASE ROOT FILL LOCK
259 
260  }
261  else{
262  jerr << " FDC time-to-distance table not available... bailing..." << endl;
263  exit(0);
264  }
265 
266 
267 
268  return NOERROR;
269 }
270 
271 //------------------
272 // erun
273 //------------------
275 {
276 
277 
278  return NOERROR;
279 }
280 
281 //------------------
282 // fini
283 //------------------
285 {
286 
287  return NOERROR;
288 }
289 
290 //------------------
291 // evnt
292 //------------------
293 jerror_t DEventProcessor_fdc_hists::evnt(JEventLoop *loop, uint64_t eventnumber){
294  myevt++;
295 
296  // Get BCAL showers, FCAL showers and FDC space points
297  vector<const DFCALShower*>fcalshowers;
298  loop->Get(fcalshowers);
299  vector<const DBCALShower*>bcalshowers;
300  //loop->Get(bcalshowers);
301  vector<const DFDCPseudo*>pseudos;
302  loop->Get(pseudos);
303 
304  for (unsigned int i=0;i<pseudos.size();i++){
305  vector<const DFDCCathodeCluster*>cathode_clusters;
306  pseudos[i]->GetT(cathode_clusters);
307 
308  unsigned int wire_number
309  =96*(pseudos[i]->wire->layer-1)+pseudos[i]->wire->wire;
310  double q1=0,q2=0;
311  for (unsigned int j=0;j<cathode_clusters[0]->members.size();j++){
312  double q=cathode_clusters[0]->members[j]->q;
313  if (q>q1) q1=q;
314  }
315  for (unsigned int j=0;j<cathode_clusters[1]->members.size();j++){
316  double q=cathode_clusters[1]->members[j]->q;
317  if (q>q2) q2=q;
318  }
319 
320  double ratio=q2/q1;
321 
322  // FILL HISTOGRAMS
323  // Since we are filling histograms local to this plugin, it will not interfere with other ROOT operations: can use plugin-wide ROOT fill lock
324  japp->RootFillLock(this); //ACQUIRE ROOT FILL LOCK
325 
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); //RELEASE ROOT FILL LOCK
329  }
330 
331  if (pseudos.size()>2 && (bcalshowers.size()>0 || fcalshowers.size()>0)){
332  // Group FDC hits by package
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]);
336  }
337 
338  // Link hits in each package together into track segments
339  vector<segment_t>segments[4];
340  for (unsigned int i=0;i<4;i++){
341  FindSegments(packages[i],segments[i]);
342  }
343 
344  // Link the segments together to form track candidadates
345  vector<vector<const DFDCPseudo *> >LinkedSegments;
346  LinkSegments(segments,LinkedSegments);
347 
348  // Loop over linked segments
349  for (unsigned int k=0;k<LinkedSegments.size();k++){
350  vector<const DFDCPseudo *>hits=LinkedSegments[k];
351 
352  // Perform preliminary line fits for current set of linked segments
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,
356  var_ty,chi2y);
357  double probx=TMath::Prob(chi2x,myndf);
358  double proby=TMath::Prob(chi2y,myndf);
359 
360  // FILL HISTOGRAMS
361  // Since we are filling histograms local to this plugin, it will not interfere with other ROOT operations: can use plugin-wide ROOT fill lock
362  japp->RootFillLock(this); //ACQUIRE ROOT FILL LOCK
363 
364  Hxcand_prob->Fill(probx);
365  Hycand_prob->Fill(proby);
366 
367  japp->RootFillUnLock(this); //RELEASE ROOT FILL LOCK
368 
369  // Match to outer detectors
370  bool got_match=false;
371  // First match to FCAL
372  double dz=0.,outer_time=0.;
373  double drmin=1000.;
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);
381 
382  if (dr<drmin){
383  drmin=dr;
384  dz=fcal_z-endplate_z;
385  outer_time=fcalshowers[i]->getTime();
386  }
387 
388  }
389 
390  // FILL HISTOGRAMS
391  // Since we are filling histograms local to this plugin, it will not interfere with other ROOT operations: can use plugin-wide ROOT fill lock
392  japp->RootFillLock(this); //ACQUIRE ROOT FILL LOCK
393 
394  Hfcal_match->Fill(drmin);
395 
396  japp->RootFillUnLock(this); //RELEASE ROOT FILL LOCK
397 
398  if (drmin<4.){
399  got_match=true;
400  // outer_time-=2.218; //empirical correction
401  }
402  else{
403  // Match to BCAL
404  drmin=1000.;
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)))
415  /(tx*tx+ty*ty);
416  double x=x0+myz*tx;
417  double y=y0+myz*ty;
418  double dx=bcalshowers[i]->x-x;
419  double dy=bcalshowers[i]->y-y;
420  double dr=sqrt(dx*dx+dy*dy);
421 
422  if (dr<drmin){
423  drmin=dr;
424  dz=bcal_z-endplate_z;
425  outer_time=bcalshowers[i]->t;
426  }
427  }
428 
429  // FILL HISTOGRAMS
430  // Since we are filling histograms local to this plugin, it will not interfere with other ROOT operations: can use plugin-wide ROOT fill lock
431  japp->RootFillLock(this); //ACQUIRE ROOT FILL LOCK
432 
433  Hbcal_match->Fill(drmin);
434 
435  japp->RootFillUnLock(this); //RELEASE ROOT FILL LOCK
436 
437  if (drmin<4.)
438  got_match=true;
439  }
440  if (got_match){
441 
442  // FILL HISTOGRAMS
443  // Since we are filling histograms local to this plugin, it will not interfere with other ROOT operations: can use plugin-wide ROOT fill lock
444  japp->RootFillLock(this); //ACQUIRE ROOT FILL LOCK
445 
446  Hcand_ty_vs_tx->Fill(S(state_tx),S(state_ty));
447 
448  //compute tangent of dip angle and related angular quantities
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));
452 
453  Htheta->Fill(theta);
454 
455  // Compute dEdx
456  double dEdx=0.;
457  vector<double>dElist;
458  for (unsigned int i=0;i<hits.size();i++){
459  dElist.push_back(hits[i]->dE);
460  }
461  sort(dElist.begin(),dElist.end(),dEdxSort);
462  unsigned int N=dElist.size()/2;
463  for (unsigned int i=0;i<N;i++){
464  dEdx+=dElist[i];
465  }
466  dEdx/=N*sinl;
467 
468  HdEdx->Fill(dEdx);
469 
470  // Estimate z-position of "vertex" in target
471  double one_over_tx2=1./(S(state_tx)*S(state_tx));
472  double varz_x
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));
476  double varz_y
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);
479  // Weighted average from x and y line fits, ignoring x-y correlations
480  double z_target
481  =(-S(state_x)/(S(state_tx)*varz_x)-S(state_y)/(S(state_ty)*varz_y))/
482  (1./varz_x+1/varz_y);
483 
484  Hz_target->Fill(z_target);
485 
486  japp->RootFillUnLock(this); //RELEASE ROOT FILL LOCK
487 
488  // Estimate for t0 at the beginning of track assuming particle is
489  // moving at the speed of light
490  mT0=outer_time-dz/(29.98*sinl);
491 
492  // Run the Kalman Filter algorithm
493  DoFilter(S,hits);
494  }
495  }
496  }
497 
498  return NOERROR;
499 }
500 
501 // Steering routine for the kalman filter
502 jerror_t
504  vector<const DFDCPseudo *>&hits){
505  unsigned int num_hits=hits.size();
506  // vector<strip_update_t>strip_updates(num_hits);
507  //vector<strip_update_t>smoothed_strip_updates(num_hits);
508  vector<update_t>updates(num_hits);
509  vector<update_t>smoothed_updates(num_hits);
510 
511  int NEVENTS=200000;
512  double anneal_factor=1.;
513  if (DoAlign){
514  anneal_factor=pow(1e6,(double(NEVENTS-myevt))/(NEVENTS-1.));
515  if (myevt>NEVENTS) anneal_factor=1.;
516  }
517  //anneal_factor=1.;
518  //anneal_factor=10000.;
519 
520  // Best guess for state vector at "vertex"
521  DMatrix4x1 Sbest;
522 
523  // Use the result from the initial line fit to form a reference trajectory
524  // for the track.
525  deque<trajectory_t>trajectory;
526  // double start_z=hits[0]->wire->origin.z()-1.;
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);
530 
531  // Intial guess for covariance matrix
532  DMatrix4x4 C,C0,Cbest;
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;
535 
536  // Chi-squared and degrees of freedom
537  double chi2=1e16,chi2_old=1e16;
538  unsigned int ndof=0,ndof_old=0;
539  unsigned iter=0;
540  for(;;){
541  iter++;
542  chi2_old=chi2;
543  ndof_old=ndof;
544  C=C0;
545  if (KalmanFilter(anneal_factor,S,C,hits,trajectory,updates,chi2,ndof)
546  !=NOERROR) break;
547 
548  //printf("=======chi2 %f\n",chi2);
549  if (chi2>chi2_old || fabs(chi2_old-chi2)<0.1 || iter==ITER_MAX) break;
550 
551  // Save the current state and covariance matrixes
552  Cbest=C;
553  Sbest=S;
554 
555  // run the smoother (opposite direction to filter)
556  Smooth(S,C,trajectory,hits,updates,smoothed_updates);
557  }
558 
559  if (iter>1){
560 
561  japp->RootFillLock(this); //ACQUIRE ROOT FILL LOCK
562 
563  double reduced_chi2=chi2_old/double(ndof_old);
564  double prob=TMath::Prob(chi2_old,ndof_old);
565  Hprob->Fill(prob);
566  Hreduced_chi2->Fill(reduced_chi2);
567 
568  if (prob>0.05){
569  Hty_vs_tx->Fill(Sbest(state_tx),Sbest(state_ty));
570 
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);
575  if (layer==1){
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);
580 
581  Hdv_vs_dE->Fill(hits[i]->dE,smoothed_updates[i].vres);
582  }
583  }
584  }
585 
586  japp->RootFillUnLock(this); //RELEASE ROOT FILL LOCK
587 
588  if (DoAlign){
589  FindOffsets(hits,smoothed_updates);
590 
591  // Although we are only filling objects local to this plugin, TTree::Fill() periodically writes to file: Global ROOT lock
592  japp->RootWriteLock(); //ACQUIRE ROOT LOCK
593 
594  for (unsigned int layer=0;layer<24;layer++){
595  // Set up to fill tree
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;
603  fdc.dX=dx;
604  fdc.dY=dy;
605  fdc.layer=layer;
606  fdc.N=myevt;
607 
608 
609  fdctree->Fill();
610 
611  }
612 
613  japp->RootUnLock(); //RELEASE ROOT LOCK
614 
615  }
616 
617 
618 
619  /*
620  printf("-------Event %d\n",myevt);
621  for (unsigned int i=0;i<24;i++)
622  printf("A %f %f %f\n",alignments[i].A(kDx),alignments[i].A(kDy),
623  alignments[i].A(kDPhi));
624  */
625 
626  return NOERROR;
627  }
628 
629  return VALUE_OUT_OF_RANGE;
630 }
631 
632 
633 // Link segments from package to package by doing straight-line projections
634 jerror_t
635 DEventProcessor_fdc_hists::LinkSegments(vector<segment_t>segments[4],
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());
642 
643  unsigned i_plus_1=i+1;
644  if (i_plus_1<4){
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);
649 
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();
653  DVector2 proj(x0+tx*z,y0+ty*z);
654 
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());
659 
660  unsigned int i_plus_2=i_plus_1+1;
661  if (i_plus_2<4){
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);
666 
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);
671 
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());
676 
677  unsigned int i_plus_3=i_plus_2+1;
678  if (i_plus_3<4){
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);
683 
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);
688 
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());
693 
694  break;
695  } // matched a segment
696  }
697  } // loop over last set of segments
698  } // if we have another package to loop over
699  break;
700  } // matched a segment
701  }
702  } // loop over second-to-last set of segments
703  }
704  break;
705  } // matched a segment
706  }
707  } // loop over third-to-last set of segments
708  }
709  LinkedSegments.push_back(myhits);
710  myhits.clear();
711  } // check if we have already used this segment
712  } // loop over first set of segments
713  } // loop over packages
714 
715  return NOERROR;
716 }
717 
718 // Find segments by associating adjacent hits within a package together.
719 jerror_t DEventProcessor_fdc_hists::FindSegments(vector<const DFDCPseudo*>&points,
720  vector<segment_t>&segments){
721  if (points.size()==0) return RESOURCE_UNAVAILABLE;
722  vector<int>used(points.size());
723 
724  // Put indices for the first point in each plane before the most downstream
725  // plane in the vector x_list.
726  double old_z=points[0]->wire->origin.z();
727  vector<unsigned int>x_list;
728  x_list.push_back(0);
729  for (unsigned int i=0;i<points.size();i++){
730  used.push_back(false);
731  if (points[i]->wire->origin.z()!=old_z){
732  x_list.push_back(i);
733  }
734  old_z=points[i]->wire->origin.z();
735  }
736  x_list.push_back(points.size());
737 
738  unsigned int start=0;
739  // loop over the start indices, starting with the first plane
740  while (start<x_list.size()-1){
741  // Now loop over the list of track segment start points
742  for (unsigned int i=x_list[start];i<x_list[start+1];i++){
743  if (used[i]==false){
744  used[i]=true;
745 
746  // Point in the current plane in the package
747  DVector2 XY=points[i]->xy;
748 
749  // Create list of nearest neighbors
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++){
755  delta_min=1000.;
756  match=0;
757  for (unsigned int m=x_list[k];m<x_list[k+1];m++){
758  delta=(XY-points[m]->xy).Mod();
759  if (delta<delta_min && delta<MATCH_RADIUS){
760  delta_min=delta;
761  match=m;
762  }
763  }
764  if (match!=0
765  && used[match]==false
766  ){
767  XY=points[match]->xy;
768  used[match]=true;
769  neighbors.push_back(points[match]);
770  }
771  }
772  unsigned int num_neighbors=neighbors.size();
773 
774  bool do_sort=false;
775  // Look for hits adjacent to the ones we have in our segment candidate
776  for (unsigned int k=0;k<points.size();k++){
777  if (!used[k]){
778  for (unsigned int j=0;j<num_neighbors;j++){
779  delta=(points[k]->xy-neighbors[j]->xy).Mod();
780 
781  if (delta<ADJACENT_MATCH_RADIUS &&
782  abs(neighbors[j]->wire->wire-points[k]->wire->wire)<=1
783  && neighbors[j]->wire->origin.z()==points[k]->wire->origin.z()){
784  used[k]=true;
785  neighbors.push_back(points[k]);
786  do_sort=true;
787  }
788  }
789  }
790  } // loop looking for hits adjacent to hits on segment
791 
792  if (neighbors.size()>4){
793  segment_t mysegment;
794  mysegment.matched=false;
795  mysegment.S=FitLine(neighbors);
796  mysegment.hits=neighbors;
797  segments.push_back(mysegment);
798  }
799  }
800  }// loop over start points within a plane
801 
802  // Look for a new plane to start looking for a segment
803  while (start<x_list.size()-1){
804  if (used[x_list[start]]==false) break;
805  start++;
806  }
807 
808  }
809 
810  return NOERROR;
811 }
812 
813 
814 DMatrix4x1
815 DEventProcessor_fdc_hists::FitLine(vector<const DFDCPseudo*> &fdchits){
816  double var_x,var_tx,cov_x_tx,chi2x;
817  double var_y,var_ty,cov_y_ty,chi2y;
818 
819  return FitLine(fdchits,var_x,cov_x_tx,var_tx,chi2x,var_y,cov_y_ty,var_ty,
820  chi2y);
821 }
822 
823 
824 // Use linear regression on the hits to obtain a first guess for the state
825 // vector. Method taken from Numerical Recipes in C.
826 DMatrix4x1
827 DEventProcessor_fdc_hists::FitLine(vector<const DFDCPseudo*> &fdchits,
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){
832  double S1=0.;
833  double S1z=0.;
834  double S1y=0.;
835  double S1zz=0.;
836  double S1zy=0.;
837  double S2=0.;
838  double S2z=0.;
839  double S2x=0.;
840  double S2zz=0.;
841  double S2zx=0.;
842 
843  double sig2v=0.01; // rough guess;
844 
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;
855 
856  S1+=one_over_var1;
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;
861 
862  S2+=one_over_var2;
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;
867  }
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;
874 
875  // Covariance matrix for x
876  var_x=S2zz/D2;
877  var_tx=S2/D2;
878  cov_x_tx=-S2z/D2;
879 
880  // Covariance matrix for y
881  var_y=S1zz/D1;
882  var_ty=S1/D1;
883  cov_y_ty=-S1z/D1;
884 
885  // Compute chi2 for the line fits, ignoring correlations between x and y
886  chi2x=0;
887  chi2y=0;
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;
895 
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);
899 
900  chi2x+=dx*dx*one_over_var2;
901  chi2y+=dy*dy*one_over_var1;
902  }
903  return DMatrix4x1(x_intercept,y_intercept,x_slope,y_slope);
904 
905 }
906 
907 
908 // Kalman smoother
910  deque<trajectory_t>&trajectory,
911  vector<const DFDCPseudo *>&hits,
912  vector<strip_update_t>updates,
913  vector<strip_update_t>&smoothed_updates
914  ){
915  DMatrix4x1 S;
916  DMatrix4x4 C,dC;
917  DMatrix4x4 JT,A;
918 
919  unsigned int max=trajectory.size()-1;
920  S=(trajectory[max].Skk);
921  C=(trajectory[max].Ckk);
922  JT=(trajectory[max].J.Transpose());
923  //Ss=S;
924  //Cs=C;
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();
930  }
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();
934  dC=A*(Cs-C)*A.Transpose();
935  Ss=updates[id].S+A*(Ss-S);
936  Cs=updates[id].C+dC;
937  /*
938  printf("-------\n");
939  updates[id].C;
940  Cs.Print();
941  */
942 
943  // Nominal rotation of wire planes
944  double cosa=hits[id]->wire->udir.y();
945  double sina=hits[id]->wire->udir.x();
946 
947  // State vector
948  double x=Ss(state_x);
949  double y=Ss(state_y);
950  double tx=Ss(state_tx);
951  double ty=Ss(state_ty);
952 
953  // Get the aligment vector and error matrix for this layer
954  unsigned int layer=hits[id]->wire->layer-1;
955  DMatrix3x3 E=alignments[layer].E;
956  DMatrix3x1 A=alignments[layer].A;
957  double dx=A(kDx);
958  double dy=A(kDy);
959  double sindphi=sin(A(kDPhi));
960  double cosdphi=cos(A(kDPhi));
961 
962  // Components of rotation matrix for converting global to local coords.
963  double cospsi=cosa*cosdphi+sina*sindphi;
964  double sinpsi=sina*cosdphi-cosa*sindphi;
965 
966  // x,y and tx,ty in local coordinate system
967  // To transform from (x,y) to (u,v), need to do a rotation:
968  // u = x*cosa-y*sina
969  // v = y*cosa+x*sina
970  // (without alignment offsets)
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;
975 
976  // Variables for angle of incidence with respect to the z-direction in
977  // the u-z plane
978  double alpha=atan(tu);
979  double cosalpha=cos(alpha);
980  double sinalpha=sin(alpha);
981 
982  // Smoothed residuals
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;
988 
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;
995  }
996  S=trajectory[m].Skk;
997  C=trajectory[m].Ckk;
998  JT=trajectory[m].J.Transpose();
999  }
1000  /*
1001  printf("-----end\n");
1002  Ss.Print();
1003  Cs.Print();
1004  printf("--ckk \n");
1005  C.Print();
1006  */
1007 
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();
1011 
1012  return NOERROR;
1013 }
1014 
1015 // Kalman smoother
1017  deque<trajectory_t>&trajectory,
1018  vector<const DFDCPseudo *>&hits,
1019  vector<update_t>updates,
1020  vector<update_t>&smoothed_updates
1021  ){
1022  DMatrix4x1 S;
1023  DMatrix4x4 C,dC;
1024  DMatrix4x4 JT,A;
1025 
1026  unsigned int max=trajectory.size()-1;
1027  S=(trajectory[max].Skk);
1028  C=(trajectory[max].Ckk);
1029  JT=(trajectory[max].J.Transpose());
1030  //Ss=S;
1031  //Cs=C;
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();
1037  }
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();
1041  dC=A*(Cs-C)*A.Transpose();
1042  Ss=updates[id].S+A*(Ss-S);
1043  Cs=updates[id].C+dC;
1044  /*
1045  printf("-------\n");
1046  updates[id].C;
1047  Cs.Print();
1048  */
1049 
1050  // Nominal rotation of wire planes
1051  double cosa=hits[id]->wire->udir.y();
1052  double sina=hits[id]->wire->udir.x();
1053 
1054  // State vector
1055  double x=Ss(state_x);
1056  double y=Ss(state_y);
1057  double tx=Ss(state_tx);
1058  double ty=Ss(state_ty);
1059 
1060  // Get the aligment vector and error matrix for this layer
1061  unsigned int layer=hits[id]->wire->layer-1;
1062  DMatrix3x3 E=alignments[layer].E;
1063  DMatrix3x1 A=alignments[layer].A;
1064  double dx=A(kDx);
1065  double dy=A(kDy);
1066  double sindphi=sin(A(kDPhi));
1067  double cosdphi=cos(A(kDPhi));
1068 
1069  // Components of rotation matrix for converting global to local coords.
1070  double cospsi=cosa*cosdphi+sina*sindphi;
1071  double sinpsi=sina*cosdphi-cosa*sindphi;
1072 
1073  // x,y and tx,ty in local coordinate system
1074  // To transform from (x,y) to (u,v), need to do a rotation:
1075  // u = x*cosa-y*sina
1076  // v = y*cosa+x*sina
1077  // (without alignment offsets)
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;
1082 
1083  // Variables for angle of incidence with respect to the z-direction in
1084  // the u-z plane
1085  double alpha=atan(tu);
1086  double cosalpha=cos(alpha);
1087  double sinalpha=sin(alpha);
1088 
1089  // Smoothed residuals
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);
1096 
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;
1103  }
1104  S=trajectory[m].Skk;
1105  C=trajectory[m].Ckk;
1106  JT=trajectory[m].J.Transpose();
1107  }
1108  /*
1109  printf("-----end\n");
1110  Ss.Print();
1111  Cs.Print();
1112  printf("--ckk \n");
1113  C.Print();
1114  */
1115 
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();
1119 
1120  return NOERROR;
1121 }
1122 
1123 // Perform Kalman Filter for the current trajectory
1124 jerror_t
1126  DMatrix4x1 &S,DMatrix4x4 &C,
1127  vector<const DFDCPseudo *>&hits,
1128  deque<trajectory_t>&trajectory,
1129  vector<strip_update_t>&updates,
1130  double &chi2,unsigned int &ndof){
1131  DMatrix1x4 H; // Track projection matrix
1132  DMatrix4x1 H_T; // Transpose of track projection matrix
1133  DMatrix4x1 K; // Kalman gain matrix
1134  double V=0; // Measurement covariance
1135  double Mdiff=0.;
1136  DMatrix4x4 I; // identity matrix
1137  DMatrix4x4 J; // Jacobian matrix
1138  DMatrix4x1 S0; // State vector from reference trajectory
1139 
1140  //Initialize chi2 and ndof
1141  chi2=0.;
1142  ndof=0;
1143 
1144  // Loop over all steps in the trajectory
1145  S0=trajectory[0].S;
1146  J=trajectory[0].J;
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;
1152 
1153  // Propagate the state and covariance matrix forward in z
1154  S=trajectory[k].S+J*(S-S0);
1155  C=J*C*J.Transpose();
1156 
1157  // Save the current state and covariance matrix
1158  trajectory[k].Skk=S;
1159  trajectory[k].Ckk=C;
1160 
1161  // Save S and J for the next step
1162  S0=trajectory[k].S;
1163  J=trajectory[k].J;
1164 
1165  // Correct S and C for the hit
1166  if (trajectory[k].h_id>0){
1167  unsigned int id=trajectory[k].h_id-1;
1168 
1169  double cosa=hits[id]->wire->udir.y();
1170  double sina=hits[id]->wire->udir.x();
1171 
1172  // State vector
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;
1178 
1179  // Get the aligment vector and error matrix for this layer
1180  unsigned int layer=hits[id]->wire->layer-1;
1181  DMatrix3x3 E=alignments[layer].E;
1182  DMatrix3x1 A=alignments[layer].A;
1183  double dx=A(kDx);
1184  double dy=A(kDy);
1185  double sindphi=sin(A(kDPhi));
1186  double cosdphi=cos(A(kDPhi));
1187 
1188  // Components of rotation matrix for converting global to local coords.
1189  double cospsi=cosa*cosdphi+sina*sindphi;
1190  double sinpsi=sina*cosdphi-cosa*sindphi;
1191 
1192  // x,y and tx,ty in local coordinate system
1193  // To transform from (x,y) to (u,v), need to do a rotation:
1194  // u = x*cosa-y*sina
1195  // v = y*cosa+x*sina
1196  // (without alignment offsets)
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;
1201 
1202  // Variables for angle of incidence with respect to the z-direction in
1203  // the u-z plane
1204  double alpha=atan(tu);
1205  double cosalpha=cos(alpha);
1206  double sinalpha=sin(alpha);
1207 
1208  // Difference between measurement and projection
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;
1213 
1214  //printf("cospsi %f sinpsi %f\n",cospsi,sinpsi);
1215  //printf("u %f %f v %f %f\n",uwire,upred,v,vpred);
1216 
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;
1221 
1222  double du=upred-uwire;
1223  double d=du*cosalpha;
1224 
1225  Mdiff=v-vpred+tv*d*sinalpha;
1226 
1227  //printf("tdrift %f d %f %f\n",drift_time,drift,d);
1228  //printf("dv %f ddv %f\n",v-vpred,tv*d*sinalpha);
1229 
1230  // Variance of measurement error
1231  double sigma=3.7e-8/hits[my_id]->dE+0.0062;
1232  V=anneal_factor*sigma*sigma;
1233 
1234  // Matrix for transforming from state-vector space to measurement space
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;
1238 
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);
1242 
1243  // H-matrix transpose
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);
1248 
1249  updates[my_id].H=H;
1250  updates[my_id].H_T=H_T;
1251 
1252  double Vtemp=V+H*C*H_T;
1253 
1254  if (DoAlign){
1255  // Matrices to rotate alignment error matrix into measurement space
1256  DMatrix1x3 G;
1257  DMatrix3x1 G_T;
1258 
1259  G_T(kDx)=-sina+tv*cosa*sinalpha_cosalpha;
1260  G_T(kDy)=-cosa-tv*sina*sinalpha_cosalpha;
1261 
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)
1265  ;
1266 
1267  // G-matrix transpose
1268  G(kDx)=G_T(kDx);
1269  G(kDy)=G_T(kDy);
1270  G(kDPhi)=G_T(kDPhi);
1271 
1272  Vtemp=Vtemp+G*E*G_T;
1273  }
1274 
1275  // Compute Kalman gain matrix
1276  K=(1./Vtemp)*(C*H_T);
1277 
1278  // Update the state vector
1279  S+=Mdiff*K;
1280  updates[my_id].S=S;
1281 
1282  // Update state vector covariance matrix
1283  C=C-K*(H*C);
1284  updates[my_id].C=C;
1285 
1286  // Update chi2 for this trajectory
1287  double scale=1-H*K;
1288  updates[my_id].R=scale*V;
1289 
1290  //printf("chi2 %f for %d\n",RC.Chi2(R),my_id);
1291 
1292  chi2+=scale*Mdiff*Mdiff/V;
1293  ndof++;
1294  }
1295 
1296  }
1297 
1298  }
1299  chi2*=anneal_factor;
1300 
1301  ndof-=4;
1302 
1303  return NOERROR;
1304 }
1305 
1306 
1307 
1308 
1309 // Perform Kalman Filter for the current trajectory
1310 jerror_t
1312  DMatrix4x1 &S,DMatrix4x4 &C,
1313  vector<const DFDCPseudo *>&hits,
1314  deque<trajectory_t>&trajectory,
1315  vector<update_t>&updates,
1316  double &chi2,unsigned int &ndof){
1317  DMatrix2x4 H; // Track projection matrix
1318  DMatrix4x2 H_T; // Transpose of track projection matrix
1319  DMatrix4x2 K; // Kalman gain matrix
1320  DMatrix2x2 V(0.020833,0,0,0); // Measurement covariance
1321  DMatrix2x2 Vtemp;
1322  DMatrix2x1 Mdiff;
1323  DMatrix2x2 InvV; // Inverse of error matrix
1324  DMatrix4x4 I; // identity matrix
1325  DMatrix4x4 J; // Jacobian matrix
1326  DMatrix4x1 S0; // State vector from reference trajectory
1327 
1328  //Initialize chi2 and ndof
1329  chi2=0.;
1330  ndof=0;
1331 
1332  // Loop over all steps in the trajectory
1333  S0=trajectory[0].S;
1334  J=trajectory[0].J;
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;
1340 
1341  // Propagate the state and covariance matrix forward in z
1342  S=trajectory[k].S+J*(S-S0);
1343  C=J*C*J.Transpose();
1344 
1345  // Save the current state and covariance matrix
1346  trajectory[k].Skk=S;
1347  trajectory[k].Ckk=C;
1348 
1349  // Save S and J for the next step
1350  S0=trajectory[k].S;
1351  J=trajectory[k].J;
1352 
1353  // Correct S and C for the hit
1354  if (trajectory[k].h_id>0){
1355  unsigned int id=trajectory[k].h_id-1;
1356 
1357  double cosa=hits[id]->wire->udir.y();
1358  double sina=hits[id]->wire->udir.x();
1359 
1360  // State vector
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;
1366 
1367  // Get the aligment vector and error matrix for this layer
1368  unsigned int layer=hits[id]->wire->layer-1;
1369  DMatrix3x3 E=alignments[layer].E;
1370  DMatrix3x1 A=alignments[layer].A;
1371  double dx=A(kDx);
1372  double dy=A(kDy);
1373  double sindphi=sin(A(kDPhi));
1374  double cosdphi=cos(A(kDPhi));
1375 
1376  // Components of rotation matrix for converting global to local coords.
1377  double cospsi=cosa*cosdphi+sina*sindphi;
1378  double sinpsi=sina*cosdphi-cosa*sindphi;
1379 
1380  // x,y and tx,ty in local coordinate system
1381  // To transform from (x,y) to (u,v), need to do a rotation:
1382  // u = x*cosa-y*sina
1383  // v = y*cosa+x*sina
1384  // (without alignment offsets)
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;
1389 
1390  // Variables for angle of incidence with respect to the z-direction in
1391  // the u-z plane
1392  double alpha=atan(tu);
1393  double cosalpha=cos(alpha);
1394  double sinalpha=sin(alpha);
1395 
1396  // Difference between measurement and projection
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;
1401 
1402  //printf("cospsi %f sinpsi %f\n",cospsi,sinpsi);
1403  //printf("u %f %f v %f %f\n",uwire,upred,v,vpred);
1404 
1405 
1406  // Find drift distance
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;
1411 
1412  double du=upred-uwire;
1413  double d=du*cosalpha;
1414 
1415  double sign=(du>0)?1.:-1.;
1416 
1417  // Difference between measured and predicted vectors
1418  //drift=0.25;
1419  Mdiff(0)=sign*drift-d;
1420  Mdiff(1)=v-vpred+tv*d*sinalpha;
1421 
1422  //printf("tdrift %f d %f %f\n",drift_time,drift,d);
1423  //printf("dv %f ddv %f\n",v-vpred,tv*d*sinalpha);
1424 
1425  // Variance of measurement error
1426  V(0,0)=anneal_factor*GetDriftVariance(drift_time);
1427 
1428  double sigma=3.7e-8/hits[my_id]->dE+0.0062;
1429  V(1,1)=anneal_factor*sigma*sigma;
1430 
1431  // Matrix for transforming from state-vector space to measurement space
1432  double sinalpha_cosalpha=sinalpha*cosalpha;
1433  H_T(state_x,0)=cospsi*cosalpha;
1434  H_T(state_y,0)=-sinpsi*cosalpha;
1435 
1436  H_T(state_x,1)=sinpsi-tv*cospsi*sinalpha_cosalpha;
1437  H_T(state_y,1)=cospsi+tv*sinpsi*sinalpha_cosalpha;
1438 
1439  double temp=d*sinalpha_cosalpha;
1440  H_T(state_tx,0)=-temp*cospsi;
1441  H_T(state_ty,0)=+temp*sinpsi;
1442 
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);
1446 
1447  // H-matrix transpose
1448  H(0,state_x)=H_T(state_x,0);
1449  H(0,state_y)=H_T(state_y,0);
1450 
1451  H(1,state_x)=H_T(state_x,1);
1452  H(1,state_y)=H_T(state_y,1);
1453 
1454  H(0,state_tx)=H_T(state_tx,0);
1455  H(0,state_ty)=H_T(state_ty,0);
1456 
1457  updates[my_id].H=H;
1458  updates[my_id].H_T=H_T;
1459 
1460  DMatrix2x2 Vtemp=V+H*C*H_T;
1461 
1462  if (DoAlign){
1463  // Matrices to rotate alignment error matrix into measurement space
1464  DMatrix2x3 G;
1465  DMatrix3x2 G_T;
1466 
1467  G_T(kDx,0)=-cosa*cosalpha;
1468  G_T(kDy,0)=+sina*cosalpha;
1469 
1470  G_T(kDx,1)=-sina+tv*cosa*sinalpha_cosalpha;
1471  G_T(kDy,1)=-cosa-tv*sina*sinalpha_cosalpha;
1472 
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)
1477  ;
1478 
1479  // G-matrix transpose
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);
1486 
1487  Vtemp=Vtemp+G*E*G_T;
1488  }
1489 
1490  // Variance for this hit
1491  InvV=Vtemp.Invert();
1492 
1493  // Compute Kalman gain matrix
1494  K=(C*H_T)*InvV;
1495 
1496  // Update the state vector
1497  S+=K*Mdiff;
1498  updates[my_id].S=S;
1499 
1500  // Update state vector covariance matrix
1501  C=C-K*(H*C);
1502  updates[my_id].C=C;
1503 
1504  // Update chi2 for this trajectory
1505  DMatrix2x1 R=Mdiff-H*K*Mdiff;
1506  DMatrix2x2 RC=V-H*K*V;
1507  updates[my_id].R=RC;
1508 
1509  //printf("chi2 %f for %d\n",RC.Chi2(R),my_id);
1510 
1511  chi2+=RC.Chi2(R);
1512  ndof+=2;
1513  }
1514 
1515  }
1516 
1517  }
1518  chi2*=anneal_factor;
1519  ndof-=4;
1520 
1521  return NOERROR;
1522 }
1523 
1524 // Reference trajectory for the track
1526 ::SetReferenceTrajectory(double z,DMatrix4x1 &S,deque<trajectory_t>&trajectory,
1527  vector<const DFDCPseudo *>&pseudos){
1528  // Jacobian matrix
1529  DMatrix4x4 J(1.,0.,1.,0., 0.,1.,0.,1., 0.,0.,1.,0., 0.,0.,0.,1.);
1530 
1531  double dz=1.1;
1532  double t=0.;
1534  temp.S=S;
1535  temp.J=J;
1536  temp.Skk=Zero4x1;
1537  temp.Ckk=Zero4x4;
1538  temp.h_id=0;
1539  temp.num_hits=0;
1540  temp.z=z;
1541  temp.t=0.;
1542  trajectory.push_front(temp);
1543  double zhit=z;
1544  double old_zhit=z;
1545  unsigned int itrajectory=0;
1546  for (unsigned int i=0;i<pseudos.size();i++){
1547  zhit=pseudos[i]->wire->origin.z();
1548 
1549  if (fabs(zhit-old_zhit)<EPS){
1551  temp.J=J;
1552  temp.S=trajectory[0].S;
1553  temp.t=trajectory[0].t;
1554  temp.h_id=i+1;
1555  temp.z=trajectory[0].z;
1556  temp.Skk=Zero4x1;
1557  temp.Ckk=Zero4x4;
1558  trajectory.push_front(temp);
1559 
1560  continue;
1561  }
1562  bool done=false;
1563  while (!done){
1564  double new_z=z+dz;
1566  temp.J=J;
1567  temp.J(state_x,state_tx)=dz;
1568  temp.J(state_y,state_ty)=dz;
1569  // Flight time: assume particle is moving at the speed of light
1570  temp.t=(t+=dz*sqrt(1+S(state_tx)*S(state_tx)+S(state_ty)*S(state_ty))
1571  /29.98);
1572  //propagate the state to the next z position
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);
1577  temp.Skk=Zero4x1;
1578  temp.Ckk=Zero4x4;
1579  temp.h_id=0;
1580  temp.num_hits=0;
1581  if (new_z>zhit){
1582  new_z=zhit;
1583  temp.h_id=i+1;
1584  temp.num_hits=1;
1585  done=true;
1586  }
1587  temp.z=new_z;
1588  trajectory.push_front(temp);
1589  S=temp.S;
1590  itrajectory++;
1591  z=new_z;
1592  }
1593  old_zhit=zhit;
1594  }
1595  temp.Skk=Zero4x1;
1596  temp.Ckk=Zero4x4;
1597  temp.h_id=0;
1598  temp.z=z+dz;
1599  temp.J=J;
1600  temp.J(state_x,state_tx)=-dz;
1601  temp.J(state_y,state_ty)=-dz;
1602  // Flight time: assume particle is moving at the speed of light
1603  temp.t=(t+=dz*sqrt(1+S(state_tx)*S(state_tx)+S(state_ty)*S(state_ty))
1604  /29.98);
1605  //propagate the state to the next z position
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);
1610  S=temp.S;
1611  trajectory.push_front(temp);
1612 
1613  if (false){
1614  printf("Trajectory:\n");
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);
1618  }
1619  }
1620 
1621  return NOERROR;
1622 }
1623 
1624 // Crude approximation for the variance in drift distance due to smearing
1626  if (t<0) t=0;
1627  if (t>115)t=115;
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);
1631 
1632  //double sigma=0.0159-0.000281*t+3.82e-6*t*t;
1633 
1634  return sigma*sigma;
1635 }
1636 
1637 #define FDC_T0_OFFSET 20.
1638 // Interpolate on a table to convert time to distance for the fdc
1640  int id=int((t+FDC_T0_OFFSET)/2.);
1641  if (id<0) id=0;
1642  if (id>138) id=138;
1643  double d=fdc_drift_table[id];
1644 
1645  if (id!=138){
1646  double frac=0.5*(t+FDC_T0_OFFSET-2.*double(id));
1647  double dd=fdc_drift_table[id+1]-fdc_drift_table[id];
1648  d+=frac*dd;
1649  }
1650  return d;
1651 }
1652 
1653 jerror_t
1654 DEventProcessor_fdc_hists::FindOffsets(vector<const DFDCPseudo *>&hits,
1655  vector<update_t>smoothed_updates){
1656  DMatrix2x3 G;//matrix relating alignment vector to measurement coords
1657  DMatrix3x2 G_T; // .. and its transpose
1658 
1659  unsigned int num_hits=hits.size();
1660 
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);
1668 
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;
1673 
1674  // Get the aligment vector and error matrix for this layer
1675  unsigned int layer=hits[id]->wire->layer-1;
1676  DMatrix3x1 A=alignments[layer].A;
1677  DMatrix3x3 E=alignments[layer].E;
1678  double dx=A(kDx);
1679  double dy=A(kDy);
1680  double sindphi=sin(A(kDPhi));
1681  double cosdphi=cos(A(kDPhi));
1682 
1683  // Components of rotation matrix for converting global to local coords.
1684  double cospsi=cosa*cosdphi+sina*sindphi;
1685  double sinpsi=sina*cosdphi-cosa*sindphi;
1686 
1687  // x,y and tx,ty in local coordinate system
1688  // To transform from (x,y) to (u,v), need to do a rotation:
1689  // u = x*cosa-y*sina
1690  // v = y*cosa+x*sina
1691  // (without alignment offsets)
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.;
1698 
1699  // Variables for angle of incidence with respect to the z-direction in
1700  // the u-z plane
1701  double alpha=atan(tu);
1702  double cosalpha=cos(alpha);
1703  double sinalpha=sin(alpha);
1704 
1705  // Transform from alignment vector coords to measurement coords
1706  G_T(kDx,0)=-cosa*cosalpha;
1707  G_T(kDy,0)=+sina*cosalpha;
1708 
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;
1712 
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)
1718  ;
1719 
1720  // G-matrix transpose
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);
1727 
1728  // Inverse of error matrix
1729  DMatrix2x2 InvV=(smoothed_updates[i].R+G*E*G_T).Invert();
1730 
1731  // Difference between measurement and projection
1732  DMatrix2x1 Mdiff;
1733  Mdiff(0)=sign*smoothed_updates[i].drift-du*cosalpha;
1734  Mdiff(1)=v-vpred+tv*du*cosalpha*sinalpha;
1735 
1736  // update the alignment vector and covariance
1737  DMatrix3x2 Ka=(E*G_T)*InvV;
1738  DMatrix3x1 dA=Ka*Mdiff;
1739  DMatrix3x3 Etemp=E-Ka*G*E;
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;
1743  }
1744  }
1745  }
1746 
1747  return NOERROR;
1748 }
1749 
1750 jerror_t
1751 DEventProcessor_fdc_hists::FindOffsets(vector<const DFDCPseudo *>&hits,
1752  vector<strip_update_t>smoothed_updates){
1753  DMatrix1x3 G;//matrix relating alignment vector to measurement coords
1754  DMatrix3x1 G_T; // .. and its transpose
1755 
1756  unsigned int num_hits=hits.size();
1757 
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);
1765 
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;
1770 
1771  // Get the aligment vector and error matrix for this layer
1772  unsigned int layer=hits[id]->wire->layer-1;
1773  DMatrix3x1 A=alignments[layer].A;
1774  DMatrix3x3 E=alignments[layer].E;
1775  double dx=A(kDx);
1776  double dy=A(kDy);
1777  double sindphi=sin(A(kDPhi));
1778  double cosdphi=cos(A(kDPhi));
1779 
1780  // Components of rotation matrix for converting global to local coords.
1781  double cospsi=cosa*cosdphi+sina*sindphi;
1782  double sinpsi=sina*cosdphi-cosa*sindphi;
1783 
1784  // x,y and tx,ty in local coordinate system
1785  // To transform from (x,y) to (u,v), need to do a rotation:
1786  // u = x*cosa-y*sina
1787  // v = y*cosa+x*sina
1788  // (without alignment offsets)
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;
1794 
1795  // Variables for angle of incidence with respect to the z-direction
1796  // in the u-z plane
1797  double alpha=atan(tu);
1798  double cosalpha=cos(alpha);
1799  double sinalpha=sin(alpha);
1800 
1801  // Transform from alignment vector coords to measurement coords
1802  double sinalpha_cosalpha=sinalpha*cosalpha;
1803  G_T(kDx)=-sina+tv*cosa*sinalpha_cosalpha;
1804  G_T(kDy)=-cosa-tv*sina*sinalpha_cosalpha;
1805 
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)
1810  ;
1811 
1812  // G-matrix transpose
1813  G(kDx)=G_T(kDx);
1814  G(kDy)=G_T(kDy);
1815  G(kDPhi)=G_T(kDPhi);
1816 
1817  // Inverse of variance
1818  double InvV=1./(smoothed_updates[i].R+G*E*G_T);
1819 
1820  // Difference between measurement and projection
1821  double Mdiff=v-vpred+tv*d*sinalpha;
1822 
1823  // update the alignment vector and covariance
1824  DMatrix3x1 Ka=InvV*(E*G_T);
1825  DMatrix3x1 dA=Mdiff*Ka;
1826  DMatrix3x3 Etemp=E-Ka*G*E;
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;
1830  }
1831  }
1832  }
1833 
1834  return NOERROR;
1835 }
jerror_t erun(void)
Invoked via DEventProcessor virtual method.
#define MATCH_RADIUS
#define ADJACENT_MATCH_RADIUS
#define FDC_T0_OFFSET
DApplication * dapp
DMatrix2x2 Invert()
Definition: DMatrix2x2.h:64
Int_t layer
TVector2 DVector2
Definition: DVector2.h:9
Double_t x[NCHANNELS]
Definition: st_tw_resols.C:39
#define EPS
#define C0
Definition: src/md5.cpp:24
#define y
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)
jerror_t evnt(JEventLoop *loop, uint64_t eventnumber)
Invoked via DEventProcessor virtual method.
jerror_t FindSegments(vector< const DFDCPseudo * > &pseudos, vector< segment_t > &segments)
DMatrix4x4 Transpose()
Definition: DMatrix4x4.h:174
JApplication * japp
double Chi2(const DMatrix2x1 &R) const
Definition: DMatrix2x2.h:78
vector< const DFDCPseudo * > hits
TEllipse * e
const double alpha
#define H(x, y, z)
TLatex tx
InitPlugin_t InitPlugin
DMatrix4x4 Invert()
Definition: DMatrix4x4.h:161
DGeometry * GetDGeometry(unsigned int run_number)
Double_t sigma[NCHANNELS]
Definition: st_tw_resols.C:37
jerror_t fini(void)
Invoked via DEventProcessor virtual method.
#define ITER_MAX
jerror_t LinkSegments(vector< segment_t >segments[4], vector< vector< const DFDCPseudo * > > &LinkedSegments)
double sqrt(double)
double sin(double)
#define S(str)
Definition: hddm-c.cpp:84
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
Definition: DGeometry.cc:1078
jerror_t DoFilter(DMatrix4x1 &S, vector< const DFDCPseudo * > &fdchits)
bool GetCDCEndplate(double &z, double &dz, double &rmin, double &rmax) const
Definition: DGeometry.cc:1497
TDirectory * dir
Definition: bcal_hist_eff.C:31
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)
TH2F * temp
#define I(x, y, z)
#define G(x, y, z)
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)