InttSeedTrackPerformance-Jingyu/src/JYInttSeedTracking.cxx

0001 // InttSeedTracking::InttSeedTracking
0002 //   ├── HitMatching
0003 //   │     ├── TempINTTIOMatching (INTT 内外层匹配)
0004 //   │     ├── TempInttCalMatch (INTT 和 EMCal 匹配)
0005 //   │     │     ├── TempCalcdPhidR (计算 dPhi/dR)
0006 //   │     │     └── AddHCalE (加入 HCal 能量)
0007 //   │     └── tracKuma 构造轨迹
0008 //   ├── TrackPropertiesEstimation
0009 //   │     ├── TrackPtEstimation (估算横向动量 pT)
0010 //   │     │     ├── Set3PointsXY (设置 3 点坐标)
0011 //   │     │     ├── RoughEstiSagittaCenter3Point (粗略估算 Sagitta 中心与半径)
0012 //   │     │     ├── AccuratePtEstimation (精确估算动量 pT)
0013 //   │     │     │     ├── SagittaRByCircleFit (弧线拟合，精确估算 Sagitta)
0014 //   │     │     │     ├── AddMvtxHits (匹配 MVTX 命中点)
0015 //   │     │     │     └── CalcSagittaPt (基于 Sagitta 半径计算横向动量)
0016 //   │     │     └── FitFunctionPt (基于 dPhi 校正动量)
0017 //   │     ├── EstimateRecoTheta (估算重建角度 θ)
0018 //   │     │     └── 不同探测器组合加权计算 θ
0019 //   │     └── CrossLineCircle (计算顶点位置)
0020 //   └── 数据清理与存储
0021 //         └── clear() 和 shrink_to_fit() 清理命中点数据
0022
0023 #include <TStyle.h>
0024 #include <TCanvas.h>
0025
0026 #include "JYInttSeedTracking.h"
0027 #include "SPHTracKuma.h"
0028
0029 InttSeedTracking::InttSeedTracking(std::vector<tracKuma>& tracks,\
0030    std::vector<hitStruct > vFMvtxHits,\
0031    std::vector<hitStruct > vSMvtxHits, std::vector<hitStruct > vTMvtxHits,\
0032    std::vector<hitStruct > vIInttHits, std::vector<hitStruct > vOInttHits,\
0033    std::vector<hitStruct > vEmcalHits,\
0034    std::vector<hitStruct > vIHCalHits, std::vector<hitStruct > vOHCalHits)
0035 {
0036     HitMatching(tracks, vFMvtxHits, vSMvtxHits, vTMvtxHits, vIInttHits, vOInttHits, vEmcalHits, vIHCalHits, vOHCalHits);
0037     for(Int_t iTrk = 0; iTrk < tracks.size(); iTrk++)
0038     {
0039         TrackPropertiesEstimation(tracks.at(iTrk), vFMvtxHits, vSMvtxHits, vTMvtxHits);
0040     }
0041
0042     // Vertex estimation using all tracks
0043     // VertexFinder(); //
0044
0045     // Re-finding tracks for single INTT hits or only mvtx using estimated vertex
0046     // ReFindTrack(); //
0047 }
0048
0049 InttSeedTracking::~InttSeedTracking()
0050 {
0051     m_fMvtxHits.clear();
0052     m_sMvtxHits.clear();
0053     m_tMvtxHits.clear();
0054     m_iInttHits.clear();
0055     m_oInttHits.clear();
0056     m_emcalHits.clear();
0057     m_iHCalHits.clear();
0058     m_oHCalHits.clear();
0059
0060     m_fMvtxHits.shrink_to_fit();
0061     m_sMvtxHits.shrink_to_fit();
0062     m_tMvtxHits.shrink_to_fit();
0063     m_iInttHits.shrink_to_fit();
0064     m_oInttHits.shrink_to_fit();
0065     m_emcalHits.shrink_to_fit();
0066     m_iHCalHits.shrink_to_fit();
0067     m_oHCalHits.shrink_to_fit();
0068 }
0069
0070
0071 // Hit Matching to get which hit point from one track
0072 void InttSeedTracking::HitMatching(std::vector<tracKuma>& tracks,\
0073    std::vector<hitStruct> vFMvtxHits,\
0074    std::vector<hitStruct> vSMvtxHits, std::vector<hitStruct> vTMvtxHits,\
0075    std::vector<hitStruct> vIInttHits, std::vector<hitStruct> vOInttHits,\
0076    std::vector<hitStruct> vEmcalHits,\
0077    std::vector<hitStruct> vIHCalHits, std::vector<hitStruct> vOHcalHits)
0078 {
0079     Int_t numOfInttClus = vOInttHits.size();
0080     if(numOfInttClus == 0) return;
0081
0082     for(Int_t iOInttClus = 0; iOInttClus < numOfInttClus; iOInttClus++)
0083     {
0084         Double_t oInttPhi = vOInttHits.at(iOInttClus).phi;
0085         // Double_t refPhi = oInttPhi;
0086         // refPhi = TMath::Pi()/2 - refPhi;
0087
0088         // oInttPhi += refPhi;
0089
0090         // Int_t closetIINTTID = TempINTTIOMatching(oInttPhi, refPhi, vIInttHits);
0091         // if(closetIINTTID == 9999) continue;
0092         // Int_t matchEcalID = TempInttCalMatch(closetIINTTID, iOInttClus, refPhi, vIInttHits, vOInttHits, vEmcalHits);
0093         Int_t closetIINTTID = TempINTTIOMatching(oInttPhi, vIInttHits);
0094         if(closetIINTTID == 9999) continue;
0095         Int_t matchEcalID = TempInttCalMatch(closetIINTTID, iOInttClus, vIInttHits, vOInttHits, vEmcalHits);
0096         if(matchEcalID == 9999) continue;
0097
0098         tracKuma trk;
0099         // trk.setRefPhi(refPhi);
0100         Double_t tempTheta = 0.;
0101
0102         trk.setHitIs(0, true); // 0 vertex, 1-3 mvtx, 4-5 ioINTT, 6 emcal,
0103         trk.setHitR(0, 0.);
0104         trk.setHitZ(0, 0.);
0105         trk.setHitPhi(0, 0.);
0106         trk.setHitTheta(0, 0.);
0107
0108         trk.setHitIs(4, true);
0109         trk.setHitR(4, vIInttHits.at(closetIINTTID).r);
0110         trk.setHitZ(4, vIInttHits.at(closetIINTTID).z);
0111         trk.setHitPhi(4, vIInttHits.at(closetIINTTID).phi);
0112         tempTheta = 2*atan(std::exp(-vIInttHits.at(closetIINTTID).eta));
0113         trk.setHitTheta(4, tempTheta);
0114
0115         trk.setHitIs(5, true);
0116         trk.setHitR(5, vOInttHits.at(iOInttClus).r);
0117         trk.setHitZ(5, vOInttHits.at(iOInttClus).z);
0118         trk.setHitPhi(5, vOInttHits.at(iOInttClus).phi);
0119         tempTheta = 2*atan(std::exp(-vOInttHits.at(iOInttClus).eta));
0120         trk.setHitTheta(5, tempTheta);
0121
0122         trk.setHitIs(6, true);
0123         trk.setHitR(6, vEmcalHits.at(matchEcalID).r);
0124         trk.setHitZ(6, vEmcalHits.at(matchEcalID).z);
0125         trk.setHitPhi(6, vEmcalHits.at(matchEcalID).phi);
0126         tempTheta = 2*atan(std::exp(-vEmcalHits.at(matchEcalID).eta));
0127         trk.setHitTheta(6, tempTheta);
0128
0129         Double_t calE = vEmcalHits.at(matchEcalID).energy;
0130         calE = AddHCalE(vEmcalHits.at(matchEcalID).phi, calE, vIHCalHits, vOHcalHits); // kuma???
0131         trk.setTrackE(calE);
0132
0133         tracks.push_back(trk);
0134     }
0135 }
0136
0137 // match the inner INTT cluster and outer INTT cluster, return closest iINTTCluID
0138 Int_t InttSeedTracking::TempINTTIOMatching(Double_t oINTTPhi, std::vector<hitStruct> vIInttHits)
0139 {
0140     Double_t minDeltaPhi = 9999.;
0141     Int_t closestIINTTCluID = 9999;
0142     Int_t numOfIInttClu = vIInttHits.size();
0143
0144     if(numOfIInttClu == 0) return 9999;
0145     for(Int_t iINTTClu = 0; iINTTClu < numOfIInttClu; iINTTClu++)
0146     {
0147         Double_t iInttPhi =  vIInttHits.at(iINTTClu).phi;
0148         Double_t dPhi = std::abs(iInttPhi - oINTTPhi);
0149         if(std::abs(dPhi) > TMath::Pi()) dPhi += (-1) * (dPhi/std::abs(dPhi)) * 2*TMath::Pi(); // turn to -PI ~ PI
0150
0151         // checkumaDAYO!!! need to optimize the matching range.
0152         if((minDeltaPhi > dPhi)&&(dPhi > m_InttMatchPhiMin)&&(dPhi < m_InttMatchPhiMax))
0153         {
0154             minDeltaPhi = dPhi;
0155             closestIINTTCluID = iINTTClu;
0156         }
0157     }
0158     return closestIINTTCluID;
0159 }
0160
0161 // calculate dPhi_intt dR_intt, return dPhi/dr, applied to TempInttCalMatch
0162 Double_t InttSeedTracking::TempCalcdPhidR(Int_t iInttID, Int_t oInttID, std::vector<hitStruct> vIInttHits, std::vector<hitStruct> vOInttHits)
0163 {
0164     Double_t dPhidR = 0.;
0165     Double_t iInttPhi = vIInttHits.at(iInttID).phi;
0166     Double_t oInttPhi = vOInttHits.at(oInttID).phi;
0167
0168     Double_t dPhi = oInttPhi - iInttPhi;
0169     if(std::abs(dPhi) > TMath::Pi()) dPhi += (-1) * (dPhi/std::abs(dPhi)) * 2*TMath::Pi();
0170     Double_t dR = vOInttHits.at(oInttID).r - vIInttHits.at(iInttID).r;
0171
0172     dPhidR = dPhi/dR;
0173
0174     return dPhidR;
0175 }
0176
0177 // calculate the dphi_cal, return all satisify emcal cluster ID
0178 Int_t InttSeedTracking::TempInttCalMatch(Int_t iInttID, Int_t oInttID, std::vector<hitStruct> vIInttHits, std::vector<hitStruct> vOInttHits, std::vector<hitStruct> vEmcalHits)
0179 {
0180    Int_t matchiECalID = 9999;
0181    Double_t matchiECalE = 0.;
0182    Int_t numOfECalT = vEmcalHits.size();
0183
0184    for(Int_t iECalT = 0; iECalT < numOfECalT ;iECalT++)
0185    {
0186         Double_t ecalE = vEmcalHits.at(iECalT).energy;
0187         // checkumaDAYO!!! need to optimize the calorimeter threshold energy.
0188         if(ecalE < m_EcalEThre) continue;
0189
0190         Double_t dPhidR = TempCalcdPhidR(iInttID, oInttID, vIInttHits, vOInttHits);
0191         Double_t dRInttEmCal = m_ECalR - vOInttHits.at(oInttID).r;
0192         Double_t phiRangeMin = 0.;
0193         Double_t phiRangeMax = 0.;
0194
0195         // checkumaDAYO!!! need to optimize the matching range. ?TMath::Pi()/2 should be Phi_oINTT ?dPhi_cal = dPhidR * dRInttEmCal
0196         if(dPhidR < 0)
0197         {
0198             phiRangeMin = vOInttHits.at(oInttID).phi + (dPhidR * dRInttEmCal - 0.177);
0199             phiRangeMax = vOInttHits.at(oInttID).phi + 0.087;
0200         }
0201         else
0202         {
0203             phiRangeMin = vOInttHits.at(oInttID).phi - 0.087;
0204             phiRangeMax = vOInttHits.at(oInttID).phi + (dPhidR * dRInttEmCal + 0.177);
0205         }
0206
0207         Double_t emcalPhi = vEmcalHits.at(iECalT).phi;
0208         if((phiRangeMin < - TMath::Pi())&&(phiRangeMin > - TMath::Pi())&&(emcalPhi > 0)) emcalPhi -= 2*TMath::Pi();
0209         else if((phiRangeMin < TMath::Pi())&&(phiRangeMax > TMath::Pi())&&(emcalPhi < 0)) emcalPhi += 2*TMath::Pi();
0210
0211         if((phiRangeMin < emcalPhi)&&(emcalPhi < phiRangeMax))
0212         {
0213             if(matchiECalE < ecalE)
0214             {
0215                 matchiECalID = iECalT;
0216                 matchiECalE = ecalE;
0217             }
0218         }
0219         // std::cout << "ecalE = " << ecalE << std::endl;
0220    }
0221    // std::cout << "matchiECalE = " << matchiECalE << std::endl;
0222    return matchiECalID;
0223 }
0224
0225 // Add HCalE, return calE = emcalE + iHcalE + oHcalE
0226 Double_t InttSeedTracking::AddHCalE(Double_t emcalPhi, Double_t emcalE, std::vector<hitStruct > vIHCalHits, std::vector<hitStruct > vOHCalHits)
0227 {
0228     // inner HCal match
0229     Int_t matchIhcalId = 99999;
0230     Double_t closePhiIhcal = 99999;
0231
0232     for(Int_t iIhcal = 0; iIhcal < vIHCalHits.size(); iIhcal++)
0233     {
0234         if(closePhiIhcal > std::abs(vIHCalHits.at(iIhcal).phi - emcalPhi))
0235         {
0236             closePhiIhcal = std::abs(vIHCalHits.at(iIhcal).phi - emcalPhi);
0237             matchIhcalId = iIhcal;
0238         }
0239     }
0240
0241     // outer HCal match
0242     Int_t matchOhcalId = 99999;
0243     Double_t closePhiOhcal = 99999;
0244
0245     for(Int_t iOhcal = 0; iOhcal < vOHCalHits.size(); iOhcal++)
0246     {
0247         if(closePhiOhcal > std::abs(vOHCalHits.at(iOhcal).phi - emcalPhi))
0248         {
0249             closePhiOhcal = std::abs(vOHCalHits.at(iOhcal).phi - emcalPhi);
0250             matchOhcalId = iOhcal;
0251         }
0252     }
0253
0254     Double_t iHcalE = 0.;
0255     Double_t oHcalE = 0.;
0256     if(matchIhcalId != 99999) iHcalE = vIHCalHits.at(matchIhcalId).energy;
0257     if(matchOhcalId != 99999) oHcalE = vOHCalHits.at(matchOhcalId).energy;
0258
0259     Double_t calE = emcalE + iHcalE + oHcalE;
0260
0261     return calE;
0262 }
0263
0264
0265 // ├────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
0266 // ├────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
0267 // ├────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
0268 // ├────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
0269 // ├────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
0270
0271
0272 // TrackPropertiesEstimation, get the track pt
0273 void InttSeedTracking::TrackPropertiesEstimation(tracKuma& trk, std::vector<hitStruct> vFMvtxHits, std::vector<hitStruct> vSMvtxHits, std::vector<hitStruct> vTMvtxHits)
0274 {
0275     Double_t recoPt = TrackPtEstimation(trk, vFMvtxHits, vSMvtxHits, vTMvtxHits);
0276     trk.setTrackPt(recoPt);
0277
0278     Double_t recoTheta = EstimateRecoTheta(trk, 1);
0279     trk.setTrackTheta(recoTheta);
0280
0281     Double_t recoP = recoPt/sin(recoTheta);
0282     Double_t recoE = trk.getTrackE();
0283     Double_t EOverP = recoE/recoP;
0284     trk.setTrackP(recoP);
0285
0286     // Double_t vtxX = 0.;
0287     // Double_t vtxY = 0.;
0288     EstiVertex(trk);
0289     ParticleIdentify(trk);
0290
0291     // // EstiVertex(vtxX, vtxY, trk.getTrackSagR(), trk.getTrackSagX(), trk.getTrackSagY());   1227
0292     // CrossLineCircle(vtxX, vtxY, trk.getTrackSagX(), trk.getTrackSagY(), trk.getTrackSagR());
0293
0294     // Double_t vtxR = std::sqrt(vtxX*vtxX + vtxY*vtxY);
0295     // Double_t vtxPhi = atan(vtxY/vtxX);
0296     // if((vtxPhi < 0)&&(vtxX < 0)) vtxPhi += TMath::Pi();
0297     // else if((vtxPhi > 0)&&(vtxY < 0)) vtxPhi -= TMath::Pi();
0298
0299     // trk.setHitR(0, vtxR);
0300     // trk.setHitPhi(0, vtxPhi);
0301
0302     // Double_t trkPhi = -1/std::tan(vtxPhi);
0303     // Double_t vtxZ = 0.;
0304     // if(trk.getHitIs(6))
0305     // {
0306     //     vtxZ = trk.getHitZ(6) - trk.getHitZ(6)/std::tan(recoTheta);
0307     // }
0308     // trk.setHitZ(0, vtxZ);
0309 }
0310
0311 // Track Pt Estimation
0312 Double_t InttSeedTracking::TrackPtEstimation(tracKuma& trk, std::vector<hitStruct> vFMvtxHits, std::vector<hitStruct > vSMvtxHits, std::vector<hitStruct> vTMvtxHits)
0313 {
0314     Double_t sagittaR = 9999.;
0315     Double_t sagittaCenterX = 0.;
0316     Double_t sagittaCenterY = 0.;
0317     Double_t HitsXY[3][2];
0318     Set3PointsXY(HitsXY, trk, 2);
0319     RoughEstiSagittaCenter3Point(sagittaR, sagittaCenterX, sagittaCenterY, HitsXY);
0320
0321     Double_t trackPt = AccuratePtEstimation(sagittaR, sagittaCenterX, sagittaCenterY, trk, vFMvtxHits, vSMvtxHits, vTMvtxHits);
0322
0323     Double_t dPhiOInttEmcal = dPhiOInttEmcalEsti(trk);
0324     trackPt = FitFunctionPt(dPhiOInttEmcal);
0325
0326     return trackPt;
0327 }
0328
0329 // a RoughEstiSagittaCenter3Point to get initial value for fit
0330 void InttSeedTracking::RoughEstiSagittaCenter3Point(Double_t& sagittaR, Double_t& sagittaCenterX, Double_t& sagittaCenterY, Double_t HitsXY[3][2])
0331 {
0332     Double_t iSlope = 0.;
0333     Double_t iSection = 0.;
0334     CalcPerpBis(iSlope, iSection, HitsXY[0], HitsXY[2]);
0335
0336     Double_t oSlope = 0.;
0337     Double_t oSection = 0.;
0338     CalcPerpBis(oSlope, oSection, HitsXY[1], HitsXY[2]);
0339
0340     sagittaCenterX = - (oSection - iSection)/(oSlope - iSlope);
0341     sagittaCenterY = iSlope * sagittaCenterX + iSection;
0342
0343     sagittaR = std::sqrt((HitsXY[0][0] - sagittaCenterX)*(HitsXY[0][0] - sagittaCenterX) + (HitsXY[0][1] - sagittaCenterY)*(HitsXY[0][1] - sagittaCenterY));
0344 }
0345
0346 // AccuratePtEstimation, return track pt from sagitta_Mvtx_Intt_Emcal
0347 Double_t InttSeedTracking::AccuratePtEstimation(Double_t sagittaR, Double_t centerX, Double_t centerY, tracKuma& trk, std::vector<hitStruct> vFMvtxHits, std::vector<hitStruct> vSMvtxHits, std::vector<hitStruct> vTMvtxHits)
0348 {
0349     Double_t sagittaVtxInttPt = 99999.;
0350     Double_t sagittaInttEmcalPt = 99999.;
0351     Double_t sagittaVtxInttEmcalPt = 99999.;
0352     Double_t sagittaMvtxInttEmcalPt = 99999.;
0353     Double_t sagittaVtxMvtxInttEmcalPt = 99999.;
0354
0355     // checkumaDaYo!!!
0356     std::vector<Int_t> subDetIds_InttEmcal = {4, 5, 6};
0357     std::vector<Double_t> hitsR_InttEmcal = {};
0358     std::vector<Double_t> hitsPhi_InttEmcal = {};
0359     if(ReturnHitsRPhiVect(hitsR_InttEmcal, hitsPhi_InttEmcal, subDetIds_InttEmcal, trk))
0360     {
0361         Double_t tempHitOInttPhi = trk.getHitPhi(4);
0362         Double_t tempHitEmcalPhi = trk.getHitPhi(6);
0363         SagittaRByCircleFit(centerX, centerY, sagittaR, hitsR_InttEmcal, hitsPhi_InttEmcal, tempHitOInttPhi, tempHitEmcalPhi);
0364
0365         sagittaMvtxInttEmcalPt = CalcSagittaPt(sagittaR);
0366     }
0367
0368     Double_t dRThre_Mvtx = 0.01;
0369     AddMvtxHits(trk, vFMvtxHits, vSMvtxHits, vTMvtxHits, dRThre_Mvtx);
0370
0371     std::vector<Int_t> subDetIds_MvtxInttEmcal = {1, 2, 3, 4, 5, 6};
0372     std::vector<Double_t> hitsR_MvtxInttEmcal = {};
0373     std::vector<Double_t> hitsPhi_MvtxInttEmcal = {};
0374     if(ReturnHitsRPhiVect(hitsR_MvtxInttEmcal, hitsPhi_MvtxInttEmcal, subDetIds_MvtxInttEmcal, trk))
0375     {
0376         Double_t tempHitOInttPhi = trk.getHitPhi(4);
0377         Double_t tempHitEmcalPhi = trk.getHitPhi(6);
0378         SagittaRByCircleFit(centerX, centerY, sagittaR, hitsR_MvtxInttEmcal, hitsPhi_MvtxInttEmcal,
0379         tempHitOInttPhi, tempHitEmcalPhi);
0380         sagittaMvtxInttEmcalPt = CalcSagittaPt(sagittaR);
0381     }
0382
0383     trk.setTrackSagR(sagittaR);
0384     trk.setTrackSagX(centerX);
0385     trk.setTrackSagY(centerY);
0386
0387     return sagittaMvtxInttEmcalPt;
0388 }
0389
0390 // calculate the SagittaR by CircleFit ???
0391 void InttSeedTracking::SagittaRByCircleFit(Double_t& centerX, Double_t& centerY, Double_t& sagittaR, std::vector<Double_t > r, std::vector<Double_t > phi, Double_t oInttPhi, Double_t emcalPhi)
0392 {
0393     Double_t basePhi = TMath::Pi()/4 - oInttPhi;
0394     bool bFlip = false;
0395     if((emcalPhi + basePhi)  > TMath::Pi()/4) bFlip = true;
0396
0397     TH2D* hHitMap = new TH2D("hHitMap", "hHitMap; x [cm]; y [cm]", 2000, -100., 100., 2000, -100., 100.);
0398
0399     Int_t numOfHits = r.size();
0400     for(Int_t iHit = 0; iHit < numOfHits; iHit++)
0401     {
0402         Double_t tempPhi = phi.at(iHit) + basePhi;
0403         if(bFlip) tempPhi = TMath::Pi()/2 - tempPhi;
0404         hHitMap->Fill(r.at(iHit)*cos(tempPhi), r.at(iHit)*sin(tempPhi));
0405     }
0406
0407     Double_t tempCPhi = std::atan(centerY/centerX);
0408     if((tempCPhi < 0)&&(centerX < 0)) tempCPhi += TMath::Pi();
0409     else if((tempCPhi > 0)&&(centerX < 0)) tempCPhi -= TMath::Pi();
0410     tempCPhi = tempCPhi + basePhi;
0411     if(bFlip) tempCPhi = TMath::Pi()/2 - tempCPhi;
0412     Double_t tempCR = std::sqrt(centerX*centerX + centerY*centerY);
0413     Double_t tempCX = tempCR*std::cos(tempCPhi);
0414     Double_t tempCY = tempCR*std::sin(tempCPhi);
0415
0416     // [0]: radius, [1]: x center, [2]: y center
0417     auto fCircle = new TF1("fCircle", "std::sqrt([0]*[0]-(x-[1])*(x-[1]))+[2]", -300, 300);
0418     fCircle->SetParameter(0, sagittaR);
0419     fCircle->SetParameter(1, tempCX);
0420     fCircle->SetParameter(2, tempCY);
0421
0422     hHitMap->Fit("fCircle", "QM");
0423     sagittaR = fCircle->GetParameter(0);
0424     centerX = fCircle->GetParameter(1);
0425     centerY = fCircle->GetParameter(2);
0426
0427     delete fCircle;
0428     delete hHitMap;
0429
0430     Double_t tempRetCPhi = std::atan(centerY/centerX);
0431     if((tempRetCPhi < 0)&&(centerX < 0)) tempRetCPhi += TMath::Pi();
0432     else if((tempRetCPhi > 0)&&(centerX < 0)) tempRetCPhi -= TMath::Pi();
0433     if(bFlip) tempRetCPhi = TMath::Pi()/2 - tempRetCPhi;
0434     tempRetCPhi = tempRetCPhi - basePhi;
0435     centerX = sagittaR*std::cos(tempRetCPhi);
0436     centerY = sagittaR*std::sin(tempRetCPhi);
0437
0438     return;
0439 }
0440
0441 //  AddMvtxHits
0442 void InttSeedTracking::AddMvtxHits(tracKuma& trk, std::vector<hitStruct> vFMvtxHits,\
0443    std::vector<hitStruct> vSMvtxHits, std::vector<hitStruct> vTMvtxHits, Double_t dRThre)
0444 {
0445     Double_t sagiR = trk.getTrackSagR();
0446     Double_t sagiX = trk.getTrackSagX();
0447     Double_t sagiY = trk.getTrackSagY();
0448
0449     // 逐层检查，匹配一层，就在trk里set一层
0450     for(Int_t iMvtxLay = 0; iMvtxLay < 3; iMvtxLay++)
0451     {
0452         Int_t mvtxMatchId = 9999;
0453         Double_t mvtxClusX = 9999.;
0454         Double_t mvtxClusY = 9999.;
0455         if(iMvtxLay == 0)
0456         {
0457             if(FindHitsOnCircle(mvtxMatchId, mvtxClusX, mvtxClusY, sagiR, sagiX, sagiY, vFMvtxHits, dRThre))
0458             {
0459                 trk.setHitIs(1, true);
0460                 trk.setHitR(1, vFMvtxHits.at(mvtxMatchId).r);
0461                 trk.setHitPhi(1, vFMvtxHits.at(mvtxMatchId).phi);
0462                 Double_t tempTheta = 2*atan(std::exp(-vFMvtxHits.at(mvtxMatchId).eta));
0463                 trk.setHitTheta(1, tempTheta);
0464             }
0465         }
0466         else if(iMvtxLay == 1)
0467         {
0468             if(FindHitsOnCircle(mvtxMatchId, mvtxClusX, mvtxClusY, sagiR, sagiX, sagiY, vSMvtxHits, dRThre))
0469             {
0470                 trk.setHitIs(2, true);
0471                 trk.setHitR(2, vSMvtxHits.at(mvtxMatchId).r);
0472                 trk.setHitPhi(2, vSMvtxHits.at(mvtxMatchId).phi);
0473                 Double_t tempTheta = 2*atan(std::exp(-vSMvtxHits.at(mvtxMatchId).eta));
0474                 trk.setHitTheta(2, tempTheta);
0475             }
0476         }
0477         else if(iMvtxLay == 2)
0478         {
0479             if(FindHitsOnCircle(mvtxMatchId, mvtxClusX, mvtxClusY, sagiR, sagiX, sagiY, vTMvtxHits, dRThre))
0480             {
0481                 trk.setHitIs(3, true);
0482                 trk.setHitR(3, vTMvtxHits.at(mvtxMatchId).r);
0483                 trk.setHitPhi(3, vTMvtxHits.at(mvtxMatchId).phi);
0484                 Double_t tempTheta = 2*atan(std::exp(-vTMvtxHits.at(mvtxMatchId).eta));
0485                 trk.setHitTheta(3, tempTheta);
0486             }
0487         }
0488     }
0489 }
0490
0491 // loop mvtx ilayer cluster,use delta_R = dR-sagittaR to check is match, get mvtxClusX and mvtxClusY ->FindHitsOnCircle
0492 bool InttSeedTracking::FindHitsOnCircle(Int_t& hitMatchId,\
0493    Double_t& hitX, Double_t& hitY,\
0494    Double_t sagittaR, Double_t centerX, Double_t centerY, std::vector<hitStruct > vHits,\
0495    Double_t dRThre)
0496 {
0497     Double_t minDeltaR = 99999.;
0498
0499     if(vHits.size() == 0) return false;
0500
0501     for(Int_t iHit=0; iHit< vHits.size(); iHit++)
0502     {
0503         Double_t tempHitX = vHits.at(iHit).r*std::cos(vHits.at(iHit).phi);
0504         Double_t tempHitY = vHits.at(iHit).r*std::sin(vHits.at(iHit).phi);
0505
0506         Double_t dX = tempHitX - centerX;
0507         Double_t dY = tempHitY - centerY;
0508         Double_t dR = std::sqrt(dX*dX + dY*dY);
0509
0510         Double_t deltaR = std::abs(dR - sagittaR);
0511         if(deltaR > dRThre) continue;
0512         if(minDeltaR > deltaR)
0513         {
0514             hitMatchId = iHit;
0515             minDeltaR = deltaR;
0516         }
0517     }
0518
0519     if(hitMatchId == 9999) return false;  // no mvtxcluster to match
0520     hitX = vHits.at(hitMatchId).r * std::cos(vHits.at(hitMatchId).phi);
0521     hitY = vHits.at(hitMatchId).r * std::sin(vHits.at(hitMatchId).phi);
0522
0523     return true;
0524 }
0525
0526 // ReturnHitsRPhiVect,
0527 bool InttSeedTracking::ReturnHitsRPhiVect(std::vector<Double_t>& hitR, std::vector<Double_t>& hitPhi, std::vector<Int_t> subDetSet, tracKuma trk)
0528 {
0529     for(Int_t iDet = 0; iDet < subDetSet.size(); iDet++)
0530     {
0531         Int_t detID = subDetSet.at(iDet);
0532         if(!trk.getHitIs(detID)) continue;
0533
0534         hitR.push_back(trk.getHitR(detID));
0535         hitPhi.push_back(trk.getHitPhi(detID));
0536     }
0537     if(hitR.size() < 3) return false;
0538     return true;
0539 }
0540
0541 // Track Theta(eta?) Estimation
0542 Double_t InttSeedTracking::EstimateRecoTheta(tracKuma trk, Int_t type)
0543 {
0544     // fitting weight for eta
0545     Double_t mvtxEtaW = 1.;
0546     Double_t inttEtaW = 1.;
0547     Double_t ecalEtaW = 1.;
0548
0549     Double_t recoTheta = 9999.;
0550     if(type == 0) // intt + EMcal
0551     {
0552         recoTheta = (inttEtaW*trk.getHitTheta(4)\
0553           + inttEtaW*trk.getHitTheta(5)\
0554           + ecalEtaW*trk.getHitTheta(6))/(2*inttEtaW+ecalEtaW);
0555     }
0556     else if(type == 1) // mvtx + intt + EMcal, 3 layers mvtx 1-3 bool getHitIs
0557     {
0558         if((trk.getHitIs(1))&&(trk.getHitIs(2))&&(trk.getHitIs(3)))
0559         {
0560             recoTheta = (mvtxEtaW*trk.getHitTheta(1) \
0561               + mvtxEtaW*trk.getHitTheta(2)\
0562               + mvtxEtaW*trk.getHitTheta(3)\
0563               + inttEtaW*trk.getHitTheta(4)\
0564               + inttEtaW*trk.getHitTheta(5)\
0565               + ecalEtaW*trk.getHitTheta(6))\
0566               /(3*mvtxEtaW+2*inttEtaW+ecalEtaW);
0567         }
0568         else if((trk.getHitIs(2))&&(trk.getHitIs(3)))
0569         {
0570             recoTheta = (mvtxEtaW*trk.getHitTheta(2)\
0571               + mvtxEtaW*trk.getHitTheta(3)\
0572               + inttEtaW*trk.getHitTheta(4)\
0573               + inttEtaW*trk.getHitTheta(5)\
0574               + ecalEtaW*trk.getHitTheta(6))/(2*mvtxEtaW+2*inttEtaW+ecalEtaW);
0575         }
0576         else if((trk.getHitIs(1))&&(trk.getHitIs(3)))
0577         {
0578             recoTheta = (mvtxEtaW*trk.getHitTheta(1)\
0579               + mvtxEtaW*trk.getHitTheta(3)\
0580               + inttEtaW*trk.getHitTheta(4)\
0581               + inttEtaW*trk.getHitTheta(5)\
0582               + ecalEtaW*trk.getHitTheta(6))/(2*mvtxEtaW+2*inttEtaW+ecalEtaW);
0583         }
0584         else if((trk.getHitIs(1))&&(trk.getHitIs(2)))
0585         {
0586             recoTheta = (mvtxEtaW*trk.getHitTheta(1)\
0587               + mvtxEtaW*trk.getHitTheta(2)\
0588               + inttEtaW*trk.getHitTheta(4)\
0589               + inttEtaW*trk.getHitTheta(5)\
0590               + ecalEtaW*trk.getHitTheta(6))/(2*mvtxEtaW+2*inttEtaW+ecalEtaW);
0591         }
0592         else if(trk.getHitIs(1))
0593         {
0594             recoTheta = (mvtxEtaW*trk.getHitTheta(1)
0595               + inttEtaW*trk.getHitTheta(4)\
0596               + inttEtaW*trk.getHitTheta(5)\
0597               + ecalEtaW*trk.getHitTheta(6))/(mvtxEtaW+2*inttEtaW+ecalEtaW);
0598         }
0599         else if(trk.getHitIs(2))
0600         {
0601             recoTheta = (mvtxEtaW*trk.getHitTheta(2)\
0602               + inttEtaW*trk.getHitTheta(4)\
0603               + inttEtaW*trk.getHitTheta(5)\
0604               + ecalEtaW*trk.getHitTheta(6))/(mvtxEtaW+2*inttEtaW+ecalEtaW);
0605         }
0606         else if(trk.getHitIs(3))
0607         {
0608             recoTheta = (mvtxEtaW*trk.getHitTheta(3)\
0609               + inttEtaW*trk.getHitTheta(4)
0610               + inttEtaW*trk.getHitTheta(5)\
0611               + ecalEtaW*trk.getHitTheta(6))/(mvtxEtaW+2*inttEtaW+ecalEtaW);
0612         }
0613     }
0614
0615     return recoTheta;
0616 }
0617
0618 // track curve vertex Estimation, get vX and vY, what method be used?
0619 void InttSeedTracking::EstiVertex(tracKuma trk)
0620 {
0621     Double_t vtxX = 0.;
0622     Double_t vtxY = 0.;
0623
0624     CrossLineCircle(vtxX, vtxY, trk.getTrackSagX(), trk.getTrackSagY(), trk.getTrackSagR());
0625     Double_t vtxR = std::sqrt(vtxX*vtxX + vtxY*vtxY);
0626     Double_t vtxPhi = atan(vtxY/vtxX);
0627     if((vtxPhi < 0)&&(vtxX < 0)) vtxPhi += TMath::Pi();
0628     else if((vtxPhi > 0)&&(vtxY < 0)) vtxPhi -= TMath::Pi();
0629     trk.setHitR(0, vtxR);
0630     trk.setHitPhi(0, vtxPhi);
0631
0632     Double_t trkPhi = -1/std::tan(vtxPhi);
0633     Double_t vtxZ = 0.;
0634     if(trk.getHitIs(6))
0635     {
0636         vtxZ = trk.getHitZ(6) - trk.getHitZ(6)/std::tan(trk.getTrackTheta());
0637     }
0638     trk.setHitZ(0, vtxZ);
0639 }
0640
0641 void InttSeedTracking::ParticleIdentify(tracKuma trk)
0642 {
0643     Double_t tempInttPhi = 0.;
0644     if(trk.getHitIs(4)) tempInttPhi = trk.getHitPhi(4);
0645     else tempInttPhi = trk.getHitPhi(5);
0646     Double_t tempCalPhi = trk.getHitPhi(6);
0647
0648     if(tempInttPhi < 0) tempInttPhi += 2*TMath::Pi();
0649     if(tempCalPhi < 0)  tempCalPhi += 2*TMath::Pi();
0650     Int_t charge = (tempCalPhi - tempCalPhi)/std::abs(tempCalPhi - tempCalPhi);
0651     trk.setTrackCharge(charge);
0652
0653     bool tempElectronIs = false;
0654     if(trk.getTrackE()/trk.getTrackP() > m_EOverPElectron) tempElectronIs = true;
0655     trk.setTrackElectronIs(tempElectronIs);
0656 }
0657
0658 void InttSeedTracking::CrossLineCircle(Double_t& xv, Double_t& yv, Double_t xc, Double_t yc, Double_t R)
0659 {
0660     // line (0,0) -> (xc, yc): y = slope*x
0661     // circle: (y-yc)^2 = R^2 - (x-xc)^2
0662     // x^2 + 2bx + c = 0
0663     // -> x = -b pm sqrt(b^2 - c)
0664     Double_t slope = yc/xc;
0665
0666     Double_t b = (xc - slope*yc)/(slope*slope + 1);
0667     Double_t c = (yc*yc - R*R + xc*xc)/(slope*slope + 1);
0668
0669     Double_t xv1 = -b + std::sqrt(b*b - c);
0670     Double_t xv2 = -b - std::sqrt(b*b - c);
0671     Double_t yv1 = slope*xv1;
0672     Double_t yv2 = slope*xv2;
0673
0674     if((xv1*xv1+yv1*yv1) < (xv2*xv2+yv2*yv2)) xv = xv1;
0675     else xv = xv2;
0676
0677     yv = slope*xv;
0678
0679     return;
0680 }
0681
0682
0683
0684
0685
0686
0687
0688
0689
0690
0691
0692
0693
0694
0695
0696
0697
0698
0699
0700
0701
0702
0703
0704
0705
0706
0707
0708
0709
0710
0711
0712
0713
0714
0715
0716
0717
0718
0719
0720
0721
0722
0723
0724
0725
0726
0727
0728
0729
0730
0731
0732