diff --git a/analysis/common/analysis_tree_converter/CbmConverterManager.cxx b/analysis/common/analysis_tree_converter/CbmConverterManager.cxx index 256d5ffde3e65c1dca3553726013f06cc154abd0..0ed6c89ae03e70542b7d780900bdfa7a847a2804 100644 --- a/analysis/common/analysis_tree_converter/CbmConverterManager.cxx +++ b/analysis/common/analysis_tree_converter/CbmConverterManager.cxx @@ -127,7 +127,8 @@ void CbmConverterManager::FillDataHeader() auto* module = specDet_mod_pos.AddChannel(); module->SetPosition(x, y, frontFaceGlobal[2]); - LOG(info) << "Module " << nSpecDetModules << " : " << Form("(%.3f, %.3f)", x, y); + LOG(info) << "Module " << nSpecDetModules << " : " << Form("(%.3f, %.3f)", x, y) << " id " + << curNode->GetNumber(); } } diff --git a/analysis/common/analysis_tree_converter/CbmFsdHitsConverter.cxx b/analysis/common/analysis_tree_converter/CbmFsdHitsConverter.cxx index 02c7f5994dbecac89806349328de30a508e5fa63..7c1b9c7ce6b3d449abc7d4ab48ecaa3bd0c2f05e 100644 --- a/analysis/common/analysis_tree_converter/CbmFsdHitsConverter.cxx +++ b/analysis/common/analysis_tree_converter/CbmFsdHitsConverter.cxx @@ -87,6 +87,15 @@ void CbmFsdHitsConverter::Init() fsd_branch.AddField<float>("dtHitPoint", "Time difference between FSD hit and matched MC point [ps]"); fsd_branch.AddFields<float>({"dxHitPoint", "dyHitPoint", "dzHitPoint"}, "Component of a 3D distance between FSD hit and matched MC point [cm]"); + fsd_branch.AddFields<float>({"xPoint", "yPoint", "zPoint"}, "MC point distribution [cm]"); + fsd_branch.AddFields<float>({"pxPoint", "pyPoint", "pzPoint"}, "MC point momentum"); + fsd_branch.AddField<float>({"phiPoint"}, "Angle of the point"); + fsd_branch.AddField<float>({"lengthPoint"}, "Lenght of the point"); + fsd_branch.AddField<float>({"tPoint"}, "Time of the point"); + fsd_branch.AddField<float>({"elossPoint"}, "Energy loss of the point"); + fsd_branch.AddField<float>({"dist_middle_x"}, "Absolute value of the hit distance x from zero"); + fsd_branch.AddField<float>({"dist_middle_y"}, "Absolute value of the hit distance y from zero"); + i_edep_ = fsd_branch.GetFieldId("dEdx"); i_t_ = fsd_branch.GetFieldId("t"); @@ -100,6 +109,15 @@ void CbmFsdHitsConverter::Init() i_bestMatchedGTrack_ = fsd_branch.GetFieldId("bestMatchedGtrack2HitId"); i_multMC_ = fsd_branch.GetFieldId("multMCtracks"); i_topW_ = fsd_branch.GetFieldId("maxWeightMCtrack"); + i_xpoint_ = fsd_branch.GetFieldId("xPoint"); + i_pxpoint_ = fsd_branch.GetFieldId("pxPoint"); + i_phipoint_ = fsd_branch.GetFieldId("phiPoint"); + i_lengthpoint_ = fsd_branch.GetFieldId("lengthPoint"); + i_tpoint_ = fsd_branch.GetFieldId("tPoint"); + i_eloss_ = fsd_branch.GetFieldId("elossPoint"); + i_dist_middle_x_ = fsd_branch.GetFieldId("dist_middle_x"); + i_dist_middle_y_ = fsd_branch.GetFieldId("dist_middle_y"); + if (fsdgtrack_minChi2_ == 0) fsdgtrack_minChi2_ = -1.; if (fsdgtrack_maxChi2_ == 0) fsdgtrack_maxChi2_ = 10000.; @@ -175,6 +193,7 @@ void CbmFsdHitsConverter::ProcessData(CbmEvent* event) fsd_hits_->Reserve(n_fsd_hits); + // first loop over all FSD hits for (Int_t ifh = 0; ifh < n_fsd_hits; ifh++) { const auto fsdHI = event ? event->GetIndex(ECbmDataType::kFsdHit, ifh) : ifh; @@ -188,10 +207,22 @@ void CbmFsdHitsConverter::ProcessData(CbmEvent* event) const Float_t eLoss = fsdHit->GetEdep(); const Float_t time = fsdHit->GetTime(); + const Float_t dist_x = std::fabs(fsdHit->GetX()); + const Float_t dist_y = std::fabs(fsdHit->GetY()); + + hit.SetField(dist_x, i_dist_middle_x_); + hit.SetField(dist_y, i_dist_middle_y_); + hit.SetPosition(hitX, hitY, hitZ); hit.SetSignal(time); hit.SetField(eLoss, i_edep_); + Float_t phi_hit = atan2(fsdHit->GetY(), fsdHit->GetX()); + if (phi_hit < 0) { + phi_hit = phi_hit + 2 * 3.1415; + } + + Int_t multMC = 0; Float_t highestWeight = 0.; const auto* fsdHitMatch = dynamic_cast<CbmMatch*>(cbm_fsd_hitmatch_->At(fsdHI)); @@ -213,9 +244,9 @@ void CbmFsdHitsConverter::ProcessData(CbmEvent* event) } multMC++; // increase counter of MC multiplicity - // add matching between FSD hit and MC track - const auto* fsdPoint = dynamic_cast<FairMCPoint*>(cbm_fsd_points_new_->Get(pointLink)); + //const auto* fsdPoint = dynamic_cast<FairMCPoint*>(cbm_fsd_points_new_->Get(pointLink)); + const auto* fsdPoint = dynamic_cast<FairMCPoint*>(cbm_fsd_points_new_->Get(digiLink)); Int_t mc_track_id = fsdPoint->GetTrackID(); if (mc_track_id >= 0) { auto it = sim_tracks_map.find(mc_track_id); @@ -223,15 +254,32 @@ void CbmFsdHitsConverter::ProcessData(CbmEvent* event) fsd_hits_2_mc_tracks_->AddMatch(hit.GetId(), it->second); } } - // for the "matched" links store the difference between Hit and Point - if (fsdHitMatch->GetMatchedLink().GetIndex() == ilDigi - && fsdDigiMatch->GetMatchedLink().GetIndex() == ilPoint) { - hit.SetField(float(fsdPoint->GetTime() - time), i_dtHP_); - hit.SetField(float(fsdPoint->GetX() - hitX), i_dxHP_); - hit.SetField(float(fsdPoint->GetY() - hitY), i_dxHP_ + 1); - hit.SetField(float(fsdPoint->GetZ() - hitZ), i_dxHP_ + 2); + /////if (fsdHitMatch->GetMatchedLink().GetIndex() == ilDigi + ///// && fsdDigiMatch->GetMatchedLink().GetIndex() == ilPoint) { + //if(fsdPoint->GetLength() > 0.1){ + hit.SetField(float(fsdPoint->GetTime() - time), i_dtHP_); + hit.SetField(float(fsdPoint->GetX() - hitX), i_dxHP_); + hit.SetField(float(fsdPoint->GetY() - hitY), i_dxHP_ + 1); + hit.SetField(float(fsdPoint->GetZ() - hitZ), i_dxHP_ + 2); + hit.SetField(float(fsdPoint->GetEnergyLoss()), i_eloss_); + + hit.SetField(float(fsdPoint->GetX()), i_xpoint_); + hit.SetField(float(fsdPoint->GetY()), i_xpoint_ + 1); + hit.SetField(float(fsdPoint->GetZ()), i_xpoint_ + 2); + hit.SetField(float(fsdPoint->GetLength()), i_lengthpoint_); + hit.SetField(float(fsdPoint->GetPx()), i_pxpoint_); + hit.SetField(float(fsdPoint->GetPy()), i_pxpoint_ + 1); + hit.SetField(float(fsdPoint->GetPz()), i_pxpoint_ + 2); + + Float_t phi_point = atan2(fsdPoint->GetY(), fsdPoint->GetX()); + if (phi_point < 0) { + phi_point = phi_point + 2 * 3.1415; } + hit.SetField(float(phi_point), i_phipoint_); + hit.SetField(float(fsdPoint->GetTime()), i_tpoint_); + /////} + //} } // end of loop over links to points } @@ -259,19 +307,11 @@ void CbmFsdHitsConverter::ProcessData(CbmEvent* event) } } // end of loop over GlobalTracks - const Int_t storedIndex = bestMatchedIndex; - const Float_t storedChi2 = static_cast<Float_t>(bestChi2); - hit.SetField(storedChi2, i_chi2_); - hit.SetField(storedIndex, i_bestMatchedGTrack_); + hit.SetField(static_cast<Float_t>(bestChi2), i_chi2_); + hit.SetField(static_cast<Int_t>(bestMatchedIndex), i_bestMatchedGTrack_); if ((bestMatchedIndex >= 0) && (bestChi2 > fsdgtrack_minChi2_) && (bestChi2 < fsdgtrack_maxChi2_)) { - const Int_t trackIndex = event ? event->GetIndex(ECbmDataType::kGlobalTrack, bestMatchedIndex) : bestMatchedIndex; - if (trackIndex > cbm_global_tracks_->GetEntriesFast()) { - LOG(error) << "Trying to access element " << trackIndex << " which is out ot the global track array of size" - << cbm_global_tracks_->GetEntriesFast(); - continue; - } - const auto* globalTrack = static_cast<const CbmGlobalTrack*>(cbm_global_tracks_->At(trackIndex)); + const auto* globalTrack = static_cast<const CbmGlobalTrack*>(cbm_global_tracks_->At(bestMatchedIndex)); FairTrackParam param_last = *(globalTrack->GetParamLast()); FairTrackParam param_fsd = ExtrapolateGtrack(hitZ, param_last); @@ -309,3 +349,4 @@ CbmFsdHitsConverter::~CbmFsdHitsConverter() delete fsd_hits_; delete vtx_tracks_2_fsd_; }; +� diff --git a/analysis/common/analysis_tree_converter/CbmFsdHitsConverter.h b/analysis/common/analysis_tree_converter/CbmFsdHitsConverter.h index 2e3ca7abb36386b5ccd476085360c957c186c223..134838262fb4bc615069fc2305806429ba7d979d 100644 --- a/analysis/common/analysis_tree_converter/CbmFsdHitsConverter.h +++ b/analysis/common/analysis_tree_converter/CbmFsdHitsConverter.h @@ -84,6 +84,15 @@ class CbmFsdHitsConverter final : public CbmConverterTask { int i_dxHP_{AnalysisTree::UndefValueInt}; int i_dtHP_{AnalysisTree::UndefValueInt}; + int i_xpoint_{AnalysisTree::UndefValueInt}; + int i_pxpoint_{AnalysisTree::UndefValueInt}; + int i_phipoint_{AnalysisTree::UndefValueInt}; + int i_lengthpoint_{AnalysisTree::UndefValueInt}; + int i_tpoint_{AnalysisTree::UndefValueInt}; + int i_eloss_{AnalysisTree::UndefValueInt}; + int i_dist_middle_x_{AnalysisTree::UndefValueInt}; + int i_dist_middle_y_{AnalysisTree::UndefValueInt}; + Double_t fsdgtrack_minChi2_{0.}; Double_t fsdgtrack_maxChi2_{0.}; diff --git a/core/base/CbmMatchRecoToMC.cxx b/core/base/CbmMatchRecoToMC.cxx index 528651ebe9df7ff36c8630e0516d4404d24d6610..ab1604a2aefafc61e22141afb8e051c2abf90804 100644 --- a/core/base/CbmMatchRecoToMC.cxx +++ b/core/base/CbmMatchRecoToMC.cxx @@ -10,8 +10,10 @@ */ #include "CbmMatchRecoToMC.h" -#include "CbmCluster.h" // for CbmCluster -#include "CbmDigiManager.h" // for CbmDigiManager +#include "CbmCluster.h" // for CbmCluster +#include "CbmDigiManager.h" // for CbmDigiManager +#include "CbmFsdDigi.h" +#include "CbmFsdHit.h" // for CbmFsdHit #include "CbmHit.h" // for CbmHit, kMUCHPIXELHIT #include "CbmLink.h" // for CbmLink #include "CbmMCDataArray.h" // for CbmMCDataArray @@ -557,8 +559,15 @@ void CbmMatchRecoToMC::MatchHitsFsd(const TClonesArray* hits, TClonesArray* hitM for (Int_t iHit = 0; iHit < nofHits; iHit++) { CbmPixelHit* hit = static_cast<CbmPixelHit*>(hits->At(iHit)); CbmMatch* hitMatch = new ((*hitMatches)[iHit]) CbmMatch(); - const CbmMatch* digiMatch = fDigiManager->GetMatch(ECbmModuleId::kFsd, hit->GetRefId()); - hitMatch->AddLinks(*digiMatch); + Int_t nofDigis = fDigiManager->GetNofDigis(ECbmModuleId::kFsd); + for (Int_t iDigi = 0; iDigi < nofDigis; iDigi++) { + const CbmFsdDigi* digi = fDigiManager->Get<const CbmFsdDigi>(iDigi); + Int_t adresa_digi = digi->GetAddress(); + if (adresa_digi == hit->GetAddress()) { + const CbmMatch* digiMatch = fDigiManager->GetMatch(ECbmModuleId::kFsd, iDigi); + hitMatch->AddLinks(*digiMatch); + } + } } } diff --git a/core/detectors/fsd/CbmFsdGeoHandler.cxx b/core/detectors/fsd/CbmFsdGeoHandler.cxx index 1535e5e477a4edacf3c94ebafd1ca35ce3100c7d..8106efb7f284ce8e4c5b36eeb66d183443e773c4 100644 --- a/core/detectors/fsd/CbmFsdGeoHandler.cxx +++ b/core/detectors/fsd/CbmFsdGeoHandler.cxx @@ -51,7 +51,6 @@ void CbmFsdGeoHandler::InitMaps() TString nodeName(curNode->GetName()); if (nodeName.Contains(fUnitStr)) { currentUnitId = curNode->GetNumber(); - const TGeoMatrix* curUnitMatrix = geoIterator.GetCurrentMatrix(); unitToGlobalMatrix = new TGeoCombiTrans(*(curUnitMatrix)); @@ -64,6 +63,7 @@ void CbmFsdGeoHandler::InitMaps() if (nodeName.Contains(fModuleStr)) { currentModuleId = curNode->GetNumber(); + TGeoMatrix* moduleToUnitMatrix = curNode->GetMatrix(); TVector3 localModuleCoord(moduleToUnitMatrix->GetTranslation()); TVector3 globalModuleCoord; diff --git a/reco/detectors/fsd/CbmFsdHitProducer.cxx b/reco/detectors/fsd/CbmFsdHitProducer.cxx index 8b50cef22308cecdec580115904a10b20a4b9122..9b9a6f5a53bc0f9a6e6dcb7c87266e3dd636f349 100644 --- a/reco/detectors/fsd/CbmFsdHitProducer.cxx +++ b/reco/detectors/fsd/CbmFsdHitProducer.cxx @@ -224,7 +224,6 @@ pair<Int_t, Int_t> CbmFsdHitProducer::ProcessData(CbmEvent* event) Double_t sumEdep = std::get<0>(entry.second); Double_t fastestDigiTime = std::get<1>(entry.second); Int_t fastestDigiIndex = std::get<2>(entry.second); - CbmFsdModuleSpecs* fsdModSpecs = CbmFsdGeoHandler::GetInstance().GetModuleSpecsById( CbmFsdAddress::GetElementId(address, static_cast<int32_t>(CbmFsdAddress::Level::Module))); if (!fsdModSpecs) { diff --git a/sim/transport/base/CbmEventGenerator.cxx b/sim/transport/base/CbmEventGenerator.cxx index de4ac566502e8b9ff25f53103e0587036dd0ad58..aeb944a7d1a95325656ca6fc1dce1ba62aefe61a 100644 --- a/sim/transport/base/CbmEventGenerator.cxx +++ b/sim/transport/base/CbmEventGenerator.cxx @@ -100,12 +100,11 @@ Bool_t CbmEventGenerator::GenerateEvent(FairGenericStack* stack) fEvent->SetRotY(fBeamAngleY); Double_t phiGen = fEvent->GetRotZ(); // from concrete generators fEvent->SetRotZ(phiGen + fPhi); // total event plane angle - // Event info to screen LOG(info) << GetName() << ": Event ID " << fEvent->GetEventID() << ", tracks " << fNTracks << ", vertex (" << fVertex.X() << ", " << fVertex.Y() << ", " << fVertex.Z() << ") cm"; LOG(info) << GetName() << ": Beam angle (" << fBeamAngleX << ", " << fBeamAngleY << ") rad, event plane angle " - << fPhi << " rad (total " << fEvent->GetRotZ() << " rad)"; + << fPhi << " rad (total " << fEvent->GetRotZ() << " rad) change"; // Update global track counter (who knows what it's good for ...) fTotPrim += fNTracks; diff --git a/sim/transport/generators/CbmUnigenGenerator.cxx b/sim/transport/generators/CbmUnigenGenerator.cxx index 6be11ca78e347220ab9fa1909466195326233cc9..e21a942a02e8dfe00456bf6c8b235d7c1a1cef56 100644 --- a/sim/transport/generators/CbmUnigenGenerator.cxx +++ b/sim/transport/generators/CbmUnigenGenerator.cxx @@ -263,8 +263,8 @@ Bool_t CbmUnigenGenerator::ReadEvent(FairPrimaryGenerator* primGen) // Lorentz transformation into lab (target) frame Double_t pz = fGammaCM * (particle->Pz() - fBetaCM * particle->E()); - TVector3 momentum(particle->Px(), particle->Py(), pz); + TVector3 momentum(particle->Px(), particle->Py(), pz); // Apply event rotation if needed if (phi2 != 0.) momentum.RotateZ(phi2);