Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Implementation of the Tprev/Tnext task. #892

Merged
merged 1 commit into from
Nov 21, 2023
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions r3bbase/BaseLinkDef.h
Original file line number Diff line number Diff line change
Expand Up @@ -34,4 +34,5 @@
#pragma link C++ class R3BMSOffsetContFact+;
#pragma link C++ class R3BMSOffsetPar+;
#pragma link C++ class R3BMSOffsetFinder+;
#pragma link C++ class R3BTprevTnext+;
#endif
2 changes: 2 additions & 0 deletions r3bbase/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,7 @@ set(SRCS
R3BTcutPar.cxx
R3BTsplinePar.cxx
R3BWhiterabbitPropagator.cxx
R3BTprevTnext.cxx
./pars/R3BMSOffsetPar.cxx
./pars/R3BMSOffsetContFact.cxx
./pars/R3BMSOffsetFinder.cxx)
Expand All @@ -62,6 +63,7 @@ set(HEADERS
R3BTcutPar.h
R3BTsplinePar.h
R3BWhiterabbitPropagator.h
R3BTprevTnext.h
./pars/R3BMSOffsetPar.h
./pars/R3BMSOffsetContFact.h
./pars/R3BMSOffsetFinder.h)
Expand Down
156 changes: 156 additions & 0 deletions r3bbase/R3BTprevTnext.cxx
Original file line number Diff line number Diff line change
@@ -0,0 +1,156 @@
/******************************************************************************
* Copyright (C) 2019 GSI Helmholtzzentrum für Schwerionenforschung GmbH *
* Copyright (C) 2019-2023 Members of R3B Collaboration *
* *
* This software is distributed under the terms of the *
* GNU General Public Licence (GPL) version 3, *
* copied verbatim in the file "LICENSE". *
* *
* In applying this license GSI does not waive the privileges and immunities *
* granted to it by virtue of its status as an Intergovernmental Organization *
* or submit itself to any jurisdiction. *
******************************************************************************/

#include "R3BTprevTnext.h"
#include "FairLogger.h"
#include "FairRootManager.h"
#include "FairRunAna.h"
#include "FairRuntimeDb.h"
#include "R3BEventHeader.h"
#include "R3BMSOffsetPar.h"
#include "R3BSamplerMappedData.h"
#include "R3BShared.h"
#include "TClonesArray.h"
#include "TMath.h"
#include "TObjArray.h"
#include "TRandom.h"
#include <cstdlib>
#include <iostream>
R3BTprevTnext::R3BTprevTnext()
: R3BTprevTnext("R3B Tprev/Tnext", 1)
{
}

R3BTprevTnext::R3BTprevTnext(const TString& name, Int_t iVerbose)
: FairTask(name, iVerbose)
{
}

InitStatus R3BTprevTnext::Init()
{

FairRootManager* rootManager = FairRootManager::Instance();
if (rootManager == nullptr)
{
LOG(error) << "FairRootManager not opened!";
return kFATAL;
}

FairRuntimeDb* rtdb = FairRuntimeDb::instance();
if (rtdb == nullptr)
{
LOG(error) << "FairRuntimeDb not opened!";
return kFATAL;
}

fR3BEventHeader = dynamic_cast<R3BEventHeader*>(rootManager->GetObject("EventHeader."));
if (fR3BEventHeader == nullptr)
{
LOG(error) << "R3BTprevTnext::Init() EventHeader. not found";
return kFATAL;
}

fSamplerMapped = dynamic_cast<TClonesArray*>(rootManager->GetObject("SamplerMapped"));
if (fSamplerMapped == nullptr)
{
LOG(error) << "SamplerMapped not found or wrong type.";
return kFATAL;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would prefer if you did one of the following:

  • assert(fSamplerMapped)
  • if () { LOG(fatal)<< "SamplerMapped not found or wrong type.";}
  • if () {LOG(error) << "..."; return kFATAL; }

If you just return kFATAL, FairRoot is (hopefully) going to report which initialization of a Task failed, but not why exactly.

}

fSamplerMSMapped = dynamic_cast<TClonesArray*>(rootManager->GetObject("SamplerMSMapped"));
if (fSamplerMSMapped == nullptr)
{
LOG(error) << "SamplerMSMapped not found or wrong type.";
return kFATAL;
}

fMSOffsetPar = dynamic_cast<R3BMSOffsetPar*>(rtdb->getContainer("MSOffsetPar"));
if (fMSOffsetPar == nullptr)
{
LOG(error) << "Could not find MSOffsetPar container!";
return kFATAL;
}

return kSUCCESS;
}

void R3BTprevTnext::Exec(Option_t* /*opt*/)
{
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

constexpr auto NOTFOUND = -10.0;
constexpr double WRAP_PERIOD = (1L<<32)*CLOCK_PERIOD; //42.9 seconds
constexpr auto ONLYHIT = -20.0;

If you ever return these, I would make the constants a public part of the class interface. If you only use them internally, put them with CLOCK_PERIOD.

Double_t fTprev = INVALID_TPTN; // Initialization with invalid values.
Double_t fTnext = INVALID_TPTN;
Int_t dMScounter = 0;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

First, please add // number of sampler hits in coincidence with MS
Second, if we use this, sawMS should be removed and replaced by dMScounter>0.
Third, do we actually need to handle the special case where LOS saw two ions within 20ns? My feeling is no. In that case, running the algorithm as usual should result in Tprev or Tnext being really small, which will in turn make the users reject the hit as an obvious pileup.

const Int_t sampHits = fSamplerMapped->GetEntriesFast();
const Int_t sampmsHits = fSamplerMSMapped->GetEntriesFast();
const Double_t MSOffset = fMSOffsetPar->GetMSOffset();
Double_t SAMPTime = 0;
Double_t SAMPMSTime = 0;
if (sampmsHits == 1 && sampHits > 0)
{
R3BSamplerMappedData* SAMPMapped = nullptr;
R3BSamplerMappedData* SAMPMSMapped = nullptr;
SAMPMSMapped = dynamic_cast<R3BSamplerMappedData*>(fSamplerMSMapped->At(0));
SAMPMSTime = SAMPMSMapped->GetTime();
for (Int_t i = 0; i < sampHits; ++i)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe put this for-loop into another member function and give it a nice name indicating the meaning?

{
SAMPMapped = dynamic_cast<R3BSamplerMappedData*>(fSamplerMapped->At(i));
SAMPTime = SAMPMapped->GetTime();
auto tpn = SAMPTime - SAMPMSTime - MSOffset;
if (tpn < -fDelta_clk)
{
fTprev = -CLOCK_PERIOD * tpn;
}
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

While I am not sure that the above code is correct as is, I am sure that the code could be simplified.
First, initialize both fTnext and fTprev to invalid. bool sawMS=false;
Then, loop.
Calculate auto tpn = SAMPTime - SAMPMSTime - MSOffset; (so it increases).
As long as tpn is smaller than -1, upate fTprev to -10*tpn.
If abs(tpn)<=1, set a bool sawMS to one.
If tpn>1, set fTnext and break.

This way, iff fTPrev or fTNext or both are never set, they automatically contain the invalid value.

Also, if sawMS is false for more than a small fraction of events, then something terrible has happened (likely the offset is wrong). We should complain periodically with LOG(error).

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok, seems like a good idea, but again, I need some time to understand completely. But I will come with something

if (abs(tpn) <= fDelta_clk)
{
++dMScounter;
}
if (tpn > fDelta_clk)
{
fTnext = CLOCK_PERIOD * tpn;
break;
}
}
}
else
{

This comment was marked as resolved.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You are right

fR3BEventHeader->SetTprev(INVALID_EVENT);
fR3BEventHeader->SetTnext(INVALID_EVENT);
return;
}

if (dMScounter ==
0) // This if statement checks if there is no Hit correspondent to the MSHit even after the first if statement
{
fTprev = ERROR_NO_MS;
fTnext = ERROR_NO_MS;
}
else if (dMScounter > 1)
{
fTprev = ERROR_MULTI_MS;
fTnext = ERROR_MULTI_MS;
}
else if (dMScounter == 1)
{
if (fTprev == INVALID_TPTN)
{
fTprev = fNo_record_time;
}
if (fTnext == INVALID_TPTN)
{
fTnext = fNo_record_time;
}
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would replace the last if block with:

else if (dMScounter == 1)
{
  if (fTprev == INVALID_TPTN)
     fTprev = fOne_hit;
  if (fTnext == INVALID_TPTN)
     fTnext = fOne_hit;
}
fR3BEventHeader->SetTprev(fTprev);
fR3BEventHeader->SetTnext(fTnext);

That way, whenever there is no previous/next hit, the "NoRecord" value will be used. For the special case of one hit (no previous and no next hit), both of these will automatically be set to that value.

}
fR3BEventHeader->SetTprev(fTprev);
fR3BEventHeader->SetTnext(fTnext);
}

ClassImp(R3BTprevTnext) // NOLINT
57 changes: 57 additions & 0 deletions r3bbase/R3BTprevTnext.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,57 @@
/******************************************************************************
* Copyright (C) 2019 GSI Helmholtzzentrum für Schwerionenforschung GmbH *
* Copyright (C) 2019-2023 Members of R3B Collaboration *
* *
* This software is distributed under the terms of the *
* GNU General Public Licence (GPL) version 3, *
* copied verbatim in the file "LICENSE". *
* *
* In applying this license GSI does not waive the privileges and immunities *
* granted to it by virtue of its status as an Intergovernmental Organization *
* or submit itself to any jurisdiction. *
******************************************************************************/

#pragma once

#include "FairTask.h"
#include "R3BMSOffsetPar.h"
class TClonesArray;
class R3BMSOffsetPar;
class R3BEventHeader;
class R3BSamplerMappedData;
class R3BTprevTnext : public FairTask
{
public:
/** Default constructor **/
R3BTprevTnext();

/** Standard constructor **/
explicit R3BTprevTnext(const TString& name, Int_t iVerbose = 1);

/** Virtual method Init **/
InitStatus Init() override;

/** Virtual method Exec **/
void Exec(Option_t* opt) override;

void SetNoRecordTimeGap(Double_t No_Record_Time_Gap)
{
fNo_record_time = No_Record_Time_Gap;
} // Setter for a value of Tprev/Tnext if there is no recorded hit before and/or after the MS
void SetDelta_clk(Double_t clock) { fDelta_clk = clock; }

static constexpr auto CLOCK_PERIOD = 10; // ns
static constexpr auto ERROR_NO_MS = -30; // Assigned value for events where there is no MS
static constexpr auto INVALID_TPTN = -10; // Assigned value for invalid or nonexistent TPrev or TNext
static constexpr auto ERROR_MULTI_MS = -20; // Assigned value for events where there are multiple MS
static constexpr auto INVALID_EVENT = -40; // Assigned value for events where the MS was not correctly written
private:
R3BMSOffsetPar* fMSOffsetPar = nullptr;
TClonesArray* fSamplerMapped = nullptr;
TClonesArray* fSamplerMSMapped = nullptr;
R3BEventHeader* fR3BEventHeader = nullptr; /**< Event header - input data. */
Double_t fDelta_clk = 1.0;
Double_t fNo_record_time = 2e3;

ClassDefOverride(R3BTprevTnext, 1); // NOLINT
};
6 changes: 0 additions & 6 deletions r3bbase/pars/R3BMSOffsetFinder.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -43,12 +43,6 @@ R3BMSOffsetFinder::R3BMSOffsetFinder()
// R3BMSOffsetFinder: Standard Constructor --------------------------
R3BMSOffsetFinder::R3BMSOffsetFinder(const TString& name, Int_t iVerbose)
: FairTask(name, iVerbose)
, fMSOffset(0.0)
, fMinStatistics(1)
, fMSOffsetPar(nullptr)
, fSamplerMapped(nullptr)
, fSamplerMSMapped(nullptr)
, fh_Offset_Finder(nullptr)
{
}

Expand Down
12 changes: 6 additions & 6 deletions r3bbase/pars/R3BMSOffsetFinder.h
Original file line number Diff line number Diff line change
Expand Up @@ -69,16 +69,16 @@ class R3BMSOffsetFinder : public FairTask // NOLINT

private:
// Number of histograms, limits and bining
Double_t fMSOffset;
Double_t fMSOffset = 0.0;

// Minimum statistics and parameters
Int_t fMinStatistics{};
Int_t fMinStatistics = 1;

R3BMSOffsetPar* fMSOffsetPar; /**< Parameter container. >*/
TClonesArray* fSamplerMapped; /**< Array with SAMP Mapped input data. >*/
TClonesArray* fSamplerMSMapped; /**< Array with SAMPMS MApped input data. >*/
R3BMSOffsetPar* fMSOffsetPar = nullptr; /**< Parameter container. >*/
TClonesArray* fSamplerMapped = nullptr; /**< Array with SAMP Mapped input data. >*/
TClonesArray* fSamplerMSMapped = nullptr; /**< Array with SAMPMS MApped input data. >*/

TH1F* fh_Offset_Finder;
TH1F* fh_Offset_Finder = nullptr;

public:
ClassDefOverride(R3BMSOffsetFinder, 1); // NOLINT
Expand Down