@@ -73,6 +73,8 @@ class TrackingStudySpec : public Task
7373 std::unique_ptr<o2::utils::TreeStreamRedirector> mDBGOut ;
7474 std::unique_ptr<o2::utils::TreeStreamRedirector> mDBGOutVtx ;
7575 float mITSROFrameLengthMUS = 0 .;
76+ int mMaxNeighbours = 3 ;
77+ float mMaxVTTimeDiff = 25 .; // \mus
7678 GTrackID::mask_t mTracksSrc {};
7779 o2::steer::MCKinematicsReader mcReader; // reader of MC information
7880};
@@ -82,6 +84,9 @@ void TrackingStudySpec::init(InitContext& ic)
8284 o2::base::GRPGeomHelper::instance ().setRequest (mGGCCDBRequest );
8385 mDBGOut = std::make_unique<o2::utils::TreeStreamRedirector>(" trackStudy.root" , " recreate" );
8486 mDBGOutVtx = std::make_unique<o2::utils::TreeStreamRedirector>(" trackStudyVtx.root" , " recreate" );
87+
88+ mMaxVTTimeDiff = ic.options ().get <float >(" max-vtx-timediff" );
89+ mMaxNeighbours = ic.options ().get <int >(" max-vtx-neighbours" );
8590}
8691
8792void TrackingStudySpec::run (ProcessingContext& pc)
@@ -116,6 +121,7 @@ void TrackingStudySpec::process(o2::globaltracking::RecoContainer& recoData)
116121 auto vtxRefs = recoData.getPrimaryVertexMatchedTrackRefs (); // references from vertex to these track IDs
117122 auto prop = o2::base::Propagator::Instance ();
118123 auto FITInfo = recoData.getFT0RecPoints ();
124+ static int TFCount = 0 ;
119125 int nv = vtxRefs.size () - 1 ;
120126 o2::dataformats::PrimaryVertexExt pveDummy;
121127 std::vector<o2::dataformats::PrimaryVertexExt> pveVec;
@@ -128,35 +134,40 @@ void TrackingStudySpec::process(o2::globaltracking::RecoContainer& recoData)
128134 auto & pve = pveVec.emplace_back ();
129135 static_cast <o2::dataformats::PrimaryVertex&>(pve) = pv;
130136
131- float bestTimeDiff = 1000 , bestTime = -999 , bestAmp = -1 ;
137+ float bestTimeDiff = 1000 , bestTime = -999 ;
138+ int bestFTID = -1 ;
132139 if (mTracksSrc [GTrackID::FT0]) {
133140 for (int ift0 = vtref.getFirstEntryOfSource (GTrackID::FT0); ift0 < vtref.getFirstEntryOfSource (GTrackID::FT0) + vtref.getEntriesOfSource (GTrackID::FT0); ift0++) {
134141 const auto & ft0 = FITInfo[trackIndex[ift0]];
135142 if (ft0.isValidTime (o2::ft0::RecPoints::TimeMean) && ft0.getTrigger ().getAmplA () + ft0.getTrigger ().getAmplC () > 20 ) {
136143 auto fitTime = ft0.getInteractionRecord ().differenceInBCMUS (recoData.startIR );
137144 if (std::abs (fitTime - pv.getTimeStamp ().getTimeStamp ()) < bestTimeDiff) {
138145 bestTimeDiff = fitTime - pv.getTimeStamp ().getTimeStamp ();
139- bestTime = fitTime;
140- bestAmp = ft0.getTrigger ().getAmplA () + ft0.getTrigger ().getAmplC ();
146+ bestFTID = trackIndex[ift0];
141147 }
142148 }
143149 }
144150 } else {
145151 LOGP (warn, " FT0 is not requested, cannot set complete vertex info" );
146152 }
147- pve.FT0Amp = bestAmp;
148- pve.FT0Time = bestTime;
153+ if (bestFTID >= 0 ) {
154+ pve.FT0A = FITInfo[bestFTID].getTrigger ().getAmplA ();
155+ pve.FT0Amp = pve.FT0A + FITInfo[bestFTID].getTrigger ().getAmplC ();
156+ pve.FT0Time = FITInfo[bestFTID].getInteractionRecord ().differenceInBCMUS (recoData.startIR );
157+ }
158+ pve.VtxID = iv;
149159 for (int is = 0 ; is < GTrackID::NSources; is++) {
150- if (!GTrackID::getSourceDetectorsMask (is)[GTrackID::ITS] || !mTracksSrc [is] || !recoData.isTrackSourceLoaded (is)) {
151- continue ;
152- }
160+ bool skipTracks = (!GTrackID::getSourceDetectorsMask (is)[GTrackID::ITS] || !mTracksSrc [is] || !recoData.isTrackSourceLoaded (is));
153161 int idMin = vtref.getFirstEntryOfSource (is), idMax = idMin + vtref.getEntriesOfSource (is);
154162 for (int i = idMin; i < idMax; i++) {
155163 auto vid = trackIndex[i];
156164 bool pvCont = vid.isPVContributor ();
157165 if (pvCont) {
158166 pve.nSrc [is]++;
159167 }
168+ if (skipTracks) {
169+ continue ;
170+ }
160171 bool ambig = vid.isAmbiguous ();
161172 auto trc = recoData.getTrackParam (vid);
162173 float xmin = trc.getX ();
@@ -169,17 +180,87 @@ void TrackingStudySpec::process(o2::globaltracking::RecoContainer& recoData)
169180 }
170181 }
171182 }
172-
173- for (size_t cnt = 0 ; cnt < pveVec. size () ; cnt++) {
183+ int nvtot = mMaxNeighbours < 0 ? - 1 : ( int )pveVec. size ();
184+ for (int cnt = 0 ; cnt < nvtot ; cnt++) {
174185 const auto & pve = pveVec[cnt];
175- const auto & pvePrev = cnt > 0 ? pveVec[cnt - 1 ] : pveDummy;
176- const auto & pveNext = cnt == pveVec.size () - 1 ? pveDummy : pveVec[cnt + 1 ];
186+ float tv = pve.getTimeStamp ().getTimeStamp ();
187+ std::vector<o2::dataformats::PrimaryVertexExt> pveT (mMaxNeighbours ); // neighbours in time
188+ std::vector<o2::dataformats::PrimaryVertexExt> pveZ (mMaxNeighbours ); // neighbours in Z
189+ std::vector<int > idT (mMaxNeighbours ), idZ (mMaxNeighbours );
190+ std::vector<float > dT (mMaxNeighbours ), dZ (mMaxNeighbours );
191+ for (int i = 0 ; i < mMaxNeighbours ; i++) {
192+ idT[i] = idZ[i] = -1 ;
193+ dT[i] = mMaxVTTimeDiff ;
194+ dZ[i] = 1e9 ;
195+ }
196+ int cntM = cnt - 1 , cntP = cnt + 1 ;
197+ for (; cntM >= 0 ; cntM--) { // backward
198+ const auto & vt = pveVec[cntM];
199+ auto dtime = std::abs (tv - vt.getTimeStamp ().getTimeStamp ());
200+ if (dtime > mMaxVTTimeDiff ) {
201+ continue ;
202+ }
203+ for (int i = 0 ; i < mMaxNeighbours ; i++) {
204+ if (dT[i] > dtime) {
205+ dT[i] = dtime;
206+ idT[i] = cntM;
207+ break ;
208+ }
209+ }
210+ auto dz = std::abs (pve.getZ () - vt.getZ ());
211+ for (int i = 0 ; i < mMaxNeighbours ; i++) {
212+ if (dZ[i] > dz) {
213+ dZ[i] = dz;
214+ idZ[i] = cntM;
215+ break ;
216+ }
217+ }
218+ }
219+ for (; cntP < nvtot; cntP++) { // forward
220+ const auto & vt = pveVec[cntP];
221+ auto dtime = std::abs (tv - vt.getTimeStamp ().getTimeStamp ());
222+ if (dtime > mMaxVTTimeDiff ) {
223+ continue ;
224+ }
225+ for (int i = 0 ; i < mMaxNeighbours ; i++) {
226+ if (dT[i] > dtime) {
227+ dT[i] = dtime;
228+ idT[i] = cntP;
229+ break ;
230+ }
231+ }
232+ auto dz = std::abs (pve.getZ () - vt.getZ ());
233+ for (int i = 0 ; i < mMaxNeighbours ; i++) {
234+ if (dZ[i] > dz) {
235+ dZ[i] = dz;
236+ idZ[i] = cntP;
237+ break ;
238+ }
239+ }
240+ }
241+ for (int i = 0 ; i < mMaxNeighbours ; i++) {
242+ if (idT[i] != -1 ) {
243+ pveT[i] = pveVec[idT[i]];
244+ } else {
245+ break ;
246+ }
247+ }
248+ for (int i = 0 ; i < mMaxNeighbours ; i++) {
249+ if (idZ[i] != -1 ) {
250+ pveZ[i] = pveVec[idZ[i]];
251+ } else {
252+ break ;
253+ }
254+ }
177255 (*mDBGOutVtx ) << " pvExt"
178256 << " pve=" << pve
179- << " pveP=" << pvePrev
180- << " pveN=" << pveNext
257+ << " pveT=" << pveT
258+ << " pveZ=" << pveZ
259+ << " tfID=" << TFCount
181260 << " \n " ;
182261 }
262+
263+ TFCount++;
183264}
184265
185266void TrackingStudySpec::endOfStream (EndOfStreamContext& ec)
@@ -217,7 +298,10 @@ DataProcessorSpec getTrackingStudySpec(GTrackID::mask_t srcTracks, GTrackID::mas
217298 dataRequest->inputs ,
218299 outputs,
219300 AlgorithmSpec{adaptFromTask<TrackingStudySpec>(dataRequest, ggRequest, srcTracks, useMC)},
220- Options{}};
301+ Options{
302+ {" max-vtx-neighbours" , VariantType::Int, 2 , {" Max PV neighbours fill, no PV study if < 0" }},
303+ {" max-vtx-timediff" , VariantType::Float, 25 .f , {" Max PV time difference to consider" }},
304+ }};
221305}
222306
223307} // namespace o2::trackstudy
0 commit comments