From 70fc0126afc6c6b0384bb80e10fd56919add9c81 Mon Sep 17 00:00:00 2001 From: delico18 Date: Thu, 12 Sep 2024 17:21:32 +0900 Subject: [PATCH 1/6] ITS alignment monitoring (check with QCv20240903) --- Modules/ITS/include/ITS/ITSTrackTask.h | 148 +++++++ Modules/ITS/itsTrack.json | 19 +- Modules/ITS/src/ITSTrackTask.cxx | 561 +++++++++++++++++++++++++ 3 files changed, 726 insertions(+), 2 deletions(-) diff --git a/Modules/ITS/include/ITS/ITSTrackTask.h b/Modules/ITS/include/ITS/ITSTrackTask.h index 34c2d06f01..1aa0896cd1 100644 --- a/Modules/ITS/include/ITS/ITSTrackTask.h +++ b/Modules/ITS/include/ITS/ITSTrackTask.h @@ -25,6 +25,13 @@ #include #include "Common/TH1Ratio.h" #include "Common/TH2Ratio.h" +#include +#include +#include +#include +#include "TMath.h" +#include +#include "ITS/ITSHelpers.h" class TH1; class TH2; @@ -118,6 +125,147 @@ class ITSTrackTask : public TaskInterface double ptBins[141]; // pt bins o2::itsmft::TopologyDictionary* mDict; + + private: + //analysis for its-only residual + o2::its::GeometryTGeo* mGeom; + + double FitStepSize0 = 0.3; + double FitStepSize1 = 1.0e-5; + double FitStepSize2 = 1.0e-5; + double FitStepSize3 = 1.0e-5; + + double FitTolerance = 1.0e-8; + double ITS_AbsBz = 0.5; //T + + double inputG_C[3*NLayer]; + double FitparXY[6]; + double FitparDZ[2]; + ROOT::Fit::Fitter fitterA; + ROOT::Fit::Fitter fitterB; + Int_t mDefaultMomResPar = 0; + float mResMonTrackMinPt = 0; + + TH1D* hResidualRealTimePerTrackwFinerBin; + TH1D* hResidualCPUTimePerTrackwFinerBin; + std::array, NLayer> hdxySensor{};//[NLayer]; + std::array, NLayer> hdzSensor{};//[NLayer]; + + void circlefit_XY(double* input, double* par, double &MSEvalue, std::vector hitUpdate, int step = 0); + + //default setting + // function Object to be minimized + struct se_circlefitXY { + // the TGraph is a data member of the object + std::vector fHits; + double thetaR; + double Bz; + double sigma_meas[2][NLayer] = {{45,45,45,55,55,55,55}, + {40,40,40,40,40,40,40}}; //um unit + double sigma_msc[2][NLayer] = {{30,30,30,110,110,110,110}, + {25,25,25,75,75,75,75}}; //um unit + + se_circlefitXY() = default; + se_circlefitXY(std::vector &h, double tR, double bz) : fHits(h), thetaR(tR), Bz(bz) { }; + + void loadparameters(double arrpar_meas[][NLayer], double arrpar_msc[][NLayer]){ + for(int a = 0; a < 2; a++){ + for(int l = 0; l < NLayer; l++){ + sigma_meas[a][l] = arrpar_meas[a][l]; + sigma_msc[a][l] = arrpar_msc[a][l]; + } + } + }; + + void init(){ + fHits.clear(); + thetaR = 0; + Bz = 0; + }; + + void set(std::vector &h, double tR, double bz){ + fHits = h; + thetaR = tR; + Bz = bz; + }; + + double getsigma(double R, int L, double B, int axis){ + //R : cm + //B : T + if(L<0 || L>=NLayer) return 1; + double aL = sigma_meas[axis][L]*1e-4; //um -> cm + double bL = sigma_msc[axis][L]*1e-4; //um -> cm + double Beff = 0.3*B; + double sigma = std::sqrt( std::pow(aL,2) + std::pow(bL,2)/(std::pow(Beff,2)*std::pow(R*1e-2,2))); + + return sigma; + }; + + // calculate distance line-point + double distance2(double x, double y, const double *p, double tR, int charge) { + + double R = std::abs(1/p[0]); + + double Xc = p[0] > 0 ? R*std::cos(p[1] + tR + 0.5*TMath::Pi()) : R*std::cos(p[1] + tR - 0.5*TMath::Pi()); + double Yc = p[0] > 0 ? R*std::sin(p[1] + tR + 0.5*TMath::Pi()) : R*std::sin(p[1] + tR - 0.5*TMath::Pi()); + + double dx = x - (Xc + p[2]); + double dy = y - (Yc + p[3]); + double dxy = R - std::sqrt(dx*dx + dy*dy); + + double d2 = dxy*dxy; + return d2; + }; + + // implementation of the function to be minimized + double operator() (const double *par) { //const double -> double + assert(fHits != 0); + + int nhits = fHits.size(); + double sum = 0; + + double charge = par[0]>0 ? +1 : -1; + double RecRadius = std::abs(1/par[0]); + + double Sigma_tot[7]; + double sum_Sigma_tot = 0; + for(int l = 0; l < nhits; l++){ + Sigma_tot[l] = getsigma(RecRadius, fHits[l].Z(), Bz, 1); + sum_Sigma_tot += 1/(std::pow(Sigma_tot[l],2)); + } + + for (int l = 0; l < nhits; l++) { + double d = distance2(fHits[l].X(), fHits[l].Y(), par, thetaR, charge)/(std::pow(Sigma_tot[l],2)); + sum += d; + } + sum = sum / sum_Sigma_tot; + + return sum; + }; + + }; + + se_circlefitXY fitfuncXY; + + void linefit_DZ(double* zIn, double* betaIn, double* parz, double Radius, bool vertex, std::vector hitUpdate); + + double sigma_meas[2][NLayer] = {{45,45,45,55,55,55,55}, + {40,40,40,40,40,40,40}}; //um unit + double sigma_msc[2][NLayer] = {{30,30,30,110,110,110,110}, + {25,25,25,75,75,75,75}}; //um unit + + double getsigma(double R, int L, double B, int axis){ + //R : cm + //B : T + if(L<0 || L>=7) return 1; + double aL = sigma_meas[axis][L]*1e-4; //um -> cm + double bL = sigma_msc[axis][L]*1e-4; //um -> cm + double Beff = 0.3*B; + double sigma = std::sqrt( std::pow(aL,2) + std::pow(bL,2)/(std::pow(Beff,2)*std::pow(R*1e-2,2))); + + return sigma; + }; + }; } // namespace o2::quality_control_modules::its diff --git a/Modules/ITS/itsTrack.json b/Modules/ITS/itsTrack.json index fb48c9f02d..a65a115a6c 100644 --- a/Modules/ITS/itsTrack.json +++ b/Modules/ITS/itsTrack.json @@ -45,9 +45,24 @@ "nBCbins": "103", "dicttimestamp" : "0", "doNorm" : "1", - "InvMasses" : "0" + "InvMasses" : "0", + "UseDefaultMomResPar" : "0", + "MomResParMEAS1": "45, 45, 45, 55, 55, 55, 55", + "MomResParMEAS2": "40, 40, 40, 40, 40, 40, 40", + "MomResParMSC1": "30, 30, 30, 110, 110, 110, 110", + "MomResParMSC2": "25, 25, 25, 75, 75, 75, 75", + "ResidualMonitorTrackMinPt" : "0.00" + }, + "grpGeomRequest" : { + "geomRequest": "Aligned", + "askGRPECS": "false", + "askGRPLHCIF": "false", + "askGRPMagField": "false", + "askMatLUT": "false", + "askTime": "false", + "askOnceAllButField": "true", + "needPropagatorD": "false" } - } }, "checks" : { diff --git a/Modules/ITS/src/ITSTrackTask.cxx b/Modules/ITS/src/ITSTrackTask.cxx index bc8aa44ba1..378cf69f3d 100644 --- a/Modules/ITS/src/ITSTrackTask.cxx +++ b/Modules/ITS/src/ITSTrackTask.cxx @@ -67,6 +67,25 @@ void ITSTrackTask::initialize(o2::framework::InitContext& /*ctx*/) mInvMasses = o2::quality_control_modules::common::getFromConfig(mCustomParameters, "InvMasses", mInvMasses); mPublishMore = o2::quality_control_modules::common::getFromConfig(mCustomParameters, "publishMore", mPublishMore); + //analysis for its-only residual + mDefaultMomResPar = o2::quality_control_modules::common::getFromConfig(mCustomParameters, "UseDefaultMomResPar", mDefaultMomResPar); + if(mDefaultMomResPar == 0){ + std::vector vMomResParMEAS1 = convertToArray(o2::quality_control_modules::common::getFromConfig(mCustomParameters, "MomResParMEAS1", "")); + std::vector vMomResParMEAS2 = convertToArray(o2::quality_control_modules::common::getFromConfig(mCustomParameters, "MomResParMEAS2", "")); + std::vector vMomResParMSC1 = convertToArray(o2::quality_control_modules::common::getFromConfig(mCustomParameters, "MomResParMSC1", "")); + std::vector vMomResParMSC2 = convertToArray(o2::quality_control_modules::common::getFromConfig(mCustomParameters, "MomResParMSC2", "")); + + for(int l = 0; l < NLayer; l++){ + sigma_meas[0][l] = (double)vMomResParMEAS1[l]; + sigma_meas[1][l] = (double)vMomResParMEAS2[l]; + sigma_msc[0][l] = (double)vMomResParMSC1[l]; + sigma_msc[1][l] = (double)vMomResParMSC2[l]; + } + } + mResMonTrackMinPt = o2::quality_control_modules::common::getFromConfig(mCustomParameters, "ResidualMonitorTrackMinPt", mResMonTrackMinPt); + + fitfuncXY.loadparameters(sigma_meas, sigma_msc); + // pt bins definition: 20 MeV/c width up to 1 GeV/c, 100 MeV/c afterwards ptBins[0] = 0.; for (int i = 1; i < 141; i++) { @@ -92,6 +111,10 @@ void ITSTrackTask::monitorData(o2::framework::ProcessingContext& ctx) ILOG(Debug, Devel) << "START DOING QC General" << ENDM; + //o2::base::GeometryManager::loadGeometry(); + mGeom = o2::its::GeometryTGeo::Instance(); + mGeom->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::L2G)); + if (mTimestamp == -1) { // get dict from ccdb mTimestamp = std::stol(o2::quality_control_modules::common::getFromConfig(mCustomParameters, "dicttimestamp", "0")); long int ts = mTimestamp ? mTimestamp : ctx.services().get().creation; @@ -226,6 +249,13 @@ void ITSTrackTask::monitorData(o2::framework::ProcessingContext& ctx) } } nClusterCntTrack += track.getNumberOfClusters(); + + std::vector hitUpdate; + int chipIDs[7] = {-1, -1, -1, -1, -1, -1, -1}; + if(out.getPt()>mResMonTrackMinPt) { + for (int iLayer = 0; iLayer < NLayer; iLayer++) hitUpdate.push_back(false); + } + for (int icluster = 0; icluster < track.getNumberOfClusters(); icluster++) { const int index = clusIdx[track.getFirstClusterEntry() + icluster]; auto& cluster = clusArr[index]; @@ -240,6 +270,66 @@ void ITSTrackTask::monitorData(o2::framework::ProcessingContext& ctx) if (mPublishMore) { hNClusterVsChipITS->Fill(ChipID + 1, clusterSizeWithCorrection); } + + //Residual Monitoring + if(out.getPt()>mResMonTrackMinPt) { + + if(layer < 0 || layer >= NLayer) continue; + hitUpdate[layer] = true; + + auto locC = mDict->getClusterCoordinates(cluster); + auto gloC = mGeom->getMatrixL2G(ChipID) * locC; + inputG_C[3*layer + 0] = gloC.X(); + inputG_C[3*layer + 1] = gloC.Y(); + inputG_C[3*layer + 2] = gloC.Z(); + + chipIDs[layer] = ChipID; + } + } + + //Residual Monitoring + if(out.getPt()>mResMonTrackMinPt) { + for (int ipar = 0; ipar < 6; ipar++) FitparXY[ipar] = 0; + + double Cost_Fit = 0; + circlefit_XY(inputG_C, FitparXY, Cost_Fit, hitUpdate, 0); + + double RecRadius = std::abs(1/FitparXY[0]); + double CircleXc = FitparXY[0] >0 ? RecRadius*std::cos(FitparXY[1]+FitparXY[4] + 0.5*TMath::Pi()) : RecRadius*std::cos(FitparXY[1]+FitparXY[4] - 0.5*TMath::Pi()); + double CircleYc = FitparXY[0] >0 ? RecRadius*std::sin(FitparXY[1]+FitparXY[4] + 0.5*TMath::Pi()) : RecRadius*std::sin(FitparXY[1]+FitparXY[4] - 0.5*TMath::Pi()); + + CircleXc += FitparXY[2]; + CircleYc += FitparXY[3]; + + for (int ipar = 0; ipar < 2; ipar++) FitparDZ[ipar] = 0; + double z_meas[NLayer]; + double beta[NLayer]; + for(int l = 0; l < NLayer; l++) { + z_meas[l] = inputG_C[3*l + 2]; + beta[l] = std::atan2(inputG_C[3*l + 1] - CircleYc, inputG_C[3*l + 0] - CircleXc); + } + + linefit_DZ(z_meas, beta, FitparDZ, RecRadius, false, hitUpdate); + + for(int l=0; l0) sign = +1; + else sign = -1; + double dxy = sign*std::sqrt(std::pow(proj_GXc-meas_GXc,2) + std::pow(proj_GYc-meas_GYc,2)); + double dz = proj_GZc - meas_GZc; + + hdxySensor[l]->Fill(dxy,chipIDs[l]); + hdzSensor[l]->Fill(dz,chipIDs[l]); + } } // Find V0s @@ -441,6 +531,11 @@ void ITSTrackTask::reset() hTrackPtVsEta->Reset(); hTrackPtVsPhi->Reset(); + + for (int l = 0; l < NLayer; l++) { + hdxySensor[l]->Reset(); + hdzSensor[l]->Reset(); + } } void ITSTrackTask::createAllHistos() @@ -628,6 +723,21 @@ void ITSTrackTask::createAllHistos() addObject(hTrackPtVsPhi); formatAxes(hTrackPtVsPhi, "#it{p}_{T} (GeV/#it{c})", "#phi", 1, 1.10); hTrackPtVsPhi->SetStats(0); + + for (int l = 0; l < NLayer; l++) { + //sensor + hdxySensor[l] = std::make_unique(Form("hdxy%dSensor",l), Form("ChipID vs dxy, Layer %d",l), + 1000, -0.05, 0.05, ChipBoundary[l+1] - ChipBoundary[l], ChipBoundary[l], ChipBoundary[l+1]); + addObject(hdxySensor[l].get()); + formatAxes(hdxySensor[l].get(), "dxy(cm)", "chipID", 1, 1.10); + hdxySensor[l]->SetStats(0); + + hdzSensor[l] = std::make_unique(Form("hdz%dSensor",l), Form("ChipID vs dz, Layer %d",l), + 1000, -0.05, 0.05, ChipBoundary[l+1] - ChipBoundary[l], ChipBoundary[l], ChipBoundary[l+1]); + addObject(hdzSensor[l].get()); + formatAxes(hdzSensor[l].get(), "dz(cm)", "chipID", 1, 1.10); + hdzSensor[l]->SetStats(0); + } } void ITSTrackTask::addObject(TObject* aObject) @@ -648,4 +758,455 @@ void ITSTrackTask::publishHistos() } } + +void ITSTrackTask::circlefit_XY(double* input, double* par, double &MSEvalue, std::vector hitUpdate, int step){ + + int hitentries = 0; + for(int a=0; a hr; + + double frphiX = 0; + double frphiY = 0; + int nfrdet = 0; + for(int a=0; a i, j, k; //[hitentries] + std::vector irot, jrot; //[hitentries] + + int fa=0; + int index[7] = {-1, -1, -1, -1, -1, -1, -1}; + for(int a=0; a initR; + + //standard seeding + //012 + (2)34(5) + 56 + int hit1[] = {0, 1, 2}; + int hit2[] = {3, 4, 5, 2}; bool hit_mid = false; + int hit3[] = {6, 5}; + for(int i1 = 0; i1 < 3; i1++){ + // i1 -> hit1[i1] + for(int i2 = 0; i2 < 4; i2++){ + // i2 -> hit2[i2] + if(hit_mid==true && i2>= 2) continue; + for(int i3 = 0; i3 < 2; i3++){ + // i3 -> hit3[i3] + if(hit1[i1]==hit2[i2] || hit2[i2]==hit3[i3]) continue; + //if(hit_mid==true) continue; + if(hitUpdate[hit1[i1]]==false) continue; + if(hitUpdate[hit2[i2]]==false) continue; + if(hitUpdate[hit3[i3]]==false) continue; + double hitX[3] = {i[index[hit1[i1]]], i[index[hit2[i2]]], i[index[hit3[i3]]]}; + double hitY[3] = {j[index[hit1[i1]]], j[index[hit2[i2]]], j[index[hit3[i3]]]}; + + double d12 = -(hitX[1] - hitX[0])/(hitY[1] - hitY[0]); + double d23 = -(hitX[2] - hitX[1])/(hitY[2] - hitY[1]); + + double x12 = 0.5*(hitX[0] + hitX[1]); + double x23 = 0.5*(hitX[1] + hitX[2]); + double y12 = 0.5*(hitY[0] + hitY[1]); + double y23 = 0.5*(hitY[1] + hitY[2]); + + double CenterX = ((-d23*x23 + d12*x12) + (y23 - y12))/(-d23 + d12); + double CenterY = d12*(CenterX - x12) + y12; + + double temp_R = std::sqrt(std::pow(CenterX - hitX[0],2) + std::pow(CenterY - hitY[0],2)); + initR.push_back(TVector3(CenterX,CenterY,temp_R)); + if(i2<2) { + hit_mid=true; // mid hit is successfully used. Do not find inner or outer hits for initial radius searching + } + cntR[0]++; + } + } + } + + /* + { + // fast seeding + //int index[7] = {-1, -1, -1, -1, -1, -1, -1, -1}; + int hit1[] = {0, 1, 2}; + int hit2[] = {3, 4, 5, 2}; bool hit_mid = false; + int hit3[] = {6, 5}; + + int nhit[3] = {0, 0, 0}; + double hitX[3] = {0, 0, 0}; + double hitY[3] = {0, 0, 0}; + + for(int i1 = 0; i1 < 3; i1++){ + // i1 -> hit1[i1] + if(hitUpdate[hit1[i1]]==true) { + hitX[0] += i[index[hit1[i1]]]; + hitY[0] += j[index[hit1[i1]]]; + nhit[0]++; + } + } + for(int i2 = 0; i2 < 4; i2++){ + // i2 -> hit2[i2] + if(hitUpdate[hit2[i2]]==true) { + hitX[1] += i[index[hit2[i2]]]; + hitY[1] += j[index[hit2[i2]]]; + nhit[1]++; + } + } + for(int i3 = 0; i3 < 2; i3++){ + // i3 -> hit3[i3] + if(hitUpdate[hit3[i3]]==true) { + hitX[2] += i[index[hit3[i3]]]; + hitY[2] += j[index[hit3[i3]]]; + nhit[2]++; + } + } + + hitX[0] /= (double)nhit[0]; hitY[0] /= (double)nhit[0]; + hitX[1] /= (double)nhit[1]; hitY[1] /= (double)nhit[1]; + hitX[2] /= (double)nhit[2]; hitY[2] /= (double)nhit[2]; + + double d12 = -(hitX[1] - hitX[0])/(hitY[1] - hitY[0]); + double d23 = -(hitX[2] - hitX[1])/(hitY[2] - hitY[1]); + + double x12 = 0.5*(hitX[0] + hitX[1]); + double x23 = 0.5*(hitX[1] + hitX[2]); + double y12 = 0.5*(hitY[0] + hitY[1]); + double y23 = 0.5*(hitY[1] + hitY[2]); + + double CenterX = ((-d23*x23 + d12*x12) + (y23 - y12))/(-d23 + d12); + double CenterY = d12*(CenterX - x12) + y12; + + double temp_R = std::sqrt(std::pow(CenterX - hitX[0],2) + std::pow(CenterY - hitY[0],2)); + initR.push_back(TVector3(CenterX,CenterY,temp_R)); + cntR[0]++; + } + */ + + if(initR.size()==0) { + initR.push_back(TVector3(0,0,10000)); + cntR[0]++; + } + + double mean_X[] = {0, 0}; + double mean_Y[] = {0, 0}; + double mean_R[] = {0, 0}; + for(int i = 0; i < cntR[0]; i++) { + mean_X[0] += initR[i].X()/(double)cntR[0]; + mean_Y[0] += initR[i].Y()/(double)cntR[0]; + mean_R[0] += initR[i](2)/(double)cntR[0]; + } + + for(int i = 0; i < cntR[0]; i++) { + if(std::abs(mean_R[0] - initR[i](2)) < mean_R[0]) { + mean_X[1] += initR[i].X(); + mean_Y[1] += initR[i].Y(); + mean_R[1] += initR[i](2); + cntR[1]++; + } + } + mean_R[1] /= cntR[1]; + + mean_R[1] *= std::pow(sqrt(10),step); + if(mean_R[1]<1.0e+1) mean_R[1] = 1.0e+1; + if(mean_R[1]>1.0e+6) mean_R[1] = 1.0e+6; + + double thetaR = std::atan2( jrot[0], irot[0]); + + double temp_parA[4]; + temp_parA[0] = + 1/mean_R[1]; + temp_parA[1] = 0; + + double temp_parB[4]; + temp_parB[0] = - 1/mean_R[1]; + temp_parB[1] = 0; + + // make the functor object + fitfuncXY.init(); + fitfuncXY.set(hr, thetaR, ITS_AbsBz); + ROOT::Math::Functor fcn(fitfuncXY, 4); + + double pStartA[4] = {temp_parA[0],temp_parA[1], 0, 0}; + fitterA.SetFCN(fcn,pStartA); + fitterA.Config().ParSettings(0).SetStepSize(FitStepSize0); + fitterA.Config().ParSettings(1).SetStepSize(FitStepSize1); + fitterA.Config().ParSettings(2).SetStepSize(FitStepSize2); + fitterA.Config().ParSettings(3).SetStepSize(FitStepSize3); + fitterA.Config().ParSettings(0).SetLimits(+1.0e-10, +1.0e-1); // + side + + double pStartB[4] = {temp_parB[0],temp_parB[1], 0, 0}; + fitterB.SetFCN(fcn,pStartB); + fitterB.Config().ParSettings(0).SetStepSize(FitStepSize0); + fitterB.Config().ParSettings(1).SetStepSize(FitStepSize1); + fitterB.Config().ParSettings(2).SetStepSize(FitStepSize2); + fitterB.Config().ParSettings(3).SetStepSize(FitStepSize3); + fitterB.Config().ParSettings(0).SetLimits(-1.0e-1, -1.0e-10); // - side + + ROOT::Math::MinimizerOptions minOpt; + for(int iTol = 0; iTol < 4; iTol++){ + + minOpt.SetTolerance(std::pow(10,iTol)*FitTolerance); + + fitterA.Config().SetMinimizerOptions(minOpt); + fitterB.Config().SetMinimizerOptions(minOpt); + bool okA = fitterA.FitFCN(); + bool okB = fitterB.FitFCN(); + + if (!okA) { + if(!okB) { + const ROOT::Fit::FitResult & resultA = fitterA.Result(); + const double * parFitA = resultA.GetParams(); + double MSEvalueA = resultA.MinFcnValue(); + const ROOT::Fit::FitResult & resultB = fitterB.Result(); + const double * parFitB = resultB.GetParams(); + double MSEvalueB = resultB.MinFcnValue(); + + TMatrixD vxyA[2]; + vxyA[0].ResizeTo(1,2); + vxyA[0][0] = { parFitA[2], parFitA[3]}; + vxyA[0].T(); + vxyA[1].ResizeTo(2,1); + vxyA[1] = RotFInv * vxyA[0]; + vxyA[1].T(); + + TMatrixD vxyB[2]; + vxyB[0].ResizeTo(1,2); + vxyB[0][0] = { parFitB[2], parFitB[3]}; + vxyB[0].T(); + vxyB[1].ResizeTo(2,1); + vxyB[1] = RotFInv * vxyB[0]; + vxyB[1].T(); + + if(MSEvalueA hitUpdate){ + + int hitentries = 0; + for(int a=0; a z, beta; + std::vector index; + + int fa=0; + for(int a=0; a Date: Thu, 19 Sep 2024 15:21:26 +0900 Subject: [PATCH 2/6] ITS alignment monitoring (check with QCv20240903), revised 20240919(5) --- Modules/ITS/include/ITS/ITSTrackTask.h | 6 +- Modules/ITS/include/ITS/README.md | 31 +++++++++ Modules/ITS/itsTrack.json | 12 +--- Modules/ITS/src/ITSTrackTask.cxx | 92 ++++++-------------------- 4 files changed, 56 insertions(+), 85 deletions(-) create mode 100644 Modules/ITS/include/ITS/README.md diff --git a/Modules/ITS/include/ITS/ITSTrackTask.h b/Modules/ITS/include/ITS/ITSTrackTask.h index 1aa0896cd1..39ccf26559 100644 --- a/Modules/ITS/include/ITS/ITSTrackTask.h +++ b/Modules/ITS/include/ITS/ITSTrackTask.h @@ -143,7 +143,9 @@ class ITSTrackTask : public TaskInterface double FitparDZ[2]; ROOT::Fit::Fitter fitterA; ROOT::Fit::Fitter fitterB; + Int_t mAlignmentMonitor = 0; Int_t mDefaultMomResPar = 0; + Int_t mResMonNclMin = 0; float mResMonTrackMinPt = 0; TH1D* hResidualRealTimePerTrackwFinerBin; @@ -151,7 +153,7 @@ class ITSTrackTask : public TaskInterface std::array, NLayer> hdxySensor{};//[NLayer]; std::array, NLayer> hdzSensor{};//[NLayer]; - void circlefit_XY(double* input, double* par, double &MSEvalue, std::vector hitUpdate, int step = 0); + void circleFitXY(double* input, double* par, double &MSEvalue, std::vector hitUpdate, int step = 0); //default setting // function Object to be minimized @@ -247,7 +249,7 @@ class ITSTrackTask : public TaskInterface se_circlefitXY fitfuncXY; - void linefit_DZ(double* zIn, double* betaIn, double* parz, double Radius, bool vertex, std::vector hitUpdate); + void lineFitDZ(double* zIn, double* betaIn, double* parz, double Radius, bool vertex, std::vector hitUpdate); double sigma_meas[2][NLayer] = {{45,45,45,55,55,55,55}, {40,40,40,40,40,40,40}}; //um unit diff --git a/Modules/ITS/include/ITS/README.md b/Modules/ITS/include/ITS/README.md new file mode 100644 index 0000000000..ea3f3ce79c --- /dev/null +++ b/Modules/ITS/include/ITS/README.md @@ -0,0 +1,31 @@ +# Analysis for ITS-only residual + +We use the Weighted Least Square(WLS) fit and calculate the residuals from each cluster positions to the fitted trajectory. + +## Circular Fit in the XY plane : circleFitXY + +In the XY plane, the trajectory of particle under the magnetic field in z-axis can be fitted with circular model. +There are 4 fit parameters : radius(R), direction at origin($$\theta_R$$), estimated vertex X(vx), and estimated vertex Y(vy). + +$$ x_c = R \cdot cos(\theta_R + \frac{R}{|R|} \cdot \pi/2) $$ + +$$ y_c = R \cdot sin(\theta_R + \frac{R}{|R|} \cdot \pi/2) $$ + +$$ (x - (x_c + v_x))^2 + (y - (y_c + v_y))^2 = |R|^2 $$ + +For the fit stability, +1. A simple algorithm to determine seeds of fit parameters(radius and direction at origin) is applied. +2. Circular fit is performed for both (+) and (-) signs of the circular trajectory (even if we have approximated seeds from 1.), the one with the lower chi2(sum of weighted residuals) is finally adopted. This procedure ensures the stability in case of hard(high pT) tracks. + +## Linear Fit in the DZ plane : lineFitDZ + +After 4 fit parameters of circular trajectory is determined, we can describe the fit positions of associated clusters as well as the global positions. +This makes it possible to represent each cluster position in the poloar coordinates $$(r, \beta)$$, whose origin is $$(x_c, y_c)$$. + +One of the simplest residual in z direction can be calculated in the RZ plane, where $$r = \sqrt(x^2 + y^2)$$. +However, in this case, the fitting stability cannot be guaranteed if the collision position is not located in the same quadrant. + +On the other hand, another parameter beta determined above the equation of the circle determined (in the XY plane) corresponds to the distance of the particle moving(=D) in the XY plane. +In this case, the fitting stability can be guaranteed even if the collision position is not located in the same quadrant. + +So we proceed with the linear fit in the DZ plane. diff --git a/Modules/ITS/itsTrack.json b/Modules/ITS/itsTrack.json index a65a115a6c..969ab6b4eb 100644 --- a/Modules/ITS/itsTrack.json +++ b/Modules/ITS/itsTrack.json @@ -46,22 +46,14 @@ "dicttimestamp" : "0", "doNorm" : "1", "InvMasses" : "0", + "doAlignmentMonitor" : "1", "UseDefaultMomResPar" : "0", "MomResParMEAS1": "45, 45, 45, 55, 55, 55, 55", "MomResParMEAS2": "40, 40, 40, 40, 40, 40, 40", "MomResParMSC1": "30, 30, 30, 110, 110, 110, 110", "MomResParMSC2": "25, 25, 25, 75, 75, 75, 75", + "ResidualMonitorNclMin" : "5", "ResidualMonitorTrackMinPt" : "0.00" - }, - "grpGeomRequest" : { - "geomRequest": "Aligned", - "askGRPECS": "false", - "askGRPLHCIF": "false", - "askGRPMagField": "false", - "askMatLUT": "false", - "askTime": "false", - "askOnceAllButField": "true", - "needPropagatorD": "false" } } }, diff --git a/Modules/ITS/src/ITSTrackTask.cxx b/Modules/ITS/src/ITSTrackTask.cxx index 378cf69f3d..6532e426a8 100644 --- a/Modules/ITS/src/ITSTrackTask.cxx +++ b/Modules/ITS/src/ITSTrackTask.cxx @@ -68,8 +68,9 @@ void ITSTrackTask::initialize(o2::framework::InitContext& /*ctx*/) mPublishMore = o2::quality_control_modules::common::getFromConfig(mCustomParameters, "publishMore", mPublishMore); //analysis for its-only residual + mAlignmentMonitor = o2::quality_control_modules::common::getFromConfig(mCustomParameters, "doAlignmentMonitor", mAlignmentMonitor); mDefaultMomResPar = o2::quality_control_modules::common::getFromConfig(mCustomParameters, "UseDefaultMomResPar", mDefaultMomResPar); - if(mDefaultMomResPar == 0){ + if(mAlignmentMonitor == 1 && mDefaultMomResPar == 0){ std::vector vMomResParMEAS1 = convertToArray(o2::quality_control_modules::common::getFromConfig(mCustomParameters, "MomResParMEAS1", "")); std::vector vMomResParMEAS2 = convertToArray(o2::quality_control_modules::common::getFromConfig(mCustomParameters, "MomResParMEAS2", "")); std::vector vMomResParMSC1 = convertToArray(o2::quality_control_modules::common::getFromConfig(mCustomParameters, "MomResParMSC1", "")); @@ -82,6 +83,7 @@ void ITSTrackTask::initialize(o2::framework::InitContext& /*ctx*/) sigma_msc[1][l] = (double)vMomResParMSC2[l]; } } + mResMonNclMin = o2::quality_control_modules::common::getFromConfig(mCustomParameters, "ResidualMonitorNclMin", mResMonNclMin); mResMonTrackMinPt = o2::quality_control_modules::common::getFromConfig(mCustomParameters, "ResidualMonitorTrackMinPt", mResMonTrackMinPt); fitfuncXY.loadparameters(sigma_meas, sigma_msc); @@ -111,10 +113,6 @@ void ITSTrackTask::monitorData(o2::framework::ProcessingContext& ctx) ILOG(Debug, Devel) << "START DOING QC General" << ENDM; - //o2::base::GeometryManager::loadGeometry(); - mGeom = o2::its::GeometryTGeo::Instance(); - mGeom->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::L2G)); - if (mTimestamp == -1) { // get dict from ccdb mTimestamp = std::stol(o2::quality_control_modules::common::getFromConfig(mCustomParameters, "dicttimestamp", "0")); long int ts = mTimestamp ? mTimestamp : ctx.services().get().creation; @@ -123,6 +121,11 @@ void ITSTrackTask::monitorData(o2::framework::ProcessingContext& ctx) std::map metadata; mDict = TaskInterface::retrieveConditionAny("ITS/Calib/ClusterDictionary", metadata, ts); ILOG(Debug, Devel) << "Dictionary size: " << mDict->getSize() << ENDM; + + o2::its::GeometryTGeo::adopt(TaskInterface::retrieveConditionAny("ITS/Config/Geometry", metadata, ts)); + mGeom = o2::its::GeometryTGeo::Instance(); + //mGeom->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::L2G)); + ILOG(Debug, Devel) << "Loaded new instance of mGeom" << ENDM; } auto trackArr = ctx.inputs().get>("tracks"); @@ -252,7 +255,7 @@ void ITSTrackTask::monitorData(o2::framework::ProcessingContext& ctx) std::vector hitUpdate; int chipIDs[7] = {-1, -1, -1, -1, -1, -1, -1}; - if(out.getPt()>mResMonTrackMinPt) { + if(mAlignmentMonitor == 1 && out.getPt()>mResMonTrackMinPt && track.getNumberOfClusters() >= mResMonNclMin) { for (int iLayer = 0; iLayer < NLayer; iLayer++) hitUpdate.push_back(false); } @@ -262,8 +265,9 @@ void ITSTrackTask::monitorData(o2::framework::ProcessingContext& ctx) auto ChipID = cluster.getSensorID(); int layer = 0; - while (ChipID > ChipBoundary[layer]) + while (ChipID >= ChipBoundary[layer] && layer <= NLayer){ layer++; + } layer--; double clusterSizeWithCorrection = (double)clSize[index] * (std::cos(std::atan(out.getTgl()))); @@ -272,7 +276,7 @@ void ITSTrackTask::monitorData(o2::framework::ProcessingContext& ctx) } //Residual Monitoring - if(out.getPt()>mResMonTrackMinPt) { + if(mAlignmentMonitor == 1 && out.getPt()>mResMonTrackMinPt && track.getNumberOfClusters() >= mResMonNclMin) { if(layer < 0 || layer >= NLayer) continue; hitUpdate[layer] = true; @@ -288,11 +292,11 @@ void ITSTrackTask::monitorData(o2::framework::ProcessingContext& ctx) } //Residual Monitoring - if(out.getPt()>mResMonTrackMinPt) { + if(mAlignmentMonitor == 1 && out.getPt()>mResMonTrackMinPt && track.getNumberOfClusters() >= mResMonNclMin) { for (int ipar = 0; ipar < 6; ipar++) FitparXY[ipar] = 0; double Cost_Fit = 0; - circlefit_XY(inputG_C, FitparXY, Cost_Fit, hitUpdate, 0); + circleFitXY(inputG_C, FitparXY, Cost_Fit, hitUpdate, 0); double RecRadius = std::abs(1/FitparXY[0]); double CircleXc = FitparXY[0] >0 ? RecRadius*std::cos(FitparXY[1]+FitparXY[4] + 0.5*TMath::Pi()) : RecRadius*std::cos(FitparXY[1]+FitparXY[4] - 0.5*TMath::Pi()); @@ -309,7 +313,7 @@ void ITSTrackTask::monitorData(o2::framework::ProcessingContext& ctx) beta[l] = std::atan2(inputG_C[3*l + 1] - CircleYc, inputG_C[3*l + 0] - CircleXc); } - linefit_DZ(z_meas, beta, FitparDZ, RecRadius, false, hitUpdate); + lineFitDZ(z_meas, beta, FitparDZ, RecRadius, false, hitUpdate); for(int l=0; l(Form("hdxy%dSensor",l), Form("ChipID vs dxy, Layer %d",l), - 1000, -0.05, 0.05, ChipBoundary[l+1] - ChipBoundary[l], ChipBoundary[l], ChipBoundary[l+1]); + 500, -0.05, 0.05, ChipBoundary[l+1] - ChipBoundary[l], ChipBoundary[l] - 0.5, ChipBoundary[l+1] - 0.5); addObject(hdxySensor[l].get()); formatAxes(hdxySensor[l].get(), "dxy(cm)", "chipID", 1, 1.10); hdxySensor[l]->SetStats(0); hdzSensor[l] = std::make_unique(Form("hdz%dSensor",l), Form("ChipID vs dz, Layer %d",l), - 1000, -0.05, 0.05, ChipBoundary[l+1] - ChipBoundary[l], ChipBoundary[l], ChipBoundary[l+1]); + 500, -0.05, 0.05, ChipBoundary[l+1] - ChipBoundary[l], ChipBoundary[l] - 0.5, ChipBoundary[l+1] - 0.5); addObject(hdzSensor[l].get()); formatAxes(hdzSensor[l].get(), "dz(cm)", "chipID", 1, 1.10); hdzSensor[l]->SetStats(0); @@ -759,7 +763,7 @@ void ITSTrackTask::publishHistos() } -void ITSTrackTask::circlefit_XY(double* input, double* par, double &MSEvalue, std::vector hitUpdate, int step){ +void ITSTrackTask::circleFitXY(double* input, double* par, double &MSEvalue, std::vector hitUpdate, int step){ int hitentries = 0; for(int a=0; a hit1[i1] - if(hitUpdate[hit1[i1]]==true) { - hitX[0] += i[index[hit1[i1]]]; - hitY[0] += j[index[hit1[i1]]]; - nhit[0]++; - } - } - for(int i2 = 0; i2 < 4; i2++){ - // i2 -> hit2[i2] - if(hitUpdate[hit2[i2]]==true) { - hitX[1] += i[index[hit2[i2]]]; - hitY[1] += j[index[hit2[i2]]]; - nhit[1]++; - } - } - for(int i3 = 0; i3 < 2; i3++){ - // i3 -> hit3[i3] - if(hitUpdate[hit3[i3]]==true) { - hitX[2] += i[index[hit3[i3]]]; - hitY[2] += j[index[hit3[i3]]]; - nhit[2]++; - } - } - - hitX[0] /= (double)nhit[0]; hitY[0] /= (double)nhit[0]; - hitX[1] /= (double)nhit[1]; hitY[1] /= (double)nhit[1]; - hitX[2] /= (double)nhit[2]; hitY[2] /= (double)nhit[2]; - - double d12 = -(hitX[1] - hitX[0])/(hitY[1] - hitY[0]); - double d23 = -(hitX[2] - hitX[1])/(hitY[2] - hitY[1]); - - double x12 = 0.5*(hitX[0] + hitX[1]); - double x23 = 0.5*(hitX[1] + hitX[2]); - double y12 = 0.5*(hitY[0] + hitY[1]); - double y23 = 0.5*(hitY[1] + hitY[2]); - - double CenterX = ((-d23*x23 + d12*x12) + (y23 - y12))/(-d23 + d12); - double CenterY = d12*(CenterX - x12) + y12; - - double temp_R = std::sqrt(std::pow(CenterX - hitX[0],2) + std::pow(CenterY - hitY[0],2)); - initR.push_back(TVector3(CenterX,CenterY,temp_R)); - cntR[0]++; - } - */ if(initR.size()==0) { initR.push_back(TVector3(0,0,10000)); @@ -1121,7 +1067,7 @@ void ITSTrackTask::circlefit_XY(double* input, double* par, double &MSEvalue, st } } -void ITSTrackTask::linefit_DZ(double* zIn, double* betaIn, double* parz, double Radius, bool vertex, std::vector hitUpdate){ +void ITSTrackTask::lineFitDZ(double* zIn, double* betaIn, double* parz, double Radius, bool vertex, std::vector hitUpdate){ int hitentries = 0; for(int a=0; a Date: Thu, 19 Sep 2024 15:27:31 +0900 Subject: [PATCH 3/6] ITS alignment monitoring (check with QCv20240903), revised 20240919(6) --- Modules/ITS/include/ITS/ITSTrackTask.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Modules/ITS/include/ITS/ITSTrackTask.h b/Modules/ITS/include/ITS/ITSTrackTask.h index 39ccf26559..c96026deff 100644 --- a/Modules/ITS/include/ITS/ITSTrackTask.h +++ b/Modules/ITS/include/ITS/ITSTrackTask.h @@ -229,7 +229,7 @@ class ITSTrackTask : public TaskInterface double charge = par[0]>0 ? +1 : -1; double RecRadius = std::abs(1/par[0]); - double Sigma_tot[7]; + double Sigma_tot[NLayer]; double sum_Sigma_tot = 0; for(int l = 0; l < nhits; l++){ Sigma_tot[l] = getsigma(RecRadius, fHits[l].Z(), Bz, 1); @@ -259,7 +259,7 @@ class ITSTrackTask : public TaskInterface double getsigma(double R, int L, double B, int axis){ //R : cm //B : T - if(L<0 || L>=7) return 1; + if(L<0 || L>=NLayer) return 1; double aL = sigma_meas[axis][L]*1e-4; //um -> cm double bL = sigma_msc[axis][L]*1e-4; //um -> cm double Beff = 0.3*B; From a56d41e254b3752c0dc20254b4e8a6cd40cf5669 Mon Sep 17 00:00:00 2001 From: delico18 Date: Mon, 21 Oct 2024 14:38:59 +0900 Subject: [PATCH 4/6] ITS alignment monitoring (check with QCv20240903), revised 20241021(1) --- Modules/ITS/include/ITS/ITSTrackTask.h | 23 ++--- Modules/ITS/src/ITSTrackTask.cxx | 127 +++++++++++++++---------- 2 files changed, 87 insertions(+), 63 deletions(-) diff --git a/Modules/ITS/include/ITS/ITSTrackTask.h b/Modules/ITS/include/ITS/ITSTrackTask.h index c96026deff..81a912029a 100644 --- a/Modules/ITS/include/ITS/ITSTrackTask.h +++ b/Modules/ITS/include/ITS/ITSTrackTask.h @@ -130,10 +130,7 @@ class ITSTrackTask : public TaskInterface //analysis for its-only residual o2::its::GeometryTGeo* mGeom; - double FitStepSize0 = 0.3; - double FitStepSize1 = 1.0e-5; - double FitStepSize2 = 1.0e-5; - double FitStepSize3 = 1.0e-5; + std::vector FitStepSize{0.3, 1.0e-5, 1.0e-5, 1.0e-5}; double FitTolerance = 1.0e-8; double ITS_AbsBz = 0.5; //T @@ -148,10 +145,8 @@ class ITSTrackTask : public TaskInterface Int_t mResMonNclMin = 0; float mResMonTrackMinPt = 0; - TH1D* hResidualRealTimePerTrackwFinerBin; - TH1D* hResidualCPUTimePerTrackwFinerBin; - std::array, NLayer> hdxySensor{};//[NLayer]; - std::array, NLayer> hdzSensor{};//[NLayer]; + std::array, NLayer> hResidualXY{};//[NLayer]; + std::array, NLayer> hResidualZD{};//[NLayer]; void circleFitXY(double* input, double* par, double &MSEvalue, std::vector hitUpdate, int step = 0); @@ -197,7 +192,7 @@ class ITSTrackTask : public TaskInterface if(L<0 || L>=NLayer) return 1; double aL = sigma_meas[axis][L]*1e-4; //um -> cm double bL = sigma_msc[axis][L]*1e-4; //um -> cm - double Beff = 0.3*B; + double Beff = 0.299792458*B; double sigma = std::sqrt( std::pow(aL,2) + std::pow(bL,2)/(std::pow(Beff,2)*std::pow(R*1e-2,2))); return sigma; @@ -251,18 +246,18 @@ class ITSTrackTask : public TaskInterface void lineFitDZ(double* zIn, double* betaIn, double* parz, double Radius, bool vertex, std::vector hitUpdate); - double sigma_meas[2][NLayer] = {{45,45,45,55,55,55,55}, + double mSigmaMeas[2][NLayer] = {{45,45,45,55,55,55,55}, {40,40,40,40,40,40,40}}; //um unit - double sigma_msc[2][NLayer] = {{30,30,30,110,110,110,110}, + double mSigmaMsc[2][NLayer] = {{30,30,30,110,110,110,110}, {25,25,25,75,75,75,75}}; //um unit double getsigma(double R, int L, double B, int axis){ //R : cm //B : T if(L<0 || L>=NLayer) return 1; - double aL = sigma_meas[axis][L]*1e-4; //um -> cm - double bL = sigma_msc[axis][L]*1e-4; //um -> cm - double Beff = 0.3*B; + double aL = mSigmaMeas[axis][L]*1e-4; //um -> cm + double bL = mSigmaMsc[axis][L]*1e-4; //um -> cm + double Beff = 0.299792458*B; double sigma = std::sqrt( std::pow(aL,2) + std::pow(bL,2)/(std::pow(Beff,2)*std::pow(R*1e-2,2))); return sigma; diff --git a/Modules/ITS/src/ITSTrackTask.cxx b/Modules/ITS/src/ITSTrackTask.cxx index 6532e426a8..91be439bcc 100644 --- a/Modules/ITS/src/ITSTrackTask.cxx +++ b/Modules/ITS/src/ITSTrackTask.cxx @@ -77,16 +77,16 @@ void ITSTrackTask::initialize(o2::framework::InitContext& /*ctx*/) std::vector vMomResParMSC2 = convertToArray(o2::quality_control_modules::common::getFromConfig(mCustomParameters, "MomResParMSC2", "")); for(int l = 0; l < NLayer; l++){ - sigma_meas[0][l] = (double)vMomResParMEAS1[l]; - sigma_meas[1][l] = (double)vMomResParMEAS2[l]; - sigma_msc[0][l] = (double)vMomResParMSC1[l]; - sigma_msc[1][l] = (double)vMomResParMSC2[l]; + mSigmaMeas[0][l] = (double)vMomResParMEAS1[l]; + mSigmaMeas[1][l] = (double)vMomResParMEAS2[l]; + mSigmaMsc[0][l] = (double)vMomResParMSC1[l]; + mSigmaMsc[1][l] = (double)vMomResParMSC2[l]; } } mResMonNclMin = o2::quality_control_modules::common::getFromConfig(mCustomParameters, "ResidualMonitorNclMin", mResMonNclMin); mResMonTrackMinPt = o2::quality_control_modules::common::getFromConfig(mCustomParameters, "ResidualMonitorTrackMinPt", mResMonTrackMinPt); - fitfuncXY.loadparameters(sigma_meas, sigma_msc); + fitfuncXY.loadparameters(mSigmaMeas, mSigmaMsc); // pt bins definition: 20 MeV/c width up to 1 GeV/c, 100 MeV/c afterwards ptBins[0] = 0.; @@ -122,10 +122,15 @@ void ITSTrackTask::monitorData(o2::framework::ProcessingContext& ctx) mDict = TaskInterface::retrieveConditionAny("ITS/Calib/ClusterDictionary", metadata, ts); ILOG(Debug, Devel) << "Dictionary size: " << mDict->getSize() << ENDM; - o2::its::GeometryTGeo::adopt(TaskInterface::retrieveConditionAny("ITS/Config/Geometry", metadata, ts)); - mGeom = o2::its::GeometryTGeo::Instance(); - //mGeom->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::L2G)); - ILOG(Debug, Devel) << "Loaded new instance of mGeom" << ENDM; + if(mAlignmentMonitor == 1) { + o2::its::GeometryTGeo::adopt(TaskInterface::retrieveConditionAny("ITS/Config/Geometry", metadata, ts)); + mGeom = o2::its::GeometryTGeo::Instance(); + if (!mGeom) { + ILOG(Fatal, Support) << "Can't retrive ITS geometry from ccdb - timestamp: " << ts << ENDM; + throw std::runtime_error("Can't retrive ITS geometry from ccdb!"); + } + ILOG(Debug, Devel) << "Loaded new instance of mGeom (for ITS alignment monitoring)" << ENDM; + } } auto trackArr = ctx.inputs().get>("tracks"); @@ -263,6 +268,7 @@ void ITSTrackTask::monitorData(o2::framework::ProcessingContext& ctx) const int index = clusIdx[track.getFirstClusterEntry() + icluster]; auto& cluster = clusArr[index]; auto ChipID = cluster.getSensorID(); + int ClusterID = cluster.getPatternID(); // used for normal (frequent) cluster shapes int layer = 0; while (ChipID >= ChipBoundary[layer] && layer <= NLayer){ @@ -270,8 +276,8 @@ void ITSTrackTask::monitorData(o2::framework::ProcessingContext& ctx) } layer--; - double clusterSizeWithCorrection = (double)clSize[index] * (std::cos(std::atan(out.getTgl()))); if (mPublishMore) { + double clusterSizeWithCorrection = (double)clSize[index] * (std::cos(std::atan(out.getTgl()))); hNClusterVsChipITS->Fill(ChipID + 1, clusterSizeWithCorrection); } @@ -279,9 +285,20 @@ void ITSTrackTask::monitorData(o2::framework::ProcessingContext& ctx) if(mAlignmentMonitor == 1 && out.getPt()>mResMonTrackMinPt && track.getNumberOfClusters() >= mResMonNclMin) { if(layer < 0 || layer >= NLayer) continue; - hitUpdate[layer] = true; - auto locC = mDict->getClusterCoordinates(cluster); + o2::math_utils::Point3D locC; // local coordinates + if (ClusterID != o2::itsmft::CompCluster::InvalidPatternID) { // Normal (frequent) cluster shapes + if (!mDict->isGroup(ClusterID)) { + locC = mDict->getClusterCoordinates(cluster); + } else { + o2::itsmft::ClusterPattern patt(pattIt); + locC = mDict->getClusterCoordinates(cluster, patt, true); + } + } else { // invalid pattern + continue; + } + + hitUpdate[layer] = true; auto gloC = mGeom->getMatrixL2G(ChipID) * locC; inputG_C[3*layer + 0] = gloC.X(); inputG_C[3*layer + 1] = gloC.Y(); @@ -293,6 +310,12 @@ void ITSTrackTask::monitorData(o2::framework::ProcessingContext& ctx) //Residual Monitoring if(mAlignmentMonitor == 1 && out.getPt()>mResMonTrackMinPt && track.getNumberOfClusters() >= mResMonNclMin) { + int NclValid = 0; + for (int iLayer = 0; iLayer < NLayer; iLayer++) { + if(hitUpdate[iLayer]) NclValid++; + } + if(NclValidFill(dxy,chipIDs[l]); - hdzSensor[l]->Fill(dz,chipIDs[l]); + hResidualXY[iLayer]->Fill(dxy,chipIDs[iLayer]); + hResidualZD[iLayer]->Fill(dz,chipIDs[iLayer]); } } @@ -536,9 +559,11 @@ void ITSTrackTask::reset() hTrackPtVsEta->Reset(); hTrackPtVsPhi->Reset(); - for (int l = 0; l < NLayer; l++) { - hdxySensor[l]->Reset(); - hdzSensor[l]->Reset(); + if(mAlignmentMonitor == 1) { + for (int l = 0; l < NLayer; l++) { + hResidualXY[l]->Reset(); + hResidualZD[l]->Reset(); + } } } @@ -728,19 +753,23 @@ void ITSTrackTask::createAllHistos() formatAxes(hTrackPtVsPhi, "#it{p}_{T} (GeV/#it{c})", "#phi", 1, 1.10); hTrackPtVsPhi->SetStats(0); - for (int l = 0; l < NLayer; l++) { - //sensor - hdxySensor[l] = std::make_unique(Form("hdxy%dSensor",l), Form("ChipID vs dxy, Layer %d",l), - 500, -0.05, 0.05, ChipBoundary[l+1] - ChipBoundary[l], ChipBoundary[l] - 0.5, ChipBoundary[l+1] - 0.5); - addObject(hdxySensor[l].get()); - formatAxes(hdxySensor[l].get(), "dxy(cm)", "chipID", 1, 1.10); - hdxySensor[l]->SetStats(0); - - hdzSensor[l] = std::make_unique(Form("hdz%dSensor",l), Form("ChipID vs dz, Layer %d",l), - 500, -0.05, 0.05, ChipBoundary[l+1] - ChipBoundary[l], ChipBoundary[l] - 0.5, ChipBoundary[l+1] - 0.5); - addObject(hdzSensor[l].get()); - formatAxes(hdzSensor[l].get(), "dz(cm)", "chipID", 1, 1.10); - hdzSensor[l]->SetStats(0); + if(mAlignmentMonitor == 1) { + for (int l = 0; l < NLayer; l++) { + //sensor + int NBinsChipID = (l < NLayerIB) ? (ChipBoundary[l+1] - ChipBoundary[l]) : (ChipBoundary[l+1] - ChipBoundary[l])/14; + + hResidualXY[l] = std::make_unique(Form("hResidualXY%d",l), Form("ChipID vs dxy, Layer %d",l), + 500, -0.05, 0.05, NBinsChipID, ChipBoundary[l] - 0.5, ChipBoundary[l+1] - 0.5); + addObject(hResidualXY[l].get()); + formatAxes(hResidualXY[l].get(), "dxy(cm)", (l < NLayerIB) ? "ChipID(Sensor Unit)" : "ChipID(HIC Unit)", 1, 1.10); + hResidualXY[l]->SetStats(0); + + hResidualZD[l] = std::make_unique(Form("hResidualZD%d",l), Form("ChipID vs dz, Layer %d",l), + 500, -0.05, 0.05, NBinsChipID, ChipBoundary[l] - 0.5, ChipBoundary[l+1] - 0.5); + addObject(hResidualZD[l].get()); + formatAxes(hResidualZD[l].get(), "dz(cm)", (l < NLayerIB) ? "ChipID(Sensor Unit)" : "ChipID(HIC Unit)", 1, 1.10); + hResidualZD[l]->SetStats(0); + } } } @@ -910,18 +939,18 @@ void ITSTrackTask::circleFitXY(double* input, double* par, double &MSEvalue, std double pStartA[4] = {temp_parA[0],temp_parA[1], 0, 0}; fitterA.SetFCN(fcn,pStartA); - fitterA.Config().ParSettings(0).SetStepSize(FitStepSize0); - fitterA.Config().ParSettings(1).SetStepSize(FitStepSize1); - fitterA.Config().ParSettings(2).SetStepSize(FitStepSize2); - fitterA.Config().ParSettings(3).SetStepSize(FitStepSize3); + fitterA.Config().ParSettings(0).SetStepSize((float)FitStepSize[0]); + fitterA.Config().ParSettings(1).SetStepSize((float)FitStepSize[1]); + fitterA.Config().ParSettings(2).SetStepSize((float)FitStepSize[2]); + fitterA.Config().ParSettings(3).SetStepSize((float)FitStepSize[3]); fitterA.Config().ParSettings(0).SetLimits(+1.0e-10, +1.0e-1); // + side double pStartB[4] = {temp_parB[0],temp_parB[1], 0, 0}; fitterB.SetFCN(fcn,pStartB); - fitterB.Config().ParSettings(0).SetStepSize(FitStepSize0); - fitterB.Config().ParSettings(1).SetStepSize(FitStepSize1); - fitterB.Config().ParSettings(2).SetStepSize(FitStepSize2); - fitterB.Config().ParSettings(3).SetStepSize(FitStepSize3); + fitterB.Config().ParSettings(0).SetStepSize((float)FitStepSize[0]); + fitterB.Config().ParSettings(1).SetStepSize((float)FitStepSize[1]); + fitterB.Config().ParSettings(2).SetStepSize((float)FitStepSize[2]); + fitterB.Config().ParSettings(3).SetStepSize((float)FitStepSize[3]); fitterB.Config().ParSettings(0).SetLimits(-1.0e-1, -1.0e-10); // - side ROOT::Math::MinimizerOptions minOpt; From bf91c26eb2b00a5e470ccf352448df4fec63e1af Mon Sep 17 00:00:00 2001 From: delico18 Date: Thu, 24 Oct 2024 10:54:41 +0900 Subject: [PATCH 5/6] ITS alignment monitoring (check with QCv20240903), revised 20241024(1), clang-format applied, TH2I->TH2F --- Modules/ITS/include/ITS/ITSTrackTask.h | 133 ++--- Modules/ITS/src/ITSTrackTask.cxx | 653 +++++++++++++------------ 2 files changed, 408 insertions(+), 378 deletions(-) diff --git a/Modules/ITS/include/ITS/ITSTrackTask.h b/Modules/ITS/include/ITS/ITSTrackTask.h index 81a912029a..9841c5194b 100644 --- a/Modules/ITS/include/ITS/ITSTrackTask.h +++ b/Modules/ITS/include/ITS/ITSTrackTask.h @@ -127,15 +127,15 @@ class ITSTrackTask : public TaskInterface o2::itsmft::TopologyDictionary* mDict; private: - //analysis for its-only residual + // analysis for its-only residual o2::its::GeometryTGeo* mGeom; - std::vector FitStepSize{0.3, 1.0e-5, 1.0e-5, 1.0e-5}; + std::vector FitStepSize{ 0.3, 1.0e-5, 1.0e-5, 1.0e-5 }; double FitTolerance = 1.0e-8; - double ITS_AbsBz = 0.5; //T + double ITS_AbsBz = 0.5; // T - double inputG_C[3*NLayer]; + double inputG_C[3 * NLayer]; double FitparXY[6]; double FitparDZ[2]; ROOT::Fit::Fitter fitterA; @@ -145,124 +145,131 @@ class ITSTrackTask : public TaskInterface Int_t mResMonNclMin = 0; float mResMonTrackMinPt = 0; - std::array, NLayer> hResidualXY{};//[NLayer]; - std::array, NLayer> hResidualZD{};//[NLayer]; + std::array, NLayer> hResidualXY{}; //[NLayer]; + std::array, NLayer> hResidualZD{}; //[NLayer]; - void circleFitXY(double* input, double* par, double &MSEvalue, std::vector hitUpdate, int step = 0); + void circleFitXY(double* input, double* par, double& MSEvalue, std::vector hitUpdate, int step = 0); - //default setting - // function Object to be minimized + // default setting + // function Object to be minimized struct se_circlefitXY { // the TGraph is a data member of the object std::vector fHits; double thetaR; double Bz; - double sigma_meas[2][NLayer] = {{45,45,45,55,55,55,55}, - {40,40,40,40,40,40,40}}; //um unit - double sigma_msc[2][NLayer] = {{30,30,30,110,110,110,110}, - {25,25,25,75,75,75,75}}; //um unit + double sigma_meas[2][NLayer] = { { 45, 45, 45, 55, 55, 55, 55 }, + { 40, 40, 40, 40, 40, 40, 40 } }; // um unit + double sigma_msc[2][NLayer] = { { 30, 30, 30, 110, 110, 110, 110 }, + { 25, 25, 25, 75, 75, 75, 75 } }; // um unit se_circlefitXY() = default; - se_circlefitXY(std::vector &h, double tR, double bz) : fHits(h), thetaR(tR), Bz(bz) { }; + se_circlefitXY(std::vector& h, double tR, double bz) : fHits(h), thetaR(tR), Bz(bz) {}; - void loadparameters(double arrpar_meas[][NLayer], double arrpar_msc[][NLayer]){ - for(int a = 0; a < 2; a++){ - for(int l = 0; l < NLayer; l++){ + void loadparameters(double arrpar_meas[][NLayer], double arrpar_msc[][NLayer]) + { + for (int a = 0; a < 2; a++) { + for (int l = 0; l < NLayer; l++) { sigma_meas[a][l] = arrpar_meas[a][l]; sigma_msc[a][l] = arrpar_msc[a][l]; } } }; - void init(){ + void init() + { fHits.clear(); thetaR = 0; Bz = 0; }; - void set(std::vector &h, double tR, double bz){ + void set(std::vector& h, double tR, double bz) + { fHits = h; thetaR = tR; Bz = bz; }; - double getsigma(double R, int L, double B, int axis){ - //R : cm - //B : T - if(L<0 || L>=NLayer) return 1; - double aL = sigma_meas[axis][L]*1e-4; //um -> cm - double bL = sigma_msc[axis][L]*1e-4; //um -> cm - double Beff = 0.299792458*B; - double sigma = std::sqrt( std::pow(aL,2) + std::pow(bL,2)/(std::pow(Beff,2)*std::pow(R*1e-2,2))); + double getsigma(double R, int L, double B, int axis) + { + // R : cm + // B : T + if (L < 0 || L >= NLayer) + return 1; + double aL = sigma_meas[axis][L] * 1e-4; // um -> cm + double bL = sigma_msc[axis][L] * 1e-4; // um -> cm + double Beff = 0.299792458 * B; + double sigma = std::sqrt(std::pow(aL, 2) + std::pow(bL, 2) / (std::pow(Beff, 2) * std::pow(R * 1e-2, 2))); return sigma; }; // calculate distance line-point - double distance2(double x, double y, const double *p, double tR, int charge) { + double distance2(double x, double y, const double* p, double tR, int charge) + { + + double R = std::abs(1 / p[0]); + + double Xc = p[0] > 0 ? R * std::cos(p[1] + tR + 0.5 * TMath::Pi()) : R * std::cos(p[1] + tR - 0.5 * TMath::Pi()); + double Yc = p[0] > 0 ? R * std::sin(p[1] + tR + 0.5 * TMath::Pi()) : R * std::sin(p[1] + tR - 0.5 * TMath::Pi()); - double R = std::abs(1/p[0]); - - double Xc = p[0] > 0 ? R*std::cos(p[1] + tR + 0.5*TMath::Pi()) : R*std::cos(p[1] + tR - 0.5*TMath::Pi()); - double Yc = p[0] > 0 ? R*std::sin(p[1] + tR + 0.5*TMath::Pi()) : R*std::sin(p[1] + tR - 0.5*TMath::Pi()); - double dx = x - (Xc + p[2]); - double dy = y - (Yc + p[3]); - double dxy = R - std::sqrt(dx*dx + dy*dy); - - double d2 = dxy*dxy; + double dy = y - (Yc + p[3]); + double dxy = R - std::sqrt(dx * dx + dy * dy); + + double d2 = dxy * dxy; return d2; }; // implementation of the function to be minimized - double operator() (const double *par) { //const double -> double - assert(fHits != 0); + double operator()(const double* par) + { // const double -> double + assert(fHits != 0); - int nhits = fHits.size(); - double sum = 0; + int nhits = fHits.size(); + double sum = 0; - double charge = par[0]>0 ? +1 : -1; - double RecRadius = std::abs(1/par[0]); + double charge = par[0] > 0 ? +1 : -1; + double RecRadius = std::abs(1 / par[0]); double Sigma_tot[NLayer]; double sum_Sigma_tot = 0; - for(int l = 0; l < nhits; l++){ + for (int l = 0; l < nhits; l++) { Sigma_tot[l] = getsigma(RecRadius, fHits[l].Z(), Bz, 1); - sum_Sigma_tot += 1/(std::pow(Sigma_tot[l],2)); - } - + sum_Sigma_tot += 1 / (std::pow(Sigma_tot[l], 2)); + } + for (int l = 0; l < nhits; l++) { - double d = distance2(fHits[l].X(), fHits[l].Y(), par, thetaR, charge)/(std::pow(Sigma_tot[l],2)); - sum += d; + double d = distance2(fHits[l].X(), fHits[l].Y(), par, thetaR, charge) / (std::pow(Sigma_tot[l], 2)); + sum += d; } sum = sum / sum_Sigma_tot; return sum; }; - }; se_circlefitXY fitfuncXY; void lineFitDZ(double* zIn, double* betaIn, double* parz, double Radius, bool vertex, std::vector hitUpdate); - double mSigmaMeas[2][NLayer] = {{45,45,45,55,55,55,55}, - {40,40,40,40,40,40,40}}; //um unit - double mSigmaMsc[2][NLayer] = {{30,30,30,110,110,110,110}, - {25,25,25,75,75,75,75}}; //um unit + double mSigmaMeas[2][NLayer] = { { 45, 45, 45, 55, 55, 55, 55 }, + { 40, 40, 40, 40, 40, 40, 40 } }; // um unit + double mSigmaMsc[2][NLayer] = { { 30, 30, 30, 110, 110, 110, 110 }, + { 25, 25, 25, 75, 75, 75, 75 } }; // um unit - double getsigma(double R, int L, double B, int axis){ - //R : cm - //B : T - if(L<0 || L>=NLayer) return 1; - double aL = mSigmaMeas[axis][L]*1e-4; //um -> cm - double bL = mSigmaMsc[axis][L]*1e-4; //um -> cm - double Beff = 0.299792458*B; - double sigma = std::sqrt( std::pow(aL,2) + std::pow(bL,2)/(std::pow(Beff,2)*std::pow(R*1e-2,2))); + double getsigma(double R, int L, double B, int axis) + { + // R : cm + // B : T + if (L < 0 || L >= NLayer) + return 1; + double aL = mSigmaMeas[axis][L] * 1e-4; // um -> cm + double bL = mSigmaMsc[axis][L] * 1e-4; // um -> cm + double Beff = 0.299792458 * B; + double sigma = std::sqrt(std::pow(aL, 2) + std::pow(bL, 2) / (std::pow(Beff, 2) * std::pow(R * 1e-2, 2))); return sigma; }; - }; } // namespace o2::quality_control_modules::its diff --git a/Modules/ITS/src/ITSTrackTask.cxx b/Modules/ITS/src/ITSTrackTask.cxx index 91be439bcc..af949e9a1f 100644 --- a/Modules/ITS/src/ITSTrackTask.cxx +++ b/Modules/ITS/src/ITSTrackTask.cxx @@ -67,16 +67,16 @@ void ITSTrackTask::initialize(o2::framework::InitContext& /*ctx*/) mInvMasses = o2::quality_control_modules::common::getFromConfig(mCustomParameters, "InvMasses", mInvMasses); mPublishMore = o2::quality_control_modules::common::getFromConfig(mCustomParameters, "publishMore", mPublishMore); - //analysis for its-only residual + // analysis for its-only residual mAlignmentMonitor = o2::quality_control_modules::common::getFromConfig(mCustomParameters, "doAlignmentMonitor", mAlignmentMonitor); mDefaultMomResPar = o2::quality_control_modules::common::getFromConfig(mCustomParameters, "UseDefaultMomResPar", mDefaultMomResPar); - if(mAlignmentMonitor == 1 && mDefaultMomResPar == 0){ + if (mAlignmentMonitor == 1 && mDefaultMomResPar == 0) { std::vector vMomResParMEAS1 = convertToArray(o2::quality_control_modules::common::getFromConfig(mCustomParameters, "MomResParMEAS1", "")); std::vector vMomResParMEAS2 = convertToArray(o2::quality_control_modules::common::getFromConfig(mCustomParameters, "MomResParMEAS2", "")); std::vector vMomResParMSC1 = convertToArray(o2::quality_control_modules::common::getFromConfig(mCustomParameters, "MomResParMSC1", "")); std::vector vMomResParMSC2 = convertToArray(o2::quality_control_modules::common::getFromConfig(mCustomParameters, "MomResParMSC2", "")); - for(int l = 0; l < NLayer; l++){ + for (int l = 0; l < NLayer; l++) { mSigmaMeas[0][l] = (double)vMomResParMEAS1[l]; mSigmaMeas[1][l] = (double)vMomResParMEAS2[l]; mSigmaMsc[0][l] = (double)vMomResParMSC1[l]; @@ -122,7 +122,7 @@ void ITSTrackTask::monitorData(o2::framework::ProcessingContext& ctx) mDict = TaskInterface::retrieveConditionAny("ITS/Calib/ClusterDictionary", metadata, ts); ILOG(Debug, Devel) << "Dictionary size: " << mDict->getSize() << ENDM; - if(mAlignmentMonitor == 1) { + if (mAlignmentMonitor == 1) { o2::its::GeometryTGeo::adopt(TaskInterface::retrieveConditionAny("ITS/Config/Geometry", metadata, ts)); mGeom = o2::its::GeometryTGeo::Instance(); if (!mGeom) { @@ -259,9 +259,10 @@ void ITSTrackTask::monitorData(o2::framework::ProcessingContext& ctx) nClusterCntTrack += track.getNumberOfClusters(); std::vector hitUpdate; - int chipIDs[7] = {-1, -1, -1, -1, -1, -1, -1}; - if(mAlignmentMonitor == 1 && out.getPt()>mResMonTrackMinPt && track.getNumberOfClusters() >= mResMonNclMin) { - for (int iLayer = 0; iLayer < NLayer; iLayer++) hitUpdate.push_back(false); + int chipIDs[7] = { -1, -1, -1, -1, -1, -1, -1 }; + if (mAlignmentMonitor == 1 && out.getPt() > mResMonTrackMinPt && track.getNumberOfClusters() >= mResMonNclMin) { + for (int iLayer = 0; iLayer < NLayer; iLayer++) + hitUpdate.push_back(false); } for (int icluster = 0; icluster < track.getNumberOfClusters(); icluster++) { @@ -271,7 +272,7 @@ void ITSTrackTask::monitorData(o2::framework::ProcessingContext& ctx) int ClusterID = cluster.getPatternID(); // used for normal (frequent) cluster shapes int layer = 0; - while (ChipID >= ChipBoundary[layer] && layer <= NLayer){ + while (ChipID >= ChipBoundary[layer] && layer <= NLayer) { layer++; } layer--; @@ -281,12 +282,13 @@ void ITSTrackTask::monitorData(o2::framework::ProcessingContext& ctx) hNClusterVsChipITS->Fill(ChipID + 1, clusterSizeWithCorrection); } - //Residual Monitoring - if(mAlignmentMonitor == 1 && out.getPt()>mResMonTrackMinPt && track.getNumberOfClusters() >= mResMonNclMin) { + // Residual Monitoring + if (mAlignmentMonitor == 1 && out.getPt() > mResMonTrackMinPt && track.getNumberOfClusters() >= mResMonNclMin) { - if(layer < 0 || layer >= NLayer) continue; + if (layer < 0 || layer >= NLayer) + continue; - o2::math_utils::Point3D locC; // local coordinates + o2::math_utils::Point3D locC; // local coordinates if (ClusterID != o2::itsmft::CompCluster::InvalidPatternID) { // Normal (frequent) cluster shapes if (!mDict->isGroup(ClusterID)) { locC = mDict->getClusterCoordinates(cluster); @@ -295,67 +297,74 @@ void ITSTrackTask::monitorData(o2::framework::ProcessingContext& ctx) locC = mDict->getClusterCoordinates(cluster, patt, true); } } else { // invalid pattern - continue; + continue; } hitUpdate[layer] = true; - auto gloC = mGeom->getMatrixL2G(ChipID) * locC; - inputG_C[3*layer + 0] = gloC.X(); - inputG_C[3*layer + 1] = gloC.Y(); - inputG_C[3*layer + 2] = gloC.Z(); + auto gloC = mGeom->getMatrixL2G(ChipID) * locC; + inputG_C[3 * layer + 0] = gloC.X(); + inputG_C[3 * layer + 1] = gloC.Y(); + inputG_C[3 * layer + 2] = gloC.Z(); chipIDs[layer] = ChipID; } } - //Residual Monitoring - if(mAlignmentMonitor == 1 && out.getPt()>mResMonTrackMinPt && track.getNumberOfClusters() >= mResMonNclMin) { + // Residual Monitoring + if (mAlignmentMonitor == 1 && out.getPt() > mResMonTrackMinPt && track.getNumberOfClusters() >= mResMonNclMin) { int NclValid = 0; for (int iLayer = 0; iLayer < NLayer; iLayer++) { - if(hitUpdate[iLayer]) NclValid++; + if (hitUpdate[iLayer]) + NclValid++; } - if(NclValid0 ? RecRadius*std::cos(FitparXY[1]+FitparXY[4] + 0.5*TMath::Pi()) : RecRadius*std::cos(FitparXY[1]+FitparXY[4] - 0.5*TMath::Pi()); - double CircleYc = FitparXY[0] >0 ? RecRadius*std::sin(FitparXY[1]+FitparXY[4] + 0.5*TMath::Pi()) : RecRadius*std::sin(FitparXY[1]+FitparXY[4] - 0.5*TMath::Pi()); - + double RecRadius = std::abs(1 / FitparXY[0]); + double CircleXc = FitparXY[0] > 0 ? RecRadius * std::cos(FitparXY[1] + FitparXY[4] + 0.5 * TMath::Pi()) : RecRadius * std::cos(FitparXY[1] + FitparXY[4] - 0.5 * TMath::Pi()); + double CircleYc = FitparXY[0] > 0 ? RecRadius * std::sin(FitparXY[1] + FitparXY[4] + 0.5 * TMath::Pi()) : RecRadius * std::sin(FitparXY[1] + FitparXY[4] - 0.5 * TMath::Pi()); + CircleXc += FitparXY[2]; CircleYc += FitparXY[3]; - for (int ipar = 0; ipar < 2; ipar++) FitparDZ[ipar] = 0; + for (int ipar = 0; ipar < 2; ipar++) + FitparDZ[ipar] = 0; double z_meas[NLayer]; double beta[NLayer]; for (int iLayer = 0; iLayer < NLayer; iLayer++) { - z_meas[iLayer] = inputG_C[3*iLayer + 2]; - beta[iLayer] = std::atan2(inputG_C[3*iLayer + 1] - CircleYc, inputG_C[3*iLayer + 0] - CircleXc); + z_meas[iLayer] = inputG_C[3 * iLayer + 2]; + beta[iLayer] = std::atan2(inputG_C[3 * iLayer + 1] - CircleYc, inputG_C[3 * iLayer + 0] - CircleXc); } lineFitDZ(z_meas, beta, FitparDZ, RecRadius, false, hitUpdate); - + for (int iLayer = 0; iLayer < NLayer; iLayer++) { - if(chipIDs[iLayer]<0) continue; - double meas_GXc = inputG_C[(3*iLayer)+0]; //alpha - double meas_GYc = inputG_C[(3*iLayer)+1]; //beta - double meas_GZc = inputG_C[(3*iLayer)+2]; //gamma - double proj_GXc = RecRadius*std::cos(beta[iLayer]) + CircleXc; - double proj_GYc = RecRadius*std::sin(beta[iLayer]) + CircleYc; - double proj_GZc = (FitparDZ[0])*(beta[iLayer]) + (FitparDZ[1]); + if (chipIDs[iLayer] < 0) + continue; + double meas_GXc = inputG_C[(3 * iLayer) + 0]; // alpha + double meas_GYc = inputG_C[(3 * iLayer) + 1]; // beta + double meas_GZc = inputG_C[(3 * iLayer) + 2]; // gamma + double proj_GXc = RecRadius * std::cos(beta[iLayer]) + CircleXc; + double proj_GYc = RecRadius * std::sin(beta[iLayer]) + CircleYc; + double proj_GZc = (FitparDZ[0]) * (beta[iLayer]) + (FitparDZ[1]); TVector3 measXY(meas_GXc, meas_GYc, 0); TVector3 projXY(proj_GXc, proj_GYc, 0); double sign = +1; - if(measXY.Cross(projXY).Z()>0) sign = +1; - else sign = -1; - double dxy = sign*std::sqrt(std::pow(proj_GXc-meas_GXc,2) + std::pow(proj_GYc-meas_GYc,2)); - double dz = proj_GZc - meas_GZc; - - hResidualXY[iLayer]->Fill(dxy,chipIDs[iLayer]); - hResidualZD[iLayer]->Fill(dz,chipIDs[iLayer]); + if (measXY.Cross(projXY).Z() > 0) + sign = +1; + else + sign = -1; + double dxy = sign * std::sqrt(std::pow(proj_GXc - meas_GXc, 2) + std::pow(proj_GYc - meas_GYc, 2)); + double dz = proj_GZc - meas_GZc; + + hResidualXY[iLayer]->Fill(dxy, chipIDs[iLayer]); + hResidualZD[iLayer]->Fill(dz, chipIDs[iLayer]); } } @@ -559,7 +568,7 @@ void ITSTrackTask::reset() hTrackPtVsEta->Reset(); hTrackPtVsPhi->Reset(); - if(mAlignmentMonitor == 1) { + if (mAlignmentMonitor == 1) { for (int l = 0; l < NLayer; l++) { hResidualXY[l]->Reset(); hResidualZD[l]->Reset(); @@ -753,19 +762,19 @@ void ITSTrackTask::createAllHistos() formatAxes(hTrackPtVsPhi, "#it{p}_{T} (GeV/#it{c})", "#phi", 1, 1.10); hTrackPtVsPhi->SetStats(0); - if(mAlignmentMonitor == 1) { + if (mAlignmentMonitor == 1) { for (int l = 0; l < NLayer; l++) { - //sensor - int NBinsChipID = (l < NLayerIB) ? (ChipBoundary[l+1] - ChipBoundary[l]) : (ChipBoundary[l+1] - ChipBoundary[l])/14; + // sensor + int NBinsChipID = (l < NLayerIB) ? (ChipBoundary[l + 1] - ChipBoundary[l]) : (ChipBoundary[l + 1] - ChipBoundary[l]) / 14; - hResidualXY[l] = std::make_unique(Form("hResidualXY%d",l), Form("ChipID vs dxy, Layer %d",l), - 500, -0.05, 0.05, NBinsChipID, ChipBoundary[l] - 0.5, ChipBoundary[l+1] - 0.5); + hResidualXY[l] = std::make_unique(Form("hResidualXY%d", l), Form("ChipID vs dxy, Layer %d", l), + 500, -0.05, 0.05, NBinsChipID, ChipBoundary[l] - 0.5, ChipBoundary[l + 1] - 0.5); addObject(hResidualXY[l].get()); formatAxes(hResidualXY[l].get(), "dxy(cm)", (l < NLayerIB) ? "ChipID(Sensor Unit)" : "ChipID(HIC Unit)", 1, 1.10); hResidualXY[l]->SetStats(0); - hResidualZD[l] = std::make_unique(Form("hResidualZD%d",l), Form("ChipID vs dz, Layer %d",l), - 500, -0.05, 0.05, NBinsChipID, ChipBoundary[l] - 0.5, ChipBoundary[l+1] - 0.5); + hResidualZD[l] = std::make_unique(Form("hResidualZD%d", l), Form("ChipID vs dz, Layer %d", l), + 500, -0.05, 0.05, NBinsChipID, ChipBoundary[l] - 0.5, ChipBoundary[l + 1] - 0.5); addObject(hResidualZD[l].get()); formatAxes(hResidualZD[l].get(), "dz(cm)", (l < NLayerIB) ? "ChipID(Sensor Unit)" : "ChipID(HIC Unit)", 1, 1.10); hResidualZD[l]->SetStats(0); @@ -791,125 +800,135 @@ void ITSTrackTask::publishHistos() } } - -void ITSTrackTask::circleFitXY(double* input, double* par, double &MSEvalue, std::vector hitUpdate, int step){ +void ITSTrackTask::circleFitXY(double* input, double* par, double& MSEvalue, std::vector hitUpdate, int step) +{ int hitentries = 0; - for(int a=0; a hr; - + std::vector hr; + double frphiX = 0; double frphiY = 0; - int nfrdet = 0; - for(int a=0; a i, j, k; //[hitentries] - std::vector irot, jrot; //[hitentries] - - int fa=0; - int index[7] = {-1, -1, -1, -1, -1, -1, -1}; - for(int a=0; a i, j, k; //[hitentries] + std::vector irot, jrot; //[hitentries] + + int fa = 0; + int index[7] = { -1, -1, -1, -1, -1, -1, -1 }; + for (int a = 0; a < hitUpdate.size(); a++) { + if (hitUpdate[a] == false) + continue; + i.push_back(input[(3 * a) + 0]); + j.push_back(input[(3 * a) + 1]); + k.push_back(input[(3 * a) + 2]); + TMatrixD gloX[2]; - gloX[0].ResizeTo(1,2); - gloX[0][0] = { i[fa], j[fa]}; + gloX[0].ResizeTo(1, 2); + gloX[0][0] = { i[fa], j[fa] }; gloX[0].T(); - gloX[1].ResizeTo(2,1); + gloX[1].ResizeTo(2, 1); gloX[1] = RotF * gloX[0]; - gloX[1].T(); - + gloX[1].T(); + irot.push_back(gloX[1][0][0]); jrot.push_back(gloX[1][0][1]); - hr.push_back(TVector3(irot[fa], jrot[fa], a)); - index[a]=fa; - fa++; + hr.push_back(TVector3(irot[fa], jrot[fa], a)); + index[a] = fa; + fa++; } - int cntR[] = {0, 0}; + int cntR[] = { 0, 0 }; std::vector initR; - //standard seeding - //012 + (2)34(5) + 56 - int hit1[] = {0, 1, 2}; - int hit2[] = {3, 4, 5, 2}; bool hit_mid = false; - int hit3[] = {6, 5}; - for(int i1 = 0; i1 < 3; i1++){ + // standard seeding + // 012 + (2)34(5) + 56 + int hit1[] = { 0, 1, 2 }; + int hit2[] = { 3, 4, 5, 2 }; + bool hit_mid = false; + int hit3[] = { 6, 5 }; + for (int i1 = 0; i1 < 3; i1++) { // i1 -> hit1[i1] - for(int i2 = 0; i2 < 4; i2++){ + for (int i2 = 0; i2 < 4; i2++) { // i2 -> hit2[i2] - if(hit_mid==true && i2>= 2) continue; - for(int i3 = 0; i3 < 2; i3++){ + if (hit_mid == true && i2 >= 2) + continue; + for (int i3 = 0; i3 < 2; i3++) { // i3 -> hit3[i3] - if(hit1[i1]==hit2[i2] || hit2[i2]==hit3[i3]) continue; - //if(hit_mid==true) continue; - if(hitUpdate[hit1[i1]]==false) continue; - if(hitUpdate[hit2[i2]]==false) continue; - if(hitUpdate[hit3[i3]]==false) continue; - double hitX[3] = {i[index[hit1[i1]]], i[index[hit2[i2]]], i[index[hit3[i3]]]}; - double hitY[3] = {j[index[hit1[i1]]], j[index[hit2[i2]]], j[index[hit3[i3]]]}; - - double d12 = -(hitX[1] - hitX[0])/(hitY[1] - hitY[0]); - double d23 = -(hitX[2] - hitX[1])/(hitY[2] - hitY[1]); - - double x12 = 0.5*(hitX[0] + hitX[1]); - double x23 = 0.5*(hitX[1] + hitX[2]); - double y12 = 0.5*(hitY[0] + hitY[1]); - double y23 = 0.5*(hitY[1] + hitY[2]); - - double CenterX = ((-d23*x23 + d12*x12) + (y23 - y12))/(-d23 + d12); - double CenterY = d12*(CenterX - x12) + y12; - - double temp_R = std::sqrt(std::pow(CenterX - hitX[0],2) + std::pow(CenterY - hitY[0],2)); - initR.push_back(TVector3(CenterX,CenterY,temp_R)); - if(i2<2) { - hit_mid=true; // mid hit is successfully used. Do not find inner or outer hits for initial radius searching + if (hit1[i1] == hit2[i2] || hit2[i2] == hit3[i3]) + continue; + // if(hit_mid==true) continue; + if (hitUpdate[hit1[i1]] == false) + continue; + if (hitUpdate[hit2[i2]] == false) + continue; + if (hitUpdate[hit3[i3]] == false) + continue; + double hitX[3] = { i[index[hit1[i1]]], i[index[hit2[i2]]], i[index[hit3[i3]]] }; + double hitY[3] = { j[index[hit1[i1]]], j[index[hit2[i2]]], j[index[hit3[i3]]] }; + + double d12 = -(hitX[1] - hitX[0]) / (hitY[1] - hitY[0]); + double d23 = -(hitX[2] - hitX[1]) / (hitY[2] - hitY[1]); + + double x12 = 0.5 * (hitX[0] + hitX[1]); + double x23 = 0.5 * (hitX[1] + hitX[2]); + double y12 = 0.5 * (hitY[0] + hitY[1]); + double y23 = 0.5 * (hitY[1] + hitY[2]); + + double CenterX = ((-d23 * x23 + d12 * x12) + (y23 - y12)) / (-d23 + d12); + double CenterY = d12 * (CenterX - x12) + y12; + + double temp_R = std::sqrt(std::pow(CenterX - hitX[0], 2) + std::pow(CenterY - hitY[0], 2)); + initR.push_back(TVector3(CenterX, CenterY, temp_R)); + if (i2 < 2) { + hit_mid = true; // mid hit is successfully used. Do not find inner or outer hits for initial radius searching } cntR[0]++; - } + } } } - if(initR.size()==0) { - initR.push_back(TVector3(0,0,10000)); + if (initR.size() == 0) { + initR.push_back(TVector3(0, 0, 10000)); cntR[0]++; } - double mean_X[] = {0, 0}; - double mean_Y[] = {0, 0}; - double mean_R[] = {0, 0}; - for(int i = 0; i < cntR[0]; i++) { - mean_X[0] += initR[i].X()/(double)cntR[0]; - mean_Y[0] += initR[i].Y()/(double)cntR[0]; - mean_R[0] += initR[i](2)/(double)cntR[0]; + double mean_X[] = { 0, 0 }; + double mean_Y[] = { 0, 0 }; + double mean_R[] = { 0, 0 }; + for (int i = 0; i < cntR[0]; i++) { + mean_X[0] += initR[i].X() / (double)cntR[0]; + mean_Y[0] += initR[i].Y() / (double)cntR[0]; + mean_R[0] += initR[i](2) / (double)cntR[0]; } - for(int i = 0; i < cntR[0]; i++) { - if(std::abs(mean_R[0] - initR[i](2)) < mean_R[0]) { + for (int i = 0; i < cntR[0]; i++) { + if (std::abs(mean_R[0] - initR[i](2)) < mean_R[0]) { mean_X[1] += initR[i].X(); mean_Y[1] += initR[i].Y(); mean_R[1] += initR[i](2); @@ -918,270 +937,274 @@ void ITSTrackTask::circleFitXY(double* input, double* par, double &MSEvalue, std } mean_R[1] /= cntR[1]; - mean_R[1] *= std::pow(sqrt(10),step); - if(mean_R[1]<1.0e+1) mean_R[1] = 1.0e+1; - if(mean_R[1]>1.0e+6) mean_R[1] = 1.0e+6; + mean_R[1] *= std::pow(sqrt(10), step); + if (mean_R[1] < 1.0e+1) + mean_R[1] = 1.0e+1; + if (mean_R[1] > 1.0e+6) + mean_R[1] = 1.0e+6; + + double thetaR = std::atan2(jrot[0], irot[0]); - double thetaR = std::atan2( jrot[0], irot[0]); - double temp_parA[4]; - temp_parA[0] = + 1/mean_R[1]; - temp_parA[1] = 0; + temp_parA[0] = +1 / mean_R[1]; + temp_parA[1] = 0; double temp_parB[4]; - temp_parB[0] = - 1/mean_R[1]; - temp_parB[1] = 0; + temp_parB[0] = -1 / mean_R[1]; + temp_parB[1] = 0; // make the functor object fitfuncXY.init(); fitfuncXY.set(hr, thetaR, ITS_AbsBz); ROOT::Math::Functor fcn(fitfuncXY, 4); - double pStartA[4] = {temp_parA[0],temp_parA[1], 0, 0}; - fitterA.SetFCN(fcn,pStartA); + double pStartA[4] = { temp_parA[0], temp_parA[1], 0, 0 }; + fitterA.SetFCN(fcn, pStartA); fitterA.Config().ParSettings(0).SetStepSize((float)FitStepSize[0]); fitterA.Config().ParSettings(1).SetStepSize((float)FitStepSize[1]); fitterA.Config().ParSettings(2).SetStepSize((float)FitStepSize[2]); - fitterA.Config().ParSettings(3).SetStepSize((float)FitStepSize[3]); + fitterA.Config().ParSettings(3).SetStepSize((float)FitStepSize[3]); fitterA.Config().ParSettings(0).SetLimits(+1.0e-10, +1.0e-1); // + side - double pStartB[4] = {temp_parB[0],temp_parB[1], 0, 0}; - fitterB.SetFCN(fcn,pStartB); + double pStartB[4] = { temp_parB[0], temp_parB[1], 0, 0 }; + fitterB.SetFCN(fcn, pStartB); fitterB.Config().ParSettings(0).SetStepSize((float)FitStepSize[0]); fitterB.Config().ParSettings(1).SetStepSize((float)FitStepSize[1]); fitterB.Config().ParSettings(2).SetStepSize((float)FitStepSize[2]); - fitterB.Config().ParSettings(3).SetStepSize((float)FitStepSize[3]); + fitterB.Config().ParSettings(3).SetStepSize((float)FitStepSize[3]); fitterB.Config().ParSettings(0).SetLimits(-1.0e-1, -1.0e-10); // - side ROOT::Math::MinimizerOptions minOpt; - for(int iTol = 0; iTol < 4; iTol++){ + for (int iTol = 0; iTol < 4; iTol++) { - minOpt.SetTolerance(std::pow(10,iTol)*FitTolerance); + minOpt.SetTolerance(std::pow(10, iTol) * FitTolerance); fitterA.Config().SetMinimizerOptions(minOpt); - fitterB.Config().SetMinimizerOptions(minOpt); + fitterB.Config().SetMinimizerOptions(minOpt); bool okA = fitterA.FitFCN(); bool okB = fitterB.FitFCN(); - + if (!okA) { - if(!okB) { - const ROOT::Fit::FitResult & resultA = fitterA.Result(); - const double * parFitA = resultA.GetParams(); - double MSEvalueA = resultA.MinFcnValue(); - const ROOT::Fit::FitResult & resultB = fitterB.Result(); - const double * parFitB = resultB.GetParams(); - double MSEvalueB = resultB.MinFcnValue(); - + if (!okB) { + const ROOT::Fit::FitResult& resultA = fitterA.Result(); + const double* parFitA = resultA.GetParams(); + double MSEvalueA = resultA.MinFcnValue(); + const ROOT::Fit::FitResult& resultB = fitterB.Result(); + const double* parFitB = resultB.GetParams(); + double MSEvalueB = resultB.MinFcnValue(); + TMatrixD vxyA[2]; - vxyA[0].ResizeTo(1,2); - vxyA[0][0] = { parFitA[2], parFitA[3]}; + vxyA[0].ResizeTo(1, 2); + vxyA[0][0] = { parFitA[2], parFitA[3] }; vxyA[0].T(); - vxyA[1].ResizeTo(2,1); + vxyA[1].ResizeTo(2, 1); vxyA[1] = RotFInv * vxyA[0]; - vxyA[1].T(); - + vxyA[1].T(); + TMatrixD vxyB[2]; - vxyB[0].ResizeTo(1,2); - vxyB[0][0] = { parFitB[2], parFitB[3]}; + vxyB[0].ResizeTo(1, 2); + vxyB[0][0] = { parFitB[2], parFitB[3] }; vxyB[0].T(); - vxyB[1].ResizeTo(2,1); - vxyB[1] = RotFInv * vxyB[0]; + vxyB[1].ResizeTo(2, 1); + vxyB[1] = RotFInv * vxyB[0]; vxyB[1].T(); - if(MSEvalueA hitUpdate){ +void ITSTrackTask::lineFitDZ(double* zIn, double* betaIn, double* parz, double Radius, bool vertex, std::vector hitUpdate) +{ int hitentries = 0; - for(int a=0; a z, beta; + if (vertex == true) + hitentries++; + + std::vector z, beta; std::vector index; - int fa=0; - for(int a=0; a Date: Mon, 28 Oct 2024 12:23:08 +0900 Subject: [PATCH 6/6] ITS alignment monitoring (check with QCv20240903), revised 20241028(1), clang-format applied, TH2I->TH2F --- Modules/ITS/include/ITS/ITSTrackTask.h | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/Modules/ITS/include/ITS/ITSTrackTask.h b/Modules/ITS/include/ITS/ITSTrackTask.h index 9841c5194b..01fe7fbb12 100644 --- a/Modules/ITS/include/ITS/ITSTrackTask.h +++ b/Modules/ITS/include/ITS/ITSTrackTask.h @@ -157,13 +157,16 @@ class ITSTrackTask : public TaskInterface std::vector fHits; double thetaR; double Bz; - double sigma_meas[2][NLayer] = { { 45, 45, 45, 55, 55, 55, 55 }, - { 40, 40, 40, 40, 40, 40, 40 } }; // um unit - double sigma_msc[2][NLayer] = { { 30, 30, 30, 110, 110, 110, 110 }, - { 25, 25, 25, 75, 75, 75, 75 } }; // um unit + double sigma_meas[2][NLayer] = { { 45, 45, 45, 55, 55, 55, 55 }, { 40, 40, 40, 40, 40, 40, 40 } }; // um unit + double sigma_msc[2][NLayer] = { { 30, 30, 30, 110, 110, 110, 110 }, { 25, 25, 25, 75, 75, 75, 75 } }; // um unit se_circlefitXY() = default; - se_circlefitXY(std::vector& h, double tR, double bz) : fHits(h), thetaR(tR), Bz(bz) {}; + se_circlefitXY(std::vector& h, double tR, double bz) + { + fHits = h; + thetaR = tR; + Bz = bz; + }; void loadparameters(double arrpar_meas[][NLayer], double arrpar_msc[][NLayer]) { @@ -252,10 +255,8 @@ class ITSTrackTask : public TaskInterface void lineFitDZ(double* zIn, double* betaIn, double* parz, double Radius, bool vertex, std::vector hitUpdate); - double mSigmaMeas[2][NLayer] = { { 45, 45, 45, 55, 55, 55, 55 }, - { 40, 40, 40, 40, 40, 40, 40 } }; // um unit - double mSigmaMsc[2][NLayer] = { { 30, 30, 30, 110, 110, 110, 110 }, - { 25, 25, 25, 75, 75, 75, 75 } }; // um unit + double mSigmaMeas[2][NLayer] = { { 45, 45, 45, 55, 55, 55, 55 }, { 40, 40, 40, 40, 40, 40, 40 } }; // um unit + double mSigmaMsc[2][NLayer] = { { 30, 30, 30, 110, 110, 110, 110 }, { 25, 25, 25, 75, 75, 75, 75 } }; // um unit double getsigma(double R, int L, double B, int axis) {