Project
Loading...
Searching...
No Matches
AODProducerWorkflowSpec.cxx
Go to the documentation of this file.
1// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
2// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3// All rights not expressly granted are reserved.
4//
5// This software is distributed under the terms of the GNU General Public
6// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7//
8// In applying this license CERN does not waive the privileges and immunities
9// granted to it by virtue of its status as an Intergovernmental Organization
10// or submit itself to any jurisdiction.
11
13
41#include "MathUtils/Utils.h"
51#include "Framework/DataTypes.h"
55#include "FT0Base/Geometry.h"
80#include "ZDCBase/Constants.h"
83#include "TOFBase/Utils.h"
84#include "O2Version.h"
85#include "TMath.h"
86#include "MathUtils/Utils.h"
87#include "Math/SMatrix.h"
88#include "TString.h"
89#include <limits>
90#include <map>
91#include <numeric>
92#include <type_traits>
93#include <unordered_map>
94#include <set>
95#include <string>
96#include <vector>
97#include <thread>
98#include "TLorentzVector.h"
99#include "TVector3.h"
100#include "MathUtils/Tsallis.h"
101#include <random>
102#ifdef WITH_OPENMP
103#include <omp.h>
104#endif
105#include <filesystem>
106#include <nlohmann/json.hpp>
107
108using namespace o2::framework;
109using namespace o2::math_utils::detail;
116
117namespace o2::aodproducer
118{
119
120void AODProducerWorkflowDPL::createCTPReadout(const o2::globaltracking::RecoContainer& recoData, std::vector<o2::ctp::CTPDigit>& ctpDigits, ProcessingContext& pc)
121{
122 // Extraxt CTP Config from CCDB
123 const auto ctpcfg = pc.inputs().get<o2::ctp::CTPConfiguration*>("ctpconfig");
124 ctpcfg->printStream(std::cout);
125 // o2::ctp::CTPConfiguration ctpcfg = o2::ctp::CTPRunManager::getConfigFromCCDB(-1, std::to_string(runNumber)); // how to get run
126 // Extract inputs from recoData
127 uint64_t classMaskEMCAL = 0, classMaskTRD = 0, classMaskPHOSCPV = 0;
128 for (const auto& trgclass : ctpcfg->getCTPClasses()) {
129 if (trgclass.cluster->getClusterDetNames().find("EMC") != std::string::npos) {
130 classMaskEMCAL = trgclass.classMask;
131 }
132 if (trgclass.cluster->getClusterDetNames().find("PHS") != std::string::npos) {
133 classMaskPHOSCPV = trgclass.classMask;
134 }
135 if (trgclass.cluster->getClusterDetNames().find("TRD") != std::string::npos) {
136 classMaskTRD = trgclass.classMask;
137 }
138 }
139 LOG(info) << "createCTPReadout: Class Mask EMCAL -> " << classMaskEMCAL;
140 LOG(info) << "createCTPReadout: Class Mask PHOS/CPV -> " << classMaskPHOSCPV;
141 LOG(info) << "createCTPReadout: Class Mask TRD -> " << classMaskTRD;
142
143 // const auto& fddRecPoints = recoData.getFDDRecPoints();
144 // const auto& fv0RecPoints = recoData.getFV0RecPoints();
145 const auto& triggerrecordEMCAL = recoData.getEMCALTriggers();
146 const auto& triggerrecordPHOSCPV = recoData.getPHOSTriggers();
147 const auto& triggerrecordTRD = recoData.getTRDTriggerRecords();
148 // For EMCAL filter remove calibration triggers
149 std::vector<o2::emcal::TriggerRecord> triggerRecordEMCALPhys;
150 for (const auto& trg : triggerrecordEMCAL) {
151 if (trg.getTriggerBits() & o2::trigger::Cal) {
152 continue;
153 }
154 triggerRecordEMCALPhys.push_back(trg);
155 }
156 // const auto& triggerrecordTRD =recoData.getITSTPCTRDTriggers()
157 //
158
159 // Find TVX triggers, only TRD/EMCAL/PHOS/CPV triggers in coincidence will be accepted
160 std::set<uint64_t> bcsMapT0triggers;
161 const auto& ft0RecPoints = recoData.getFT0RecPoints();
162 for (auto& ft0RecPoint : ft0RecPoints) {
163 auto t0triggers = ft0RecPoint.getTrigger();
164 if (t0triggers.getVertex()) {
165 uint64_t globalBC = ft0RecPoint.getInteractionRecord().toLong();
166 bcsMapT0triggers.insert(globalBC);
167 }
168 }
169
170 auto genericCTPDigitizer = [&bcsMapT0triggers, &ctpDigits](auto triggerrecords, uint64_t classmask) -> int {
171 // Strategy:
172 // find detector trigger based on trigger record from readout and add CTPDigit if trigger there
173 int cntwarnings = 0;
174 uint32_t orbitPrev = 0;
175 uint16_t bcPrev = 0;
176 for (auto& trigger : triggerrecords) {
177 auto orbitPrevT = orbitPrev;
178 auto bcPrevT = bcPrev;
179 bcPrev = trigger.getBCData().bc;
180 orbitPrev = trigger.getBCData().orbit;
181 // dedicated for TRD: remove bogus triggers
182 if (orbitPrev < orbitPrevT || bcPrev >= o2::constants::lhc::LHCMaxBunches || (orbitPrev == orbitPrevT && bcPrev < bcPrevT)) {
183 cntwarnings++;
184 // LOGP(warning, "Bogus TRD trigger at bc:{}/orbit:{} (previous was {}/{}), with {} tracklets and {} digits",bcPrev, orbitPrev, bcPrevT, orbitPrevT, trig.getNumberOfTracklets(), trig.getNumberOfDigits());
185 } else {
186 uint64_t globalBC = trigger.getBCData().toLong();
187 auto t0entry = bcsMapT0triggers.find(globalBC);
188 if (t0entry != bcsMapT0triggers.end()) {
189 auto ctpdig = std::find_if(ctpDigits.begin(), ctpDigits.end(), [globalBC](const o2::ctp::CTPDigit& dig) { return static_cast<uint64_t>(dig.intRecord.toLong()) == globalBC; });
190 if (ctpdig != ctpDigits.end()) {
191 // CTP digit existing from other trigger, merge detector class mask
192 ctpdig->CTPClassMask |= std::bitset<64>(classmask);
193 LOG(debug) << "createCTPReadout: Merging " << classmask << " CTP digits with existing digit, CTP mask " << ctpdig->CTPClassMask;
194 } else {
195 // New CTP digit needed
196 LOG(debug) << "createCTPReadout: New CTP digit needed for class " << classmask << std::endl;
197 auto& ctpdigNew = ctpDigits.emplace_back();
198 ctpdigNew.intRecord.setFromLong(globalBC);
199 ctpdigNew.CTPClassMask = classmask;
200 }
201 } else {
202 LOG(warning) << "createCTPReadout: Found " << classmask << " and no MTVX:" << globalBC;
203 }
204 }
205 }
206 return cntwarnings;
207 };
208
209 auto warningsTRD = genericCTPDigitizer(triggerrecordTRD, classMaskTRD);
210 auto warningsEMCAL = genericCTPDigitizer(triggerRecordEMCALPhys, classMaskEMCAL);
211 auto warningsPHOSCPV = genericCTPDigitizer(triggerrecordPHOSCPV, classMaskPHOSCPV);
212
213 LOG(info) << "createCTPReadout:# of TRD bogus triggers:" << warningsTRD;
214 LOG(info) << "createCTPReadout:# of EMCAL bogus triggers:" << warningsEMCAL;
215 LOG(info) << "createCTPReadout:# of PHOS/CPV bogus triggers:" << warningsPHOSCPV;
216}
217
218void AODProducerWorkflowDPL::collectBCs(const o2::globaltracking::RecoContainer& data,
219 const std::vector<o2::InteractionTimeRecord>& mcRecords,
220 std::map<uint64_t, int>& bcsMap)
221{
222 const auto& primVertices = data.getPrimaryVertices();
223 const auto& fddRecPoints = data.getFDDRecPoints();
224 const auto& ft0RecPoints = data.getFT0RecPoints();
225 const auto& fv0RecPoints = data.getFV0RecPoints();
226 const auto& caloEMCCellsTRGR = data.getEMCALTriggers();
227 const auto& caloPHOSCellsTRGR = data.getPHOSTriggers();
228 const auto& cpvTRGR = data.getCPVTriggers();
229 const auto& ctpDigits = data.getCTPDigits();
230 const auto& zdcBCRecData = data.getZDCBCRecData();
231
232 bcsMap[mStartIR.toLong()] = 1; // store the start of TF
233
234 // collecting non-empty BCs and enumerating them
235 for (auto& rec : mcRecords) {
236 uint64_t globalBC = rec.toLong();
237 bcsMap[globalBC] = 1;
238 }
239
240 for (auto& fddRecPoint : fddRecPoints) {
241 uint64_t globalBC = fddRecPoint.getInteractionRecord().toLong();
242 bcsMap[globalBC] = 1;
243 }
244
245 for (auto& ft0RecPoint : ft0RecPoints) {
246 uint64_t globalBC = ft0RecPoint.getInteractionRecord().toLong();
247 bcsMap[globalBC] = 1;
248 }
249
250 for (auto& fv0RecPoint : fv0RecPoints) {
251 uint64_t globalBC = fv0RecPoint.getInteractionRecord().toLong();
252 bcsMap[globalBC] = 1;
253 }
254
255 for (auto& zdcRecData : zdcBCRecData) {
256 uint64_t globalBC = zdcRecData.ir.toLong();
257 bcsMap[globalBC] = 1;
258 }
259
260 for (auto& vertex : primVertices) {
261 auto& timeStamp = vertex.getTimeStamp();
262 double tsTimeStamp = timeStamp.getTimeStamp() * 1E3; // mus to ns
263 uint64_t globalBC = relativeTime_to_GlobalBC(tsTimeStamp);
264 bcsMap[globalBC] = 1;
265 }
266
267 for (auto& emcaltrg : caloEMCCellsTRGR) {
268 uint64_t globalBC = emcaltrg.getBCData().toLong();
269 bcsMap[globalBC] = 1;
270 }
271
272 for (auto& phostrg : caloPHOSCellsTRGR) {
273 uint64_t globalBC = phostrg.getBCData().toLong();
274 bcsMap[globalBC] = 1;
275 }
276
277 for (auto& cpvtrg : cpvTRGR) {
278 uint64_t globalBC = cpvtrg.getBCData().toLong();
279 bcsMap[globalBC] = 1;
280 }
281
282 for (auto& ctpDigit : ctpDigits) {
283 uint64_t globalBC = ctpDigit.intRecord.toLong();
284 bcsMap[globalBC] = 1;
285 }
286
287 int bcID = 0;
288 for (auto& item : bcsMap) {
289 item.second = bcID;
290 bcID++;
291 }
292}
293
294template <typename TracksCursorType, typename TracksCovCursorType>
295void AODProducerWorkflowDPL::addToTracksTable(TracksCursorType& tracksCursor, TracksCovCursorType& tracksCovCursor,
296 const o2::track::TrackParCov& track, int collisionID, aod::track::TrackTypeEnum type)
297{
298 tracksCursor(collisionID,
299 type,
300 truncateFloatFraction(track.getX(), mTrackX),
301 truncateFloatFraction(track.getAlpha(), mTrackAlpha),
302 track.getY(),
303 track.getZ(),
304 truncateFloatFraction(track.getSnp(), mTrackSnp),
305 truncateFloatFraction(track.getTgl(), mTrackTgl),
306 truncateFloatFraction(track.getQ2Pt(), mTrack1Pt));
307 // trackscov
308 float sY = TMath::Sqrt(track.getSigmaY2()), sZ = TMath::Sqrt(track.getSigmaZ2()), sSnp = TMath::Sqrt(track.getSigmaSnp2()),
309 sTgl = TMath::Sqrt(track.getSigmaTgl2()), sQ2Pt = TMath::Sqrt(track.getSigma1Pt2());
310 tracksCovCursor(truncateFloatFraction(sY, mTrackCovDiag),
311 truncateFloatFraction(sZ, mTrackCovDiag),
312 truncateFloatFraction(sSnp, mTrackCovDiag),
313 truncateFloatFraction(sTgl, mTrackCovDiag),
314 truncateFloatFraction(sQ2Pt, mTrackCovDiag),
315 (Char_t)(128. * track.getSigmaZY() / (sZ * sY)),
316 (Char_t)(128. * track.getSigmaSnpY() / (sSnp * sY)),
317 (Char_t)(128. * track.getSigmaSnpZ() / (sSnp * sZ)),
318 (Char_t)(128. * track.getSigmaTglY() / (sTgl * sY)),
319 (Char_t)(128. * track.getSigmaTglZ() / (sTgl * sZ)),
320 (Char_t)(128. * track.getSigmaTglSnp() / (sTgl * sSnp)),
321 (Char_t)(128. * track.getSigma1PtY() / (sQ2Pt * sY)),
322 (Char_t)(128. * track.getSigma1PtZ() / (sQ2Pt * sZ)),
323 (Char_t)(128. * track.getSigma1PtSnp() / (sQ2Pt * sSnp)),
324 (Char_t)(128. * track.getSigma1PtTgl() / (sQ2Pt * sTgl)));
325}
326
327template <typename TracksExtraCursorType>
328void AODProducerWorkflowDPL::addToTracksExtraTable(TracksExtraCursorType& tracksExtraCursor, TrackExtraInfo& extraInfoHolder)
329{
330 // In case of TPC-only tracks, do not truncate the time error since we encapsulate there a special encoding of
331 // the deltaFwd/Bwd times
332 auto trackTimeRes = extraInfoHolder.trackTimeRes;
333 if (!extraInfoHolder.isTPConly) {
334 trackTimeRes = truncateFloatFraction(trackTimeRes, mTrackTimeError);
335 }
336
337 // extra
338 tracksExtraCursor(truncateFloatFraction(extraInfoHolder.tpcInnerParam, mTrack1Pt),
339 extraInfoHolder.flags,
340 extraInfoHolder.itsClusterSizes,
341 extraInfoHolder.tpcNClsFindable,
342 extraInfoHolder.tpcNClsFindableMinusFound,
343 extraInfoHolder.tpcNClsFindableMinusPID,
344 extraInfoHolder.tpcNClsFindableMinusCrossedRows,
345 extraInfoHolder.tpcNClsShared,
346 extraInfoHolder.trdPattern,
347 truncateFloatFraction(extraInfoHolder.itsChi2NCl, mTrackChi2),
348 truncateFloatFraction(extraInfoHolder.tpcChi2NCl, mTrackChi2),
349 truncateFloatFraction(extraInfoHolder.trdChi2, mTrackChi2),
350 truncateFloatFraction(extraInfoHolder.tofChi2, mTrackChi2),
351 truncateFloatFraction(extraInfoHolder.tpcSignal, mTrackSignal),
352 truncateFloatFraction(extraInfoHolder.trdSignal, mTrackSignal),
353 truncateFloatFraction(extraInfoHolder.length, mTrackSignal),
354 truncateFloatFraction(extraInfoHolder.tofExpMom, mTrack1Pt),
355 truncateFloatFraction(extraInfoHolder.trackEtaEMCAL, mTrackPosEMCAL),
356 truncateFloatFraction(extraInfoHolder.trackPhiEMCAL, mTrackPosEMCAL),
357 truncateFloatFraction(extraInfoHolder.trackTime, mTrackTime),
358 trackTimeRes);
359}
360
361template <typename TracksQACursorType>
362void AODProducerWorkflowDPL::addToTracksQATable(TracksQACursorType& tracksQACursor, TrackQA& trackQAInfoHolder)
363{
364 tracksQACursor(
365 trackQAInfoHolder.trackID,
366 truncateFloatFraction(trackQAInfoHolder.tpcTime0, mTPCTime0),
367 truncateFloatFraction(trackQAInfoHolder.tpcdEdxNorm, mTrackSignal),
368 trackQAInfoHolder.tpcdcaR,
369 trackQAInfoHolder.tpcdcaZ,
370 trackQAInfoHolder.tpcClusterByteMask,
371 trackQAInfoHolder.tpcdEdxMax0R,
372 trackQAInfoHolder.tpcdEdxMax1R,
373 trackQAInfoHolder.tpcdEdxMax2R,
374 trackQAInfoHolder.tpcdEdxMax3R,
375 trackQAInfoHolder.tpcdEdxTot0R,
376 trackQAInfoHolder.tpcdEdxTot1R,
377 trackQAInfoHolder.tpcdEdxTot2R,
378 trackQAInfoHolder.tpcdEdxTot3R,
379 trackQAInfoHolder.dRefContY,
380 trackQAInfoHolder.dRefContZ,
381 trackQAInfoHolder.dRefContSnp,
382 trackQAInfoHolder.dRefContTgl,
383 trackQAInfoHolder.dRefContQ2Pt,
384 trackQAInfoHolder.dRefGloY,
385 trackQAInfoHolder.dRefGloZ,
386 trackQAInfoHolder.dRefGloSnp,
387 trackQAInfoHolder.dRefGloTgl,
388 trackQAInfoHolder.dRefGloQ2Pt,
389 trackQAInfoHolder.dTofdX,
390 trackQAInfoHolder.dTofdZ);
391}
392
393template <typename mftTracksCursorType, typename AmbigMFTTracksCursorType>
394void AODProducerWorkflowDPL::addToMFTTracksTable(mftTracksCursorType& mftTracksCursor, AmbigMFTTracksCursorType& ambigMFTTracksCursor,
395 GIndex trackID, const o2::globaltracking::RecoContainer& data, int collisionID,
396 std::uint64_t collisionBC, const std::map<uint64_t, int>& bcsMap)
397{
398 // mft tracks
399 int bcSlice[2] = {-1, -1};
400 const auto& track = data.getMFTTrack(trackID);
401 const auto& rof = data.getMFTTracksROFRecords()[mMFTROFs[trackID.getIndex()]];
402 float trackTime = rof.getBCData().differenceInBC(mStartIR) * o2::constants::lhc::LHCBunchSpacingNS + mMFTROFrameHalfLengthNS + mMFTROFBiasNS;
403 float trackTimeRes = mMFTROFrameHalfLengthNS;
404 bool needBCSlice = collisionID < 0;
405 std::uint64_t bcOfTimeRef;
406 if (needBCSlice) {
407 double error = mTimeMarginTrackTime + trackTimeRes;
408 bcOfTimeRef = fillBCSlice(bcSlice, trackTime - error, trackTime + error, bcsMap);
409 } else {
410 bcOfTimeRef = collisionBC - mStartIR.toLong(); // by default (unambiguous) track time is wrt collision BC
411 }
412 trackTime -= bcOfTimeRef * o2::constants::lhc::LHCBunchSpacingNS;
413
414 // the Cellular Automaton track-finding algorithm flag is stored in first of the 4 bits not used for the cluster size
415 uint64_t mftClusterSizesAndTrackFlags = track.getClusterSizes();
416 mftClusterSizesAndTrackFlags |= (track.isCA()) ? (1ULL << (60)) : 0;
417
418 mftTracksCursor(collisionID,
419 track.getX(),
420 track.getY(),
421 truncateFloatFraction(track.getZ(), mTrackX), // for the forward tracks Z has the same role as X in barrel
422 truncateFloatFraction(track.getPhi(), mTrackAlpha),
423 truncateFloatFraction(track.getTanl(), mTrackTgl),
424 truncateFloatFraction(track.getInvQPt(), mTrack1Pt),
425 mftClusterSizesAndTrackFlags,
426 truncateFloatFraction(track.getTrackChi2(), mTrackChi2),
427 truncateFloatFraction(trackTime, mTrackTime),
428 truncateFloatFraction(trackTimeRes, mTrackTimeError));
429 if (needBCSlice) {
430 ambigMFTTracksCursor(mTableTrMFTID, bcSlice);
431 }
432}
433
434template <typename TracksCursorType, typename TracksCovCursorType, typename TracksExtraCursorType, typename TracksQACursorType, typename AmbigTracksCursorType,
435 typename MFTTracksCursorType, typename MFTTracksCovCursorType, typename AmbigMFTTracksCursorType,
436 typename FwdTracksCursorType, typename FwdTracksCovCursorType, typename AmbigFwdTracksCursorType, typename FwdTrkClsCursorType>
437void AODProducerWorkflowDPL::fillTrackTablesPerCollision(int collisionID,
438 std::uint64_t collisionBC,
439 const o2::dataformats::VtxTrackRef& trackRef,
440 const gsl::span<const GIndex>& GIndices,
442 TracksCursorType& tracksCursor,
443 TracksCovCursorType& tracksCovCursor,
444 TracksExtraCursorType& tracksExtraCursor,
445 TracksQACursorType& tracksQACursor,
446 AmbigTracksCursorType& ambigTracksCursor,
447 MFTTracksCursorType& mftTracksCursor,
448 MFTTracksCovCursorType& mftTracksCovCursor,
449 AmbigMFTTracksCursorType& ambigMFTTracksCursor,
450 FwdTracksCursorType& fwdTracksCursor,
451 FwdTracksCovCursorType& fwdTracksCovCursor,
452 AmbigFwdTracksCursorType& ambigFwdTracksCursor,
453 FwdTrkClsCursorType& fwdTrkClsCursor,
454 const std::map<uint64_t, int>& bcsMap)
455{
456 for (int src = GIndex::NSources; src--;) {
457 if (!GIndex::isTrackSource(src)) {
458 continue;
459 }
460 int start = trackRef.getFirstEntryOfSource(src);
461 int end = start + trackRef.getEntriesOfSource(src);
462 int nToReserve = end - start; // + last index for a given table
463 if (src == GIndex::Source::MFT) {
464 mftTracksCursor.reserve(nToReserve + mftTracksCursor.lastIndex());
466 fwdTracksCursor.reserve(nToReserve + fwdTracksCursor.lastIndex());
467 fwdTracksCovCursor.reserve(nToReserve + fwdTracksCovCursor.lastIndex());
469 mftTracksCovCursor.reserve(nToReserve + mftTracksCovCursor.lastIndex());
470 }
471 } else {
472 tracksCursor.reserve(nToReserve + tracksCursor.lastIndex());
473 tracksCovCursor.reserve(nToReserve + tracksCovCursor.lastIndex());
474 tracksExtraCursor.reserve(nToReserve + tracksExtraCursor.lastIndex());
475 }
476 for (int ti = start; ti < end; ti++) {
477 const auto& trackIndex = GIndices[ti];
478 if (GIndex::includesSource(src, mInputSources)) {
479 if (src == GIndex::Source::MFT) { // MFT tracks are treated separately since they are stored in a different table
480 if (trackIndex.isAmbiguous() && mGIDToTableMFTID.find(trackIndex) != mGIDToTableMFTID.end()) { // was it already stored ?
481 continue;
482 }
483 addToMFTTracksTable(mftTracksCursor, ambigMFTTracksCursor, trackIndex, data, collisionID, collisionBC, bcsMap);
484 mGIDToTableMFTID.emplace(trackIndex, mTableTrMFTID);
485 mTableTrMFTID++;
486 } else if (src == GIndex::Source::MCH || src == GIndex::Source::MFTMCH || src == GIndex::Source::MCHMID) { // FwdTracks tracks are treated separately since they are stored in a different table
487 if (trackIndex.isAmbiguous() && mGIDToTableFwdID.find(trackIndex) != mGIDToTableFwdID.end()) { // was it already stored ?
488 continue;
489 }
490 addToFwdTracksTable(fwdTracksCursor, fwdTracksCovCursor, ambigFwdTracksCursor, mftTracksCovCursor, trackIndex, data, collisionID, collisionBC, bcsMap);
491 mGIDToTableFwdID.emplace(trackIndex, mTableTrFwdID);
492 addClustersToFwdTrkClsTable(data, fwdTrkClsCursor, trackIndex, mTableTrFwdID);
493 mTableTrFwdID++;
494 } else {
495 // barrel track: normal tracks table
496 if (trackIndex.isAmbiguous() && mGIDToTableID.find(trackIndex) != mGIDToTableID.end()) { // was it already stored ?
497 continue;
498 }
499
500 float weight = 0;
501 static std::uniform_real_distribution<> distr(0., 1.);
502 bool writeQAData = o2::math_utils::Tsallis::downsampleTsallisCharged(data.getTrackParam(trackIndex).getPt(), mTrackQCFraction, mSqrtS, weight, distr(mGenerator));
503 auto extraInfoHolder = processBarrelTrack(collisionID, collisionBC, trackIndex, data, bcsMap);
504
505 if (writeQAData) {
506 auto trackQAInfoHolder = processBarrelTrackQA(collisionID, collisionBC, trackIndex, data, bcsMap);
507 if (std::bitset<8>(trackQAInfoHolder.tpcClusterByteMask).count() >= mTrackQCNTrCut) {
508 trackQAInfoHolder.trackID = mTableTrID;
509 // LOGP(info, "orig time0 in bc: {} diffBCRef: {}, ttime: {} -> {}", trackQAInfoHolder.tpcTime0*8, extraInfoHolder.diffBCRef, extraInfoHolder.trackTime, (trackQAInfoHolder.tpcTime0 * 8 - extraInfoHolder.diffBCRef) * o2::constants::lhc::LHCBunchSpacingNS - extraInfoHolder.trackTime);
510 trackQAInfoHolder.tpcTime0 = (trackQAInfoHolder.tpcTime0 * 8 - extraInfoHolder.diffBCRef) * o2::constants::lhc::LHCBunchSpacingNS - extraInfoHolder.trackTime;
511 // difference between TPC track time0 and stored track nominal time in ns instead of TF start
512 addToTracksQATable(tracksQACursor, trackQAInfoHolder);
513 } else {
514 writeQAData = false;
515 }
516 }
517
518 // include specific selection of tpc standalone tracks if thinning is active
519 if (mThinTracks && extraInfoHolder.isTPConly && !writeQAData) { // if trackQA is written then no check has to be done
520 auto trk = data.getTPCTrack(trackIndex);
521 if (trk.getNClusters() >= mTrackQCNCls && trk.getPt() >= mTrackQCPt) {
522 o2::dataformats::DCA dcaInfo{999.f, 999.f, 999.f, 999.f, 999.f};
523 o2::dataformats::VertexBase v = mVtx.getMeanVertex(collisionID < 0 ? 0.f : data.getPrimaryVertex(collisionID).getZ());
524 if (o2::base::Propagator::Instance()->propagateToDCABxByBz(v, trk, 2., mMatCorr, &dcaInfo) && std::abs(dcaInfo.getY()) < mTrackQCDCAxy) {
525 writeQAData = true; // just setting this to not thin the track
526 }
527 }
528 }
529
530 // Skip thinning if not enabled or track is not tpc standalone or assoc. to a V0 or qa'ed
531 if (mThinTracks && src == GIndex::Source::TPC && mGIDUsedBySVtx.find(trackIndex) == mGIDUsedBySVtx.end() && mGIDUsedByStr.find(trackIndex) == mGIDUsedByStr.end() && !writeQAData) {
532 mGIDToTableID.emplace(trackIndex, -1); // pretend skipped tracks are stored; this is safe since they are are not written to disk and -1 indicates to all users to not use this track
533 continue;
534 }
535
536 if (!extraInfoHolder.isTPConly && extraInfoHolder.trackTimeRes < 0.f) { // failed or rejected?
537 LOG(warning) << "Barrel track " << trackIndex << " has no time set, rejection is not expected : time=" << extraInfoHolder.trackTime
538 << " timeErr=" << extraInfoHolder.trackTimeRes << " BCSlice: " << extraInfoHolder.bcSlice[0] << ":" << extraInfoHolder.bcSlice[1];
539 continue;
540 }
541 const auto& trOrig = data.getTrackParam(trackIndex);
542 bool isProp = false;
543 if (mPropTracks && trOrig.getX() < mMaxPropXiu &&
544 mGIDUsedBySVtx.find(trackIndex) == mGIDUsedBySVtx.end() &&
545 mGIDUsedByStr.find(trackIndex) == mGIDUsedByStr.end()) { // Do not propagate track assoc. to V0s and str. tracking
546 auto trackPar(trOrig);
547 isProp = propagateTrackToPV(trackPar, data, collisionID);
548 if (isProp) {
549 addToTracksTable(tracksCursor, tracksCovCursor, trackPar, collisionID, aod::track::Track);
550 }
551 }
552 if (!isProp) {
553 addToTracksTable(tracksCursor, tracksCovCursor, trOrig, collisionID, aod::track::TrackIU);
554 }
555 addToTracksExtraTable(tracksExtraCursor, extraInfoHolder);
556
557 // collecting table indices of barrel tracks for V0s table
558 if (extraInfoHolder.bcSlice[0] >= 0 && collisionID < 0) {
559 ambigTracksCursor(mTableTrID, extraInfoHolder.bcSlice);
560 }
561 mGIDToTableID.emplace(trackIndex, mTableTrID);
562 mTableTrID++;
563 }
564 }
565 }
566 }
567 if (collisionID < 0) {
568 return;
569 }
571 auto sTracks = data.getStrangeTracks();
572 tracksCursor.reserve(mVertexStrLUT[collisionID + 1] + tracksCursor.lastIndex());
573 tracksCovCursor.reserve(mVertexStrLUT[collisionID + 1] + tracksCovCursor.lastIndex());
574 tracksExtraCursor.reserve(mVertexStrLUT[collisionID + 1] + tracksExtraCursor.lastIndex());
575 for (int iS{mVertexStrLUT[collisionID]}; iS < mVertexStrLUT[collisionID + 1]; ++iS) {
576 auto& collStrTrk = mCollisionStrTrk[iS];
577 auto& sTrk = sTracks[collStrTrk.second];
578 TrackExtraInfo extraInfo;
579 extraInfo.itsChi2NCl = sTrk.mTopoChi2; // TODO: this is the total chi2 of adding the ITS clusters, the topology chi2 meaning might change in the future
580 extraInfo.itsClusterSizes = sTrk.getClusterSizes();
581 addToTracksTable(tracksCursor, tracksCovCursor, sTrk.mMother, collisionID, aod::track::StrangeTrack);
582 addToTracksExtraTable(tracksExtraCursor, extraInfo);
583 mStrTrkIndices[collStrTrk.second] = mTableTrID;
584 mTableTrID++;
585 }
586}
587
588void AODProducerWorkflowDPL::fillIndexTablesPerCollision(const o2::dataformats::VtxTrackRef& trackRef, const gsl::span<const GIndex>& GIndices, const o2::globaltracking::RecoContainer& data)
589{
590 const auto& mchmidMatches = data.getMCHMIDMatches();
591
593 int start = trackRef.getFirstEntryOfSource(src);
594 int end = start + trackRef.getEntriesOfSource(src);
595 for (int ti = start; ti < end; ti++) {
596 auto& trackIndex = GIndices[ti];
597 if (GIndex::includesSource(src, mInputSources)) {
598 if (src == GIndex::Source::MFT) {
599 if (trackIndex.isAmbiguous() && mGIDToTableMFTID.find(trackIndex) != mGIDToTableMFTID.end()) {
600 continue;
601 }
602 mGIDToTableMFTID.emplace(trackIndex, mIndexMFTID);
603 mIndexTableMFT[trackIndex.getIndex()] = mIndexMFTID;
604 mIndexMFTID++;
606 if (trackIndex.isAmbiguous() && mGIDToTableFwdID.find(trackIndex) != mGIDToTableFwdID.end()) {
607 continue;
608 }
609 mGIDToTableFwdID.emplace(trackIndex, mIndexFwdID);
610 if (src == GIndex::Source::MCH) {
611 mIndexTableFwd[trackIndex.getIndex()] = mIndexFwdID;
612 } else if (src == GIndex::Source::MCHMID) {
613 const auto& mchmidMatch = mchmidMatches[trackIndex.getIndex()];
614 const auto mchTrackID = mchmidMatch.getMCHRef().getIndex();
615 mIndexTableFwd[mchTrackID] = mIndexFwdID;
616 }
617 mIndexFwdID++;
618 }
619 }
620 }
621 }
622}
623
624template <typename FwdTracksCursorType, typename FwdTracksCovCursorType, typename AmbigFwdTracksCursorType, typename mftTracksCovCursorType>
625void AODProducerWorkflowDPL::addToFwdTracksTable(FwdTracksCursorType& fwdTracksCursor, FwdTracksCovCursorType& fwdTracksCovCursor,
626 AmbigFwdTracksCursorType& ambigFwdTracksCursor, mftTracksCovCursorType& mftTracksCovCursor, GIndex trackID,
627 const o2::globaltracking::RecoContainer& data, int collisionID, std::uint64_t collisionBC,
628 const std::map<uint64_t, int>& bcsMap)
629{
630 const auto& mchTracks = data.getMCHTracks();
631 const auto& midTracks = data.getMIDTracks();
632 const auto& mchmidMatches = data.getMCHMIDMatches();
633 const auto& mchClusters = data.getMCHTrackClusters();
634
635 FwdTrackInfo fwdInfo;
636 FwdTrackCovInfo fwdCovInfo;
637 int bcSlice[2] = {-1, -1};
638
639 // helper lambda for mch bitmap -- common for global and standalone tracks
640 auto getMCHBitMap = [&](int mchTrackID) {
641 if (mchTrackID != -1) { // check matching just in case
642 const auto& mchTrack = mchTracks[mchTrackID];
643 int first = mchTrack.getFirstClusterIdx();
644 int last = mchTrack.getLastClusterIdx();
645 for (int i = first; i <= last; i++) { // check chamberIds of all clusters
646 const auto& cluster = mchClusters[i];
647 int chamberId = cluster.getChamberId();
648 fwdInfo.mchBitMap |= 1 << chamberId;
649 }
650 }
651 };
652
653 auto getMIDBitMapBoards = [&](int midTrackID) {
654 if (midTrackID != -1) { // check matching just in case
655 const auto& midTrack = midTracks[midTrackID];
656 fwdInfo.midBitMap = midTrack.getHitMap();
657 fwdInfo.midBoards = midTrack.getEfficiencyWord();
658 }
659 };
660
661 auto extrapMCHTrack = [&](int mchTrackID) {
662 const auto& track = mchTracks[mchTrackID];
663
664 // mch standalone tracks extrapolated to vertex
665 // compute 3 sets of tracks parameters :
666 // - at vertex
667 // - at DCA
668 // - at the end of the absorber
669 // extrapolate to vertex
670 float vx = 0, vy = 0, vz = 0;
671 if (collisionID >= 0) {
672 const auto& v = data.getPrimaryVertex(collisionID);
673 vx = v.getX();
674 vy = v.getY();
675 vz = v.getZ();
676 }
677
678 o2::mch::TrackParam trackParamAtVertex(track.getZ(), track.getParameters(), track.getCovariances());
679 if (mPropMuons) {
680 double errVtx{0.0}; // FIXME: get errors associated with vertex if available
681 double errVty{0.0};
682 if (!o2::mch::TrackExtrap::extrapToVertex(trackParamAtVertex, vx, vy, vz, errVtx, errVty)) {
683 return false;
684 }
685 }
686
687 // extrapolate to DCA
688 o2::mch::TrackParam trackParamAtDCA(track.getZ(), track.getParameters());
689 if (!o2::mch::TrackExtrap::extrapToVertexWithoutBranson(trackParamAtDCA, vz)) {
690 return false;
691 }
692
693 // extrapolate to the end of the absorber
694 o2::mch::TrackParam trackParamAtRAbs(track.getZ(), track.getParameters());
695 if (!o2::mch::TrackExtrap::extrapToZ(trackParamAtRAbs, -505.)) { // FIXME: replace hardcoded 505
696 return false;
697 }
698
699 double dcaX = trackParamAtDCA.getNonBendingCoor() - vx;
700 double dcaY = trackParamAtDCA.getBendingCoor() - vy;
701 double dca = std::sqrt(dcaX * dcaX + dcaY * dcaY);
702
703 double xAbs = trackParamAtRAbs.getNonBendingCoor();
704 double yAbs = trackParamAtRAbs.getBendingCoor();
705
706 double dpdca = track.getP() * dca;
707 double dchi2 = track.getChi2OverNDF();
708
709 auto fwdmuon = mMatching.MCHtoFwd(trackParamAtVertex);
710
711 fwdInfo.x = fwdmuon.getX();
712 fwdInfo.y = fwdmuon.getY();
713 fwdInfo.z = fwdmuon.getZ();
714 fwdInfo.phi = fwdmuon.getPhi();
715 fwdInfo.tanl = fwdmuon.getTgl();
716 fwdInfo.invqpt = fwdmuon.getInvQPt();
717 fwdInfo.rabs = std::sqrt(xAbs * xAbs + yAbs * yAbs);
718 fwdInfo.chi2 = dchi2;
719 fwdInfo.pdca = dpdca;
720 fwdInfo.nClusters = track.getNClusters();
721
722 fwdCovInfo.sigX = TMath::Sqrt(fwdmuon.getCovariances()(0, 0));
723 fwdCovInfo.sigY = TMath::Sqrt(fwdmuon.getCovariances()(1, 1));
724 fwdCovInfo.sigPhi = TMath::Sqrt(fwdmuon.getCovariances()(2, 2));
725 fwdCovInfo.sigTgl = TMath::Sqrt(fwdmuon.getCovariances()(3, 3));
726 fwdCovInfo.sig1Pt = TMath::Sqrt(fwdmuon.getCovariances()(4, 4));
727 fwdCovInfo.rhoXY = (Char_t)(128. * fwdmuon.getCovariances()(0, 1) / (fwdCovInfo.sigX * fwdCovInfo.sigY));
728 fwdCovInfo.rhoPhiX = (Char_t)(128. * fwdmuon.getCovariances()(0, 2) / (fwdCovInfo.sigPhi * fwdCovInfo.sigX));
729 fwdCovInfo.rhoPhiY = (Char_t)(128. * fwdmuon.getCovariances()(1, 2) / (fwdCovInfo.sigPhi * fwdCovInfo.sigY));
730 fwdCovInfo.rhoTglX = (Char_t)(128. * fwdmuon.getCovariances()(0, 3) / (fwdCovInfo.sigTgl * fwdCovInfo.sigX));
731 fwdCovInfo.rhoTglY = (Char_t)(128. * fwdmuon.getCovariances()(1, 3) / (fwdCovInfo.sigTgl * fwdCovInfo.sigY));
732 fwdCovInfo.rhoTglPhi = (Char_t)(128. * fwdmuon.getCovariances()(2, 3) / (fwdCovInfo.sigTgl * fwdCovInfo.sigPhi));
733 fwdCovInfo.rho1PtX = (Char_t)(128. * fwdmuon.getCovariances()(0, 4) / (fwdCovInfo.sig1Pt * fwdCovInfo.sigX));
734 fwdCovInfo.rho1PtY = (Char_t)(128. * fwdmuon.getCovariances()(1, 4) / (fwdCovInfo.sig1Pt * fwdCovInfo.sigY));
735 fwdCovInfo.rho1PtPhi = (Char_t)(128. * fwdmuon.getCovariances()(2, 4) / (fwdCovInfo.sig1Pt * fwdCovInfo.sigPhi));
736 fwdCovInfo.rho1PtTgl = (Char_t)(128. * fwdmuon.getCovariances()(3, 4) / (fwdCovInfo.sig1Pt * fwdCovInfo.sigTgl));
737
738 return true;
739 };
740
741 if (trackID.getSource() == GIndex::MCH) { // This is an MCH track
742 int mchTrackID = trackID.getIndex();
743 getMCHBitMap(mchTrackID);
744 if (!extrapMCHTrack(mchTrackID)) {
745 LOGF(warn, "Unable to extrapolate MCH track with ID %d! Dummy parameters will be used", mchTrackID);
746 }
747 fwdInfo.trackTypeId = o2::aod::fwdtrack::MCHStandaloneTrack;
748 const auto& rof = data.getMCHTracksROFRecords()[mMCHROFs[mchTrackID]];
749 auto time = rof.getTimeMUS(mStartIR).first;
750 fwdInfo.trackTime = time.getTimeStamp() * 1.e3;
751 fwdInfo.trackTimeRes = time.getTimeStampError() * 1.e3;
752 } else if (trackID.getSource() == GIndex::MCHMID) { // This is an MCH-MID track
753 fwdInfo.trackTypeId = o2::aod::fwdtrack::MuonStandaloneTrack;
754 auto mchmidMatch = mchmidMatches[trackID.getIndex()];
755 auto mchTrackID = mchmidMatch.getMCHRef().getIndex();
756 if (!extrapMCHTrack(mchTrackID)) {
757 LOGF(warn, "Unable to extrapolate MCH track with ID %d! Dummy parameters will be used", mchTrackID);
758 }
759 auto midTrackID = mchmidMatch.getMIDRef().getIndex();
760 fwdInfo.chi2matchmchmid = mchmidMatch.getMatchChi2OverNDF();
761 getMCHBitMap(mchTrackID);
762 getMIDBitMapBoards(midTrackID);
763 auto time = mchmidMatch.getTimeMUS(mStartIR).first;
764 fwdInfo.trackTime = time.getTimeStamp() * 1.e3;
765 fwdInfo.trackTimeRes = time.getTimeStampError() * 1.e3;
766 } else { // This is a GlobalMuonTrack or a GlobalForwardTrack
767 const auto& track = data.getGlobalFwdTrack(trackID);
768 const auto& mftTracks = data.getMFTTracks();
769 const auto& mfttrack = mftTracks[track.getMFTTrackID()];
770 if (!extrapMCHTrack(track.getMCHTrackID())) {
771 LOGF(warn, "Unable to extrapolate MCH track with ID %d! Dummy parameters will be used", track.getMCHTrackID());
772 }
773 fwdInfo.x = track.getX();
774 fwdInfo.y = track.getY();
775 fwdInfo.z = track.getZ();
776 fwdInfo.phi = track.getPhi();
777 fwdInfo.tanl = track.getTanl();
778 fwdInfo.invqpt = track.getInvQPt();
779 fwdInfo.chi2 = track.getTrackChi2();
780 // fwdInfo.nClusters = track.getNumberOfPoints();
781 fwdInfo.chi2matchmchmid = track.getMIDMatchingChi2();
782 fwdInfo.chi2matchmchmft = track.getMFTMCHMatchingChi2();
783 fwdInfo.matchscoremchmft = track.getMFTMCHMatchingScore();
784 fwdInfo.matchmfttrackid = mIndexTableMFT[track.getMFTTrackID()];
785 fwdInfo.matchmchtrackid = mIndexTableFwd[track.getMCHTrackID()];
786 fwdInfo.trackTime = track.getTimeMUS().getTimeStamp() * 1.e3;
787 fwdInfo.trackTimeRes = track.getTimeMUS().getTimeStampError() * 1.e3;
788
789 getMCHBitMap(track.getMCHTrackID());
790 getMIDBitMapBoards(track.getMIDTrackID());
791
792 fwdCovInfo.sigX = TMath::Sqrt(track.getCovariances()(0, 0));
793 fwdCovInfo.sigY = TMath::Sqrt(track.getCovariances()(1, 1));
794 fwdCovInfo.sigPhi = TMath::Sqrt(track.getCovariances()(2, 2));
795 fwdCovInfo.sigTgl = TMath::Sqrt(track.getCovariances()(3, 3));
796 fwdCovInfo.sig1Pt = TMath::Sqrt(track.getCovariances()(4, 4));
797 fwdCovInfo.rhoXY = (Char_t)(128. * track.getCovariances()(0, 1) / (fwdCovInfo.sigX * fwdCovInfo.sigY));
798 fwdCovInfo.rhoPhiX = (Char_t)(128. * track.getCovariances()(0, 2) / (fwdCovInfo.sigPhi * fwdCovInfo.sigX));
799 fwdCovInfo.rhoPhiY = (Char_t)(128. * track.getCovariances()(1, 2) / (fwdCovInfo.sigPhi * fwdCovInfo.sigY));
800 fwdCovInfo.rhoTglX = (Char_t)(128. * track.getCovariances()(0, 3) / (fwdCovInfo.sigTgl * fwdCovInfo.sigX));
801 fwdCovInfo.rhoTglY = (Char_t)(128. * track.getCovariances()(1, 3) / (fwdCovInfo.sigTgl * fwdCovInfo.sigY));
802 fwdCovInfo.rhoTglPhi = (Char_t)(128. * track.getCovariances()(2, 3) / (fwdCovInfo.sigTgl * fwdCovInfo.sigPhi));
803 fwdCovInfo.rho1PtX = (Char_t)(128. * track.getCovariances()(0, 4) / (fwdCovInfo.sig1Pt * fwdCovInfo.sigX));
804 fwdCovInfo.rho1PtY = (Char_t)(128. * track.getCovariances()(1, 4) / (fwdCovInfo.sig1Pt * fwdCovInfo.sigY));
805 fwdCovInfo.rho1PtPhi = (Char_t)(128. * track.getCovariances()(2, 4) / (fwdCovInfo.sig1Pt * fwdCovInfo.sigPhi));
806 fwdCovInfo.rho1PtTgl = (Char_t)(128. * track.getCovariances()(3, 4) / (fwdCovInfo.sig1Pt * fwdCovInfo.sigTgl));
807
808 fwdInfo.trackTypeId = (fwdInfo.chi2matchmchmid >= 0) ? o2::aod::fwdtrack::GlobalMuonTrack : o2::aod::fwdtrack::GlobalForwardTrack;
809
810 float sX = TMath::Sqrt(mfttrack.getSigma2X()), sY = TMath::Sqrt(mfttrack.getSigma2Y()), sPhi = TMath::Sqrt(mfttrack.getSigma2Phi()),
811 sTgl = TMath::Sqrt(mfttrack.getSigma2Tanl()), sQ2Pt = TMath::Sqrt(mfttrack.getSigma2InvQPt());
812
813 mftTracksCovCursor(fwdInfo.matchmfttrackid,
814 truncateFloatFraction(sX, mTrackCovDiag),
815 truncateFloatFraction(sY, mTrackCovDiag),
816 truncateFloatFraction(sPhi, mTrackCovDiag),
817 truncateFloatFraction(sTgl, mTrackCovDiag),
818 truncateFloatFraction(sQ2Pt, mTrackCovDiag),
819 (Char_t)(128. * mfttrack.getCovariances()(0, 1) / (sX * sY)),
820 (Char_t)(128. * mfttrack.getCovariances()(0, 2) / (sPhi * sX)),
821 (Char_t)(128. * mfttrack.getCovariances()(1, 2) / (sPhi * sY)),
822 (Char_t)(128. * mfttrack.getCovariances()(0, 3) / (sTgl * sX)),
823 (Char_t)(128. * mfttrack.getCovariances()(1, 3) / (sTgl * sY)),
824 (Char_t)(128. * mfttrack.getCovariances()(2, 3) / (sTgl * sPhi)),
825 (Char_t)(128. * mfttrack.getCovariances()(0, 4) / (sQ2Pt * sX)),
826 (Char_t)(128. * mfttrack.getCovariances()(1, 4) / (sQ2Pt * sY)),
827 (Char_t)(128. * mfttrack.getCovariances()(2, 4) / (sQ2Pt * sPhi)),
828 (Char_t)(128. * mfttrack.getCovariances()(3, 4) / (sQ2Pt * sTgl)));
829 }
830
831 std::uint64_t bcOfTimeRef;
832 bool needBCSlice = collisionID < 0;
833 if (needBCSlice) { // need to store BC slice
834 float err = mTimeMarginTrackTime + fwdInfo.trackTimeRes;
835 bcOfTimeRef = fillBCSlice(bcSlice, fwdInfo.trackTime - err, fwdInfo.trackTime + err, bcsMap);
836 } else {
837 bcOfTimeRef = collisionBC - mStartIR.toLong(); // by default track time is wrt collision BC (unless no collision assigned)
838 }
839 fwdInfo.trackTime -= bcOfTimeRef * o2::constants::lhc::LHCBunchSpacingNS;
840
841 fwdTracksCursor(collisionID,
842 fwdInfo.trackTypeId,
843 fwdInfo.x,
844 fwdInfo.y,
845 truncateFloatFraction(fwdInfo.z, mTrackX), // for the forward tracks Z has the same role as X in the barrel
846 truncateFloatFraction(fwdInfo.phi, mTrackAlpha),
847 truncateFloatFraction(fwdInfo.tanl, mTrackTgl),
848 truncateFloatFraction(fwdInfo.invqpt, mTrack1Pt),
849 fwdInfo.nClusters,
850 truncateFloatFraction(fwdInfo.pdca, mTrackX),
851 truncateFloatFraction(fwdInfo.rabs, mTrackX),
852 truncateFloatFraction(fwdInfo.chi2, mTrackChi2),
853 truncateFloatFraction(fwdInfo.chi2matchmchmid, mTrackChi2),
854 truncateFloatFraction(fwdInfo.chi2matchmchmft, mTrackChi2),
855 truncateFloatFraction(fwdInfo.matchscoremchmft, mTrackChi2),
856 fwdInfo.matchmfttrackid,
857 fwdInfo.matchmchtrackid,
858 fwdInfo.mchBitMap,
859 fwdInfo.midBitMap,
860 fwdInfo.midBoards,
861 truncateFloatFraction(fwdInfo.trackTime, mTrackTime),
862 truncateFloatFraction(fwdInfo.trackTimeRes, mTrackTimeError));
863
864 fwdTracksCovCursor(truncateFloatFraction(fwdCovInfo.sigX, mTrackCovDiag),
865 truncateFloatFraction(fwdCovInfo.sigY, mTrackCovDiag),
866 truncateFloatFraction(fwdCovInfo.sigPhi, mTrackCovDiag),
867 truncateFloatFraction(fwdCovInfo.sigTgl, mTrackCovDiag),
868 truncateFloatFraction(fwdCovInfo.sig1Pt, mTrackCovDiag),
869 fwdCovInfo.rhoXY,
870 fwdCovInfo.rhoPhiX,
871 fwdCovInfo.rhoPhiY,
872 fwdCovInfo.rhoTglX,
873 fwdCovInfo.rhoTglY,
874 fwdCovInfo.rhoTglPhi,
875 fwdCovInfo.rho1PtX,
876 fwdCovInfo.rho1PtY,
877 fwdCovInfo.rho1PtPhi,
878 fwdCovInfo.rho1PtTgl);
879
880 if (needBCSlice) {
881 ambigFwdTracksCursor(mTableTrFwdID, bcSlice);
882 }
883}
884
885//------------------------------------------------------------------
886void AODProducerWorkflowDPL::updateMCHeader(MCCollisionCursor& collisionCursor,
887 XSectionCursor& xSectionCursor,
888 PdfInfoCursor& pdfInfoCursor,
889 HeavyIonCursor& heavyIonCursor,
890 const MCEventHeader& header,
891 int collisionID,
892 int bcID,
893 float time,
894 short generatorID,
895 int sourceID)
896{
901
902 auto genID = updateMCCollisions(collisionCursor,
903 bcID,
904 time,
905 header,
906 generatorID,
907 sourceID,
908 mCollisionPosition);
909 mXSectionUpdate = (updateHepMCXSection(xSectionCursor, //
910 collisionID, //
911 genID, //
912 header, //
913 mXSectionUpdate)
914 ? HepMCUpdate::always
915 : HepMCUpdate::never);
916 mPdfInfoUpdate = (updateHepMCPdfInfo(pdfInfoCursor, //
917 collisionID, //
918 genID, //
919 header, //
920 mPdfInfoUpdate)
921 ? HepMCUpdate::always
922 : HepMCUpdate::never);
923 mHeavyIonUpdate = (updateHepMCHeavyIon(heavyIonCursor, //
924 collisionID, //
925 genID, //
926 header,
927 mHeavyIonUpdate)
928 ? HepMCUpdate::always
929 : HepMCUpdate::never);
930}
931
932void dimensionMCKeepStore(std::vector<std::vector<std::unordered_map<int, int>>>& store, int Nsources, int NEvents)
933{
934 store.resize(Nsources);
935 for (int s = 0; s < Nsources; ++s) {
936 store[s].resize(NEvents);
937 }
938}
939
940void clearMCKeepStore(std::vector<std::vector<std::unordered_map<int, int>>>& store)
941{
942 for (auto s = 0U; s < store.size(); ++s) {
943 for (auto e = 0U; e < store[s].size(); ++e) {
944 store[s][e].clear();
945 }
946 }
947}
948
949// helper function to add a particle/track to the MC keep store
950void keepMCParticle(std::vector<std::vector<std::unordered_map<int, int>>>& store, int source, int event, int track, int value = 1)
951{
952 if (track < 0) {
953 LOG(warn) << "trackID is smaller than 0. Neglecting";
954 return;
955 }
956 store[source][event][track] = value;
957}
958
959void AODProducerWorkflowDPL::fillMCParticlesTable(o2::steer::MCKinematicsReader& mcReader,
960 MCParticlesCursor& mcParticlesCursor,
961 const gsl::span<const o2::dataformats::VtxTrackRef>& primVer2TRefs,
962 const gsl::span<const GIndex>& GIndices,
964 const std::vector<std::vector<int>>& mcColToEvSrc)
965{
966 int NSources = 0;
967 int NEvents = 0;
968 for (auto& p : mcColToEvSrc) {
969 NSources = std::max(p[1], NSources);
970 NEvents = std::max(p[2], NEvents);
971 }
972 NSources++; // 0 - indexed
973 NEvents++;
974 LOG(info) << " number of events " << NEvents;
975 LOG(info) << " number of sources " << NSources;
976 dimensionMCKeepStore(mToStore, NSources, NEvents);
977
978 std::vector<int> particleIDsToKeep;
979
980 auto markMCTrackForSrc = [&](std::array<GID, GID::NSources>& contributorsGID, uint8_t src) {
981 auto mcLabel = data.getTrackMCLabel(contributorsGID[src]);
982 if (!mcLabel.isValid()) {
983 return;
984 }
985 keepMCParticle(mToStore, mcLabel.getSourceID(), mcLabel.getEventID(), mcLabel.getTrackID());
986 };
987
988 // mark reconstructed MC particles to store them into the table
989 for (auto& trackRef : primVer2TRefs) {
990 for (int src = GIndex::NSources; src--;) {
991 int start = trackRef.getFirstEntryOfSource(src);
992 int end = start + trackRef.getEntriesOfSource(src);
993 for (int ti = start; ti < end; ti++) {
994 auto& trackIndex = GIndices[ti];
995 if (GIndex::includesSource(src, mInputSources)) {
996 auto mcTruth = data.getTrackMCLabel(trackIndex);
997 if (!mcTruth.isValid()) {
998 continue;
999 }
1000 keepMCParticle(mToStore, mcTruth.getSourceID(), mcTruth.getEventID(), mcTruth.getTrackID());
1001 // treating contributors of global tracks
1002 auto contributorsGID = data.getSingleDetectorRefs(trackIndex);
1003 if (contributorsGID[GIndex::Source::TPC].isIndexSet()) {
1004 markMCTrackForSrc(contributorsGID, GIndex::Source::TPC);
1005 }
1006 if (contributorsGID[GIndex::Source::ITS].isIndexSet()) {
1007 markMCTrackForSrc(contributorsGID, GIndex::Source::ITS);
1008 }
1009 if (contributorsGID[GIndex::Source::TOF].isIndexSet()) {
1010 const auto& labelsTOF = data.getTOFClustersMCLabels()->getLabels(contributorsGID[GIndex::Source::TOF]);
1011 for (auto& mcLabel : labelsTOF) {
1012 if (!mcLabel.isValid()) {
1013 continue;
1014 }
1015 keepMCParticle(mToStore, mcLabel.getSourceID(), mcLabel.getEventID(), mcLabel.getTrackID());
1016 }
1017 }
1018 }
1019 }
1020 }
1021 }
1022 // mark calorimeter signals as reconstructed particles
1023 if (mInputSources[GIndex::EMC]) {
1024 auto& mcCaloEMCCellLabels = data.getEMCALCellsMCLabels()->getTruthArray();
1025 for (auto& mcTruth : mcCaloEMCCellLabels) {
1026 if (!mcTruth.isValid()) {
1027 continue;
1028 }
1029 keepMCParticle(mToStore, mcTruth.getSourceID(), mcTruth.getEventID(), mcTruth.getTrackID());
1030 }
1031 }
1032 if (mInputSources[GIndex::PHS]) {
1033 auto& mcCaloPHOSCellLabels = data.getPHOSCellsMCLabels()->getTruthArray();
1034 for (auto& mcTruth : mcCaloPHOSCellLabels) {
1035 if (!mcTruth.isValid()) {
1036 continue;
1037 }
1038 keepMCParticle(mToStore, mcTruth.getSourceID(), mcTruth.getEventID(), mcTruth.getTrackID());
1039 }
1040 }
1041 using namespace aodmchelpers;
1042 using MCTrackNavigator = o2::mcutils::MCTrackNavigator;
1043
1044 size_t offset = 0;
1045 for (auto& colInfo : mcColToEvSrc) { // loop over "<eventID, sourceID> <-> combined MC col. ID" key pairs
1046 int event = colInfo[2];
1047 int source = colInfo[1];
1048 int mcColId = colInfo[0];
1049 std::vector<MCTrack> const& mcParticles = mcReader.getTracks(source, event);
1050 LOG(debug) << "Event=" << event << " source=" << source << " collision=" << mcColId;
1051 auto& preselect = mToStore[source][event];
1052
1053 offset = updateParticles(mcParticlesCursor,
1054 mcColId,
1055 mcParticles,
1056 preselect,
1057 offset,
1058 mRecoOnly,
1059 source == 0, // background
1060 mMcParticleW,
1061 mMcParticleMom,
1062 mMcParticlePos);
1063
1065 }
1066}
1067
1068template <typename MCTrackLabelCursorType, typename MCMFTTrackLabelCursorType, typename MCFwdTrackLabelCursorType>
1069void AODProducerWorkflowDPL::fillMCTrackLabelsTable(MCTrackLabelCursorType& mcTrackLabelCursor,
1070 MCMFTTrackLabelCursorType& mcMFTTrackLabelCursor,
1071 MCFwdTrackLabelCursorType& mcFwdTrackLabelCursor,
1072 const o2::dataformats::VtxTrackRef& trackRef,
1073 const gsl::span<const GIndex>& primVerGIs,
1075 int vertexId)
1076{
1077 // labelMask (temporary) usage:
1078 // bit 13 -- ITS/TPC with ITS label (track of AB tracklet) different from TPC
1079 // bit 14 -- isNoise() == true
1080 // bit 15 -- isFake() == true (defined by the fakeness of the top level global track, i.e. if TOF is present, fake means that the track of the TPC label does not contribute to TOF cluster)
1081 // labelID = -1 -- label is not set
1082
1083 for (int src = GIndex::NSources; src--;) {
1084 int start = trackRef.getFirstEntryOfSource(src);
1085 int end = start + trackRef.getEntriesOfSource(src);
1086 mcMFTTrackLabelCursor.reserve(end - start + mcMFTTrackLabelCursor.lastIndex());
1087 mcFwdTrackLabelCursor.reserve(end - start + mcFwdTrackLabelCursor.lastIndex());
1088 mcTrackLabelCursor.reserve(end - start + mcTrackLabelCursor.lastIndex());
1089 for (int ti = start; ti < end; ti++) {
1090 const auto trackIndex = primVerGIs[ti];
1091
1092 // check if the label was already stored (or the track was rejected for some reason in the fillTrackTablesPerCollision)
1093 auto needToStore = [trackIndex](std::unordered_map<GIndex, int>& mp) {
1094 auto entry = mp.find(trackIndex);
1095 if (entry == mp.end() || entry->second == -1) {
1096 return false;
1097 }
1098 entry->second = -1;
1099 return true;
1100 };
1101
1102 if (GIndex::includesSource(src, mInputSources)) {
1103 auto mcTruth = data.getTrackMCLabel(trackIndex);
1104 MCLabels labelHolder{};
1105 if ((src == GIndex::Source::MFT) || (src == GIndex::Source::MFTMCH) || (src == GIndex::Source::MCH) || (src == GIndex::Source::MCHMID)) { // treating mft and fwd labels separately
1106 if (!needToStore(src == GIndex::Source::MFT ? mGIDToTableMFTID : mGIDToTableFwdID)) {
1107 continue;
1108 }
1109 if (mcTruth.isValid()) { // if not set, -1 will be stored
1110 labelHolder.labelID = (mToStore[mcTruth.getSourceID()][mcTruth.getEventID()])[mcTruth.getTrackID()];
1111 }
1112 if (mcTruth.isFake()) {
1113 labelHolder.fwdLabelMask |= (0x1 << 7);
1114 }
1115 if (mcTruth.isNoise()) {
1116 labelHolder.fwdLabelMask |= (0x1 << 6);
1117 }
1118 if (src == GIndex::Source::MFT) {
1119 mcMFTTrackLabelCursor(labelHolder.labelID,
1120 labelHolder.fwdLabelMask);
1121 } else {
1122 mcFwdTrackLabelCursor(labelHolder.labelID,
1123 labelHolder.fwdLabelMask);
1124 }
1125 } else {
1126 if (!needToStore(mGIDToTableID)) {
1127 continue;
1128 }
1129 if (mcTruth.isValid()) { // if not set, -1 will be stored
1130 labelHolder.labelID = (mToStore[mcTruth.getSourceID()][mcTruth.getEventID()])[mcTruth.getTrackID()]; // defined by TPC if it contributes, otherwise: by ITS
1131 if (mcTruth.isFake()) {
1132 labelHolder.labelMask |= (0x1 << 15);
1133 }
1134 if (trackIndex.includesDet(DetID::TPC) && trackIndex.getSource() != GIndex::Source::TPC) { // this is global track
1135 auto contributorsGID = data.getSingleDetectorRefs(trackIndex);
1136 if (contributorsGID[GIndex::Source::ITSTPC].isIndexSet()) { // there is a match to ITS tracks or ITSAB tracklet!
1137 if (data.getTrackMCLabel(contributorsGID[GIndex::Source::ITSTPC]).isFake()) {
1138 labelHolder.labelMask |= (0x1 << 13);
1139 }
1140 }
1141 }
1142 } else if (mcTruth.isNoise()) {
1143 labelHolder.labelMask |= (0x1 << 14);
1144 }
1145 mcTrackLabelCursor(labelHolder.labelID, labelHolder.labelMask);
1146 }
1147 }
1148 }
1149 }
1150
1151 // filling the tables with the strangeness tracking labels
1152 auto sTrackLabels = data.getStrangeTracksMCLabels();
1153 // check if vertexId and vertexId + 1 maps into mVertexStrLUT
1154 if (!(vertexId < 0 || vertexId >= mVertexStrLUT.size() - 1)) {
1155 mcTrackLabelCursor.reserve(mVertexStrLUT[vertexId + 1] + mcTrackLabelCursor.lastIndex());
1156 for (int iS{mVertexStrLUT[vertexId]}; iS < mVertexStrLUT[vertexId + 1]; ++iS) {
1157 auto& collStrTrk = mCollisionStrTrk[iS];
1158 auto& label = sTrackLabels[collStrTrk.second];
1159 MCLabels labelHolder;
1160 labelHolder.labelID = label.isValid() ? (mToStore[label.getSourceID()][label.getEventID()])[label.getTrackID()] : -1;
1161 labelHolder.labelMask = (label.isFake() << 15) | (label.isNoise() << 14);
1162 mcTrackLabelCursor(labelHolder.labelID, labelHolder.labelMask);
1163 }
1164 }
1165}
1166
1167template <typename V0CursorType, typename CascadeCursorType, typename Decay3BodyCursorType>
1168void AODProducerWorkflowDPL::fillSecondaryVertices(const o2::globaltracking::RecoContainer& recoData, V0CursorType& v0Cursor, CascadeCursorType& cascadeCursor, Decay3BodyCursorType& decay3BodyCursor)
1169{
1170
1171 auto v0s = recoData.getV0sIdx();
1172 auto cascades = recoData.getCascadesIdx();
1173 auto decays3Body = recoData.getDecays3BodyIdx();
1174
1175 v0Cursor.reserve(v0s.size());
1176 // filling v0s table
1177 for (size_t iv0 = 0; iv0 < v0s.size(); iv0++) {
1178 const auto& v0 = v0s[iv0];
1179 auto trPosID = v0.getProngID(0);
1180 auto trNegID = v0.getProngID(1);
1181 uint8_t v0flags = v0.getBits();
1182 int posTableIdx = -1, negTableIdx = -1, collID = -1;
1183 auto item = mGIDToTableID.find(trPosID);
1184 if (item != mGIDToTableID.end()) {
1185 posTableIdx = item->second;
1186 } else {
1187 LOG(warn) << "Could not find a positive track index for prong ID " << trPosID;
1188 }
1189 item = mGIDToTableID.find(trNegID);
1190 if (item != mGIDToTableID.end()) {
1191 negTableIdx = item->second;
1192 } else {
1193 LOG(warn) << "Could not find a negative track index for prong ID " << trNegID;
1194 }
1195 auto itemV = mVtxToTableCollID.find(v0.getVertexID());
1196 if (itemV == mVtxToTableCollID.end()) {
1197 LOG(warn) << "Could not find V0 collisionID for the vertex ID " << v0.getVertexID();
1198 } else {
1199 collID = itemV->second;
1200 }
1201 if (posTableIdx != -1 and negTableIdx != -1 and collID != -1) {
1202 v0Cursor(collID, posTableIdx, negTableIdx, v0flags);
1203 mV0ToTableID[int(iv0)] = mTableV0ID++;
1204 }
1205 }
1206
1207 // filling cascades table
1208 cascadeCursor.reserve(cascades.size());
1209 for (auto& cascade : cascades) {
1210 auto itemV0 = mV0ToTableID.find(cascade.getV0ID());
1211 if (itemV0 == mV0ToTableID.end()) {
1212 continue;
1213 }
1214 int v0tableID = itemV0->second, bachTableIdx = -1, collID = -1;
1215 auto bachelorID = cascade.getBachelorID();
1216 auto item = mGIDToTableID.find(bachelorID);
1217 if (item != mGIDToTableID.end()) {
1218 bachTableIdx = item->second;
1219 } else {
1220 LOG(warn) << "Could not find a bachelor track index";
1221 continue;
1222 }
1223 auto itemV = mVtxToTableCollID.find(cascade.getVertexID());
1224 if (itemV != mVtxToTableCollID.end()) {
1225 collID = itemV->second;
1226 } else {
1227 LOG(warn) << "Could not find cascade collisionID for the vertex ID " << cascade.getVertexID();
1228 continue;
1229 }
1230 cascadeCursor(collID, v0tableID, bachTableIdx);
1231 }
1232
1233 // filling 3 body decays table
1234 decay3BodyCursor.reserve(decays3Body.size());
1235 for (size_t i3Body = 0; i3Body < decays3Body.size(); i3Body++) {
1236 const auto& decay3Body = decays3Body[i3Body];
1237 GIndex trIDs[3]{
1238 decay3Body.getProngID(0),
1239 decay3Body.getProngID(1),
1240 decay3Body.getProngID(2)};
1241 int tableIdx[3]{-1, -1, -1}, collID = -1;
1242 bool missing{false};
1243 for (int i{0}; i < 3; ++i) {
1244 auto item = mGIDToTableID.find(trIDs[i]);
1245 if (item != mGIDToTableID.end()) {
1246 tableIdx[i] = item->second;
1247 } else {
1248 LOG(warn) << fmt::format("Could not find a track index for prong ID {}", (int)trIDs[i]);
1249 missing = true;
1250 }
1251 }
1252 auto itemV = mVtxToTableCollID.find(decay3Body.getVertexID());
1253 if (itemV == mVtxToTableCollID.end()) {
1254 LOG(warn) << "Could not find 3 body collisionID for the vertex ID " << decay3Body.getVertexID();
1255 missing = true;
1256 } else {
1257 collID = itemV->second;
1258 }
1259 if (missing) {
1260 continue;
1261 }
1262 decay3BodyCursor(collID, tableIdx[0], tableIdx[1], tableIdx[2]);
1263 }
1264}
1265
1266template <typename FwdTrkClsCursorType>
1267void AODProducerWorkflowDPL::addClustersToFwdTrkClsTable(const o2::globaltracking::RecoContainer& recoData, FwdTrkClsCursorType& fwdTrkClsCursor, GIndex trackID, int fwdTrackId)
1268{
1269 const auto& mchTracks = recoData.getMCHTracks();
1270 const auto& mchmidMatches = recoData.getMCHMIDMatches();
1271 const auto& mchClusters = recoData.getMCHTrackClusters();
1272
1273 int mchTrackID = -1;
1274 if (trackID.getSource() == GIndex::MCH) { // This is an MCH track
1275 mchTrackID = trackID.getIndex();
1276 } else if (trackID.getSource() == GIndex::MCHMID) { // This is an MCH-MID track
1277 auto mchmidMatch = mchmidMatches[trackID.getIndex()];
1278 mchTrackID = mchmidMatch.getMCHRef().getIndex();
1279 } // Others are Global Forward Tracks, their clusters will be or were added with the corresponding MCH track
1280
1281 if (mchTrackID > -1 && mchTrackID < mchTracks.size()) {
1282 const auto& mchTrack = mchTracks[mchTrackID];
1283 fwdTrkClsCursor.reserve(mchTrack.getNClusters() + fwdTrkClsCursor.lastIndex());
1284 int first = mchTrack.getFirstClusterIdx();
1285 int last = mchTrack.getLastClusterIdx();
1286 for (int i = first; i <= last; i++) {
1287 const auto& cluster = mchClusters[i];
1288 fwdTrkClsCursor(fwdTrackId,
1289 truncateFloatFraction(cluster.x, mMuonCl),
1290 truncateFloatFraction(cluster.y, mMuonCl),
1291 truncateFloatFraction(cluster.z, mMuonCl),
1292 (((cluster.ey < 5.) & 0x1) << 12) | (((cluster.ex < 5.) & 0x1) << 11) | cluster.getDEId());
1293 }
1294 }
1295}
1296
1297template <typename HMPCursorType>
1298void AODProducerWorkflowDPL::fillHMPID(const o2::globaltracking::RecoContainer& recoData, HMPCursorType& hmpCursor)
1299{
1300 auto hmpMatches = recoData.getHMPMatches();
1301
1302 hmpCursor.reserve(hmpMatches.size());
1303
1304 // filling HMPs table
1305 for (size_t iHmp = 0; iHmp < hmpMatches.size(); iHmp++) {
1306
1307 const auto& match = hmpMatches[iHmp];
1308
1309 float xTrk, yTrk, theta, phi;
1310 float xMip, yMip;
1311 int charge, nph;
1312
1313 match.getHMPIDtrk(xTrk, yTrk, theta, phi);
1314 match.getHMPIDmip(xMip, yMip, charge, nph);
1315
1316 auto photChargeVec = match.getPhotCharge();
1317
1318 float photChargeVec2[10]; // = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
1319
1320 for (Int_t i = 0; i < 10; i++) {
1321 photChargeVec2[i] = photChargeVec[i];
1322 }
1323 auto tref = mGIDToTableID.find(match.getTrackRef());
1324 if (tref != mGIDToTableID.end()) {
1325 hmpCursor(tref->second, match.getHMPsignal(), xTrk, yTrk, xMip, yMip, nph, charge, match.getIdxHMPClus(), match.getHmpMom(), photChargeVec2);
1326 } else {
1327 LOG(error) << "Could not find AOD track table entry for HMP-matched track " << match.getTrackRef().asString();
1328 }
1329 }
1330}
1331
1332void AODProducerWorkflowDPL::prepareStrangenessTracking(const o2::globaltracking::RecoContainer& recoData)
1333{
1334 auto v0s = recoData.getV0sIdx();
1335 auto cascades = recoData.getCascadesIdx();
1336 auto decays3Body = recoData.getDecays3BodyIdx();
1337
1338 int sTrkID = 0;
1339 mCollisionStrTrk.clear();
1340 mCollisionStrTrk.reserve(recoData.getStrangeTracks().size());
1341 mVertexStrLUT.clear();
1342 mVertexStrLUT.resize(recoData.getPrimaryVertices().size() + 1, 0);
1343 for (auto& sTrk : recoData.getStrangeTracks()) {
1344 auto ITSIndex = GIndex{sTrk.mITSRef, GIndex::ITS};
1345 int vtxId{0};
1346 if (sTrk.mPartType == dataformats::kStrkV0) {
1347 vtxId = v0s[sTrk.mDecayRef].getVertexID();
1348 } else if (sTrk.mPartType == dataformats::kStrkCascade) {
1349 vtxId = cascades[sTrk.mDecayRef].getVertexID();
1350 } else {
1351 vtxId = decays3Body[sTrk.mDecayRef].getVertexID();
1352 }
1353 mCollisionStrTrk.emplace_back(vtxId, sTrkID++);
1354 mVertexStrLUT[vtxId]++;
1355 }
1356 std::exclusive_scan(mVertexStrLUT.begin(), mVertexStrLUT.end(), mVertexStrLUT.begin(), 0);
1357
1358 // sort by collision ID
1359 std::sort(mCollisionStrTrk.begin(), mCollisionStrTrk.end(), [](const auto& a, const auto& b) { return a.first < b.first; });
1360 mStrTrkIndices.clear();
1361 mStrTrkIndices.resize(mCollisionStrTrk.size(), -1);
1362}
1363
1364template <typename V0C, typename CC, typename D3BC>
1365void AODProducerWorkflowDPL::fillStrangenessTrackingTables(const o2::globaltracking::RecoContainer& recoData, V0C& v0Curs, CC& cascCurs, D3BC& d3BodyCurs)
1366{
1367 int itsTableIdx = -1;
1368 int sTrkID = 0;
1369 int nV0 = 0;
1370 int nCasc = 0;
1371 int nD3Body = 0;
1372
1373 for (const auto& sTrk : recoData.getStrangeTracks()) {
1374 if (sTrk.mPartType == dataformats::kStrkV0) {
1375 nV0++;
1376 } else if (sTrk.mPartType == dataformats::kStrkCascade) {
1377 nCasc++;
1378 } else {
1379 nD3Body++;
1380 }
1381 }
1382
1383 v0Curs.reserve(nV0);
1384 cascCurs.reserve(nCasc);
1385 d3BodyCurs.reserve(nD3Body);
1386
1387 for (const auto& sTrk : recoData.getStrangeTracks()) {
1388 auto ITSIndex = GIndex{sTrk.mITSRef, GIndex::ITS};
1389 auto item = mGIDToTableID.find(ITSIndex);
1390 if (item != mGIDToTableID.end()) {
1391 itsTableIdx = item->second;
1392 } else {
1393 LOG(warn) << "Could not find a ITS strange track index " << ITSIndex;
1394 continue;
1395 }
1396 if (sTrk.mPartType == dataformats::kStrkV0) {
1397 v0Curs(mStrTrkIndices[sTrkID++],
1398 itsTableIdx,
1399 sTrk.mDecayRef,
1400 sTrk.mDecayVtx[0],
1401 sTrk.mDecayVtx[1],
1402 sTrk.mDecayVtx[2],
1403 sTrk.mMasses[0],
1404 sTrk.mMasses[1],
1405 sTrk.mMatchChi2,
1406 sTrk.mTopoChi2,
1407 sTrk.getAverageClusterSize());
1408 } else if (sTrk.mPartType == dataformats::kStrkCascade) {
1409 cascCurs(mStrTrkIndices[sTrkID++],
1410 itsTableIdx,
1411 sTrk.mDecayRef,
1412 sTrk.mDecayVtx[0],
1413 sTrk.mDecayVtx[1],
1414 sTrk.mDecayVtx[2],
1415 sTrk.mMasses[0],
1416 sTrk.mMasses[1],
1417 sTrk.mMatchChi2,
1418 sTrk.mTopoChi2,
1419 sTrk.getAverageClusterSize());
1420 } else {
1421 d3BodyCurs(mStrTrkIndices[sTrkID++],
1422 itsTableIdx,
1423 sTrk.mDecayRef,
1424 sTrk.mDecayVtx[0],
1425 sTrk.mDecayVtx[1],
1426 sTrk.mDecayVtx[2],
1427 sTrk.mMasses[0],
1428 sTrk.mMasses[1],
1429 sTrk.mMatchChi2,
1430 sTrk.mTopoChi2,
1431 sTrk.getAverageClusterSize());
1432 }
1433 }
1434}
1435
1436void AODProducerWorkflowDPL::countTPCClusters(const o2::globaltracking::RecoContainer& data)
1437{
1438 const auto& tpcTracks = data.getTPCTracks();
1439 const auto& tpcClusRefs = data.getTPCTracksClusterRefs();
1440 const auto& tpcClusShMap = data.clusterShMapTPC;
1441 const auto& tpcClusAcc = data.getTPCClusters();
1442 constexpr int maxRows = 152;
1443 constexpr int neighbour = 2;
1444 int ntr = tpcTracks.size();
1445 mTPCCounters.clear();
1446 mTPCCounters.resize(ntr);
1447#ifdef WITH_OPENMP
1448 int ngroup = std::min(50, std::max(1, ntr / mNThreads));
1449#pragma omp parallel for schedule(dynamic, ngroup) num_threads(mNThreads)
1450#endif
1451 for (int itr = 0; itr < ntr; itr++) {
1452 std::array<bool, maxRows> clMap{}, shMap{};
1453 uint8_t sectorIndex, rowIndex;
1454 uint32_t clusterIndex;
1455 auto& counters = mTPCCounters[itr];
1456 const auto& track = tpcTracks[itr];
1457 for (int i = 0; i < track.getNClusterReferences(); i++) {
1458 o2::tpc::TrackTPC::getClusterReference(tpcClusRefs, i, sectorIndex, rowIndex, clusterIndex, track.getClusterRef());
1459 unsigned int absoluteIndex = tpcClusAcc.clusterOffset[sectorIndex][rowIndex] + clusterIndex;
1460 clMap[rowIndex] = true;
1461 if (tpcClusShMap[absoluteIndex] & o2::gpu::GPUTPCGMMergedTrackHit::flagShared) {
1462 if (!shMap[rowIndex]) {
1463 counters.shared++;
1464 }
1465 shMap[rowIndex] = true;
1466 }
1467 }
1468 int last = -1;
1469 for (int i = 0; i < maxRows; i++) {
1470 if (clMap[i]) {
1471 counters.crossed++;
1472 counters.found++;
1473 last = i;
1474 } else if ((i - last) <= neighbour) {
1475 counters.crossed++;
1476 } else {
1477 int lim = std::min(i + 1 + neighbour, maxRows);
1478 for (int j = i + 1; j < lim; j++) {
1479 if (clMap[j]) {
1480 counters.crossed++;
1481 break;
1482 }
1483 }
1484 }
1485 }
1486 }
1487}
1488
1489uint8_t AODProducerWorkflowDPL::getTRDPattern(const o2::trd::TrackTRD& track)
1490{
1491 uint8_t pattern = 0;
1492 for (int il = o2::trd::TrackTRD::EGPUTRDTrack::kNLayers - 1; il >= 0; il--) {
1493 if (track.getTrackletIndex(il) != -1) {
1494 pattern |= 0x1 << il;
1495 }
1496 }
1497 if (track.getHasNeighbor()) {
1498 pattern |= 0x1 << 6;
1499 }
1500 if (track.getHasPadrowCrossing()) {
1501 pattern |= 0x1 << 7;
1502 }
1503 return pattern;
1504}
1505
1506template <typename TCaloHandler, typename TCaloCursor, typename TCaloTRGCursor, typename TMCCaloLabelCursor>
1507void AODProducerWorkflowDPL::addToCaloTable(TCaloHandler& caloHandler, TCaloCursor& caloCellCursor, TCaloTRGCursor& caloTRGCursor,
1508 TMCCaloLabelCursor& mcCaloCellLabelCursor, int eventID, int bcID, int8_t caloType)
1509{
1510 auto inputEvent = caloHandler.buildEvent(eventID);
1511 auto cellsInEvent = inputEvent.mCells; // get cells belonging to current event
1512 auto cellMClabels = inputEvent.mMCCellLabels; // get MC labels belonging to current event (only implemented for EMCal currently!)
1513 caloCellCursor.reserve(cellsInEvent.size() + caloCellCursor.lastIndex());
1514 caloTRGCursor.reserve(cellsInEvent.size() + caloTRGCursor.lastIndex());
1515 if (mUseMC) {
1516 mcCaloCellLabelCursor.reserve(cellsInEvent.size() + mcCaloCellLabelCursor.lastIndex());
1517 }
1518 for (auto iCell = 0U; iCell < cellsInEvent.size(); iCell++) {
1519 caloCellCursor(bcID,
1520 CellHelper::getCellNumber(cellsInEvent[iCell]),
1521 truncateFloatFraction(CellHelper::getAmplitude(cellsInEvent[iCell]), mCaloAmp),
1522 truncateFloatFraction(CellHelper::getTimeStamp(cellsInEvent[iCell]), mCaloTime),
1523 cellsInEvent[iCell].getType(),
1524 caloType); // 1 = emcal, -1 = undefined, 0 = phos
1525
1526 // todo: fix dummy values in CellHelper when it is clear what is filled for trigger information
1527 if (CellHelper::isTRU(cellsInEvent[iCell])) { // Write only trigger cells into this table
1528 caloTRGCursor(bcID,
1529 CellHelper::getFastOrAbsID(cellsInEvent[iCell]),
1530 CellHelper::getLnAmplitude(cellsInEvent[iCell]),
1531 CellHelper::getTriggerBits(cellsInEvent[iCell]),
1532 caloType);
1533 }
1534 if (mUseMC) {
1535 // Common for PHOS and EMCAL
1536 // loop over all MC Labels for the current cell
1537 std::vector<int32_t> particleIds;
1538 std::vector<float> amplitudeFraction;
1539 if (!mEMCselectLeading) {
1540 particleIds.reserve(cellMClabels.size());
1541 amplitudeFraction.reserve(cellMClabels.size());
1542 }
1543 float tmpMaxAmplitude = 0;
1544 int32_t tmpindex = 0;
1545 for (auto& mclabel : cellMClabels[iCell]) {
1546 // do not fill noise lables!
1547 if (mclabel.isValid()) {
1548 if (mEMCselectLeading) {
1549 if (mclabel.getAmplitudeFraction() > tmpMaxAmplitude) {
1550 // Check if this MCparticle added to be kept?
1551 if (mToStore.at(mclabel.getSourceID()).at(mclabel.getEventID()).find(mclabel.getTrackID()) !=
1552 mToStore.at(mclabel.getSourceID()).at(mclabel.getEventID()).end()) {
1553 tmpMaxAmplitude = mclabel.getAmplitudeFraction();
1554 tmpindex = (mToStore.at(mclabel.getSourceID()).at(mclabel.getEventID())).at(mclabel.getTrackID());
1555 }
1556 }
1557 } else {
1558 auto trackStore = mToStore.at(mclabel.getSourceID()).at(mclabel.getEventID());
1559 auto iter = trackStore.find(mclabel.getTrackID());
1560 if (iter != trackStore.end()) {
1561 amplitudeFraction.emplace_back(mclabel.getAmplitudeFraction());
1562 particleIds.emplace_back(iter->second);
1563 } else {
1564 particleIds.emplace_back(-1); // should the mc particle not be in mToStore make sure something (e.g. -1) is saved in particleIds so the length of particleIds is the same es amplitudeFraction!
1565 amplitudeFraction.emplace_back(0.f);
1566 LOG(warn) << "CaloTable: Could not find track for mclabel (" << mclabel.getSourceID() << "," << mclabel.getEventID() << "," << mclabel.getTrackID() << ") in the AOD MC store";
1567 if (mMCKineReader) {
1568 auto mctrack = mMCKineReader->getTrack(mclabel);
1569 TVector3 vec;
1570 mctrack->GetStartVertex(vec);
1571 LOG(warn) << " ... this track is of PDG " << mctrack->GetPdgCode() << " produced by " << mctrack->getProdProcessAsString() << " at (" << vec.X() << "," << vec.Y() << "," << vec.Z() << ")";
1572 }
1573 }
1574 }
1575 }
1576 } // end of loop over all MC Labels for the current cell
1577 if (mEMCselectLeading) {
1578 amplitudeFraction.emplace_back(tmpMaxAmplitude);
1579 particleIds.emplace_back(tmpindex);
1580 }
1581 if (particleIds.size() == 0) {
1582 particleIds.emplace_back(-1);
1583 amplitudeFraction.emplace_back(0.f);
1584 }
1585 mcCaloCellLabelCursor(particleIds,
1586 amplitudeFraction);
1587 }
1588 } // end of loop over cells in current event
1589}
1590
1591// fill calo related tables (cells and calotrigger table)
1592template <typename TCaloCursor, typename TCaloTRGCursor, typename TMCCaloLabelCursor>
1593void AODProducerWorkflowDPL::fillCaloTable(TCaloCursor& caloCellCursor, TCaloTRGCursor& caloTRGCursor,
1594 TMCCaloLabelCursor& mcCaloCellLabelCursor, const std::map<uint64_t, int>& bcsMap,
1596{
1597 // get calo information
1598 auto caloEMCCells = data.getEMCALCells();
1599 auto caloEMCCellsTRGR = data.getEMCALTriggers();
1600 auto mcCaloEMCCellLabels = data.getEMCALCellsMCLabels();
1601
1602 auto caloPHOSCells = data.getPHOSCells();
1603 auto caloPHOSCellsTRGR = data.getPHOSTriggers();
1604 auto mcCaloPHOSCellLabels = data.getPHOSCellsMCLabels();
1605
1606 if (!mInputSources[GIndex::PHS]) {
1607 caloPHOSCells = {};
1608 caloPHOSCellsTRGR = {};
1609 mcCaloPHOSCellLabels = {};
1610 }
1611
1612 if (!mInputSources[GIndex::EMC]) {
1613 caloEMCCells = {};
1614 caloEMCCellsTRGR = {};
1615 mcCaloEMCCellLabels = {};
1616 }
1617
1620
1621 // get cell belonging to an eveffillnt instead of timeframe
1622 emcEventHandler.reset();
1623 emcEventHandler.setCellData(caloEMCCells, caloEMCCellsTRGR);
1624 emcEventHandler.setCellMCTruthContainer(mcCaloEMCCellLabels);
1625
1626 phsEventHandler.reset();
1627 phsEventHandler.setCellData(caloPHOSCells, caloPHOSCellsTRGR);
1628 phsEventHandler.setCellMCTruthContainer(mcCaloPHOSCellLabels);
1629
1630 int emcNEvents = emcEventHandler.getNumberOfEvents();
1631 int phsNEvents = phsEventHandler.getNumberOfEvents();
1632
1633 std::vector<std::tuple<uint64_t, int8_t, int>> caloEvents; // <bc, caloType, eventID>
1634
1635 caloEvents.reserve(emcNEvents + phsNEvents);
1636
1637 for (int iev = 0; iev < emcNEvents; ++iev) {
1638 uint64_t bc = emcEventHandler.getInteractionRecordForEvent(iev).toLong();
1639 caloEvents.emplace_back(std::make_tuple(bc, 1, iev));
1640 }
1641
1642 for (int iev = 0; iev < phsNEvents; ++iev) {
1643 uint64_t bc = phsEventHandler.getInteractionRecordForEvent(iev).toLong();
1644 caloEvents.emplace_back(std::make_tuple(bc, 0, iev));
1645 }
1646
1647 std::sort(caloEvents.begin(), caloEvents.end(),
1648 [](const auto& left, const auto& right) { return std::get<0>(left) < std::get<0>(right); });
1649
1650 // loop over events
1651 for (int i = 0; i < emcNEvents + phsNEvents; ++i) {
1652 uint64_t globalBC = std::get<0>(caloEvents[i]);
1653 int8_t caloType = std::get<1>(caloEvents[i]);
1654 int eventID = std::get<2>(caloEvents[i]);
1655 auto item = bcsMap.find(globalBC);
1656 int bcID = -1;
1657 if (item != bcsMap.end()) {
1658 bcID = item->second;
1659 } else {
1660 LOG(warn) << "Error: could not find a corresponding BC ID for a calo point; globalBC = " << globalBC << ", caloType = " << (int)caloType;
1661 }
1662 if (caloType == 0) { // phos
1663 addToCaloTable(phsEventHandler, caloCellCursor, caloTRGCursor, mcCaloCellLabelCursor, eventID, bcID, caloType);
1664 }
1665 if (caloType == 1) { // emc
1666 addToCaloTable(emcEventHandler, caloCellCursor, caloTRGCursor, mcCaloCellLabelCursor, eventID, bcID, caloType);
1667 }
1668 }
1669
1670 caloEvents.clear();
1671}
1672
1674{
1675 mTimer.Stop();
1677 mLPMProdTag = ic.options().get<std::string>("lpmp-prod-tag");
1678 mAnchorPass = ic.options().get<std::string>("anchor-pass");
1679 mAnchorProd = ic.options().get<std::string>("anchor-prod");
1680 mUser = ic.options().get<std::string>("created-by");
1681 mRecoPass = ic.options().get<std::string>("reco-pass");
1682 mTFNumber = ic.options().get<int64_t>("aod-timeframe-id");
1683 mRecoOnly = ic.options().get<int>("reco-mctracks-only");
1684 mTruncate = ic.options().get<int>("enable-truncation");
1685 mRunNumber = ic.options().get<int>("run-number");
1686 mCTPReadout = ic.options().get<int>("ctpreadout-create");
1687 mNThreads = std::max(1, ic.options().get<int>("nthreads"));
1688 mEMCselectLeading = ic.options().get<bool>("emc-select-leading");
1689 mThinTracks = ic.options().get<bool>("thin-tracks");
1690 mPropTracks = ic.options().get<bool>("propagate-tracks");
1691 mMaxPropXiu = ic.options().get<float>("propagate-tracks-max-xiu");
1692 mPropMuons = ic.options().get<bool>("propagate-muons");
1693 if (auto s = ic.options().get<std::string>("with-streamers"); !s.empty()) {
1694 mStreamerFlags.set(s);
1695 if (mStreamerFlags) {
1696 LOGP(info, "Writing streamer data with mask:");
1697 LOG(info) << mStreamerFlags;
1698 } else {
1699 LOGP(warn, "Specified non-default empty streamer mask!");
1700 }
1701 }
1702 mTrackQCFraction = ic.options().get<float>("trackqc-fraction");
1703 mTrackQCNTrCut = ic.options().get<int64_t>("trackqc-NTrCut");
1704 mTrackQCDCAxy = ic.options().get<float>("trackqc-tpc-dca");
1705 mTrackQCPt = ic.options().get<float>("trackqc-tpc-pt");
1706 mTrackQCNCls = ic.options().get<int>("trackqc-tpc-cls");
1707 if (auto seed = ic.options().get<int>("seed"); seed == 0) {
1708 LOGP(info, "Using random device for seeding");
1709 std::random_device rd;
1710 std::array<int, std::mt19937::state_size> seed_data{};
1711 std::generate(std::begin(seed_data), std::end(seed_data), std::ref(rd));
1712 std::seed_seq seq(std::begin(seed_data), std::end(seed_data));
1713 mGenerator = std::mt19937(seq);
1714 } else {
1715 LOGP(info, "Using seed {} for sampling", seed);
1716 mGenerator.seed(seed);
1717 }
1718#ifdef WITH_OPENMP
1719 LOGP(info, "Multi-threaded parts will run with {} OpenMP threads", mNThreads);
1720#else
1721 mNThreads = 1;
1722 LOG(info) << "OpenMP is disabled";
1723#endif
1724 if (mTFNumber == -1L) {
1725 LOG(info) << "TFNumber will be obtained from CCDB";
1726 }
1727 if (mRunNumber == -1L) {
1728 LOG(info) << "The Run number will be obtained from DPL headers";
1729 }
1730
1731 // set no truncation if selected by user
1732 if (mTruncate != 1) {
1733 LOG(info) << "Truncation is not used!";
1734 mCollisionPosition = 0xFFFFFFFF;
1735 mCollisionPositionCov = 0xFFFFFFFF;
1736 mTrackX = 0xFFFFFFFF;
1737 mTrackAlpha = 0xFFFFFFFF;
1738 mTrackSnp = 0xFFFFFFFF;
1739 mTrackTgl = 0xFFFFFFFF;
1740 mTrack1Pt = 0xFFFFFFFF;
1741 mTrackChi2 = 0xFFFFFFFF;
1742 mTrackCovDiag = 0xFFFFFFFF;
1743 mTrackCovOffDiag = 0xFFFFFFFF;
1744 mTrackSignal = 0xFFFFFFFF;
1745 mTrackTime = 0xFFFFFFFF;
1746 mTPCTime0 = 0xFFFFFFFF;
1747 mTrackTimeError = 0xFFFFFFFF;
1748 mTrackPosEMCAL = 0xFFFFFFFF;
1749 mTracklets = 0xFFFFFFFF;
1750 mMcParticleW = 0xFFFFFFFF;
1751 mMcParticlePos = 0xFFFFFFFF;
1752 mMcParticleMom = 0xFFFFFFFF;
1753 mCaloAmp = 0xFFFFFFFF;
1754 mCaloTime = 0xFFFFFFFF;
1755 mCPVPos = 0xFFFFFFFF;
1756 mCPVAmpl = 0xFFFFFFFF;
1757 mMuonTr1P = 0xFFFFFFFF;
1758 mMuonTrThetaX = 0xFFFFFFFF;
1759 mMuonTrThetaY = 0xFFFFFFFF;
1760 mMuonTrZmu = 0xFFFFFFFF;
1761 mMuonTrBend = 0xFFFFFFFF;
1762 mMuonTrNonBend = 0xFFFFFFFF;
1763 mMuonTrCov = 0xFFFFFFFF;
1764 mMuonCl = 0xFFFFFFFF;
1765 mMuonClErr = 0xFFFFFFFF;
1766 mV0Time = 0xFFFFFFFF;
1767 mV0ChannelTime = 0xFFFFFFFF;
1768 mFDDTime = 0xFFFFFFFF;
1769 mFDDChannelTime = 0xFFFFFFFF;
1770 mT0Time = 0xFFFFFFFF;
1771 mT0ChannelTime = 0xFFFFFFFF;
1772 mV0Amplitude = 0xFFFFFFFF;
1773 mFDDAmplitude = 0xFFFFFFFF;
1774 mT0Amplitude = 0xFFFFFFFF;
1775 }
1776 // Initialize ZDC helper maps
1777 for (int ic = 0; ic < o2::zdc::NChannels; ic++) {
1778 mZDCEnergyMap[ic] = -std::numeric_limits<float>::infinity();
1779 }
1780 for (int ic = 0; ic < o2::zdc::NTDCChannels; ic++) {
1781 mZDCTDCMap[ic] = -std::numeric_limits<float>::infinity();
1782 }
1783
1784 std::string hepmcUpdate = ic.options().get<std::string>("hepmc-update");
1785 HepMCUpdate when = (hepmcUpdate == "never" ? HepMCUpdate::never : hepmcUpdate == "always" ? HepMCUpdate::always
1786 : hepmcUpdate == "all" ? HepMCUpdate::allKeys
1787 : HepMCUpdate::anyKey);
1788 mXSectionUpdate = when;
1789 mPdfInfoUpdate = when;
1790 mHeavyIonUpdate = when;
1791
1792 mTimer.Reset();
1793
1794 if (mStreamerFlags) {
1795 mStreamer = std::make_unique<o2::utils::TreeStreamRedirector>("AO2DStreamer.root", "RECREATE");
1796 }
1797}
1798
1799namespace
1800{
1801void add_additional_meta_info(std::vector<TString>& keys, std::vector<TString>& values)
1802{
1803 // see if we should put additional meta info (e.g. from MC)
1804 auto aod_external_meta_info_file = getenv("AOD_ADDITIONAL_METADATA_FILE");
1805 if (aod_external_meta_info_file != nullptr) {
1806 LOG(info) << "Trying to inject additional AOD meta-data from " << aod_external_meta_info_file;
1807 if (std::filesystem::exists(aod_external_meta_info_file)) {
1808 std::ifstream input_file(aod_external_meta_info_file);
1809 if (input_file) {
1810 nlohmann::json json_data;
1811 try {
1812 input_file >> json_data;
1813 } catch (nlohmann::json::parse_error& e) {
1814 std::cerr << "JSON Parse Error: " << e.what() << "\n";
1815 std::cerr << "Exception ID: " << e.id << "\n";
1816 std::cerr << "Byte position: " << e.byte << "\n";
1817 return;
1818 }
1819 // If parsing succeeds, iterate over key-value pairs
1820 for (const auto& [key, value] : json_data.items()) {
1821 LOG(info) << "Adding AOD MetaData" << key << " : " << value;
1822 keys.push_back(key.c_str());
1823 values.push_back(value.get<std::string>());
1824 }
1825 }
1826 }
1827 }
1828}
1829} // namespace
1830
1832{
1833 mTimer.Start(false);
1835 recoData.collectData(pc, *mDataRequest);
1836 updateTimeDependentParams(pc); // Make sure that this is called after the RecoContainer collect data, since some condition objects are fetched there
1837
1838 mStartIR = recoData.startIR;
1839
1840 auto primVertices = recoData.getPrimaryVertices();
1841 auto primVer2TRefs = recoData.getPrimaryVertexMatchedTrackRefs();
1842 auto primVerGIs = recoData.getPrimaryVertexMatchedTracks();
1843 auto primVerLabels = recoData.getPrimaryVertexMCLabels();
1844
1845 auto fddChData = recoData.getFDDChannelsData();
1846 auto fddRecPoints = recoData.getFDDRecPoints();
1847 auto ft0ChData = recoData.getFT0ChannelsData();
1848 auto ft0RecPoints = recoData.getFT0RecPoints();
1849 auto fv0ChData = recoData.getFV0ChannelsData();
1850 auto fv0RecPoints = recoData.getFV0RecPoints();
1851
1852 auto zdcEnergies = recoData.getZDCEnergy();
1853 auto zdcBCRecData = recoData.getZDCBCRecData();
1854 auto zdcTDCData = recoData.getZDCTDCData();
1855
1856 auto cpvClusters = recoData.getCPVClusters();
1857 auto cpvTrigRecs = recoData.getCPVTriggers();
1858
1859 auto ctpDigits = recoData.getCTPDigits();
1860 const auto& tinfo = pc.services().get<o2::framework::TimingInfo>();
1861 std::vector<o2::ctp::CTPDigit> ctpDigitsCreated;
1862 if (mCTPReadout == 1) {
1863 LOG(info) << "CTP : creating ctpreadout in AOD producer";
1864 createCTPReadout(recoData, ctpDigitsCreated, pc);
1865 LOG(info) << "CTP : ctpreadout created from AOD";
1866 ctpDigits = gsl::span<o2::ctp::CTPDigit>(ctpDigitsCreated);
1867 }
1868 LOG(debug) << "FOUND " << primVertices.size() << " primary vertices";
1869 LOG(debug) << "FOUND " << ft0RecPoints.size() << " FT0 rec. points";
1870 LOG(debug) << "FOUND " << fv0RecPoints.size() << " FV0 rec. points";
1871 LOG(debug) << "FOUND " << fddRecPoints.size() << " FDD rec. points";
1872 LOG(debug) << "FOUND " << cpvClusters.size() << " CPV clusters";
1873 LOG(debug) << "FOUND " << cpvTrigRecs.size() << " CPV trigger records";
1874
1875 LOG(info) << "FOUND " << primVertices.size() << " primary vertices";
1876
1877 using namespace o2::aodhelpers;
1878
1879 auto bcCursor = createTableCursor<o2::aod::BCs>(pc);
1880 auto bcFlagsCursor = createTableCursor<o2::aod::BCFlags>(pc);
1881 auto cascadesCursor = createTableCursor<o2::aod::Cascades>(pc);
1882 auto collisionsCursor = createTableCursor<o2::aod::Collisions>(pc);
1883 auto decay3BodyCursor = createTableCursor<o2::aod::Decay3Bodys>(pc);
1884 auto trackedCascadeCursor = createTableCursor<o2::aod::TrackedCascades>(pc);
1885 auto trackedV0Cursor = createTableCursor<o2::aod::TrackedV0s>(pc);
1886 auto tracked3BodyCurs = createTableCursor<o2::aod::Tracked3Bodys>(pc);
1887 auto fddCursor = createTableCursor<o2::aod::FDDs>(pc);
1888 auto fddExtraCursor = createTableCursor<o2::aod::FDDsExtra>(pc);
1889 auto ft0Cursor = createTableCursor<o2::aod::FT0s>(pc);
1890 auto ft0ExtraCursor = createTableCursor<o2::aod::FT0sExtra>(pc);
1891 auto fv0aCursor = createTableCursor<o2::aod::FV0As>(pc);
1892 auto fv0aExtraCursor = createTableCursor<o2::aod::FV0AsExtra>(pc);
1893 auto fwdTracksCursor = createTableCursor<o2::aod::StoredFwdTracks>(pc);
1894 auto fwdTracksCovCursor = createTableCursor<o2::aod::StoredFwdTracksCov>(pc);
1895 auto fwdTrkClsCursor = createTableCursor<o2::aod::FwdTrkCls>(pc);
1896 auto mftTracksCursor = createTableCursor<o2::aod::StoredMFTTracks>(pc);
1897 auto mftTracksCovCursor = createTableCursor<o2::aod::StoredMFTTracksCov>(pc);
1898 auto tracksCursor = createTableCursor<o2::aod::StoredTracksIU>(pc);
1899 auto tracksCovCursor = createTableCursor<o2::aod::StoredTracksCovIU>(pc);
1900 auto tracksExtraCursor = createTableCursor<o2::aod::StoredTracksExtra>(pc);
1901 auto tracksQACursor = createTableCursor<o2::aod::TracksQAVersion>(pc);
1902 auto ambigTracksCursor = createTableCursor<o2::aod::AmbiguousTracks>(pc);
1903 auto ambigMFTTracksCursor = createTableCursor<o2::aod::AmbiguousMFTTracks>(pc);
1904 auto ambigFwdTracksCursor = createTableCursor<o2::aod::AmbiguousFwdTracks>(pc);
1905 auto v0sCursor = createTableCursor<o2::aod::V0s>(pc);
1906 auto zdcCursor = createTableCursor<o2::aod::Zdcs>(pc);
1907 auto hmpCursor = createTableCursor<o2::aod::HMPIDs>(pc);
1908 auto caloCellsCursor = createTableCursor<o2::aod::Calos>(pc);
1909 auto caloCellsTRGTableCursor = createTableCursor<o2::aod::CaloTriggers>(pc);
1910 auto cpvClustersCursor = createTableCursor<o2::aod::CPVClusters>(pc);
1911 auto originCursor = createTableCursor<o2::aod::Origins>(pc);
1912
1913 // Declare MC cursors type without adding the output for a table
1924 if (mUseMC) { // This creates the actual writercursor
1925 mcColLabelsCursor = createTableCursor<o2::aod::McCollisionLabels>(pc);
1926 mcCollisionsCursor = createTableCursor<o2::aod::McCollisions>(pc);
1927 hepmcXSectionsCursor = createTableCursor<o2::aod::HepMCXSections>(pc);
1928 hepmcPdfInfosCursor = createTableCursor<o2::aod::HepMCPdfInfos>(pc);
1929 hepmcHeavyIonsCursor = createTableCursor<o2::aod::HepMCHeavyIons>(pc);
1930 mcMFTTrackLabelCursor = createTableCursor<o2::aod::McMFTTrackLabels>(pc);
1931 mcFwdTrackLabelCursor = createTableCursor<o2::aod::McFwdTrackLabels>(pc);
1932 mcParticlesCursor = createTableCursor<o2::aod::StoredMcParticles_001>(pc);
1933 mcTrackLabelCursor = createTableCursor<o2::aod::McTrackLabels>(pc);
1934 mcCaloLabelsCursor = createTableCursor<o2::aod::McCaloLabels_001>(pc);
1935 }
1936
1937 std::unique_ptr<o2::steer::MCKinematicsReader> mcReader;
1938 if (mUseMC) {
1939 mcReader = std::make_unique<o2::steer::MCKinematicsReader>("collisioncontext.root");
1940 }
1941 mMCKineReader = mcReader.get(); // for use in different functions
1942 std::map<uint64_t, int> bcsMap;
1943 collectBCs(recoData, mUseMC ? mcReader->getDigitizationContext()->getEventRecords() : std::vector<o2::InteractionTimeRecord>{}, bcsMap);
1944 if (!primVer2TRefs.empty()) { // if the vertexing was done, the last slot refers to orphan tracks
1945 addRefGlobalBCsForTOF(primVer2TRefs.back(), primVerGIs, recoData, bcsMap);
1946 }
1947 // initialize the bunch crossing container for further use below
1948 mBCLookup.init(bcsMap);
1949
1950 uint64_t tfNumber;
1951 const int runNumber = (mRunNumber == -1) ? int(tinfo.runNumber) : mRunNumber;
1952 if (mTFNumber == -1L) {
1953 // TODO has to use absolute time of TF
1954 tfNumber = uint64_t(tinfo.firstTForbit) + (uint64_t(tinfo.runNumber) << 32); // getTFNumber(mStartIR, runNumber);
1955 } else {
1956 tfNumber = mTFNumber;
1957 }
1958
1959 std::vector<float> aAmplitudes, aTimes;
1960 std::vector<uint8_t> aChannels;
1961 fv0aCursor.reserve(fv0RecPoints.size());
1962 for (auto& fv0RecPoint : fv0RecPoints) {
1963 aAmplitudes.clear();
1964 aChannels.clear();
1965 aTimes.clear();
1966 const auto channelData = fv0RecPoint.getBunchChannelData(fv0ChData);
1967 for (auto& channel : channelData) {
1968 if (channel.charge > 0) {
1969 aAmplitudes.push_back(truncateFloatFraction(channel.charge, mV0Amplitude));
1970 aTimes.push_back(truncateFloatFraction(channel.time * 1.E-3, mV0ChannelTime));
1971 aChannels.push_back(channel.channel);
1972 }
1973 }
1974 uint64_t bc = fv0RecPoint.getInteractionRecord().toLong();
1975 auto item = bcsMap.find(bc);
1976 int bcID = -1;
1977 if (item != bcsMap.end()) {
1978 bcID = item->second;
1979 } else {
1980 LOG(fatal) << "Error: could not find a corresponding BC ID for a FV0 rec. point; BC = " << bc;
1981 }
1982 fv0aCursor(bcID,
1983 aAmplitudes,
1984 aChannels,
1985 truncateFloatFraction(fv0RecPoint.getCollisionGlobalMeanTime() * 1E-3, mV0Time), // ps to ns
1986 fv0RecPoint.getTrigger().getTriggersignals());
1987
1988 if (mEnableFITextra) {
1989 fv0aExtraCursor(bcID,
1990 aTimes);
1991 }
1992 }
1993
1994 std::vector<float> zdcEnergy, zdcAmplitudes, zdcTime;
1995 std::vector<uint8_t> zdcChannelsE, zdcChannelsT;
1996 zdcCursor.reserve(zdcBCRecData.size());
1997 for (auto zdcRecData : zdcBCRecData) {
1998 uint64_t bc = zdcRecData.ir.toLong();
1999 auto item = bcsMap.find(bc);
2000 int bcID = -1;
2001 if (item != bcsMap.end()) {
2002 bcID = item->second;
2003 } else {
2004 LOG(fatal) << "Error: could not find a corresponding BC ID for a ZDC rec. point; BC = " << bc;
2005 }
2006 int fe, ne, ft, nt, fi, ni;
2007 zdcRecData.getRef(fe, ne, ft, nt, fi, ni);
2008 zdcEnergy.clear();
2009 zdcChannelsE.clear();
2010 zdcAmplitudes.clear();
2011 zdcTime.clear();
2012 zdcChannelsT.clear();
2013 for (int ie = 0; ie < ne; ie++) {
2014 auto& zdcEnergyData = zdcEnergies[fe + ie];
2015 zdcEnergy.emplace_back(zdcEnergyData.energy());
2016 zdcChannelsE.emplace_back(zdcEnergyData.ch());
2017 }
2018 for (int it = 0; it < nt; it++) {
2019 auto& tdc = zdcTDCData[ft + it];
2020 zdcAmplitudes.emplace_back(tdc.amplitude());
2021 zdcTime.emplace_back(tdc.value());
2022 zdcChannelsT.emplace_back(o2::zdc::TDCSignal[tdc.ch()]);
2023 }
2024 zdcCursor(bcID,
2025 zdcEnergy,
2026 zdcChannelsE,
2027 zdcAmplitudes,
2028 zdcTime,
2029 zdcChannelsT);
2030 }
2031
2032 // keep track event/source id for each mc-collision
2033 // using map and not unordered_map to ensure
2034 // correct ordering when iterating over container elements
2035 std::vector<std::vector<int>> mcColToEvSrc;
2036
2037 if (mUseMC) {
2038 using namespace o2::aodmchelpers;
2039
2040 // filling mcCollision table
2041 int nMCCollisions = mcReader->getDigitizationContext()->getNCollisions();
2042 const auto& mcRecords = mcReader->getDigitizationContext()->getEventRecords();
2043 const auto& mcParts = mcReader->getDigitizationContext()->getEventParts();
2044
2045 // count all parts
2046 int totalNParts = 0;
2047 for (int iCol = 0; iCol < nMCCollisions; iCol++) {
2048 totalNParts += mcParts[iCol].size();
2049 }
2050 mcCollisionsCursor.reserve(totalNParts);
2051
2052 for (int iCol = 0; iCol < nMCCollisions; iCol++) {
2053 const auto time = mcRecords[iCol].getTimeOffsetWrtBC();
2054 auto globalBC = mcRecords[iCol].toLong();
2055 auto item = bcsMap.find(globalBC);
2056 int bcID = -1;
2057 if (item != bcsMap.end()) {
2058 bcID = item->second;
2059 } else {
2060 LOG(fatal) << "Error: could not find a corresponding BC ID "
2061 << "for MC collision; BC = " << globalBC
2062 << ", mc collision = " << iCol;
2063 }
2064 auto& colParts = mcParts[iCol];
2065 auto nParts = colParts.size();
2066 for (auto colPart : colParts) {
2067 auto eventID = colPart.entryID;
2068 auto sourceID = colPart.sourceID;
2069 // enable embedding: if several colParts exist, then they are
2070 // saved as one collision
2071 if (nParts == 1 || sourceID == 0) {
2072 // FIXME:
2073 // use generators' names for generatorIDs (?)
2074 auto& header = mcReader->getMCEventHeader(sourceID, eventID);
2075 updateMCHeader(mcCollisionsCursor.cursor,
2076 hepmcXSectionsCursor.cursor,
2077 hepmcPdfInfosCursor.cursor,
2078 hepmcHeavyIonsCursor.cursor,
2079 header,
2080 iCol,
2081 bcID,
2082 time,
2083 0,
2084 sourceID);
2085 }
2086 mcColToEvSrc.emplace_back(std::vector<int>{iCol, sourceID, eventID}); // point background and injected signal events to one collision
2087 }
2088 }
2089 }
2090
2091 std::sort(mcColToEvSrc.begin(), mcColToEvSrc.end(),
2092 [](const std::vector<int>& left, const std::vector<int>& right) { return (left[0] < right[0]); });
2093
2094 // vector of FDD amplitudes
2095 int16_t aFDDAmplitudesA[8] = {0u}, aFDDAmplitudesC[8] = {0u};
2096 float aFDDTimesA[8] = {0.f}, aFDDTimesC[8] = {0.f};
2097 // filling FDD table
2098 fddCursor.reserve(fddRecPoints.size());
2099 for (const auto& fddRecPoint : fddRecPoints) {
2100 for (int i = 0; i < 8; i++) {
2101 aFDDAmplitudesA[i] = 0;
2102 aFDDAmplitudesC[i] = 0;
2103 aFDDTimesA[i] = 0.f;
2104 aFDDTimesC[i] = 0.f;
2105 }
2106 uint64_t globalBC = fddRecPoint.getInteractionRecord().toLong();
2107 uint64_t bc = globalBC;
2108 auto item = bcsMap.find(bc);
2109 int bcID = -1;
2110 if (item != bcsMap.end()) {
2111 bcID = item->second;
2112 } else {
2113 LOG(fatal) << "Error: could not find a corresponding BC ID for a FDD rec. point; BC = " << bc;
2114 }
2115 const auto channelData = fddRecPoint.getBunchChannelData(fddChData);
2116 for (const auto& channel : channelData) {
2117 if (channel.mPMNumber < 8) {
2118 aFDDAmplitudesC[channel.mPMNumber] = channel.mChargeADC; // amplitude
2119 aFDDTimesC[channel.mPMNumber] = truncateFloatFraction(channel.mTime * 1E-3, mFDDChannelTime); // time
2120 } else {
2121 aFDDAmplitudesA[channel.mPMNumber - 8] = channel.mChargeADC; // amplitude
2122 aFDDTimesA[channel.mPMNumber - 8] = truncateFloatFraction(channel.mTime * 1E-3, mFDDChannelTime); // time
2123 }
2124 }
2125
2126 fddCursor(bcID,
2127 aFDDAmplitudesA,
2128 aFDDAmplitudesC,
2129 truncateFloatFraction(fddRecPoint.getCollisionTimeA() * 1E-3, mFDDTime), // ps to ns
2130 truncateFloatFraction(fddRecPoint.getCollisionTimeC() * 1E-3, mFDDTime), // ps to ns
2131 fddRecPoint.getTrigger().getTriggersignals());
2132 if (mEnableFITextra) {
2133 fddExtraCursor(bcID,
2134 aFDDTimesA,
2135 aFDDTimesC);
2136 }
2137 }
2138
2139 // filling FT0 table
2140 std::vector<float> aAmplitudesA, aAmplitudesC, aTimesA, aTimesC;
2141 std::vector<uint8_t> aChannelsA, aChannelsC;
2142 ft0Cursor.reserve(ft0RecPoints.size());
2143 for (auto& ft0RecPoint : ft0RecPoints) {
2144 aAmplitudesA.clear();
2145 aAmplitudesC.clear();
2146 aTimesA.clear();
2147 aTimesC.clear();
2148 aChannelsA.clear();
2149 aChannelsC.clear();
2150 const auto channelData = ft0RecPoint.getBunchChannelData(ft0ChData);
2151 for (auto& channel : channelData) {
2152 // TODO: switch to calibrated amplitude
2153 if (channel.QTCAmpl > 0) {
2154 constexpr int nFT0ChannelsAside = o2::ft0::Geometry::NCellsA * 4;
2155 if (channel.ChId < nFT0ChannelsAside) {
2156 aChannelsA.push_back(channel.ChId);
2157 aAmplitudesA.push_back(truncateFloatFraction(channel.QTCAmpl, mT0Amplitude));
2158 aTimesA.push_back(truncateFloatFraction(channel.CFDTime * 1E-3, mT0ChannelTime));
2159 } else {
2160 aChannelsC.push_back(channel.ChId - nFT0ChannelsAside);
2161 aAmplitudesC.push_back(truncateFloatFraction(channel.QTCAmpl, mT0Amplitude));
2162 aTimesC.push_back(truncateFloatFraction(channel.CFDTime * 1E-3, mT0ChannelTime));
2163 }
2164 }
2165 }
2166 uint64_t globalBC = ft0RecPoint.getInteractionRecord().toLong();
2167 uint64_t bc = globalBC;
2168 auto item = bcsMap.find(bc);
2169 int bcID = -1;
2170 if (item != bcsMap.end()) {
2171 bcID = item->second;
2172 } else {
2173 LOG(fatal) << "Error: could not find a corresponding BC ID for a FT0 rec. point; BC = " << bc;
2174 }
2175 ft0Cursor(bcID,
2176 aAmplitudesA,
2177 aChannelsA,
2178 aAmplitudesC,
2179 aChannelsC,
2180 truncateFloatFraction(ft0RecPoint.getCollisionTimeA() * 1E-3, mT0Time), // ps to ns
2181 truncateFloatFraction(ft0RecPoint.getCollisionTimeC() * 1E-3, mT0Time), // ps to ns
2182 ft0RecPoint.getTrigger().getTriggersignals());
2183 if (mEnableFITextra) {
2184 ft0ExtraCursor(bcID,
2185 aTimesA,
2186 aTimesC);
2187 }
2188 }
2189
2190 if (mUseMC) {
2191 // filling MC collision labels
2192 mcColLabelsCursor.reserve(primVerLabels.size());
2193 for (auto& label : primVerLabels) {
2194 auto it = std::find_if(mcColToEvSrc.begin(), mcColToEvSrc.end(),
2195 [&label](const auto& mcColInfo) { return mcColInfo[1] == label.getSourceID() && mcColInfo[2] == label.getEventID(); });
2196 int32_t mcCollisionID = -1;
2197 if (it != mcColToEvSrc.end()) {
2198 mcCollisionID = it->at(0);
2199 }
2200 uint16_t mcMask = 0; // todo: set mask using normalized weights?
2201 mcColLabelsCursor(mcCollisionID, mcMask);
2202 }
2203 }
2204
2205 cacheTriggers(recoData);
2206 countTPCClusters(recoData);
2207
2208 int collisionID = 0;
2209 mIndexTableMFT.resize(recoData.getMFTTracks().size());
2210 mIndexTableFwd.resize(recoData.getMCHTracks().size());
2211
2212 auto& trackReffwd = primVer2TRefs.back();
2213 fillIndexTablesPerCollision(trackReffwd, primVerGIs, recoData);
2214 collisionID = 0;
2215 for (auto& vertex : primVertices) {
2216 auto& trackReffwd = primVer2TRefs[collisionID];
2217 fillIndexTablesPerCollision(trackReffwd, primVerGIs, recoData); // this function must follow the same track order as 'fillTrackTablesPerCollision' to fill the map of track indices
2218 collisionID++;
2219 }
2220
2222 prepareStrangenessTracking(recoData);
2223
2224 mGIDToTableFwdID.clear(); // reset the tables to be used by 'fillTrackTablesPerCollision'
2225 mGIDToTableMFTID.clear();
2226
2227 if (mPropTracks || mThinTracks) {
2228 auto v0s = recoData.getV0sIdx();
2229 auto cascades = recoData.getCascadesIdx();
2230 auto decays3Body = recoData.getDecays3BodyIdx();
2231 mGIDUsedBySVtx.reserve(v0s.size() * 2 + cascades.size() + decays3Body.size() * 3);
2232 for (const auto& v0 : v0s) {
2233 mGIDUsedBySVtx.insert(v0.getProngID(0));
2234 mGIDUsedBySVtx.insert(v0.getProngID(1));
2235 }
2236 for (const auto& cascade : cascades) {
2237 mGIDUsedBySVtx.insert(cascade.getBachelorID());
2238 }
2239 for (const auto& id3Body : decays3Body) {
2240 mGIDUsedBySVtx.insert(id3Body.getProngID(0));
2241 mGIDUsedBySVtx.insert(id3Body.getProngID(1));
2242 mGIDUsedBySVtx.insert(id3Body.getProngID(2));
2243 }
2244
2245 mGIDUsedByStr.reserve(recoData.getStrangeTracks().size());
2246 for (const auto& sTrk : recoData.getStrangeTracks()) {
2247 mGIDUsedByStr.emplace(sTrk.mITSRef, GIndex::ITS);
2248 }
2249 }
2250
2251 // filling unassigned tracks first
2252 // so that all unassigned tracks are stored in the beginning of the table together
2253 auto& trackRef = primVer2TRefs.back(); // references to unassigned tracks are at the end
2254 // fixme: interaction time is undefined for unassigned tracks (?)
2255 fillTrackTablesPerCollision(-1, std::uint64_t(-1), trackRef, primVerGIs, recoData, tracksCursor, tracksCovCursor, tracksExtraCursor, tracksQACursor,
2256 ambigTracksCursor, mftTracksCursor, mftTracksCovCursor, ambigMFTTracksCursor,
2257 fwdTracksCursor, fwdTracksCovCursor, ambigFwdTracksCursor, fwdTrkClsCursor, bcsMap);
2258
2259 // filling collisions and tracks into tables
2260 collisionID = 0;
2261 collisionsCursor.reserve(primVertices.size());
2262 for (auto& vertex : primVertices) {
2263 auto& cov = vertex.getCov();
2264 auto& timeStamp = vertex.getTimeStamp(); // this is a relative time
2265 const double interactionTime = timeStamp.getTimeStamp() * 1E3; // mus to ns
2266 uint64_t globalBC = relativeTime_to_GlobalBC(interactionTime);
2267 uint64_t localBC = relativeTime_to_LocalBC(interactionTime);
2268 LOG(debug) << "global BC " << globalBC << " local BC " << localBC << " relative interaction time " << interactionTime;
2269 // collision timestamp in ns wrt the beginning of collision BC
2270 const float relInteractionTime = static_cast<float>(localBC * o2::constants::lhc::LHCBunchSpacingNS - interactionTime);
2271 auto item = bcsMap.find(globalBC);
2272 int bcID = -1;
2273 if (item != bcsMap.end()) {
2274 bcID = item->second;
2275 } else {
2276 LOG(fatal) << "Error: could not find a corresponding BC ID for a collision; BC = " << globalBC << ", collisionID = " << collisionID;
2277 }
2278 collisionsCursor(bcID,
2279 truncateFloatFraction(vertex.getX(), mCollisionPosition),
2280 truncateFloatFraction(vertex.getY(), mCollisionPosition),
2281 truncateFloatFraction(vertex.getZ(), mCollisionPosition),
2282 truncateFloatFraction(cov[0], mCollisionPositionCov),
2283 truncateFloatFraction(cov[1], mCollisionPositionCov),
2284 truncateFloatFraction(cov[2], mCollisionPositionCov),
2285 truncateFloatFraction(cov[3], mCollisionPositionCov),
2286 truncateFloatFraction(cov[4], mCollisionPositionCov),
2287 truncateFloatFraction(cov[5], mCollisionPositionCov),
2288 vertex.getFlags(),
2289 truncateFloatFraction(vertex.getChi2(), mCollisionPositionCov),
2290 vertex.getNContributors(),
2291 truncateFloatFraction(relInteractionTime, mCollisionPosition),
2292 truncateFloatFraction(timeStamp.getTimeStampError() * 1E3, mCollisionPositionCov));
2293 mVtxToTableCollID[collisionID] = mTableCollID++;
2294
2295 auto& trackRef = primVer2TRefs[collisionID];
2296 // passing interaction time in [ps]
2297 fillTrackTablesPerCollision(collisionID, globalBC, trackRef, primVerGIs, recoData, tracksCursor, tracksCovCursor, tracksExtraCursor, tracksQACursor, ambigTracksCursor,
2298 mftTracksCursor, mftTracksCovCursor, ambigMFTTracksCursor,
2299 fwdTracksCursor, fwdTracksCovCursor, ambigFwdTracksCursor, fwdTrkClsCursor, bcsMap);
2300 collisionID++;
2301 }
2302
2303 fillSecondaryVertices(recoData, v0sCursor, cascadesCursor, decay3BodyCursor);
2304 fillHMPID(recoData, hmpCursor);
2305 fillStrangenessTrackingTables(recoData, trackedV0Cursor, trackedCascadeCursor, tracked3BodyCurs);
2306
2307 // helper map for fast search of a corresponding class mask for a bc
2308 auto emcalIncomplete = filterEMCALIncomplete(recoData.getEMCALTriggers());
2309 std::unordered_map<uint64_t, std::pair<uint64_t, uint64_t>> bcToClassMask;
2310 if (mInputSources[GID::CTP]) {
2311 LOG(debug) << "CTP input available";
2312 for (auto& ctpDigit : ctpDigits) {
2313 uint64_t bc = ctpDigit.intRecord.toLong();
2314 uint64_t classMask = ctpDigit.CTPClassMask.to_ulong();
2315 uint64_t inputMask = ctpDigit.CTPInputMask.to_ulong();
2316 if (emcalIncomplete.find(bc) != emcalIncomplete.end()) {
2317 // reject EMCAL triggers as BC was rejected as incomplete at readout level
2318 auto classMaskOrig = classMask;
2319 classMask = classMask & ~mEMCALTrgClassMask;
2320 LOG(debug) << "Found EMCAL incomplete event, mask before " << std::bitset<64>(classMaskOrig) << ", after " << std::bitset<64>(classMask);
2321 }
2322 bcToClassMask[bc] = {classMask, inputMask};
2323 // LOG(debug) << Form("classmask:0x%llx", classMask);
2324 }
2325 }
2326
2327 // filling BC table
2328 bcCursor.reserve(bcsMap.size());
2329 for (auto& item : bcsMap) {
2330 uint64_t bc = item.first;
2331 std::pair<uint64_t, uint64_t> masks{0, 0};
2332 if (mInputSources[GID::CTP]) {
2333 auto bcClassPair = bcToClassMask.find(bc);
2334 if (bcClassPair != bcToClassMask.end()) {
2335 masks = bcClassPair->second;
2336 }
2337 }
2338 bcCursor(runNumber,
2339 bc,
2340 masks.first,
2341 masks.second);
2342 }
2343
2344 bcToClassMask.clear();
2345
2346 // filling BC flags table:
2347 auto bcFlags = fillBCFlags(recoData, bcsMap);
2348 bcFlagsCursor.reserve(bcFlags.size());
2349 for (auto f : bcFlags) {
2350 bcFlagsCursor(f);
2351 }
2352
2353 // fill cpvcluster table
2354 if (mInputSources[GIndex::CPV]) {
2355 float posX, posZ;
2356 cpvClustersCursor.reserve(cpvTrigRecs.size());
2357 for (auto& cpvEvent : cpvTrigRecs) {
2358 uint64_t bc = cpvEvent.getBCData().toLong();
2359 auto item = bcsMap.find(bc);
2360 int bcID = -1;
2361 if (item != bcsMap.end()) {
2362 bcID = item->second;
2363 } else {
2364 LOG(fatal) << "Error: could not find a corresponding BC ID for a CPV Trigger Record; BC = " << bc;
2365 }
2366 for (int iClu = cpvEvent.getFirstEntry(); iClu < cpvEvent.getFirstEntry() + cpvEvent.getNumberOfObjects(); iClu++) {
2367 auto& clu = cpvClusters[iClu];
2368 clu.getLocalPosition(posX, posZ);
2369 cpvClustersCursor(bcID,
2370 truncateFloatFraction(posX, mCPVPos),
2371 truncateFloatFraction(posZ, mCPVPos),
2372 truncateFloatFraction(clu.getEnergy(), mCPVAmpl),
2373 clu.getPackedClusterStatus());
2374 }
2375 }
2376 }
2377
2378 if (mUseMC) {
2379 TStopwatch timer;
2380 timer.Start();
2381 // filling mc particles table
2382 fillMCParticlesTable(*mcReader,
2383 mcParticlesCursor.cursor,
2384 primVer2TRefs,
2385 primVerGIs,
2386 recoData,
2387 mcColToEvSrc);
2388 timer.Stop();
2389 LOG(info) << "FILL MC took " << timer.RealTime() << " s";
2390 mcColToEvSrc.clear();
2391
2392 // ------------------------------------------------------
2393 // filling track labels
2394
2395 // need to go through labels in the same order as for tracks
2396 fillMCTrackLabelsTable(mcTrackLabelCursor, mcMFTTrackLabelCursor, mcFwdTrackLabelCursor, primVer2TRefs.back(), primVerGIs, recoData);
2397 for (auto iref = 0U; iref < primVer2TRefs.size() - 1; iref++) {
2398 auto& trackRef = primVer2TRefs[iref];
2399 fillMCTrackLabelsTable(mcTrackLabelCursor, mcMFTTrackLabelCursor, mcFwdTrackLabelCursor, trackRef, primVerGIs, recoData, iref);
2400 }
2401 }
2402
2403 // Fill calo tables and if MC also the MCCaloTable, therefore, has to be after fillMCParticlesTable call!
2404 if (mInputSources[GIndex::PHS] || mInputSources[GIndex::EMC]) {
2405 fillCaloTable(caloCellsCursor, caloCellsTRGTableCursor, mcCaloLabelsCursor, bcsMap, recoData);
2406 }
2407
2408 bcsMap.clear();
2409 clearMCKeepStore(mToStore);
2410 mGIDToTableID.clear();
2411 mTableTrID = 0;
2412 mGIDToTableFwdID.clear();
2413 mTableTrFwdID = 0;
2414 mGIDToTableMFTID.clear();
2415 mTableTrMFTID = 0;
2416 mVtxToTableCollID.clear();
2417 mTableCollID = 0;
2418 mV0ToTableID.clear();
2419 mTableV0ID = 0;
2420
2421 mIndexTableFwd.clear();
2422 mIndexFwdID = 0;
2423 mIndexTableMFT.clear();
2424 mIndexMFTID = 0;
2425
2426 mBCLookup.clear();
2427
2428 mGIDUsedBySVtx.clear();
2429 mGIDUsedByStr.clear();
2430
2431 originCursor(tfNumber);
2432
2433 // sending metadata to writer
2434 TString dataType = mUseMC ? "MC" : "RAW";
2435 TString O2Version = o2::fullVersion();
2436 TString ROOTVersion = ROOT_RELEASE;
2437 mMetaDataKeys = {"DataType", "Run", "O2Version", "ROOTVersion", "RecoPassName", "AnchorProduction", "AnchorPassName", "LPMProductionTag", "CreatedBy"};
2438 mMetaDataVals = {dataType, "3", O2Version, ROOTVersion, mRecoPass, mAnchorProd, mAnchorPass, mLPMProdTag, mUser};
2439 add_additional_meta_info(mMetaDataKeys, mMetaDataVals);
2440
2441 pc.outputs().snapshot(Output{"AMD", "AODMetadataKeys", 0}, mMetaDataKeys);
2442 pc.outputs().snapshot(Output{"AMD", "AODMetadataVals", 0}, mMetaDataVals);
2443
2444 pc.outputs().snapshot(Output{"TFN", "TFNumber", 0}, tfNumber);
2445 pc.outputs().snapshot(Output{"TFF", "TFFilename", 0}, "");
2446
2447 mTimer.Stop();
2448}
2449
2450void AODProducerWorkflowDPL::cacheTriggers(const o2::globaltracking::RecoContainer& recoData)
2451{
2452 // ITS tracks->ROF
2453 {
2454 mITSROFs.clear();
2455 const auto& rofs = recoData.getITSTracksROFRecords();
2456 uint16_t count = 0;
2457 for (const auto& rof : rofs) {
2458 int first = rof.getFirstEntry(), last = first + rof.getNEntries();
2459 for (int i = first; i < last; i++) {
2460 mITSROFs.push_back(count);
2461 }
2462 count++;
2463 }
2464 }
2465 // MFT tracks->ROF
2466 {
2467 mMFTROFs.clear();
2468 const auto& rofs = recoData.getMFTTracksROFRecords();
2469 uint16_t count = 0;
2470 for (const auto& rof : rofs) {
2471 int first = rof.getFirstEntry(), last = first + rof.getNEntries();
2472 for (int i = first; i < last; i++) {
2473 mMFTROFs.push_back(count);
2474 }
2475 count++;
2476 }
2477 }
2478 // ITSTPCTRD tracks -> TRD trigger
2479 {
2480 mITSTPCTRDTriggers.clear();
2481 const auto& itstpctrigs = recoData.getITSTPCTRDTriggers();
2482 int count = 0;
2483 for (const auto& trig : itstpctrigs) {
2484 int first = trig.getFirstTrack(), last = first + trig.getNumberOfTracks();
2485 for (int i = first; i < last; i++) {
2486 mITSTPCTRDTriggers.push_back(count);
2487 }
2488 count++;
2489 }
2490 }
2491 // TPCTRD tracks -> TRD trigger
2492 {
2493 mTPCTRDTriggers.clear();
2494 const auto& tpctrigs = recoData.getTPCTRDTriggers();
2495 int count = 0;
2496 for (const auto& trig : tpctrigs) {
2497 int first = trig.getFirstTrack(), last = first + trig.getNumberOfTracks();
2498 for (int i = first; i < last; i++) {
2499 mTPCTRDTriggers.push_back(count);
2500 }
2501 count++;
2502 }
2503 }
2504 // MCH tracks->ROF
2505 {
2506 mMCHROFs.clear();
2507 const auto& rofs = recoData.getMCHTracksROFRecords();
2508 uint16_t count = 0;
2509 for (const auto& rof : rofs) {
2510 int first = rof.getFirstIdx(), last = first + rof.getNEntries();
2511 for (int i = first; i < last; i++) {
2512 mMCHROFs.push_back(count);
2513 }
2514 count++;
2515 }
2516 }
2517}
2518
2519AODProducerWorkflowDPL::TrackExtraInfo AODProducerWorkflowDPL::processBarrelTrack(int collisionID, std::uint64_t collisionBC, GIndex trackIndex,
2520 const o2::globaltracking::RecoContainer& data, const std::map<uint64_t, int>& bcsMap)
2521{
2522 TrackExtraInfo extraInfoHolder;
2523 if (collisionID < 0) {
2524 extraInfoHolder.flags |= o2::aod::track::OrphanTrack;
2525 }
2526 bool needBCSlice = collisionID < 0; // track is associated to multiple vertices
2527 uint64_t bcOfTimeRef = collisionBC - mStartIR.toLong(); // by default track time is wrt collision BC (unless no collision assigned)
2528
2529 auto setTrackTime = [&](double t, double terr, bool gaussian) {
2530 // set track time and error, for ambiguous tracks define the bcSlice as it was used in vertex-track association
2531 // provided track time (wrt TF start) and its error should be in ns, gaussian flag tells if the error is assumed to be gaussin or half-interval
2532 if (!gaussian) {
2533 extraInfoHolder.flags |= o2::aod::track::TrackTimeResIsRange;
2534 }
2535 extraInfoHolder.trackTimeRes = terr;
2536 if (needBCSlice) { // need to define BC slice
2537 double error = this->mTimeMarginTrackTime + (gaussian ? extraInfoHolder.trackTimeRes * this->mNSigmaTimeTrack : extraInfoHolder.trackTimeRes);
2538 bcOfTimeRef = fillBCSlice(extraInfoHolder.bcSlice, t - error, t + error, bcsMap);
2539 }
2540 extraInfoHolder.trackTime = float(t - bcOfTimeRef * o2::constants::lhc::LHCBunchSpacingNS);
2541 extraInfoHolder.diffBCRef = int(bcOfTimeRef);
2542 LOGP(debug, "time : {}/{} -> {}/{} -> trunc: {}/{} CollID: {} Amb: {}", t, terr, t - bcOfTimeRef * o2::constants::lhc::LHCBunchSpacingNS, terr,
2543 truncateFloatFraction(extraInfoHolder.trackTime, mTrackTime), truncateFloatFraction(extraInfoHolder.trackTimeRes, mTrackTimeError),
2544 collisionID, trackIndex.isAmbiguous());
2545 };
2546 auto contributorsGID = data.getSingleDetectorRefs(trackIndex);
2547 const auto& trackPar = data.getTrackParam(trackIndex);
2548 extraInfoHolder.flags |= trackPar.getPID() << 28;
2549 auto src = trackIndex.getSource();
2550 if (contributorsGID[GIndex::Source::TOF].isIndexSet()) { // ITS-TPC-TRD-TOF, ITS-TPC-TOF, TPC-TRD-TOF, TPC-TOF
2551 const auto& tofMatch = data.getTOFMatch(trackIndex);
2552 extraInfoHolder.tofChi2 = tofMatch.getChi2();
2553 const auto& tofInt = tofMatch.getLTIntegralOut();
2554 float intLen = tofInt.getL();
2555 extraInfoHolder.length = intLen;
2556 const float mass = o2::constants::physics::MassPionCharged; // default pid = pion
2557 if (tofInt.getTOF(o2::track::PID::Pion) > 0.f) {
2558 float expBeta = (intLen / (tofInt.getTOF(o2::track::PID::Pion) * cSpeed));
2559 if (expBeta > o2::constants::math::Almost1) {
2561 }
2562 extraInfoHolder.tofExpMom = mass * expBeta / std::sqrt(1.f - expBeta * expBeta);
2563 }
2564 // correct the time of the track
2565 const double massZ = o2::track::PID::getMass2Z(trackPar.getPID());
2566 const double energy = sqrt((massZ * massZ) + (extraInfoHolder.tofExpMom * extraInfoHolder.tofExpMom));
2567 const double exp = extraInfoHolder.length * energy / (cSpeed * extraInfoHolder.tofExpMom);
2568 auto tofSignal = (tofMatch.getSignal() - exp) * 1e-3; // time in ns wrt TF start
2569 setTrackTime(tofSignal, 0.2, true); // FIXME: calculate actual resolution (if possible?)
2570 }
2571 if (contributorsGID[GIndex::Source::TRD].isIndexSet()) { // ITS-TPC-TRD-TOF, TPC-TRD-TOF, TPC-TRD, ITS-TPC-TRD
2572 const auto& trdOrig = data.getTrack<o2::trd::TrackTRD>(contributorsGID[GIndex::Source::TRD]); // refitted TRD trac
2573 extraInfoHolder.trdChi2 = trdOrig.getChi2();
2574 extraInfoHolder.trdSignal = trdOrig.getSignal();
2575 extraInfoHolder.trdPattern = getTRDPattern(trdOrig);
2576 if (extraInfoHolder.trackTimeRes < 0.) { // time is not set yet, this is possible only for TPC-TRD and ITS-TPC-TRD tracks, since those with TOF are set upstream
2577 // TRD is triggered: time uncertainty is within a BC
2578 const auto& trdTrig = (src == GIndex::Source::TPCTRD) ? data.getTPCTRDTriggers()[mTPCTRDTriggers[trackIndex.getIndex()]] : data.getITSTPCTRDTriggers()[mITSTPCTRDTriggers[trackIndex.getIndex()]];
2579 double ttrig = trdTrig.getBCData().differenceInBC(mStartIR) * o2::constants::lhc::LHCBunchSpacingNS; // 1st get time wrt TF start
2580 setTrackTime(ttrig, 1., true); // FIXME: calculate actual resolution (if possible?)
2581 }
2582 }
2583 if (contributorsGID[GIndex::Source::ITS].isIndexSet()) {
2584 const auto& itsTrack = data.getITSTrack(contributorsGID[GIndex::ITS]);
2585 int nClusters = itsTrack.getNClusters();
2586 float chi2 = itsTrack.getChi2();
2587 extraInfoHolder.itsChi2NCl = nClusters != 0 ? chi2 / (float)nClusters : 0;
2588 extraInfoHolder.itsClusterSizes = itsTrack.getClusterSizes();
2589 if (src == GIndex::ITS) { // standalone ITS track should set its time from the ROF
2590 const auto& rof = data.getITSTracksROFRecords()[mITSROFs[trackIndex.getIndex()]];
2591 double t = rof.getBCData().differenceInBC(mStartIR) * o2::constants::lhc::LHCBunchSpacingNS + mITSROFrameHalfLengthNS + mITSROFBiasNS;
2592 setTrackTime(t, mITSROFrameHalfLengthNS, false);
2593 }
2594 } else if (contributorsGID[GIndex::Source::ITSAB].isIndexSet()) { // this is an ITS-TPC afterburner contributor
2595 extraInfoHolder.itsClusterSizes = data.getITSABRefs()[contributorsGID[GIndex::Source::ITSAB].getIndex()].getClusterSizes();
2596 }
2597 if (contributorsGID[GIndex::Source::TPC].isIndexSet()) {
2598 const auto& tpcOrig = data.getTPCTrack(contributorsGID[GIndex::TPC]);
2599 const auto& tpcClData = mTPCCounters[contributorsGID[GIndex::TPC]];
2600 const auto& dEdx = tpcOrig.getdEdx().dEdxTotTPC > 0 ? tpcOrig.getdEdx() : tpcOrig.getdEdxAlt();
2601 if (tpcOrig.getdEdx().dEdxTotTPC == 0) {
2602 extraInfoHolder.flags |= o2::aod::track::TPCdEdxAlt;
2603 }
2604 extraInfoHolder.tpcInnerParam = tpcOrig.getP() / tpcOrig.getAbsCharge();
2605 extraInfoHolder.tpcChi2NCl = tpcOrig.getNClusters() ? tpcOrig.getChi2() / tpcOrig.getNClusters() : 0;
2606 extraInfoHolder.tpcSignal = dEdx.dEdxTotTPC;
2607 extraInfoHolder.tpcNClsFindable = tpcOrig.getNClusters();
2608 extraInfoHolder.tpcNClsFindableMinusFound = tpcOrig.getNClusters() - tpcClData.found;
2609 extraInfoHolder.tpcNClsFindableMinusCrossedRows = tpcOrig.getNClusters() - tpcClData.crossed;
2610 extraInfoHolder.tpcNClsShared = tpcClData.shared;
2611 uint32_t clsUsedForPID = dEdx.NHitsIROC + dEdx.NHitsOROC1 + dEdx.NHitsOROC2 + dEdx.NHitsOROC3;
2612 extraInfoHolder.tpcNClsFindableMinusPID = tpcOrig.getNClusters() - clsUsedForPID;
2613 if (src == GIndex::TPC) { // standalone TPC track should set its time from their timebins range
2614 if (needBCSlice) {
2615 double t = (tpcOrig.getTime0() + 0.5 * (tpcOrig.getDeltaTFwd() - tpcOrig.getDeltaTBwd())) * mTPCBinNS; // central value
2616 double terr = 0.5 * (tpcOrig.getDeltaTFwd() + tpcOrig.getDeltaTBwd()) * mTPCBinNS;
2617 double err = mTimeMarginTrackTime + terr;
2618 bcOfTimeRef = fillBCSlice(extraInfoHolder.bcSlice, t - err, t + err, bcsMap);
2619 }
2621 p.setDeltaTFwd(tpcOrig.getDeltaTFwd());
2622 p.setDeltaTBwd(tpcOrig.getDeltaTBwd());
2623 extraInfoHolder.trackTimeRes = p.getTimeErr();
2624 extraInfoHolder.trackTime = float(tpcOrig.getTime0() * mTPCBinNS - bcOfTimeRef * o2::constants::lhc::LHCBunchSpacingNS);
2625 extraInfoHolder.diffBCRef = int(bcOfTimeRef);
2626 extraInfoHolder.isTPConly = true; // no truncation
2627 extraInfoHolder.flags |= o2::aod::track::TrackTimeAsym;
2628 } else if (src == GIndex::ITSTPC) { // its-tpc matched tracks have gaussian time error and the time was not set above
2629 const auto& trITSTPC = data.getTPCITSTrack(trackIndex);
2630 auto ts = trITSTPC.getTimeMUS();
2631 setTrackTime(ts.getTimeStamp() * 1.e3, ts.getTimeStampError() * 1.e3, true);
2632 }
2633 }
2634
2635 extrapolateToCalorimeters(extraInfoHolder, data.getTrackParamOut(trackIndex));
2636 // set bit encoding for PVContributor property as part of the flag field
2637 if (trackIndex.isPVContributor()) {
2638 extraInfoHolder.flags |= o2::aod::track::PVContributor;
2639 }
2640 return extraInfoHolder;
2641}
2642
2643AODProducerWorkflowDPL::TrackQA AODProducerWorkflowDPL::processBarrelTrackQA(int collisionID, std::uint64_t collisionBC, GIndex trackIndex,
2644 const o2::globaltracking::RecoContainer& data, const std::map<uint64_t, int>& bcsMap)
2645{
2646 TrackQA trackQAHolder;
2647 auto contributorsGID = data.getTPCContributorGID(trackIndex);
2648 const auto& trackPar = data.getTrackParam(trackIndex);
2649 if (contributorsGID.isIndexSet()) {
2650 auto prop = o2::base::Propagator::Instance();
2651 const auto& tpcOrig = data.getTPCTrack(contributorsGID);
2654 const o2::base::Propagator::MatCorrType mMatType = o2::base::Propagator::MatCorrType::USEMatCorrLUT;
2655 const o2::dataformats::VertexBase v = mVtx.getMeanVertex(collisionID < 0 ? 0.f : data.getPrimaryVertex(collisionID).getZ());
2656 std::array<float, 2> dcaInfo{-999., -999.};
2657 if (prop->propagateToDCABxByBz({v.getX(), v.getY(), v.getZ()}, tpcTMP, 2.f, mMatType, &dcaInfo)) {
2658 trackQAHolder.tpcdcaR = 100. * dcaInfo[0] / sqrt(1. + trackPar.getQ2Pt() * trackPar.getQ2Pt());
2659 trackQAHolder.tpcdcaZ = 100. * dcaInfo[1] / sqrt(1. + trackPar.getQ2Pt() * trackPar.getQ2Pt());
2660 }
2661 // This allows to safely clamp any float to one byte, using the
2662 // minmal/maximum values as under-/overflow borders and rounding to the nearest integer
2663 auto safeInt8Clamp = [](auto value) -> int8_t {
2664 using ValType = decltype(value);
2665 return static_cast<int8_t>(TMath::Nint(std::clamp(value, static_cast<ValType>(std::numeric_limits<int8_t>::min()), static_cast<ValType>(std::numeric_limits<int8_t>::max()))));
2666 };
2667 auto safeUInt8Clamp = [](auto value) -> uint8_t {
2668 using ValType = decltype(value);
2669 return static_cast<uint8_t>(TMath::Nint(std::clamp(value, static_cast<ValType>(std::numeric_limits<uint8_t>::min()), static_cast<ValType>(std::numeric_limits<uint8_t>::max()))));
2670 };
2671
2673 uint8_t clusterCounters[8] = {0};
2674 {
2675 uint8_t sectorIndex, rowIndex;
2676 uint32_t clusterIndex;
2677 const auto& tpcClusRefs = data.getTPCTracksClusterRefs();
2678 for (int i = 0; i < tpcOrig.getNClusterReferences(); i++) {
2679 o2::tpc::TrackTPC::getClusterReference(tpcClusRefs, i, sectorIndex, rowIndex, clusterIndex, tpcOrig.getClusterRef());
2680 char indexTracklet = (rowIndex % 152) / 19;
2681 clusterCounters[indexTracklet]++;
2682 }
2683 }
2684 uint8_t byteMask = 0;
2685 for (int i = 0; i < 8; i++) {
2686 if (clusterCounters[i] > 5) {
2687 byteMask |= (1 << i);
2688 }
2689 }
2690 trackQAHolder.tpcTime0 = tpcOrig.getTime0();
2691 trackQAHolder.tpcClusterByteMask = byteMask;
2692 const auto& dEdxInfoAlt = tpcOrig.getdEdxAlt(); // tpcOrig.getdEdx()
2693 const float dEdxNorm = (dEdxInfoAlt.dEdxTotTPC > 0) ? 100. / dEdxInfoAlt.dEdxTotTPC : 0;
2694 trackQAHolder.tpcdEdxNorm = dEdxInfoAlt.dEdxTotTPC;
2695 trackQAHolder.tpcdEdxMax0R = safeUInt8Clamp(dEdxInfoAlt.dEdxMaxIROC * dEdxNorm);
2696 trackQAHolder.tpcdEdxMax1R = safeUInt8Clamp(dEdxInfoAlt.dEdxMaxOROC1 * dEdxNorm);
2697 trackQAHolder.tpcdEdxMax2R = safeUInt8Clamp(dEdxInfoAlt.dEdxMaxOROC2 * dEdxNorm);
2698 trackQAHolder.tpcdEdxMax3R = safeUInt8Clamp(dEdxInfoAlt.dEdxMaxOROC3 * dEdxNorm);
2699 //
2700 trackQAHolder.tpcdEdxTot0R = safeUInt8Clamp(dEdxInfoAlt.dEdxTotIROC * dEdxNorm);
2701 trackQAHolder.tpcdEdxTot1R = safeUInt8Clamp(dEdxInfoAlt.dEdxTotOROC1 * dEdxNorm);
2702 trackQAHolder.tpcdEdxTot2R = safeUInt8Clamp(dEdxInfoAlt.dEdxTotOROC2 * dEdxNorm);
2703 trackQAHolder.tpcdEdxTot3R = safeUInt8Clamp(dEdxInfoAlt.dEdxTotOROC3 * dEdxNorm);
2705 float scaleTOF{0};
2706 auto contributorsGIDA = data.getSingleDetectorRefs(trackIndex);
2707 if (contributorsGIDA[GIndex::Source::TOF].isIndexSet()) { // ITS-TPC-TRD-TOF, ITS-TPC-TOF, TPC-TRD-TOF, TPC-TOF
2708 const auto& tofMatch = data.getTOFMatch(trackIndex);
2709 const float qpt = trackPar.getQ2Pt();
2711 trackQAHolder.dTofdX = safeInt8Clamp(tofMatch.getDXatTOF() / scaleTOF);
2712 trackQAHolder.dTofdZ = safeInt8Clamp(tofMatch.getDZatTOF() / scaleTOF);
2713 }
2714
2715 // Add matching information at a reference point (defined by
2716 // o2::aod::track::trackQARefRadius) in the same frame as the global track
2717 // without material corrections and error propagation
2718 if (auto itsContGID = data.getITSContributorGID(trackIndex); itsContGID.isIndexSet() && itsContGID.getSource() != GIndex::ITSAB) {
2719 const auto& itsOrig = data.getITSTrack(itsContGID);
2720 o2::track::TrackPar gloCopy = trackPar;
2721 o2::track::TrackPar itsCopy = itsOrig.getParamOut();
2722 o2::track::TrackPar tpcCopy = tpcOrig;
2723 if (prop->propagateToX(gloCopy, o2::aod::track::trackQARefRadius, prop->getNominalBz(), o2::base::Propagator::MAX_SIN_PHI, o2::base::Propagator::MAX_STEP, mMatCorr) &&
2724 prop->propagateToAlphaX(tpcCopy, gloCopy.getAlpha(), o2::aod::track::trackQARefRadius, false, o2::base::Propagator::MAX_SIN_PHI, o2::base::Propagator::MAX_STEP, 1, mMatCorr) &&
2725 prop->propagateToAlphaX(itsCopy, gloCopy.getAlpha(), o2::aod::track::trackQARefRadius, false, o2::base::Propagator::MAX_SIN_PHI, o2::base::Propagator::MAX_STEP, 1, mMatCorr)) {
2726 // All tracks are now at the same radius and in the same frame and we can calculate the deltas wrt. to the global track
2727 // The scale is defined by the global track scaling depending on beta0
2728 const float beta0 = std::sqrt(std::min(50.f / tpcOrig.getdEdx().dEdxMaxTPC, 1.f));
2729 const float qpt = gloCopy.getQ2Pt();
2730 const float x = qpt / beta0;
2731 // scaling is defined as sigmaBins/sqrt(p0^2 + (p1 * q/pt / beta)^2)
2732 auto scaleCont = [&x](int i) -> float {
2734 };
2735 auto scaleGlo = [&x](int i) -> float {
2737 };
2738
2739 // Calculate deltas for contributors
2740 trackQAHolder.dRefContY = safeInt8Clamp((itsCopy.getY() - tpcCopy.getY()) * scaleCont(0));
2741 trackQAHolder.dRefContZ = safeInt8Clamp((itsCopy.getZ() - tpcCopy.getZ()) * scaleCont(1));
2742 trackQAHolder.dRefContSnp = safeInt8Clamp((itsCopy.getSnp() - tpcCopy.getSnp()) * scaleCont(2));
2743 trackQAHolder.dRefContTgl = safeInt8Clamp((itsCopy.getTgl() - tpcCopy.getTgl()) * scaleCont(3));
2744 trackQAHolder.dRefContQ2Pt = safeInt8Clamp((itsCopy.getQ2Pt() - tpcCopy.getQ2Pt()) * scaleCont(4));
2745 // Calculate deltas for global track against averaged contributors
2746 trackQAHolder.dRefGloY = safeInt8Clamp(((itsCopy.getY() + tpcCopy.getY()) * 0.5f - gloCopy.getY()) * scaleGlo(0));
2747 trackQAHolder.dRefGloZ = safeInt8Clamp(((itsCopy.getZ() + tpcCopy.getZ()) * 0.5f - gloCopy.getZ()) * scaleGlo(1));
2748 trackQAHolder.dRefGloSnp = safeInt8Clamp(((itsCopy.getSnp() + tpcCopy.getSnp()) * 0.5f - gloCopy.getSnp()) * scaleGlo(2));
2749 trackQAHolder.dRefGloTgl = safeInt8Clamp(((itsCopy.getTgl() + tpcCopy.getTgl()) * 0.5f - gloCopy.getTgl()) * scaleGlo(3));
2750 trackQAHolder.dRefGloQ2Pt = safeInt8Clamp(((itsCopy.getQ2Pt() + tpcCopy.getQ2Pt()) * 0.5f - gloCopy.getQ2Pt()) * scaleGlo(4));
2751 //
2752
2753 if (mStreamerFlags[AODProducerStreamerFlags::TrackQA]) {
2754 (*mStreamer) << "trackQA"
2755 << "trackITSOrig=" << itsOrig
2756 << "trackTPCOrig=" << tpcOrig
2757 << "trackITSTPCOrig=" << trackPar
2758 << "trackITSProp=" << itsCopy
2759 << "trackTPCProp=" << tpcCopy
2760 << "trackITSTPCProp=" << gloCopy
2761 << "refRadius=" << o2::aod::track::trackQARefRadius
2762 << "scaleBins=" << o2::aod::track::trackQAScaleBins
2763 << "scaleCont0=" << scaleCont(0)
2764 << "scaleCont1=" << scaleCont(1)
2765 << "scaleCont2=" << scaleCont(2)
2766 << "scaleCont3=" << scaleCont(3)
2767 << "scaleCont4=" << scaleCont(4)
2768 << "scaleGlo0=" << scaleGlo(0)
2769 << "scaleGlo1=" << scaleGlo(1)
2770 << "scaleGlo2=" << scaleGlo(2)
2771 << "scaleGlo3=" << scaleGlo(3)
2772 << "scaleGlo4=" << scaleGlo(4)
2773 << "trackQAHolder.tpcTime0=" << trackQAHolder.tpcTime0
2774 << "trackQAHolder.tpcdEdxNorm=" << trackQAHolder.tpcdEdxNorm
2775 << "trackQAHolder.tpcdcaR=" << trackQAHolder.tpcdcaR
2776 << "trackQAHolder.tpcdcaZ=" << trackQAHolder.tpcdcaZ
2777 << "trackQAHolder.tpcdcaClusterByteMask=" << trackQAHolder.tpcClusterByteMask
2778 << "trackQAHolder.tpcdEdxMax0R=" << trackQAHolder.tpcdEdxMax0R
2779 << "trackQAHolder.tpcdEdxMax1R=" << trackQAHolder.tpcdEdxMax1R
2780 << "trackQAHolder.tpcdEdxMax2R=" << trackQAHolder.tpcdEdxMax2R
2781 << "trackQAHolder.tpcdEdxMax3R=" << trackQAHolder.tpcdEdxMax3R
2782 << "trackQAHolder.tpcdEdxTot0R=" << trackQAHolder.tpcdEdxTot0R
2783 << "trackQAHolder.tpcdEdxTot1R=" << trackQAHolder.tpcdEdxTot1R
2784 << "trackQAHolder.tpcdEdxTot2R=" << trackQAHolder.tpcdEdxTot2R
2785 << "trackQAHolder.tpcdEdxTot3R=" << trackQAHolder.tpcdEdxTot3R
2786 << "trackQAHolder.dRefContY=" << trackQAHolder.dRefContY
2787 << "trackQAHolder.dRefContZ=" << trackQAHolder.dRefContZ
2788 << "trackQAHolder.dRefContSnp=" << trackQAHolder.dRefContSnp
2789 << "trackQAHolder.dRefContTgl=" << trackQAHolder.dRefContTgl
2790 << "trackQAHolder.dRefContQ2Pt=" << trackQAHolder.dRefContQ2Pt
2791 << "trackQAHolder.dRefGloY=" << trackQAHolder.dRefGloY
2792 << "trackQAHolder.dRefGloZ=" << trackQAHolder.dRefGloZ
2793 << "trackQAHolder.dRefGloSnp=" << trackQAHolder.dRefGloSnp
2794 << "trackQAHolder.dRefGloTgl=" << trackQAHolder.dRefGloTgl
2795 << "trackQAHolder.dRefGloQ2Pt=" << trackQAHolder.dRefGloQ2Pt
2796 << "trackQAHolder.dTofdX=" << trackQAHolder.dTofdX
2797 << "trackQAHolder.dTofdZ=" << trackQAHolder.dTofdZ
2798 << "scaleTOF=" << scaleTOF
2799 << "\n";
2800 }
2801 }
2802 }
2803 }
2804
2805 return trackQAHolder;
2806}
2807
2808bool AODProducerWorkflowDPL::propagateTrackToPV(o2::track::TrackParametrizationWithError<float>& trackPar,
2810 int colID)
2811{
2812 o2::dataformats::DCA dcaInfo;
2813 dcaInfo.set(999.f, 999.f, 999.f, 999.f, 999.f);
2814 o2::dataformats::VertexBase v = mVtx.getMeanVertex(colID < 0 ? 0.f : data.getPrimaryVertex(colID).getZ());
2815 return o2::base::Propagator::Instance()->propagateToDCABxByBz(v, trackPar, 2.f, mMatCorr, &dcaInfo);
2816}
2817
2818void AODProducerWorkflowDPL::extrapolateToCalorimeters(TrackExtraInfo& extraInfoHolder, const o2::track::TrackPar& track)
2819{
2820 constexpr float XEMCAL = 440.f, XPHOS = 460.f, XEMCAL2 = XEMCAL * XEMCAL;
2821 constexpr float ETAEMCAL = 0.75; // eta of EMCAL/DCAL with margin
2822 constexpr float ZEMCALFastCheck = 460.; // Max Z (with margin to check with straightline extrapolarion)
2823 constexpr float ETADCALINNER = 0.22; // eta of the DCAL PHOS Hole (at XEMCAL)
2824 constexpr float ETAPHOS = 0.13653194; // nominal eta of the PHOS acceptance (at XPHOS): -log(tan((TMath::Pi()/2 - atan2(63, 460))/2))
2825 constexpr float ETAPHOSMARGIN = 0.17946979; // etat of the PHOS acceptance with 20 cm margin (at XPHOS): -log(tan((TMath::Pi()/2 + atan2(63+20., 460))/2)), not used, for the ref only
2826 constexpr float ETADCALPHOSSWITCH = (ETADCALINNER + ETAPHOS) / 2; // switch to DCAL to PHOS check if eta < this value
2827 constexpr short SNONE = 0, SEMCAL = 0x1, SPHOS = 0x2;
2828 constexpr short SECTORTYPE[18] = {
2829 SNONE, SNONE, SNONE, SNONE, // 0:3
2830 SEMCAL, SEMCAL, SEMCAL, SEMCAL, SEMCAL, SEMCAL, // 3:9 EMCAL only
2831 SNONE, SNONE, // 10:11
2832 SPHOS, // 12 PHOS only
2833 SPHOS | SEMCAL, SPHOS | SEMCAL, SPHOS | SEMCAL, // 13:15 PHOS & DCAL
2834 SEMCAL, // 16 DCAL only
2835 SNONE // 17
2836 };
2837
2838 o2::track::TrackPar outTr{track};
2839 auto prop = o2::base::Propagator::Instance();
2840 // 1st propagate to EMCAL nominal radius
2841 float xtrg = 0;
2842 // quick check with straight line propagtion
2843 if (!outTr.getXatLabR(XEMCAL, xtrg, prop->getNominalBz(), o2::track::DirType::DirOutward) ||
2844 (std::abs(outTr.getZAt(xtrg, 0)) > ZEMCALFastCheck) ||
2845 !prop->PropagateToXBxByBz(outTr, xtrg, 0.95, 10, o2::base::Propagator::MatCorrType::USEMatCorrLUT)) {
2846 LOGP(debug, "preliminary step: does not reach R={} {}", XEMCAL, outTr.asString());
2847 return;
2848 }
2849 // we do not necessarilly reach wanted radius in a single propagation
2850 if ((outTr.getX() * outTr.getX() + outTr.getY() * outTr.getY() < XEMCAL2) &&
2851 (!outTr.rotateParam(outTr.getPhi()) ||
2852 !outTr.getXatLabR(XEMCAL, xtrg, prop->getNominalBz(), o2::track::DirType::DirOutward) ||
2853 !prop->PropagateToXBxByBz(outTr, xtrg, 0.95, 10, o2::base::Propagator::MatCorrType::USEMatCorrLUT))) {
2854 LOGP(debug, "does not reach R={} {}", XEMCAL, outTr.asString());
2855 return;
2856 }
2857 // rotate to proper sector
2858 int sector = o2::math_utils::angle2Sector(outTr.getPhiPos());
2859
2860 auto propExactSector = [&outTr, &sector, prop](float xprop) -> bool { // propagate exactly to xprop in the proper sector frame
2861 int ntri = 0;
2862 while (ntri < 2) {
2863 auto outTrTmp = outTr;
2864 float alpha = o2::math_utils::sector2Angle(sector);
2865 if ((std::abs(outTr.getZ()) > ZEMCALFastCheck) || !outTrTmp.rotateParam(alpha) ||
2866 !prop->PropagateToXBxByBz(outTrTmp, xprop, 0.95, 10, o2::base::Propagator::MatCorrType::USEMatCorrLUT)) {
2867 LOGP(debug, "failed on rotation to {} (sector {}) or propagation to X={} {}", alpha, sector, xprop, outTrTmp.asString());
2868 return false;
2869 }
2870 // make sure we are still in the target sector
2871 int sectorTmp = o2::math_utils::angle2Sector(outTrTmp.getPhiPos());
2872 if (sectorTmp == sector) {
2873 outTr = outTrTmp;
2874 break;
2875 }
2876 sector = sectorTmp;
2877 ntri++;
2878 }
2879 if (ntri == 2) {
2880 LOGP(debug, "failed to rotate to sector, {}", outTr.asString());
2881 return false;
2882 }
2883 return true;
2884 };
2885
2886 // we are at the EMCAL X, check if we are in the good sector
2887 if (!propExactSector(XEMCAL) || SECTORTYPE[sector] == SNONE) { // propagation failed or neither EMCAL not DCAL/PHOS
2888 return;
2889 }
2890
2891 // check if we are in a good eta range
2892 float r = std::sqrt(outTr.getX() * outTr.getX() + outTr.getY() * outTr.getY()), tg = std::atan2(r, outTr.getZ());
2893 float eta = -std::log(std::tan(0.5f * tg)), etaAbs = std::abs(eta);
2894 if (etaAbs > ETAEMCAL) {
2895 LOGP(debug, "eta = {} is off at EMCAL radius", eta, outTr.asString());
2896 return;
2897 }
2898 // are we in the PHOS hole (with margin)?
2899 if ((SECTORTYPE[sector] & SPHOS) && etaAbs < ETADCALPHOSSWITCH) { // propagate to PHOS radius
2900 if (!propExactSector(XPHOS)) {
2901 return;
2902 }
2903 r = std::sqrt(outTr.getX() * outTr.getX() + outTr.getY() * outTr.getY());
2904 tg = std::atan2(r, outTr.getZ());
2905 eta = -std::log(std::tan(0.5f * tg));
2906 } else if (!(SECTORTYPE[sector] & SEMCAL)) { // are in the sector with PHOS only
2907 return;
2908 }
2909 extraInfoHolder.trackPhiEMCAL = outTr.getPhiPos();
2910 extraInfoHolder.trackEtaEMCAL = eta;
2911 LOGP(debug, "eta = {} phi = {} sector {} for {}", extraInfoHolder.trackEtaEMCAL, extraInfoHolder.trackPhiEMCAL, sector, outTr.asString());
2912 //
2913}
2914
2915std::set<uint64_t> AODProducerWorkflowDPL::filterEMCALIncomplete(const gsl::span<const o2::emcal::TriggerRecord> triggers)
2916{
2917 std::set<uint64_t> emcalIncompletes;
2918 for (const auto& trg : triggers) {
2919 if (trg.getTriggerBits() & o2::emcal::triggerbits::Inc) {
2920 // trigger record masked at incomplete at readout level
2921 emcalIncompletes.insert(trg.getBCData().toLong());
2922 }
2923 }
2924 return emcalIncompletes;
2925}
2926
2927void AODProducerWorkflowDPL::updateTimeDependentParams(ProcessingContext& pc)
2928{
2930 static bool initOnceDone = false;
2931 if (!initOnceDone) { // this params need to be queried only once
2932 initOnceDone = true;
2933 // Note: DPLAlpideParam for ITS and MFT will be loaded by the RecoContainer
2934 mSqrtS = o2::base::GRPGeomHelper::instance().getGRPLHCIF()->getSqrtS();
2935 // apply settings
2938 std::bitset<3564> bs = bcf.getBCPattern();
2939 for (auto i = 0U; i < bs.size(); i++) {
2940 if (bs.test(i)) {
2942 }
2943 }
2944
2946 mITSROFrameHalfLengthNS = 0.5 * (grpECS->isDetContinuousReadOut(o2::detectors::DetID::ITS) ? alpParamsITS.roFrameLengthInBC * o2::constants::lhc::LHCBunchSpacingNS : alpParamsITS.roFrameLengthTrig);
2949 mMFTROFrameHalfLengthNS = 0.5 * (grpECS->isDetContinuousReadOut(o2::detectors::DetID::MFT) ? alpParamsMFT.roFrameLengthInBC * o2::constants::lhc::LHCBunchSpacingNS : alpParamsMFT.roFrameLengthTrig);
2951
2952 // RS FIXME: this is not yet fetched from the CCDB
2954 mTPCBinNS = elParam.ZbinWidth * 1.e3;
2955
2956 const auto& pvParams = o2::vertexing::PVertexerParams::Instance();
2957 mNSigmaTimeTrack = pvParams.nSigmaTimeTrack;
2958 mTimeMarginTrackTime = pvParams.timeMarginTrackTime * 1.e3;
2959 mFieldON = std::abs(o2::base::Propagator::Instance()->getNominalBz()) > 0.01;
2960
2961 pc.inputs().get<o2::ctp::CTPConfiguration*>("ctpconfig");
2962 }
2963 if (mPropTracks) {
2965 }
2966}
2967
2968//_______________________________________
2970{
2971 // Note: strictly speaking, for Configurable params we don't need finaliseCCDB check, the singletons are updated at the CCDB fetcher level
2972 if (o2::base::GRPGeomHelper::instance().finaliseCCDB(matcher, obj)) {
2973 if (matcher == ConcreteDataMatcher("GLO", "GRPMAGFIELD", 0)) {
2975 }
2976 return;
2977 }
2978 if (matcher == ConcreteDataMatcher("ITS", "ALPIDEPARAM", 0)) {
2979 LOG(info) << "ITS Alpide param updated";
2981 par.printKeyValues();
2982 return;
2983 }
2984 if (matcher == ConcreteDataMatcher("MFT", "ALPIDEPARAM", 0)) {
2985 LOG(info) << "MFT Alpide param updated";
2987 par.printKeyValues();
2988 return;
2989 }
2990 if (matcher == ConcreteDataMatcher("GLO", "MEANVERTEX", 0)) {
2991 LOG(info) << "Imposing new MeanVertex: " << ((const o2::dataformats::MeanVertexObject*)obj)->asString();
2992 mVtx = *(const o2::dataformats::MeanVertexObject*)obj;
2993 return;
2994 }
2995 if (matcher == ConcreteDataMatcher("CTP", "CTPCONFIG", 0)) {
2996 // construct mask with EMCAL trigger classes for rejection of incomplete triggers
2997 auto ctpconfig = *(const o2::ctp::CTPConfiguration*)obj;
2998 mEMCALTrgClassMask = 0;
2999 for (const auto& trgclass : ctpconfig.getCTPClasses()) {
3000 if (trgclass.cluster->maskCluster[o2::detectors::DetID::EMC]) {
3001 mEMCALTrgClassMask |= trgclass.classMask;
3002 }
3003 }
3004 LOG(info) << "Loaded EMCAL trigger class mask: " << std::bitset<64>(mEMCALTrgClassMask);
3005 }
3006}
3007
3008void AODProducerWorkflowDPL::addRefGlobalBCsForTOF(const o2::dataformats::VtxTrackRef& trackRef, const gsl::span<const GIndex>& GIndices,
3009 const o2::globaltracking::RecoContainer& data, std::map<uint64_t, int>& bcsMap)
3010{
3011 // Orphan tracks need to refer to some globalBC and for tracks with TOF this BC should be whithin an orbit
3012 // from the track abs time (to guarantee time precision). Therefore, we may need to insert some dummy globalBCs
3013 // to guarantee proper reference.
3014 // complete globalBCs by dummy entries necessary to provide BC references for TOF tracks with requested precision
3015 // to provide a reference for the time of the orphan tracks we should make sure that there are no gaps longer
3016 // than needed to store the time with sufficient precision
3017
3018 // estimate max distance in BCs between TOF time and eventual reference BC
3019 int nbitsFrac = 24 - (32 - o2::math_utils::popcount(mTrackTime)); // number of bits used to encode the fractional part of float truncated by the mask
3020 int nbitsLoss = std::max(0, int(std::log2(TOFTimePrecPS))); // allowed bit loss guaranteeing needed precision in PS
3021 assert(nbitsFrac > 1);
3022 std::uint64_t maxRangePS = std::uint64_t(0x1) << (nbitsFrac + nbitsLoss);
3023 int maxGapBC = maxRangePS / (o2::constants::lhc::LHCBunchSpacingNS * 1e3); // max gap in BCs allowing to store time with required precision
3024 LOG(info) << "Max gap of " << maxGapBC << " BCs to closest globalBC reference is needed for TOF tracks to provide precision of "
3025 << TOFTimePrecPS << " ps";
3026
3027 // check if there are tracks orphan tracks at all
3028 if (!trackRef.getEntries()) {
3029 return;
3030 }
3031 // the bscMap has at least TF start BC
3032 std::uint64_t maxBC = mStartIR.toLong();
3033 const auto& tofClus = data.getTOFClusters();
3034 for (int src = GIndex::NSources; src--;) {
3035 if (!GIndex::getSourceDetectorsMask(src)[o2::detectors::DetID::TOF]) { // check only tracks with TOF contribution
3036 continue;
3037 }
3038 int start = trackRef.getFirstEntryOfSource(src);
3039 int end = start + trackRef.getEntriesOfSource(src);
3040 for (int ti = start; ti < end; ti++) {
3041 auto& trackIndex = GIndices[ti];
3042 const auto& tofMatch = data.getTOFMatch(trackIndex);
3043 const auto& tofInt = tofMatch.getLTIntegralOut();
3044 float intLen = tofInt.getL();
3045 float tofExpMom = 0.;
3046 if (tofInt.getTOF(o2::track::PID::Pion) > 0.f) {
3047 float expBeta = (intLen / (tofInt.getTOF(o2::track::PID::Pion) * cSpeed));
3048 if (expBeta > o2::constants::math::Almost1) {
3050 }
3051 tofExpMom = o2::constants::physics::MassPionCharged * expBeta / std::sqrt(1.f - expBeta * expBeta);
3052 } else {
3053 continue;
3054 }
3055 double massZ = o2::track::PID::getMass2Z(data.getTrackParam(trackIndex).getPID());
3056 double energy = sqrt((massZ * massZ) + (tofExpMom * tofExpMom));
3057 double exp = intLen * energy / (cSpeed * tofExpMom);
3058 auto tofSignal = (tofMatch.getSignal() - exp) * 1e-3; // time in ns wrt TF start
3059 auto bc = relativeTime_to_GlobalBC(tofSignal);
3060
3061 auto it = bcsMap.lower_bound(bc);
3062 if (it == bcsMap.end() || it->first > bc + maxGapBC) {
3063 bcsMap.emplace_hint(it, bc, 1);
3064 LOG(debug) << "adding dummy BC " << bc;
3065 }
3066 if (bc > maxBC) {
3067 maxBC = bc;
3068 }
3069 }
3070 }
3071 // make sure there is a globalBC exceeding the max encountered bc
3072 if ((--bcsMap.end())->first <= maxBC) {
3073 bcsMap.emplace_hint(bcsMap.end(), maxBC + 1, 1);
3074 }
3075 // renumber BCs
3076 int bcID = 0;
3077 for (auto& item : bcsMap) {
3078 item.second = bcID;
3079 bcID++;
3080 }
3081}
3082
3083std::uint64_t AODProducerWorkflowDPL::fillBCSlice(int (&slice)[2], double tmin, double tmax, const std::map<uint64_t, int>& bcsMap) const
3084{
3085 // for ambiguous tracks (no or multiple vertices) we store the BC slice corresponding to track time window used for track-vertex matching,
3086 // see VertexTrackMatcher::extractTracks creator method, i.e. central time estimated +- uncertainty defined as:
3087 // 1) for tracks having a gaussian time error: PVertexerParams.nSigmaTimeTrack * trackSigma + PVertexerParams.timeMarginTrackTime
3088 // 2) for tracks having time uncertainty in a fixed time interval (TPC,ITS,MFT..): half of the interval + PVertexerParams.timeMarginTrackTime
3089 // The track time in the TrackExtraInfo is stored in ns wrt the collision BC for unambigous tracks and wrt bcSlice[0] for ambiguous ones,
3090 // with convention for errors: trackSigma in case (1) and half of the time interval for case (2) above.
3091
3092 // find indices of widest slice of global BCs in the map compatible with provided BC range. bcsMap is guaranteed to be non-empty.
3093 // We also assume that tmax >= tmin.
3094
3095 uint64_t bcMin = relativeTime_to_GlobalBC(tmin), bcMax = relativeTime_to_GlobalBC(tmax);
3096
3097 /*
3098 // brute force way of searching bcs via direct binary search in the map
3099 auto lower = bcsMap.lower_bound(bcMin), upper = bcsMap.upper_bound(bcMax);
3100
3101 if (lower == bcsMap.end()) {
3102 --lower;
3103 }
3104 if (upper != lower) {
3105 --upper;
3106 }
3107 slice[0] = std::distance(bcsMap.begin(), lower);
3108 slice[1] = std::distance(bcsMap.begin(), upper);
3109 */
3110
3111 // faster way to search in bunch crossing via the accelerated bunch crossing lookup structure
3112 auto p = mBCLookup.lower_bound(bcMin);
3113 // assuming that bcMax will be >= bcMin and close to bcMin; we can find
3114 // the upper bound quickly by lineary iterating from p.first to the point where
3115 // the time becomes larger than bcMax.
3116 // (if this is not the case we could determine it with a similar call to mBCLookup)
3117 auto& bcvector = mBCLookup.getBCTimeVector();
3118 auto upperindex = p.first;
3119 while (upperindex < bcvector.size() && bcvector[upperindex] <= bcMax) {
3120 upperindex++;
3121 }
3122 if (upperindex != p.first) {
3123 upperindex--;
3124 }
3125 slice[0] = p.first;
3126 slice[1] = upperindex;
3127
3128 auto bcOfTimeRef = p.second - this->mStartIR.toLong();
3129 LOG(debug) << "BC slice t:" << tmin << " " << slice[0]
3130 << " t: " << tmax << " " << slice[1]
3131 << " bcref: " << bcOfTimeRef;
3132 return bcOfTimeRef;
3133}
3134
3135std::vector<uint8_t> AODProducerWorkflowDPL::fillBCFlags(const o2::globaltracking::RecoContainer& data, std::map<uint64_t, int>& bcsMap) const
3136{
3137 std::vector<uint8_t> flags(bcsMap.size());
3138
3139 // flag BCs belonging to UPC mode ITS ROFs
3140 auto bcIt = bcsMap.begin();
3141 auto itsrofs = data.getITSTracksROFRecords();
3144 for (auto& rof : itsrofs) {
3145 if (!rof.getFlag(o2::itsmft::ROFRecord::VtxUPCMode)) {
3146 continue;
3147 }
3148 uint64_t globalBC0 = rof.getBCData().toLong() + bROF, globalBC1 = globalBC0 + lROF - 1;
3149 // BCs are sorted, iterate until the start of ROF
3150 while (bcIt != bcsMap.end()) {
3151 if (bcIt->first < globalBC0) {
3152 ++bcIt;
3153 continue;
3154 }
3155 if (bcIt->first > globalBC1) {
3156 break;
3157 }
3158 flags[bcIt->second] |= o2::aod::bc::ITSUPCMode;
3159 ++bcIt;
3160 }
3161 }
3162 return flags;
3163}
3164
3166{
3167 LOGF(info, "aod producer dpl total timing: Cpu: %.3e Real: %.3e s in %d slots",
3168 mTimer.CpuTime(), mTimer.RealTime(), mTimer.Counter() - 1);
3169
3170 mStreamer.reset();
3171}
3172
3173DataProcessorSpec getAODProducerWorkflowSpec(GID::mask_t src, bool enableSV, bool enableStrangenessTracking, bool useMC, bool CTPConfigPerRun, bool enableFITextra)
3174{
3175 auto dataRequest = std::make_shared<DataRequest>();
3176 dataRequest->inputs.emplace_back("ctpconfig", "CTP", "CTPCONFIG", 0, Lifetime::Condition, ccdbParamSpec("CTP/Config/Config", CTPConfigPerRun));
3177
3178 dataRequest->requestTracks(src, useMC);
3179 dataRequest->requestPrimaryVertices(useMC);
3180 if (src[GID::CTP]) {
3181 dataRequest->requestCTPDigits(useMC);
3182 }
3183 if (enableSV) {
3184 dataRequest->requestSecondaryVertices(useMC);
3185 }
3186 if (enableStrangenessTracking) {
3187 dataRequest->requestStrangeTracks(useMC);
3188 LOGF(info, "requestStrangeTracks Finish");
3189 }
3190 if (src[GID::ITS]) {
3191 dataRequest->requestClusters(GIndex::getSourcesMask("ITS"), false);
3192 }
3193 if (src[GID::TPC]) {
3194 dataRequest->requestClusters(GIndex::getSourcesMask("TPC"), false); // no need to ask for TOF clusters as they are requested with TOF tracks
3195 }
3196 if (src[GID::TOF]) {
3197 dataRequest->requestTOFClusters(useMC);
3198 }
3199 if (src[GID::PHS]) {
3200 dataRequest->requestPHOSCells(useMC);
3201 }
3202 if (src[GID::TRD]) {
3203 dataRequest->requestTRDTracklets(false);
3204 }
3205 if (src[GID::EMC]) {
3206 dataRequest->requestEMCALCells(useMC);
3207 }
3208 if (src[GID::CPV]) {
3209 dataRequest->requestCPVClusters(useMC);
3210 }
3211
3212 auto ggRequest = std::make_shared<o2::base::GRPGeomRequest>(true, // orbitResetTime
3213 true, // GRPECS=true
3214 true, // GRPLHCIF
3215 true, // GRPMagField
3216 true, // askMatLUT
3218 dataRequest->inputs,
3219 true); // query only once all objects except mag.field
3220
3221 dataRequest->inputs.emplace_back("meanvtx", "GLO", "MEANVERTEX", 0, Lifetime::Condition, ccdbParamSpec("GLO/Calib/MeanVertex", {}, 1));
3222
3223 using namespace o2::aod;
3224 using namespace o2::aodproducer;
3225
3226 std::vector<OutputSpec> outputs{
3260 OutputSpec{"TFN", "TFNumber"},
3261 OutputSpec{"TFF", "TFFilename"},
3262 OutputSpec{"AMD", "AODMetadataKeys"},
3263 OutputSpec{"AMD", "AODMetadataVals"}};
3264
3265 if (useMC) {
3266 outputs.insert(outputs.end(),
3267 {OutputForTable<McCollisions>::spec(),
3268 OutputForTable<HepMCXSections>::spec(),
3269 OutputForTable<HepMCPdfInfos>::spec(),
3270 OutputForTable<HepMCHeavyIons>::spec(),
3271 OutputForTable<McMFTTrackLabels>::spec(),
3272 OutputForTable<McFwdTrackLabels>::spec(),
3273 OutputForTable<StoredMcParticles_001>::spec(),
3274 OutputForTable<McTrackLabels>::spec(),
3275 OutputForTable<McCaloLabels_001>::spec(),
3276 // todo: use addTableToOuput helper?
3277 // currently the description is MCCOLLISLABEL, so
3278 // the name in AO2D would be O2mccollislabel
3279 // addTableToOutput<McCollisionLabels>(outputs);
3280 {OutputLabel{"McCollisionLabels"}, "AOD", "MCCOLLISIONLABEL", 0, Lifetime::Timeframe}});
3281 }
3282
3283 return DataProcessorSpec{
3284 "aod-producer-workflow",
3285 dataRequest->inputs,
3286 outputs,
3287 AlgorithmSpec{adaptFromTask<AODProducerWorkflowDPL>(src, dataRequest, ggRequest, enableSV, useMC, enableFITextra)},
3288 Options{
3289 ConfigParamSpec{"run-number", VariantType::Int64, -1L, {"The run-number. If left default we try to get it from DPL header."}},
3290 ConfigParamSpec{"aod-timeframe-id", VariantType::Int64, -1L, {"Set timeframe number"}},
3291 ConfigParamSpec{"fill-calo-cells", VariantType::Int, 1, {"Fill calo cells into cell table"}},
3292 ConfigParamSpec{"enable-truncation", VariantType::Int, 1, {"Truncation parameter: 1 -- on, != 1 -- off"}},
3293 ConfigParamSpec{"lpmp-prod-tag", VariantType::String, "", {"LPMProductionTag"}},
3294 ConfigParamSpec{"anchor-pass", VariantType::String, "", {"AnchorPassName"}},
3295 ConfigParamSpec{"anchor-prod", VariantType::String, "", {"AnchorProduction"}},
3296 ConfigParamSpec{"reco-pass", VariantType::String, "", {"RecoPassName"}},
3297 ConfigParamSpec{"created-by", VariantType::String, "", {"Who created this AO2D"}},
3298 ConfigParamSpec{"nthreads", VariantType::Int, std::max(1, int(std::thread::hardware_concurrency() / 2)), {"Number of threads"}},
3299 ConfigParamSpec{"reco-mctracks-only", VariantType::Int, 0, {"Store only reconstructed MC tracks and their mothers/daughters. 0 -- off, != 0 -- on"}},
3300 ConfigParamSpec{"ctpreadout-create", VariantType::Int, 0, {"Create CTP digits from detector readout and CTP inputs. !=1 -- off, 1 -- on"}},
3301 ConfigParamSpec{"emc-select-leading", VariantType::Bool, false, {"Flag to select if only the leading contributing particle for an EMCal cell should be stored"}},
3302 ConfigParamSpec{"propagate-tracks", VariantType::Bool, false, {"Propagate tracks (not used for secondary vertices) to IP"}},
3303 ConfigParamSpec{"propagate-tracks-max-xiu", VariantType::Float, 5.0f, {"Propagate tracks to IP if X_IU smaller than this value (and if propagate tracks enabled)"}},
3304 ConfigParamSpec{"hepmc-update", VariantType::String, "always", {"When to update HepMC Aux tables: always - force update, never - never update, all - if all keys are present, any - when any key is present (not valid yet)"}},
3305 ConfigParamSpec{"propagate-muons", VariantType::Bool, false, {"Propagate muons to IP"}},
3306 ConfigParamSpec{"thin-tracks", VariantType::Bool, false, {"Produce thinned track tables"}},
3307 ConfigParamSpec{"trackqc-fraction", VariantType::Float, float(0.1), {"Fraction of tracks to QC"}},
3308 ConfigParamSpec{"trackqc-NTrCut", VariantType::Int64, 4L, {"Minimal length of the track - in amount of tracklets"}},
3309 ConfigParamSpec{"trackqc-tpc-dca", VariantType::Float, 3.f, {"Keep TPC standalone track with this DCAxy to the PV"}},
3310 ConfigParamSpec{"trackqc-tpc-cls", VariantType::Int, 80, {"Keep TPC standalone track with this #clusters"}},
3311 ConfigParamSpec{"trackqc-tpc-pt", VariantType::Float, 0.2f, {"Keep TPC standalone track with this pt"}},
3312 ConfigParamSpec{"with-streamers", VariantType::String, "", {"Bit-mask to steer writing of intermediate streamer files"}},
3313 ConfigParamSpec{"seed", VariantType::Int, 0, {"Set seed for random generator used for sampling (0 (default) means using a random_device)"}},
3314 }};
3315}
3316
3317} // namespace o2::aodproducer
Class to refer to the reconstructed information.
Definition of the 32 Central Trigger System (CTS) Trigger Types defined in https://twiki....
General auxilliary methods.
definition of CTPConfiguration and related CTP structures
Wrapper container for different reconstructed object types.
Definition of the MCH cluster minimal structure.
Reconstructed MID track.
Base track model for the Barrel, params only, w/o covariance.
uint64_t exp(uint64_t base, uint8_t exp) noexcept
definition of CTPDigit, CTPInputDigit
int16_t charge
Definition RawEventData.h:5
uint64_t vertex
Definition RawEventData.h:9
uint64_t bc
Definition RawEventData.h:5
int16_t time
Definition RawEventData.h:4
Definition of the FIT RecPoints class.
Definition of the FV0 RecPoints class.
int32_t i
std::string getType()
Definition Utils.h:138
Header of the AggregatedRunInfo struct.
Global Forward Muon tracks.
Global index for barrel track: provides provenance (detectors combination), index in respective array...
Definition of the MCTrack class.
Definition of a container to keep Monte Carlo truth external to simulation objects.
Utility functions for MC particles.
Definition of the MCH ROFrame record.
Class to perform MFT MCH (and MID) matching.
Class to store the output of the matching to HMPID.
Class to perform TOF matching to global tracks.
Definition of the parameter class for the detector electronics.
Header to collect physics constants.
uint32_t j
Definition RawData.h:0
Definition of the FDD RecPoint class.
Definition of tools for track extrapolation.
Definition of the ITS track.
Definition of the MUON track.
Definition of the MCH track.
Definition of the MCH track parameters for internal use.
Result of refitting TPC-ITS matched track.
std::ostringstream debug
Extention of GlobalTrackID by flags relevant for verter-track association.
Referenc on track indices contributing to the vertex, with possibility chose tracks from specific sou...
Container class to store energy released in the ZDC.
Container class to store a TDC hit in a ZDC channel.
int nClusters
StringRef key
const auto & getBCPattern() const
void GetStartVertex(TVector3 &vertex) const
Definition MCTrack.h:326
void endOfStream(framework::EndOfStreamContext &ec) final
This is invoked whenever we have an EndOfStream event.
void finaliseCCDB(ConcreteDataMatcher &matcher, void *obj) final
std::pair< size_t, uint64_t > lower_bound(uint64_t timestamp) const
void init(std::map< uint64_t, int > const &bcs)
initialize this container (to be ready for lookup/search queries)
void clear()
clear/reset this container
std::vector< uint64_t > const & getBCTimeVector() const
return the sorted vector of increaing BC times
static float getAmplitude(const o2::emcal::Cell &cell)
static int16_t getLnAmplitude(const o2::emcal::Cell &)
static int8_t getTriggerBits(const o2::emcal::Cell &)
static int16_t getCellNumber(const o2::emcal::Cell &cell)
static int16_t getFastOrAbsID(const o2::emcal::Cell &)
static bool isTRU(const o2::emcal::Cell &cell)
static float getTimeStamp(const o2::emcal::Cell &cell)
void checkUpdates(o2::framework::ProcessingContext &pc)
static GRPGeomHelper & instance()
void setRequest(std::shared_ptr< GRPGeomRequest > req)
static constexpr float MAX_SIN_PHI
Definition Propagator.h:72
static constexpr float MAX_STEP
Definition Propagator.h:73
GPUd() value_type estimateLTFast(o2 static GPUd() float estimateLTIncrement(const o2 PropagatorImpl * Instance(bool uninitialized=false)
Definition Propagator.h:143
void printStream(std::ostream &stream) const
static mask_t getSourcesMask(const std::string_view srcList)
VertexBase getMeanVertex(float z) const
int getEntriesOfSource(int s) const
Definition VtxTrackRef.h:62
int getFirstEntryOfSource(int s) const
Definition VtxTrackRef.h:55
Static class with identifiers, bitmasks and names for ALICE detectors.
Definition DetID.h:58
static constexpr ID ITS
Definition DetID.h:63
static constexpr ID MFT
Definition DetID.h:71
static constexpr ID TPC
Definition DetID.h:64
static constexpr ID EMC
Definition DetID.h:69
static constexpr ID TOF
Definition DetID.h:66
Handler for EMCAL event data.
void reset()
Reset containers with empty ranges.
void setCellData(CellRange cells, TriggerRange triggers)
Setting the data at cell level.
void setCellMCTruthContainer(const o2::dataformats::MCTruthContainer< o2::emcal::MCLabel > *mclabels)
Setting the pointer for the MCTruthContainer for cells.
InteractionRecord getInteractionRecordForEvent(int eventID) const
Get the interaction record for the given event.
int getNumberOfEvents() const
Get the number of events handled by the event handler.
void snapshot(const Output &spec, T const &object)
ConfigParamRegistry const & options()
Definition InitContext.h:33
decltype(auto) get(R binding, int part=0) const
DataAllocator & outputs()
The data allocator is used to allocate memory for the output data.
InputRecord & inputs()
The inputs associated with this processing context.
ServiceRegistryRef services()
The services registry associated with this processing context.
static constexpr int NCellsA
Definition Geometry.h:52
o2::dataformats::GlobalFwdTrack MCHtoFwd(const o2::mch::TrackParam &mchTrack)
Converts mchTrack parameters to Forward coordinate system.
static bool extrapToZ(TrackParam &trackParam, double zEnd)
static void setField()
static bool extrapToVertex(TrackParam &trackParam, double xVtx, double yVtx, double zVtx, double errXVtx, double errYVtx)
Definition TrackExtrap.h:62
static bool extrapToVertexWithoutBranson(TrackParam &trackParam, double zVtx, double xUpstream=0., double yUpstream=0., std::optional< double > zUpstream=std::nullopt)
Definition TrackExtrap.h:74
track parameters for internal use
Definition TrackParam.h:34
Double_t getNonBendingCoor() const
return non bending coordinate (cm)
Definition TrackParam.h:51
Double_t getBendingCoor() const
return bending coordinate (cm)
Definition TrackParam.h:59
void setCellMCTruthContainer(const o2::dataformats::MCTruthContainer< o2::phos::MCLabel > *mclabels)
Setting the pointer for the MCTruthContainer for cells.
void setCellData(CellRange cells, TriggerRange triggers)
Setting the data at cell level.
void reset()
Reset containers with empty ranges.
InteractionRecord getInteractionRecordForEvent(int eventID) const
int getNumberOfEvents() const
MCTrack const * getTrack(o2::MCCompLabel const &) const
void releaseTracksForSourceAndEvent(int source, int event)
API to ask releasing tracks (freeing memory) for source + event.
std::vector< MCTrack > const & getTracks(int source, int event) const
variant returning all tracks for source and event at once
static void addInteractionBC(int bc, bool fromCollisonCotext=false)
Definition Utils.cxx:52
static constexpr ID Pion
Definition PID.h:96
Double_t getX() const
Definition TrackFwd.h:49
void set(const std::string &s="", int base=2)
Definition EnumFlags.h:401
bool match(const std::vector< std::string > &queries, const char *pattern)
Definition dcs-ccdb.cxx:229
struct _cl_event * event
Definition glcorearb.h:2982
GLfloat GLfloat GLfloat alpha
Definition glcorearb.h:279
GLint GLenum GLint x
Definition glcorearb.h:403
GLenum src
Definition glcorearb.h:1767
GLint GLsizei count
Definition glcorearb.h:399
GLuint entry
Definition glcorearb.h:5735
GLint GLint GLsizei GLuint * counters
Definition glcorearb.h:3985
GLuint GLuint end
Definition glcorearb.h:469
const GLdouble * v
Definition glcorearb.h:832
GLdouble GLdouble right
Definition glcorearb.h:4077
GLdouble f
Definition glcorearb.h:310
GLuint GLuint GLfloat weight
Definition glcorearb.h:5477
GLboolean GLboolean GLboolean b
Definition glcorearb.h:1233
GLsizei GLsizei GLchar * source
Definition glcorearb.h:798
GLsizei const GLfloat * value
Definition glcorearb.h:819
GLint GLint GLsizei GLint GLenum GLenum type
Definition glcorearb.h:275
GLenum GLsizei GLsizei GLint * values
Definition glcorearb.h:1576
GLboolean * data
Definition glcorearb.h:298
GLintptr offset
Definition glcorearb.h:660
GLuint GLsizei const GLchar * label
Definition glcorearb.h:2519
GLfloat v0
Definition glcorearb.h:811
GLbitfield flags
Definition glcorearb.h:1570
GLboolean r
Definition glcorearb.h:1233
GLuint start
Definition glcorearb.h:469
GLboolean GLboolean GLboolean GLboolean a
Definition glcorearb.h:1233
constexpr std::array< float, 2 > trackQAScaledTOF
Definition DataTypes.h:135
constexpr std::array< float, 5 > trackQAScaleContP1
Definition DataTypes.h:132
uint8_t itsSharedClusterMap uint8_t
constexpr std::array< float, 5 > trackQAScaleContP0
Definition DataTypes.h:131
constexpr std::array< float, 5 > trackQAScaleGloP0
Definition DataTypes.h:133
constexpr std::array< float, 5 > trackQAScaleGloP1
Definition DataTypes.h:134
constexpr float trackQAScaleBins
Definition DataTypes.h:129
constexpr float trackQARefRadius
Definition DataTypes.h:128
bool updateHepMCHeavyIon(const HeavyIonCursor &cursor, int collisionID, short generatorID, o2::dataformats::MCEventHeader const &header, HepMCUpdate when=HepMCUpdate::anyKey)
short updateMCCollisions(const CollisionCursor &cursor, int bcId, float time, o2::dataformats::MCEventHeader const &header, short generatorId=0, int sourceId=0, unsigned int mask=0xFFFFFFF0)
bool updateHepMCPdfInfo(const PdfInfoCursor &cursor, int collisionID, short generatorID, o2::dataformats::MCEventHeader const &header, HepMCUpdate when=HepMCUpdate::anyKey)
uint32_t updateParticles(const ParticleCursor &cursor, int collisionID, std::vector< MCTrack > const &tracks, TrackToIndex &preselect, uint32_t offset=0, bool filter=false, bool background=false, uint32_t weightMask=0xFFFFFFF0, uint32_t momentumMask=0xFFFFFFF0, uint32_t positionMask=0xFFFFFFF0)
bool updateHepMCXSection(const XSectionCursor &cursor, int collisionID, short generatorID, o2::dataformats::MCEventHeader const &header, HepMCUpdate when=HepMCUpdate::anyKey)
void keepMCParticle(std::vector< std::vector< std::unordered_map< int, int > > > &store, int source, int event, int track, int value=1)
framework::DataProcessorSpec getAODProducerWorkflowSpec(GID::mask_t src, bool enableSV, bool enableST, bool useMC, bool CTPConfigPerRun, bool enableFITextra)
create a processor spec
void dimensionMCKeepStore(std::vector< std::vector< std::unordered_map< int, int > > > &store, int Nsources, int NEvents)
void clearMCKeepStore(std::vector< std::vector< std::unordered_map< int, int > > > &store)
constexpr int LHCMaxBunches
constexpr double LHCBunchSpacingNS
constexpr float Almost1
constexpr double MassPionCharged
Defining PrimaryVertex explicitly as messageable.
Definition TFIDInfo.h:20
std::vector< ConfigParamSpec > ccdbParamSpec(std::string const &path, int runDependent, std::vector< CCDBMetadata > metadata={}, int qrate=0)
int angle2Sector(float phi)
Definition Utils.h:183
float sector2Angle(int sect)
Definition Utils.h:193
TrackParCovF TrackParCov
Definition Track.h:33
constexpr uint32_t Cal
Definition Triggers.h:32
const int TDCSignal[NTDCChannels]
Definition Constants.h:181
constexpr int NTDCChannels
Definition Constants.h:90
constexpr int NChannels
Definition Constants.h:65
std::string fullVersion()
get full version information (official O2 release and git commit)
GPUReconstruction * rec
static OutputSpec const spec()
decltype(FFL(std::declval< cursor_t >())) cursor
gsl::span< const o2::trd::TriggerRecord > getTRDTriggerRecords() const
void collectData(o2::framework::ProcessingContext &pc, const DataRequest &request)
static bool downsampleTsallisCharged(float pt, float factorPt, float sqrts, float &weight, float rnd, float mass=0.13957)
Definition Tsallis.cxx:31
std::vector< o2::ctf::BufferType > vec
LOG(info)<< "Compressed in "<< sw.CpuTime()<< " s"
Cluster clu
std::uniform_int_distribution< unsigned long long > distr
std::random_device rd
std::array< uint16_t, 5 > pattern