From a56d41e254b3752c0dc20254b4e8a6cd40cf5669 Mon Sep 17 00:00:00 2001 From: delico18 Date: Mon, 21 Oct 2024 14:38:59 +0900 Subject: [PATCH] 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;