126#include "TRestDetectorSignalToHitsProcess.h"
128#include "TRestDetectorSetup.h"
190 fGas = GetMetadata<TRestDetectorGas>();
191 if (
fGas !=
nullptr) {
193 RESTError <<
"A TRestDetectorGas definition was found but REST was not linked to Garfield libraries."
196 <<
"Please, remove the TRestDetectorGas definition, and add gas parameters inside the process "
197 "TRestDetectorSignalToHitsProcess"
200 if (!this->
GetError()) this->
SetError(
"Attempt to use TRestDetectorGas without Garfield");
219 this->
SetError(
"Drift velocity is negative.");
224 fReadout = GetMetadata<TRestDetectorReadout>();
228 this->
SetError(
"The readout was not properly initialized.");
248 RESTDebug <<
"TRestDetectorSignalToHitsProcess. Event id : " <<
fHitsEvent->GetID() <<
RESTendl;
253 Int_t numberOfSignals =
fSignalEvent->GetNumberOfSignals();
255 if (numberOfSignals == 0)
return nullptr;
257 Int_t planeID, readoutChannel = -1, readoutModule;
258 for (
int i = 0; i < numberOfSignals; i++) {
260 Int_t signalID = signal->GetSignalID();
263 cout <<
"Searching readout coordinates for signal ID : " << signalID << endl;
265 fReadout->GetPlaneModuleChannel(signalID, planeID, readoutModule, readoutChannel);
267 if (readoutChannel == -1) {
277 Double_t fieldZDirection = plane->
GetNormal().Z();
280 Double_t x = plane->
GetX(readoutModule, readoutChannel);
281 Double_t y = plane->
GetY(readoutModule, readoutChannel);
283 REST_HitType type = XYZ;
285 if (TMath::IsNaN(x)) {
288 }
else if (TMath::IsNaN(y)) {
294 Double_t hitTime = signal->GetMaxPeakTime();
298 cout <<
"Distance to plane : " << distanceToPlane << endl;
300 Double_t z = zPosition + fieldZDirection * distanceToPlane;
302 Double_t energy = signal->GetMaxPeakValue();
305 cout <<
"Adding hit. Time : " << hitTime <<
" x : " << x <<
" y : " << y <<
" z : " << z
306 <<
" Energy : " << energy << endl;
309 }
else if (
fMethod ==
"tripleMax") {
310 Int_t bin = signal->GetMaxIndex();
311 int binprev = (bin - 1) < 0 ? bin : bin - 1;
312 int binnext = (bin + 1) > signal->GetNumberOfPoints() - 1 ? bin : bin + 1;
314 Double_t hitTime = signal->GetTime(bin);
315 Double_t energy = signal->GetData(bin);
318 Double_t z = zPosition + fieldZDirection * distanceToPlane;
322 hitTime = signal->GetTime(binprev);
323 energy = signal->GetData(binprev);
326 z = zPosition + fieldZDirection * distanceToPlane;
330 hitTime = signal->GetTime(binnext);
331 energy = signal->GetData(binnext);
334 z = zPosition + fieldZDirection * distanceToPlane;
339 cout <<
"Distance to plane : " << distanceToPlane << endl;
340 cout <<
"Adding hit. Time : " << hitTime <<
" x : " << x <<
" y : " << y <<
" z : " << z
341 <<
" Energy : " << energy << endl;
343 }
else if (
fMethod ==
"tripleMaxAverage") {
344 Int_t bin = signal->GetMaxIndex();
345 int binprev = (bin - 1) < 0 ? bin : bin - 1;
346 int binnext = (bin + 1) > signal->GetNumberOfPoints() - 1 ? bin : bin + 1;
348 Double_t hitTime = signal->GetTime(bin);
349 Double_t energy1 = signal->GetData(bin);
352 Double_t z1 = zPosition + fieldZDirection * distanceToPlane;
354 hitTime = signal->GetTime(binprev);
355 Double_t energy2 = signal->GetData(binprev);
358 Double_t z2 = zPosition + fieldZDirection * distanceToPlane;
360 hitTime = signal->GetTime(binnext);
361 Double_t energy3 = signal->GetData(binnext);
364 Double_t z3 = zPosition + fieldZDirection * distanceToPlane;
366 Double_t eTot = energy1 + energy2 + energy3;
368 Double_t zAvg = ((z1 * energy1) + (z2 * energy2) + (z3 * energy3)) / eTot;
369 Double_t eAvg = eTot / 3.0;
374 cout <<
"Distance to plane: " << distanceToPlane << endl;
375 cout <<
"Adding hit. Time : " << hitTime <<
" x : " << x <<
" y : " << y <<
" z : " << zAvg
376 <<
" Energy : " << eAvg << endl;
377 cout <<
"z1, z2, z3 = " << z1 <<
", " << z2 <<
", " << z3 << endl;
378 cout <<
"E1, E2, E3 = " << energy1 <<
", " << energy2 <<
", " << energy3 << endl;
381 }
else if (
fMethod ==
"gaussFit") {
382 TVector2 gaussFit = signal->GetMaxGauss();
384 Double_t hitTime = 0;
386 if (gaussFit.X() >= 0.0) {
387 hitTime = gaussFit.X();
389 z = zPosition + fieldZDirection * distanceToPlane;
391 Double_t energy = -1.0;
392 if (gaussFit.Y() >= 0.0) {
393 energy = gaussFit.Y();
397 cout <<
"Signal event : " << signal->GetSignalID()
398 <<
"--------------------------------------------------------" << endl;
399 cout <<
"GausFit : time bin " << gaussFit.X() <<
" and energy : " << gaussFit.Y() << endl;
400 cout <<
"Signal to hit info : zPosition : " << zPosition
401 <<
"; fieldZDirection : " << fieldZDirection <<
" and driftV : " <<
fDriftVelocity
403 cout <<
"Adding hit. Time : " << hitTime <<
" x : " << x <<
" y : " << y <<
" z : " << z
404 <<
" Energy : " << energy << endl;
409 }
else if (
fMethod ==
"landauFit") {
410 TVector2 landauFit = signal->GetMaxLandau();
412 Double_t hitTime = 0;
414 if (landauFit.X() >= 0.0) {
415 hitTime = landauFit.X();
417 z = zPosition + fieldZDirection * distanceToPlane;
419 Double_t energy = -1.0;
420 if (landauFit.Y() >= 0.0) {
421 energy = landauFit.Y();
425 cout <<
"Signal event : " << signal->GetSignalID()
426 <<
"--------------------------------------------------------" << endl;
427 cout <<
"landauFit : time bin " << landauFit.X() <<
" and energy : " << landauFit.Y() << endl;
428 cout <<
"Signal to hit info : zPosition : " << zPosition
429 <<
"; fieldZDirection : " << fieldZDirection <<
" and driftV : " <<
fDriftVelocity
431 cout <<
"Adding hit. Time : " << hitTime <<
" x : " << x <<
" y : " << y <<
" z : " << z
432 <<
" Energy : " << energy << endl;
437 }
else if (
fMethod ==
"agetFit") {
438 TVector2 agetFit = signal->GetMaxAget();
440 Double_t hitTime = 0;
442 if (agetFit.X() >= 0.0) {
443 hitTime = agetFit.X();
445 z = zPosition + fieldZDirection * distanceToPlane;
447 Double_t energy = -1.0;
448 if (agetFit.Y() >= 0.0) {
449 energy = agetFit.Y();
453 cout <<
"Signal event : " << signal->GetSignalID()
454 <<
"--------------------------------------------------------" << endl;
455 cout <<
"agetFit : time bin " << agetFit.X() <<
" and energy : " << agetFit.Y() << endl;
456 cout <<
"Signal to hit info : zPosition : " << zPosition
457 <<
"; fieldZDirection : " << fieldZDirection <<
" and driftV : " <<
fDriftVelocity
459 cout <<
"Adding hit. Time : " << hitTime <<
" x : " << x <<
" y : " << y <<
" z : " << z
460 <<
" Energy : " << energy << endl;
465 }
else if (
fMethod ==
"qCenter") {
466 Double_t energy_signal = 0;
467 Double_t distanceToPlane = 0;
469 for (
int j = 0; j < signal->GetNumberOfPoints(); j++) {
470 Double_t energy_point = signal->GetData(j);
471 energy_signal += energy_point;
472 distanceToPlane += signal->GetTime(j) *
fDriftVelocity * energy_point;
474 Double_t energy = energy_signal / signal->GetNumberOfPoints();
476 Double_t z = zPosition + fieldZDirection * (distanceToPlane / energy_signal);
479 for (
int j = 0; j < signal->GetNumberOfPoints(); j++) {
480 Double_t energy = signal->GetData(j);
485 cout <<
"Time : " << signal->GetTime(j) <<
" Drift velocity : " <<
fDriftVelocity << endl;
486 cout <<
"Distance to plane : " << distanceToPlane << endl;
489 Double_t z = zPosition + fieldZDirection * distanceToPlane;
492 cout <<
"Adding hit. Time : " << signal->GetTime(j) <<
" x : " << x <<
" y : " << y
493 <<
" z : " << z << endl;
497 }
else if (
fMethod ==
"intwindow") {
498 Int_t nPoints = signal->GetNumberOfPoints();
499 std::map<int, std::pair<int, double> > windowMap;
500 for (
int j = 0; j < nPoints; j++) {
501 int index = signal->GetTime(j) / fIntWindow;
502 auto it = windowMap.find(index);
503 if (it != windowMap.end()) {
505 it->second.second += signal->GetData(j);
507 windowMap[index] = std::make_pair(1, signal->GetData(j));
511 for (
const auto& [index, pair] : windowMap) {
512 Double_t hitTime = index * fIntWindow + fIntWindow / 2.;
513 Double_t energy = pair.second / pair.first;
514 if (energy < fThreshold) {
517 RESTDebug <<
"TimeBin " << index <<
" Time " << hitTime <<
" Charge: " << energy
518 <<
" Thr: " << (fThreshold) <<
RESTendl;
520 Double_t z = zPosition + fieldZDirection * distanceToPlane;
522 RESTDebug <<
"Time : " << hitTime <<
" Drift velocity : " <<
fDriftVelocity
523 <<
"\nDistance to plane : " << distanceToPlane <<
RESTendl;
524 RESTDebug <<
"Adding hit. Time : " << hitTime <<
" x : " << x <<
" y : " << y <<
" z : " << z
530 string errMsg =
"The method " + (string)
fMethod +
" is not implemented!";
535 RESTDebug <<
"TRestDetectorSignalToHitsProcess. Hits added : " <<
fHitsEvent->GetNumberOfHits()
537 RESTDebug <<
"TRestDetectorSignalToHitsProcess. Hits total energy : " <<
fHitsEvent->GetTotalEnergy()
548 ". Failed to find readout positions in channel to hit conversion.";
virtual void SetElectricField(double value)
Sets the electric field. Must be given in V/mm.
virtual Double_t GetElectricField() const
Returns the electric field in V/mm.
virtual void PrintEvent() const
void AddHit(Double_t x, Double_t y, Double_t z, Double_t en, Double_t t=0, REST_HitType type=XYZ)
Adds a new hit to this event.
TVector2 GetSize() const
Returns the module size (x, y) in mm.
TVector2 GetPlaneCoordinates(const TVector2 &p)
TRestDetectorReadoutModule * GetModuleByID(Int_t modID)
Returns a pointer to a module using its internal module id.
TVector3 GetNormal() const
Returns a TVector3 with a vector normal to the readout plane.
TVector3 GetPosition() const
Returns a TVector3 with the readout plane position.
Double_t GetX(Int_t modID, Int_t chID)
Returns the X coordinate of a given channel in a given module using their internal module and channel...
Double_t GetY(Int_t modID, Int_t chID)
Returns the Y coordinate of a given channel in a given module using their internal module and channel...
TRestDetectorReadoutPlane * GetReadoutPlaneWithID(int id)
Returns a pointer to the readout plane by ID.
A process to transform a daq channel and physical time to spatial coordinates.
Double_t fDriftVelocity
The drift velocity in standard REST units (mm/us).
void InitProcess() override
Process initialization. Observable names are interpreted and auxiliar observable members,...
TRestDetectorHitsEvent * fHitsEvent
A pointer to the specific TRestDetectorHitsEvent output.
TRestEvent * ProcessEvent(TRestEvent *inputEvent) override
The main processing event function.
~TRestDetectorSignalToHitsProcess() override
Default destructor.
Double_t fGasPressure
The gas pressure in atm. Only relevant if TRestDetectorGas is used.
TRestDetectorGas * fGas
A pointer to the detector gas definition accessible to TRestRun.
TRestDetectorSignalToHitsProcess()
Default constructor.
void Initialize() override
Function to initialize input/output event members and define the section name.
void LoadDefaultConfig()
Function to load the default config in absence of RML input.
Double_t fElectricField
The electric field in standard REST units (V/mm). Only relevant if TRestDetectorGas is used.
TString fMethod
The method used to transform the signal points to hits.
TRestDetectorSignalEvent * fSignalEvent
A pointer to the specific TRestDetectorHitsEvent input.
TRestDetectorReadout * fReadout
A pointer to the detector readout definition accessible to TRestRun.
A base class for any REST event.
@ REST_Extreme
show everything
@ REST_Debug
+show the defined debug messages
std::string IntegerToString(Int_t n, std::string format="%d")
Gets a string from an integer.