From 6d1c2c0e316294e0f9372ecb4a748fd5da3f0391 Mon Sep 17 00:00:00 2001 From: Sean Murray Date: Fri, 7 Jun 2024 17:55:41 +0200 Subject: [PATCH] key mistmatches --- .../qc/include/TRDQC/CoordinateTransformer.h | 25 +- .../TRD/qc/include/TRDQC/RawDataManager.h | 63 ++- Detectors/TRD/qc/include/TRDQC/RawDisplay.h | 1 + Detectors/TRD/qc/macros/DrawMCMs.C | 1 - Detectors/TRD/qc/macros/PlotTrackRef.C | 4 - .../TRD/qc/src/CoordinateTransformer.cxx | 15 +- Detectors/TRD/qc/src/RawDataManager.cxx | 512 +++++++++++++++++- Detectors/TRD/qc/src/RawDisplay.cxx | 44 +- 8 files changed, 615 insertions(+), 50 deletions(-) diff --git a/Detectors/TRD/qc/include/TRDQC/CoordinateTransformer.h b/Detectors/TRD/qc/include/TRDQC/CoordinateTransformer.h index 14cb9d6b9bfc4..a76d85ac9378c 100644 --- a/Detectors/TRD/qc/include/TRDQC/CoordinateTransformer.h +++ b/Detectors/TRD/qc/include/TRDQC/CoordinateTransformer.h @@ -18,6 +18,8 @@ /// #include "DataFormatsTRD/Hit.h" +#include "DataFormatsTRD/HelperMethods.h" +#include "SimulationDataFormat/TrackReference.h" #include @@ -36,7 +38,7 @@ class ChamberSpacePoint ChamberSpacePoint(int det = -999) : mDetector(det){}; ChamberSpacePoint(int id, int detector, float x, float y, float z, std::array rct, bool inDrift) : mID(id), mDetector(detector), mX(x), mY(y), mZ(z), mPadrow(rct[0]), mPadcol(rct[1]), mTimebin(rct[2]), mInDrift(inDrift){}; - + //ChamberSpacePoint(o2::TrackReference &ref): mID(ref->getTrackID()), mDetector(ref->getDetectorId()), mX(ref->X()),mY(ref->Y()), mZ(ref->Z()), /// check if the space point has been initialized bool isValid() const { return mDetector >= 0; } @@ -74,10 +76,10 @@ class ChamberSpacePoint bool isInMCM(int detector, int padrow, int mcmcol) const; /// calculate MCM corresponding to pad row/column - // int getMCM() const { return o2::trd::HelperMethods::getMCMfromPad(mPadrow, mPadcol); } + int getMCM() const { return o2::trd::HelperMethods::getMCMfromPad(mPadrow, mPadcol); } /// calculate readout board corresponding to pad row/column - // int getROB() const { return o2::trd::HelperMethods::getROBfromPad(mPadrow, mPadcol); } + int getROB() const { return o2::trd::HelperMethods::getROBfromPad(mPadrow, mPadcol); } protected: float mX, mY, mZ; @@ -98,17 +100,23 @@ std::ostream& operator<<(std::ostream& os, const ChamberSpacePoint& p); class HitPoint : public ChamberSpacePoint { public: - HitPoint(ChamberSpacePoint position, float charge) - : ChamberSpacePoint(position), mCharge(charge) + HitPoint(ChamberSpacePoint position, float charge=0., int trackid=0) + : ChamberSpacePoint(position), mCharge(charge), mTrackID(trackid) { } HitPoint(){}; float getCharge() { return mCharge; } + const uint32_t getUserId()const { return mUserId; } + void setUserId(uint32_t userid){ mUserId=userid; } + const uint32_t getTrackID()const { return mTrackID; } + void setTrackID(uint32_t trackid){ mTrackID=trackid; } private: float mCharge{0.0}; + uint32_t mUserId{0}; // pulled in from a trackref. + uint32_t mTrackID{0}; }; /// A track segment: a straight line connecting two points. The points are generally given in spatial coordinates, @@ -124,6 +132,11 @@ class TrackSegment { assert(start.getDetector() == end.getDetector()); } + /* TrackSegment(o2::TrackReference start, o2::TrackReference end, int id, int det) + : mStartPoint(start), mEndPoint(end), mTrackID(id) + { + assert(start.getDetector() == end.getDetector()); + }*/ /// check if the space point has been initialized bool isValid() const { return mStartPoint.isValid(); } @@ -187,6 +200,8 @@ class CoordinateTransformer o2::trd::ChamberSpacePoint MakeSpacePoint(o2::trd::Hit& hit); + o2::trd::ChamberSpacePoint MakeSpacePoint(o2::TrackReference& ref); + /// Legacy, less accurate method to convert local spatial to row/column/time coordinate. /// This method is only included for comparision, and should be removed in the future. std::array OrigLocal2RCT(int det, float x, float y, float z); diff --git a/Detectors/TRD/qc/include/TRDQC/RawDataManager.h b/Detectors/TRD/qc/include/TRDQC/RawDataManager.h index ee8ce96128f95..118c53ac8aad1 100644 --- a/Detectors/TRD/qc/include/TRDQC/RawDataManager.h +++ b/Detectors/TRD/qc/include/TRDQC/RawDataManager.h @@ -21,6 +21,7 @@ #include "DataFormatsTRD/Tracklet64.h" #include "DataFormatsTRD/TriggerRecord.h" #include "DataFormatsTRD/Hit.h" +#include "TRDBase/Geometry.h" #include "TRDQC/CoordinateTransformer.h" @@ -44,6 +45,40 @@ class TTree; namespace o2::trd { +class MCTrackletSegmentInfo { + // simply a place to store the geant entry and exit points and the associated trd tracklet +public: + MCTrackletSegmentInfo(){}; + MCTrackletSegmentInfo(o2::TrackReference &in, o2::TrackReference &out, int trackid, int det) + : mEnter(in), mExit(out), mTrackId(trackid), mDet(det) {}; + void setEntry(o2::TrackReference& entry){mEnter=entry;mHasEnter=true;} + void setExit(o2::TrackReference& exit){mExit=exit;mHasExit=true;} + void setMidPoint(){/*std::cout << "setting midpoint " << mEnter << " :: " << mExit << "\n";*/ mMidPoint[0]=(mEnter.X()+mExit.X())/2.; mMidPoint[1]=(mEnter.Y()+mExit.Y())/2.; mMidPoint[2]=(mEnter.Z()+mExit.Z())/2.;/* std::cout << "set midpoint\n";*/} // contract the enter and exit points to guarantee both are with in UJ UK; + o2::TrackReference mEnter; + ChamberSpacePoint mEnterSpace; + o2::TrackReference mExit; + ChamberSpacePoint mExitSpace; + std::array mMidPoint; // store the mid point of the segment for purposes of figuring out the position in the geometry + o2::trd::Tracklet64 mTracklet; // enables us to match against a tracklet as last resort. + int mTrackId; + int mDet; + int mcmid; + bool mHasEnter{false}; + bool mHasExit{false}; + //methods for satisfying templates using the exit point (pad side) + int getDetector()const{return mDet;} + int getPadRow()const;//{return 1;} + int getPadCol()const;//{return 1;} + float getROBf()const;//{return 1;} + float getMCMf()const;//{return 1;} + int getROB()const;//{return 1;} + int getMCM()const;//{return 1;} + bool isGood(){ if (mEnter.getLength() >0.1 && mExit.getLength()>0.1) return true; else return false;} +}; + + +std::ostream& operator<<(std::ostream& os, const MCTrackletSegmentInfo& p); + /// RawDataSpan holds ranges pointing to parts of other containers /// This class helps with restricting the view of digits/trackets/... to part /// of a timeframe, e.g. data from one event, detector, padrow or MCM. @@ -54,7 +89,8 @@ struct RawDataSpan { boost::iterator_range::iterator> digits; boost::iterator_range::iterator> tracklets; boost::iterator_range::iterator> hits; - boost::iterator_range::iterator> trackrefs; +// boost::iterator_range::iterator> trackrefs; + boost::iterator_range::iterator> trackrefsegments; /// Sort digits, tracklets and space points by detector, pad row, column /// The digits, tracklets, hits and other future data members must be sorted @@ -77,9 +113,19 @@ struct RawDataSpan { std::vector iterateByPadRow(); // std::vector iterateDetector(); - std::vector makeMCTrackSegments(); + std::vector makeMCTrackSegments(); // this is based on hits +// std::vector makeMCTrackSegmentsEntryExit(); + //this one is based on TrackReferences. + std::vector makeMCTrackSegmentsEntryExit(); + //std::vector makeMCTrackSegmentsEntryExit(); std::vector makeMCTrackReferences(); + // std::vector makeMCTrackSegmentsGeant(); + std::vector makeMCTrackSegmentsGeant(); + + //convert a track ref to a hit so that the global and local coordinates are kept togther. + Hit convertTrackReferenceToHit(o2::TrackReference &ref, int trackid, int detector); + // pair getMaxADCsumAndChannel(); // int getMaxADCsum(){ return getMaxADCsumAndChannel().first; } @@ -162,16 +208,25 @@ class RawDataManager TFile* mMCFile{0}; TTree* mMCTree{0}; // TTreeReader* mMCReader{0}; - std::vector* mMCEventHeader{0}; + //std::vector* mMCEventHeader{0}; + o2::dataformats::MCEventHeader* mMCEventHeader{0}; std::vector>* mMCTracks{0}; std::vector* mMCTrackReferences{0}; + std::vector mFilteredMCTrackReferences{0}; + std::vector mMCTrackletSegmentInfo{0}; std::vector* mHits{0}; // MC hits, converted to chamber coordinates std::vector mHitPoints; + // process the MC Trackreferences into MCTrackletSegmentInfo + void processTrackReferences(); + // MC track references, converted to chamber coordinates - std::vector mTrackReferences; + std::vector mTrackReferences; + + // MC track references, converted to chamber coordinates, for those matching entry and exit of a chamber. + std::vector mMatchedTrackReferences; // time frame information (for data only) std::vector* mTFIDs{0}; diff --git a/Detectors/TRD/qc/include/TRDQC/RawDisplay.h b/Detectors/TRD/qc/include/TRDQC/RawDisplay.h index fb7f7014d0342..9ae0967943739 100644 --- a/Detectors/TRD/qc/include/TRDQC/RawDisplay.h +++ b/Detectors/TRD/qc/include/TRDQC/RawDisplay.h @@ -39,6 +39,7 @@ class RawDisplay void drawMCTrackSegments(); void drawMCTrackReferences(); void drawTrackRef(); + void drawMCTrackRefEntryExit(); void draw() { diff --git a/Detectors/TRD/qc/macros/DrawMCMs.C b/Detectors/TRD/qc/macros/DrawMCMs.C index 9d7bb2b948907..5ce317af45288 100644 --- a/Detectors/TRD/qc/macros/DrawMCMs.C +++ b/Detectors/TRD/qc/macros/DrawMCMs.C @@ -86,7 +86,6 @@ void DrawMCMs(std::string dirname = ".") disp.drawDigits("text,same"); disp.drawClusters(); disp.drawTracklets(); - disp.drawHits(); disp.drawMCTrackSegments(); } diff --git a/Detectors/TRD/qc/macros/PlotTrackRef.C b/Detectors/TRD/qc/macros/PlotTrackRef.C index a141c511d0d17..3a3dd3fe07443 100644 --- a/Detectors/TRD/qc/macros/PlotTrackRef.C +++ b/Detectors/TRD/qc/macros/PlotTrackRef.C @@ -57,10 +57,6 @@ using namespace o2::detectors; -/*std::vector transformtracklet(o2::trd::Tracklet64 tracklet) { - // this is ported kind of from tracklettransformer but that needs ccdb - // connections and other stuff. -}*/ class MCTracklet { public: MCTracklet(o2::TrackReference &in, o2::TrackReference &out, int trackid, diff --git a/Detectors/TRD/qc/src/CoordinateTransformer.cxx b/Detectors/TRD/qc/src/CoordinateTransformer.cxx index 02206db10ca30..f7217d0e8992c 100644 --- a/Detectors/TRD/qc/src/CoordinateTransformer.cxx +++ b/Detectors/TRD/qc/src/CoordinateTransformer.cxx @@ -96,7 +96,7 @@ std::array CoordinateTransformer::OrigLocal2RCT(int det, float x, floa double lengthCorr = padPlane->getLengthOPad() / padPlane->getLengthIPad(); // calculate position based on inner pad length - rct[0] = -z / padPlane->getLengthIPad() + padPlane->getNrows() / 2; + rct[0] = -z / padPlane->getLengthIPad() + padPlane->getNrows() / 2.0; // correct row for outer pad rows if (rct[0] <= 1.0) { @@ -125,7 +125,7 @@ std::array CoordinateTransformer::OrigLocal2RCT(int det, float x, floa // anode region: very rough guess rct[2] = mT0 - 1.0 + fabs(x); } - +LOGP(info," rct after OrigLocal2RCT : x:y:z {}:{}:{} => {}:{}:{}",x,y,z,rct[0], rct[1], rct[2]); return rct; } @@ -138,6 +138,15 @@ o2::trd::ChamberSpacePoint CoordinateTransformer::MakeSpacePoint(o2::trd::Hit& h return o2::trd::ChamberSpacePoint(hit.GetTrackID(), hit.GetDetectorID(), x, y, y, rct, hit.isFromDriftRegion()); } +o2::trd::ChamberSpacePoint CoordinateTransformer::MakeSpacePoint(o2::TrackReference& ref) +{ + float x = ref.X(); + float y = ref.Y(); + float z = ref.Z(); + auto rct = Local2RCT(ref.getDetectorId(), x, y, z); + return o2::trd::ChamberSpacePoint(ref.getTrackID(), ref.getDetectorId(), x, y, z, rct, true); +} + namespace o2::trd { std::ostream& operator<<(std::ostream& os, const ChamberSpacePoint& p) @@ -154,4 +163,4 @@ std::ostream& operator<<(std::ostream& os, const ChamberSpacePoint& p) << " pad " << std::setprecision(5) << p.getPadCol() << "]"; return os; } -}; // namespace o2::trd \ No newline at end of file +}; // namespace o2::trd diff --git a/Detectors/TRD/qc/src/RawDataManager.cxx b/Detectors/TRD/qc/src/RawDataManager.cxx index 462a44d5995eb..187fcdb726046 100644 --- a/Detectors/TRD/qc/src/RawDataManager.cxx +++ b/Detectors/TRD/qc/src/RawDataManager.cxx @@ -10,6 +10,12 @@ // or submit itself to any jurisdiction. #include "TRDQC/RawDataManager.h" +#include "DataFormatsTRD/HelperMethods.h" +#include "TRDQC/CoordinateTransformer.h" + +#include "Framework/Logger.h" +#include +#include #include #include @@ -18,9 +24,6 @@ #include #include #include -#include "TRDQC/CoordinateTransformer.h" -#include "Framework/Logger.h" - #include #include @@ -92,11 +95,60 @@ bool comp_spacepoint(const ChamberSpacePoint& a, const ChamberSpacePoint& b) return true; } +bool comp_trackrefs(const o2::TrackReference& a, const o2::TrackReference& b) +{ + if (a.getTrackID() != b.getTrackID()) { + return a.getTrackID() < b.getTrackID(); + } + + if (a.getUserId()>>2 != b.getUserId()>>2) { + return a.getUserId()>>2 < b.getUserId()>>2; // this gets the detector/stack/layer + } + + if (a.getTime() != b.getTime()) { + return a.getTime() < b.getTime(); // time of flight, so order from interaction region + } + + return true; +} + +bool comp_trackrefshits(const MCTrackletSegmentInfo& a, const MCTrackletSegmentInfo& b) +{ + if (a.mEnter.getTime() != b.mEnter.getTime()) { + return a.mEnter.getTime() < b.mEnter.getTime(); // time of flight, so order from interaction region + } + if (a.getDetector() != b.getDetector()) { + return a.getDetector() < b.getDetector(); + } + + if (a.getPadRow() != b.getPadRow()) { + return a.getPadRow() < b.getPadRow(); + } + + if (a.getPadCol() != b.getPadCol()) { + return a.getPadCol() < b.getPadCol(); + } + /*if (a.mTrackId != b.mTrackId) { + return a.mTrackId < b.mTrackId; + } + + if (a.mDet>>2 != b.mDet>>2) { + return (a.mDet>>2) < (b.mDet>>2); // this gets the detector/stack/layer + } + + if (a.mEnter.getTime() != b.mEnter.getTime()) { + return a.mEnter.getTime() < b.mEnter.getTime(); // time of flight, so order from interaction region + } +*/ + return true; +} + void RawDataSpan::sort() { std::stable_sort(std::begin(digits), std::end(digits), comp_digit); std::stable_sort(std::begin(tracklets), std::end(tracklets), comp_tracklet); std::stable_sort(std::begin(hits), std::end(hits), comp_spacepoint); + std::stable_sort(std::begin(trackrefsegments), std::end(trackrefsegments), comp_trackrefshits); } template @@ -104,6 +156,10 @@ std::vector RawDataSpan::iterateBy() { // an map for keeping track which ranges correspond to which key std::map spanmap; + std::vector foundkeys; + + //TODO come back and try assign the trackrefhits to particular keys. + //for now we just add all of them to each key :-( // sort digits and tracklets sort(); @@ -111,7 +167,8 @@ std::vector RawDataSpan::iterateBy() // add all the digits to a map for (auto cur = digits.begin(); cur != digits.end(); /* noop */) { // calculate the key of the current (first unprocessed) digit - auto key = keyfunc::key(*cur); + auto key = keyfunc::key(*cur,true); + foundkeys.push_back(key); // find the first digit with a different key auto nxt = std::find_if(cur, digits.end(), [key](auto x) { return keyfunc::key(x) != key; }); // store the range cur:nxt in the map @@ -123,6 +180,7 @@ std::vector RawDataSpan::iterateBy() // add tracklets to the map for (auto cur = tracklets.begin(); cur != tracklets.end(); /* noop */) { auto key = keyfunc::key(*cur); + foundkeys.push_back(key); auto nxt = std::find_if(cur, tracklets.end(), [key](auto x) { return keyfunc::key(x) != key; }); spanmap[key].tracklets = boost::make_iterator_range(cur, nxt); cur = nxt; @@ -134,12 +192,12 @@ std::vector RawDataSpan::iterateBy() std::map::iterator> firsthit; for (auto cur = hits.begin(); cur != hits.end(); ++cur) { // calculate the keys for this hit - auto keys = keyfunc::keys(*cur); + auto keys = keyfunc::keys(*cur,true); // if we are not yet aware of this key, register the current hit as the first hit for (auto key : keys) { firsthit.insert({key, cur}); } - // remote the keys from the firsthit map that are no longer found in the hits + // remove the keys from the firsthit map that are no longer found in the hits for (auto it = firsthit.cbegin(); it != firsthit.cend(); /* no increment */) { if (keys.find(it->first) == keys.end()) { spanmap[it->first].hits = boost::make_iterator_range(it->second, cur); @@ -149,6 +207,40 @@ std::vector RawDataSpan::iterateBy() } } } + // trackreferences underpinning the tracksegments are in time order, match the times with the + // given a key, find the first and last tracklet, take the first and last time. + // advance to the first time in the trackrefsegments + auto uniquekeys=std::unique(foundkeys.begin(),foundkeys.end()); + foundkeys.erase(uniquekeys,foundkeys.end()); + for( auto &key : foundkeys) { + //spanmap[key].trackrefsegments= boost::make_iterator_range(trackrefsegments.begin(),trackrefsegments.end()); + // std::cout << "key : " << key << "\n"; + } + int count=0; + for (auto cur = trackrefsegments.begin(); cur != trackrefsegments.end(); /* noop */) { + //LOGP(info,"********************** tracksegment start: ***********************************"); + //std::cout << "index: " << count++ << "\n"; + + + auto key = keyfunc::key(*cur,true); + auto it= std::find(foundkeys.begin(),foundkeys.end(),key); + if(it == foundkeys.end()){ + LOGP(info, "tracksegment key is not present so error key:{}", key); + + } + else{ + LOGP(info, "tracksegment key is present so error key:{}", key); + } + + //LOGP(info,"!!!!!!!!!!!!!!!!!!!!!! tracksegment find_if: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"); + //std::cout << "!!!!!!!!!!!!!!!!!!!!!! cur " << *cur << "\n"; + auto nxt = std::find_if(cur, trackrefsegments.end(), [key](auto x) { return keyfunc::key(x) != key; }); + //std::cout << "!!!!!!!!!!!!!!!!!!!!!! nxt " << *nxt << "\n"; + //LOGP(info,"!!!!!!!!!!!!!!!!!!!!!! tracksegment find_if: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"); + spanmap[key].trackrefsegments = boost::make_iterator_range(cur, nxt); + cur = nxt; + //LOGP(info,"********************** tracksegment end: ***********************************"); + } // convert the map of spans into a vector, as we do not need the access by key // and longer, and having a vector makes the looping by the user easier. @@ -164,12 +256,12 @@ std::vector RawDataSpan::iterateBy() struct PadRowID { /// The static `key` method calculates a padrow ID for digits and tracklets template - static uint32_t key(const T& x) + static uint32_t key(const T& x, bool print=false) { return 100 * x.getDetector() + x.getPadRow(); } - static std::set keys(const o2::trd::ChamberSpacePoint& x) + static std::set keys(const o2::trd::ChamberSpacePoint& x, bool print=false) { uint32_t key = 100 * x.getDetector() + x.getPadRow(); return {key}; @@ -191,24 +283,31 @@ std::vector RawDataSpan::iterateByPadRow() { return iterateBy - static uint32_t key(const T& x) + static uint32_t key(const T& x, bool print=false) { + if(print)std::cout << " calculate key : " << 1000 * x.getDetector() + 8 * x.getPadRow() + 4 * (x.getROB() % 2) + x.getMCM() % 4 << " for det:" << x.getDetector() << " padrow:" << x.getPadRow() << " rob:" << x.getROB() << " mcm:" << x.getMCM() << std::endl; return 1000 * x.getDetector() + 8 * x.getPadRow() + 4 * (x.getROB() % 2) + x.getMCM() % 4; } - static std::set keys(const o2::trd::ChamberSpacePoint& x) + static std::set keys(const o2::trd::ChamberSpacePoint& x, bool print=false) { uint32_t detrow = 1000 * x.getDetector() + 8 * x.getPadRow(); uint32_t mcmcol = uint32_t(x.getPadCol() / float(o2::trd::constants::NCOLMCM)); // float c = x.getPadCol() - float(mcmcol * o2::trd::constants::NCOLMCM); float c = x.getMCMChannel(mcmcol); + if(print)std::cout << " calculate key : " << 1000 * x.getDetector() + 8 * x.getPadRow() + 4 * (x.getROB() % 2) + x.getMCM() % 4 << " for det:" << x.getDetector() << " padrow:" << x.getPadRow() << " rob:" << x.getROB() << " mcm:" << x.getMCM() << std::endl; + if(print)std::cout << " calculates key : detrow : " << detrow << " mcmcol : " << mcmcol << std::endl; if (c >= 19.0 && mcmcol >= 1) { + + if(print)std::cout << "returning detrow+mcmcol-1 : " << detrow+mcmcol-1 << " detrow+mcmcol : " << detrow+mcmcol << std::endl; return {detrow + mcmcol - 1, detrow + mcmcol}; } else if (c <= 1.0 && mcmcol <= 6) { + if(print)std::cout << "returning detrow+mcmcol : " << detrow+mcmcol << " detrow+mcmcol+1 : " << detrow+mcmcol+1 << std::endl; return {detrow + mcmcol, detrow + mcmcol + 1}; } else { + if(print)std::cout << "returning detrow+mcmcol : " << detrow+mcmcol << std::endl; return {detrow + mcmcol}; } } @@ -261,7 +360,7 @@ std::vector RawDataSpan::makeMCTrackReferences() // the case of processing a whole event, the distinction by detector will be needed. std::map, TrackReferencesInfo> trackReferenceInfo; - for (int iTrackRef = 0; iTrackRef < trackrefs.size(); ++iTrackRef) { +/* for (int iTrackRef = 0; iTrackRef < trackrefs.size(); ++iTrackRef) { auto ref = trackrefs[iTrackRef]; // we look for track references classified as entering the drift region @@ -279,7 +378,7 @@ std::vector RawDataSpan::makeMCTrackReferences() } } } // trackreference loop - +*/ std::vector trackReferences; for (auto x : trackReferenceInfo) { auto trackid = x.first.first; @@ -311,7 +410,7 @@ std::vector RawDataSpan::makeMCTrackSegments() auto hit = hits[iHit]; // in the following, we will look for track segments using hits in the drift region - if (hit.isFromDriftRegion()) { +// if (hit.isFromDriftRegion()) { // The first hit is the hit closest to the anode region, i.e. with the largest x coordinate. auto id = std::make_pair(hit.getID(), hit.getDetector()); if (hit.getX() > trackSegmentInfo[id].start) { @@ -323,9 +422,8 @@ std::vector RawDataSpan::makeMCTrackSegments() trackSegmentInfo[id].lasthit = iHit; trackSegmentInfo[id].end = hit.getX(); } - } + // } } // hit loop - std::vector trackSegments; for (auto x : trackSegmentInfo) { auto trackid = x.first.first; @@ -337,6 +435,202 @@ std::vector RawDataSpan::makeMCTrackSegments() return trackSegments; } +std::vector RawDataSpan::makeMCTrackSegmentsGeant() +{ + std::vector trackSegments; +// go through the TrackletSegmentInfo and match on time and return the relevant segemnts to this time period. + std::cout << " in " << __func__ << " trackrefhits has : " << trackrefsegments.size() << " entries " << std::endl; + for(auto &segmentinfo : trackrefsegments ){ + //convert these to TrackSegments + + } + return trackSegments; + +} +std::vector RawDataSpan::makeMCTrackSegmentsEntryExit() +{ + // define a struct to keep track of the first and last MC hit of one track in one chamber + struct SegmentInfo { + //enering and exiting trackreferences. + size_t entertrackref{0}; // + size_t exittrackref{0}; + float start{-999.9}; // local x cordinate of the first trackref. + float end{999.9}; // local x cordinate of the last trackref. + }; + // Keep information about found track segments in a map indexed by track ID and detector number. + // If the span only covers (part of) a detector, the detector information is redundant, but in + // the case of processing a whole event, the distinction by detector will be needed. + std::map, SegmentInfo> trackRefSegmentInfo; + int trackrefcounter=0; + // std::cout << " trackrefs has : " << trackrefs.size() << " entries " << std::endl; + std::cout << " trackrefs has : " << trackrefsegments.size() << " entries " << std::endl; + + for (auto &trackref : trackrefsegments) { + // LOGP(info, "trackrefs det id : {} ", trackref.getDetectorId()); + // only copy the trd trackrefs. + // a TrackID is for a specific geant track, so it and detectorId uniquely identifies the trackref entry and exit points, defined below. + auto id = std::make_pair(trackref.mEnter.getTrackID(), trackref.mEnter.getUserId()>>2); + //case 0x1: + // direction entering + trackRefSegmentInfo[id].entertrackref=trackrefcounter; + // std::cout << " Now to look for corresponding hit of trackref entering" << std::endl; + for (int iHit = 0; iHit < hits.size(); ++iHit) { + auto hit = hits[iHit]; + float distance = sqrt(pow((hit.getX() - trackref.mEnter.X()),2) + pow(( hit.getY()- trackref.mEnter.Y()),2)+ pow( (hit.getZ() - trackref.mEnter.Z()),2)); + //std::cout << " distance of " << distance << " for " < rct; + Hit a=convertTrackReferenceToHit(trackref.mEnter,trackref.mEnter.getTrackID(),trackref.mEnter.getUserId()>>2); + //ChamberSpacePoint(; + // std::cout << " distance of " << distance << " for " <( + // trackref.X(), trackref.Y(), trackref.Z())); + // direction exiting + //std::cout << " Now to look for corresponding hit of trackref exiting" << std::endl; + for (int iHit = 0; iHit < hits.size(); ++iHit) { + auto hit = hits[iHit]; + float distance = sqrt(pow((hit.getX() - trackref.mExit.X()),2) + pow(( hit.getY()- trackref.mExit.Y()),2)+ pow( (hit.getZ() - trackref.mExit.Z()),2)); + // std::cout << " distance of " << distance << " for " < rct; + Hit a=convertTrackReferenceToHit(trackref.mEnter,trackref.mEnter.getTrackID(),trackref.mEnter.getUserId()>>2); + //std::cout << " distance of " << distance << " for " <( +// // trackref.X(), trackref.Y(), trackref.Z())); +// std::cout << " Now to look for corresponding hit of trackref entering" << std::endl; +// for (int iHit = 0; iHit < hits.size(); ++iHit) { +// auto hit = hits[iHit]; +// float distance = sqrt(pow((hit.getX() - trackref.getX()),2) + pow(( hit.getY()- trackref.getY()),2)+ pow( (hit.getZ() - trackref.getZ()),2)); +// std::cout << " distance of " << distance << " for " <( +// // trackref.X(), trackref.Y(), trackref.Z())); +// // direction exiting +// std::cout << " Now to look for corresponding hit of trackref exiting" << std::endl; +// for (int iHit = 0; iHit < hits.size(); ++iHit) { +// auto hit = hits[iHit]; +// float distance = sqrt(pow((hit.getX() - trackref.getX()),2) + pow(( hit.getY()- trackref.getY()),2)+ pow( (hit.getZ() - trackref.getZ()),2)); +// std::cout << " distance of " << distance << " for " < trackSegments; + + auto ctrans = o2::trd::CoordinateTransformer::instance(); + for (auto x : trackRefSegmentInfo) { + HitPoint enter,exit; + // enter=ctrans->MakeSpacePoint(trackrefs[x.second.entertrackref]); + // exit=ctrans->MakeSpacePoint(trackrefs[x.second.exittrackref]); + auto trackid = x.first.first; + auto detector = x.first.second; +// o2::trd::Hit entering= convertTrackReferenceToHit(enter,trackid,detector); +// o2::trd::Hit leaving= convertTrackReferenceToHit(exit,trackid,detector); + //find the corresponding hit to each. + //// this is massively expensive but will fix later. + auto firsthit = hits[x.second.entertrackref]; + auto lasthit = hits[x.second.exittrackref]; + trackSegments.emplace_back(firsthit, lasthit, trackid); + } + return trackSegments; +} + +Hit RawDataSpan::convertTrackReferenceToHit(o2::TrackReference &ref, + int trackid, int detector) +{ + //std::cout << __func__ << " " << ref << " trackid: " << trackid << " det: " << detector << std::endl; + std::cout << std::dec; + float tof = ref.getTime() * 1e6; // The time of flight in micro-seconds + double pos[3] = {ref.X(),ref.Y(),ref.Z()}; + double loc[3] = {-99, -99, -99}; + // gGeoManager->Export("geometry.root"); // is there a corresponding import ? + if(gGeoManager == nullptr){ + gGeoManager = new TGeoManager(); + gGeoManager->Import("o2sim_geometry.root"); + } + auto node = gGeoManager->FindNode(ref.X(), ref.Y(), ref.Z()); + gGeoManager->MasterToLocal(pos, loc); // Go to the local coordinate system (locR, locC, locT) + // if (!node) { + // std::cout << "node name : unknown\n"; + // } + auto userid=ref.getUserId(); + auto enterleave=userid&0x3; + // std::cout << "node name : " << node->GetVolume()->GetName() << " userid: " << userid << " ref enter/leave : " << enterleave << "\n"; + //std::cout << "node name : " << node->GetVolume()->GetName() << " ref enter/leave : " << (ref.getUserId()&0x3) << "\n"; + // find node + float locC = loc[0], locR = loc[1], locT = loc[2]; + locT = locT - 0.5 * (Geometry::drThick() + Geometry::amThick()); // Relative to middle of amplification region + float charge=0.;// nominal error charge as this is not a hit. + //std::cout << "converTrackReferenceToHit : " << pos[0] <<"," << pos[1] <<"," << pos[2] <<" :: " << locC << "," << locR << "," << locT << " local : " << ref.LocalX() << ":" << ref.LocalY() << std::endl; + /********/ + /*int sector = Geometry::getSector(ref.getUserId()>>2); + int layer = Geometry::getLayer(detector); + int stack = Geometry::getStack(detector); + float phi = 2.0 * TMath::Pi() / (float)o2::trd::constants::NSECTOR * ((float)sector + 0.5); + double loc1[3] = {-99, -99, -99}; + double pos1[3] = {ref.X(),ref.Y(),ref.Z()}; + loc1[0] = pos1[0] * TMath::Cos(phi) - pos1[1] * TMath::Sin(phi); + loc1[1] = -1.0*pos1[0] * TMath::Sin(phi) + pos1[1] * TMath::Cos(phi); + loc1[2] = pos1[2];*/ + // std::cout << fmt::format("converTrackReferenceToHit manual phi:{} : {:.5}:{:.5}:{:.5} :: {:.5}:{:.5}:{:.5}",phi,pos1[0],pos1[1],pos1[2],loc1[0],loc1[1],loc1[2]); +/* double rmin = mgeo->getTime0(layer); + double rmax = mgeo->getTime0(layer) - mgeo->drThick() - mgeo->amThick(); + double ymin = -mgeo->getChamberWidth(layer) / 2; + double ymax = mgeo->getChamberWidth(layer) / 2; + double zmin = -mgeo->getChamberLength(layer, stack) / 2; + double zmax = mgeo->getChamberLength(layer, stack) / 2; + */ + /********/ + + + o2::trd::Hit thehit(ref.X(),ref.Y(),ref.Z(),locC,locR,locT, tof, charge, trackid, detector,true); + return thehit; +} + /// The RawDataManager constructor: connects all data files and sets up trees, readers etc. RawDataManager::RawDataManager(std::filesystem::path dir) { @@ -380,7 +674,7 @@ RawDataManager::RawDataManager(std::filesystem::path dir) if (std::filesystem::exists(dir / "collisioncontext.root")) { TFile fInCollCtx((dir / "collisioncontext.root").c_str()); mCollisionContext = (o2::steer::DigitizationContext*)fInCollCtx.Get("DigitizationContext"); - // mCollisionContext->printCollisionSummary(); + mCollisionContext->printCollisionSummary(); } // We create the MC TTree using event header and tracks from the kinematics file @@ -389,9 +683,10 @@ RawDataManager::RawDataManager(std::filesystem::path dir) mMCFile->GetObject("o2sim", mMCTree); mMCTree->SetBranchAddress("MCEventHeader.", &mMCEventHeader); mMCTree->SetBranchAddress("MCTrack", &mMCTracks); - //mMCTree->SetBranchAddress("TrackRefs", &mMCTrackReferences); + mMCTree->SetBranchAddress("TrackRefs", &mMCTrackReferences); } - + //process trackreferences into MCTrackletSegmentInfo + processTrackReferences(); // We then add the TRD hits to the MC tree if (mMCFile && std::filesystem::exists(dir / "o2sim_HitsTRD.root")) { mMCTree->AddFriend("o2sim", (dir / "o2sim_HitsTRD.root").c_str()); @@ -399,6 +694,74 @@ RawDataManager::RawDataManager(std::filesystem::path dir) } } +void RawDataManager::processTrackReferences() +{ + //Trackreferences come in ordered by : trackid, then time. + std::map, MCTrackletSegmentInfo> trackreferences; + int nev = mMCTree->GetEntries(); + // take the incoming trackreferences, remove those not applicable to TRD, and order them by det, trackid, and time. + for (int iev = 0; iev < nev; ++iev) { + mMCTree->GetEvent(iev); // all trackreferences are in a single branch. + // sort track ref by time and enter before exit. + for (const auto &trackref : *mMCTrackReferences) { + if (trackref.getDetectorId() == 2) { +// LOGP(info, "trackrefs det id : {} trackid: {} time : {} ", trackref.getDetectorId(), trackref.getTrackID(),trackref.getTime()); + uint32_t trackid=trackref.getTrackID(); + uint32_t det; + auto key = std::make_pair(trackref.getUserId()>>2,trackref.getTrackID()); + // only copy the trd trackrefs. + // std::map inner; + // inner.insert(std::make_pair(det, trackref)); + // trackreferences.insert(trackref.getTrackID(), inner); + if(!trackreferences.contains(key)){ + trackreferences[key]=MCTrackletSegmentInfo(); // map[key]. trackrefsvectorwithextra.push_back(trackref); + if((trackref.getUserId() & 0x3) == 0x2) { + // we have an exit but no preceding enter ?? + std::cout << "Exit but no entry key : " << key.first << ":"<< key.second << " Setting entry point with entry point of : " << trackreferences[key].mEnter << " and corresponding exit of : " << trackreferences[key].mExit << std::endl; + } + } + switch (trackref.getUserId() & 0x3) { + case 0x1: + // direction entering + //trackreferences[key].setEntry(trackref.X(),trackref.Y(),trackref.Z(),trackid, det); + trackreferences[key].mEnter=trackref;//setEntry(trackref.X(),trackref.Y(),trackref.Z(),trackid, det); + //std::cout << " key : " << key.first << ":"<< key.second << " Setting entry point with entry point of : " << trackreferences[key].mEnter << " and corresponding exit of : " << trackreferences[key].mExit << std::endl; + break; + case 0x2: + // trackreferences[key].setExit(trackref.X(),trackref.Y(),trackref.Z(), trackid,det); + trackreferences[key].mExit=trackref; + //std::cout << " key : " << key.first << ":"<< key.second << " Setting exit point with exit point of : " << trackreferences[key].mExit << " and corresponding enter of : " << trackreferences[key].mEnter << std::endl; + // direction exiting + break; + } + } + } + } +// now walk through the map and insert the entry/exit ref points into std::vector mMCTrackletSegmentInfo{0}; + // Iterate using C++17 facilities + for (const auto& [key, value] : trackreferences){ + //std::cout << '[' << key << "] = " << value << "; "; + auto det = key.first; + auto trackid = key.second; + if(trackreferences[key].mEnter.getTrackID() == 0){ + std::cout << " something wrong, trackid for key.1 "<< key.first << " key.2" << key.second << " is : " << trackreferences[key].mEnter.getTrackID() << " " << trackreferences[key].mExit.getTrackID() << std::endl; + } + trackreferences[key].setMidPoint(); + //mMCTrackletSegmentInfo.emplace_back(trackreferences[key].mEnter,trackreferences[key].mExit,det,trackid); + if(trackreferences[key].isGood()){ + // remove those that dont have a enter or exit + mMCTrackletSegmentInfo.emplace_back(value); + // std::cout << "GOOD: " << trackreferences[key] << "\n"; + } + //else { + // std::cout << "BAD: " << trackreferences[key] << "\n"; + // } + //mMCTrackletSegmentInfo.emplace_back(value.mEnter,value.mExit,det,trackid); +// std::cout << "trackid for key.1 "<< key.first << " key.2" << key.second << "\n Enter: " << trackreferences[key].mEnter << "\n Exit: " << trackreferences[key].mExit << "\n"; +// std::cout << " det:" << trackreferences[key].getDetector() << " padrow:"<< trackreferences[key].getPadRow() <<" padcol:"<< trackreferences[key].getPadRow()<<" rob:"<< trackreferences[key].getROB()<<" mcm:"<< trackreferences[key].getMCM() << "\n"; + } +} + bool RawDataManager::nextTimeFrame() { if (!mDataTree->GetEntry(mTimeFrameNo)) { @@ -438,14 +801,21 @@ bool RawDataManager::nextEvent() mMCTree->GetEntry(i); // } - O2INFO("Loaded matching MC event #%d with time offset %f ns and %d hits", - i, mTriggerRecord.getBCData().differenceInBCNS(evrec), mHits->size()); + O2INFO("Loaded matching MC event #%d with time offset %f ns and %d hits and %d trackreferences ", + i, mTriggerRecord.getBCData().differenceInBCNS(evrec), mHits->size(), mMCTrackReferences->size()); + O2INFO("Loaded matching MC event #%d with time offset %f ns and %d hits and %d trackreferences ", + i, mTriggerRecord.getBCData().differenceInBCNS(evrec), mHits->size(), mMCTrackReferences->size()); // convert hits to spacepoints auto ctrans = o2::trd::CoordinateTransformer::instance(); for (auto& hit : *mHits) { - mHitPoints.emplace_back(ctrans->MakeSpacePoint(hit), hit.GetCharge()); + mHitPoints.emplace_back(ctrans->MakeSpacePoint(hit), hit.GetCharge(),hit.GetTrackID()); + // std::cout << " building mHitPoints trackid : " << hit.GetTrackID() << std::endl; } + /* for (auto& trackref : *mMCTrackReferences) { + std::cout << "added track ref at : " << trackref.X()<<":"<MakeSpacePoint(trackref)); + }*/ } } } @@ -462,8 +832,24 @@ RawDataSpan RawDataManager::getEvent() ev.tracklets = boost::make_iterator_range_n(mTracklets->begin() + mTriggerRecord.getFirstTracklet(), mTriggerRecord.getNumberOfTracklets()); ev.hits = boost::make_iterator_range(mHitPoints.begin(), mHitPoints.end()); - + + //ev.trackrefs = boost::make_iterator_range(mTrackReferences.begin(), mTrackReferences.end()); + //ev.trackrefhits = boost::make_iterator_range(mMCTrackReferences.begin(), mMCTrackReferences.end()); + ev.trackrefsegments = boost::make_iterator_range(mMCTrackletSegmentInfo.begin(), mMCTrackletSegmentInfo.end()); + auto evtime = getTriggerTime(); + std::cout << "event time : " << evtime << std::endl; + std::vector trackids; + for(auto hit : ev.hits){ + trackids.push_back(hit.getTrackID()); + + } + int idcountpre=trackids.size(); + sort(trackids.begin(),trackids.end()); + auto uniqueids=std::unique(trackids.begin(),trackids.end()); + trackids.erase(uniqueids,trackids.end()); + int idcountpost=trackids.size(); + std::cout << " pre trackids:" << idcountpre << " post trackids:" << idcountpost << std::endl; // if (tpctracks) { // for (auto &track : *mTpcTracks) { @@ -494,7 +880,7 @@ RawDataSpan RawDataManager::getEvent() // ev.trackpoints.begin() = ev.evtrackpoints.begin(); // ev.trackpoints.end() = ev.evtrackpoints.end(); - + std::cout << " returning event" << std::endl; return ev; } @@ -561,3 +947,83 @@ std::string RawDataManager::describeEvent() << mTriggerRecord.getNumberOfTracklets() << " tracklets"; return out.str(); } + + + + +int MCTrackletSegmentInfo::getPadRow() const +{ + auto enterlocalx=mEnter.LocalX(); + auto exitlocalx=mExit.LocalX(); + //: mID(id), mDetector(detector), mX(x), mY(y), mZ(z), mPadrow(rct[0]), mPadcol(rct[1]), mTimebin(rct[2]), mInDrift(inDrift){}; + double pos[3] = {mEnter.X(),mEnter.Y(),mEnter.Z()}; + double posexit[3] = {mExit.X(),mExit.Y(),mExit.Z()}; + double loc1[3] = {-99, -99, -99}; + double loc[3] = {-99, -99, -99}; + // gGeoManager->Export("geometry.root"); // is there a corresponding import ? + if(gGeoManager == nullptr){ + gGeoManager = new TGeoManager(); + gGeoManager->Import("o2sim_geometry.root"); + } + auto node = gGeoManager->FindNode(mMidPoint[0],mMidPoint[1],mMidPoint[2]); + gGeoManager->MasterToLocal(posexit, loc); // Go to the local coordinate system (locR, locC, locT) + //std::cout <<__func__ <<"aapos : " << pos[0] << ":"<GetVolume()->GetName() << "\n"; + + return (int)loc[1]; +}; + +int MCTrackletSegmentInfo::getPadCol() const +{ + double pos[3] = {mEnter.X(),mEnter.Y(),mEnter.Z()}; + double posexit[3] = {mExit.X(),mExit.Y(),mExit.Z()}; + double loc[3] = {-99, -99, -99}; + // gGeoManager->Export("geometry.root"); // is there a corresponding import ? + if(gGeoManager == nullptr){ + gGeoManager = new TGeoManager(); + gGeoManager->Import("o2sim_geometry.root"); + } + auto node = gGeoManager->FindNode(mMidPoint[0],mMidPoint[1],mMidPoint[2]); + gGeoManager->MasterToLocal(posexit, loc); // Go to the local coordinate system (locR, locC, locT) + float locC = loc[0], locR = loc[1], locT = loc[2]; + //std::cout <<__func__ <<"pos : " << pos[0] << ":"<GetVolume()->GetName() << "\n"; + return (int)locC; +} + +float MCTrackletSegmentInfo::getMCMf() const +{ + return o2::trd::HelperMethods::getMCMfromPad(getPadRow(),getPadCol()); +}; + +float MCTrackletSegmentInfo::getROBf() const +{ + return o2::trd::HelperMethods::getROBfromPad(getPadRow(),getPadCol()); +}; + +int MCTrackletSegmentInfo::getMCM() const +{ + auto padrow = getPadRow(); + auto padcol = getPadCol(); + //std::cout << __func__ << " padrow:" << padrow << "[0:"<< constants::NROWC1<<"] padcol:" << padcol << "[0:"<< constants::NCOLUMN << "] \n"; + //std::cout <<__func__ << " midpoint:" << mMidPoint[0] << ":" << mMidPoint[1] <<":" << mMidPoint[2] <<" " << gGeoManager->FindNode(mMidPoint[0],mMidPoint[1],mMidPoint[2])->GetVolume()->GetName() << "\n"; + return o2::trd::HelperMethods::getMCMfromPad(getPadRow(),getPadCol()); +}; + +int MCTrackletSegmentInfo::getROB() const +{ + return o2::trd::HelperMethods::getROBfromPad(getPadRow(),getPadCol()); +}; + + +namespace o2::trd +{ + +std::ostream& operator<<(std::ostream& os, const MCTrackletSegmentInfo& p) +{ + int sector = p.getDetector() / 30; + int stack = (p.getDetector() % 30) / 6; + int layer = p.getDetector() % 6; + os << fmt::format("mctrackletsegmentinfo : ({:.3}:{:.3}:{:.3}--{:.3}:{:.3}:{:.3} : mid:{:.3}:{:.3}:{:.3})",p.mEnter.X(),p.mEnter.Y(),p.mEnter.Z(),p.mExit.X(),p.mExit.Y(),p.mExit.Z(), p.mMidPoint[0],p.mMidPoint[1],p.mMidPoint[2] ); + return os; +} + +} diff --git a/Detectors/TRD/qc/src/RawDisplay.cxx b/Detectors/TRD/qc/src/RawDisplay.cxx index da1d0b115be45..491561ade7ea3 100644 --- a/Detectors/TRD/qc/src/RawDisplay.cxx +++ b/Detectors/TRD/qc/src/RawDisplay.cxx @@ -33,11 +33,11 @@ float PadColF(o2::trd::Tracklet64& tracklet) float padLocal = tracklet.getPositionBinSigned() * constants::GRANULARITYTRKLPOS; // MCM number in column direction (0..7) int mcmCol = (tracklet.getMCM() % constants::NMCMROBINCOL) + constants::NMCMROBINCOL * (tracklet.getROB() % 2); - // original calculation // FIXME: understand why the offset seems to be 6 pads and not nChannels / 2 = 10.5 // return CAMath::Round(6.f + mcmCol * ((float)constants::NCOLMCM) + padLocal); + //std::cout << "tracklet64 cal: " << std::dec << tracklet.getPadCol() << " :: " << float((mcmCol + 1) * constants::NCOLMCM) + padLocal - 10.0 << std::endl; // my calculation return float((mcmCol + 1) * constants::NCOLMCM) + padLocal - 10.0; } @@ -87,7 +87,7 @@ MCMDisplay::MCMDisplay(RawDataSpan& mcmdata, TVirtualPad* pad) mPad->SetTitle(mDesc.c_str()); } - mDigitsHisto = new TH2F(mName.c_str(), (mDesc + ";pad;time bin").c_str(), (mLastPad - mFirstPad), mFirstPad, mLastPad, 30, 0., 30.); + mDigitsHisto = new TH2F(mName.c_str(), (mDesc + ";pad;time bin").c_str(), (mLastPad - mFirstPad), mFirstPad, mLastPad, 36, -6., 30.); for (auto digit : mDataSpan.digits) { auto adc = digit.getADC(); @@ -150,7 +150,7 @@ void RawDisplay::drawClusters() // << " pos = " << pos << " ~ " << clpos // << endl; clustermarker.DrawMarker(clpos, t - 0.5); - // cogmarker.DrawMarker(cog, t - 0.5); + cogmarker.DrawMarker(cog, t - 0.5); } } } @@ -168,8 +168,8 @@ void RawDisplay::drawHits() } } } -/* -void RawDisplay::drawMCTrackReferences() + +/*void RawDisplay::drawMCTrackReferences() { TLine line; line.SetLineColor(kMagenta); @@ -179,15 +179,39 @@ void RawDisplay::drawMCTrackReferences() line.DrawLine(trkref.getStartPoint().getPadCol(), trkref.getStartPoint().getTimeBin(), trkref.getEndPoint().getPadCol(), trkref.getEndPoint().getTimeBin()); } } - +*/ void RawDisplay::drawMCTrackSegments() { TLine line; - line.SetLineColor(kBlue); + line.SetLineColor(kBlack); line.SetLineWidth(2.0); - for (auto& trkl : mDataSpan.makeMCTrackSegments()) { - line.DrawLine(trkl.getStartPoint().getPadCol(), trkl.getStartPoint().getTimeBin(), trkl.getEndPoint().getPadCol(), trkl.getEndPoint().getTimeBin()); + for (auto& trkl : mDataSpan.trackrefsegments) { + line.DrawLine(trkl.mEnter.X(), trkl.mEnter.Y(),trkl.mExit.X(),trkl.mExit.Y()); + } +} + + +void RawDisplay::drawMCTrackRefEntryExit() +{ + TLine line; + line.SetLineColor(kBlack); + line.SetLineWidth(2.0); + TMarker hitmarker; + hitmarker.SetMarkerColor(kBlack); + hitmarker.SetMarkerStyle(38); + std::cout << "About to display the geant entry and exit point !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl; + int count=0; + for (auto& geanttrkl : mDataSpan.makeMCTrackSegmentsGeant()) { + count++; + } + std::cout << "makeMCTrackSegmentsGeant created " << count << " items" << std::endl; + + for (auto& geanttrkl : mDataSpan.makeMCTrackSegmentsEntryExit()) { + hitmarker.SetMarkerSize(log10(10)); + hitmarker.DrawMarker(geanttrkl.getStartPoint().getPadCol(), geanttrkl.getStartPoint().getTimeBin()); + //std::cout << "entry at :"<< geanttrkl.getStartPoint().getPadCol() <<":"<< geanttrkl.getStartPoint().getTimeBin() << std::endl; + //std::cout << "exit at :" <