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