From 054460e57256ed72a4c094ee3f796cd75e48d284 Mon Sep 17 00:00:00 2001 From: SEUNGHWAN LEE Date: Wed, 20 Nov 2024 05:09:31 -0500 Subject: [PATCH 1/2] RHICf Filter for event generator --- StRoot/StarGenerator/FILT/RHICfFilter.cxx | 209 ++++++++++++++++++++++ StRoot/StarGenerator/FILT/RHICfFilter.h | 51 ++++++ 2 files changed, 260 insertions(+) create mode 100644 StRoot/StarGenerator/FILT/RHICfFilter.cxx create mode 100644 StRoot/StarGenerator/FILT/RHICfFilter.h diff --git a/StRoot/StarGenerator/FILT/RHICfFilter.cxx b/StRoot/StarGenerator/FILT/RHICfFilter.cxx new file mode 100644 index 00000000000..ae53b43c9ec --- /dev/null +++ b/StRoot/StarGenerator/FILT/RHICfFilter.cxx @@ -0,0 +1,209 @@ +#include "RHICfFilter.h" + +#include "StarGenerator/EVENT/StarGenEvent.h" +#include "StarGenerator/EVENT/StarGenPPEvent.h" +#include "StarGenerator/EVENT/StarGenParticle.h" + +RHICfFilter::RHICfFilter( const char* name ) +: StarFilterMaker(name), mRHICfRunType(-1), mIsOnlyPi0(0), mIsOnlyNeutron(0), mIsPi0Event(0), mIsNeuEvent(0), mEnergyCut(20.) +{ + mRHICfPoly = 0; + memset(mRHICfTowerBoundary, 0, sizeof(mRHICfTowerBoundary)); + memset(mRHICfTowerCenterPos, 0, sizeof(mRHICfTowerCenterPos)); +} + +RHICfFilter::~RHICfFilter() +{ + if(mRHICfPoly){ + delete mRHICfPoly; + mRHICfPoly = 0; + } +} + +void RHICfFilter::SetRHICfRunType(int type) +{ + if(0 <= type && type <= 2){ + mRHICfRunType = type; + } +} + +void RHICfFilter::SetEnergyCut(double val){mEnergyCut = val;} +void RHICfFilter::SetOnlyPi0(){mIsOnlyPi0 = true;} +void RHICfFilter::SetOnlyNeutron(){mIsOnlyNeutron = true;} + +int RHICfFilter::Init() +{ + if(!InitRHICfGeometry()){return kStErr;} + if(mIsOnlyPi0 && mIsOnlyNeutron){return kStErr;} + return kStOk; +} + +int RHICfFilter::Filter( StarGenEvent *_event ) +{ + // Get a reference to the current event + StarGenEvent& event = (_event)? *_event : *mEvent; + + ClearEvent(); + StarGenParticle *part; + + // Loop over tracks to find particles of interest + int npart = event.GetNumberOfParticles(); + for ( int ipart=1; ipartGetId()) < 10){continue;} + if(part -> GetStatus() != StarGenParticle::kFinal){continue;} + + int pid = part -> GetId(); + double posX = part -> GetVx(); + double posY = part -> GetVy(); + double posZ = part -> GetVz(); + double px = part -> GetPx(); + double py = part -> GetPy(); + double pz = part -> GetPz(); + double e = part -> GetEnergy(); + + int hit = GetRHICfGeoHit(posX, posY, posZ, px, py, pz, e); + if(hit < 0){continue;} + + if(pid == 22){ // for RHICf gamma + mRHICfGammaIdx.push_back(ipart); + } + if(pid == 2112){ // for RHICf neutron + mIsNeuEvent = true; + } + } + + // Find a RHICf pi0 events + int RHICfgammaNum = mRHICfGammaIdx.size(); + for(int i=0; i GetFirstMother(); + + part = event[gamma1MotherIdx]; + int gamma1MotherPID = part -> GetId(); + int gamma1MotherNumber = part -> GetIndex(); + + if(gamma1MotherPID != 111){continue;} + + for(int j=i+1; j GetFirstMother(); + + part = event[gamma2MotherIdx]; + int gamma2MotherPID = part -> GetId(); + int gamma2MotherNumber = part -> GetIndex(); + + if(gamma2MotherPID != 111){continue;} + + if(gamma1MotherNumber == gamma2MotherNumber){ + mIsPi0Event = true; + } + } + } + + if(!mIsOnlyPi0 && !mIsOnlyNeutron){ + if(mIsPi0Event || mIsNeuEvent){ + return StarGenEvent::kAccept; + } + } + else if(mIsOnlyPi0 && mIsPi0Event){ + return StarGenEvent::kAccept; + } + else if(mIsOnlyNeutron && mIsNeuEvent){ + return StarGenEvent::kAccept; + } + return StarGenEvent::kReject; +} + +int RHICfFilter::InitRHICfGeometry() +{ + double tsDetSize = 20.; // [mm] + double tlDetSize = 40.; // [mm] + double detBoundCut = 2.; // [mm] + double distTStoTL = 47.4; // [mm] + + double detBeamCenter = 0.; + + if(mRHICfRunType < 0 || mRHICfRunType > 2){ + cout << "RHICfFilter::InitRHICfGeometry() warning!!! RHICf run type is not setted!!!" << endl; + return 0; + } + if(mRHICfRunType == 0){detBeamCenter = 0.;} // TS + if(mRHICfRunType == 1){detBeamCenter = -47.4;} // TL + if(mRHICfRunType == 2){detBeamCenter = 21.6;} // TOP + + mRHICfTowerBoundary[0][0][0] = sqrt(2)*((tsDetSize - detBoundCut*2.)/2.); + mRHICfTowerBoundary[0][0][1] = 0.; + mRHICfTowerBoundary[0][1][0] = 0.; + mRHICfTowerBoundary[0][1][1] = sqrt(2)*((tsDetSize - detBoundCut*2.)/2.); + mRHICfTowerBoundary[0][2][0] = -1.*sqrt(2)*((tsDetSize - detBoundCut*2.)/2.); + mRHICfTowerBoundary[0][2][1] = 0.; + mRHICfTowerBoundary[0][3][0] = 0.; + mRHICfTowerBoundary[0][3][1] = -1.*sqrt(2)*((tsDetSize - detBoundCut*2.)/2.); + + mRHICfTowerBoundary[1][0][0] = sqrt(2)*((tlDetSize - detBoundCut*2.)/2.); + mRHICfTowerBoundary[1][0][1] = 0.; + mRHICfTowerBoundary[1][1][0] = 0.; + mRHICfTowerBoundary[1][1][1] = sqrt(2)*((tlDetSize - detBoundCut*2.)/2.); + mRHICfTowerBoundary[1][2][0] = -1.*sqrt(2)*((tlDetSize - detBoundCut*2.)/2.); + mRHICfTowerBoundary[1][2][1] = 0.; + mRHICfTowerBoundary[1][3][0] = 0.; + mRHICfTowerBoundary[1][3][1] = -1.*sqrt(2)*((tlDetSize - detBoundCut*2.)/2.); + + mRHICfTowerCenterPos[0] = detBeamCenter; + mRHICfTowerCenterPos[1] = distTStoTL + detBeamCenter; + + mRHICfPoly = new TH2Poly(); + mRHICfPoly -> SetName("RHICfPoly"); + mRHICfPoly -> SetStats(0); + + double x[4]; + double y[4]; + + for(int t=0; t<2; t++){ + for(int i=0; i<4; i++){ + double xPos = mRHICfTowerBoundary[t][i][0]; + double yPos = mRHICfTowerCenterPos[t] + mRHICfTowerBoundary[t][i][1]; + x[i] = xPos; + y[i] = yPos; + } + mRHICfPoly -> AddBin(4, x, y); + } + if(!mRHICfPoly){return 0;} + return 1; +} + +int RHICfFilter::GetRHICfGeoHit(double posX, double posY, double posZ, double px, double py, double pz, double e) +{ + if(e < mEnergyCut){return -1;} // energy cut + + double momMag = sqrt(px*px + py*py + pz*pz); + double unitVecX = px/momMag; + double unitVecY = py/momMag; + double unitVecZ = pz/momMag; + + if(unitVecZ < 0){return -1;} // opposite side cut + + double z = mRHICfDetZ - posZ; + if(z < 0.){return -1;} // create z-position cut + + double x = z * (unitVecX/unitVecZ) + posX; + double y = z * (unitVecY/unitVecZ) + posY; + + int type = mRHICfPoly -> FindBin(x, y); + if(type < 1 || type > 2){return -1;} // RHICf geometrical hit cut + + return type; +} + +void RHICfFilter::ClearEvent() +{ + mRHICfGammaIdx.clear(); + mIsPi0Event = false; + mIsNeuEvent = false; +} \ No newline at end of file diff --git a/StRoot/StarGenerator/FILT/RHICfFilter.h b/StRoot/StarGenerator/FILT/RHICfFilter.h new file mode 100644 index 00000000000..02de06c96c6 --- /dev/null +++ b/StRoot/StarGenerator/FILT/RHICfFilter.h @@ -0,0 +1,51 @@ +#ifndef __RHICfFilter_h__ +#define __RHICfFilter_h__ + +#include +#include "TH2Poly.h" + +/*! + \class RHICfFilter + \brief RHICf pi0 and Neutron events filter as RHICf Run type + */ + +#include "StarFilterMaker.h" + +class RHICfFilter : public StarFilterMaker +{ +public: + RHICfFilter( const char* name = "rhicffilt" ); + virtual ~RHICfFilter(); + + void SetRHICfRunType(int type); // [0=TS, 1=TL, 2=TOP] + void SetEnergyCut(double val); + void SetOnlyPi0(); + void SetOnlyNeutron(); + + int Init(); + int Filter( StarGenEvent *event = 0 ); + +private: + int InitRHICfGeometry(); + int GetRHICfGeoHit(double posX, double posY, double posZ, double px, double py, double pz, double e); + void ClearEvent(); // clear the temporary array in events + + int mRHICfRunType; // [0=TS, 1=TL, 2=TOP] + bool mIsOnlyPi0; + bool mIsOnlyNeutron; + + TH2Poly* mRHICfPoly; // only west + double mRHICfTowerBoundary[2][4][2]; // [TS, TL][bound square][x, y] + double mRHICfTowerCenterPos[2]; // [TS, TL] y pos + + vector mRHICfGammaIdx; + bool mIsPi0Event; + bool mIsNeuEvent; + + double mEnergyCut; // [GeV] + double mRHICfDetZ = 17800.; // [mm] + + ClassDef(RHICfFilter,0); +}; + +#endif \ No newline at end of file From b9a58cbea16566f525d8ad39dd9379ab4912e825 Mon Sep 17 00:00:00 2001 From: Seunghwan Lee Date: Fri, 13 Dec 2024 03:46:59 -0500 Subject: [PATCH 2/2] Change the unit for RHICf geometry mm to cm --- StRoot/StarGenerator/FILT/RHICfFilter.cxx | 14 +++++++------- StRoot/StarGenerator/FILT/RHICfFilter.h | 2 +- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/StRoot/StarGenerator/FILT/RHICfFilter.cxx b/StRoot/StarGenerator/FILT/RHICfFilter.cxx index ae53b43c9ec..12cbc89e80a 100644 --- a/StRoot/StarGenerator/FILT/RHICfFilter.cxx +++ b/StRoot/StarGenerator/FILT/RHICfFilter.cxx @@ -122,20 +122,20 @@ int RHICfFilter::Filter( StarGenEvent *_event ) int RHICfFilter::InitRHICfGeometry() { - double tsDetSize = 20.; // [mm] - double tlDetSize = 40.; // [mm] - double detBoundCut = 2.; // [mm] - double distTStoTL = 47.4; // [mm] + double tsDetSize = 2.; // [cm] + double tlDetSize = 4.; // [cm] + double detBoundCut = 0.2; // [cm] + double distTStoTL = 4.74; // [cm] - double detBeamCenter = 0.; + double detBeamCenter = 0.; // [cm] if(mRHICfRunType < 0 || mRHICfRunType > 2){ cout << "RHICfFilter::InitRHICfGeometry() warning!!! RHICf run type is not setted!!!" << endl; return 0; } if(mRHICfRunType == 0){detBeamCenter = 0.;} // TS - if(mRHICfRunType == 1){detBeamCenter = -47.4;} // TL - if(mRHICfRunType == 2){detBeamCenter = 21.6;} // TOP + if(mRHICfRunType == 1){detBeamCenter = -4.74;} // TL + if(mRHICfRunType == 2){detBeamCenter = 2.16;} // TOP mRHICfTowerBoundary[0][0][0] = sqrt(2)*((tsDetSize - detBoundCut*2.)/2.); mRHICfTowerBoundary[0][0][1] = 0.; diff --git a/StRoot/StarGenerator/FILT/RHICfFilter.h b/StRoot/StarGenerator/FILT/RHICfFilter.h index 02de06c96c6..fa1a575b50a 100644 --- a/StRoot/StarGenerator/FILT/RHICfFilter.h +++ b/StRoot/StarGenerator/FILT/RHICfFilter.h @@ -43,7 +43,7 @@ class RHICfFilter : public StarFilterMaker bool mIsNeuEvent; double mEnergyCut; // [GeV] - double mRHICfDetZ = 17800.; // [mm] + double mRHICfDetZ = 1780.; // [cm] ClassDef(RHICfFilter,0); };