From 43a05d13a9d62309d9ed081559530f1ed7836894 Mon Sep 17 00:00:00 2001 From: Daniel Brandenburg Date: Fri, 28 Jun 2024 17:46:39 -0400 Subject: [PATCH 1/3] StFwdQAMaker for optional QA of fwd systems and tracking --- StRoot/StFwdTrackMaker/StFwdQAMaker.cxx | 450 ++++++++++++++++++++++++ StRoot/StFwdTrackMaker/StFwdQAMaker.h | 396 +++++++++++++++++++++ 2 files changed, 846 insertions(+) create mode 100644 StRoot/StFwdTrackMaker/StFwdQAMaker.cxx create mode 100644 StRoot/StFwdTrackMaker/StFwdQAMaker.h diff --git a/StRoot/StFwdTrackMaker/StFwdQAMaker.cxx b/StRoot/StFwdTrackMaker/StFwdQAMaker.cxx new file mode 100644 index 00000000000..b4e13b6f5e9 --- /dev/null +++ b/StRoot/StFwdTrackMaker/StFwdQAMaker.cxx @@ -0,0 +1,450 @@ +#include "StFwdTrackMaker/StFwdQAMaker.h" +#include "StFwdQAMaker.h" +#include "St_base/StMessMgr.h" +#include "StBFChain/StBFChain.h" +#include "StFwdTrackMaker/StFwdTrackMaker.h" + +#include "StFwdTrackMaker/include/Tracker/FwdTracker.h" +#include "StFwdTrackMaker/include/Tracker/ObjExporter.h" +// StEvent includes +#include "StEvent/StBTofCollection.h" +#include "StEvent/StBTofHeader.h" +#include "StEvent/StEvent.h" +#include "StEvent/StFttCluster.h" +#include "StEvent/StFttCollection.h" +#include "StEvent/StFcsCluster.h" +#include "StEvent/StFcsCollection.h" +#include "StFcsDbMaker/StFcsDb.h" +#include "StRoot/StEpdUtil/StEpdGeom.h" +#include "StEvent/StFwdTrackCollection.h" +#include "StEvent/StFwdTrack.h" + + +#include "StMuDSTMaker/COMMON/StMuDstMaker.h" +#include "StMuDSTMaker/COMMON/StMuDst.h" +#include "StMuDSTMaker/COMMON/StMuFstCollection.h" +#include "StMuDSTMaker/COMMON/StMuFstHit.h" +#include "StMuDSTMaker/COMMON/StMuPrimaryVertex.h" +#include "StMuDSTMaker/COMMON/StMuFwdTrack.h" +#include "StMuDSTMaker/COMMON/StMuFwdTrackCollection.h" +#include "StMuDSTMaker/COMMON/StMuFcsCollection.h" +#include "StMuDSTMaker/COMMON/StMuFcsCluster.h" +#include "StMuDSTMaker/COMMON/StMuFcsHit.h" + + +/** Clear the FwdTreeData from one event to next */ +void FwdTreeData::clear(){ + header.clear(); + ftt.reset(); + fttClusters.reset(); + fst.reset(); + reco.reset(); + seeds.reset(); + wcal.reset(); + hcal.reset(); +} + +/** @brief Sets the Tree cluster from StFcsCluster +* @param clu - StFcsCluster +* @param fcsDb - StFcsDb +* @return void +*/ +void FwdTreeFcsCluster::set( StFcsCluster *clu, StFcsDb * fcsDb ) { + if (!clu) { + mId = -1; + return; + } + if ( fcsDb ){ + StThreeVectorD xyz = fcsDb->getStarXYZfromColumnRow(clu->detectorId(),clu->x(),clu->y()); + pos.SetXYZ( xyz.x(), xyz.y(), xyz.z() ); + } else { + pos.SetXYZ( 0, 0, 0 ); + } + + mDetectorId = clu->detectorId(); + + mEnergy = clu->energy(); + mNTowers = clu->nTowers(); + mId = clu->id(); + mCategory = clu->category(); + mX = clu->x(); + mY = clu->y(); + mSigmaMin = clu->sigmaMin(); + mSigmaMax = clu->sigmaMax(); + mTheta = clu->theta(); + mChi2Ndf1Photon = clu->chi2Ndf1Photon(); + mChi2Ndf2Photon = clu->chi2Ndf2Photon(); +} + +void FwdTreeFcsCluster::set( StMuFcsCluster *clu, StFcsDb * fcsDb ) { + if (!clu) { + mId = -1; + return; + } + if ( fcsDb ){ + StThreeVectorD xyz = fcsDb->getStarXYZfromColumnRow(clu->detectorId(),clu->x(),clu->y()); + pos.SetXYZ( xyz.x(), xyz.y(), xyz.z() ); + } else { + pos.SetXYZ( 0, 0, 0 ); + } + + mDetectorId = clu->detectorId(); + mEnergy = clu->energy(); + mNTowers = clu->nTowers(); + mId = clu->id(); + mCategory = clu->category(); + mX = clu->x(); + mY = clu->y(); + mSigmaMin = clu->sigmaMin(); + mSigmaMax = clu->sigmaMax(); + mTheta = clu->theta(); + mChi2Ndf1Photon = clu->chi2Ndf1Photon(); + mChi2Ndf2Photon = clu->chi2Ndf2Photon(); +} + +void FwdTreeRecoTrack::set( StMuFwdTrack *muFwd ){ + if ( !muFwd ){ + LOG_DEBUG << "Invalid muFwdTrack found, skipping FwdTreeRecoTrack::set" << endm; + return; + } + + if (muFwd->didFitConverge()) status = 2; + if (muFwd->didFitConvergeFully()) status = 1; + nFailedPoints = muFwd->numberOfFailedPoints(); + q = muFwd->charge(); + mom = muFwd->momentum(); + projs.clear(); + seeds.clear(); + for ( auto p : muFwd->mProjections ){ + FwdTreeTrackProjection tp; + tp.set( p.mDetId, p.mXYZ, p.mMom, p.mCov ); + projs.push_back( tp ); + } + for ( auto s : muFwd->mFSTPoints ){ + FwdTreeHit h; + h.set( s.mXYZ.X(), s.mXYZ.Y(), s.mXYZ.Z(), 1, 1 ); + seeds.push_back( h ); + } + for ( auto s : muFwd->mFTTPoints ){ + FwdTreeHit h; + h.set( s.mXYZ.X(), s.mXYZ.Y(), s.mXYZ.Z(), 1, 1 ); + seeds.push_back( h ); + } +} + +void FwdTreeRecoTrack::set( StFwdTrack *fwd ){ + if ( !fwd ){ + LOG_DEBUG << "Invalid fwdTrack found, skipping FwdTreeRecoTrack::set" << endm; + return; + } + + if (fwd->didFitConverge()) status = 2; + if (fwd->didFitConvergeFully()) status = 1; + nFailedPoints = fwd->numberOfFailedPoints(); + mChi2 = fwd->chi2(); + q = fwd->charge(); + mom.SetXYZ( fwd->momentum().x(), fwd->momentum().y(), fwd->momentum().z() ); + projs.clear(); + seeds.clear(); + for ( auto p : fwd->mProjections ){ + FwdTreeTrackProjection tp; + tp.set( p.mDetId, TVector3(p.mXYZ.x(), p.mXYZ.y(), p.mXYZ.z()), TVector3(p.mMom.x(), p.mMom.y(), p.mMom.z()), p.mCov ); + projs.push_back( tp ); + } + for ( auto s : fwd->mFSTPoints ){ + FwdTreeHit h; + h.set( s.mXYZ.x(), s.mXYZ.y(), s.mXYZ.z(), 1, 1 ); + seeds.push_back( h ); + } + for ( auto s : fwd->mFTTPoints ){ + FwdTreeHit h; + h.set( s.mXYZ.x(), s.mXYZ.y(), s.mXYZ.z(), 1, 1 ); + seeds.push_back( h ); + } +} + + +StFwdQAMaker::StFwdQAMaker() : StMaker("fwdQAMaker"), mTreeFile(nullptr), mTree(nullptr) { + +} + +int StFwdQAMaker::Init() { + + mTreeFile = new TFile("fwdtree.root", "RECREATE"); + mTree = new TTree("fwd", "fwd tracking tree"); + + mTree->Branch("header", &mTreeData. header, 3200, 99 ); + mTree->Branch("nSeedTracks", &mTreeData.nSeedTracks, "nSeedTracks/I"); + mTreeData.fst.createBranch(mTree, "fst"); + mTreeData.ftt.createBranch(mTree, "ftt"); + mTreeData.fttClusters.createBranch(mTree, "fttClusters"); + mTreeData.wcal.createBranch(mTree, "wcalClusters"); + mTreeData.hcal.createBranch(mTree, "hcalClusters"); + + mTreeData.reco.createBranch(mTree, "reco"); + mTreeData.seeds.createBranch(mTree, "seeds"); + return kStOk; +} +int StFwdQAMaker::Finish() { + + if ( mTreeFile && mTree ){ + mTreeFile->cd(); + mTree->Write(); + mTreeFile->Write(); + LOG_DEBUG << "StFwdQA File written" << endm; + } + return kStOk; +} +int StFwdQAMaker::Make() { + LOG_DEBUG << "SETUP START" << endm; + // setup the datasets / makers + mMuDstMaker = (StMuDstMaker *)GetMaker("MuDst"); + if(mMuDstMaker) { + mMuDst = mMuDstMaker->muDst(); + mMuForwardTrackCollection = mMuDst->muFwdTrackCollection(); + mMuFcsCollection = mMuDst->muFcsCollection(); + if (mMuForwardTrackCollection){ + LOG_DEBUG << "Number of StMuFwdTracks: " << mMuForwardTrackCollection->numberOfFwdTracks() << endm; + } + } else { + LOG_DEBUG << "No StMuDstMaker found: " << mMuDstMaker << endm; + } + mFcsDb = static_cast(GetDataSet("fcsDb")); + + mFwdTrackMaker = (StFwdTrackMaker*) GetMaker( "fwdTrack" ); + if (!mFwdTrackMaker) { + LOG_WARN << "No StFwdTrackMaker found, skipping StFwdQAMaker" << endm; + // return kStOk; + } + + // Event header info from stevent + StEvent *mStEvent = static_cast(GetInputDS("StEvent")); + LOG_DEBUG << "SETUP COMPLETE" << endm; + // Get the primary vertex used by the FWD Tracker + auto eventPV = mFwdTrackMaker->GetEventPrimaryVertex(); + mTreeData.header.set( 0, 0, 0, eventPV ); + + // Fill Header from StEvent + if ( mMuDstMaker ){ + + } + else if ( mStEvent ){ + StBTofCollection *btofC = mStEvent->btofCollection(); + if (btofC) { + StBTofHeader * btofHeader = btofC->tofHeader(); + if (btofHeader){ + + int nEast = btofHeader->numberOfVpdHits( east ); + int nWest = btofHeader->numberOfVpdHits( west ); + int nTof = btofC->tofHits().size(); + float vpdVz = btofHeader->vpdVz(); + + mTreeData.header.set( 1, GetEventNumber(), nTof, eventPV ); + mTreeData.header.vpdVz = vpdVz; + } + } // btofC != nullptr + } + LOG_DEBUG << "HEADER COMPLETE" << endm; + + FwdTreeHit fh; + + auto fttHits = mFwdTrackMaker -> GetFttHits(); + for ( auto h : fttHits ){ + fh.set( h.getX(), h.getY(), h.getZ(), 1, 2 ); + mTreeData.ftt.add( fh ); + } + LOG_DEBUG << "FTT COMPLETE" << endm; + auto fstHits = mFwdTrackMaker -> GetFstHits(); + for ( auto h : fstHits ){ + fh.set( h.getX(), h.getY(), h.getZ(), 1, 1 ); + mTreeData.fst.add( fh ); + } + LOG_DEBUG << "FST COMPLETE" << endm; + auto seeds = mFwdTrackMaker -> getTrackSeeds(); + LOG_DEBUG << "Number of Seeds to save: " << seeds.size() << endm; + if (seeds.size() > 5000){ + LOG_WARN << "Number of seeds is too large, truncating to 5000" << endm; + // seeds.resize(1000); + } + size_t iSeed = 0; + for ( auto s : seeds ){ + for ( auto h : s ){ + fh.set( h->getX(), h->getY(), h->getZ(), 1, 1 ); + fh.trackId = iSeed; + mTreeData.seeds.add( fh ); + } + iSeed++; + if (iSeed > 5000) break; + } // for s in seeds + LOG_DEBUG << "SEEDS COMPLETE, saved " << iSeed << " seed hits" << endm; + mTreeData.nSeedTracks = (int)iSeed; + + size_t iTrack = 0; + FwdTreeRecoTrack frt; + FwdTreeHit fth; + if ( mMuForwardTrackCollection ){ + LOG_DEBUG << "Adding " << mMuForwardTrackCollection->numberOfFwdTracks() << " FwdTracks (MuDst)" << endm; + for ( size_t iTrack = 0; iTrack < mMuForwardTrackCollection->numberOfFwdTracks(); iTrack++ ){ + auto muTrack = mMuForwardTrackCollection->getFwdTrack(iTrack); + frt.set( muTrack ); + mTreeData.reco.add( frt ); + if ( iTrack > 5000 ) break; + } + } + + if ( mStEvent && mStEvent->fwdTrackCollection() ){ + iTrack = 0; + auto fwdTracks = mStEvent->fwdTrackCollection()->tracks(); + LOG_DEBUG << "Adding " << fwdTracks.size() << " FwdTracks (StEvent)" << endm; + for ( auto ft : fwdTracks ){ + frt.set( ft ); + mTreeData.reco.add( frt ); + if ( iTrack > 5000 ) break; + iTrack++; + } + } + LOG_DEBUG << "TRACKS COMPLETE" << endm; + + + FillFttClusters(); + if ( mMuDstMaker ){ + FillFcsStMuDst(); + } else { + FillFcsStEvent(); + } + + mTree->Fill(); + return kStOk; +} +void StFwdQAMaker::Clear(const Option_t *opts) { + mTreeData.clear(); + return; +} + +void StFwdQAMaker::FillFcsStMuDst( ) { + + if ( !mMuDst ){ + LOG_DEBUG << "No mMuDst found, skipping StFwdQAMaker::FillFcsStEvent" << endm; + return; + } + StMuFcsCollection* fcs = mMuDst->muFcsCollection(); + if ( !fcs ){ + LOG_DEBUG << "No muFcsCollection found, skipping StFwdQAMaker::FillFcsStEvent" << endm; + return; + } + + StEpdGeom epdgeo; + FwdTreeFcsCluster fcsclu; + // LOAD ECAL / HCAL CLUSTERS + for( size_t i = 0; i < fcs->numberOfClusters(); i++){ + StMuFcsCluster * clu = fcs->getCluster(i); + LOG_DEBUG << "FCS CLUSTERS: " << fcs->numberOfClusters() << endm; + fcsclu.set( clu, mFcsDb ); + + if ( clu->detectorId() == kFcsEcalNorthDetId || clu->detectorId() == kFcsEcalSouthDetId ){ + LOG_DEBUG << "Adding WCAL Cluster to FwdTree" << endm; + mTreeData.wcal.add( fcsclu ); + } else if ( clu->detectorId() == kFcsHcalNorthDetId || clu->detectorId() == kFcsHcalSouthDetId ){ + LOG_DEBUG << "Adding HCAL Cluster to FwdTree" << endm; + mTreeData.hcal.add( fcsclu ); + } + } +} + +void StFwdQAMaker::FillFcsStEvent( ) { + if ( !mStEvent ){ + LOG_DEBUG << "No StEvent found, skipping StFwdQAMaker::FillFcsStEvent" << endm; + return; + } + StFcsCollection* fcsCol = mStEvent->fcsCollection(); + if ( !fcsCol ){ + LOG_DEBUG << "No StFcsCollection found, skipping StFwdQAMaker::FillFcsStEvent" << endm; + return; + } + + StEpdGeom epdgeo; + + FwdTreeFcsCluster fcsclu; + // LOAD ECAL / HCAL CLUSTERS + for ( int idet = 0; idet < 4; idet++ ){ + StSPtrVecFcsCluster& clusters = fcsCol->clusters(idet); + LOG_DEBUG << "FCS CLUSTERS [" << idet << "]: " << clusters.size() << endm; + int nc=fcsCol->numberOfClusters(idet); + for ( int i = 0; i < nc; i++ ){ + StFcsCluster* clu = clusters[i]; + + fcsclu.set( clu, mFcsDb ); + + if ( clu->detectorId() == kFcsEcalNorthDetId || clu->detectorId() == kFcsEcalSouthDetId ){ + LOG_DEBUG << "Adding WCAL Cluster to FwdTree" << endm; + mTreeData.wcal.add( fcsclu ); + } else if ( clu->detectorId() == kFcsHcalNorthDetId || clu->detectorId() == kFcsHcalSouthDetId ){ + LOG_DEBUG << "Adding HCAL Cluster to FwdTree" << endm; + mTreeData.hcal.add( fcsclu ); + } + + } // i + } // idet + +} + +void StFwdQAMaker::FillFttClusters(){ + float pz[] = {280.90499, 303.70498, 326.60501, 349.40499}; + float SCALE = 1.0; + TVector3 cp; + FwdTreeHit fh; + FwdTreeFttCluster ftc; + StEvent *event = static_cast(GetInputDS("StEvent")); + if ( !event || !event->fttCollection() ) return; + + LOG_DEBUG << "FTT RawHits: " << event->fttCollection()->numberOfRawHits() << endm; + LOG_DEBUG << "FTT Clusters: " << event->fttCollection()->numberOfClusters() << endm; + + for ( size_t i = 0; i < event->fttCollection()->numberOfClusters(); i++ ){ + StFttCluster* c = event->fttCollection()->clusters()[i]; + if ( c->nStrips() < 1 ) continue; + float dw = 0.05, dlh = 60.0, dlv = 60.0; + float mx = 0.0, my = 0.0; + float sx = 1.0, sy = 1.0; + + + if ( c->quadrant() == kFttQuadrantA ){ + mx = 0; my = 0; + sx = 1.0; sy = 1.0; + } else if ( c->quadrant() == kFttQuadrantB ){ + mx = 10.16*SCALE; my = 0.0*SCALE; + sy = -1; + dlv = -dlv; + + } else if ( c->quadrant() == kFttQuadrantC ){ + mx = -10.16*SCALE ; my = -00.0*SCALE; + sx = -1.0; sy = -1.0; + dlh = -dlh; dlv = -dlv; + + } else if ( c->quadrant() == kFttQuadrantD ){ + sx = -1; + dlh = -dlh; + } + + cp.SetZ( -pz[ c->plane() ] * SCALE ); + if ( c->orientation() == kFttHorizontal ){ + cp.SetY( my + sy * c->x()/10.0 * SCALE ); + cp.SetX( mx ); + } else if ( c->orientation() == kFttVertical ){ + cp.SetX( mx + sx * c->x()/10.0 * SCALE ); + cp.SetY( my ); + } + // fh.set( cp.X(), cp.Y(), cp.Z(), c->nStrips(), c->quadrant() ); + ftc.pos.SetXYZ( cp.X(), cp.Y(), cp.Z() ); + ftc.mQuadrant = c->quadrant(); + ftc.mRow = c->row(); + ftc.mOrientation = c->orientation(); + ftc.mPlane = c->plane(); + ftc.mId = c->id(); + ftc.mNStrips = c->nStrips(); + ftc.mX = c->x(); + ftc.mSumAdc = c->sumAdc(); + + mTreeData.fttClusters.add( ftc ); + } +} diff --git a/StRoot/StFwdTrackMaker/StFwdQAMaker.h b/StRoot/StFwdTrackMaker/StFwdQAMaker.h new file mode 100644 index 00000000000..8d6c7415801 --- /dev/null +++ b/StRoot/StFwdTrackMaker/StFwdQAMaker.h @@ -0,0 +1,396 @@ +#ifndef ST_FWD_TREE_MAKER_H +#define ST_FWD_TREE_MAKER_H + +#include "TClonesArray.h" +#ifndef __CINT__ +#include "GenFit/Track.h" +#include "StFwdTrackMaker/include/Tracker/FwdHit.h" +#include "StMuDSTMaker/COMMON/StMuFwdTrack.h" +#endif + +#include "StChain/StMaker.h" +#include "TTree.h" +#include "TVector3.h" +#include "TLorentzVector.h" +#include "StEvent/StEnumerations.h" + +class StMuFwdTrack; +class StMuFwdTrackProjection; +class ForwardTracker; +class StFwdTrack; + +/** @brief TClonesArray writer + * Helper class for writing TClonesArrays to TTree of custom class type + */ +template +class TClonesArrayWriter { + public: + TClonesArrayWriter() {} + ~TClonesArrayWriter() {} + + void createBranch( TTree *tree, const char* name, int buffSize = 256000, int splitLevel = 99){ + _tca = new TClonesArray( BranchType().classname() ); + tree->Branch( name, &this->_tca, buffSize, splitLevel ); + } + + void add( BranchType &branch ){ + if ( nullptr == this->_tca ) return; + BranchType *new_branch = new ((*this->_tca)[this->_n]) BranchType( ); + new_branch->copy( &branch ); + this->_n++; + } + + void add( BranchType *branch ){ + if ( nullptr == this->_tca || nullptr == branch) return; + BranchType *new_branch = new ((*this->_tca)[this->_n]) BranchType( ); + new_branch->copy( branch ); + this->_n++; + } + + void reset(){ + this->_n = 0; + if( nullptr != this->_tca ) + this->_tca->Clear(); + } + + UInt_t N() const { return _n; } + BranchType *at( UInt_t i ){ + if ( nullptr == _tca ) + return nullptr; + return (BranchType*)_tca->At( i ); + } + + protected: + TClonesArray * _tca = nullptr; + UInt_t _n = 0; +}; + +class FwdTreeHeader : public TObject { + public: + FwdTreeHeader() : TObject() { + run = 0; + event = 0; + tofmult = 0; + vpdVz = -999; + pv.SetXYZ(0, 0, 0); + } + + void set( int r, int e, int t, TVector3 &p ){ + run = r; + event = e; + tofmult = t; + pv = p; + } + + void clear() { + run = 0; + event = 0; + tofmult = 0; + TVector3 pv(-999, -999, -999); + vpdVz = -999; + } + + TVector3 pv; + int run, event, tofmult; + float vpdVz; + + ClassDef(FwdTreeHeader, 1) +}; + +class FwdTreeFttCluster : public TObject { + public: + TString classname() { return "FwdTreeFttCluster"; } + void clear(){ + TObject::Clear(); + pos.SetXYZ(-1, -1, -1); + mPlane = 0; + mQuadrant = 0; + mRow = 0; + mOrientation = kFttUnknownOrientation; + mNStrips = 0; + mSumAdc = 0.0; + mX = 0.0; + mSigma = 0.0; + } + + void copy( FwdTreeFttCluster *cluster ){ + pos = cluster->pos; + mPlane = cluster->mPlane; + mQuadrant = cluster->mQuadrant; + mRow = cluster->mRow; + mOrientation = cluster->mOrientation; + mNStrips = cluster->mNStrips; + mSumAdc = cluster->mSumAdc; + mX = cluster->mX; + mSigma = cluster->mSigma; + } + + TVector3 pos; + UChar_t mId = 0; + UChar_t mPlane = 0; + UChar_t mQuadrant = 0; + UChar_t mRow = 0; + UChar_t mOrientation = kFttUnknownOrientation; // Orientation of cluster + Int_t mNStrips=0; // Number of strips + Float_t mSumAdc=0.0; // Total ADC (0th moment) + Float_t mX=0.0; // Mean x ("center of gravity") in local grid coordinate (1st moment) + Float_t mSigma=0.0; + + ClassDef(FwdTreeFttCluster, 1); +}; + +class FwdTreeHit : public TObject { + public: + TString classname() { return "FwdTreeHit"; } + FwdTreeHit() : TObject() { + pos.SetXYZ(-1, -1, -1); + id = 0; + vol = 0; + det = 0; + trackId = -1; + } + FwdTreeHit(float x, float y, float z, int v, int d, int trkId = -1) : TObject() { + pos.SetXYZ(x, y, z); + id = 0; + vol = v; + det = d; + trackId = trkId; + } + + void set(float x, float y, float z, int v, int d, int trkId = -1) { + pos.SetXYZ(x, y, z); + id = 0; + vol = v; + det = d; + trackId = trkId; + } + + int id, vol, det, trackId; + TVector3 pos; + + void copy( FwdTreeHit *hit ){ + id = hit->id; + vol = hit->vol; + det = hit->det; + pos = hit->pos; + trackId = hit->trackId; + } + + ClassDef(FwdTreeHit, 2) +}; + +class FwdTreeTrackProjection : public TObject { + public: + FwdTreeTrackProjection() {} + FwdTreeTrackProjection( unsigned short detId, + TVector3 xyz, + TVector3 mom, + float c[9] ) { + set( detId, xyz, mom, c ); + } + + void set( unsigned short detId, + TVector3 xyz, + TVector3 mom, + float c[9]) { + mDetId = detId; + mXYZ = xyz; + mMom = mom; + memcpy( mCov, c, sizeof(mCov) ); + } + void copy( FwdTreeTrackProjection *other ){ + mDetId = other->mDetId; + mXYZ = other->mXYZ; + mMom = other->mMom; + memcpy( mCov, other->mCov, sizeof(mCov) ); + } + TVector3 mXYZ; + TVector3 mMom; + unsigned char mDetId; + float mCov[9]; + + float dx(){ + return sqrt( mCov[0] ); + } + float dy(){ + return sqrt( mCov[4] ); + } + float dz(){ + return sqrt( mCov[8] ); + } + + ClassDef(FwdTreeTrackProjection, 1) +}; + +class FwdTreeRecoTrack : public TObject { + public: + TString classname() { return "FwdTreeRecoTrack"; } + FwdTreeRecoTrack() : TObject() { + id = 0; + q = 0; + status = 0; + mom.SetXYZ(0, 0, 0); + projs.clear(); + } + + virtual void Clear(){ + TObject::Clear(); + seeds.clear(); + projs.clear(); + } + + void set( StMuFwdTrack *muFwdTrack ); + void set( StFwdTrack *fwdTrack ); + + /** + * Copies the values of the given FwdTreeRecoTrack object to the current object. + * + * @param hit The FwdTreeRecoTrack object to copy from. + */ + void copy( FwdTreeRecoTrack *hit ){ + id = hit->id; + q = hit->q; + status = hit->status; + mom = hit->mom; + seeds = hit->seeds; + mChi2 = hit->mChi2; + projs = hit->projs; + nFailedPoints = hit->nFailedPoints; + } + + int id, q, status; + int nFailedPoints; + float mChi2; + TVector3 mom; + vector projs; + vector seeds; + + ClassDef(FwdTreeRecoTrack, 4); +}; + +class StFcsDb; +class StFcsCluster; +class StMuFcsCluster; +class FwdTreeFcsCluster : public TObject { + public: + TString classname() { return "FwdTreeFcsCluster"; } + void copy( FwdTreeFcsCluster * clu ){ + mId = clu->mId; + mDetectorId = clu->mDetectorId; + mCategory = clu->mCategory; + mNTowers = clu->mNTowers; + mEnergy = clu->mEnergy; + mX = clu->mX; + mY = clu->mY; + mSigmaMin = clu->mSigmaMin; + mSigmaMax = clu->mSigmaMax; + mTheta = clu->mTheta; + mChi2Ndf1Photon = clu->mChi2Ndf1Photon; + mChi2Ndf2Photon = clu->mChi2Ndf2Photon; + mFourMomentum = clu->mFourMomentum; + pos = clu->pos; + } + + void set( StFcsCluster *clu, StFcsDb* fcsDb ); + void set( StMuFcsCluster *clu, StFcsDb* fcsDb ); + + Int_t mId=-1; // Eventwise cluster ID + UShort_t mDetectorId=0; // Detector starts from 1 + Int_t mCategory=0; // Category of cluster (see StFcsClusterCategory) + Int_t mNTowers=0; // Number of non-zero-energy tower hits in the cluster + Float_t mEnergy=0.0; // Total energy contained in this cluster (0th moment) + Float_t mX=0.0; // Mean x ("center of gravity") in local grid coordinate (1st moment) + Float_t mY=0.0; // Mean y ("center of gravity") in local grid coordinate (1st moment) + Float_t mSigmaMin=0.0; // Minimum 2nd moment + Float_t mSigmaMax=0.0; // Maximum 2nd moment (along major axis) + Float_t mTheta=0.0; //Angle in x-y plane that defines the direction of least-2nd-sigma + Float_t mChi2Ndf1Photon=0.0; // χ2 / ndf for 1-photon fit + Float_t mChi2Ndf2Photon=0.0; // χ2 / ndf for 2-photon fit + TLorentzVector mFourMomentum; // Cluster four momentum + TVector3 pos; // STAR XYZ position + + ClassDef(FwdTreeFcsCluster, 1); +}; + +class FwdTreeMonteCarloTrack : public TObject { + public: + + FwdTreeMonteCarloTrack() : TObject() { + id = 0; + q = 0; + status = 0; + mom.SetXYZ(0, 0, 0); + } + + int id, q, status; + TVector3 mom; + + ClassDef(FwdTreeMonteCarloTrack, 1); +}; + +/** @brief +* This class is a container for the data that will be written to the output tree. +*/ +struct FwdTreeData { + + /** @brief Primary event vertex*/ + FwdTreeHeader header; + TClonesArrayWriter ftt; + // TClonesArrayWriter fttClusters; + TClonesArrayWriter fttClusters; + TClonesArrayWriter fst; + + TClonesArrayWriter wcal; + TClonesArrayWriter hcal; + + TClonesArrayWriter reco; + + int nSeedTracks; + TClonesArrayWriter seeds; + + + void clear(); +}; + + +class StMuDstMaker; +class StMuDst; +class StMuFwdTrackCollection; +class StMuFcsCollection; +class StFwdTrackMaker; +class StEvent; + +class StFwdQAMaker : public StMaker { + + ClassDef(StFwdQAMaker, 0); + + public: + StFwdQAMaker(); + ~StFwdQAMaker(){/* nada */}; + + int Init(); + int Finish(); + int Make(); + void Clear(const Option_t *opts = ""); + + void FillFttClusters(); + void FillFcsStEvent(); + void FillFcsStMuDst(); + + protected: + TFile *mTreeFile = nullptr; + TTree *mTree = nullptr; + FwdTreeData mTreeData; + + StEvent *mStEvent = nullptr; + StMuDstMaker *mMuDstMaker = nullptr; + StMuDst *mMuDst = nullptr; + StMuFwdTrackCollection * mMuForwardTrackCollection = nullptr; + StMuFcsCollection *mMuFcsCollection = nullptr; + StFwdTrackMaker *mFwdTrackMaker = nullptr; + StFcsDb *mFcsDb = nullptr; + +}; + + +#endif From 660f552ea8a6369c0a789be77e45873292f9f44b Mon Sep 17 00:00:00 2001 From: Daniel Brandenburg Date: Mon, 29 Jul 2024 14:30:34 -0400 Subject: [PATCH 2/3] Use MuDst IO containers everywhere possible. --- StRoot/StFwdTrackMaker/StFwdQAMaker.cxx | 493 +++++++++--------------- StRoot/StFwdTrackMaker/StFwdQAMaker.h | 289 +++----------- 2 files changed, 243 insertions(+), 539 deletions(-) diff --git a/StRoot/StFwdTrackMaker/StFwdQAMaker.cxx b/StRoot/StFwdTrackMaker/StFwdQAMaker.cxx index b4e13b6f5e9..b7ce1cd225e 100644 --- a/StRoot/StFwdTrackMaker/StFwdQAMaker.cxx +++ b/StRoot/StFwdTrackMaker/StFwdQAMaker.cxx @@ -22,6 +22,7 @@ #include "StMuDSTMaker/COMMON/StMuDstMaker.h" #include "StMuDSTMaker/COMMON/StMuDst.h" +#include "StMuDSTMaker/COMMON/StMuEvent.h" #include "StMuDSTMaker/COMMON/StMuFstCollection.h" #include "StMuDSTMaker/COMMON/StMuFstHit.h" #include "StMuDSTMaker/COMMON/StMuPrimaryVertex.h" @@ -30,140 +31,67 @@ #include "StMuDSTMaker/COMMON/StMuFcsCollection.h" #include "StMuDSTMaker/COMMON/StMuFcsCluster.h" #include "StMuDSTMaker/COMMON/StMuFcsHit.h" +#include "StMuDSTMaker/COMMON/StMuFttCluster.h" +#include "StMuDSTMaker/COMMON/StMuFttPoint.h" +#include "StMuDSTMaker/COMMON/StMuMcTrack.h" +#include "StMuDSTMaker/COMMON/StMuFstHit.h" +// ClassImp(FcsClusterWithStarXYZ); -/** Clear the FwdTreeData from one event to next */ -void FwdTreeData::clear(){ +/** Clear the FwdQATreeData from one event to next */ +void FwdQATreeData::clear(){ header.clear(); - ftt.reset(); + mcTracks.reset(); + fttPoints.reset(); fttClusters.reset(); - fst.reset(); + fstPoints.reset(); reco.reset(); seeds.reset(); wcal.reset(); hcal.reset(); + wcalHits.reset(); + hcalHits.reset(); + epdHits.reset(); } -/** @brief Sets the Tree cluster from StFcsCluster -* @param clu - StFcsCluster -* @param fcsDb - StFcsDb -* @return void -*/ -void FwdTreeFcsCluster::set( StFcsCluster *clu, StFcsDb * fcsDb ) { - if (!clu) { - mId = -1; - return; - } - if ( fcsDb ){ - StThreeVectorD xyz = fcsDb->getStarXYZfromColumnRow(clu->detectorId(),clu->x(),clu->y()); - pos.SetXYZ( xyz.x(), xyz.y(), xyz.z() ); - } else { - pos.SetXYZ( 0, 0, 0 ); - } - - mDetectorId = clu->detectorId(); - - mEnergy = clu->energy(); - mNTowers = clu->nTowers(); - mId = clu->id(); - mCategory = clu->category(); - mX = clu->x(); - mY = clu->y(); - mSigmaMin = clu->sigmaMin(); - mSigmaMax = clu->sigmaMax(); - mTheta = clu->theta(); - mChi2Ndf1Photon = clu->chi2Ndf1Photon(); - mChi2Ndf2Photon = clu->chi2Ndf2Photon(); -} - -void FwdTreeFcsCluster::set( StMuFcsCluster *clu, StFcsDb * fcsDb ) { - if (!clu) { - mId = -1; - return; +FcsClusterWithStarXYZ::FcsClusterWithStarXYZ( StMuFcsCluster *clu, StFcsDb *fcsDb ) { + if ( nullptr == clu ) return; + StThreeVectorD xyz = fcsDb->getStarXYZfromColumnRow(clu->detectorId(),clu->x(),clu->y()); + float detOffset = 0.0; + if ( clu->detectorId() == kFcsEcalNorthDetId || clu->detectorId() == kFcsEcalSouthDetId ){ + detOffset = 715.0; // cm from IP + } else if ( clu->detectorId() == kFcsHcalNorthDetId || clu->detectorId() == kFcsHcalSouthDetId ){ + detOffset = 807.0; // cm from IP } - if ( fcsDb ){ - StThreeVectorD xyz = fcsDb->getStarXYZfromColumnRow(clu->detectorId(),clu->x(),clu->y()); - pos.SetXYZ( xyz.x(), xyz.y(), xyz.z() ); - } else { - pos.SetXYZ( 0, 0, 0 ); - } - - mDetectorId = clu->detectorId(); - mEnergy = clu->energy(); - mNTowers = clu->nTowers(); - mId = clu->id(); - mCategory = clu->category(); - mX = clu->x(); - mY = clu->y(); - mSigmaMin = clu->sigmaMin(); - mSigmaMax = clu->sigmaMax(); - mTheta = clu->theta(); - mChi2Ndf1Photon = clu->chi2Ndf1Photon(); - mChi2Ndf2Photon = clu->chi2Ndf2Photon(); + mXYZ.SetXYZ( xyz.x(), xyz.y(), xyz.z() + detOffset ); + mClu = clu; } -void FwdTreeRecoTrack::set( StMuFwdTrack *muFwd ){ - if ( !muFwd ){ - LOG_DEBUG << "Invalid muFwdTrack found, skipping FwdTreeRecoTrack::set" << endm; - return; - } - - if (muFwd->didFitConverge()) status = 2; - if (muFwd->didFitConvergeFully()) status = 1; - nFailedPoints = muFwd->numberOfFailedPoints(); - q = muFwd->charge(); - mom = muFwd->momentum(); - projs.clear(); - seeds.clear(); - for ( auto p : muFwd->mProjections ){ - FwdTreeTrackProjection tp; - tp.set( p.mDetId, p.mXYZ, p.mMom, p.mCov ); - projs.push_back( tp ); - } - for ( auto s : muFwd->mFSTPoints ){ - FwdTreeHit h; - h.set( s.mXYZ.X(), s.mXYZ.Y(), s.mXYZ.Z(), 1, 1 ); - seeds.push_back( h ); - } - for ( auto s : muFwd->mFTTPoints ){ - FwdTreeHit h; - h.set( s.mXYZ.X(), s.mXYZ.Y(), s.mXYZ.Z(), 1, 1 ); - seeds.push_back( h ); - } -} - -void FwdTreeRecoTrack::set( StFwdTrack *fwd ){ - if ( !fwd ){ - LOG_DEBUG << "Invalid fwdTrack found, skipping FwdTreeRecoTrack::set" << endm; - return; - } - - if (fwd->didFitConverge()) status = 2; - if (fwd->didFitConvergeFully()) status = 1; - nFailedPoints = fwd->numberOfFailedPoints(); - mChi2 = fwd->chi2(); - q = fwd->charge(); - mom.SetXYZ( fwd->momentum().x(), fwd->momentum().y(), fwd->momentum().z() ); - projs.clear(); - seeds.clear(); - for ( auto p : fwd->mProjections ){ - FwdTreeTrackProjection tp; - tp.set( p.mDetId, TVector3(p.mXYZ.x(), p.mXYZ.y(), p.mXYZ.z()), TVector3(p.mMom.x(), p.mMom.y(), p.mMom.z()), p.mCov ); - projs.push_back( tp ); - } - for ( auto s : fwd->mFSTPoints ){ - FwdTreeHit h; - h.set( s.mXYZ.x(), s.mXYZ.y(), s.mXYZ.z(), 1, 1 ); - seeds.push_back( h ); - } - for ( auto s : fwd->mFTTPoints ){ - FwdTreeHit h; - h.set( s.mXYZ.x(), s.mXYZ.y(), s.mXYZ.z(), 1, 1 ); - seeds.push_back( h ); +FcsHitWithStarXYZ::FcsHitWithStarXYZ( StMuFcsHit *hit, StFcsDb *fcsDb ) { + if ( nullptr == hit ) return; + StThreeVectorD xyz = fcsDb->getStarXYZ(hit->detectorId(),hit->id()); + float detOffset = 0.0; + if ( hit->detectorId() == kFcsEcalNorthDetId || hit->detectorId() == kFcsEcalSouthDetId ){ + detOffset = 715.0; // cm from IP + } else if ( hit->detectorId() == kFcsHcalNorthDetId || hit->detectorId() == kFcsHcalSouthDetId ){ + detOffset = 807.0; // cm from IP + } else if ( hit->detectorId() == kFcsPresNorthDetId || hit->detectorId() == kFcsPresSouthDetId ){ + StEpdGeom epdgeo; + double zepd=375.0; + int pp,tt,n; + double x[5],y[5]; + fcsDb->getEPDfromId(hit->detectorId(),hit->id(),pp,tt); + epdgeo.GetCorners(100*pp+tt,&n,x,y); + double x0 = (x[0] + x[1] + x[2] + x[3]) / 4.0; + double y0 = (y[0] + y[1] + y[2] + y[3]) / 4.0; + xyz.setX(x0); + xyz.setY(y0); + xyz.setZ(zepd); } + mXYZ.SetXYZ( xyz.x(), xyz.y(), xyz.z() + detOffset ); + mHit = hit; } - StFwdQAMaker::StFwdQAMaker() : StMaker("fwdQAMaker"), mTreeFile(nullptr), mTree(nullptr) { } @@ -174,13 +102,19 @@ int StFwdQAMaker::Init() { mTree = new TTree("fwd", "fwd tracking tree"); mTree->Branch("header", &mTreeData. header, 3200, 99 ); + mTreeData.mcTracks.createBranch(mTree, "mcTracks"); mTree->Branch("nSeedTracks", &mTreeData.nSeedTracks, "nSeedTracks/I"); - mTreeData.fst.createBranch(mTree, "fst"); - mTreeData.ftt.createBranch(mTree, "ftt"); + mTreeData.fstPoints.createBranch(mTree, "fstHits"); + mTreeData.fttPoints.createBranch(mTree, "fttPoints"); mTreeData.fttClusters.createBranch(mTree, "fttClusters"); + mTreeData.fstPoints.createBranch(mTree, "fstPoints"); mTreeData.wcal.createBranch(mTree, "wcalClusters"); mTreeData.hcal.createBranch(mTree, "hcalClusters"); + mTreeData.wcalHits.createBranch(mTree, "wcalHits"); + mTreeData.hcalHits.createBranch(mTree, "hcalHits"); + mTreeData.epdHits.createBranch(mTree, "epdHits"); + mTreeData.reco.createBranch(mTree, "reco"); mTreeData.seeds.createBranch(mTree, "seeds"); return kStOk; @@ -196,7 +130,15 @@ int StFwdQAMaker::Finish() { return kStOk; } int StFwdQAMaker::Make() { - LOG_DEBUG << "SETUP START" << endm; + LOG_INFO << "FWD Report:" << endm; + StEvent *mStEvent = static_cast(GetInputDS("StEvent")); + if ( mStEvent ){ + // report number of fwd tracks + auto fwdTracks = mStEvent->fwdTrackCollection(); + LOG_INFO << "Number of FwdTracks (StFwdTrackCollection): " << fwdTracks->tracks().size() << endm; + LOG_INFO << "Number of Ftt Points (StEvent)" << mStEvent->fttCollection()->points().size() << endm; + } + LOG_INFO << "SETUP START" << endm; // setup the datasets / makers mMuDstMaker = (StMuDstMaker *)GetMaker("MuDst"); if(mMuDstMaker) { @@ -217,108 +159,93 @@ int StFwdQAMaker::Make() { // return kStOk; } - // Event header info from stevent - StEvent *mStEvent = static_cast(GetInputDS("StEvent")); + mTreeData.header.run = mMuDst->event()->runNumber(); LOG_DEBUG << "SETUP COMPLETE" << endm; - // Get the primary vertex used by the FWD Tracker - auto eventPV = mFwdTrackMaker->GetEventPrimaryVertex(); - mTreeData.header.set( 0, 0, 0, eventPV ); - - // Fill Header from StEvent - if ( mMuDstMaker ){ + auto muFstCollection = mMuDst->muFstCollection(); + if ( muFstCollection ){ + LOG_DEBUG << "MuDst has #fst hits: " << muFstCollection->numberOfHits() << endm; + for ( size_t i = 0; i < muFstCollection->numberOfHits(); i++ ){ + StMuFstHit * h = muFstCollection->getHit(i); + mTreeData.fstPoints.add( h ); + } } - else if ( mStEvent ){ - StBTofCollection *btofC = mStEvent->btofCollection(); - if (btofC) { - StBTofHeader * btofHeader = btofC->tofHeader(); - if (btofHeader){ - - int nEast = btofHeader->numberOfVpdHits( east ); - int nWest = btofHeader->numberOfVpdHits( west ); - int nTof = btofC->tofHits().size(); - float vpdVz = btofHeader->vpdVz(); - - mTreeData.header.set( 1, GetEventNumber(), nTof, eventPV ); - mTreeData.header.vpdVz = vpdVz; - } - } // btofC != nullptr + FillMcTracks(); + FillTracks(); + FillFstPoints(); + FillFttClusters(); + FillFcsStMuDst(); + mTree->Fill(); + return kStOk; +} +void StFwdQAMaker::Clear(const Option_t *opts) { + mTreeData.clear(); + return; +} + +void StFwdQAMaker::FillFstPoints(){ + StMuFstCollection * fst = mMuDst->muFstCollection(); + if (!fst) { + LOG_WARN << "No StMuFstCollection ... bye-bye" << endm; + return; } - LOG_DEBUG << "HEADER COMPLETE" << endm; - FwdTreeHit fh; + // size_t numFwdHitsPrior = mFwdHitsFst.size(); + LOG_INFO << "Loading " << fst->numberOfHits() << " StMuFstHits" << endm; + // TMatrixDSym hitCov3(3); + for ( unsigned int index = 0; index < fst->numberOfHits(); index++){ + StMuFstHit * muFstHit = fst->getHit( index ); + mTreeData.fstPoints.add( muFstHit ); - auto fttHits = mFwdTrackMaker -> GetFttHits(); - for ( auto h : fttHits ){ - fh.set( h.getX(), h.getY(), h.getZ(), 1, 2 ); - mTreeData.ftt.add( fh ); - } - LOG_DEBUG << "FTT COMPLETE" << endm; - auto fstHits = mFwdTrackMaker -> GetFstHits(); - for ( auto h : fstHits ){ - fh.set( h.getX(), h.getY(), h.getZ(), 1, 1 ); - mTreeData.fst.add( fh ); - } - LOG_DEBUG << "FST COMPLETE" << endm; - auto seeds = mFwdTrackMaker -> getTrackSeeds(); - LOG_DEBUG << "Number of Seeds to save: " << seeds.size() << endm; - if (seeds.size() > 5000){ - LOG_WARN << "Number of seeds is too large, truncating to 5000" << endm; - // seeds.resize(1000); - } - size_t iSeed = 0; - for ( auto s : seeds ){ - for ( auto h : s ){ - fh.set( h->getX(), h->getY(), h->getZ(), 1, 1 ); - fh.trackId = iSeed; - mTreeData.seeds.add( fh ); - } - iSeed++; - if (iSeed > 5000) break; - } // for s in seeds - LOG_DEBUG << "SEEDS COMPLETE, saved " << iSeed << " seed hits" << endm; - mTreeData.nSeedTracks = (int)iSeed; - - size_t iTrack = 0; - FwdTreeRecoTrack frt; - FwdTreeHit fth; + + // float vR = muFstHit->localPosition(0); + // float vPhi = muFstHit->localPosition(1); + // float vZ = muFstHit->localPosition(2); + + // const float dz0 = fabs( vZ - mFstZFromGeom[0] ); + // const float dz1 = fabs( vZ - mFstZFromGeom[1] ); + // const float dz2 = fabs( vZ - mFstZFromGeom[2] ); + // static const float fstThickness = 2.0; // thickness in cm between inner and outer on sigle plane + + // // assign disk according to which z value the hit has, within the z-plane thickness + // int d = 0 * ( dz0 < fstThickness ) + 1 * ( dz1 < fstThickness ) + 2 * ( dz2 < fstThickness ); + + // float x0 = vR * cos( vPhi ); + // float y0 = vR * sin( vPhi ); + // hitCov3 = makeSiCovMat( TVector3( x0, y0, vZ ), mFwdConfig ); + + // LOG_DEBUG << "FST HIT: d = " << d << ", x=" << x0 << ", y=" << y0 << ", z=" << vZ << endm; + // mFstHits.push_back( TVector3( x0, y0, vZ) ); + + // // we use d+4 so that both FTT and FST start at 4 + // mFwdHitsFst.push_back(FwdHit(count++, x0, y0, vZ, d+4, 0, hitCov3, nullptr)); + // count++; + } // index +} + +void StFwdQAMaker::FillTracks() { + mTreeData.nSeedTracks = 0; if ( mMuForwardTrackCollection ){ LOG_DEBUG << "Adding " << mMuForwardTrackCollection->numberOfFwdTracks() << " FwdTracks (MuDst)" << endm; for ( size_t iTrack = 0; iTrack < mMuForwardTrackCollection->numberOfFwdTracks(); iTrack++ ){ auto muTrack = mMuForwardTrackCollection->getFwdTrack(iTrack); - frt.set( muTrack ); - mTreeData.reco.add( frt ); - if ( iTrack > 5000 ) break; - } - } + mTreeData.reco.add( muTrack ); - if ( mStEvent && mStEvent->fwdTrackCollection() ){ - iTrack = 0; - auto fwdTracks = mStEvent->fwdTrackCollection()->tracks(); - LOG_DEBUG << "Adding " << fwdTracks.size() << " FwdTracks (StEvent)" << endm; - for ( auto ft : fwdTracks ){ - frt.set( ft ); - mTreeData.reco.add( frt ); - if ( iTrack > 5000 ) break; - iTrack++; + for (auto fsth : muTrack->mFSTPoints){ + mTreeData.seeds.add( fsth ); + mTreeData.nSeedTracks++; + } + for (auto ftth : muTrack->mFTTPoints){ + mTreeData.seeds.add( ftth ); + mTreeData.nSeedTracks++; + } + if ( iTrack > 5000 ) { + LOG_WARN << "Truncating to 5000 tracks" << endm; + break; + } } } LOG_DEBUG << "TRACKS COMPLETE" << endm; - - - FillFttClusters(); - if ( mMuDstMaker ){ - FillFcsStMuDst(); - } else { - FillFcsStEvent(); - } - - mTree->Fill(); - return kStOk; -} -void StFwdQAMaker::Clear(const Option_t *opts) { - mTreeData.clear(); - return; } void StFwdQAMaker::FillFcsStMuDst( ) { @@ -333,118 +260,74 @@ void StFwdQAMaker::FillFcsStMuDst( ) { return; } - StEpdGeom epdgeo; - FwdTreeFcsCluster fcsclu; + StFcsDb* fcsDb=static_cast(GetDataSet("fcsDb")); + // LOAD ECAL / HCAL CLUSTERS + LOG_INFO << "MuDst has #fcs clusters: " << fcs->numberOfClusters() << endm; for( size_t i = 0; i < fcs->numberOfClusters(); i++){ StMuFcsCluster * clu = fcs->getCluster(i); - LOG_DEBUG << "FCS CLUSTERS: " << fcs->numberOfClusters() << endm; - fcsclu.set( clu, mFcsDb ); - + FcsClusterWithStarXYZ *cluSTAR = new FcsClusterWithStarXYZ(clu, fcsDb); if ( clu->detectorId() == kFcsEcalNorthDetId || clu->detectorId() == kFcsEcalSouthDetId ){ - LOG_DEBUG << "Adding WCAL Cluster to FwdTree" << endm; - mTreeData.wcal.add( fcsclu ); + LOG_INFO << "Adding WCAL Cluster to FwdTree" << endm; + mTreeData.wcal.add( cluSTAR ); } else if ( clu->detectorId() == kFcsHcalNorthDetId || clu->detectorId() == kFcsHcalSouthDetId ){ - LOG_DEBUG << "Adding HCAL Cluster to FwdTree" << endm; - mTreeData.hcal.add( fcsclu ); - } - } -} + LOG_INFO << "Adding HCAL Cluster to FwdTree" << endm; + mTreeData.hcal.add( cluSTAR ); + } -void StFwdQAMaker::FillFcsStEvent( ) { - if ( !mStEvent ){ - LOG_DEBUG << "No StEvent found, skipping StFwdQAMaker::FillFcsStEvent" << endm; - return; + delete cluSTAR; } - StFcsCollection* fcsCol = mStEvent->fcsCollection(); - if ( !fcsCol ){ - LOG_DEBUG << "No StFcsCollection found, skipping StFwdQAMaker::FillFcsStEvent" << endm; - return; - } - - StEpdGeom epdgeo; - FwdTreeFcsCluster fcsclu; // LOAD ECAL / HCAL CLUSTERS - for ( int idet = 0; idet < 4; idet++ ){ - StSPtrVecFcsCluster& clusters = fcsCol->clusters(idet); - LOG_DEBUG << "FCS CLUSTERS [" << idet << "]: " << clusters.size() << endm; - int nc=fcsCol->numberOfClusters(idet); - for ( int i = 0; i < nc; i++ ){ - StFcsCluster* clu = clusters[i]; - - fcsclu.set( clu, mFcsDb ); - - if ( clu->detectorId() == kFcsEcalNorthDetId || clu->detectorId() == kFcsEcalSouthDetId ){ - LOG_DEBUG << "Adding WCAL Cluster to FwdTree" << endm; - mTreeData.wcal.add( fcsclu ); - } else if ( clu->detectorId() == kFcsHcalNorthDetId || clu->detectorId() == kFcsHcalSouthDetId ){ - LOG_DEBUG << "Adding HCAL Cluster to FwdTree" << endm; - mTreeData.hcal.add( fcsclu ); - } + LOG_INFO << "MuDst has #fcs hits: " << fcs->numberOfHits() << endm; + for( size_t i = 0; i < fcs->numberOfHits(); i++){ + StMuFcsHit * hit = fcs->getHit(i); + FcsHitWithStarXYZ *hitSTAR = new FcsHitWithStarXYZ(hit, fcsDb); + if ( hit->detectorId() == kFcsEcalNorthDetId || hit->detectorId() == kFcsEcalSouthDetId ){ + LOG_DEBUG << "Adding WCAL Cluster to FwdTree" << endm; + mTreeData.wcalHits.add( hitSTAR ); + } else if ( hit->detectorId() == kFcsHcalNorthDetId || hit->detectorId() == kFcsHcalSouthDetId ){ + LOG_DEBUG << "Adding HCAL Cluster to FwdTree" << endm; + mTreeData.hcalHits.add( hitSTAR ); + } else if ( hit->detectorId() == kFcsPresNorthDetId || hit->detectorId() == kFcsPresSouthDetId ){ + LOG_DEBUG << "Adding PRES hit to FwdTree" << endm; + mTreeData.epdHits.add( hitSTAR ); + } + delete hitSTAR; + } - } // i - } // idet + // TODO, cleanup? +} +void StFwdQAMaker::FillMcTracks(){ + // Retrieve pointer to MC tracks + TClonesArray *mcTracks = mMuDst->mcArray(1); + LOG_INFO << "MuDst has #mc tracks: " << mcTracks->GetEntriesFast() << endm; + // Loop over MC vertices + for (Int_t iVtx=0; iVtxGetEntriesFast(); iVtx++) { + // Retrieve i-th MC vertex from MuDst + StMuMcTrack *mcTrack = (StMuMcTrack*)mcTracks->UncheckedAt(iVtx); + if ( !mcTrack ) continue; + + // Add MC track to the tree + mTreeData.mcTracks.add( mcTrack ); + } } + void StFwdQAMaker::FillFttClusters(){ - float pz[] = {280.90499, 303.70498, 326.60501, 349.40499}; - float SCALE = 1.0; - TVector3 cp; - FwdTreeHit fh; - FwdTreeFttCluster ftc; - StEvent *event = static_cast(GetInputDS("StEvent")); - if ( !event || !event->fttCollection() ) return; - - LOG_DEBUG << "FTT RawHits: " << event->fttCollection()->numberOfRawHits() << endm; - LOG_DEBUG << "FTT Clusters: " << event->fttCollection()->numberOfClusters() << endm; - - for ( size_t i = 0; i < event->fttCollection()->numberOfClusters(); i++ ){ - StFttCluster* c = event->fttCollection()->clusters()[i]; - if ( c->nStrips() < 1 ) continue; - float dw = 0.05, dlh = 60.0, dlv = 60.0; - float mx = 0.0, my = 0.0; - float sx = 1.0, sy = 1.0; - - - if ( c->quadrant() == kFttQuadrantA ){ - mx = 0; my = 0; - sx = 1.0; sy = 1.0; - } else if ( c->quadrant() == kFttQuadrantB ){ - mx = 10.16*SCALE; my = 0.0*SCALE; - sy = -1; - dlv = -dlv; - - } else if ( c->quadrant() == kFttQuadrantC ){ - mx = -10.16*SCALE ; my = -00.0*SCALE; - sx = -1.0; sy = -1.0; - dlh = -dlh; dlv = -dlv; - - } else if ( c->quadrant() == kFttQuadrantD ){ - sx = -1; - dlh = -dlh; + + auto muFttCollection = mMuDst->muFttCollection(); + if ( muFttCollection ){ + LOG_DEBUG << "MuDst has #ftt clusters: " << muFttCollection->numberOfClusters() << endm; + for ( size_t i = 0; i < muFttCollection->numberOfClusters(); i++ ){ + StMuFttCluster * c = muFttCollection->getCluster(i); + mTreeData.fttClusters.add( c ); } - cp.SetZ( -pz[ c->plane() ] * SCALE ); - if ( c->orientation() == kFttHorizontal ){ - cp.SetY( my + sy * c->x()/10.0 * SCALE ); - cp.SetX( mx ); - } else if ( c->orientation() == kFttVertical ){ - cp.SetX( mx + sx * c->x()/10.0 * SCALE ); - cp.SetY( my ); + for ( size_t i = 0; i < muFttCollection->numberOfPoints(); i++ ){ + StMuFttPoint * c = muFttCollection->getPoint(i); + mTreeData.fttPoints.add( c ); } - // fh.set( cp.X(), cp.Y(), cp.Z(), c->nStrips(), c->quadrant() ); - ftc.pos.SetXYZ( cp.X(), cp.Y(), cp.Z() ); - ftc.mQuadrant = c->quadrant(); - ftc.mRow = c->row(); - ftc.mOrientation = c->orientation(); - ftc.mPlane = c->plane(); - ftc.mId = c->id(); - ftc.mNStrips = c->nStrips(); - ftc.mX = c->x(); - ftc.mSumAdc = c->sumAdc(); - - mTreeData.fttClusters.add( ftc ); } } diff --git a/StRoot/StFwdTrackMaker/StFwdQAMaker.h b/StRoot/StFwdTrackMaker/StFwdQAMaker.h index 8d6c7415801..bc60f98d5b4 100644 --- a/StRoot/StFwdTrackMaker/StFwdQAMaker.h +++ b/StRoot/StFwdTrackMaker/StFwdQAMaker.h @@ -13,11 +13,13 @@ #include "TVector3.h" #include "TLorentzVector.h" #include "StEvent/StEnumerations.h" +#include "StThreeVectorD.hh" class StMuFwdTrack; class StMuFwdTrackProjection; class ForwardTracker; class StFwdTrack; +class StMuMcTrack; /** @brief TClonesArray writer * Helper class for writing TClonesArrays to TTree of custom class type @@ -29,21 +31,21 @@ class TClonesArrayWriter { ~TClonesArrayWriter() {} void createBranch( TTree *tree, const char* name, int buffSize = 256000, int splitLevel = 99){ - _tca = new TClonesArray( BranchType().classname() ); + _tca = new TClonesArray( BranchType::Class_Name() ); tree->Branch( name, &this->_tca, buffSize, splitLevel ); } void add( BranchType &branch ){ if ( nullptr == this->_tca ) return; BranchType *new_branch = new ((*this->_tca)[this->_n]) BranchType( ); - new_branch->copy( &branch ); + *new_branch = branch; this->_n++; } void add( BranchType *branch ){ if ( nullptr == this->_tca || nullptr == branch) return; BranchType *new_branch = new ((*this->_tca)[this->_n]) BranchType( ); - new_branch->copy( branch ); + *new_branch = *branch; this->_n++; } @@ -97,256 +99,72 @@ class FwdTreeHeader : public TObject { ClassDef(FwdTreeHeader, 1) }; -class FwdTreeFttCluster : public TObject { - public: - TString classname() { return "FwdTreeFttCluster"; } - void clear(){ - TObject::Clear(); - pos.SetXYZ(-1, -1, -1); - mPlane = 0; - mQuadrant = 0; - mRow = 0; - mOrientation = kFttUnknownOrientation; - mNStrips = 0; - mSumAdc = 0.0; - mX = 0.0; - mSigma = 0.0; - } - - void copy( FwdTreeFttCluster *cluster ){ - pos = cluster->pos; - mPlane = cluster->mPlane; - mQuadrant = cluster->mQuadrant; - mRow = cluster->mRow; - mOrientation = cluster->mOrientation; - mNStrips = cluster->mNStrips; - mSumAdc = cluster->mSumAdc; - mX = cluster->mX; - mSigma = cluster->mSigma; - } - - TVector3 pos; - UChar_t mId = 0; - UChar_t mPlane = 0; - UChar_t mQuadrant = 0; - UChar_t mRow = 0; - UChar_t mOrientation = kFttUnknownOrientation; // Orientation of cluster - Int_t mNStrips=0; // Number of strips - Float_t mSumAdc=0.0; // Total ADC (0th moment) - Float_t mX=0.0; // Mean x ("center of gravity") in local grid coordinate (1st moment) - Float_t mSigma=0.0; - - ClassDef(FwdTreeFttCluster, 1); -}; - -class FwdTreeHit : public TObject { - public: - TString classname() { return "FwdTreeHit"; } - FwdTreeHit() : TObject() { - pos.SetXYZ(-1, -1, -1); - id = 0; - vol = 0; - det = 0; - trackId = -1; - } - FwdTreeHit(float x, float y, float z, int v, int d, int trkId = -1) : TObject() { - pos.SetXYZ(x, y, z); - id = 0; - vol = v; - det = d; - trackId = trkId; - } - - void set(float x, float y, float z, int v, int d, int trkId = -1) { - pos.SetXYZ(x, y, z); - id = 0; - vol = v; - det = d; - trackId = trkId; - } - - int id, vol, det, trackId; - TVector3 pos; - - void copy( FwdTreeHit *hit ){ - id = hit->id; - vol = hit->vol; - det = hit->det; - pos = hit->pos; - trackId = hit->trackId; - } +class StFcsDb; +class StFcsCluster; +class StFcsHit; +class StMuFcsCluster; +class StMuFcsHit; +class StMuFttCluster; +class StMuFttPoint; +class StMuFstHit; +class StMuFwdTrackSeedPoint; - ClassDef(FwdTreeHit, 2) -}; -class FwdTreeTrackProjection : public TObject { +/** + * @brief Store Cluster with STAR XYZ position + * + */ +class FcsClusterWithStarXYZ: public TObject { public: - FwdTreeTrackProjection() {} - FwdTreeTrackProjection( unsigned short detId, - TVector3 xyz, - TVector3 mom, - float c[9] ) { - set( detId, xyz, mom, c ); - } - - void set( unsigned short detId, - TVector3 xyz, - TVector3 mom, - float c[9]) { - mDetId = detId; - mXYZ = xyz; - mMom = mom; - memcpy( mCov, c, sizeof(mCov) ); - } - void copy( FwdTreeTrackProjection *other ){ - mDetId = other->mDetId; - mXYZ = other->mXYZ; - mMom = other->mMom; - memcpy( mCov, other->mCov, sizeof(mCov) ); + FcsClusterWithStarXYZ() { + mXYZ.SetXYZ(0, 0, 0); + mClu = nullptr; } + FcsClusterWithStarXYZ( StMuFcsCluster *clu, StFcsDb *fcsDb ); TVector3 mXYZ; - TVector3 mMom; - unsigned char mDetId; - float mCov[9]; - - float dx(){ - return sqrt( mCov[0] ); - } - float dy(){ - return sqrt( mCov[4] ); - } - float dz(){ - return sqrt( mCov[8] ); - } - - ClassDef(FwdTreeTrackProjection, 1) + StMuFcsCluster *mClu; + ClassDef(FcsClusterWithStarXYZ, 1); }; -class FwdTreeRecoTrack : public TObject { - public: - TString classname() { return "FwdTreeRecoTrack"; } - FwdTreeRecoTrack() : TObject() { - id = 0; - q = 0; - status = 0; - mom.SetXYZ(0, 0, 0); - projs.clear(); - } - - virtual void Clear(){ - TObject::Clear(); - seeds.clear(); - projs.clear(); - } - - void set( StMuFwdTrack *muFwdTrack ); - void set( StFwdTrack *fwdTrack ); - - /** - * Copies the values of the given FwdTreeRecoTrack object to the current object. - * - * @param hit The FwdTreeRecoTrack object to copy from. - */ - void copy( FwdTreeRecoTrack *hit ){ - id = hit->id; - q = hit->q; - status = hit->status; - mom = hit->mom; - seeds = hit->seeds; - mChi2 = hit->mChi2; - projs = hit->projs; - nFailedPoints = hit->nFailedPoints; - } - - int id, q, status; - int nFailedPoints; - float mChi2; - TVector3 mom; - vector projs; - vector seeds; - - ClassDef(FwdTreeRecoTrack, 4); -}; - -class StFcsDb; -class StFcsCluster; -class StMuFcsCluster; -class FwdTreeFcsCluster : public TObject { +/** + * @brief Store Hit with STAR XYZ position + * + */ +class FcsHitWithStarXYZ: public TObject { public: - TString classname() { return "FwdTreeFcsCluster"; } - void copy( FwdTreeFcsCluster * clu ){ - mId = clu->mId; - mDetectorId = clu->mDetectorId; - mCategory = clu->mCategory; - mNTowers = clu->mNTowers; - mEnergy = clu->mEnergy; - mX = clu->mX; - mY = clu->mY; - mSigmaMin = clu->mSigmaMin; - mSigmaMax = clu->mSigmaMax; - mTheta = clu->mTheta; - mChi2Ndf1Photon = clu->mChi2Ndf1Photon; - mChi2Ndf2Photon = clu->mChi2Ndf2Photon; - mFourMomentum = clu->mFourMomentum; - pos = clu->pos; + FcsHitWithStarXYZ() { + mXYZ.SetXYZ(0, 0, 0); + mHit = nullptr; } - - void set( StFcsCluster *clu, StFcsDb* fcsDb ); - void set( StMuFcsCluster *clu, StFcsDb* fcsDb ); - - Int_t mId=-1; // Eventwise cluster ID - UShort_t mDetectorId=0; // Detector starts from 1 - Int_t mCategory=0; // Category of cluster (see StFcsClusterCategory) - Int_t mNTowers=0; // Number of non-zero-energy tower hits in the cluster - Float_t mEnergy=0.0; // Total energy contained in this cluster (0th moment) - Float_t mX=0.0; // Mean x ("center of gravity") in local grid coordinate (1st moment) - Float_t mY=0.0; // Mean y ("center of gravity") in local grid coordinate (1st moment) - Float_t mSigmaMin=0.0; // Minimum 2nd moment - Float_t mSigmaMax=0.0; // Maximum 2nd moment (along major axis) - Float_t mTheta=0.0; //Angle in x-y plane that defines the direction of least-2nd-sigma - Float_t mChi2Ndf1Photon=0.0; // χ2 / ndf for 1-photon fit - Float_t mChi2Ndf2Photon=0.0; // χ2 / ndf for 2-photon fit - TLorentzVector mFourMomentum; // Cluster four momentum - TVector3 pos; // STAR XYZ position - - ClassDef(FwdTreeFcsCluster, 1); + FcsHitWithStarXYZ( StMuFcsHit *hit, StFcsDb *fcsDb ); + TVector3 mXYZ; + StMuFcsHit *mHit; + ClassDef(FcsHitWithStarXYZ, 1); }; -class FwdTreeMonteCarloTrack : public TObject { - public: - - FwdTreeMonteCarloTrack() : TObject() { - id = 0; - q = 0; - status = 0; - mom.SetXYZ(0, 0, 0); - } - - int id, q, status; - TVector3 mom; - - ClassDef(FwdTreeMonteCarloTrack, 1); -}; /** @brief * This class is a container for the data that will be written to the output tree. */ -struct FwdTreeData { +struct FwdQATreeData { /** @brief Primary event vertex*/ FwdTreeHeader header; - TClonesArrayWriter ftt; - // TClonesArrayWriter fttClusters; - TClonesArrayWriter fttClusters; - TClonesArrayWriter fst; - - TClonesArrayWriter wcal; - TClonesArrayWriter hcal; - - TClonesArrayWriter reco; + /** @brief MC tracks */ + TClonesArrayWriter mcTracks; + TClonesArrayWriter fttPoints; + TClonesArrayWriter fttClusters; + TClonesArrayWriter fstPoints; + + TClonesArrayWriter wcal; + TClonesArrayWriter wcalHits; + TClonesArrayWriter hcal; + TClonesArrayWriter hcalHits; + TClonesArrayWriter epdHits; + TClonesArrayWriter reco; int nSeedTracks; - TClonesArrayWriter seeds; + TClonesArrayWriter seeds; void clear(); @@ -374,13 +192,16 @@ class StFwdQAMaker : public StMaker { void Clear(const Option_t *opts = ""); void FillFttClusters(); + void FillFstPoints(); void FillFcsStEvent(); void FillFcsStMuDst(); + void FillTracks(); + void FillMcTracks(); protected: TFile *mTreeFile = nullptr; TTree *mTree = nullptr; - FwdTreeData mTreeData; + FwdQATreeData mTreeData; StEvent *mStEvent = nullptr; StMuDstMaker *mMuDstMaker = nullptr; From d86c4d369bd77e264757f4babd0f4df0e75d9405 Mon Sep 17 00:00:00 2001 From: Daniel Brandenburg Date: Mon, 29 Jul 2024 14:30:56 -0400 Subject: [PATCH 3/3] fix global consts -> move to local variables to avoid multiple definitions --- StRoot/StFwdTrackMaker/include/Tracker/ObjExporter.h | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/StRoot/StFwdTrackMaker/include/Tracker/ObjExporter.h b/StRoot/StFwdTrackMaker/include/Tracker/ObjExporter.h index 08ffda49a8e..623b0172452 100644 --- a/StRoot/StFwdTrackMaker/include/Tracker/ObjExporter.h +++ b/StRoot/StFwdTrackMaker/include/Tracker/ObjExporter.h @@ -9,8 +9,6 @@ class StEvent; -float DEGS_TO_RAD = 3.14159f/180.0f; -float SCALE = 0.1; class ObjExporter { public: @@ -33,6 +31,7 @@ class ObjExporter { int p, s, i, j; float x, y, z, out; int nPitch = nLongitude + 1; + const float DEGS_TO_RAD = 3.14159f/180.0f; float pitchInc = (180. / (float)nPitch) * DEGS_TO_RAD; float rotInc = (360. / (float)nLatitude) * DEGS_TO_RAD; @@ -209,6 +208,7 @@ class ObjExporter { fout << "usemtl stgc_hits\n" << endl; float pz[] = {280.90499, 303.70498, 326.60501, 349.40499}; TVector3 cp; + const float SCALE = 0.1; for ( size_t i = 0; i < event->fttCollection()->numberOfClusters(); i++ ){ @@ -263,6 +263,7 @@ class ObjExporter { std::vector &fcsClusterEnergy ){ + const float SCALE = 0.1; LOG_INFO << "Writing: " << filename << endm; numVertices = 0; // OPEN output