From d993d8e70481ea28b0307eea3dc8315f0cb40c3f Mon Sep 17 00:00:00 2001 From: Daniel Brandenburg Date: Wed, 14 Aug 2024 10:38:55 -0400 Subject: [PATCH] Update StEvent structs for StFttPoint StFttCluster and update StFttFastSimMaker (#698) 1. Updates StEvent containers for Ftt data (clusters and points) to include idTruth (mirror FST and FCS) 2. Update StFttFastSimMaker to copy Ftt GEANT hits to StFttCollection as StFttPoints The Fast sim just copies GEANT hits from front modules and writes them into the StFttCollection as StFttPoints (skips the entire reco chain). The "slow" sim (still WIP) will write clusters and therefore idTruth is needed there. Removed old code based on square modules (RnD geometry) and different strip layout. This also simplifies the Fwd tracking because it no longer needs to read directly from GEANT --- StRoot/StEvent/StFttCluster.h | 5 +- StRoot/StEvent/StFttPoint.h | 6 +- StRoot/StFttSimMaker/StFttFastSimMaker.cxx | 432 ++------------------- StRoot/StFttSimMaker/StFttFastSimMaker.h | 96 +---- 4 files changed, 53 insertions(+), 486 deletions(-) diff --git a/StRoot/StEvent/StFttCluster.h b/StRoot/StEvent/StFttCluster.h index 4b1b8ce6515..d478c798f32 100644 --- a/StRoot/StEvent/StFttCluster.h +++ b/StRoot/StEvent/StFttCluster.h @@ -27,6 +27,7 @@ class StFttCluster : public StObject { float sumAdc() const; float x() const; // Mean x ("center of gravity") in local grid coordinate (1st moment). float sigma() const; // Maximum 2nd moment (along major axis). + UShort_t idTruth() const { return mIdTruth; } // Get the truth ID void setId(int cluid); void setPlane(UChar_t plane); @@ -37,6 +38,7 @@ class StFttCluster : public StObject { void setSumAdc(int theSumAdc); void setX(float x0); void setSigma(float sigma); + void setIdTruth(UShort_t id) { mIdTruth = id; } StPtrVecFttRawHit& rawHits(); const StPtrVecFttRawHit& rawHits() const; @@ -62,8 +64,9 @@ class StFttCluster : public StObject { StPtrVecFttRawHit mRawHits; // Tower hits of the current cluster StPtrVecFttCluster mNeighbors; // Neighbor clusters StPtrVecFttPoint mPoints; // Fitted points (photons) in the cluster + UShort_t mIdTruth=0; // Truth ID - ClassDef(StFttCluster, 2) + ClassDef(StFttCluster, 3) }; std::ostream& operator << ( std::ostream&, const StFttCluster& clu ); // Printing operator diff --git a/StRoot/StEvent/StFttPoint.h b/StRoot/StEvent/StFttPoint.h index 71bcbd8cd79..bab29a1971f 100644 --- a/StRoot/StEvent/StFttPoint.h +++ b/StRoot/StEvent/StFttPoint.h @@ -31,6 +31,7 @@ class StFttPoint : public StObject { int nClusters() const; // Number of points in the parent cluster. StFttCluster* cluster( size_t i); // Parent cluster of the photon. const StThreeVectorD& xyz() const; // XYZ position in global STAR coordinate + UShort_t idTruth() const { return mIdTruth; } // Get the truth ID void setPlane(UChar_t plane); void setQuadrant(UChar_t quad); @@ -38,7 +39,7 @@ class StFttPoint : public StObject { void setY(float y); void addCluster(StFttCluster* cluster, UChar_t dir); void setXYZ(const StThreeVectorD& p3); - + void setIdTruth(UShort_t id) { mIdTruth = id; } void print(int option=0); @@ -49,8 +50,9 @@ class StFttPoint : public StObject { Float_t mY=0.0; // y-position in local coordinate StFttCluster *mClusters[4]; StThreeVectorD mXYZ; // Photon position in STAR coordinate + UShort_t mIdTruth=0; // Truth ID - ClassDef(StFttPoint, 1) + ClassDef(StFttPoint, 2) }; inline UChar_t StFttPoint::plane() const { return mPlane; } diff --git a/StRoot/StFttSimMaker/StFttFastSimMaker.cxx b/StRoot/StFttSimMaker/StFttFastSimMaker.cxx index da9b723f741..607a0333cfd 100644 --- a/StRoot/StFttSimMaker/StFttFastSimMaker.cxx +++ b/StRoot/StFttSimMaker/StFttFastSimMaker.cxx @@ -4,46 +4,19 @@ #include "StEvent/StEvent.h" #include "St_base/StMessMgr.h" -#include "StEvent/StRnDHit.h" -#include "StEvent/StRnDHitCollection.h" #include "StThreeVectorF.hh" -#include "TCanvas.h" -#include "TCernLib.h" -#include "TH2F.h" -#include "TLine.h" -#include "TString.h" -#include "TVector3.h" #include "tables/St_g2t_fts_hit_Table.h" -#include "tables/St_g2t_track_Table.h" -#include +#include "StEvent/StFttCollection.h" +#include "StEvent/StFttPoint.h" #include "StarGenerator/UTIL/StarRandom.h" -namespace FttGlobal { - const bool verbose = false; -} StFttFastSimMaker::StFttFastSimMaker(const Char_t *name) - : StMaker{name}, - hGlobalYX(0), - hOctantYX(0), - hOctantWireYX(0), - hOctantStripYX(0), - hWireDeltasX(0), - hWireDeltasY(0), - hStripDeltasX(0), - hStripDeltasY(0), - hWirePullsX(0), - hWirePullsY(0), - hStripPullsX(0), - hStripPullsY(0), - hPointsPullsX(0), - hPointsPullsY(0) {} + : StMaker{name} {} int StFttFastSimMaker::Init() { - iEvent = 0; - return StMaker::Init(); } @@ -58,373 +31,46 @@ Int_t StFttFastSimMaker::Make() { LOG_DEBUG << "Creating StEvent" << endm; } - if (0 == event->rndHitCollection()) { - event->setRnDHitCollection(new StRnDHitCollection()); - LOG_DEBUG << "Creating StRnDHitCollection for FTS" << endm; - } - - FillThinGapChambers(event); - iEvent++; - - return kStOk; -} - -/** - * Maps a global hit to a local coordinate system for a given quadrant - * The quadrants are numbered clockwise as: - * 0 1 - * 3 2 - * Does NOT support rotations. Maybe added later if we need to - * The coordinate system is x positive to the right and y positive up (top coords are greater than bottom coords) - */ -void StFttFastSimMaker::GlobalToLocal(float x, float y, int disk, int &quad, float &localX, float &localY) { - // quad RECT - float qr = -1; - float ql = -1; - float qb = -1; - float qt = -1; - - QuadBottomLeft(disk, 0, qb, ql); - qr = ql + STGC_QUAD_WIDTH; - qt = qb + STGC_QUAD_HEIGHT; - if (x >= ql && x < qr && y >= qb && y < qt) { - quad = 0; - localX = x - ql; - localY = y - qb; - } - - QuadBottomLeft(disk, 1, qb, ql); - qr = ql + STGC_QUAD_WIDTH; - qt = qb + STGC_QUAD_HEIGHT; - if (x >= ql && x < qr && y >= qb && y < qt) { - quad = 1; - localX = x - ql; - localY = y - qb; + if ( event->fttCollection() == nullptr ){ + LOG_INFO << "Creating FttCollection" << endm; + StFttCollection *fttcollection = new StFttCollection(); + event->setFttCollection(fttcollection); } - QuadBottomLeft(disk, 2, qb, ql); - qr = ql + STGC_QUAD_WIDTH; - qt = qb + STGC_QUAD_HEIGHT; - if (x >= ql && x < qr && y >= qb && y < qt) { - quad = 2; - localX = x - ql; - localY = y - qb; - } - QuadBottomLeft(disk, 3, qb, ql); - qr = ql + STGC_QUAD_WIDTH; - qt = qb + STGC_QUAD_HEIGHT; - if (x >= ql && x < qr && y >= qb && y < qt) { - quad = 3; - localX = x - ql; - localY = y - qb; - } - - return; -} - -float StFttFastSimMaker::DiskOffset(int disk) { - assert(disk >= 9 && disk <= 12); - if (disk == 9) - return 10; - if (disk == 10) - return 11; - if (disk == 11) - return 12; - if (disk == 12) - return 13; - return 10; -} - -float StFttFastSimMaker::DiskRotation(int disk) { - assert(disk >= 9 && disk <= 12); - // these are - if (disk == 9) - return this->sTGC_disk9_theta; - if (disk == 10) - return this->sTGC_disk10_theta; - if (disk == 11) - return this->sTGC_disk11_theta; - if (disk == 12) - return this->sTGC_disk12_theta; - return 0; -} - -void StFttFastSimMaker::QuadBottomLeft(int disk, int quad, float &bottom, float &left) { - float hbp = DiskOffset(disk); - if ( FttGlobal::verbose ) {LOG_INFO << "disk: " << disk << ", offset = " << hbp << endm;} - - // quad 0 RECT - float q0l = hbp - STGC_QUAD_WIDTH; - float q0b = hbp; - - float q1l = hbp; - float q1b = -hbp; - float q2l = -hbp; - float q2b = -hbp - STGC_QUAD_HEIGHT; - - float q3l = -hbp - STGC_QUAD_WIDTH; - float q3b = hbp - STGC_QUAD_HEIGHT; - - if (0 == quad) { - bottom = q0b; - left = q0l; - } - if (1 == quad) { - bottom = q1b; - left = q1l; + St_g2t_fts_hit *g2t_stg_hits = (St_g2t_fts_hit *)GetDataSet("geant/g2t_stg_hit"); + // size_t numFwdHitsPrior = mFwdHitsFtt.size(); + if (!g2t_stg_hits){ + LOG_WARN << "geant/g2t_stg_hit is empty" << endm; + return kStOk; } - if (2 == quad) { - bottom = q2b; - left = q2l; - } - if (3 == quad) { - bottom = q3b; - left = q3l; - } - return; -} - -/** - * Map a local coordinate back to global coords - */ -void StFttFastSimMaker::LocalToGlobal(float localX, float localY, int disk, int quad, float &globalX, float &globalY) { - // quad RECT - float qb = -1; - float ql = -1; - QuadBottomLeft(disk, quad, qb, ql); - - globalX = localX + ql; - globalY = localY + qb; - return; -} - -/** - * Checks if a vertical (face=0) and horizontal (face=1) wires are overlapping - * used for determining if ghost hits can be created from the overlapped wires - * - */ -bool StFttFastSimMaker::Overlaps(StRnDHit *hitA, StRnDHit *hitB) { - // require that they are in the same disk! - if (hitA->layer() != hitB->layer()) - return false; - int disk = hitA->layer(); - // require that they are in the same quadrant of the detector - if (hitA->wafer() != hitB->wafer()) - return false; - int quad = hitA->wafer(); - - float x1 = hitA->double2(); - float y1 = hitA->double3(); - float x2 = hitB->double2(); - float y2 = hitB->double3(); - - float b = -1, l = -1; - QuadBottomLeft(disk, quad, b, l); - - float lx1 = x1 - l; - float ly1 = y1 - b; - float lx2 = x2 - l; - float ly2 = y2 - b; - - int chunkx1 = lx1 / STGC_WIRE_LENGTH; - int chunky1 = ly1 / STGC_WIRE_LENGTH; - - int chunkx2 = lx2 / STGC_WIRE_LENGTH; - int chunky2 = ly2 / STGC_WIRE_LENGTH; - - if (chunkx1 != chunkx2) - return false; - if (chunky1 != chunky2) - return false; - return true; -} - -void StFttFastSimMaker::FillThinGapChambers(StEvent *event) { - // Read the g2t table - St_g2t_fts_hit *hitTable = static_cast(GetDataSet("g2t_stg_hit")); - if (!hitTable) { - LOG_INFO << "g2t_stg_hit table is empty" << endm; - return; - } // if !hitTable - - StRnDHitCollection *ftscollection = event->rndHitCollection(); - - std::vector hits; - - // Prepare to loop over all hits - const int nhits = hitTable->GetNRows(); - const g2t_fts_hit_st *hit = hitTable->GetTable(); - - StarRandom &rand = StarRandom::Instance(); - float dx = STGC_SIGMA_X; - float dy = STGC_SIGMA_Y; - float dz = STGC_SIGMA_Z; - - int nSTGCHits = 0; - sTGCNRealPoints = 0; - sTGCNGhostPoints = 0; - - for (int i = 0; i < nhits; i++) { - hit = (g2t_fts_hit_st *)hitTable->At(i); - if (0 == hit) + const double sigXY = 0.01; + int nstg = g2t_stg_hits->GetNRows(); + LOG_DEBUG << "This event has " << nstg << " stg hits in geant/g2t_stg_hit " << endm; + for (int i = 0; i < nstg; i++) { + g2t_fts_hit_st *git = (g2t_fts_hit_st *)g2t_stg_hits->At(i); + if (0 == git) + continue; // geant hit + int track_id = git->track_p; + int volume_id = git->volume_id; + int plane_id = (volume_id - 1) / 100; // from 1 - 16. four chambers per station + + // only use the hits on the front modules + if ( volume_id % 2 ==0 ) continue; - float xhit = hit->x[0]; - float yhit = hit->x[1]; - float zhit = hit->x[2]; - int volume_id = hit->volume_id; - // volume_id = (1 front | 2 back) + 10 * (quadrant 0-3) + 100 * (station 0-4) - int disk = ((volume_id - 1) / 100) + 9 ; - LOG_DEBUG << "sTGC hit: volume_id = " << volume_id << " disk = " << disk << endm; + float x = git->x[0] + gRandom->Gaus(0, sigXY); // 100 micron blur according to approx sTGC reso + float y = git->x[1] + gRandom->Gaus(0, sigXY); // 100 micron blur according to approx sTGC reso + float z = git->x[2]; - // Now that geometry has a front and back, we skip points on the back module for fast sim - if (disk < 9 || volume_id % 2 == 0) - continue; - - float theta = DiskRotation(disk); - - float x_blurred = xhit + rand.gauss( STGC_SIGMA_X); - float y_blurred = yhit + rand.gauss( STGC_SIGMA_Y); - - float x_rot = -999, y_rot = -999; - this->rot(-theta, x_blurred, y_blurred, x_rot, y_rot); - - int quad = -1; - float localX = -999, localY = -999; - GlobalToLocal(x_rot, y_rot, disk, quad, localX, localY); - - // not in the active region - if (quad < 0 || quad > 3) - continue; - nSTGCHits++; - - StRnDHit *ahit = new StRnDHit(); - - ahit->setPosition({x_blurred, y_blurred, zhit}); - ahit->setPositionError({dx, dy, 0.1}); - - ahit->setDouble0(xhit); - ahit->setDouble1(yhit); - ahit->setDouble2(x_rot); - ahit->setDouble3(y_rot); - - ahit->setLayer(disk); // disk mapped to layer - ahit->setLadder(2); // indicates a point - ahit->setWafer(quad); // quadrant number - - ahit->setIdTruth(hit->track_p, 0); - ahit->setDetectorId(kFtsId); // TODO: use dedicated ID for Ftt when StEvent is updated - - float Ematrix[] = { - dx * dx, 0.f, 0.f, - 0.f, dy * dy, 0.f, - 0.f, 0, 0.f, dz * dz}; - - ahit->setErrorMatrix(Ematrix); - hits.push_back(ahit); - - if (!STGC_MAKE_GHOST_HITS) { - // Make this "REAL" hit. - if (FttGlobal::verbose){ - ahit->Print(); - } - ftscollection->addHit(ahit); - sTGCNRealPoints++; - } - } - - if (STGC_MAKE_GHOST_HITS) { - for (auto &hit0 : hits) { // first loop on hits - float hit0_x = hit0->double2(); - float hit0_y = hit0->double3(); - int disk0 = hit0->layer(); - int quad0 = hit0->wafer(); - float theta = DiskRotation(disk0); - - for (auto &hit1 : hits) { // second loop on hits - float hit1_x = hit1->double2(); - float hit1_y = hit1->double3(); - int disk1 = hit1->layer(); - int quad1 = hit1->wafer(); - - if (disk0 != disk1) - continue; - if (quad0 != quad1) - continue; - - // check on overlapping segments - if (false == Overlaps(hit0, hit1)) - continue; - - float x = hit0_x; - float y = hit1_y; - - int qaTruth = 0; - int idTruth = 0; - if (hit1_x == hit0_x && hit1_y == hit0_y) { - sTGCNRealPoints++; - qaTruth = 1; - idTruth = hit0->idTruth(); - } else { - sTGCNGhostPoints++; - qaTruth = 0; - } - - float rx = -999, ry = -999; - this->rot(theta, x, y, rx, ry); - // the trick here is that rotations (in 2D) will commute - // so the earlier -theta rotation and this +theta - // rotation cancel for real hits - // but not so for ghost hits - - StRnDHit *ahit = new StRnDHit(); - - ahit->setPosition({rx, ry, hit0->position().z()}); - ahit->setPositionError({dx, dy, 0.1}); - - ahit->setDouble0(hit0_x); - ahit->setDouble1(hit0_y); - ahit->setDouble2(hit1_x); - ahit->setDouble3(hit1_y); - - ahit->setLayer(disk0); // disk mapped to layer - ahit->setLadder(2); // indicates a point - ahit->setWafer(quad0); // quadrant number - - ahit->setIdTruth(idTruth, qaTruth); - ahit->setDetectorId(kFtsId); // TODO: use dedicated ID for Ftt when StEvent is updated - - float Ematrix[] = { - dx * dx, 0.f, 0.f, - 0.f, dy * dy, 0.f, - 0.f, 0, 0.f, dz * dz}; - ahit->setErrorMatrix(Ematrix); - ftscollection->addHit(ahit); - - } // loop hit1 - } // loop hit0 - - - // in this case the hits used in the original array were not saved, but copied so we need to delete them - - for ( size_t i = 0; i < hits.size(); i++ ){ - delete hits[i]; - hits[i] = nullptr; - } - } // make Ghost Hits - - if (FttGlobal::verbose) { - LOG_INFO << "nHits (all FTS) = " << nhits << endm; - } - if (FttGlobal::verbose) { - LOG_INFO << "nSTGC = " << nSTGCHits << endm; - } - if (FttGlobal::verbose) { - LOG_INFO << "nReal = " << sTGCNRealPoints << endm; - } - if (FttGlobal::verbose) { - LOG_INFO << "nGhost = " << sTGCNGhostPoints << endm; - } - -} // fillThinGap + StFttPoint *point = new StFttPoint(); + point->setPlane(plane_id); + point->setQuadrant(0); // TODO this could be improved, but it is not used in the current implementation + StThreeVectorD xyz; + xyz.set(x, y, z); + point->setXYZ( xyz ); + point->setIdTruth( track_id ); + event->fttCollection()->addPoint(point); + } // loop on hits + return kStOk; +} diff --git a/StRoot/StFttSimMaker/StFttFastSimMaker.h b/StRoot/StFttSimMaker/StFttFastSimMaker.h index 5a5fe10a7d4..7b10bfc9ccd 100644 --- a/StRoot/StFttSimMaker/StFttFastSimMaker.h +++ b/StRoot/StFttSimMaker/StFttFastSimMaker.h @@ -1,10 +1,6 @@ -#ifndef ST_FTT_FAST_SIM_MAKER_H -#define ST_FTT_FAST_SIM_MAKER_H - -class g2t_emc_hit_st; -class StFtsHit; -class StEvent; +#ifndef ST_FTT_FASTER_SIM_MAKER_H +#define ST_FTT_FASTER_SIM_MAKER_H #include "StChain/StMaker.h" #include @@ -14,99 +10,19 @@ class StEvent; #include "TH2F.h" #include "TNtuple.h" -class StRnDHit; class StFttFastSimMaker : public StMaker { public: - explicit StFttFastSimMaker(const Char_t *name = "fttSim"); - virtual ~StFttFastSimMaker() {} + StFttFastSimMaker(const Char_t *name = "fttSim"); + ~StFttFastSimMaker() {} Int_t Make(); int Init(); int Finish() { return kStOk; } - virtual const char *GetCVS() const; - - void SetDiskRotation(int disk, float degrees) { - - const float deg_to_radians = 0.017453292f; // = 3.1415926 / 180.0; - if (9 == disk) - sTGC_disk9_theta = degrees * deg_to_radians; - else if (10 == disk) - sTGC_disk10_theta = degrees * deg_to_radians; - else if (11 == disk) - sTGC_disk11_theta = degrees * deg_to_radians; - else if (12 == disk) - sTGC_disk12_theta = degrees * deg_to_radians; - return; - } - - private: - void FillThinGapChambers(StEvent *event); - - int iEvent; - - TH2F *hGlobalYX; - TH2F *hOctantYX; - - TH2F *hOctantWireYX; - TH2F *hOctantStripYX; - - TH2F *hWireDeltasX; - TH2F *hWireDeltasY; - TH2F *hStripDeltasX; - TH2F *hStripDeltasY; - - TH2F *hWirePullsX; - TH2F *hWirePullsY; - TH2F *hStripPullsX; - TH2F *hStripPullsY; - - TH2F *hPointsPullsX; - TH2F *hPointsPullsY; - - //table to keep pointer to hit for each disc, r & phi strips - - // convert x, y to quandrant and local X, Y - // quadrants are - // 0 1 - // 3 2 - void GlobalToLocal(float x, float y, int disk, int &quad, float &localX, float &localY); - void LocalToGlobal(float localX, float localY, int disk, int quad, float &globalX, float &globalY); - bool Overlaps(StRnDHit *hitA, StRnDHit *hitB); - void QuadBottomLeft(int disk, int quad, float &bottom, float &left); - float DiskOffset(int disk); - float DiskRotation(int disk); - - void rot(float theta, float x, float y, float &xp, float &yp) { - xp = x * cos(theta) - y * sin(theta); - yp = x * sin(theta) + y * cos(theta); - } - - const double STGC_BEAM_CUT_OUT = 6.0; // cm - const double STGC_QUAD_WIDTH = 60.0; // cm - const double STGC_QUAD_HEIGHT = 60.0; // cm - const double STGC_WIRE_WIDTH = 0.32; // cm - const double STGC_SIGMA_X = 0.01; // 100 microns - const double STGC_SIGMA_Y = 0.01; // 100 microns - const double STGC_SIGMA_Z = 0.001; // 10 microns - const double STGC_WIRE_LENGTH = 15.0; // cm - const bool STGC_MAKE_GHOST_HITS = true; //should be moved to run-time opt - - float sTGC_disk9_theta = 0.0f; - float sTGC_disk10_theta = 0.0f; - float sTGC_disk11_theta = 0.0f; - float sTGC_disk12_theta = 0.0f; - - int sTGCNRealPoints = 0; - int sTGCNGhostPoints = 0; - ClassDef(StFttFastSimMaker, 0) + ClassDef(StFttFastSimMaker, 0); }; -inline const char *StFttFastSimMaker::GetCVS() const { - static const char cvs[] = "Tag $Name: $ $Id: StFttFastSimMaker.h,v 1.1 2021/03/26 14:11:40 jdb Exp $ built " __DATE__ " " __TIME__; - return cvs; -} -#endif +#endif \ No newline at end of file