Skip to content
Snippets Groups Projects
Commit 51c16fe3 authored by Axel Puntke's avatar Axel Puntke
Browse files

Add TOF beta information

parent caf1c68e
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
analysis/detectors/trd/CbmTrdTrackTreeProducerTask.cxx
+ 108
− 28
View file @ 51c16fe3
......@@ -7,7 +7,7 @@
/// @author Sergey Gorbunov
/// @author Nora Bluhme
/// @date Feb 05 2025
/// @brief Task class for crreating trees containing tracks and TRD hits
/// @brief Task class for creating trees containing tracks and TRD hits
///
// Cbm Headers ----------------------
......@@ -42,6 +42,7 @@
#include "CbmTrdTrackingInterface.h"
#include "TCanvas.h"
#include "TGeoPhysicalNode.h"
#include "CbmKfTarget.h"
#include <cmath>
#include <iostream>
......@@ -93,6 +94,7 @@ CbmTrdTrackTreeProducerTask::CbmTrdTrackTreeProducerTask(const char* name, Int_t
fTree->Branch("trackX", &fCurRow.trackX, "trackX/F");
fTree->Branch("trackY", &fCurRow.trackY, "trackY/F");
fTree->Branch("trackZ", &fCurRow.trackZ, "trackZ/F");
fTree->Branch("tofBeta", &fCurRow.tofBeta, "tofBeta/D");
gFile = curFile;
gDirectory = curDirectory;
......@@ -123,32 +125,17 @@ InitStatus CbmTrdTrackTreeProducerTask::Init()
LOG(error) << "CbmTrdTrackTreeProducerTask::Init: Global track array not found!";
return kERROR;
}
fInputGlobalTrackMatches = static_cast<TClonesArray*>(ioman->GetObject("GlobalTrackMatch"));
if (!fInputGlobalTrackMatches) {
LOG(error) << "CbmTrdTrackTreeProducerTask::Init: Global track matches array not found!";
//return kERROR;
}
fInputStsTracks = static_cast<TClonesArray*>(ioman->GetObject("StsTrack"));
if (!fInputStsTracks) {
LOG(error) << "CbmTrdTrackTreeProducerTask::Init: Sts track array not found!";
fInputTofHits = static_cast<TClonesArray*>(ioman->GetObject("TofHit"));
if (!fInputTofHits) {
LOG(error) << "CbmTrdTrackTreeProducerTask::Init: Tof hit array not found!";
return kERROR;
}
fInputStsTrackMatches = static_cast<TClonesArray*>(ioman->GetObject("StsTrackMatch"));
if (!fInputStsTrackMatches) {
LOG(error) << "CbmTrdTrackTreeProducerTask::Init: Sts track matches array not found!";
//return kERROR;
}
// MC tracks
fInputMcTracks = static_cast<TClonesArray*>(ioman->GetObject("MCTrack"));
if (!fInputMcTracks) {
Warning("CbmTrdTrackTreeProducerTask::Init", "mc track array not found!");
//return kERROR;
fInputTrdHits = static_cast<TClonesArray*>(ioman->GetObject("TrdHit"));
if (!fInputTrdHits) {
LOG(error) << "CbmTrdTrackTreeProducerTask::Init: Trd hit array not found!";
return kERROR;
}
fTracks.clear();
......@@ -169,6 +156,16 @@ void CbmTrdTrackTreeProducerTask::Exec(Option_t* /*opt*/)
LOG(info) << "CbmTrdTrackTreeProducerTask: exec event N " << fNEvents;
fNEvents++;
// Obtain the correct timing information for each TOF hit (for TOF beta calculation)
TClonesArray* tcEvents = (TClonesArray*) FairRootManager::Instance()->GetObject("CbmEvent");
unsigned int nEvents = tcEvents->GetEntriesFast();
for (unsigned int i = 0; i < nEvents; ++i)
{
CbmEvent* event = static_cast<CbmEvent*>(tcEvents->At(i));
ShiftTofHitsToTzero(event);
}
// END TOF beta calculation initialization
ca::Framework* l1 = CbmL1::Instance()->fpAlgo;
const ca::Parameters<ca::fvec>& l1Par = l1->GetParameters();
......@@ -306,7 +303,6 @@ void CbmTrdTrackTreeProducerTask::Exec(Option_t* /*opt*/)
<< fTracks.size() << " tracks selected for alignment, " << active_tracks << " tracks active";
// fill the tree with the remaining hit-track pairs
TClonesArray* tcTrdHits = (TClonesArray*) FairRootManager::Instance()->GetObject("TrdHit");
for (const auto& t : fTracks) {
if (!t.fIsActive) continue;
......@@ -328,8 +324,20 @@ void CbmTrdTrackTreeProducerTask::Exec(Option_t* /*opt*/)
if (!node.fIsFitted) {
continue;
}
CbmTrdHit* trdHit = static_cast<CbmTrdHit*>(tcTrdHits->At(node.fHitIndex));
if (node.fHitIndex < 0 || node.fHitIndex >= fInputTrdHits->GetEntriesFast())
{
std::cout << "ERROR: node.fHitIndex is invalid! (idx " << node.fHitIndex << ", max " << fInputTrdHits->GetEntriesFast() << ")" << std::endl;
continue;
}
CbmTrdHit* trdHit = static_cast<CbmTrdHit*>(fInputTrdHits->At(node.fHitIndex));
if (!trdHit)
{
std::cout << "ERROR: trdHit is null!" << std::endl;
continue;
}
double x_residual = node.fMxy.X() - node.fParamDn.X(); // this should be the difference vector
double y_residual = node.fMxy.Y() - node.fParamDn.Y();
......@@ -356,6 +364,7 @@ void CbmTrdTrackTreeProducerTask::Exec(Option_t* /*opt*/)
fCurRow.trackX = node.fParamDn.X();
fCurRow.trackY = node.fParamDn.Y();
fCurRow.trackZ = node.fParamDn.Z();
fCurRow.tofBeta = CalculateTofBeta(tr);
//printf("fCurRow: hitpos (%f, %f, %f), module %d, charge %f, trackpos (%f, %f, %f)\n", fCurRow.hitX, fCurRow.hitY, fCurRow.hitZ, fCurRow.moduleID, fCurRow.charge, fCurRow.trackX, fCurRow.trackY, fCurRow.trackZ);
......@@ -388,4 +397,75 @@ void CbmTrdTrackTreeProducerTask::Finish()
std::cout << "-I- CbmTrdTrackTreeProducerTask: successfully wrote tree to " << fTreeFileName << std::endl;
}
double CbmTrdTrackTreeProducerTask::ShiftTofHitsToTzero(const CbmEvent* pEvent)
{
// Define t0 from the Bmon hits
double t0{std::numeric_limits<double>::signaling_NaN()};
int nTofHits = pEvent->GetNofData(ECbmDataType::kTofHit);
for (int iTofHitEvt{0}; iTofHitEvt < nTofHits; ++iTofHitEvt) {
int iTofHit = pEvent->GetIndex(ECbmDataType::kTofHit, iTofHitEvt);
const auto* pTofHit{static_cast<const CbmTofHit*>(fInputTofHits->At(iTofHit))};
//if (5 == CbmTofAddress::GetSmType(pTofHit->GetAddress())) { // selection by the supermodule type
if (pTofHit->GetZ() == 0) { // Take some small z-coordinate to identify BMon (TODO: provide a flag by some task)
t0 = pTofHit->GetTime();
}
}
// NOTE: t0 must be defined for each event, since the Bmon digis are used to seed the digi event builder. Basically,
// the tZero must be defined as a field of CbmEvent, e.g. in a separate FairTask, or directly
// in CbmAlgoBuildRawEvent
if (std::isnan(t0)) {
return t0;
}
// Shift TOF times to found t0:
for (int iTofHitEvt{0}; iTofHitEvt < nTofHits; ++iTofHitEvt) {
int iTofHit = pEvent->GetIndex(ECbmDataType::kTofHit, iTofHitEvt);
auto* pTofHit{static_cast<CbmTofHit*>(fInputTofHits->At(iTofHit))};
pTofHit->SetTime(pTofHit->GetTime() - t0 - fTzeroOffset);
}
return t0;
}
double CbmTrdTrackTreeProducerTask::CalculateTofBeta(CbmKfTrackFitter::Trajectory tr)
{
//Find Index of last TOF hit
int iLstTofHit = -1;
double tofHighestZ = 0;
for (unsigned int in = 0; in < tr.fNodes.size(); in++)
{
CbmKfTrackFitter::TrajectoryNode& node = tr.fNodes[in];
if (node.fHitSystemId == ECbmModuleId::kTof && node.fParamDn.Z() >= tofHighestZ)
{
tofHighestZ = node.fParamDn.Z();
iLstTofHit = node.fHitIndex;
}
}
if (iLstTofHit < 0 || iLstTofHit >= fInputTofHits->GetEntriesFast())
return 0;
const CbmTofHit* pTofHit{static_cast<CbmTofHit*>(fInputTofHits->At(iLstTofHit))};
if (!pTofHit)
return 0;
if (pTofHit->GetTime() < 0) { //< wrongly calibrated T0
return 0;
}
//calculate beta from last TOF hit
const auto* pTarget = cbm::kf::Target::Instance(); // CBM target info
double x{pTofHit->GetX() - pTarget->GetX()};
double y{pTofHit->GetY() - pTarget->GetY()};
double z{pTofHit->GetZ() - pTarget->GetZ()};
double x2{x * x};
double y2{y * y};
double z2{z * z};
double r2{x2 + y2 + z2};
double beta = std::sqrt(r2) / (pTofHit->GetTime() * kSpeedOfLight);
return beta;
}
ClassImp(CbmTrdTrackTreeProducerTask);
analysis/detectors/trd/CbmTrdTrackTreeProducerTask.h
+ 11
− 50
View file @ 51c16fe3
......@@ -7,7 +7,7 @@
/// @author Sergey Gorbunov
/// @author Nora Bluhme
/// @date Feb 05 2025
/// @brief Task class for crreating trees containing tracks and TRD hits
/// @brief Task class for creating trees containing tracks and TRD hits
///
#ifndef CbmTrdTrackTreeProducerTask_H
......@@ -28,16 +28,9 @@ class TClonesArray;
class TFile;
class TDirectory;
class TH1F;
class CbmEvent;
///
/// an example of alignment using BBA package
///
/// you need to switch to the double precision in /algo/ca/CaSimdVc.h
/// by uncommenting this line there:
///
/// typedef Vc::double_v fvec;
///
class CbmTrdTrackTreeProducerTask : public FairTask {
public:
// Constructors/Destructors ---------
......@@ -45,8 +38,6 @@ class CbmTrdTrackTreeProducerTask : public FairTask {
TString treeFileName = "CbmTrdTrackTree.root");
~CbmTrdTrackTreeProducerTask();
Int_t GetZtoNStation(Double_t getZ);
InitStatus Init();
void Exec(Option_t* opt);
void Finish();
......@@ -74,55 +65,25 @@ class CbmTrdTrackTreeProducerTask : public FairTask {
float charge; // Charge of the hit
int moduleID; // Module ID
float trackX, trackY, trackZ; // Fitted track position
};
struct AlignmentBody {
int fTrackingStation{-1};
};
struct Sensor {
ECbmModuleId fSystemId{0};
int fSensorId{-1};
int fTrackingStation{-1};
int fAlignmentBody{-1};
std::string fNodePath{""};
bool operator<(const Sensor& other) const
{
if (fSystemId < other.fSystemId) return true;
if (fSystemId > other.fSystemId) return false;
if (fTrackingStation < other.fTrackingStation) return true;
if (fTrackingStation > other.fTrackingStation) return false;
return (fSensorId < other.fSensorId);
};
bool operator==(const Sensor& other) const
{
return (fSystemId == other.fSystemId) && (fSensorId == other.fSensorId);
};
double tofBeta; // TOF beta information
};
private:
const CbmTrdTrackTreeProducerTask& operator=(const CbmTrdTrackTreeProducerTask&);
CbmTrdTrackTreeProducerTask(const CbmTrdTrackTreeProducerTask&);
double ShiftTofHitsToTzero(const CbmEvent* pEvent);
double CalculateTofBeta(CbmKfTrackFitter::Trajectory tr);
void WriteHistosCurFile(TObject* obj);
void ApplyAlignment(const std::vector<double>& par);
double CostFunction(const std::vector<double>& par);
void ApplyConstraints(std::vector<double>& par);
void ConstrainStation(std::vector<double>& par, int iSta, int ixyz);
static constexpr double kSpeedOfLight{29.9792458}; ///< Speed of light [cm/ns]
double fTzeroOffset{0.}; ///< Offset for T0
// input data arrays
TClonesArray* fInputGlobalTracks{nullptr};
TClonesArray* fInputStsTracks{nullptr};
TClonesArray* fInputMcTracks{nullptr}; // Mc info for debugging
TClonesArray* fInputGlobalTrackMatches{nullptr}; // Mc info for debugging
TClonesArray* fInputStsTrackMatches{nullptr}; // Mc info for debugging
TClonesArray* fInputTofHits{nullptr};
TClonesArray* fInputTrdHits{nullptr};
CbmKfTrackFitter fFitter;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment