Belle II Software development
CDCGeometryPar.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9#include <framework/gearbox/GearDir.h>
10#include <framework/logging/Logger.h>
11#include <framework/utilities/FileSystem.h>
12
13#include <cdc/geometry/CDCGeometryPar.h>
14#include <cdc/geometry/CDCGeoControlPar.h>
15#include <cdc/simulation/CDCSimControlPar.h>
16#include <cdc/utilities/OpenFile.h>
17
18#include <cmath>
19#include <iomanip>
20
21#include <boost/iostreams/filtering_stream.hpp>
22
23#include <Math/ChebyshevPol.h>
24
25using namespace std;
26using namespace Belle2;
27using namespace CDC;
28
30
36
38{
39
41
42 if (gcp.getT0InputType()) {
44 if ((*m_t0FromDB).isValid()) {
45 (*m_t0FromDB).addCallback(this, &CDCGeometryPar::setT0);
46 }
47 }
48
49 if (gcp.getBwInputType()) {
51 if ((*m_badWireFromDB).isValid()) {
52 (*m_badWireFromDB).addCallback(this, &CDCGeometryPar::setBadWire);
53 }
54 }
55 // Bad boards are on DB only:
57 if ((*m_badBoardsFromDB).isValid()) {
58 (*m_badBoardsFromDB).addCallback(this, &CDCGeometryPar::setBadBoard);
59 }
60
61 if (gcp.getPropSpeedInputType()) {
63 if ((*m_propSpeedFromDB).isValid()) {
64 (*m_propSpeedFromDB).addCallback(this, &CDCGeometryPar::setPropSpeed);
65 }
66 }
67
68 if (gcp.getTwInputType()) {
70 if ((*m_timeWalkFromDB).isValid()) {
71 (*m_timeWalkFromDB).addCallback(this, &CDCGeometryPar::setTW);
72 }
73 }
74
75 if (gcp.getXtInputType()) {
77 if ((*m_xtRelFromDB).isValid()) {
78 (*m_xtRelFromDB).addCallback(this, &CDCGeometryPar::setXtRel);
79 }
80 }
81
82 if (gcp.getSigmaInputType()) {
84 if ((*m_sResolFromDB).isValid()) {
85 (*m_sResolFromDB).addCallback(this, &CDCGeometryPar::setSResol);
86 }
87 }
88
89 if (gcp.getFFactorInputType()) {
91 if ((*m_fFactorFromDB).isValid()) {
92 (*m_fFactorFromDB).addCallback(this, &CDCGeometryPar::setFFactor);
93 }
94 }
95
96 if (gcp.getChMapInputType()) {
98 if ((*m_chMapFromDB).isValid()) {
99 (*m_chMapFromDB).addCallback(this, &CDCGeometryPar::setChMap);
100 }
101 }
102
103 if (gcp.getDisplacementInputType()) {
105 if ((*m_displacementFromDB).isValid()) {
106 (*m_displacementFromDB).addCallback(this, &CDCGeometryPar::setDisplacement);
107 }
108 }
109
110 if (gcp.getAlignmentInputType()) {
112 if ((*m_alignmentFromDB).isValid()) {
113 (*m_alignmentFromDB).addCallback(this, &CDCGeometryPar::setWirPosAlignParams);
114 }
115 }
116
117 if (gcp.getMisalignment()) {
118 if (gcp.getMisalignmentInputType()) {
120 if ((*m_misalignmentFromDB).isValid()) {
121 (*m_misalignmentFromDB).addCallback(this, &CDCGeometryPar::setWirPosMisalignParams);
122 }
123 }
124 }
125
126 //TODO in future: make a new (singleton?) class and move all EDepToADC things there.
127 if (gcp.getEDepToADCInputType()) {
129 if ((*m_eDepToADCConversionsFromDB).isValid()) {
130 (*m_eDepToADCConversionsFromDB).addCallback(this, &CDCGeometryPar::setEDepToADCConversions);
131 }
132 }
133
134 clear();
135 if (geom) {
136 // B2INFO("CDCGeometryPar: Read Geometry object");
137 readFromDB(*geom);
138 } else {
139 // std::cout <<"readcalled" << std::endl;
140 // read();
141 // B2FATAL("CDCGeometryPar: Strange that readFromDB is not called !");
142 B2WARNING("CDCGeometryPar: Strange that readFromDB is not called! Please make sure that CDC is included in Geometry.");
143 }
144}
145
147{
148 // B2INFO("CDCGeometryPar: destructor called");
149 // if (m_t0FromDB) delete m_t0FromDB;
150 // if (m_badWireFromDB) delete m_badWireFromDB;
151 // if (m_propSpeedFromDB) delete m_propSpeedFromDB;
152 // if (m_timeWalkFromDB) delete m_timeWalkFromDB;
153 // if (m_xtRelFromDB) delete m_xtRelFromDB;
154 // if (m_sResolFromDB) delete m_sResolFromDB;
155 // if (m_chMapFromDB) delete [] m_chMapFromDB;
156 // if (m_displacementFromDB) delete [] m_displacementFromDB;
157 // if (m_chMapFromDB) delete m_chMapFromDB;
158 // if (m_displacementFromDB) delete m_displacementFromDB;
159 // if (m_alignmentFromDB) delete m_alignmentFromDB;
160 // if (m_misalignmentFromDB) delete m_misalignmentFromDB;
161 // B2INFO("CDCGeometryPar: destructor ended");
162}
163
165{
166 m_version = "unknown";
167 m_nSLayer = 0;
168 m_nFLayer = 0;
170 m_senseWireTension = 0.0;
171 m_senseWireDensity = 0.0;
173
174 m_tdcOffset = 0; //not used; to be removed later
175 m_clockFreq4TDC = 0.0;
176 m_tdcBinWidth = 0.0;
177 m_nominalDriftV = 0.0;
178 m_nominalPropSpeed = 0.0;
180
181 for (unsigned i = 0; i < 4; ++i) {
182 m_rWall[i] = 0;
183 for (unsigned j = 0; j < 2; ++j)
184 m_zWall[i][j] = 0;
185 }
186 for (unsigned i = 0; i < c_maxNSenseLayers; ++i) {
187 m_rSLayer[i] = 0;
188 m_zSForwardLayer[i] = 0;
189 m_dzSForwardLayer[i] = 0;
190 m_zSBackwardLayer[i] = 0;
191 m_dzSBackwardLayer[i] = 0;
192 m_cellSize[i] = 0;
193 m_nWires[i] = 0;
194 m_offSet[i] = 0;
195 m_nShifts[i] = 0;
196 m_propSpeedInv[i] = 0.;
197 }
198 for (unsigned i = 0; i < c_maxNFieldLayers; ++i) {
199 m_rFLayer[i] = 0;
200 m_zFForwardLayer[i] = 0;
201 m_zFBackwardLayer[i] = 0;
202 }
203
204 for (unsigned L = 0; L < c_maxNSenseLayers; ++L) {
205 for (unsigned C = 0; C < c_maxNDriftCells; ++C) {
206 for (unsigned i = 0; i < 3; ++i) {
207 m_FWirPos [L][C][i] = 0.;
208 m_BWirPos [L][C][i] = 0.;
209 m_FWirPosMisalign[L][C][i] = 0.;
210 m_BWirPosMisalign[L][C][i] = 0.;
211 m_FWirPosAlign [L][C][i] = 0.;
212 m_BWirPosAlign [L][C][i] = 0.;
213 }
214 for (unsigned i = 0; i < 7; ++i) {
215 m_eDepToADCParams[L][C][i] = 0.;
216 }
217 m_WireSagCoef [L][C] = 0.;
218 m_WireSagCoefMisalign[L][C] = 0.;
219 m_WireSagCoefAlign [L][C] = 0.;
220 m_t0 [L][C] = 0.;
221 }
222 }
223
224 for (unsigned L = 0; L < c_maxNSenseLayers; ++L) {
225 for (unsigned i = 0; i < 2; ++i) {
226 for (unsigned alpha = 0; alpha < c_maxNAlphaPoints; ++alpha) {
227 for (unsigned theta = 0; theta < c_maxNThetaPoints; ++theta) {
228 for (unsigned xtparam = 0; xtparam < c_nXTParams; ++xtparam) {
229 m_XT[L][i][alpha][theta][xtparam] = 0.;
230 }
231
232 for (unsigned sigmaparam = 0; sigmaparam < c_nSigmaParams; ++sigmaparam) {
233 m_Sigma[L][i][alpha][theta][sigmaparam] = 0.;
234 }
235 }
236 }
237 }
238 }
239
240 for (unsigned board = 0; board < c_nBoards; ++board) {
241 for (unsigned i = 0; i < 2; ++i) {
242 m_timeWalkCoef[board][i] = 0.;
243 }
244 for (unsigned channel = 0; channel < 48; ++channel) {
245 m_boardAndChannelToWire[board][channel] = 0.;
246 }
247 }
248
249 for (unsigned superLayer = 0; superLayer < c_nSuperLayers; ++superLayer) {
250 for (unsigned layer = 0; layer < 8; ++layer) {
251 m_shiftInSuperLayer[superLayer][layer] = 0;
252 }
253 }
254
255}
256
258{
259 m_globalPhiRotation = geom.getGlobalPhiRotation();
260 // m_globalPhiRotation = geom.getGlobalOffsetC();
261
262 // Get inner wall parameters
263 m_rWall[0] = geom.getInnerWall(2).getRmin();
264 m_zWall[0][0] = geom.getInnerWall(0).getZbwd();
265 m_zWall[0][1] = geom.getInnerWall(0).getZfwd();
266
267 m_rWall[1] = geom.getInnerWall(0).getRmax();
268 m_zWall[1][0] = geom.getInnerWall(0).getZbwd();
269 m_zWall[1][1] = geom.getInnerWall(0).getZbwd();
270
271 // Get outer wall parameters
272 m_rWall[2] = geom.getOuterWall(0).getRmin();
273 m_zWall[2][0] = geom.getOuterWall(0).getZbwd();
274 m_zWall[2][1] = geom.getOuterWall(0).getZfwd();
275
276 m_rWall[3] = geom.getOuterWall(1).getRmax();
277 m_zWall[3][0] = geom.getOuterWall(0).getZbwd();
278 m_zWall[3][1] = geom.getOuterWall(0).getZfwd();
279
280 // Get sense layers parameters
282 m_nSLayer = geom.getNSenseLayers();
283
285 // std::cout << m_materialDefinitionMode << std::endl;
286 if (m_materialDefinitionMode == 0) {
287 B2DEBUG(100, "CDCGeometryPar: Define a mixture of gases and wires in the tracking volume.");
288 } else if (m_materialDefinitionMode == 2) {
289 // B2INFO("CDCGeometryPar: Define all sense and field wires explicitly in the tracking volume.");
290 B2FATAL("CDCGeometryPar: Materialdefinition=2 is disabled for now.");
291 } else {
292 B2FATAL("CDCGeometryPar: Materialdefinition mode you specify is invalid.");
293 }
294
295 // Get mode for wire z-position
297 //Set z corrections (from input data)
298 B2DEBUG(100, "CDCGeometryPar: Sense wire z mode:" << m_senseWireZposMode);
299
300 //
301 // The DB version should be implemented ASAP.
302 //
303 GearDir content = GearDir("/Detector/DetectorComponent[@name=\"CDC\"]/Content/");
304 GearDir gbxParams(content);
305 // if (m_senseWireZposMode == 1) readDeltaz(gbxParams);
306
307
308 //
309 // Sense wires.
310 //
311 for (const auto& sense : geom.getSenseLayers()) {
312 uint layerId = sense.getId();
313
314 m_rSLayer[layerId] = sense.getR();
315 m_zSBackwardLayer[layerId] = sense.getZbwd();
316 m_zSForwardLayer[layerId] = sense.getZfwd();
317 m_nWires[layerId] = sense.getNWires();
318 m_nShifts[layerId] = sense.getNShifts();
319 m_offSet[layerId] = sense.getOffset();
320 m_cellSize[layerId] = 2 * M_PI * m_rSLayer[layerId] / (double) m_nWires[layerId];
321 m_dzSBackwardLayer[layerId] = sense.getDZbwd();
322 m_dzSForwardLayer[layerId] = sense.getDZfwd();
323
324 //correction to z-position
325 if (m_senseWireZposMode == 0) {
326 } else if (m_senseWireZposMode == 1) {
327 // m_zSBackwardLayer[layerId] += m_bwdDz[layerId];
328 // m_zSForwardLayer [layerId] += m_fwdDz[layerId];
329 m_zSBackwardLayer[layerId] += m_dzSBackwardLayer[layerId];
330 m_zSForwardLayer [layerId] -= m_dzSForwardLayer [layerId];
331 } else {
332 B2FATAL("CDCGeometryPar: invalid wire z definition mode specified");
333 }
334
335 //Set design sense-wire related params.
336 const int nWires = m_nWires[layerId];
337 for (int iCell = 0; iCell < nWires; ++iCell) {
338 setDesignWirParam(layerId, iCell);
339 }
340
341 }
342
343 // Get field layers parameters
344 for (const auto& field : geom.getFieldLayers()) {
345 uint layerId = field.getId();
346
347 m_rFLayer[layerId] = field.getR();
348 m_zFBackwardLayer[layerId] = field.getZbwd();
349 m_zFForwardLayer[layerId] = field.getZfwd();
350 }
351
352 // Get sense wire diameter
353 m_senseWireDiameter = geom.getSenseDiameter();
354
355 // Get sense wire tension
356 m_senseWireTension = geom.getSenseTension();
357
358 // // Get sense wire density
359 m_senseWireDensity = 19.3; // g/cm3 <- tentatively hard-coded here
360
361 // Get field wire diameter
362 m_fieldWireDiameter = geom.getFieldDiameter();
363
364 // Get information on the number of (super) layers etc.
365 m_nSenseWires = geom.getNSenseWires();
366 m_nFieldWires = geom.getNFieldWires();
367 m_maxNSenseLayers = geom.getNumberOfSenseLayers();
368 m_maxNFieldLayers = geom.getNumberOfFieldLayers();
369 m_maxNSuperLayers = geom.getMaxNumberOfSuperLayers();
370 m_firstLayerOffset = geom.getOffsetOfFirstLayer();
371 m_firstSuperLayerOffset = geom.getOffsetOfFirstSuperLayer();
372 m_maxNCellsPerLayer = geom.getMaxNumberOfCellsPerLayer();
373
374 //Set various quantities (should be moved to CDC.xml later...)
375 m_clockFreq4TDC = geom.getClockFrequency();
376 if (not m_clockSettings.isValid())
377 B2FATAL("HardwareClockSettings payloads are not valid.");
378 const double officialClockFreq4TDC = 2 * m_clockSettings->getAcceleratorRF(); // in GHz
379 if (abs(m_clockFreq4TDC - officialClockFreq4TDC) / m_clockFreq4TDC > 1.e-4) {
380 B2WARNING("ClockFreq4TDC changed from cdclocal " << scientific << setprecision(6) << m_clockFreq4TDC << " to official " <<
381 officialClockFreq4TDC << " (GHz) (difference larger than 0.01%)");
382 m_clockFreq4TDC = officialClockFreq4TDC;
383 }
384 B2DEBUG(100, "CDCGeometryPar: Clock freq. for TDC= " << m_clockFreq4TDC << " (GHz).");
385 m_tdcBinWidth = 1. / m_clockFreq4TDC; //in ns
386 B2DEBUG(100, "CDCGeometryPar: TDC bin width= " << m_tdcBinWidth << " (ns).");
387
388 m_nominalDriftV = 4.e-3; //in cm/ns
389 m_nominalDriftVInv = 1. / m_nominalDriftV; //in ns/cm
390 m_nominalPropSpeed = 27.25; //in cm/nsec (Belle's result, provided by iwasaki san)
391
392 m_nominalSpaceResol = geom.getNominalSpaceResolution();
394
395 //Set displacement params. (from input data)
397 B2DEBUG(100, "CDCGeometryPar: Load displacement params. (=1); not load (=0):" << m_displacement);
398 if (m_displacement) {
399 if (gcp.getDisplacementInputType()) {
400 B2DEBUG(100, "CDCGeometryPar: Read displacement from DB");
402 } else {
403 readWirePositionParams(c_Base, &geom);
404 }
405 }
406
407 //Set alignment params. (from input data)
409 B2DEBUG(100, "CDCGeometryPar: Load alignment params. (=1); not load (=0):" <<
411 if (m_alignment) {
412 if (gcp.getAlignmentInputType()) {
413 B2DEBUG(100, "CDCGeometryPar: Read alignment from DB");
415 } else {
416 readWirePositionParams(c_Aligned, &geom);
417 }
418 }
419
420 //Set misalignment params. (from input data)
422 B2DEBUG(100, "CDCGeometryPar: Load misalignment params. (=1); not load (=0):" <<
424 if (m_misalignment) {
425 if (gcp.getMisalignmentInputType()) {
426 B2DEBUG(100, "CDCGeometryPar: Read misalignment from DB");
428 } else {
429 readWirePositionParams(c_Misaligned, &geom);
430 }
431 }
432
433 // Get control params. for CDC FullSim
438 if (m_modLeftRightFlag) {
439 B2FATAL("ModifiedLeftRightFlag = true is disabled for now; need to update a G4-related code in framework...");
440 }
441 //N.B. The following two lines are hard-coded since only =1 are used now.
442 m_xtFileFormat = 1;
444
445 m_XTetc = true;
446 if (m_XTetc) {
447 if (gcp.getXtInputType()) {
448 B2DEBUG(100, "CDCGeometryPar: Read xt from DB");
449 setXtRel(); //Set xt param. (from DB)
450 } else {
451 readXT(gbxParams); //Read xt params. (from file)
452 }
453
454 if (gcp.getSigmaInputType()) {
455 B2DEBUG(100, "CDCGeometryPar: Read sigma from DB");
456 setSResol(); //Set sigma param. (from DB)
457 } else {
458 readSigma(gbxParams); //Read sigma params. (from file)
459 }
460
461 if (gcp.getFFactorInputType()) {
462 B2DEBUG(100, "CDCGeometryPar: Read fudge factors from DB");
463 setFFactor(); //Set fudge factors (from DB)
464 } else {
465 readFFactor(gbxParams); //Read fudge factors (from file)
466 }
467
468 if (gcp.getPropSpeedInputType()) {
469 B2DEBUG(100, "CDCGeometryPar: Read prop-speed from DB");
470 setPropSpeed(); //Set prop-speed (from DB)
471 } else {
472 readPropSpeed(gbxParams); //Read propagation speed
473 }
474
475 if (gcp.getT0InputType()) {
476 B2DEBUG(100, "CDCGeometryPar: Read t0 from DB");
477 setT0(); //Set t0 (from DB)
478 } else {
479 readT0(gbxParams); //Read t0 (from file)
480 }
481
482 if (gcp.getBwInputType()) {
483 B2DEBUG(100, "CDCGeometryPar: Read badwire from DB");
484 setBadWire(); //Set bad-wire (from DB)
485 } else {
486 // readBadWire(gbxParams); //Read bad-wire (from file)
487 B2FATAL("Text file input mode for bdwires is disabled now!");
488 }
489
490 if (gcp.getChMapInputType()) {
491 B2DEBUG(100, "CDCGeometryPar: Read ch-map from DB");
492 setChMap(); //Set ch-map (from DB)
493 } else {
494 readChMap(); //Read ch-map
495 }
496
497 if (gcp.getTwInputType()) {
498 B2DEBUG(100, "CDCGeometryPar: Read time-walk from DB");
499 setTW(); //Set time-walk coeffs. (from DB)
500 } else {
501 readTW(gbxParams); //Read time-walk coeffs. (from file)
502 }
503 B2DEBUG(100, "CDCGeometryPar: Time-walk param. mode= " << m_twParamMode);
504
505 if (gcp.getEDepToADCInputType()) {
506 B2DEBUG(29, "CDCGeometryPar: Read EDepToADC from DB");
507 if ((*m_eDepToADCConversionsFromDB).isValid()) {
508 setEDepToADCConversions(); //Set edep-to-adc (from DB)
509 }
510 } else {
511 readEDepToADC(gbxParams); //Read edep-to-adc params. (from file)
512 }
513 }
514
515 m_XTetc4Recon = 0;
516 // B2INFO("CDCGeometryPar: Load x-t etc. params. for reconstruction (=1); not load and use the same ones for digitization (=0):" <<
517 // B2INFO("CDCGeometryPar: Use the same x-t etc. for reconstruction as those used for digitization");
518 if (m_XTetc4Recon) {
519 readXT(gbxParams, 1);
520 readSigma(gbxParams, 1);
521 readPropSpeed(gbxParams, 1);
522 readT0(gbxParams, 1);
523 readTW(gbxParams, 1);
524 }
525
526 //calculate and save shifts in super-layers
528
529}
530
531/*
532//TODO: move the following two functions to cdc/utilities <-done
533// Open a file
534void CDCGeometryPar::openFile(std::ifstream& ifs, const std::string& fileName0) const
535{
536 std::string fileName1 = "/data/cdc/" + fileName0;
537 std::string fileName = FileSystem::findFile(fileName1);
538
539 if (fileName == "") {
540 fileName = FileSystem::findFile(fileName0);
541 }
542
543 if (fileName == "") {
544 B2FATAL("CDCGeometryPar: " << fileName1 << " not exist!");
545 } else {
546 B2INFO("CDCGeometryPar: open " << fileName1);
547 ifs.open(fileName.c_str());
548 if (!ifs) B2FATAL("CDCGeometryPar: cannot open " << fileName1 << " !");
549 }
550}
551// Open a file using boost (to be able to read a gzipped file)
552void CDCGeometryPar::openFile(boost::iostreams::filtering_istream& ifs, const std::string& fileName0) const
553{
554 std::string fileName1 = "/data/cdc/" + fileName0;
555 std::string fileName = FileSystem::findFile(fileName1);
556
557 if (fileName == "") {
558 fileName = FileSystem::findFile(fileName0);
559 }
560
561 if (fileName == "") {
562 B2FATAL("CDCGeometryPar: " << fileName1 << " not exist!");
563 } else {
564 B2INFO("CDCGeometryPar: open " << fileName1);
565 if ((fileName.rfind(".gz") != string::npos) && (fileName.length() - fileName.rfind(".gz") == 3)) {
566 ifs.push(boost::iostreams::gzip_decompressor());
567 }
568 ifs.push(boost::iostreams::file_source(fileName));
569 if (!ifs) B2FATAL("CDCGeometryPar: cannot open " << fileName1 << " !");
570 }
571}
572*/
573
574// Read displacement or (mis)alignment params.
575//void CDCGeometryPar::readWirePositionParams(EWirePosition set, const CDCGeometry* geom, const GearDir gbxParams)
577{
578
579 std::string fileName0;
581 if (geom) {
582 if (set == c_Base) {
583 fileName0 = gcp.getDisplacementFile();
584 } else if (set == c_Misaligned) {
585 fileName0 = gcp.getMisalignmentFile();
586 } else if (set == c_Aligned) {
587 fileName0 = gcp.getAlignmentFile();
588 }
589 } else {
590 if (set == c_Base) {
591 fileName0 = gcp.getDisplacementFile();
592 } else if (set == c_Misaligned) {
593 fileName0 = gcp.getMisalignmentFile();
594 } else if (set == c_Aligned) {
595 fileName0 = gcp.getAlignmentFile();
596 }
597 }
598
599 // ifstream ifs;
600 // openFile(ifs, fileName0);
601 boost::iostreams::filtering_istream ifs;
602 openFileB(ifs, fileName0);
603
604 uint iL(0), iC(0);
605 const int np = 3;
606 double back[np], fwrd[np], tension;
607 unsigned nRead = 0;
608
609 while (true) {
610 ifs >> iL >> iC;
611 for (int i = 0; i < np; ++i) {
612 ifs >> back[i];
613 }
614 for (int i = 0; i < np; ++i) {
615 ifs >> fwrd[i];
616 }
617 // if (set != c_Base) ifs >> tension;
618 ifs >> tension;
619
620 if (ifs.eof()) break;
621
622 if (iL < m_firstLayerOffset) {
623 continue;
624 }
625
626 ++nRead;
627
628 for (int i = 0; i < np; ++i) {
629 if (set == c_Base) {
630 m_BWirPos[iL][iC][i] += (iL < m_firstLayerOffset) ? 0 : back[i];
631 m_FWirPos[iL][iC][i] += (iL < m_firstLayerOffset) ? 0 : fwrd[i];
632 } else if (set == c_Misaligned) {
633 m_BWirPosMisalign[iL][iC][i] = m_BWirPos[iL][iC][i] + ((iL < m_firstLayerOffset) ? 0 : back[i]);
634 m_FWirPosMisalign[iL][iC][i] = m_FWirPos[iL][iC][i] + ((iL < m_firstLayerOffset) ? 0 : fwrd[i]);
635 } else if (set == c_Aligned) {
636 m_BWirPosAlign[iL][iC][i] = m_BWirPos[iL][iC][i] + ((iL < m_firstLayerOffset) ? 0 : back[i]);
637 m_FWirPosAlign[iL][iC][i] = m_FWirPos[iL][iC][i] + ((iL < m_firstLayerOffset) ? 0 : fwrd[i]);
638 }
639 }
640
641 // double baseTension = 0.;
642
643 if (set == c_Base) {
645 (8.*(m_senseWireTension + tension));
646 // std::cout <<"base iL, iC, m_senseWireTension, tension= " << iL <<" " << iC <<" "<< m_senseWireTension <<" "<< tension << std::endl;
647 } else if (set == c_Misaligned) {
648 double baseTension = M_PI * m_senseWireDensity * m_senseWireDiameter * m_senseWireDiameter / (8.* m_WireSagCoef[iL][iC]);
650 m_senseWireDiameter / (8.* (baseTension + tension));
651 // std::cout <<"misa iL, iC,basetension, tension= " << iL <<" " << iC <<" "<< baseTension <<" "<< tension << std::endl;
652 } else if (set == c_Aligned) {
653 double baseTension = M_PI * m_senseWireDensity * m_senseWireDiameter * m_senseWireDiameter / (8.* m_WireSagCoef[iL][iC]);
655 m_senseWireDiameter / (8.*(baseTension + tension));
656 // std::cout <<"algn iL, iC,basetension, tension= " << iL <<" " << iC <<" "<< baseTension <<" "<< tension << std::endl;
657 }
658 // std::cout << "baseTension,tension= " << baseTension <<" "<< tension << std::endl;
659
660 /*
661 if (m_debug) {
662 std::cout << iL << " " << iC;
663 for (int i = 0; i < np; ++i) cout << " " << back[i];
664 for (int i = 0; i < np; ++i) cout << " " << fwrd[i];
665 std::cout << " " << tension << std::endl;
666 }
667 */
668
669 }
670
671 if (nRead != m_nSenseWires) B2FATAL("CDCGeometryPar::readWirePositionParams: #lines read-in (=" << nRead <<
672 ") is inconsistent with total #sense wires (=" << m_nSenseWires << ") !");
673
674 // ifs.close();
675 boost::iostreams::close(ifs);
676}
677
678
679// Set alignment wire positions
681{
682 // Layer alignment
683 for (unsigned iL = 0; iL < c_maxNSenseLayers; ++iL) {
684
685 if (iL < m_firstLayerOffset) {
686 continue;
687 }
688
689 // wire number 511 = no wire
690 auto layerID = WireID(iL, 511);
691
692 // Alignment parameters for layer iL
693 double d_layerXbwd = (*m_alignmentFromDB)->get(layerID, CDCAlignment::layerX);
694 double d_layerYbwd = (*m_alignmentFromDB)->get(layerID, CDCAlignment::layerY);
695 double d_layerPhiBwd = (*m_alignmentFromDB)->get(layerID, CDCAlignment::layerPhi);
696
697 double d_layerXfwd = (*m_alignmentFromDB)->get(layerID, CDCAlignment::layerDx) + d_layerXbwd;
698 double d_layerYfwd = (*m_alignmentFromDB)->get(layerID, CDCAlignment::layerDy) + d_layerYbwd;
699 double d_layerPhiFwd = (*m_alignmentFromDB)->get(layerID, CDCAlignment::layerDPhi) + d_layerPhiBwd;
700
701 for (unsigned iC = 0; iC < m_nWires[iL]; ++iC) {
702 // Positions (nominal+displacement) of wire-ends of wire iC in layer iL
703 double wireXbwd = m_BWirPos[iL][iC][0];
704 double wireYbwd = m_BWirPos[iL][iC][1];
705 double wireZbwd = m_BWirPos[iL][iC][2];
706
707 double wireXfwd = m_FWirPos[iL][iC][0];
708 double wireYfwd = m_FWirPos[iL][iC][1];
709 double wireZfwd = m_FWirPos[iL][iC][2];
710
711 // Aligned positions of wire-ends are obtained by rotating "nominal+displacement" positions and shifting them using
712 // common parameters for layer rotation and shifts (at corresponding end-caps)
713 m_BWirPosAlign[iL][iC][0] = d_layerXbwd + cos(d_layerPhiBwd) * wireXbwd + sin(d_layerPhiBwd) * wireYbwd;
714 m_BWirPosAlign[iL][iC][1] = d_layerYbwd - sin(d_layerPhiBwd) * wireXbwd + cos(d_layerPhiBwd) * wireYbwd;
715 m_BWirPosAlign[iL][iC][2] = wireZbwd;
716
717 m_FWirPosAlign[iL][iC][0] = d_layerXfwd + cos(d_layerPhiFwd) * wireXfwd + sin(d_layerPhiFwd) * wireYfwd;
718 m_FWirPosAlign[iL][iC][1] = d_layerYfwd - sin(d_layerPhiFwd) * wireXfwd + cos(d_layerPhiFwd) * wireYfwd;
719 m_FWirPosAlign[iL][iC][2] = wireZfwd;
720 } //end of cell loop
721 } //end of layer loop
722
723 const int np = 3;
724 double back[np], fwrd[np];
725
726 for (unsigned iL = 0; iL < c_maxNSenseLayers; ++iL) {
727
728 if (iL < m_firstLayerOffset) {
729 continue;
730 }
731
732 for (unsigned iC = 0; iC < m_nWires[iL]; ++iC) {
733 // std::cout << "iLiC= " << iL <<" "<< iC << std::endl;
734 WireID wire(iL, iC);
735 back[0] = (*m_alignmentFromDB)->get(wire, CDCAlignment::wireBwdX);
736 back[1] = (*m_alignmentFromDB)->get(wire, CDCAlignment::wireBwdY);
737 back[2] = (*m_alignmentFromDB)->get(wire, CDCAlignment::wireBwdZ);
738
739 fwrd[0] = (*m_alignmentFromDB)->get(wire, CDCAlignment::wireFwdX);
740 fwrd[1] = (*m_alignmentFromDB)->get(wire, CDCAlignment::wireFwdY);
741 fwrd[2] = (*m_alignmentFromDB)->get(wire, CDCAlignment::wireFwdZ);
742
743 for (int i = 0; i < np; ++i) {
744 // On top of the wire-end positions corrected for layer alignment, we apply possible
745 // fine corrections per wire
746 m_BWirPosAlign[iL][iC][i] += back[i];
747 m_FWirPosAlign[iL][iC][i] += fwrd[i];
748 }
749
750 // double baseTension = 0.;
751 double baseTension = M_PI * m_senseWireDensity * m_senseWireDiameter * m_senseWireDiameter / (8.* m_WireSagCoef[iL][iC]);
752 double tension = (*m_alignmentFromDB)->get(wire, CDCAlignment::wireTension);
753 // std::cout << back[0] <<" "<< back[1] <<" "<< back[2] <<" "<< fwrd[0] <<" "<< fwrd[1] <<" "<< fwrd[2] <<" "<< tension << std::endl;
754 m_WireSagCoefAlign[iL][iC] = M_PI * m_senseWireDensity *
755 m_senseWireDiameter * m_senseWireDiameter / (8.*(baseTension + tension));
756 // std::cout << "baseTension,tension= " << baseTension <<" "<< tension << std::endl;
757 } //end of layer loop
758 } //end of cell loop
759}
760
761
762// Set misalignment wire positions
763//TODO: merge this and setWirPosAlignParam() somehow
765{
766 const int np = 3;
767 double back[np], fwrd[np];
768
769 for (unsigned iL = 0; iL < c_maxNSenseLayers; ++iL) {
770
771 if (iL < m_firstLayerOffset) {
772 continue;
773 }
774
775 for (unsigned iC = 0; iC < m_nWires[iL]; ++iC) {
776 // std::cout << "iLiC= " << iL <<" "<< iC << std::endl;
777 WireID wire(iL, iC);
778 back[0] = (*m_misalignmentFromDB)->get(wire, CDCMisalignment::wireBwdX);
779 back[1] = (*m_misalignmentFromDB)->get(wire, CDCMisalignment::wireBwdY);
780 back[2] = (*m_misalignmentFromDB)->get(wire, CDCMisalignment::wireBwdZ);
781
782 fwrd[0] = (*m_misalignmentFromDB)->get(wire, CDCMisalignment::wireFwdX);
783 fwrd[1] = (*m_misalignmentFromDB)->get(wire, CDCMisalignment::wireFwdY);
784 fwrd[2] = (*m_misalignmentFromDB)->get(wire, CDCMisalignment::wireFwdZ);
785
786 for (int i = 0; i < np; ++i) {
787 m_BWirPosMisalign[iL][iC][i] = m_BWirPos[iL][iC][i] + back[i];
788 m_FWirPosMisalign[iL][iC][i] = m_FWirPos[iL][iC][i] + fwrd[i];
789 }
790
791 // double baseTension = 0.;
792 double baseTension = M_PI * m_senseWireDensity * m_senseWireDiameter * m_senseWireDiameter / (8.* m_WireSagCoef[iL][iC]);
793 double tension = (*m_misalignmentFromDB)->get(wire, CDCMisalignment::wireTension);
794 // std::cout << back[0] <<" "<< back[1] <<" "<< back[2] <<" "<< fwrd[0] <<" "<< fwrd[1] <<" "<< fwrd[2] <<" "<< tension << std::endl;
796 m_senseWireDiameter * m_senseWireDiameter / (8.*(baseTension + tension));
797 // std::cout << "baseTension,tension= " << baseTension <<" "<< tension << std::endl;
798 } //end of layer loop
799 } //end of cell loop
800}
801
802
803// Read x-t params.
804void CDCGeometryPar::readXT(const GearDir& gbxParams, const int mode)
805{
806 if (m_xtFileFormat == 0) {
807 // oldReadXT(gbxParams, mode);
808 } else {
809 newReadXT(gbxParams, mode);
810 }
811}
812
813
814// Read x-t params. (new)
815void CDCGeometryPar::newReadXT(const GearDir& gbxParams, const int mode)
816{
817 m_linearInterpolationOfXT = true; //must be true now
818
819 std::string fileName0 = CDCGeoControlPar::getInstance().getXtFile();
820 if (mode == 1) {
821 fileName0 = gbxParams.getString("xt4ReconFileName");
822 }
823
824 boost::iostreams::filtering_istream ifs;
825 openFileB(ifs, fileName0);
826 //TODO: use openFile() in cdc/utilities instead of the following 18 lines <- done
827 /*
828 std::string fileName1 = "/data/cdc/" + fileName0;
829 std::string fileName = FileSystem::findFile(fileName1);
830
831 if (fileName == "") {
832 fileName = FileSystem::findFile(fileName0);
833 }
834
835 if (fileName == "") {
836 B2FATAL("CDCGeometryPar: " << fileName1 << " not exist!");
837 } else {
838 B2INFO("CDCGeometryPar: open " << fileName1);
839 if ((fileName.rfind(".gz") != string::npos) && (fileName.length() - fileName.rfind(".gz") == 3)) {
840 ifs.push(boost::iostreams::gzip_decompressor());
841 }
842 ifs.push(boost::iostreams::file_source(fileName));
843 if (!ifs) B2FATAL("CDCGeometryPar: cannot open " << fileName1 << " !");
844 }
845 */
846
847 //read alpha bin info.
848 unsigned short nAlphaBins = 0;
849 if (ifs >> nAlphaBins) {
850 if (nAlphaBins == 0 || nAlphaBins > c_maxNAlphaPoints) B2FATAL("Fail to read alpha bins !");
851 } else {
852 B2FATAL("Fail to read alpha bins !");
853 }
854 m_nAlphaPoints = nAlphaBins;
855 double alpha0, alpha1, alpha2;
856 for (unsigned short i = 0; i < nAlphaBins; ++i) {
857 ifs >> alpha0 >> alpha1 >> alpha2;
858 m_alphaPoints[i] = alpha2;
859 }
860
861 //read theta bin info.
862 unsigned short nThetaBins = 0;
863 if (ifs >> nThetaBins) {
864 if (nThetaBins == 0 || nThetaBins > c_maxNThetaPoints) B2FATAL("CDCGeometryPar: fail to read theta bins !");
865 } else {
866 B2FATAL("CDCGeometryPar: fail to read theta bins !");
867 }
868 m_nThetaPoints = nThetaBins;
869 double theta0, theta1, theta2;
870
871 for (unsigned short i = 0; i < nThetaBins; ++i) {
872 ifs >> theta0 >> theta1 >> theta2;
873 m_thetaPoints[i] = theta2;
874 }
875
876 short np = 0;
877 unsigned short iCL, iLR;
878 const unsigned short npx = c_nXTParams - 1;
879 double xtc[npx];
880 double theta, alpha, dummy1;
881
882 ifs >> m_xtParamMode >> np;
883 if (m_xtParamMode < 0 || m_xtParamMode > 3) B2FATAL("CDCGeometryPar: invalid xt-parameterization mode read !");
884
885 if (np <= 0 || np > npx) B2FATAL("CDCGeometryPar: no. of xt-params. outside limits !");
886
887 const double epsi = 0.1;
888
889 while (ifs >> iCL) {
890
891 if (iCL < m_firstLayerOffset) {
892 continue;
893 }
894
895 ifs >> theta >> alpha >> dummy1 >> iLR;
896 for (int i = 0; i < np; ++i) {
897 ifs >> xtc[i];
898 }
899
900 int itheta = -99;
901 for (unsigned short i = 0; i < nThetaBins; ++i) {
902 if (fabs(theta - m_thetaPoints[i]) < epsi) {
903 itheta = i;
904 break;
905 }
906 }
907 if (itheta < 0) B2FATAL("CDCGeometryPar: thetas in xt.dat are inconsistent !");
908
909 int ialpha = -99;
910 for (unsigned short i = 0; i < nAlphaBins; ++i) {
911 if (fabs(alpha - m_alphaPoints[i]) < epsi) {
912 ialpha = i;
913 break;
914 }
915 }
916 if (ialpha < 0) B2FATAL("CDCGeometryPar: alphas in xt.dat are inconsistent !");
917
918 for (int i = 0; i < np; ++i) {
919 m_XT[iCL][iLR][ialpha][itheta][i] = xtc[i];
920 }
921
922 double boundT = xtc[6];
923 if (m_xtParamMode == 1) {
924 m_XT[iCL][iLR][ialpha][itheta][np] = ROOT::Math::Chebyshev5(boundT, xtc[0], xtc[1], xtc[2], xtc[3], xtc[4], xtc[5]);
925 } else {
926 m_XT[iCL][iLR][ialpha][itheta][np] =
927 xtc[0] + boundT
928 * (xtc[1] + boundT
929 * (xtc[2] + boundT
930 * (xtc[3] + boundT
931 * (xtc[4] + boundT
932 * (xtc[5])))));
933 }
934 } //end of while loop
935
936 // ifs.close();
937 boost::iostreams::close(ifs);
938
939 //convert unit
940 const double degrad = M_PI / 180.;
941 for (unsigned i = 0; i < nAlphaBins; ++i) {
942 m_alphaPoints[i] *= degrad;
943 }
944 for (unsigned i = 0; i < nThetaBins; ++i) {
945 m_thetaPoints[i] *= degrad;
946 }
947
948}
949
950
951// Read space resol. params.
952void CDCGeometryPar::readSigma(const GearDir& gbxParams, const int mode)
953{
954 if (m_sigmaFileFormat == 0) {
955 // oldReadSigma(gbxParams, mode);
956 } else {
957 newReadSigma(gbxParams, mode);
958 }
959}
960
961void CDCGeometryPar::newReadSigma(const GearDir& gbxParams, const int mode)
962{
963 m_linearInterpolationOfSgm = true; //must be true now
964
965 std::string fileName0 = CDCGeoControlPar::getInstance().getSigmaFile();
966 if (mode == 1) {
967 fileName0 = gbxParams.getString("sigma4ReconFileName");
968 }
969
970 ifstream ifs;
971 // openFile(ifs, fileName0);
972 openFileA(ifs, fileName0);
973
974 //read alpha bin info.
975 unsigned short nAlphaBins = 0;
976 if (ifs >> nAlphaBins) {
977 if (nAlphaBins == 0 || nAlphaBins > c_maxNAlphaPoints) B2FATAL("Fail to read alpha bins !");
978 } else {
979 B2FATAL("Fail to read alpha bins !");
980 }
981 m_nAlphaPoints4Sgm = nAlphaBins;
982 // std:: cout << nAlphaBins << std::endl;
983 double alpha0, alpha1, alpha2;
984 for (unsigned short i = 0; i < nAlphaBins; ++i) {
985 ifs >> alpha0 >> alpha1 >> alpha2;
986 m_alphaPoints4Sgm[i] = alpha2;
987 // std:: cout << alpha2 << std::endl;
988 }
989
990 //read theta bin info.
991 unsigned short nThetaBins = 0;
992 if (ifs >> nThetaBins) {
993 if (nThetaBins == 0 || nThetaBins > c_maxNThetaPoints) B2FATAL("CDCGeometryPar: fail to read theta bins !");
994 } else {
995 B2FATAL("CDCGeometryPar: fail to read theta bins !");
996 }
997 m_nThetaPoints4Sgm = nThetaBins;
998 // std:: cout << nThetaBins << std::endl;
999 double theta0, theta1, theta2;
1000
1001 for (unsigned short i = 0; i < nThetaBins; ++i) {
1002 ifs >> theta0 >> theta1 >> theta2;
1003 m_thetaPoints4Sgm[i] = theta2;
1004 // std:: cout << theta2 << std::endl;
1005 }
1006
1007 unsigned short np = 0;
1008 unsigned short iCL, iLR;
1009 double sigma[c_nSigmaParams]; // cppcheck-suppress constVariable
1010 double theta, alpha;
1011
1012 ifs >> m_sigmaParamMode >> np;
1013 // std:: cout << m_sigmaParamMode <<" "<< np << std::endl;
1014 if (m_sigmaParamMode < 0 || m_sigmaParamMode > 4) B2FATAL("CDCGeometryPar: invalid sigma-parameterization mode read !");
1015
1016 if (np > c_nSigmaParams) B2FATAL("CDCGeometryPar: no. of sigma-params. outside limits !");
1017
1018 ifs >> m_maxSpaceResol;
1019
1020 const double epsi = 0.1;
1021
1022 while (ifs >> iCL) {
1023
1024 if (iCL < m_firstLayerOffset) {
1025 continue;
1026 }
1027
1028 ifs >> theta >> alpha >> iLR;
1029 // std::cout << iCL <<" "<< theta <<" "<< alpha <<" "<< iLR << std::endl;
1030 for (int i = 0; i < np; ++i) {
1031 ifs >> sigma[i];
1032 }
1033
1034 int itheta = -99;
1035 for (unsigned short i = 0; i < nThetaBins; ++i) {
1036 if (fabs(theta - m_thetaPoints4Sgm[i]) < epsi) {
1037 itheta = i;
1038 break;
1039 }
1040 }
1041 if (itheta < 0) B2FATAL("CDCGeometryPar: thetas in sigma.dat are inconsistent !");
1042
1043 int ialpha = -99;
1044 for (unsigned short i = 0; i < nAlphaBins; ++i) {
1045 if (fabs(alpha - m_alphaPoints4Sgm[i]) < epsi) {
1046 ialpha = i;
1047 break;
1048 }
1049 }
1050 if (ialpha < 0) B2FATAL("CDCGeometryPar: alphas in sigma.dat are inconsistent !");
1051
1052 for (int i = 0; i < np; ++i) {
1053 m_Sigma[iCL][iLR][ialpha][itheta][i] = sigma[i];
1054 }
1055 } //end of while loop
1056
1057 ifs.close();
1058
1059 //convert unit
1060 const double degrad = M_PI / 180.;
1061 for (unsigned i = 0; i < nAlphaBins; ++i) {
1062 m_alphaPoints4Sgm[i] *= degrad;
1063 }
1064 for (unsigned i = 0; i < nThetaBins; ++i) {
1065 m_thetaPoints4Sgm[i] *= degrad;
1066 }
1067
1068 // std::cout << "end of newreadsigma " << std::endl;
1069}
1070
1071
1072// Read fudge factors
1073void CDCGeometryPar::readFFactor(const GearDir& gbxParams, const int mode)
1074{
1075 std::string fileName0 = CDCGeoControlPar::getInstance().getFFactorFile();
1076 if (mode == 1) {
1077 fileName0 = gbxParams.getString("fudgeFactorFileName");
1078 }
1079 B2WARNING("readFFactor is not ready! " << fileName0);
1080 //TODO; implement the following part.
1081}
1082
1083
1084// Read propagation speed param.
1085void CDCGeometryPar::readPropSpeed(const GearDir& gbxParams, const int mode)
1086{
1087 std::string fileName0 = CDCGeoControlPar::getInstance().getPropSpeedFile();
1088 if (mode == 1) {
1089 fileName0 = gbxParams.getString("propSpeed4ReconFileName");
1090 }
1091
1092 ifstream ifs;
1093 // openFile(ifs, fileName0);
1094 openFileA(ifs, fileName0);
1095
1096 uint iL;
1097 double speed;
1098 unsigned nRead = 0;
1099
1100 while (true) {
1101 ifs >> iL >> speed;
1102 if (ifs.eof()) break;
1103
1104 ++nRead;
1105
1106 m_propSpeedInv[iL] = (iL < m_firstLayerOffset) ? 0 : 1. / speed;
1107
1108 if (m_debug) B2DEBUG(150, iL << " " << speed);
1109 }
1110
1111 if (nRead != c_maxNSenseLayers) B2FATAL("CDCGeometryPar::readPropSpeed: #lines read-in (=" << nRead <<
1112 ") is inconsistent with total #layers (=" << c_maxNSenseLayers << ") !");
1113
1114 ifs.close();
1115}
1116
1117/*
1118// Read deltaz params.
1119void CDCGeometryPar::readDeltaz(const GearDir gbxParams)
1120{
1121 std::string fileName0 = gbxParams.getString("deltazFileName");
1122 fileName0 = "/data/cdc/" + fileName0;
1123 std::string fileName = FileSystem::findFile(fileName0);
1124
1125 ifstream ifs;
1126
1127 if (fileName == "") {
1128 B2FATAL("CDCGeometryPar: " << fileName0 << " not exist!");
1129 } else {
1130 B2INFO("CDCGeometryPar: " << fileName0 << " exists.");
1131 ifs.open(fileName.c_str());
1132 if (!ifs) B2FATAL("CDCGeometryPar: cannot open " << fileName0 << " !");
1133 }
1134
1135 int iL;
1136 unsigned nRead = 0;
1137
1138 while (ifs >> iL) {
1139 ifs >> m_bwdDz[iL] >> m_fwdDz[iL];
1140 ++nRead;
1141 if (m_debug) cout << iL << " " << m_bwdDz[iL] << " " << m_fwdDz[iL] << endl;
1142 }
1143
1144 if (nRead != c_maxNSenseLayers) B2FATAL("CDCGeometryPar::readDeltaz: #lines read-in (=" << nRead <<
1145 ") is inconsistent with total #layers (=" << c_maxNSenseLayers << ") !");
1146
1147 ifs.close();
1148}
1149*/
1150
1151
1152// Read t0 params.
1153void CDCGeometryPar::readT0(const GearDir& gbxParams, int mode)
1154{
1155 std::string fileName0 = CDCGeoControlPar::getInstance().getT0File();
1156 if (mode == 1) {
1157 fileName0 = gbxParams.getString("t04ReconFileName");
1158 }
1159
1160 ifstream ifs;
1161 // openFile(ifs, fileName0);
1162 openFileA(ifs, fileName0);
1163
1164 uint iL(0), iC(0);
1165 float t0(0);
1166 unsigned nRead = 0;
1167
1168 while (true) {
1169 ifs >> iL >> iC >> t0;
1170
1171 if (iL < m_firstLayerOffset) {
1172 continue;
1173 }
1174
1175 if (ifs.eof()) break;
1176
1177 ++nRead;
1178
1179 m_t0[iL][iC] = (iL < m_firstLayerOffset) ? 0. : t0;
1180
1181 if (m_debug) {
1182 B2DEBUG(150, iL << " " << iC << " " << t0);
1183 }
1184 }
1185
1186 if (nRead != m_nSenseWires) B2FATAL("CDCGeometryPar::readT0: #lines read-in (=" << nRead <<
1187 ") is inconsistent with total #sense wires (=" << m_nSenseWires << ") !");
1188
1189 ifs.close();
1190
1191 calcMeanT0();
1192}
1193
1194
1195/*
1196// Read bad-wires.
1197void CDCGeometryPar::readBadWire(const GearDir gbxParams, int mode)
1198{
1199 std::string fileName0 = CDCGeoControlPar::getInstance().getBwFile();
1200 if (mode == 1) {
1201 fileName0 = gbxParams.getString("bw4ReconFileName");
1202 }
1203
1204 ifstream ifs;
1205 // openFile(ifs, fileName0);
1206 openFileA(ifs, fileName0);
1207
1208 int iCL(0), iW(0);
1209 unsigned nRead = 0;
1210
1211 while (true) {
1212 ifs >> iCL >> iW;
1213
1214 if (ifs.eof()) break;
1215
1216 ++nRead;
1217
1218 m_badWire.push_back(WireID(iCL, iW));
1219
1220 if (m_debug) {
1221 B2DEBUG(150, iCL << " " << iW);
1222 }
1223 }
1224
1225 if (nRead > m_nSenseWires) B2FATAL("CDCGeometryPar::readBadWire: #lines read-in (=" << nRead <<
1226 ") is larger than the total #sense wires (=" << m_nSenseWires << ") !");
1227
1228 ifs.close();
1229}
1230*/
1231
1232
1233// Read time-walk parameters
1234void CDCGeometryPar::readTW(const GearDir& gbxParams, const int mode)
1235{
1236 std::string fileName0 = CDCGeoControlPar::getInstance().getTwFile();
1237 if (mode == 1) {
1238 fileName0 = gbxParams.getString("tw4ReconFileName");
1239 }
1240
1241 ifstream ifs;
1242 // openFile(ifs, fileName0);
1243 openFileA(ifs, fileName0);
1244
1245 unsigned short nPars(0);
1246 ifs >> m_twParamMode >> nPars;
1247 if (m_twParamMode > 1) {
1248 B2FATAL("CDCGeometryPar::readTW: invalid mode specified!");
1249 }
1250 if (nPars > 2) {
1251 B2FATAL("CDCGeometryPar::readTW: invalid #params specified!");
1252 }
1253
1254 unsigned iBoard = 0;
1255 unsigned nRead = 0;
1256 // Read board id and coefficients
1257 while (ifs >> iBoard) {
1258 for (unsigned short i = 0; i < nPars; ++i) {
1259 ifs >> m_timeWalkCoef[iBoard][i];
1260 }
1261 ++nRead;
1262 }
1263
1264 if (nRead != c_nBoards) B2FATAL("CDCGeometryPar::readTW: #lines read-in (=" << nRead << ") is inconsistent with #boards (=" <<
1265 c_nBoards
1266 << ") !");
1267
1268 ifs.close();
1269}
1270
1271
1272// Read ch-map
1273//void CDCGeometryPar::readChMap(const GearDir gbxParams)
1275{
1276 std::string fileName0 = CDCGeoControlPar::getInstance().getChMapFile();
1277
1278 ifstream ifs;
1279 // openFile(ifs, fileName0);
1280 openFileA(ifs, fileName0);
1281
1282 unsigned short iSL, iL, iW, iB, iC;
1283 unsigned nRead = 0;
1284
1285 while (true) {
1286 // Read a relation
1287 ifs >> iSL >> iL >> iW >> iB >> iC;
1288 if (ifs.eof()) break;
1289 if (iSL >= c_nSuperLayers or iSL < m_firstSuperLayerOffset) continue;
1290
1291 ++nRead;
1292 WireID wID(iSL, iL, iW);
1293 m_wireToBoard.insert(pair<WireID, unsigned short>(wID, iB));
1294 }
1295
1296 if (nRead != m_nSenseWires) B2FATAL("CDCGeometryPar::readChMap: #lines read-in (=" << nRead <<
1297 ") is inconsistent with #sense-wires (="
1298 << m_nSenseWires << ") !");
1299
1300 ifs.close();
1301}
1302
1303
1304// Read edep-to-adc
1305void CDCGeometryPar::readEDepToADC(const GearDir& gbxParams, const int mode)
1306{
1307 // B2WARNING("CDCGeometryPar: readEDepToADC is not ready!");
1308 std::string fileName0 = CDCGeoControlPar::getInstance().getEDepToADCFile();
1309 if (mode == 1) {
1310 fileName0 = gbxParams.getString("fudgeFactorFileName");
1311 }
1312
1313 ifstream ifs;
1314 // openFileA(ifs, fileName0);
1315 std::string fileName1 = "/data/cdc/" + fileName0;
1316 std::string fileName = FileSystem::findFile(fileName1, true);
1317
1318 if (fileName == "") {
1319 fileName = FileSystem::findFile(fileName0, true);
1320 }
1321
1322 if (fileName == "") {
1323 B2FATAL("CDC::openFile: " << fileName0 << " not exist!");
1324 } else {
1325 // B2INFO("CDC::openFile: open " << fileName);
1326 B2DEBUG(29, "CDC::openFile: open " << fileName);
1327 ifs.open(fileName.c_str());
1328 if (!ifs) B2FATAL("CDC::openFile: cannot open " << fileName << " !");
1329 }
1330
1331 unsigned short paramMode(0), nParams(0);
1332 ifs >> paramMode >> nParams;
1333 if (paramMode > 1) B2FATAL("Param mode > 1!");
1334 if (nParams > 7) B2FATAL("No. of params. > 7!");
1335 unsigned short groupId(0);
1336 ifs >> groupId;
1337 B2DEBUG(29, paramMode << " " << nParams << " " << groupId);
1338 if (groupId > 0) B2FATAL("GgroupId > 0!");
1339
1340
1341 unsigned short cLMin[c_nSuperLayers], cLMax[c_nSuperLayers]; //min and max clayer per super-layer
1342 cLMin[0] = 0;
1343 cLMax[0] = 7;
1344 for (unsigned int sl = 1; sl < c_nSuperLayers; ++sl) {
1345 cLMin[sl] = cLMax[0] + 6 * sl - 5;
1346 cLMax[sl] = cLMax[0] + 6 * sl;
1347 }
1348
1349 unsigned short id = 0;
1350 double coef = 0.;
1351 unsigned short nRead = 0;
1352 while (ifs >> id) {
1353 for (unsigned short i = 0; i < nParams; ++i) {
1354 ifs >> coef;
1355 for (unsigned short cL = cLMin[id]; cL <= cLMax[id]; ++cL) { //clayer loop
1356 for (unsigned short cell = 0; cell < m_nWires[cL]; ++cell) { //cell loop
1357 m_eDepToADCParams[cL][cell][i] = (cL < m_firstLayerOffset) ? 0 : coef;
1358 // B2DEBUG(29, "cL,cell,i,coef= "<< cL <<" "<< cell <<" "<< i <<" "<< coef);
1359 }
1360 }
1361 }
1362 ++nRead;
1363 if (nRead > c_nSuperLayers) B2FATAL("No. of read in lines > " << c_nSuperLayers << " !");
1364 }
1365
1366 ifs.close();
1367}
1368
1369
1370// Set t0 (from DB)
1372{
1373 for (unsigned short iCL = 0; iCL < c_maxNSenseLayers; ++iCL) {
1374 for (unsigned short iW = 0; iW < c_maxNDriftCells; ++iW) {
1375 m_t0[iCL][iW] = 0.;
1376 }
1377 }
1378
1379 for (auto const& ent : (*m_t0FromDB)->getT0s()) {
1380 const WireID wid = WireID(ent.first);
1381 const unsigned short iCL = wid.getICLayer();
1382 const unsigned short iW = wid.getIWire();
1383 m_t0[iCL][iW] = (iCL < m_firstLayerOffset) ? 0. : ent.second;
1384 }
1385
1386 calcMeanT0();
1387}
1388
1389
1390// Calculate mean t0
1391void CDCGeometryPar::calcMeanT0(double minT0, double maxT0, int maxIt, double nStdv, double epsi)
1392{
1393 double oldMeanT0 = 0;
1394 unsigned short it1 = 0;
1395 for (unsigned short it = 0; it < maxIt; ++it) {
1396 it1 = it;
1397 double effiSum = 0.;
1398 m_meanT0 = 0.;
1399 double stdvT0 = 0;
1400 for (unsigned short iCL = 0; iCL < c_maxNSenseLayers; ++iCL) {
1401 for (unsigned short iW = 0; iW < m_nWires[iCL]; ++iW) {
1402 if (m_t0[iCL][iW] < minT0 || m_t0[iCL][iW] > maxT0) continue;
1403 const WireID wid = WireID(iCL, iW);
1404 if (isHotWire(wid)) continue;
1405 if (isBadWire(wid)) continue;
1406 double effi = 1.;
1407 isDeadWire(wid, effi);
1408 effiSum += effi;
1409 m_meanT0 += (iCL < m_firstLayerOffset) ? 0. : effi * m_t0[iCL][iW];
1410 stdvT0 += (iCL < m_firstLayerOffset) ? 0. : effi * m_t0[iCL][iW] * m_t0[iCL][iW];
1411 }
1412 }
1413 if (effiSum > 0.) {
1414 m_meanT0 /= effiSum;
1415 stdvT0 /= effiSum;
1416 stdvT0 = sqrt(fabs(stdvT0 - m_meanT0 * m_meanT0));
1417 B2DEBUG(29, it << " " << effiSum << " " << m_meanT0 << " " << stdvT0);
1418 if (fabs(m_meanT0 - oldMeanT0) < epsi) break;
1419 oldMeanT0 = m_meanT0;
1420 minT0 = m_meanT0 - nStdv * stdvT0;
1421 maxT0 = m_meanT0 + nStdv * stdvT0;
1422 } else {
1423 B2FATAL("Wire efficiency sum <= 0!");
1424 }
1425 }
1426 if (it1 == maxIt - 1) B2WARNING("Max. iterations(=" << maxIt << ") needed to calculate the mean t0. Strange.");
1427}
1428
1429
1430// Set bad-wire (from DB)
1432{
1433 // m_badWire = (*m_badWireFromDB)->getWires();
1434 calcMeanT0();
1435}
1436
1437// Set bad-boards
1439{
1440 calcMeanT0();
1441}
1442
1443// Set prop.-speed (from DB)
1445{
1446 for (unsigned short iCL = 0; iCL < (*m_propSpeedFromDB)->getEntries(); ++iCL) {
1447 m_propSpeedInv[iCL] = (iCL < m_firstLayerOffset) ? 0. : 1. / (*m_propSpeedFromDB)->getSpeed(iCL);
1448 }
1449}
1450
1451
1452// Set time-walk coefficient (from DB)
1454{
1455 // (*m_timeWalkFromDB)->dump();
1456 m_twParamMode = (*m_timeWalkFromDB)->getTwParamMode();
1457
1458 for (unsigned short iBd = 0; iBd < (*m_timeWalkFromDB)->getEntries(); ++iBd) {
1459 int np = ((*m_timeWalkFromDB)->getTimeWalkParams(iBd)).size();
1460 for (int i = 0; i < np; ++i) {
1461 m_timeWalkCoef[iBd][i] = ((*m_timeWalkFromDB)->getTimeWalkParams(iBd))[i];
1462 }
1463 }
1464}
1465
1466
1467// Set xt params. (from DB)
1469{
1470 m_linearInterpolationOfXT = true; //must be true now
1471
1472 // std::cout <<"setXtRelation called" << std::endl;
1473 m_nAlphaPoints = (*m_xtRelFromDB)->getNoOfAlphaBins();
1474 for (unsigned short i = 0; i < m_nAlphaPoints; ++i) {
1475 m_alphaPoints[i] = (*m_xtRelFromDB)->getAlphaPoint(i);
1476 // std::cout << m_alphaPoints[i]*180./M_PI << std::endl;
1477 }
1478
1479 m_nThetaPoints = (*m_xtRelFromDB)->getNoOfThetaBins();
1480 for (unsigned short i = 0; i < m_nThetaPoints; ++i) {
1481 m_thetaPoints[i] = (*m_xtRelFromDB)->getThetaPoint(i);
1482 // std::cout << m_thetaPoints[i]*180./M_PI << std::endl;
1483 }
1484
1485 m_xtParamMode = (*m_xtRelFromDB)->getXtParamMode();
1486
1487 for (unsigned short iCL = 0; iCL < c_maxNSenseLayers; ++iCL) {
1488 if (iCL < m_firstLayerOffset) {
1489 // m_XT is initialized to 0, but reading
1490 // (*m_xtRelFromDB)->getXtParams(iCL, iLR, iA, iT)
1491 // could fail if iCL < m_firstLayerOffset, thus continue as m_XT would be set to 0 in this case anyway
1492 continue;
1493 }
1494 for (unsigned short iLR = 0; iLR < 2; ++iLR) {
1495 for (unsigned short iA = 0; iA < m_nAlphaPoints; ++iA) {
1496 for (unsigned short iT = 0; iT < m_nThetaPoints; ++iT) {
1497 const std::vector<float> params = (*m_xtRelFromDB)->getXtParams(iCL, iLR, iA, iT);
1498 unsigned short np = params.size();
1499 // std::cout <<"np4xt= " << np << std::endl;
1500 for (unsigned short i = 0; i < np; ++i) {
1501 m_XT[iCL][iLR][iA][iT][i] = params[i];
1502 }
1503
1504 double boundT = m_XT[iCL][iLR][iA][iT][6];
1505 if (m_xtParamMode == 1) {
1506 m_XT[iCL][iLR][iA][iT][np] = ROOT::Math::Chebyshev5(boundT, m_XT[iCL][iLR][iA][iT][0], m_XT[iCL][iLR][iA][iT][1],
1507 m_XT[iCL][iLR][iA][iT][2], m_XT[iCL][iLR][iA][iT][3], m_XT[iCL][iLR][iA][iT][4], m_XT[iCL][iLR][iA][iT][5]);
1508 } else {
1509 m_XT[iCL][iLR][iA][iT][np] =
1510 m_XT[iCL][iLR][iA][iT][0] + boundT
1511 * (m_XT[iCL][iLR][iA][iT][1] + boundT
1512 * (m_XT[iCL][iLR][iA][iT][2] + boundT
1513 * (m_XT[iCL][iLR][iA][iT][3] + boundT
1514 * (m_XT[iCL][iLR][iA][iT][4] + boundT
1515 * (m_XT[iCL][iLR][iA][iT][5])))));
1516 }
1517 }
1518 }
1519 }
1520 }
1521
1522}
1523
1524
1525// Set sigma params. (from DB)
1527{
1528 m_linearInterpolationOfSgm = true; //must be true now
1529
1530 // std::cout <<"setSResol called" << std::endl;
1531 m_nAlphaPoints4Sgm = (*m_sResolFromDB)->getNoOfAlphaBins();
1532 for (unsigned short i = 0; i < m_nAlphaPoints4Sgm; ++i) {
1533 m_alphaPoints4Sgm[i] = (*m_sResolFromDB)->getAlphaPoint(i);
1534 // std::cout << m_alphaPoints4Sgm[i]*180./M_PI << std::endl;
1535 }
1536
1537 m_nThetaPoints4Sgm = (*m_sResolFromDB)->getNoOfThetaBins();
1538 for (unsigned short i = 0; i < m_nThetaPoints4Sgm; ++i) {
1539 m_thetaPoints4Sgm[i] = (*m_sResolFromDB)->getThetaPoint(i);
1540 // std::cout << m_thetaPoints4Sgm[i]*180./M_PI << std::endl;
1541 }
1542
1543 // std::cout << "m_nAlphaPoints4Sgm= " << m_nAlphaPoints4Sgm << std::endl;
1544 // std::cout << "m_nThetaPoints4Sgm= " << m_nThetaPoints4Sgm << std::endl;
1545
1546 m_sigmaParamMode = (*m_sResolFromDB)->getSigmaParamMode();
1547
1548 m_maxSpaceResol = (*m_sResolFromDB)->getMaxSpaceResol();
1549
1550 for (unsigned short iCL = 0; iCL < c_maxNSenseLayers; ++iCL) {
1551 for (unsigned short iLR = 0; iLR < 2; ++iLR) {
1552 for (unsigned short iA = 0; iA < m_nAlphaPoints4Sgm; ++iA) {
1553 for (unsigned short iT = 0; iT < m_nThetaPoints4Sgm; ++iT) {
1554 const std::vector<float> params = (*m_sResolFromDB)->getSigmaParams(iCL, iLR, iA, iT);
1555 unsigned short np = params.size();
1556 // std::cout <<"np4sigma= " << np << std::endl;
1557 for (unsigned short i = 0; i < np; ++i) {
1558 m_Sigma[iCL][iLR][iA][iT][i] = (iCL < m_firstLayerOffset) ? 0. : params[i];
1559 }
1560 }
1561 }
1562 }
1563 }
1564
1565}
1566
1567
1568// Set fudge factors (from DB)
1570{
1571 unsigned short groupId = (*m_fFactorFromDB)->getGroupID();
1572 unsigned short nEnt = (*m_fFactorFromDB)->getEntries();
1573 B2DEBUG(29, "setFFactor called: groupId,nEnt= " << groupId << " " << nEnt);
1574
1575 if (groupId == 0) { //per all-layers mode
1576 } else {
1577 B2FATAL("CDCGeometryPar:: Invalid group-id " << groupId << " specified!");
1578 }
1579
1580 for (unsigned short id = 0; id < nEnt; ++id) {
1581 unsigned short np = ((*m_fFactorFromDB)->getFactors(id)).size();
1582 if (np != 3) B2FATAL("CDCGeometryPar:: No. of fudge factors != 3!");
1583 for (unsigned short i = 0; i < np; ++i) {
1584 m_fudgeFactorForSigma[i] = ((*m_fFactorFromDB)->getFactors(id))[i];
1585 B2DEBUG(29, i << " " << m_fudgeFactorForSigma[i]);
1586 }
1587 }
1588
1592 B2DEBUG(29, "fudge factors= " << m_fudgeFactorForSigma[0] << " " << m_fudgeFactorForSigma[1] << " " << m_fudgeFactorForSigma[2]);
1593}
1594
1595
1596// Set ch-map (from DB)
1598{
1599 for (const auto& cm : (*m_chMapFromDB)) {
1600 const unsigned short isl = cm.getISuperLayer();
1601 if (isl >= c_nSuperLayers or isl < m_firstSuperLayerOffset) continue;
1602 const uint il = cm.getILayer();
1603 const int iw = cm.getIWire();
1604 const int iBd = cm.getBoardID();
1605 const WireID wID(isl, il, iw);
1606 m_wireToBoard.insert(pair<WireID, unsigned short>(wID, iBd));
1607 const int iCh = cm.getBoardChannel();
1608 m_wireToChannel.insert(pair<WireID, unsigned short>(wID, iCh));
1609 m_boardAndChannelToWire[iBd][iCh] = wID.getEWire();
1610 }
1611}
1612
1613// Set edep-to-ADC conversion params. (from DB)
1615{
1616 unsigned short groupId = (*m_eDepToADCConversionsFromDB)->getGroupID();
1617 unsigned short nEnt = (*m_eDepToADCConversionsFromDB)->getEntries();
1618 if (groupId == 0) { //per super-layer mode
1619 if (nEnt > c_nSuperLayers) B2FATAL("CDCGeometryPar:: group-id " << groupId << " and #entries " << nEnt << " are inconsistent!");
1620 } else if (groupId == 1) { //per layer mode
1621 if (nEnt > c_maxNSenseLayers) B2FATAL("CDCGeometryPar:: group-id " << groupId << " and #entries " << nEnt << " are inconsistent!");
1622 } else {
1623 B2FATAL("CDCGeometryPar:: Invalid group-id " << groupId << " specified !");
1624 }
1625
1626 unsigned short cLMin[c_nSuperLayers], cLMax[c_nSuperLayers]; //min and max clayer per super-layer
1627 cLMin[0] = 0;
1628 cLMax[0] = 7;
1629 for (unsigned int sl = 1; sl < c_nSuperLayers; ++sl) {
1630 cLMin[sl] = cLMax[0] + 6 * sl - 5;
1631 cLMax[sl] = cLMax[0] + 6 * sl;
1632 }
1633
1634 for (unsigned short id = 0; id < nEnt; ++id) {
1635 unsigned short np = ((*m_eDepToADCConversionsFromDB)->getParams(id)).size();
1636 if (np > 7) B2FATAL("CDCGeometryPar:: No. of edep-to-ADC conversion params. > 7");
1637 if (groupId == 0) { //per super-layer; id=super-layer
1638 for (unsigned short cL = cLMin[id]; cL <= cLMax[id]; ++cL) { //clayer loop
1639 for (unsigned short cell = 0; cell < m_nWires[cL]; ++cell) { //cell loop
1640 for (unsigned short i = 0; i < np; ++i) {
1641 m_eDepToADCParams[cL][cell][i] = (cL < m_firstLayerOffset) ? 0. : ((*m_eDepToADCConversionsFromDB)->getParams(id))[i];
1642 }
1643 }
1644 }
1645 } else if (groupId == 1) { //per clayer; id=clayer
1646 for (unsigned short cell = 0; cell < m_nWires[id]; ++cell) { //cell loop
1647 for (unsigned short i = 0; i < np; ++i) {
1648 m_eDepToADCParams[id][cell][i] = (id < m_firstLayerOffset) ? 0. : ((*m_eDepToADCConversionsFromDB)->getParams(id))[i];
1649 }
1650 }
1651 } else if (groupId == 2) { //per wire
1652 //not ready
1653 B2FATAL("CDCGeometryPar::setEDepToADCConversions(): groupId=2 not ready!");
1654 }
1655 }
1656}
1657
1658
1659double CDCGeometryPar::getEDepToADCConvFactor(unsigned short iCL, unsigned short iW, double edep, double dx, double costh)
1660{
1661 // double convF = (100.0 / 3.2); //keV -> count
1662 //Model assumed here is from CLEO-c:
1663 //Igen = Imea * [1 + alf*Imea/cth] / [1 + gam*Imea/cth];
1664 //cth = |costh| + dlt;
1665 //Igen: original dE/dx; Imea: measured dE/dx with space-charge effect
1666 const double mainF = (iCL < m_firstLayerOffset) ? 0. : m_eDepToADCParams[iCL][iW][0];
1667 const double& alf = (iCL < m_firstLayerOffset) ? 0. : m_eDepToADCParams[iCL][iW][1];
1668 const double& gam = (iCL < m_firstLayerOffset) ? 0. : m_eDepToADCParams[iCL][iW][2];
1669 const double& dlt = (iCL < m_firstLayerOffset) ? 0. : m_eDepToADCParams[iCL][iW][3];
1670 const double& a = (iCL < m_firstLayerOffset) ? 0. : m_eDepToADCParams[iCL][iW][4];
1671 const double& b = (iCL < m_firstLayerOffset) ? 0. : m_eDepToADCParams[iCL][iW][5];
1672 const double cth = fabs(costh) + dlt;
1673 const double iGen = edep / dx; // keV/cm
1674 const double tmp = cth - gam * iGen;
1675 const double disc = tmp * tmp + 4.*alf * cth * iGen;
1676
1677 double iMea = 0.;
1678 if (alf == 0.) {
1679 iMea = cth * iGen / tmp;
1680 } else if (disc >= 0.) {
1681 iMea = (-tmp + sqrt(disc)) / (2.*alf);
1682 // if (alf < 0.) {
1683 // B2INFO("alf<0:CHECK");
1684 // B2INFO(-tmp + sqrt(disc)) / (2.*alf));
1685 // B2INFO(-tmp - sqrt(disc)) / (2.*alf));
1686 // }
1687 }
1688
1689 double convF = mainF;
1690 if (iMea > 0.) {
1691 convF = mainF * std::min(iMea / iGen, 1.);
1692 // if (iMea > iGen) B2DEBUG(29, "iMea > iGen " << iMea <<" "<< iGen);
1693 } else {
1694 //TODO: check the following issue more
1695 // B2WARNING("CDCGeometryPar: Measured dE/dx <= 0!");
1696 // B2DEBUG(29, "CDCGeometryPar: Measured dE/dx <= 0!");
1697 // B2DEBUG(29, "iGen,iMea= " << std::setw(15) << std::scientific << std::setprecision(8) << iGen <<" "<< iMea);
1698 // B2DEBUG(29, "dx,mainF,alf,gam,dlt,cth,tmp,disc= " << dx <<" "<< mainF <<" "<< alf <<" "<< gam <<" "<< dlt <<" "<<" "<< tmp <<" "<< disc);
1699 }
1700 convF *= 1. + a * (costh - b);
1701 return convF;
1702}
1703
1704
1706{}
1707
1708const B2Vector3D CDCGeometryPar::wireForwardPosition(uint layerID, int cellID, EWirePosition set) const
1709{
1710 // return early in case of empty layer, i.e. layerID < m_firstLayerOffset
1711 if (layerID < m_firstLayerOffset) {
1712 return B2Vector3D(0, 0, 0);
1713 }
1714
1715 // std::cout <<"cdcgeopar::fwdpos set= " << set << std::endl;
1716 B2Vector3D wPos(m_FWirPosAlign[layerID][cellID][0],
1717 m_FWirPosAlign[layerID][cellID][1],
1718 m_FWirPosAlign[layerID][cellID][2]);
1719
1720 if (set == c_Misaligned) {
1721 wPos.SetX(m_FWirPosMisalign[layerID][cellID][0]);
1722 wPos.SetY(m_FWirPosMisalign[layerID][cellID][1]);
1723 wPos.SetZ(m_FWirPosMisalign[layerID][cellID][2]);
1724 } else if (set == c_Base) {
1725 wPos.SetX(m_FWirPos [layerID][cellID][0]);
1726 wPos.SetY(m_FWirPos [layerID][cellID][1]);
1727 wPos.SetZ(m_FWirPos [layerID][cellID][2]);
1728 }
1729
1730 return wPos;
1731}
1732
1733const B2Vector3D CDCGeometryPar::wireForwardPosition(uint layerID, int cellID, double z, EWirePosition set) const
1734{
1735 // return early in case of empty layer, i.e. layerID < m_firstLayerOffset
1736 if (layerID < m_firstLayerOffset) {
1737 return B2Vector3D(0, 0, 0);
1738 }
1739
1740 double yb_sag = 0.;
1741 double yf_sag = 0.;
1742 getWireSagEffect(set, layerID, cellID, z, yb_sag, yf_sag);
1743
1744 B2Vector3D wPos(m_FWirPosAlign[layerID][cellID][0], yf_sag,
1745 m_FWirPosAlign[layerID][cellID][2]);
1746 if (set == c_Misaligned) {
1747 wPos.SetX(m_FWirPosMisalign[layerID][cellID][0]);
1748 wPos.SetZ(m_FWirPosMisalign[layerID][cellID][2]);
1749 } else if (set == c_Base) {
1750 wPos.SetX(m_FWirPos [layerID][cellID][0]);
1751 wPos.SetZ(m_FWirPos [layerID][cellID][2]);
1752 }
1753
1754 return wPos;
1755}
1756
1757const B2Vector3D CDCGeometryPar::wireBackwardPosition(uint layerID, int cellID, EWirePosition set) const
1758{
1759 // return early in case of empty layer, i.e. layerID < m_firstLayerOffset
1760 if (layerID < m_firstLayerOffset) {
1761 return B2Vector3D(0, 0, 0);
1762 }
1763
1764 B2Vector3D wPos(m_BWirPosAlign[layerID][cellID][0],
1765 m_BWirPosAlign[layerID][cellID][1],
1766 m_BWirPosAlign[layerID][cellID][2]);
1767
1768 if (set == c_Misaligned) {
1769 wPos.SetX(m_BWirPosMisalign[layerID][cellID][0]);
1770 wPos.SetY(m_BWirPosMisalign[layerID][cellID][1]);
1771 wPos.SetZ(m_BWirPosMisalign[layerID][cellID][2]);
1772 } else if (set == c_Base) {
1773 wPos.SetX(m_BWirPos [layerID][cellID][0]);
1774 wPos.SetY(m_BWirPos [layerID][cellID][1]);
1775 wPos.SetZ(m_BWirPos [layerID][cellID][2]);
1776 }
1777
1778 return wPos;
1779}
1780
1781const B2Vector3D CDCGeometryPar::wireBackwardPosition(uint layerID, int cellID, double z, EWirePosition set) const
1782{
1783 // return early in case of empty layer, i.e. layerID < m_firstLayerOffset
1784 if (layerID < m_firstLayerOffset) {
1785 return B2Vector3D(0, 0, 0);
1786 }
1787
1788 double yb_sag = 0.;
1789 double yf_sag = 0.;
1790 getWireSagEffect(set, layerID, cellID, z, yb_sag, yf_sag);
1791
1792 B2Vector3D wPos(m_BWirPosAlign[layerID][cellID][0], yb_sag,
1793 m_BWirPosAlign[layerID][cellID][2]);
1794 if (set == c_Misaligned) {
1795 wPos.SetX(m_BWirPosMisalign[layerID][cellID][0]);
1796 wPos.SetZ(m_BWirPosMisalign[layerID][cellID][2]);
1797 } else if (set == c_Base) {
1798 wPos.SetX(m_BWirPos [layerID][cellID][0]);
1799 wPos.SetZ(m_BWirPos [layerID][cellID][2]);
1800 }
1801
1802 return wPos;
1803}
1804
1805double CDCGeometryPar::getWireSagCoef(EWirePosition set, uint layerID, int cellID) const
1806{
1807 double coef = m_WireSagCoef[layerID][cellID];
1808 if (set == c_Misaligned) {
1809 coef = m_WireSagCoefMisalign[layerID][cellID];
1810 } else if (set == c_Aligned) {
1811 coef = m_WireSagCoefAlign [layerID][cellID];
1812 }
1813 return coef;
1814}
1815
1817{
1818 static double IRWL[c_maxNSenseLayers] = {0};
1819
1820 IRWL[0] = outerRadiusInnerWall();
1821 for (unsigned i = 1; i < nWireLayers(); i++)
1822 //IRWL[i] = (m_rSLayer[i - 1] + m_rSLayer[i]) / 2.;
1823 IRWL[i] = (i == m_firstLayerOffset) ? outerRadiusInnerWall() : m_rFLayer[i - 1];
1824
1825 return IRWL;
1826}
1827
1829{
1830 static double ORWL[c_maxNSenseLayers] = {0};
1831
1832 ORWL[nWireLayers() - 1] = innerRadiusOuterWall();
1833 for (unsigned i = 0; i < nWireLayers() - 1; i++)
1834 //ORWL[i] = (m_rSLayer[i] + m_rSLayer[i + 1]) / 2.;
1835 ORWL[i] = m_rFLayer[i];
1836
1837 return ORWL;
1838}
1839
1840unsigned CDCGeometryPar::cellId(unsigned layerId, const B2Vector3D& position) const
1841{
1842 if (layerId < m_firstLayerOffset) {
1843 return 0;
1844 }
1845
1846 const unsigned nWires = m_nWires[layerId];
1847
1848 double offset = m_offSet[layerId];
1849 //...Offset modification to be aligned to axial at z=0...
1850 const double phiSize = 2 * M_PI / double(nWires);
1851 /*{
1852 const double phiF = phiSize * offset
1853 + phiSize * 0.5 * double(m_nShifts[layerId]);
1854 const double phiB = phiSize * offset;
1855 const B2Vector3D f(m_rSLayer[layerId] * cos(phiF), m_rSLayer[layerId] * sin(phiF), m_zSForwardLayer[layerId]);
1856 const B2Vector3D b(m_rSLayer[layerId] * cos(phiB), m_rSLayer[layerId] * sin(phiB), m_zSBackwardLayer[layerId]);
1857
1858 const B2Vector3D v = f - b;
1859 const B2Vector3D u = v.Unit();
1860 const double beta = (0 - b.Z()) / u.Z();
1861 const B2Vector3D p = b + beta * u;
1862 double phi0 = - atan2(p.Y(), p.X());
1863 offset += phi0 / (2 * M_PI / double(nWires));
1864 }*/
1865
1866 unsigned j = 0;
1867 for (unsigned i = 0; i < 1; ++i) {
1868 const double phiF = phiSize * (double(i) + offset)
1869 + phiSize * 0.5 * double(m_nShifts[layerId]) + m_globalPhiRotation;
1870 const double phiB = phiSize * (double(i) + offset) + m_globalPhiRotation;
1871 const B2Vector3D f(m_rSLayer[layerId] * cos(phiF), m_rSLayer[layerId] * sin(phiF), m_zSForwardLayer[layerId]);
1872 const B2Vector3D b(m_rSLayer[layerId] * cos(phiB), m_rSLayer[layerId] * sin(phiB), m_zSBackwardLayer[layerId]);
1873 const B2Vector3D v = f - b;
1874 const B2Vector3D u = v.Unit();
1875 const double beta = (position.Z() - b.Z()) / u.Z();
1876 const B2Vector3D p = b + beta * u;
1877 double dPhi = std::atan2(position.Y(), position.X())
1878 - std::atan2(p.Y(), p.X())
1879 + phiSize / 2.;
1880 while (dPhi < 0) dPhi += (2. * M_PI);
1881 j = int(dPhi / phiSize);
1882 while (j >= nWires) j -= nWires;
1883 }
1884
1885 return j;
1886}
1887
1888void CDCGeometryPar::generateXML(const string& of)
1889{
1890 //...Open xml file...
1891 std::ofstream ofs(of.c_str(), std::ios::out);
1892 if (! ofs) {
1893 B2ERROR("CDCGeometryPar::read !!! can not open file : "
1894 << of);
1895 }
1896 ofs << "<?xml version=\"1.0\" encoding=\"UTF-8\"?>"
1897 << endl
1898 << "<Subdetector type=\"CDC\">"
1899 << endl
1900 << " <Name>CDC BelleII </Name>"
1901 << endl
1902 << " <Description>CDC geometry parameters</Description>"
1903 << endl
1904 << " <Version>0</Version>"
1905 << endl
1906 << " <GeoCreator>CDCBelleII</GeoCreator>"
1907 << endl
1908 << " <Content>"
1909 << endl
1910 << " <Rotation desc=\"Rotation of the whole cdc detector (should be the same as beampipe)\" unit=\"mrad\">0.0</Rotation>"
1911 << endl
1912 << " <OffsetZ desc=\"The offset of the whole cdc in z with respect to the IP (should be the same as beampipe)\" unit=\"mm\">0.0</OffsetZ>"
1913 << endl
1914 << " <Material>CDCGas</Material>"
1915 << endl
1916 << endl;
1917
1918 ofs << " <SLayers>" << endl;
1919
1920 for (int i = 0; i < m_nSLayer; i++) {
1921 ofs << " <SLayer id=\"" << i << "\">" << endl;
1922 ofs << " <Radius desc=\"Radius of wires in this layer\" unit=\"mm\">" << senseWireR(i) << "</Radius>" << endl;
1923 ofs << " <BackwardZ desc=\"z position of this wire layer at backward endplate\" unit=\"mm\">" << senseWireBZ(
1924 i) << "</BackwardZ>" << endl;
1925 ofs << " <ForwardZ desc=\"z position of this wire layer at forward endplate\" unit=\"mm\">" << senseWireFZ(
1926 i) << "</ForwardZ>" << endl;
1927// ofs << " <BackwardPhi desc=\"azimuth angle of the first wire in this layer at backward endplate\" unit=\"rad\">" << wireBackwardPosition(i).phi() << "</BackwardPhi>" << endl;
1928// ofs << " <ForwardPhi desc=\"azimuth angle of the first wire in this layer at forward endplate\" unit=\"rad\">" << wireForwardPosition(i).phi() << "</ForwardPhi>" << endl;
1929 ofs << " <NHoles desc=\"the number of holes in this layer, 2*(cell number)\">" << nWiresInLayer(
1930 i) * 2 << "</NHoles>" << endl;
1931 ofs << " <NShift desc=\"the shifted hole number of each wire in this layer\">" << nShifts(i) << "</NShift>" << endl;
1932 ofs << " <Offset desc=\"wire offset in phi direction at endplate\">" << m_offSet[i] << "</Offset>" << endl;
1933 ofs << " </SLayer>" << endl;
1934 }
1935
1936 ofs << " </SLayers>" << endl;
1937 ofs << " <FLayers>" << endl;
1938
1939 for (int i = 0; i < m_nFLayer; i++) {
1940 ofs << " <FLayer id=\"" << i << "\">" << endl;
1941 ofs << " <Radius desc=\"Radius of field wires in this layer\" unit=\"mm\">" << fieldWireR(i) << "</Radius>" << endl;
1942 ofs << " <BackwardZ desc=\"z position of this field wire layer at backward endplate\" unit=\"mm\">" << fieldWireBZ(
1943 i) << "</BackwardZ>" << endl;
1944 ofs << " <ForwardZ desc=\"z position of this field wire layer at forward endplate\" unit=\"mm\">" << fieldWireFZ(
1945 i) << "</ForwardZ>" << endl;
1946 ofs << " </FLayer>" << endl;
1947 }
1948
1949 ofs << " </FLayers>" << endl;
1950
1951 ofs << " <InnerWall name=\"InnerWall\">" << endl;
1952 ofs << " <InnerR desc=\"Inner radius\" unit=\"mm\">" << innerRadiusInnerWall() << "</InnerR>" << endl;
1953 ofs << " <OuterR desc=\"Outer radius\" unit=\"mm\">" << outerRadiusInnerWall() << "</OuterR>" << endl;
1954 ofs << " <BackwardZ desc=\"z position at backward endplate\" unit=\"mm\">" << m_zWall[0][0] << "</BackwardZ>" << endl;
1955 ofs << " <ForwardZ desc=\"z position at forward endplate\" unit=\"mm\">" << m_zWall[0][1] << "</ForwardZ>" << endl;
1956 ofs << " </InnerWall>" << endl;
1957
1958 ofs << " <OuterWall name=\"OuterWall\">" << endl;
1959 ofs << " <InnerR desc=\"Inner radius\" unit=\"mm\">" << innerRadiusOuterWall() << "</InnerR>" << endl;
1960 ofs << " <OuterR desc=\"Outer radius\" unit=\"mm\">" << outerRadiusOuterWall() << "</OuterR>" << endl;
1961 ofs << " <BackwardZ desc=\"z position at backward endplate\" unit=\"mm\">" << m_zWall[2][0] << "</BackwardZ>" << endl;
1962 ofs << " <ForwardZ desc=\"z position at forward endplate\" unit=\"mm\">" << m_zWall[2][1] << "</ForwardZ>" << endl;
1963 ofs << " </OuterWall>" << endl;
1964
1965 ofs << " </Content>" << endl
1966 << "</Subdetector>" << endl;
1967}
1968
1969void CDCGeometryPar::getWireSagEffect(const EWirePosition set, const unsigned layerID, const unsigned cellID, const double Z,
1970 double& Yb_sag, double& Yf_sag) const
1971{
1972 //Input
1973 // set : c_Base, c_Misaligned or c_Aligned
1974 // layerID: layer id (0 - 55);
1975 // cellID: cell id in the layer;
1976 // Z: Z-coord. (cm) at which sense wire sag is computed.
1977 //
1978 //Output Yb_sag: Y-corrd. (cm) of intersection of a tangent and the backward endplate.
1979 // Here the tangent is computed from the 1'st derivative of
1980 // a paraboric wire (due to gravity) defined at Z.
1981 // Yf_sag: ibid. but for forward.
1982 //
1983 //N.B.- Maybe replaced with a bit more accurate formula.
1984 // - The electrostatic force effect is not included.
1985
1986 // return early in case of empty layer, i.e. layerID < m_firstLayerOffset
1987 if (layerID < m_firstLayerOffset) {
1988 Yb_sag = 0.;
1989 Yf_sag = 0.;
1990 return;
1991 }
1992
1993 double Xb = 0.;
1994 double Xf = 0.;
1995 double Yb = 0.;
1996 double Yf = 0.;
1997 double Zb = 0.;
1998 double Zf = 0.;
1999 double Coef = 0.;
2000
2001 if (set == c_Aligned) {
2002 Coef = m_WireSagCoefAlign[layerID][cellID];
2003 Yb = m_BWirPosAlign[layerID][cellID][1];
2004 Yf = m_FWirPosAlign[layerID][cellID][1];
2005 if (Coef == 0.) {
2006 Yb_sag = Yb;
2007 Yf_sag = Yf;
2008 return;
2009 }
2010 Xb = m_BWirPosAlign[layerID][cellID][0];
2011 Xf = m_FWirPosAlign[layerID][cellID][0];
2012 Zb = m_BWirPosAlign[layerID][cellID][2];
2013 Zf = m_FWirPosAlign[layerID][cellID][2];
2014
2015 } else if (set == c_Misaligned) {
2016 Coef = m_WireSagCoefMisalign[layerID][cellID];
2017 Yb = m_BWirPosMisalign[layerID][cellID][1];
2018 Yf = m_FWirPosMisalign[layerID][cellID][1];
2019 if (Coef == 0.) {
2020 Yb_sag = Yb;
2021 Yf_sag = Yf;
2022 return;
2023 }
2024 Xb = m_BWirPosMisalign[layerID][cellID][0];
2025 Xf = m_FWirPosMisalign[layerID][cellID][0];
2026 Zb = m_BWirPosMisalign[layerID][cellID][2];
2027 Zf = m_FWirPosMisalign[layerID][cellID][2];
2028
2029 } else if (set == c_Base) {
2030 Coef = m_WireSagCoef[layerID][cellID];
2031 Yb = m_BWirPos[layerID][cellID][1];
2032 Yf = m_FWirPos[layerID][cellID][1];
2033 if (Coef == 0.) {
2034 Yb_sag = Yb;
2035 Yf_sag = Yf;
2036 return;
2037 }
2038 Xb = m_BWirPos[layerID][cellID][0];
2039 Xf = m_FWirPos[layerID][cellID][0];
2040 Zb = m_BWirPos[layerID][cellID][2];
2041 Zf = m_FWirPos[layerID][cellID][2];
2042
2043 } else {
2044 B2FATAL("CDCGeometryPar::getWireSagEffect: called with an invalid set: " << " " << set);
2045 }
2046
2047 const double dx = Xf - Xb;
2048 const double dy = Yf - Yb;
2049 const double dz = Zf - Zb;
2050
2051 const double Zfp = sqrt(dz * dz + dx * dx); // Wire length in z-x plane since Zbp==0
2052 const double Zp = (Z - Zb) * Zfp / dz;
2053
2054 const double Y_sag = (Coef * (Zp - Zfp) + dy / Zfp) * Zp + Yb;
2055 const double dydz = (Coef * (2.*Zp - Zfp) * Zfp + dy) / dz;
2056
2057 Yb_sag = Y_sag + dydz * (Zb - Z);
2058 Yf_sag = Y_sag + dydz * (Zf - Z);
2059
2060}
2061
2062void CDCGeometryPar::setDesignWirParam(const unsigned layerID, const unsigned cellID)
2063{
2064 const unsigned L = layerID;
2065 const unsigned C = cellID;
2066
2067 const double offset = m_offSet[L];
2068 //...Offset modification to be aligned to axial at z=0...
2069 const double phiSize = 2 * M_PI / double(m_nWires[L]);
2070
2071 const double phiF = phiSize * (double(C) + offset)
2072 + phiSize * 0.5 * double(m_nShifts[L]) + m_globalPhiRotation;
2073
2074 m_FWirPos[L][C][0] = (L < m_firstLayerOffset) ? 0. : m_rSLayer[L] * cos(phiF);
2075 m_FWirPos[L][C][1] = (L < m_firstLayerOffset) ? 0. : m_rSLayer[L] * sin(phiF);
2076 m_FWirPos[L][C][2] = (L < m_firstLayerOffset) ? 0. : m_zSForwardLayer[L];
2077
2078 const double phiB = phiSize * (double(C) + offset) + m_globalPhiRotation;
2079
2080 m_BWirPos[L][C][0] = (L < m_firstLayerOffset) ? 0. : m_rSLayer[L] * cos(phiB);
2081 m_BWirPos[L][C][1] = (L < m_firstLayerOffset) ? 0. : m_rSLayer[L] * sin(phiB);
2082 m_BWirPos[L][C][2] = (L < m_firstLayerOffset) ? 0. : m_zSBackwardLayer[L];
2083
2084 for (int i = 0; i < 3; ++i) {
2085 m_FWirPosMisalign[L][C][i] = (L < m_firstLayerOffset) ? 0. : m_FWirPos[L][C][i];
2086 m_BWirPosMisalign[L][C][i] = (L < m_firstLayerOffset) ? 0. : m_BWirPos[L][C][i];
2087 m_FWirPosAlign [L][C][i] = (L < m_firstLayerOffset) ? 0. : m_FWirPos[L][C][i];
2088 m_BWirPosAlign [L][C][i] = (L < m_firstLayerOffset) ? 0. : m_BWirPos[L][C][i];
2089 }
2090
2092 (8. * m_senseWireTension);
2093 // m_WireSagCoef [L][C] = 0.;
2094 m_WireSagCoefMisalign[L][C] = (L < m_firstLayerOffset) ? 0. : m_WireSagCoef[L][C];
2095 m_WireSagCoefAlign [L][C] = (L < m_firstLayerOffset) ? 0. : m_WireSagCoef[L][C];
2096
2097}
2098
2099void CDCGeometryPar::outputDesignWirParam(const unsigned layerID, const unsigned cellID) const
2100{
2101
2102 const unsigned L = layerID;
2103 const unsigned C = cellID;
2104
2105 static bool first = true;
2106 static ofstream ofs;
2107 if (first) {
2108 first = false;
2109 ofs.open("alignment.dat");
2110 }
2111
2112 ofs << L << " " << C;
2113
2114 ofs << setiosflags(ios::showpoint | ios::uppercase);
2115
2116 for (int i = 0; i < 3; ++i) ofs << " " << setw(15) << setprecision(8) << m_BWirPos[L][C][i];
2117
2118 for (int i = 0; i < 3; ++i) ofs << " " << setw(15) << setprecision(8) << m_FWirPos[L][C][i];
2119 ofs << setiosflags(ios::fixed);
2120 ofs << " " << setw(4) << setprecision(1) << m_senseWireTension;
2121
2122 ofs << endl;
2123}
2124
2125double CDCGeometryPar::getDriftV(const double time, const unsigned short iCLayer, const unsigned short lr, const double alpha,
2126 const double theta) const
2127{
2128 if (iCLayer < m_firstLayerOffset || iCLayer >= c_maxNSenseLayers) {
2129 return 0.;
2130 }
2131
2132 double dDdt = 0.;
2133
2134 //calculate min. drift time
2135 double minTime = getMinDriftTime(iCLayer, lr, alpha, theta);
2136 double delta = time - minTime;
2137
2138 //convert incoming- to outgoing-lr
2139 unsigned short lro = getOutgoingLR(lr, alpha);
2140
2142 B2FATAL("linearInterpolationOfXT = false is not allowed now !");
2143 } else {
2144 double wal(0.);
2145 unsigned short ial[2] = {0};
2146 unsigned short ilr[2] = {lro, lro};
2147 getClosestAlphaPoints(alpha, wal, ial, ilr);
2148 double wth(0.);
2149 unsigned short ith[2] = {0};
2150 getClosestThetaPoints(alpha, theta, wth, ith);
2151
2152 unsigned short jal(0), jlr(0), jth(0);
2153 double w = 0.;
2154
2155 //use xt reversed at (x=0,t=tmin) for delta<0 ("negative drifttime")
2156 double timep = delta < 0. ? minTime - delta : time;
2157
2158 //compute linear interpolation (=weithed average over 4 points) in (alpha-theta) space
2159 for (unsigned k = 0; k < 4; ++k) {
2160 if (k == 0) {
2161 jal = ial[0];
2162 jlr = ilr[0];
2163 jth = ith[0];
2164 w = (1. - wal) * (1. - wth);
2165 } else if (k == 1) {
2166 jal = ial[0];
2167 jlr = ilr[0];
2168 jth = ith[1];
2169 w = (1. - wal) * wth;
2170 } else if (k == 2) {
2171 jal = ial[1];
2172 jlr = ilr[1];
2173 jth = ith[0];
2174 w = wal * (1. - wth);
2175 } else if (k == 3) {
2176 jal = ial[1];
2177 jlr = ilr[1];
2178 jth = ith[1];
2179 w = wal * wth;
2180 }
2181
2182 double boundary = m_XT[iCLayer][jlr][jal][jth][6];
2183
2184 if (timep < boundary) {
2185 if (m_xtParamMode == 1) {
2186 const double& c1 = m_XT[iCLayer][jlr][jal][jth][1];
2187 const double& c2 = m_XT[iCLayer][jlr][jal][jth][2];
2188 const double& c3 = m_XT[iCLayer][jlr][jal][jth][3];
2189 const double& c4 = m_XT[iCLayer][jlr][jal][jth][4];
2190 const double& c5 = m_XT[iCLayer][jlr][jal][jth][5];
2191 dDdt += w * ROOT::Math::Chebyshev4(timep, c1 + 3.*c3 + 5.*c5, 4.*c2 + 8.*c4, 6.*c3 + 10.*c5, 8.*c4, 10.*c5);
2192 } else {
2193 dDdt += w * (m_XT[iCLayer][jlr][jal][jth][1] + timep
2194 * (2.*m_XT[iCLayer][jlr][jal][jth][2] + timep
2195 * (3.*m_XT[iCLayer][jlr][jal][jth][3] + timep
2196 * (4.*m_XT[iCLayer][jlr][jal][jth][4] + timep
2197 * (5.*m_XT[iCLayer][jlr][jal][jth][5])))));
2198 }
2199 } else {
2200 dDdt += w * m_XT[iCLayer][jlr][jal][jth][7];
2201 }
2202 } //end of weighted mean loop
2203 }
2204
2205 dDdt = fabs(dDdt);
2206 //n.b. following line not needed since dDdt > 0 even for delta < 0
2207 // if (delta < 0.) dDdt *= -1.;
2208 return dDdt;
2209
2210}
2211
2212//TODO: mearge with getDriftLength
2213double CDCGeometryPar::getDriftLength0(const double time, const unsigned short iCLayer, const unsigned short lr, const double alpha,
2214 const double theta) const
2215{
2216 if (iCLayer < m_firstLayerOffset || iCLayer >= c_maxNSenseLayers) {
2217 return 0.;
2218 }
2219
2220 double dist = 0.;
2221
2222 //convert incoming- to outgoing-lr
2223 unsigned short lro = getOutgoingLR(lr, alpha);
2224
2225 // std::cout << m_linearInterpolationOfXT << std::endl;
2226 // exit(-1);
2228 B2FATAL("linearInterpolationOfXT = false is not allowed now !");
2229 } else {
2230 double wal(0.);
2231 unsigned short ial[2] = {0};
2232 unsigned short ilr[2] = {lro, lro};
2233 getClosestAlphaPoints(alpha, wal, ial, ilr);
2234 double wth(0.);
2235 unsigned short ith[2] = {0};
2236 getClosestThetaPoints(alpha, theta, wth, ith);
2237
2238 unsigned short jal(0), jlr(0), jth(0);
2239 double w = 0.;
2240
2241 //use xt reversed at (x=0,t=tmin) for delta<0 ("negative drifttime")
2242 double timep = time;
2243 // std::cout << "iCLayer,alpha,theta,lro= " << iCLayer <<" "<< (180./M_PI)*alpha <<" "<< (180./M_PI)*theta <<" "<< lro << std::endl;
2244
2245 //compute linear interpolation (=weithed average over 4 points) in (alpha-theta) space
2246 for (unsigned k = 0; k < 4; ++k) {
2247 if (k == 0) {
2248 jal = ial[0];
2249 jlr = ilr[0];
2250 jth = ith[0];
2251 w = (1. - wal) * (1. - wth);
2252 } else if (k == 1) {
2253 jal = ial[0];
2254 jlr = ilr[0];
2255 jth = ith[1];
2256 w = (1. - wal) * wth;
2257 } else if (k == 2) {
2258 jal = ial[1];
2259 jlr = ilr[1];
2260 jth = ith[0];
2261 w = wal * (1. - wth);
2262 } else if (k == 3) {
2263 jal = ial[1];
2264 jlr = ilr[1];
2265 jth = ith[1];
2266 w = wal * wth;
2267 }
2268
2269 double boundary = m_XT[iCLayer][jlr][jal][jth][6];
2270
2271 if (timep < boundary) {
2272 if (m_xtParamMode == 1) {
2273 dist += w * ROOT::Math::Chebyshev5(timep, m_XT[iCLayer][jlr][jal][jth][0], m_XT[iCLayer][jlr][jal][jth][1],
2274 m_XT[iCLayer][jlr][jal][jth][2], m_XT[iCLayer][jlr][jal][jth][3], m_XT[iCLayer][jlr][jal][jth][4], m_XT[iCLayer][jlr][jal][jth][5]);
2275 } else {
2276 dist += w * (m_XT[iCLayer][jlr][jal][jth][0] + timep
2277 * (m_XT[iCLayer][jlr][jal][jth][1] + timep
2278 * (m_XT[iCLayer][jlr][jal][jth][2] + timep
2279 * (m_XT[iCLayer][jlr][jal][jth][3] + timep
2280 * (m_XT[iCLayer][jlr][jal][jth][4] + timep
2281 * (m_XT[iCLayer][jlr][jal][jth][5]))))));
2282 }
2283 } else {
2284 dist += w * (m_XT[iCLayer][jlr][jal][jth][7] * (timep - boundary) + m_XT[iCLayer][jlr][jal][jth][8]);
2285 }
2286 // std::cout <<"k,w,dist= " << k <<" "<< w <<" "<< dist << std::endl;
2287 } //end of weighted mean loop
2288 }
2289
2290 // dist = fabs(dist);
2291 return dist;
2292
2293}
2294
2295double CDCGeometryPar::getDriftLength(const double time, const unsigned short iCLayer, const unsigned short lr, const double alpha,
2296 const double theta,
2297 const bool calculateMinTime,
2298 const double inputMinTime) const
2299{
2300 if (iCLayer < m_firstLayerOffset || iCLayer >= c_maxNSenseLayers) {
2301 return 0.;
2302 }
2303
2304 double dist = 0.;
2305
2306 //calculate min. drift time
2307 double minTime = calculateMinTime ? getMinDriftTime(iCLayer, lr, alpha, theta) : inputMinTime;
2308 double delta = time - minTime;
2309
2310 //convert incoming- to outgoing-lr
2311 unsigned short lro = getOutgoingLR(lr, alpha);
2312
2313 // std::cout << m_linearInterpolationOfXT << std::endl;
2314 // exit(-1);
2316 B2FATAL("linearInterpolationOfXT = false is not allowed now !");
2317 } else {
2318 double wal(0.);
2319 unsigned short ial[2] = {0};
2320 unsigned short ilr[2] = {lro, lro};
2321 getClosestAlphaPoints(alpha, wal, ial, ilr);
2322 double wth(0.);
2323 unsigned short ith[2] = {0};
2324 getClosestThetaPoints(alpha, theta, wth, ith);
2325
2326 unsigned short jal(0), jlr(0), jth(0);
2327 double w = 0.;
2328
2329 //use xt reversed at (x=0,t=tmin) for delta<0 ("negative drifttime")
2330 double timep = delta < 0. ? minTime - delta : time;
2331
2332 // std::cout << "iCLayer,alpha,theta,lro= " << iCLayer <<" "<< (180./M_PI)*alpha <<" "<< (180./M_PI)*theta <<" "<< lro << std::endl;
2333
2334 //compute linear interpolation (=weithed average over 4 points) in (alpha-theta) space
2335 for (unsigned k = 0; k < 4; ++k) {
2336 if (k == 0) {
2337 jal = ial[0];
2338 jlr = ilr[0];
2339 jth = ith[0];
2340 w = (1. - wal) * (1. - wth);
2341 } else if (k == 1) {
2342 jal = ial[0];
2343 jlr = ilr[0];
2344 jth = ith[1];
2345 w = (1. - wal) * wth;
2346 } else if (k == 2) {
2347 jal = ial[1];
2348 jlr = ilr[1];
2349 jth = ith[0];
2350 w = wal * (1. - wth);
2351 } else if (k == 3) {
2352 jal = ial[1];
2353 jlr = ilr[1];
2354 jth = ith[1];
2355 w = wal * wth;
2356 }
2357
2358 // std::cout << "k,w= " << k <<" "<< w << std::endl;
2359 // std::cout << "ial[0],[1],jal,wal= " << ial[0] <<" "<< ial[1] <<" "<< jal <<" "<< wal << std::endl;
2360 // std::cout << "ith[0],[1],jth,wth= " << ith[0] <<" "<< ith[1] <<" "<< jth <<" "<< wth << std::endl;
2361
2362 /*
2363 std::cout <<"iCLayer= " << iCLayer << std::endl;
2364 std::cout <<"lr= " << lr << std::endl;
2365 std::cout <<"jal,jth= " << jal <<" "<< jth << std::endl;
2366 std::cout <<"wal,wth= " << wal <<" "<< wth << std::endl;
2367 for (int i=0; i<9; ++i) {
2368 std::cout <<"a= "<< i <<" "<< m_XT[iCLayer][lro][jal][jth][i] << std::endl;
2369 }
2370 */
2371 double boundary = m_XT[iCLayer][jlr][jal][jth][6];
2372
2373 if (timep < boundary) {
2374 if (m_xtParamMode == 1) {
2375 dist += w * ROOT::Math::Chebyshev5(timep, m_XT[iCLayer][jlr][jal][jth][0], m_XT[iCLayer][jlr][jal][jth][1],
2376 m_XT[iCLayer][jlr][jal][jth][2], m_XT[iCLayer][jlr][jal][jth][3], m_XT[iCLayer][jlr][jal][jth][4], m_XT[iCLayer][jlr][jal][jth][5]);
2377 } else {
2378 dist += w * (m_XT[iCLayer][jlr][jal][jth][0] + timep
2379 * (m_XT[iCLayer][jlr][jal][jth][1] + timep
2380 * (m_XT[iCLayer][jlr][jal][jth][2] + timep
2381 * (m_XT[iCLayer][jlr][jal][jth][3] + timep
2382 * (m_XT[iCLayer][jlr][jal][jth][4] + timep
2383 * (m_XT[iCLayer][jlr][jal][jth][5]))))));
2384 }
2385 } else {
2386 dist += w * (m_XT[iCLayer][jlr][jal][jth][7] * (timep - boundary) + m_XT[iCLayer][jlr][jal][jth][8]);
2387 }
2388 // std::cout <<"k,w,dist= " << k <<" "<< w <<" "<< dist << std::endl;
2389 } //end of weighted mean loop
2390 }
2391
2392 dist = fabs(dist);
2393 if (delta < 0.) dist *= -1.;
2394 return dist;
2395
2396}
2397
2398double CDCGeometryPar::getMinDriftTime(const unsigned short iCLayer, const unsigned short lr, const double alpha,
2399 const double theta) const
2400{
2401 if (iCLayer < m_firstLayerOffset || iCLayer >= c_maxNSenseLayers) {
2402 return 0.;
2403 }
2404
2405 double minTime = 0.;
2406
2407 //convert incoming- to outgoing-lr
2408 unsigned short lro = getOutgoingLR(lr, alpha);
2409
2411 B2FATAL("linearInterpolationOfXT = false is not allowed now !");
2412 } else {
2413 double wal(0.);
2414 unsigned short ial[2] = {0};
2415 unsigned short ilr[2] = {lro, lro};
2416 getClosestAlphaPoints(alpha, wal, ial, ilr);
2417 double wth(0.);
2418 unsigned short ith[2] = {0};
2419 getClosestThetaPoints(alpha, theta, wth, ith);
2420
2421 unsigned short jal(0), jlr(0), jth(0);
2422 double w = 0.;
2423
2424 double c[6] = {0.}, a[6] = {0.};
2425 for (unsigned k = 0; k < 4; ++k) {
2426 if (k == 0) {
2427 jal = ial[0];
2428 jlr = ilr[0];
2429 jth = ith[0];
2430 w = (1. - wal) * (1. - wth);
2431 } else if (k == 1) {
2432 jal = ial[0];
2433 jlr = ilr[0];
2434 jth = ith[1];
2435 w = (1. - wal) * wth;
2436 } else if (k == 2) {
2437 jal = ial[1];
2438 jlr = ilr[1];
2439 jth = ith[0];
2440 w = wal * (1. - wth);
2441 } else if (k == 3) {
2442 jal = ial[1];
2443 jlr = ilr[1];
2444 jth = ith[1];
2445 w = wal * wth;
2446 }
2447
2448 for (int i = 0; i < 5; ++i) {
2449 c[i] += w * m_XT[iCLayer][jlr][jal][jth][i];
2450 }
2451 }
2452
2453 if (m_xtParamMode == 1) { //convert c to coeff for normal-poly if Chebyshev
2454 a[0] = c[0] - c[2] + c[4];
2455 a[1] = c[1] - 3.*c[3] + 5.*c[5];
2456 a[2] = 2.*c[2] - 8.*c[4];
2457 a[3] = 4.*c[3] - 20.*c[5];
2458 a[4] = 8.*c[4];
2459 a[5] = 16.*c[5];
2460 } else { //normal-poly
2461 for (int i = 0; i < 5; ++i) a[i] = c[i];
2462 }
2463
2464 //estimate an initial value
2465 // bool check = false;
2466 // bool nointersection = false;
2467 if (a[2] != 0.) { //2nd-order approx. near t=0
2468 const double det = a[1] * a[1] - 4.*a[2] * a[0];
2469 if (det >= 0.) {
2470 //Choose the solution with dx/dt > 0 which gives x=0
2471 minTime = (-a[1] + sqrt(det)) / (2.*a[2]);
2472 } else {
2473 //Choose the solution with smallest x
2474 // nointersection = true;
2475 minTime = -a[1] / (2.*a[2]);
2476 // cout <<"smallest-x solution= " << minTime << endl;
2477 }
2478 } else if (a[1] != 0.) {
2479 minTime = -a[0] / a[1]; //1st-order approx.
2480 } else {
2481 B2WARNING("CDCGeometryPar::getMinDriftTime: minDriftTime not determined; assume zero.\n" << "layer(#0-55),lr,alpha(rad),theta= " <<
2482 iCLayer << " " << lr << " " << alpha << " " << theta);
2483 return minTime;
2484 }
2485
2486 // double minTime0 = minTime;
2487 //higher-order corr. using Newton method; trial to minimize x^2
2488 double edm; // = 10.; //(cm) (SG: fix to avoid cpp-check warning)
2489 // const double epsi4t = 0.01; //(ns)
2490 // const double epsi4x = 1.e-5; //(cm)
2491 const double epsi4x = 5.e-6; //(cm)
2492 // const unsigned short maxIter = 4;
2493 const unsigned short maxIter = 8;
2494 const double maxDt = 20.; //(ns)
2495 unsigned short nIter = 0;
2496 double minXsq = 1.e10; //(cm^2)
2497 double minMinTime = minTime;
2498 // double told = minTime + 1000.*epsi4t;
2499 // while (fabs(minTime - told) > epsi && nIter <= maxIter) {
2500 for (nIter = 0; nIter <= maxIter; ++nIter) {
2501 // told = minTime;
2502 double t = minTime;
2503 double x = a[0] + t * (a[1] + t * (a[2] + t * (a[3] + t * (a[4] + t * a[5]))));
2504 double x2 = x * x;
2505 if (x2 < minXsq) {
2506 minXsq = x2;
2507 minMinTime = t;
2508 }
2509 double xp = a[1] + t * (2 * a[2] + t * (3 * a[3] + t * (4 * a[4] + t * 5 * a[5])));
2510 double xpp = 2 * a[2] + t * (6 * a[3] + t * (12 * a[4] + t * 20 * a[5]));
2511 double den = xp * xp + x * xpp;
2512 if (den <= 0.) {
2513 den = xp * xp;
2514 }
2515
2516 if (den > 0.) {
2517 //estimated distance to min.
2518 edm = fabs(x * xp) / sqrt(den); //not in distance^2 but in distance
2519 if (edm < epsi4x) break; //converged
2520 }
2521
2522 double dt = 1.; //dt for den=0 (ns)
2523 if (den != 0.) {
2524 dt = x * xp / den;
2525 if (dt >= 0.) {
2526 dt = std::min(dt, maxDt);
2527 } else {
2528 dt = std::max(dt, -maxDt);
2529 }
2530 } else {
2531 B2WARNING("CDCGeometryPar::getMinDriftTime: den = 0\n" << "layer(#0-55),lr,alpha(rad),theta= " <<
2532 iCLayer << " "
2533 << lr <<
2534 " " << alpha << " " << theta);
2535 }
2536 minTime -= dt;
2537 } //end of iteration loop
2538
2539 //choose minMinTime for not-converged case
2540 if (nIter == (maxIter + 1)) minTime = minMinTime;
2541
2542 /*
2543 if (fabs(minTime) > 20.) {
2544 B2WARNING("CDCGeometryPar::getMinDriftTime: |minDriftTime| > 20ns. Ok ?\n" << "layer(#0-55),lr,alpha(rad),theta,minTime(ns)= " <<
2545 iCLayer << " "
2546 << lr <<
2547 " " << alpha << " " << theta << " " << minTime);
2548 }
2549 if (nointersection) {
2550 cout <<"final minTime= " << minTime << endl;
2551 cout <<"final minx = " << a[0] + a[1] * minTime + a[2] *pow(minTime,2) + a[3] *pow(minTime,3) + a[4] *pow(minTime,4) + a[5] *pow(minTime,5) << endl;
2552 }
2553 */
2554 /*
2555 if (check) {
2556 double dmin = 999.;
2557 double tmin = 999.;
2558 for (int i = -10000; i < 10000; ++i) {
2559 double ti = 0.01 * i;
2560 double dl = fabs(getDriftLength0(ti, iCLayer, lr, alpha, theta));
2561 if (dl < dmin) {
2562 dmin = dl;
2563 tmin = ti;
2564 }
2565 }
2566
2567 double smartd = getDriftLength0(minTime, iCLayer, lr, alpha, theta);
2568 if (check) {
2569 // if (fabs(smartd) > dmin && minTime < tmin && fabs(minTime - tmin) > 0.1) {
2570 B2WARNING("CDCGeometryPar::getMinDriftTime \n" << "layer(#0-55),lr,alpha(rad),theta= " <<
2571 iCLayer << " "
2572 << lr <<
2573 " " << alpha << " " << theta);
2574 B2INFO("det, minTime0= " << det << " " << minTime0);
2575 B2INFO("direct search n,tmin,dmin= " << nIter << " " << tmin << " " << dmin);
2576 B2INFO(" smart search n,tmin,dmin= " << nIter << " " << minTime << " " << getDriftLength0(minTime, iCLayer, lr, alpha, theta));
2577
2578 for (int i=-200; i < 200; ++i) {
2579 double ti = 0.25*i;
2580 double dl = getDriftLength0(ti, iCLayer, lr, alpha, theta);
2581 std::cout << ti <<" "<< dl << std::endl;
2582 }
2583 exit(-1);
2584 }
2585 }
2586 */
2587 }
2588
2589 return minTime;
2590}
2591
2592double CDCGeometryPar::getDriftTime(const double dist, const unsigned short iCLayer, const unsigned short lr, const double alpha,
2593 const double theta) const
2594{
2595 if (iCLayer < m_firstLayerOffset || iCLayer >= c_maxNSenseLayers) {
2596 return 0.;
2597 }
2598
2599 //to be replaced with a smarter algorithm...
2600
2601 const double eps = 2.5e-1;
2602 const double maxTrials = 100;
2603
2604 // int ialpha = getAlphaBin(alpha);
2605 // int itheta = getThetaBin(theta);
2606
2607 //convert incoming- to outgoing-lr
2608 // unsigned short lrp = getOutgoingLR(lr, alpha);
2609
2610 double maxTime = 2000.; //in ns (n.b. further reduction, 2->1us could be ok)
2611 // if (m_XT[iCLayer][lrp][ialpha][itheta][7] == 0.) {
2612 // maxTime = m_XT[iCLayer][lrp][ialpha][itheta][6];
2613 // }
2614
2615 double minTime = getMinDriftTime(iCLayer, lr, alpha, theta);
2616 double t0 = minTime;
2617 // std::cout << "minTime,x= " << t0 <<" "<< getDriftLength(t0, iCLayer, lr, alpha, theta) << std::endl;
2618 const bool calMinTime = false;
2619 // double d0 = getDriftLength(t0, iCLayer, lr, alpha, theta, calMinTime, minTime) - dist;
2620 double d0 = - dist;
2621
2622 unsigned i = 0;
2623 double t1 = maxTime;
2624 double time = dist * m_nominalDriftVInv;
2625 while (((t1 - t0) > eps) && (i < maxTrials)) {
2626 time = 0.5 * (t0 + t1);
2627 double d1 = getDriftLength(time, iCLayer, lr, alpha, theta, calMinTime, minTime) - dist;
2628 // std::cout <<"i,dist,t0,t1,d0,d1= " << i <<" "<< dist <<" "<< t0 <<" "<< t1 <<" "<< d0 <<" "<< d1 << std::endl;
2629 if (d0 * d1 > 0.) {
2630 t0 = time;
2631 } else {
2632 t1 = time;
2633 }
2634 ++i;
2635 }
2636
2637 if (i >= maxTrials - 1 || time > maxTime) {
2638 B2WARNING("CDCGeometryPar::getDriftTime " << dist << " " << iCLayer << " " << alpha << " " << lr << " " << t0 << " " << t1 << " " <<
2639 time << " " << d0);
2640 }
2641
2642 // std::cout <<"dist0,dist1= " << dist <<" "<< getDriftLength(time, iCLayer, lr, alpha, theta) <<" "<< getDriftLength(time, iCLayer, lr, alpha, theta) - dist << std::endl;
2643 // std::cout <<"dist0,dist1= " << dist <<" "<< getDriftLength(time, iCLayer, lr, 0., theta) <<" "<< getDriftLength(time, iCLayer, lr, 0., theta) - dist << std::endl;
2644 return time;
2645
2646}
2647
2648double CDCGeometryPar::getSigma(const double DriftL0, const unsigned short iCLayer, const unsigned short lr, const double alpha,
2649 const double theta) const
2650{
2651 if (iCLayer < m_firstLayerOffset || iCLayer >= c_maxNSenseLayers) {
2652 return 0.;
2653 }
2654
2655
2656 double sigma = 0.;
2657 //DriftL0 < 0 for the hit w/driftTime < 0; use |DriftL0| to avoid sigma=nan
2658 const double driftL = fabs(DriftL0);
2659
2660 //convert incoming- to outgoing-lr
2661 unsigned short lro = getOutgoingLR(lr, alpha);
2662
2664 B2FATAL("linearInterpolationOfSgm = false is not allowed now !");
2665 }
2667 double wal(0.);
2668 unsigned short ial[2] = {0};
2669 unsigned short ilr[2] = {lro, lro};
2670 getClosestAlphaPoints4Sgm(alpha, wal, ial, ilr);
2671 double wth(0.);
2672 unsigned short ith[2] = {0};
2673 getClosestThetaPoints4Sgm(alpha, theta, wth, ith);
2674
2675 //compute linear interpolation (=weithed average over 4 points) in (alpha-theta) space
2676 unsigned short jal(0), jlr(0), jth(0);
2677 double w = 0.;
2678 for (unsigned k = 0; k < 4; ++k) {
2679 if (k == 0) {
2680 jal = ial[0];
2681 jlr = ilr[0];
2682 jth = ith[0];
2683 w = (1. - wal) * (1. - wth);
2684 } else if (k == 1) {
2685 jal = ial[0];
2686 jlr = ilr[0];
2687 jth = ith[1];
2688 w = (1. - wal) * wth;
2689 } else if (k == 2) {
2690 jal = ial[1];
2691 jlr = ilr[1];
2692 jth = ith[0];
2693 w = wal * (1. - wth);
2694 } else if (k == 3) {
2695 jal = ial[1];
2696 jlr = ilr[1];
2697 jth = ith[1];
2698 w = wal * wth;
2699 }
2700 // std::cout << "k,w= " << k <<" "<< w << std::endl;
2701 // std::cout << "ial[0],[1],jal,wal= " << ial[0] <<" "<< ial[1] <<" "<< jal <<" "<< wal << std::endl;
2702 // std::cout << "ith[0],[1],jth,wth= " << ith[0] <<" "<< ith[1] <<" "<< jth <<" "<< wth << std::endl;
2703 /*
2704 std::cout <<"iCLayer= " << iCLayer << std::endl;
2705 std::cout <<"lr= " << lr << std::endl;
2706 std::cout <<"jal,jth= " << jal <<" "<< jth << std::endl;
2707 std::cout <<"wal,wth= " << wal <<" "<< wth << std::endl;
2708 for (int i=0; i<9; ++i) {
2709 std::cout <<"a= "<< i <<" "<< m_XT[iCLayer][lro][jal][jth][i] << std::endl;
2710 }
2711 */
2712 const double& P0 = m_Sigma[iCLayer][jlr][jal][jth][0];
2713 const double& P1 = m_Sigma[iCLayer][jlr][jal][jth][1];
2714 const double& P2 = m_Sigma[iCLayer][jlr][jal][jth][2];
2715 const double& P3 = m_Sigma[iCLayer][jlr][jal][jth][3];
2716 const double& P4 = m_Sigma[iCLayer][jlr][jal][jth][4];
2717 const double& P5 = m_Sigma[iCLayer][jlr][jal][jth][5];
2718 const double& P6 = m_Sigma[iCLayer][jlr][jal][jth][6];
2719
2720#if defined(CDC_DEBUG)
2721 cout << "driftL= " << driftL << endl;
2722 cout << "iCLayer= " << iCLayer << " " << jlr << " " << jal << " " << jth << endl;
2723 cout << "P0= " << P0 << endl;
2724 cout << "P1= " << P1 << endl;
2725 cout << "P2= " << P2 << endl;
2726 cout << "P3= " << P3 << endl;
2727 cout << "P4= " << P4 << endl;
2728 cout << "P5= " << P5 << endl;
2729 cout << "P6= " << P6 << endl;
2730#endif
2731 const double P7 = m_sigmaParamMode == 0 ? DBL_MAX : m_Sigma[iCLayer][jlr][jal][jth][7];
2732
2733 if (driftL < P7) {
2734 sigma += w * sqrt(P0 / (driftL * driftL + P1) + P2 * driftL + P3 +
2735 P4 * exp(P5 * (driftL - P6) * (driftL - P6)));
2736 } else {
2737 double forthTermAtP7 = P4 * exp(P5 * (P7 - P6) * (P7 - P6));
2738 const double& P8 = m_Sigma[iCLayer][jlr][jal][jth][8];
2739 if (m_sigmaParamMode == 1) {
2740 double sigmaAtP7 = sqrt(P0 / (P7 * P7 + P1) + P2 * P7 + P3 + forthTermAtP7);
2741 sigma += w * (P8 * (driftL - P7) + sigmaAtP7);
2742 } else if (m_sigmaParamMode == 2) {
2743 double onePls4AtP7 = sqrt(P0 / (P7 * P7 + P1) + forthTermAtP7);
2744 const double onePls4 = P8 * (driftL - P7) + onePls4AtP7;
2745 sigma += w * sqrt(P2 * driftL + P3 + onePls4 * onePls4);
2746 } else if (m_sigmaParamMode == 3) {
2747 forthTermAtP7 = sqrt(forthTermAtP7);
2748 const double forthTerm = P8 * (driftL - P7) + forthTermAtP7;
2749 sigma += w * sqrt(P0 / (driftL * driftL + P1) + P2 * driftL + P3 +
2750 forthTerm * forthTerm);
2751 } //end of mode
2752 } // end of driftL
2753 } //end of for loop
2754 }
2755
2756 sigma = std::min(sigma, m_maxSpaceResol);
2757 return sigma;
2758}
2759
2760unsigned short CDCGeometryPar::getOldLeftRight(const B2Vector3D& posOnWire, const B2Vector3D& posOnTrack,
2761 const B2Vector3D& momentum) const
2762{
2763 unsigned short lr = 0;
2764 double wCrossT = (posOnWire.Cross(posOnTrack)).Z();
2765
2766 if (wCrossT < 0.) {
2767 lr = 0;
2768 } else if (wCrossT > 0.) {
2769 lr = 1;
2770 } else {
2771 if ((posOnTrack - posOnWire).Perp() != 0.) {
2772 double wCrossP = (posOnWire.Cross(momentum)).Z();
2773 if (wCrossP > 0.) {
2774 if (posOnTrack.Perp() > posOnWire.Perp()) {
2775 lr = 0;
2776 } else {
2777 lr = 1;
2778 }
2779 } else if (wCrossP < 0.) {
2780 if (posOnTrack.Perp() < posOnWire.Perp()) {
2781 lr = 0;
2782 } else {
2783 lr = 1;
2784 }
2785 } else {
2786 lr = 0;
2787 }
2788 } else {
2789 lr = 0;
2790 }
2791 }
2792 return lr;
2793}
2794
2795unsigned short CDCGeometryPar::getNewLeftRightRaw(const B2Vector3D& posOnWire, const B2Vector3D& posOnTrack,
2796 const B2Vector3D& momentum) const
2797{
2798 const double distanceCrossP = ((posOnWire - posOnTrack).Cross(momentum)).Z();
2799 unsigned short int lr = (distanceCrossP > 0.) ? 1 : 0;
2800 return lr;
2801}
2802
2803//N.B. The following alpha and theta calculations are directly implemented in CDCRecoHit.cc tentatively to avoid a circular dependence betw cdc_dataobjects and cdclib. So be careful when changing the calculations !
2804double CDCGeometryPar::getAlpha(const B2Vector3D& posOnWire, const B2Vector3D& momentum) const
2805{
2806 const double wx = posOnWire.X();
2807 const double wy = posOnWire.Y();
2808 const double px = momentum.X();
2809 const double py = momentum.Y();
2810
2811 const double cross = wx * py - wy * px;
2812 const double dot = wx * px + wy * py;
2813
2814 return atan2(cross, dot);
2815}
2816
2817double CDCGeometryPar::getTheta(const B2Vector3D& momentum) const
2818{
2819 return atan2(momentum.Perp(), momentum.Z());
2820}
2821
2822
2823unsigned short CDCGeometryPar::getOutgoingLR(const unsigned short lr, const double alpha) const
2824{
2825 unsigned short lro = (fabs(alpha) <= 0.5 * M_PI) ? lr : abs(lr - 1);
2826 return lro;
2827}
2828
2829
2830double CDCGeometryPar::getOutgoingAlpha(const double alpha) const
2831{
2832 //convert incoming- to outgoing-alpha
2833 double alphao = alpha;
2834 if (alpha > 0.5 * M_PI) {
2835 alphao -= M_PI;
2836 } else if (alpha < -0.5 * M_PI) {
2837 alphao += M_PI;
2838 }
2839
2840 return alphao;
2841}
2842
2843double CDCGeometryPar::getOutgoingTheta(const double alpha, const double theta) const
2844{
2845 //convert incoming- to outgoing-theta
2846 double thetao = fabs(alpha) > 0.5 * M_PI ? M_PI - theta : theta;
2847 // std::cout << alpha <<" "<< thetao << std::endl;
2848 return thetao;
2849}
2850
2851void CDCGeometryPar::getClosestAlphaPoints(const double alpha, double& weight, unsigned short points[2],
2852 unsigned short lrs[2]) const
2853{
2854 double alphao = getOutgoingAlpha(alpha);
2855 weight = 1.;
2856
2857 if (alphao < m_alphaPoints[0]) {
2858 points[0] = m_nAlphaPoints - 1;
2859 points[1] = 0;
2860 if (m_nAlphaPoints > 1) {
2861 lrs[0] = abs(lrs[0] - 1); //flip lr
2862 weight = (alphao - (m_alphaPoints[points[0]] - M_PI)) / (m_alphaPoints[points[1]] - (m_alphaPoints[points[0]] - M_PI));
2863 }
2864 } else if (m_alphaPoints[m_nAlphaPoints - 1] <= alphao) {
2865 points[0] = m_nAlphaPoints - 1;
2866 points[1] = 0;
2867 if (m_nAlphaPoints > 1) {
2868 lrs[1] = abs(lrs[1] - 1); //flip lr
2869 weight = (alphao - m_alphaPoints[points[0]]) / (m_alphaPoints[points[1]] + M_PI - m_alphaPoints[points[0]]);
2870 }
2871 } else {
2872 const float* upper = std::upper_bound(m_alphaPoints,
2873 m_alphaPoints + m_nAlphaPoints, alphao);
2874 points[1] = upper - m_alphaPoints;
2875 points[0] = points[1] - 1;
2876 weight = (alphao - m_alphaPoints[points[0]]) / (m_alphaPoints[points[1]] - m_alphaPoints[points[0]]);
2877 }
2878}
2879
2880
2881void CDCGeometryPar::getClosestAlphaPoints4Sgm(const double alpha, double& weight, unsigned short points[2],
2882 unsigned short lrs[2]) const
2883{
2884 double alphao = getOutgoingAlpha(alpha);
2885 weight = 1.;
2886
2887 if (alphao < m_alphaPoints4Sgm[0]) {
2888 points[0] = m_nAlphaPoints4Sgm - 1;
2889 points[1] = 0;
2890 if (m_nAlphaPoints4Sgm > 1) {
2891 lrs[0] = abs(lrs[0] - 1); //flip lr
2892 weight = (alphao - (m_alphaPoints4Sgm[points[0]] - M_PI)) / (m_alphaPoints4Sgm[points[1]] - (m_alphaPoints4Sgm[points[0]] - M_PI));
2893 }
2894 } else if (m_alphaPoints4Sgm[m_nAlphaPoints4Sgm - 1] <= alphao) {
2895 points[0] = m_nAlphaPoints4Sgm - 1;
2896 points[1] = 0;
2897 if (m_nAlphaPoints4Sgm > 1) {
2898 lrs[1] = abs(lrs[1] - 1); //flip lr
2899 weight = (alphao - m_alphaPoints4Sgm[points[0]]) / (m_alphaPoints4Sgm[points[1]] + M_PI - m_alphaPoints4Sgm[points[0]]);
2900 }
2901 } else {
2902 const float* upper = std::upper_bound(m_alphaPoints4Sgm,
2904 points[1] = upper - m_alphaPoints4Sgm;
2905 points[0] = points[1] - 1;
2906 weight = (alphao - m_alphaPoints4Sgm[points[0]]) / (m_alphaPoints4Sgm[points[1]] - m_alphaPoints4Sgm[points[0]]);
2907 }
2908}
2909
2910
2911void CDCGeometryPar::getClosestThetaPoints(const double alpha, const double theta, double& weight, unsigned short points[2]) const
2912{
2913 const double thetao = getOutgoingTheta(alpha, theta);
2914
2915 if (thetao < m_thetaPoints[0]) {
2916 // points[0] = 0;
2917 // points[1] = 1;
2918 points[0] = 0;
2919 points[1] = 0;
2920 weight = 1.;
2921 } else if (m_thetaPoints[m_nThetaPoints - 1] <= thetao) {
2922 // points[0] = m_nThetaPoints - 2;
2923 // points[1] = m_nThetaPoints - 1;
2924 points[0] = m_nThetaPoints - 1;
2925 points[1] = m_nThetaPoints - 1;
2926 weight = 1.;
2927 } else {
2928 const float* upper = std::upper_bound(m_thetaPoints,
2929 m_thetaPoints + m_nThetaPoints, thetao);
2930 points[1] = upper - m_thetaPoints;
2931 points[0] = points[1] - 1;
2932 weight = (thetao - m_thetaPoints[points[0]]) / (m_thetaPoints[points[1]] - m_thetaPoints[points[0]]);
2933 }
2934}
2935
2936
2937void CDCGeometryPar::getClosestThetaPoints4Sgm(const double alpha, const double theta, double& weight,
2938 unsigned short points[2]) const
2939{
2940 const double thetao = getOutgoingTheta(alpha, theta);
2941
2942 if (thetao < m_thetaPoints4Sgm[0]) {
2943 points[0] = 0;
2944 points[1] = 0;
2945 weight = 1.;
2946 } else if (m_thetaPoints4Sgm[m_nThetaPoints4Sgm - 1] <= thetao) {
2947 points[0] = m_nThetaPoints4Sgm - 1;
2948 points[1] = m_nThetaPoints4Sgm - 1;
2949 weight = 1.;
2950 } else {
2951 const float* upper = std::upper_bound(m_thetaPoints4Sgm,
2953 points[1] = upper - m_thetaPoints4Sgm;
2954 points[0] = points[1] - 1;
2955 weight = (thetao - m_thetaPoints4Sgm[points[0]]) / (m_thetaPoints4Sgm[points[1]] - m_thetaPoints4Sgm[points[0]]);
2956 }
2957}
2958
2959
2961{
2962 // std::cout <<"setDisplacement called" << std::endl;
2963 for (const auto& disp : (*m_displacementFromDB)) {
2964 // const int iLayer0 = disp.getICLayer();
2965 // const int iWire0 = disp.getIWire();
2966 const int iLayer = WireID(disp.getEWire()).getICLayer();
2967 const int iWire = WireID(disp.getEWire()).getIWire();
2968 // if (iLayer0 != iLayer) B2FATAL("Layer0 != Layer");
2969 // if (iWire0 != iWire) B2FATAL("Wire0 != Wire");
2970 m_FWirPos[iLayer][iWire][0] += (iLayer < m_firstLayerOffset) ? 0. : disp.getXFwd();
2971 m_FWirPos[iLayer][iWire][1] += (iLayer < m_firstLayerOffset) ? 0. : disp.getYFwd();
2972 m_FWirPos[iLayer][iWire][2] += (iLayer < m_firstLayerOffset) ? 0. : disp.getZFwd();
2973 m_BWirPos[iLayer][iWire][0] += (iLayer < m_firstLayerOffset) ? 0. : disp.getXBwd();
2974 m_BWirPos[iLayer][iWire][1] += (iLayer < m_firstLayerOffset) ? 0. : disp.getYBwd();
2975 m_BWirPos[iLayer][iWire][2] += (iLayer < m_firstLayerOffset) ? 0. : disp.getZBwd();
2976 m_WireSagCoef[iLayer][iWire] = (iLayer < m_firstLayerOffset) ? 0. : M_PI * m_senseWireDensity * m_senseWireDiameter *
2977 m_senseWireDiameter / (8.*
2978 (m_senseWireTension + disp.getTension()));
2979 // std::cout <<"setdisp iL, iC, nominaltension, tension= " << iLayer <<" " << iWire <<" "<< m_senseWireTension <<" "<< disp.getTension() << std::endl;
2980 }
2981}
2982
2983
2985{
2986 const unsigned short nLayers[c_nSuperLayers] = {8, 6, 6, 6, 6, 6, 6, 6, 6}; //tentative
2987
2988 for (unsigned short SLayer = 0; SLayer < c_nSuperLayers; ++SLayer) {
2989 unsigned short firstCLayer = 0;
2990 for (unsigned short i = 0; i < SLayer; ++i) {
2991 firstCLayer += nLayers[i];
2992 }
2993 // std::cout <<"SLayer,firstCLayer= " << SLayer <<" "<< firstCLayer << std::endl;
2994
2995 B2Vector3D firstBPos = wireBackwardPosition(firstCLayer, 0);
2996 for (unsigned short Layer = 0; Layer < nLayers[SLayer]; ++Layer) {
2997 unsigned short CLayer = firstCLayer + Layer;
2998
2999 if (CLayer == firstCLayer) {
3000 m_shiftInSuperLayer[SLayer][Layer] = 0;
3001
3002 } else if (CLayer == firstCLayer + 1) {
3003 B2Vector3D BPos = wireBackwardPosition(CLayer, 0);
3004 m_shiftInSuperLayer[SLayer][Layer] = (BPos.Cross(firstBPos)).Z() > 0. ? -1 : 1;
3005 // std::cout <<"CLayer,Layer,shift= " << CLayer <<" "<< Layer <<" "<< m_shiftInSuperLayer[SLayer][Layer] <<" "<< (BPos.Cross(firstBPos)).Z() << std::endl;
3006
3007 } else {
3008 if (Layer % 2 == 0) {
3009 m_shiftInSuperLayer[SLayer][Layer] = 0;
3010 } else {
3011 m_shiftInSuperLayer[SLayer][Layer] = m_shiftInSuperLayer[SLayer][1];
3012 }
3013 }
3014 // std::cout <<"CLayer,Layer,shift= " << CLayer <<" "<< Layer <<" "<< m_shiftInSuperLayer[SLayer][Layer] << std::endl;
3015 }
3016 }
3017}
3018
3019signed short CDCGeometryPar::getShiftInSuperLayer(unsigned short iSuperLayer, unsigned short iLayer) const
3020{
3021 return m_shiftInSuperLayer[iSuperLayer][iLayer];
3022}
DataType Z() const
access variable Z (= .at(2) without boundary check)
Definition B2Vector3.h:435
void SetX(DataType x)
set X/1st-coordinate
Definition B2Vector3.h:457
B2Vector3< DataType > Cross(const B2Vector3< DataType > &p) const
Cross product.
Definition B2Vector3.h:296
DataType X() const
access variable X (= .at(0) without boundary check)
Definition B2Vector3.h:431
DataType Y() const
access variable Y (= .at(1) without boundary check)
Definition B2Vector3.h:433
void SetZ(DataType z)
set Z/3rd-coordinate
Definition B2Vector3.h:461
void SetY(DataType y)
set Y/2nd-coordinate
Definition B2Vector3.h:459
DataType Perp() const
The transverse component (R in cylindrical coordinate system).
Definition B2Vector3.h:200
B2Vector3< DataType > Unit() const
Unit vector parallel to this.
Definition B2Vector3.h:269
static const baseType layerDPhi
Layer rotation in global X-Y plane (gamma) dPhi = forward - backward endplate.
static const baseType layerDy
Layer shift in global Y dY = forward - backward endplate.
static const baseType wireBwdZ
Wire Z position w.r.t. nominal on backward endplate.
static const baseType layerDx
Layer shift in global X dX = forward - backward endplate.
static const baseType wireBwdY
Wire Y position w.r.t. nominal on backward endplate.
static const baseType wireFwdZ
Wire Z position w.r.t. nominal on forward endplate.
static const baseType wireFwdY
Wire Y position w.r.t. nominal on forward endplate.
static const baseType wireFwdX
Wire X position w.r.t. nominal on forward endplate.
static const baseType wireBwdX
Wire X position w.r.t. nominal on backward endplate.
static const baseType layerY
Layer shift in global Y at backward endplate.
static const baseType layerX
Layer shift in global X at backward endplate.
static const baseType layerPhi
Layer rotation in global X-Y plane (gamma) at backward endplate.
static const baseType wireTension
Wire tension w.r.t. nominal (=50. ?)
The Class for CDC geometry.
Definition CDCGeometry.h:27
static const baseType wireBwdZ
Wire Z position w.r.t. nominal on backward endplate.
static const baseType wireBwdY
Wire Y position w.r.t. nominal on backward endplate.
static const baseType wireFwdZ
Wire Z position w.r.t. nominal on forward endplate.
static const baseType wireFwdY
Wire Y position w.r.t. nominal on forward endplate.
static const baseType wireFwdX
Wire X position w.r.t. nominal on forward endplate.
static const baseType wireBwdX
Wire X position w.r.t. nominal on backward endplate.
static const baseType wireTension
Wire tension w.r.t. nominal (=50. ?)
The Class for CDC Geometry Control Parameters.
bool getDebug() const
Get debug flag.
bool getSigmaInputType()
Get input type for sigma.
bool getMisalignmentInputType()
Get input type for wire misalignment.
std::string getT0File() const
Get input file name for t0.
bool getDisplacementInputType()
Get input type for wire displacement.
double getAddFudgeFactorForSigmaForMC() const
Get additional fudge factor for space resol for MC.
std::string getEDepToADCFile() const
Get input file name for edeptoadc.
std::string getDisplacementFile() const
Get input file name for wire displacement.
std::string getMisalignmentFile() const
Get input file name for wire misalignment.
std::string getAlignmentFile() const
Get input file name for wire alignment.
bool getAlignmentInputType()
Get input type for wire alignment.
double getMaterialDefinitionMode() const
Get material definition mode.
std::string getPropSpeedFile() const
Get input file name for prop-speed.
bool getT0InputType()
Get input type for t0.
bool getEDepToADCInputType()
Get input type for edeptoadc.
std::string getSigmaFile() const
Get input file name for sigma.
bool getAlignment() const
Get alignment switch.
bool getMisalignment() const
Get misalignment switch.
bool getDisplacement() const
Get displacement switch.
int getSenseWireZposMode() const
Get sense wire z position mode.
double getAddFudgeFactorForSigmaForData() const
Get additional fudge factor for space resol for data.
std::string getXtFile() const
Get input file name for xt-relation.
bool getTwInputType()
Get input type for time-walk.
std::string getFFactorFile() const
Get input file name for fudge factor.
bool getChMapInputType()
Get input type for channel map.
std::string getTwFile() const
Get input file name for time-walk.
bool getFFactorInputType()
Get input type for fuge factor.
bool getXtInputType()
Get input type for xt.
bool getBwInputType()
Get input type for bad wire.
bool getPropSpeedInputType()
Get input type for prop.
static CDCGeoControlPar & getInstance()
Static method to get a reference to the CDCGeoControlPar instance.
std::string getChMapFile() const
Get input file name for channel map.
The Class for CDC Geometry Parameters.
void outputDesignWirParam(unsigned layerID, unsigned cellID) const
Write the designed wire parameters to the alignment.dat (default).
std::map< WireID, unsigned short > m_wireToChannel
map relating wire-id and channel-id.
void setWirPosAlignParams()
Set wire alignment params.
int m_materialDefinitionMode
Control switch for gas and wire material definition.
unsigned short m_nAlphaPoints4Sgm
No.
void readSigma(const GearDir &gbxParams, int mode=0)
Read spatial resolution table.
void readTW(const GearDir &gbxParams, int mode=0)
Read time-walk parameter.
void setXtRel()
Set XT-relation table (from DB) (new).
float m_BWirPosMisalign[c_maxNSenseLayers][c_maxNDriftCells][3]
Wire position incl.
DBObjPtr< CDCBadWires > * m_badWireFromDB
bad-wires retrieved from DB.
double getTheta(const B2Vector3D &momentum) const
Returns track incident angle (theta in rad.).
void readT0(const GearDir &gbxParams, int mode=0)
Read t0 parameters (from a file).
void setT0()
Set t0 parameters (from DB)
float m_alphaPoints4Sgm[c_maxNAlphaPoints]
alpha sampling points for sigma (rad)
virtual ~CDCGeometryPar()
Destructor.
ushort m_maxNSuperLayers
Maximum number of Super Layers.
double m_fudgeFactorForSigma[3]
Fuge factor for space resol.
double outerRadiusInnerWall() const
Returns the outer radius of the inner wall.
int m_sigmaParamMode
Mode for sigma parameterization.
bool m_alignment
Switch for alignment.
float m_alphaPoints[c_maxNAlphaPoints]
alpha sampling points for xt (rad)
void getClosestAlphaPoints(const double alpha, double &wal, unsigned short points[2], unsigned short lrs[2]) const
Returns the two closest alpha points for the input track incident angle (alpha).
void readEDepToADC(const GearDir &gbxParams, int mode=0)
Read spatial edep-to-adc conv.
double m_globalPhiRotation
Global ratation in phi (rad.); only for sense wires now.
EWirePosition
Wire position set.
double innerRadiusOuterWall() const
Returns the inner radius of the outer wall.
unsigned cellId(unsigned layerId, const B2Vector3D &position) const
The method to get cell id based on given layer id and the position.
void newReadSigma(const GearDir &gbxParams, int mode=0)
Read spatial resolution table in new format.
void setEDepToADCConversions()
Set edep-to-ADC conversion params.
double m_nominalPropSpeed
Nominal propagation speed of the sense wire (27.25 cm/nsec).
int nShifts(int layerId) const
Returns number shift.
float m_thetaPoints[c_maxNThetaPoints]
theta sampling points for xt (rad)
void setDesignWirParam(unsigned layerID, unsigned cellID)
Set the desizend wire parameters.
void getWireSagEffect(EWirePosition set, unsigned layerID, unsigned cellID, double zw, double &ywb_sag, double &ywf_sag) const
Compute effects of the sense wire sag.
DBArray< CDCDisplacement > * m_displacementFromDB
displacement params.
bool m_XTetc
Switch for reading x-t etc.
void setShiftInSuperLayer()
Calculates and saves shifts in super-layers (to be used in searching hits in neighboring cells)
int m_nShifts[c_maxNSenseLayers]
The array to store shifted cell number in each sense wire layer.
bool m_wireSag
Switch for sense wire sag.
bool m_XTetc4Recon
Switch for selecting xt etc.
double getDriftLength0(double dt, unsigned short layer, unsigned short lr, double alpha=0., double theta=0.5 *M_PI) const
Return the drift dength to the sense wire; tentative ver.
unsigned short m_boardAndChannelToWire[c_nBoards][48]
array relating board-channel-id and wire-id.
unsigned short m_nThetaPoints4Sgm
No.
float m_WireSagCoef[c_maxNSenseLayers][c_maxNDriftCells]
Wire sag coefficient for each cell; ibid.
void generateXML(const std::string &of)
Generate an xml file used in gearbox.
void getClosestAlphaPoints4Sgm(const double alpha, double &wal, unsigned short points[2], unsigned short lrs[2]) const
Returns the two closest alpha points for sigma for the input track incident angle (alpha).
DBObjPtr< HardwareClockSettings > m_clockSettings
hardware clock settings
float m_eDepToADCParams[c_maxNSenseLayers][c_maxNDriftCells][7]
edep-to-ADC conv.
double m_minTrackLength
Minimum track length for G4 step.
double getAlpha(const B2Vector3D &posOnWire, const B2Vector3D &momentum) const
Returns track incident angle in rphi plane (alpha in rad.).
double fieldWireR(int layerId) const
Returns radius of field wire in each layer.
float m_FWirPosMisalign[c_maxNSenseLayers][c_maxNDriftCells][3]
Wire position incl.
const B2Vector3D wireForwardPosition(uint layerId, int cellId, EWirePosition set=c_Base) const
Returns the forward position of the input sense wire.
bool isDeadWire(const WireID &wid, double &eff)
Inquire if the wire is dead.
unsigned m_nWires[c_maxNSenseLayers]
The array to store the wire number in each sense wire layre.
double getDriftV(double dt, unsigned short layer, unsigned short lr, double alpha=0., double theta=0.5 *M_PI) const
Get the realistic drift velocity.
double m_clockFreq4TDC
Clock frequency used for TDC (GHz).
int m_nFLayer
The number of field wire layer.
static CDCGeometryPar * m_B4CDCGeometryParDB
Pointer that saves the instance of this class.
DBObjPtr< CDCMisalignment > * m_misalignmentFromDB
misalignment params.
float m_XT[c_maxNSenseLayers][2][c_maxNAlphaPoints][c_maxNThetaPoints][c_nXTParams]
XT-relation coefficients for each layer, Left/Right, entrance angle and polar angle.
signed short m_shiftInSuperLayer[c_nSuperLayers][8]
shift in phi-direction wrt the 1st layer in each super layer
std::string m_version
The version of geometry parameters.
const B2Vector3D wireBackwardPosition(uint layerId, int cellId, EWirePosition set=c_Base) const
Returns the backward position of the input sense wire.
bool m_linearInterpolationOfXT
Switch for linear interpolation of xt.
DBObjPtr< CDCSpaceResols > * m_sResolFromDB
sigma params.
void setDisplacement()
Set displacement of sense wire.
double getOutgoingAlpha(const double alpha) const
Converts incoming- to outgoing-alpha.
bool isHotWire(const WireID &wid)
Inquire if the wire is hot.
double m_meanT0
mean t0 over all wires; should be double.
bool m_debug
Switch for debug printing.
DBObjPtr< CDCTimeWalks > * m_timeWalkFromDB
time-walk coeffs.
ushort m_nSenseWires
Maximum number of Sense Wires.
unsigned short m_tdcOffset
Not used; to be removed later.
double getSigma(double dist, unsigned short layer, unsigned short lr, double alpha=0., double theta=0.5 *M_PI) const
Return the basic resolution of drift length (cm).
CDCGeometryPar(const CDCGeometry *=nullptr)
Singleton class.
unsigned short getOutgoingLR(const unsigned short lr, const double alpha) const
Converts incoming-lr to outgoing-lr.
ushort m_maxNSenseLayers
Maximum number of Sense Wire Layers.
void readWirePositionParams(EWirePosition set, const CDCGeometry *geom)
Read displacement or (mis)alignment params from text file.
float m_thetaPoints4Sgm[c_maxNThetaPoints]
theta sampling points for sigma (rad)
std::map< WireID, unsigned short > m_wireToBoard
map relating wire-id and board-id.
double getOutgoingTheta(const double alpha, const double theta) const
Converts incoming- to outgoing-theta.
double m_senseWireDiameter
The diameter of sense wires.
double offset(int layerID) const
Return wire offset in phi direction at endplate.
bool m_displacement
Switch for displacement.
int m_twParamMode
Mode for tw parameterization.
double senseWireBZ(int layerId) const
Returns backward z position of sense wire in each layer.
double m_dzSBackwardLayer[c_maxNSenseLayers]
Corrections for backward z position of sense wire layers.
double getMinDriftTime(unsigned short layer, unsigned short lr, double alpha=0., double theta=0.5 *M_PI) const
Return the min.
DBObjPtr< CDCBadBoards > * m_badBoardsFromDB
bad-boards retrieved from DB.
unsigned nWiresInLayer(int layerId) const
Returns wire numbers in a layer.
double fieldWireBZ(int layerId) const
Returns backward z position of field wire in each layer.
void getClosestThetaPoints(const double alpha, const double theta, double &wth, unsigned short points[2]) const
Returns the two closest theta points for the input track incident angle (theta).
double senseWireFZ(int layerId) const
Returns forward z position of sense wire in each layer.
void readPropSpeed(const GearDir &gbxParams, int mode=0)
Read the propagation speed along the sense wire.
double innerRadiusInnerWall() const
Returns the inner radius of the inner wall.
DBObjPtr< CDCTimeZeros > * m_t0FromDB
t0s retrieved from DB.
int m_xtFileFormat
Format of xt input file.
int m_sigmaFileFormat
Format of sigma input file.
ushort m_firstSuperLayerOffset
Offset of the first super layer (for reduced CDC studies)
void setChMap()
Set channel map (from DB)
DBObjPtr< CDCFudgeFactorsForSigma > * m_fFactorFromDB
fudge factors retrieved from DB.
int m_nSLayer
The number of sense wire layer.
DBArray< CDCChannelMap > * m_chMapFromDB
channel map retrieved from DB.
void calcMeanT0(double minT0=3800, double maxT0=5800, int maxIt=10, double nStdv=3, double epsi=0.1)
Calculate mean t0 in ns (over all good wires)
void readFromDB(const CDCGeometry &)
Gets geometry parameters from database.
DBObjPtr< CDCXtRelations > * m_xtRelFromDB
xt params.
double m_nominalDriftV
Nominal drift velocity (4.0x10^-3 cm/nsec).
float m_FWirPos[c_maxNSenseLayers][c_maxNDriftCells][3]
Wire position incl.
float m_WireSagCoefAlign[c_maxNSenseLayers][c_maxNDriftCells]
Wire sag coefficient incl.
double m_zSForwardLayer[c_maxNSenseLayers]
The array to store forward z position of sense wire layers.
bool m_linearInterpolationOfSgm
Switch for linear interpolation of sigma.
DBObjPtr< CDCPropSpeeds > * m_propSpeedFromDB
prop.
void setBadWire()
Set bad-wires (from DB)
int m_xtParamMode
Mode for xt parameterization.
float m_FWirPosAlign[c_maxNSenseLayers][c_maxNDriftCells][3]
Wire position incl.
double fieldWireFZ(int layerId) const
Returns forward z position of field wire in each layer.
double m_dzSForwardLayer[c_maxNSenseLayers]
Corrections for forward z position of sense wire layers.
double getEDepToADCConvFactor(unsigned short layer, unsigned short cell, double edep, double dx, double costh)
Return edep-to-ADC conversion factor.
double m_nominalDriftVInv
Inverse of the nominal drift velocity.
const double * innerRadiusWireLayer() const
Returns an array of inner radius of wire layers.
double getDriftTime(double dist, unsigned short layer, unsigned short lr, double alpha, double theta) const
Return the drift time to the sense wire.
signed short getShiftInSuperLayer(unsigned short iSuperLayer, unsigned short iLayer) const
Returns shift in the super-layer.
void setBadBoard()
Set bad-boards (from DB)
double getWireSagCoef(EWirePosition set, uint layerId, int cellId) const
Returns coefficient for the sense wire sag.
void readChMap()
Read channel map between wire-id and electronics-id.
double m_zFForwardLayer[c_maxNFieldLayers]
The array to store forward z position of field wire layers.
double m_rSLayer[c_maxNSenseLayers]
The array to store radius of sense wire layers.
double m_fieldWireDiameter
The diameter of field wires.
double m_zFBackwardLayer[c_maxNFieldLayers]
The array to store backward z position of field wire layers.
ushort m_nFieldWires
Maximum number of Field Wires.
float m_WireSagCoefMisalign[c_maxNSenseLayers][c_maxNDriftCells]
Wire sag coefficient incl.
double getDriftLength(double dt, unsigned short layer, unsigned short lr, double alpha=0., double theta=0.5 *M_PI, bool calculateMinTime=true, double minTime=0.) const
Return the drift dength to the sense wire.
double m_cellSize[c_maxNSenseLayers]
The array to store cell size in each sense wire layer.
float m_BWirPosAlign[c_maxNSenseLayers][c_maxNDriftCells][3]
Wire position incl.
double outerRadiusOuterWall() const
Returns the outer radius of the outer wall.
float m_t0[c_maxNSenseLayers][c_maxNDriftCells]
t0 for each sense-wire (in nsec).
void setWirPosMisalignParams()
Set wire misalignment params.
double m_offSet[c_maxNSenseLayers]
The array to store z offset of sense wire layers.
void setTW()
Set time-walk parameters.
void setSResol()
Set spatial resolution (from DB).
double m_zSBackwardLayer[c_maxNSenseLayers]
The array to store backward z position of sense wire layers.
unsigned short getOldLeftRight(const B2Vector3D &posOnWire, const B2Vector3D &posOnTrack, const B2Vector3D &momentum) const
Returns old left/right.
float m_propSpeedInv[c_maxNSenseLayers]
Inverse of propagation speed of the sense wire.
ushort m_maxNCellsPerLayer
Maximum number wires within a layer.
double m_zWall[4][2]
The array to store z position of inner wall and outer wall.
bool isBadWire(const WireID &wid)
Inquire if the wire is totally-dead.
DBObjPtr< CDCAlignment > * m_alignmentFromDB
alignment params.
ushort m_maxNFieldLayers
Maximum number of Field Wire Layers.
bool m_misalignment
Switch for misalignment.
double m_senseWireDensity
The density of sense wires.
static CDCGeometryPar & Instance(const CDCGeometry *=nullptr)
Static method to get a reference to the CDCGeometryPar instance.
ushort m_firstLayerOffset
Offset of the first layer (for reduced CDC studies)
void Print() const
Print some debug information.
float m_timeWalkCoef[c_nBoards][2]
coefficients for time walk.
void readFFactor(const GearDir &gbxParams, int mode=0)
Read fudge factors.
unsigned short m_nThetaPoints
No.
double m_maxSpaceResol
max space resolution allowed (cm).
double m_thresholdEnergyDeposit
Energy thresh.
float m_BWirPos[c_maxNSenseLayers][c_maxNDriftCells][3]
Wire position incl.
float m_Sigma[c_maxNSenseLayers][2][c_maxNAlphaPoints][c_maxNThetaPoints][c_nSigmaParams]
position resolution for each layer.
void readXT(const GearDir &gbxParams, int mode=0)
Read XT-relation table.
double m_rFLayer[c_maxNFieldLayers]
The array to store radius of field wire layers.
void newReadXT(const GearDir &gbxParams, int mode=0)
Read XT-relation table in new format.
double m_nominalSpaceResol
Nominal spatial resolution (0.0130 cm).
unsigned nWireLayers() const
Returns a number of wire layers.
unsigned short m_nAlphaPoints
No.
unsigned short getNewLeftRightRaw(const B2Vector3D &posOnWire, const B2Vector3D &posOnTrack, const B2Vector3D &momentum) const
Returns new left/right_raw.
int m_senseWireZposMode
Mode for sense wire z position corr.
double m_rWall[4]
The array to store radius of inner wall and outer wall.
double m_tdcBinWidth
TDC bin width (nsec/bin).
bool m_modLeftRightFlag
Switch for modified left/right flag.
void getClosestThetaPoints4Sgm(const double alpha, const double theta, double &wth, unsigned short points[2]) const
Returns the two closest theta points for sigma for the input track incident angle (theta).
DBObjPtr< CDCEDepToADCConversions > * m_eDepToADCConversionsFromDB
Pointer to edep-to-ADC conv.
const double * outerRadiusWireLayer() const
Returns an array of outer radius of wire layers.
double senseWireR(int layerId) const
Returns radius of sense wire in each layer.
double m_senseWireTension
The tension of sense wires.
void setFFactor()
Set fudge factors (from DB).
bool getModLeftRightFlag() const
Get modified left/right flag.
double getMinTrackLength() const
Get minimum track length.
double getThresholdEnergyDeposit() const
Get threshold for Energy Deposit;.
static CDCSimControlPar & getInstance()
Static method to get a reference to the CDCSimControlPar instance.
bool getWireSag() const
Get wiresag flag.
Class for accessing arrays of objects in the database.
Definition DBArray.h:26
Class for accessing objects in the database.
Definition DBObjPtr.h:21
static std::string findFile(const std::string &path, bool silent=false)
Search for given file or directory in local or central release directory, and return absolute path if...
GearDir is the basic class used for accessing the parameter store.
Definition GearDir.h:31
virtual std::string getString(const std::string &path="") const noexcept(false) override
Get the parameter path as a string.
Definition GearDir.h:69
Optional DBObjPtr: This class behaves the same as the DBObjPtr except that it will not raise errors w...
Definition DBObjPtr.h:48
Class to identify a wire inside the CDC.
Definition WireID.h:34
unsigned short getICLayer() const
Getter for continuous layer numbering.
Definition WireID.cc:24
unsigned short getIWire() const
Getter for wire within the layer.
Definition WireID.h:145
unsigned short getEWire() const
Getter for encoded wire number.
Definition WireID.h:154
B2Vector3< double > B2Vector3D
typedef for common usage with double
Definition B2Vector3.h:516
double sqrt(double a)
sqrt for double
Definition beamHelpers.h:28
T dot(GeneralVector< T > a, GeneralVector< T > b)
dot product of two general vectors
void openFileB(boost::iostreams::filtering_istream &ifs, const std::string &fileName0)
Open a file using boost (to be able to read a gzipped file)
Definition OpenFile.cc:44
void openFileA(std::ifstream &ifs, const std::string &fileName0)
Open a file.
Definition OpenFile.cc:25
Abstract base class for different kinds of events.
STL namespace.