153 double* refractiveInd,
double* z,
int n,
int opt)
164 double angmir = 0;
int section[2] = {0, 0};
166 unsigned tileID =
m_arichgp->getAerogelPlane().getAerogelTileID(r.X(), r.Y());
168 if (tileID == 0 && opt == 1)
return ROOT::Math::XYZVector();
170 int nmir =
m_arichgp->getMirrors().getNMirrors();
172 double dangle = 2 * M_PI / nmir;
173 angmir =
m_arichgp->getMirrors().getStartAngle() - dangle / 2.;
175 double trkangle = r.Phi() - angmir;
176 if (trkangle < 0) trkangle += 2 * M_PI;
177 if (trkangle > 2 * M_PI) trkangle -= 2 * M_PI;
179 section[1] = int(trkangle / dangle) + 1;
182 bool reflok =
false;
bool refl =
false;
183 double path = (z[0] - r.Z()) / dirf.Z();
185 for (
int a = 1; a <= n + 1 ; a++) {
186 double rind = refractiveInd[a] / refractiveInd[a - 1];
187 dirf = Refraction(dirf, rind);
188 if (dirf.R() == 0)
return ROOT::Math::XYZVector();
189 path = (z[a] - r.Z()) / dirf.Z();
190 ROOT::Math::XYZVector r0 = r;
191 if (a == n && opt == 1) {
192 if (
m_arichgp->getAerogelPlane().getAerogelTileID(r.X(), r.Y()) != tileID)
return ROOT::Math::XYZVector();
195 ROOT::Math::XYVector rxy(r);
197 if (a != n || nmir == 0)
continue;
198 double angle = rxy.Phi() - angmir;
199 if (angle < 0) angle += 2 * M_PI;
200 if (angle > 2 * M_PI) angle -= 2 * M_PI;
201 double dangle = 2 * M_PI / nmir;
202 section[0] = int(angle / dangle) + 1;
206 if (section[0] == section[1]) nrefl = 1;
207 for (
int k = 0; k < nrefl; k++) {
208 if (!
HitsMirror(r0, dirf, section[k]))
continue;
210 ROOT::Math::XYZVector mirnorm =
m_mirrorNorms[section[k] - 1];
211 double s = dirf.Dot(mirnorm);
212 double s1 = (mirpoint - r0).Dot(mirnorm);
213 r = r0 + s1 / s * dirf;
214 dirf = dirf - 2 * (dirf.Dot(mirnorm)) * mirnorm;
215 path = (z[a] - r.Z()) / dirf.Z();
223 if (!reflok && refl)
return ROOT::Math::XYZVector();
255 ROOT::Math::XYZVector& rf, ROOT::Math::XYZVector& dirf,
256 double* refractiveInd,
double* z,
int n,
int mirrorID)
275 static ROOT::Math::XYZVector norm(0, 0, 1);
277 double dwin =
m_arichgp->getHAPDGeometry().getWinThickness();
278 double refractiveInd0 =
m_arichgp->getHAPDGeometry().getWinRefIndex();
282 const double dmin = 0.0000001;
283 const int niter = 100;
284 ROOT::Math::XYZVector dirw;
285 ROOT::Math::XYZVector rh1 = rh - dwin * norm;
287 std::vector <ROOT::Math::XYZVector > rf0(n + 2);
291 std::vector <ROOT::Math::XYZVector > dirf0(n + 2);
304 for (
int iter = 0; iter < niter; iter++) {
308 if (iter == 0) dirf0[n] = (rh1 - rf0[1]).
Unit();
309 else dirf0[n] = (rh1 - rf0[n]).
Unit();
313 for (
int a = n - 1; a >= 0 ; a--) {
314 double rind = refractiveInd[a] / refractiveInd[a + 1];
315 dirf0[a] = Refraction(dirf0[a + 1], rind);
318 double path = (z[0] - r.Z()) / dirf0[0].Z();
319 double x1 = rf0[1].X();
320 double y1 = rf0[1].Y();
321 for (
int a = 0; a < n ; a++) {
322 rf0[a + 1] = rf0[a] + dirf0[a] * path;
323 path = (z[a + 1] - rf0[a + 1].Z()) / dirf0[a + 1].Z();
326 Refraction(dirf0[n], norm, refractiveInd0, dirw);
330 path = dwin / (dirw.Dot(norm));
331 rh1 = rh - dirw * path;
333 double d2 = (rf0[1].X() - x1) * (rf0[1].X() - x1) + (rf0[1].Y() - y1) * (rf0[1].Y() - y1);
335 if ((d2 < dmin) && (iter)) {
339 if (!
HitsMirror(rf0[n], dirf0[n], mirrorID))
return -1;
354 const unsigned int nPhotonHits = arichHits.getEntries();
367 double padSize =
m_arichgp->getHAPDGeometry().getPadSize();
368 int nMirSeg =
m_arichgp->getMirrors().getNMirrors();
369 double angmir =
m_arichgp->getMirrors().getStartAngle();
379 if (edir.Z() < 0)
return 0;
382 double thcResolution =
m_recPars->getThcResolution(momentum);
383 if (thcResolution < 0) thcResolution = 0.01;
385 double wideGaussFract = (
m_recPars->getParameters())[0];
386 double wideGaussSigma = (
m_recPars->getParameters())[1];
396 float nphot_scaling = 20.;
408 thetaCh[iHyp][iAerogel] = ExpectedCherenkovAngle(momentum,
414 if (pathLengthRadiator) pathLengthRadiator =
m_thickness[iAerogel] / pathLengthRadiator;
417 double dxx = pathLengthRadiator / double(nStep);
419 double nPhot =
m_n0[iAerogel] * sin(thetaCh[iHyp][iAerogel]) * sin(thetaCh[iHyp][iAerogel]) * dxx * nphot_scaling;
423 for (
int iepoint = 0; iepoint < nStep; iepoint++) {
425 ROOT::Math::XYZVector epoint = exit_point - (0.5 + iepoint) * dxx * edir;
427 unsigned int genPhot = nPhot * abs;
430 for (
unsigned int iPhoton = 0; iPhoton < genPhot; iPhoton++) {
431 double fi = 2 * M_PI * iPhoton / float(genPhot);
432 ROOT::Math::XYZVector adirf = setThetaPhi(thetaCh[iHyp][iAerogel], fi);
433 adirf = TransformFromFixed(edir) * adirf;
434 int ifi = int (fi * 20 / 2. / M_PI);
438 if (dposition.R() > 1.0) {nSig_wo_acc[iHyp][iAerogel][ifi] += 1; nSig_wo_accInt[iHyp][iAerogel] += 1;}
440 unsigned copyno =
m_arichgp->getDetectorPlane().pointSlotID(dposition.X(), dposition.Y());
441 if (!copyno)
continue;
443 if (
InsideDetector(dposition, copyno)) nSig_w_acc[iHyp][iAerogel] += 1;
448 for (
int ik = 0; ik < 20; ik++) {
449 nSig_wo_acc[iHyp][iAerogel][ik] /= nphot_scaling;
451 nSig_w_acc[iHyp][iAerogel] /= nphot_scaling;
452 nSig_wo_accInt[iHyp][iAerogel] /= nphot_scaling;
460 for (
int ik = 0; ik < 20; ik++) {
466 std::vector<double> bkgPars = {momentum /
sqrt(
p_mass[iHyp]*
p_mass[iHyp] + momentum * momentum), double(arichTrack.
hitsWindow())};
467 nBgr[iHyp] =
m_recPars->getExpectedBackgroundHits(bkgPars);
481 if (track_at_detector.Rho() > 85.0) {
482 double trackang = track_at_detector.Phi() - angmir;
483 if (trackang < 0) trackang += 2 * M_PI;
484 if (trackang > 2 * M_PI) trackang -= 2 * M_PI;
485 int section1 = int(trackang * nMirSeg / 2. / M_PI) + 1;
486 int section2 = section1 + 1;
487 if (section1 == nMirSeg) section2 = 1;
488 mirrors[1] = section1; mirrors[2] = section2;
494 for (
unsigned int iPhoton = 0; iPhoton < nPhotonHits; iPhoton++) {
497 int modID = h->getModule();
498 int channel = h->getChannel();
499 ROOT::Math::XYZVector hitpos =
m_arichgp->getMasterVolume().pointToLocal(h->getPosition());
500 bool bkgAdded =
false;
501 int nfoo = nDetPhotons;
502 for (
int iHyp = 0; iHyp <
c_noOfHypotheses; iHyp++) { esigi[iHyp] = 0; ebgri[iHyp] = 0;}
507 for (
int mirr = 0; mirr < refl; mirr++) {
516 if ((track_at_detector - virthitpos).
R() > 25.0)
continue;
525 double fi_cer_trk = 0.;
532 ROOT::Math::XYZVector edirr = arichTrack.
getDirection();
533 ROOT::Math::XYZVector photonDirection;
538 ROOT::Math::XYZVector dirch = TransformToFixed(edirr) * photonDirection;
539 double fi_cer = dirch.Phi();
540 double th_cer = dirch.Theta();
543 th_cer_all[iAerogel] = th_cer;
544 fi_cer_all[iAerogel] = fi_cer;
545 auto deltaPhi = dirch.Phi() - edir.Phi();
547 deltaPhi -= 2 * M_PI;
548 if (deltaPhi < -M_PI)
549 deltaPhi += 2 * M_PI;
550 fi_cer_trk = deltaPhi;
552 if (mirr == 0 && th_cer < 0.1) reflOK =
false;
554 if (th_cer > 0.5 || th_cer < 0.1)
continue;
557 if (nfoo == nDetPhotons) nDetPhotons++;
559 if (fi_cer < 0) fi_cer += 2 * M_PI;
563 fii = 2 * fi_mir - fi_cer - M_PI;
571 ROOT::Math::XYZVector photonDirection1 = setThetaPhi(thetaCh[iHyp][iAerogel], fi_cer);
572 photonDirection1 = TransformFromFixed(edirr) * photonDirection1;
573 int ifi = int (fi_cer * 20 / 2. / M_PI);
574 ROOT::Math::XYZVector detector_position;
579 ROOT::Math::XYZVector meanr = detector_position - epoint;
580 double path = meanr.R();
581 meanr = meanr.Unit();
583 double meanpath = (
m_recPars->getParameters())[2];
584 if (iAerogel == 1) meanpath = meanpath -
m_thickness[iAerogel];
586 double detector_sigma = thcResolution * meanpath / meanr.Z();
587 double wide_sigma = wideGaussSigma * path / meanr.Z();
589 double modphi =
m_arichgp->getDetectorPlane().getSlotPhi(modID);
591 double pad_fi = fii - modphi;
592 double dx = (detector_position - hitpos).
R();
593 double dr = (track_at_detector - detector_position).
R();
596 double normalizacija = nSig_wo_acc[iHyp][iAerogel][ifi] * padSize / (0.1 * M_PI * dr * meanr.Z());
597 weight[iHyp][iAerogel] = normalizacija;
598 weight_sum[iHyp] += weight[iHyp][iAerogel];
599 double integralMain = SquareInt(padSize, pad_fi, dx, detector_sigma) /
sqrt(2.);
600 double integralWide = SquareInt(padSize, pad_fi, dx, wide_sigma) /
sqrt(2.);
602 double exp = normalizacija * ((1 - wideGaussFract) * integralMain + wideGaussFract * integralWide);
604 sigExpArr[iHyp] += exp;
615 std::vector<double> pars = {momentum /
sqrt(
p_mass[iHyp]*
p_mass[iHyp] + momentum * momentum), double(arichTrack.
hitsWindow())};
616 ebgri[iHyp] +=
m_recPars->getBackgroundPerPad(th_cer_all[1], pars);
622 if (
m_storePhot && th_cer_all[1] > 0 && th_cer_all[1] < 0.6) {
626 if (weight_sum[iHyp] > 0) {
628 double emission_prob = weight[iHyp][iAerogel] / weight_sum[iHyp];
629 n_cos_theta_ch[iHyp] += emission_prob *
m_refractiveInd[iAerogel] * cos(th_cer_all[iAerogel]);
630 phi_ch[iHyp] += emission_prob * fi_cer_all[iAerogel];
634 n_cos_theta_ch[iHyp] = -99999.;
635 phi_ch[iHyp] = -99999.;
638 ARICHPhoton phot(iPhoton, th_cer_all[1], fi_cer_all[1], mirrors[mirr]);
644 phot.
setXY(hitpos.X(), hitpos.Y());
658 double expected = esigi[iHyp] + ebgri[iHyp];
659 if (bkgAdded) logL[iHyp] += expected + log(1 - exp(-expected));
669 exppho[iHyp] = nSig_w_acc[iHyp][
m_nAerogelLayers] * (1 - wideGaussFract) + wideGaussFract * 0.7 *
671 logL[iHyp] -= exppho[iHyp];
672 if (std::isnan(logL[iHyp]) || std::isinf(logL[iHyp])) {
673 B2WARNING(
"ARICHReconstruction: log likelihood value infinite! Flat background hit probability is " << ebgri[iHyp] <<
"!");
683 if ((thetaCh[0][0] > 0 || thetaCh[0][1] > 0) && nSig_w_acc[0][
m_nAerogelLayers] == 0) flag = 0;
686 arichLikelihood.
setValues(flag, logL, nDetPhotons, exppho);