TRestDetectorHitsEvent.cxx

/*************************************************************************
 * This file is part of the REST software framework.                     *
 *                                                                       *
 * Copyright (C) 2016 GIFNA/TREX (University of Zaragoza)                *
 * For more information see http://gifna.unizar.es/trex                  *
 *                                                                       *
 * REST is free software: you can redistribute it and/or modify          *
 * it under the terms of the GNU General Public License as published by  *
 * the Free Software Foundation, either version 3 of the License, or     *
 * (at your option) any later version.                                   *
 *                                                                       *
 * REST is distributed in the hope that it will be useful,               *
 * but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the          *
 * GNU General Public License for more details.                          *
 *                                                                       *
 * You should have a copy of the GNU General Public License along with   *
 * REST in $REST_PATH/LICENSE.                                           *
 * If not, see http://www.gnu.org/licenses/.                             *
 * For the list of contributors see $REST_PATH/CREDITS.                  *
 *************************************************************************/
 
#include "TRestDetectorHitsEvent.h"
 
#include <array>
 
#include "TCanvas.h"
#include "TColor.h"
#include "TH3F.h"
#include "TPolyLine3D.h"
#include "TPolyMarker3D.h"
#include "TRandom.h"
#include "TRestStringHelper.h"
#include "TRestTools.h"
#include "TStyle.h"
 
using namespace std;
using namespace TMath;
 
ClassImp(TRestDetectorHitsEvent);
 
TRestDetectorHitsEvent::TRestDetectorHitsEvent() {
    fHits = new TRestHits();
 
    fPad = nullptr;
 
    fXYHitGraph = nullptr;
    fXZHitGraph = nullptr;
    fYZHitGraph = nullptr;
 
    fXYHisto = nullptr;
    fXZHisto = nullptr;
    fYZHisto = nullptr;
 
    fXZHits = nullptr;
    fYZHits = nullptr;
    fXYZHits = nullptr;
 
    fXHisto = nullptr;
    fYHisto = nullptr;
    fZHisto = nullptr;
}
 
TRestDetectorHitsEvent::~TRestDetectorHitsEvent() {
    delete fHits;
    delete fXZHits;
    delete fYZHits;
    delete fXYZHits;
 
    delete fXYHitGraph;
    delete fXZHitGraph;
    delete fYZHitGraph;
 
    delete fXYHisto;
    delete fXZHisto;
    delete fYZHisto;
 
    delete fXHisto;
    delete fYHisto;
    delete fZHisto;
}
 
void TRestDetectorHitsEvent::AddHit(Double_t x, Double_t y, Double_t z, Double_t en, Double_t t,
                                    REST_HitType type) {
    fHits->AddHit({x, y, z}, en, t, type);
}
 
void TRestDetectorHitsEvent::AddHit(const TVector3& position, Double_t energy, Double_t time,
                                    REST_HitType type) {
    fHits->AddHit(position, energy, time, type);
}
 
void TRestDetectorHitsEvent::Initialize() {
    TRestEvent::Initialize();
    fHits->RemoveHits();
 
    if (fXZHits) {
        delete fXZHits;
        fXZHits = nullptr;
    }
    if (fYZHits) {
        delete fYZHits;
        fYZHits = nullptr;
    }
    if (fXYZHits) {
        delete fXYZHits;
        fXYZHits = nullptr;
    }
 
    fXZHits = new TRestHits();
    fYZHits = new TRestHits();
    fXYZHits = new TRestHits();
}
 
void TRestDetectorHitsEvent::Sort(bool(compareCondition)(const TRestHits::iterator& hit1,
                                                         const TRestHits::iterator& hit2)) {
#ifndef __APPLE__
    if (compareCondition == 0) {
        // default sort logic: z from smaller to greater
        std::sort(fHits->begin(), fHits->end(),
                  [](const TRestHits::iterator& hit1, const TRestHits::iterator& hit2) -> bool {
                      return hit1.z() < hit2.z();
                  });
    } else {
        std::sort(fHits->begin(), fHits->end(), compareCondition);
    }
#else
    std::cout << "TRestDetectorHitsEvent::Sort is not implemented on MacOs!!" << std::endl;
    std::cout << "This method implementation should be reviewed for proper operation in Mac" << std::endl;
#endif
}
 
void TRestDetectorHitsEvent::Shuffle(int NLoop) {
    Int_t nHits = fHits->GetNumberOfHits();
    if (nHits >= 2) {
        for (int n = 0; n < NLoop; n++) {
            Int_t hit1 = (Int_t)(nHits * gRandom->Uniform(0, 1));
            Int_t hit2 = (Int_t)(nHits * gRandom->Uniform(0, 1));
 
            fHits->SwapHits(hit1, hit2);
        }
    }
}
 
TRestHits* TRestDetectorHitsEvent::GetXZHits() {
    fXZHits->RemoveHits();
 
    for (unsigned int i = 0; i < this->GetNumberOfHits(); i++) {
        if (GetType(i) == XZ) {
            fXZHits->AddHit({this->GetX(i), this->GetY(i), this->GetZ(i)}, this->GetEnergy(i),
                            this->GetTime(i), XZ);
        }
    }
 
    return fXZHits;
}
 
TRestHits* TRestDetectorHitsEvent::GetYZHits() {
    fYZHits->RemoveHits();
 
    for (unsigned int i = 0; i < this->GetNumberOfHits(); i++) {
        if (GetType(i) == YZ) {
            fYZHits->AddHit({this->GetX(i), this->GetY(i), this->GetZ(i)}, this->GetEnergy(i),
                            this->GetTime(i), YZ);
        }
    }
 
    return fYZHits;
}
 
TRestHits* TRestDetectorHitsEvent::GetXYZHits() {
    fXYZHits->RemoveHits();
 
    for (unsigned int i = 0; i < this->GetNumberOfHits(); i++) {
        if (GetType(i) == XYZ) {
            {
                fXYZHits->AddHit({this->GetX(i), this->GetY(i), this->GetZ(i)}, this->GetEnergy(i),
                                 this->GetTime(i), XYZ);
            }
        }
    }
 
    return fXYZHits;
}
 
Bool_t TRestDetectorHitsEvent::anyHitInsideCylinder(TVector3 x0, TVector3 x1, Double_t radius) {
    if (fHits->GetNumberOfHitsInsideCylinder(x0, x1, radius) > 0) {
        return true;
    }
 
    return false;
}
 
Bool_t TRestDetectorHitsEvent::allHitsInsideCylinder(TVector3 x0, TVector3 x1, Double_t radius) {
    if ((size_t)fHits->GetNumberOfHitsInsideCylinder(x0, x1, radius) == GetNumberOfHits()) return true;
 
    return false;
}
 
Double_t TRestDetectorHitsEvent::GetEnergyInCylinder(TVector3 x0, TVector3 x1, Double_t radius) {
    return fHits->GetEnergyInCylinder(x0, x1, radius);
}
 
Int_t TRestDetectorHitsEvent::GetNumberOfHitsInsideCylinder(TVector3 x0, TVector3 x1, Double_t radius) {
    return fHits->GetNumberOfHitsInsideCylinder(x0, x1, radius);
}
 
TVector3 TRestDetectorHitsEvent::GetMeanPositionInCylinder(TVector3 x0, TVector3 x1, Double_t radius) {
    return fHits->GetMeanPositionInCylinder(x0, x1, radius);
}
 
Bool_t TRestDetectorHitsEvent::anyHitInsidePrism(TVector3 x0, TVector3 x1, Double_t sizeX, Double_t sizeY,
                                                 Double_t theta) {
    if (fHits->GetNumberOfHitsInsidePrism(x0, x1, sizeX, sizeY, theta) > 0) return true;
 
    return false;
}
 
Bool_t TRestDetectorHitsEvent::allHitsInsidePrism(TVector3 x0, TVector3 x1, Double_t sizeX, Double_t sizeY,
                                                  Double_t theta) {
    if ((size_t)fHits->GetNumberOfHitsInsidePrism(x0, x1, sizeX, sizeY, theta) == GetNumberOfHits())
        return true;
 
    return false;
}
 
Double_t TRestDetectorHitsEvent::GetEnergyInPrism(TVector3 x0, TVector3 x1, Double_t sizeX, Double_t sizeY,
                                                  Double_t theta) {
    return fHits->GetEnergyInPrism(x0, x1, sizeX, sizeY, theta);
}
 
Int_t TRestDetectorHitsEvent::GetNumberOfHitsInsidePrism(TVector3 x0, TVector3 x1, Double_t sizeX,
                                                         Double_t sizeY, Double_t theta) {
    return fHits->GetNumberOfHitsInsidePrism(x0, x1, sizeX, sizeY, theta);
}
 
TVector3 TRestDetectorHitsEvent::GetMeanPositionInPrism(TVector3 x0, TVector3 x1, Double_t sizeX,
                                                        Double_t sizeY, Double_t theta) {
    return fHits->GetMeanPositionInPrism(x0, x1, sizeX, sizeY, theta);
}
 
Double_t TRestDetectorHitsEvent::GetClosestHitInsideDistanceToCylinderWall(TVector3 x0, TVector3 x1,
                                                                           Double_t radius) {
    Double_t rad2 = radius * radius;
    Double_t hitDistance = rad2;
 
    Double_t d2, l;
 
    TVector3 axis = x1 - x0;
    Double_t cylLength = axis.Mag();
 
    Int_t nhits = 0;
    for (unsigned int n = 0; n < GetNumberOfHits(); n++) {
        if (fHits->isHitNInsideCylinder(n, x0, x1, radius)) {
            l = axis.Dot(this->GetPosition(n) - x0) / cylLength;
 
            d2 = rad2 - (this->GetPosition(n) - x0).Mag2() + l * l;
 
            if (d2 < hitDistance) hitDistance = d2;
 
            nhits++;
        }
    }
 
    if (nhits == 0) return -1;
 
    return TMath::Sqrt(hitDistance);
}
 
Double_t TRestDetectorHitsEvent::GetClosestHitInsideDistanceToCylinderTop(TVector3 x0, TVector3 x1,
                                                                          Double_t radius) {
    TVector3 axis = x1 - x0;
    Double_t cylLength = axis.Mag();
 
    Double_t hitDistance = cylLength;
    Double_t d = cylLength;
 
    Int_t nhits = 0;
    for (unsigned int n = 0; n < GetNumberOfHits(); n++) {
        if (fHits->isHitNInsideCylinder(n, x0, x1, radius)) {
            d = cylLength - axis.Dot(this->GetPosition(n) - x0) / cylLength;
 
            if (d < hitDistance) hitDistance = d;
 
            nhits++;
        }
    }
 
    if (nhits == 0) return -1;
 
    return hitDistance;
}
 
Double_t TRestDetectorHitsEvent::GetClosestHitInsideDistanceToCylinderBottom(TVector3 x0, TVector3 x1,
                                                                             Double_t radius) {
    TVector3 axis = x1 - x0;
    Double_t cylLength = axis.Mag();
 
    Double_t hitDistance = cylLength;
    Double_t d = cylLength;
 
    Int_t nhits = 0;
    for (unsigned int n = 0; n < GetNumberOfHits(); n++) {
        if (fHits->isHitNInsideCylinder(n, x0, x1, radius)) {
            d = axis.Dot(this->GetPosition(n) - x0) / cylLength;
 
            if (d < hitDistance) hitDistance = d;
 
            nhits++;
        }
    }
 
    if (nhits == 0) return -1;
 
    return hitDistance;
}
 
Double_t TRestDetectorHitsEvent::GetClosestHitInsideDistanceToPrismWall(TVector3 x0, TVector3 x1,
                                                                        Double_t sizeX, Double_t sizeY,
                                                                        Double_t theta) {
    Double_t dX = sizeX / 2;
    Double_t dY = sizeY / 2;
 
    Double_t hitDistance = max(dX, dY);
 
    Double_t d;
    Int_t nhits = 0;
    for (unsigned int n = 0; n < GetNumberOfHits(); n++) {
        if (fHits->isHitNInsidePrism(n, x0, x1, sizeX, sizeY, theta)) {
            dX = sizeX / 2 - TMath::Abs((this->GetPosition(n) - x0).X());
            dY = sizeY / 2 - TMath::Abs((this->GetPosition(n) - x0).Y());
 
            d = min(dX, dY);
 
            if (d < hitDistance) hitDistance = d;
 
            nhits++;
        }
    }
 
    if (nhits == 0) return -1;
 
    return hitDistance;
}
 
Double_t TRestDetectorHitsEvent::GetClosestHitInsideDistanceToPrismTop(TVector3 x0, TVector3 x1,
                                                                       Double_t sizeX, Double_t sizeY,
                                                                       Double_t theta) {
    TVector3 axis = x1 - x0;
    Double_t prismLength = axis.Mag();
 
    Double_t hitDistance = prismLength;
 
    Double_t d;
    Int_t nhits = 0;
    for (unsigned int n = 0; n < GetNumberOfHits(); n++) {
        if (fHits->isHitNInsidePrism(n, x0, x1, sizeX, sizeY, theta)) {
            d = prismLength - axis.Dot(this->GetPosition(n) - x0) / prismLength;
 
            if (d < hitDistance) hitDistance = d;
 
            nhits++;
        }
    }
 
    if (nhits == 0) return -1;
 
    return hitDistance;
}
 
Double_t TRestDetectorHitsEvent::GetClosestHitInsideDistanceToPrismBottom(const TVector3& x0,
                                                                          const TVector3& x1, Double_t sizeX,
                                                                          Double_t sizeY, Double_t theta) {
    TVector3 axis = x1 - x0;
    Double_t prismLength = axis.Mag();
 
    Double_t hitDistance = prismLength;
 
    Double_t d;
    Int_t nhits = 0;
    for (unsigned int n = 0; n < GetNumberOfHits(); n++) {
        if (fHits->isHitNInsidePrism(n, x0, x1, sizeX, sizeY, theta)) {
            d = axis.Dot(this->GetPosition(n) - x0) / prismLength;
 
            if (d < hitDistance) hitDistance = d;
 
            nhits++;
        }
    }
 
    if (nhits == 0) return -1;
 
    return hitDistance;
}
 
TPad* TRestDetectorHitsEvent::DrawEvent(const TString& option) {
    vector<TString> optList = Vector_cast<string, TString>(TRestTools::GetOptions((string)option));
 
    for (unsigned int n = 0; n < optList.size(); n++) {
        if (optList[n] == "print") this->PrintEvent();
    }
 
    optList.erase(std::remove(optList.begin(), optList.end(), "print"), optList.end());
 
    //  which means that it should be extracted from the hit array
    if (optList.size() == 0) optList.push_back("hist(Cont1,col)");
 
    if (fPad != nullptr) {
        delete fPad;
        fPad = nullptr;
    }
 
    fPad = new TPad(this->GetName(), " ", 0, 0, 1, 1);
    fPad->Divide(3, 2 * optList.size());
    fPad->Draw();
 
    // Give each sub-pad enough margin so the axis titles (in particular the bottom-row
    // x-axis titles and the wider 2-digit y-axis labels) are not clipped at the edges.
    for (int i = 1; i <= 6 * (int)optList.size(); i++) {
        fPad->cd(i);
        gPad->SetBottomMargin(0.13);
        gPad->SetLeftMargin(0.15);
    }
 
    Int_t column = 0;
    for (unsigned int n = 0; n < optList.size(); n++) {
        string optionStr = (string)optList[n];
        TString drawEventOption = optList[n];
 
        // Generating drawOption argument
        size_t startPos = optionStr.find("(");
        if (startPos != string::npos) drawEventOption = optList[n](0, startPos);
 
        // Generating histogram option argument
        string histOption = "";
        size_t endPos = optionStr.find(")");
        if (endPos != string::npos) {
            TString histOptionTmp = optList[n](startPos + 1, endPos - startPos - 1);
 
            histOption = (string)histOptionTmp;
            size_t pos = 0;
            while ((pos = histOption.find(",", pos)) != string::npos) {
                histOption.replace(pos, 1, ":");
                pos = pos + 1;
            }
        }
 
        // Generating histogram pitch argument
        string pitchOption = "";
 
        startPos = optionStr.find("[");
        endPos = optionStr.find("]");
        Double_t pitch = 0;
        if (endPos != string::npos) {
            TString pitchOption = optList[n](startPos + 1, endPos - startPos - 1);
            pitch = stod((string)pitchOption);
        }
 
        if (drawEventOption == "graph") this->DrawGraphs(column);
 
        if (drawEventOption == "hist") this->DrawHistograms(column, histOption, pitch);
    }
 
    return fPad;
}
 
void TRestDetectorHitsEvent::DrawGraphs(Int_t& column) {
    if (fXYHitGraph != nullptr) {
        delete fXYHitGraph;
        fXYHitGraph = nullptr;
    }
    if (fXZHitGraph != nullptr) {
        delete fXZHitGraph;
        fXZHitGraph = nullptr;
    }
    if (fYZHitGraph != nullptr) {
        delete fYZHitGraph;
        fYZHitGraph = nullptr;
    }
 
    vector<vector<Double_t>> xz(2, vector<Double_t>(this->GetNumberOfHits()));
    vector<vector<Double_t>> yz(2, vector<Double_t>(this->GetNumberOfHits()));
    vector<vector<Double_t>> xy(2, vector<Double_t>(this->GetNumberOfHits()));
 
    /* {{{ Creating xz, yz, and xy arrays and initializing graphs */
    Int_t nXZ = 0;
    Int_t nYZ = 0;
    Int_t nXY = 0;
 
    for (unsigned int nhit = 0; nhit < this->GetNumberOfHits(); nhit++) {
        Double_t x = fHits->GetX(nhit);
        Double_t y = fHits->GetY(nhit);
        Double_t z = fHits->GetZ(nhit);
        int type = fHits->GetType(nhit);
 
        if (type % XZ == 0) {
            xz[0][nXZ] = x;
            xz[1][nXZ] = z;
            nXZ++;
        }
 
        if (type % YZ == 0) {
            yz[0][nYZ] = y;
            yz[1][nYZ] = z;
            nYZ++;
        }
 
        if (type % XY == 0) {
            xy[0][nXY] = x;
            xy[1][nXY] = y;
            nXY++;
        }
    }
 
    fXZHitGraph = new TGraph(nXZ, &xz[0][0], &xz[1][0]);
    fXZHitGraph->SetMarkerColor(kBlue);
    fXZHitGraph->SetMarkerSize(0.3);
    fXZHitGraph->SetMarkerStyle(20);
 
    fYZHitGraph = new TGraph(nYZ, &yz[0][0], &yz[1][0]);
    fYZHitGraph->SetMarkerColor(kRed);
    fYZHitGraph->SetMarkerSize(0.3);
    fYZHitGraph->SetMarkerStyle(20);
 
    fXYHitGraph = new TGraph(nXY, &xy[0][0], &xy[1][0]);
    fXYHitGraph->SetMarkerColor(kBlack);
    fXYHitGraph->SetMarkerSize(0.3);
    fXYHitGraph->SetMarkerStyle(20);
    /* }}} */
 
    char title[256];
    sprintf(title, "Event ID %d", this->GetID());
 
    if (nXZ > 0) {
        fPad->cd(1 + 3 * column);
        fXZHitGraph->SetTitle(title);
        fXZHitGraph->Draw("AP*");
 
        fXZHitGraph->GetXaxis()->SetTitle("X-axis (mm)");
        fXZHitGraph->GetYaxis()->SetTitle("Z-axis (mm)");
    }
 
    if (nYZ > 0) {
        fPad->cd(2 + 3 * column);
        fYZHitGraph->SetTitle(title);
        fYZHitGraph->Draw("AP");
 
        fYZHitGraph->GetXaxis()->SetTitle("Y-axis (mm)");
        fYZHitGraph->GetYaxis()->SetTitle("Z-axis (mm)");
    }
 
    if (nXY > 0) {
        fPad->cd(3 + 3 * column);
        fXYHitGraph->SetTitle(title);
        fXYHitGraph->Draw("AP");
 
        fXYHitGraph->GetXaxis()->SetTitle("X-axis (mm)");
        fXYHitGraph->GetYaxis()->SetTitle("Y-axis (mm)");
    }
 
    column++;
}
 
void TRestDetectorHitsEvent::DrawHistograms(Int_t& column, const TString& histOption, double pitch) {
    std::vector<double> fX, fY, fZ;
    for (unsigned int i = 0; i < GetNumberOfHits(); i++) {
        if (GetType(i) % X == 0) fX.emplace_back(GetX(i));
        if (GetType(i) % Y == 0) fY.emplace_back(GetY(i));
        if (GetType(i) % Z == 0) fZ.emplace_back(GetZ(i));
    }
 
    double maxX, minX, maxY, minY, maxZ, minZ;
    int nBinsX, nBinsY, nBinsZ;
    // GetBoundaries reads dist.front()/back(), which is undefined for an empty projection
    // (e.g. an event with no hits of a given type). Guard it so DrawEvent never crashes.
    auto getBoundaries = [](std::vector<double>& v, double& mx, double& mn, int& nb) {
        if (v.empty()) {
            mn = -1;
            mx = 1;
            nb = 1;
            return;
        }
        TRestHits::GetBoundaries(v, mx, mn, nb);
    };
    getBoundaries(fX, maxX, minX, nBinsX);
    getBoundaries(fY, maxY, minY, nBinsY);
    getBoundaries(fZ, maxZ, minZ, nBinsZ);
 
    if (pitch > 0) {
        nBinsX = std::round((maxX - minX) / pitch);
        nBinsY = std::round((maxY - minY) / pitch);
        nBinsZ = std::round((maxZ - minZ) / pitch);
    }
 
    // A histogram needs at least one bin; ensure we never request 0 bins for a degenerate
    // event (all hits sharing a coordinate, or a sub-pitch extent when a pitch is given).
    nBinsX = std::max(1, nBinsX);
    nBinsY = std::max(1, nBinsY);
    nBinsZ = std::max(1, nBinsZ);
 
    delete fXYHisto;
    fXYHisto = nullptr;
    delete fXZHisto;
    fXZHisto = nullptr;
    delete fYZHisto;
    fYZHisto = nullptr;
 
    delete fXHisto;
    fXHisto = nullptr;
    delete fYHisto;
    fYHisto = nullptr;
    delete fZHisto;
    fZHisto = nullptr;
 
    fXYHisto = new TH2F("XY", "", nBinsX, minX, maxX, nBinsY, minY, maxY);
    fXYHisto->SetDirectory(nullptr);
    fXZHisto = new TH2F("XZ", "", nBinsX, minX, maxX, nBinsZ, minZ, maxZ);
    fXZHisto->SetDirectory(nullptr);
    fYZHisto = new TH2F("YZ", "", nBinsY, minY, maxY, nBinsZ, minZ, maxZ);
    fYZHisto->SetDirectory(nullptr);
 
    fXHisto = new TH1F("X", "", nBinsX, minX, maxX);
    fXHisto->SetDirectory(nullptr);
    fYHisto = new TH1F("Y", "", nBinsY, minY, maxY);
    fYHisto->SetDirectory(nullptr);
    fZHisto = new TH1F("Z", "", nBinsZ, minZ, maxZ);
    fZHisto->SetDirectory(nullptr);
 
    fXYHisto->SetStats(false);
    fXZHisto->SetStats(false);
    fYZHisto->SetStats(false);
 
    fXHisto->SetStats(false);
    fYHisto->SetStats(false);
    fZHisto->SetStats(false);
 
    Int_t nYZ = 0, nXZ = 0, nXY = 0;
    Int_t nX = 0, nY = 0, nZ = 0;
 
    for (unsigned int nhit = 0; nhit < this->GetNumberOfHits(); nhit++) {
        Double_t x = fHits->GetX(nhit);
        Double_t y = fHits->GetY(nhit);
        Double_t z = fHits->GetZ(nhit);
        Double_t en = fHits->GetEnergy(nhit);
        int type = fHits->GetType(nhit);
 
        if (type % XZ == 0) {
            fXZHisto->Fill(x, z, en);
            nXZ++;
        }
 
        if (type % YZ == 0) {
            fYZHisto->Fill(y, z, en);
            nYZ++;
        }
 
        if (type % XY == 0) {
            fXYHisto->Fill(x, y, en);
            nXY++;
        }
 
        if (type % X == 0) {
            fXHisto->Fill(x);
            nX++;
        }
 
        if (type % Y == 0) {
            fYHisto->Fill(y);
            nY++;
        }
 
        if (type % Z == 0) {
            fZHisto->Fill(z);
            nZ++;
        }
    }
 
    TStyle style;
    style.SetPalette(1);
 
    if (nXZ > 0) {
        fPad->cd(1 + 3 * column);
        fXZHisto->Draw(histOption);
        fXZHisto->GetXaxis()->SetTitle("X-axis (mm)");
        fXZHisto->GetYaxis()->SetTitle("Z-axis (mm)");
        fXZHisto->GetYaxis()->SetTitleSize(1.4 * fXZHisto->GetYaxis()->GetTitleSize());
        fXZHisto->GetXaxis()->SetTitleSize(1.4 * fXZHisto->GetXaxis()->GetTitleSize());
        fXZHisto->GetYaxis()->SetLabelSize(1.25 * fXZHisto->GetYaxis()->GetLabelSize());
        fXZHisto->GetXaxis()->SetLabelSize(1.25 * fXZHisto->GetXaxis()->GetLabelSize());
        fXZHisto->GetYaxis()->SetTitleOffset(1);
    }
 
    if (nYZ > 0) {
        fPad->cd(2 + 3 * column);
        fYZHisto->Draw(histOption);
        fYZHisto->GetXaxis()->SetTitle("Y-axis (mm)");
        fYZHisto->GetYaxis()->SetTitle("Z-axis (mm)");
        fYZHisto->GetYaxis()->SetTitleSize(1.4 * fYZHisto->GetYaxis()->GetTitleSize());
        fYZHisto->GetXaxis()->SetTitleSize(1.4 * fYZHisto->GetXaxis()->GetTitleSize());
        fYZHisto->GetYaxis()->SetLabelSize(1.25 * fYZHisto->GetYaxis()->GetLabelSize());
        fYZHisto->GetXaxis()->SetLabelSize(1.25 * fYZHisto->GetXaxis()->GetLabelSize());
        fYZHisto->GetYaxis()->SetTitleOffset(1);
    }
 
    if (nXY > 0) {
        fPad->cd(3 + 3 * column);
        fXYHisto->Draw(histOption);
        fXYHisto->GetXaxis()->SetTitle("X-axis (mm)");
        fXYHisto->GetYaxis()->SetTitle("Y-axis (mm)");
        fXYHisto->GetYaxis()->SetTitleSize(1.4 * fXYHisto->GetYaxis()->GetTitleSize());
        fXYHisto->GetXaxis()->SetTitleSize(1.4 * fXYHisto->GetXaxis()->GetTitleSize());
        fXYHisto->GetYaxis()->SetLabelSize(1.25 * fXYHisto->GetYaxis()->GetLabelSize());
        fXYHisto->GetXaxis()->SetLabelSize(1.25 * fXYHisto->GetXaxis()->GetLabelSize());
        fXYHisto->GetYaxis()->SetTitleOffset(1);
    } else {
        // When there is no XY projection (a strip readout, where hits are XZ/YZ) the pad is
        // otherwise empty, so we use it for a 3D view of the event. A strip hit has one
        // undetermined coordinate, so it is drawn as a line spanning that axis (XZ -> line
        // along Y, YZ -> along X, XY -> along Z), colored by energy. An XZ and a YZ line
        // cross in 3D where they share the same Z; that crossing (x_from_XZ, y_from_YZ, z)
        // is the reconstructed hit position, so we highlight it with a marker, which sketches
        // the event in 3D. (With several strips at the same Z these crossings also include
        // ambiguous "ghost" combinations, inherent to strip readouts.) All objects are owned
        // by the pad (kCanDelete) so they are freed when fPad is recreated on next DrawEvent.
        fPad->cd(3 + 3 * column);
 
        TH3F* frame = new TH3F("XYZframe", "3D hits;X-axis (mm);Y-axis (mm);Z-axis (mm)", 1, minX, maxX, 1,
                               minY, maxY, 1, minZ, maxZ);
        frame->SetDirectory(nullptr);
        frame->SetStats(false);
        frame->SetBit(kCanDelete);
        frame->Draw();
 
        Double_t eMin = 1e9, eMax = -1e9;
        for (unsigned int nhit = 0; nhit < GetNumberOfHits(); nhit++) {
            Double_t e = fHits->GetEnergy(nhit);
            if (e < eMin) eMin = e;
            if (e > eMax) eMax = e;
        }
        const Int_t nColors = TColor::GetNumberOfColors();
        auto energyColor = [&](Double_t e) {
            Int_t ci = (eMax > eMin) ? (Int_t)((e - eMin) / (eMax - eMin) * (nColors - 1)) : 0;
            return TColor::GetColorPalette(ci);
        };
 
        std::vector<std::array<Double_t, 3>> xzHits;  // {x, z, energy}
        std::vector<std::array<Double_t, 3>> yzHits;  // {y, z, energy}
        for (unsigned int nhit = 0; nhit < GetNumberOfHits(); nhit++) {
            Double_t x = fHits->GetX(nhit);
            Double_t y = fHits->GetY(nhit);
            Double_t z = fHits->GetZ(nhit);
            Double_t en = fHits->GetEnergy(nhit);
            int type = fHits->GetType(nhit);
            bool hasX = (type % X == 0);
            bool hasY = (type % Y == 0);
 
            TPolyLine3D* line = new TPolyLine3D(2);
            if (!hasX) {  // undetermined X (e.g. YZ hit): line along X
                line->SetPoint(0, minX, y, z);
                line->SetPoint(1, maxX, y, z);
                yzHits.push_back({y, z, en});
            } else if (!hasY) {  // undetermined Y (e.g. XZ hit): line along Y
                line->SetPoint(0, x, minY, z);
                line->SetPoint(1, x, maxY, z);
                xzHits.push_back({x, z, en});
            } else {  // undetermined Z (e.g. XY hit): line along Z
                line->SetPoint(0, x, y, minZ);
                line->SetPoint(1, x, y, maxZ);
            }
            line->SetLineColor(energyColor(en));
            line->SetBit(kCanDelete);
            line->Draw();
        }
 
        // Highlight the XZ/YZ crossings (same Z) -> reconstructed 3D positions of the event.
        const Double_t zTol = 0.05 * (maxZ - minZ);
        for (const auto& xz : xzHits) {
            for (const auto& yz : yzHits) {
                if (std::abs(xz[1] - yz[1]) > zTol) continue;
                TPolyMarker3D* marker = new TPolyMarker3D(1);
                marker->SetPoint(0, xz[0], yz[0], 0.5 * (xz[1] + yz[1]));
                marker->SetMarkerStyle(20);
                marker->SetMarkerSize(1.2);
                marker->SetMarkerColor(energyColor(0.5 * (xz[2] + yz[2])));
                marker->SetBit(kCanDelete);
                marker->Draw();
            }
        }
    }
 
    column++;
 
    if (nX > 0) {
        fPad->cd(1 + 3 * column);
        fXHisto->Draw(histOption);
        fXHisto->GetXaxis()->SetTitle("X-axis (mm)");
        fXHisto->GetYaxis()->SetTitle("Number of events");
        fXHisto->GetYaxis()->SetTitleSize(1.4 * fXHisto->GetYaxis()->GetTitleSize());
        fXHisto->GetXaxis()->SetTitleSize(1.4 * fXHisto->GetXaxis()->GetTitleSize());
        fXHisto->GetYaxis()->SetLabelSize(1.25 * fXHisto->GetYaxis()->GetLabelSize());
        fXHisto->GetXaxis()->SetLabelSize(1.25 * fXHisto->GetXaxis()->GetLabelSize());
        fXHisto->GetYaxis()->SetTitleOffset(1);
    }
 
    if (nY > 0) {
        fPad->cd(2 + 3 * column);
        fYHisto->Draw(histOption);
        fYHisto->GetXaxis()->SetTitle("Y-axis (mm)");
        fYHisto->GetYaxis()->SetTitle("Number of events");
        fYHisto->GetYaxis()->SetTitleSize(1.4 * fYHisto->GetYaxis()->GetTitleSize());
        fYHisto->GetXaxis()->SetTitleSize(1.4 * fYHisto->GetXaxis()->GetTitleSize());
        fYHisto->GetYaxis()->SetLabelSize(1.25 * fYHisto->GetYaxis()->GetLabelSize());
        fYHisto->GetXaxis()->SetLabelSize(1.25 * fYHisto->GetXaxis()->GetLabelSize());
        fYHisto->GetYaxis()->SetTitleOffset(1);
    }
 
    if (nZ > 0) {
        fPad->cd(3 + 3 * column);
        fZHisto->Draw(histOption);
        fZHisto->GetXaxis()->SetTitle("Z-axis (mm)");
        fZHisto->GetYaxis()->SetTitle("Number of events");
        fZHisto->GetYaxis()->SetTitleSize(1.4 * fYHisto->GetYaxis()->GetTitleSize());
        fZHisto->GetXaxis()->SetTitleSize(1.4 * fYHisto->GetXaxis()->GetTitleSize());
        fZHisto->GetYaxis()->SetLabelSize(1.25 * fYHisto->GetYaxis()->GetLabelSize());
        fZHisto->GetXaxis()->SetLabelSize(1.25 * fYHisto->GetXaxis()->GetLabelSize());
        fZHisto->GetYaxis()->SetTitleOffset(1);
    }
 
    column++;
}
 
void TRestDetectorHitsEvent::PrintEvent(Int_t nHits) const {
    TRestEvent::PrintEvent();
 
    cout << "Total energy : " << GetTotalEnergy() << endl;
    cout << "Mean position : ( " << GetMeanPositionX() << " , " << GetMeanPositionY() << " , "
         << GetMeanPositionZ() << " ) " << endl;
    cout << "Number of hits : " << fHits->GetNumberOfHits() << endl;
    if (nHits != -1) {
        cout << "+++++++++++++++++++++++" << endl;
        cout << "Printing only the first " << nHits << " hits" << endl;
    }
 
    fHits->PrintHits(nHits);
}
 
TH2F* TRestDetectorHitsEvent::GetXYHistogram(std::vector<float> ranges, Double_t pitch, Double_t border) {
    double maxX, minX, maxY, minY;
    int nBinsX, nBinsY;
 
    if (ranges.size() == 6) {
        nBinsX = ranges[0];
        minX = ranges[1];
        maxX = ranges[2];
 
        nBinsY = ranges[3];
        minY = ranges[4];
        maxY = ranges[5];
    } else {
        std::vector<double> x, y, en;
 
        for (unsigned int i = 0; i < GetNumberOfHits(); i++) {
            if (GetType(i) % X == 0) x.emplace_back(GetX(i));
            if (GetType(i) % Y == 0) y.emplace_back(GetY(i));
        }
 
        if (x.empty() && y.empty()) {
            std::cout << "TRestDetectorHitsEvent::GetXYHistogram. Empty histogram!" << std::endl;
            return nullptr;
        }
 
        TRestHits::GetBoundaries(x, maxX, minX, nBinsX);
        TRestHits::GetBoundaries(y, maxY, minY, nBinsY);
 
        maxX += border;
        minX -= border;
 
        maxY += border;
        minY -= border;
 
        if (pitch > 0) {
            nBinsX = std::round((maxX - minX) / pitch);
            nBinsY = std::round((maxY - minY) / pitch);
        }
    }
 
    delete fXYHisto;
    fXYHisto = nullptr;
    fXYHisto = new TH2F("XY", "", nBinsX, minX, maxX, nBinsY, minY, maxY);
    fXYHisto->SetDirectory(nullptr);
    fXYHisto->SetStats(false);
 
    Int_t nXY = 0;
 
    for (unsigned int nhit = 0; nhit < this->GetNumberOfHits(); nhit++) {
        Double_t x = fHits->GetX(nhit);
        Double_t y = fHits->GetY(nhit);
        Double_t en = fHits->GetEnergy(nhit);
        int type = fHits->GetType(nhit);
 
        if (type % XY == 0 || type % XYZ == 0) {
            fXYHisto->Fill(x, y, en);
            nXY++;
        }
    }
 
    TStyle style;
    style.SetPalette(1);
 
    if (nXY > 0) {
        fXYHisto->Draw("col");
        fXYHisto->Draw("cont3 same");
        fXYHisto->GetXaxis()->SetTitle("X-axis (mm)");
        fXYHisto->GetYaxis()->SetTitle("Y-axis (mm)");
        fXYHisto->GetYaxis()->SetTitleSize(1.4 * fXYHisto->GetYaxis()->GetTitleSize());
        fXYHisto->GetXaxis()->SetTitleSize(1.4 * fXYHisto->GetXaxis()->GetTitleSize());
        fXYHisto->GetYaxis()->SetLabelSize(1.25 * fXYHisto->GetYaxis()->GetLabelSize());
        fXYHisto->GetXaxis()->SetLabelSize(1.25 * fXYHisto->GetXaxis()->GetLabelSize());
        fXYHisto->GetYaxis()->SetTitleOffset(1);
    }
 
    return fXYHisto;
}
 
TH2F* TRestDetectorHitsEvent::GetXZHistogram(std::vector<float> ranges, Double_t pitch, Double_t border) {
    double maxX, minX, maxZ, minZ;
    int nBinsX, nBinsZ;
 
    if (ranges.size() == 6) {
        nBinsX = ranges[0];
        minX = ranges[1];
        maxX = ranges[2];
 
        nBinsZ = ranges[3];
        minZ = ranges[4];
        maxZ = ranges[5];
    } else {
        std::vector<double> x, z;
 
        for (unsigned int i = 0; i < GetNumberOfHits(); i++) {
            if (GetType(i) % X == 0) x.emplace_back(GetX(i));
            if (GetType(i) % Z == 0) z.emplace_back(GetZ(i));
        }
 
        if (x.empty() || z.empty()) {
            std::cout << "TRestDetectorHitsEvent::GetXZHistogram. Empty histogram!" << std::endl;
            return nullptr;
        }
 
        TRestHits::GetBoundaries(x, maxX, minX, nBinsX);
        TRestHits::GetBoundaries(z, maxZ, minZ, nBinsZ);
 
        maxX += border;
        minX -= border;
 
        maxZ += border;
        minZ -= border;
 
        if (pitch > 0) {
            nBinsX = std::round((maxX - minX) / pitch);
            nBinsZ = std::round((maxZ - minZ) / pitch);
        }
    }
 
    delete fXZHisto;
    fXZHisto = nullptr;
    fXZHisto = new TH2F("XZ", "", nBinsX, minX, maxX, nBinsZ, minZ, maxZ);
    fXZHisto->SetDirectory(nullptr);
    fXZHisto->SetStats(false);
 
    Int_t nXZ = 0;
 
    for (unsigned int nhit = 0; nhit < this->GetNumberOfHits(); nhit++) {
        Double_t x = fHits->GetX(nhit);
        Double_t z = fHits->GetZ(nhit);
        Double_t en = fHits->GetEnergy(nhit);
        int type = fHits->GetType(nhit);
 
        if (type % XZ == 0 || type % XYZ == 0) {
            fXZHisto->Fill(x, z, en);
            nXZ++;
        }
    }
 
    TStyle style;
    style.SetPalette(1);
 
    if (nXZ > 0) {
        fXZHisto->Draw("col");
        fXZHisto->Draw("cont3 same");
        fXZHisto->GetXaxis()->SetTitle("X-axis (mm)");
        fXZHisto->GetYaxis()->SetTitle("Z-axis (mm)");
        fXZHisto->GetYaxis()->SetTitleSize(1.4 * fXZHisto->GetYaxis()->GetTitleSize());
        fXZHisto->GetXaxis()->SetTitleSize(1.4 * fXYHisto->GetXaxis()->GetTitleSize());
        fXZHisto->GetYaxis()->SetLabelSize(1.25 * fXYHisto->GetYaxis()->GetLabelSize());
        fXZHisto->GetXaxis()->SetLabelSize(1.25 * fXYHisto->GetXaxis()->GetLabelSize());
        fXZHisto->GetYaxis()->SetTitleOffset(1);
    }
 
    return fXZHisto;
}
 
TH2F* TRestDetectorHitsEvent::GetYZHistogram(std::vector<float> ranges, Double_t pitch, Double_t border) {
    double maxY, minY, maxZ, minZ;
    int nBinsY, nBinsZ;
 
    if (ranges.size() == 6) {
        nBinsY = ranges[0];
        minY = ranges[1];
        maxY = ranges[2];
 
        nBinsZ = ranges[3];
        minZ = ranges[4];
        maxZ = ranges[5];
    } else {
        std::vector<double> y, z, en;
 
        for (unsigned int i = 0; i < GetNumberOfHits(); i++) {
            if (GetType(i) % Y == 0) y.emplace_back(GetY(i));
            if (GetType(i) % Z == 0) z.emplace_back(GetZ(i));
        }
 
        if (y.empty() || z.empty()) {
            std::cout << "TRestDetectorHitsEvent::GetYZHistogram. Empty histogram!" << std::endl;
            return nullptr;
        }
 
        TRestHits::GetBoundaries(y, maxY, minY, nBinsY);
        TRestHits::GetBoundaries(z, maxZ, minZ, nBinsZ);
 
        maxY += border;
        minY -= border;
 
        maxZ += border;
        minZ -= border;
 
        if (pitch > 0) {
            nBinsY = std::round((maxY - minZ) / pitch);
            nBinsZ = std::round((maxZ - minZ) / pitch);
        }
    }
 
    delete fYZHisto;
    fYZHisto = nullptr;
    fYZHisto = new TH2F("YZ", "", nBinsY, minY, maxY, nBinsZ, minZ, maxZ);
    fYZHisto->SetDirectory(nullptr);
    fYZHisto->SetStats(false);
 
    Int_t nYZ = 0;
 
    for (unsigned int nhit = 0; nhit < this->GetNumberOfHits(); nhit++) {
        Double_t y = fHits->GetY(nhit);
        Double_t z = fHits->GetZ(nhit);
        Double_t en = fHits->GetEnergy(nhit);
        int type = fHits->GetType(nhit);
 
        if (type % YZ == 0 || type % XYZ == 0) {
            fYZHisto->Fill(y, z, en);
            nYZ++;
        }
    }
 
    TStyle style;
    style.SetPalette(1);
 
    if (nYZ > 0) {
        fYZHisto->Draw("col");
        fYZHisto->Draw("cont3 same");
        fYZHisto->GetXaxis()->SetTitle("Y-axis (mm)");
        fYZHisto->GetYaxis()->SetTitle("Z-axis (mm)");
        fYZHisto->GetYaxis()->SetTitleSize(1.4 * fYZHisto->GetYaxis()->GetTitleSize());
        fYZHisto->GetXaxis()->SetTitleSize(1.4 * fYZHisto->GetXaxis()->GetTitleSize());
        fYZHisto->GetYaxis()->SetLabelSize(1.25 * fYZHisto->GetYaxis()->GetLabelSize());
        fYZHisto->GetXaxis()->SetLabelSize(1.25 * fYZHisto->GetXaxis()->GetLabelSize());
        fYZHisto->GetYaxis()->SetTitleOffset(1);
    }
 
    return fYZHisto;
}