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