diff --git a/images/rotation.png b/images/rotation.png new file mode 100644 index 00000000..af69855f Binary files /dev/null and b/images/rotation.png differ diff --git a/images/specular.png b/images/specular.png new file mode 100644 index 00000000..cf0b6b4d Binary files /dev/null and b/images/specular.png differ diff --git a/images/transformedHitmap.png b/images/transformedHitmap.png new file mode 100644 index 00000000..1f18695b Binary files /dev/null and b/images/transformedHitmap.png differ diff --git a/images/translation.png b/images/translation.png new file mode 100644 index 00000000..40d77130 Binary files /dev/null and b/images/translation.png differ diff --git a/inc/TRestDetectorHitmapAnalysisProcess.h b/inc/TRestDetectorHitmapAnalysisProcess.h new file mode 100644 index 00000000..8912088c --- /dev/null +++ b/inc/TRestDetectorHitmapAnalysisProcess.h @@ -0,0 +1,88 @@ +/************************************************************************* + * 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. * + *************************************************************************/ + +#ifndef RestCore_TRestDetectorHitmapAnalysisProcess +#define RestCore_TRestDetectorHitmapAnalysisProcess + +#include + +#include "TRestEventProcess.h" + +/// A structure to define a given hit transformation +struct HitTransformation { + /// The center of rotation or a point contained in the specular plane + TVector3 position = {0, 0, 0}; + + /// The translation vector, the axis of rotation or the specular plane normal + TVector3 vector = {0, 0, 1}; + + /// The angle of rotation + Double_t angle = 0; + + /// The type of transformation (specular/rotation/translation) + std::string type = ""; + + /// A given name that allows to place an order to the transformations + std::string name = ""; +}; + +//! An analysis process to apply rotation/translation/specular to mean hit positions +class TRestDetectorHitmapAnalysisProcess : public TRestEventProcess { + private: + /// A pointer to the specific TRestDetectorHitsEvent input + TRestDetectorHitsEvent* fHitsEvent; //! + + void InitFromConfigFile() override; + + void Initialize() override; + + protected: + /// A list of transformations that can be applied to the mean positions + std::vector fTransDefinitions; + + /// An ordered list with the names of transformations that will be applied + std::vector fTransformations; + + /// It returns the transformation with a given name + HitTransformation GetTransformation(const std::string& name); + + TVector3 Specular(const TVector3& pos, const HitTransformation& tr); + TVector3 Rotation(const TVector3& pos, const HitTransformation& tr); + TVector3 Translation(const TVector3& pos, const HitTransformation& tr); + + public: + any GetInputEvent() const override { return fHitsEvent; } + any GetOutputEvent() const override { return fHitsEvent; } + + TRestEvent* ProcessEvent(TRestEvent* inputEvent) override; + + void PrintMetadata() override; + + const char* GetProcessName() const override { return "hitsTransformation"; } + + TRestDetectorHitmapAnalysisProcess(); + TRestDetectorHitmapAnalysisProcess(const char* configFilename); + ~TRestDetectorHitmapAnalysisProcess(); + + ClassDefOverride(TRestDetectorHitmapAnalysisProcess, 1); +}; +#endif diff --git a/inc/TRestDetectorHitsEvent.h b/inc/TRestDetectorHitsEvent.h index 6b5c5c7d..15f171d7 100644 --- a/inc/TRestDetectorHitsEvent.h +++ b/inc/TRestDetectorHitsEvent.h @@ -161,6 +161,10 @@ class TRestDetectorHitsEvent : public TRestEvent { TPad* DrawEvent(const TString& option = ""); + TH2F* GetXYHistogram(std::vector ranges, Double_t pitch = 3, Double_t border = 5); + TH2F* GetXZHistogram(std::vector ranges, Double_t pitch = 3, Double_t border = 5); + TH2F* GetYZHistogram(std::vector ranges, Double_t pitch = 3, Double_t border = 5); + void DrawHistograms(Int_t& column, const TString& histOption = "", double pitch = 0); void DrawGraphs(Int_t& column); diff --git a/inc/TRestDetectorHitsRotateAndTranslateProcess.h b/inc/TRestDetectorHitsRotateAndTranslateProcess.h deleted file mode 100644 index cc156110..00000000 --- a/inc/TRestDetectorHitsRotateAndTranslateProcess.h +++ /dev/null @@ -1,79 +0,0 @@ -///______________________________________________________________________________ -///______________________________________________________________________________ -///______________________________________________________________________________ -/// -/// -/// RESTSoft : Software for Rare Event Searches with TPCs -/// -/// TRestDetectorHitsRotateAndTranslateProcess.h -/// -/// march 2016: First concept -/// Created as part of the conceptualization of existing REST -/// software. -/// Javier G. Garza -///_______________________________________________________________________________ - -#ifndef RestCore_TRestDetectorHitsRotateAndTranslateProcess -#define RestCore_TRestDetectorHitsRotateAndTranslateProcess - -#include - -#include "TRestDetectorGas.h" -#include "TRestDetectorHitsEvent.h" - -class TRestDetectorHitsRotateAndTranslateProcess : public TRestEventProcess { - private: - TRestDetectorHitsEvent* fInputHitsEvent; //! - TRestDetectorHitsEvent* fOutputHitsEvent; //! - - TVector3 fTranslationVector = {0, 0, 0}; ///< translation vector (x,y,z) - TVector3 fRotationVector = {0, 0, 0}; ///< rotation vector around axis (x,y,z) in radians - TVector3 fRotationCenter = {0, 0, 0}; ///< rotation center (x,y,z) - - void InitFromConfigFile() override; - void Initialize() override; - void LoadDefaultConfig(); - - protected: - // add here the members of your event process - - public: - any GetInputEvent() const override { return fInputHitsEvent; } - any GetOutputEvent() const override { return fOutputHitsEvent; } - - void InitProcess() override; - TRestEvent* ProcessEvent(TRestEvent* inputEvent) override; - void EndProcess() override; - - void LoadConfig(const std::string& configFilename); - - void PrintMetadata() override { - BeginPrintProcess(); - - RESTMetadata << "Translation vector (mm): (" << fTranslationVector.X() << ", " - << fTranslationVector.Y() << ", " << fTranslationVector.Z() << ")" << RESTendl; - TVector3 rotationVectorDegrees = fRotationVector * (180. / TMath::Pi()); - RESTMetadata << "Rotation vector (degrees): (" << rotationVectorDegrees.X() << ", " - << rotationVectorDegrees.Y() << ", " << rotationVectorDegrees.Z() << ")" << RESTendl; - - EndPrintProcess(); - } - - const char* GetProcessName() const override { return "rotateAndTranslate"; } - - inline Double_t GetTranslationX() const { return fTranslationVector.X(); } - inline Double_t GetTranslationY() const { return fTranslationVector.Y(); } - inline Double_t GetTranslationZ() const { return fTranslationVector.Z(); } - - inline Double_t GetRotationX() const { return fRotationVector.X(); } - inline Double_t GetRotationY() const { return fRotationVector.Y(); } - inline Double_t GetRotationZ() const { return fRotationVector.Z(); } - - TRestDetectorHitsRotateAndTranslateProcess(); - explicit TRestDetectorHitsRotateAndTranslateProcess(const char* configFilename); - ~TRestDetectorHitsRotateAndTranslateProcess(); - - ClassDefOverride(TRestDetectorHitsRotateAndTranslateProcess, 2); -}; - -#endif diff --git a/inc/TRestDetectorHitsRotationProcess.h b/inc/TRestDetectorHitsRotationProcess.h new file mode 100644 index 00000000..efa38dc1 --- /dev/null +++ b/inc/TRestDetectorHitsRotationProcess.h @@ -0,0 +1,71 @@ +/************************************************************************* + * 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. * + *************************************************************************/ + +#ifndef RestCore_TRestDetectorHitsRotationProcess +#define RestCore_TRestDetectorHitsRotationProcess + +#include +#include + +/// A process to rotate hits from a given center and axis +class TRestDetectorHitsRotationProcess : public TRestEventProcess { + private: + /// A pointer to the process input event + TRestDetectorHitsEvent* fInputEvent; //! + + /// A pointer to the process output event + TRestDetectorHitsEvent* fOutputEvent; //! + + void InitFromConfigFile() override; + void Initialize() override; + + protected: + /// Angle of rotation respect to the given axis + Double_t fAngle = 0; //< + + /// Center of rotation + TVector3 fCenter = {0, 0, 0}; //< + + /// Axis of rotation + TVector3 fAxis = {0, 0, 1}; //< + + public: + any GetInputEvent() const override { return fInputEvent; } + any GetOutputEvent() const override { return fOutputEvent; } + + TRestEvent* ProcessEvent(TRestEvent* inputEvent) override; + + void PrintMetadata() override; + + const char* GetProcessName() const override { return "hitsRotation"; } + + inline Double_t GetAngle() const { return fAngle; } + inline TVector3 GetAxis() const { return fAxis; } + inline TVector3 GetCenter() const { return fCenter; } + + TRestDetectorHitsRotationProcess(); + TRestDetectorHitsRotationProcess(const char* configFilename); + ~TRestDetectorHitsRotationProcess(); + + ClassDefOverride(TRestDetectorHitsRotationProcess, 1); +}; +#endif diff --git a/inc/TRestDetectorHitsSpecularProcess.h b/inc/TRestDetectorHitsSpecularProcess.h new file mode 100644 index 00000000..dcbb3270 --- /dev/null +++ b/inc/TRestDetectorHitsSpecularProcess.h @@ -0,0 +1,71 @@ +/************************************************************************* + * 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. * + *************************************************************************/ + +#ifndef RestCore_TRestDetectorHitsSpecularProcess +#define RestCore_TRestDetectorHitsSpecularProcess + +#include +#include + +/// A process to create specular hit images using a plane definition +class TRestDetectorHitsSpecularProcess : public TRestEventProcess { + private: + /// A pointer to the process input event + TRestDetectorHitsEvent* fInputEvent; //! + + /// A pointer to the process output event + TRestDetectorHitsEvent* fOutputEvent; //! + + void InitFromConfigFile() override; + void Initialize() override; + + protected: + /// Defines the position of the plane. A point inside the plane. + TVector3 fPosition = {0, 0, 0}; //< + + /// Defines the normal to the plane we use to generate the specular image + TVector3 fNormal = {0, 0, 1}; //< + + public: + any GetInputEvent() const override { return fInputEvent; } + any GetOutputEvent() const override { return fOutputEvent; } + + TRestEvent* ProcessEvent(TRestEvent* inputEvent) override; + + void PrintMetadata() override; + + /// It returns the given process name + const char* GetProcessName() const override { return "hitsSpecular"; } + + /// It returns the plane normal + inline TVector3 GetNormal() const { return fNormal; } + + /// It returns the position of the point contained in the plane + inline TVector3 GetPosition() const { return fPosition; } + + TRestDetectorHitsSpecularProcess(); + TRestDetectorHitsSpecularProcess(const char* configFilename); + ~TRestDetectorHitsSpecularProcess(); + + ClassDefOverride(TRestDetectorHitsSpecularProcess, 1); +}; +#endif diff --git a/inc/TRestDetectorHitsTranslationProcess.h b/inc/TRestDetectorHitsTranslationProcess.h new file mode 100644 index 00000000..21a465f2 --- /dev/null +++ b/inc/TRestDetectorHitsTranslationProcess.h @@ -0,0 +1,63 @@ +/************************************************************************* + * 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. * + *************************************************************************/ + +#ifndef RestCore_TRestDetectorHitsTranslationProcess +#define RestCore_TRestDetectorHitsTranslationProcess + +#include +#include + +/// A process to translate hits by a given user amount +class TRestDetectorHitsTranslationProcess : public TRestEventProcess { + private: + /// A pointer to the process input event + TRestDetectorHitsEvent* fInputEvent; //! + + /// A pointer to the process output event + TRestDetectorHitsEvent* fOutputEvent; //! + + void InitFromConfigFile() override; + void Initialize() override; + + protected: + /// The amount to be translated every hit + TVector3 fTranslation = {0, 0, 0}; //< + + public: + any GetInputEvent() const override { return fInputEvent; } + any GetOutputEvent() const override { return fOutputEvent; } + + TRestEvent* ProcessEvent(TRestEvent* inputEvent) override; + + void PrintMetadata() override; + + const char* GetProcessName() const override { return "hitsTranslation"; } + + inline TVector3 GetTranslation() const { return fTranslation; } + + TRestDetectorHitsTranslationProcess(); + TRestDetectorHitsTranslationProcess(const char* configFilename); + ~TRestDetectorHitsTranslationProcess(); + + ClassDefOverride(TRestDetectorHitsTranslationProcess, 1); +}; +#endif diff --git a/pipeline/analysis/hitmap.C b/pipeline/analysis/hitmap.C new file mode 100644 index 00000000..ddb56d90 --- /dev/null +++ b/pipeline/analysis/hitmap.C @@ -0,0 +1,25 @@ + +Int_t hitmap(Bool_t draw = false) { + TRestDetectorHitsEvent* ev = new TRestDetectorHitsEvent(); + + Double_t xCenter = 5; + Double_t yCenter = 10; + + TRandom3* random = new TRandom3(137); + for (int n = 0; n < 10000; n++) { + Double_t x = random->Gaus(xCenter, 3); + Double_t y = random->Gaus(yCenter, 3); + + ev->AddHit(x, y, 0, 1); + } + + // Initializing processes + string cfgFile1 = "transform.rml"; + TRestDetectorHitmapAnalysisProcess* t = new TRestDetectorHitmapAnalysisProcess((char*)cfgFile1.c_str()); + t->PrintMetadata(); + + TRestDetectorHitsEvent* tEv = (TRestDetectorHitsEvent*)t->ProcessEvent(ev); + + cout << "[\033[92m OK \x1b[0m]" << endl; + return 0; +} diff --git a/pipeline/analysis/hitmap.rml b/pipeline/analysis/hitmap.rml new file mode 100644 index 00000000..09c84a39 --- /dev/null +++ b/pipeline/analysis/hitmap.rml @@ -0,0 +1,10 @@ + + + + + + + + + + diff --git a/pipeline/hits/rotation/rotation.C b/pipeline/hits/rotation/rotation.C new file mode 100644 index 00000000..02ec3573 --- /dev/null +++ b/pipeline/hits/rotation/rotation.C @@ -0,0 +1,53 @@ + +Int_t rotation(Bool_t draw = false) { + TRestDetectorHitsEvent* ev = new TRestDetectorHitsEvent(); + + Double_t xCenter = 5; + Double_t yCenter = 10; + + TRandom3* random = new TRandom3(137); + for (int n = 0; n < 10000; n++) { + Double_t x = random->Gaus(xCenter, 3); + Double_t y = random->Gaus(yCenter, 3); + + ev->AddHit(x, y, 0, 1); + } + + // Initializing processes + string cfgFile1 = "rotation90.rml"; + TRestDetectorHitsRotationProcess* rotation1 = + new TRestDetectorHitsRotationProcess((char*)cfgFile1.c_str()); + rotation1->PrintMetadata(); + + cfgFile1 = "rotation45off.rml"; + TRestDetectorHitsRotationProcess* rotation2 = + new TRestDetectorHitsRotationProcess((char*)cfgFile1.c_str()); + rotation2->PrintMetadata(); + + TRestDetectorHitsEvent* rotEvent1 = (TRestDetectorHitsEvent*)rotation1->ProcessEvent(ev); + TRestDetectorHitsEvent* rotEvent2 = (TRestDetectorHitsEvent*)rotation2->ProcessEvent(ev); + + /* Used to generate a combined plot */ + if (draw) { + TCanvas* c = new TCanvas("c0", "", 800, 600); + std::vector ranges{50, -25, 25, 50, -25, 25}; + TH2F* h = ev->GetXYHistogram(ranges); + h->Draw("col"); + h->Draw("cont3 same"); + c->Print("original.png"); + + h = rotEvent1->GetXYHistogram(ranges); + // pad->cd(); + h->Draw("col"); + h->Draw("cont3 same"); + c->Print("rotation1.png"); + + h = rotEvent2->GetXYHistogram(ranges); + h->Draw("col"); + h->Draw("cont3 same"); + c->Print("rotation2.png"); + } + + cout << "[\033[92m OK \x1b[0m]" << endl; + return 0; +} diff --git a/pipeline/hits/rotation/rotation45off.rml b/pipeline/hits/rotation/rotation45off.rml new file mode 100644 index 00000000..f24072ed --- /dev/null +++ b/pipeline/hits/rotation/rotation45off.rml @@ -0,0 +1,8 @@ + + + + + + + + diff --git a/pipeline/hits/rotation/rotation90.rml b/pipeline/hits/rotation/rotation90.rml new file mode 100644 index 00000000..193e0072 --- /dev/null +++ b/pipeline/hits/rotation/rotation90.rml @@ -0,0 +1,8 @@ + + + + + + + + diff --git a/pipeline/hits/rotation/specularXY.rml b/pipeline/hits/rotation/specularXY.rml new file mode 100644 index 00000000..0e0d7351 --- /dev/null +++ b/pipeline/hits/rotation/specularXY.rml @@ -0,0 +1,7 @@ + + + + + + + diff --git a/pipeline/hits/specular/specular.C b/pipeline/hits/specular/specular.C new file mode 100644 index 00000000..317a6ea9 --- /dev/null +++ b/pipeline/hits/specular/specular.C @@ -0,0 +1,67 @@ + +Int_t specular(Bool_t draw = false) { + TRestDetectorHitsEvent* ev = new TRestDetectorHitsEvent(); + + Double_t xCenter = 5; + Double_t yCenter = 10; + + TRandom3* random = new TRandom3(137); + for (int n = 0; n < 10000; n++) { + Double_t x = random->Gaus(xCenter, 3); + Double_t y = random->Gaus(yCenter, 3); + + ev->AddHit(x, y, 0, 1); + } + + // Initializing processes + string cfgFile1 = "specularXY.rml"; + TRestDetectorHitsSpecularProcess* specular1 = + new TRestDetectorHitsSpecularProcess((char*)cfgFile1.c_str()); + specular1->PrintMetadata(); + + cfgFile1 = "specularY.rml"; + TRestDetectorHitsSpecularProcess* specular2 = + new TRestDetectorHitsSpecularProcess((char*)cfgFile1.c_str()); + specular2->PrintMetadata(); + + cfgFile1 = "specularYoff.rml"; + TRestDetectorHitsSpecularProcess* specular3 = + new TRestDetectorHitsSpecularProcess((char*)cfgFile1.c_str()); + specular3->PrintMetadata(); + + TRestDetectorHitsEvent* specEvent1 = (TRestDetectorHitsEvent*)specular1->ProcessEvent(ev); + TRestDetectorHitsEvent* specEvent2 = (TRestDetectorHitsEvent*)specular2->ProcessEvent(ev); + TRestDetectorHitsEvent* specEvent3 = (TRestDetectorHitsEvent*)specular3->ProcessEvent(ev); + + /* Used to generate a combined plot */ + if (draw) { + TCanvas* c = new TCanvas("c0", "", 800, 600); + c->SetBorderSize(0.1); + // TPad *pad = new TPad("pad","The pad with the function", 0.05,0.50,0.95,0.95,21); + // pad->SetLeftMargin(0.15); + std::vector ranges{50, -25, 25, 50, -25, 25}; + TH2F* h = ev->GetXYHistogram(ranges); + h->Draw("col"); + h->Draw("cont3 same"); + c->Print("original.png"); + + h = specEvent1->GetXYHistogram(ranges); + // pad->cd(); + h->Draw("col"); + h->Draw("cont3 same"); + c->Print("specular1.png"); + + h = specEvent2->GetXYHistogram(ranges); + h->Draw("col"); + h->Draw("cont3 same"); + c->Print("specular2.png"); + + h = specEvent3->GetXYHistogram(ranges); + h->Draw("col"); + h->Draw("cont3 same"); + c->Print("specular3.png"); + } + + cout << "[\033[92m OK \x1b[0m]" << endl; + return 0; +} diff --git a/pipeline/hits/specular/specularXY.rml b/pipeline/hits/specular/specularXY.rml new file mode 100644 index 00000000..0e0d7351 --- /dev/null +++ b/pipeline/hits/specular/specularXY.rml @@ -0,0 +1,7 @@ + + + + + + + diff --git a/pipeline/hits/specular/specularY.rml b/pipeline/hits/specular/specularY.rml new file mode 100644 index 00000000..ce00b57e --- /dev/null +++ b/pipeline/hits/specular/specularY.rml @@ -0,0 +1,7 @@ + + + + + + + diff --git a/pipeline/hits/specular/specularYoff.rml b/pipeline/hits/specular/specularYoff.rml new file mode 100644 index 00000000..457b45a9 --- /dev/null +++ b/pipeline/hits/specular/specularYoff.rml @@ -0,0 +1,7 @@ + + + + + + + diff --git a/pipeline/hits/translation/specularXY.rml b/pipeline/hits/translation/specularXY.rml new file mode 100644 index 00000000..0e0d7351 --- /dev/null +++ b/pipeline/hits/translation/specularXY.rml @@ -0,0 +1,7 @@ + + + + + + + diff --git a/pipeline/hits/translation/translation.C b/pipeline/hits/translation/translation.C new file mode 100644 index 00000000..cb2d1554 --- /dev/null +++ b/pipeline/hits/translation/translation.C @@ -0,0 +1,42 @@ + +Int_t translation(Bool_t draw = false) { + TRestDetectorHitsEvent* ev = new TRestDetectorHitsEvent(); + + Double_t xCenter = 5; + Double_t yCenter = 10; + + TRandom3* random = new TRandom3(137); + for (int n = 0; n < 10000; n++) { + Double_t x = random->Gaus(xCenter, 3); + Double_t y = random->Gaus(yCenter, 3); + + ev->AddHit(x, y, 0, 1); + } + + // Initializing processes + string cfgFile1 = "translation.rml"; + TRestDetectorHitsTranslationProcess* translation = + new TRestDetectorHitsTranslationProcess((char*)cfgFile1.c_str()); + translation->PrintMetadata(); + + TRestDetectorHitsEvent* trEvent = (TRestDetectorHitsEvent*)translation->ProcessEvent(ev); + + /* Used to generate a combined plot */ + if (draw) { + TCanvas* c = new TCanvas("c0", "", 800, 600); + std::vector ranges{50, -25, 25, 50, -25, 25}; + TH2F* h = ev->GetXYHistogram(ranges); + h->Draw("col"); + h->Draw("cont3 same"); + c->Print("original.png"); + + h = trEvent->GetXYHistogram(ranges); + // pad->cd(); + h->Draw("col"); + h->Draw("cont3 same"); + c->Print("translation.png"); + } + + cout << "[\033[92m OK \x1b[0m]" << endl; + return 0; +} diff --git a/pipeline/hits/translation/translation.rml b/pipeline/hits/translation/translation.rml new file mode 100644 index 00000000..05f7d332 --- /dev/null +++ b/pipeline/hits/translation/translation.rml @@ -0,0 +1,4 @@ + + + + diff --git a/src/TRestDetectorHitmapAnalysisProcess.cxx b/src/TRestDetectorHitmapAnalysisProcess.cxx new file mode 100644 index 00000000..30948ff8 --- /dev/null +++ b/src/TRestDetectorHitmapAnalysisProcess.cxx @@ -0,0 +1,358 @@ +/************************************************************************* + * 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. * + *************************************************************************/ + +////////////////////////////////////////////////////////////////////////// +/// +/// This process allows to transform the detector hits mean positions using +/// rotation, translation and specular transformations. +/// +/// This process will not modify the detector hits inside the event. It +/// will just create new observables with the resulting transformed mean +/// positions. +/// +/// The process allows to include different definitions: +/// - **specular**: It requires a `position` and a `vector` that emplace the plane +/// position and orientation used for the specular transformation. +/// - **rotation**: It requires a `position` that defines the rotation center, +/// a `vector` that defines the rotation axis, and an `angle` that defines the +/// rotation amount. +/// - **translation**: It only requires to define the `vector` field with the +/// translation to be applied. +/// +/// The different transformations can be defined inside the RML section as +/// illustrated in the following example: +/// +/// \code +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// \endcode +/// +/// \htmlonly \endhtmlonly +/// ![The effect of different transfromation processes](transformedHitmap.png) +/// +/// The transformations to be applied, and the order in which those transformations will +/// be applied are defined inside the parameter `transformations` that contains a list +/// with the names from the transformation definitions. +/// +/// ### Observables added by this process +/// +/// * **xMean**: The mean X-position after the transformation +/// * **yMean**: The mean Y-position after the transformation +/// * **zMean**: The mean Z-position after the transformation +/// +///______________________________________________________________________________ +/// +/// RESTsoft - Software for Rare Event Searches with TPCs +/// +/// History of developments: +/// +/// 2023-July: First implementation +/// +/// \class TRestDetectorHitmapAnalysisProcess +/// \author Javier Galan +/// +///______________________________________________________________________________ +/// +////////////////////////////////////////////////////////////////////////// + +#include "TRestDetectorHitmapAnalysisProcess.h" + +using namespace std; + +ClassImp(TRestDetectorHitmapAnalysisProcess); + +/////////////////////////////////////////////// +/// \brief Default constructor +/// +TRestDetectorHitmapAnalysisProcess::TRestDetectorHitmapAnalysisProcess() { Initialize(); } + +/////////////////////////////////////////////// +/// \brief Default destructor +/// +TRestDetectorHitmapAnalysisProcess::~TRestDetectorHitmapAnalysisProcess() {} + +/////////////////////////////////////////////// +/// \brief Constructor loading data from a config file +/// +/// If no configuration path is defined using TRestMetadata::SetConfigFilePath +/// the path to the config file must be specified using full path, absolute or +/// relative. +/// +/// The default behaviour is that the config file must be specified with +/// full path, absolute or relative. +/// +/// \param configFilename A const char* giving the path to an RML file. +/// +TRestDetectorHitmapAnalysisProcess::TRestDetectorHitmapAnalysisProcess(const char* configFilename) { + Initialize(); + LoadConfigFromFile(configFilename); +} + +/////////////////////////////////////////////// +/// \brief Function to initialize input/output event members and define the section name +/// +void TRestDetectorHitmapAnalysisProcess::Initialize() { + SetSectionName(this->ClassName()); + SetLibraryVersion(LIBRARY_VERSION); + + fHitsEvent = nullptr; +} + +/////////////////////////////////////////////// +/// \brief The main processing event function +/// +TRestEvent* TRestDetectorHitmapAnalysisProcess::ProcessEvent(TRestEvent* inputEvent) { + fHitsEvent = (TRestDetectorHitsEvent*)inputEvent; + + TVector3 hitMean = fHitsEvent->GetMeanPosition(); + + TVector3 transformedMean = hitMean; + + for (const auto& t : fTransformations) { + HitTransformation tr; + tr = GetTransformation(t); + + if (tr.type == "specular") transformedMean = Specular(transformedMean, tr); + if (tr.type == "rotation") transformedMean = Rotation(transformedMean, tr); + if (tr.type == "translation") transformedMean = Translation(transformedMean, tr); + } + + SetObservableValue("xMean", transformedMean.X()); + SetObservableValue("yMean", transformedMean.Y()); + SetObservableValue("zMean", transformedMean.Z()); + + return fHitsEvent; +} + +/////////////////////////////////////////////// +/// \brief It performs a specular transformation of the given `pos` in the argument and the transformation +/// properties stored in the HitTransformation. +/// +TVector3 TRestDetectorHitmapAnalysisProcess::Specular(const TVector3& pos, const HitTransformation& tr) { + if (tr.type != "specular") return {0, 0, 0}; + + TVector3 V = pos - tr.position; + Double_t vByN = V.Dot(tr.vector); + TVector3 reflectionVector = V - 2 * vByN * tr.vector; + return tr.position + reflectionVector; +} + +/////////////////////////////////////////////// +/// \brief It performs a rotation of the given `pos` in the argument and the transformation properties +/// stored in the HitTransformation. +/// +TVector3 TRestDetectorHitmapAnalysisProcess::Rotation(const TVector3& pos, const HitTransformation& tr) { + if (tr.type != "rotation") return {0, 0, 0}; + + TVector3 position = pos - tr.position; + position.Rotate(tr.angle, tr.vector); + position += tr.position; + + return position; +} + +/////////////////////////////////////////////// +/// \brief It performs a translation of the given `pos` in the argument and the transformation properties +/// stored in the HitTransformation. +/// +TVector3 TRestDetectorHitmapAnalysisProcess::Translation(const TVector3& pos, const HitTransformation& tr) { + if (tr.type != "translation") return {0, 0, 0}; + + TVector3 position = pos + tr.vector; + return position; +} + +/////////////////////////////////////////////// +/// \brief It returns the transformation structure that matches the name in the argument +/// +HitTransformation TRestDetectorHitmapAnalysisProcess::GetTransformation(const std::string& name) { + for (const auto& t : fTransDefinitions) { + if (t.name == name) return t; + } + + HitTransformation Tr; + return Tr; +} + +/////////////////////////////////////////////// +/// \brief A custom initalization from a config file +/// +void TRestDetectorHitmapAnalysisProcess::InitFromConfigFile() { + TRestEventProcess::InitFromConfigFile(); + + auto specularDef = GetElement("specular"); + while (specularDef) { + HitTransformation trans; + trans.type = "specular"; + + string name = GetFieldValue("name", specularDef); + if (name == "Not defined") { + RESTError << "The `name` field was not found in specular definition!" << RESTendl; + exit(1); + } + trans.name = name; + + TVector3 position = Get3DVectorParameterWithUnits("position", specularDef); + if (position == TVector3(-1, -1, -1)) { + RESTError << "The `position` field was not found in specular definition!" << RESTendl; + exit(1); + } + trans.position = position; + + TVector3 vect = Get3DVectorParameterWithUnits("vector", specularDef); + if (vect == TVector3(-1, -1, -1)) { + RESTError << "The `vector` field was not found in specular definition!" << RESTendl; + exit(1); + } + trans.vector = vect; + + fTransDefinitions.push_back(trans); + specularDef = GetNextElement(specularDef); + } + + auto rotationDef = GetElement("rotation"); + while (rotationDef) { + HitTransformation trans; + trans.type = "rotation"; + + string name = GetFieldValue("name", rotationDef); + if (name == "Not defined") { + RESTError << "The `name` field was not found in rotation definition!" << RESTendl; + exit(1); + } + trans.name = name; + + TVector3 position = Get3DVectorParameterWithUnits("position", rotationDef); + if (position == TVector3(-1, -1, -1)) { + RESTError << "The `position` field was not found in rotation definition!" << RESTendl; + exit(1); + } + trans.position = position; + + TVector3 vect = Get3DVectorParameterWithUnits("vector", rotationDef); + if (vect == TVector3(-1, -1, -1)) { + RESTError << "The `vector` field was not found in rotation definition!" << RESTendl; + exit(1); + } + trans.vector = vect; + + Double_t angle = GetDblParameterWithUnits("angle", rotationDef); + if (angle == PARAMETER_NOT_FOUND_DBL) { + RESTError << "The `angle` field was not found in rotation definition!" << RESTendl; + exit(1); + } + trans.angle = angle; + + fTransDefinitions.push_back(trans); + rotationDef = GetNextElement(rotationDef); + } + + auto translationDef = GetElement("translation"); + while (translationDef) { + HitTransformation trans; + trans.type = "translation"; + + string name = GetFieldValue("name", translationDef); + if (name == "Not defined") { + RESTError << "The `name` field was not found in translation definition!" << RESTendl; + exit(1); + } + trans.name = name; + + TVector3 vect = Get3DVectorParameterWithUnits("vector", translationDef); + if (vect == TVector3(-1, -1, -1)) { + RESTError << "The `vector` field was not found in translation definition!" << RESTendl; + exit(1); + } + trans.vector = vect; + + fTransDefinitions.push_back(trans); + translationDef = GetNextElement(translationDef); + } +} + +/////////////////////////////////////////////// +/// \brief Prints out the metadata member values +/// +void TRestDetectorHitmapAnalysisProcess::PrintMetadata() { + BeginPrintProcess(); + + if (!fTransDefinitions.empty()) { + RESTMetadata << " List of transformation definitions: " << RESTendl; + RESTMetadata << " ----------------------------------- " << RESTendl; + RESTMetadata << " " << RESTendl; + } else + RESTMetadata << " No transformations found! " << RESTendl; + + for (const auto& t : fTransDefinitions) { + if (t.type == "specular") { + RESTMetadata << " - Specular transformation: " << RESTendl; + RESTMetadata << " + Plane position : (" << t.position.X() << ", " << t.position.Y() << ", " + << t.position.Z() << ") mm" << RESTendl; + RESTMetadata << " + Plane normal : (" << t.vector.X() << ", " << t.vector.Y() << ", " + << t.vector.Z() << ")" << RESTendl; + RESTMetadata << " " << RESTendl; + } + + if (t.type == "rotation") { + RESTMetadata << " - Rotation transformation: " << RESTendl; + RESTMetadata << " + Rotation center : (" << t.position.X() << ", " << t.position.Y() << ", " + << t.position.Z() << ") mm" << RESTendl; + RESTMetadata << " + Rotation axis : (" << t.vector.X() << ", " << t.vector.Y() << ", " + << t.vector.Z() << ")" << RESTendl; + RESTMetadata << " + Rotation angle : " << t.angle * units("degrees") << " degrees" << RESTendl; + RESTMetadata << " " << RESTendl; + } + + if (t.type == "translation") { + RESTMetadata << " - Translation transformation: " << RESTendl; + RESTMetadata << " + Translation vector : (" << t.vector.X() << ", " << t.vector.Y() << ", " + << t.vector.Z() << ") mm" << RESTendl; + RESTMetadata << " " << RESTendl; + } + } + + RESTMetadata << " Order of transformations: " << RESTendl; + RESTMetadata << " ------------------------- " << RESTendl; + RESTMetadata << " " << RESTendl; + + int n = 0; + for (const auto& t : fTransformations) { + n++; + RESTMetadata << n << ". " << t << RESTendl; + } + + EndPrintProcess(); +} diff --git a/src/TRestDetectorHitsEvent.cxx b/src/TRestDetectorHitsEvent.cxx index 0367c099..08f1340a 100644 --- a/src/TRestDetectorHitsEvent.cxx +++ b/src/TRestDetectorHitsEvent.cxx @@ -202,6 +202,7 @@ TRestHits* TRestDetectorHitsEvent::GetYZHits() { /// A XYZ hit compatible are those hits that have valid X, Y and Z coordinates. /// /// \return It returns back a TRestHits structure with the hits fulfilling the XYZ condition. +/// TRestHits* TRestDetectorHitsEvent::GetXYZHits() { fXYZHits->RemoveHits(); @@ -876,7 +877,7 @@ void TRestDetectorHitsEvent::DrawHistograms(Int_t& column, const TString& histOp 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.75); + fXZHisto->GetYaxis()->SetTitleOffset(1); } if (nYZ > 0) { @@ -888,7 +889,7 @@ void TRestDetectorHitsEvent::DrawHistograms(Int_t& column, const TString& histOp 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.75); + fYZHisto->GetYaxis()->SetTitleOffset(1); } if (nXY > 0) { @@ -900,7 +901,7 @@ void TRestDetectorHitsEvent::DrawHistograms(Int_t& column, const TString& histOp 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.75); + fXYHisto->GetYaxis()->SetTitleOffset(1); } column++; @@ -914,7 +915,7 @@ void TRestDetectorHitsEvent::DrawHistograms(Int_t& column, const TString& histOp 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.75); + fXHisto->GetYaxis()->SetTitleOffset(1); } if (nY > 0) { @@ -926,7 +927,7 @@ void TRestDetectorHitsEvent::DrawHistograms(Int_t& column, const TString& histOp 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.75); + fYHisto->GetYaxis()->SetTitleOffset(1); } if (nZ > 0) { @@ -938,7 +939,7 @@ void TRestDetectorHitsEvent::DrawHistograms(Int_t& column, const TString& histOp 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.75); + fZHisto->GetYaxis()->SetTitleOffset(1); } column++; @@ -958,3 +959,264 @@ void TRestDetectorHitsEvent::PrintEvent(Int_t nHits) const { fHits->PrintHits(nHits); } + +/////////////////////////////////////////////// +/// \brief This method draws the hits found on XY as a TH2F and it returns the +/// generated histogram. +/// +/// The first argument allows to define the ranges in the format: +/// {binsX, minX, maxX, binsY, minY, maxY}. +/// +/// The pitch allows to adjust the binning in automatic mode (when ranges vector is empty). +/// +/// The border gives an additional margin to the automatic minX, maxX, minY, maxY found. +/// +TH2F* TRestDetectorHitsEvent::GetXYHistogram(std::vector 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 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; +} + +/////////////////////////////////////////////// +/// \brief This method draws the hits found on XY as a TH2F and it returns the +/// generated histogram. +/// +/// The first argument allows to define the ranges in the format: +/// {binsX, minX, maxX, binsZ, minZ, maxZ}. +/// +/// The pitch allows to adjust the binning in automatic mode (when ranges vector is empty). +/// +/// The border gives an additional margin to the automatic minX, maxX, minZ, maxZ found. +/// +TH2F* TRestDetectorHitsEvent::GetXZHistogram(std::vector 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 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; +} + +/////////////////////////////////////////////// +/// \brief This method draws the hits found on YZ as a TH2F and it returns the +/// generated histogram. +/// +/// The first argument allows to define the ranges in the format: +/// {binsY, minY, maxY, binsZ, minZ, maxZ}. +/// +/// The pitch allows to adjust the binning in automatic mode (when ranges vector is empty). +/// +/// The border gives an additional margin to the automatic minY, maxY, minZ, maxZ found. +/// +TH2F* TRestDetectorHitsEvent::GetYZHistogram(std::vector 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 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; +} diff --git a/src/TRestDetectorHitsRotateAndTranslateProcess.cxx b/src/TRestDetectorHitsRotateAndTranslateProcess.cxx deleted file mode 100644 index 6509c46a..00000000 --- a/src/TRestDetectorHitsRotateAndTranslateProcess.cxx +++ /dev/null @@ -1,129 +0,0 @@ -///______________________________________________________________________________ -///______________________________________________________________________________ -///______________________________________________________________________________ -/// -/// -/// RESTSoft : Software for Rare Event Searches with TPCs -/// -/// TRestDetectorHitsRotateAndTranslateProcess.cxx -/// -/// Template to use to design "event process" classes inherited from -/// TRestDetectorHitsRotateAndTranslateProcess -/// How to use: replace TRestDetectorHitsRotateAndTranslateProcess by your -/// name, fill the required functions following instructions and add -/// all needed additional members and funcionality -/// -/// march 2016: First concept -/// Created as part of the conceptualization of existing REST -/// software. -/// Javier G. Garza -///_______________________________________________________________________________ - -////////////////////////////////////////////////////////////////////////// -/// -/// This process rotates and translates hits in a TRestDetectorHitsEvent. -/// Hits of type "VETO" are not affected by this process (we typically only want to rotate TPC hits). -/// -/// Usage: -/// -/// -/// -/// -/// -/// -/// - -#include "TRestDetectorHitsRotateAndTranslateProcess.h" - -using namespace std; - -ClassImp(TRestDetectorHitsRotateAndTranslateProcess); - -TRestDetectorHitsRotateAndTranslateProcess::TRestDetectorHitsRotateAndTranslateProcess() { Initialize(); } - -TRestDetectorHitsRotateAndTranslateProcess::TRestDetectorHitsRotateAndTranslateProcess( - const char* configFilename) { - Initialize(); - - if (LoadConfigFromFile(configFilename)) { - LoadDefaultConfig(); - } - - PrintMetadata(); -} - -TRestDetectorHitsRotateAndTranslateProcess::~TRestDetectorHitsRotateAndTranslateProcess() = default; - -void TRestDetectorHitsRotateAndTranslateProcess::LoadDefaultConfig() { SetTitle("Default config"); } - -void TRestDetectorHitsRotateAndTranslateProcess::Initialize() { - SetSectionName(this->ClassName()); - SetLibraryVersion(LIBRARY_VERSION); - - fInputHitsEvent = nullptr; - fOutputHitsEvent = nullptr; - - // Get volume name from parameter -} - -void TRestDetectorHitsRotateAndTranslateProcess::LoadConfig(const string& configFilename) { - if (LoadConfigFromFile(configFilename)) { - LoadDefaultConfig(); - } - - PrintMetadata(); -} - -void TRestDetectorHitsRotateAndTranslateProcess::InitProcess() {} - -TRestEvent* TRestDetectorHitsRotateAndTranslateProcess::ProcessEvent(TRestEvent* inputEvent) { - fInputHitsEvent = (TRestDetectorHitsEvent*)inputEvent; - fOutputHitsEvent = fInputHitsEvent; - - if (fOutputHitsEvent->GetNumberOfHits() == 0) { - return nullptr; - } - - for (unsigned int hit = 0; hit < fOutputHitsEvent->GetNumberOfHits(); hit++) { - const auto& type = fOutputHitsEvent->GetHits()->GetType(hit); - if (type == VETO) { - // Do not rotate VETO hits (we typically only rotate TPC hits) - continue; - } - fOutputHitsEvent->GetHits()->RotateIn3D(hit, fRotationVector.X(), fRotationVector.Y(), - fRotationVector.Z(), fRotationCenter); - fOutputHitsEvent->GetHits()->Translate(hit, fTranslationVector.X(), fTranslationVector.Y(), - fTranslationVector.Z()); - } - - return fOutputHitsEvent; -} - -void TRestDetectorHitsRotateAndTranslateProcess::EndProcess() {} - -void TRestDetectorHitsRotateAndTranslateProcess::InitFromConfigFile() { - fRotationCenter = Get3DVectorParameterWithUnits("rotationCenter", fRotationCenter); - fTranslationVector = Get3DVectorParameterWithUnits("translation", fTranslationVector); - - // rotation units should be specified in the rml (e.g. "90deg") - double rotationX = GetDblParameterWithUnits("rotationX", fRotationVector.X()); - double rotationY = GetDblParameterWithUnits("rotationY", fRotationVector.Y()); - double rotationZ = GetDblParameterWithUnits("rotationZ", fRotationVector.Z()); - fRotationVector = {rotationX, rotationY, rotationZ}; - - // legacy (maybe deprecated soon) - if (fTranslationVector.Mag2() == 0.0) { - fTranslationVector = { - GetDblParameterWithUnits("deltaX", 0.0), // - GetDblParameterWithUnits("deltaY", 0.0), // - GetDblParameterWithUnits("deltaZ", 0.0), // - }; - } - if (fRotationVector.Mag2() == 0.0) { - fRotationVector = { - GetDblParameterWithUnits("alpha", 0.0), // - GetDblParameterWithUnits("beta", 0.0), // - GetDblParameterWithUnits("gamma", 0.0), // - }; - } -} diff --git a/src/TRestDetectorHitsRotationProcess.cxx b/src/TRestDetectorHitsRotationProcess.cxx new file mode 100644 index 00000000..4f9bdd0b --- /dev/null +++ b/src/TRestDetectorHitsRotationProcess.cxx @@ -0,0 +1,146 @@ +/************************************************************************* + * 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. * + *************************************************************************/ + +////////////////////////////////////////////////////////////////////////// +/// TRestDetectorHitsRotationProcess will add a rotation of the x,y,z hits +/// in the space. The rotation will be done respect to a center and an +/// axis provided by the user. +/// +/// \code +/// +/// +/// +/// +/// +/// \endcode +/// +/// \warning This process will only be valid for those detector hits +/// of type XYZ, where the (x,y,z) coordinates take simultaneously physical +/// values. When the hits contain a non-valid (undetermined) coordinate +/// component this transformation will have non-sense and this process +/// will do nothing. +/// +/// The following figure has been produced using the `rotation.C` defined +/// under `detector/pipeline/hits/rotation/`. On the left we have the +/// original hits event distribution. On the middle the hits +/// with a 90 degrees clock-wise rotation with center at (0,0,0). On the +/// right the original hits after a clock-wise rotation of 45 degrees with +/// center at (-25,25,0). The axis for both rotations is along the Z-axis. +/// +/// \htmlonly \endhtmlonly +/// ![The effect of different rotation processes](rotation.png) +/// +///-------------------------------------------------------------------------- +/// +/// RESTsoft - Software for Rare Event Searches with TPCs +/// +/// History of developments: +/// +/// 2023-July: First implementation of TRestDetectorHitsRotationProcess +/// \author Javier Galan +/// +/// \class TRestDetectorHitsRotationProcess +/// +///
+/// +#include "TRestDetectorHitsRotationProcess.h" + +using namespace std; + +ClassImp(TRestDetectorHitsRotationProcess); + +TRestDetectorHitsRotationProcess::TRestDetectorHitsRotationProcess() { Initialize(); } + +/////////////////////////////////////////////// +/// \brief Constructor loading data from a config file +/// +/// If no configuration path is defined using TRestMetadata::SetConfigFilePath +/// the path to the config file must be specified using full path, absolute or +/// relative. +/// +/// The default behaviour is that the config file must be specified with +/// full path, absolute or relative. +/// +/// \param configFilename A const char* giving the path to an RML file. +/// +TRestDetectorHitsRotationProcess::TRestDetectorHitsRotationProcess(const char* configFilename) { + Initialize(); + LoadConfigFromFile(configFilename); +} + +TRestDetectorHitsRotationProcess::~TRestDetectorHitsRotationProcess() {} + +void TRestDetectorHitsRotationProcess::Initialize() { + SetSectionName(this->ClassName()); + SetLibraryVersion(LIBRARY_VERSION); + + fInputEvent = nullptr; + fOutputEvent = new TRestDetectorHitsEvent(); +} + +TRestEvent* TRestDetectorHitsRotationProcess::ProcessEvent(TRestEvent* inputEvent) { + fInputEvent = (TRestDetectorHitsEvent*)inputEvent; + fOutputEvent->SetEventInfo(fInputEvent); + + for (unsigned int hit = 0; hit < fInputEvent->GetNumberOfHits(); hit++) { + TVector3 position = {fInputEvent->GetX(hit), fInputEvent->GetY(hit), fInputEvent->GetZ(hit)}; + const auto type = fInputEvent->GetType(hit); + const auto energy = fInputEvent->GetEnergy(hit); + const auto time = fInputEvent->GetTime(hit); + + if (type != XYZ) { + fOutputEvent->AddHit(position, energy, time, type); + continue; + } + + // veto hits won't be rotated because they are not XYZ + + position -= fCenter; + position.Rotate(fAngle, fAxis); + position += fCenter; + + fOutputEvent->AddHit(position, energy, time, type); + } + + return fOutputEvent; +} + +void TRestDetectorHitsRotationProcess::InitFromConfigFile() { + TRestEventProcess::InitFromConfigFile(); + + fAxis = Get3DVectorParameterWithUnits("axis", fAxis); + fAxis.Unit(); + fCenter = Get3DVectorParameterWithUnits("center", fCenter); +} + +void TRestDetectorHitsRotationProcess::PrintMetadata() { + BeginPrintProcess(); + + RESTMetadata << " - Rotation center : ( " << fCenter.X() << ", " << fCenter.Y() << ", " << fCenter.Z() + << ") mm" << RESTendl; + RESTMetadata << " - Rotation axis : ( " << fAxis.X() << ", " << fAxis.Y() << ", " << fAxis.Z() << ")" + << RESTendl; + RESTMetadata << " - Rotation angle : " << fAngle * units("degrees") << " degrees" << RESTendl; + + EndPrintProcess(); +} diff --git a/src/TRestDetectorHitsSpecularProcess.cxx b/src/TRestDetectorHitsSpecularProcess.cxx new file mode 100644 index 00000000..4bc93ac8 --- /dev/null +++ b/src/TRestDetectorHitsSpecularProcess.cxx @@ -0,0 +1,144 @@ +/************************************************************************* + * 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. * + *************************************************************************/ + +////////////////////////////////////////////////////////////////////////// +/// TRestDetectorHitsSpecularProcess will update the hits position using +/// a specular image respect to a plane defined using its `normal` and +/// a `position` that is contained in the plane. +/// +/// \code +/// +/// +/// +/// +/// \endcode +/// +/// \warning This process will only be valid for those detector hits +/// of type XYZ, where the (x,y,z) coordinates take simultaneously physical +/// values. When the hits contain a non-valid (undetermined) coordinate +/// component this transformation will have non-sense and this process +/// will do nothing. +/// +/// The following figure has been produced using the `specular.C` defined +/// under `detector/pipeline/hits/specular/`. On the left-top corner we +/// have the original hits event distribution. On the right-top the hits +/// using a specular plane vector defined as (1,1,0) (or y=-x). On the +/// left-bottom a specular image respect to the X-axis, using the normal +/// plane vector (0,1,0). On the right-bottom image the same axis is used +/// but the plane position has been shifted by 5mm to the bottom. Black +/// lines identify the specular planes. +/// +/// \htmlonly \endhtmlonly +/// ![The effect of different specular processes](specular.png) +/// +///-------------------------------------------------------------------------- +/// +/// RESTsoft - Software for Rare Event Searches with TPCs +/// +/// History of developments: +/// +/// 2023-July: First implementation of TRestDetectorHitsSpecularProcess +/// \author Javier Galan +/// +/// \class TRestDetectorHitsSpecularProcess +/// +///
+/// +#include "TRestDetectorHitsSpecularProcess.h" + +using namespace std; + +ClassImp(TRestDetectorHitsSpecularProcess); + +TRestDetectorHitsSpecularProcess::TRestDetectorHitsSpecularProcess() { Initialize(); } + +/////////////////////////////////////////////// +/// \brief Constructor loading data from a config file +/// +/// If no configuration path is defined using TRestMetadata::SetConfigFilePath +/// the path to the config file must be specified using full path, absolute or +/// relative. +/// +/// The default behaviour is that the config file must be specified with +/// full path, absolute or relative. +/// +/// \param configFilename A const char* giving the path to an RML file. +/// +TRestDetectorHitsSpecularProcess::TRestDetectorHitsSpecularProcess(const char* configFilename) { + Initialize(); + LoadConfigFromFile(configFilename); +} + +TRestDetectorHitsSpecularProcess::~TRestDetectorHitsSpecularProcess() {} + +void TRestDetectorHitsSpecularProcess::Initialize() { + SetSectionName(this->ClassName()); + SetLibraryVersion(LIBRARY_VERSION); + + fInputEvent = nullptr; + fOutputEvent = new TRestDetectorHitsEvent(); +} + +TRestEvent* TRestDetectorHitsSpecularProcess::ProcessEvent(TRestEvent* inputEvent) { + fInputEvent = (TRestDetectorHitsEvent*)inputEvent; + fOutputEvent->SetEventInfo(fInputEvent); + + for (unsigned int hit = 0; hit < fInputEvent->GetNumberOfHits(); hit++) { + TVector3 position(fInputEvent->GetX(hit), fInputEvent->GetY(hit), fInputEvent->GetZ(hit)); + const auto type = fInputEvent->GetType(hit); + const auto energy = fInputEvent->GetEnergy(hit); + const auto time = fInputEvent->GetTime(hit); + + if (type != XYZ) { + fOutputEvent->AddHit(position, energy, time, type); + continue; + } + + const TVector3 relativePosition = position - fPosition; + const TVector3 reflectionVector = relativePosition - 2 * relativePosition.Dot(fNormal) * fNormal; + + position = fPosition + reflectionVector; + + fOutputEvent->AddHit(position, energy, time, type); + } + return fOutputEvent; +} + +void TRestDetectorHitsSpecularProcess::InitFromConfigFile() { + TRestEventProcess::InitFromConfigFile(); + + fNormal = Get3DVectorParameterWithUnits("normal", {0, 0, 1}); + fNormal = fNormal.Unit(); + fPosition = Get3DVectorParameterWithUnits("position", {0, 0, 1}); +} + +void TRestDetectorHitsSpecularProcess::PrintMetadata() { + BeginPrintProcess(); + + RESTMetadata << " - Plane normal : ( " << fNormal.X() << ", " << fNormal.Y() << ", " << fNormal.Z() + << ") mm" << RESTendl; + RESTMetadata << " - Plane position : ( " << fPosition.X() << ", " << fPosition.Y() << ", " + << fPosition.Z() << ")" << RESTendl; + + EndPrintProcess(); +} diff --git a/src/TRestDetectorHitsTranslationProcess.cxx b/src/TRestDetectorHitsTranslationProcess.cxx new file mode 100644 index 00000000..a5d1e1d3 --- /dev/null +++ b/src/TRestDetectorHitsTranslationProcess.cxx @@ -0,0 +1,128 @@ +/************************************************************************* + * 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. * + *************************************************************************/ + +////////////////////////////////////////////////////////////////////////// +/// TRestDetectorHitsTranslationProcess will move every hit inside a +/// TRestDetectorHitsEvent by an amount given at the `translation` vector. +/// +/// \code +/// +/// \endcode +/// +/// \warning This process will only be valid for those detector hits +/// of type XYZ, where the (x,y,z) coordinates take simultaneously physical +/// values. When the hits contain a non-valid (undetermined) coordinate +/// component this transformation will have non-sense and this process +/// will do nothing. +/// +/// The following figure has been produced using the `translation.C` +/// defined under `detector/pipeline/hits/translation/`. On the left we +/// have the original hits event distribution. On the right the hits +/// distribution after a translation by a vector of (-1.5,0.5,0)cm. +/// +/// \htmlonly \endhtmlonly +/// ![The effect of translation process](translation.png) +/// +///-------------------------------------------------------------------------- +/// +/// RESTsoft - Software for Rare Event Searches with TPCs +/// +/// History of developments: +/// +/// 2023-July: First implementation of TRestDetectorHitsTranslationProcess +/// \author Javier Galan +/// +/// \class TRestDetectorHitsTranslationProcess +/// +///
+/// +#include "TRestDetectorHitsTranslationProcess.h" + +using namespace std; + +ClassImp(TRestDetectorHitsTranslationProcess); + +TRestDetectorHitsTranslationProcess::TRestDetectorHitsTranslationProcess() { Initialize(); } + +/////////////////////////////////////////////// +/// \brief Constructor loading data from a config file +/// +/// If no configuration path is defined using TRestMetadata::SetConfigFilePath +/// the path to the config file must be specified using full path, absolute or +/// relative. +/// +/// The default behaviour is that the config file must be specified with +/// full path, absolute or relative. +/// +/// \param configFilename A const char* giving the path to an RML file. +/// +TRestDetectorHitsTranslationProcess::TRestDetectorHitsTranslationProcess(const char* configFilename) { + Initialize(); + LoadConfigFromFile(configFilename); +} + +TRestDetectorHitsTranslationProcess::~TRestDetectorHitsTranslationProcess() {} + +void TRestDetectorHitsTranslationProcess::Initialize() { + SetSectionName(this->ClassName()); + SetLibraryVersion(LIBRARY_VERSION); + + fInputEvent = nullptr; + fOutputEvent = new TRestDetectorHitsEvent(); +} + +TRestEvent* TRestDetectorHitsTranslationProcess::ProcessEvent(TRestEvent* inputEvent) { + fInputEvent = (TRestDetectorHitsEvent*)inputEvent; + fOutputEvent->SetEventInfo(fInputEvent); + + for (unsigned int hit = 0; hit < fInputEvent->GetNumberOfHits(); hit++) { + TVector3 position(fInputEvent->GetX(hit), fInputEvent->GetY(hit), fInputEvent->GetZ(hit)); + const auto type = fInputEvent->GetType(hit); + const auto energy = fInputEvent->GetEnergy(hit); + const auto time = fInputEvent->GetTime(hit); + + if (type != XYZ) { + fOutputEvent->AddHit(position, energy, time, type); + continue; + } + + position += fTranslation; + fOutputEvent->AddHit(position, energy, time, type); + } + return fOutputEvent; +} + +void TRestDetectorHitsTranslationProcess::InitFromConfigFile() { + TRestEventProcess::InitFromConfigFile(); + + fTranslation = Get3DVectorParameterWithUnits("translation", {0, 0, 1}); +} + +void TRestDetectorHitsTranslationProcess::PrintMetadata() { + BeginPrintProcess(); + + RESTMetadata << " - Translation vector : ( " << fTranslation.X() << ", " << fTranslation.Y() << ", " + << fTranslation.Z() << ") mm" << RESTendl; + + EndPrintProcess(); +}