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 "TCanvas.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; }
 
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)[3]");
 
    if (fPad != nullptr) {
        delete fPad;
        fPad = nullptr;
    }
 
    fPad = new TPad(this->GetName(), " ", 0, 0, 1, 1);
    fPad->Divide(3, 2 * optList.size());
    fPad->Draw();
 
    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;
    TRestHits::GetBoundaries(fX, maxX, minX, nBinsX);
    TRestHits::GetBoundaries(fY, maxY, minY, nBinsY);
    TRestHits::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);
    }
 
    delete fXYHisto;
    delete fXZHisto;
    delete fYZHisto;
 
    delete fXHisto;
    delete fYHisto;
    delete fZHisto;
 
    fXYHisto = new TH2F("XY", "", nBinsX, minX, maxX, nBinsY, minY, maxY);
    fXZHisto = new TH2F("XZ", "", nBinsX, minX, maxX, nBinsZ, minZ, maxZ);
    fYZHisto = new TH2F("YZ", "", nBinsY, minY, maxY, nBinsZ, minZ, maxZ);
 
    fXHisto = new TH1F("X", "", nBinsX, minX, maxX);
    fYHisto = new TH1F("Y", "", nBinsY, minY, maxY);
    fZHisto = new TH1F("Z", "", nBinsZ, minZ, maxZ);
 
    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);
    }
 
    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 = new TH2F("XY", "", nBinsX, minX, maxX, nBinsY, minY, maxY);
    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 = new TH2F("XZ", "", nBinsX, minX, maxX, nBinsZ, minZ, maxZ);
    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 = new TH2F("YZ", "", nBinsY, minY, maxY, nBinsZ, minZ, maxZ);
    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;
}