1010// or submit itself to any jurisdiction.
1111// /
1212
13+ #include < format>
14+ #include < cmath>
15+ #include < limits>
16+
1317#include " ITStracking/TrackerTraitsMC.h"
18+ #include " ITStracking/TuneExt.h"
19+ #include " CommonUtils/TreeStreamRedirector.h"
20+ #include " Steer/MCKinematicsReader.h"
1421
1522#include " Framework/Logger.h"
1623
1724namespace o2 ::its
1825{
1926
27+ static o2::utils::TreeStreamRedirector* gDBG {nullptr };
28+ static o2::steer::MCKinematicsReader* gMCReader {nullptr };
29+
2030template <int NLayers>
2131void TrackerTraitsMC<NLayers>::computeLayerTracklets(const int iteration, int iVertex)
2232{
33+ if (!gDBG ) {
34+ gDBG = new o2::utils::TreeStreamRedirector (" its_tune.root" );
35+ }
36+ if (!gMCReader ) {
37+ gMCReader = new o2::steer::MCKinematicsReader (" collisioncontext.root" );
38+ }
39+
40+ // Create all tracklets we find in this iteration and dump them.
41+ const std::string treeName = std::format (" trklt_{}" , iteration);
2342 TrackerTraits<NLayers>::computeLayerTracklets (iteration, iVertex);
43+ this ->createTrackletMC ();
44+ const auto topology = this ->mTimeFrame ->getTrackingTopologyView ();
45+ for (int transitionId{0 }; transitionId < topology.nTransitions ; ++transitionId) {
46+ const auto & transition = topology.getTransition (transitionId);
47+ for (int iTrklt{0 }; iTrklt < this ->mTimeFrame ->getTracklets ()[transitionId].size (); ++iTrklt) {
48+ const auto & lbl = this ->mTimeFrame ->getTrackletsLabel (transitionId)[iTrklt];
49+ const auto trklt = this ->mTimeFrame ->getTracklets ()[transitionId][iTrklt];
50+ TrackletMC trkltMC{
51+ .tgl = trklt.tanLambda ,
52+ .phi = trklt.phi ,
53+ .ok = lbl.isValid (),
54+ };
55+ float dcaXY = std::numeric_limits<float >::max (), dcaZ = std::numeric_limits<float >::max ();
56+ if (lbl.isValid ()) {
57+ const auto & eve = gMCReader ->getMCEventHeader (lbl.getSourceID (), lbl.getEventID ());
58+ const auto & firstCluster = this ->mTimeFrame ->getClusters ()[transition.fromLayer ][trklt.firstClusterIndex ];
59+ const auto & secondCluster = this ->mTimeFrame ->getClusters ()[transition.toLayer ][trklt.secondClusterIndex ];
60+ const float dx = secondCluster.xCoordinate - firstCluster.xCoordinate ;
61+ const float dy = secondCluster.yCoordinate - firstCluster.yCoordinate ;
62+ const float dz = secondCluster.zCoordinate - firstCluster.zCoordinate ;
63+ const float dxy2 = math_utils::hypot (dx, dy);
64+ if (dxy2 > constants::Tolerance) {
65+ const float t = ((eve.GetX () - firstCluster.xCoordinate ) * dx + (eve.GetY () - firstCluster.yCoordinate ) * dy) / dxy2;
66+ const float xAtDCA = firstCluster.xCoordinate + t * dx;
67+ const float yAtDCA = firstCluster.yCoordinate + t * dy;
68+ const float zAtDCA = firstCluster.zCoordinate + t * dz;
69+ const float curDCAx = xAtDCA - eve.GetX ();
70+ const float curDCAy = yAtDCA - eve.GetY ();
71+ dcaXY = math_utils::hypot (curDCAx, curDCAy);
72+ dcaZ = zAtDCA - eve.GetZ ();
73+ trkltMC.dXY = dcaXY;
74+ trkltMC.dZ = dcaZ;
75+ }
76+ const auto * mcTrk = gMCReader ->getTrack (lbl);
77+ if (mcTrk) {
78+ trkltMC.prim = mcTrk->isPrimary ();
79+ }
80+ }
81+ (*gDBG ) << treeName.c_str ()
82+ << " from=" << transition.fromLayer
83+ << " to=" << transition.toLayer
84+ << " trklt=" << trkltMC
85+ << " \n " ;
86+ }
87+ }
2488}
2589
2690template <int NLayers>
@@ -39,6 +103,10 @@ template <int NLayers>
39103void TrackerTraitsMC<NLayers>::findRoads(const int iteration)
40104{
41105 TrackerTraits<NLayers>::findRoads (iteration);
106+
107+ if (this ->mTrkParams .size () - 1 == iteration) {
108+ gDBG ->Close ();
109+ }
42110}
43111
44112template class TrackerTraitsMC <7 >;
0 commit comments