From b9c228ec9964ec626e2a34ab94e444496f3a8c54 Mon Sep 17 00:00:00 2001
From: Alexander Kiselev <ayk@bnl.gov>
Date: Fri, 19 Nov 2021 02:11:13 -0500
Subject: [PATCH] Substantial clean up in the IRT codes
---
drich-testIRT.py | 4 +-
evaluation/source/reader.cc | 2 +
include/ChargedParticle.h | 2 +-
include/ChargedParticleStep.h | 5 +-
include/CherenkovDetector.h | 50 ++++++++++-
include/CherenkovRadiator.h | 9 +-
include/OpticalPhoton.h | 6 +-
include/ParametricSurface.h | 7 +-
include/RadiatorHistory.h | 35 ++++++++
scripts/evaluation.C | 41 ++++++---
source/ChargedParticle.cc | 164 ++++++++++++----------------------
source/FlatSurface.cc | 5 +-
source/SphericalSurface.cc | 5 +-
13 files changed, 195 insertions(+), 140 deletions(-)
diff --git a/drich-testIRT.py b/drich-testIRT.py
index cba1e4a..12baf5d 100644
--- a/drich-testIRT.py
+++ b/drich-testIRT.py
@@ -32,8 +32,8 @@ qe_data = [
]
radiators = [
- "Aerogel zbins=1 smearing=gaussian 1mrad rindex=1.0189",
- "GasVolume zbins=1 smearing=gaussian 1mrad rindex=1.00080"
+ "Aerogel zbins=1 smearing=gaussian 2mrad rindex=1.0190",
+ "GasVolume zbins=10 smearing=gaussian 1mrad rindex=1.00076"
]
podioinput = PodioInput(
diff --git a/evaluation/source/reader.cc b/evaluation/source/reader.cc
index f6089f6..82a6e8b 100644
--- a/evaluation/source/reader.cc
+++ b/evaluation/source/reader.cc
@@ -20,6 +20,7 @@
int main(int argc, char** argv)
{
+#if _BACK_
// Command line "parser";
if (argc != 3 && argc != 4) {
printf("usage: %s <root-data-file> <root-config-file> [DNAME]\n", argv[0]);
@@ -188,6 +189,7 @@ int main(int argc, char** argv)
//auto cv = new TCanvas("cv", "", 800, 600);
//cv->cd(1); np->Draw();
+#endif
return 0;
} // main()
diff --git a/include/ChargedParticle.h b/include/ChargedParticle.h
index cd96d34..9f41fbf 100644
--- a/include/ChargedParticle.h
+++ b/include/ChargedParticle.h
@@ -64,7 +64,7 @@ class ChargedParticle: public TransientParticle {
};
// Single particle case for now;
- void PIDReconstruction(CherenkovPID &pid, std::vector<OpticalPhoton*> *photons = 0);
+ void PIDReconstruction(CherenkovPID &pid);
private:
// Optical photons produced elsewhere;
diff --git a/include/ChargedParticleStep.h b/include/ChargedParticleStep.h
index 48540ae..2eeebbd 100644
--- a/include/ChargedParticleStep.h
+++ b/include/ChargedParticleStep.h
@@ -12,8 +12,9 @@ class ChargedParticleStep: public TObject {
m_Position(position), m_Momentum(momentum)/*, m_Length(length)*/ {};
~ChargedParticleStep() {};
- inline const TVector3 &GetPosition( void ) const { return m_Position; };
- inline const TVector3 &GetMomentum( void ) const { return m_Momentum; };
+ inline const TVector3 &GetPosition ( void ) const { return m_Position; };
+ inline const TVector3 &GetMomentum ( void ) const { return m_Momentum; };
+ inline TVector3 GetDirection( void ) const { return m_Momentum.Unit(); };
// Nodes along the particle trajectory;
private:
diff --git a/include/CherenkovDetector.h b/include/CherenkovDetector.h
index 941e112..9d7a505 100644
--- a/include/CherenkovDetector.h
+++ b/include/CherenkovDetector.h
@@ -16,7 +16,7 @@ class G4LogicalVolume;
class CherenkovDetector: public TObject {
public:
CherenkovDetector(const char *name = 0): /*m_ContainerVolume(0),*/ m_Name(name ? name : ""),
- m_ReadoutCellMask(0x0) {};
+ m_ReadoutCellMask(0x0)/*, m_SectorBoundaryOffset(0.0)*/ {};
~CherenkovDetector() {};
void AddOpticalBoundary(unsigned sector, OpticalBoundary *boundary) {
@@ -69,15 +69,61 @@ class CherenkovDetector: public TObject {
void SetReadoutCellMask(uint64_t mask) { m_ReadoutCellMask = mask; };
inline uint64_t GetReadoutCellMask( void ) const { return m_ReadoutCellMask; };
+ // FIXME: not at all clean (uses implicit phase assumptions);
+ unsigned GetSector(const TVector3 &pt) {
+ // FIXME: may require tuning for a dual-mirror setup;
+ unsigned nSectors = _m_OpticalBoundaries.size();
+
+ // Either a single "sector" or sector structure not define yet -> return 0;
+ if (nSectors <= 1) return 0;
+
+ // FIXME: this offset is only defined by the way Chris positions sector #0;
+ double bin = 2*M_PI/nSectors, offset = -bin/2;
+
+ return (unsigned)floor((pt.Phi() + 4*M_PI - offset)/bin) % nSectors;
+ //return 0;
+ };
+
+ // FIXME: get rid of the second argument here;
+ CherenkovRadiator *GuessRadiator(const TVector3 &x0, const TVector3 &n0) {
+ // FIXME: may want to do a better check;
+ if (_m_OpticalBoundaries.empty()) return 0;
+
+ // Determine sector (in EIC DRICH terminology);
+ unsigned isec = GetSector(x0);
+
+ // Now loop through all radiators, and check boundaries in this sector;
+ for(auto rptr: _m_Radiators) {
+ const auto radiator = rptr.second;
+
+ // Front and rear surfaces for this particular sector;
+ auto s1 = radiator->GetFrontSide(isec);
+ auto s2 = radiator->GetRearSide (isec);
+
+ TVector3 from, to;
+ // Go backwards and ignore surface orientation mismatch;
+ bool b1 = s1->GetCrossing(x0, -1*n0, &from, false);
+ bool b2 = s2->GetCrossing(x0, n0, &to);
+ if (!b1 || !b2) continue;
+
+ if ((x0 - from).Dot(to - x0) > 0.0) return radiator;
+ } //for radiator
+
+ // Seemingly this 3D point does not belong to any radiator;
+ return 0;
+ };
+
private:
TString m_Name;
// This is needed for dd4hep cell index decoding;
uint64_t m_ReadoutCellMask;
+ //+double m_SectorBoundaryOffset;
+
std::map<TString, CherenkovRadiator*> _m_Radiators;
//+std::vector<CherenkovMirrorGroup*> m_MirrorGroups;
- ClassDef(CherenkovDetector, 4);
+ ClassDef(CherenkovDetector, 5);
};
#endif
diff --git a/include/CherenkovRadiator.h b/include/CherenkovRadiator.h
index ef1c7ab..8ee9279 100644
--- a/include/CherenkovRadiator.h
+++ b/include/CherenkovRadiator.h
@@ -70,9 +70,9 @@ class CherenkovRadiator: public TObject {
bool UseGaussianSmearing( void ) const { return m_GaussianSmearing; };
// FIXME: memory leak;
- void ResetLocations( void ) { _m_Locations.clear(); }
- void AddLocation(unsigned sector, const TVector3 &x, const TVector3 &p) {
- _m_Locations[sector].push_back(std::make_pair(x, p));
+ void ResetLocations( void ) { m_Locations.clear(); }
+ void AddLocation(/*unsigned sector,*/ const TVector3 &x, const TVector3 &n) {
+ m_Locations/*[sector]*/.push_back(std::make_pair(x, n));
};
// Transient variables for the analysis script convenience;
@@ -85,7 +85,8 @@ class CherenkovRadiator: public TObject {
// This is a hack for now;
double m_Smearing; //!
bool m_GaussianSmearing; //!
- std::map<unsigned, std::vector<std::pair<TVector3, TVector3>>> _m_Locations; //!
+ //std::map<unsigned, std::vector<std::pair<TVector3, TVector3>>> _m_Locations; //!
+ std::vector<std::pair<TVector3, TVector3>> m_Locations; //!
std::vector<std::pair<double, double>> m_ri_lookup_table; //!
diff --git a/include/OpticalPhoton.h b/include/OpticalPhoton.h
index 10d806d..eb595b8 100644
--- a/include/OpticalPhoton.h
+++ b/include/OpticalPhoton.h
@@ -1,4 +1,5 @@
+#include <set>
#include <map>
#include <vector>
@@ -17,7 +18,7 @@ class OpticalPhoton: public TransientParticle {
public:
OpticalPhoton(): TransientParticle(0), m_VertexRefractiveIndex(0.0), m_PhotonDetector(0),
m_VolumeCopy(0), m_DetectionTime(0.0),
- m_Detected(false), m_Selected(false) {};
+ m_Detected(false) {};//, m_Selected(false) {};
~OpticalPhoton() {};
inline bool IsCharged( void ) const { return false; };
@@ -84,7 +85,8 @@ class OpticalPhoton: public TransientParticle {
// Transient variables for some convenience in an analysis script;
public:
//inline bool WasSelected( void ) const { return m_Selected; };
- bool m_Selected; //!
+ //bool m_Selected; //!
+ std::set<std::pair<unsigned, CherenkovRadiator*>> _m_Selected; //!
std::map<CherenkovRadiator*, VectorPDF> _m_PDF; //!
ClassDef(OpticalPhoton, 3);
diff --git a/include/ParametricSurface.h b/include/ParametricSurface.h
index 3b72cce..00b6112 100644
--- a/include/ParametricSurface.h
+++ b/include/ParametricSurface.h
@@ -43,7 +43,8 @@ class ParametricSurface: public TObject {
virtual double GetDistance(const TVector3 &xx) const = 0;
// Crossing with the straight line defined by {x0,n0};
- virtual bool GetCrossing(const TVector3 &x0, const TVector3 &n0, TVector3 *crs) const = 0;
+ virtual bool GetCrossing(const TVector3 &x0, const TVector3 &n0, TVector3 *crs,
+ bool check_normal = true) const = 0;
// Introduce only the transformations needed for the current task;
virtual ParametricSurface *_Clone(double angle, const TVector3 &axis) const = 0;
@@ -80,7 +81,7 @@ class SphericalSurface: public ParametricSurface {
};
double GetRadius( void ) const { return m_Radius; };
- bool GetCrossing(const TVector3 &x0, const TVector3 &n0, TVector3 *crs) const;
+ bool GetCrossing(const TVector3 &x0, const TVector3 &n0, TVector3 *crs, bool check_normal = true) const;
double GetDistance(const TVector3 &xx) const;
void SetConvex( void ) { m_Concave = false; };
@@ -139,7 +140,7 @@ class FlatSurface: public ParametricSurface, public LocalCoordinatesXY {
return (xx - GetCenter()).Dot(m_Ny);
};
- bool GetCrossing(const TVector3 &x0, const TVector3 &n0, TVector3 *crs) const;
+ bool GetCrossing(const TVector3 &x0, const TVector3 &n0, TVector3 *crs, bool check_normal = true) const;
double GetDistance(const TVector3 &xx) const;
ParametricSurface *_Clone(double angle, const TVector3 &axis) const {
auto copy = new FlatSurface(*this);
diff --git a/include/RadiatorHistory.h b/include/RadiatorHistory.h
index 33ad9cc..bc74b83 100644
--- a/include/RadiatorHistory.h
+++ b/include/RadiatorHistory.h
@@ -1,4 +1,5 @@
+#include <map>
#include <vector>
#include <TObject.h>
@@ -32,11 +33,45 @@ class RadiatorHistory: public TObject {
};
inline unsigned StepCount( void ) { return m_Steps.size(); };
+ void AddStepBufferPoint(double time, const TVector3 &x) {
+ // Will be in ascending order of time;
+ m_StepBuffer[time] = x;
+ };
+
+ void CalculateSteps(bool dirty = true) {
+ // FIXME: memory leak;
+ m_Steps.clear();
+
+ // Add origin to the top of the buffer; FIXME: can times be negative
+ // as provided by GEANT?;
+ if (m_StepBuffer.size() == 1) {
+ if (!dirty) return;
+
+ m_StepBuffer[0.0] = TVector3(0,0,0);
+ } //if
+
+ // FIXME: efficiency sucks here;
+ std::vector<TVector3> buffer;
+ for(auto entry: m_StepBuffer)
+ buffer.push_back(entry.second);
+
+ for(unsigned iq=1; iq<buffer.size(); iq++) {
+ m_Steps.push_back(new ChargedParticleStep(buffer[iq-1],
+ (buffer[iq] - buffer[iq-1]).Unit()));
+ // Well, would not hurt to add the last point in the same direction;
+ if (iq == buffer.size()-1)
+ m_Steps.push_back(new ChargedParticleStep(buffer[iq],
+ (buffer[iq] - buffer[iq-1]).Unit()));
+ } //for iq
+ };
+
// Charged particle trajectory and optical photons generated in this radiator;
private:
std::vector<OpticalPhoton*> m_Photons;
std::vector<ChargedParticleStep*> m_Steps;
+ std::map<double, TVector3> m_StepBuffer; //!
+
ClassDef(RadiatorHistory, 1);
};
diff --git a/scripts/evaluation.C b/scripts/evaluation.C
index b872b8a..372d973 100644
--- a/scripts/evaluation.C
+++ b/scripts/evaluation.C
@@ -1,11 +1,11 @@
-
-void evaluation(const char *dfname)
-{
#define _DETECTOR_ "DRICH"
+//#define _DETECTOR_ "ERICH"
#define _AEROGEL_
- //#define _USE_RECONSTRUCTED_TRACKS_
+//#define _USE_RECONSTRUCTED_TRACKS_
+void evaluation(const char *dfname)
+{
// .root file with event tree;
auto fdata = new TFile(dfname);
if (!fdata) {
@@ -21,13 +21,15 @@ void evaluation(const char *dfname)
auto np = new TH1D("np", "Photon count", 50, 0, 50);
#ifdef _AEROGEL_
unsigned id = 0;
- auto th = new TH1D("th", "Cherenkov theta", 100, 180, 200);
- auto ri = new TH1D("ri", "Refractive Index - 1.0", 100, 0.017, 0.022);
+ auto th = new TH1D("th", "Cherenkov theta", 50, 180, 200);
+ auto ri = new TH1D("ri", "Refractive Index - 1.0", 50, 0.018, 0.020);
+ auto dt = new TH1D("dt", "Cherenkov theta diff", 50, -5, 5);
+ //auto dt = new TH1D("dt", "Cherenkov theta diff", 50, -1, 1);
#else
unsigned id = 1;
- auto th = new TH1D("th", "Cherenkov theta", 100, 33, 43);
- auto ri = new TH1D("ri", "Refractive Index - 1.0", 100, 0.00074, 0.00079);
- //auto th = new TH1D("th", "Cherenkov theta", 100, 33, 53);
+ auto th = new TH1D("th", "Cherenkov theta", 50, 35, 39);
+ auto ri = new TH1D("ri", "Refractive Index - 1.0", 50, 0.00075, 0.00077);
+ auto dt = new TH1D("dt", "Cherenkov theta diff", 50, -2, 2);
#endif
// Use MC truth particles for a "main" loop;
@@ -82,6 +84,9 @@ void evaluation(const char *dfname)
#endif
if (!cherenkov) continue;
+ double pp = mctrack.ps.mag(), m = mctrack.mass;
+ //printf("%f %f\n", mctrack.mass, mctrack.ps.mag());
+
// Loop through all of the mass hypotheses available for this reconstructed track;
{
const eic::CherenkovPdgHypothesis *best = 0;
@@ -108,18 +113,26 @@ void evaluation(const char *dfname)
}
// This assumes of course that at least one radiator was requested in juggler;
- th->Fill(1000* (*angles)[id].theta);
- ri->Fill( (*angles)[id].rindex - 1.0);
- printf("%f\n", (*angles)[id].rindex - 1.0);
+ {
+ double rindex = (*angles)[id].rindex, theta = (*angles)[id].theta;
+ double argument = sqrt(pp*pp + m*m)/(rindex*pp);
+ double thp = fabs(argument) <= 1.0 ? acos(argument) : theta;//sqrt(pp*pp + m*m)/(rindex*pp));
+
+ th->Fill(1000* theta);
+ dt->Fill(1000* (theta - thp));
+ ri->Fill( rindex - 1.0);
+ printf("%f\n", rindex - 1.0);
+ }
} //for track
} //for ev
printf("%3d (%3d) false out of %lld\n", false_assignment_stat[0],
false_assignment_stat[1], t->GetEntries());
- auto cv = new TCanvas("cv", "", 1500, 600);
- cv->Divide(3, 1);
+ auto cv = new TCanvas("cv", "", 1500, 500);
+ cv->Divide(4, 1);
cv->cd(1); np->Draw(); np->Fit("gaus");
cv->cd(2); th->Draw(); th->Fit("gaus");
cv->cd(3); ri->Draw(); ri->Fit("gaus");
+ cv->cd(4); dt->Draw(); dt->Fit("gaus");
} // evaluation()
diff --git a/source/ChargedParticle.cc b/source/ChargedParticle.cc
index e2e7ca2..6f7df0e 100644
--- a/source/ChargedParticle.cc
+++ b/source/ChargedParticle.cc
@@ -11,30 +11,17 @@
// -------------------------------------------------------------------------------------
-void ChargedParticle::PIDReconstruction(CherenkovPID &pid, std::vector<OpticalPhoton*> *photons)
+void ChargedParticle::PIDReconstruction(CherenkovPID &pid)
{
-#if 0
- if (!photons) {
- photons = new std::vector<OpticalPhoton*>();
+ std::vector<OpticalPhoton*> photons;
- for(auto rhistory: GetRadiatorHistory())
- for(auto photon: GetHistory(rhistory)->Photons())
- photons->push_back(photon);
- } //if
-#endif
-#if 1//_OLD_
- // Reset individual photon "selected" flags;
- if (photons)
- for(auto photon: *photons)
- photon->m_Selected = false;
- else
- for(auto rhistory: GetRadiatorHistory())
- for(auto photon: GetHistory(rhistory)->Photons())
- photon->m_Selected = false;
-#endif
+ // Mix all photons together;
+ for(auto rhistory: GetRadiatorHistory())
+ for(auto photon: GetHistory(rhistory)->Photons())
+ photons.push_back(photon);
- unsigned qzdim = 0;
- TVector3 momentum;
+ for(auto photon: photons)
+ photon->_m_Selected.clear();// = false;
// Loop through all of the photons recorded in all radiators; apply IRT on a fixed grid
// of emission vertex locations and build kind of a PDF out of that to calculate weights;
@@ -43,74 +30,49 @@ void ChargedParticle::PIDReconstruction(CherenkovPID &pid, std::vector<OpticalPh
// hypothesis would evaluate each of the detected photons on the aerogel and gas grids
// independently and build a sum of their PDFs without much thinking about outlier photon
// rejection, average theta calculation and such;
- for(auto rhistory: GetRadiatorHistory()) {
- auto radiator = GetRadiator(rhistory);
- // FIXME: error message;
- //if (radiator->m_Locations.size() < 2) continue;//return;
- //if (radiator->_m_Locations.empty() || radiator->_m_Locationssize() < 2) continue;//return;
+ for(auto photon: photons) {
+ if (!photon->WasDetected()) continue;
- //const unsigned zdim = radiator->m_Locations.size()-1;//radiator->GetTrajectoryBinCount();
-#if 1//_TODAY_
- //if (radiator->m_Locations.size()) {
- for(auto photon: (photons ? *photons : GetHistory(rhistory)->Photons())) {
- if (!photon->WasDetected()) continue;
-
- auto pd = photon->GetPhotonDetector();
- const auto irt = pd->GetIRT(photon->GetVolumeCopy());//vcopy);
- if (!irt) {//pd->GetIRT(photon->GetVolumeCopy())) {
- printf("No photosensor with this cellID found!\n");
- continue;
- } //if
-
- unsigned isec = irt->GetSector();
-
- if (radiator->_m_Locations.find(isec) == radiator->_m_Locations.end() ||
- radiator->_m_Locations[isec].size() < 2)
- continue;
-
- const unsigned zdim = radiator->_m_Locations[isec].size()-1;//radiator->GetTrajectoryBinCount();
- if (zdim) {
- qzdim = zdim;
- momentum = 0.5*(radiator->_m_Locations[isec][0].second +
- radiator->_m_Locations[isec][radiator->_m_Locations[isec].size()-1].second);
- } //if
-
- TVector3 phx = photon->GetDetectionPosition();
+ auto pd = photon->GetPhotonDetector();
+ const auto irt = pd->GetIRT(photon->GetVolumeCopy());
+ if (!irt) {
+ printf("No photosensor with this cellID found!\n");
+ continue;
+ } //if
- {
- bool all_converged = true;
- IRTSolution solutions[zdim+1];
-
- for(unsigned iq=0; iq<zdim+1; iq++) {
- //printf("--> %ld %ld\n", photon->GetVolumeCopy(), (unsigned long)(pd->GetIRT(photon->GetVolumeCopy())));
- //printf("--> %ld %ld\n", photon->GetVolumeCopy(), (photon->GetVolumeCopy() >> 8) & 0x7);
- auto &solution = solutions[iq] =
- //pd->GetIRT(photon->GetVolumeCopy())->Solve(radiator->m_Locations[iq].first,
- irt->Solve(radiator->_m_Locations[isec][iq].first,
- radiator->_m_Locations[isec][iq].second.Unit(),
- // FIXME: give beam line as a parameter;
- phx, TVector3(0,0,1), false);
- if (!solution.Converged()) {
- all_converged = false;
- break;
- } //if
- } //for iq
+ for(auto rhistory: GetRadiatorHistory()) {
+ auto radiator = GetRadiator(rhistory);
+ unsigned zdim = radiator->GetTrajectoryBinCount();
+ if (radiator->m_Locations.size() != zdim+1) continue;
+
+ TVector3 phx = photon->GetDetectionPosition();
+
+ {
+ bool all_converged = true;
+ IRTSolution solutions[zdim+1];
+
+ for(unsigned iq=0; iq<zdim+1; iq++) {
+ auto &solution = solutions[iq] = irt->Solve(radiator->m_Locations[iq].first,
+ // FIXME: give beam line as a parameter;
+ radiator->m_Locations[iq].second.Unit(), phx, TVector3(0,0,1), false);
+ if (!solution.Converged()) {
+ all_converged = false;
+ break;
+ } //if
+ } //for iq
+
+ if (!all_converged) continue;
+
+ for(unsigned iq=0; iq<zdim; iq++) {
+ auto &s0 = solutions[iq], &s1 = solutions[iq+1];
- if (!all_converged) continue;
-
- for(unsigned iq=0; iq<zdim; iq++) {
- auto &s0 = solutions[iq], &s1 = solutions[iq+1];
-
- //printf(" converged: %d %d\n", s0.Converged(), s1.Converged());
- // NB: y0 & y1 values do not matter; what matters is that they were equidistant
- // in the previous loop; FIXME: add some smearing later;
- photon->_m_PDF[radiator].AddMember(new UniformPDF(s0.GetTheta(), s1.GetTheta(), 1.0));
- } //for iq
- }
- } //for photon
- //} //if
-#endif
- } //for rhistory
+ // NB: y0 & y1 values do not matter; what matters is that they were equidistant
+ // in the previous loop; FIXME: add some smearing later;
+ photon->_m_PDF[radiator].AddMember(new UniformPDF(s0.GetTheta(), s1.GetTheta(), 1.0));
+ } //for iq
+ }
+ } //for rhistory
+ } //for photon
// And now that IRT is performed on a grid of possible emission vertex locations,
// populate the requested mass hypothesis array;
@@ -123,20 +85,12 @@ void ChargedParticle::PIDReconstruction(CherenkovPID &pid, std::vector<OpticalPh
for(auto rhistory: GetRadiatorHistory()) {
auto radiator = GetRadiator(rhistory);
- //const unsigned zdim = radiator->GetTrajectoryBinCount();
-#if 1//_TODAY_
- //unsigned isec = irt->GetSector();
-
- //auto mloc = radiator->_m_Locations.begin();
- //+++const unsigned zdim = radiator->m_Locations.size()-1;//radiator->GetTrajectoryBinCount();
- //const unsigned zdim = (*mloc).size()-1;//radiator->GetTrajectoryBinCount();
-
- //if (!radiator->m_Locations.size()) continue;
-
- //+TVector3 p = 0.5*(radiator->m_Locations[0].second + radiator->m_Locations[radiator->m_Locations.size()-1].second);
- //TVector3 p = 0.5*((*radiator->_m_Locations.begin())[0].second + radiator->m_Locations[radiator->m_Locations.size()-1].second);
+ unsigned zdim = radiator->GetTrajectoryBinCount();
+ if (radiator->m_Locations.size() != zdim+1) continue;
- double pp = momentum.Mag(), arg = sqrt(pp*pp + m*m)/(radiator->m_AverageRefractiveIndex*pp);
+ // Asume this estimate is good enough;
+ double pp = radiator->m_Locations[0].second.Mag();
+ double arg = sqrt(pp*pp + m*m)/(radiator->m_AverageRefractiveIndex*pp);
// Threshold check; FIXME: do it better?;
if (fabs(arg) > 1.0) continue;
@@ -146,24 +100,22 @@ void ChargedParticle::PIDReconstruction(CherenkovPID &pid, std::vector<OpticalPh
// +/-3 sigma; FIXME: do it better later;
if (gaussian) dth *= 3;
- for(auto photon: (photons ? *photons : GetHistory(rhistory)->Photons())) {
+ for(auto photon: photons) {
if (!photon->WasDetected()) continue;
auto pdf = &photon->_m_PDF[radiator];
// FIXME: unreadable;
- //printf("-> %7.2f\n", pdf->GetRangeIntegral(theta - dth, theta + dth));
auto within_range = float(pdf->GetWithinRangeCount(theta, dth));
- if (within_range) photon->m_Selected = true;
+ //if (within_range) photon->m_Selected = true;
+ if (within_range) photon->_m_Selected.insert(std::make_pair(ih, radiator));// = true;
- hypothesis->IncrementWeight(radiator, within_range/qzdim,
+ hypothesis->IncrementWeight(radiator, within_range/zdim,
(dth ? (gaussian ? pdf->GetGaussianIntegral(theta,dth) :
pdf->GetRangeIntegral(theta - dth, theta + dth)) :
- pdf->GetValue(theta))/qzdim);
- //printf("@W@ %2d -> %7.2f %7.2f\n", ih, m, 1000*theta);
+ pdf->GetValue(theta))/zdim);
} //for photon
}
-#endif
} //for rhistory
} //for ih
} // ChargedParticle::PIDReconstruction()
diff --git a/source/FlatSurface.cc b/source/FlatSurface.cc
index 492a377..edf204a 100644
--- a/source/FlatSurface.cc
+++ b/source/FlatSurface.cc
@@ -3,7 +3,8 @@
// -------------------------------------------------------------------------------------
-bool FlatSurface::GetCrossing(const TVector3 &x0, const TVector3 &n0, TVector3 *crs) const
+bool FlatSurface::GetCrossing(const TVector3 &x0, const TVector3 &n0, TVector3 *crs,
+ bool check_normal) const
{
double npl = n0.Dot(m_Nz);
@@ -11,7 +12,7 @@ bool FlatSurface::GetCrossing(const TVector3 &x0, const TVector3 &n0, TVector3 *
if (!crs || !npl) return false;
// May want to check that approach the surface from a right direction;
- if (/*check_normal &&*/ n0.Dot(m_Nz) >= 0.0) return false;
+ if (check_normal && n0.Dot(m_Nz) >= 0.0) return false;
{
TVector3 dx = x0 - GetCenter();
diff --git a/source/SphericalSurface.cc b/source/SphericalSurface.cc
index b8d323f..e78e755 100644
--- a/source/SphericalSurface.cc
+++ b/source/SphericalSurface.cc
@@ -3,7 +3,8 @@
// -------------------------------------------------------------------------------------
-bool SphericalSurface::GetCrossing(const TVector3 &x0, const TVector3 &n0, TVector3 *crs) const
+bool SphericalSurface::GetCrossing(const TVector3 &x0, const TVector3 &n0, TVector3 *crs,
+ bool check_normal) const
{
TVector3 dx = x0 - GetCenter();
double a = 1.0, b = 2*dx.Dot(n0), c = dx.Mag2() - pow(m_Radius, 2);
@@ -20,7 +21,7 @@ bool SphericalSurface::GetCrossing(const TVector3 &x0, const TVector3 &n0, TVect
if (t < 0.0) continue;
*crs = x0 + t*n0;
- if (n0.Dot(GetNormal(*crs)) >= 0.0) continue;
+ if (check_normal && n0.Dot(GetNormal(*crs)) >= 0.0) continue;
return true;
} //for iq
--
GitLab