Skip to content

Commit

Permalink
Merge pull request #162 from DUNE/kleykamp_issue_158
Browse files Browse the repository at this point in the history 
Kleykamp issue 158
  • Loading branch information
@jdkio
jdkio authored Sep 25, 2024
2 parents 8ce130c + eff42e0 commit 8fc2e62
Show file tree
Hide file tree
Showing 5 changed files with 54 additions and 31 deletions.
2 changes: 2 additions & 0 deletions src/Material.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -71,6 +71,8 @@ class Material {
fMaterialType = kAir;
} else {
fMaterialType = kUnknown;
std::cout<<"Made material with unknown fMaterialType, given density "<<density<<std::endl;
throw std::invalid_argument("Material error: do not know how to deal with density");
}

SetProperties();
Expand Down
13 changes: 6 additions & 7 deletions src/TMS_Event.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,6 @@ TMS_Event::TMS_Event() {
LightWeight = true;
}


static bool TMS_TrueParticle_NotWorthSaving(TMS_TrueParticle tp) {
if (tp.GetTrueVisibleEnergy() == 0 && !tp.IsPrimary()) return true;
// Don't worry about really low visible energy
Expand Down Expand Up @@ -1016,16 +1015,16 @@ double TMS_Event::GetMuonTrueTrackLength() {

std::vector<TLorentzVector> pos = (*it).GetPositionPoints();
int num = 0;
for (auto pnt = pos.begin(); pnt != pos.end(); ++pnt, ++num) {
for (auto pnt = pos.begin(); (pnt+1) != pos.end(); ++pnt, ++num) {
auto nextpnt = *(pnt+1);
TVector3 point1((*pnt).X(), (*pnt).Y(), (*pnt).Z()); //-200
TVector3 point2(nextpnt.X(), nextpnt.Y(), nextpnt.Z()); //-200
if (TMS_Geom::GetInstance().IsInsideTMS(point1) && TMS_Geom::GetInstance().IsInsideTMS(point2)) {
if ((point2-point1).Mag() > 100) {
continue;
}
double tracklength = TMS_Geom::GetInstance().GetTrackLength(point1, point2);
total += tracklength;
if ((point2-point1).Mag() > 100) {
continue;
}
double tracklength = TMS_Geom::GetInstance().GetTrackLength(point1, point2);
total += tracklength;
}
}
}
Expand Down
38 changes: 29 additions & 9 deletions src/TMS_Kalman.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -56,7 +56,9 @@ TMS_Kalman::TMS_Kalman(std::vector<TMS_Hit> &Candidates) :
}
}

const int N_LAYER_BACK = 10;
int N_LAYER_BACK = 10;
// Can't look back further than the first element
if (Candidates.size() < (unsigned)N_LAYER_BACK) N_LAYER_BACK = Candidates.size();

AverageXSlope = (Candidates[Candidates.size() - N_LAYER_BACK].GetRecoX() - Candidates.back().GetRecoX())/(Candidates[Candidates.size() - N_LAYER_BACK].GetZ() - Candidates.back().GetZ());
AverageYSlope = (Candidates.front().GetRecoY() - Candidates.back().GetRecoY())/(Candidates.front().GetZ() - Candidates.back().GetZ());
Expand Down Expand Up @@ -166,10 +168,18 @@ void TMS_Kalman::Predict(TMS_KalmanNode &Node) {
PreviousVec[4] = PreviousState.qp;


if ( (PreviousState.x < TMS_Const::TMS_Start[0]) || (PreviousState.x > TMS_Const::TMS_End[0]) ) // point outside x region
if ( (PreviousState.x < TMS_Const::TMS_Start[0]) || (PreviousState.x > TMS_Const::TMS_End[0]) ) { // point outside x region
std::cerr << "Predicted x value outside TMS: " << PreviousState.y << "\tTMS: [" << TMS_Const::TMS_Start[0] << ", "<< TMS_Const::TMS_End[0] << "]" << std::endl;
if ( (PreviousState.y < TMS_Const::TMS_Start[1]) || (PreviousState.y > TMS_Const::TMS_End[1]) ) // point outside y region
Node.PrintTrueReco();
PreviousState.Print();
CurrentState.Print();
}
if ( (PreviousState.y < TMS_Const::TMS_Start[1]) || (PreviousState.y > TMS_Const::TMS_End[1]) ) { // point outside y region
std::cerr << "Predicted y value outside TMS: " << PreviousState.y << "\tTMS: [" << TMS_Const::TMS_Start[1] << ", "<< TMS_Const::TMS_End[1] << "]" << std::endl;
Node.PrintTrueReco();
PreviousState.Print();
CurrentState.Print();
}


TVectorD UpdateVec = Transfer*(PreviousVec);
Expand Down Expand Up @@ -231,9 +241,15 @@ void TMS_Kalman::Predict(TMS_KalmanNode &Node) {
TotalLength += thickness;

// Update the Bethe Bloch calculator to use this material
Material matter(density);

Bethe.fMaterial = matter;
try {
Material matter(density);
Bethe.fMaterial = matter;
// Set the material for the multiple scattering
MSC.fMaterial = matter;
}
catch (std::invalid_argument const& ex) {
std::cout<<"Could not make a material using density "<<density<<", is position within tms bounds?"<<std::endl;
}

// Subtract off the energy loss for this material
if (ForwardFitting) en -= Bethe.Calc_dEdx(en)*density*thickness;
Expand All @@ -244,8 +260,6 @@ void TMS_Kalman::Predict(TMS_KalmanNode &Node) {
double en_var = Bethe.Calc_dEdx_Straggling(en)*density*thickness;
total_en_var += en_var*en_var;

// Set the material for the multiple scattering
MSC.fMaterial = matter;
// Calculate this before or after the energy subtraction/addition?
MSC.Calc_MS(en, thickness*density);

Expand Down Expand Up @@ -300,7 +314,7 @@ void TMS_Kalman::Predict(TMS_KalmanNode &Node) {
CovarianceMatrix(2,2) = 1.50;
CovarianceMatrix(3,3) = 2.50;
CovarianceMatrix(4,4) = 1.0;
std::cout << "Initialising covariance!" << std::endl;
if (Talk) std::cout << "Initialising covariance!" << std::endl;
}

Node.FillNoiseMatrix(); // Full the matrix for multiple scattering
Expand Down Expand Up @@ -345,10 +359,16 @@ void TMS_Kalman::Predict(TMS_KalmanNode &Node) {
if ( (CurrentState.x < TMS_Const::TMS_Start[0]) || (CurrentState.x > TMS_Const::TMS_End[0]) ) // point outside x region
{
std::cerr << "lol x value outside TMS: " << CurrentState.y << "\tTMS: [" << TMS_Const::TMS_Start[0] << ", "<< TMS_Const::TMS_End[0] << "]" << std::endl;
Node.PrintTrueReco();
PreviousState.Print();
CurrentState.Print();
}
if ( (CurrentState.y < TMS_Const::TMS_Start[1]) || (CurrentState.y > TMS_Const::TMS_End[1]) ) // point outside y region
{
std::cerr << "lol y value outside TMS: " << CurrentState.y << "\tTMS: [" << TMS_Const::TMS_Start[1] << ", "<< TMS_Const::TMS_End[1] << "]" << std::endl;
Node.PrintTrueReco();
PreviousState.Print();
CurrentState.Print();
}

Node.SetRecoXY(CurrentState);
Expand Down
30 changes: 16 additions & 14 deletions src/TMS_Reco.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -733,19 +733,20 @@ void TMS_TrackFinder::FindTracks(TMS_Event &event) {
KalmanFilter = TMS_Kalman(trk.Hits);
kalman_reco_mom = KalmanFilter.GetMomentum();

std::cout << "Kalman filter momentum: " << kalman_reco_mom << " MeV" << std::endl;
bool verbose_kalman = false;
if (verbose_kalman) std::cout << "Kalman filter momentum: " << kalman_reco_mom << " MeV" << std::endl;
trk.SetMomentum(kalman_reco_mom); // Fill the momentum of the TMS_Track obj

std::cout << "Kalman filter start pos : " << KalmanFilter.Start[0] << ", " << KalmanFilter.Start[1] << ", " << KalmanFilter.Start[2] << std::endl;
if (verbose_kalman) std::cout << "Kalman filter start pos : " << KalmanFilter.Start[0] << ", " << KalmanFilter.Start[1] << ", " << KalmanFilter.Start[2] << std::endl;
trk.SetStartPosition(KalmanFilter.Start[0], KalmanFilter.Start[1], KalmanFilter.Start[2]); // Fill the momentum of the TMS_Track obj

std::cout << "Kalman filter end pos : " << KalmanFilter.End[0] << ", " << KalmanFilter.End[1] << ", " << KalmanFilter.End[2] << std::endl;
if (verbose_kalman) std::cout << "Kalman filter end pos : " << KalmanFilter.End[0] << ", " << KalmanFilter.End[1] << ", " << KalmanFilter.End[2] << std::endl;
trk.SetEndPosition(KalmanFilter.End[0], KalmanFilter.End[1], KalmanFilter.End[2]); // Fill the momentum of the TMS_Track obj

std::cout << "Kalman filter start dir : " << KalmanFilter.StartDirection[0] << ", " << KalmanFilter.StartDirection[1] << ", " << KalmanFilter.StartDirection[2] << std::endl;
if (verbose_kalman) std::cout << "Kalman filter start dir : " << KalmanFilter.StartDirection[0] << ", " << KalmanFilter.StartDirection[1] << ", " << KalmanFilter.StartDirection[2] << std::endl;
trk.SetStartDirection(KalmanFilter.StartDirection[0], KalmanFilter.StartDirection[1], KalmanFilter.StartDirection[2]); // Fill the momentum of the TMS_Track obj

std::cout << "Kalman filter end dir : " << KalmanFilter.EndDirection[0] << ", " << KalmanFilter.EndDirection[1] << ", " << KalmanFilter.EndDirection[2] << std::endl;
if (verbose_kalman) std::cout << "Kalman filter end dir : " << KalmanFilter.EndDirection[0] << ", " << KalmanFilter.EndDirection[1] << ", " << KalmanFilter.EndDirection[2] << std::endl;
trk.SetEndDirection(KalmanFilter.EndDirection[0], KalmanFilter.EndDirection[1], KalmanFilter.EndDirection[2]); // Fill the momentum of the TMS_Track obj
trk.KalmanNodes = KalmanFilter.GetKalmanNodes(); // Fill the KalmanNodes of the TMS_Track

Expand Down Expand Up @@ -1260,13 +1261,13 @@ std::vector<TMS_Track> TMS_TrackFinder::TrackMatching3D() {
// Handling cases of two hits in same plane to be matched
// By adding a loop into these statements one could also take care of more than 2 hits in the same plane/with the same z coordinate

if (UTracks[itU].GetZ() == UTracks[itU - 1].GetZ()) {
if (itU > 0 && UTracks[itU].GetZ() == UTracks[itU - 1].GetZ()) {
CalculateRecoY(UTracks[itU - 1], UTracks[itU - 1], VTracks[itV]);
CalculateRecoX(UTracks[itU - 1], VTracks[itV], UTracks[itU - 1]);
(aTrack.Hits).push_back(UTracks[itU]); // This adds the original hit
if (itU > 0) --itU; // and this allows for the other hit then to be added with the next push_back statement
}
if (VTracks[itV].GetZ() == VTracks[itV - 1].GetZ()) {
if (itV > 0 && VTracks[itV].GetZ() == VTracks[itV - 1].GetZ()) {
CalculateRecoY(VTracks[itV - 1], UTracks[itU], VTracks[itV - 1]);
CalculateRecoX(UTracks[itU], VTracks[itV - 1], VTracks[itV - 1]);
(aTrack.Hits).push_back(VTracks[itV]);
Expand All @@ -1284,20 +1285,20 @@ std::vector<TMS_Track> TMS_TrackFinder::TrackMatching3D() {
CalculateRecoX(UTracks[itU], VTracks[itV], VTracks[itV]);

// Handling cases of two hits in same plane
if (UTracks[itU].GetZ() == UTracks[itU - 1].GetZ()) {
if (itU > 0 && UTracks[itU].GetZ() == UTracks[itU - 1].GetZ()) {
UTracks[itU - 1].SetRecoY(CompareY(UTracks[itU - 1], VTracks[itV], XTracks[itX]));
CalculateRecoX(UTracks[itU - 1], VTracks[itV], UTracks[itU - 1]);

(aTrack.Hits).push_back(UTracks[itU]); // This adds the original hit
if (itU > 0) --itU; // and this allows for the other hit then to be aded with the next push_back statement
}
if (VTracks[itV].GetZ() == VTracks[itV - 1].GetZ()) {
if (itV > 0 && VTracks[itV].GetZ() == VTracks[itV - 1].GetZ()) {
VTracks[itV - 1].SetRecoY(CompareY(UTracks[itU], VTracks[itV - 1], XTracks[itX]));
CalculateRecoX(UTracks[itU], VTracks[itV - 1], VTracks[itV - 1]);
(aTrack.Hits).push_back(VTracks[itV]);
if (itV > 0) --itV;
}
if (XTracks[itX].GetZ() == XTracks[itX - 1].GetZ()) {
if (itX > 0 && XTracks[itX].GetZ() == XTracks[itX - 1].GetZ()) {
XTracks[itX - 1].SetRecoY(XTracks[itX - 1].GetNotZ());

CalculateRecoX(UTracks[itU], VTracks[itV], XTracks[itX - 1]);
Expand All @@ -1322,13 +1323,13 @@ std::vector<TMS_Track> TMS_TrackFinder::TrackMatching3D() {
std::cout << "same in UV, not X" << std::endl;
std::cout << "Hit U: " << UTracks[itU].GetRecoX() << " | " << UTracks[itU].GetRecoY() << " | " << UTracks[itU].GetZ() << " / V: " << VTracks[itV].GetRecoX() << " | " << VTracks[itV].GetRecoY() << " | " << VTracks[itV].GetZ() << " / X: " << XTracks[itX].GetNotZ() << " | " << XTracks[itX].GetZ() << std::endl;
#endif
if (UTracks[itU].GetZ() == UTracks[itU - 1].GetZ()) {
if (itU > 0 && UTracks[itU].GetZ() == UTracks[itU - 1].GetZ()) {
CalculateRecoY(UTracks[itU - 1], UTracks[itU - 1], VTracks[itV]);
CalculateRecoX(UTracks[itU - 1], VTracks[itV], UTracks[itU - 1]);
(aTrack.Hits).push_back(UTracks[itU]); // This adds the original hit
if (itU > 0) --itU; // and this allows for the other hit then to be added with the next push_back statement
}
if (VTracks[itV].GetZ() == VTracks[itV - 1].GetZ()) {
if (itV > 0 && VTracks[itV].GetZ() == VTracks[itV - 1].GetZ()) {
CalculateRecoY(VTracks[itV - 1], UTracks[itU], VTracks[itV - 1]);
CalculateRecoX(UTracks[itU], VTracks[itV - 1], VTracks[itV - 1]);
(aTrack.Hits).push_back(VTracks[itV]);
Expand Down Expand Up @@ -3794,8 +3795,9 @@ void TMS_TrackFinder::Accumulate(double xhit, double zhit) {
std::cout << "cbin: " << c_bin << std::endl;
}
*/
// Fill the accumulator
Accumulator[i][c_bin]++;
// Fill the accumulator, but only within bounds
// We don't care about lines outside of bounds
if (i >= 0 && c_bin >= 0 && i < nSlope && c_bin < nIntercept) Accumulator[i][c_bin]++;
}
}

2 changes: 1 addition & 1 deletion src/TMS_TreeWriter.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -234,7 +234,7 @@ void TMS_TreeWriter::MakeBranches() {
Reco_Tree->Branch("KalmanPos", RecoTrackKalmanPos, "TrackHitPos[nTracks][200][3]/F");
Reco_Tree->Branch("KalmanTruePos", RecoTrackKalmanTruePos, "TrackHitTruePos[nTracks][200][3]/F");
Reco_Tree->Branch("StartDirection", RecoTrackStartDirection, "StartDirection[nTracks][3]/F");
Reco_Tree->Branch("EndDirectioN", RecoTrackEndDirection, "EndDirection[nTracks][3]/F");
Reco_Tree->Branch("EndDirection", RecoTrackEndDirection, "EndDirection[nTracks][3]/F");
Reco_Tree->Branch("StartPos", RecoTrackStartPos, "StartPos[nTracks][3]/F");
Reco_Tree->Branch("EndPos", RecoTrackEndPos, "EndPos[nTracks][3]/F");
Reco_Tree->Branch("EnergyRange", RecoTrackEnergyRange, "EnergyRange[nTracks]/F");
Expand Down

0 comments on commit 8fc2e62

Please sign in to comment.