Back to home page

sPhenix code displayed by LXR
 
 

    

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

0001 // InttSeedTracking
0002 //   ├── InttSeedTracking::InttSeedTracking (构造函数)
0003 //   │     ├── HitMatching (未使用，代码注释)
0004 //   │     ├── RecoTracksInttSeed2 (主重建流程)
0005 //   │     └── TrackPropertiesEstimation2 (轨迹属性估算)
0006 //   ├── InttSeedTracking::~InttSeedTracking (析构函数)
0007 //   │     ├── 清理数据
0008 //   │     └── 收缩内存
0009 //   ├── HitMatching (轨迹匹配流程)
0010 //   ├── RecoTracksInttSeed2
0011 //   │     ├── 遍历 oINTT 数据并进行匹配
0012 //   │     └── 遍历 iINTT 数据并进行匹配
0013 //   ├── InttSeedMatching (匹配算法核心)
0014 //   │     ├── Phi 和 Theta 范围检查
0015 //   │     ├── dR 和 Chi² 优化
0016 //   │     └── 校验轨迹有效性
0017 //   ├── TrackPropertiesEstimation2 (轨迹属性估算)
0018 //   │     ├── 估算 pT (CalcSagittaPt)
0019 //   │     ├── 重建角度 (EstimateRecoTheta)
0020 //   │     ├── 顶点估算 (EstiVertex)
0021 //   │     └── 粒子鉴别 (ParticleIdentify)
0022 //   ├── 其他辅助功能
0023 //   │     ├── TempINTTIOMatching (内外层匹配)
0024 //   │     ├── TempCalcdPhidR (计算 ΔPhi/ΔR)
0025 //   │     ├── AddMvtxHits (添加 MVTX 命中点)
0026 //   │     ├── RoughEstiSagittaCenter3Point (粗略圆拟合)
0027 //   │     ├── AccuratePtEstimation (精确动量估算)
0028 //   │     ├── ReturnHitsRPhiVect (命中点极坐标返回)
0029 //   │     └── ParticleIdentify (粒子属性)
0030 
0031 #ifndef InttSeedTracking_cxx
0032 #define InttSeedTracking_cxx
0033 
0034 
0035 #include <TStyle.h>
0036 #include <TCanvas.h>
0037 
0038 #include "SPHTracKuma.h"
0039 
0040 InttSeedTracking::InttSeedTracking(std::vector<tracKuma>& tracks,\
0041    std::vector<hitStruct > vFMvtxHits,\
0042    std::vector<hitStruct > vSMvtxHits, std::vector<hitStruct > vTMvtxHits,\
0043    std::vector<hitStruct > vIInttHits, std::vector<hitStruct > vOInttHits,\
0044    std::vector<hitStruct > vEmcalHits,\
0045    std::vector<hitStruct > vIHCalHits, std::vector<hitStruct > vOHCalHits)
0046 {
0047    // reco way1
0048    // HitMatching(tracks, vFMvtxHits, vSMvtxHits, vTMvtxHits, vIInttHits, vOInttHits,\
0049    //    vEmcalHits, vIHCalHits, vOHCalHits);
0050    // for(Int_t iTrk = 0; iTrk < tracks.size(); iTrk++){
0051    //    TrackPropertiesEstimation(tracks.at(iTrk), vFMvtxHits, vSMvtxHits, vTMvtxHits);
0052    // }
0053 
0054    // reco way2
0055    RecoTracksInttSeed2(tracks, vFMvtxHits, vSMvtxHits, vTMvtxHits,\
0056       vIInttHits, vOInttHits, vEmcalHits, vIHCalHits, vOHCalHits);
0057    for(Int_t iTrk = 0; iTrk < tracks.size(); iTrk++){
0058       TrackPropertiesEstimation2(tracks.at(iTrk));
0059    }
0060 
0061    // Vertex estimation using all tracks
0062    // VertexFinder(); // 
0063 
0064    // Re-finding tracks for single INTT hits or only mvtx using estimated vertex
0065    // ReFindTrack(); // 
0066 
0067 }
0068 
0069 InttSeedTracking::~InttSeedTracking()
0070 {
0071    m_fMvtxHits.clear();
0072    m_sMvtxHits.clear();
0073    m_tMvtxHits.clear();
0074    m_iInttHits.clear();
0075    m_oInttHits.clear();
0076    m_emcalHits.clear();
0077    m_iHCalHits.clear();
0078    m_oHCalHits.clear();
0079 
0080    m_fMvtxHits.shrink_to_fit();
0081    m_sMvtxHits.shrink_to_fit();
0082    m_tMvtxHits.shrink_to_fit();
0083    m_iInttHits.shrink_to_fit();
0084    m_oInttHits.shrink_to_fit();
0085    m_emcalHits.shrink_to_fit();
0086    m_iHCalHits.shrink_to_fit();
0087    m_oHCalHits.shrink_to_fit();
0088 }
0089 
0090 
0091 void InttSeedTracking::HitMatching(std::vector<tracKuma>& tracks,\
0092    std::vector<hitStruct > vFMvtxHits, std::vector<hitStruct > vSMvtxHits, std::vector<hitStruct > vTMvtxHits,\
0093    std::vector<hitStruct > vIInttHits, std::vector<hitStruct > vOInttHits,\
0094    std::vector<hitStruct > vEmcalHits,\
0095    std::vector<hitStruct > vIHCalHits, std::vector<hitStruct > vOHcalHits)
0096 {
0097     Int_t numOfInttClus = vOInttHits.size();
0098     if(numOfInttClus == 0) return;
0099 
0100     for(Int_t iOInttClus = 0; iOInttClus < numOfInttClus; iOInttClus++)
0101     {
0102         Double_t oInttPhi = vOInttHits.at(iOInttClus).phi;
0103 
0104         Int_t closetIINTTID = TempINTTIOMatching(oInttPhi, vIInttHits);
0105         if(closetIINTTID == 9999) continue;      
0106         Int_t matchEcalID = TempInttCalMatch(closetIINTTID, iOInttClus, vIInttHits, vOInttHits, vEmcalHits);
0107         if(matchEcalID == 9999) continue;
0108 
0109         tracKuma trk;
0110         Double_t tempTheta = 0.;
0111 
0112         trk.setHitIs(0, true);
0113         trk.setHitR(0, 0.);
0114         trk.setHitZ(0, 0.);
0115         trk.setHitPhi(0, 0.);
0116         trk.setHitTheta(0, 0.);
0117 
0118         trk.setHitIs(4, true);
0119         trk.setHitR(4, vIInttHits.at(closetIINTTID).r);
0120         trk.setHitZ(4, vIInttHits.at(closetIINTTID).z);
0121         trk.setHitPhi(4, vIInttHits.at(closetIINTTID).phi);
0122         tempTheta = 2*atan(std::exp(-vIInttHits.at(closetIINTTID).eta));
0123         trk.setHitTheta(4, tempTheta);
0124 
0125         trk.setHitIs(5, true);
0126         trk.setHitR(5, vOInttHits.at(iOInttClus).r);
0127         trk.setHitZ(5, vOInttHits.at(iOInttClus).z);
0128         trk.setHitPhi(5, vOInttHits.at(iOInttClus).phi);
0129         tempTheta = 2*atan(std::exp(-vOInttHits.at(iOInttClus).eta));
0130         trk.setHitTheta(5, tempTheta);
0131 
0132         trk.setHitIs(6, true);
0133         trk.setHitR(6, vEmcalHits.at(matchEcalID).r);
0134         trk.setHitZ(6, vEmcalHits.at(matchEcalID).z);
0135         trk.setHitPhi(6, vEmcalHits.at(matchEcalID).phi);
0136         tempTheta = 2*atan(std::exp(-vEmcalHits.at(matchEcalID).eta));
0137         trk.setHitTheta(6, tempTheta);
0138 
0139         Double_t calE = vEmcalHits.at(matchEcalID).energy;
0140         calE = AddHCalE(vEmcalHits.at(matchEcalID).phi, calE, vIHCalHits, vOHcalHits); // kuma???
0141         trk.setTrackE(calE);
0142 
0143         tracks.push_back(trk);
0144         // m_tracks.push_back(trk);
0145     }
0146 }
0147 
0148 // == s == new tracking algorithm =====================
0149 bool InttSeedTracking::RecoTracksInttSeed2(std::vector<tracKuma>& tracks,\
0150    std::vector<hitStruct > vFMvtxHits,\
0151    std::vector<hitStruct > vSMvtxHits, std::vector<hitStruct > vTMvtxHits,\
0152    std::vector<hitStruct > vIInttHits, std::vector<hitStruct > vOInttHits,\
0153    std::vector<hitStruct > vEmcalHits,\
0154    std::vector<hitStruct > vIHCalHits, std::vector<hitStruct > vOHCalHits){
0155    
0156    std::vector<hitStruct > vRemainFMvtxHits;
0157    std::vector<hitStruct > vRemainSMvtxHits;
0158    std::vector<hitStruct > vRemainTMvtxHits;
0159    std::vector<hitStruct > vRemainIInttHits;
0160    std::copy(vFMvtxHits.begin(), vFMvtxHits.end(), std::back_inserter(vRemainFMvtxHits));
0161    std::copy(vSMvtxHits.begin(), vSMvtxHits.end(), std::back_inserter(vRemainSMvtxHits));
0162    std::copy(vTMvtxHits.begin(), vTMvtxHits.end(), std::back_inserter(vRemainTMvtxHits));
0163    std::copy(vIInttHits.begin(), vIInttHits.end(), std::back_inserter(vRemainIInttHits));
0164    
0165    for(Int_t iOInttClus = 0; iOInttClus < vOInttHits.size(); iOInttClus++){
0166       InttSeedMatching(tracks, 5, vOInttHits.at(iOInttClus),\
0167          vRemainFMvtxHits, vRemainSMvtxHits, vRemainTMvtxHits, vRemainIInttHits,\
0168             vEmcalHits, vIHCalHits, vOHCalHits);
0169    }// == e == oIntt Loop
0170    
0171    // ChecKumaDaYo!!!! I cannot use this function because it makes too many tracks.
0172    // This function is to recover tracks without oINTT 
0173    // I think to use it, you need to remove used MVTX hits for upper algorithm.
0174    // On the other hand, it will use long time...
0175    for(Int_t iIInttClus = 0; iIInttClus < vRemainIInttHits.size(); iIInttClus++){
0176       InttSeedMatching(tracks, 4, vRemainIInttHits.at(iIInttClus),\
0177          vRemainFMvtxHits, vRemainSMvtxHits, vRemainTMvtxHits,\
0178             vRemainIInttHits, vEmcalHits, vIHCalHits, vOHCalHits);
0179    }
0180 
0181    return true;
0182 }
0183 
0184 bool InttSeedTracking::InttSeedMatching(std::vector<tracKuma>& tracks, Int_t inttId,\
0185    hitStruct baseInttHit,\
0186    std::vector<hitStruct >& vFMvtxHits,\
0187    std::vector<hitStruct >& vSMvtxHits, std::vector<hitStruct >& vTMvtxHits,\
0188    std::vector<hitStruct >& vIInttHits, std::vector<hitStruct > vEmcalHits,\
0189    std::vector<hitStruct > vIHCalHits, std::vector<hitStruct > vOHCalHits)
0190 {
0191     Double_t inttTheta = 2*atan(std::exp(-baseInttHit.eta));
0192     Double_t inttPhi = baseInttHit.phi;
0193     Double_t inttR = baseInttHit.r;
0194     
0195     std::vector<tracKuma > v_tempTracks;
0196     tracKuma bestTrk;
0197     Int_t matchFMvtxId = 99999;
0198     Int_t matchSMvtxId = 99999;
0199     Int_t matchTMvtxId = 99999;
0200     Int_t matchIInttId = 99999;
0201     Double_t smallestChi = 9999.;
0202     Int_t matchiECalID = 9999;
0203     for(Int_t iECalT = 0; iECalT < vEmcalHits.size(); iECalT++)
0204     {
0205         Double_t ecalE = vEmcalHits.at(iECalT).energy;
0206         if(ecalE < m_EcalEThre) continue; // 能量阈值的筛选，太低能量的hit不要
0207 
0208         Double_t ecalTheta = 2*atan(std::exp(-vEmcalHits.at(iECalT).eta));
0209         // ChecKumaDaYo!!! you need to optimize the search range
0210         if((inttTheta - TMath::Pi()/10 < ecalTheta)&&(ecalTheta < inttTheta + TMath::Pi()/10)) // emcal eta match
0211         {
0212             Int_t minChiSqrtEMCalTrkId = 9999;
0213             Double_t minChiSqrt = 9999.;
0214             Double_t ecalPhi = vEmcalHits.at(iECalT).phi;
0215             // ChecKumaDaYo!!! you need to optimize the search range
0216             if((inttPhi - TMath::Pi()/10 < ecalPhi)&&(ecalPhi < inttPhi + TMath::Pi()/10)) // emcal phi match
0217             {
0218                 tracKuma tempTrack;
0219 
0220                 tempTrack.setHitIs(0, true); // vertex set
0221                 tempTrack.setHitR(0, 0.);
0222                 tempTrack.setHitZ(0, 0.);
0223                 tempTrack.setHitPhi(0, 0.);
0224                 tempTrack.setHitTheta(0, 0.);
0225 
0226                 SetHitParaInTrack(tempTrack, inttId, baseInttHit); // set INTT parameters, inttid 4-iintt, 5-oIntt
0227                 SetHitParaInTrack(tempTrack, 6, vEmcalHits.at(iECalT)); // set ECal parameters
0228 
0229                 Double_t HitsXY[3][2];
0230                 Set3PointsXY(HitsXY, tempTrack, 2); // convert r, phi -> x, y for the three points.
0231                 Double_t sagittaR = 9999.;
0232                 Double_t cX = 0.;
0233                 Double_t cY = 0.;
0234                 // estimate the circle using three points roughly
0235                 RoughEstiSagittaCenter3Point(sagittaR, cX, cY, HitsXY);
0236                 tempTrack.setTrackSagR(sagittaR);
0237                 tempTrack.setTrackSagX(cX);
0238                 tempTrack.setTrackSagY(cY);
0239 
0240                 // 先match emcal，再match iintt吗？
0241                 Int_t iInttMatchId = 9999;
0242                 if(inttId == 5)
0243                 {
0244                     Double_t iInttClusX = 9999.;
0245                     Double_t iInttClusY = 9999.;
0246                     Double_t iInttDRThre = 1.0;
0247                     // searching for the best iINTT hit satisfying dR threshold
0248                     // ChecKumaDaYo!!! you should optimize the dR threshold iInttDRThre
0249                     if(FindHitsOnCircle(iInttMatchId, iInttClusX, iInttClusY, sagittaR, cX, cY, vIInttHits, iInttDRThre))
0250                     {         
0251                         SetHitParaInTrack(tempTrack, 4, vIInttHits.at(iInttMatchId));
0252                     }      
0253                 }
0254 
0255                 // searching for the best MVTX hits satisfying dR threshold
0256                 // ChecKumaDaYo!!! you should optimize the dR threshold dRThre_Mvtx
0257                 Int_t seleFMvtxId = 99999;
0258                 Int_t seleSMvtxId = 99999;
0259                 Int_t seleTMvtxId = 99999;
0260                 Double_t dRThre_Mvtx = 5.;
0261                 AddMvtxHits(tempTrack, vFMvtxHits, vSMvtxHits, vTMvtxHits, dRThre_Mvtx, seleFMvtxId, seleSMvtxId, seleTMvtxId);
0262                 Double_t dRChi = EstiChiTrkOnCircle(tempTrack, cX, cY, sagittaR);
0263                 if(minChiSqrt > dRChi)
0264                 {
0265                     minChiSqrt = dRChi;
0266                     bestTrk = tempTrack;
0267                     matchFMvtxId = seleFMvtxId;
0268                     matchSMvtxId = seleSMvtxId;
0269                     matchTMvtxId = seleTMvtxId;
0270                     matchIInttId = iInttMatchId;
0271                 }
0272             }// == e == Ecal Phi range if
0273         }// == e == Ecal Theta range if
0274     }// == e == Ecal Loop
0275     
0276     // Judge the reconstructed track is available or not.
0277     // The requirements: 
0278     // (1) It has both iINTT and oINTT hits
0279     // (2) Single INTT and two MVTX hits
0280     if(!CheckTrkRequirements(bestTrk)) return false;
0281     RefindCalHit(bestTrk, vEmcalHits, vIHCalHits, vOHCalHits);
0282     tracks.push_back(bestTrk);
0283     if(matchFMvtxId!=99999) vFMvtxHits.erase(vFMvtxHits.begin() + matchFMvtxId);
0284     if(matchSMvtxId!=99999) vSMvtxHits.erase(vSMvtxHits.begin() + matchSMvtxId);
0285     if(matchTMvtxId!=99999) vTMvtxHits.erase(vTMvtxHits.begin() + matchTMvtxId);
0286     if((inttId == 5)&&(matchIInttId!=9999)) vIInttHits.erase(vIInttHits.begin() + matchIInttId);
0287 
0288     return true;
0289 }
0290 
0291 
0292 Double_t InttSeedTracking::EstiChiTrkOnCircle(tracKuma trk,\
0293    Double_t cX, Double_t cY, Double_t sagittaR){
0294    Double_t chi = 0.;
0295    Double_t numOfHits = 0.;
0296    for(Int_t iHit = 0; iHit < 7; iHit++){
0297       Double_t hitX = trk.getHitR(iHit)*cos(trk.getHitPhi(iHit));
0298       Double_t hitY = trk.getHitR(iHit)*sin(trk.getHitPhi(iHit));
0299 
0300       if(!trk.getHitIs(iHit)) continue;
0301       chi += std::sqrt((hitX - cX)*(hitX - cX) + (hitY - cY)*(hitY - cY)) - sagittaR;
0302       numOfHits++;
0303    }
0304    chi /= numOfHits;
0305    
0306    return chi;
0307 }
0308 
0309 // ointt + emcal + iintt/mvtx 3points
0310 bool InttSeedTracking::CheckTrkRequirements(tracKuma trk)
0311 {
0312     if(trk.getHitIs(4)) return true; // iintt
0313     Int_t numMvtx = 0;
0314     for(Int_t iMvtx = 1; iMvtx < 4; iMvtx++) if(trk.getHitIs(iMvtx)) numMvtx += 1; // mvtx
0315 
0316     if(numMvtx > 1) return true;
0317     else return false; 
0318 }
0319 
0320 void InttSeedTracking::RefindCalHit(tracKuma trk, std::vector<hitStruct> vEmcalHits, std::vector<hitStruct> vIHCalHits, std::vector<hitStruct> vOHcalHits)
0321 {
0322     std::vector<Int_t> subDetIds = {1, 2, 3, 4, 5};
0323     std::vector<Double_t> vHitR = {};
0324     std::vector<Double_t> vHitsPhi = {};
0325     if(!ReturnHitsRPhiVect(vHitR, vHitsPhi, subDetIds, trk)) return;
0326     Double_t cX = 0.;
0327     Double_t cY = 0.;
0328     Double_t sagittaR = 0.;
0329     Double_t tempHitIInttPhi = trk.getHitPhi(4);
0330     Double_t tempHitEmcalPhi = trk.getHitPhi(6);
0331     // ChecKumaDaYo!!! it does not work in only sPHENIX server?????
0332     // SagittaRByCircleFit(cX, cY, sagittaR, vHitR, vHitsPhi, trk.getHitPhi(4), trk.getHitPhi(6));
0333 
0334     Double_t targetCalX = 0.;
0335     Double_t targetCalY = 0.;
0336     CrossCircleCircle(targetCalX, targetCalY, cX, cY, sagittaR, trk.getHitPhi(4));
0337     Double_t calHighestE = 0.;
0338     Int_t highestECalId = 99999;
0339     for(Int_t iECalT = 0; iECalT < vEmcalHits.size(); iECalT++)
0340     {
0341         Double_t ecalE = vEmcalHits.at(iECalT).energy;
0342         if(ecalE < m_EcalEThre) continue;
0343         if(ecalE < calHighestE) continue;
0344 
0345         Double_t tempEcalPhi = vEmcalHits.at(iECalT).phi;
0346         Double_t targetECalPhi = std::tan(targetCalY/targetCalX);
0347         if((targetECalPhi < 0)&&(cX < 0)) targetECalPhi += TMath::Pi();
0348         else if((targetECalPhi > 0)&&(cX < 0)) targetECalPhi -= TMath::Pi();
0349         if((targetECalPhi - 0.05 < tempEcalPhi)&&(tempEcalPhi < targetECalPhi + 0.05))
0350         {
0351             calHighestE = ecalE;
0352             highestECalId = iECalT;
0353         }
0354     }
0355     
0356     if(highestECalId == 99999) return;
0357     trk.setHitIs(6, true);
0358     trk.setHitR(6, vEmcalHits.at(highestECalId).r);
0359     trk.setHitZ(6, vEmcalHits.at(highestECalId).z);
0360     trk.setHitPhi(6, vEmcalHits.at(highestECalId).phi);
0361     Double_t ecalTheta = 2*atan(std::exp(-vEmcalHits.at(highestECalId).eta));
0362     trk.setHitTheta(6, ecalTheta);
0363     
0364     // ChecKumaDaYo!!!
0365     CalESumAndCorrPosi(trk, vEmcalHits, vIHCalHits, vOHcalHits);
0366 
0367     // Double_t calE = vEmcalHits.at(highestECalId).energy;
0368     // // ChecKumaDaYo!!!
0369     // calE = AddHCalE(vEmcalHits.at(highestECalId).phi, calE, vIHCalHits, vOHcalHits);
0370     // trk.setTrackE(calE);
0371 
0372 }
0373 
0374 // 对emcalhit取加权平均(w-energy)，energy加上ihcal和ohcal
0375 void InttSeedTracking::CalESumAndCorrPosi(tracKuma trk, std::vector<hitStruct> vEmcalHits,\
0376     std::vector<hitStruct> vIHCalHits, std::vector<hitStruct> vOHCalHits)
0377 {
0378     Double_t refCalPhi = trk.getHitPhi(6);
0379     Double_t refCalTheta = trk.getHitTheta(6);
0380     Double_t TotEMCalE = 0.;
0381     Double_t ModifEMCalPhi = 0.;
0382     Double_t ModifEMCalTheta = 0.;
0383 
0384     // 加权平均算phi，theta
0385     for(Int_t iEmcal = 0; iEmcal < vEmcalHits.size(); iEmcal++)
0386     {
0387         Double_t hitTheta = 2*atan(std::exp(-vEmcalHits.at(iEmcal).eta));
0388         if((hitTheta < refCalTheta - TMath::Pi()/20)&&(refCalTheta + TMath::Pi()/20 < hitTheta)) continue;
0389         if((vEmcalHits.at(iEmcal).phi < refCalPhi - TMath::Pi()/20)&&(refCalPhi + TMath::Pi()/20 < vEmcalHits.at(iEmcal).phi)) continue;
0390         
0391         TotEMCalE += vEmcalHits.at(iEmcal).energy;
0392         ModifEMCalPhi += vEmcalHits.at(iEmcal).energy*vEmcalHits.at(iEmcal).phi;
0393         ModifEMCalTheta += vEmcalHits.at(iEmcal).energy*hitTheta;
0394     }
0395     ModifEMCalPhi /= TotEMCalE;
0396     ModifEMCalTheta /= TotEMCalE;
0397 
0398     refCalPhi = ModifEMCalPhi;
0399     for(Int_t iIHcal = 0; iIHcal < vIHCalHits.size(); iIHcal++)
0400     {
0401         Double_t hitTheta = 2*atan(std::exp(-vIHCalHits.at(iIHcal).eta));
0402         if((hitTheta < refCalTheta - TMath::Pi()/20)&&(refCalTheta + TMath::Pi()/20 < hitTheta)) continue;
0403         if((vIHCalHits.at(iIHcal).phi < refCalPhi - TMath::Pi()/20)&&(refCalPhi + TMath::Pi()/20 < vIHCalHits.at(iIHcal).phi)) continue;
0404         
0405         TotEMCalE += vIHCalHits.at(iIHcal).energy;
0406     }
0407 
0408     for(Int_t iOHcal = 0; iOHcal < vOHCalHits.size(); iOHcal++)
0409     {
0410         Double_t hitTheta = 2*atan(std::exp(-vOHCalHits.at(iOHcal).eta));
0411         if((hitTheta < refCalTheta - TMath::Pi()/20)&&(refCalTheta + TMath::Pi()/20 < hitTheta)) continue;
0412         if((vOHCalHits.at(iOHcal).phi < refCalPhi - TMath::Pi()/20)&&(refCalPhi + TMath::Pi()/20 < vIHCalHits.at(iOHcal).phi)) continue;
0413         TotEMCalE += vIHCalHits.at(iOHcal).energy;
0414     }
0415     trk.setHitPhi(6, ModifEMCalPhi);
0416     trk.setHitTheta(6, ModifEMCalTheta);
0417     trk.setTrackE(TotEMCalE);
0418 }
0419 
0420 void InttSeedTracking::TrackPropertiesEstimation2(tracKuma& trk){
0421    Double_t dPhiOInttEmcal = dPhiOInttEmcalEsti(trk);
0422    Double_t recoPt = CalcSagittaPt(trk.getTrackSagR());
0423    // Double_t recoPt = FitFunctionPt(dPhiOInttEmcal);
0424    trk.setTrackPt(recoPt);
0425    
0426    Double_t recoTheta = EstimateRecoTheta(trk, 1);
0427    trk.setTrackTheta(recoTheta);
0428 
0429    Double_t recoP = recoPt/sin(recoTheta);
0430    trk.setTrackP(recoP);
0431 
0432    EstiVertex(trk);
0433 
0434    ParticleIdentify(trk);
0435 }
0436 // == e == new tracking algorithm =====================
0437 
0438 
0439 Int_t InttSeedTracking::TempINTTIOMatching(Double_t oINTTPhi, std::vector<hitStruct > vIInttHits){
0440    Double_t minDeltaPhi = 9999.;
0441    Int_t closestIINTTCluID = 9999;
0442    Int_t numOfIInttClu = vIInttHits.size();
0443 
0444    if(numOfIInttClu == 0) return 9999;
0445    for(Int_t iINTTClu = 0; iINTTClu < numOfIInttClu; iINTTClu++){
0446       Double_t iInttPhi =  vIInttHits.at(iINTTClu).phi;
0447       Double_t dPhi = std::abs(iInttPhi - oINTTPhi);
0448       if(std::abs(dPhi) > TMath::Pi()) dPhi += (-1) * (dPhi/std::abs(dPhi)) * 2*TMath::Pi();
0449 
0450       // checkumaDAYO!!! need to optimize the matching range.
0451       if((minDeltaPhi > dPhi)&&(dPhi > m_InttMatchPhiMin)&&(dPhi < m_InttMatchPhiMax)){
0452          minDeltaPhi = dPhi;
0453          closestIINTTCluID = iINTTClu;
0454       }
0455    }
0456    return closestIINTTCluID;
0457 }
0458 
0459 
0460 Double_t InttSeedTracking::TempCalcdPhidR(Int_t iInttID, Int_t oInttID,\
0461    std::vector<hitStruct > vIInttHits, std::vector<hitStruct > vOInttHits){
0462    Double_t dPhidR = 0.;
0463    Double_t iInttPhi = vIInttHits.at(iInttID).phi;
0464    Double_t oInttPhi = vOInttHits.at(oInttID).phi;
0465 
0466    Double_t dPhi = oInttPhi - iInttPhi;
0467    if(std::abs(dPhi) > TMath::Pi()) dPhi += (-1) * (dPhi/std::abs(dPhi)) * 2*TMath::Pi();
0468    Double_t dR = vOInttHits.at(oInttID).r - vIInttHits.at(iInttID).r;
0469 
0470    dPhidR = dPhi/dR;
0471 
0472    return dPhidR;
0473 }
0474 
0475 
0476 Int_t InttSeedTracking::TempInttCalMatch(Int_t iInttID, Int_t oInttID,
0477    std::vector<hitStruct > vIInttHits, std::vector<hitStruct > vOInttHits,\
0478    std::vector<hitStruct > vEmcalHits){
0479    Int_t matchiECalID = 9999;
0480    Double_t matchiECalE = 0.;
0481    Int_t numOfECalT = vEmcalHits.size();
0482    for(Int_t iECalT = 0; iECalT < numOfECalT ;iECalT++){
0483       Double_t ecalE = vEmcalHits.at(iECalT).energy;
0484       // checkumaDAYO!!! need to optimize the calorimeter threshold energy.
0485       if(ecalE < m_EcalEThre) continue;
0486 
0487       Double_t dPhidR = TempCalcdPhidR(iInttID, oInttID, vIInttHits, vOInttHits);
0488       Double_t dRInttEmCal = m_ECalR - vOInttHits.at(oInttID).r;
0489       Double_t phiRangeMin = 0.;
0490       Double_t phiRangeMax = 0.;
0491 
0492       // checkumaDAYO!!! need to optimize the matching range.
0493       if(dPhidR < 0){
0494          phiRangeMin = vOInttHits.at(oInttID).phi + (dPhidR * dRInttEmCal - 0.177);
0495          phiRangeMax = vOInttHits.at(oInttID).phi + 0.087;
0496       }else{
0497          phiRangeMin = vOInttHits.at(oInttID).phi - 0.087;
0498          phiRangeMax = vOInttHits.at(oInttID).phi + (dPhidR * dRInttEmCal + 0.177);
0499       }
0500       
0501       Double_t emcalPhi = vEmcalHits.at(iECalT).phi;
0502       if((phiRangeMin < - TMath::Pi())&&(phiRangeMin > - TMath::Pi())&&(emcalPhi > 0)) emcalPhi -= 2*TMath::Pi();
0503       else if((phiRangeMin < TMath::Pi())&&(phiRangeMax > TMath::Pi())&&(emcalPhi < 0)) emcalPhi += 2*TMath::Pi();
0504 
0505       if((phiRangeMin < emcalPhi)&&(emcalPhi < phiRangeMax)){
0506          if(matchiECalE < ecalE){
0507             matchiECalID = iECalT;
0508             matchiECalE = ecalE;
0509          }
0510       }
0511       // std::cout << "ecalE = " << ecalE << std::endl;
0512    }
0513    // std::cout << "matchiECalE = " << matchiECalE << std::endl;
0514    return matchiECalID;
0515 
0516 }
0517  
0518 
0519 Double_t InttSeedTracking::AddHCalE(Double_t emcalPhi, Double_t emcalE,\
0520    std::vector<hitStruct > vIHCalHits, std::vector<hitStruct > vOHCalHits){
0521    Int_t matchIhcalId = 99999;
0522    Double_t closePhiIhcal = 99999;
0523    for(Int_t iIhcal = 0; iIhcal < vIHCalHits.size(); iIhcal++){
0524       if(closePhiIhcal > std::abs(vIHCalHits.at(iIhcal).phi - emcalPhi)){
0525          closePhiIhcal = std::abs(vIHCalHits.at(iIhcal).phi - emcalPhi);
0526          matchIhcalId = iIhcal;
0527       }
0528    }
0529    Int_t matchOhcalId = 99999;
0530    Double_t closePhiOhcal = 99999;
0531    for(Int_t iOhcal = 0; iOhcal < vOHCalHits.size(); iOhcal++){
0532       if(closePhiOhcal > std::abs(vOHCalHits.at(iOhcal).phi - emcalPhi)){
0533          closePhiOhcal = std::abs(vOHCalHits.at(iOhcal).phi - emcalPhi);
0534          matchOhcalId = iOhcal;
0535       }
0536    }
0537    
0538    Double_t iHcalE = 0.;
0539    Double_t oHcalE = 0.;
0540    if(matchIhcalId != 99999) iHcalE = vIHCalHits.at(matchIhcalId).energy;
0541    if(matchOhcalId != 99999) oHcalE = vOHCalHits.at(matchOhcalId).energy;
0542 
0543    Double_t calE = emcalE + iHcalE + oHcalE;
0544 
0545    return calE;
0546 }
0547 
0548 void InttSeedTracking::AddMvtxHits(tracKuma& trk, std::vector<hitStruct > vFMvtxHits,\
0549    std::vector<hitStruct > vSMvtxHits, std::vector<hitStruct > vTMvtxHits, Double_t dRThre,\
0550    Int_t& seleFMvtxId, Int_t& seleSMvtxId, Int_t& seleTMvtxId){
0551    Double_t sagiR = trk.getTrackSagR();
0552    Double_t sagiX = trk.getTrackSagX();
0553    Double_t sagiY = trk.getTrackSagY();
0554 
0555    for(Int_t iMvtxLay = 0; iMvtxLay < 3; iMvtxLay++){
0556       Int_t mvtxMatchId = 9999;
0557       Double_t mvtxClusX = 9999.;
0558       Double_t mvtxClusY = 9999.;
0559       if(iMvtxLay == 0){
0560          if(FindHitsOnCircle(mvtxMatchId, mvtxClusX, mvtxClusY, sagiR, sagiX, sagiY, \
0561             vFMvtxHits, dRThre)){
0562             SetHitParaInTrack(trk, 1, vFMvtxHits.at(mvtxMatchId));
0563             seleFMvtxId = mvtxMatchId;
0564          }
0565       }else if(iMvtxLay == 1){
0566          if(FindHitsOnCircle(mvtxMatchId, mvtxClusX, mvtxClusY, sagiR, sagiX, sagiY,\
0567             vSMvtxHits, dRThre)){
0568             SetHitParaInTrack(trk, 2, vSMvtxHits.at(mvtxMatchId));
0569             seleSMvtxId = mvtxMatchId;
0570          }
0571       }else if(iMvtxLay == 2){
0572          if(FindHitsOnCircle(mvtxMatchId, mvtxClusX, mvtxClusY, sagiR, sagiX, sagiY,\
0573             vTMvtxHits, dRThre)){
0574             SetHitParaInTrack(trk, 3, vTMvtxHits.at(mvtxMatchId));
0575             seleTMvtxId = mvtxMatchId;
0576          }
0577       }
0578    }
0579 }
0580 
0581 void InttSeedTracking::TrackPropertiesEstimation(tracKuma& trk, std::vector<hitStruct > vFMvtxHits,\
0582    std::vector<hitStruct > vSMvtxHits, std::vector<hitStruct > vTMvtxHits){
0583    Double_t recoPt = TrackPtEstimation(trk, vFMvtxHits, vSMvtxHits, vTMvtxHits);
0584    trk.setTrackPt(recoPt);
0585    
0586    Double_t recoTheta = EstimateRecoTheta(trk, 1);
0587    trk.setTrackTheta(recoTheta);
0588 
0589    Double_t recoP = recoPt/sin(recoTheta);
0590    trk.setTrackP(recoP);
0591    // Double_t recoE = trk.getTrackE();
0592    // Double_t EOverP = recoE/recoP;
0593 
0594    EstiVertex(trk);
0595 
0596    ParticleIdentify(trk);
0597 }
0598 
0599 Double_t InttSeedTracking::TrackPtEstimation(tracKuma& trk, std::vector<hitStruct > vFMvtxHits,\
0600    std::vector<hitStruct > vSMvtxHits, std::vector<hitStruct > vTMvtxHits){
0601    Double_t sagittaR = 9999.;
0602    Double_t sagittaCenterX = 0.;
0603    Double_t sagittaCenterY = 0.;
0604    Double_t HitsXY[3][2];
0605    Set3PointsXY(HitsXY, trk, 2);
0606    RoughEstiSagittaCenter3Point(sagittaR, sagittaCenterX, sagittaCenterY, HitsXY);
0607 
0608    Double_t trackPt = AccuratePtEstimation(sagittaR, sagittaCenterX, sagittaCenterY, trk, \
0609    vFMvtxHits, vSMvtxHits, vTMvtxHits);
0610 
0611    Double_t dPhiOInttEmcal = dPhiOInttEmcalEsti(trk);
0612    trackPt = FitFunctionPt(dPhiOInttEmcal);
0613 
0614    return trackPt;
0615 }
0616 
0617 
0618 void InttSeedTracking::RoughEstiSagittaCenter3Point(Double_t& sagittaR,\
0619    Double_t& sagittaCenterX, Double_t& sagittaCenterY, Double_t HitsXY[3][2]){
0620    Double_t iSlope = 0.;
0621    Double_t iSection = 0.;
0622    CalcPerpBis(iSlope, iSection, HitsXY[0], HitsXY[2]);
0623 
0624    Double_t oSlope = 0.;
0625    Double_t oSection = 0.;
0626    CalcPerpBis(oSlope, oSection, HitsXY[1], HitsXY[2]);
0627 
0628    sagittaCenterX = - (oSection - iSection)/(oSlope - iSlope);
0629    sagittaCenterY = iSlope * sagittaCenterX + iSection;
0630 
0631    sagittaR = std::sqrt((HitsXY[0][0] - sagittaCenterX)*(HitsXY[0][0] - sagittaCenterX) \
0632       + (HitsXY[0][1] - sagittaCenterY)*(HitsXY[0][1] - sagittaCenterY));
0633 
0634 }
0635 
0636 Double_t InttSeedTracking::AccuratePtEstimation(\
0637    Double_t sagittaR, Double_t centerX, Double_t centerY, tracKuma& trk,
0638    std::vector<hitStruct > vFMvtxHits,\
0639    std::vector<hitStruct > vSMvtxHits, std::vector<hitStruct > vTMvtxHits){
0640 
0641    Double_t sagittaVtxInttPt = 99999.;
0642    Double_t sagittaInttEmcalPt = 99999.;
0643    Double_t sagittaVtxInttEmcalPt = 99999.;
0644    Double_t sagittaMvtxInttEmcalPt = 99999.;
0645    Double_t sagittaVtxMvtxInttEmcalPt = 99999.;
0646    
0647    // checkumaDaYo!!!
0648    std::vector<Int_t > subDetIds_InttEmcal = {4, 5, 6};
0649    std::vector<Double_t > hitsR_InttEmcal = {};
0650    std::vector<Double_t > hitsPhi_InttEmcal = {};
0651    if(ReturnHitsRPhiVect(hitsR_InttEmcal, hitsPhi_InttEmcal, subDetIds_InttEmcal, trk)){
0652       Double_t tempHitOInttPhi = trk.getHitPhi(4);
0653       Double_t tempHitEmcalPhi = trk.getHitPhi(6);
0654       SagittaRByCircleFit(centerX, centerY, sagittaR, hitsR_InttEmcal, hitsPhi_InttEmcal, tempHitOInttPhi, tempHitEmcalPhi);
0655       
0656       sagittaMvtxInttEmcalPt = CalcSagittaPt(sagittaR);
0657    }
0658 
0659    Int_t seleFMvtxId = 99999;
0660    Int_t seleSMvtxId = 99999;
0661    Int_t seleTMvtxId = 99999;
0662    Double_t dRThre_Mvtx = 5.;
0663    AddMvtxHits(trk, vFMvtxHits, vSMvtxHits, vTMvtxHits, dRThre_Mvtx,\
0664       seleFMvtxId, seleSMvtxId, seleTMvtxId);
0665    
0666    std::vector<Int_t > subDetIds_MvtxInttEmcal = {1, 2, 3, 4, 5, 6};
0667    std::vector<Double_t > hitsR_MvtxInttEmcal = {};
0668    std::vector<Double_t > hitsPhi_MvtxInttEmcal = {};
0669    if(ReturnHitsRPhiVect(hitsR_MvtxInttEmcal, hitsPhi_MvtxInttEmcal, subDetIds_MvtxInttEmcal, trk)){
0670       Double_t tempHitOInttPhi = trk.getHitPhi(4);
0671       Double_t tempHitEmcalPhi = trk.getHitPhi(6);
0672       SagittaRByCircleFit(centerX, centerY, sagittaR, hitsR_MvtxInttEmcal, hitsPhi_MvtxInttEmcal,\
0673       tempHitOInttPhi, tempHitEmcalPhi);
0674       sagittaMvtxInttEmcalPt = CalcSagittaPt(sagittaR);
0675    }
0676 
0677    trk.setTrackSagR(sagittaR);
0678    trk.setTrackSagX(centerX);
0679    trk.setTrackSagY(centerY);
0680 
0681    return sagittaMvtxInttEmcalPt;
0682 }
0683 
0684 void InttSeedTracking::SagittaRByCircleFit(Double_t& centerX, Double_t& centerY, Double_t& sagittaR,
0685    std::vector<Double_t > r, std::vector<Double_t > phi, Double_t oInttPhi, Double_t emcalPhi){
0686    Double_t basePhi = TMath::Pi()/4 - oInttPhi;
0687    bool bFlip = false;
0688    if((emcalPhi + basePhi)  > TMath::Pi()/4) bFlip = true;
0689    
0690    TH2D* hHitMap = new TH2D("hHitMap", "hHitMap; x [cm]; y [cm]",\
0691       2000, -100., 100., 2000, -100., 100.);
0692 
0693    Int_t numOfHits = r.size();
0694    for(Int_t iHit = 0; iHit < numOfHits; iHit++){
0695       Double_t tempPhi = phi.at(iHit) + basePhi;
0696       if(bFlip) tempPhi = TMath::Pi()/2 - tempPhi;
0697       hHitMap->Fill(r.at(iHit)*cos(tempPhi), r.at(iHit)*sin(tempPhi));
0698    }
0699 
0700    Double_t tempCPhi = std::atan(centerY/centerX);
0701    if((tempCPhi < 0)&&(centerX < 0)) tempCPhi += TMath::Pi();
0702    else if((tempCPhi > 0)&&(centerX < 0)) tempCPhi -= TMath::Pi();
0703    tempCPhi = tempCPhi + basePhi;
0704    if(bFlip) tempCPhi = TMath::Pi()/2 - tempCPhi;
0705    Double_t tempCR = std::sqrt(centerX*centerX + centerY*centerY);
0706    Double_t tempCX = tempCR*std::cos(tempCPhi);
0707    Double_t tempCY = tempCR*std::sin(tempCPhi);
0708       
0709    // [0]: radius, [1]: x center, [2]: y center 
0710    auto fCircle = new TF1("fCircle", "std::sqrt([0]*[0]-(x-[1])*(x-[1]))+[2]", -300, 300);
0711    fCircle->SetParameter(0, sagittaR);
0712    fCircle->SetParameter(1, tempCX);
0713    fCircle->SetParameter(2, tempCY);
0714    
0715    hHitMap->Fit("fCircle", "Q");
0716    sagittaR = fCircle->GetParameter(0);
0717    centerX = fCircle->GetParameter(1);
0718    centerY = fCircle->GetParameter(2);
0719 
0720    delete fCircle;
0721    delete hHitMap;
0722 
0723    Double_t tempRetCPhi = std::atan(centerY/centerX);
0724    if((tempRetCPhi < 0)&&(centerX < 0)) tempRetCPhi += TMath::Pi();
0725    else if((tempRetCPhi > 0)&&(centerX < 0)) tempRetCPhi -= TMath::Pi();
0726    if(bFlip) tempRetCPhi = TMath::Pi()/2 - tempRetCPhi;
0727    tempRetCPhi = tempRetCPhi - basePhi;
0728    centerX = sagittaR*std::cos(tempRetCPhi);
0729    centerY = sagittaR*std::sin(tempRetCPhi);
0730 
0731    return;
0732 }
0733 
0734 bool InttSeedTracking::FindHitsOnCircle(Int_t& hitMatchId,\
0735    Double_t& hitX, Double_t& hitY,\
0736    Double_t sagittaR, Double_t centerX, Double_t centerY, std::vector<hitStruct > vHits,\
0737    Double_t dRThre){
0738    Double_t minDeltaR = 99999.;
0739 
0740    if(vHits.size() == 0) return false;
0741    for(Int_t iHit=0; iHit< vHits.size(); iHit++){
0742       Double_t tempHitX = vHits.at(iHit).r*std::cos(vHits.at(iHit).phi);
0743       Double_t tempHitY = vHits.at(iHit).r*std::sin(vHits.at(iHit).phi);
0744 
0745       Double_t dX = tempHitX - centerX;
0746       Double_t dY = tempHitY - centerY;
0747       Double_t dR = std::sqrt(dX*dX + dY*dY);
0748       
0749       Double_t deltaR = std::abs(dR - sagittaR);
0750       if(deltaR > dRThre) continue;
0751       if(minDeltaR > deltaR){
0752          hitMatchId = iHit;
0753          minDeltaR = deltaR;
0754       }
0755    }
0756    
0757    if(hitMatchId == 9999) return false;
0758    hitX = vHits.at(hitMatchId).r * std::cos(vHits.at(hitMatchId).phi);
0759    hitY = vHits.at(hitMatchId).r * std::sin(vHits.at(hitMatchId).phi);
0760 
0761    return true;
0762 }
0763 
0764 
0765 void InttSeedTracking::EstiVertex(tracKuma trk){
0766    Double_t vtxX = 0.;
0767    Double_t vtxY = 0.;
0768 
0769    CrossLineCircle(vtxX, vtxY, trk.getTrackSagX(), trk.getTrackSagY(), trk.getTrackSagR());
0770    Double_t vtxR = std::sqrt(vtxX*vtxX + vtxY*vtxY);
0771    Double_t vtxPhi = atan(vtxY/vtxX);
0772    if((vtxPhi < 0)&&(vtxX < 0)) vtxPhi += TMath::Pi();
0773    else if((vtxPhi > 0)&&(vtxY < 0)) vtxPhi -= TMath::Pi();
0774    trk.setHitR(0, vtxR);
0775    trk.setHitPhi(0, vtxPhi);
0776 
0777    Double_t trkPhi = -1/std::tan(vtxPhi);
0778    Double_t vtxZ = 0.;
0779    if(trk.getHitIs(6)){
0780       vtxZ = trk.getHitZ(6) - trk.getHitZ(6)/std::tan(trk.getTrackTheta());
0781    }
0782    trk.setHitZ(0, vtxZ);
0783 }
0784 
0785 Double_t InttSeedTracking::EstimateRecoTheta(tracKuma trk, Int_t type){
0786    // fitting weight for eta
0787    Double_t mvtxEtaW = 1.;
0788    Double_t inttEtaW = 1.;
0789    Double_t ecalEtaW = 1.;
0790    
0791    Double_t recoTheta = 9999.;
0792    if(type == 0){ // intt + EMcal
0793       recoTheta = (inttEtaW*trk.getHitTheta(4)\
0794          + inttEtaW*trk.getHitTheta(5)\
0795          + ecalEtaW*trk.getHitTheta(6))/(2*inttEtaW+ecalEtaW);
0796    }else if(type == 1){ // mvtx + intt + EMcal
0797       if((trk.getHitIs(1))&&(trk.getHitIs(2))\
0798       &&(trk.getHitIs(3))){
0799          recoTheta = (mvtxEtaW*trk.getHitTheta(1) \
0800             + mvtxEtaW*trk.getHitTheta(2)\
0801             + mvtxEtaW*trk.getHitTheta(3)\
0802             + inttEtaW*trk.getHitTheta(4)\
0803             + inttEtaW*trk.getHitTheta(5)\
0804             + ecalEtaW*trk.getHitTheta(6))\
0805             /(3*mvtxEtaW+2*inttEtaW+ecalEtaW);
0806       }else if((trk.getHitIs(2))&&(trk.getHitIs(3))){
0807          recoTheta = (mvtxEtaW*trk.getHitTheta(2)\
0808             + mvtxEtaW*trk.getHitTheta(3)\
0809             + inttEtaW*trk.getHitTheta(4)\
0810             + inttEtaW*trk.getHitTheta(5)\
0811             + ecalEtaW*trk.getHitTheta(6))/(2*mvtxEtaW+2*inttEtaW+ecalEtaW);
0812       }else if((trk.getHitIs(1))&&(trk.getHitIs(3))){
0813          recoTheta = (mvtxEtaW*trk.getHitTheta(1)\
0814             + mvtxEtaW*trk.getHitTheta(3)\
0815             + inttEtaW*trk.getHitTheta(4)\
0816             + inttEtaW*trk.getHitTheta(5)\
0817             + ecalEtaW*trk.getHitTheta(6))/(2*mvtxEtaW+2*inttEtaW+ecalEtaW);
0818       }else if((trk.getHitIs(1))&&(trk.getHitIs(2))){
0819          recoTheta = (mvtxEtaW*trk.getHitTheta(1)\
0820             + mvtxEtaW*trk.getHitTheta(2)\
0821             + inttEtaW*trk.getHitTheta(4)\
0822             + inttEtaW*trk.getHitTheta(5)\
0823             + ecalEtaW*trk.getHitTheta(6))/(2*mvtxEtaW+2*inttEtaW+ecalEtaW);
0824       }else if(trk.getHitIs(1)){
0825          recoTheta = (mvtxEtaW*trk.getHitTheta(1)
0826             + inttEtaW*trk.getHitTheta(4)\
0827             + inttEtaW*trk.getHitTheta(5)\
0828             + ecalEtaW*trk.getHitTheta(6))/(mvtxEtaW+2*inttEtaW+ecalEtaW);
0829       }else if(trk.getHitIs(2)){
0830          recoTheta = (mvtxEtaW*trk.getHitTheta(2)\
0831             + inttEtaW*trk.getHitTheta(4)\
0832             + inttEtaW*trk.getHitTheta(5)\
0833             + ecalEtaW*trk.getHitTheta(6))/(mvtxEtaW+2*inttEtaW+ecalEtaW);
0834       }else if(trk.getHitIs(3)){
0835          recoTheta = (mvtxEtaW*trk.getHitTheta(3)\
0836             + inttEtaW*trk.getHitTheta(4)
0837             + inttEtaW*trk.getHitTheta(5)\
0838             + ecalEtaW*trk.getHitTheta(6))/(mvtxEtaW+2*inttEtaW+ecalEtaW);
0839       }
0840 
0841    }
0842 
0843    return recoTheta;
0844 }
0845 
0846 void InttSeedTracking::ParticleIdentify(tracKuma trk){
0847    Double_t tempInttPhi = 0.;
0848    if(trk.getHitIs(4)) tempInttPhi = trk.getHitPhi(4);
0849    else tempInttPhi = trk.getHitPhi(5);
0850    Double_t tempCalPhi = trk.getHitPhi(6);
0851    
0852    if(tempInttPhi < 0) tempInttPhi += 2*TMath::Pi();
0853    if(tempCalPhi < 0)  tempCalPhi += 2*TMath::Pi();
0854 
0855    Int_t charge = (tempCalPhi - tempCalPhi)/std::abs(tempCalPhi - tempCalPhi);
0856 
0857    trk.setTrackCharge(charge);
0858 
0859    bool tempElectronIs = false;
0860    // ChecKumaDaYo!!! You need to optimize the threshold to identify electron
0861    if(trk.getTrackE()/trk.getTrackP() > m_EOverPElectron) tempElectronIs = true;
0862    trk.setTrackElectronIs(tempElectronIs);
0863 
0864 }
0865 
0866 // 获取track在各个子探测器的R,Phi，命中点极坐标返回，以及检查是否超过3个hit点
0867 bool InttSeedTracking::ReturnHitsRPhiVect(std::vector<Double_t>& hitR, std::vector<Double_t>& hitPhi,\
0868     std::vector<Int_t> subDetSet, tracKuma trk)
0869 {
0870     for(Int_t iDet = 0; iDet < subDetSet.size(); iDet++)
0871     {
0872         Int_t detID = subDetSet.at(iDet);
0873         if(!trk.getHitIs(detID)) continue;
0874 
0875         hitR.push_back(trk.getHitR(detID)); 
0876         hitPhi.push_back(trk.getHitPhi(detID));
0877     }
0878     if(hitR.size() < 3) return false;
0879     return true;
0880 }
0881 
0882 
0883 void InttSeedTracking::SetHitParaInTrack(tracKuma& trk, Int_t detId, hitStruct hitSt)
0884 {
0885    trk.setHitIs(detId, true);
0886    trk.setHitR(detId, hitSt.r);
0887    trk.setHitZ(detId, hitSt.z);
0888    trk.setHitPhi(detId, hitSt.phi);
0889    Double_t hitTheta = 2*atan(std::exp(-hitSt.eta));
0890    trk.setHitTheta(detId, hitTheta);
0891 }
0892 
0893 #endif // #ifdef InttSeedTracking_cxx
0894 
0895 
0896 
0897
[ 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!