Skip to content

Commit

Permalink
TwiM drift time calib fixed, suppress unnecessary cout
Browse files Browse the repository at this point in the history 
The drift time calibration for TwiM in s467/s444 (2020) has to be independent of the standard configuration since two MDPP16 were used for section 0 for those cases. The mapping of the time reference was wrong previously, and the R3BCoarseTimeStitch has been removed since it was causing fake reconstructions. Clang tidy for the TwimMapped2Cal and R3BCalifa applied to satisfy the requirements by the static-analysis in Github... New structures for handling Energy and Drift Time in Twim have been defined.
  • Loading branch information
@ryotani @jose-luis-rs
ryotani authored and jose-luis-rs committed Jul 16, 2023
1 parent 059771a commit b61b772
Show file tree
Hide file tree
Showing 6 changed files with 178 additions and 110 deletions.
13 changes: 9 additions & 4 deletions califa/sim/R3BCalifa.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -133,8 +133,11 @@ Bool_t R3BCalifa::ProcessHits(FairVolume* vol)
// stack when they crash (e.g. because of SIGFPE, which GEANT4 helpfully
// activates on Tuesdays and Debug builds (->G4FPE_debug).

auto makeCompFun = [](std::array<double, 5> p)
{ return [p](double x) { return (x > 0.0) ? 1. / (p[0] + p[1] * pow(x, p[2]) + p[3] / pow(x, p[4])) : 0.0; }; };
auto makeCompFun = [](std::array<double, 5> par)
{
return [par](double val)
{ return (val > 0.0) ? 1. / (par[0] + par[1] * pow(val, par[2]) + par[3] / pow(val, par[4])) : 0.0; };
};
auto tf_dNf_dE = makeCompFun({ -1.79, 1.36e-2, 7.84e-1, 4.97, 1.75e-1 });
auto tf_dNs_dE = makeCompFun({ -1.24e2, 6.3e-3, 1.27, 1.262e2, 2.3e-3 });

Expand Down Expand Up @@ -220,8 +223,10 @@ Bool_t R3BCalifa::ProcessHits(FairVolume* vol)

void R3BCalifa::EndOfEvent()
{
// if (fVerboseLevel > 1)
Print();
if (fVerboseLevel > 1)
{
Print();
}
fCalifaCollection->Clear();
ResetParameters();
}
Expand Down
6 changes: 3 additions & 3 deletions r3bdata/R3BMCStack.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -384,11 +384,11 @@ void R3BStack::Register() { FairRootManager::Instance()->Register("MCTrack", "St
// ----- Public method Print --------------------------------------------
void R3BStack::PrintStack(Int_t iVerbose) const
{
LOG(info) << "R3BStack: Number of primaries = " << fNPrimaries;
LOG(info) << " Total number of particles = " << fNParticles;
LOG(info) << " Number of tracks in output = " << fNTracks;
if (1 == iVerbose)
{
LOG(info) << "R3BStack: Number of primaries = " << fNPrimaries;
LOG(info) << " Total number of particles = " << fNParticles;
LOG(info) << " Number of tracks in output = " << fNTracks;
for (Int_t iTrack = 0; iTrack < fNTracks; iTrack++)
{
char str[100];
Expand Down
2 changes: 1 addition & 1 deletion r3bsource/trloii/R3BTrloiiTpatReader.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,7 +43,7 @@ R3BTrloiiTpatReader::R3BTrloiiTpatReader(EXT_STR_h101_TPAT* data, size_t offset)
{
}

R3BTrloiiTpatReader::~R3BTrloiiTpatReader() { R3BLOG(info, ""); }
R3BTrloiiTpatReader::~R3BTrloiiTpatReader() { R3BLOG(debug, ""); }

Bool_t R3BTrloiiTpatReader::Init(ext_data_struct_info* a_struct_info)
{
Expand Down
5 changes: 3 additions & 2 deletions r3bsource/twim/R3BTwimReader.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,13 +47,14 @@ class R3BTwimReader : public R3BReader
// Accessor to select events without pileup
void SetSkipPileup() { fPileup = kTRUE; }

private:
Bool_t ReadData(EXT_STR_h101_SOFTWIM_onion*, UShort_t);
void SetNumSections(Int_t num) { fSections = num; }
void SetNumAnodes(Int_t num) { fAnodes = num; }
void SetNumTref(Int_t num) { fTref = num; }
void SetNumTtrig(Int_t num) { fTtrig = num; }

private:
Bool_t ReadData(EXT_STR_h101_SOFTWIM_onion*, UShort_t);

/* Reader specific data structure from ucesb */
EXT_STR_h101_SOFTWIM* fData;
/* Data offset */
Expand Down
208 changes: 113 additions & 95 deletions twim/calibration/R3BTwimMapped2Cal.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,14 +28,16 @@
#include "FairRuntimeDb.h"

// TWIM headers
#include "R3BCoarseTimeStitch.h"
#include "R3BEventHeader.h"
#include "R3BLogger.h"
#include "R3BTwimCalData.h"
#include "R3BTwimCalPar.h"
#include "R3BTwimMapped2Cal.h"
#include "R3BTwimMappedData.h"

constexpr int S444 = 444;
constexpr int S467 = 467;

// R3BTwimMapped2Cal::Default Constructor --------------------------
R3BTwimMapped2Cal::R3BTwimMapped2Cal()
: R3BTwimMapped2Cal("R3BTwimMapped2Cal", 1)
Expand All @@ -60,13 +62,9 @@ R3BTwimMapped2Cal::R3BTwimMapped2Cal(const TString& name, Int_t iVerbose)
, fExpId(0)
, fOnline(kFALSE)
{
CalEParams.resize(fNumSec);
PosParams.resize(fNumSec);
for (Int_t s = 0; s < fNumSec; s++)
{
CalEParams[s] = NULL;
PosParams[s] = NULL;
}
CalEParams.resize(fNumSec, nullptr);
PosParams.resize(fNumSec, nullptr);
fTwimCal.resize(fNumSec);
}

void R3BTwimMapped2Cal::SetParContainers()
Expand Down Expand Up @@ -127,13 +125,18 @@ void R3BTwimMapped2Cal::SetParameter()
}

// Count the number of dead anodes or not used
for (Int_t s = 0; s < fNumSec; s++)
for (Int_t secId = 0; secId < fNumSec; secId++)
{
Int_t numdeadanodes = 0;
for (Int_t i = 0; i < fNumAnodes; i++)
if (CalEParams[s]->GetAt(fNumEParams * i + 1) == -1 || fCal_Par->GetInUse(s + 1, i + 1) == 0)
for (Int_t anodeId = 0; anodeId < fNumAnodes; anodeId++)
{
if (CalEParams[secId]->GetAt(fNumEParams * anodeId + 1) == -1 ||
fCal_Par->GetInUse(secId + 1, anodeId + 1) == 0)
{
numdeadanodes++;
R3BLOG(info, "Dead (or not used) anodes in section " << s + 1 << ": " << numdeadanodes);
}
}
R3BLOG(info, "Dead (or not used) anodes in section " << secId + 1 << ": " << numdeadanodes);
}
return;
}
Expand Down Expand Up @@ -164,9 +167,6 @@ InitStatus R3BTwimMapped2Cal::Init()
frm->Register("TwimCalData", "TWIM_Cal", fTwimCalDataCA, !fOnline);
Reset();

// Definition of a time stich object to correlate VFTX times
fTimeStitch = new R3BCoarseTimeStitch();

SetParameter();
return kSUCCESS;
}
Expand All @@ -179,8 +179,44 @@ InitStatus R3BTwimMapped2Cal::ReInit()
return kSUCCESS;
}

void R3BTwimMapped2Cal::CalAnode::Init()
{
fmult = 0;
fE.clear();
fDT.clear();
fE.resize(MAXNUMANODE, NAN);
fDT.resize(MAXNUMANODE, 0);
}

R3BTwimMapped2Cal::CalSection::CalSection()
{
fAnode.clear();
fAnode.resize(MAXNUMANODE);

fTref.clear();
fTref.resize(MAXNUMTREF);

fTrig.clear();
fTrig.resize(MAXNUMTRIG);
}
void R3BTwimMapped2Cal::CalSection::Init()
{
for (auto& val : fAnode)
{
val.Init();
}
for (auto& val : fTref)
{
val.Init();
}
for (auto& val : fTrig)
{
val.Init();
}
}

// ----- Public method Execution --------------------------------------------
void R3BTwimMapped2Cal::Exec(Option_t*)
void R3BTwimMapped2Cal::Exec(Option_t* /*option*/) // NOLINT(readability-function-cognitive-complexity)
{
// Reset entries in output arrays, local arrays
Reset();
Expand All @@ -190,101 +226,83 @@ void R3BTwimMapped2Cal::Exec(Option_t*)
if (nHits == 0)
return;

auto** mappedData = new R3BTwimMappedData*[nHits];
Int_t secId = 0;
Int_t anodeId = 0;
Double_t pedestal = 0.;
Double_t slope = 1.;

for (Int_t s = 0; s < fNumSec; s++)
for (Int_t i = 0; i < (fNumAnodes + fNumAnodesRef + fNumAnodesTrig); i++)
for (Int_t sec = 0; sec < fNumSec; sec++)
{
fTwimCal.at(sec).Init();
}
for (Int_t i = 0; i < nHits; i++)
{
auto* mappedData = dynamic_cast<R3BTwimMappedData*>(fTwimMappedDataCA->At(i));
auto secId = static_cast<Int_t>(mappedData->GetSecID()) - 1;
auto anodeId = static_cast<Int_t>(mappedData->GetAnodeID()) - 1;
if (anodeId < fNumAnodes && fCal_Par->GetInUse(secId + 1, anodeId + 1) == 1)
{
auto pedestal = CalEParams[secId]->GetAt(fNumEParams * anodeId);
auto slope = CalEParams[secId]->GetAt(fNumEParams * anodeId + 1);
fTwimCal.at(secId).GetAnode(anodeId).SetVal(pedestal + slope * mappedData->GetEnergy(),
mappedData->GetTime());
}
else if (anodeId >= fNumAnodes && anodeId < fNumAnodes + fNumAnodesRef)
{
mulanode[s][i] = 0;
for (Int_t j = 0; j < fMaxMult; j++)
auto TrefId = anodeId - fNumAnodes;
if ((fExpId == S444 || fExpId == S467) && secId != 0)
{
fE[s][j][i] = 0.;
fDT[s][j][i] = 0.;
secId = 0;
TrefId++;
}
fTwimCal.at(secId).GetTref(TrefId).SetVal(NAN, mappedData->GetTime());
}
}

for (Int_t i = 0; i < nHits; i++)
auto calc_dtime = [&](Int_t sec, Int_t mul_anode, Int_t id_anode, Int_t mul_tref)
{
mappedData[i] = dynamic_cast<R3BTwimMappedData*>(fTwimMappedDataCA->At(i));
secId = mappedData[i]->GetSecID() - 1;
anodeId = mappedData[i]->GetAnodeID() - 1;

if (anodeId < fNumAnodes && fCal_Par->GetInUse(secId + 1, anodeId + 1) == 1)
Int_t id_tref = 0;
if ((fExpId == S444 || fExpId == S467) && id_anode >= fNumAnodes / 2)
{
pedestal = CalEParams[secId]->GetAt(fNumEParams * anodeId);
slope = CalEParams[secId]->GetAt(fNumEParams * anodeId + 1);

fE[secId][mulanode[secId][anodeId]][anodeId] = pedestal + slope * mappedData[i]->GetEnergy();
fDT[secId][mulanode[secId][anodeId]][anodeId] = mappedData[i]->GetTime();
mulanode[secId][anodeId]++;
// s444 and s467: 2020
// anodes 1 to 16 : energy and time
// anode 17 and 18 : reference time
// will be cabled anode 19 and 20 : trigger time
id_tref = 1;
}
else if (anodeId >= fNumAnodes)
const auto tdiff = static_cast<Double_t>(fTwimCal.at(sec).GetAnode(id_anode).GetDT(mul_anode) -
fTwimCal.at(sec).GetTref(id_tref).GetDT(mul_tref));
//
auto i_param = fNumPosParams * id_anode;
auto par0 = PosParams[sec]->GetAt(i_param);
auto par1 = PosParams[sec]->GetAt(i_param + 1);
auto par2 = 0.0;
if (fNumPosParams > 2)
{
fDT[secId][mulanode[secId][anodeId]][anodeId] = mappedData[i]->GetTime(); // Ref. Time
mulanode[secId][anodeId]++;
par2 = PosParams[sec]->GetAt(i_param + 2);
}
}
return par0 + par1 * tdiff + par2 * pow(tdiff, 2);
};

// Fill data only if there is TREF signal
for (Int_t s = 0; s < fNumSec; s++)
if (mulanode[s][fNumAnodes] == 1)
for (Int_t sec = 0; sec < fNumSec; sec++)
{
if (fTwimCal.at(sec).GetTref().GetMult() != 1)
{
continue;
}
for (Int_t anodeId = 0; anodeId < fNumAnodes; anodeId++)
{
for (Int_t i = 0; i < fNumAnodes; i++)
for (Int_t mult = 0; mult < fTwimCal.at(sec).GetAnode(anodeId).GetMult(); mult++)
{
Int_t ii = fNumPosParams * i;
Float_t a0 = PosParams[s]->GetAt(ii);
Float_t a1 = PosParams[s]->GetAt(ii + 1);
Float_t a2 = 0.0;
if (fNumPosParams > 2)
a2 = PosParams[s]->GetAt(ii + 2);
for (Int_t j = 0; j < mulanode[s][fNumAnodes]; j++)
for (Int_t k = 0; k < mulanode[s][i]; k++)
{
if (fExpId == 444 || fExpId == 467)
{
// s444 and s467: 2020
// anodes 1 to 16 : energy and time
// anode 17 and 18 : reference time --> will be changed to 17 only when the full Twin-MUSIC
// will be cabled anode 19 and 20 : trigger time --> will be changed to 18 only when the
// full Twin-MUSIC will be cabled
if (i < 8) // Check if fNumAnodes+1 and +2 values
{
if (fE[s][k][i] > 0.)
AddCalData(s + 1,
i + 1,
a0 + a1 * fTimeStitch->GetTime(
fDT[s][k][i] - fDT[s][j][fNumAnodes + 1], "vftx", "vftx"),
fE[s][k][i]);
}
else
{
if (fE[s][k][i] > 0.)
AddCalData(s + 1,
i + 1,
a0 + a1 * fTimeStitch->GetTime(
fDT[s][k][i] - fDT[s][j][fNumAnodes + 2], "vftx", "vftx"),
fE[s][k][i]);
}
}
else if (fExpId == 455 && fE[s][k][i] > 0.)
{
auto dtime = a0 + a1 * (fDT[s][k][i] - fDT[s][j][fNumAnodes]) +
a2 * pow((fDT[s][k][i] - fDT[s][j][fNumAnodes]), 2);
if (dtime > fMinDT && dtime < fMaxDT)
{
AddCalData(s + 1, i + 1, dtime, fE[s][k][i]);
}
}
}
if (fTwimCal.at(sec).GetAnode(anodeId).GetE(mult) <= 0.)
{
continue;
}
const Int_t mult_Tref = 0; // mulanode != 1 hits are rejected already.
auto dtime = calc_dtime(sec, mult, anodeId, mult_Tref);
if (dtime > fMinDT && dtime < fMaxDT)
{
AddCalData(sec + 1, anodeId + 1, dtime, fTwimCal.at(sec).GetAnode(anodeId).GetE(mult));
}
}
}

if (mappedData)
delete[] mappedData;
}
return;
}

Expand Down
Loading

0 comments on commit b61b772

Please sign in to comment.