Belle II Software  release-05-01-25
ARICHReconstruction.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2010 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Luka Santelj, Rok Pestotnik *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include "arich/modules/arichReconstruction/ARICHReconstruction.h"
12 #include "arich/dbobjects/ARICHGeometryConfig.h"
13 #include "arich/modules/arichReconstruction/Utility.h"
14 #include "arich/dataobjects/ARICHHit.h"
15 #include "arich/dataobjects/ARICHTrack.h"
16 #include "arich/dataobjects/ARICHPhoton.h"
17 
18 // DataStore
19 #include <framework/datastore/StoreArray.h>
20 
21 // framework aux
22 #include <framework/logging/Logger.h>
23 #include <framework/gearbox/Const.h>
24 
25 #include <vector>
26 #include <TRotation.h>
27 #include <TRandom3.h>
28 
29 using namespace std;
30 using namespace boost;
31 
32 namespace Belle2 {
38  using namespace arich;
39 
40  ARICHReconstruction::ARICHReconstruction(int storePhot):
41  m_arichgp(),
42  m_recPars(),
43  m_trackPosRes(0),
44  m_trackAngRes(0),
45  m_alignMirrors(true),
46  m_nAerogelLayers(0),
47  m_storePhot(storePhot)
48  {
49  for (unsigned i = 0; i < c_noOfHypotheses; i++) {p_mass[i] = 0;}
50  for (unsigned i = 0; i < c_noOfAerogels; i++) {
51  m_refractiveInd[i] = 0;
52  m_zaero[i] = 0;
53  m_thickness[i] = 0;
54  m_transmissionLen[i] = 0;
55  m_n0[i] = 0;
56  m_anorm[i] = TVector3();
57  }
58 
59  }
60 
61  void ARICHReconstruction::initialize()
62  {
63 
64  for (const auto& part : Const::chargedStableSet) {
65  p_mass[part.getIndex()] = part.getMass();
66  }
67 
68  m_nAerogelLayers = m_arichgp->getAerogelPlane().getNLayers();
69  m_thickness[m_nAerogelLayers] = 0;
70  for (unsigned int i = 0; i < m_nAerogelLayers; i++) {
71  m_refractiveInd[i] = m_arichgp->getAerogelPlane().getLayerRefIndex(i + 1);
72  m_anorm[i] = TVector3(0, 0, 1);
73  m_thickness[i] = m_arichgp->getAerogelPlane().getLayerThickness(i + 1);
74  if (i == 0) m_zaero[i] = m_arichgp->getAerogelPlane().getAerogelZPosition() + m_thickness[i];
75  else m_zaero[i] = m_zaero[i - 1] + m_thickness[i];
76 
77  m_transmissionLen[i] = m_arichgp->getAerogelPlane().getLayerTrLength(i + 1) ; // aerogel transmission length;
78 
79  // measured FOM; aeroMerit is number of detected photons for beam of beta=1 and perpedicular incidence to aerogel tile
80  // (corrected for geometrical acceptance). n0[i] is then calculated from
81  // aeroMerit[i]=n0[i]*sin2(thc)*transmissionLen[i] * (1 - exp(-thickness[i] / transmissionLen[i])
82  m_n0[i] = m_recPars->getAerogelFOM(i) / ((1. - 1. / m_refractiveInd[i] / m_refractiveInd[i]) * m_transmissionLen[i] * (1 - exp(
83  -m_thickness[i] / m_transmissionLen[i])));
84  m_thickness[m_nAerogelLayers] += m_thickness[i];
85  }
86  m_refractiveInd[m_nAerogelLayers ] = 1.0;
87  m_refractiveInd[m_nAerogelLayers + 1] = m_arichgp->getHAPDGeometry().getWinRefIndex();
88  m_zaero[m_nAerogelLayers ] = m_arichgp->getDetectorZPosition();
89  m_zaero[m_nAerogelLayers + 1] = m_zaero[m_nAerogelLayers] + m_arichgp->getHAPDGeometry().getWinThickness();
90 
91  if (m_mirrAlign.hasChanged()) {
92  m_mirrorNorms.clear();
93  m_mirrorPoints.clear();
94  for (unsigned i = 1; i < m_arichgp->getMirrors().getNMirrors() + 1; i++) {
95  m_mirrorNorms.push_back(getMirrorNorm(i));
96  m_mirrorPoints.push_back(getMirrorPoint(i));
97  }
98  }
99 
100  if (m_tileAlign) {
101  if (m_tileAlign.hasChanged()) {
102  for (int iTile = 1; iTile < 125; iTile++) {
103  m_tilePars[iTile - 1][0] = m_tileAlign->getAlignmentElement(iTile).getAlpha();
104  m_tilePars[iTile - 1][1] = m_tileAlign->getAlignmentElement(iTile).getBeta();
105  }
106  }
107  }
108  }
109 
110 
111  int ARICHReconstruction::InsideDetector(TVector3 a, int copyno)
112  {
113  if (copyno == -1) return 0;
114  TVector2 origin;
115  origin.SetMagPhi(m_arichgp->getDetectorPlane().getSlotR(copyno), m_arichgp->getDetectorPlane().getSlotPhi(copyno));
116  TVector2 a2(a.X(), a.Y());
117  double phi = m_arichgp->getDetectorPlane().getSlotPhi(copyno);
118  TVector2 diff = a2 - origin;
119  diff = diff.Rotate(-phi);
120  const double size = m_arichgp->getHAPDGeometry().getAPDSizeX();
121  if (fabs(diff.X()) < size / 2. && fabs(diff.Y()) < size / 2.) {
122  int chX, chY;
123  m_arichgp->getHAPDGeometry().getXYChannel(diff.X(), diff.Y(), chX, chY);
124  if (chX < 0 || chY < 0) return 0;
125  int asicChannel = m_chnMap->getAsicFromXY(chX, chY);
126  // eliminate un-active channels
127  if (asicChannel < 0 || !m_chnMask->isActive(copyno, asicChannel)) return 0;
128  return 1;
129  }
130  return 0;
131  }
132 
133 
134  int ARICHReconstruction::smearTrack(ARICHTrack& arichTrack)
135  {
136  double a = gRandom->Gaus(0, m_trackAngRes);
137  double b = gRandom->Gaus(0, m_trackAngRes);
138  TVector3 dirf(a, b, sqrt(1 - a * a - b * b));
139  double dx = gRandom->Gaus(0, m_trackPosRes);
140  double dy = gRandom->Gaus(0, m_trackPosRes);
141  TVector3 mod = TVector3(dx, dy, 0);
142  TVector3 rpoint = arichTrack.getPosition() + mod;
143  TVector3 odir = arichTrack.getDirection();
144  double omomentum = arichTrack.getMomentum();
145  TVector3 rdir = TransformFromFixed(odir) * dirf; // global system
146  double rmomentum = omomentum;
147  arichTrack.setReconstructedValues(rpoint, rdir, rmomentum);
148  return 1;
149  }
150 
151 
152  TVector3 ARICHReconstruction::FastTracking(TVector3 dirf, TVector3 r, double* refractiveInd, double* z, int n, int opt)
153  {
154  //
155  // Description:
156  // The method calculates the intersection of the cherenkov photon
157  // with the detector plane
158 
159  // z[n+1]
160  // z[0] .. 1st aerogel exit
161  // z[n-1] .. 2nd aerogel exit
162 
163  double angmir = 0; int section[2] = {0, 0};
164 
165  unsigned tileID = m_arichgp->getAerogelPlane().getAerogelTileID(r.X(), r.Y());
166 
167  if (tileID == 0 && opt == 1) return TVector3();
168 
169  int nmir = m_arichgp->getMirrors().getNMirrors();
170  if (nmir > 0) {
171  double dangle = 2 * M_PI / nmir;
172  angmir = m_arichgp->getMirrors().getStartAngle() - dangle / 2.;
173 
174  double trkangle = r.XYvector().Phi() - angmir;
175  if (trkangle < 0) trkangle += 2 * M_PI;
176  if (trkangle > 2 * M_PI) trkangle -= 2 * M_PI;
177 
178  section[1] = int(trkangle / dangle) + 1;
179  }
180 
181  bool reflok = false; bool refl = false;
182  double path = (z[0] - r.z()) / dirf.z();
183  r += dirf * path;
184  for (int a = 1; a <= n + 1 ; a++) {
185  double rind = refractiveInd[a] / refractiveInd[a - 1];
186  dirf = Refraction(dirf, rind);
187  if (dirf.Mag() == 0) return TVector3();
188  path = (z[a] - r.z()) / dirf.z();
189  TVector3 r0 = r;
190  if (a == n && opt == 1) {
191  if (m_arichgp->getAerogelPlane().getAerogelTileID(r.X(), r.Y()) != tileID) return TVector3();
192  }
193  r += dirf * path;
194  TVector2 rxy = r.XYvector();
195  // check for possible reflections
196  if (a != n || nmir == 0) continue;
197  double angle = rxy.Phi() - angmir;
198  if (angle < 0) angle += 2 * M_PI;
199  if (angle > 2 * M_PI) angle -= 2 * M_PI;
200  double dangle = 2 * M_PI / nmir;
201  section[0] = int(angle / dangle) + 1;
202  if (r.Mag() > (r - 2 * m_mirrorPoints[section[0] - 1]).Mag()) {
203  refl = true;
204  int nrefl = 2;
205  if (section[0] == section[1]) nrefl = 1;
206  for (int k = 0; k < nrefl; k++) {
207  if (!HitsMirror(r0, dirf, section[k])) continue;
208  TVector3 mirpoint = m_mirrorPoints[section[k] - 1];
209  TVector3 mirnorm = m_mirrorNorms[section[k] - 1];
210  double s = dirf * mirnorm;
211  double s1 = (mirpoint - r0) * mirnorm;
212  r = r0 + s1 / s * dirf;
213  dirf = dirf - 2 * (dirf * mirnorm) * mirnorm;
214  path = (z[a] - r.z()) / dirf.z();
215  r += dirf * path;
216  reflok = true;
217  break;
218  }
219  }
220  }
221 
222  if (!reflok && refl) return TVector3();
223  return r;
224  }
225 
226  TVector3 ARICHReconstruction::HitVirtualPosition(const TVector3& hitpos, int mirrorID)
227  {
228 
229  if (mirrorID == 0) return hitpos;
230  TVector3 mirpoint = m_mirrorPoints[mirrorID - 1];
231  TVector3 mirnorm = m_mirrorNorms[mirrorID - 1];
232  return hitpos - 2 * ((hitpos - mirpoint) * mirnorm) * mirnorm;
233  }
234 
235 
236  bool ARICHReconstruction::HitsMirror(const TVector3& pos, const TVector3& dir, int mirrorID)
237  {
238 
239  TVector3 mirnorm = m_mirrorNorms[mirrorID - 1];
240  TVector3 mirpoint = m_mirrorPoints[mirrorID - 1];
241  TRotation rot = TransformToFixed(mirnorm);
242  TVector3 dirTr = rot * dir;
243  if (dirTr.Z() < 0) return 0; // comes from outter side
244  TVector3 posTr = rot * (pos - mirpoint);
245  TVector3 pointOnMirr = posTr - (posTr.Z() / dirTr.Z()) * dirTr;
246  if (fabs(pointOnMirr.Y()) < m_arichgp->getMirrors().getPlateLength() / 2.
247  && fabs(pointOnMirr.X()) < m_arichgp->getMirrors().getPlateWidth() / 2.) return 1;
248 
249  return 0;
250  }
251 
252 
253  int ARICHReconstruction::CherenkovPhoton(TVector3 r, TVector3 rh,
254  TVector3& rf, TVector3& dirf,
255  double* refractiveInd, double* z, int n, int mirrorID)
256  {
257  //
258  // Description:
259  // The method calculates the direction of the cherenkov photon
260  // and intersection with the aerogel exit point
261  //
262  // Called by: CerenkovAngle
263 
264 
265  // Arguments:
266  // Input:
267  // dir, r track position
268  // rh photon hit
269 
270 
271  // Output:
272  // rf aerogel exit from which the photon was emitted
273  // dirf photon direction in aerogel
274  static TVector3 norm(0, 0, 1); // detector plane normal vector
275 
276  double dwin = m_arichgp->getHAPDGeometry().getWinThickness();
277  double refractiveInd0 = m_arichgp->getHAPDGeometry().getWinRefIndex();
278 
279  // iteration is stoped when the difference of photon positions on first aerogel exit
280  // between two iterations is smaller than this value.
281  const double dmin = 0.0000001;
282  const int niter = 100; // maximal number of iterations
283  TVector3 dirw;
284  TVector3 rh1 = rh - dwin * norm;
285 
286  std::vector <TVector3> rf0(n + 2);
287  // rf0[0] .. track point
288  // rf0[1] 1. 1st aerogel exit
289  // rf0[n] n. aerogel exit ...
290  std::vector <TVector3> dirf0(n + 2);
291  // dirf0[0] .. 1st aerogel direction
292  // dirf0[1] .. 2nd aerogel direction
293  // dirf0[n] .. direction after aerogels
294 
295  // z[0] .. 1st aerogel exit
296  // z[n-1] .. 2nd aerogel exit
297  // z[n] .. detector position
298  // z[n+1] .. detector + window
299 
300  rf0[0] = r;
301  rf0[1] = rf;
302 
303  for (int iter = 0; iter < niter; iter++) {
304 
305  // direction in the space between aerogels and detector
306  // *************************************
307  if (iter == 0) dirf0[n] = (rh1 - rf0[1]).Unit();
308  else dirf0[n] = (rh1 - rf0[n]).Unit();
309 
310  // *************************************
311  // n-layers of aerogel // refractiveInd relative refractive index
312  for (int a = n - 1; a >= 0 ; a--) {
313  double rind = refractiveInd[a] / refractiveInd[a + 1];
314  dirf0[a] = Refraction(dirf0[a + 1], rind);
315  }
316 
317  double path = (z[0] - r.z()) / dirf0[0].z();
318  double x1 = rf0[1].x();
319  double y1 = rf0[1].y();
320  for (int a = 0; a < n ; a++) {
321  rf0[a + 1] = rf0[a] + dirf0[a] * path;
322  path = (z[a + 1] - rf0[a + 1].z()) / dirf0[a + 1].z();
323  }
324 
325  Refraction(dirf0[n], norm, refractiveInd0, dirw);
326 
327  // *************************************
328 
329  path = dwin / (dirw * norm);
330  rh1 = rh - dirw * path;
331 
332  double d2 = (rf0[1].x() - x1) * (rf0[1].x() - x1) + (rf0[1].y() - y1) * (rf0[1].y() - y1);
333 
334  if ((d2 < dmin) && (iter)) {
335 
336  // if mirror hypothesis check if reflection point lies on mirror plate
337  if (mirrorID) {
338  if (!HitsMirror(rf0[n], dirf0[n], mirrorID)) return -1;
339  }
340 
341  rf = rf0[1];
342  dirf = dirf0[0];
343  return iter;
344  }
345  }
346  return -1;
347  }
348 
349  int ARICHReconstruction::likelihood2(ARICHTrack& arichTrack, const StoreArray<ARICHHit>& arichHits,
350  ARICHLikelihood& arichLikelihood)
351  {
352 
353  const unsigned int nPhotonHits = arichHits.getEntries(); // number of detected photons
354 
355  if (m_nAerogelLayers + 1 > c_noOfAerogels) B2ERROR("ARICHReconstrucion: number of aerogel layers defined in the xml file exceeds "
356  << c_noOfAerogels);
357 
358  double logL[c_noOfHypotheses] = {0.0};
359  double nSig_w_acc[c_noOfHypotheses][c_noOfAerogels] = { {0.0} }; // expected no. of signal photons, including geometrical acceptance
360  double nSig_wo_acc[c_noOfHypotheses][c_noOfAerogels][20] = { { {0.0} } }; // expected no. of signal photons, without geometrical acceptance, divided in 20 phi bins (used for PDF normalization)
361  double nSig_wo_accInt[c_noOfHypotheses][c_noOfAerogels] = { {0.0} }; // expected no. of signal photons, without geometrical acceptance, integrated over phi
362  double esigi[c_noOfHypotheses] = {0.0}; // expected number of signal photons in hit pixel
363  double thetaCh[c_noOfHypotheses][c_noOfAerogels] = { {0.0} }; // expected Cherenkov angle
364 
365  // read some geometry parameters
366  double padSize = m_arichgp->getHAPDGeometry().getPadSize();
367  int nMirSeg = m_arichgp->getMirrors().getNMirrors();
368  double angmir = m_arichgp->getMirrors().getStartAngle();
369 
370  // Detected photons in cherenkov ring (integrated over 0.1<theta<0.5)
371  int nDetPhotons = 0;
372 
373  double ebgri[c_noOfHypotheses] = {0.0}; // number of expected background photons per pad
374  double nBgr[c_noOfHypotheses] = {0.0}; // total number of expected background photons (in 0.1-0.5 rad ring)
375 
376  // reconstructed track direction
377  TVector3 edir = arichTrack.getDirection();
378  if (edir.z() < 0) return 0;
379  double momentum = arichTrack.getMomentum();
380 
381  double thcResolution = m_recPars->getThcResolution(momentum);
382  if (thcResolution < 0) thcResolution = 0.01; // happens for spurious tracks with 100 GeV momentum!
383 
384  double wideGaussFract = (m_recPars->getParameters())[0];
385  double wideGaussSigma = (m_recPars->getParameters())[1];
386 
387  unsigned tileID = m_arichgp->getAerogelPlane().getAerogelTileID(arichTrack.getPosition().X(), arichTrack.getPosition().Y());
388  double r = arichTrack.getPosition().XYvector().Mod();
389  if (tileID > 0) correctEmissionPoint(tileID, r);
390 
391  //------------------------------------------------------
392  // Calculate number of expected detected photons (emmited x geometrical acceptance).
393  // -----------------------------------------------------
394 
395  float nphot_scaling = 20.; // number of photons to be traced is (expected number of emitted photons * nphot_scaling)
396  int nStep = 5; // number of steps in one aerogel layer
397 
398  // loop over all particle hypotheses
399  for (int iHyp = 0; iHyp < c_noOfHypotheses; iHyp++) {
400 
401  // absorbtion factor (scattering)
402  double abs = 1;
403 
404  // loop over aerogel layers
405  for (int iAerogel = m_nAerogelLayers - 1; iAerogel >= 0; iAerogel--) {
406 
407  thetaCh[iHyp][iAerogel] = ExpectedCherenkovAngle(momentum,
408  p_mass[iHyp],
409  m_refractiveInd[iAerogel]);
410 
411  // track length in the radiator
412  double pathLengthRadiator = arichTrack.getDirection() * m_anorm[iAerogel];
413  if (pathLengthRadiator) pathLengthRadiator = m_thickness[iAerogel] / pathLengthRadiator;
414 
415  // step length
416  double dxx = pathLengthRadiator / double(nStep);
417  // number of photons to be emmited per step (number of expected photons * nphot_scaling)
418  double nPhot = m_n0[iAerogel] * sin(thetaCh[iHyp][iAerogel]) * sin(thetaCh[iHyp][iAerogel]) * dxx * nphot_scaling;
419  TVector3 exit_point = getTrackPositionAtZ(arichTrack, m_zaero[iAerogel]);
420 
421  // loop over emmision point steps
422  for (int iepoint = 0; iepoint < nStep; iepoint++) {
423 
424  TVector3 epoint = exit_point - (0.5 + iepoint) * dxx * edir;
425  abs *= exp(-dxx / m_transmissionLen[iAerogel]);
426  unsigned int genPhot = nPhot * abs; // number of photons to emmit in current step, including scattering correction
427 
428  // loop over emmited "photons"
429  for (unsigned int iPhoton = 0; iPhoton < genPhot; iPhoton++) {
430  double fi = 2 * M_PI * iPhoton / float(genPhot); // uniformly distributed in phi
431  TVector3 adirf = setThetaPhi(thetaCh[iHyp][iAerogel], fi); // photon direction in track system
432  adirf = TransformFromFixed(edir) * adirf; // photon direction in global system
433  int ifi = int (fi * 20 / 2. / M_PI); // phi bin
434  // track photon from emission point to the detector plane
435  TVector3 dposition = FastTracking(adirf, epoint, &m_refractiveInd[iAerogel], &m_zaero[iAerogel], m_nAerogelLayers - iAerogel, 1);
436  if (dposition.Mag() > 1.0) {nSig_wo_acc[iHyp][iAerogel][ifi] += 1; nSig_wo_accInt[iHyp][iAerogel] += 1;}
437  else continue;
438  unsigned copyno = m_arichgp->getDetectorPlane().pointSlotID(dposition.X(), dposition.Y());
439  if (!copyno) continue;
440  // check if photon fell on photosensitive area
441  if (InsideDetector(dposition, copyno)) nSig_w_acc[iHyp][iAerogel] += 1;
442  }
443  }
444 
445  // scale the obtained numbers
446  for (int ik = 0; ik < 20; ik++) {
447  nSig_wo_acc[iHyp][iAerogel][ik] /= nphot_scaling;
448  }
449  nSig_w_acc[iHyp][iAerogel] /= nphot_scaling;
450  nSig_wo_accInt[iHyp][iAerogel] /= nphot_scaling;
451  } // for (unsigned int iAerogel=0;iAerogel<m_nAerogelLayers;iAerogel++)
452 
453  // sum up contribution from all aerogel layers
454  for (unsigned int iAerogel = 0; iAerogel < m_nAerogelLayers; iAerogel++) {
455  nSig_w_acc[iHyp][m_nAerogelLayers] += nSig_w_acc[iHyp][iAerogel];
456  nSig_wo_accInt[iHyp][m_nAerogelLayers] += nSig_wo_accInt[iHyp][iAerogel];
457 
458  for (int ik = 0; ik < 20; ik++) {
459  nSig_wo_acc[iHyp][m_nAerogelLayers][ik] += nSig_wo_acc[iHyp][iAerogel][ik];
460  }
461  }
462 
463  // get number of expected background photons in ring (integrated from 0.1 to 0.5 rad)
464  std::vector<double> bkgPars = {momentum / sqrt(p_mass[iHyp]*p_mass[iHyp] + momentum * momentum), double(arichTrack.hitsWindow())};
465  nBgr[iHyp] = m_recPars->getExpectedBackgroundHits(bkgPars);
466 
467  } // for (int iHyp=0;iHyp < c_noOfHypotheses; iHyp++ )
468  //#####################################################
469 
470  TVector3 track_at_detector = getTrackPositionAtZ(arichTrack, m_zaero[m_nAerogelLayers + 1]);
471 
472  // the id numbers of mirrors from which the photons could possibly reflect are calculated
473  int mirrors[3];
474  mirrors[0] = 0; // for no reflection
475  int refl = 1;
476 
477  // only if particle track on detector is at radius larger than 850mm (for now hardcoded)
478  // possible reflections are taken into account.
479  if (track_at_detector.XYvector().Mod() > 85.0) {
480  double trackang = track_at_detector.Phi() - angmir;
481  if (trackang < 0) trackang += 2 * M_PI;
482  if (trackang > 2 * M_PI) trackang -= 2 * M_PI;
483  int section1 = int(trackang * nMirSeg / 2. / M_PI) + 1;
484  int section2 = section1 + 1;
485  if (section1 == nMirSeg) section2 = 1;
486  mirrors[1] = section1; mirrors[2] = section2;
487  refl = 3;
488  }
489 
490  // loop over all detected photon hits
491 
492  for (unsigned int iPhoton = 0; iPhoton < nPhotonHits; iPhoton++) {
493 
494  ARICHHit* h = arichHits[iPhoton];
495  int modID = h->getModule();
496  int channel = h->getChannel();
497  TVector3 hitpos = m_arichgp->getMasterVolume().pointToLocal(h->getPosition());
498  bool bkgAdded = false;
499  int nfoo = nDetPhotons;
500  for (int iHyp = 0; iHyp < c_noOfHypotheses; iHyp++) { esigi[iHyp] = 0; ebgri[iHyp] = 0;}
501 
502  bool reflOK = true; // remove window photons from reflected hypothesis
503 
504  // loop over possible mirror reflections
505  for (int mirr = 0; mirr < refl; mirr++) {
506 
507  if (!reflOK) break; // photon from window so break
508 
509  // calculate fi_ch for a given track refl
510  TVector3 virthitpos = HitVirtualPosition(hitpos, mirrors[mirr]);
511 
512  // if hit is more than 25cm from the track position on the detector plane, skip it.
513  // (not reconstructing hits with irrelevantly large Cherenkov angle)
514  if ((track_at_detector - virthitpos).Mag() > 25.0) continue;
515 
516  double sigExpArr[c_noOfHypotheses] = {0.0}; // esigi for given mirror hypothesis only
517  double th_cer_all[c_noOfAerogels] = {0.0};
518  double fi_cer_all[c_noOfAerogels] = {0.0};
519 
520  double weight[c_noOfHypotheses][c_noOfAerogels] = { {0.0} };
521  double weight_sum[c_noOfHypotheses] = {0.0};
522  int proc = 0;
523  double fi_cer_trk = 0.;
524 
525  // loop over all aerogel layers
526  for (unsigned int iAerogel = 0; iAerogel < m_nAerogelLayers; iAerogel++) {
527 
528  TVector3 initialrf = getTrackPositionAtZ(arichTrack, m_zaero[iAerogel]);
529  TVector3 epoint = getTrackMeanEmissionPosition(arichTrack, iAerogel);
530  TVector3 edirr = arichTrack.getDirection();
531  TVector3 photonDirection; // calculated photon direction
532 
533  if (CherenkovPhoton(epoint, virthitpos, initialrf, photonDirection, &m_refractiveInd[iAerogel], &m_zaero[iAerogel],
534  m_nAerogelLayers - iAerogel, mirrors[mirr]) < 0) break;
535 
536  TVector3 dirch = TransformToFixed(edirr) * photonDirection;
537  double fi_cer = dirch.Phi();
538  double th_cer = dirch.Theta();
539 
540 
541  th_cer_all[iAerogel] = th_cer;
542  fi_cer_all[iAerogel] = fi_cer;
543  fi_cer_trk = dirch.XYvector().DeltaPhi(edirr.XYvector());
544 
545  if (mirr == 0 && th_cer < 0.1) reflOK = false;
546  // skip photons with irrelevantly large/small Cherenkov angle
547  if (th_cer > 0.5 || th_cer < 0.1) continue;
548 
549  // count photons with 0.1<thc<0.5
550  if (nfoo == nDetPhotons) nDetPhotons++;
551 
552  if (fi_cer < 0) fi_cer += 2 * M_PI;
553  double fii = fi_cer;
554  if (mirr > 0) {
555  double fi_mir = m_mirrorNorms[mirrors[mirr] - 1].XYvector().Phi();
556  fii = 2 * fi_mir - fi_cer - M_PI;
557  }
558 
559 
560  // loop over all particle hypotheses
561  for (int iHyp = 0; iHyp < c_noOfHypotheses; iHyp++) {
562 
563  // track a photon from the mean emission point to the detector surface
564  TVector3 photonDirection1 = setThetaPhi(thetaCh[iHyp][iAerogel], fi_cer); // particle system
565  photonDirection1 = TransformFromFixed(edirr) * photonDirection1; // global system
566  int ifi = int (fi_cer * 20 / 2. / M_PI);
567  TVector3 photonAtAerogelExit = photonDirection1 * (m_thickness[iAerogel] / photonDirection1.z());
568  TVector3 trackAtAerogelExit = edirr * (m_thickness[iAerogel] / edirr.z());
569  TVector3 dtrackphoton = photonAtAerogelExit - trackAtAerogelExit;
570  TVector3 detector_position;
571 
572  detector_position = FastTracking(photonDirection1, epoint, &m_refractiveInd[iAerogel], &m_zaero[iAerogel],
573  m_nAerogelLayers - iAerogel, 0);
574 
575  TVector3 meanr = detector_position - epoint;
576  double path = meanr.Mag();
577  meanr = meanr.Unit();
578 
579  double detector_sigma = thcResolution * path / meanr.z();
580  double wide_sigma = wideGaussSigma * path / meanr.z();
581  // calculate pad orientation and distance relative to that photon
582  double modphi = m_arichgp->getDetectorPlane().getSlotPhi(modID);
583 
584  double pad_fi = fii - modphi;
585  double dx = (detector_position - hitpos).Mag();
586  double dr = (track_at_detector - detector_position).Mag();
587 
588  if (dr > 0.01) {
589  double normalizacija = nSig_wo_acc[iHyp][iAerogel][ifi] * padSize / (0.1 * M_PI * dr);
590  weight[iHyp][iAerogel] = normalizacija;
591  weight_sum[iHyp] += weight[iHyp][iAerogel];
592  double integralMain = SquareInt(padSize, pad_fi, dx, detector_sigma) / sqrt(2.);
593  double integralWide = SquareInt(padSize, pad_fi, dx, wide_sigma) / sqrt(2.);
594  // expected number of signal photons in each pixel
595  double exp = normalizacija * ((1 - wideGaussFract) * integralMain + wideGaussFract * integralWide);
596  esigi[iHyp] += exp;
597  sigExpArr[iHyp] += exp;
598  } // if (dr>0 && thetaCh[iHyp][iAerogel])
599 
600  }// for (int iHyp=0;iHyp< c_noOfHypotheses; iHyp++)
601  if (iAerogel == m_nAerogelLayers - 1) proc = 1; // successfully processed for all layers
602  }// for (unsigned int iAerogel=0; iAerogel<m_nAerogelLayers;iAerogel++)
603 
604  if (proc) {
605  // add background contribution if not yet (add only once)
606  if (!bkgAdded) {
607  for (int iHyp = 0; iHyp < c_noOfHypotheses; iHyp++) {
608  std::vector<double> pars = {momentum / sqrt(p_mass[iHyp]*p_mass[iHyp] + momentum * momentum), double(arichTrack.hitsWindow())};
609  ebgri[iHyp] += m_recPars->getBackgroundPerPad(th_cer_all[1], pars);
610  }
611  bkgAdded = true;
612  }
613  }
614  // create ARICHPhoton if desired
615  if (m_storePhot && th_cer_all[1] > 0 && th_cer_all[1] < 0.6) {
616  double n_cos_theta_ch[c_noOfHypotheses] = {0.0};
617  double phi_ch[c_noOfHypotheses] = {0.0};
618  for (int iHyp = 0; iHyp < c_noOfHypotheses; iHyp++) {
619  if (weight_sum[iHyp] > 0) {
620  for (unsigned int iAerogel = 0; iAerogel < m_nAerogelLayers; iAerogel++) {
621  double emission_prob = weight[iHyp][iAerogel] / weight_sum[iHyp];
622  n_cos_theta_ch[iHyp] += emission_prob * m_refractiveInd[iAerogel] * cos(th_cer_all[iAerogel]);
623  phi_ch[iHyp] += emission_prob * fi_cer_all[iAerogel];
624  }
625  //std::cout << iHyp << " " << n_cos_theta_ch[iHyp] << " " << phi_ch[iHyp] << std::endl;
626  } else {
627  n_cos_theta_ch[iHyp] = -99999.;
628  phi_ch[iHyp] = -99999.;
629  }
630  }
631  ARICHPhoton phot(iPhoton, th_cer_all[1], fi_cer_all[1], mirrors[mirr]); // th_cer of the first aerogel layer assumption is stored
632  phot.setBkgExp(ebgri); // store expected number of background hits
633  phot.setSigExp(sigExpArr); // store expected number of signal hits
634  phot.setPhiCerTrk(fi_cer_trk); // store phi angle in track coordinates
635  phot.setNCosThetaCh(n_cos_theta_ch); // store n cos(theta_th) for all particle hypotheses
636  phot.setPhiCh(phi_ch); // store phi_ch for all particle hypotheses
637  phot.setXY(hitpos.X(), hitpos.Y()); // store x-y hit position
638  phot.setModuleID(modID); // store module id
639  phot.setChannel(channel); // store channel
640  arichTrack.addPhoton(phot);
641  }
642 
643 
644  }// for (int mirr = 0; mirr < refl; mirr++)
645 
646  //******************************************
647  // LIKELIHOOD construction
648  //*******************************************
649 
650  for (int iHyp = 0; iHyp < c_noOfHypotheses; iHyp++) {
651  double expected = esigi[iHyp] + ebgri[iHyp];
652  if (bkgAdded) logL[iHyp] += expected + log(1 - exp(-expected));
653  }
654 
655  } // for (unsigned int iPhoton=0; iPhoton< nPhotonHits; iPhoton++)
656 
657  //*********************************************
658  // add constant term to the LIKELIHOOD function
659  //*********************************************
660  double exppho[c_noOfHypotheses] = {0.0};
661  for (int iHyp = 0; iHyp < c_noOfHypotheses; iHyp++) {
662  exppho[iHyp] = nSig_w_acc[iHyp][m_nAerogelLayers] * (1 - wideGaussFract) + wideGaussFract * 0.7 *
663  nSig_wo_accInt[iHyp][m_nAerogelLayers] + nBgr[iHyp];
664  logL[iHyp] -= exppho[iHyp];
665  if (isnan(logL[iHyp]) || isinf(logL[iHyp])) {
666  B2WARNING("ARICHReconstruction: log likelihood value infinite! Flat background hit probability is " << ebgri[iHyp] << "!");
667  logL[iHyp] = 0;
668  }
669  }
670 
671  //******************************************
672  // store LikeliHOOD info
673  //******************************************
674 
675  int flag = 1;
676  if ((thetaCh[0][0] > 0 || thetaCh[0][1] > 0) && nSig_w_acc[0][m_nAerogelLayers] == 0) flag = 0;
677 
678  // set values of ARICHLikelihood
679  arichLikelihood.setValues(flag, logL, nDetPhotons, exppho);
680 
681  return 1;
682  }
683 
684  void ARICHReconstruction::setTrackPositionResolution(double pRes)
685  {
686  m_trackPosRes = pRes;
687  }
688  void ARICHReconstruction::setTrackAngleResolution(double aRes)
689  {
690  m_trackAngRes = aRes;
691  }
692 
693  TVector3 ARICHReconstruction::getTrackMeanEmissionPosition(const ARICHTrack& track, int iAero)
694  {
695  // Emission length measured from aerogel exit
696 
697  TVector3 dir = track.getDirection();
698  if (dir.Z() == 0) return TVector3();
699  double d = m_thickness[iAero] / dir.Z() / m_transmissionLen[iAero];
700  double dmean = 1 - d / expm1(d);
701  //double dmean = -log((1 + exp(-d)) / 2.);
702  double mel = dmean * m_transmissionLen[iAero];
703 
704  return (getTrackPositionAtZ(track, m_zaero[iAero]) - mel * dir);
705  }
706 
707  TVector3 ARICHReconstruction::getTrackPositionAtZ(const ARICHTrack& track, double zout)
708  {
709  TVector3 dir = track.getDirection();
710  TVector3 pos = track.getPosition();
711  if (dir.Z() == 0) return TVector3(0, 0, 0);
712  double path = (zout - pos.Z()) / dir.Z();
713  return pos + dir * path;
714  }
715 
716  void ARICHReconstruction::transformTrackToLocal(ARICHTrack& arichTrack, bool align)
717  {
718  // tranform track from Belle II to local ARICH frame
719  TVector3 locPos = m_arichgp->getMasterVolume().pointToLocal(arichTrack.getPosition());
720  TVector3 locDir = m_arichgp->getMasterVolume().momentumToLocal(arichTrack.getDirection());
721 
722  // apply the alignment correction
723  if (align && m_alignp.isValid()) {
724  // apply global alignment correction
725  locPos = m_alignp->pointToLocal(locPos);
726  locDir = m_alignp->momentumToLocal(locDir);
727  }
728 
729  // set parameters and return
730  // is it needed to extrapolate to z of aerogel in local frame?? tabun
731  arichTrack.setReconstructedValues(locPos, locDir, arichTrack.getMomentum());
732  return;
733  }
734 
735 
736  TVector3 ARICHReconstruction::getMirrorPoint(int mirrorID)
737  {
738 
739  TVector3 mirpoint = m_arichgp->getMirrors().getPoint(mirrorID);
740  if (m_alignMirrors && m_mirrAlign.isValid()) mirpoint += m_mirrAlign->getAlignmentElement(mirrorID).getTranslation();
741  return mirpoint;
742 
743  }
744 
745  TVector3 ARICHReconstruction::getMirrorNorm(int mirrorID)
746  {
747  if (m_alignMirrors && m_mirrAlign.isValid()) {
748 
749  TVector3 mirnorm = m_arichgp->getMirrors().getNormVector(mirrorID);
750  mirnorm.SetTheta(mirnorm.Theta() + m_mirrAlign->getAlignmentElement(mirrorID).getAlpha());
751  mirnorm.SetPhi(mirnorm.Phi() + m_mirrAlign->getAlignmentElement(mirrorID).getBeta());
752  return mirnorm;
753 
754  }
755  return m_arichgp->getMirrors().getNormVector(mirrorID);
756  }
757 
758  void ARICHReconstruction::correctEmissionPoint(int tileID, double r)
759  {
760 
761  double ang = m_tilePars[tileID - 1][0] + m_tilePars[tileID - 1][1] * r;
762  m_zaero[0] = m_arichgp->getAerogelPlane().getAerogelZPosition() + m_thickness[0] - ang * 50.;
763  m_zaero[1] = m_zaero[0] + m_thickness[1];
764 
765  }
766 
768 }
Belle2::ARICHReconstruction::m_trackPosRes
double m_trackPosRes
track position resolution (from tracking)
Definition: ARICHReconstruction.h:113
Belle2::ARICHReconstruction::m_tileAlign
OptionalDBObjPtr< ARICHAeroTilesAlignment > m_tileAlign
alignment of aerogel tiles from DB
Definition: ARICHReconstruction.h:108
Belle2::ARICHPhoton::setPhiCerTrk
void setPhiCerTrk(float phi)
Set hit phi angle in track coordinates.
Definition: ARICHPhoton.h:58
Belle2::ARICHTrack::getMomentum
double getMomentum() const
returns track momentum
Definition: ARICHTrack.h:155
Belle2::ARICHReconstruction::InsideDetector
int InsideDetector(TVector3 a, int copyno)
Returns 1 if vector "a" lies on "copyno"-th detector active surface of detector and 0 else.
Definition: ARICHReconstruction.cc:111
Belle2::ARICHReconstruction::FastTracking
TVector3 FastTracking(TVector3 dirf, TVector3 r, double *refind, double *z, int n, int opt)
Calculates the intersection of the Cherenkov photon emitted from point "r" in "dirf" direction with t...
Definition: ARICHReconstruction.cc:152
Belle2::ARICHReconstruction::m_transmissionLen
double m_transmissionLen[c_noOfAerogels]
transmission lengths of aerogel layers
Definition: ARICHReconstruction.h:121
Belle2::ARICHLikelihood
This is a class to store ARICH likelihoods in the datastore.
Definition: ARICHLikelihood.h:38
Belle2::ARICHReconstruction::m_alignMirrors
bool m_alignMirrors
if set to true mirror alignment constants from DB are used
Definition: ARICHReconstruction.h:115
Belle2::ARICHReconstruction::m_recPars
DBObjPtr< ARICHReconstructionPar > m_recPars
reconstruction parameters from the DB
Definition: ARICHReconstruction.h:103
Belle2::ARICHPhoton::setPhiCh
void setPhiCh(const double *phi_ch)
Set phi_ch.
Definition: ARICHPhoton.h:132
Belle2::ARICHReconstruction::m_mirrorNorms
std::vector< TVector3 > m_mirrorNorms
vector of nomal vectors of all mirror plates
Definition: ARICHReconstruction.h:111
Belle2::ARICHTrack::getDirection
TVector3 getDirection() const
returns track direction vector
Definition: ARICHTrack.h:149
Belle2::ARICHReconstruction::p_mass
double p_mass[c_noOfHypotheses]
particle masses
Definition: ARICHReconstruction.h:100
Belle2::ARICHPhoton::setModuleID
void setModuleID(int modID)
Set id of hit module.
Definition: ARICHPhoton.h:75
Belle2::ARICHReconstruction::m_anorm
TVector3 m_anorm[c_noOfAerogels]
normal vector of the aerogle plane
Definition: ARICHReconstruction.h:123
Belle2::Const::chargedStableSet
static const ParticleSet chargedStableSet
set of charged stable particles
Definition: Const.h:494
Belle2::ARICHReconstruction::HitVirtualPosition
TVector3 HitVirtualPosition(const TVector3 &hitpos, int mirrorID)
Returns the hit virtual position, assuming that it was reflected from mirror.
Definition: ARICHReconstruction.cc:226
Belle2::ARICHReconstruction::m_mirrAlign
DBObjPtr< ARICHMirrorAlignment > m_mirrAlign
global alignment parameters from the DB
Definition: ARICHReconstruction.h:107
Belle2::ARICHReconstruction::m_n0
double m_n0[c_noOfAerogels]
number of emmited photons per unit length
Definition: ARICHReconstruction.h:122
Belle2::ARICHReconstruction::m_storePhot
int m_storePhot
set to 1 to store individual reconstructed photon information
Definition: ARICHReconstruction.h:124
Belle2::ARICHReconstruction::m_thickness
double m_thickness[c_noOfAerogels]
thicknesses of areogel layers
Definition: ARICHReconstruction.h:120
Belle2::ARICHTrack::setReconstructedValues
void setReconstructedValues(TVector3 r, TVector3 dir, double p)
Sets the reconstructed value of track parameters.
Definition: ARICHTrack.h:176
Belle2::ARICHReconstruction::m_alignp
DBObjPtr< ARICHGlobalAlignment > m_alignp
global alignment parameters from the DB
Definition: ARICHReconstruction.h:106
Belle2::ARICHReconstruction::smearTrack
int smearTrack(ARICHTrack &arichTrack)
Smeares track parameters ("simulate" the uncertainties of tracking).
Definition: ARICHReconstruction.cc:134
Belle2::ARICHLikelihood::setValues
void setValues(int flag, const double *logL, int detPhot, const double *expPhot)
Set values.
Definition: ARICHLikelihood.h:58
Belle2::ARICHPhoton::setBkgExp
void setBkgExp(const double *bkgExp)
Set expected background contribution.
Definition: ARICHPhoton.h:106
Belle2::ARICHTrack::hitsWindow
bool hitsWindow() const
Returns true if track hits HAPD window.
Definition: ARICHTrack.h:161
Belle2::ARICHReconstruction::initialize
void initialize()
read geomerty parameters from xml and initialize class memebers
Definition: ARICHReconstruction.cc:61
Belle2::ARICHReconstruction::getTrackPositionAtZ
TVector3 getTrackPositionAtZ(const ARICHTrack &track, double zout)
Returns track direction at point with z coordinate "zout" (assumes straight track).
Definition: ARICHReconstruction.cc:707
Belle2::ARICHReconstruction::HitsMirror
bool HitsMirror(const TVector3 &pos, const TVector3 &dir, int mirrorID)
returns true if photon at position pos with direction dir hits mirror plate with ID mirrorID
Definition: ARICHReconstruction.cc:236
Belle2::ARICHHit
Datastore class that holds photon hits. Input to the reconstruction.
Definition: ARICHHit.h:33
Belle2::ARICHPhoton::setXY
void setXY(float x, float y)
Set X-Y position of hit.
Definition: ARICHPhoton.h:66
Belle2::ARICHReconstruction::setTrackPositionResolution
void setTrackPositionResolution(double pRes)
Sets track position resolution (from tracking)
Definition: ARICHReconstruction.cc:684
Belle2::ARICHReconstruction::getMirrorPoint
TVector3 getMirrorPoint(int mirrorID)
Returns point on the mirror plate with id mirrorID.
Definition: ARICHReconstruction.cc:736
Belle2::ARICHPhoton::setNCosThetaCh
void setNCosThetaCh(const double *n_cos_theta_ch)
Set n cos(theta_ch)
Definition: ARICHPhoton.h:119
Belle2::ARICHTrack::getPosition
TVector3 getPosition() const
returns track position vector
Definition: ARICHTrack.h:143
Belle2::ARICHReconstruction::transformTrackToLocal
void transformTrackToLocal(ARICHTrack &arichTrack, bool align)
Transforms track parameters from global Belle2 to ARICH local frame.
Definition: ARICHReconstruction.cc:716
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::asicChannel
record to be used to store ASIC info
Definition: CDCCrossTalkClasses.h:23
Belle2::ARICHTrack
Datastore class that holds position and momentum information of tracks that hit ARICH.
Definition: ARICHTrack.h:42
Belle2::ARICHPhoton::setSigExp
void setSigExp(const double *sigExp)
Set expected signal contribution.
Definition: ARICHPhoton.h:93
Belle2::ARICHReconstruction::m_mirrorPoints
std::vector< TVector3 > m_mirrorPoints
vector of points on all mirror plates
Definition: ARICHReconstruction.h:110
Belle2::ARICHReconstruction::getTrackMeanEmissionPosition
TVector3 getTrackMeanEmissionPosition(const ARICHTrack &track, int iAero)
Returns mean emission position of Cherenkov photons from i-th aerogel layer.
Definition: ARICHReconstruction.cc:693
Belle2::ARICHReconstruction::m_zaero
double m_zaero[c_noOfAerogels]
z-positions of aerogel layers
Definition: ARICHReconstruction.h:119
Belle2::ARICHTrack::addPhoton
void addPhoton(ARICHPhoton photon)
Add ARICHPhoton to collection of photons.
Definition: ARICHTrack.h:190
Belle2::ARICHReconstruction::CherenkovPhoton
int CherenkovPhoton(TVector3 r, TVector3 rh, TVector3 &rf, TVector3 &dirf, double *refind, double *z, int n, int mirrorID=0)
Calculates the direction of photon emission.
Definition: ARICHReconstruction.cc:253
Belle2::ARICHPhoton::setChannel
void setChannel(int chn)
set channel (asic) of hit
Definition: ARICHPhoton.h:83
Belle2::ARICHReconstruction::likelihood2
int likelihood2(ARICHTrack &arichTrack, const StoreArray< ARICHHit > &arichHits, ARICHLikelihood &arichLikelihood)
Computes the value of identity likelihood function for different particle hypotheses.
Definition: ARICHReconstruction.cc:349
Belle2::ARICHReconstruction::correctEmissionPoint
void correctEmissionPoint(int tileID, double r)
correct mean emission point z position
Definition: ARICHReconstruction.cc:758
Belle2::ARICHReconstruction::setTrackAngleResolution
void setTrackAngleResolution(double aRes)
Sets track direction resolution (from tracking)
Definition: ARICHReconstruction.cc:688
Belle2::ARICHReconstruction::m_refractiveInd
double m_refractiveInd[c_noOfAerogels]
refractive indices of aerogel layers
Definition: ARICHReconstruction.h:118
Belle2::ARICHReconstruction::c_noOfHypotheses
static const int c_noOfHypotheses
Number of hypotheses to loop over.
Definition: ARICHReconstruction.h:98
Belle2::StoreArray
Accessor to arrays stored in the data store.
Definition: ECLMatchingPerformanceExpertModule.h:33
Belle2::ARICHReconstruction::m_nAerogelLayers
unsigned int m_nAerogelLayers
number of aerogel layers
Definition: ARICHReconstruction.h:117
Belle2::ARICHReconstruction::getMirrorNorm
TVector3 getMirrorNorm(int mirrorID)
Returns normal vector of the mirror plate with id mirrorID.
Definition: ARICHReconstruction.cc:745
Belle2::ARICHReconstruction::m_trackAngRes
double m_trackAngRes
track direction resolution (from tracking)
Definition: ARICHReconstruction.h:114
Belle2::Unit
The Unit class.
Definition: Unit.h:50
Belle2::ARICHReconstruction::m_arichgp
DBObjPtr< ARICHGeometryConfig > m_arichgp
geometry configuration parameters from the DB
Definition: ARICHReconstruction.h:102
Belle2::ARICHReconstruction::c_noOfAerogels
static const int c_noOfAerogels
Maximal number of aerogel layers to loop over.
Definition: ARICHReconstruction.h:99
Belle2::ARICHReconstruction::m_chnMap
DBObjPtr< ARICHChannelMapping > m_chnMap
map x,y channels to asic channels from the DB
Definition: ARICHReconstruction.h:105
Belle2::ARICHPhoton
Struct for ARICH reconstructed photon (hit related to track) information.
Definition: ARICHPhoton.h:27