diff --git a/CODEOWNERS b/CODEOWNERS index 3c31287f6727b..6b51b6cfc630e 100644 --- a/CODEOWNERS +++ b/CODEOWNERS @@ -73,6 +73,7 @@ /Detectors/ZDC @coppedis /Detectors/CTF @shahor02 /Detectors/Raw @shahor02 +/Detectors/StrangenessTracking @mconcas @mpuccio @fmazzasc /EventVisualisation @jmyrcha #/EventVisualisation/Base diff --git a/DataFormats/Reconstruction/src/ReconstructionDataFormatsLinkDef.h b/DataFormats/Reconstruction/src/ReconstructionDataFormatsLinkDef.h index ad235797d79ff..e26510b5723d5 100644 --- a/DataFormats/Reconstruction/src/ReconstructionDataFormatsLinkDef.h +++ b/DataFormats/Reconstruction/src/ReconstructionDataFormatsLinkDef.h @@ -24,6 +24,8 @@ #pragma link C++ class o2::track::TrackParCovD + ; #pragma link C++ class o2::track::TrackParCov + ; #pragma link C++ class o2::track::TrackParametrizationWithError < float> + ; +#pragma link C++ class std::vector < o2::track::TrackParametrizationWithError < float>> + ; + #pragma link C++ class o2::track::TrackParametrizationWithError < double> + ; #pragma link C++ class o2::track::TrackParFwd + ; #pragma link C++ class o2::track::PID + ; diff --git a/Detectors/CMakeLists.txt b/Detectors/CMakeLists.txt index 255a155136e0e..473f6dffa5cd6 100644 --- a/Detectors/CMakeLists.txt +++ b/Detectors/CMakeLists.txt @@ -36,7 +36,7 @@ add_subdirectory(FOCAL) add_subdirectory(GlobalTracking) add_subdirectory(GlobalTrackingWorkflow) add_subdirectory(Vertexing) - +add_subdirectory(StrangenessTracking) if(BUILD_ANALYSIS) add_subdirectory(AOD) endif() diff --git a/Detectors/StrangenessTracking/CMakeLists.txt b/Detectors/StrangenessTracking/CMakeLists.txt new file mode 100644 index 0000000000000..723a7693e0c6f --- /dev/null +++ b/Detectors/StrangenessTracking/CMakeLists.txt @@ -0,0 +1,14 @@ +# Copyright 2019-2020 CERN and copyright holders of ALICE O2. +# See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +# All rights not expressly granted are reserved. +# +# This software is distributed under the terms of the GNU General Public +# License v3 (GPL Version 3), copied verbatim in the file "COPYING". +# +# In applying this license CERN 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. + +add_subdirectory(workflow) +add_subdirectory(tracking) +add_subdirectory(macros) \ No newline at end of file diff --git a/Detectors/StrangenessTracking/macros/CMakeLists.txt b/Detectors/StrangenessTracking/macros/CMakeLists.txt new file mode 100644 index 0000000000000..5de6df71aa8fb --- /dev/null +++ b/Detectors/StrangenessTracking/macros/CMakeLists.txt @@ -0,0 +1,10 @@ +o2_add_test_root_macro(XiTrackingStudy.C + PUBLIC_LINK_LIBRARIES O2::MathUtils + O2::ITSBase + O2::ITSMFTReconstruction + O2::ITSMFTSimulation + O2::DataFormatsITSMFT + O2::DataFormatsITS + O2::SimulationDataFormat + O2::StrangenessTracking + LABELS strangeness-tracking) \ No newline at end of file diff --git a/Detectors/StrangenessTracking/macros/XiTrackingStudy.C b/Detectors/StrangenessTracking/macros/XiTrackingStudy.C new file mode 100644 index 0000000000000..6e666df11177c --- /dev/null +++ b/Detectors/StrangenessTracking/macros/XiTrackingStudy.C @@ -0,0 +1,644 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN 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. +/// \file XiTrackingStudy.C +/// \brief Simple macro to check Xi strange tracks + +#if !defined(__CLING__) || defined(__ROOTCLING__) +#include "ReconstructionDataFormats/PID.h" +#include "ReconstructionDataFormats/V0.h" +#include "SimulationDataFormat/MCTruthContainer.h" +#include "SimulationDataFormat/MCCompLabel.h" +#include "SimulationDataFormat/MCTrack.h" +#include "ITSMFTSimulation/Hit.h" + +#include "DataFormatsITSMFT/TopologyDictionary.h" +#include "DetectorsCommonDataFormats/DetectorNameConf.h" +#include "ITSBase/GeometryTGeo.h" +#include "DataFormatsITS/TrackITS.h" +#include "DataFormatsITSMFT/CompCluster.h" +#include "DataFormatsITSMFT/ROFRecord.h" +#include "ITStracking/IOUtils.h" + +#include +#include +#include "TCanvas.h" +#include "TFile.h" +#include "TH1F.h" +#include "TH2D.h" +#include "TSystemDirectory.h" +#include "TMath.h" +#include "TString.h" +#include "TTree.h" +#include "TLegend.h" +#include "CommonDataFormat/RangeReference.h" +#include "DetectorsVertexing/DCAFitterN.h" +#include "StrangenessTracking/StrangenessTracker.h" + +#endif + +using GIndex = o2::dataformats::VtxTrackIndex; +using V0 = o2::dataformats::V0; +using Cascade = o2::dataformats::Cascade; + +using MCTrack = o2::MCTrack; +using VBracket = o2::math_utils::Bracket; +using namespace o2::itsmft; +using CompClusterExt = o2::itsmft::CompClusterExt; +using ITSCluster = o2::BaseCluster; +using Vec3 = ROOT::Math::SVector; +using StrangeTrack = o2::strangeness_tracking::StrangeTrack; + +const int motherPDG = 3312; +const int v0PDG = 3122; +const int bachPDG = 211; +const int firstV0dauPDG = 2212; +const int secondV0dauPDG = 211; + +const float XiMass = 1.32171; +const int kDauPdgs[3] = {2212, 211, 211}; +const float kDauMasses[3] = {0.938272, 0.13957, 0.13957}; +const float kFirstDauMasses[2] = {1.115683, 0.13957}; + +std::array matchCascToMC(const std::vector>& mcTracksMatrix, std::map*>& map, std::vector* v0vec, Cascade& casc, bool& isV0reco); +std::array checkV0mother(const std::vector>& mcTracksMatrix, std::map*>& map, V0& v0); +double calcDecLength(std::vector* MCTracks, const MCTrack& motherTrack); +std::array matchITStracktoMC(const std::vector>& mcTracksMatrix, o2::MCCompLabel ITSlabel); +double calcMass(std::vector tracks); +std::array matchCompLabelToMC(const std::vector>& mcTracksMatrix, o2::MCCompLabel compLabel); + +void XiTrackingStudy(std::string path) +{ + TSystemDirectory dir("MyDir", path.data()); + auto files = dir.GetListOfFiles(); + std::vector dirs; + std::vector kine_files; + + for (auto fileObj : *files) { + std::string file = ((TSystemFile*)fileObj)->GetName(); + if (file.substr(0, 2) == "tf") { + dirs.push_back(path + "/" + file); + auto innerdir = (TSystemDirectory*)fileObj; + auto innerfiles = innerdir->GetListOfFiles(); + for (auto innerfileObj : *innerfiles) { + TString innerfile = ((TSystemFile*)innerfileObj)->GetName(); + if (innerfile.EndsWith("Kine.root") && innerfile.Contains("sgn")) { + kine_files.push_back(innerfile); + } + } + } + } + + TH2D* hHyperhisto = new TH2D("histo hyperV0", ";#it{p}_{T} (GeV/#it{c}); Radius^2 (cm) ; Counts", 20, 1, 10, 30, 1, 900); + TH1D* hResV0histo = new TH1D("pT resolution before hypertracking", ";(#it{p}_{T}^{gen} - #it{p}_{T}^{rec})/#it{p}_{T}^{gen}; Counts", 20, -0.2, 0.2); + TH1D* hResHyperhisto = new TH1D("pT resolution after hypertracking", ";(#it{p}_{T}^{gen} - #it{p}_{T}^{rec})/#it{p}_{T}^{gen}; Counts", 20, -0.2, 0.2); + TH1D* hResCascR2 = new TH1D("R2 resolution before tracking", ";(R2^{gen} - R2^{rec})/R2^{gen}; Counts", 200, -0.1, 0.1); + TH1D* hResCascTrackedR2 = new TH1D("R2 resolution after tracking", ";(R2^{gen} - R2^{rec})/R2^{gen}; Counts", 200, -0.1, 0.1); + + TH1D* hCascCounter = new TH1D("Casc counter", ";Casc counter; Counts", 1, 0.5, 1.5); + TH1D* hGenXiCounter = new TH1D("Gen Casc counter", ";Casc counter; Counts", 1, 0.5, 1.5); + + TH1D* hFakeAssocCounter = new TH1D("Fake assoc counter", ";Fake assoc counter; Counts", 1, 0.5, 1.5); + + TH1D* hRecXiCounter = new TH1D("Rec Xi counter", "; ; Counts", 1, 0.5, 1.5); + TH1D* hCascMomsInV0s = new TH1D("Rec V0s from Xi counter", "; ; Counts", 1, 0.5, 1.5); + TH1D* hFindableBachfromXiCounter = new TH1D("Rec bachelors from Xi counter", "; ; Counts", 1, 0.5, 1.5); + + TH1D* hRecCascInvMass = new TH1D("Rec Casc InvMass", "; M (GeV/c^{2}); Counts", 150, 1.2, 1.5); + TH1D* hStrangeTrackInvMass = new TH1D("Strange track inv mass", ";M (GeV/c^{2}); Counts", 150, 1.2, 1.5); + TH1D* hXiStats = new TH1D("cascade_stats", "; ; Counts", 3, 0.5, 3.5); + + TH1D* hGenXiRadius = new TH1D("gen_casc_r", "; Radius (cm); Counts", 40, 0., 40.); + TH1D* hRecXiRad = new TH1D("rec_casc_r", ";Radius_{trackable} (cm); Counts", 40, 0, 40); + TH1D* hTrackedXiRad = new TH1D("r_dist_aft_track", ";R2^{rec} (cm); Counts", 40, 0, 40); + + TH1D* hGenXiMom = new TH1D("gen_casc_pt", "; #it{p}_{T}^{gen} (GeV/#it{c}); Counts", 40, 1, 6); + TH1D* hRecXiMom = new TH1D("rec_casc_pt", "; #it{p}_{T}^{gen} (GeV/#it{c}); Counts", 40, 1, 6); + TH1D* hTrackedXiMom = new TH1D("pT_dist_aft_track", ";#it{p}_{T}^{rec} (GeV/#it{c}); Counts", 40, 1, 6); + + // Geometry + o2::base::GeometryManager::loadGeometry(dirs[0] + "/o2sim_geometry.root"); + auto gman = o2::its::GeometryTGeo::Instance(); + gman->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::T2L, o2::math_utils::TransformType::L2G)); + + int counter = 0; + for (unsigned int i = 0; i < dirs.size(); i++) { + auto& dir = dirs[i]; + auto& kine_file = kine_files[i]; + LOG(info) << "Processing " << dir; + LOG(info) << "kine file " << kine_file; + // Files + auto fMCTracks = TFile::Open((TString(dir + "/") + kine_file)); + auto fStrangeTracks = TFile::Open((dir + "/o2_strange_tracks.root").data()); + auto fSecondaries = TFile::Open((dir + "/o2_secondary_vertex.root").data()); + auto fITSTPC = TFile::Open((dir + "/o2match_itstpc.root").data()); + auto fTPCTOF = TFile::Open((dir + "/o2match_tof_tpc.root").data()); + auto fTPCTRD = TFile::Open((dir + "/trdmatches_tpc.root").data()); + auto fITSTPCTOF = TFile::Open((dir + "/o2match_tof_itstpc.root").data()); + auto fITS = TFile::Open((dir + "/o2trac_its.root").data()); + auto fClusITS = TFile::Open((dir + "/o2clus_its.root").data()); + auto fTPC = TFile::Open((dir + "/tpctracks.root").data()); + + // Trees + auto treeMCTracks = (TTree*)fMCTracks->Get("o2sim"); + auto treeStrangeTracks = (TTree*)fStrangeTracks->Get("o2sim"); + auto treeSecondaries = (TTree*)fSecondaries->Get("o2sim"); + auto treeITSTPC = (TTree*)fITSTPC->Get("matchTPCITS"); + auto treeITSTPCTOF = (TTree*)fITSTPCTOF->Get("matchTOF"); + auto treeTPCTOF = (TTree*)fTPCTOF->Get("matchTOF"); + auto treeTPCTRD = (TTree*)fTPCTRD->Get("tracksTRD"); + + auto treeITS = (TTree*)fITS->Get("o2sim"); + auto treeITSclus = (TTree*)fClusITS->Get("o2sim"); + auto treeTPC = (TTree*)fTPC->Get("tpcrec"); + + // MC Tracks + std::vector* MCtracks = nullptr; + std::vector* ITSHits = nullptr; + + // Hypertracks + std::vector* strangeTrackVec = nullptr; + std::vector>* hypertrackVec = nullptr; + std::vector* nAttachments = nullptr; + + // Secondary Vertices + std::vector* cascVec = nullptr; + std::vector* v0Vec = nullptr; + + // ITS tracks + std::vector* ITStracks = nullptr; + + // Labels + std::vector* labITSvec = nullptr; + std::vector* labTPCvec = nullptr; + std::vector* labITSTPCvec = nullptr; + std::vector* labITSTPCTOFvec = nullptr; + std::vector* labTPCTOFvec = nullptr; + std::vector* labTPCTRDvec = nullptr; + + // Clusters + std::vector* ITSclus = nullptr; + o2::dataformats::MCTruthContainer* clusLabArr = nullptr; + std::vector* ITSTrackClusIdx = nullptr; + std::vector* ITSpatt = nullptr; + + // Setting branches + treeStrangeTracks->SetBranchAddress("StrangeTracks", &strangeTrackVec); + treeStrangeTracks->SetBranchAddress("ClusUpdates", &nAttachments); + + treeSecondaries->SetBranchAddress("Cascades", &cascVec); + treeSecondaries->SetBranchAddress("V0s", &v0Vec); + + treeMCTracks->SetBranchAddress("MCTrack", &MCtracks); + + treeITS->SetBranchAddress("ITSTrackMCTruth", &labITSvec); + treeITS->SetBranchAddress("ITSTrack", &ITStracks); + treeTPC->SetBranchAddress("TPCTracksMCTruth", &labTPCvec); + treeITSTPC->SetBranchAddress("MatchMCTruth", &labITSTPCvec); + treeTPCTOF->SetBranchAddress("MatchTOFMCTruth", &labTPCTOFvec); + treeTPCTRD->SetBranchAddress("labels", &labTPCTRDvec); + + treeITSTPCTOF->SetBranchAddress("MatchTOFMCTruth", &labITSTPCTOFvec); + + treeITS->SetBranchAddress("ITSTrackClusIdx", &ITSTrackClusIdx); + treeITSclus->SetBranchAddress("ITSClusterComp", &ITSclus); + treeITSclus->SetBranchAddress("ITSClusterMCTruth", &clusLabArr); + + // define detector map + std::map*> map{{"ITS", labITSvec}, {"TPC", labTPCvec}, {"ITS-TPC", labITSTPCvec}, {"TPC-TOF", labTPCTOFvec}, {"TPC-TRD", labTPCTRDvec}, {"ITS-TPC-TOF", labITSTPCTOFvec}}; + + // fill MC matrix + std::vector> mcTracksMatrix; + auto nev = treeMCTracks->GetEntriesFast(); + mcTracksMatrix.resize(nev); + for (int n = 0; n < nev; n++) { // loop over MC events + treeMCTracks->GetEvent(n); + mcTracksMatrix[n].resize(MCtracks->size()); + for (unsigned int mcI{0}; mcI < MCtracks->size(); ++mcI) { + mcTracksMatrix[n][mcI] = MCtracks->at(mcI); + if (abs(MCtracks->at(mcI).GetPdgCode()) == motherPDG) { + auto& motherTrack = mcTracksMatrix[n][mcI]; + hGenXiCounter->Fill(1); + hGenXiMom->Fill(motherTrack.GetPt()); + hGenXiRadius->Fill(calcDecLength(MCtracks, motherTrack)); + } + } + } + + // Starting matching Cascades and ITS tracks + int counterV0 = 0; + for (int frame = 0; frame < treeITS->GetEntriesFast(); frame++) { + if (!treeITS->GetEvent(frame) || !treeITS->GetEvent(frame) || !treeSecondaries->GetEvent(frame) || !treeITSTPC->GetEvent(frame) || !treeTPC->GetEvent(frame) || + !treeITSTPCTOF->GetEvent(frame) || !treeTPCTOF->GetEvent(frame) || !treeITSclus->GetEvent(frame) || !treeTPCTRD->GetEvent(frame) || !treeStrangeTracks->GetEvent(frame)) + continue; + + for (unsigned int iCascVec = 0; iCascVec < cascVec->size(); iCascVec++) { + hCascCounter->Fill(1); + auto& casc = cascVec->at(iCascVec); + + hRecCascInvMass->Fill(sqrt(casc.calcMass2())); + + bool isV0reco = false; + auto cascMCref = matchCascToMC(mcTracksMatrix, map, v0Vec, casc, isV0reco); + + if (cascMCref[0] == -1 || cascMCref[1] == -1) + continue; + + hRecXiCounter->Fill(1); + hRecXiRad->Fill(sqrt(casc.calcR2())); + hRecXiMom->Fill(casc.getPt()); + auto& mcCasc = mcTracksMatrix[cascMCref[0]][cascMCref[1]]; + + // Matching ITS tracks to MC tracks and V0 + std::array ITSref = {-1, 1}; + o2::its::TrackITS ITStrack; + std::array, 7> clsRef; + + int iTrack = -1; + bool isMatched = false; + + for (unsigned int iITStrack = 0; iITStrack < ITStracks->size(); iITStrack++) { + auto& labITS = (*labITSvec)[iITStrack]; + auto& trackIdx = (*ITSTrackClusIdx)[iITStrack]; + + ITSref = matchITStracktoMC(mcTracksMatrix, labITS); + ITStrack = (*ITStracks)[iITStrack]; + + if (ITSref[0] == cascMCref[0] && ITSref[1] == cascMCref[1]) { + LOG(info) << "++++++++++++++++"; + LOG(info) << "Cascade + ITS track found! "; + LOG(info) << "CascV0: " << casc.getV0ID() << ", Bach ID: " << casc.getBachelorID() << ", ITS track ref: " << iITStrack; + LOG(info) << "Casc Momentum: " << casc.getPt() << ", ITS track momentum: " << ITStrack.getPt(); + LOG(info) << "Casc Cov Y: " << casc.getSigmaY2() << ", Cov Z: " << casc.getSigmaZ2(); + + auto& MCTrack = mcTracksMatrix[ITSref[0]][ITSref[1]]; + auto& MCTracks = mcTracksMatrix[ITSref[0]]; + auto MCrad = calcDecLength(&MCTracks, MCTrack); + LOG(info) << "Casc rad: " << casc.calcR2() << ", MC radius: " << MCrad; + hXiStats->Fill(1); + + auto firstClus = ITStrack.getFirstClusterEntry(); + auto ncl = ITStrack.getNumberOfClusters(); + + // check whether the corresponding strange track is found + bool isStrangeTrackFound = false; + for (unsigned int iStTr = 0; iStTr < strangeTrackVec->size(); iStTr++) { + + auto& strangeTrack = strangeTrackVec->at(iStTr); + if (!(strangeTrack.mPartType == o2::strangeness_tracking::kCascade)) { + continue; + } + + auto& itsRef = strangeTrack.mITSRef; + auto& cascRef = strangeTrack.mDecayRef; + if (itsRef == int(iITStrack) && cascRef == int(iCascVec)) { + isStrangeTrackFound = true; + break; + } + } + + LOG(info) << "is Strange Track Found? " << isStrangeTrackFound; + + for (int icl = 0; icl < ncl; icl++) { + auto& labCls = (clusLabArr->getLabels(ITSTrackClusIdx->at(firstClus + icl)))[0]; + auto& clus = (*ITSclus)[(*ITSTrackClusIdx)[firstClus + icl]]; + auto layer = gman->getLayer(clus.getSensorID()); + clsRef[layer] = matchCompLabelToMC(mcTracksMatrix, labCls); + if (clsRef[layer][0] > -1 && clsRef[layer][1] > -1) { + auto& mcTrack = mcTracksMatrix[clsRef[layer][0]][clsRef[layer][1]]; + auto pdg = mcTrack.GetPdgCode(); + auto pdgMom1 = -1, pdgMom2 = -1; + if (!mcTrack.isPrimary()) { + auto& mcTrackMom1 = mcTracksMatrix[clsRef[layer][0]][mcTrack.getMotherTrackId()]; + pdgMom1 = mcTrackMom1.GetPdgCode(); + if (!mcTrackMom1.isPrimary()) { + pdgMom2 = mcTracksMatrix[clsRef[layer][0]][mcTrackMom1.getMotherTrackId()].GetPdgCode(); + } + } + LOG(info) << "Layer: " << layer << ", pdg " << pdg << ", pdgMom1 " << pdgMom1 << ", pdgMom2 " << pdgMom2 << ", refs: " << clsRef[layer][0] << ", " << clsRef[layer][1]; + } else + LOG(info) << "Layer: " << layer << ", No valid cluster ref" + << ", isNoise: " << labCls.isNoise() << ", isFake: " << labCls.isFake(); + } + } + } + } + + LOG(info) << "+++++++++++++++++++++++++++++++++++++++++++++"; + + for (unsigned int iHyperVec = 0; iHyperVec < strangeTrackVec->size(); iHyperVec++) { + + auto& strangeTrack = strangeTrackVec->at(iHyperVec); + if (!(strangeTrack.mPartType == o2::strangeness_tracking::kCascade)) { + continue; + } + + auto& hyperChi2 = strangeTrack.mMatchChi2; + auto& clusAttachments = nAttachments->at(iHyperVec); + auto& ITStrack = ITStracks->at(strangeTrack.mITSRef); + auto& ITStrackLab = labITSvec->at(strangeTrack.mITSRef); + + auto clusAttArr = clusAttachments.arr; + auto& casc = cascVec->at(strangeTrack.mDecayRef); + auto& v0Casc = casc.getV0Track(); + auto& bach = casc.getBachelorTrack(); + + auto ITStrackRef = matchCompLabelToMC(mcTracksMatrix, ITStrackLab); + bool isV0reco = false; + auto cascMCref = matchCascToMC(mcTracksMatrix, map, v0Vec, casc, isV0reco); + + std::vector dauTracks = {v0Casc, bach}; + + hStrangeTrackInvMass->Fill(calcMass(dauTracks)); + + if (cascMCref[0] == -1 || cascMCref[1] == -1 || ITStrackRef[0] == -1 || ITStrackRef[1] == -1) + continue; + + if (ITStrackRef[0] == cascMCref[0] && ITStrackRef[1] == cascMCref[1]) { + hXiStats->Fill(2); + hTrackedXiRad->Fill(sqrt(casc.calcR2())); + hTrackedXiMom->Fill(casc.getPt()); + + auto& mcTracks = mcTracksMatrix[cascMCref[0]]; + auto& mcTrack = mcTracksMatrix[cascMCref[0]][cascMCref[1]]; + hResCascR2->Fill((sqrt(casc.calcR2()) - calcDecLength(&mcTracks, mcTrack)) / calcDecLength(&mcTracks, mcTrack)); + + auto trackedR = sqrt(strangeTrack.decayVtx[0] * strangeTrack.decayVtx[0] + strangeTrack.decayVtx[1] * strangeTrack.decayVtx[1]); + hResCascTrackedR2->Fill((trackedR - calcDecLength(&mcTracks, mcTrack)) / calcDecLength(&mcTracks, mcTrack)); + } else { + hFakeAssocCounter->Fill(1); + } + } + } + } + + auto outFile = TFile("cascade_study.root", "recreate"); + + hResV0histo->Write(); + hHyperhisto->Write(); + hResHyperhisto->Write(); + hResCascR2->Write(); + hResCascTrackedR2->Write(); + + hGenXiCounter->Write(); + hCascCounter->Write(); + hXiStats->Write(); + hRecXiCounter->Write(); + hFakeAssocCounter->Write(); + hCascMomsInV0s->Write(); + hRecCascInvMass->Write(); + hStrangeTrackInvMass->Write(); + + TH1D* hMomEffFound = (TH1D*)hTrackedXiMom->Clone("efficiency_found_mom"); + hMomEffFound->Divide(hRecXiMom); + auto cv1 = TCanvas("found_efficiency_pt", "", 1000, 1000); + hMomEffFound->GetYaxis()->SetTitle("Tracked/Found"); + hMomEffFound->SetLineColor(kBlue); + hMomEffFound->Draw(); + cv1.Write(); + + TH1D* hRadEffFound = (TH1D*)hTrackedXiRad->Clone("efficiency_found_rad"); + hRadEffFound->Divide(hRecXiRad); + auto cv2 = TCanvas("found_efficiency_rad", "", 1000, 1000); + hRadEffFound->GetXaxis()->SetTitle("Radius (cm)"); + hRadEffFound->GetYaxis()->SetTitle("Tracked/Found"); + hRadEffFound->SetLineColor(kBlue); + hRadEffFound->Draw(); + cv2.Write(); + + TH1D* hMomEffGen = (TH1D*)hTrackedXiMom->Clone("efficiency_gen_mom"); + hMomEffGen->Divide(hGenXiMom); + auto cv3 = TCanvas("gen_efficiency_pt", "", 1000, 1000); + hMomEffGen->GetYaxis()->SetTitle("Efficiency"); + hMomEffGen->SetLineColor(kBlue); + hMomEffGen->Draw(); + cv3.Write(); + + TH1D* hRadEffGen = (TH1D*)hTrackedXiRad->Clone("efficiency_gen_rad"); + hRadEffGen->Divide(hGenXiRadius); + auto cv4 = TCanvas("gen_efficiency_rad", "", 1000, 1000); + hRadEffGen->GetXaxis()->SetTitle("Radius (cm)"); + hRadEffGen->GetYaxis()->SetTitle("Efficiency"); + hRadEffGen->SetLineColor(kBlue); + hRadEffGen->Draw(); + cv4.Write(); + + auto cv5 = TCanvas("inv_mass_res_study", "", 1000, 1000); + hStrangeTrackInvMass->GetXaxis()->SetTitle("M(GeV/#it{c}^{2})"); + hStrangeTrackInvMass->GetYaxis()->SetTitle("Normalised Counts"); + hStrangeTrackInvMass->SetLineColor(kRed); + hStrangeTrackInvMass->SetLineWidth(2); + hRecCascInvMass->SetLineWidth(2); + hStrangeTrackInvMass->DrawNormalized(); + hRecCascInvMass->DrawNormalized("same"); + auto leg2 = new TLegend(0.5, 0.5, 0.8, 0.8); + leg2->AddEntry(hStrangeTrackInvMass, "After tracking"); + leg2->AddEntry(hRecCascInvMass, "Before tracking"); + leg2->Draw(); + cv5.Write(); + + hRecXiRad->Sumw2(); + hRecXiRad->Divide(hGenXiRadius); + auto cv6 = TCanvas("reco_efficiency_r", "", 1000, 1000); + hRecXiRad->GetXaxis()->SetTitle("Radius (cm)"); + hRecXiRad->GetYaxis()->SetTitle("Efficiency"); + hRecXiRad->SetLineColor(kBlue); + hRecXiRad->Draw(); + cv6.Write(); + + hRecXiMom->Sumw2(); + hRecXiMom->Divide(hGenXiMom); + auto cv7 = TCanvas("reco_efficiency_pt", "", 1000, 1000); + hRecXiMom->GetYaxis()->SetTitle("Efficiency"); + hRecXiMom->SetLineColor(kBlue); + hRecXiMom->Draw("pe"); + cv7.Write(); + + outFile.Close(); +} +std::array matchCascToMC(const std::vector>& mcTracksMatrix, std::map*>& map, std::vector* v0vec, Cascade& casc, bool& isV0reco) +{ + std::array motherVec{-1, -1}; + std::array, 2> v0DauRefs; + std::array bachRef; + + auto v0Idx = casc.getV0ID(); + auto& v0 = v0vec->at(v0Idx); + + auto bachID = casc.getBachelorID(); + if (!map[bachID.getSourceName()]) + return motherVec; + auto& bachLab = map[bachID.getSourceName()]->at(bachID.getIndex()); + if (bachLab.isValid()) + bachRef = {bachLab.getEventID(), bachLab.getTrackID()}; + else + return motherVec; + + for (unsigned int iV0 = 0; iV0 < 2; iV0++) { + v0DauRefs[iV0] = {-1, -1}; + if (map[v0.getProngID(iV0).getSourceName()]) { + auto labTrackType = map[v0.getProngID(iV0).getSourceName()]; + auto lab = labTrackType->at(v0.getProngID(iV0).getIndex()); + + int trackID, evID, srcID; + bool fake; + lab.get(trackID, evID, srcID, fake); + if (lab.isValid()) { + v0DauRefs[iV0] = {lab.getEventID(), lab.getTrackID()}; + } + } + } + + if (v0DauRefs[0][1] == -1 || v0DauRefs[1][1] == -1) + return motherVec; + + auto& dau1MC = mcTracksMatrix[v0DauRefs[0][0]][v0DauRefs[0][1]]; + auto& dau2MC = mcTracksMatrix[v0DauRefs[1][0]][v0DauRefs[1][1]]; + + if (!(std::abs(dau1MC.GetPdgCode()) == firstV0dauPDG && std::abs(dau2MC.GetPdgCode()) == secondV0dauPDG) && !(std::abs(dau1MC.GetPdgCode()) == secondV0dauPDG && std::abs(dau2MC.GetPdgCode()) == firstV0dauPDG)) + return motherVec; + + if (!dau1MC.isSecondary() || !dau2MC.isSecondary() || dau1MC.getMotherTrackId() != dau2MC.getMotherTrackId()) + return motherVec; + + auto v0MC = mcTracksMatrix[v0DauRefs[0][0]][dau1MC.getMotherTrackId()]; + auto& bachMC = mcTracksMatrix[bachRef[0]][bachRef[1]]; + + if (v0MC.getMotherTrackId() >= 0) { + + if (std::abs(mcTracksMatrix[v0DauRefs[0][0]][v0MC.getMotherTrackId()].GetPdgCode()) == motherPDG) + isV0reco = true; + } + + if (std::abs(v0MC.GetPdgCode()) != v0PDG || !v0MC.isSecondary() || !bachMC.isSecondary()) + return motherVec; + + auto cascMC = mcTracksMatrix[v0DauRefs[0][0]][v0MC.getMotherTrackId()]; + if (v0MC.getMotherTrackId() != bachMC.getMotherTrackId() || std::abs(cascMC.GetPdgCode()) != motherPDG) + return motherVec; + + motherVec = {v0DauRefs[0][0], v0MC.getMotherTrackId()}; + return motherVec; +} + +std::array checkV0mother(const std::vector>& mcTracksMatrix, std::map*>& map, V0& v0) +{ + std::array motherVec{-1, -1}; + std::array, 2> v0DauRefs; + + for (unsigned int iV0 = 0; iV0 < 2; iV0++) { + v0DauRefs[iV0] = {-1, -1}; + if (map[v0.getProngID(iV0).getSourceName()]) { + auto labTrackType = map[v0.getProngID(iV0).getSourceName()]; + auto lab = labTrackType->at(v0.getProngID(iV0).getIndex()); + + int trackID, evID, srcID; + bool fake; + lab.get(trackID, evID, srcID, fake); + if (lab.isValid()) { + v0DauRefs[iV0] = {lab.getEventID(), lab.getTrackID()}; + } + } + } + + if (v0DauRefs[0][1] == -1 || v0DauRefs[1][1] == -1) + return motherVec; + + auto& dau1MC = mcTracksMatrix[v0DauRefs[0][0]][v0DauRefs[0][1]]; + auto& dau2MC = mcTracksMatrix[v0DauRefs[1][0]][v0DauRefs[1][1]]; + + if (!(std::abs(dau1MC.GetPdgCode()) == firstV0dauPDG && std::abs(dau2MC.GetPdgCode()) == secondV0dauPDG) && !(std::abs(dau1MC.GetPdgCode()) == secondV0dauPDG && std::abs(dau2MC.GetPdgCode()) == firstV0dauPDG)) + return motherVec; + + if (!dau1MC.isSecondary() || !dau2MC.isSecondary() || dau1MC.getMotherTrackId() != dau2MC.getMotherTrackId()) + return motherVec; + + LOG(info) << "-----------------------------------------"; + LOG(info) << "V0 daughters: " << dau1MC.GetPdgCode() << " " << dau2MC.GetPdgCode(); + LOG(info) << " Dau refs: " << dau1MC.getMotherTrackId() << " " << dau2MC.getMotherTrackId(); + + auto v0MC = mcTracksMatrix[v0DauRefs[0][0]][dau1MC.getMotherTrackId()]; + LOG(info) << "V0 mother: " << v0MC.GetPdgCode() << ", Ref: " << v0MC.getMotherTrackId(); + + auto cascMC = mcTracksMatrix[v0DauRefs[0][0]][v0MC.getMotherTrackId()]; + if (std::abs(cascMC.GetPdgCode()) != motherPDG) + return motherVec; + + LOG(info) << "Casc from V0 PDG: " << mcTracksMatrix[v0DauRefs[0][0]][v0MC.getMotherTrackId()].GetPdgCode(); + + motherVec = {v0DauRefs[0][0], v0MC.getMotherTrackId()}; + return motherVec; +} + +double calcDecLength(std::vector* MCTracks, const MCTrack& motherTrack) +{ + auto idStart = motherTrack.getFirstDaughterTrackId(); + auto idStop = motherTrack.getLastDaughterTrackId(); + + if (idStart == -1 || idStop == -1) + return -1; + for (auto iD{idStart}; iD <= idStop; ++iD) { + auto dauTrack = MCTracks->at(iD); + if (std::abs(dauTrack.GetPdgCode()) == v0PDG) { + auto decLength = (dauTrack.GetStartVertexCoordinatesX() - + motherTrack.GetStartVertexCoordinatesX()) * + (dauTrack.GetStartVertexCoordinatesX() - + motherTrack.GetStartVertexCoordinatesX()) + + (dauTrack.GetStartVertexCoordinatesY() - + motherTrack.GetStartVertexCoordinatesY()) * + (dauTrack.GetStartVertexCoordinatesY() - + motherTrack.GetStartVertexCoordinatesY()); + return sqrt(decLength); + } + } + return -1; +} + +std::array matchCompLabelToMC(const std::vector>& mcTracksMatrix, o2::MCCompLabel compLabel) +{ + std::array compRef = {-1, -1}; + int trackID, evID, srcID; + bool fake; + compLabel.get(trackID, evID, srcID, fake); + if (compLabel.isValid()) { + compRef = {evID, trackID}; + } + return compRef; +} + +std::array matchITStracktoMC(const std::vector>& mcTracksMatrix, + o2::MCCompLabel ITSlabel) + +{ + std::array outArray = {-1, -1}; + int trackID, evID, srcID; + bool fake; + ITSlabel.get(trackID, evID, srcID, fake); + if (ITSlabel.isValid() && + std::abs(mcTracksMatrix[evID][trackID].GetPdgCode()) == motherPDG) { + outArray = {evID, trackID}; + } + + return outArray; +} + +double calcMass(std::vector tracks) +{ + TLorentzVector moth, prong; + std::array p; + for (unsigned int i = 0; i < tracks.size(); i++) { + auto& track = tracks[i]; + auto mass = tracks.size() == 2 ? kFirstDauMasses[i] : kDauMasses[i]; + track.getPxPyPzGlo(p); + prong.SetVectM({p[0], p[1], p[2]}, mass); + moth += prong; + } + return moth.M(); +} \ No newline at end of file diff --git a/Detectors/StrangenessTracking/tracking/CMakeLists.txt b/Detectors/StrangenessTracking/tracking/CMakeLists.txt new file mode 100644 index 0000000000000..c5fceb6888f29 --- /dev/null +++ b/Detectors/StrangenessTracking/tracking/CMakeLists.txt @@ -0,0 +1,32 @@ +# Copyright 2019-2020 CERN and copyright holders of ALICE O2. +# See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +# All rights not expressly granted are reserved. +# +# This software is distributed under the terms of the GNU General Public +# License v3 (GPL Version 3), copied verbatim in the file "COPYING". +# +# In applying this license CERN 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. + +o2_add_library(StrangenessTracking + SOURCES src/StrangenessTracker.cxx + src/StrangenessTrackingConfigParam.cxx + PUBLIC_LINK_LIBRARIES O2::MathUtils + O2::ITSBase + O2::ITSMFTReconstruction + O2::ITSMFTSimulation + O2::DataFormatsITSMFT + O2::DataFormatsITS + O2::SimulationDataFormat + O2::DetectorsVertexing + O2::ReconstructionDataFormats + ) + +o2_target_root_dictionary(StrangenessTracking + HEADERS include/StrangenessTracking/StrangenessTrackingConfigParam.h + include/StrangenessTracking/IndexTableUtils.h + include/StrangenessTracking/StrangenessTracker.h + + + LINKDEF src/StrangenessTrackingLinkDef.h) \ No newline at end of file diff --git a/Detectors/StrangenessTracking/tracking/include/StrangenessTracking/IndexTableUtils.h b/Detectors/StrangenessTracking/tracking/include/StrangenessTracking/IndexTableUtils.h new file mode 100644 index 0000000000000..bfa2362155983 --- /dev/null +++ b/Detectors/StrangenessTracking/tracking/include/StrangenessTracking/IndexTableUtils.h @@ -0,0 +1,88 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN 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. + +/// \file IndexTableUtils.h +/// \brief +/// +#ifndef STRTRACKING_INCLUDE_INDEXTABLEUTILS_H_ +#define STRTRACKING_INCLUDE_INDEXTABLEUTILS_H_ + +#include "TMath.h" + +namespace o2 +{ +namespace strangeness_tracking +{ + +struct IndexTableUtils { + int getEtaBin(float eta); + int getPhiBin(float phi); + int getBinIndex(float eta, float phi); + std::vector getBinRect(float eta, float phi, float deltaEta, float deltaPhi); + int mEtaBins = 64, mPhiBins = 64; + float minEta = -1.5, maxEta = 1.5; + float minPhi = 0., maxPhi = 2 * TMath::Pi(); +}; + +inline int IndexTableUtils::getEtaBin(float eta) +{ + float deltaEta = (maxEta - minEta) / (mEtaBins); + int bEta = (eta - minEta) / deltaEta; // bins recentered to 0 + return bEta; +}; + +inline int IndexTableUtils::getPhiBin(float phi) +{ + float deltaPhi = (maxPhi - minPhi) / (mPhiBins); + int bPhi = (phi - minPhi) / deltaPhi; // bin recentered to 0 + return bPhi; +} + +inline int IndexTableUtils::getBinIndex(float eta, float phi) +{ + float deltaPhi = (maxPhi - minPhi) / (mPhiBins); + float deltaEta = (maxEta - minEta) / (mEtaBins); + int bEta = getEtaBin(eta); + int bPhi = getPhiBin(phi); + return (bEta >= mEtaBins || bPhi >= mPhiBins || bEta < 0 || bPhi < 0) ? mEtaBins * mPhiBins : bEta + mEtaBins * bPhi; +} + +inline std::vector IndexTableUtils::getBinRect(float eta, float phi, float deltaEta, float deltaPhi) +{ + std::vector idxVec; + int centralBin = getBinIndex(eta, phi); + if (centralBin == mPhiBins * mEtaBins) { // condition for overflows + idxVec.push_back(centralBin); + return idxVec; + } + int minEtaBin = TMath::Max(0, getEtaBin(eta - deltaEta)); + int maxEtaBin = getEtaBin(eta + deltaEta); + int minPhiBin = TMath::Max(0, getPhiBin(phi - deltaPhi)); + int maxPhiBin = getPhiBin(phi + deltaPhi); + + for (int iPhi{minPhiBin}; iPhi <= maxPhiBin; iPhi++) { + if (iPhi >= mPhiBins) { + break; + } + for (int iEta{minEtaBin}; iEta <= maxEtaBin; iEta++) { + if (iEta >= mEtaBins) { + break; + } + idxVec.push_back(iEta + mEtaBins * iPhi); + } + } + return idxVec; +}; + +} // namespace strangeness_tracking +} // namespace o2 + +#endif \ No newline at end of file diff --git a/Detectors/StrangenessTracking/tracking/include/StrangenessTracking/StrangenessTracker.h b/Detectors/StrangenessTracking/tracking/include/StrangenessTracking/StrangenessTracker.h new file mode 100644 index 0000000000000..2ba048d8c028c --- /dev/null +++ b/Detectors/StrangenessTracking/tracking/include/StrangenessTracking/StrangenessTracker.h @@ -0,0 +1,140 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN 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. + +/// \file StrangenessTracker.h +/// \brief +/// + +#ifndef _ALICEO2_STRANGENESS_TRACKER_ +#define _ALICEO2_STRANGENESS_TRACKER_ + +#include +#include +#include "TMath.h" +#include "StrangenessTracking/IndexTableUtils.h" +#include "StrangenessTracking/StrangenessTrackingConfigParam.h" +#include "ReconstructionDataFormats/PID.h" +#include "ReconstructionDataFormats/V0.h" +#include "ReconstructionDataFormats/Cascade.h" + +#include "DataFormatsITS/TrackITS.h" +#include "ITSBase/GeometryTGeo.h" +#include "ReconstructionDataFormats/Track.h" +#include "DataFormatsITSMFT/CompCluster.h" + +#include "DetectorsVertexing/DCAFitterN.h" +#include "DetectorsBase/Propagator.h" + +namespace o2 +{ +namespace strangeness_tracking +{ + +enum kPartType { kV0, + kCascade, + kThreeBody }; + +struct ClusAttachments { + + std::array arr; +}; + +struct StrangeTrack { + kPartType mPartType; + o2::track::TrackParCovF mMother; + int mITSRef = -1; + int mDecayRef = -1; + std::array decayVtx; + std::array decayMom; + float mInvMass; + float mMatchChi2; + float mTopoChi2; +}; + +class StrangenessTracker +{ + public: + using PID = o2::track::PID; + using TrackITS = o2::its::TrackITS; + using ITSCluster = o2::BaseCluster; + using V0 = o2::dataformats::V0; + using Cascade = o2::dataformats::Cascade; + using GIndex = o2::dataformats::VtxTrackIndex; + using DCAFitter2 = o2::vertexing::DCAFitterN<2>; + using DCAFitter3 = o2::vertexing::DCAFitterN<3>; + + StrangenessTracker() = default; + ~StrangenessTracker() = default; + + void initialise(); + void process(); + + std::vector& getClusAttachments() { return mClusAttachments; }; + std::vector& getStrangeTrackVec() { return mStrangeTrackVec; }; + + float getBz() const { return mBz; } + void setBz(float d) { mBz = d; } + void setCorrType(const o2::base::PropagatorImpl::MatCorrType& type) { mCorrType = type; } + + void setupFitters() + { + mFitterV0.setBz(mBz); + mFitter3Body.setBz(mBz); + mFitterV0.setUseAbsDCA(true); + mFitter3Body.setUseAbsDCA(true); + } + + bool loadData(gsl::span InputITStracks, std::vector& InputITSclusters, gsl::span InputITSidxs, gsl::span InputV0tracks, gsl::span InputCascadeTracks, o2::its::GeometryTGeo* geomITS); + double calcV0alpha(const V0& v0); + std::vector getTrackClusters(); + float getMatchingChi2(o2::track::TrackParCovF, const TrackITS ITSTrack, ITSCluster matchingClus); + bool recreateV0(const o2::track::TrackParCov& posTrack, const o2::track::TrackParCov& negTrack, V0& newV0); + + bool updateTrack(const ITSCluster& clus, o2::track::TrackParCov& track); + bool matchDecayToITStrack(float decayR); + + protected: + gsl::span mInputITStracks; // input ITS tracks + std::vector mTracksIdxTable; // index table for ITS tracks + std::vector mInputITSclusters; // input ITS clusters + gsl::span mInputITSidxs; // input ITS track-cluster indexes + gsl::span mInputV0tracks; // input V0 of decay daughters + gsl::span mInputCascadeTracks; // input V0 of decay daughters + + std::vector mSortedITStracks; // sorted ITS tracks + std::vector mSortedITSindexes; // indexes of sorted ITS tracks + IndexTableUtils mUtils; // structure for computing eta/phi matching selections + + std::vector mStrangeTrackVec; // structure containing updated mother and daughter tracks + std::vector mClusAttachments; // # of attached tracks, 1 for mother, 2 for daughter + + const StrangenessTrackingParamConfig* mStrParams = nullptr; + float mBz = -5; // Magnetic field + + DCAFitter2 mFitterV0; // optional DCA Fitter for recreating V0 with hypertriton mass hypothesis + DCAFitter3 mFitter3Body; // optional DCA Fitter for final 3 Body refit + + o2::base::PropagatorImpl::MatCorrType mCorrType = o2::base::PropagatorImpl::MatCorrType::USEMatCorrNONE; // use mat correction + o2::its::GeometryTGeo* mGeomITS; // ITS geometry + + std::vector mDaughterTracks; // vector of daughter tracks + StrangeTrack mStrangeTrack; // structure containing updated mother and daughter track refs + ClusAttachments mStructClus; // # of attached tracks, 1 for mother, 2 for daughter + o2::its::TrackITS mITStrack; // ITS track + std::array mV0dauIDs; // V0 daughter IDs + + ClassDefNV(StrangenessTracker, 1); +}; + +} // namespace strangeness_tracking +} // namespace o2 + +#endif // _ALICEO2_STRANGENESS_TRACKER_ diff --git a/Detectors/StrangenessTracking/tracking/include/StrangenessTracking/StrangenessTrackingConfigParam.h b/Detectors/StrangenessTracking/tracking/include/StrangenessTracking/StrangenessTrackingConfigParam.h new file mode 100644 index 0000000000000..5af773122691b --- /dev/null +++ b/Detectors/StrangenessTracking/tracking/include/StrangenessTracking/StrangenessTrackingConfigParam.h @@ -0,0 +1,40 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN 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. +/// \file StrangenessTrackingConfigParam.h +/// \brief + +#ifndef ALICEO2_STRANGENESS_TRACKING_PARAM_H_ +#define ALICEO2_STRANGENESS_TRACKING_PARAM_H_ + +#include "CommonUtils/ConfigurableParam.h" +#include "CommonUtils/ConfigurableParamHelper.h" + +namespace o2 +{ +namespace strangeness_tracking +{ + +struct StrangenessTrackingParamConfig : public o2::conf::ConfigurableParamHelper { + + // parameters + float mRadiusTolIB = .3; // Radius tolerance for matching V0s in the IB + float mRadiusTolOB = .1; // Radius tolerance for matching V0s in the OB + float mPhiBinSize = 0.1; // Phi bin size for the matching grid + float mEtaBinSize = 0.1; // Eta bin size for the matching grid + float mMinMotherClus = 3.; // minimum number of cluster to be attached to the mother + float mMaxChi2 = 50; // Maximum matching chi2 + + O2ParamDef(StrangenessTrackingParamConfig, "strtracker"); +}; + +} // namespace strangeness_tracking +} // namespace o2 +#endif \ No newline at end of file diff --git a/Detectors/StrangenessTracking/tracking/src/StrangenessTracker.cxx b/Detectors/StrangenessTracking/tracking/src/StrangenessTracker.cxx new file mode 100644 index 0000000000000..6e01694bdaa01 --- /dev/null +++ b/Detectors/StrangenessTracking/tracking/src/StrangenessTracker.cxx @@ -0,0 +1,374 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN 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. +/// \file StrangenessTracker.cxx +/// \brief + +#include +#include "StrangenessTracking/StrangenessTracker.h" + +namespace o2 +{ +namespace strangeness_tracking +{ + +bool StrangenessTracker::loadData(gsl::span InputITStracks, std::vector& InputITSclusters, gsl::span InputITSidxs, gsl::span InputV0tracks, gsl::span InputCascadeTracks, o2::its::GeometryTGeo* geomITS) +{ + mInputV0tracks = InputV0tracks; + mInputCascadeTracks = InputCascadeTracks; + mInputITStracks = InputITStracks; + mInputITSclusters = InputITSclusters; + mInputITSidxs = InputITSidxs; + LOG(info) << "all tracks loaded"; + LOG(info) << "V0 tracks size: " << mInputV0tracks.size(); + LOG(info) << "Cascade tracks size: " << mInputCascadeTracks.size(); + LOG(info) << "ITS tracks size: " << mInputITStracks.size(); + LOG(info) << "ITS clusters size: " << mInputITSclusters.size(); + LOG(info) << "ITS idxs size: " << mInputITSidxs.size(); + mGeomITS = geomITS; + setupFitters(); + return true; +} + +void StrangenessTracker::initialise() +{ + mTracksIdxTable.clear(); + mSortedITStracks.clear(); + mSortedITSindexes.clear(); + + for (int iTrack{0}; iTrack < mInputITStracks.size(); iTrack++) { + mSortedITStracks.push_back(mInputITStracks[iTrack]); + mSortedITSindexes.push_back(iTrack); + } + + mTracksIdxTable.resize(mUtils.mPhiBins * mUtils.mEtaBins + 1); + std::sort(mSortedITStracks.begin(), mSortedITStracks.end(), [&](o2::its::TrackITS& a, o2::its::TrackITS& b) { return mUtils.getBinIndex(a.getEta(), a.getPhi()) < mUtils.getBinIndex(b.getEta(), b.getPhi()); }); + std::sort(mSortedITSindexes.begin(), mSortedITSindexes.end(), [&](int i, int j) { return mUtils.getBinIndex(mInputITStracks[i].getEta(), mInputITStracks[i].getPhi()) < mUtils.getBinIndex(mInputITStracks[j].getEta(), mInputITStracks[j].getPhi()); }); + + for (auto& track : mSortedITStracks) { + mTracksIdxTable[mUtils.getBinIndex(track.getEta(), track.getPhi())]++; + } + std::exclusive_scan(mTracksIdxTable.begin(), mTracksIdxTable.begin() + mUtils.mPhiBins * mUtils.mEtaBins, mTracksIdxTable.begin(), 0); + mTracksIdxTable[mUtils.mPhiBins * mUtils.mEtaBins] = mSortedITStracks.size(); + + // create config param instance + mStrParams = &StrangenessTrackingParamConfig::Instance(); +} + +void StrangenessTracker::process() +{ + // Loop over V0s + mDaughterTracks.resize(2); // resize to 2 prongs + + for (int iV0{0}; iV0 < mInputV0tracks.size(); iV0++) { + LOG(debug) << "Analysing V0: " << iV0 + 1 << "/" << mInputV0tracks.size(); + auto& DecIndexRef = iV0; + auto& v0 = mInputV0tracks[iV0]; + mV0dauIDs[0] = v0.getProngID(0), mV0dauIDs[1] = v0.getProngID(1); + auto posTrack = v0.getProng(0); + auto negTrack = v0.getProng(1); + auto alphaV0 = calcV0alpha(v0); + alphaV0 > 0 ? posTrack.setAbsCharge(2) : negTrack.setAbsCharge(2); + V0 correctedV0; // recompute V0 for Hypertriton + + if (!recreateV0(posTrack, negTrack, correctedV0)) { + continue; + } + + mStrangeTrack.mPartType = kV0; + + auto v0R2 = v0.calcR2(); + auto iBinsV0 = mUtils.getBinRect(correctedV0.getEta(), correctedV0.getPhi(), mStrParams->mEtaBinSize, mStrParams->mPhiBinSize); + for (int& iBinV0 : iBinsV0) { + for (int iTrack{mTracksIdxTable[iBinV0]}; iTrack < TMath::Min(mTracksIdxTable[iBinV0 + 1], int(mSortedITStracks.size())); iTrack++) { + mStrangeTrack.mMother = (o2::track::TrackParCovF)correctedV0; + mDaughterTracks[0] = correctedV0.getProng(0); + mDaughterTracks[1] = correctedV0.getProng(1); + mITStrack = mSortedITStracks[iTrack]; + + LOG(debug) << "V0 pos: " << v0.getProngID(0) << " V0 neg: " << v0.getProngID(1) << ", ITS track ref: " << mSortedITSindexes[iTrack]; + + if (matchDecayToITStrack(sqrt(v0R2))) { + LOG(debug) << "ITS Track matched with a V0 decay topology ...."; + LOG(debug) << "Number of ITS track clusters attached: " << mITStrack.getNumberOfClusters(); + mStrangeTrack.mDecayRef = iV0; + mStrangeTrack.mITSRef = mSortedITSindexes[iTrack]; + mStrangeTrackVec.push_back(mStrangeTrack); + mClusAttachments.push_back(mStructClus); + } + } + } + } + + // Loop over Cascades + mDaughterTracks.resize(3); // resize to 3 prongs + + for (int iCasc{0}; iCasc < mInputCascadeTracks.size(); iCasc++) { + LOG(debug) << "Analysing Cascade: " << iCasc + 1 << "/" << mInputCascadeTracks.size(); + auto& DecIndexRef = iCasc; + auto& casc = mInputCascadeTracks[iCasc]; + auto& cascV0 = mInputV0tracks[casc.getV0ID()]; + mV0dauIDs[0] = cascV0.getProngID(0), mV0dauIDs[1] = cascV0.getProngID(1); + + mStrangeTrack.mPartType = kCascade; + // first: bachelor, second: V0 pos, third: V0 neg + auto cascR2 = casc.calcR2(); + auto iBinsCasc = mUtils.getBinRect(casc.getEta(), casc.getPhi(), mStrParams->mEtaBinSize, mStrParams->mPhiBinSize); + for (int& iBinCasc : iBinsCasc) { + for (int iTrack{mTracksIdxTable[iBinCasc]}; iTrack < TMath::Min(mTracksIdxTable[iBinCasc + 1], int(mSortedITStracks.size())); iTrack++) { + mStrangeTrack.mMother = (o2::track::TrackParCovF)casc; + mDaughterTracks[0] = casc.getBachelorTrack(), mDaughterTracks[1] = cascV0.getProng(0), mDaughterTracks[2] = cascV0.getProng(1); + mITStrack = mSortedITStracks[iTrack]; + auto& ITSindexRef = mSortedITSindexes[iTrack]; + LOG(debug) << "----------------------"; + LOG(debug) << "CascV0: " << casc.getV0ID() << ", Bach ID: " << casc.getBachelorID() << ", ITS track ref: " << mSortedITSindexes[iTrack]; + if (matchDecayToITStrack(sqrt(cascR2))) { + LOG(debug) << "ITS Track matched with a Cascade decay topology ...."; + LOG(debug) << "Number of ITS track clusters attached: " << mITStrack.getNumberOfClusters(); + mStrangeTrack.mDecayRef = iCasc; + mStrangeTrack.mITSRef = mSortedITSindexes[iTrack]; + mStrangeTrackVec.push_back(mStrangeTrack); + mClusAttachments.push_back(mStructClus); + } + } + } + } +} + +bool StrangenessTracker::matchDecayToITStrack(float decayR) +{ + + auto trackClusters = getTrackClusters(); + auto& lastClus = trackClusters[0]; + mStrangeTrack.mMatchChi2 = getMatchingChi2(mStrangeTrack.mMother, mITStrack, lastClus); + + auto radTol = decayR < 4 ? mStrParams->mRadiusTolIB : mStrParams->mRadiusTolOB; + auto nMinClusMother = trackClusters.size() < 4 ? 2 : mStrParams->mMinMotherClus; + + std::vector motherClusters; + std::array nAttachments; + + int nUpdates = 0; + bool isMotherUpdated = false; + + for (auto& clus : trackClusters) { + int nUpdOld = nUpdates; + double clusRad = sqrt(clus.getX() * clus.getX() - clus.getY() * clus.getY()); + auto diffR = decayR - clusRad; + auto relDiffR = diffR / decayR; + // Look for the Mother if the Decay radius allows for it, within a tolerance + LOG(debug) << "++++++++"; + LOG(debug) << "decayR: " << decayR << ", diffR: " << diffR << ", clus rad: " << clusRad << ", radTol: " << radTol; + if (relDiffR > -radTol) { + LOG(debug) << "Try to attach cluster to Mother, layer: " << mGeomITS->getLayer(clus.getSensorID()); + if (updateTrack(clus, mStrangeTrack.mMother)) { + motherClusters.push_back(clus); + nAttachments[mGeomITS->getLayer(clus.getSensorID())] = 0; + isMotherUpdated = true; + nUpdates++; + LOG(debug) << "Cluster attached to Mother"; + continue; // if the cluster is attached to the mother, skip the rest of the loop + } + } + + // if Mother is not found, check for V0 daughters compatibility + if (relDiffR < radTol && !isMotherUpdated) { + bool isDauUpdated = false; + LOG(debug) << "Try to attach cluster to Daughters, layer: " << mGeomITS->getLayer(clus.getSensorID()); + for (int iDau{0}; iDau < mDaughterTracks.size(); iDau++) { + auto& dauTrack = mDaughterTracks[iDau]; + if (updateTrack(clus, dauTrack)) { + nAttachments[mGeomITS->getLayer(clus.getSensorID())] = iDau + 1; + isDauUpdated = true; + break; + } + } + if (!isDauUpdated) { + break; // no daughter track updated, stop the loop + } + nUpdates++; + } + if (nUpdates == nUpdOld) { + break; // no track updated, stop the loop + } + } + + if (nUpdates < trackClusters.size() || motherClusters.size() < nMinClusMother) { + return false; + } + + o2::track::TrackParCov motherTrackClone = mStrangeTrack.mMother; // clone and reset covariance for final topology refit + motherTrackClone.resetCovariance(); + + LOG(debug) << "Clusters attached, starting inward-outward refit"; + + std::reverse(motherClusters.begin(), motherClusters.end()); + for (auto& clus : motherClusters) { + if (!updateTrack(clus, motherTrackClone)) { + break; + } + } + LOG(debug) << "Inward-outward refit finished, starting final topology refit"; + + // final Topology refit + + int cand = 0; // best V0 candidate + int nCand; + + // refit cascade + if (mStrangeTrack.mPartType == kCascade) { + V0 cascV0Upd; + if (!recreateV0(mDaughterTracks[1], mDaughterTracks[2], cascV0Upd)) { + LOG(debug) << "Cascade V0 refit failed"; + return false; + } + try { + nCand = mFitter3Body.process(cascV0Upd, mDaughterTracks[0], motherTrackClone); + } catch (std::runtime_error& e) { + LOG(debug) << "Fitter3Body failed: " << e.what(); + return false; + } + if (!nCand || !mFitter3Body.propagateTracksToVertex()) { + LOG(debug) << "Fitter3Body failed: propagation to vertex failed"; + return false; + } + } + + // refit V0 + else if (mStrangeTrack.mPartType == kV0) { + try { + nCand = mFitter3Body.process(mDaughterTracks[0], mDaughterTracks[1], motherTrackClone); + } catch (std::runtime_error& e) { + LOG(debug) << "Fitter3Body failed: " << e.what(); + return false; + } + if (!nCand || !mFitter3Body.propagateTracksToVertex()) { + LOG(debug) << "Fitter3Body failed: propagation to vertex failed"; + return false; + } + } + + mStrangeTrack.decayVtx = mFitter3Body.getPCACandidatePos(); + mStrangeTrack.mTopoChi2 = mFitter3Body.getChi2AtPCACandidate(); + mStructClus.arr = nAttachments; + + return true; +} + +bool StrangenessTracker::updateTrack(const ITSCluster& clus, o2::track::TrackParCov& track) +{ + auto propInstance = o2::base::Propagator::Instance(); + float alpha = mGeomITS->getSensorRefAlpha(clus.getSensorID()), x = clus.getX(); + int layer{mGeomITS->getLayer(clus.getSensorID())}; + + auto stringOld = track.asString(); + LOG(debug) << "Track before update, Y2: " << track.getSigmaY2() << ", Z2: " << track.getSigmaZ2(); + + if (!track.rotate(alpha)) { + return false; + } + + LOG(debug) << "Track rotated, Y2: " << track.getSigmaY2() << ", Z2: " << track.getSigmaZ2(); + + if (!propInstance->propagateToX(track, x, getBz(), o2::base::PropagatorImpl::MAX_SIN_PHI, o2::base::PropagatorImpl::MAX_STEP, mCorrType)) { + return false; + } + + auto stringNew = track.asString(); + LOG(debug) << "Track after propagation, Y2: " << track.getSigmaY2() << ", Z2: " << track.getSigmaZ2(); + + if (mCorrType == o2::base::PropagatorF::MatCorrType::USEMatCorrNONE) { + float thick = layer < 3 ? 0.005 : 0.01; + constexpr float radl = 9.36f; // Radiation length of Si [cm] + constexpr float rho = 2.33f; // Density of Si [g/cm^3] + if (!track.correctForMaterial(thick, thick * rho * radl)) { + return false; + } + } + auto chi2 = std::abs(track.getPredictedChi2(clus)); // abs to be understood + LOG(debug) << "Chi2: " << chi2; + if (chi2 > mStrParams->mMaxChi2 || chi2 < 0) { + return false; + } + + if (!track.update(clus)) { + return false; + } + + return true; +} + +bool StrangenessTracker::recreateV0(const o2::track::TrackParCov& posTrack, const o2::track::TrackParCov& negTrack, V0& newV0) +{ + + int nCand; + try { + nCand = mFitterV0.process(posTrack, negTrack); + } catch (std::runtime_error& e) { + return false; + } + if (!nCand || !mFitterV0.propagateTracksToVertex()) { + return false; + } + + const auto& v0XYZ = mFitterV0.getPCACandidatePos(); + + auto& propPos = mFitterV0.getTrack(0, 0); + auto& propNeg = mFitterV0.getTrack(1, 0); + + std::array pP, pN; + propPos.getPxPyPzGlo(pP); + propNeg.getPxPyPzGlo(pN); + std::array pV0 = {pP[0] + pN[0], pP[1] + pN[1], pP[2] + pN[2]}; + newV0 = V0(v0XYZ, pV0, mFitterV0.calcPCACovMatrixFlat(0), propPos, propNeg, mV0dauIDs[0], mV0dauIDs[1], o2::track::PID::HyperTriton); + return true; +}; + +std::vector StrangenessTracker::getTrackClusters() +{ + std::vector outVec; + auto firstClus = mITStrack.getFirstClusterEntry(); + auto ncl = mITStrack.getNumberOfClusters(); + for (int icl = 0; icl < ncl; icl++) { + outVec.push_back(mInputITSclusters[mInputITSidxs[firstClus + icl]]); + } + return outVec; +}; + +float StrangenessTracker::getMatchingChi2(o2::track::TrackParCovF v0, const TrackITS ITStrack, ITSCluster matchingClus) +{ + float alpha = mGeomITS->getSensorRefAlpha(matchingClus.getSensorID()), x = matchingClus.getX(); + if (v0.rotate(alpha)) { + if (v0.propagateTo(x, mBz)) { + return v0.getPredictedChi2(ITStrack.getParamOut()); + } + } + return -100; +}; + +double StrangenessTracker::calcV0alpha(const V0& v0) +{ + std::array fV0mom, fPmom, fNmom = {0, 0, 0}; + v0.getProng(0).getPxPyPzGlo(fPmom); + v0.getProng(1).getPxPyPzGlo(fNmom); + v0.getPxPyPzGlo(fV0mom); + + TVector3 momNeg(fNmom[0], fNmom[1], fNmom[2]); + TVector3 momPos(fPmom[0], fPmom[1], fPmom[2]); + TVector3 momTot(fV0mom[0], fV0mom[1], fV0mom[2]); + + Double_t lQlNeg = momNeg.Dot(momTot) / momTot.Mag(); + Double_t lQlPos = momPos.Dot(momTot) / momTot.Mag(); + + return (lQlPos - lQlNeg) / (lQlPos + lQlNeg); +}; + +} // namespace strangeness_tracking +} // namespace o2 \ No newline at end of file diff --git a/Detectors/StrangenessTracking/tracking/src/StrangenessTrackingConfigParam.cxx b/Detectors/StrangenessTracking/tracking/src/StrangenessTrackingConfigParam.cxx new file mode 100644 index 0000000000000..9f324b9f5db0c --- /dev/null +++ b/Detectors/StrangenessTracking/tracking/src/StrangenessTrackingConfigParam.cxx @@ -0,0 +1,24 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN 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. +/// \file StrangenessTrackingConfigParam.cxx +/// \brief + +#include "StrangenessTracking/StrangenessTrackingConfigParam.h" + +namespace o2 +{ +namespace strangeness_tracking +{ +static auto& sStrangenessTrackingParam = o2::strangeness_tracking::StrangenessTrackingParamConfig::Instance(); + +O2ParamImpl(o2::strangeness_tracking::StrangenessTrackingParamConfig); +} // namespace strangeness_tracking +} // namespace o2 diff --git a/Detectors/StrangenessTracking/tracking/src/StrangenessTrackingLinkDef.h b/Detectors/StrangenessTracking/tracking/src/StrangenessTrackingLinkDef.h new file mode 100644 index 0000000000000..25221e0dd6ea9 --- /dev/null +++ b/Detectors/StrangenessTracking/tracking/src/StrangenessTrackingLinkDef.h @@ -0,0 +1,28 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN 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. + +#ifdef __CLING__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class o2::strangeness_tracking::StrangenessTrackingParamConfig + ; +#pragma link C++ class o2::conf::ConfigurableParamHelper < o2::strangeness_tracking::StrangenessTrackingParamConfig> + ; +#pragma link C++ struct o2::strangeness_tracking::ClusAttachments + ; +#pragma link C++ struct o2::strangeness_tracking::StrangeTrack + ; +#pragma link C++ class o2::strangeness_tracking::IndexTableUtils + ; + +#pragma link C++ class std::vector < o2::strangeness_tracking::ClusAttachments> + ; +#pragma link C++ class std::vector < o2::strangeness_tracking::StrangeTrack> + ; +#pragma link C++ class std::vector < o2::strangeness_tracking::IndexTableUtils> + ; + +#endif diff --git a/Detectors/StrangenessTracking/workflow/CMakeLists.txt b/Detectors/StrangenessTracking/workflow/CMakeLists.txt new file mode 100644 index 0000000000000..1954e8ea7548a --- /dev/null +++ b/Detectors/StrangenessTracking/workflow/CMakeLists.txt @@ -0,0 +1,38 @@ +# Copyright 2019-2020 CERN and copyright holders of ALICE O2. +# See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +# All rights not expressly granted are reserved. +# +# This software is distributed under the terms of the GNU General Public +# License v3 (GPL Version 3), copied verbatim in the file "COPYING". +# +# In applying this license CERN 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. + +o2_add_library(StrangenessTrackingWorkflow + TARGETVARNAME targetName + SOURCES src/StrangenessTrackingSpec.cxx + SOURCES src/StrangenessTrackingWriterSpec.cxx + PUBLIC_LINK_LIBRARIES O2::GlobalTrackingWorkflowReaders + O2::GlobalTrackingWorkflow + O2::ITSWorkflow + O2::StrangenessTracking + O2::SimulationDataFormat + O2::ReconstructionDataFormats + O2::Framework + O2::DetectorsRaw + O2::SimConfig + O2::DataFormatsITS + O2::DataFormatsDCS + O2::ITStracking + O2::ITSReconstruction + O2::ITSMFTReconstruction + O2::ITSMFTWorkflow + O2::DetectorsCalibration + O2::GlobalTrackingWorkflowWriters + O2::CCDB) + +o2_add_executable(strangeness-tracking-workflow + SOURCES src/strangeness-tracking-workflow.cxx + PUBLIC_LINK_LIBRARIES O2::StrangenessTrackingWorkflow + O2::GlobalTrackingWorkflowReaders) \ No newline at end of file diff --git a/Detectors/StrangenessTracking/workflow/include/StrangenessTrackingWorkflow/StrangenessTrackingSpec.h b/Detectors/StrangenessTracking/workflow/include/StrangenessTrackingWorkflow/StrangenessTrackingSpec.h new file mode 100644 index 0000000000000..9b109d2d8706b --- /dev/null +++ b/Detectors/StrangenessTracking/workflow/include/StrangenessTrackingWorkflow/StrangenessTrackingSpec.h @@ -0,0 +1,64 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN 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. +/// \file StrangenessTrackingSpec.h +/// \brief + +#ifndef O2_STRANGENESS_SPEC_H +#define O2_STRANGENESS_SPEC_H + +#include "Framework/WorkflowSpec.h" +#include "Framework/DataProcessorSpec.h" +#include "Framework/Task.h" +#include "DataFormatsITSMFT/CompCluster.h" +#include "DataFormatsParameters/GRPObject.h" +#include "DetectorsBase/GRPGeomHelper.h" +#include "ITSMFTReconstruction/ClustererParam.h" +#include "DataFormatsITSMFT/TopologyDictionary.h" + +#include "TStopwatch.h" + +#include "StrangenessTracking/StrangenessTracker.h" + +namespace o2 +{ +namespace strangeness_tracking +{ +class StrangenessTrackerSpec : public framework::Task +{ + public: + using ITSCluster = o2::BaseCluster; + + StrangenessTrackerSpec(std::shared_ptr gr, bool isMC); + ~StrangenessTrackerSpec() override = default; + + void init(framework::InitContext& ic) final; + void run(framework::ProcessingContext& pc) final; + void endOfStream(framework::EndOfStreamContext& ec) final; + void finaliseCCDB(framework::ConcreteDataMatcher& matcher, void* obj) final; + void setClusterDictionary(const o2::itsmft::TopologyDictionary* d) { mDict = d; } + + private: + void updateTimeDependentParams(framework::ProcessingContext& pc); + + bool mIsMC = false; + TStopwatch mTimer; + StrangenessTracker mTracker; + std::shared_ptr mGGCCDBRequest; + std::unique_ptr mGRP = nullptr; + const o2::itsmft::TopologyDictionary* mDict = nullptr; +}; + +o2::framework::DataProcessorSpec getStrangenessTrackerSpec(); +o2::framework::WorkflowSpec getWorkflow(bool upstreamClusters = false, bool upstreamV0s = false); + +} // namespace strangeness_tracking +} // namespace o2 +#endif \ No newline at end of file diff --git a/Detectors/StrangenessTracking/workflow/include/StrangenessTrackingWorkflow/StrangenessTrackingWriterSpec.h b/Detectors/StrangenessTracking/workflow/include/StrangenessTrackingWorkflow/StrangenessTrackingWriterSpec.h new file mode 100644 index 0000000000000..2cd3262370eb5 --- /dev/null +++ b/Detectors/StrangenessTracking/workflow/include/StrangenessTrackingWorkflow/StrangenessTrackingWriterSpec.h @@ -0,0 +1,33 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN 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. + +/// \file StrangenessTrackingWriterSpec.h +/// \brief +/// + +#ifndef O2_STRANGENESSTRACKINGWRITER +#define O2_STRANGENESSTRACKINGWRITER + +#include "Framework/DataProcessorSpec.h" + +namespace o2 +{ +namespace strangeness_tracking +{ + +/// create a processor spec +/// write ITS tracks to ROOT file +o2::framework::DataProcessorSpec getStrangenessTrackingWriterSpec(); + +} // namespace strangeness_tracking +} // namespace o2 + +#endif /* O2_STRANGENESSTRACKINGWRITER */ diff --git a/Detectors/StrangenessTracking/workflow/src/StrangenessTrackingSpec.cxx b/Detectors/StrangenessTracking/workflow/src/StrangenessTrackingSpec.cxx new file mode 100644 index 0000000000000..9974a9f55d1fd --- /dev/null +++ b/Detectors/StrangenessTracking/workflow/src/StrangenessTrackingSpec.cxx @@ -0,0 +1,216 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN 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. +/// \file StrangenessTrackingSpec.cxx +/// \brief + +#include "TGeoGlobalMagField.h" +#include "Framework/ConfigParamRegistry.h" +#include "Field/MagneticField.h" + +#include "StrangenessTrackingWorkflow/StrangenessTrackingSpec.h" +#include "ITSWorkflow/ClusterWriterSpec.h" +#include "ITSWorkflow/TrackerSpec.h" +#include "ITSWorkflow/TrackReaderSpec.h" +#include "ITSMFTWorkflow/ClusterReaderSpec.h" +#include "GlobalTrackingWorkflowReaders/SecondaryVertexReaderSpec.h" +#include "GlobalTrackingWorkflowReaders/TrackTPCITSReaderSpec.h" +#include "GlobalTrackingWorkflow/TOFMatcherSpec.h" +#include "Framework/CCDBParamSpec.h" +#include "DataFormatsParameters/GRPObject.h" + +#include "DataFormatsITSMFT/ROFRecord.h" +#include "DataFormatsITS/TrackITS.h" +#include "ReconstructionDataFormats/TrackTPCITS.h" +#include "SimulationDataFormat/MCCompLabel.h" +#include "ITStracking/IOUtils.h" +#include "DetectorsCommonDataFormats/DetectorNameConf.h" + +#include + +namespace o2 +{ +using namespace o2::framework; +namespace strangeness_tracking +{ + +framework::WorkflowSpec getWorkflow(bool useMC, bool useRootInput) +{ + framework::WorkflowSpec specs; + if (useRootInput) { + specs.emplace_back(o2::itsmft::getITSClusterReaderSpec(useMC, true)); + specs.emplace_back(o2::its::getITSTrackReaderSpec(useMC)); + specs.emplace_back(o2::vertexing::getSecondaryVertexReaderSpec()); + specs.emplace_back(o2::globaltracking::getTrackTPCITSReaderSpec(true)); + // auto src = o2::dataformats::GlobalTrackID::Source::ITSTPCTOF | o2::dataformats::GlobalTrackID::Source::ITSTPC | o2::dataformats::GlobalTrackID::Source::TPCTOF; + // specs.emplace_back(o2::globaltracking::getTOFMatcherSpec(src, true, false, false, 0)); + } + specs.emplace_back(getStrangenessTrackerSpec()); + return specs; +} + +StrangenessTrackerSpec::StrangenessTrackerSpec(std::shared_ptr gr, bool isMC) : mGGCCDBRequest(gr), mIsMC{isMC} +{ + // no ops +} + +void StrangenessTrackerSpec::init(framework::InitContext& ic) +{ + mTimer.Stop(); + mTimer.Reset(); + + // load propagator + + o2::base::GRPGeomHelper::instance().setRequest(mGGCCDBRequest); + + LOG(info) << "Initialized strangeness tracker..."; +} + +void StrangenessTrackerSpec::run(framework::ProcessingContext& pc) +{ + mTimer.Start(false); + LOG(info) << "Running strangeness tracker..."; + updateTimeDependentParams(pc); + // ITS + auto ITSclus = pc.inputs().get>("compClusters"); + auto ITSpatt = pc.inputs().get>("patterns"); + auto ITStracks = pc.inputs().get>("ITSTrack"); + auto ROFsInput = pc.inputs().get>("ROframes"); + auto ITSTrackClusIdx = pc.inputs().get>("trackITSClIdx"); + + // V0 + auto v0Vec = pc.inputs().get>("v0s"); + auto cascadeVec = pc.inputs().get>("cascs"); + auto tpcITSTracks = pc.inputs().get>("trackTPCITS"); + + // Monte Carlo + auto labITSTPC = pc.inputs().get>("trackITSTPCMCTR"); + // auto labTPCTOF = pc.inputs().get>("clsTOF_TPC_MCTR"); + // auto labITSTPCTOF = pc.inputs().get>("clsTOF_GLO_MCTR"); + auto labITS = pc.inputs().get>("trackITSMCTR"); + LOGF(debug, "ITSclus: %d \nITSpatt: %d \nITStracks: %d \nROFsInput: %d \nITSTrackClusIdx: %d \nTPCITStracks: %d \nv0s: %d \nlabITSTPC: %d\nlabITS: %d", + ITSclus.size(), + ITSpatt.size(), + ITStracks.size(), + ROFsInput.size(), + ITSTrackClusIdx.size(), + tpcITSTracks.size(), + v0Vec.size(), + cascadeVec.size(), + labITSTPC.size(), + // labTPCTOF.size(), + // labITSTPCTOF.size(), + labITS.size()); + // \nlabTPCTOF: %d\nlabITSTPCTOF: %d + + mTracker.setCorrType(o2::base::PropagatorImpl::MatCorrType::USEMatCorrLUT); + LOG(debug) << "Bz: " << o2::base::Propagator::Instance()->getNominalBz(); + mTracker.setBz(o2::base::Propagator::Instance()->getNominalBz()); + + auto pattIt = ITSpatt.begin(); + std::vector ITSClustersArray; + ITSClustersArray.reserve(ITSclus.size()); + o2::its::ioutils::convertCompactClusters(ITSclus, pattIt, ITSClustersArray, mDict); + auto geom = o2::its::GeometryTGeo::Instance(); + + mTracker.loadData(ITStracks, ITSClustersArray, ITSTrackClusIdx, v0Vec, cascadeVec, geom); + mTracker.initialise(); + mTracker.process(); + + pc.outputs().snapshot(Output{"STK", "STRTRACKS", 0, Lifetime::Timeframe}, mTracker.getStrangeTrackVec()); + pc.outputs().snapshot(Output{"STK", "CLUSUPDATES", 0, Lifetime::Timeframe}, mTracker.getClusAttachments()); + + mTimer.Stop(); +} + +///_______________________________________ +void StrangenessTrackerSpec::updateTimeDependentParams(ProcessingContext& pc) +{ + o2::base::GRPGeomHelper::instance().checkUpdates(pc); + static bool initOnceDone = false; + if (!initOnceDone) { // this params need to be queried only once + initOnceDone = true; + pc.inputs().get("cldict"); // just to trigger the finaliseCCDB + o2::its::GeometryTGeo* geom = o2::its::GeometryTGeo::Instance(); + geom->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::T2L, o2::math_utils::TransformType::T2GRot, o2::math_utils::TransformType::T2G)); + } +} + +///_______________________________________ +void StrangenessTrackerSpec::finaliseCCDB(ConcreteDataMatcher& matcher, void* obj) +{ + if (o2::base::GRPGeomHelper::instance().finaliseCCDB(matcher, obj)) { + return; + } + if (matcher == ConcreteDataMatcher("ITS", "CLUSDICT", 0)) { + LOG(info) << "cluster dictionary updated"; + setClusterDictionary((const o2::itsmft::TopologyDictionary*)obj); + return; + } +} + +void StrangenessTrackerSpec::endOfStream(framework::EndOfStreamContext& ec) +{ + LOGF(info, "Strangeness tracking total timing: Cpu: %.3e Real: %.3e s in %d slots", + mTimer.CpuTime(), mTimer.RealTime(), mTimer.Counter() - 1); +} + +DataProcessorSpec getStrangenessTrackerSpec() +{ + std::vector inputs; + + // ITS + inputs.emplace_back("compClusters", "ITS", "COMPCLUSTERS", 0, Lifetime::Timeframe); + inputs.emplace_back("patterns", "ITS", "PATTERNS", 0, Lifetime::Timeframe); + inputs.emplace_back("ROframes", "ITS", "CLUSTERSROF", 0, Lifetime::Timeframe); + inputs.emplace_back("trackITSClIdx", "ITS", "TRACKCLSID", 0, Lifetime::Timeframe); + inputs.emplace_back("ITSTrack", "ITS", "TRACKS", 0, Lifetime::Timeframe); + inputs.emplace_back("cldict", "ITS", "CLUSDICT", 0, Lifetime::Condition, ccdbParamSpec("ITS/Calib/ClusterDictionary")); + auto ggRequest = std::make_shared(false, // orbitResetTime + true, // GRPECS=true + false, // GRPLHCIF + true, // GRPMagField + true, // askMatLUT + o2::base::GRPGeomRequest::Aligned, // geometry + inputs, + true); + + // V0 + inputs.emplace_back("v0s", "GLO", "V0S", 0, Lifetime::Timeframe); // found V0s + inputs.emplace_back("v02pvrf", "GLO", "PVTX_V0REFS", 0, Lifetime::Timeframe); // prim.vertex -> V0s refs + inputs.emplace_back("cascs", "GLO", "CASCS", 0, Lifetime::Timeframe); // found Cascades + inputs.emplace_back("cas2pvrf", "GLO", "PVTX_CASCREFS", 0, Lifetime::Timeframe); // prim.vertex -> Cascades refs + + // TPC-ITS + inputs.emplace_back("trackTPCITS", "GLO", "TPCITS", 0, Lifetime::Timeframe); + + // TPC-TOF + // inputs.emplace_back("matchTPCTOF", "TOF", "MTC_TPC", 0, Lifetime::Timeframe); // Matching input type manually set to 0 + + // Monte Carlo + inputs.emplace_back("trackITSTPCMCTR", "GLO", "TPCITS_MC", 0, Lifetime::Timeframe); // MC truth + // inputs.emplace_back("clsTOF_GLO_MCTR", "TOF", "MCMTC_ITSTPC", 0, Lifetime::Timeframe); // MC truth + inputs.emplace_back("trackITSMCTR", "ITS", "TRACKSMCTR", 0, Lifetime::Timeframe); // MC truth + // inputs.emplace_back("clsTOF_TPC_MCTR", "TOF", "MCMTC_TPC", 0, Lifetime::Timeframe); // MC truth, // Matching input type manually set to 0 + + std::vector outputs; + outputs.emplace_back("STK", "STRTRACKS", 0, Lifetime::Timeframe); + outputs.emplace_back("STK", "CLUSUPDATES", 0, Lifetime::Timeframe); + + return DataProcessorSpec{ + "strangeness-tracker", + inputs, + outputs, + AlgorithmSpec{adaptFromTask(ggRequest, false)}, + Options{}}; +} + +} // namespace strangeness_tracking +} // namespace o2 \ No newline at end of file diff --git a/Detectors/StrangenessTracking/workflow/src/StrangenessTrackingWriterSpec.cxx b/Detectors/StrangenessTracking/workflow/src/StrangenessTrackingWriterSpec.cxx new file mode 100644 index 0000000000000..a7409b1cbb413 --- /dev/null +++ b/Detectors/StrangenessTracking/workflow/src/StrangenessTrackingWriterSpec.cxx @@ -0,0 +1,54 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN 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. + +/// \file StrangenessWriterSpec.cxx +/// \brief + +#include +#include "StrangenessTrackingWorkflow/StrangenessTrackingWriterSpec.h" +#include "DPLUtils/MakeRootTreeWriterSpec.h" +#include "CommonDataFormat/TimeStamp.h" +#include "CommonDataFormat/RangeReference.h" +#include "ReconstructionDataFormats/V0.h" +#include "ReconstructionDataFormats/Track.h" +#include "StrangenessTracking/StrangenessTracker.h" + +using namespace o2::framework; + +namespace o2 +{ +namespace strangeness_tracking +{ +using V0 = o2::dataformats::V0; + +template +using BranchDefinition = MakeRootTreeWriterSpec::BranchDefinition; + +using namespace o2::header; + +DataProcessorSpec getStrangenessTrackingWriterSpec() +{ + auto loggerV = [](std::vector const& v) { + LOG(info) << "StrangenessTracker writer pulled " << v.size() << " strange tracks"; + }; + + auto inpStTrkID = InputSpec{"strangetracks", "STK", "STRTRACKS", 0}; + auto inpClusAtt = InputSpec{"clusupdates", "STK", "CLUSUPDATES", 0}; + + return MakeRootTreeWriterSpec("strangenesstracking-writer", + "o2_strange_tracks.root", + MakeRootTreeWriterSpec::TreeAttributes{"o2sim", "Tree with Strange Tracks"}, + BranchDefinition>{inpStTrkID, "StrangeTracks", loggerV}, + BranchDefinition>{inpClusAtt, "ClusUpdates"})(); +} + +} // namespace strangeness_tracking +} // namespace o2 \ No newline at end of file diff --git a/Detectors/StrangenessTracking/workflow/src/strangeness-tracking-workflow.cxx b/Detectors/StrangenessTracking/workflow/src/strangeness-tracking-workflow.cxx new file mode 100644 index 0000000000000..ddfae84841129 --- /dev/null +++ b/Detectors/StrangenessTracking/workflow/src/strangeness-tracking-workflow.cxx @@ -0,0 +1,66 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN 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 "StrangenessTrackingWorkflow/StrangenessTrackingSpec.h" +#include "StrangenessTrackingWorkflow/StrangenessTrackingWriterSpec.h" + +#include "CommonUtils/ConfigurableParam.h" +#include "StrangenessTracking/StrangenessTrackingConfigParam.h" + +#include "DetectorsRaw/HBFUtilsInitializer.h" +#include "Framework/CallbacksPolicy.h" +#include "Framework/ConfigContext.h" + +#include + +using namespace o2::framework; +using namespace o2::strangeness_tracking; + +void customize(std::vector& policies) +{ + o2::raw::HBFUtilsInitializer::addNewTimeSliceCallback(policies); +} + +void customize(std::vector& workflowOptions) +{ + // option allowing to set parameters + std::vector options{ + {"disable-root-input", o2::framework::VariantType::Bool, false, {"disable root-files input reader"}}, + {"disable-mc", o2::framework::VariantType::Bool, false, {"disable MC"}}, + {"configKeyValues", VariantType::String, "", {"Semicolon separated key=value strings"}}}; + o2::raw::HBFUtilsInitializer::addConfigOption(options); + std::swap(workflowOptions, options); +} + +// ------------------------------------------------------------------ + +#include "Framework/runDataProcessing.h" +#include "Framework/Logger.h" + +WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) +{ + // Update the (declared) parameters if changed from the command line + auto useMC = !configcontext.options().get("disable-mc"); + auto useRootInput = !configcontext.options().get("disable-root-input"); + + o2::conf::ConfigurableParam::updateFromString(configcontext.options().get("configKeyValues")); + + auto wf = o2::strangeness_tracking::getWorkflow(useMC, useRootInput); + wf.emplace_back(getStrangenessTrackingWriterSpec()); + + // configure dpl timer to inject correct firstTFOrbit: start from the 1st orbit of TF containing 1st sampled orbit + o2::raw::HBFUtilsInitializer hbfIni(configcontext, wf); + + // write the configuration used for the reco workflow + o2::conf::ConfigurableParam::writeINI("o2strangeness_tracking_workflow_configuration.ini"); + + return std::move(wf); +}