From dfc419c1e0075caa8d89c64664dd8eb8052d52cf Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Tue, 23 Jan 2024 15:33:02 +0000 Subject: [PATCH] Please consider the following formatting changes --- .../Aligner/include/MCHAlign/MinResSolve.h | 1 - .../MCH/Align/Aligner/src/AlignmentSpec.cxx | 1661 ++++++++--------- .../Align/Aligner/src/alignment-workflow.cxx | 6 +- .../MCH/Tracking/include/MCHTracking/Track.h | 1 - 4 files changed, 813 insertions(+), 856 deletions(-) diff --git a/Detectors/MUON/MCH/Align/Aligner/include/MCHAlign/MinResSolve.h b/Detectors/MUON/MCH/Align/Aligner/include/MCHAlign/MinResSolve.h index 8b6ab97833566..15ddaf5ed28a6 100644 --- a/Detectors/MUON/MCH/Align/Aligner/include/MCHAlign/MinResSolve.h +++ b/Detectors/MUON/MCH/Align/Aligner/include/MCHAlign/MinResSolve.h @@ -13,7 +13,6 @@ /// \brief General class (from AliROOT) for solving large system of linear equations /// \author ruben.shahoyan@cern.ch - #ifndef ALICEO2_MCH_MINRESSOLVE_H #define ALICEO2_MCH_MINRESSOLVE_H diff --git a/Detectors/MUON/MCH/Align/Aligner/src/AlignmentSpec.cxx b/Detectors/MUON/MCH/Align/Aligner/src/AlignmentSpec.cxx index 308f3575d177a..37fcedea7df60 100644 --- a/Detectors/MUON/MCH/Align/Aligner/src/AlignmentSpec.cxx +++ b/Detectors/MUON/MCH/Align/Aligner/src/AlignmentSpec.cxx @@ -13,7 +13,6 @@ /// \brief Implementation of alignment process for muon spectrometer /// \author Chi ZHANG(CEA-Saclay) - #include "MCHAlign/AlignmentSpec.h" #include @@ -79,758 +78,729 @@ #include "DetectorsCommonDataFormats/DetectorNameConf.h" #include "MathUtils/Cartesian.h" - - -namespace o2{ -namespace mch{ +namespace o2 +{ +namespace mch +{ using namespace std; using namespace o2::framework; -using namespace o2; +using namespace o2; class AlignmentTask { -public: - double Reso_X{0.4}; - double Reso_Y{0.4}; - double ImproveCut{6.0}; - - int tracksGood = 0; - int tracksGoodwithoutFit = 0; - int tracksAll = 0; - int trackMCHMID = 0; + public: + double Reso_X{0.4}; + double Reso_Y{0.4}; + double ImproveCut{6.0}; - const int fgNDetElemCh[10] = {4, 4, 4, 4, 18, 18, 26, 26, 26, 26}; - const int fgSNDetElemCh[11] = {0, 4, 8, 12, 16, 34, 52, 78, 104, 130, 156}; - const int fgNDetElemHalfCh[20] = {2, 2, 2, 2, 2, 2, 2, 2, 9, - 9, 9, 9, 13, 13, 13, 13, 13, 13, 13, 13}; - const int fgSNDetElemHalfCh[21] = {0, 3, 6, 9, 12, 15, 18, 21, 24, 34, 44, 54, 64, - 78, 92, 106, 120, 134, 148, 162, 176}; - const int fgDetElemHalfCh[20][13] = - { - {100, 103, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, - {101, 102, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + int tracksGood = 0; + int tracksGoodwithoutFit = 0; + int tracksAll = 0; + int trackMCHMID = 0; - {200, 203, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, - {201, 202, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + const int fgNDetElemCh[10] = {4, 4, 4, 4, 18, 18, 26, 26, 26, 26}; + const int fgSNDetElemCh[11] = {0, 4, 8, 12, 16, 34, 52, 78, 104, 130, 156}; + const int fgNDetElemHalfCh[20] = {2, 2, 2, 2, 2, 2, 2, 2, 9, + 9, 9, 9, 13, 13, 13, 13, 13, 13, 13, 13}; + const int fgSNDetElemHalfCh[21] = {0, 3, 6, 9, 12, 15, 18, 21, 24, 34, 44, 54, 64, + 78, 92, 106, 120, 134, 148, 162, 176}; + const int fgDetElemHalfCh[20][13] = + { + {100, 103, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {101, 102, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, - {300, 303, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, - {301, 302, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {200, 203, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {201, 202, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, - {400, 403, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, - {401, 402, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {300, 303, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {301, 302, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, - {500, 501, 502, 503, 504, 514, 515, 516, 517, 0, 0, 0, 0}, - {505, 506, 507, 508, 509, 510, 511, 512, 513, 0, 0, 0, 0}, + {400, 403, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {401, 402, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, - {600, 601, 602, 603, 604, 614, 615, 616, 617, 0, 0, 0, 0}, - {605, 606, 607, 608, 609, 610, 611, 612, 613, 0, 0, 0, 0}, + {500, 501, 502, 503, 504, 514, 515, 516, 517, 0, 0, 0, 0}, + {505, 506, 507, 508, 509, 510, 511, 512, 513, 0, 0, 0, 0}, - {700, 701, 702, 703, 704, 705, 706, 720, 721, 722, 723, 724, 725}, - {707, 708, 709, 710, 711, 712, 713, 714, 715, 716, 717, 718, 719}, + {600, 601, 602, 603, 604, 614, 615, 616, 617, 0, 0, 0, 0}, + {605, 606, 607, 608, 609, 610, 611, 612, 613, 0, 0, 0, 0}, - {800, 801, 802, 803, 804, 805, 806, 820, 821, 822, 823, 824, 825}, - {807, 808, 809, 810, 811, 812, 813, 814, 815, 816, 817, 818, 819}, + {700, 701, 702, 703, 704, 705, 706, 720, 721, 722, 723, 724, 725}, + {707, 708, 709, 710, 711, 712, 713, 714, 715, 716, 717, 718, 719}, - {900, 901, 902, 903, 904, 905, 906, 920, 921, 922, 923, 924, 925}, - {907, 908, 909, 910, 911, 912, 913, 914, 915, 916, 917, 918, 919}, + {800, 801, 802, 803, 804, 805, 806, 820, 821, 822, 823, 824, 825}, + {807, 808, 809, 810, 811, 812, 813, 814, 815, 816, 817, 818, 819}, - {1000, 1001, 1002, 1003, 1004, 1005, 1006, 1020, 1021, 1022, 1023, 1024, 1025}, - {1007, 1008, 1009, 1010, 1011, 1012, 1013, 1014, 1015, 1016, 1017, 1018, 1019} + {900, 901, 902, 903, 904, 905, 906, 920, 921, 922, 923, 924, 925}, + {907, 908, 909, 910, 911, 912, 913, 914, 915, 916, 917, 918, 919}, - }; + {1000, 1001, 1002, 1003, 1004, 1005, 1006, 1020, 1021, 1022, 1023, 1024, 1025}, + {1007, 1008, 1009, 1010, 1011, 1012, 1013, 1014, 1015, 1016, 1017, 1018, 1019} + }; - double params[624]; - double errors[624]; - double pulls[624]; + double params[624]; + double errors[624]; + double pulls[624]; - constexpr double pi() { return 3.14159265358979323846; } + constexpr double pi() { return 3.14159265358979323846; } - //_________________________________________________________________________________________________ - AlignmentTask(shared_ptr req):mCCDBRequest(req){} + //_________________________________________________________________________________________________ + AlignmentTask(shared_ptr req) : mCCDBRequest(req) {} - //_________________________________________________________________________________________________ - void init(framework::InitContext& ic){ + //_________________________________________________________________________________________________ + void init(framework::InitContext& ic) + { - LOG(info) << "Initializing aligner"; - // Initialize alignment algorithm + LOG(info) << "Initializing aligner"; + // Initialize alignment algorithm - doAlign = ic.options().get("do-align"); - if(doAlign){ - LOG(info) << "Alignment mode"; - }else{ - LOG(info) << "No alignment mode, only residuals will be stored"; - } + doAlign = ic.options().get("do-align"); + if (doAlign) { + LOG(info) << "Alignment mode"; + } else { + LOG(info) << "No alignment mode, only residuals will be stored"; + } - doReAlign = ic.options().get("do-realign"); + doReAlign = ic.options().get("do-realign"); if (doReAlign) { LOG(info) << "Re-alignment mode"; } auto param_config = ic.options().get("fitter-config"); - if(param_config == "PbPb"){ - Reso_X = 0.2; - Reso_Y = 0.2; - ImproveCut = 4.0; - LOG(info) << "Using PbPb parameter set for TrackFitter"; - }else if(param_config == "pp"){ - Reso_X = 0.4; - Reso_Y = 0.4; - ImproveCut = 6.0; - LOG(info) << "Using pp parameter set for TrackFitter"; - }else{ - LOG(fatal) << "Please enter a correct parameter configuration option"; - exit(-1); - } - - - if (mCCDBRequest) { - LOG(info) << "Loading magnetic field and reference geometry from CCDB"; - base::GRPGeomHelper::instance().setRequest(mCCDBRequest); - } else { - - LOG(info) << "Loading magnetic field and reference geometry from input files"; - - auto grpFile = ic.options().get("grp-file"); - if (filesystem::exists(grpFile)) { - const auto grp = parameters::GRPObject::loadFrom(grpFile); - base::Propagator::initFieldFromGRP(grp); - TrackExtrap::setField(); - mAlign.SetBFieldOn(TrackExtrap::isFieldON()); - TrackExtrap::useExtrapV2(); - trackFitter.initField(grp->getL3Current(), grp->getDipoleCurrent()); - }else{ - LOG(fatal) << "No GRP file"; - } - - auto geoIdealFile = ic.options().get("geo-file-ideal"); - if (filesystem::exists(geoIdealFile)) { - base::GeometryManager::loadGeometry(geoIdealFile.c_str()); - transformation = geo::transformationFromTGeoManager(*gGeoManager); - for (int i = 0; i < 156; i++) { - int iDEN = GetDetElemId(i); - transformIdeal[iDEN] = transformation(iDEN); - } - } else { - LOG(fatal) << "No ideal geometry"; - } - - - auto geoRefFile = ic.options().get("geo-file-ref"); - if (filesystem::exists(geoRefFile)) { - base::GeometryManager::loadGeometry(geoRefFile.c_str()); - transformation = geo::transformationFromTGeoManager(*gGeoManager); - for (int i = 0; i < 156; i++) { - int iDEN = GetDetElemId(i); - transformRef[iDEN] = transformation(iDEN); - } - } else { - LOG(fatal) << "No reference geometry"; - } - } - - trackFitter.smoothTracks(true); - trackFitter.setChamberResolution(Reso_X, Reso_Y); - trackFitter.useChamberResolution(); - - mAlign.SetDoEvaluation(true); - // Variation range for parameters - mAlign.SetAllowedVariation(0, 2.0); - mAlign.SetAllowedVariation(1, 0.3); - mAlign.SetAllowedVariation(2, 0.002); - mAlign.SetAllowedVariation(3, 2.0); - - // Fix chambers - auto chambers = ic.options().get("fix-chamber"); - for (int i = 0; i < chambers.length(); ++i) { + if (param_config == "PbPb") { + Reso_X = 0.2; + Reso_Y = 0.2; + ImproveCut = 4.0; + LOG(info) << "Using PbPb parameter set for TrackFitter"; + } else if (param_config == "pp") { + Reso_X = 0.4; + Reso_Y = 0.4; + ImproveCut = 6.0; + LOG(info) << "Using pp parameter set for TrackFitter"; + } else { + LOG(fatal) << "Please enter a correct parameter configuration option"; + exit(-1); + } + + if (mCCDBRequest) { + LOG(info) << "Loading magnetic field and reference geometry from CCDB"; + base::GRPGeomHelper::instance().setRequest(mCCDBRequest); + } else { + + LOG(info) << "Loading magnetic field and reference geometry from input files"; + + auto grpFile = ic.options().get("grp-file"); + if (filesystem::exists(grpFile)) { + const auto grp = parameters::GRPObject::loadFrom(grpFile); + base::Propagator::initFieldFromGRP(grp); + TrackExtrap::setField(); + mAlign.SetBFieldOn(TrackExtrap::isFieldON()); + TrackExtrap::useExtrapV2(); + trackFitter.initField(grp->getL3Current(), grp->getDipoleCurrent()); + } else { + LOG(fatal) << "No GRP file"; + } + + auto geoIdealFile = ic.options().get("geo-file-ideal"); + if (filesystem::exists(geoIdealFile)) { + base::GeometryManager::loadGeometry(geoIdealFile.c_str()); + transformation = geo::transformationFromTGeoManager(*gGeoManager); + for (int i = 0; i < 156; i++) { + int iDEN = GetDetElemId(i); + transformIdeal[iDEN] = transformation(iDEN); + } + } else { + LOG(fatal) << "No ideal geometry"; + } + + auto geoRefFile = ic.options().get("geo-file-ref"); + if (filesystem::exists(geoRefFile)) { + base::GeometryManager::loadGeometry(geoRefFile.c_str()); + transformation = geo::transformationFromTGeoManager(*gGeoManager); + for (int i = 0; i < 156; i++) { + int iDEN = GetDetElemId(i); + transformRef[iDEN] = transformation(iDEN); + } + } else { + LOG(fatal) << "No reference geometry"; + } + } + + trackFitter.smoothTracks(true); + trackFitter.setChamberResolution(Reso_X, Reso_Y); + trackFitter.useChamberResolution(); + + mAlign.SetDoEvaluation(true); + // Variation range for parameters + mAlign.SetAllowedVariation(0, 2.0); + mAlign.SetAllowedVariation(1, 0.3); + mAlign.SetAllowedVariation(2, 0.002); + mAlign.SetAllowedVariation(3, 2.0); + + // Fix chambers + auto chambers = ic.options().get("fix-chamber"); + for (int i = 0; i < chambers.length(); ++i) { if (chambers[i] == ',') { continue; } int chamber = chambers[i] - '0'; - LOG(info) << Form("%s%d","Fixing chamber: ",chamber); - mAlign.FixChamber(chamber); - } - + LOG(info) << Form("%s%d", "Fixing chamber: ", chamber); + mAlign.FixChamber(chamber); + } - doMatched = ic.options().get("matched"); - outFileName = ic.options().get("output"); - readFromRec = ic.options().get("use-record"); + doMatched = ic.options().get("matched"); + outFileName = ic.options().get("output"); + readFromRec = ic.options().get("use-record"); if (readFromRec) { LOG(info) << "Reading records as input"; } mAlign.init("recDataFile.root", "recConsFile.root", readFromRec); - ic.services().get().set([this](){ - LOG(info) << "Alignment duration = " << mElapsedTime.count() << " s"; - }); - - - } - - //_________________________________________________________________________________________________ - void finaliseCCDB(framework::ConcreteDataMatcher& matcher, void* obj) - { - /// finalize the track extrapolation setting - if (mCCDBRequest && base::GRPGeomHelper::instance().finaliseCCDB(matcher, obj)) { - - if (matcher == framework::ConcreteDataMatcher("GLO", "GRPMAGFIELD", 0)) { - LOG(info) << "Loading magnetic field from CCDB"; - auto grp = base::GRPGeomHelper::instance().getGRPMagField(); - base::Propagator::initFieldFromGRP(grp); - TrackExtrap::setField(); - mAlign.SetBFieldOn(TrackExtrap::isFieldON()); - TrackExtrap::useExtrapV2(); - trackFitter.initField(grp->getL3Current(), grp->getDipoleCurrent()); - } - - - if (matcher == framework::ConcreteDataMatcher("GLO", "GEOMALIGN", 0)) { - LOG(info) << "Loading reference geometry from CCDB"; - transformation = geo::transformationFromTGeoManager(*gGeoManager); - for (int i = 0; i < 156; i++) { - int iDEN = GetDetElemId(i); - transformRef[iDEN] = transformation(iDEN); - } - } - } - } - - //_________________________________________________________________________________________________ - void processWithMatching(vector &mchROFs, vector &mchTracks, vector &mchClusters, - vector &muonTracks) - { - // processing for each track - for (const auto &mchROF : mchROFs) { - - for (int iMCHTrack = mchROF.getFirstIdx(); - iMCHTrack <= mchROF.getLastIdx(); ++iMCHTrack) { - tracksAll +=1; - // MCH-MID matching + ic.services().get().set([this]() { + LOG(info) << "Alignment duration = " << mElapsedTime.count() << " s"; + }); + } + + //_________________________________________________________________________________________________ + void finaliseCCDB(framework::ConcreteDataMatcher& matcher, void* obj) + { + /// finalize the track extrapolation setting + if (mCCDBRequest && base::GRPGeomHelper::instance().finaliseCCDB(matcher, obj)) { + + if (matcher == framework::ConcreteDataMatcher("GLO", "GRPMAGFIELD", 0)) { + LOG(info) << "Loading magnetic field from CCDB"; + auto grp = base::GRPGeomHelper::instance().getGRPMagField(); + base::Propagator::initFieldFromGRP(grp); + TrackExtrap::setField(); + mAlign.SetBFieldOn(TrackExtrap::isFieldON()); + TrackExtrap::useExtrapV2(); + trackFitter.initField(grp->getL3Current(), grp->getDipoleCurrent()); + } + + if (matcher == framework::ConcreteDataMatcher("GLO", "GEOMALIGN", 0)) { + LOG(info) << "Loading reference geometry from CCDB"; + transformation = geo::transformationFromTGeoManager(*gGeoManager); + for (int i = 0; i < 156; i++) { + int iDEN = GetDetElemId(i); + transformRef[iDEN] = transformation(iDEN); + } + } + } + } + + //_________________________________________________________________________________________________ + void processWithMatching(vector& mchROFs, vector& mchTracks, vector& mchClusters, + vector& muonTracks) + { + // processing for each track + for (const auto& mchROF : mchROFs) { + + for (int iMCHTrack = mchROF.getFirstIdx(); + iMCHTrack <= mchROF.getLastIdx(); ++iMCHTrack) { + tracksAll += 1; + // MCH-MID matching if (!FindMuon(iMCHTrack, muonTracks)) { continue; } - trackMCHMID += 1; + trackMCHMID += 1; - auto mchTrack = mchTracks.at(iMCHTrack); - int id_track = iMCHTrack; - int nb_clusters = mchTrack.getNClusters(); + auto mchTrack = mchTracks.at(iMCHTrack); + int id_track = iMCHTrack; + int nb_clusters = mchTrack.getNClusters(); - // Track selection, saving only tracks having exactly 10 clusters + // Track selection, saving only tracks having exactly 10 clusters if (nb_clusters <= 9) { continue; } tracksGoodwithoutFit += 1; - // Format conversion from TrackMCH to Track(MCH internal use) - mch::Track convertedTrack = MCHFormatConvert(mchTrack, mchClusters, doReAlign); + // Format conversion from TrackMCH to Track(MCH internal use) + mch::Track convertedTrack = MCHFormatConvert(mchTrack, mchClusters, doReAlign); - // Erase removable track - if(RemoveTrack(convertedTrack, ImproveCut)){ - continue; - }else{ - tracksGood += 1; - } + // Erase removable track + if (RemoveTrack(convertedTrack, ImproveCut)) { + continue; + } else { + tracksGood += 1; + } - // Track processing, saving residuals + // Track processing, saving residuals mAlign.ProcessTrack(convertedTrack, transformation, doAlign, weightRecord); } - } - } + } + } - //_________________________________________________________________________________________________ - void processWithOutMatching(vector &mchROFs, vector &mchTracks, vector &mchClusters) - { + //_________________________________________________________________________________________________ + void processWithOutMatching(vector& mchROFs, vector& mchTracks, vector& mchClusters) + { - // processing for each track - for (const auto &mchROF : mchROFs) { + // processing for each track + for (const auto& mchROF : mchROFs) { - for (int iMCHTrack = mchROF.getFirstIdx(); - iMCHTrack <= mchROF.getLastIdx(); ++iMCHTrack) { + for (int iMCHTrack = mchROF.getFirstIdx(); + iMCHTrack <= mchROF.getLastIdx(); ++iMCHTrack) { - auto mchTrack = mchTracks.at(iMCHTrack); - int id_track = iMCHTrack; - int nb_clusters = mchTrack.getNClusters(); - tracksAll +=1; + auto mchTrack = mchTracks.at(iMCHTrack); + int id_track = iMCHTrack; + int nb_clusters = mchTrack.getNClusters(); + tracksAll += 1; - // Track selection, saving only tracks having exactly 10 clusters + // Track selection, saving only tracks having exactly 10 clusters if (nb_clusters <= 9) { continue; } tracksGoodwithoutFit += 1; - // Format conversion from TrackMCH to Track(MCH internal use) - Track convertedTrack = MCHFormatConvert(mchTrack, mchClusters, doReAlign); + // Format conversion from TrackMCH to Track(MCH internal use) + Track convertedTrack = MCHFormatConvert(mchTrack, mchClusters, doReAlign); - // Erase removable track - if(RemoveTrack(convertedTrack, ImproveCut)){ - continue; - }else{ - tracksGood += 1; - } + // Erase removable track + if (RemoveTrack(convertedTrack, ImproveCut)) { + continue; + } else { + tracksGood += 1; + } - // Track processing, saving residuals + // Track processing, saving residuals mAlign.ProcessTrack(convertedTrack, transformation, doAlign, weightRecord); } - } - } - - //_________________________________________________________________________________________________ - void run(framework::ProcessingContext& pc) - { - auto tStart = std::chrono::high_resolution_clock::now(); - LOG(info) << "Starting alignment process"; + } + } + + //_________________________________________________________________________________________________ + void run(framework::ProcessingContext& pc) + { + auto tStart = std::chrono::high_resolution_clock::now(); + LOG(info) << "Starting alignment process"; + if (doMatched) { + LOG(info) << "Using MCH-MID matched tracks"; + } + if (mCCDBRequest) { + + LOG(info) << "Checking CCDB updates with processing context"; + base::GRPGeomHelper::instance().checkUpdates(pc); + + auto geoIdeal = pc.inputs().get("geomIdeal"); + LOG(info) << "Loading ideal geometry from CCDB"; + transformation = geo::transformationFromTGeoManager(*geoIdeal); + for (int i = 0; i < 156; i++) { + int iDEN = GetDetElemId(i); + transformIdeal[iDEN] = transformation(iDEN); + } + } + + // Load new geometry if we need to do re-align + if (doReAlign) { + if (NewGeoFileName != "") { + LOG(info) << "Loading re-alignment geometry"; + base::GeometryManager::loadGeometry(NewGeoFileName.c_str()); + transformation = geo::transformationFromTGeoManager(*gGeoManager); + for (int i = 0; i < 156; i++) { + int iDEN = GetDetElemId(i); + transformNew[iDEN] = transformation(iDEN); + } + } else { + LOG(fatal) << "No re-alignment geometry"; + } + } + + if (!readFromRec) { + // Loading input data + LOG(info) << "Loading MCH tracks"; + auto [fMCH, mchReader] = LoadData(mchFileName, "o2sim"); + TTreeReaderValue> mchROFs = {*mchReader, "trackrofs"}; + TTreeReaderValue> mchTracks = {*mchReader, "tracks"}; + TTreeReaderValue> mchClusters = {*mchReader, "trackclusters"}; + if (doMatched) { - LOG(info) << "Using MCH-MID matched tracks"; + LOG(info) << "Loading MID tracks"; + auto [fMUON, muonReader] = LoadData(muonFileName.c_str(), "o2sim"); + TTreeReaderValue> muonTracks = {*muonReader, "tracks"}; + int nTF = muonReader->GetEntries(false); + if (mchReader->GetEntries(false) != nTF) { + LOG(error) << mchFileName << " and " << muonFileName << " do not contain the same number of TF"; + exit(-1); + } + + LOG(info) << "Starting track processing"; + while (mchReader->Next() && muonReader->Next()) { + int id_event = mchReader->GetCurrentEntry(); + processWithMatching(*mchROFs, *mchTracks, *mchClusters, *muonTracks); + } + } else { + LOG(info) << "Starting track processing"; + while (mchReader->Next()) { + int id_event = mchReader->GetCurrentEntry(); + processWithOutMatching(*mchROFs, *mchTracks, *mchClusters); + } } - if (mCCDBRequest) { - - LOG(info) << "Checking CCDB updates with processing context"; - base::GRPGeomHelper::instance().checkUpdates(pc); - - auto geoIdeal = pc.inputs().get("geomIdeal"); - LOG(info) << "Loading ideal geometry from CCDB"; - transformation = geo::transformationFromTGeoManager(*geoIdeal); - for (int i = 0; i < 156; i++) { - int iDEN = GetDetElemId(i); - transformIdeal[iDEN] = transformation(iDEN); - } - - } - - // Load new geometry if we need to do re-align - if(doReAlign){ - if(NewGeoFileName!=""){ - LOG(info) << "Loading re-alignment geometry"; - base::GeometryManager::loadGeometry(NewGeoFileName.c_str()); - transformation = geo::transformationFromTGeoManager(*gGeoManager); - for (int i = 0; i < 156; i++) { - int iDEN = GetDetElemId(i); - transformNew[iDEN] = transformation(iDEN); - } - }else{ - LOG(fatal) << "No re-alignment geometry"; - } - } - - if(!readFromRec){ - // Loading input data - LOG(info) << "Loading MCH tracks"; - auto [fMCH, mchReader] = LoadData(mchFileName, "o2sim"); - TTreeReaderValue> mchROFs = {*mchReader, "trackrofs"}; - TTreeReaderValue> mchTracks = {*mchReader,"tracks"}; - TTreeReaderValue> mchClusters = {*mchReader,"trackclusters"}; - - - if(doMatched){ - LOG(info) << "Loading MID tracks"; - auto [fMUON, muonReader] = LoadData(muonFileName.c_str(), "o2sim"); - TTreeReaderValue> muonTracks = {*muonReader, "tracks"}; - int nTF = muonReader->GetEntries(false); - if (mchReader->GetEntries(false) != nTF) { - LOG(error) << mchFileName << " and " << muonFileName << " do not contain the same number of TF"; - exit(-1); - } - - LOG(info) << "Starting track processing"; - while(mchReader->Next() && muonReader->Next()){ - int id_event = mchReader->GetCurrentEntry(); - processWithMatching(*mchROFs, *mchTracks, *mchClusters, *muonTracks); - } - }else{ - LOG(info) << "Starting track processing"; - while(mchReader->Next()){ - int id_event = mchReader->GetCurrentEntry(); - processWithOutMatching(*mchROFs, *mchTracks, *mchClusters); - } - } - - } - - - // Global fit + } + + // Global fit if (doAlign) { mAlign.GlobalFit(params, errors, pulls); } auto tEnd = std::chrono::high_resolution_clock::now(); - mElapsedTime = tEnd - tStart; - // Evaluation for track removing and selection - LOG(info) << Form("%s%d", "Number of good tracks used in alignment process: ",tracksGood); - LOG(info) << Form("%s%d", "Number of good tracks without fit processing: ",tracksGoodwithoutFit); - LOG(info) << Form("%s%d", "Number of MCH-MID tracks: ",trackMCHMID); - LOG(info) << Form("%s%d","Total number of tracks loaded: ", tracksAll); - LOG(info) << Form("%s%f","Ratio of MCH-MID track: ", double(trackMCHMID)/tracksAll); - LOG(info) << Form("%s%f","Ratio before fit: ", double(tracksGoodwithoutFit)/tracksAll); - LOG(info) << Form("%s%f","Ratio after fit: ", double(tracksGood)/tracksAll); - - // Generate new geometry w.r.t alignment results - if(doAlign){ - - LOG(info) << "Generating aligned geometry using global parameters"; - vector ParamAligned; - mAlign.ReAlign(ParamAligned, params); - - TFile *FileAlign = TFile::Open("AlignParam.root", "RECREATE"); - FileAlign->cd(); - FileAlign->WriteObjectAny(&ParamAligned, "std::vector", "alignment"); - FileAlign->Close(); - - string Geo_file; - - if(doReAlign){ - Geo_file = Form("%s%s", "o2sim_geometry_ReAlign", ".root"); - }else{ - Geo_file = Form("%s%s", "o2sim_geometry_Align", ".root"); - } + mElapsedTime = tEnd - tStart; + // Evaluation for track removing and selection + LOG(info) << Form("%s%d", "Number of good tracks used in alignment process: ", tracksGood); + LOG(info) << Form("%s%d", "Number of good tracks without fit processing: ", tracksGoodwithoutFit); + LOG(info) << Form("%s%d", "Number of MCH-MID tracks: ", trackMCHMID); + LOG(info) << Form("%s%d", "Total number of tracks loaded: ", tracksAll); + LOG(info) << Form("%s%f", "Ratio of MCH-MID track: ", double(trackMCHMID) / tracksAll); + LOG(info) << Form("%s%f", "Ratio before fit: ", double(tracksGoodwithoutFit) / tracksAll); + LOG(info) << Form("%s%f", "Ratio after fit: ", double(tracksGood) / tracksAll); + + // Generate new geometry w.r.t alignment results + if (doAlign) { - // Store aligned geometry - gGeoManager->Export(Geo_file.c_str()); + LOG(info) << "Generating aligned geometry using global parameters"; + vector ParamAligned; + mAlign.ReAlign(ParamAligned, params); - // Store param plots - drawHisto(params, errors, pulls, *(mAlign.GetResTree()), outFileName); + TFile* FileAlign = TFile::Open("AlignParam.root", "RECREATE"); + FileAlign->cd(); + FileAlign->WriteObjectAny(&ParamAligned, "std::vector", "alignment"); + FileAlign->Close(); - // Export align params in ideal frame - TransRef(ParamAligned); + string Geo_file; + if (doReAlign) { + Geo_file = Form("%s%s", "o2sim_geometry_ReAlign", ".root"); + } else { + Geo_file = Form("%s%s", "o2sim_geometry_Align", ".root"); + } - } + // Store aligned geometry + gGeoManager->Export(Geo_file.c_str()); - mAlign.terminate(); + // Store param plots + drawHisto(params, errors, pulls, *(mAlign.GetResTree()), outFileName); - pc.services().get().endOfStream(); - pc.services().get().readyToQuit(QuitRequest::Me); + // Export align params in ideal frame + TransRef(ParamAligned); + } + mAlign.terminate(); - } + pc.services().get().endOfStream(); + pc.services().get().readyToQuit(QuitRequest::Me); + } + + private: + //_________________________________________________________________________________________________ + void TransRef(vector& ParamsTrack) + { + LOG(info) << "Transforming align params to the frame of ideal geometry"; + vector ParamsRef; + o2::detectors::AlignParam param_Ref; + for (int hc = 0; hc < 20; hc++) { -private: - //_________________________________________________________________________________________________ - void TransRef(vector &ParamsTrack) - { - LOG(info) << "Transforming align params to the frame of ideal geometry"; - vector ParamsRef; - o2::detectors::AlignParam param_Ref; + ParamsRef.emplace_back(ParamsTrack.at(fgSNDetElemHalfCh[hc])); - for (int hc = 0; hc < 20; hc++) { + for (int de = 0; de < fgNDetElemHalfCh[hc]; de++) { - ParamsRef.emplace_back(ParamsTrack.at(fgSNDetElemHalfCh[hc])); + int iDEN = fgDetElemHalfCh[hc][de]; + o2::detectors::AlignParam param_Track = ParamsTrack.at(fgSNDetElemHalfCh[hc] + 1 + de); - for (int de = 0; de < fgNDetElemHalfCh[hc]; de++) { + LOG(debug) << Form("%s%s", "Processing DET Elem: ", (param_Track.getSymName()).c_str()); - int iDEN = fgDetElemHalfCh[hc][de]; - o2::detectors::AlignParam param_Track = ParamsTrack.at(fgSNDetElemHalfCh[hc]+1+de); + TGeoHMatrix delta_track; + TGeoRotation r("Rotation/Track", param_Track.getPsi() / pi() * 180.0, param_Track.getTheta() / pi() * 180.0, param_Track.getPhi() / pi() * 180.0); + delta_track.SetRotation(r.GetRotationMatrix()); + delta_track.SetDx(param_Track.getX()); + delta_track.SetDy(param_Track.getY()); + delta_track.SetDz(param_Track.getZ()); - LOG(debug) << Form("%s%s","Processing DET Elem: ", (param_Track.getSymName()).c_str()); + TGeoHMatrix transRef = transformIdeal[iDEN]; + TGeoHMatrix transTrack = doReAlign ? transformNew[iDEN] : transformRef[iDEN]; + TGeoHMatrix transRefTrack = transTrack * transRef.Inverse(); + TGeoHMatrix delta_ref = delta_track * transRefTrack; - TGeoHMatrix delta_track; - TGeoRotation r("Rotation/Track", param_Track.getPsi()/pi()*180.0, param_Track.getTheta()/pi()*180.0, param_Track.getPhi()/pi()*180.0); - delta_track.SetRotation(r.GetRotationMatrix()); - delta_track.SetDx(param_Track.getX()); - delta_track.SetDy(param_Track.getY()); - delta_track.SetDz(param_Track.getZ()); + param_Ref.setSymName((param_Track.getSymName()).c_str()); + param_Ref.setGlobalParams(delta_ref); + param_Ref.applyToGeometry(); + ParamsRef.emplace_back(param_Ref); + } + } - TGeoHMatrix transRef = transformIdeal[iDEN]; - TGeoHMatrix transTrack = doReAlign ? transformNew[iDEN] : transformRef[iDEN]; - TGeoHMatrix transRefTrack = transTrack * transRef.Inverse(); - TGeoHMatrix delta_ref = delta_track * transRefTrack; + TFile* fOut = TFile::Open("AlignParam@ideal.root", "RECREATE"); + fOut->WriteObjectAny(&ParamsRef, "std::vector", "alignment"); + fOut->Close(); + } - param_Ref.setSymName((param_Track.getSymName()).c_str()); - param_Ref.setGlobalParams(delta_ref); - param_Ref.applyToGeometry(); - ParamsRef.emplace_back(param_Ref); + //_________________________________________________________________________________________________ + Track MCHFormatConvert(TrackMCH& mchTrack, vector& mchClusters, bool doReAlign) + { - } - } + Track convertedTrack = Track(); - TFile *fOut = TFile::Open("AlignParam@ideal.root","RECREATE"); - fOut->WriteObjectAny(&ParamsRef, "std::vector", "alignment"); - fOut->Close(); - } + // Get clusters for current track + int id_cluster_first = mchTrack.getFirstClusterIdx(); + int id_cluster_last = mchTrack.getLastClusterIdx(); + for (int id_cluster = id_cluster_first; + id_cluster < id_cluster_last + 1; ++id_cluster) { + Cluster* cluster = &(mchClusters.at(id_cluster)); + const int DEId_cluster = cluster->getDEId(); + const int CId_cluster = cluster->getChamberId(); + const int ind_cluster = cluster->getClusterIndex(); - //_________________________________________________________________________________________________ - Track MCHFormatConvert(TrackMCH &mchTrack, vector &mchClusters, bool doReAlign) - { - - Track convertedTrack = Track(); + // Transformations to new geometry from reference geometry + if (doReAlign) { - // Get clusters for current track - int id_cluster_first = mchTrack.getFirstClusterIdx(); - int id_cluster_last = mchTrack.getLastClusterIdx(); + math_utils::Point3D local; + math_utils::Point3D master; - for (int id_cluster = id_cluster_first; - id_cluster < id_cluster_last + 1; ++id_cluster) { + master.SetXYZ(cluster->getX(), cluster->getY(), cluster->getZ()); + transformRef[cluster->getDEId()].MasterToLocal(master, local); + transformNew[cluster->getDEId()].LocalToMaster(local, master); - Cluster *cluster = &(mchClusters.at(id_cluster)); - const int DEId_cluster = cluster->getDEId(); - const int CId_cluster = cluster->getChamberId(); - const int ind_cluster = cluster->getClusterIndex(); + cluster->x = master.x(); + cluster->y = master.y(); + cluster->z = master.z(); + } + convertedTrack.createParamAtCluster(*cluster); + } - // Transformations to new geometry from reference geometry - if(doReAlign){ + return Track(convertedTrack); + } - math_utils::Point3D local; - math_utils::Point3D master; + //_________________________________________________________________________________________________ + bool RemoveTrack(Track& track, double ImproveCut) + { - - master.SetXYZ(cluster->getX(), cluster->getY(), cluster->getZ()); + double maxChi2Cluster = 2 * ImproveCut * ImproveCut; + bool removeTrack = false; - transformRef[cluster->getDEId()].MasterToLocal(master, local); - transformNew[cluster->getDEId()].LocalToMaster(local, master); - - cluster->x = master.x(); - cluster->y = master.y(); - cluster->z = master.z(); + try { + trackFitter.fit(track, false); + } catch (exception const& e) { + removeTrack = true; + return removeTrack; + } - } - convertedTrack.createParamAtCluster(*cluster); + auto itStartingParam = std::prev(track.rend()); - } + while (true) { - return Track(convertedTrack); + try { + trackFitter.fit(track, true, false, (itStartingParam == track.rbegin()) ? nullptr : &itStartingParam); + } catch (exception const&) { + removeTrack = true; + break; + } - } + double worstLocalChi2 = -1.0; - //_________________________________________________________________________________________________ - bool RemoveTrack(Track &track, double ImproveCut) - { + track.tagRemovableClusters(0x1F, false); - double maxChi2Cluster = 2*ImproveCut*ImproveCut; - bool removeTrack = false; + auto itWorstParam = track.end(); - try{ - trackFitter.fit(track, false); - }catch(exception const& e){ - removeTrack = true; - return removeTrack; - } + for (auto itParam = track.begin(); itParam != track.end(); ++itParam) { + if (itParam->getLocalChi2() > worstLocalChi2) { + worstLocalChi2 = itParam->getLocalChi2(); + itWorstParam = itParam; + } + } - auto itStartingParam = std::prev(track.rend()); + if (worstLocalChi2 < maxChi2Cluster) { + break; + } - while(true){ + if (!itWorstParam->isRemovable()) { + removeTrack = true; + track.removable(); + break; + } - try { - trackFitter.fit(track, true, false, (itStartingParam == track.rbegin()) ? nullptr : &itStartingParam); - } catch (exception const&) { - removeTrack = true; - break; - } + auto itNextParam = track.removeParamAtCluster(itWorstParam); + auto itNextToNextParam = (itNextParam == track.end()) ? itNextParam : std::next(itNextParam); + itStartingParam = track.rbegin(); - double worstLocalChi2 = -1.0; + if (track.getNClusters() < 10) { + removeTrack = true; + break; + } else { + while (itNextToNextParam != track.end()) { + if (itNextToNextParam->getClusterPtr()->getChamberId() != itNextParam->getClusterPtr()->getChamberId()) { + itStartingParam = std::make_reverse_iterator(++itNextParam); + break; + } + ++itNextToNextParam; + } + } + } - track.tagRemovableClusters(0x1F, false); + if (!removeTrack) { + for (auto& param : track) { + param.setParameters(param.getSmoothParameters()); + param.setCovariances(param.getSmoothCovariances()); + } + } - auto itWorstParam = track.end(); + return removeTrack; + } - for(auto itParam = track.begin(); itParam != track.end(); ++itParam){ - if(itParam->getLocalChi2() > worstLocalChi2){ - worstLocalChi2 = itParam->getLocalChi2(); - itWorstParam = itParam; - } - } + //_________________________________________________________________________________________________ + void drawHisto(double* params, double* errors, double* pulls, TTree& Res_Tree, string outFileName) + { - if (worstLocalChi2 < maxChi2Cluster) { - break; + TH1F* hPullX = new TH1F("hPullX", "hPullX", 201, -10, 10); + TH1F* hPullY = new TH1F("hPullY", "hPullY", 201, -10, 10); + TH1F* hPullZ = new TH1F("hPullZ", "hPullZ", 201, -10, 10); + TH1F* hPullPhi = new TH1F("hPullPhi", "hPullPhi", 201, -10, 10); + + double deNumber[156]; + + double alignX[156]; + double alignY[156]; + double alignZ[156]; + double alignPhi[156]; + double pullX[156]; + double pullY[156]; + double pullZ[156]; + double pullPhi[156]; + + for (int iDEN = 0; iDEN < 156; iDEN++) { + deNumber[iDEN] = iDEN + 0.5; + alignX[iDEN] = params[iDEN * 4]; + alignY[iDEN] = params[iDEN * 4 + 1]; + alignZ[iDEN] = params[iDEN * 4 + 3]; + alignPhi[iDEN] = params[iDEN * 4 + 2]; + pullX[iDEN] = pulls[iDEN * 4]; + pullY[iDEN] = pulls[iDEN * 4 + 1]; + pullZ[iDEN] = pulls[iDEN * 4 + 3]; + pullPhi[iDEN] = pulls[iDEN * 4 + 2]; + if (params[iDEN * 4]) { + + hPullX->Fill(pulls[iDEN * 4]); + hPullY->Fill(pulls[iDEN * 4 + 1]); + hPullZ->Fill(pulls[iDEN * 4 + 3]); + hPullPhi->Fill(pulls[iDEN * 4 + 2]); } + } - if(!itWorstParam->isRemovable()){ - removeTrack = true; - track.removable(); - break; - } - - - auto itNextParam = track.removeParamAtCluster(itWorstParam); - auto itNextToNextParam = (itNextParam == track.end()) ? itNextParam : std::next(itNextParam); - itStartingParam = track.rbegin(); - - if(track.getNClusters()<10){ - removeTrack = true; - break; - }else{ - while (itNextToNextParam != track.end()) { - if (itNextToNextParam->getClusterPtr()->getChamberId() != itNextParam->getClusterPtr()->getChamberId()) { - itStartingParam = std::make_reverse_iterator(++itNextParam); - break; - } - ++itNextToNextParam; - } - } - - - } - - if(!removeTrack){ - for (auto& param : track) { - param.setParameters(param.getSmoothParameters()); - param.setCovariances(param.getSmoothCovariances()); - } - } - - return removeTrack; - - } - - //_________________________________________________________________________________________________ - void drawHisto(double *params, double *errors, double *pulls, TTree &Res_Tree, string outFileName) - { - - TH1F *hPullX = new TH1F("hPullX", "hPullX", 201, -10, 10); - TH1F *hPullY = new TH1F("hPullY", "hPullY", 201, -10, 10); - TH1F *hPullZ = new TH1F("hPullZ", "hPullZ", 201, -10, 10); - TH1F *hPullPhi = new TH1F("hPullPhi", "hPullPhi", 201, -10, 10); - - double deNumber[156]; - - double alignX[156]; - double alignY[156]; - double alignZ[156]; - double alignPhi[156]; - double pullX[156]; - double pullY[156]; - double pullZ[156]; - double pullPhi[156]; - - for (int iDEN = 0; iDEN < 156; iDEN++) { - deNumber[iDEN] = iDEN + 0.5; - alignX[iDEN] = params[iDEN * 4]; - alignY[iDEN] = params[iDEN * 4 + 1]; - alignZ[iDEN] = params[iDEN * 4 + 3]; - alignPhi[iDEN] = params[iDEN * 4 + 2]; - pullX[iDEN] = pulls[iDEN * 4]; - pullY[iDEN] = pulls[iDEN * 4 + 1]; - pullZ[iDEN] = pulls[iDEN * 4 + 3]; - pullPhi[iDEN] = pulls[iDEN * 4 + 2]; - if (params[iDEN * 4]) { - - hPullX->Fill(pulls[iDEN * 4]); - hPullY->Fill(pulls[iDEN * 4 + 1]); - hPullZ->Fill(pulls[iDEN * 4 + 3]); - hPullPhi->Fill(pulls[iDEN * 4 + 2]); - } - } - - TGraph *graphAlignX = new TGraph(156, deNumber, alignX); - TGraph *graphAlignY = new TGraph(156, deNumber, alignY); - TGraph *graphAlignZ = new TGraph(156, deNumber, alignZ); - TGraph *graphAlignPhi = new TGraph(156, deNumber, alignPhi); - TGraph *graphAlignYZ = new TGraph(156, alignY, alignZ); - - TGraph *graphPullX = new TGraph(156, deNumber, pullX); - TGraph *graphPullY = new TGraph(156, deNumber, pullY); - TGraph *graphPullZ = new TGraph(156, deNumber, pullZ); - TGraph *graphPullPhi = new TGraph(156, deNumber, pullPhi); - - - graphAlignX->SetMarkerStyle(24); - graphPullX->SetMarkerStyle(25); - - // graphAlignX->Draw("AP"); - - graphAlignY->SetMarkerStyle(24); - graphPullY->SetMarkerStyle(25); - - // graphAlignY->Draw("Psame"); - - graphAlignZ->SetMarkerStyle(24); - graphPullZ->SetMarkerStyle(25); - - // graphAlignZ->Draw("AP"); - graphAlignPhi->SetMarkerStyle(24); - graphPullPhi->SetMarkerStyle(25); - - graphAlignYZ->SetMarkerStyle(24); - // graphAlignYZ->Draw("AP"); - - // Saving plots - string PlotFiles_name = Form("%s%s",outFileName.c_str(),"_results.root"); - TFile *PlotFiles = TFile::Open(PlotFiles_name.c_str(),"RECREATE"); - PlotFiles->WriteObjectAny(hPullX,"TH1F","hPullX"); - PlotFiles->WriteObjectAny(hPullY,"TH1F","hPullY"); - PlotFiles->WriteObjectAny(hPullZ,"TH1F","hPullZ"); - PlotFiles->WriteObjectAny(hPullPhi,"TH1F","hPullPhi"); - PlotFiles->WriteObjectAny(graphAlignX,"TGraph","graphAlignX"); - PlotFiles->WriteObjectAny(graphAlignY,"TGraph","graphAlignY"); - PlotFiles->WriteObjectAny(graphAlignZ,"TGraph","graphAlignZ"); - PlotFiles->WriteObjectAny(graphAlignYZ,"TGraph","graphAlignYZ"); - - - TCanvas *cvn1 = new TCanvas("cvn1", "cvn1", 1200, 1600); - //cvn1->Draw(); - cvn1->Divide(1, 4); - TLine limLine(4, -5, 4, 5); - TH1F *aHisto = new TH1F("aHisto", "AlignParam", 161, 0, 160); - aHisto->SetXTitle("Det. Elem. Number"); - for (int i = 1; i < 5; i++) { - cvn1->cd(i); - double Range[4] = {5.0, 1.0, 5.0, 0.01}; - switch (i) { - case 1: - aHisto->SetYTitle("#delta_{#X} (cm)"); - aHisto->GetYaxis()->SetRangeUser(-5.0, 5.0); - aHisto->DrawCopy("goff"); - graphAlignX->Draw("Psame goff"); - limLine.DrawLine(4, -Range[i-1], 4, Range[i-1]); - limLine.DrawLine(8, -Range[i-1], 8, Range[i-1]); - limLine.DrawLine(12, -Range[i-1], 12, Range[i-1]); - limLine.DrawLine(16, -Range[i-1], 16, Range[i-1]); - limLine.DrawLine(16 + 18, -Range[i-1], 16 + 18, Range[i-1]); - limLine.DrawLine(16 + 2 * 18, -Range[i-1], 16 + 2 * 18, Range[i-1]); - limLine.DrawLine(16 + 2 * 18 + 26, -Range[i-1], 16 + 2 * 18 + 26, Range[i-1]); - limLine.DrawLine(16 + 2 * 18 + 2 * 26, -Range[i-1], 16 + 2 * 18 + 2 * 26, Range[i-1]); - limLine.DrawLine(16 + 2 * 18 + 3 * 26, -Range[i-1], 16 + 2 * 18 + 3 * 26, Range[i-1]); - break; - case 2: - aHisto->SetYTitle("#delta_{#Y} (cm)"); - aHisto->GetYaxis()->SetRangeUser(-1.0, 1.0); - aHisto->DrawCopy("goff"); - graphAlignY->Draw("Psame goff"); - limLine.DrawLine(4, -Range[i-1], 4, Range[i-1]); - limLine.DrawLine(8, -Range[i-1], 8, Range[i-1]); - limLine.DrawLine(12, -Range[i-1], 12, Range[i-1]); - limLine.DrawLine(16, -Range[i-1], 16, Range[i-1]); - limLine.DrawLine(16 + 18, -Range[i-1], 16 + 18, Range[i-1]); - limLine.DrawLine(16 + 2 * 18, -Range[i-1], 16 + 2 * 18, Range[i-1]); - limLine.DrawLine(16 + 2 * 18 + 26, -Range[i-1], 16 + 2 * 18 + 26, Range[i-1]); - limLine.DrawLine(16 + 2 * 18 + 2 * 26, -Range[i-1], 16 + 2 * 18 + 2 * 26, Range[i-1]); - limLine.DrawLine(16 + 2 * 18 + 3 * 26, -Range[i-1], 16 + 2 * 18 + 3 * 26, Range[i-1]); - break; - case 3: - aHisto->SetYTitle("#delta_{#Z} (cm)"); - aHisto->GetYaxis()->SetRangeUser(-5.0, 5.0); - aHisto->DrawCopy("goff"); - graphAlignZ->Draw("Psame goff"); - limLine.DrawLine(4, -Range[i-1], 4, Range[i-1]); - limLine.DrawLine(8, -Range[i-1], 8, Range[i-1]); - limLine.DrawLine(12, -Range[i-1], 12, Range[i-1]); - limLine.DrawLine(16, -Range[i-1], 16, Range[i-1]); - limLine.DrawLine(16 + 18, -Range[i-1], 16 + 18, Range[i-1]); - limLine.DrawLine(16 + 2 * 18, -Range[i-1], 16 + 2 * 18, Range[i-1]); - limLine.DrawLine(16 + 2 * 18 + 26, -Range[i-1], 16 + 2 * 18 + 26, Range[i-1]); - limLine.DrawLine(16 + 2 * 18 + 2 * 26, -Range[i-1], 16 + 2 * 18 + 2 * 26, Range[i-1]); - limLine.DrawLine(16 + 2 * 18 + 3 * 26, -Range[i-1], 16 + 2 * 18 + 3 * 26, Range[i-1]); - break; - case 4: - aHisto->SetYTitle("#delta_{#varphi} (cm)"); - aHisto->GetYaxis()->SetRangeUser(-0.01, 0.01); - aHisto->DrawCopy("goff"); - graphAlignPhi->Draw("Psame goff"); - limLine.DrawLine(4, -Range[i-1], 4, Range[i-1]); - limLine.DrawLine(8, -Range[i-1], 8, Range[i-1]); - limLine.DrawLine(12, -Range[i-1], 12, Range[i-1]); - limLine.DrawLine(16, -Range[i-1], 16, Range[i-1]); - limLine.DrawLine(16 + 18, -Range[i-1], 16 + 18, Range[i-1]); - limLine.DrawLine(16 + 2 * 18, -Range[i-1], 16 + 2 * 18, Range[i-1]); - limLine.DrawLine(16 + 2 * 18 + 26, -Range[i-1], 16 + 2 * 18 + 26, Range[i-1]); - limLine.DrawLine(16 + 2 * 18 + 2 * 26, -Range[i-1], 16 + 2 * 18 + 2 * 26, Range[i-1]); - limLine.DrawLine(16 + 2 * 18 + 3 * 26, -Range[i-1], 16 + 2 * 18 + 3 * 26, Range[i-1]); - break; - } - } + TGraph* graphAlignX = new TGraph(156, deNumber, alignX); + TGraph* graphAlignY = new TGraph(156, deNumber, alignY); + TGraph* graphAlignZ = new TGraph(156, deNumber, alignZ); + TGraph* graphAlignPhi = new TGraph(156, deNumber, alignPhi); + TGraph* graphAlignYZ = new TGraph(156, alignY, alignZ); + + TGraph* graphPullX = new TGraph(156, deNumber, pullX); + TGraph* graphPullY = new TGraph(156, deNumber, pullY); + TGraph* graphPullZ = new TGraph(156, deNumber, pullZ); + TGraph* graphPullPhi = new TGraph(156, deNumber, pullPhi); + + graphAlignX->SetMarkerStyle(24); + graphPullX->SetMarkerStyle(25); + + // graphAlignX->Draw("AP"); + + graphAlignY->SetMarkerStyle(24); + graphPullY->SetMarkerStyle(25); + + // graphAlignY->Draw("Psame"); + + graphAlignZ->SetMarkerStyle(24); + graphPullZ->SetMarkerStyle(25); + + // graphAlignZ->Draw("AP"); + graphAlignPhi->SetMarkerStyle(24); + graphPullPhi->SetMarkerStyle(25); + + graphAlignYZ->SetMarkerStyle(24); + // graphAlignYZ->Draw("AP"); + + // Saving plots + string PlotFiles_name = Form("%s%s", outFileName.c_str(), "_results.root"); + TFile* PlotFiles = TFile::Open(PlotFiles_name.c_str(), "RECREATE"); + PlotFiles->WriteObjectAny(hPullX, "TH1F", "hPullX"); + PlotFiles->WriteObjectAny(hPullY, "TH1F", "hPullY"); + PlotFiles->WriteObjectAny(hPullZ, "TH1F", "hPullZ"); + PlotFiles->WriteObjectAny(hPullPhi, "TH1F", "hPullPhi"); + PlotFiles->WriteObjectAny(graphAlignX, "TGraph", "graphAlignX"); + PlotFiles->WriteObjectAny(graphAlignY, "TGraph", "graphAlignY"); + PlotFiles->WriteObjectAny(graphAlignZ, "TGraph", "graphAlignZ"); + PlotFiles->WriteObjectAny(graphAlignYZ, "TGraph", "graphAlignYZ"); + + TCanvas* cvn1 = new TCanvas("cvn1", "cvn1", 1200, 1600); + // cvn1->Draw(); + cvn1->Divide(1, 4); + TLine limLine(4, -5, 4, 5); + TH1F* aHisto = new TH1F("aHisto", "AlignParam", 161, 0, 160); + aHisto->SetXTitle("Det. Elem. Number"); + for (int i = 1; i < 5; i++) { + cvn1->cd(i); + double Range[4] = {5.0, 1.0, 5.0, 0.01}; + switch (i) { + case 1: + aHisto->SetYTitle("#delta_{#X} (cm)"); + aHisto->GetYaxis()->SetRangeUser(-5.0, 5.0); + aHisto->DrawCopy("goff"); + graphAlignX->Draw("Psame goff"); + limLine.DrawLine(4, -Range[i - 1], 4, Range[i - 1]); + limLine.DrawLine(8, -Range[i - 1], 8, Range[i - 1]); + limLine.DrawLine(12, -Range[i - 1], 12, Range[i - 1]); + limLine.DrawLine(16, -Range[i - 1], 16, Range[i - 1]); + limLine.DrawLine(16 + 18, -Range[i - 1], 16 + 18, Range[i - 1]); + limLine.DrawLine(16 + 2 * 18, -Range[i - 1], 16 + 2 * 18, Range[i - 1]); + limLine.DrawLine(16 + 2 * 18 + 26, -Range[i - 1], 16 + 2 * 18 + 26, Range[i - 1]); + limLine.DrawLine(16 + 2 * 18 + 2 * 26, -Range[i - 1], 16 + 2 * 18 + 2 * 26, Range[i - 1]); + limLine.DrawLine(16 + 2 * 18 + 3 * 26, -Range[i - 1], 16 + 2 * 18 + 3 * 26, Range[i - 1]); + break; + case 2: + aHisto->SetYTitle("#delta_{#Y} (cm)"); + aHisto->GetYaxis()->SetRangeUser(-1.0, 1.0); + aHisto->DrawCopy("goff"); + graphAlignY->Draw("Psame goff"); + limLine.DrawLine(4, -Range[i - 1], 4, Range[i - 1]); + limLine.DrawLine(8, -Range[i - 1], 8, Range[i - 1]); + limLine.DrawLine(12, -Range[i - 1], 12, Range[i - 1]); + limLine.DrawLine(16, -Range[i - 1], 16, Range[i - 1]); + limLine.DrawLine(16 + 18, -Range[i - 1], 16 + 18, Range[i - 1]); + limLine.DrawLine(16 + 2 * 18, -Range[i - 1], 16 + 2 * 18, Range[i - 1]); + limLine.DrawLine(16 + 2 * 18 + 26, -Range[i - 1], 16 + 2 * 18 + 26, Range[i - 1]); + limLine.DrawLine(16 + 2 * 18 + 2 * 26, -Range[i - 1], 16 + 2 * 18 + 2 * 26, Range[i - 1]); + limLine.DrawLine(16 + 2 * 18 + 3 * 26, -Range[i - 1], 16 + 2 * 18 + 3 * 26, Range[i - 1]); + break; + case 3: + aHisto->SetYTitle("#delta_{#Z} (cm)"); + aHisto->GetYaxis()->SetRangeUser(-5.0, 5.0); + aHisto->DrawCopy("goff"); + graphAlignZ->Draw("Psame goff"); + limLine.DrawLine(4, -Range[i - 1], 4, Range[i - 1]); + limLine.DrawLine(8, -Range[i - 1], 8, Range[i - 1]); + limLine.DrawLine(12, -Range[i - 1], 12, Range[i - 1]); + limLine.DrawLine(16, -Range[i - 1], 16, Range[i - 1]); + limLine.DrawLine(16 + 18, -Range[i - 1], 16 + 18, Range[i - 1]); + limLine.DrawLine(16 + 2 * 18, -Range[i - 1], 16 + 2 * 18, Range[i - 1]); + limLine.DrawLine(16 + 2 * 18 + 26, -Range[i - 1], 16 + 2 * 18 + 26, Range[i - 1]); + limLine.DrawLine(16 + 2 * 18 + 2 * 26, -Range[i - 1], 16 + 2 * 18 + 2 * 26, Range[i - 1]); + limLine.DrawLine(16 + 2 * 18 + 3 * 26, -Range[i - 1], 16 + 2 * 18 + 3 * 26, Range[i - 1]); + break; + case 4: + aHisto->SetYTitle("#delta_{#varphi} (cm)"); + aHisto->GetYaxis()->SetRangeUser(-0.01, 0.01); + aHisto->DrawCopy("goff"); + graphAlignPhi->Draw("Psame goff"); + limLine.DrawLine(4, -Range[i - 1], 4, Range[i - 1]); + limLine.DrawLine(8, -Range[i - 1], 8, Range[i - 1]); + limLine.DrawLine(12, -Range[i - 1], 12, Range[i - 1]); + limLine.DrawLine(16, -Range[i - 1], 16, Range[i - 1]); + limLine.DrawLine(16 + 18, -Range[i - 1], 16 + 18, Range[i - 1]); + limLine.DrawLine(16 + 2 * 18, -Range[i - 1], 16 + 2 * 18, Range[i - 1]); + limLine.DrawLine(16 + 2 * 18 + 26, -Range[i - 1], 16 + 2 * 18 + 26, Range[i - 1]); + limLine.DrawLine(16 + 2 * 18 + 2 * 26, -Range[i - 1], 16 + 2 * 18 + 2 * 26, Range[i - 1]); + limLine.DrawLine(16 + 2 * 18 + 3 * 26, -Range[i - 1], 16 + 2 * 18 + 3 * 26, Range[i - 1]); + break; + } + } int RefClDetElem; int RefClDetElemNumber; @@ -840,32 +810,31 @@ class AlignmentTask float RefTrackY; float RefTrackSlopeX; float RefTrackSlopeY; - Res_Tree.SetBranchAddress("fClusterX",&RefClusterX); - Res_Tree.SetBranchAddress("fClusterY",&RefClusterY); - Res_Tree.SetBranchAddress("fTrackX",&RefTrackX); - Res_Tree.SetBranchAddress("fTrackY",&RefTrackY); - Res_Tree.SetBranchAddress("fTrackSlopeX",&RefTrackSlopeX); - Res_Tree.SetBranchAddress("fTrackSlopeY",&RefTrackSlopeY); - Res_Tree.SetBranchAddress("fClDetElem",&RefClDetElem); - Res_Tree.SetBranchAddress("fClDetElemNumber",&RefClDetElemNumber); - - TH1F *Histos_Res[2][11]; - for(int i=0;i<2;i++){ - for(int j=0;j<11;j++){ - if(i==0){ - Histos_Res[i][j] = new TH1F(Form("%s%d","Residual_X_Ch",j),Form("%s%d","Residual_x_Ch",j),200,-5,5); - } - - if(i==1){ - Histos_Res[i][j] = new TH1F(Form("%s%d","Residual_Y_Ch",j),Form("%s%d","Residual_y_Ch",j),200,-5,5); - } - } - } - - - TH1F *Histos_DETRes[2][156]; - for(int i=0;i<2;i++){ - for(int j=0;j<156;j++){ + Res_Tree.SetBranchAddress("fClusterX", &RefClusterX); + Res_Tree.SetBranchAddress("fClusterY", &RefClusterY); + Res_Tree.SetBranchAddress("fTrackX", &RefTrackX); + Res_Tree.SetBranchAddress("fTrackY", &RefTrackY); + Res_Tree.SetBranchAddress("fTrackSlopeX", &RefTrackSlopeX); + Res_Tree.SetBranchAddress("fTrackSlopeY", &RefTrackSlopeY); + Res_Tree.SetBranchAddress("fClDetElem", &RefClDetElem); + Res_Tree.SetBranchAddress("fClDetElemNumber", &RefClDetElemNumber); + + TH1F* Histos_Res[2][11]; + for (int i = 0; i < 2; i++) { + for (int j = 0; j < 11; j++) { + if (i == 0) { + Histos_Res[i][j] = new TH1F(Form("%s%d", "Residual_X_Ch", j), Form("%s%d", "Residual_x_Ch", j), 200, -5, 5); + } + + if (i == 1) { + Histos_Res[i][j] = new TH1F(Form("%s%d", "Residual_Y_Ch", j), Form("%s%d", "Residual_y_Ch", j), 200, -5, 5); + } + } + } + + TH1F* Histos_DETRes[2][156]; + for (int i = 0; i < 2; i++) { + for (int j = 0; j < 156; j++) { if (i == 0) { Histos_DETRes[i][j] = new TH1F(Form("%s%d", "Hist_x_DET", j + 1), Form("%s%d", "Hist_x_DET", j + 1), 200, -5, 5); } @@ -873,39 +842,37 @@ class AlignmentTask Histos_DETRes[i][j] = new TH1F(Form("%s%d", "Hist_y_DET", j + 1), Form("%s%d", "Hist_y_DET", j + 1), 200, -5, 5); } } - } - - int Ref_NbEntries = Res_Tree.GetEntries(); - for(int i=0; i < Ref_NbEntries;i++){ - - Res_Tree.GetEntry(i); - double Res_X = RefClusterX - RefTrackX; - double Res_Y = RefClusterY - RefTrackY; - - for(int iCh=0; iCh<11;iCh++){ - if(iCh==0){ - Histos_Res[0][iCh]->Fill(Res_X); - Histos_Res[1][iCh]->Fill(Res_Y); - }else{ - if(iCh == int(RefClDetElem/100)){ - Histos_Res[0][iCh]->Fill(Res_X); - Histos_Res[1][iCh]->Fill(Res_Y); - } - } - } - - for(int iDEN = 0; iDEN < 156; iDEN++){ - if(RefClDetElemNumber == iDEN){ - Histos_DETRes[0][iDEN]->Fill(Res_X); - Histos_DETRes[1][iDEN]->Fill(Res_Y); - } - - } - - } - - for(int i=0;i<2;i++){ - for(int j=0;j<11;j++){ + } + + int Ref_NbEntries = Res_Tree.GetEntries(); + for (int i = 0; i < Ref_NbEntries; i++) { + + Res_Tree.GetEntry(i); + double Res_X = RefClusterX - RefTrackX; + double Res_Y = RefClusterY - RefTrackY; + + for (int iCh = 0; iCh < 11; iCh++) { + if (iCh == 0) { + Histos_Res[0][iCh]->Fill(Res_X); + Histos_Res[1][iCh]->Fill(Res_Y); + } else { + if (iCh == int(RefClDetElem / 100)) { + Histos_Res[0][iCh]->Fill(Res_X); + Histos_Res[1][iCh]->Fill(Res_Y); + } + } + } + + for (int iDEN = 0; iDEN < 156; iDEN++) { + if (RefClDetElemNumber == iDEN) { + Histos_DETRes[0][iDEN]->Fill(Res_X); + Histos_DETRes[1][iDEN]->Fill(Res_Y); + } + } + } + + for (int i = 0; i < 2; i++) { + for (int j = 0; j < 11; j++) { if (i == 0) { PlotFiles->WriteObjectAny(Histos_Res[i][j], "TH1F", Form("%s%d", "Residual_X_Ch", j)); } @@ -913,193 +880,188 @@ class AlignmentTask PlotFiles->WriteObjectAny(Histos_Res[i][j], "TH1F", Form("%s%d", "Residual_Y_Ch", j)); } } - } - - double ResX_mean[156]={}; - double ResX_err[156]={}; - double ResY_mean[156]={}; - double ResY_err[156]={}; - for(int iDEN = 0; iDEN < 156; iDEN++){ - - ResX_mean[iDEN]=Histos_DETRes[0][iDEN]->GetMean(); - ResX_err[iDEN]=Histos_DETRes[0][iDEN]->GetRMS(); - - ResY_mean[iDEN]=Histos_DETRes[1][iDEN]->GetMean(); - ResY_err[iDEN]=Histos_DETRes[1][iDEN]->GetRMS(); - } - - TGraphErrors *graphResX = new TGraphErrors(156, deNumber, ResX_mean, nullptr, ResX_err); - TGraphErrors *graphResY = new TGraphErrors(156, deNumber, ResY_mean, nullptr, ResY_err); - - TCanvas *graphRes = new TCanvas("graphRes","graphRes",1200, 1600); - TH1F *gHisto = new TH1F("gHisto", "Residuals", 161, 0, 160); - gHisto->SetXTitle("Det. Elem. Number"); - - graphRes->Divide(1,2); - - graphRes->cd(1); - gHisto->SetYTitle("TrackX - ClusterX (cm)"); - gHisto->GetYaxis()->SetRangeUser(-5.0, 5.0); - gHisto->DrawCopy("goff"); - graphResX->SetMarkerStyle(8); - graphResX->SetMarkerSize(0.7); - graphResX->SetLineColor(kBlue); - graphResX->Draw("PZsame goff"); - limLine.DrawLine(4, -5, 4, 5); - limLine.DrawLine(8, -5, 8, 5); - limLine.DrawLine(12, -5, 12, 5); - limLine.DrawLine(16, -5, 16, 5); - limLine.DrawLine(16 + 18, -5, 16 + 18, 5); - limLine.DrawLine(16 + 2 * 18, -5, 16 + 2 * 18, 5); - limLine.DrawLine(16 + 2 * 18 + 26, -5, 16 + 2 * 18 + 26, 5); - limLine.DrawLine(16 + 2 * 18 + 2 * 26, -5, 16 + 2 * 18 + 2 * 26, 5); - limLine.DrawLine(16 + 2 * 18 + 3 * 26, -5, 16 + 2 * 18 + 3 * 26, 5); - - graphRes->cd(2); - gHisto->SetYTitle("TrackY - ClusterY (cm)"); - gHisto->GetYaxis()->SetRangeUser(-5.0, 5.0); - gHisto->DrawCopy("goff"); - graphResY->SetMarkerStyle(8); - graphResY->SetMarkerSize(0.7); - graphResY->SetLineColor(kBlue); - graphResY->Draw("PZsame goff"); - limLine.DrawLine(4, -5, 4, 5); - limLine.DrawLine(8, -5, 8, 5); - limLine.DrawLine(12, -5, 12, 5); - limLine.DrawLine(16, -5, 16, 5); - limLine.DrawLine(16 + 18, -5, 16 + 18, 5); - limLine.DrawLine(16 + 2 * 18, -5, 16 + 2 * 18, 5); - limLine.DrawLine(16 + 2 * 18 + 26, -5, 16 + 2 * 18 + 26, 5); - limLine.DrawLine(16 + 2 * 18 + 2 * 26, -5, 16 + 2 * 18 + 2 * 26, 5); - limLine.DrawLine(16 + 2 * 18 + 3 * 26, -5, 16 + 2 * 18 + 3 * 26, 5); - - - PlotFiles->WriteObjectAny(cvn1,"TCanvas","AlignParam"); - PlotFiles->WriteObjectAny(graphRes,"TCanvas","ResGraph"); - PlotFiles->Close(); - - } - - //_________________________________________________________________________________________________ - tuple LoadData(const string fileName, const string treeName) - { - /// open the input file and get the intput tree - - TFile *f = TFile::Open(fileName.c_str(), "READ"); - if (!f || f->IsZombie()) { - LOG(error) << "Opening file " << fileName << " failed"; - exit(-1); - } - - TTreeReader *r = new TTreeReader(treeName.c_str(), f); - if (r->IsZombie()) { - LOG(error) << "Tree " << treeName << " not found"; - exit(-1); - } - - return std::make_tuple(f, r); - } - - //_________________________________________________________________________________________________ - bool FindMuon(int iMCHTrack, vector &muonTracks) { - /// find the MCH-MID matched track corresponding to this MCH track - for (const auto &muon : muonTracks) { - // cout << "Muon track index: " << muon.getMCHRef().getIndex()<GetMean(); + ResX_err[iDEN] = Histos_DETRes[0][iDEN]->GetRMS(); + + ResY_mean[iDEN] = Histos_DETRes[1][iDEN]->GetMean(); + ResY_err[iDEN] = Histos_DETRes[1][iDEN]->GetRMS(); + } + + TGraphErrors* graphResX = new TGraphErrors(156, deNumber, ResX_mean, nullptr, ResX_err); + TGraphErrors* graphResY = new TGraphErrors(156, deNumber, ResY_mean, nullptr, ResY_err); + + TCanvas* graphRes = new TCanvas("graphRes", "graphRes", 1200, 1600); + TH1F* gHisto = new TH1F("gHisto", "Residuals", 161, 0, 160); + gHisto->SetXTitle("Det. Elem. Number"); + + graphRes->Divide(1, 2); + + graphRes->cd(1); + gHisto->SetYTitle("TrackX - ClusterX (cm)"); + gHisto->GetYaxis()->SetRangeUser(-5.0, 5.0); + gHisto->DrawCopy("goff"); + graphResX->SetMarkerStyle(8); + graphResX->SetMarkerSize(0.7); + graphResX->SetLineColor(kBlue); + graphResX->Draw("PZsame goff"); + limLine.DrawLine(4, -5, 4, 5); + limLine.DrawLine(8, -5, 8, 5); + limLine.DrawLine(12, -5, 12, 5); + limLine.DrawLine(16, -5, 16, 5); + limLine.DrawLine(16 + 18, -5, 16 + 18, 5); + limLine.DrawLine(16 + 2 * 18, -5, 16 + 2 * 18, 5); + limLine.DrawLine(16 + 2 * 18 + 26, -5, 16 + 2 * 18 + 26, 5); + limLine.DrawLine(16 + 2 * 18 + 2 * 26, -5, 16 + 2 * 18 + 2 * 26, 5); + limLine.DrawLine(16 + 2 * 18 + 3 * 26, -5, 16 + 2 * 18 + 3 * 26, 5); + + graphRes->cd(2); + gHisto->SetYTitle("TrackY - ClusterY (cm)"); + gHisto->GetYaxis()->SetRangeUser(-5.0, 5.0); + gHisto->DrawCopy("goff"); + graphResY->SetMarkerStyle(8); + graphResY->SetMarkerSize(0.7); + graphResY->SetLineColor(kBlue); + graphResY->Draw("PZsame goff"); + limLine.DrawLine(4, -5, 4, 5); + limLine.DrawLine(8, -5, 8, 5); + limLine.DrawLine(12, -5, 12, 5); + limLine.DrawLine(16, -5, 16, 5); + limLine.DrawLine(16 + 18, -5, 16 + 18, 5); + limLine.DrawLine(16 + 2 * 18, -5, 16 + 2 * 18, 5); + limLine.DrawLine(16 + 2 * 18 + 26, -5, 16 + 2 * 18 + 26, 5); + limLine.DrawLine(16 + 2 * 18 + 2 * 26, -5, 16 + 2 * 18 + 2 * 26, 5); + limLine.DrawLine(16 + 2 * 18 + 3 * 26, -5, 16 + 2 * 18 + 3 * 26, 5); + + PlotFiles->WriteObjectAny(cvn1, "TCanvas", "AlignParam"); + PlotFiles->WriteObjectAny(graphRes, "TCanvas", "ResGraph"); + PlotFiles->Close(); + } + + //_________________________________________________________________________________________________ + tuple LoadData(const string fileName, const string treeName) + { + /// open the input file and get the intput tree + + TFile* f = TFile::Open(fileName.c_str(), "READ"); + if (!f || f->IsZombie()) { + LOG(error) << "Opening file " << fileName << " failed"; + exit(-1); + } + + TTreeReader* r = new TTreeReader(treeName.c_str(), f); + if (r->IsZombie()) { + LOG(error) << "Tree " << treeName << " not found"; + exit(-1); + } + + return std::make_tuple(f, r); + } + + //_________________________________________________________________________________________________ + bool FindMuon(int iMCHTrack, vector& muonTracks) + { + /// find the MCH-MID matched track corresponding to this MCH track + for (const auto& muon : muonTracks) { + // cout << "Muon track index: " << muon.getMCHRef().getIndex()< 0 && iCh <= 10 && iDet < fgNDetElemCh[iCh - 1])) { - LOG(fatal) << "Invalid detector element id: " << iDetElemId; - } + if (!(iCh > 0 && iCh <= 10 && iDet < fgNDetElemCh[iCh - 1])) { + LOG(fatal) << "Invalid detector element id: " << iDetElemId; + } - // add number of detectors up to this chamber - return iDet + fgSNDetElemCh[iCh - 1]; + // add number of detectors up to this chamber + return iDet + fgSNDetElemCh[iCh - 1]; } //_________________________________________________________________________________________________ int GetDetElemId(int iDetElemNumber) { // make sure detector number is valid - if (!(iDetElemNumber >= fgSNDetElemCh[0] && - iDetElemNumber < fgSNDetElemCh[10])) { - LOG(fatal) << "Invalid detector element number: " << iDetElemNumber; - } - /// get det element number from ID - // get chamber and element number in chamber - int iCh = 0; - int iDet = 0; - for (int i = 1; i <= 10; i++) { - if (iDetElemNumber < fgSNDetElemCh[i]) { - iCh = i; - iDet = iDetElemNumber - fgSNDetElemCh[i - 1]; - break; - } - } - - // make sure detector index is valid - if (!(iCh > 0 && iCh <= 10 && iDet < fgNDetElemCh[iCh - 1])) { - LOG(fatal) << "Invalid detector element id: " << 100 * iCh + iDet; - } - - // add number of detectors up to this chamber - return 100 * iCh + iDet; + if (!(iDetElemNumber >= fgSNDetElemCh[0] && + iDetElemNumber < fgSNDetElemCh[10])) { + LOG(fatal) << "Invalid detector element number: " << iDetElemNumber; + } + /// get det element number from ID + // get chamber and element number in chamber + int iCh = 0; + int iDet = 0; + for (int i = 1; i <= 10; i++) { + if (iDetElemNumber < fgSNDetElemCh[i]) { + iCh = i; + iDet = iDetElemNumber - fgSNDetElemCh[i - 1]; + break; + } + } + + // make sure detector index is valid + if (!(iCh > 0 && iCh <= 10 && iDet < fgNDetElemCh[iCh - 1])) { + LOG(fatal) << "Invalid detector element id: " << 100 * iCh + iDet; + } + + // add number of detectors up to this chamber + return 100 * iCh + iDet; } const string mchFileName{"mchtracks.root"}; - const string muonFileName{"muontracks.root"}; - string outFileName{"Alignment"}; - string RefGeoFileName{""}; - string NewGeoFileName{""}; - bool doAlign{false}; - bool doReAlign{false}; - bool doMatched{false}; - bool readFromRec{false}; + const string muonFileName{"muontracks.root"}; + string outFileName{"Alignment"}; + string RefGeoFileName{""}; + string NewGeoFileName{""}; + bool doAlign{false}; + bool doReAlign{false}; + bool doMatched{false}; + bool readFromRec{false}; const double weightRecord{1.0}; Aligner mAlign{}; shared_ptr mCCDBRequest{}; - map transformRef{}; // reference geometry w.r.t track data - map transformNew{}; // new geometry - map transformIdeal{}; // Ideal geometry - - geo::TransformationCreator transformation{}; - TrackFitter trackFitter{}; - - std::chrono::duration mElapsedTime{}; - + map transformRef{}; // reference geometry w.r.t track data + map transformNew{}; // new geometry + map transformIdeal{}; // Ideal geometry + geo::TransformationCreator transformation{}; + TrackFitter trackFitter{}; + std::chrono::duration mElapsedTime{}; }; //_________________________________________________________________________________________________ o2::framework::DataProcessorSpec getAlignmentSpec(bool disableCCDB) { - vector inputSpecs{}; - inputSpecs.emplace_back("STFDist", "FLP", "DISTSUBTIMEFRAME", 0); - inputSpecs.emplace_back("geomIdeal", "GLO", "GEOMIDEAL", 0, Lifetime::Condition, framework::ccdbParamSpec("GLO/Config/Geometry")); - - vector outputSpecs{}; - auto ccdbRequest = disableCCDB ? nullptr : std::make_shared( false, // orbitResetTime - false, // GRPECS=true - false, // GRPLHCIF - true, // GRPMagField - false, // askMatLUT - base::GRPGeomRequest::Aligned, // geometry - inputSpecs); - + vector inputSpecs{}; + inputSpecs.emplace_back("STFDist", "FLP", "DISTSUBTIMEFRAME", 0); + inputSpecs.emplace_back("geomIdeal", "GLO", "GEOMIDEAL", 0, Lifetime::Condition, framework::ccdbParamSpec("GLO/Config/Geometry")); + + vector outputSpecs{}; + auto ccdbRequest = disableCCDB ? nullptr : std::make_shared(false, // orbitResetTime + false, // GRPECS=true + false, // GRPLHCIF + true, // GRPMagField + false, // askMatLUT + base::GRPGeomRequest::Aligned, // geometry + inputSpecs); return DataProcessorSpec{ "mch-alignment", @@ -1107,7 +1069,7 @@ o2::framework::DataProcessorSpec getAlignmentSpec(bool disableCCDB) outputSpecs, AlgorithmSpec{o2::framework::adaptFromTask(ccdbRequest)}, Options{{"geo-file-ref", VariantType::String, o2::base::NameConf::getAlignedGeomFileName(), {"Name of the reference geometry file"}}, - {"geo-file-ideal", VariantType::String, o2::base::NameConf::getGeomFileName(), {"Name of the ideal geometry file"}}, + {"geo-file-ideal", VariantType::String, o2::base::NameConf::getGeomFileName(), {"Name of the ideal geometry file"}}, {"grp-file", VariantType::String, o2::base::NameConf::getGRPFileName(), {"Name of the grp file"}}, {"fitter-config", VariantType::String, "", {"Option of parameter set for TrackFitter, pp or PbPb"}}, {"do-align", VariantType::Bool, false, {"Switch for alignment, otherwise only residuals will be stored"}}, @@ -1115,9 +1077,8 @@ o2::framework::DataProcessorSpec getAlignmentSpec(bool disableCCDB) {"matched", VariantType::Bool, false, {"Switch for using MCH-MID matched tracks"}}, {"fix-chamber", VariantType::String, "", {"Chamber fixing, ex 1,2,3"}}, {"use-record", VariantType::Bool, false, {"Option for directly using record in alignment if provided"}}, - {"output", VariantType::String, "Alignment", {"Option for name of output file"}}}} ; + {"output", VariantType::String, "Alignment", {"Option for name of output file"}}}}; } - } // namespace mch } // namespace o2 \ No newline at end of file diff --git a/Detectors/MUON/MCH/Align/Aligner/src/alignment-workflow.cxx b/Detectors/MUON/MCH/Align/Aligner/src/alignment-workflow.cxx index 037d5ee46261d..94c858ef18a18 100644 --- a/Detectors/MUON/MCH/Align/Aligner/src/alignment-workflow.cxx +++ b/Detectors/MUON/MCH/Align/Aligner/src/alignment-workflow.cxx @@ -10,7 +10,7 @@ // or submit itself to any jurisdiction. /// \file alignment-workflow.cxx -/// \brief Implementation of a DPL device to perform alignment for muon spectrometer +/// \brief Implementation of a DPL device to perform alignment for muon spectrometer /// \author Chi ZHANG(CEA-Saclay) #include "MCHAlign/AlignmentSpec.h" @@ -32,7 +32,6 @@ #include "Headers/STFHeader.h" #include "DetectorsRaw/HBFUtils.h" - using namespace o2::framework; using namespace std; @@ -62,7 +61,7 @@ class SeederTask : public Task } }; -} +} // namespace o2::mch o2::framework::DataProcessorSpec getSeederSpec() { @@ -74,7 +73,6 @@ o2::framework::DataProcessorSpec getSeederSpec() Options{}}; } - // we need to add workflow options before including Framework/runDataProcessing void customize(vector& workflowOptions) { diff --git a/Detectors/MUON/MCH/Tracking/include/MCHTracking/Track.h b/Detectors/MUON/MCH/Tracking/include/MCHTracking/Track.h index ae02d3a1c4d17..e9b78c176c8e9 100644 --- a/Detectors/MUON/MCH/Tracking/include/MCHTracking/Track.h +++ b/Detectors/MUON/MCH/Tracking/include/MCHTracking/Track.h @@ -104,7 +104,6 @@ class Track int mCurrentChamber = -1; ///< current chamber on which the current parameters are given bool mConnected = false; ///< flag telling if this track shares cluster(s) with another bool mRemovable = false; ///< flag telling if this track should be deleted - }; } // namespace mch