Home » Blogs » slee5's blog

RHICf pi0 and Neutron event selection code

Updated on Wed, 2024-09-04 07:09. Originally created by slee5 on 2024-09-04 06:51.
// RHICf geometrical acceptance selection code for RHICf Pi0 event Event
 
#include <TH2Poly.h>
 
// Initialization RHICf geometry acceptance
TH2Poly* InitRHICfGeo(int runType);
 
// RHICfHitPos function convert the particle position at RHICf detector plane (18 m)
double RHICfPlanePos(int xy, double xProd, double yProd, double zProd, double px, double py, double pz)
 
 
// main selection code
void RHICfEventSelection()
{
    int RHICfRunType = 0; // [TS=0, TL=1, TOP=2]
    TH2Poly* RHICfgeo = InitRHICfGeo(RHICfRunType)
 
    // temporary parameters for your convenience
    enum tempParameters{
        kGamma = 22; // pdg 
        kPi0 = 111; // pdg
    };
 
    // init temporary RHICf gamma array
    int tmpRHICfHitGammaIdx[20];
    memset(tmpRHICfHitGammaIdx, 0, sizeof(tmpRHICfHitGammaIdx));
    int tmpRHICfHitGammaNum = 0;
    
 
    // main particle loop in events
    int particleNumInEvents = 0;
    for (int par = 0; par < particleNumInEvents; par++) { 
 
        // event[par] ptr will be PYTHIA event data structrue
        
        // temporary variables 
        bool isParticle_final = event[par].isFinal(); // just only final state
        int particle_pid = event[par].id();           // pdg encoding
        double particle_posX = event[par].xProd();    // particle originate position x
        double particle_posY = event[par].yProd();    // particle originate position y
        double particle_posZ = event[par].zProd();    // particle originate position z
        double particle_px = event[par].px();         // particle momentum px
        double particle_py = event[par].py();         // particle momentum py
        double particle_pz = event[par].pz();         // particle momentum pz
        double particle_e = event[par].e();           // particle energy
 
        if( !isParticle_final ){continue;} // if not a final state, just skip
        if( particle_pid != kGamma ){continue;} // it particle not a gamma, just skip
 
        // Convert the Gamma position on transverse plane at RHICf 18 m 
        double gamma_posX = RHICfPlanePos(0, particle_posX, particle_posY, particle_posZ, particle_px, particle_py, particle_pz);
        double gamma_posY = RHICfPlanePos(1, particle_posX, particle_posY, particle_posZ, particle_px, particle_py, particle_pz);
 
        if(gamma_posX < -900 || gamma_posY < -900){continue;} // gamma position is significantly out of RHICf detector plane
        
        int RHICftowerHitIdx = RHICfgeo -> FindBin(gamma_posX, gamma_posY);
 
        int gammaType = -1;
        if(bin == 1){gammaType = 0;} // hitted RHICf TS tower
        if(bin == 2){gammaType = 1;} // hitted RHICf TL tower
        if(gammaType == -1){continue;}
 
        // save the RHICf hit gamma index
        tmpRHICfHitGammaIdx[tmpRHICfHitGammaNum] = par;
        tmpRHICfHitGammaNum++;
    }
 
    bool isRHICfPi0Event = false;
 
    // Find a two RHICf gamma from truth pi0
    for(int i=0; i<tmpRHICfHitGammaNum; i++){
        int gamma1Idx = tmpRHICfHitGammaIdx[i];
 
        int gamma1MotherIdx = event[gamma1Idx].mother1();
        int gamma1MotherPID = event[gamma1MotherIdx].id();
        if(gamma1MotherPID != kPi0){continue;} // if gamma1 is not from pi0, skip
 
        for(int j=i+1; j<tmpRHICfHitGammaNum; j++){
            int gamma2Idx = tmpRHICfGammaIdx[j];
 
            int gamma2MotherIdx = event[gamma2Idx].mother1();
            int gamma2MotherPID = event[gamma2MotherIdx].id();
            if(gamma2MotherPID != kPi0){continue;} // if gamma2 is not from pi0, skip
 
            if(gamma1MotherIdx == gamma2MotherIdx){
                isRHICfPi0Event = true; // RHICf pi0 event has found
            }
        }
    }
 
    if(isRHICfPi0Event == false){continue;} // if not a RHICf pi0 event, skip
 
    // it is RHICf pi0 Events !!
}
 
TH2Poly* InitRHICfGeo(int runType)
{
    double mRHICfTowerBoundary[2][4][2]; // [TS, TL][bound square][x, y]
    double mRHICfTowerCenterPos[2]; // [TS, TL] y pos
 
    double tsDetSize = 20.; // [mm]
    double tlDetSize = 40.; // [mm]
    double detBoundCut = 2.; // [mm]
    double distTStoTL = 47.4; // [mm]
 
    double detBeamCenter = 0.;
    if(runType == 0){detBeamCenter = 0.;} // TS
    if(runType == 1){detBeamCenter = -47.4;} // TL
    if(runType == 2){detBeamCenter = 21.6;} // TOP
 
    mRHICfTowerBoundary[0][0][0] = sqrt(2)*((tsDetSize - detBoundCut*2.)/2.); // right point
    mRHICfTowerBoundary[0][0][1] = 0.; // right point
    mRHICfTowerBoundary[0][1][0] = 0.; // up point
    mRHICfTowerBoundary[0][1][1] = sqrt(2)*((tsDetSize - detBoundCut*2.)/2.); // up point
    mRHICfTowerBoundary[0][2][0] = -1.*sqrt(2)*((tsDetSize - detBoundCut*2.)/2.); // left point
    mRHICfTowerBoundary[0][2][1] = 0.; // left point
    mRHICfTowerBoundary[0][3][0] = 0.; // down point
    mRHICfTowerBoundary[0][3][1] = -1.*sqrt(2)*((tsDetSize - detBoundCut*2.)/2.); // down point
 
    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;
 
    TH2Poly* rhicfGeo = new TH2Poly();
 
    double mX[4];
    double mY[4];
 
    for(int t=0; t<2; t++){ // tower
        for(int i=0; i<4; i++){ // edge points
            double xPos = mRHICfTowerBoundary[t][i][0];
            double yPos = mRHICfTowerCenterPos[t] + mRHICfTowerBoundary[t][i][1];
            mX[i] = xPos;
            mY[i] = yPos;
        }
        rhicfGeo -> AddBin(4, mX, mY);
    }
 
    return rhicfGeo;
}
 
double RHICfPlanePos(int xy, double xProd, double yProd, double zProd, double px, double py, double pz)
{
    double RHICfZ = 17800.0; // [mm] RHICf detector installed at about 18 m from IP
 
    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 -999;} // cut opposit-side (East-side region cut, RHICf detector is in only west-side)
 
    double z = RHICfZ - zProd;
    if(z < 0.){return -999;}
 
    double x = z * (unitVecX/unitVecZ) + xProd;
    double y = z * (unitVecY/unitVecZ) + yProd;
 
    if(xy == 0){return x;}
    if(xy == 1){return y;}
    return -999; 
}
 
 
»