diff --git a/StRoot/StFwdTrackMaker/StFwdQAMaker.cxx b/StRoot/StFwdTrackMaker/StFwdQAMaker.cxx new file mode 100644 index 00000000000..b7ce1cd225e --- /dev/null +++ b/StRoot/StFwdTrackMaker/StFwdQAMaker.cxx @@ -0,0 +1,333 @@ +#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/StMuEvent.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" +#include "StMuDSTMaker/COMMON/StMuFttCluster.h" +#include "StMuDSTMaker/COMMON/StMuFttPoint.h" +#include "StMuDSTMaker/COMMON/StMuMcTrack.h" +#include "StMuDSTMaker/COMMON/StMuFstHit.h" + +// ClassImp(FcsClusterWithStarXYZ); + +/** Clear the FwdQATreeData from one event to next */ +void FwdQATreeData::clear(){ + header.clear(); + mcTracks.reset(); + fttPoints.reset(); + fttClusters.reset(); + fstPoints.reset(); + reco.reset(); + seeds.reset(); + wcal.reset(); + hcal.reset(); + wcalHits.reset(); + hcalHits.reset(); + epdHits.reset(); +} + +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 + } + mXYZ.SetXYZ( xyz.x(), xyz.y(), xyz.z() + detOffset ); + mClu = clu; +} + +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) { + +} + +int StFwdQAMaker::Init() { + + mTreeFile = new TFile("fwdtree.root", "RECREATE"); + 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.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; +} +int StFwdQAMaker::Finish() { + + if ( mTreeFile && mTree ){ + mTreeFile->cd(); + mTree->Write(); + mTreeFile->Write(); + LOG_DEBUG << "StFwdQA File written" << endm; + } + return kStOk; +} +int StFwdQAMaker::Make() { + 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) { + 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; + } + + mTreeData.header.run = mMuDst->event()->runNumber(); + LOG_DEBUG << "SETUP COMPLETE" << endm; + + 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 ); + } + } + 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; + } + + // 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 ); + + + // 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); + mTreeData.reco.add( muTrack ); + + 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; +} + +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; + } + + 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); + FcsClusterWithStarXYZ *cluSTAR = new FcsClusterWithStarXYZ(clu, fcsDb); + if ( clu->detectorId() == kFcsEcalNorthDetId || clu->detectorId() == kFcsEcalSouthDetId ){ + LOG_INFO << "Adding WCAL Cluster to FwdTree" << endm; + mTreeData.wcal.add( cluSTAR ); + } else if ( clu->detectorId() == kFcsHcalNorthDetId || clu->detectorId() == kFcsHcalSouthDetId ){ + LOG_INFO << "Adding HCAL Cluster to FwdTree" << endm; + mTreeData.hcal.add( cluSTAR ); + } + + delete cluSTAR; + } + + // LOAD ECAL / HCAL CLUSTERS + 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; + } + + // 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(){ + + 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 ); + } + + for ( size_t i = 0; i < muFttCollection->numberOfPoints(); i++ ){ + StMuFttPoint * c = muFttCollection->getPoint(i); + mTreeData.fttPoints.add( c ); + } + } +} diff --git a/StRoot/StFwdTrackMaker/StFwdQAMaker.h b/StRoot/StFwdTrackMaker/StFwdQAMaker.h new file mode 100644 index 00000000000..bc60f98d5b4 --- /dev/null +++ b/StRoot/StFwdTrackMaker/StFwdQAMaker.h @@ -0,0 +1,217 @@ +#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" +#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 + */ +template +class TClonesArrayWriter { + public: + TClonesArrayWriter() {} + ~TClonesArrayWriter() {} + + void createBranch( TTree *tree, const char* name, int buffSize = 256000, int splitLevel = 99){ + _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 = 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 = *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 StFcsDb; +class StFcsCluster; +class StFcsHit; +class StMuFcsCluster; +class StMuFcsHit; +class StMuFttCluster; +class StMuFttPoint; +class StMuFstHit; +class StMuFwdTrackSeedPoint; + + +/** + * @brief Store Cluster with STAR XYZ position + * + */ +class FcsClusterWithStarXYZ: public TObject { + public: + FcsClusterWithStarXYZ() { + mXYZ.SetXYZ(0, 0, 0); + mClu = nullptr; + } + FcsClusterWithStarXYZ( StMuFcsCluster *clu, StFcsDb *fcsDb ); + TVector3 mXYZ; + StMuFcsCluster *mClu; + ClassDef(FcsClusterWithStarXYZ, 1); +}; + +/** + * @brief Store Hit with STAR XYZ position + * + */ +class FcsHitWithStarXYZ: public TObject { + public: + FcsHitWithStarXYZ() { + mXYZ.SetXYZ(0, 0, 0); + mHit = nullptr; + } + FcsHitWithStarXYZ( StMuFcsHit *hit, StFcsDb *fcsDb ); + TVector3 mXYZ; + StMuFcsHit *mHit; + ClassDef(FcsHitWithStarXYZ, 1); +}; + + +/** @brief +* This class is a container for the data that will be written to the output tree. +*/ +struct FwdQATreeData { + + /** @brief Primary event vertex*/ + FwdTreeHeader header; + /** @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; + + + 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 FillFstPoints(); + void FillFcsStEvent(); + void FillFcsStMuDst(); + void FillTracks(); + void FillMcTracks(); + + protected: + TFile *mTreeFile = nullptr; + TTree *mTree = nullptr; + FwdQATreeData mTreeData; + + StEvent *mStEvent = nullptr; + StMuDstMaker *mMuDstMaker = nullptr; + StMuDst *mMuDst = nullptr; + StMuFwdTrackCollection * mMuForwardTrackCollection = nullptr; + StMuFcsCollection *mMuFcsCollection = nullptr; + StFwdTrackMaker *mFwdTrackMaker = nullptr; + StFcsDb *mFcsDb = nullptr; + +}; + + +#endif 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