Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-12-16 09:18:04 

0001 #ifndef InttSeedTracking_h
0002 #define InttSeedTracking_h
0003 
0004 #include <TROOT.h>
0005 #include <TChain.h>
0006 #include <TFile.h>
0007 
0008 #include <TMath.h>
0009 #include "Fit/Fitter.h"
0010 #include <Math/Functor.h>
0011 
0012 #include <TH1.h>
0013 #include <TH2.h>
0014 #include <TH3.h>
0015 #include <TF1.h>
0016 
0017 #include "TArc.h"
0018 
0019 #include "SPHTracKuma.h"
0020 
0021 #include <vector>
0022 
0023 class InttSeedTracking 
0024 {
0025 public :
0026     InttSeedTracking(std::vector<tracKuma>& tracks,\
0027         std::vector<hitStruct > vFMvtxHits,\
0028         std::vector<hitStruct > vSMvtxHits, std::vector<hitStruct > vTMvtxHits,\
0029         std::vector<hitStruct > vIInttHits, std::vector<hitStruct > vOInttHits,\
0030         std::vector<hitStruct > vEmcalHits,\
0031         std::vector<hitStruct > vIHCalHits, std::vector<hitStruct > vOHcalHits);
0032     InttSeedTracking(){};
0033     virtual ~InttSeedTracking();
0034 
0035     void HitMatching(std::vector<tracKuma>& tracks, std::vector<hitStruct > vFMvtxHits,\
0036     std::vector<hitStruct > vSMvtxHits, std::vector<hitStruct > vTMvtxHits,\
0037     std::vector<hitStruct > vIInttHits, std::vector<hitStruct > vOInttHits,\
0038     std::vector<hitStruct > vEmcalHits,\
0039     std::vector<hitStruct > vIHCalHits, std::vector<hitStruct > vOHcalHits);
0040 
0041     Double_t AddHCalE(Double_t emcalPhi, Double_t emcalE,\
0042         std::vector<hitStruct > vIHCalHits, std::vector<hitStruct > vOHcalHits);
0043 
0044     void AddMvtxHits(tracKuma& trk, std::vector<hitStruct> vFMvtxHits,\
0045         std::vector<hitStruct> vSMvtxHits, std::vector<hitStruct> vTMvtxHits, Double_t dRThre);
0046 
0047     void TrackPropertiesEstimation(tracKuma& trk, std::vector<hitStruct > vFMvtxHits,\
0048     std::vector<hitStruct > vSMvtxHits, std::vector<hitStruct > vTMvtxHits);
0049     Double_t TrackPtEstimation(tracKuma& trk, std::vector<hitStruct > vFMvtxHits,\
0050     std::vector<hitStruct > vSMvtxHits, std::vector<hitStruct > vTMvtxHits);
0051 
0052     void RoughEstiSagittaCenter3Point(Double_t& sagittaR, \
0053         Double_t& sagittaCenterX, Double_t& sagittaCenterY, Double_t HitsXY[3][2]);
0054 
0055     Double_t AccuratePtEstimation(Double_t sagittaR, Double_t centerX, Double_t centerY, tracKuma& trk, std::vector<hitStruct > vFMvtxHits,\
0056         std::vector<hitStruct > vSMvtxHits, std::vector<hitStruct > vTMvtxHits);
0057 
0058     Int_t TempINTTIOMatching(Double_t oINTTPhi, std::vector<hitStruct > vIInttHits); //????
0059     Double_t TempCalcdPhidR(Int_t iInttID, Int_t oInttID, \
0060         std::vector<hitStruct> vIInttHits, std::vector<hitStruct> vOInttHits); //????
0061     Int_t TempInttCalMatch(Int_t iInttID, Int_t oInttID,\
0062         std::vector<hitStruct > vIInttHits, std::vector<hitStruct > vOInttHits,\
0063         std::vector<hitStruct > vEmcalHits); //????
0064 
0065     void SagittaRByCircleFit(Double_t& cernterX, Double_t& centerY, Double_t& sagittaR,
0066         std::vector<Double_t > r, std::vector<Double_t > phi , Double_t oInttPhi, Double_t emcalPhi);
0067 
0068     // bool MatchingMvtxHits(Int_t& mvtxMatchId, Double_t& mvtxClusX, Double_t& mvtxClusY,\            1227
0069     //     Double_t sagittaR, Double_t centerX, Double_t centerY, std::vector<hitStruct > vMvtxHits);
0070     bool FindHitsOnCircle(Int_t& hitMatchId, Double_t& hitX, Double_t& hitY, Double_t sagittaR, \
0071         Double_t centerX, Double_t centerY, std::vector<hitStruct > vHits, Double_t dRThre);
0072 
0073     Double_t EstimateRecoTheta(tracKuma trk, Int_t type);
0074     // void EstiVertex(Double_t& vX, Double_t& vY, Double_t sagittaR, Double_t centerX, Double_t centerY);
0075     void EstiVertex(tracKuma trk);
0076     
0077     void ParticleIdentify(tracKuma trk);
0078 
0079     bool ReturnHitsRPhiVect(std::vector<Double_t >& hitR, std::vector<Double_t >& hitPhi,\
0080     std::vector<Int_t > subDetSet, tracKuma trk);
0081 
0082     void CrossLineCircle(Double_t& xv, Double_t& yv, Double_t xc, Double_t yc, Double_t R);
0083 
0084     // == s == Simple math equations  ##############################################
0085 
0086     inline Double_t CalcSagittaPt(Double_t sagittaR)
0087     {
0088         Double_t pT = 0.3*1.4*sagittaR*0.01; //0.3 * 1.4 T *(R [cm -> m])
0089         return pT;
0090     }
0091 
0092     // From vtx, mvtx(f,s,t), intt(i,o), emcal. there are some method to describe dPhi: 
0093     // dPhiOInttEmcalEsti, dPhiVtxIInttEmcalEsti, dPhiVtxOInttEmcalEsti, dPhiVtxInttEmcalEsti, 
0094     // dPhiMvtxIInttEmcalEsti, dPhiMvtxOInttEmcalEsti, dPhiMvtxInttEmcalEsti
0095     inline Double_t dPhiOInttEmcalEsti(tracKuma trk)
0096     {
0097         Double_t xIIntt = trk.getHitR(4)*std::cos(trk.getHitPhi(4));
0098         Double_t yIIntt = trk.getHitR(4)*std::sin(trk.getHitPhi(4));
0099         Double_t xOIntt = trk.getHitR(5)*std::cos(trk.getHitPhi(5));
0100         Double_t yOIntt = trk.getHitR(5)*std::sin(trk.getHitPhi(5));
0101         Double_t xEmcal = trk.getHitR(6)*std::cos(trk.getHitPhi(6));
0102         Double_t yEmcal = trk.getHitR(6)*std::sin(trk.getHitPhi(6));
0103 
0104         Double_t phiIInttOIntt = std::atan((yOIntt-yIIntt)/(xOIntt-xIIntt));
0105         Double_t phiOInttEmcal = std::atan((yEmcal-yOIntt)/(xEmcal-xOIntt));
0106 
0107         Double_t dPhiOInttEmcal = phiOInttEmcal - phiIInttOIntt;
0108 
0109         return dPhiOInttEmcal;
0110     }
0111     inline Double_t dPhiVtxIInttEmcalEsti(tracKuma trk)
0112     {
0113         Double_t xVtx = trk.getHitR(0)*std::cos(trk.getHitPhi(0));
0114         Double_t yVtx = trk.getHitR(0)*std::sin(trk.getHitPhi(0));
0115         Double_t xIIntt = trk.getHitR(4)*std::cos(trk.getHitPhi(4));
0116         Double_t yIIntt = trk.getHitR(4)*std::sin(trk.getHitPhi(4));
0117         Double_t xEmcal = trk.getHitR(6)*std::cos(trk.getHitPhi(6));
0118         Double_t yEmcal = trk.getHitR(6)*std::sin(trk.getHitPhi(6));
0119         Double_t phiVtxIIntt = std::atan((yIIntt-yVtx)/(xIIntt-xVtx));
0120         Double_t phiIInttEmcal = std::atan((yEmcal-yIIntt)/(xEmcal-xIIntt));
0121         Double_t dPhiIInttEmcal = phiIInttEmcal - phiVtxIIntt;
0122 
0123         return dPhiIInttEmcal;
0124     }
0125     inline Double_t dPhiVtxOInttEmcalEsti(tracKuma trk)
0126     {
0127         Double_t xVtx = trk.getHitR(0)*std::cos(trk.getHitPhi(0));
0128         Double_t yVtx = trk.getHitR(0)*std::sin(trk.getHitPhi(0));
0129         Double_t xOIntt = trk.getHitR(5)*std::cos(trk.getHitPhi(5));
0130         Double_t yOIntt = trk.getHitR(5)*std::sin(trk.getHitPhi(5));
0131         Double_t xEmcal = trk.getHitR(6)*std::cos(trk.getHitPhi(6));
0132         Double_t yEmcal = trk.getHitR(6)*std::sin(trk.getHitPhi(6));
0133         Double_t phiVtxOIntt = std::atan((yOIntt-yVtx)/(xOIntt-xVtx));
0134         Double_t phiOInttEmcal = std::atan((yEmcal-yOIntt)/(xEmcal-xOIntt));
0135         Double_t dPhiOInttEmcal = phiOInttEmcal - phiVtxOIntt;
0136 
0137         return dPhiOInttEmcal;
0138     }
0139     inline Double_t dPhiVtxInttEmcalEsti(tracKuma trk)
0140     {
0141         Double_t xVtx = trk.getHitR(0)*std::cos(trk.getHitPhi(0));
0142         Double_t yVtx = trk.getHitR(0)*std::sin(trk.getHitPhi(0));
0143         Double_t xIIntt = trk.getHitR(4)*std::cos(trk.getHitPhi(4));
0144         Double_t yIIntt = trk.getHitR(4)*std::sin(trk.getHitPhi(4));
0145         Double_t xOIntt = trk.getHitR(5)*std::cos(trk.getHitPhi(5));
0146         Double_t yOIntt = trk.getHitR(5)*std::sin(trk.getHitPhi(5));
0147         Double_t xEmcal = trk.getHitR(6)*std::cos(trk.getHitPhi(6));
0148         Double_t yEmcal = trk.getHitR(6)*std::sin(trk.getHitPhi(6));
0149         Double_t phiVtxIIntt = std::atan((yIIntt-yVtx)/(xIIntt-xVtx));
0150         Double_t phiVtxOIntt = std::atan((yOIntt-yVtx)/(xOIntt-xVtx));
0151         Double_t phiVtxIntt = (phiVtxOIntt - phiVtxIIntt)/2; 
0152         Double_t phiInttEmcal = std::atan((yEmcal-yOIntt)/(xEmcal-xOIntt));
0153         Double_t dPhiInttEmcal = phiInttEmcal - phiVtxIntt;
0154 
0155         return dPhiInttEmcal;
0156     }
0157     inline Double_t dPhiMvtxIInttEmcalEsti(tracKuma trk)
0158     {
0159         Int_t numOfMvtx = 0;
0160         Double_t xMvtx = 0.;
0161         Double_t yMvtx = 0.;
0162         for(Int_t i = 1; i < 4; i++)
0163         {
0164             if(trk.getHitIs(i))
0165             {
0166                 xMvtx += trk.getHitR(i)*std::cos(trk.getHitPhi(i));
0167                 yMvtx += trk.getHitR(i)*std::sin(trk.getHitPhi(i));
0168                 numOfMvtx++;
0169             }
0170         }
0171         xMvtx /= numOfMvtx;
0172         yMvtx /= numOfMvtx;
0173         Double_t xIIntt = trk.getHitR(4)*std::cos(trk.getHitPhi(4));
0174         Double_t yIIntt = trk.getHitR(4)*std::sin(trk.getHitPhi(4));
0175         Double_t xEmcal = trk.getHitR(6)*std::cos(trk.getHitPhi(6));
0176         Double_t yEmcal = trk.getHitR(6)*std::sin(trk.getHitPhi(6));
0177         Double_t phiMvtxIntt = std::atan((yIIntt-yMvtx)/(xIIntt-xMvtx));
0178         Double_t phiInttEmcal = std::atan((yEmcal-yIIntt)/(xEmcal-xIIntt));
0179         Double_t dPhiIInttEmcal = phiInttEmcal - phiMvtxIntt;
0180 
0181         return dPhiIInttEmcal;
0182     }
0183     inline Double_t dPhiMvtxOInttEmcalEsti(tracKuma trk)
0184     {
0185         Int_t numOfMvtx = 0;
0186         Double_t xMvtx = 0.;
0187         Double_t yMvtx = 0.;
0188         for(Int_t i = 1; i < 4; i++)
0189         {
0190             if(trk.getHitIs(i))
0191             {
0192                 xMvtx += trk.getHitR(i)*std::cos(trk.getHitPhi(i));
0193                 yMvtx += trk.getHitR(i)*std::sin(trk.getHitPhi(i));
0194                 numOfMvtx++;
0195             }
0196         }
0197         xMvtx /= numOfMvtx;
0198         yMvtx /= numOfMvtx;
0199         Double_t xOIntt = trk.getHitR(5)*std::cos(trk.getHitPhi(5));
0200         Double_t yOIntt = trk.getHitR(5)*std::sin(trk.getHitPhi(5));
0201         Double_t xEmcal = trk.getHitR(6)*std::cos(trk.getHitPhi(6));
0202         Double_t yEmcal = trk.getHitR(6)*std::sin(trk.getHitPhi(6));
0203         Double_t phiMvtxOIntt = std::atan((yOIntt-yMvtx)/(xOIntt-xMvtx));
0204         Double_t phiOInttEmcal = std::atan((yEmcal-yOIntt)/(xEmcal-xOIntt));
0205         Double_t dPhiOInttEmcal = phiOInttEmcal - phiMvtxOIntt;
0206 
0207         return dPhiOInttEmcal;
0208     }
0209     inline Double_t dPhiMvtxInttEmcalEsti(tracKuma trk)
0210     {
0211         Int_t numOfMvtx = 0;
0212         Double_t xMvtx = 0.;
0213         Double_t yMvtx = 0.;
0214         for(Int_t i = 1; i < 4; i++)
0215         {
0216             if(trk.getHitIs(i))
0217             {
0218                 xMvtx += trk.getHitR(i)*std::cos(trk.getHitPhi(i));
0219                 yMvtx += trk.getHitR(i)*std::sin(trk.getHitPhi(i));
0220                 numOfMvtx++;
0221             }
0222         }
0223         xMvtx /= numOfMvtx;
0224         yMvtx /= numOfMvtx;
0225         Double_t xIIntt = trk.getHitR(4)*std::cos(trk.getHitPhi(4));
0226         Double_t yIIntt = trk.getHitR(4)*std::sin(trk.getHitPhi(4));
0227         Double_t xOIntt = trk.getHitR(5)*std::cos(trk.getHitPhi(5));
0228         Double_t yOIntt = trk.getHitR(5)*std::sin(trk.getHitPhi(5));
0229         Double_t xIntt = (xIIntt + xOIntt)/2;
0230         Double_t yIntt = (yIIntt + yOIntt)/2;
0231         Double_t xEmcal = trk.getHitR(6)*std::cos(trk.getHitPhi(6));
0232         Double_t yEmcal = trk.getHitR(6)*std::sin(trk.getHitPhi(6));
0233         Double_t phiMvtxIntt = std::atan((yIntt-yMvtx)/(xIntt-xMvtx));
0234         Double_t phiInttEmcal = std::atan((yEmcal-yIntt)/(xEmcal-xIntt));
0235         Double_t dPhiInttEmcal = phiInttEmcal - phiMvtxIntt;
0236 
0237         return dPhiInttEmcal;
0238     }
0239 
0240     // fitfunction from dPhi to pt
0241     inline Double_t FitFunctionPt(Double_t dPhi)
0242     {
0243         Double_t pT = 0.0249291 + 0.232799/dPhi -0.000752825/(dPhi*dPhi);
0244         return pT;
0245     }
0246     inline Double_t FitFunctionPt_VtxIInttEmcal(Double_t dPhi)
0247     {
0248         Double_t pT = -0.075123 + 0.241657 /dPhi + 7.81643e-05/(dPhi*dPhi);
0249         return pT;
0250     }
0251     inline Double_t FitFunctionPt_VtxOInttEmcal(Double_t dPhi)
0252     {
0253         Double_t pT = 0.0592938 + 0.210859/dPhi - 0.00157569/(dPhi*dPhi);
0254         return pT;
0255     }
0256     inline Double_t FitFunctionPt_VtxInttEmcal(Double_t dPhi){
0257         Double_t pT = 0.089305 + 0.202456/dPhi + 0.0018411/(dPhi*dPhi);
0258         return pT;
0259     }
0260     inline Double_t FitFunctionPt_MVtxIInttEmcal(Double_t dPhi)
0261     {
0262         Double_t pT = -0.125412 + 0.293345/dPhi - 0.00233657/(dPhi*dPhi);
0263         return pT;
0264     }
0265     inline Double_t FitFunctionPt_MVtxOInttEmcal(Double_t dPhi)
0266     {
0267         Double_t pT = -0.0708943 + 0.269396/dPhi - 0.00158606/(dPhi*dPhi);
0268         return pT;
0269     }
0270     inline Double_t FitFunctionPt_MVtxInttEmcal(Double_t dPhi)
0271     {
0272         Double_t pT = -0.0645267 + 0.268115/dPhi -0.00156531/(dPhi*dPhi);
0273         return pT;
0274     }
0275 
0276 
0277     inline void Set3PointsXY(Double_t (&HitXY)[3][2], tracKuma trk, Int_t type)
0278     {
0279         Double_t iInttR = trk.getHitR(4);
0280         Double_t iInttPhi = trk.getHitPhi(4);
0281         Double_t oInttR = trk.getHitR(5);
0282         Double_t oInttPhi = trk.getHitPhi(5);
0283         Double_t emcalR = trk.getHitR(6);
0284         Double_t emcalPhi = trk.getHitPhi(6);
0285         if(type==0)
0286         {
0287             HitXY[0][0] = 0.;
0288             HitXY[0][1] = 0.;
0289             HitXY[1][0] = trk.getHitR(4)*cos(trk.getHitPhi(4));
0290             HitXY[1][1] = trk.getHitR(4)*sin(trk.getHitPhi(4));
0291             HitXY[2][0] = trk.getHitR(5)*cos(trk.getHitPhi(5));
0292             HitXY[2][1] = trk.getHitR(5)*sin(trk.getHitPhi(5));
0293         }
0294         else if(type==1)
0295         {
0296             HitXY[0][0] = trk.getHitR(4)*cos(trk.getHitPhi(4));
0297             HitXY[0][1] = trk.getHitR(4)*sin(trk.getHitPhi(4));
0298             HitXY[1][0] = trk.getHitR(5)*cos(trk.getHitPhi(5));
0299             HitXY[1][1] = trk.getHitR(5)*sin(trk.getHitPhi(5));
0300             HitXY[2][0] = trk.getHitR(6)*cos(trk.getHitPhi(6));
0301             HitXY[2][1] = trk.getHitR(6)*sin(trk.getHitPhi(6));
0302         }
0303         else if(type==2)
0304         {
0305             HitXY[0][0] = 0.;
0306             HitXY[0][1] = 0.;
0307             HitXY[1][0] = trk.getHitR(5)*cos(trk.getHitPhi(5));
0308             HitXY[1][1] = trk.getHitR(5)*sin(trk.getHitPhi(5));
0309             HitXY[2][0] = trk.getHitR(6)*cos(trk.getHitPhi(6));
0310             HitXY[2][1] = trk.getHitR(6)*sin(trk.getHitPhi(6));
0311         }
0312     }
0313 
0314     inline void CalcPerpBis(Double_t& slope, Double_t& sector, 
0315         Double_t HitXY1[2], Double_t HitXY2[2])
0316     {
0317         Double_t tempSlope = (HitXY2[1] - HitXY1[1])/(HitXY2[0] - HitXY1[0]);
0318         Double_t tempSector = HitXY1[1] - tempSlope*HitXY1[0];
0319 
0320         Double_t centerX = (HitXY2[0] + HitXY1[0])/2;
0321         Double_t centerY = (HitXY2[1] + HitXY1[1])/2;
0322 
0323         slope = -1/tempSlope;
0324         sector = centerY - slope*centerX;
0325 
0326         return;
0327     }
0328 
0329     inline void CalcLineEq(Double_t& slope, Double_t& section,\
0330         Double_t x1, Double_t y1, Double_t x2, Double_t y2)
0331     {
0332         slope = (y2 - y1)/(x2 - x1);
0333         section = y1 - slope*x1;
0334     }
0335 
0336     inline void CalcCrossPoint(Double_t& crossX, Double_t& crossY,\
0337         Double_t origSlope, Double_t origSection, Double_t nextPointX, Double_t nextPointY)
0338     {
0339         Double_t invSlope = - (1/origSlope);      
0340         Double_t invSection = nextPointY - invSlope*nextPointX;
0341         crossX = - (invSection - origSection)/(invSlope - origSlope);
0342         crossY = origSlope*crossX + origSection;
0343     }
0344 
0345     inline Double_t CalcSection(Double_t x, Double_t y, Double_t slope)
0346     {
0347         Double_t section = y - slope*x;
0348         return section;
0349     }
0350   
0351     // what is CrossCircleCircle?
0352     inline Double_t CrossCircleCircle(Double_t& xcal, Double_t& ycal,\
0353         Double_t xc, Double_t yc, Double_t R, Double_t phiIntt)
0354     {
0355         Double_t slope = -1*(xc/yc);
0356         Double_t section = (xc*xc+yc*yc-R*R)/(2*yc);
0357         Double_t coeff1 = (slope*section)/(slope*slope + 1);
0358         Double_t coeff2 = (section*section - m_ECalR*m_ECalR)/(slope*slope + 1);
0359         Double_t x1 = -1 * coeff1 + std::sqrt(coeff1*coeff1 - coeff2);
0360         Double_t x2 = -1 * coeff1 - std::sqrt(coeff1*coeff1 - coeff2);
0361         Double_t y1 = slope*x1 + section;
0362         Double_t y2 = slope*x2 + section;
0363         Double_t phi1 = std::atan(y1/x1);    
0364         if((phi1 < 0)&&(x1 < 0)) phi1 += TMath::Pi();
0365         else if((phi1 > 0)&&(x1 < 0)) phi1 -= TMath::Pi();
0366         
0367         Double_t phi2 = std::atan(y2/x2);        
0368         if((phi2 < 0)&&(x2 < 0)) phi2 += TMath::Pi();
0369         else if((phi2 > 0)&&(x2 < 0)) phi2 -= TMath::Pi();
0370         
0371         Double_t calPhi = 0.;
0372         if((std::abs(phi1-phiIntt)) < (std::abs(phi2-phiIntt)))
0373         {
0374             xcal = x1;
0375             ycal = y1;
0376             calPhi = phi1;
0377         }
0378         else
0379         {
0380             xcal = x2;
0381             ycal = y2;
0382             calPhi = phi2;
0383         }
0384         return calPhi;
0385    }
0386 
0387     // == e == Simple math equations  ##############################################
0388 
0389 private:
0390     Double_t m_ECalR = 93.5; //surface
0391     // const Double_t m_ECalR = 106.; // center 106
0392     Double_t m_iHCalR = 127.502;
0393     Double_t m_oHCalR = 225.87;
0394 
0395     Double_t m_InttMatchPhiMin = -0.05;
0396     Double_t m_InttMatchPhiMax = 0.05;
0397 
0398     Double_t m_EcalEThre = 0.1; // GeV
0399     Double_t m_iHcalEThre = 5; // GeV
0400     Double_t m_oHcalEThre = 5; // GeV
0401 
0402     Double_t m_EcalSearchPhiMin = 0.;
0403     Double_t m_EcalSearchPhiMax = 0.;
0404 
0405     Double_t m_EOverPElectron = 0.2; // CheckumaDayo!! not be optimize
0406 
0407     std::vector<hitStruct> m_fMvtxHits;
0408     std::vector<hitStruct> m_sMvtxHits;
0409     std::vector<hitStruct> m_tMvtxHits;
0410 
0411     std::vector<hitStruct> m_iInttHits;
0412     std::vector<hitStruct> m_oInttHits;
0413 
0414     std::vector<hitStruct> m_emcalHits;
0415     std::vector<hitStruct> m_iHCalHits;
0416     std::vector<hitStruct> m_oHCalHits;
0417 };
0418 
0419 #endif
0420 
0421 
0422 
0423 
0424 
0425 
0426 
0427
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!