4343#include " DataFormatsTPC/PIDResponse.h"
4444#include " DataFormatsITS/TrackITS.h"
4545#include " TROOT.h"
46+ #include " ReconstructionDataFormats/MatchInfoTOF.h"
47+ #include " DataFormatsTOF/Cluster.h"
4648
4749using namespace o2 ::globaltracking;
4850using GTrackID = o2::dataformats::GlobalTrackID;
@@ -170,7 +172,10 @@ class TPCTimeSeries : public Task
170172 auto primMatchedTracks = mTPCOnly ? gsl::span<o2::dataformats::VtxTrackIndex>() : recoData.getPrimaryVertexMatchedTracks (); // Global ID's for associated tracks
171173 auto primMatchedTracksRef = mTPCOnly ? gsl::span<o2::dataformats::VtxTrackRef>() : recoData.getPrimaryVertexMatchedTrackRefs (); // references from vertex to these track IDs
172174
173- LOGP (info, " Processing {} vertices, {} primary matched vertices, {} TPC tracks, {} ITS tracks, {} ITS-TPC tracks" , vertices.size (), primMatchedTracks.size (), tracksTPC.size (), tracksITS.size (), tracksITSTPC.size ());
175+ // TOF clusters
176+ const auto & tofClusters = mTPCOnly ? gsl::span<o2::tof::Cluster>() : recoData.getTOFClusters ();
177+
178+ LOGP (info, " Processing {} vertices, {} primary matched vertices, {} TPC tracks, {} ITS tracks, {} ITS-TPC tracks, {} TOF clusters" , vertices.size (), primMatchedTracks.size (), tracksTPC.size (), tracksITS.size (), tracksITSTPC.size (), tofClusters.size ());
174179
175180 // calculate mean vertex, RMS and count vertices
176181 auto indicesITSTPC_vtx = processVertices (vertices, primMatchedTracks, primMatchedTracksRef, recoData);
@@ -186,6 +191,24 @@ class TPCTimeSeries : public Task
186191 indicesITSTPC[tracksITSTPC[i].getRefTPC ().getIndex ()] = {i, idxVtx};
187192 }
188193
194+ std::vector<std::tuple<int , float , float >> idxTPCTrackToTOFCluster; // store for each tpc track index the index to the TOF cluster
195+
196+ // get matches to TOF in case skimmed data is produced
197+ if (mUnbinnedWriter ) {
198+ idxTPCTrackToTOFCluster = std::vector<std::tuple<int , float , float >>(tracksTPC.size (), {-1 , -999 , -999 });
199+ const std::vector<gsl::span<const o2::dataformats::MatchInfoTOF>> tofMatches{recoData.getTPCTOFMatches (), recoData.getTPCTRDTOFMatches (), recoData.getITSTPCTOFMatches (), recoData.getITSTPCTRDTOFMatches ()};
200+
201+ // loop over ITS-TPC-TRD-TOF and ITS-TPC-TOF tracks an store for each ITS-TPC track the TOF track index
202+ for (const auto & tofMatch : tofMatches) {
203+ for (const auto & tpctofmatch : tofMatch) {
204+ auto refTPC = recoData.getTPCContributorGID (tpctofmatch.getTrackRef ());
205+ if (refTPC.isIndexSet ()) {
206+ idxTPCTrackToTOFCluster[refTPC] = {tpctofmatch.getIdxTOFCl (), tpctofmatch.getDXatTOF (), tpctofmatch.getDZatTOF ()};
207+ }
208+ }
209+ }
210+ }
211+
189212 // find nearest vertex of tracks which have no vertex assigned
190213 findNearesVertex (tracksTPC, vertices);
191214
@@ -358,7 +381,7 @@ class TPCTimeSeries : public Task
358381 auto myThread = [&](int iThread) {
359382 for (size_t i = iThread; i < loopEnd; i += mNThreads ) {
360383 if (acceptTrack (tracksTPC[i])) {
361- fillDCA (tracksTPC, tracksITSTPC, vertices, i, iThread, indicesITSTPC, tracksITS);
384+ fillDCA (tracksTPC, tracksITSTPC, vertices, i, iThread, indicesITSTPC, tracksITS, idxTPCTrackToTOFCluster, tofClusters );
362385 }
363386 }
364387 };
@@ -375,7 +398,7 @@ class TPCTimeSeries : public Task
375398 auto myThread = [&](int iThread) {
376399 for (size_t i = iThread; i < loopEnd; i += mNThreads ) {
377400 if (acceptTrack (tracksTPC[i])) {
378- fillDCA (tracksTPC, tracksITSTPC, vertices, i, iThread, indicesITSTPC, tracksITS);
401+ fillDCA (tracksTPC, tracksITSTPC, vertices, i, iThread, indicesITSTPC, tracksITS, idxTPCTrackToTOFCluster, tofClusters );
379402 }
380403 }
381404 };
@@ -1001,9 +1024,10 @@ class TPCTimeSeries : public Task
10011024 return true ;
10021025 }
10031026
1004- void fillDCA (const gsl::span<const TrackTPC> tracksTPC, const gsl::span<const o2::dataformats::TrackTPCITS> tracksITSTPC, const gsl::span<const o2::dataformats::PrimaryVertex> vertices, const int iTrk, const int iThread, const std::unordered_map<unsigned int , std::array<int , 2 >>& indicesITSTPC, const gsl::span<const o2::its::TrackITS> tracksITS)
1027+ void fillDCA (const gsl::span<const TrackTPC> tracksTPC, const gsl::span<const o2::dataformats::TrackTPCITS> tracksITSTPC, const gsl::span<const o2::dataformats::PrimaryVertex> vertices, const int iTrk, const int iThread, const std::unordered_map<unsigned int , std::array<int , 2 >>& indicesITSTPC, const gsl::span<const o2::its::TrackITS> tracksITS, const std::vector<std::tuple< int , float , float >>& idxTPCTrackToTOFCluster, const gsl::span< const o2::tof::Cluster> tofClusters )
10051028 {
1006- TrackTPC track = tracksTPC[iTrk];
1029+ o2::track::TrackParCov track = tracksTPC[iTrk];
1030+ const auto & trackFull = tracksTPC[iTrk];
10071031
10081032 // propagate track to the DCA and fill in slice
10091033 auto propagator = o2::base::Propagator::Instance ();
@@ -1055,11 +1079,11 @@ class TPCTimeSeries : public Task
10551079 const auto vertex = (idxITSTPC.back () != -1 ) ? vertices[idxITSTPC.back ()] : ((mNearestVtxTPC [iTrk] != -1 ) ? vertices[mNearestVtxTPC [iTrk]] : o2::dataformats::PrimaryVertex{});
10561080
10571081 // calculate DCAz: (time TPC track - time vertex) * vDrift + sign_side * vertexZ
1058- const float signSide = track .hasCSideClustersOnly () ? -1 : 1 ; // invert sign for C-side
1059- const float dcaZFromDeltaTime = (vertex.getTimeStamp ().getTimeStamp () == 0 ) ? 0 : (o2::tpc::ParameterElectronics::Instance ().ZbinWidth * track .getTime0 () - vertex.getTimeStamp ().getTimeStamp ()) * mVDrift + signSide * vertex.getZ ();
1082+ const float signSide = trackFull .hasCSideClustersOnly () ? -1 : 1 ; // invert sign for C-side
1083+ const float dcaZFromDeltaTime = (vertex.getTimeStamp ().getTimeStamp () == 0 ) ? 0 : (o2::tpc::ParameterElectronics::Instance ().ZbinWidth * trackFull .getTime0 () - vertex.getTimeStamp ().getTimeStamp ()) * mVDrift + signSide * vertex.getZ ();
10601084
10611085 // for weight of DCA
1062- const float resCl = std::min (track .getNClusters (), static_cast <int >(Mapper::PADROWS)) / static_cast <float >(Mapper::PADROWS);
1086+ const float resCl = std::min (trackFull .getNClusters (), static_cast <int >(Mapper::PADROWS)) / static_cast <float >(Mapper::PADROWS);
10631087
10641088 const float div = (resCl * track.getPt ());
10651089 if (div == 0 ) {
@@ -1087,15 +1111,15 @@ class TPCTimeSeries : public Task
10871111 }
10881112
10891113 const float chi2Match = (chi2 > 0 ) ? std::sqrt (chi2) : -1 ;
1090- const float sqrtChi2TPC = (track .getChi2 () > 0 ) ? std::sqrt (track .getChi2 ()) : 0 ;
1091- const float nClTPC = track .getNClusters ();
1114+ const float sqrtChi2TPC = (trackFull .getChi2 () > 0 ) ? std::sqrt (trackFull .getChi2 ()) : 0 ;
1115+ const float nClTPC = trackFull .getNClusters ();
10921116
10931117 // const float dedx = mUseQMax ? track.getdEdx().dEdxMaxTPC : track.getdEdx().dEdxTotTPC;
1094- const float dedxRatioqTot = (track .getdEdx ().dEdxTotTPC > 0 ) ? (mMIPdEdx / track .getdEdx ().dEdxTotTPC ) : -1 ;
1095- const float dedxRatioqMax = (track .getdEdx ().dEdxMaxTPC > 0 ) ? (mMIPdEdx / track .getdEdx ().dEdxMaxTPC ) : -1 ;
1118+ const float dedxRatioqTot = (trackFull .getdEdx ().dEdxTotTPC > 0 ) ? (mMIPdEdx / trackFull .getdEdx ().dEdxTotTPC ) : -1 ;
1119+ const float dedxRatioqMax = (trackFull .getdEdx ().dEdxMaxTPC > 0 ) ? (mMIPdEdx / trackFull .getdEdx ().dEdxMaxTPC ) : -1 ;
10961120
1097- const auto dedxQTotVars = getdEdxVars (0 , track );
1098- const auto dedxQMaxVars = getdEdxVars (1 , track );
1121+ const auto dedxQTotVars = getdEdxVars (0 , trackFull );
1122+ const auto dedxQMaxVars = getdEdxVars (1 , trackFull );
10991123
11001124 // make check to avoid crash in case no or less ITS tracks have been found!
11011125 const int idxITSTrack = (hasITSTPC && (gID == o2::dataformats::GlobalTrackID::Source::ITS)) ? tracksITSTPC[idxITSTPC.front ()].getRefITS ().getIndex () : -1 ;
@@ -1108,9 +1132,9 @@ class TPCTimeSeries : public Task
11081132 }
11091133 sigmaY2 = track.getSigmaY2 ();
11101134 sigmaZ2 = track.getSigmaZ2 ();
1111- if (track .hasCSideClustersOnly ()) {
1135+ if (trackFull .hasCSideClustersOnly ()) {
11121136 mBufferVals [iThread].front ().emplace_back (Side::C, tglBin, phiBin, qPtBin, multBin, dca[0 ], dcaZFromDeltaTime, dcarW, dedxRatioqTot, dedxRatioqMax, sqrtChi2TPC, nClTPC, gID , chi2Match, hasITSTPC, nClITS, chi2ITS, dedxQTotVars, dedxQMaxVars, sigmaY2, sigmaZ2);
1113- } else if (track .hasASideClustersOnly ()) {
1137+ } else if (trackFull .hasASideClustersOnly ()) {
11141138 mBufferVals [iThread].front ().emplace_back (Side::A, tglBin, phiBin, qPtBin, multBin, dca[0 ], dcaZFromDeltaTime, dcarW, dedxRatioqTot, dedxRatioqMax, sqrtChi2TPC, nClTPC, gID , chi2Match, hasITSTPC, nClITS, chi2ITS, dedxQTotVars, dedxQMaxVars, sigmaY2, sigmaZ2);
11151139 }
11161140
@@ -1150,9 +1174,9 @@ class TPCTimeSeries : public Task
11501174 dcaITSTPCTmp[1 ] = -1 ;
11511175 }
11521176
1153- if (track .hasCSideClustersOnly ()) {
1177+ if (trackFull .hasCSideClustersOnly ()) {
11541178 mBufferVals [iThread].back ().emplace_back_ITSTPC (Side::C, tglBin, phiBin, qPtBin, multBin, dca[0 ], dcaZFromDeltaTime, dcarW, dedxRatioqTot, dedxRatioqMax, sqrtChi2TPC, nClTPC, dcaITSTPCTmp[0 ], dcaITSTPCTmp[1 ]);
1155- } else if (track .hasASideClustersOnly ()) {
1179+ } else if (trackFull .hasASideClustersOnly ()) {
11561180 mBufferVals [iThread].back ().emplace_back_ITSTPC (Side::A, tglBin, phiBin, qPtBin, multBin, dca[0 ], dcaZFromDeltaTime, dcarW, dedxRatioqTot, dedxRatioqMax, sqrtChi2TPC, nClTPC, dcaITSTPCTmp[0 ], dcaITSTPCTmp[1 ]);
11571181 }
11581182 }
@@ -1168,7 +1192,7 @@ class TPCTimeSeries : public Task
11681192 writeData = o2::math_utils::Tsallis::downsampleTsallisCharged (tracksTPC[iTrk].getPt (), factorPt, mSqrt , weight, distr (mGenerator ));
11691193 }
11701194 if (writeData) {
1171- auto clusterMask = makeClusterBitMask (track );
1195+ auto clusterMask = makeClusterBitMask (trackFull );
11721196 const auto & trkOrig = tracksTPC[iTrk];
11731197 const bool isNearestVtx = (idxITSTPC.back () == -1 ); // is nearest vertex in case no vertex was found
11741198 const float mx_ITS = hasITSTPC ? tracksITSTPC[idxITSTPC.front ()].getX () : -1 ;
@@ -1177,12 +1201,30 @@ class TPCTimeSeries : public Task
11771201 const int nClITS = idxITSCheck ? tracksITS[idxITSTrack].getNClusters () : -1 ;
11781202 const int chi2ITS = idxITSCheck ? tracksITS[idxITSTrack].getChi2 () : -1 ;
11791203 int typeSide = 2 ; // A- and C-Side cluster
1180- if (track .hasASideClustersOnly ()) {
1204+ if (trackFull .hasASideClustersOnly ()) {
11811205 typeSide = 0 ;
1182- } else if (track .hasCSideClustersOnly ()) {
1206+ } else if (trackFull .hasCSideClustersOnly ()) {
11831207 typeSide = 1 ;
11841208 }
11851209
1210+ // check for TOF and propagate TPC track to TOF cluster
1211+ bool hasTOFCluster = (std::get<0 >(idxTPCTrackToTOFCluster[iTrk]) != -1 );
1212+ auto tofCl = hasTOFCluster ? tofClusters[std::get<0 >(idxTPCTrackToTOFCluster[iTrk])] : o2::tof::Cluster ();
1213+
1214+ float tpcYDeltaAtTOF = -999 ;
1215+ float tpcZDeltaAtTOF = -999 ;
1216+ if (hasTOFCluster) {
1217+ o2::track::TrackPar trackTmpOut (tracksTPC[iTrk].getParamOut ());
1218+ if (trackTmpOut.rotate (o2::math_utils::sector2Angle (tofCl.getSector ())) && propagator->propagateTo (trackTmpOut, tofCl.getX (), false , mMaxSnp , mFineStep , mMatType )) {
1219+ tpcYDeltaAtTOF = trackTmpOut.getY () - tofCl.getY ();
1220+ tpcZDeltaAtTOF = signSide * (o2::tpc::ParameterElectronics::Instance ().ZbinWidth * trackFull.getTime0 () - vertex.getTimeStamp ().getTimeStamp ()) * mVDrift - trackTmpOut.getZ () + tofCl.getZ ();
1221+ }
1222+ }
1223+
1224+ // get delta parameter between inner and outer
1225+ float deltaTPCParamInOutTgl = trackFull.getTgl () - trackFull.getParamOut ().getTgl ();
1226+ float deltaTPCParamInOutQPt = trackFull.getQ2Pt () - trackFull.getParamOut ().getQ2Pt ();
1227+
11861228 *mStreamer [iThread] << " treeTimeSeries"
11871229 // DCAs
11881230 << " factorPt=" << factorPt
@@ -1231,6 +1273,14 @@ class TPCTimeSeries : public Task
12311273 << " mVDrift=" << mVDrift
12321274 << " its_flag=" << int (gID )
12331275 << " sqrtChi2Match=" << chi2Match
1276+ // TOF cluster
1277+ << " tpcYDeltaAtTOF=" << tpcYDeltaAtTOF
1278+ << " tpcZDeltaAtTOF=" << tpcZDeltaAtTOF
1279+ << " mDXatTOF=" << std::get<1 >(idxTPCTrackToTOFCluster[iTrk])
1280+ << " mDZatTOF=" << std::get<2 >(idxTPCTrackToTOFCluster[iTrk])
1281+ // TPC delta param
1282+ << " deltaTPCParamInOutTgl=" << deltaTPCParamInOutTgl
1283+ << " deltaTPCParamInOutQPt=" << deltaTPCParamInOutQPt
12341284 << " \n " ;
12351285 }
12361286 }
@@ -1591,7 +1641,7 @@ o2::framework::DataProcessorSpec getTPCTimeSeriesSpec(const bool disableWriter,
15911641 {" max-dedx-ratio" , VariantType::Float, 0 .3f , {" Maximum absolute log(dedx(pion)/dedx) ratio" }},
15921642 {" max-dedx-region-ratio" , VariantType::Float, 0 .5f , {" Maximum absolute log(dedx(region)/dedx) ratio" }},
15931643 {" sample-unbinned-tsallis" , VariantType::Bool, false , {" Perform sampling of unbinned data based on Tsallis function" }},
1594- {" sampling-factor" , VariantType::Float, 0 .1f , {" Sampling factor in case sample-unbinned-tsallis is used" }},
1644+ {" sampling-factor" , VariantType::Float, 0 .001f , {" Sampling factor in case sample-unbinned-tsallis is used" }},
15951645 {" out-file-unbinned" , VariantType::String, " time_series_tracks.root" , {" name of the output file for the unbinned data" }}}};
15961646}
15971647
0 commit comments