Bug Summary

File:libraries/FDC/DFDCPseudo_factory.cc
Location:line 666, column 11
Description:Assigned value is garbage or undefined

Annotated Source Code

1//********************************************************
2// DFDCPseudo_factory.cc - factory producing first-order
3// reconstructed points for the FDC
4// Author: Craig Bookwalter (craigb at jlab.org)
5// Date: March 2006
6// UVX cathode-wire matching revised by Simon Taylor, Aug 2006
7//********************************************************
8
9#include "DFDCPseudo_factory.h"
10#include "HDGEOMETRY/DGeometry.h"
11#include "DFDCGeometry.h"
12#include <TRACKING/DTrackHitSelectorTHROWN.h>
13#include <TROOT.h>
14#include <FDC/DFDCCathodeDigiHit.h>
15
16#define HALF_CELL0.5 0.5
17#define MAX_DEFLECTION0.15 0.15
18#define X00 0
19#define QA1 1
20#define K22 2
21#define ITER_MAX100 100
22#define TOLX1e-4 1e-4
23#define TOLF1e-4 1e-4
24#define A_OVER_H0.4 0.4
25#define ONE_OVER_H2.0 2.0
26#define ALPHA1e-4 1e-4 // rate parameter for Newton step backtracking algorithm
27#define W_EFF30.2e-9 30.2e-9 // GeV
28#define GAS_GAIN8e4 8e4
29#define ELECTRON_CHARGE1.6022e-4 1.6022e-4 // fC
30
31
32///
33/// DFDCAnode_gLayer_cmp():
34/// non-member function passed to std::sort() to sort DFDCHit pointers
35/// for the anode wires by their gLayer attributes.
36///
37bool DFDCAnode_gLayer_cmp(const DFDCHit* a, const DFDCHit* b) {
38 return a->gLayer < b->gLayer;
39}
40
41
42
43bool DFDCPseudo_cmp(const DFDCPseudo* a, const DFDCPseudo *b){
44 if (a->wire->wire == b->wire->wire && a->wire->layer==b->wire->layer){
45 return a->time<b->time;
46 }
47 if (a->wire->layer!=b->wire->layer) return a->wire->layer<b->wire->layer;
48
49 return a->wire->wire<b->wire->wire;
50}
51
52
53
54
55///
56/// DFDCPseudo_factory::DFDCPseudo_factory():
57/// default constructor -- initializes log file
58///
59DFDCPseudo_factory::DFDCPseudo_factory() {
60 //_log = new JStreamLog(std::cout, "FDC PSEUDO >>");
61 //*_log << "File initialized." << endMsg;
62}
63
64///
65/// DFDCPseudo_factory::~DFDCPseudo_factory():
66/// default destructor -- closes log file
67///
68DFDCPseudo_factory::~DFDCPseudo_factory() {
69 if (fdcwires.size()){
70 for (unsigned int i=0;i<fdcwires.size();i++){
71 for (unsigned int j=0;j<fdcwires[i].size();j++){
72 delete fdcwires[i][j];
73 }
74 }
75 }
76 if (fdccathodes.size()){
77 for (unsigned int i=0;i<fdccathodes.size();i++){
78 for (unsigned int j=0;j<fdccathodes[i].size();j++){
79 delete fdccathodes[i][j];
80 }
81 }
82 }
83 //delete _log;
84}
85
86//------------------
87// init
88//------------------
89jerror_t DFDCPseudo_factory::init(void)
90{
91 RIN_FIDUCIAL = 1.5;
92 ROUT_FIDUCIAL=48.0;
93 MAX_ALLOWED_FDC_HITS=20000;
94 STRIP_ANODE_TIME_CUT=10.;
95 MIDDLE_STRIP_THRESHOLD=0.;
96
97 r2_out=ROUT_FIDUCIAL*ROUT_FIDUCIAL;
98 r2_in=RIN_FIDUCIAL*RIN_FIDUCIAL;
99
100 gPARMS->SetDefaultParameter("FDC:ROUT_FIDUCIAL",ROUT_FIDUCIAL, "Outer fiducial radius of FDC in cm");
101 gPARMS->SetDefaultParameter("FDC:RIN_FIDUCIAL",RIN_FIDUCIAL, "Inner fiducial radius of FDC in cm");
102 gPARMS->SetDefaultParameter("FDC:MAX_ALLOWED_FDC_HITS",MAX_ALLOWED_FDC_HITS, "Max. number of FDC hits (includes both cathode strips and wires hits) to allow before considering event too busy to attempt FDC tracking");
103 gPARMS->SetDefaultParameter("FDC:STRIP_ANODE_TIME_CUT",STRIP_ANODE_TIME_CUT,
104 "maximum time difference between strips and wires (in ns)");
105
106 DEBUG_HISTS = false;
107 gPARMS->SetDefaultParameter("FDC:DEBUG_HISTS",DEBUG_HISTS);
108
109
110 return NOERROR;
111}
112
113
114//------------------
115// brun
116//------------------
117jerror_t DFDCPseudo_factory::brun(JEventLoop *loop, int32_t runnumber)
118{
119 // Get pointer to DGeometry object
120 DApplication* dapp=dynamic_cast<DApplication*>(eventLoop->GetJApplication());
121 const DGeometry *dgeom = dapp->GetDGeometry(runnumber);
122
123 USE_FDC=true;
124 if (!dgeom->GetFDCWires(fdcwires)){
125 _DBG_std::cerr<<"libraries/FDC/DFDCPseudo_factory.cc"<<
":"<<125<<" "<< "FDC geometry not available!" <<endl;
126 USE_FDC=false;
127 }
128 if (!dgeom->GetFDCCathodes(fdccathodes)){
129 _DBG_std::cerr<<"libraries/FDC/DFDCPseudo_factory.cc"<<
":"<<129<<" "<< "FDC geometry not available!" <<endl;
130 USE_FDC=false;
131 }
132
133 // Get offsets tweaking nominal geometry from calibration database
134 JCalibration * jcalib = dapp->GetJCalibration(runnumber);
135 vector<map<string,double> >vals;
136 if (jcalib->Get("FDC/cell_offsets",vals)==false){
137 for(unsigned int i=0; i<vals.size(); i++){
138 map<string,double> &row = vals[i];
139
140 // Get the offsets from the calibration database
141 xshifts.push_back(row["xshift"]);
142 yshifts.push_back(row["yshift"]);
143 }
144 }
145 // Get FDC resolution parameters from database
146 map<string, double> fdcparms;
147 FDC_RES_PAR1=0.;
148 FDC_RES_PAR2=0.;
149 jcalib->Get("FDC/fdc_resolution_parms",fdcparms);
150 FDC_RES_PAR1=fdcparms["res_par1"];
151 FDC_RES_PAR2=fdcparms["res_par2"];
152
153 if(DEBUG_HISTS){
154 dapp->Lock();
155
156 // Histograms may already exist. (Another thread may have created them)
157 // Try and get pointers to the existing ones.
158 v_vs_u=(TH2F*)gROOT->FindObject("v_vs_u");
159 if (!v_vs_u) v_vs_u=new TH2F("v_vs_u","v vs u",192,0.5,192.5,192,0.5,192.5);
160
161 qv_vs_qu= (TH2F*) gROOT->FindObject("qv_vs_qu");
162 if (!qv_vs_qu) qv_vs_qu=new TH2F("qv_vs_qu","Anode charge from each cathode",100,0,20000,100,0,20000);
163
164 tv_vs_tu= (TH2F*) gROOT->FindObject("tv_vs_tu");
165 if (!tv_vs_tu) tv_vs_tu=new TH2F("tv_vs_tu","t(v) vs t(u)",100,-50,250,100,-50,250);
166
167 dtv_vs_dtu= (TH2F*) gROOT->FindObject("dtv_vs_dtu");
168 if (!dtv_vs_dtu) dtv_vs_dtu=new TH2F("dtv_vs_dtu","t(wire)-t(v) vs t(wire)-t(u)",100,-50,25000,100,-50,25000);
169
170 u_wire_dt_vs_wire=(TH2F *) gROOT->FindObject("u_wire_dt_vs_wire");
171 if (!u_wire_dt_vs_wire) u_wire_dt_vs_wire=new TH2F("u_wire_dt_vs_wire","wire/u cathode time difference vs wire number",
172 96,0.5,96.5,100,-50,50);
173 v_wire_dt_vs_wire=(TH2F *) gROOT->FindObject("v_wire_dt_vs_wire");
174 if (!v_wire_dt_vs_wire) v_wire_dt_vs_wire=new TH2F("v_wire_dt_vs_wire","wire/v cathode time difference vs wire number",
175 96,0.5,96.5,100,-50,50);
176 uv_dt_vs_u=(TH2F *) gROOT->FindObject("uv_dt_vs_u");
177 if (!uv_dt_vs_u) uv_dt_vs_u=new TH2F("uv_dt_vs_u","uv time difference vs u",
178 192,0.5,192.5,100,-50,50);
179 uv_dt_vs_v=(TH2F *) gROOT->FindObject("uv_dt_vs_v");
180 if (!uv_dt_vs_v) uv_dt_vs_v=new TH2F("uv_dt_vs_v","uv time difference vs v",
181 192,0.5,192.5,100,-50,50);
182
183 ut_vs_u=(TH2F *) gROOT->FindObject("ut_vs_u");
184 if (!ut_vs_u) ut_vs_u=new TH2F("ut_vs_u","u time vs u",
185 192,0.5,192.5,100,0,6000);
186 vt_vs_v=(TH2F *) gROOT->FindObject("vt_vs_v");
187 if (!vt_vs_v) vt_vs_v=new TH2F("vt_vs_v","v time vs v",
188 192,0.5,192.5,100,0,6000);
189
190 dx_vs_dE=(TH2F*)gROOT->FindObject("dx_vs_dE");
191 if (!dx_vs_dE) dx_vs_dE=new TH2F("dx_vs_dE","dx vs dE",100,0,25.0,
192 100,-0.2,0.2);
193
194
195 for (unsigned int i=0;i<24;i++){
196 char hname[80];
197 sprintf(hname,"Hxy%d",i);
198 Hxy[i]=(TH2F *) gROOT->FindObject(hname);
199 if (!Hxy[i]){
200 Hxy[i]=new TH2F(hname,hname,4000,-50,50,200,-50,50);
201 }
202 }
203
204 dapp->Unlock();
205 }
206
207 return NOERROR;
208}
209
210jerror_t DFDCPseudo_factory::erun(void){
211 if (fdcwires.size()){
212 for (unsigned int i=0;i<fdcwires.size();i++){
213 for (unsigned int j=0;j<fdcwires[i].size();j++){
214 delete fdcwires[i][j];
215 }
216 }
217 }
218 fdcwires.clear();
219 if (fdccathodes.size()){
220 for (unsigned int i=0;i<fdccathodes.size();i++){
221 for (unsigned int j=0;j<fdccathodes[i].size();j++){
222 delete fdccathodes[i][j];
223 }
224 }
225 }
226 fdccathodes.clear();
227
228
229 return NOERROR;
230}
231///
232/// DFDCPseudo_factory::evnt():
233/// this is the place that anode hits and DFDCCathodeClusters are organized into pseudopoints.
234///
235jerror_t DFDCPseudo_factory::evnt(JEventLoop* eventLoop, uint64_t eventNo) {
236 if (!USE_FDC) return NOERROR;
237
238 // Get all FDC hits (anode and cathode)
239 vector<const DFDCHit*> fdcHits;
240 eventLoop->Get(fdcHits);
241 if (fdcHits.size()==0) return NOERROR;
242
243 // For events with a very large number of hits, assume
244 // we can't reconstruct them so bail early
245 // Feb. 8, 2008 D.L. (updated to config param. Nov. 18, 2010 D.L.)
246 if(fdcHits.size()>MAX_ALLOWED_FDC_HITS){
247 _DBG_std::cerr<<"libraries/FDC/DFDCPseudo_factory.cc"<<
":"<<247<<" "<<"Too many hits in FDC ("<<fdcHits.size()<<", max="<<MAX_ALLOWED_FDC_HITS<<")! Pseudopoint reconstruction in FDC bypassed for event "<<eventLoop->GetJEvent().GetEventNumber()<<endl;
248 return NOERROR;
249 }
250
251 // Get cathode clusters
252 vector<const DFDCCathodeCluster*> cathClus;
253 eventLoop->Get(cathClus);
254 if (cathClus.size()==0) return NOERROR;
255
256 // Sift through hits and select out anode hits.
257 vector<const DFDCHit*> xHits;
258 for (unsigned int i=0; i < fdcHits.size(); i++)
259 if (fdcHits[i]->type == 0)
260 xHits.push_back(fdcHits[i]);
261 // Make sure the wires are also in order of ascending z position
262 std::sort(xHits.begin(), xHits.end(), DFDCAnode_gLayer_cmp);
263
264 // Sift through clusters and put U and V clusters into respective vectors.
265 vector<const DFDCCathodeCluster*> uClus;
266 vector<const DFDCCathodeCluster*> vClus;
267 for (unsigned int i=0; i < cathClus.size(); i++) {
268 if (cathClus[i]->plane == 1)
269 vClus.push_back(cathClus[i]);
270 else
271 uClus.push_back(cathClus[i]);
272 }
273
274 // If this is simulated data then we want to match up the truth hit
275 // with this "real" hit. Ideally, this would be done at the
276 // DFDCHit object level, but the organization of the data in HDDM
277 // makes that difficult. Here we have the full wire definition so
278 // we make the connection here.
279 vector<const DMCTrackHit*> mctrackhits;
280 eventLoop->Get(mctrackhits);
281
282 vector<const DFDCCathodeCluster*>::iterator uIt = uClus.begin();
283 vector<const DFDCCathodeCluster*>::iterator vIt = vClus.begin();
284 vector<const DFDCHit*>::iterator xIt = xHits.begin();
285
286 // For each layer, get its sets of V, X, and U hits, and then pass them to the geometrical
287 // organization routine, DFDCPseudo_factory::makePseudo()
288 vector<const DFDCCathodeCluster*> oneLayerU;
289 vector<const DFDCCathodeCluster*> oneLayerV;
290 vector<const DFDCHit*> oneLayerX;
291 for (int iLayer=1; iLayer <= 24; iLayer++) {
292 for (; ((uIt != uClus.end() && (*uIt)->gLayer == iLayer)); uIt++)
293 oneLayerU.push_back(*uIt);
294 for (; ((vIt != vClus.end() && (*vIt)->gLayer == iLayer)); vIt++)
295 oneLayerV.push_back(*vIt);
296 for (; ((xIt != xHits.end() && (*xIt)->gLayer == iLayer)); xIt++)
297 oneLayerX.push_back(*xIt);
298 if (oneLayerU.size()>0 && oneLayerV.size()>0 && oneLayerX.size()>0)
299 makePseudo(oneLayerX, oneLayerU, oneLayerV,iLayer, mctrackhits);
300 oneLayerU.clear();
301 oneLayerV.clear();
302 oneLayerX.clear();
303 }
304 // Make sure the data are both time- and z-ordered
305 std::sort(_data.begin(),_data.end(),DFDCPseudo_cmp);
306
307 return NOERROR;
308}
309
310///
311/// DFDCPseudo_factory::makePseudo():
312/// performs UV+X matching to create pseudopoints
313///
314void DFDCPseudo_factory::makePseudo(vector<const DFDCHit*>& x,
315 vector<const DFDCCathodeCluster*>& u,
316 vector<const DFDCCathodeCluster*>& v,
317 int layer,
318 vector<const DMCTrackHit*> &mctrackhits)
319{
320 vector<const DFDCHit*>::iterator xIt;
321 vector<centroid_t>upeaks;
322 vector<centroid_t>vpeaks;
323
324
325 // printf("---------u clusters --------\n");
326 // Loop over all U and V clusters looking for peaks
327 for (unsigned int i=0;i<u.size();i++){
328 //printf("Cluster %d\n",i);
329 if (u[i]->members.size()>2)
330 {
331 for (vector<const DFDCHit*>::const_iterator strip=u[i]->members.begin();
332 strip!=u[i]->members.end();strip++){
333 //printf(" %d %f %f\n",(*strip)->element,(*strip)->q,(*strip)->t);
334 if (FindCentroid(u[i]->members,strip,upeaks)==NOERROR){
335 upeaks[upeaks.size()-1].cluster=i;
336 }
337 }
338 }
339 }
340 //printf("---------v cluster --------\n");
341 for (unsigned int i=0;i<v.size();i++){
342 //printf("Cluster %d\n",i);
343 if (v[i]->members.size()>2)
344 {
345 for (vector<const DFDCHit*>::const_iterator strip=v[i]->members.begin();
346 strip!=v[i]->members.end();strip++){
347 //printf(" %d %f %f\n",(*strip)->element,(*strip)->q,(*strip)->t);
348 if (FindCentroid(v[i]->members,strip,vpeaks)==NOERROR){
349 vpeaks[vpeaks.size()-1].cluster=i;
350 }
351 }
352 }
353 }
354 if (upeaks.size()*vpeaks.size()>0){
355 // Rotation angles for strips
356 unsigned int ilay=layer-1;
357 unsigned int ind=2*ilay;
358 float phi_u=fdccathodes[ind][0]->angle;
359 float phi_v=fdccathodes[ind+1][0]->angle;
360
361 //Loop over all u and v centroids looking for matches with wires
362 for (unsigned int i=0;i<upeaks.size();i++){
363 for (unsigned int j=0;j<vpeaks.size();j++){
364 // In the layer local coordinate system, wires are quantized
365 // in the x-direction and y is along the wire.
366 double x_from_strips=DFDCGeometry::getXLocalStrips(upeaks[i].pos,phi_u,
367 vpeaks[j].pos,phi_v);
368 double y_from_strips=DFDCGeometry::getYLocalStrips(upeaks[i].pos,phi_u,
369 vpeaks[j].pos,phi_v);
370 int old_wire_num=-1;
371 for(xIt=x.begin();xIt!=x.end();xIt++){
372 if ((*xIt)->element<=WIRES_PER_PLANE96 && (*xIt)->element>0){
373 const DFDCWire *wire=fdcwires[layer-1][(*xIt)->element-1];
374 double x_from_wire=wire->u;
375
376 //printf("xs %f xw %f\n",x_from_strips,x_from_wire);
377
378 // Test radial value for checking whether or not the hit is within
379 // the fiducial region of the detector
380 double r2test=x_from_wire*x_from_wire+y_from_strips*y_from_strips;
381 double delta_x=x_from_wire-x_from_strips;
382
383 if (old_wire_num==(*xIt)->element) continue;
384 old_wire_num=(*xIt)->element;
385
386 if (fabs(delta_x)<0.5*WIRE_SPACING1.0 && r2test<r2_out
387 && r2test>r2_in){
388 double dt1 = (*xIt)->t - upeaks[i].t;
389 double dt2 = (*xIt)->t - vpeaks[j].t;
390
391 //printf("dt1 %f dt2 %f\n",dt1,dt2);
392
393 if (DEBUG_HISTS){
394 if (layer==1){
395 dtv_vs_dtu->Fill(dt1,dt2);
396 tv_vs_tu->Fill(upeaks[i].t, vpeaks[j].t);
397 u_wire_dt_vs_wire->Fill((*xIt)->element,(*xIt)->t-upeaks[i].t);
398 v_wire_dt_vs_wire->Fill((*xIt)->element,(*xIt)->t-vpeaks[j].t);
399
400
401
402 int uid=u[upeaks[i].cluster]->members[1]->element;
403 int vid=v[vpeaks[j].cluster]->members[1]->element;
404
405 const DFDCCathodeDigiHit *vdigihit;
406 v[vpeaks[j].cluster]->members[1]->GetSingle(vdigihit);
407 const DFDCCathodeDigiHit *udigihit;
408 u[upeaks[i].cluster]->members[1]->GetSingle(udigihit);
409 if (vdigihit!=NULL__null && udigihit!=NULL__null){
410 int dt=udigihit->pulse_time-vdigihit->pulse_time;
411 // printf("%d %d\n",udigihit->pulse_time,vdigihit->pulse_time);
412 uv_dt_vs_u->Fill(uid,dt);
413 uv_dt_vs_v->Fill(vid,dt);
414 v_vs_u->Fill(uid,vid);
415 ut_vs_u->Fill(uid,udigihit->pulse_time);
416 vt_vs_v->Fill(vid,udigihit->pulse_time);
417 }
418 // Hxy->Fill(x_from_strips,y_from_strips);
419 }
420 }
421 // if (sqrt(dt1*dt1+dt2*dt2)>STRIP_ANODE_TIME_CUT) continue;
422
423 // Temporary cut until TDC timing is worked out
424 if (fabs(vpeaks[j].t-upeaks[i].t)>STRIP_ANODE_TIME_CUT) continue;
425
426 // Charge and energy loss
427 double q_cathodes=0.5*(upeaks[i].q+vpeaks[j].q);
428 double charge_to_energy=W_EFF30.2e-9/(GAS_GAIN8e4*ELECTRON_CHARGE1.6022e-4);
429 double dE=charge_to_energy*q_cathodes;
430 double q_from_pulse_height=5.0e-4*(upeaks[i].q_from_pulse_height
431 +vpeaks[j].q_from_pulse_height);
432
433 if (DEBUG_HISTS){
434 qv_vs_qu->Fill(upeaks[i].q,vpeaks[j].q);
435 }
436
437 int status=upeaks[i].numstrips+vpeaks[j].numstrips;
438 //double xres=WIRE_SPACING/2./sqrt(12.);
439
440 DFDCPseudo* newPseu = new DFDCPseudo;
441 newPseu->phi_u=phi_u;
442 newPseu->phi_v=phi_v;
443 newPseu->u = upeaks[i].pos;
444 newPseu->v = vpeaks[j].pos;
445 newPseu->w = x_from_wire-xshifts[ilay];
446 newPseu->dw = 0.; // place holder
447 newPseu->w_c = x_from_strips-xshifts[ilay];
448 newPseu->s = y_from_strips-yshifts[ilay];
449 newPseu->ds = FDC_RES_PAR1/q_from_pulse_height+FDC_RES_PAR2;
450 //newPseu->ds=0.011/q_from_pulse_height+5e-3+2.14e-10*pow(q_from_pulse_height,6);
451 newPseu->wire = wire;
452 //newPseu->time = (*xIt)->t;
453 newPseu->time=0.5*(upeaks[i].t+vpeaks[j].t);
454 newPseu->status = status;
455 newPseu->itrack = (*xIt)->itrack;
456
457 newPseu->AddAssociatedObject(v[vpeaks[j].cluster]);
458 newPseu->AddAssociatedObject(u[upeaks[i].cluster]);
459
460 newPseu->dE = dE;
461 newPseu->q = q_cathodes;
462
463 // It can occur (although rarely) that newPseu->wire is NULL
464 // which causes us to crash below. In these cases, we can't really
465 // make a psuedo point so we delete the current object
466 // and just go on to the next one.
467 if(newPseu->wire==NULL__null){
468 _DBG_std::cerr<<"libraries/FDC/DFDCPseudo_factory.cc"<<
":"<<468<<" "<<"newPseu->wire=NULL! This shouldn't happen. Complain to staylor@jlab.org"<<endl;
469 delete newPseu;
470 continue;
471 }
472 double sinangle=newPseu->wire->udir(0);
473 double cosangle=newPseu->wire->udir(1);
474
475 newPseu->xy.Set((newPseu->w)*cosangle+(newPseu->s)*sinangle,
476 -(newPseu->w)*sinangle+(newPseu->s)*cosangle);
477
478 double sigx2=HALF_CELL0.5*HALF_CELL0.5/3.;
479 double sigy2=MAX_DEFLECTION0.15*MAX_DEFLECTION0.15/3.;
480 newPseu->covxx=sigx2*cosangle*cosangle+sigy2*sinangle*sinangle;
481 newPseu->covyy=sigx2*sinangle*sinangle+sigy2*cosangle*cosangle;
482 newPseu->covxy=(sigy2-sigx2)*sinangle*cosangle;
483
484 // Try matching truth hit with this "real" hit.
485 const DMCTrackHit *mctrackhit = DTrackHitSelectorTHROWN::GetMCTrackHit(newPseu->wire, DRIFT_SPEED.0055*newPseu->time, mctrackhits);
486 if(mctrackhit)newPseu->AddAssociatedObject(mctrackhit);
487
488 _data.push_back(newPseu);
489
490 if (DEBUG_HISTS){
491 Hxy[ilay]->Fill(newPseu->w_c,newPseu->s);
492 if (ilay==6) dx_vs_dE->Fill(q_from_pulse_height,0.2588*delta_x);
493 }
494
495 } // match in x
496 } else _DBG_std::cerr<<"libraries/FDC/DFDCPseudo_factory.cc"<<
":"<<496<<" " << "Bad wire " << (*xIt)->element <<endl;
497 } // xIt loop
498 } // vpeaks loop
499 } // upeaks loop
500 } // if we have peaks in both u and v views
501
502}
503
504//
505/// DFDCPseudo_factory::CalcMeanTime()
506/// Calculate the mean and rms of the times of the hits passed in "H".
507/// The contents of H should be pointers to a single cluster in a
508/// cathode plane.
509// 1/2/2008 D.L.
510void DFDCPseudo_factory::CalcMeanTime(const vector<const DFDCHit*>& H, double &t, double &t_rms)
511{
512 // Calculate mean
513 t=0.0;
514 for(unsigned int i=0; i<H.size(); i++)t+=H[i]->t;
515 if(H.size()>0)t/=(double)H.size();
516
517 // Calculate RMS
518 t_rms=0.0;
519 for(unsigned int i=0; i<H.size(); i++)t_rms+=pow((double)(H[i]->t-t),2.0);
520 if(H.size()>0)t_rms = sqrt(t_rms/(double)H.size());
521}
522
523// Find the mean time and rms for a group of 3 hits with a maximum in the
524// center hit
525void DFDCPseudo_factory::CalcMeanTime(vector<const DFDCHit *>::const_iterator peak, double &t, double &t_rms)
526{
527 // Calculate mean
528 t=0.0;
529 for (vector<const DFDCHit*>::const_iterator j=peak-1;j<=peak+1;j++){
530 t+=(*j)->t;
531 }
532 t/=3.;
533
534 // Calculate RMS
535 t_rms=0.0;
536 for (vector<const DFDCHit*>::const_iterator j=peak-1;j<=peak+1;j++){
537 t_rms+=((*j)->t-t)*((*j)->t-t);
538 }
539 t_rms=sqrt(t_rms/3.);
540}
541
542
543//
544/// DFDCPseudo_factory::FindCentroid()
545/// Uses the Newton-Raphson method for solving the set of non-linear
546/// equations describing the charge distribution over 3 strips for the peak
547/// position x0, the anode charge qa, and the "width" parameter k2.
548/// See Numerical Recipes in C p.379-383.
549/// Updates list of centroids.
550///
551jerror_t DFDCPseudo_factory::FindCentroid(const vector<const DFDCHit*>& H,
552 vector<const DFDCHit *>::const_iterator peak,
553 vector<centroid_t>&centroids){
554 centroid_t temp;
555 // Make sure we do not exceed the range of the vector
556 if (peak>H.begin() && peak+1!=H.end()){
1
Taking true branch
→
557 double err_diff1=0.,err_diff2=0.;
558
559 // Fill in time info in temp 1/2/2008 D.L.
560 //CalcMeanTime(H, temp.t, temp.t_rms);
561
562 // Some code for checking for significance of fluctuations.
563 // Currently disabled.
564 //double dq1=(*(peak-1))->dq;
565 //double dq2=(*peak)->dq;
566 //double dq3=(*(peak+1))->dq;
567 //err_diff1=sqrt(dq1*dq1+dq2*dq2);
568 //err_diff2=sqrt(dq2*dq2+dq3*dq3);
569 if ((*peak)->pulse_height<MIDDLE_STRIP_THRESHOLD) {
2
←
Taking false branch
→
570 return VALUE_OUT_OF_RANGE;
571 }
572
573 // Check for a peak in three adjacent strips
574 if ((*peak)->pulse_height-(*(peak-1))->pulse_height > err_diff1
3
←
Taking true branch
→
575 && (*peak)->pulse_height-(*(peak+1))->pulse_height > err_diff2){
576 // Define some matrices for use in the Newton-Raphson iteration
577 DMatrix3x3 J; //Jacobean matrix
578 DMatrix3x1 F,N,X,par,newpar,f;
579 int i=0;
580 double sum=0.;
581
582 // Initialize the matrices to some suitable starting values
583 unsigned int index=2*((*peak)->gLayer-1)+(*peak)->plane/2;
584 par(K22)=1.;
585 for (vector<const DFDCHit*>::const_iterator j=peak-1;j<=peak+1;j++){
4
←
Loop condition is false. Execution continues on line 591
→
586 X(i)=fdccathodes[index][(*j)->element-1]->u;
587 N(i)=double((*j)->pulse_height);
588 sum+=N(i);
589 i++;
590 }
591 par(X00)=X(1);
592 par(QA1)=2.*sum;
593 newpar=par;
594
595 // Newton-Raphson procedure
596 double errf=0.,errx=0;
597 for (int iter=1;iter<=ITER_MAX100;iter++){
5
←
Loop condition is true. Entering loop body
→
598 errf=0.;
599 errx=0.;
600
601 // Compute Jacobian matrix: J_ij = dF_i/dx_j.
602 for (i=0;i<3;i++){
6
←
Loop condition is true. Entering loop body
→
8
←
Loop condition is true. Entering loop body
→
10
←
Loop condition is true. Entering loop body
→
12
←
Loop condition is false. Execution continues on line 622
→
603 double dx=(par(X00)-X(i))*ONE_OVER_H2.0;
604 double argp=par(K22)*(dx+A_OVER_H0.4);
605 double argm=par(K22)*(dx-A_OVER_H0.4);
606 double tanh_p=tanh(argp);
607 double tanh_m=tanh(argm);
608 double tanh2_p=tanh_p*tanh_p;
609 double tanh2_m=tanh_m*tanh_m;
610 double q_over_4=0.25*par(QA1);
611
612 f(i)=tanh_p-tanh_m;
613 J(i,QA1)=-0.25*f(i);
614 J(i,K22)=-q_over_4*(argp/par(K22)*(1.-tanh2_p)
615 -argm/par(K22)*(1.-tanh2_m));
616 J(i,X00)=-q_over_4*par(K22)*(tanh2_m-tanh2_p);
617 F(i)=N(i)-q_over_4*f(i);
618 double new_errf=fabs(F(i));
619 if (new_errf>errf) errf=new_errf;
7
←
Taking false branch
→
9
←
Taking false branch
→
11
←
Taking false branch
→
620 }
621 // Check for convergence
622 if (errf<TOLF1e-4){
13
←
Taking false branch
→
623 temp.pos=par(X00);
624 temp.q_from_pulse_height=par(QA1);
625 temp.numstrips=3;
626 temp.t=(*peak)->t;
627 temp.t_rms=0.;
628
629 // Find estimate for anode charge
630 sum=0;
631 double sum_f=0;
632 unsigned int k=0;
633 for (vector<const DFDCHit*>::const_iterator j=peak-1;j<=peak+1;j++){
634 sum_f+=f(k);
635 sum+=double((*j)->q);
636 k++;
637 }
638 temp.q=4.*sum/sum_f;
639
640 //CalcMeanTime(peak,temp.t,temp.t_rms);
641 centroids.push_back(temp);
642
643 return NOERROR;
644 }
645
646 // Find the new set of parameters
647 FindNewParmVec(N,X,F,J,par,newpar);
648
649 //Check for convergence
650 for (i=0;i<3;i++){
14
←
Loop condition is true. Entering loop body
→
16
←
Loop condition is true. Entering loop body
→
18
←
Loop condition is true. Entering loop body
→
20
←
Loop condition is false. Execution continues on line 654
→
651 double new_err=fabs(par(i)-newpar(i));
652 if (new_err>errx) errx=new_err;
15
←
Taking false branch
→
17
←
Taking false branch
→
19
←
Taking false branch
→
653 }
654 if (errx<TOLX1e-4){
21
←
Taking true branch
→
655 temp.pos=par(X00);
656 temp.numstrips=3;
657 temp.q_from_pulse_height=par(QA1);
658 temp.t=(*peak)->t;
659 temp.t_rms=0.;
660
661 // Find estimate for anode charge
662 sum=0;
663 double sum_f=0;
664 unsigned int k=0;
665 for (vector<const DFDCHit*>::const_iterator j=peak-1;j<=peak+1;j++){
22
←
Loop condition is true. Entering loop body
→
25
←
Loop condition is true. Entering loop body
→
28
←
Loop condition is true. Entering loop body
→
31
←
Loop condition is true. Entering loop body
→
666 sum_f+=f(k);
23
←
Calling 'DMatrix3x1::operator()'
→
24
←
Returning from 'DMatrix3x1::operator()'
→
26
←
Calling 'DMatrix3x1::operator()'
→
27
←
Returning from 'DMatrix3x1::operator()'
→
29
←
Calling 'DMatrix3x1::operator()'
→
30
←
Returning from 'DMatrix3x1::operator()'
→
32
←
Calling 'DMatrix3x1::operator()'
→
33
←
Returning from 'DMatrix3x1::operator()'
→
34
←
Assigned value is garbage or undefined
667 sum+=double((*j)->q);
668 k++;
669 }
670 temp.q=4.*sum/sum_f;
671
672 //CalcMeanTime(peak,temp.t,temp.t_rms);
673 centroids.push_back(temp);
674
675 return NOERROR;
676 }
677 par=newpar;
678 } // iterations
679 }
680 }
681 else{
682 return VALUE_OUT_OF_RANGE;
683 }
684 return INFINITE_RECURSION; // error placeholder
685}
686
687
688
689///
690/// DFDCPseudo_factory::FindNewParmVec()
691/// Routine used by FindCentroid to backtrack along the direction
692/// of the Newton step if the step is too large in order to avoid conditions
693/// where the iteration procedure starts to diverge.
694/// The procedure uses a scalar quantity f= 1/2 F.F for this purpose.
695/// Algorithm described in Numerical Recipes in C, pp. 383-389.
696///
697
698jerror_t DFDCPseudo_factory::FindNewParmVec(const DMatrix3x1 &N,
699 const DMatrix3x1 &X,
700 const DMatrix3x1 &F,
701 const DMatrix3x3 &J,
702 const DMatrix3x1 &par,
703 DMatrix3x1 &newpar){
704 // Invert the J matrix
705 DMatrix3x3 InvJ=J.Invert();
706
707 // Find the full Newton step
708 DMatrix3x1 fullstep=(-1.)*(InvJ*F);
709
710 // find the rate of decrease for the Newton-Raphson step
711 double slope=(-1.0)*F.Mag2(); //dot product
712
713 // This should be a negative number...
714 if (slope>=0){
715 return VALUE_OUT_OF_RANGE;
716 }
717
718 double lambda=1.0; // Start out with full Newton step
719 double lambda_temp,lambda2=lambda;
720 DMatrix3x1 newF;
721 double f2=0.,newf;
722
723 // Compute starting values for f=1/2 F.F
724 double f=-0.5*slope;
725
726 for (;;){
727 newpar=par+lambda*fullstep;
728
729 // Compute the value of the vector F and f=1/2 F.F with the current step
730 for (int i=0;i<3;i++){
731 double dx=(newpar(X00)-X(i))*ONE_OVER_H2.0;
732 double argp=newpar(K22)*(dx+A_OVER_H0.4);
733 double argm=newpar(K22)*(dx-A_OVER_H0.4);
734 newF(i)=N(i)-0.25*newpar(QA1)*(tanh(argp)-tanh(argm));
735 }
736 newf=0.5*newF.Mag2(); // dot product
737
738 if (lambda<0.1) { // make sure the step is not too small
739 newpar=par;
740 return NOERROR;
741 } // Check if we have sufficient function decrease
742 else if (newf<=f+ALPHA1e-4*lambda*slope){
743 return NOERROR;
744 }
745 else{
746 // g(lambda)=f(par+lambda*fullstep)
747 if (lambda==1.0){//first attempt: quadratic approximation for g(lambda)
748 lambda_temp=-0.5*slope/(newf-f-slope);
749 }
750 else{ // cubic approximation for g(lambda)
751 double temp1=newf-f-lambda*slope;
752 double temp2=f2-f-lambda2*slope;
753 double one_over_lambda2_sq=1./(lambda2*lambda2);
754 double one_over_lambda_sq=1./(lambda*lambda);
755 double one_over_lambda_diff=1./(lambda-lambda2);
756 double a=(temp1*one_over_lambda_sq-temp2*one_over_lambda2_sq)*one_over_lambda_diff;
757 double b=(-lambda2*temp1*one_over_lambda_sq+lambda*temp2*one_over_lambda2_sq)
758 *one_over_lambda_diff;
759 if (a==0.0) lambda_temp=-0.5*slope/b;
760 else{
761 double disc=b*b-3.0*a*slope;
762 if (disc<0.0) lambda_temp=0.5*lambda;
763 else if (b<=0.0) lambda_temp=(-b+sqrt(disc))/(3.*a);
764 else lambda_temp=-slope/(b+sqrt(disc));
765 }
766 // ensure that we are headed in the right direction...
767 if (lambda_temp>0.5*lambda) lambda_temp=0.5*lambda;
768 }
769 }
770 lambda2=lambda;
771 f2=newf;
772 // Make sure that new version of lambda is not too small
773 lambda=(lambda_temp>0.1*lambda ? lambda_temp : 0.1*lambda);
774 }
775}
776