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