diff --git a/Detectors/CMakeLists.txt b/Detectors/CMakeLists.txt index 8fd00ae6e27d3..a7cbccc7f0be8 100644 --- a/Detectors/CMakeLists.txt +++ b/Detectors/CMakeLists.txt @@ -47,6 +47,8 @@ add_subdirectory(Calibration) add_subdirectory(DCS) add_subdirectory(Align) +add_subdirectory(ForwardAlign) + if(BUILD_SIMULATION) add_subdirectory(gconfig) diff --git a/Detectors/ForwardAlign/CMakeLists.txt b/Detectors/ForwardAlign/CMakeLists.txt new file mode 100644 index 0000000000000..48acb71e2c658 --- /dev/null +++ b/Detectors/ForwardAlign/CMakeLists.txt @@ -0,0 +1,40 @@ +# Copyright 2019-2020 CERN and copyright holders of ALICE O2. +# See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +# All rights not expressly granted are reserved. +# +# This software is distributed under the terms of the GNU General Public +# License v3 (GPL Version 3), copied verbatim in the file "COPYING". +# +# In applying this license CERN does not waive the privileges and immunities +# granted to it by virtue of its status as an Intergovernmental Organization +# or submit itself to any jurisdiction. + +o2_add_library(ForwardAlign + SOURCES src/MatrixSparse.cxx + src/MatrixSq.cxx + src/MillePede2.cxx + src/MillePedeRecord.cxx + src/MilleRecordReader.cxx + src/MilleRecordWriter.cxx + src/MinResSolve.cxx + src/RectMatrix.cxx + src/SymBDMatrix.cxx + src/SymMatrix.cxx + src/VectorSparse.cxx + PUBLIC_LINK_LIBRARIES O2::CCDB + O2::Steer + ROOT::TreePlayer) + +o2_target_root_dictionary(ForwardAlign + HEADERS include/ForwardAlign/MatrixSparse.h + include/ForwardAlign/MatrixSq.h + include/ForwardAlign/MillePede2.h + include/ForwardAlign/MillePedeRecord.h + include/ForwardAlign/MilleRecordReader.h + include/ForwardAlign/MilleRecordWriter.h + include/ForwardAlign/MinResSolve.h + include/ForwardAlign/RectMatrix.h + include/ForwardAlign/SymBDMatrix.h + include/ForwardAlign/SymMatrix.h + include/ForwardAlign/VectorSparse.h + LINKDEF src/ForwardAlignLinkDef.h) diff --git a/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/MatrixSparse.h b/Detectors/ForwardAlign/include/ForwardAlign/MatrixSparse.h similarity index 95% rename from Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/MatrixSparse.h rename to Detectors/ForwardAlign/include/ForwardAlign/MatrixSparse.h index 25ea6da8928a6..6c60b4275b583 100644 --- a/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/MatrixSparse.h +++ b/Detectors/ForwardAlign/include/ForwardAlign/MatrixSparse.h @@ -13,15 +13,15 @@ /// \brief Sparse matrix class (from AliROOT), used as a global matrix for MillePede2 /// \author ruben.shahoyan@cern.ch -#ifndef ALICEO2_MFT_MATRIXSPARSE_H -#define ALICEO2_MFT_MATRIXSPARSE_H +#ifndef ALICEO2_FWDALIGN_MATRIXSPARSE_H +#define ALICEO2_FWDALIGN_MATRIXSPARSE_H -#include "MFTAlignment/MatrixSq.h" -#include "MFTAlignment/VectorSparse.h" +#include "ForwardAlign/MatrixSq.h" +#include "ForwardAlign/VectorSparse.h" namespace o2 { -namespace mft +namespace fwdalign { /// \class MatrixSparse @@ -157,7 +157,7 @@ inline Double_t& MatrixSparse::DiagElem(Int_t row) } } -} // namespace mft +} // namespace fwdalign } // namespace o2 #endif diff --git a/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/MatrixSq.h b/Detectors/ForwardAlign/include/ForwardAlign/MatrixSq.h similarity index 97% rename from Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/MatrixSq.h rename to Detectors/ForwardAlign/include/ForwardAlign/MatrixSq.h index 8c154b359f816..e79fc19004629 100644 --- a/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/MatrixSq.h +++ b/Detectors/ForwardAlign/include/ForwardAlign/MatrixSq.h @@ -13,15 +13,15 @@ /// \brief Abstract class (from AliROOT) for square matrix used for millepede2 operation /// \author ruben.shahoyan@cern.ch -#ifndef ALICEO2_MFT_MATRIXSQ_H -#define ALICEO2_MFT_MATRIXSQ_H +#ifndef ALICEO2_FWDALIGN_MATRIXSQ_H +#define ALICEO2_FWDALIGN_MATRIXSQ_H #include #include namespace o2 { -namespace mft +namespace fwdalign { /// \class MatrixSq @@ -153,7 +153,7 @@ inline void MatrixSq::MultiplyByVec(const TVectorD& vecIn, TVectorD& vecOut) con MultiplyByVec(vecIn.GetMatrixArray(), vecOut.GetMatrixArray()); } -} // namespace mft +} // namespace fwdalign } // namespace o2 #endif diff --git a/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/MillePede2.h b/Detectors/ForwardAlign/include/ForwardAlign/MillePede2.h similarity index 89% rename from Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/MillePede2.h rename to Detectors/ForwardAlign/include/ForwardAlign/MillePede2.h index 1ec7ab7bc8df7..f3707aa64dabf 100644 --- a/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/MillePede2.h +++ b/Detectors/ForwardAlign/include/ForwardAlign/MillePede2.h @@ -16,20 +16,20 @@ /// Based on the original milliped2 by Volker Blobel /// http://www.desy.de/~blobel/mptalks.html -#ifndef ALICEO2_MFT_MILLEPEDE2_H -#define ALICEO2_MFT_MILLEPEDE2_H +#ifndef ALICEO2_FWDALIGN_MILLEPEDE2_H +#define ALICEO2_FWDALIGN_MILLEPEDE2_H #include #include #include -#include "MFTAlignment/MinResSolve.h" -#include "MFTAlignment/MillePedeRecord.h" -#include "MFTAlignment/SymMatrix.h" -#include "MFTAlignment/RectMatrix.h" -#include "MFTAlignment/MatrixSparse.h" -#include "MFTAlignment/MatrixSq.h" -#include "MFTAlignment/MilleRecordWriter.h" -#include "MFTAlignment/MilleRecordReader.h" +#include "ForwardAlign/MinResSolve.h" +#include "ForwardAlign/MillePedeRecord.h" +#include "ForwardAlign/SymMatrix.h" +#include "ForwardAlign/RectMatrix.h" +#include "ForwardAlign/MatrixSparse.h" +#include "ForwardAlign/MatrixSq.h" +#include "ForwardAlign/MilleRecordWriter.h" +#include "ForwardAlign/MilleRecordReader.h" class TFile; class TStopwatch; @@ -38,7 +38,7 @@ class TArrayF; namespace o2 { -namespace mft +namespace fwdalign { class MillePede2 @@ -256,17 +256,17 @@ class MillePede2 /// \brief write tree and close file where are stored chi2 from LocalFit() void EndChi2Storage(); - o2::mft::MillePedeRecord* GetRecord() const { return fRecord; } + o2::fwdalign::MillePedeRecord* GetRecord() const { return fRecord; } long GetSelFirst() const { return fSelFirst; } long GetSelLast() const { return fSelLast; } void SetSelFirst(long v) { fSelFirst = v; } void SetSelLast(long v) { fSelLast = v; } - void SetRecord(o2::mft::MillePedeRecord* aRecord) { fRecord = aRecord; } - void SetRecordWriter(o2::mft::MilleRecordWriter* myP) { fRecordWriter = myP; } - void SetConstraintsRecWriter(o2::mft::MilleRecordWriter* myP) { fConstraintsRecWriter = myP; } - void SetRecordReader(o2::mft::MilleRecordReader* myP) { fRecordReader = myP; } - void SetConstraintsRecReader(o2::mft::MilleRecordReader* myP) { fConstraintsRecReader = myP; } + void SetRecord(o2::fwdalign::MillePedeRecord* aRecord) { fRecord = aRecord; } + void SetRecordWriter(o2::fwdalign::MilleRecordWriter* myP) { fRecordWriter = myP; } + void SetConstraintsRecWriter(o2::fwdalign::MilleRecordWriter* myP) { fConstraintsRecWriter = myP; } + void SetRecordReader(o2::fwdalign::MilleRecordReader* myP) { fRecordReader = myP; } + void SetConstraintsRecReader(o2::fwdalign::MilleRecordReader* myP) { fConstraintsRecReader = myP; } /// \brief return the limit in chi^2/nd for n sigmas stdev authorized /// @@ -340,11 +340,11 @@ class MillePede2 std::vector fCGlo2Glo; ///< [fNGloPar] compressed ID to global ID buffer // Matrices - o2::mft::SymMatrix* fMatCLoc; ///< Matrix C local - o2::mft::MatrixSq* fMatCGlo; ///< Matrix C global - o2::mft::RectMatrix* fMatCGloLoc; ///< Rectangular matrix C g*l - std::vector fFillIndex; ///< [fNGloPar] auxilary index array for fast matrix fill - std::vector fFillValue; ///< [fNGloPar] auxilary value array for fast matrix fill + o2::fwdalign::SymMatrix* fMatCLoc; ///< Matrix C local + o2::fwdalign::MatrixSq* fMatCGlo; ///< Matrix C global + o2::fwdalign::RectMatrix* fMatCGloLoc; ///< Rectangular matrix C g*l + std::vector fFillIndex; ///< [fNGloPar] auxilary index array for fast matrix fill + std::vector fFillValue; ///< [fNGloPar] auxilary value array for fast matrix fill TFile* fRecChi2File; TString fRecChi2FName; @@ -354,7 +354,7 @@ class MillePede2 bool fIsChi2BelowLimit; int fRecNDoF; - o2::mft::MillePedeRecord* fRecord; ///< Buffer of measurements records + o2::fwdalign::MillePedeRecord* fRecord; ///< Buffer of measurements records long fCurrRecDataID; ///< ID of the current data record long fCurrRecConstrID; ///< ID of the current constraint record @@ -380,15 +380,15 @@ class MillePede2 static int fgNKrylovV; ///< size of Krylov vectors buffer in FGMRES // processed data record bufferization - o2::mft::MilleRecordWriter* fRecordWriter; ///< data record writer - o2::mft::MilleRecordWriter* fConstraintsRecWriter; ///< constraints record writer - o2::mft::MilleRecordReader* fRecordReader; ///< data record reader - o2::mft::MilleRecordReader* fConstraintsRecReader; ///< constraints record reader + o2::fwdalign::MilleRecordWriter* fRecordWriter; ///< data record writer + o2::fwdalign::MilleRecordWriter* fConstraintsRecWriter; ///< constraints record writer + o2::fwdalign::MilleRecordReader* fRecordReader; ///< data record reader + o2::fwdalign::MilleRecordReader* fConstraintsRecReader; ///< constraints record reader ClassDef(MillePede2, 0); }; -} // namespace mft +} // namespace fwdalign } // namespace o2 #endif diff --git a/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/MillePedeRecord.h b/Detectors/ForwardAlign/include/ForwardAlign/MillePedeRecord.h similarity index 97% rename from Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/MillePedeRecord.h rename to Detectors/ForwardAlign/include/ForwardAlign/MillePedeRecord.h index ecc28d6f04b98..43dcbf2671c71 100644 --- a/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/MillePedeRecord.h +++ b/Detectors/ForwardAlign/include/ForwardAlign/MillePedeRecord.h @@ -16,14 +16,14 @@ /// The records for all processed tracks are stored in the temporary tree in orgder to be /// reused for multiple iterations of MillePede -#ifndef ALICEO2_MFT_MILLEPEDERECORD_H -#define ALICEO2_MFT_MILLEPEDERECORD_H +#ifndef ALICEO2_FWDALIGN_MILLEPEDERECORD_H +#define ALICEO2_FWDALIGN_MILLEPEDERECORD_H #include namespace o2 { -namespace mft +namespace fwdalign { /// \brief Store residuals and local/global deriavtives from a single track processing /// @@ -152,7 +152,7 @@ inline Bool_t MillePedeRecord::IsGroupPresent(Int_t id) const return kFALSE; } -} // namespace mft +} // namespace fwdalign } // namespace o2 #endif diff --git a/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/MilleRecordReader.h b/Detectors/ForwardAlign/include/ForwardAlign/MilleRecordReader.h similarity index 61% rename from Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/MilleRecordReader.h rename to Detectors/ForwardAlign/include/ForwardAlign/MilleRecordReader.h index 6682f547a6cc3..61c710f509f31 100644 --- a/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/MilleRecordReader.h +++ b/Detectors/ForwardAlign/include/ForwardAlign/MilleRecordReader.h @@ -13,8 +13,8 @@ /// \author arakotoz@cern.ch /// \brief Class dedicated to read MillePedeRecords from ROOT files -#ifndef ALICEO2_MFT_MILLERECORD_READER_H -#define ALICEO2_MFT_MILLERECORD_READER_H +#ifndef ALICEO2_FWDALIGN_MILLERECORD_READER_H +#define ALICEO2_FWDALIGN_MILLERECORD_READER_H #include #include @@ -22,11 +22,11 @@ #include #include -#include "MFTAlignment/MillePedeRecord.h" +#include "ForwardAlign/MillePedeRecord.h" namespace o2 { -namespace mft +namespace fwdalign { class MilleRecordReader @@ -51,7 +51,7 @@ class MilleRecordReader bool isReadEntryOk() const { return mIsReadEntryOk; } /// \brief return the record - o2::mft::MillePedeRecord* getRecord() { return mRecord; }; + o2::fwdalign::MillePedeRecord* getRecord() { return mRecord; }; /// \brief return the ID of the current record in the TTree long getCurrentDataID() const { return mCurrentDataID; } @@ -65,20 +65,23 @@ class MilleRecordReader /// \brief return the number of entries Long64_t getNEntries() const { return mNEntries; } + /// \brief return the name of record data tree + TString getDataTreeName() const { return mDataTreeName; } + protected: - TChain* mDataTree; ///< TChain container that stores the records - bool mIsSuccessfulInit; ///< boolean to monitor the success of the initialization - bool mIsConstraintsRec; ///< boolean to know if these are data records or constraints records - bool mIsReadEntryOk; ///< boolean to know if the last operation readNextEntry() was ok - TString mDataTreeName; ///< name of the record TTree/TChain - TString mDataBranchName; ///< name of the branch where records will be stored - o2::mft::MillePedeRecord* mRecord; ///< the running record - Long64_t mCurrentDataID; ///< counter indicating the ID of the current record in the tree - Long64_t mNEntries; ///< number of entries in the read TChain + TChain* mDataTree; ///< TChain container that stores the records + bool mIsSuccessfulInit; ///< boolean to monitor the success of the initialization + bool mIsConstraintsRec; ///< boolean to know if these are data records or constraints records + bool mIsReadEntryOk; ///< boolean to know if the last operation readNextEntry() was ok + TString mDataTreeName; ///< name of the record TTree/TChain + TString mDataBranchName; ///< name of the branch where records will be stored + o2::fwdalign::MillePedeRecord* mRecord; ///< the running record + Long64_t mCurrentDataID; ///< counter indicating the ID of the current record in the tree + Long64_t mNEntries; ///< number of entries in the read TChain ClassDef(MilleRecordReader, 0); }; -} // namespace mft +} // namespace fwdalign } // namespace o2 #endif diff --git a/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/MilleRecordWriter.h b/Detectors/ForwardAlign/include/ForwardAlign/MilleRecordWriter.h similarity index 60% rename from Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/MilleRecordWriter.h rename to Detectors/ForwardAlign/include/ForwardAlign/MilleRecordWriter.h index 82173111b1acf..7af35092c6ff9 100644 --- a/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/MilleRecordWriter.h +++ b/Detectors/ForwardAlign/include/ForwardAlign/MilleRecordWriter.h @@ -11,21 +11,21 @@ /// \file MilleRecordWriter.h /// \author arakotoz@cern.ch -/// \brief Class dedicated to write MillePedeRecords to output file for MFT +/// \brief Class dedicated to write MillePedeRecords to output file for FWDALIGN -#ifndef ALICEO2_MFT_MILLERECORD_WRITER_H -#define ALICEO2_MFT_MILLERECORD_WRITER_H +#ifndef ALICEO2_FWDALIGN_MILLERECORD_WRITER_H +#define ALICEO2_FWDALIGN_MILLERECORD_WRITER_H #include #include -#include "MFTAlignment/MillePedeRecord.h" +#include "ForwardAlign/MillePedeRecord.h" class TFile; class TTree; namespace o2 { -namespace mft +namespace fwdalign { class MilleRecordWriter @@ -53,7 +53,7 @@ class MilleRecordWriter bool isInitOk() const { return mIsSuccessfulInit; } /// \brief return the record - o2::mft::MillePedeRecord* getRecord() { return mRecord; }; + o2::fwdalign::MillePedeRecord* getRecord() { return mRecord; }; /// \brief return the ID of the current record in the TTree Long64_t getCurrentDataID() const { return mCurrentDataID; } @@ -71,20 +71,20 @@ class MilleRecordWriter void setRecordWeight(double wgh); protected: - TTree* mDataTree; ///< TTree container that stores the records - TFile* mDataFile; ///< output file where the records are written - bool mIsSuccessfulInit; ///< boolean to monitor the success of the initialization - bool mIsConstraintsRec; ///< boolean to know if these are data records or constraints records - long mNEntriesAutoSave; ///< max entries in the buffer after which TTree::AutoSave() is automatically used - TString mDataFileName; ///< name of the output file that will store the record TTree - TString mDataTreeName; ///< name of the record TTree - TString mDataBranchName; ///< name of the branch where records will be stored - o2::mft::MillePedeRecord* mRecord; ///< the running record - Long64_t mCurrentDataID; ///< counter increasing when adding a record to the tree + TTree* mDataTree; ///< TTree container that stores the records + TFile* mDataFile; ///< output file where the records are written + bool mIsSuccessfulInit; ///< boolean to monitor the success of the initialization + bool mIsConstraintsRec; ///< boolean to know if these are data records or constraints records + long mNEntriesAutoSave; ///< max entries in the buffer after which TTree::AutoSave() is automatically used + TString mDataFileName; ///< name of the output file that will store the record TTree + TString mDataTreeName; ///< name of the record TTree + TString mDataBranchName; ///< name of the branch where records will be stored + o2::fwdalign::MillePedeRecord* mRecord; ///< the running record + Long64_t mCurrentDataID; ///< counter increasing when adding a record to the tree ClassDef(MilleRecordWriter, 0); }; -} // namespace mft +} // namespace fwdalign } // namespace o2 #endif diff --git a/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/MinResSolve.h b/Detectors/ForwardAlign/include/ForwardAlign/MinResSolve.h similarity index 97% rename from Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/MinResSolve.h rename to Detectors/ForwardAlign/include/ForwardAlign/MinResSolve.h index 69cbe7faec624..5c0f06e93c49a 100644 --- a/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/MinResSolve.h +++ b/Detectors/ForwardAlign/include/ForwardAlign/MinResSolve.h @@ -13,8 +13,8 @@ /// \brief General class (from AliROOT) for solving large system of linear equations /// \author ruben.shahoyan@cern.ch -#ifndef ALICEO2_MFT_MINRESSOLVE_H -#define ALICEO2_MFT_MINRESSOLVE_H +#ifndef ALICEO2_FWDALIGN_MINRESSOLVE_H +#define ALICEO2_FWDALIGN_MINRESSOLVE_H #include #include @@ -22,7 +22,7 @@ namespace o2 { -namespace mft +namespace fwdalign { class MatrixSq; @@ -135,7 +135,7 @@ class MinResSolve : public TObject ClassDefOverride(MinResSolve, 0); }; -} // namespace mft +} // namespace fwdalign } // namespace o2 #endif diff --git a/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/RectMatrix.h b/Detectors/ForwardAlign/include/ForwardAlign/RectMatrix.h similarity index 94% rename from Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/RectMatrix.h rename to Detectors/ForwardAlign/include/ForwardAlign/RectMatrix.h index 125478f087b61..fae47c5b22eb4 100644 --- a/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/RectMatrix.h +++ b/Detectors/ForwardAlign/include/ForwardAlign/RectMatrix.h @@ -13,15 +13,15 @@ /// \brief Class for rectangular matrix used for millepede2 operation /// \author ruben.shahoyan@cern.ch -#ifndef ALICEO2_MFT_RECTMATRIX_H -#define ALICEO2_MFT_RECTMATRIX_H +#ifndef ALICEO2_FWDALIGN_RECTMATRIX_H +#define ALICEO2_FWDALIGN_RECTMATRIX_H #include "TObject.h" class TString; namespace o2 { -namespace mft +namespace fwdalign { /// \brief Class for rectangular matrix used for millepede2 operation @@ -82,7 +82,7 @@ inline Double_t& RectMatrix::operator()(Int_t row, Int_t col) return (Double_t&)fRows[row][col]; } -} // namespace mft +} // namespace fwdalign } // namespace o2 #endif diff --git a/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/SymBDMatrix.h b/Detectors/ForwardAlign/include/ForwardAlign/SymBDMatrix.h similarity index 97% rename from Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/SymBDMatrix.h rename to Detectors/ForwardAlign/include/ForwardAlign/SymBDMatrix.h index 8ae1056c1c184..27b7e9df0d1ff 100644 --- a/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/SymBDMatrix.h +++ b/Detectors/ForwardAlign/include/ForwardAlign/SymBDMatrix.h @@ -13,16 +13,16 @@ /// \brief Symmetric Band Diagonal matrix (from AliROOT) with half band width W (+1 for diagonal) /// \author ruben.shahoyan@cern.ch -#ifndef ALICEO2_MFT_SYMBDMATRIX_H -#define ALICEO2_MFT_SYMBDMATRIX_H +#ifndef ALICEO2_FWDALIGN_SYMBDMATRIX_H +#define ALICEO2_FWDALIGN_SYMBDMATRIX_H #include #include -#include "MFTAlignment/MatrixSq.h" +#include "ForwardAlign/MatrixSq.h" namespace o2 { -namespace mft +namespace fwdalign { /// \class SymBDMatrix @@ -170,7 +170,7 @@ inline void SymBDMatrix::MultiplyByVec(const TVectorD& vecIn, TVectorD& vecOut) MultiplyByVec(vecIn.GetMatrixArray(), vecOut.GetMatrixArray()); } -} // namespace mft +} // namespace fwdalign } // namespace o2 #endif diff --git a/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/SymMatrix.h b/Detectors/ForwardAlign/include/ForwardAlign/SymMatrix.h similarity index 98% rename from Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/SymMatrix.h rename to Detectors/ForwardAlign/include/ForwardAlign/SymMatrix.h index e86c3bb1f9886..8de3867419202 100644 --- a/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/SymMatrix.h +++ b/Detectors/ForwardAlign/include/ForwardAlign/SymMatrix.h @@ -13,17 +13,17 @@ /// \brief Fast symmetric matrix (from AliROOT) with dynamically expandable size /// \author ruben.shahoyan@cern.ch -#ifndef ALICEO2_MFT_SYMMATRIX_H -#define ALICEO2_MFT_SYMMATRIX_H +#ifndef ALICEO2_FWDALIGN_SYMMATRIX_H +#define ALICEO2_FWDALIGN_SYMMATRIX_H #include #include -#include "MFTAlignment/MatrixSq.h" +#include "ForwardAlign/MatrixSq.h" namespace o2 { -namespace mft +namespace fwdalign { /// \class SymMatrix @@ -277,7 +277,7 @@ inline void SymMatrix::AddToRow(Int_t r, Double_t* valc, Int_t* indc, Int_t n) } } -} // namespace mft +} // namespace fwdalign } // namespace o2 #endif diff --git a/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/VectorSparse.h b/Detectors/ForwardAlign/include/ForwardAlign/VectorSparse.h similarity index 95% rename from Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/VectorSparse.h rename to Detectors/ForwardAlign/include/ForwardAlign/VectorSparse.h index 993ce59014ce7..0e36af5a62c7a 100644 --- a/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/VectorSparse.h +++ b/Detectors/ForwardAlign/include/ForwardAlign/VectorSparse.h @@ -13,15 +13,15 @@ /// \brief Sparse vector class (from AliROOT), used as row of the MatrixSparse class /// \author ruben.shahoyan@cern.ch -#ifndef ALICEO2_MFT_VECTORSPARSE_H -#define ALICEO2_MFT_VECTORSPARSE_H +#ifndef ALICEO2_FWDALIGN_VECTORSPARSE_H +#define ALICEO2_FWDALIGN_VECTORSPARSE_H #include #include namespace o2 { -namespace mft +namespace fwdalign { /// \class VectorSparse @@ -97,7 +97,7 @@ inline Double_t& VectorSparse::operator()(Int_t ind) return FindIndexAdd(ind); } -} // namespace mft +} // namespace fwdalign } // namespace o2 #endif diff --git a/Detectors/ForwardAlign/src/ForwardAlignLinkDef.h b/Detectors/ForwardAlign/src/ForwardAlignLinkDef.h new file mode 100644 index 0000000000000..cfc895a101f8b --- /dev/null +++ b/Detectors/ForwardAlign/src/ForwardAlignLinkDef.h @@ -0,0 +1,30 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#ifdef __CLING__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class o2::fwdalign::MatrixSparse + ; +#pragma link C++ class o2::fwdalign::MatrixSq + ; +#pragma link C++ class o2::fwdalign::MillePede2 + ; +#pragma link C++ class o2::fwdalign::MillePedeRecord + ; +#pragma link C++ class o2::fwdalign::MilleRecordReader + ; +#pragma link C++ class o2::fwdalign::MilleRecordWriter + ; +#pragma link C++ class o2::fwdalign::MinResSolve + ; +#pragma link C++ class o2::fwdalign::RectMatrix + ; +#pragma link C++ class o2::fwdalign::SymBDMatrix + ; +#pragma link C++ class o2::fwdalign::SymMatrix + ; +#pragma link C++ class o2::fwdalign::VectorSparse + ; + +#endif diff --git a/Detectors/ITSMFT/MFT/alignment/src/MatrixSparse.cxx b/Detectors/ForwardAlign/src/MatrixSparse.cxx similarity index 98% rename from Detectors/ITSMFT/MFT/alignment/src/MatrixSparse.cxx rename to Detectors/ForwardAlign/src/MatrixSparse.cxx index acbb7ea94a334..210b0c3a77f5c 100644 --- a/Detectors/ITSMFT/MFT/alignment/src/MatrixSparse.cxx +++ b/Detectors/ForwardAlign/src/MatrixSparse.cxx @@ -15,9 +15,9 @@ #include #include "Framework/Logger.h" -#include "MFTAlignment/MatrixSparse.h" +#include "ForwardAlign/MatrixSparse.h" -using namespace o2::mft; +using namespace o2::fwdalign; ClassImp(MatrixSparse); diff --git a/Detectors/ITSMFT/MFT/alignment/src/MatrixSq.cxx b/Detectors/ForwardAlign/src/MatrixSq.cxx similarity index 96% rename from Detectors/ITSMFT/MFT/alignment/src/MatrixSq.cxx rename to Detectors/ForwardAlign/src/MatrixSq.cxx index abb794def134d..1aadc3b6563d5 100644 --- a/Detectors/ITSMFT/MFT/alignment/src/MatrixSq.cxx +++ b/Detectors/ForwardAlign/src/MatrixSq.cxx @@ -17,9 +17,9 @@ #include "TMath.h" #include "Framework/Logger.h" -#include "MFTAlignment/MatrixSq.h" +#include "ForwardAlign/MatrixSq.h" -using namespace o2::mft; +using namespace o2::fwdalign; ClassImp(MatrixSq); diff --git a/Detectors/ITSMFT/MFT/alignment/src/MillePede2.cxx b/Detectors/ForwardAlign/src/MillePede2.cxx similarity index 99% rename from Detectors/ITSMFT/MFT/alignment/src/MillePede2.cxx rename to Detectors/ForwardAlign/src/MillePede2.cxx index 6220341f20054..d66786d658c5c 100644 --- a/Detectors/ITSMFT/MFT/alignment/src/MillePede2.cxx +++ b/Detectors/ForwardAlign/src/MillePede2.cxx @@ -11,7 +11,7 @@ /// @file MillePede2.cxx -#include "MFTAlignment/MillePede2.h" +#include "ForwardAlign/MillePede2.h" #include "Framework/Logger.h" #include #include @@ -27,14 +27,14 @@ #include #include -//#define _DUMP_EQ_BEFORE_ -//#define _DUMP_EQ_AFTER_ +// #define _DUMP_EQ_BEFORE_ +// #define _DUMP_EQ_AFTER_ -//#define _DUMPEQ_BEFORE_ -//#define _DUMPEQ_AFTER_ +// #define _DUMPEQ_BEFORE_ +// #define _DUMPEQ_AFTER_ using std::ifstream; -using namespace o2::mft; +using namespace o2::fwdalign; ClassImp(MillePede2); diff --git a/Detectors/ITSMFT/MFT/alignment/src/MillePedeRecord.cxx b/Detectors/ForwardAlign/src/MillePedeRecord.cxx similarity index 99% rename from Detectors/ITSMFT/MFT/alignment/src/MillePedeRecord.cxx rename to Detectors/ForwardAlign/src/MillePedeRecord.cxx index d87a7ee9a77d6..447b01bc60035 100644 --- a/Detectors/ITSMFT/MFT/alignment/src/MillePedeRecord.cxx +++ b/Detectors/ForwardAlign/src/MillePedeRecord.cxx @@ -11,11 +11,11 @@ /// @file MillePedeRecord.cxx -#include "MFTAlignment/MillePedeRecord.h" +#include "ForwardAlign/MillePedeRecord.h" #include #include "Framework/Logger.h" -using namespace o2::mft; +using namespace o2::fwdalign; ClassImp(MillePedeRecord); diff --git a/Detectors/ITSMFT/MFT/alignment/src/MilleRecordReader.cxx b/Detectors/ForwardAlign/src/MilleRecordReader.cxx similarity index 96% rename from Detectors/ITSMFT/MFT/alignment/src/MilleRecordReader.cxx rename to Detectors/ForwardAlign/src/MilleRecordReader.cxx index 0cab0726c7e24..a169ba77b6e0d 100644 --- a/Detectors/ITSMFT/MFT/alignment/src/MilleRecordReader.cxx +++ b/Detectors/ForwardAlign/src/MilleRecordReader.cxx @@ -13,11 +13,11 @@ #include "Framework/Logger.h" -#include "MFTAlignment/MilleRecordReader.h" +#include "ForwardAlign/MilleRecordReader.h" -using namespace o2::mft; +using namespace o2::fwdalign; -ClassImp(o2::mft::MilleRecordReader); +ClassImp(o2::fwdalign::MilleRecordReader); //__________________________________________________________________________ MilleRecordReader::MilleRecordReader() diff --git a/Detectors/ITSMFT/MFT/alignment/src/MilleRecordWriter.cxx b/Detectors/ForwardAlign/src/MilleRecordWriter.cxx similarity index 97% rename from Detectors/ITSMFT/MFT/alignment/src/MilleRecordWriter.cxx rename to Detectors/ForwardAlign/src/MilleRecordWriter.cxx index 91b6438d968fd..8dbe31e91744c 100644 --- a/Detectors/ITSMFT/MFT/alignment/src/MilleRecordWriter.cxx +++ b/Detectors/ForwardAlign/src/MilleRecordWriter.cxx @@ -16,11 +16,11 @@ #include "Framework/Logger.h" -#include "MFTAlignment/MilleRecordWriter.h" +#include "ForwardAlign/MilleRecordWriter.h" -using namespace o2::mft; +using namespace o2::fwdalign; -ClassImp(o2::mft::MilleRecordWriter); +ClassImp(o2::fwdalign::MilleRecordWriter); //__________________________________________________________________________ MilleRecordWriter::MilleRecordWriter() diff --git a/Detectors/ITSMFT/MFT/alignment/src/MinResSolve.cxx b/Detectors/ForwardAlign/src/MinResSolve.cxx similarity index 99% rename from Detectors/ITSMFT/MFT/alignment/src/MinResSolve.cxx rename to Detectors/ForwardAlign/src/MinResSolve.cxx index 14eae8262c45d..6d1982307fc1e 100644 --- a/Detectors/ITSMFT/MFT/alignment/src/MinResSolve.cxx +++ b/Detectors/ForwardAlign/src/MinResSolve.cxx @@ -16,12 +16,12 @@ #include #include "Framework/Logger.h" -#include "MFTAlignment/MinResSolve.h" -#include "MFTAlignment/MatrixSq.h" -#include "MFTAlignment/MatrixSparse.h" -#include "MFTAlignment/SymBDMatrix.h" +#include "ForwardAlign/MinResSolve.h" +#include "ForwardAlign/MatrixSq.h" +#include "ForwardAlign/MatrixSparse.h" +#include "ForwardAlign/SymBDMatrix.h" -using namespace o2::mft; +using namespace o2::fwdalign; ClassImp(MinResSolve); diff --git a/Detectors/ITSMFT/MFT/alignment/src/RectMatrix.cxx b/Detectors/ForwardAlign/src/RectMatrix.cxx similarity index 97% rename from Detectors/ITSMFT/MFT/alignment/src/RectMatrix.cxx rename to Detectors/ForwardAlign/src/RectMatrix.cxx index a9b958b1fe308..94eb9fdf70359 100644 --- a/Detectors/ITSMFT/MFT/alignment/src/RectMatrix.cxx +++ b/Detectors/ForwardAlign/src/RectMatrix.cxx @@ -12,9 +12,9 @@ /// @file RectMatrix.cxx #include -#include "MFTAlignment/RectMatrix.h" +#include "ForwardAlign/RectMatrix.h" -using namespace o2::mft; +using namespace o2::fwdalign; ClassImp(RectMatrix); diff --git a/Detectors/ITSMFT/MFT/alignment/src/SymBDMatrix.cxx b/Detectors/ForwardAlign/src/SymBDMatrix.cxx similarity index 99% rename from Detectors/ITSMFT/MFT/alignment/src/SymBDMatrix.cxx rename to Detectors/ForwardAlign/src/SymBDMatrix.cxx index 6dd4129e1d7c7..84d002d0933ea 100644 --- a/Detectors/ITSMFT/MFT/alignment/src/SymBDMatrix.cxx +++ b/Detectors/ForwardAlign/src/SymBDMatrix.cxx @@ -15,9 +15,9 @@ #include "TClass.h" #include "TMath.h" -#include "MFTAlignment/SymBDMatrix.h" +#include "ForwardAlign/SymBDMatrix.h" -using namespace o2::mft; +using namespace o2::fwdalign; ClassImp(SymBDMatrix); diff --git a/Detectors/ITSMFT/MFT/alignment/src/SymMatrix.cxx b/Detectors/ForwardAlign/src/SymMatrix.cxx similarity index 99% rename from Detectors/ITSMFT/MFT/alignment/src/SymMatrix.cxx rename to Detectors/ForwardAlign/src/SymMatrix.cxx index 498cad5d080fc..c5e5267645148 100644 --- a/Detectors/ITSMFT/MFT/alignment/src/SymMatrix.cxx +++ b/Detectors/ForwardAlign/src/SymMatrix.cxx @@ -15,10 +15,10 @@ #include #include -#include "MFTAlignment/SymMatrix.h" +#include "ForwardAlign/SymMatrix.h" #include "Framework/Logger.h" -using namespace o2::mft; +using namespace o2::fwdalign; ClassImp(SymMatrix); diff --git a/Detectors/ITSMFT/MFT/alignment/src/VectorSparse.cxx b/Detectors/ForwardAlign/src/VectorSparse.cxx similarity index 98% rename from Detectors/ITSMFT/MFT/alignment/src/VectorSparse.cxx rename to Detectors/ForwardAlign/src/VectorSparse.cxx index 94463b8650297..b2e7f346b91d2 100644 --- a/Detectors/ITSMFT/MFT/alignment/src/VectorSparse.cxx +++ b/Detectors/ForwardAlign/src/VectorSparse.cxx @@ -13,9 +13,9 @@ #include -#include "MFTAlignment/VectorSparse.h" +#include "ForwardAlign/VectorSparse.h" -using namespace o2::mft; +using namespace o2::fwdalign; ClassImp(VectorSparse); diff --git a/Detectors/ITSMFT/MFT/alignment/CMakeLists.txt b/Detectors/ITSMFT/MFT/alignment/CMakeLists.txt index 7fda70e463936..7b06d70dfd413 100644 --- a/Detectors/ITSMFT/MFT/alignment/CMakeLists.txt +++ b/Detectors/ITSMFT/MFT/alignment/CMakeLists.txt @@ -15,24 +15,14 @@ o2_add_library(MFTAlignment src/AlignPointControl.cxx src/AlignPointHelper.cxx src/AlignSensorHelper.cxx - src/MatrixSparse.cxx - src/MatrixSq.cxx - src/MillePede2.cxx - src/MillePedeRecord.cxx - src/MilleRecordReader.cxx - src/MilleRecordWriter.cxx - src/MinResSolve.cxx src/RecordsToAlignParams.cxx - src/RectMatrix.cxx - src/SymBDMatrix.cxx - src/SymMatrix.cxx src/TracksToRecords.cxx - src/VectorSparse.cxx PUBLIC_LINK_LIBRARIES O2::MFTBase O2::DataFormatsMFT O2::MFTTracking O2::CCDB O2::Steer + O2::ForwardAlign ROOT::TreePlayer) o2_target_root_dictionary(MFTAlignment @@ -41,17 +31,6 @@ o2_target_root_dictionary(MFTAlignment include/MFTAlignment/AlignPointControl.h include/MFTAlignment/AlignPointHelper.h include/MFTAlignment/AlignSensorHelper.h - include/MFTAlignment/MatrixSparse.h - include/MFTAlignment/MatrixSq.h - include/MFTAlignment/MillePede2.h - include/MFTAlignment/MillePedeRecord.h - include/MFTAlignment/MilleRecordReader.h - include/MFTAlignment/MilleRecordWriter.h - include/MFTAlignment/MinResSolve.h include/MFTAlignment/RecordsToAlignParams.h - include/MFTAlignment/RectMatrix.h - include/MFTAlignment/SymBDMatrix.h - include/MFTAlignment/SymMatrix.h include/MFTAlignment/TracksToRecords.h - include/MFTAlignment/VectorSparse.h LINKDEF src/MFTAlignmentLinkDef.h) diff --git a/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/Aligner.h b/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/Aligner.h index 31fb3e7aefb5a..8f68128ece3e8 100644 --- a/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/Aligner.h +++ b/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/Aligner.h @@ -24,7 +24,7 @@ #include #include "ITSMFTReconstruction/ChipMappingMFT.h" -#include "MFTAlignment/MillePede2.h" +#include "ForwardAlign/MillePede2.h" namespace o2 { diff --git a/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/RecordsToAlignParams.h b/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/RecordsToAlignParams.h index 8a6759f491ead..b899357a2a223 100644 --- a/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/RecordsToAlignParams.h +++ b/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/RecordsToAlignParams.h @@ -19,8 +19,8 @@ #include #include -#include "MFTAlignment/MillePede2.h" -#include "MFTAlignment/MilleRecordReader.h" +#include "ForwardAlign/MillePede2.h" +#include "ForwardAlign/MilleRecordReader.h" #include "DetectorsCommonDataFormats/AlignParam.h" #include "MFTAlignment/Aligner.h" @@ -70,16 +70,16 @@ class RecordsToAlignParams : public Aligner void connectConstraintsRecReaderToChain(TChain* ch); protected: - bool mWithControl; ///< boolean to set the use of the control tree = chi2 per track filled by MillePede LocalFit() - long mNEntriesAutoSave = 10000; ///< number of entries needed to cyclically call AutoSave for the output control tree - std::vector mAlignParams; ///< vector of alignment parameters computed by MillePede simultaneous fit - o2::mft::MilleRecordReader* mRecordReader; ///< utility that handles the reading of the data records used to feed MillePede solver - bool mWithConstraintsRecReader; ///< boolean to set to true if one wants to also read constraints records - o2::mft::MilleRecordReader* mConstraintsRecReader; ///< utility that handles the reading of the constraints records - o2::mft::MillePede2* mMillepede; ///< Millepede2 implementation copied from AliROOT - std::vector mPedeOutParams; ///< Vector to store the outputs (alignment corrections) of the MillePede simulatenous fit - std::vector mPedeOutParamsErrors; ///< Vector to store the outputs (errors on the alignement corrections) of the MillePede simulatenous fit - std::vector mPedeOutParamsPulls; ///< Vector to store the outputs (pulls on the alignement corrections) of the MillePede simulatenous fit + bool mWithControl; ///< boolean to set the use of the control tree = chi2 per track filled by MillePede LocalFit() + long mNEntriesAutoSave = 10000; ///< number of entries needed to cyclically call AutoSave for the output control tree + std::vector mAlignParams; ///< vector of alignment parameters computed by MillePede simultaneous fit + o2::fwdalign::MilleRecordReader* mRecordReader; ///< utility that handles the reading of the data records used to feed MillePede solver + bool mWithConstraintsRecReader; ///< boolean to set to true if one wants to also read constraints records + o2::fwdalign::MilleRecordReader* mConstraintsRecReader; ///< utility that handles the reading of the constraints records + o2::fwdalign::MillePede2* mMillepede; ///< Millepede2 implementation copied from AliROOT + std::vector mPedeOutParams; ///< Vector to store the outputs (alignment corrections) of the MillePede simulatenous fit + std::vector mPedeOutParamsErrors; ///< Vector to store the outputs (errors on the alignement corrections) of the MillePede simulatenous fit + std::vector mPedeOutParamsPulls; ///< Vector to store the outputs (pulls on the alignement corrections) of the MillePede simulatenous fit ClassDefOverride(RecordsToAlignParams, 0); }; diff --git a/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/TracksToRecords.h b/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/TracksToRecords.h index 53a7c3ca5a3b6..4a3e3db5d1572 100644 --- a/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/TracksToRecords.h +++ b/Detectors/ITSMFT/MFT/alignment/include/MFTAlignment/TracksToRecords.h @@ -28,8 +28,8 @@ #include "DataFormatsITSMFT/ROFRecord.h" #include "DataFormatsITSMFT/CompCluster.h" #include "ReconstructionDataFormats/BaseCluster.h" -#include "MFTAlignment/MillePedeRecord.h" -#include "MFTAlignment/MilleRecordWriter.h" +#include "ForwardAlign/MillePedeRecord.h" +#include "ForwardAlign/MilleRecordWriter.h" #include "MFTAlignment/AlignPointHelper.h" #include "MFTAlignment/AlignPointControl.h" #include "MFTBase/GeometryTGeo.h" @@ -105,12 +105,12 @@ class TracksToRecords : public Aligner bool mWithControl; ///< boolean to set the use of the control tree long mNEntriesAutoSave = 10000; ///< number of entries needed to call AutoSave for the output TTrees o2::mft::AlignPointControl mPointControl; ///< AlignPointControl handles the control tree - o2::mft::MilleRecordWriter* mRecordWriter; ///< utility that handles the writing of the data records to a ROOT file + o2::fwdalign::MilleRecordWriter* mRecordWriter; ///< utility that handles the writing of the data records to a ROOT file bool mWithConstraintsRecWriter; ///< boolean to be set to true if one wants to also write constaints records - o2::mft::MilleRecordWriter* mConstraintsRecWriter; ///< utility that handles the writing of the constraints records + o2::fwdalign::MilleRecordWriter* mConstraintsRecWriter; ///< utility that handles the writing of the constraints records std::vector> mMFTClustersLocal; ///< MFT clusters in local coordinate system std::vector> mMFTClustersGlobal; ///< MFT clusters in global coordinate system - o2::mft::MillePede2* mMillepede; ///< Millepede2 implementation copied from AliROOT + o2::fwdalign::MillePede2* mMillepede; ///< Millepede2 implementation copied from AliROOT // access these data from CTFs provided uptream by the workflow diff --git a/Detectors/ITSMFT/MFT/alignment/src/Aligner.cxx b/Detectors/ITSMFT/MFT/alignment/src/Aligner.cxx index 8e553aedeb175..fd933ad1f2dce 100644 --- a/Detectors/ITSMFT/MFT/alignment/src/Aligner.cxx +++ b/Detectors/ITSMFT/MFT/alignment/src/Aligner.cxx @@ -19,7 +19,7 @@ #include "Framework/Logger.h" #include "MFTAlignment/Aligner.h" -#include "MFTAlignment/MillePede2.h" +#include "ForwardAlign/MillePede2.h" using namespace o2::mft; diff --git a/Detectors/ITSMFT/MFT/alignment/src/MFTAlignmentLinkDef.h b/Detectors/ITSMFT/MFT/alignment/src/MFTAlignmentLinkDef.h index 5a365cac251f7..e711dff895ae2 100644 --- a/Detectors/ITSMFT/MFT/alignment/src/MFTAlignmentLinkDef.h +++ b/Detectors/ITSMFT/MFT/alignment/src/MFTAlignmentLinkDef.h @@ -20,18 +20,7 @@ #pragma link C++ class o2::mft::AlignPointControl + ; #pragma link C++ class o2::mft::AlignPointHelper + ; #pragma link C++ class o2::mft::AlignSensorHelper + ; -#pragma link C++ class o2::mft::MatrixSparse + ; -#pragma link C++ class o2::mft::MatrixSq + ; -#pragma link C++ class o2::mft::MillePede2 + ; -#pragma link C++ class o2::mft::MillePedeRecord + ; -#pragma link C++ class o2::mft::MilleRecordReader + ; -#pragma link C++ class o2::mft::MilleRecordWriter + ; -#pragma link C++ class o2::mft::MinResSolve + ; #pragma link C++ class o2::mft::RecordsToAlignParams + ; -#pragma link C++ class o2::mft::RectMatrix + ; -#pragma link C++ class o2::mft::SymBDMatrix + ; -#pragma link C++ class o2::mft::SymMatrix + ; #pragma link C++ class o2::mft::TracksToRecords + ; -#pragma link C++ class o2::mft::VectorSparse + ; #endif diff --git a/Detectors/ITSMFT/MFT/alignment/src/RecordsToAlignParams.cxx b/Detectors/ITSMFT/MFT/alignment/src/RecordsToAlignParams.cxx index 7e4078ba8cddb..c632718c93a37 100644 --- a/Detectors/ITSMFT/MFT/alignment/src/RecordsToAlignParams.cxx +++ b/Detectors/ITSMFT/MFT/alignment/src/RecordsToAlignParams.cxx @@ -28,13 +28,13 @@ RecordsToAlignParams::RecordsToAlignParams() : Aligner(), mWithControl(false), mNEntriesAutoSave(10000), - mRecordReader(new MilleRecordReader()), + mRecordReader(new o2::fwdalign::MilleRecordReader()), mWithConstraintsRecReader(false), mConstraintsRecReader(nullptr), - mMillepede(new MillePede2()) + mMillepede(new o2::fwdalign::MillePede2()) { if (mWithConstraintsRecReader) { - mConstraintsRecReader = new MilleRecordReader(); + mConstraintsRecReader = new o2::fwdalign::MilleRecordReader(); } // allocate memory in vectors to store the results of the global fit diff --git a/Detectors/ITSMFT/MFT/alignment/src/TracksToRecords.cxx b/Detectors/ITSMFT/MFT/alignment/src/TracksToRecords.cxx index 0c44a4ad6c9ba..fbf2a7c8a2fee 100644 --- a/Detectors/ITSMFT/MFT/alignment/src/TracksToRecords.cxx +++ b/Detectors/ITSMFT/MFT/alignment/src/TracksToRecords.cxx @@ -18,7 +18,7 @@ #include "Framework/Logger.h" #include #include "MFTBase/Geometry.h" -#include "MFTAlignment/MillePedeRecord.h" +#include "ForwardAlign/MillePedeRecord.h" #include "MFTAlignment/TracksToRecords.h" @@ -45,13 +45,13 @@ TracksToRecords::TracksToRecords() mAlignPoint(new AlignPointHelper()), mWithControl(false), mNEntriesAutoSave(10000), - mRecordWriter(new MilleRecordWriter()), + mRecordWriter(new o2::fwdalign::MilleRecordWriter()), mWithConstraintsRecWriter(false), mConstraintsRecWriter(nullptr), - mMillepede(new MillePede2()) + mMillepede(new o2::fwdalign::MillePede2()) { if (mWithConstraintsRecWriter) { - mConstraintsRecWriter = new MilleRecordWriter(); + mConstraintsRecWriter = new o2::fwdalign::MilleRecordWriter(); } // initialise the content of each array resetGlocalDerivative(); diff --git a/Detectors/MUON/MCH/Align/CMakeLists.txt b/Detectors/MUON/MCH/Align/CMakeLists.txt new file mode 100644 index 0000000000000..4f9c3cd7f35b8 --- /dev/null +++ b/Detectors/MUON/MCH/Align/CMakeLists.txt @@ -0,0 +1,47 @@ +# Copyright 2019-2020 CERN and copyright holders of ALICE O2. +# See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +# All rights not expressly granted are reserved. +# +# This software is distributed under the terms of the GNU General Public +# License v3 (GPL Version 3), copied verbatim in the file "COPYING". +# +# In applying this license CERN does not waive the privileges and immunities +# granted to it by virtue of its status as an Intergovernmental Organization +# or submit itself to any jurisdiction. + +o2_add_library(MCHAlign + SOURCES + src/Aligner.cxx + src/AlignmentSpec.cxx + PUBLIC_LINK_LIBRARIES + O2::MathUtils + O2::CCDB + O2::DataFormatsMCH + O2::ForwardAlign + O2::MCHTracking + O2::MCHGeometryTransformer + O2::CommonUtils + O2::DataFormatsParameters + O2::DetectorsBase + O2::Framework + O2::DetectorsRaw + O2::Headers + O2::ReconstructionDataFormats + O2::DetectorsCommonDataFormats) + +o2_target_root_dictionary(MCHAlign + HEADERS + include/MCHAlign/Aligner.h + include/MCHAlign/AlignmentSpec.h) + +o2_add_executable( + alignment-workflow + SOURCES src/alignment-workflow.cxx + COMPONENT_NAME mch + PUBLIC_LINK_LIBRARIES + O2::ForwardAlign + O2::MCHAlign + Boost::program_options) + + + diff --git a/Detectors/MUON/MCH/Align/include/MCHAlign/Aligner.h b/Detectors/MUON/MCH/Align/include/MCHAlign/Aligner.h new file mode 100644 index 0000000000000..fe18c337cbe9a --- /dev/null +++ b/Detectors/MUON/MCH/Align/include/MCHAlign/Aligner.h @@ -0,0 +1,484 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/** @file CathodeSegmentation.h + * C++ Alignmnet . + * @author Javier Castillo Castellanos + */ + +#ifndef ALICEO2_MCH_ALIGNER +#define ALICEO2_MCH_ALIGNER + +#include +#include + +#include "DataFormatsMCH/Cluster.h" +#include "DetectorsCommonDataFormats/AlignParam.h" +#include "ForwardAlign/MillePede2.h" +#include "ForwardAlign/MillePedeRecord.h" +#include "MCHGeometryTransformer/Transformations.h" +#include "MCHTracking/Track.h" + +#include +#include +#include +#include +#include +#include + +namespace o2 +{ + +namespace mch +{ + +/// local track parameters, for refit +class LocalTrackParam +{ + public: + //* construction + LocalTrackParam() = default; + ~LocalTrackParam() = default; + + // private: + //* y and z + double fTrackX = 0.0; + double fTrackY = 0.0; + double fTrackZ = 0.0; + double fTrackSlopeX = 0.0; + double fTrackSlopeY = 0.0; +}; // class LocalTrackParam + +/// local track residual, for tempoarary eval +class LocalTrackClusterResidual +{ + public: + //* construction + LocalTrackClusterResidual() = default; + ~LocalTrackClusterResidual() = default; + + // private: + //* y and z + int fClDetElem = 0.0; + int fClDetElemNumber = 0.0; + + float fClusterX = 0.0; + float fClusterY = 0.0; + float fClusterZ = 0.0; + float fClusterXloc = 0.0; + float fClusterYloc = 0.0; + + float fTrackX = 0.0; + float fTrackY = 0.0; + float fTrackZ = 0.0; + float fTrackXloc = 0.0; + float fTrackYloc = 0.0; + + float fTrackSlopeX = 0.0; + float fTrackSlopeY = 0.0; + + float fBendingMomentum = 0.0; + + float fResiduXGlobal = 0.0; + float fResiduYGlobal = 0.0; + + float fResiduXLocal = 0.0; + float fResiduYLocal = 0.0; + + float fCharge = 0.0; + + float fBx = 0.0; + float fBy = 0.0; + float fBz = 0.0; + +}; // class LocalTrackClusterResidual + +class Aligner : public TObject +{ + + public: + Aligner(); + + ~Aligner() = default; + + // initialize + void init(TString DataRecFName = "recDataFile.root", TString ConsRecFName = "recConsFile.root"); + + // terminate + void terminate(void); + + // array dimendions + enum { + /// Number tracking stations + fgNSt = 5, + + /// Number tracking chambers + fgNCh = 10, + + /// Number of tracking modules + fgNTrkMod = 16, + + /// Number of half chambers + fgNHalfCh = 20, + + /// max number of detector elements per half chamber + fgNDetHalfChMax = 13, + + /// Total number of detection elements + /// (4*2 + 4*2 + 18*2 + 26*2 + 26*2) + fgNDetElem = 156, + + /// Number of local parameters + fNLocal = 4, // t_x, t_y, x0, y0 + + /// Number of degrees of freedom per chamber + fgNParCh = 4, // x,y,z,phi + + /// Number of global parameters + fNGlobal = fgNParCh * fgNDetElem + }; + + /// Number of detection elements per chamber + static const int fgNDetElemCh[fgNCh]; + + /// Sum of detection elements up to this chamber + static const int fgSNDetElemCh[fgNCh + 1]; + + /// Number of detection element per tracking module + static const int fgNDetElemHalfCh[fgNHalfCh]; + + /// list of detection elements per tracking module + static const int fgDetElemHalfCh[fgNHalfCh][fgNDetHalfChMax]; + + /// global parameter bit set, used for masks + enum ParameterMask { + ParX = 1 << 0, + ParY = 1 << 1, + ParZ = 1 << 2, + ParTZ = 1 << 3, + + ParAllTranslations = ParX | ParY | ParZ, + ParAllRotations = ParTZ, + ParAll = ParAllTranslations | ParAllRotations + + }; + + /// detector sides bit set, used for selecting sides in constrains + enum SidesMask { + SideTop = 1 << 0, + SideLeft = 1 << 1, + SideBottom = 1 << 2, + SideRight = 1 << 3, + SideTopLeft = SideTop | SideLeft, + SideTopRight = SideTop | SideRight, + SideBottomLeft = SideBottom | SideLeft, + SideBottomRight = SideBottom | SideRight, + AllSides = SideTop | SideBottom | SideLeft | SideRight + }; + + o2::fwdalign::MillePedeRecord* ProcessTrack(Track& track, const o2::mch::geo::TransformationCreator& transformation, Bool_t doAlignment, Double_t weight = 1); + + void ProcessTrack(o2::fwdalign::MillePedeRecord*); + + //@name modifiers + //@{ + + /// run number + void SetRunNumber(int id) + { + fRunNumber = id; + } + + /// Set flag for Magnetic field On/Off + void SetBFieldOn(bool value) + { + fBFieldOn = value; + } + + /// set to true to do refit evaluation + void SetDoEvaluation(bool value) + { + fDoEvaluation = value; + } + + /// set to true to refit tracks + void SetRefitStraightTracks(bool value) + { + fRefitStraightTracks = value; + } + + void SetAllowedVariation(int iPar, double value); + + void SetSigmaXY(double sigmaX, double sigmaY); + + void FixAll(unsigned int parameterMask = ParAll); + + void FixChamber(int iCh, unsigned int parameterMask = ParAll); + + void FixDetElem(int iDetElemId, unsigned int parameterMask = ParAll); + + void FixHalfSpectrometer(const bool* bChOnOff, unsigned int sidesMask = AllSides, unsigned int parameterMask = ParAll); + + void FixParameter(int iPar); + + void FixParameter(int iDetElem, int iPar) + { + FixParameter(iDetElem * fgNParCh + iPar); + } + + //@} + + //@name releasing detectors + //@{ + + void ReleaseChamber(int iCh, unsigned int parameterMask = ParAll); + + void ReleaseDetElem(int iDetElemId, unsigned int parameterMask = ParAll); + + void ReleaseParameter(int iPar); + + void ReleaseParameter(int iDetElem, int iPar) + { + ReleaseParameter(iDetElem * fgNParCh + iPar); + } + + //@} + + //@name grouping detectors + //@{ + + void GroupChamber(int iCh, unsigned int parameterMask = ParAll); + + void GroupHalfChamber(int iCh, int iHalf, unsigned int parameterMask = ParAll); + + void GroupDetElems(int detElemMin, int detElemMax, unsigned int parameterMask = ParAll); + + void GroupDetElems(const int* detElemList, int nDetElem, unsigned int parameterMask = ParAll); + + //@} + + //@name define non linearity + //@{ + + void SetChamberNonLinear(int iCh, unsigned int parameterMask); + + void SetDetElemNonLinear(int iSt, unsigned int parameterMask); + + void SetParameterNonLinear(int iPar); + + void SetParameterNonLinear(int iDetElem, int iPar) + { + SetParameterNonLinear(iDetElem * fgNParCh + iPar); + } + + //@} + + //@name constraints + //@{ + + void AddConstraints(const bool* bChOnOff, unsigned int parameterMask); + + void AddConstraints(const bool* bChOnOff, const bool* lVarXYT, unsigned int sidesMask = AllSides); + + //@} + + /// initialize global parameters to a give set of values + void InitGlobalParameters(double* par); + + /// perform global fit + void GlobalFit(std::vector& params, std::vector& errors, std::vector& pulls); + + /// print global parameters + void PrintGlobalParameters(void) const; + + /// get error on a given parameter + double GetParError(int iPar) const; + + void ReAlign(std::vector& params, std::vector& misAlignments); + + void SetAlignmentResolution(const TClonesArray* misAlignArray, int chId, double chResX, double chResY, double deResX, double deResY); + + TTree* GetResTree() + { + return fTTree; + } + + void SetReadOnly() + { + mRead = true; + } + + private: + /// Not implemented + Aligner(const Aligner& right); + + /// Not implemented + Aligner& operator=(const Aligner& right); + + /// Set array of local derivatives + void SetLocalDerivative(int index, double value) + { + fLocalDerivatives[index] = value; + } + + /// Set array of global derivatives + void SetGlobalDerivative(int index, double value) + { + fGlobalDerivatives[index] = value; + } + + /// refit track using straight track model + LocalTrackParam RefitStraightTrack(Track&, double) const; + + void FillDetElemData(const Cluster*); + + void FillRecPointData(const Cluster*); + + void FillTrackParamData(const TrackParam*); + + void LocalEquationX(const double* r); + + void LocalEquationY(const double* r); + + TGeoCombiTrans DeltaTransform(const double* detElemMisAlignment) const; + + bool isMatrixConvertedToAngles(const double* rot, double& psi, double& theta, double& phi) const; + + ///@name utilities + //@{ + + void AddConstraint(double* parameters, double value); + + int GetChamberId(int iDetElemNumber) const; + + bool DetElemIsValid(int iDetElemId) const; + + int GetDetElemNumber(int iDetElemId) const; + + TString GetParameterMaskString(unsigned int parameterMask) const; + + TString GetSidesMaskString(unsigned int sidesMask) const; + + //@} + + /// true when initialized + bool fInitialized; + + /// current run id + int fRunNumber; + + /// Flag for Magnetic filed On/Off + bool fBFieldOn; + + /// true if straight track refit is to be performed + bool fRefitStraightTracks; + + /// "Encouraged" variation for degrees of freedom + double fAllowVar[fgNParCh]; + + /// Initial value for chi2 cut + double fStartFac; + + /// Cut on residual for first iteration + double fResCutInitial; + + /// Cut on residual for other iterations + double fResCut; + + /// Detector independent alignment class + o2::fwdalign::MillePede2* fMillepede; // AliMillePede2 implementation + + /// MCH cluster class + o2::mch::Cluster* fCluster; + + /// Number of standard deviations for chi2 cut + int fNStdDev; + + /// Cluster (global) position + double fClustPos[3]; + + /// Track slope at reference point + double fTrackSlope0[2]; + + /// Track slope at current point + double fTrackSlope[2]; + + /// Track intersection at reference point + double fTrackPos0[3]; + + /// Track intersection at current point + double fTrackPos[3]; + + /// Current measurement (depend on B field On/Off) + double fMeas[2]; + + /// Estimated resolution on measurement + double fSigma[2]; + + /// degrees of freedom + enum { + kFixedParId = -1, + kFreeParId = kFixedParId - 1, + kGroupBaseId = -10 + }; + + /// Array of effective degrees of freedom + /// it is used to fix detectors, parameters, etc. + std::vector fGlobalParameterStatus; + + /// Array of global derivatives + std::vector fGlobalDerivatives; + + /// Array of local derivatives + std::vector fLocalDerivatives; + + /// current detection element number + int fDetElemNumber; + + /// running Track record + o2::fwdalign::MillePedeRecord fTrackRecord; + + /// Geometry transformation + o2::mch::geo::TransformationCreator fTransformCreator; + + /// preform evaluation + bool fDoEvaluation; + + /// original local track params + LocalTrackParam* fTrackParamOrig; + LocalTrackParam* fTrackParamNew; + + LocalTrackClusterResidual* fTrkClRes; + + /// output TFile + TFile* fTFile; + + /// output TTree + TTree* fTTree; + + /// Option switch for read/write mode + bool mRead; + + long mNEntriesAutoSave = 10000; ///< number of entries needed to call AutoSave for the output TTrees + + o2::fwdalign::MilleRecordWriter* mRecordWriter; ///< utility that handles the writing of the data records to a ROOT file + bool mWithConstraintsRecWriter; ///< boolean to be set to true if one wants to also write constaints records + o2::fwdalign::MilleRecordWriter* mConstraintsRecWriter; ///< utility that handles the writing of the constraints records + + o2::fwdalign::MilleRecordReader* mRecordReader; ///< utility that handles the reading of the data records from a ROOT file + bool mWithConstraintsRecReader = false; ///< boolean to be set to true if one wants to also read constaints records + o2::fwdalign::MilleRecordReader* mConstraintsRecReader; ///< utility that handles the reading of the constraints records + +}; // class Alignment + +} // namespace mch +} // namespace o2 +#endif // ALICEO2_MCH_ALIGNER_H_ diff --git a/Detectors/MUON/MCH/Align/include/MCHAlign/AlignmentSpec.h b/Detectors/MUON/MCH/Align/include/MCHAlign/AlignmentSpec.h new file mode 100644 index 0000000000000..283d07464bf09 --- /dev/null +++ b/Detectors/MUON/MCH/Align/include/MCHAlign/AlignmentSpec.h @@ -0,0 +1,32 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file AlignmentSpec.h +/// \brief Definition of alignment process for muon spectrometer +/// +/// \author Chi ZHANG, CEA-Saclay + +#ifndef O2_MCH_ALIGNMENTSPEC_H_ +#define O2_MCH_ALIGNMENTSPEC_H_ + +#include "Framework/DataProcessorSpec.h" + +namespace o2 +{ +namespace mch +{ + +o2::framework::DataProcessorSpec getAlignmentSpec(bool disableCCDB = false); + +} // end namespace mch +} // end namespace o2 + +#endif // O2_MCH_ALIGNMENTSPEC_H_ diff --git a/Detectors/MUON/MCH/Align/src/Aligner.cxx b/Detectors/MUON/MCH/Align/src/Aligner.cxx new file mode 100644 index 0000000000000..39a00cdd893d9 --- /dev/null +++ b/Detectors/MUON/MCH/Align/src/Aligner.cxx @@ -0,0 +1,1797 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +//----------------------------------------------------------------------------- +/// \file Aligner +/// Aligner class for the ALICE DiMuon spectrometer +/// +/// MUON specific alignment class which interface to Millepede. +/// For each track ProcessTrack calculates the local and global derivatives +/// at each cluster and fill the corresponding local equations. Provide methods +/// for fixing or constraining detection elements for best results. +/// +/// \author Javier Castillo Castellanos +//----------------------------------------------------------------------------- +#include +#include + +#include "DataFormatsMCH/Cluster.h" +#include "DetectorsCommonDataFormats/AlignParam.h" +#include "ForwardAlign/MillePede2.h" +#include "ForwardAlign/MillePedeRecord.h" +#include "Framework/Logger.h" +#include "MCHAlign/Aligner.h" +#include "MCHTracking/Track.h" +#include "MCHTracking/TrackParam.h" +#include "MCHGeometryTransformer/Transformations.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace o2 +{ +namespace mch +{ + +using namespace std; + +//_____________________________________________________________________ +// static variables +const int Aligner::fgNDetElemCh[Aligner::fgNCh] = {4, 4, 4, 4, 18, 18, 26, 26, 26, 26}; +const int Aligner::fgSNDetElemCh[Aligner::fgNCh + 1] = {0, 4, 8, 12, 16, 34, 52, 78, 104, 130, 156}; + +// number of detector elements in each half-chamber +const int Aligner::fgNDetElemHalfCh[Aligner::fgNHalfCh] = {2, 2, 2, 2, 2, 2, 2, 2, 9, 9, 9, 9, 13, 13, 13, 13, 13, 13, 13, 13}; + +// list of detector elements for each half chamber +const int Aligner::fgDetElemHalfCh[Aligner::fgNHalfCh][Aligner::fgNDetHalfChMax] = + { + {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}, + + {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}, + + {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}, + + {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}, + + {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}, + + {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}, + + {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}, + + {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}, + + {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} + +}; + +//_____________________________________________________________________ +/// self initialized array, used for adding constraints +class Array +{ + + public: + /// contructor + Array() + { + for (int i = 0; i < Aligner::fNGlobal; ++i) { + values[i] = 0; + } + } + + /// array + double values[Aligner::fNGlobal]; + + private: + /// Not implemented + Array(const Array&); + + /// Not implemented + Array& operator=(const Array&); +}; + +//________________________________________________________________________ +double Square(double x) { return x * x; } + +//_____________________________________________________________________ +Aligner::Aligner() + : TObject(), + fInitialized(false), + fRunNumber(0), + fBFieldOn(false), + fRefitStraightTracks(false), + fStartFac(65536), + fResCutInitial(1000), + fResCut(100), + fMillepede(nullptr), // to be modified + fCluster(nullptr), + fNStdDev(3), + fDetElemNumber(0), + fGlobalParameterStatus(std::vector(fNGlobal)), + fGlobalDerivatives(std::vector(fNGlobal)), + fLocalDerivatives(std::vector(fNLocal)), + fTrackRecord(), + mNEntriesAutoSave(10000), + mRecordWriter(new o2::fwdalign::MilleRecordWriter()), + mWithConstraintsRecWriter(false), + mConstraintsRecWriter(nullptr), + mRecordReader(new o2::fwdalign::MilleRecordReader()), + mWithConstraintsRecReader(false), + mConstraintsRecReader(nullptr), + fTransformCreator(), + // fGeoCombiTransInverse(), + fDoEvaluation(false), + mRead(false), + fTrackParamOrig(nullptr), + fTrackParamNew(nullptr), + fTrkClRes(nullptr), + fTFile(nullptr), + fTTree(nullptr) +{ + /// constructor + fSigma[0] = 1.5e-1; + fSigma[1] = 1.0e-2; + + // default allowed variations + fAllowVar[0] = 0.5; // x + fAllowVar[1] = 0.5; // y + fAllowVar[2] = 0.01; // phi_z + fAllowVar[3] = 5; // z + + // initialize millepede + fMillepede = new o2::fwdalign::MillePede2(); + + // initialize degrees of freedom + // by default all parameters are free + for (int iPar = 0; iPar < fNGlobal; ++iPar) { + fGlobalParameterStatus[iPar] = kFreeParId; + } + + // initialize local equations + for (int i = 0; i < fNLocal; ++i) { + fLocalDerivatives[i] = 0.0; + } + + for (int i = 0; i < fNGlobal; ++i) { + fGlobalDerivatives[i] = 0.0; + } +} + +//_____________________________________________________________________ +void Aligner::init(TString DataRecFName, TString ConsRecFName) +{ + + /// initialize + /** + initialize millipede + must be called after necessary detectors have been fixed, + but before constrains are added and before global parameters initial value are set + */ + if (fInitialized) { + LOG(fatal) << "Millepede already initialized"; + } + + if (!mRead) { + + mRecordWriter->setCyclicAutoSave(mNEntriesAutoSave); + mRecordWriter->setDataFileName(DataRecFName); + fMillepede->SetRecordWriter(mRecordWriter); + + if (mWithConstraintsRecWriter) { + mConstraintsRecWriter->setCyclicAutoSave(mNEntriesAutoSave); + mConstraintsRecWriter->setDataFileName(ConsRecFName); + fMillepede->SetConstraintsRecWriter(mConstraintsRecWriter); + } + + } else { + + TChain* ch = new TChain(mRecordReader->getDataTreeName()); + if (DataRecFName.EndsWith(".root")) { + ch->AddFile(DataRecFName); + } + int nent = ch->GetEntries(); + + if (nent < 1) { + LOG(fatal) << "Obtained chain is empty, please check your record ROOT file."; + } + + mRecordReader->connectToChain(ch); + fMillepede->SetRecordReader(mRecordReader); + + if (mConstraintsRecReader) { + + TChain* ch_cons = new TChain(mConstraintsRecReader->getDataTreeName()); + if (ConsRecFName.EndsWith(".root")) { + ch_cons->AddFile(ConsRecFName); + } + int nent_cons = ch_cons->GetEntries(); + + if (nent_cons < 1) { + LOG(fatal) << "Obtained chain is empty, please check your record ROOT file."; + } + + mConstraintsRecReader->connectToChain(ch_cons); + fMillepede->SetConstraintsRecReader(mConstraintsRecReader); + } + } + + // assign proper groupID to free parameters + int nGlobal = 0; + for (int iPar = 0; iPar < fNGlobal; ++iPar) { + + if (fGlobalParameterStatus[iPar] == kFixedParId) { + // fixed parameters are left unchanged + continue; + + } else if (fGlobalParameterStatus[iPar] == kFreeParId || fGlobalParameterStatus[iPar] == kGroupBaseId) { + + // free parameters or first element of group are assigned a new group id + fGlobalParameterStatus[iPar] = nGlobal++; + continue; + + } else if (fGlobalParameterStatus[iPar] < kGroupBaseId) { + + // get detector element id from status, get chamber parameter id + const int iDeBase(kGroupBaseId - 1 - fGlobalParameterStatus[iPar]); + const int iParBase = iPar % fgNParCh; + + // check + if (iDeBase < 0 || iDeBase >= iPar / fgNParCh) { + LOG(fatal) << "Group for parameter index " << iPar << " has wrong base detector element: " << iDeBase; + } + + // assign identical group id to current + fGlobalParameterStatus[iPar] = fGlobalParameterStatus[iDeBase * fgNParCh + iParBase]; + LOG(info) << "Parameter " << iPar << " grouped to detector " << iDeBase << " (" << GetParameterMaskString(1 << iParBase).Data() << ")"; + + } else { + LOG(fatal) << "Unrecognized parameter status for index " << iPar << ": " << fGlobalParameterStatus[iPar]; + } + } + + LOG(info) << "Free Parameters: " << nGlobal << " out of " << fNGlobal; + + // initialize millepedes + fMillepede->InitMille(fNGlobal, fNLocal, fNStdDev, fResCut, fResCutInitial, fGlobalParameterStatus); + + if (!mRead) { + mRecordWriter->init(); + } + + fInitialized = true; + + // some debug output + for (int iPar = 0; iPar < fgNParCh; ++iPar) { + LOG(info) << "fAllowVar[" << iPar << "]= " << fAllowVar[iPar]; + } + + // set allowed variations for all parameters + for (int iDet = 0; iDet < fgNDetElem; ++iDet) { + for (int iPar = 0; iPar < fgNParCh; ++iPar) { + fMillepede->SetParSigma(iDet * fgNParCh + iPar, fAllowVar[iPar]); + } + } + + // Set iterations + if (fStartFac > 1) { + fMillepede->SetIterations(fStartFac); + } + // setup monitoring TFile + if (fDoEvaluation) { + // if (fDoEvaluation && fRefitStraightTracks) { + string Path_file = Form("%s%s", "Residual", ".root"); + + fTFile = new TFile(Path_file.c_str(), "RECREATE"); + fTTree = new TTree("TreeE", "Evaluation"); + + const int kSplitlevel = 98; + const int kBufsize = 32000; + + // fTrackParamOrig = new LocalTrackParam(); + // fTTree->Branch("fTrackParamOrig", "LocalTrackParam", &fTrackParamOrig, kBufsize, kSplitlevel); + + // fTrackParamNew = new LocalTrackParam(); + // fTTree->Branch("fTrackParamNew", "LocalTrackParam", &fTrackParamNew, kBufsize, kSplitlevel); + + fTrkClRes = new o2::mch::LocalTrackClusterResidual(); + fTTree->Branch("fClDetElem", &(fTrkClRes->fClDetElem), "fClDetElem/I"); + fTTree->Branch("fClDetElemNumber", &(fTrkClRes->fClDetElemNumber), "fClDetElemNumber/I"); + fTTree->Branch("fClusterX", &(fTrkClRes->fClusterX), "fClusterX/F"); + fTTree->Branch("fClusterY", &(fTrkClRes->fClusterY), "fClusterY/F"); + fTTree->Branch("fTrackX", &(fTrkClRes->fTrackX), "fTrackX/F"); + fTTree->Branch("fTrackY", &(fTrkClRes->fTrackY), "fTrackY/F"); + fTTree->Branch("fClusterXloc", &(fTrkClRes->fClusterXloc), "fClusterXloc/F"); + fTTree->Branch("fClusterYloc", &(fTrkClRes->fClusterYloc), "fClusterYloc/F"); + fTTree->Branch("fTrackXloc", &(fTrkClRes->fTrackXloc), "fTrackXloc/F"); + fTTree->Branch("fTrackYloc", &(fTrkClRes->fTrackYloc), "fTrackYloc/F"); + fTTree->Branch("fTrackSlopeX", &(fTrkClRes->fTrackSlopeX), "fTrackSlopeX/F"); + fTTree->Branch("fTrackSlopeY", &(fTrkClRes->fTrackSlopeY), "fTrackSlopeY/F"); + fTTree->Branch("fBendingMomentum", &(fTrkClRes->fBendingMomentum), "fBendingMomentum/F"); + fTTree->Branch("fResiduXGlobal", &(fTrkClRes->fResiduXGlobal), "fResiduXGlobal/F"); + fTTree->Branch("fResiduYGlobal", &(fTrkClRes->fResiduYGlobal), "fResiduYGlobal/F"); + fTTree->Branch("fResiduXLocal", &(fTrkClRes->fResiduXLocal), "fResiduXLocal/F"); + fTTree->Branch("fResiduYLocal", &(fTrkClRes->fResiduYLocal), "fResiduYLocal/F"); + fTTree->Branch("fCharge", &(fTrkClRes->fCharge), "fCharge/F"); + fTTree->Branch("fClusterZ", &(fTrkClRes->fClusterZ), "fClusterZ/F"); + fTTree->Branch("fTrackZ", &(fTrkClRes->fTrackZ), "fTrackZ/F"); + fTTree->Branch("fBx", &(fTrkClRes->fBx), "fBx/F"); + fTTree->Branch("fBy", &(fTrkClRes->fBy), "fBy/F"); + fTTree->Branch("fBz", &(fTrkClRes->fBz), "fBz/F"); + } +} + +//_____________________________________________________ +void Aligner::terminate() +{ + mRecordWriter->terminate(); + fInitialized = kFALSE; + LOG(info) << "Closing Evaluation TFile"; + if (fTFile && fTTree) { + fTFile->cd(); + fTTree->Write(); + fTFile->Close(); + } +} + +//_____________________________________________________ +o2::fwdalign::MillePedeRecord* Aligner::ProcessTrack(Track& track, const o2::mch::geo::TransformationCreator& transformation, bool doAlignment, double weight) +{ + + /// process track for alignment minimization + + // reset track records + fTrackRecord.Reset(); + if (fMillepede->GetRecord()) { + fMillepede->GetRecord()->Reset(); + } + + // loop over clusters to get starting values + bool first(true); + + auto itTrackParam = track.begin(); + for (; itTrackParam != track.end(); ++itTrackParam) { + + // get cluster + const Cluster* cluster = itTrackParam->getClusterPtr(); + if (!cluster) { + continue; + } + // for first valid cluster, save track position as "starting" values + if (first) { + + first = false; + FillTrackParamData(&*itTrackParam); + fTrackPos0[0] = fTrackPos[0]; + fTrackPos0[1] = fTrackPos[1]; + fTrackPos0[2] = fTrackPos[2]; + fTrackSlope0[0] = fTrackSlope[0]; + fTrackSlope0[1] = fTrackSlope[1]; + + break; + } + } + + // redo straight track fit + if (fRefitStraightTracks) { + + // refit straight track + const LocalTrackParam trackParam(RefitStraightTrack(track, fTrackPos0[2])); + + /* + copy new parameters to stored ones for derivatives calculation + this is done only if BFieldOn is false, for which these parameters are used + */ + if (!fBFieldOn) { + fTrackPos0[0] = trackParam.fTrackX; + fTrackPos0[1] = trackParam.fTrackY; + fTrackPos0[2] = trackParam.fTrackZ; + fTrackSlope0[0] = trackParam.fTrackSlopeX; + fTrackSlope0[1] = trackParam.fTrackSlopeY; + } + } + + // second loop to perform alignment + itTrackParam = track.begin(); + for (; itTrackParam != track.end(); ++itTrackParam) { + + // get cluster + const Cluster* cluster = itTrackParam->getClusterPtr(); + if (!cluster) { + continue; + } + + // fill local variables for this position --> one measurement + + FillDetElemData(cluster); // function to get the transformation matrix + FillRecPointData(cluster); + FillTrackParamData(&*itTrackParam); + + // 'inverse' (GlobalToLocal) rotation matrix + // const double* r(fGeoCombiTransInverse.GetRotationMatrix()); + + auto trans = transformation(cluster->getDEId()); + // LOG(info) << Form("cluster ID: %i", cluster->getDEId()); + TMatrixD transMat(3, 4); + trans.GetTransformMatrix(transMat); + // transMat.Print(); + double r[12]; + r[0] = transMat(0, 0); + r[1] = transMat(0, 1); + r[2] = transMat(0, 2); + r[3] = transMat(1, 0); + r[4] = transMat(1, 1); + r[5] = transMat(1, 2); + r[6] = transMat(2, 0); + r[7] = transMat(2, 1); + r[8] = transMat(2, 2); + r[9] = transMat(0, 3); + r[10] = transMat(1, 3); + r[11] = transMat(2, 3); + + // calculate measurements + if (fBFieldOn) { + + // use residuals (cluster - track) for measurement + fMeas[0] = r[0] * (fClustPos[0] - fTrackPos[0]) + r[1] * (fClustPos[1] - fTrackPos[1]); + fMeas[1] = r[3] * (fClustPos[0] - fTrackPos[0]) + r[4] * (fClustPos[1] - fTrackPos[1]); + } else { + + // use cluster position for measurement + fMeas[0] = (r[0] * fClustPos[0] + r[1] * fClustPos[1]); + fMeas[1] = (r[3] * fClustPos[0] + r[4] * fClustPos[1]); + } + // printf("DE %d, X: %f %f ; Y: %f %f ; Z: %f\n", cluster->getDEId(), fClustPos[0], fTrackPos[0], fClustPos[1], fTrackPos[1], fClustPos[2]); + + if (fDoEvaluation) { + + const float InvBendingMom = itTrackParam->getInverseBendingMomentum(); + const float TrackCharge = itTrackParam->getCharge(); + + double B[3] = {0.0, 0.0, 0.0}; + double x[3] = {fTrackPos[0], fTrackPos[1], fTrackPos[2]}; + TGeoGlobalMagField::Instance()->Field(x, B); + const float Bx = B[0]; + const float By = B[1]; + const float Bz = B[2]; + + fTrkClRes->fClDetElem = cluster->getDEId(); + fTrkClRes->fClDetElemNumber = GetDetElemNumber(cluster->getDEId()); + fTrkClRes->fClusterX = fClustPos[0]; + fTrkClRes->fClusterY = fClustPos[1]; + fTrkClRes->fClusterZ = fClustPos[2]; + + // fTrkClRes->fTrackX = fTrackPos0[0] + fTrackSlope0[0] * (fTrackPos[2] - fTrackPos0[2]); // fTrackPos[0]; + // fTrkClRes->fTrackY = fTrackPos0[1] + fTrackSlope0[1] * (fTrackPos[2] - fTrackPos0[2]); // fTrackPos[1]; + // fTrkClRes->fTrackSlopeX = fTrackSlope0[0]; + // fTrkClRes->fTrackSlopeY = fTrackSlope0[1]; + + fTrkClRes->fTrackX = fTrackPos[0]; + fTrkClRes->fTrackY = fTrackPos[1]; + fTrkClRes->fTrackZ = fTrackPos[2]; + + fTrkClRes->fClusterXloc = r[0] * fClustPos[0] + r[1] * fClustPos[1]; + fTrkClRes->fClusterYloc = r[3] * fClustPos[0] + r[4] * fClustPos[1]; + + fTrkClRes->fTrackXloc = r[0] * fTrackPos[0] + r[1] * fTrackPos[1]; + fTrkClRes->fTrackYloc = r[3] * fTrackPos[0] + r[4] * fTrackPos[1]; + + fTrkClRes->fTrackSlopeX = fTrackSlope[0]; + fTrkClRes->fTrackSlopeY = fTrackSlope[1]; + + fTrkClRes->fBendingMomentum = TrackCharge / InvBendingMom; + + fTrkClRes->fResiduXGlobal = fClustPos[0] - fTrackPos[0]; + fTrkClRes->fResiduYGlobal = fClustPos[1] - fTrackPos[1]; + fTrkClRes->fResiduXLocal = r[0] * (fClustPos[0] - fTrackPos[0]) + r[1] * (fClustPos[1] - fTrackPos[1]); + fTrkClRes->fResiduYLocal = r[3] * (fClustPos[0] - fTrackPos[0]) + r[4] * (fClustPos[1] - fTrackPos[1]); + + fTrkClRes->fCharge = TrackCharge; + + fTrkClRes->fBx = Bx; + fTrkClRes->fBy = By; + fTrkClRes->fBz = Bz; + + if (fTTree) { + fTTree->Fill(); + } + } + // Set local equations + LocalEquationX(r); + LocalEquationY(r); + } + + // copy track record + mRecordWriter->setRecordRun(fRunNumber); + mRecordWriter->setRecordWeight(weight); + fTrackRecord = *fMillepede->GetRecord(); + + // save record data + if (doAlignment) { + mRecordWriter->fillRecordTree(); + } + + // return record + return &fTrackRecord; +} + +//______________________________________________________________________________ +void Aligner::ProcessTrack(o2::fwdalign::MillePedeRecord* trackRecord) +{ + LOG(fatal) << __PRETTY_FUNCTION__ << " is disabled"; + + /// process track record + if (!trackRecord) { + return; + } + + // // make sure record storage is initialized + if (!fMillepede->GetRecord()) { + mRecordWriter->init(); + } + // // copy content + *fMillepede->GetRecord() = *trackRecord; + + // save record + mRecordWriter->fillRecordTree(); + // write to local file + // mRecordWriter->terminate(); + + return; +} + +//_____________________________________________________________________ +void Aligner::FixAll(unsigned int mask) +{ + /// fix parameters matching mask, for all chambers + LOG(info) << "Fixing " << GetParameterMaskString(mask).Data() << " for all detector elements"; + + // fix all stations + for (int i = 0; i < fgNDetElem; ++i) { + if (mask & ParX) { + FixParameter(i, 0); + } + if (mask & ParY) { + FixParameter(i, 1); + } + if (mask & ParTZ) { + FixParameter(i, 2); + } + if (mask & ParZ) { + FixParameter(i, 3); + } + } +} + +//_____________________________________________________________________ +void Aligner::FixChamber(int iCh, unsigned int mask) +{ + /// fix parameters matching mask, for all detector elements in a given chamber, counting from 1 + + // check boundaries + if (iCh < 1 || iCh > 10) { + LOG(fatal) << "Invalid chamber index " << iCh; + } + + // get first and last element + const int iDetElemFirst = fgSNDetElemCh[iCh - 1]; + const int iDetElemLast = fgSNDetElemCh[iCh]; + for (int i = iDetElemFirst; i < iDetElemLast; ++i) { + + LOG(info) << "Fixing " << GetParameterMaskString(mask).Data() << " for detector element " << i; + + if (mask & ParX) { + FixParameter(i, 0); + } + if (mask & ParY) { + FixParameter(i, 1); + } + if (mask & ParTZ) { + FixParameter(i, 2); + } + if (mask & ParZ) { + FixParameter(i, 3); + } + } +} + +//_____________________________________________________________________ +void Aligner::FixDetElem(int iDetElemId, unsigned int mask) +{ + /// fix parameters matching mask, for a given detector element, counting from 0 + const int iDet(GetDetElemNumber(iDetElemId)); + if (mask & ParX) { + FixParameter(iDet, 0); + } + if (mask & ParY) { + FixParameter(iDet, 1); + } + if (mask & ParTZ) { + FixParameter(iDet, 2); + } + if (mask & ParZ) { + FixParameter(iDet, 3); + } +} + +//_____________________________________________________________________ +void Aligner::FixHalfSpectrometer(const bool* lChOnOff, unsigned int sidesMask, unsigned int mask) +{ + + /// Fix parameters matching mask for all detectors in selected chambers and selected sides of the spectrometer + for (int i = 0; i < fgNDetElem; ++i) { + + // get chamber matching detector + const int iCh(GetChamberId(i)); + if (!lChOnOff[iCh - 1]) { + continue; + } + + // get detector element in chamber + int lDetElemNumber = i - fgSNDetElemCh[iCh - 1]; + + // skip detector if its side is off + // stations 1 and 2 + if (iCh >= 1 && iCh <= 4) { + if (lDetElemNumber == 0 && !(sidesMask & SideTopRight)) { + continue; + } + if (lDetElemNumber == 1 && !(sidesMask & SideTopLeft)) { + continue; + } + if (lDetElemNumber == 2 && !(sidesMask & SideBottomLeft)) { + continue; + } + if (lDetElemNumber == 3 && !(sidesMask & SideBottomRight)) { + continue; + } + } + + // station 3 + if (iCh >= 5 && iCh <= 6) { + if (lDetElemNumber >= 0 && lDetElemNumber <= 4 && !(sidesMask & SideTopRight)) { + continue; + } + if (lDetElemNumber >= 5 && lDetElemNumber <= 10 && !(sidesMask & SideTopLeft)) { + continue; + } + if (lDetElemNumber >= 11 && lDetElemNumber <= 13 && !(sidesMask & SideBottomLeft)) { + continue; + } + if (lDetElemNumber >= 14 && lDetElemNumber <= 17 && !(sidesMask & SideBottomRight)) { + continue; + } + } + + // stations 4 and 5 + if (iCh >= 7 && iCh <= 10) { + if (lDetElemNumber >= 0 && lDetElemNumber <= 6 && !(sidesMask & SideTopRight)) { + continue; + } + if (lDetElemNumber >= 7 && lDetElemNumber <= 13 && !(sidesMask & SideTopLeft)) { + continue; + } + if (lDetElemNumber >= 14 && lDetElemNumber <= 19 && !(sidesMask & SideBottomLeft)) { + continue; + } + if (lDetElemNumber >= 20 && lDetElemNumber <= 25 && !(sidesMask & SideBottomRight)) { + continue; + } + } + + // detector is accepted, fix it + FixDetElem(i, mask); + } +} + +//______________________________________________________________________ +void Aligner::FixParameter(int iPar) +{ + + /// fix a given parameter, counting from 0 + if (fInitialized) { + LOG(fatal) << "Millepede already initialized"; + } + + fGlobalParameterStatus[iPar] = kFixedParId; +} + +//_____________________________________________________________________ +void Aligner::ReleaseChamber(int iCh, unsigned int mask) +{ + /// release parameters matching mask, for all detector elements in a given chamber, counting from 1 + + // check boundaries + if (iCh < 1 || iCh > 10) { + LOG(fatal) << "Invalid chamber index " << iCh; + } + + // get first and last element + const int iDetElemFirst = fgSNDetElemCh[iCh - 1]; + const int iDetElemLast = fgSNDetElemCh[iCh]; + for (int i = iDetElemFirst; i < iDetElemLast; ++i) { + + LOG(info) << "Releasing " << GetParameterMaskString(mask).Data() << " for detector element " << i; + + if (mask & ParX) { + ReleaseParameter(i, 0); + } + if (mask & ParY) { + ReleaseParameter(i, 1); + } + if (mask & ParTZ) { + ReleaseParameter(i, 2); + } + if (mask & ParZ) { + ReleaseParameter(i, 3); + } + } +} + +//_____________________________________________________________________ +void Aligner::ReleaseDetElem(int iDetElemId, unsigned int mask) +{ + /// release parameters matching mask, for a given detector element, counting from 0 + const int iDet(GetDetElemNumber(iDetElemId)); + if (mask & ParX) { + ReleaseParameter(iDet, 0); + } + if (mask & ParY) { + ReleaseParameter(iDet, 1); + } + if (mask & ParTZ) { + ReleaseParameter(iDet, 2); + } + if (mask & ParZ) { + ReleaseParameter(iDet, 3); + } +} + +//______________________________________________________________________ +void Aligner::ReleaseParameter(int iPar) +{ + + /// release a given parameter, counting from 0 + if (fInitialized) { + LOG(fatal) << "Millepede already initialized"; + } + + fGlobalParameterStatus[iPar] = kFreeParId; +} + +//_____________________________________________________________________ +void Aligner::GroupChamber(int iCh, unsigned int mask) +{ + /// group parameters matching mask for all detector elements in a given chamber, counting from 1 + if (iCh < 1 || iCh > fgNCh) { + LOG(fatal) << "Invalid chamber index " << iCh; + } + + const int detElemMin = 100 * iCh; + const int detElemMax = 100 * iCh + fgNDetElemCh[iCh] - 1; + GroupDetElems(detElemMin, detElemMax, mask); +} + +//_____________________________________________________________________ +void Aligner::GroupHalfChamber(int iCh, int iHalf, unsigned int mask) +{ + /// group parameters matching mask for all detector elements in a given tracking module (= half chamber), counting from 0 + if (iCh < 1 || iCh > fgNCh) { + LOG(fatal) << "Invalid chamber index " << iCh; + } + + if (iHalf < 0 || iHalf > 1) { + LOG(fatal) << "Invalid half chamber index " << iHalf; + } + + const int iHalfCh = 2 * (iCh - 1) + iHalf; + GroupDetElems(&fgDetElemHalfCh[iHalfCh][0], fgNDetElemHalfCh[iHalfCh], mask); +} + +//_____________________________________________________________________ +void Aligner::GroupDetElems(int detElemMin, int detElemMax, unsigned int mask) +{ + /// group parameters matching mask for all detector elements between min and max + // check number of detector elements + const int nDetElem = detElemMax - detElemMin + 1; + if (nDetElem < 2) { + LOG(fatal) << "Requested group of DEs " << detElemMin << "-" << detElemMax << " contains less than 2 DE's"; + } + + // create list + int* detElemList = new int[nDetElem]; + for (int i = 0; i < nDetElem; ++i) { + detElemList[i] = detElemMin + i; + } + + // group + GroupDetElems(detElemList, nDetElem, mask); + delete[] detElemList; +} + +//_____________________________________________________________________ +void Aligner::GroupDetElems(const int* detElemList, int nDetElem, unsigned int mask) +{ + /// group parameters matching mask for all detector elements in list + if (fInitialized) { + LOG(fatal) << "Millepede already initialized"; + } + + const int iDeBase(GetDetElemNumber(detElemList[0])); + for (int i = 0; i < nDetElem; ++i) { + const int iDeCurrent(GetDetElemNumber(detElemList[i])); + if (mask & ParX) { + fGlobalParameterStatus[iDeCurrent * fgNParCh + 0] = (i == 0) ? kGroupBaseId : (kGroupBaseId - iDeBase - 1); + } + if (mask & ParY) { + fGlobalParameterStatus[iDeCurrent * fgNParCh + 1] = (i == 0) ? kGroupBaseId : (kGroupBaseId - iDeBase - 1); + } + if (mask & ParTZ) { + fGlobalParameterStatus[iDeCurrent * fgNParCh + 2] = (i == 0) ? kGroupBaseId : (kGroupBaseId - iDeBase - 1); + } + if (mask & ParZ) { + fGlobalParameterStatus[iDeCurrent * fgNParCh + 3] = (i == 0) ? kGroupBaseId : (kGroupBaseId - iDeBase - 1); + } + + if (i == 0) { + LOG(info) << "Creating new group for detector " << detElemList[i] << " and variable " << GetParameterMaskString(mask).Data(); + } else { + LOG(info) << "Adding detector element " << detElemList[i] << " to current group"; + } + } +} + +//______________________________________________________________________ +void Aligner::SetChamberNonLinear(int iCh, unsigned int mask) +{ + /// Set parameters matching mask as non linear, for all detector elements in a given chamber, counting from 1 + const int iDetElemFirst = fgSNDetElemCh[iCh - 1]; + const int iDetElemLast = fgSNDetElemCh[iCh]; + for (int i = iDetElemFirst; i < iDetElemLast; ++i) { + + if (mask & ParX) { + SetParameterNonLinear(i, 0); + } + if (mask & ParY) { + SetParameterNonLinear(i, 1); + } + if (mask & ParTZ) { + SetParameterNonLinear(i, 2); + } + if (mask & ParZ) { + SetParameterNonLinear(i, 3); + } + } +} + +//_____________________________________________________________________ +void Aligner::SetDetElemNonLinear(int iDetElemId, unsigned int mask) +{ + /// Set parameters matching mask as non linear, for a given detector element, counting from 0 + const int iDet(GetDetElemNumber(iDetElemId)); + if (mask & ParX) { + SetParameterNonLinear(iDet, 0); + } + if (mask & ParY) { + SetParameterNonLinear(iDet, 1); + } + if (mask & ParTZ) { + SetParameterNonLinear(iDet, 2); + } + if (mask & ParZ) { + SetParameterNonLinear(iDet, 3); + } +} + +//______________________________________________________________________ +void Aligner::SetParameterNonLinear(int iPar) +{ + /// Set nonlinear flag for parameter iPar + if (!fInitialized) { + LOG(fatal) << "Millepede not initialized"; + } + + fMillepede->SetNonLinear(iPar); + LOG(info) << "Parameter " << iPar << " set to non linear "; +} + +//______________________________________________________________________ +void Aligner::AddConstraints(const bool* lChOnOff, unsigned int mask) +{ + /// Add constraint equations for selected chambers and degrees of freedom + + Array fConstraintX; + Array fConstraintY; + Array fConstraintTZ; + Array fConstraintZ; + + for (int i = 0; i < fgNDetElem; ++i) { + + // get chamber matching detector + const int iCh(GetChamberId(i)); + if (lChOnOff[iCh - 1]) { + + if (mask & ParX) { + fConstraintX.values[i * fgNParCh + 0] = 1.0; + } + if (mask & ParY) { + fConstraintY.values[i * fgNParCh + 1] = 1.0; + } + if (mask & ParTZ) { + fConstraintTZ.values[i * fgNParCh + 2] = 1.0; + } + if (mask & ParZ) { + fConstraintTZ.values[i * fgNParCh + 3] = 1.0; + } + } + } + + if (mask & ParX) { + AddConstraint(fConstraintX.values, 0.0); + } + if (mask & ParY) { + AddConstraint(fConstraintY.values, 0.0); + } + if (mask & ParTZ) { + AddConstraint(fConstraintTZ.values, 0.0); + } + if (mask & ParZ) { + AddConstraint(fConstraintZ.values, 0.0); + } +} + +//______________________________________________________________________ +void Aligner::AddConstraints(const bool* lChOnOff, const bool* lVarXYT, unsigned int sidesMask) +{ + /* + questions: + - is there not redundancy/inconsistency between lDetTLBR and lSpecLROnOff ? shouldn't we use only lDetTLBR ? + - why is weight ignored for ConstrainT and ConstrainB + - why is there no constrain on z + */ + + /// Add constraint equations for selected chambers, degrees of freedom and detector half + double lMeanY = 0.; + double lSigmaY = 0.; + double lMeanZ = 0.; + double lSigmaZ = 0.; + int lNDetElem = 0; + + for (int i = 0; i < fgNDetElem; ++i) { + + // get chamber matching detector + const int iCh(GetChamberId(i)); + + // skip detector if chamber is off + if (lChOnOff[iCh - 1]) { + continue; + } + + // get detector element id from detector element number + const int lDetElemNumber = i - fgSNDetElemCh[iCh - 1]; + const int lDetElemId = iCh * 100 + lDetElemNumber; + + // skip detector if its side is off + // stations 1 and 2 + if (iCh >= 1 && iCh <= 4) { + if (lDetElemNumber == 0 && !(sidesMask & SideTopRight)) { + continue; + } + if (lDetElemNumber == 1 && !(sidesMask & SideTopLeft)) { + continue; + } + if (lDetElemNumber == 2 && !(sidesMask & SideBottomLeft)) { + continue; + } + if (lDetElemNumber == 3 && !(sidesMask & SideBottomRight)) { + continue; + } + } + + // station 3 + if (iCh >= 5 && iCh <= 6) { + if (lDetElemNumber >= 0 && lDetElemNumber <= 4 && !(sidesMask & SideTopRight)) { + continue; + } + if (lDetElemNumber >= 5 && lDetElemNumber <= 10 && !(sidesMask & SideTopLeft)) { + continue; + } + if (lDetElemNumber >= 11 && lDetElemNumber <= 13 && !(sidesMask & SideBottomLeft)) { + continue; + } + if (lDetElemNumber >= 14 && lDetElemNumber <= 17 && !(sidesMask & SideBottomRight)) { + continue; + } + } + + // stations 4 and 5 + if (iCh >= 7 && iCh <= 10) { + if (lDetElemNumber >= 0 && lDetElemNumber <= 6 && !(sidesMask & SideTopRight)) { + continue; + } + if (lDetElemNumber >= 7 && lDetElemNumber <= 13 && !(sidesMask & SideTopLeft)) { + continue; + } + if (lDetElemNumber >= 14 && lDetElemNumber <= 19 && !(sidesMask & SideBottomLeft)) { + continue; + } + if (lDetElemNumber >= 20 && lDetElemNumber <= 25 && !(sidesMask & SideBottomRight)) { + continue; + } + } + + // get global x, y and z position + double lDetElemGloX = 0.; + double lDetElemGloY = 0.; + double lDetElemGloZ = 0.; + + auto fTransform = fTransformCreator(lDetElemId); + o2::math_utils::Point3D SlatPos{0.0, 0.0, 0.0}; + o2::math_utils::Point3D GlobalPos; + + fTransform.LocalToMaster(SlatPos, GlobalPos); + lDetElemGloX = GlobalPos.x(); + lDetElemGloY = GlobalPos.y(); + lDetElemGloZ = GlobalPos.z(); + // fTransform->Local2Global(lDetElemId, 0, 0, 0, lDetElemGloX, lDetElemGloY, lDetElemGloZ); + + // increment mean Y, mean Z, sigmas and number of accepted detectors + lMeanY += lDetElemGloY; + lSigmaY += lDetElemGloY * lDetElemGloY; + lMeanZ += lDetElemGloZ; + lSigmaZ += lDetElemGloZ * lDetElemGloZ; + lNDetElem++; + } + + // calculate mean values + lMeanY /= lNDetElem; + lSigmaY /= lNDetElem; + lSigmaY = TMath::Sqrt(lSigmaY - lMeanY * lMeanY); + lMeanZ /= lNDetElem; + lSigmaZ /= lNDetElem; + lSigmaZ = TMath::Sqrt(lSigmaZ - lMeanZ * lMeanZ); + LOG(info) << "Used " << lNDetElem << " DetElem, MeanZ= " << lMeanZ << ", SigmaZ= " << lSigmaZ; + + // create all possible arrays + Array fConstraintX[4]; // Array for constraint equation X + Array fConstraintY[4]; // Array for constraint equation Y + Array fConstraintP[4]; // Array for constraint equation P + Array fConstraintXZ[4]; // Array for constraint equation X vs Z + Array fConstraintYZ[4]; // Array for constraint equation Y vs Z + Array fConstraintPZ[4]; // Array for constraint equation P vs Z + + // do we really need these ? + Array fConstraintXY[4]; // Array for constraint equation X vs Y + Array fConstraintYY[4]; // Array for constraint equation Y vs Y + Array fConstraintPY[4]; // Array for constraint equation P vs Y + + // fill bool sides array based on masks, for convenience + bool lDetTLBR[4]; + lDetTLBR[0] = sidesMask & SideTop; + lDetTLBR[1] = sidesMask & SideLeft; + lDetTLBR[2] = sidesMask & SideBottom; + lDetTLBR[3] = sidesMask & SideRight; + + for (int i = 0; i < fgNDetElem; ++i) { + + // get chamber matching detector + const int iCh(GetChamberId(i)); + + // skip detector if chamber is off + if (!lChOnOff[iCh - 1]) { + continue; + } + + // get detector element id from detector element number + const int lDetElemNumber = i - fgSNDetElemCh[iCh - 1]; + const int lDetElemId = iCh * 100 + lDetElemNumber; + + // get global x, y and z position + double lDetElemGloX = 0.; + double lDetElemGloY = 0.; + double lDetElemGloZ = 0.; + + auto fTransform = fTransformCreator(lDetElemId); + o2::math_utils::Point3D SlatPos{0.0, 0.0, 0.0}; + o2::math_utils::Point3D GlobalPos; + + fTransform.LocalToMaster(SlatPos, GlobalPos); + lDetElemGloX = GlobalPos.x(); + lDetElemGloY = GlobalPos.y(); + lDetElemGloZ = GlobalPos.z(); + // fTransform->Local2Global(lDetElemId, 0, 0, 0, lDetElemGloX, lDetElemGloY, lDetElemGloZ); + + // loop over sides + for (int iSide = 0; iSide < 4; iSide++) { + + // skip if side is not selected + if (!lDetTLBR[iSide]) { + continue; + } + + // skip detector if it is not in the selected side + // stations 1 and 2 + if (iCh >= 1 && iCh <= 4) { + if (lDetElemNumber == 0 && !(iSide == 0 || iSide == 3)) { + continue; // top-right + } + if (lDetElemNumber == 1 && !(iSide == 0 || iSide == 1)) { + continue; // top-left + } + if (lDetElemNumber == 2 && !(iSide == 2 || iSide == 1)) { + continue; // bottom-left + } + if (lDetElemNumber == 3 && !(iSide == 2 || iSide == 3)) { + continue; // bottom-right + } + } + + // station 3 + if (iCh >= 5 && iCh <= 6) { + if (lDetElemNumber >= 0 && lDetElemNumber <= 4 && !(iSide == 0 || iSide == 3)) { + continue; // top-right + } + if (lDetElemNumber >= 5 && lDetElemNumber <= 9 && !(iSide == 0 || iSide == 1)) { + continue; // top-left + } + if (lDetElemNumber >= 10 && lDetElemNumber <= 13 && !(iSide == 2 || iSide == 1)) { + continue; // bottom-left + } + if (lDetElemNumber >= 14 && lDetElemNumber <= 17 && !(iSide == 2 || iSide == 3)) { + continue; // bottom-right + } + } + + // stations 4 and 5 + if (iCh >= 7 && iCh <= 10) { + if (lDetElemNumber >= 0 && lDetElemNumber <= 6 && !(iSide == 0 || iSide == 3)) { + continue; // top-right + } + if (lDetElemNumber >= 7 && lDetElemNumber <= 13 && !(iSide == 0 || iSide == 1)) { + continue; // top-left + } + if (lDetElemNumber >= 14 && lDetElemNumber <= 19 && !(iSide == 2 || iSide == 1)) { + continue; // bottom-left + } + if (lDetElemNumber >= 20 && lDetElemNumber <= 25 && !(iSide == 2 || iSide == 3)) { + continue; // bottom-right + } + } + + // constrain x + if (lVarXYT[0]) { + fConstraintX[iSide].values[i * fgNParCh + 0] = 1; + } + + // constrain y + if (lVarXYT[1]) { + fConstraintY[iSide].values[i * fgNParCh + 1] = 1; + } + + // constrain phi (rotation around z) + if (lVarXYT[2]) { + fConstraintP[iSide].values[i * fgNParCh + 2] = 1; + } + + // x-z shearing + if (lVarXYT[3]) { + fConstraintXZ[iSide].values[i * fgNParCh + 0] = (lDetElemGloZ - lMeanZ) / lSigmaZ; + } + + // y-z shearing + if (lVarXYT[4]) { + fConstraintYZ[iSide].values[i * fgNParCh + 1] = (lDetElemGloZ - lMeanZ) / lSigmaZ; + } + + // phi-z shearing + if (lVarXYT[5]) { + fConstraintPZ[iSide].values[i * fgNParCh + 2] = (lDetElemGloZ - lMeanZ) / lSigmaZ; + } + + // x-y shearing + if (lVarXYT[6]) { + fConstraintXY[iSide].values[i * fgNParCh + 0] = (lDetElemGloY - lMeanY) / lSigmaY; + } + + // y-y shearing + if (lVarXYT[7]) { + fConstraintYY[iSide].values[i * fgNParCh + 1] = (lDetElemGloY - lMeanY) / lSigmaY; + } + + // phi-y shearing + if (lVarXYT[8]) { + fConstraintPY[iSide].values[i * fgNParCh + 2] = (lDetElemGloY - lMeanY) / lSigmaY; + } + } + } + + // pass constraints to millepede + for (int iSide = 0; iSide < 4; iSide++) { + // skip if side is not selected + if (!lDetTLBR[iSide]) { + continue; + } + + if (lVarXYT[0]) { + AddConstraint(fConstraintX[iSide].values, 0.0); + } + if (lVarXYT[1]) { + AddConstraint(fConstraintY[iSide].values, 0.0); + } + if (lVarXYT[2]) { + AddConstraint(fConstraintP[iSide].values, 0.0); + } + if (lVarXYT[3]) { + AddConstraint(fConstraintXZ[iSide].values, 0.0); + } + if (lVarXYT[4]) { + AddConstraint(fConstraintYZ[iSide].values, 0.0); + } + if (lVarXYT[5]) { + AddConstraint(fConstraintPZ[iSide].values, 0.0); + } + if (lVarXYT[6]) { + AddConstraint(fConstraintXY[iSide].values, 0.0); + } + if (lVarXYT[7]) { + AddConstraint(fConstraintYY[iSide].values, 0.0); + } + if (lVarXYT[8]) { + AddConstraint(fConstraintPY[iSide].values, 0.0); + } + } +} + +//______________________________________________________________________ +void Aligner::InitGlobalParameters(double* par) +{ + /// Initialize global parameters with par array + if (!fInitialized) { + LOG(fatal) << "Millepede is not initialized"; + } + + fMillepede->SetGlobalParameters(par); +} + +//______________________________________________________________________ +void Aligner::SetAllowedVariation(int iPar, double value) +{ + /// "Encouraged" variation for degrees of freedom + // check initialization + if (fInitialized) { + LOG(fatal) << "Millepede already initialized"; + } + + // check initialization + if (!(iPar >= 0 && iPar < fgNParCh)) { + LOG(fatal) << "Invalid index: " << iPar; + } + + fAllowVar[iPar] = value; +} + +//______________________________________________________________________ +void Aligner::SetSigmaXY(double sigmaX, double sigmaY) +{ + + /// Set expected measurement resolution + fSigma[0] = sigmaX; + fSigma[1] = sigmaY; + + // print + for (int i = 0; i < 2; ++i) { + LOG(info) << "fSigma[" << i << "] =" << fSigma[i]; + } +} + +//_____________________________________________________ +void Aligner::GlobalFit(std::vector& parameters, std::vector& errors, std::vector& pulls) +{ + /// Call global fit; Global parameters are stored in parameters + fMillepede->GlobalFit(parameters, errors, pulls); + + LOG(info) << "Done fitting global parameters"; + for (int iDet = 0; iDet < fgNDetElem; ++iDet) { + LOG(info) << iDet << " " << parameters[iDet * fgNParCh + 0] << " " << parameters[iDet * fgNParCh + 1] << " " << parameters[iDet * fgNParCh + 3] << " " << parameters[iDet * fgNParCh + 2]; + } +} + +//_____________________________________________________ +void Aligner::PrintGlobalParameters() const +{ + fMillepede->PrintGlobalParameters(); +} + +//_____________________________________________________ +double Aligner::GetParError(int iPar) const +{ + return fMillepede->GetParError(iPar); +} + +//______________________________________________________________________ +void Aligner::ReAlign( + std::vector& params, + std::vector& misAlignments) +{ + + /// Returns a new AliMUONGeometryTransformer with the found misalignments + /// applied. + + // Takes the internal geometry module transformers, copies them + // and gets the Detection Elements from them. + // Takes misalignment parameters and applies these + // to the local transform of the Detection Element + // Obtains the global transform by multiplying the module transformer + // transformation with the local transformation + // Applies the global transform to a new detection element + // Adds the new detection element to a new module transformer + // Adds the new module transformer to a new geometry transformer + // Returns the new geometry transformer + + double lModuleMisAlignment[fgNParCh] = {0}; + double lDetElemMisAlignment[fgNParCh] = {0}; + + o2::detectors::AlignParam lAP; + for (int hc = 0; hc < 20; hc++) { + + TGeoCombiTrans localDeltaTransform; + localDeltaTransform = DeltaTransform(lModuleMisAlignment); + + std::string sname = fmt::format("MCH/HC{}", hc); + lAP.setSymName(sname.c_str()); + + double lPsi, lTheta, lPhi = 0.; + if (!isMatrixConvertedToAngles(localDeltaTransform.GetRotationMatrix(), + lPsi, lTheta, lPhi)) { + LOG(error) << "Problem extracting angles!"; + } + + lAP.setGlobalParams(localDeltaTransform); + + // lAP.print(); + lAP.applyToGeometry(); + params.emplace_back(lAP); + for (int de = 0; de < fgNDetElemHalfCh[hc]; de++) { + + // store detector element id and number + const int iDetElemId = fgDetElemHalfCh[hc][de]; + if (DetElemIsValid(iDetElemId)) { + + const int iDetElemNumber(GetDetElemNumber(iDetElemId)); + + for (int i = 0; i < fgNParCh; ++i) { + lDetElemMisAlignment[i] = 0.0; + if (hc < fgNHalfCh) { + lDetElemMisAlignment[i] = misAlignments[iDetElemNumber * fgNParCh + i]; + } + } + + sname = fmt::format("MCH/HC{}/DE{}", hc, fgDetElemHalfCh[hc][de]); + lAP.setSymName(sname.c_str()); + localDeltaTransform = DeltaTransform(lDetElemMisAlignment); + + if (!isMatrixConvertedToAngles(localDeltaTransform.GetRotationMatrix(), + lPsi, lTheta, lPhi)) { + LOG(error) << "Problem extracting angles for " << sname.c_str(); + } + + lAP.setGlobalParams(localDeltaTransform); + lAP.applyToGeometry(); + params.emplace_back(lAP); + + } else { + + // "invalid" detector elements come from MTR and are left unchanged + LOG(info) << fmt::format("Keeping detElement {} unchanged\n", iDetElemId); + } + } + } + + // return params; +} + +//______________________________________________________________________ +void Aligner::SetAlignmentResolution(const TClonesArray* misAlignArray, int rChId, double chResX, double chResY, double deResX, double deResY) +{ + + /// Set alignment resolution to misalign objects to be stored in CDB + /// if rChId is > 0 set parameters for this chamber only, counting from 1 + TMatrixDSym mChCorrMatrix(6); + mChCorrMatrix[0][0] = chResX * chResX; + mChCorrMatrix[1][1] = chResY * chResY; + + TMatrixDSym mDECorrMatrix(6); + mDECorrMatrix[0][0] = deResX * deResX; + mDECorrMatrix[1][1] = deResY * deResY; + + o2::detectors::AlignParam* alignMat = nullptr; + + for (int chId = 0; chId <= 9; ++chId) { + + // skip chamber if selection is valid, and does not match + if (rChId > 0 && chId + 1 != rChId) { + continue; + } + + TString chName1; + TString chName2; + if (chId < 4) { + + chName1 = Form("GM%d", chId); + chName2 = Form("GM%d", chId); + + } else { + + chName1 = Form("GM%d", 4 + (chId - 4) * 2); + chName2 = Form("GM%d", 4 + (chId - 4) * 2 + 1); + } + + for (int i = 0; i < misAlignArray->GetEntries(); ++i) { + + alignMat = (o2::detectors::AlignParam*)misAlignArray->At(i); + TString volName(alignMat->getSymName()); + if ((volName.Contains(chName1) && + ((volName.Last('/') == volName.Index(chName1) + chName1.Length()) || + (volName.Length() == volName.Index(chName1) + chName1.Length()))) || + (volName.Contains(chName2) && + ((volName.Last('/') == volName.Index(chName2) + chName2.Length()) || + (volName.Length() == volName.Index(chName2) + chName2.Length())))) { + + volName.Remove(0, volName.Last('/') + 1); + } + } + } +} + +//_____________________________________________________ +LocalTrackParam Aligner::RefitStraightTrack(Track& track, double z0) const +{ + + // initialize matrices + TMatrixD AtGASum(4, 4); + AtGASum.Zero(); + + TMatrixD AtGMSum(4, 1); + AtGMSum.Zero(); + + // loop over clusters + for (auto itTrackParam(track.begin()); itTrackParam != track.end(); ++itTrackParam) { + + // get cluster + const Cluster* cluster = itTrackParam->getClusterPtr(); + if (!cluster) { + continue; + } + + // projection matrix + TMatrixD A(2, 4); + A.Zero(); + A(0, 0) = 1; + A(0, 2) = (cluster->getZ() - z0); + A(1, 1) = 1; + A(1, 3) = (cluster->getZ() - z0); + + TMatrixD At(TMatrixD::kTransposed, A); + + // gain matrix + TMatrixD G(2, 2); + G.Zero(); + G(0, 0) = 1.0 / Square(cluster->getEx()); + G(1, 1) = 1.0 / Square(cluster->getEy()); + + const TMatrixD AtG(At, TMatrixD::kMult, G); + const TMatrixD AtGA(AtG, TMatrixD::kMult, A); + AtGASum += AtGA; + + // measurement + TMatrixD M(2, 1); + M(0, 0) = cluster->getX(); + M(1, 0) = cluster->getY(); + const TMatrixD AtGM(AtG, TMatrixD::kMult, M); + AtGMSum += AtGM; + } + + // perform inversion + TMatrixD AtGASumInv(TMatrixD::kInverted, AtGASum); + TMatrixD X(AtGASumInv, TMatrixD::kMult, AtGMSum); + + // fill output parameters + LocalTrackParam out; + out.fTrackX = X(0, 0); + out.fTrackY = X(1, 0); + out.fTrackZ = z0; + out.fTrackSlopeX = X(2, 0); + out.fTrackSlopeY = X(3, 0); + + return out; +} + +//_____________________________________________________ +void Aligner::FillDetElemData(const Cluster* cluster) +{ + + /// Get information of current detection element + // get detector element number from Alice ID + const int detElemId = cluster->getDEId(); + fDetElemNumber = GetDetElemNumber(detElemId); +} + +//______________________________________________________________________ +void Aligner::FillRecPointData(const Cluster* cluster) +{ + + /// Get information of current cluster + fClustPos[0] = cluster->getX(); + fClustPos[1] = cluster->getY(); + fClustPos[2] = cluster->getZ(); +} + +//______________________________________________________________________ +void Aligner::FillTrackParamData(const TrackParam* trackParam) +{ + + /// Get information of current track at current cluster + fTrackPos[0] = trackParam->getNonBendingCoor(); + fTrackPos[1] = trackParam->getBendingCoor(); + fTrackPos[2] = trackParam->getZ(); + fTrackSlope[0] = trackParam->getNonBendingSlope(); + fTrackSlope[1] = trackParam->getBendingSlope(); +} + +//______________________________________________________________________ +void Aligner::LocalEquationX(const double* r) +{ + /// local equation along X + + // 'inverse' (GlobalToLocal) rotation matrix + // const double* r(fGeoCombiTransInverse.GetRotationMatrix()); + + // local derivatives + SetLocalDerivative(0, r[0]); + SetLocalDerivative(1, r[0] * (fTrackPos[2] - fTrackPos0[2])); + // SetLocalDerivative(1, -r[0] * fTrackPos[2]); + + SetLocalDerivative(2, r[1]); + SetLocalDerivative(3, r[1] * (fTrackPos[2] - fTrackPos0[2])); + // SetLocalDerivative(3, -r[1] * fTrackPos[2]); + + // global derivatives + /* + alignment parameters are + 0: delta_x + 1: delta_y + 2: delta_phiz + 3: delta_z + */ + + SetGlobalDerivative(fDetElemNumber * fgNParCh + 0, -r[0]); + SetGlobalDerivative(fDetElemNumber * fgNParCh + 1, -r[1]); + + if (fBFieldOn) { + + // use local position for derivatives vs 'delta_phi_z' + SetGlobalDerivative(fDetElemNumber * fgNParCh + 2, -r[1] * fTrackPos[0] + r[0] * fTrackPos[1]); + + // use local slopes for derivatives vs 'delta_z' + SetGlobalDerivative(fDetElemNumber * fgNParCh + 3, r[0] * fTrackSlope[0] + r[1] * fTrackSlope[1]); + + } else { + + // local copy of extrapolated track positions + const double trackPosX = fTrackPos0[0] + fTrackSlope0[0] * (fTrackPos[2] - fTrackPos0[2]); + const double trackPosY = fTrackPos0[1] + fTrackSlope0[1] * (fTrackPos[2] - fTrackPos0[2]); + + // use properly extrapolated position for derivatives vs 'delta_phi_z' + SetGlobalDerivative(fDetElemNumber * fgNParCh + 2, -r[1] * trackPosX + r[0] * trackPosY); + // SetGlobalDerivative(fDetElemNumber * fgNParCh + 2, -r[1] * (trackPosX - r[9]) + r[0] * (trackPosY - r[10])); + + // use slopes at origin for derivatives vs 'delta_z' + SetGlobalDerivative(fDetElemNumber * fgNParCh + 3, r[0] * fTrackSlope0[0] + r[1] * fTrackSlope0[1]); + } + + // store local equation + fMillepede->SetLocalEquation(fGlobalDerivatives, fLocalDerivatives, fMeas[0], fSigma[0]); +} + +//______________________________________________________________________ +void Aligner::LocalEquationY(const double* r) +{ + /// local equation along Y + + // 'inverse' (GlobalToLocal) rotation matrix + // const double* r(fGeoCombiTransInverse.GetRotationMatrix()); + + // store local derivatives + SetLocalDerivative(0, r[3]); + SetLocalDerivative(1, r[3] * (fTrackPos[2] - fTrackPos0[2])); + // SetLocalDerivative(1, -r[3] * fTrackPos[2]); + + SetLocalDerivative(2, r[4]); + SetLocalDerivative(3, r[4] * (fTrackPos[2] - fTrackPos0[2])); + // SetLocalDerivative(3, -r[4] * fTrackPos[2]); + + // set global derivatives + SetGlobalDerivative(fDetElemNumber * fgNParCh + 0, -r[3]); + SetGlobalDerivative(fDetElemNumber * fgNParCh + 1, -r[4]); + + if (fBFieldOn) { + + // use local position for derivatives vs 'delta_phi' + SetGlobalDerivative(fDetElemNumber * fgNParCh + 2, -r[4] * fTrackPos[0] + r[3] * fTrackPos[1]); + + // use local slopes for derivatives vs 'delta_z' + SetGlobalDerivative(fDetElemNumber * fgNParCh + 3, r[3] * fTrackSlope[0] + r[4] * fTrackSlope[1]); + + } else { + + // local copy of extrapolated track positions + const double trackPosX = fTrackPos0[0] + fTrackSlope0[0] * (fTrackPos[2] - fTrackPos0[2]); + const double trackPosY = fTrackPos0[1] + fTrackSlope0[1] * (fTrackPos[2] - fTrackPos0[2]); + + // use properly extrapolated position for derivatives vs 'delta_phi' + SetGlobalDerivative(fDetElemNumber * fgNParCh + 2, -r[4] * trackPosX + r[3] * trackPosY); + + // use slopes at origin for derivatives vs 'delta_z' + SetGlobalDerivative(fDetElemNumber * fgNParCh + 3, r[3] * fTrackSlope0[0] + r[4] * fTrackSlope0[1]); + } + + // store local equation + fMillepede->SetLocalEquation(fGlobalDerivatives, fLocalDerivatives, fMeas[1], fSigma[1]); +} + +//_________________________________________________________________________ +TGeoCombiTrans Aligner::DeltaTransform(const double* lMisAlignment) const +{ + /// Get Delta Transformation, based on alignment parameters + + // translation + const TGeoTranslation deltaTrans(lMisAlignment[0], lMisAlignment[1], lMisAlignment[3]); + + // rotation + TGeoRotation deltaRot; + deltaRot.RotateZ(lMisAlignment[2] * 180. / TMath::Pi()); + + // combined rotation and translation. + return TGeoCombiTrans(deltaTrans, deltaRot); +} + +bool Aligner::isMatrixConvertedToAngles(const double* rot, double& psi, double& theta, double& phi) const +{ + /// Calculates the Euler angles in "x y z" notation + /// using the rotation matrix + /// Returns false in case the rotation angles can not be + /// extracted from the matrix + // + if (std::abs(rot[0]) < 1e-7 || std::abs(rot[8]) < 1e-7) { + LOG(error) << "Failed to extract roll-pitch-yall angles!"; + return false; + } + psi = std::atan2(-rot[5], rot[8]); + theta = std::asin(rot[2]); + phi = std::atan2(-rot[1], rot[0]); + return true; +} + +//______________________________________________________________________ +void Aligner::AddConstraint(double* par, double value) +{ + /// Constrain equation defined by par to value + if (!fInitialized) { + LOG(fatal) << "Millepede is not initialized"; + } + + std::vector vpar(fNGlobal); + for (int i = 0; i < fNGlobal; i++) { + vpar[i] = par[i]; + } + + fMillepede->SetGlobalConstraint(vpar, value); +} + +//______________________________________________________________________ +bool Aligner::DetElemIsValid(int iDetElemId) const +{ + /// return true if given detector element is valid (and belongs to muon tracker) + const int iCh = iDetElemId / 100; + const int iDet = iDetElemId % 100; + return (iCh > 0 && iCh <= fgNCh && iDet < fgNDetElemCh[iCh - 1]); +} + +//______________________________________________________________________ +int Aligner::GetDetElemNumber(int iDetElemId) const +{ + /// get det element number from ID + // get chamber and element number in chamber + const int iCh = iDetElemId / 100; + const int iDet = iDetElemId % 100; + + // make sure detector index is valid + if (!(iCh > 0 && iCh <= fgNCh && iDet < fgNDetElemCh[iCh - 1])) { + LOG(fatal) << "Invalid detector element id: " << iDetElemId; + } + + // add number of detectors up to this chamber + return iDet + fgSNDetElemCh[iCh - 1]; +} + +//______________________________________________________________________ +int Aligner::GetChamberId(int iDetElemNumber) const +{ + /// get chamber (counting from 1) matching a given detector element id + int iCh(0); + for (iCh = 0; iCh < fgNCh; iCh++) { + if (iDetElemNumber < fgSNDetElemCh[iCh]) { + break; + } + } + + return iCh; +} + +//______________________________________________________________________ +TString Aligner::GetParameterMaskString(unsigned int mask) const +{ + TString out; + if (mask & ParX) { + out += "X"; + } + if (mask & ParY) { + out += "Y"; + } + if (mask & ParZ) { + out += "Z"; + } + if (mask & ParTZ) { + out += "T"; + } + return out; +} + +//______________________________________________________________________ +TString Aligner::GetSidesMaskString(unsigned int mask) const +{ + TString out; + if (mask & SideTop) { + out += "T"; + } + if (mask & SideLeft) { + out += "L"; + } + if (mask & SideBottom) { + out += "B"; + } + if (mask & SideRight) { + out += "R"; + } + return out; +} + +} // namespace mch +} // namespace o2 \ No newline at end of file diff --git a/Detectors/MUON/MCH/Align/src/AlignmentSpec.cxx b/Detectors/MUON/MCH/Align/src/AlignmentSpec.cxx new file mode 100644 index 0000000000000..2d50a8235bc0c --- /dev/null +++ b/Detectors/MUON/MCH/Align/src/AlignmentSpec.cxx @@ -0,0 +1,946 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file AlignmentSpec.cxx +/// \brief Implementation of alignment process for muon spectrometer +/// +/// \author Chi ZHANG, CEA-Saclay + +#include "MCHAlign/AlignmentSpec.h" + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "Framework/CallbackService.h" +#include "Framework/ConcreteDataMatcher.h" +#include "Framework/ConfigParamRegistry.h" +#include "Framework/ControlService.h" +#include "Framework/DataProcessorSpec.h" +#include "Framework/CCDBParamSpec.h" +#include "Framework/Lifetime.h" +#include "Framework/Output.h" +#include "Framework/Task.h" +#include "Framework/Logger.h" + +#include "CCDB/BasicCCDBManager.h" +#include "CCDB/CCDBTimeStampUtils.h" +#include "CommonUtils/NameConf.h" +#include "CommonUtils/ConfigurableParam.h" +#include "DataFormatsMCH/Cluster.h" +#include "DataFormatsMCH/ROFRecord.h" +#include "DataFormatsMCH/TrackMCH.h" +#include "DataFormatsParameters/GRPObject.h" +#include "DetectorsBase/GeometryManager.h" +#include "DetectorsBase/GRPGeomHelper.h" +#include "DetectorsBase/Propagator.h" +#include "DetectorsCommonDataFormats/AlignParam.h" +#include "DetectorsCommonDataFormats/DetID.h" +#include "DetectorsCommonDataFormats/DetectorNameConf.h" +#include "MathUtils/Cartesian.h" +#include "MCHAlign/Aligner.h" +#include "MCHGeometryTransformer/Transformations.h" +#include "MCHTracking/Track.h" +#include "MCHTracking/TrackExtrap.h" +#include "MCHTracking/TrackParam.h" +#include "MCHTracking/TrackFitter.h" +#include "ReconstructionDataFormats/TrackMCHMID.h" + +namespace o2 +{ +namespace mch +{ + +using namespace std; +using namespace o2::framework; +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; + + 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}, + + {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}, + + {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}, + + {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}, + + {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}, + + {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}, + + {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}, + + {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}, + + {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} + + }; + std::vector params; + std::vector errors; + std::vector pulls; + + constexpr double pi() { return 3.14159265358979323846; } + + //_________________________________________________________________________________________________ + AlignmentTask(shared_ptr req) : mCCDBRequest(req) {} + + //_________________________________________________________________________________________________ + void init(framework::InitContext& ic) + { + + LOG(info) << "Initializing aligner"; + // Initialize alignment algorithm + + int numberOfGlobalParam = 624; + double default_value = 0; + params = std::vector(numberOfGlobalParam, default_value); + errors = std::vector(numberOfGlobalParam, default_value); + pulls = std::vector(numberOfGlobalParam, default_value); + + 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"); + 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 (std::filesystem::exists(grpFile)) { + const auto grp = parameters::GRPObject::loadFrom(grpFile); + base::Propagator::initFieldFromGRP(grp); + TrackExtrap::setField(); + mAlign.SetBFieldOn(TrackExtrap::isFieldON()); + TrackExtrap::useExtrapV2(); + } else { + LOG(fatal) << "No GRP file"; + } + + auto geoIdealFile = ic.options().get("geo-file-ideal"); + if (std::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 (std::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); + } + + doMatched = ic.options().get("matched"); + outFileName = ic.options().get("output"); + readFromRec = ic.options().get("use-record"); + + if (readFromRec) { + mAlign.SetReadOnly(); + LOG(info) << "Reading records as input"; + } + mAlign.init(); + + 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)) { + TrackExtrap::setField(); + mAlign.SetBFieldOn(TrackExtrap::isFieldON()); + TrackExtrap::useExtrapV2(); + } + + 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; + + auto mchTrack = mchTracks.at(iMCHTrack); + int id_track = iMCHTrack; + int nb_clusters = mchTrack.getNClusters(); + + // Track selection, considering only tracks having at least 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); + + // Erase removable track + if (RemoveTrack(convertedTrack, ImproveCut)) { + continue; + } else { + tracksGood += 1; + } + + // Track processing, saving residuals + o2::fwdalign::MillePedeRecord* mchRecord = mAlign.ProcessTrack(convertedTrack, transformation, + doAlign, weightRecord); + } + } + } + + //_________________________________________________________________________________________________ + void processWithOutMatching(vector& mchROFs, vector& mchTracks, vector& mchClusters) + { + + // processing for each track + for (const auto& mchROF : mchROFs) { + + 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; + + // 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); + + // Erase removable track + if (RemoveTrack(convertedTrack, ImproveCut)) { + continue; + } else { + tracksGood += 1; + } + + // Track processing, saving residuals + o2::fwdalign::MillePedeRecord* mchRecord = mAlign.ProcessTrack(convertedTrack, transformation, + doAlign, weightRecord); + } + } + } + + //_________________________________________________________________________________________________ + 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) << "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 + 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"); + } + + // Store aligned geometry + gGeoManager->Export(Geo_file.c_str()); + + // Store param plots + drawHisto(params, errors, pulls, outFileName); + + // 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++) { + + ParamsRef.emplace_back(ParamsTrack.at(fgSNDetElemHalfCh[hc])); + + for (int de = 0; de < fgNDetElemHalfCh[hc]; de++) { + + int iDEN = fgDetElemHalfCh[hc][de]; + o2::detectors::AlignParam param_Track = ParamsTrack.at(fgSNDetElemHalfCh[hc] + 1 + de); + + LOG(debug) << Form("%s%s", "Processing DET Elem: ", (param_Track.getSymName()).c_str()); + + 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()); + + TGeoHMatrix transRef = transformIdeal[iDEN]; + TGeoHMatrix transTrack = doReAlign ? transformNew[iDEN] : transformRef[iDEN]; + TGeoHMatrix transRefTrack = transTrack * transRef.Inverse(); + TGeoHMatrix delta_ref = delta_track * transRefTrack; + + param_Ref.setSymName((param_Track.getSymName()).c_str()); + param_Ref.setGlobalParams(delta_ref); + param_Ref.applyToGeometry(); + ParamsRef.emplace_back(param_Ref); + } + } + + TFile* fOut = TFile::Open("AlignParam@ideal.root", "RECREATE"); + fOut->WriteObjectAny(&ParamsRef, "std::vector", "alignment"); + fOut->Close(); + } + + //_________________________________________________________________________________________________ + Track MCHFormatConvert(TrackMCH& mchTrack, vector& mchClusters, bool doReAlign) + { + + Track convertedTrack = Track(); + + // 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(); + + // Transformations to new geometry from reference geometry + if (doReAlign) { + + math_utils::Point3D local; + math_utils::Point3D master; + + master.SetXYZ(cluster->getX(), cluster->getY(), cluster->getZ()); + + transformRef[cluster->getDEId()].MasterToLocal(master, local); + transformNew[cluster->getDEId()].LocalToMaster(local, master); + + cluster->x = master.x(); + cluster->y = master.y(); + cluster->z = master.z(); + } + convertedTrack.createParamAtCluster(*cluster); + } + + return Track(convertedTrack); + } + + //_________________________________________________________________________________________________ + bool RemoveTrack(Track& track, double ImproveCut) + { + + double maxChi2Cluster = 2 * ImproveCut * ImproveCut; + bool removeTrack = false; + + try { + trackFitter.fit(track, false); + } catch (exception const& e) { + removeTrack = true; + return removeTrack; + } + + auto itStartingParam = std::prev(track.rend()); + + while (true) { + + try { + trackFitter.fit(track, true, false, (itStartingParam == track.rbegin()) ? nullptr : &itStartingParam); + } catch (exception const&) { + removeTrack = true; + break; + } + + double worstLocalChi2 = -1.0; + + track.tagRemovableClusters(0x1F, false); + + auto itWorstParam = track.end(); + + for (auto itParam = track.begin(); itParam != track.end(); ++itParam) { + if (itParam->getLocalChi2() > worstLocalChi2) { + worstLocalChi2 = itParam->getLocalChi2(); + itWorstParam = itParam; + } + } + + if (worstLocalChi2 < maxChi2Cluster) { + break; + } + + 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(std::vector& params, std::vector& errors, std::vector& pulls, 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("P goff"); + + // 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; + } + } + + PlotFiles->WriteObjectAny(cvn1, "TCanvas", "AlignParam"); + 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; + } + + // 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; + } + + 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 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{}; +}; + +//_________________________________________________________________________________________________ +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); + + return DataProcessorSpec{ + "mch-alignment", + inputSpecs, + 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"}}, + {"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"}}, + {"do-realign", VariantType::Bool, false, {"Switch for re-alignment using another geometry"}}, + {"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"}}}}; +} + +} // namespace mch +} // namespace o2 \ No newline at end of file diff --git a/Detectors/MUON/MCH/Align/src/MCHAlignLinkDef.h b/Detectors/MUON/MCH/Align/src/MCHAlignLinkDef.h new file mode 100644 index 0000000000000..2687cb222464d --- /dev/null +++ b/Detectors/MUON/MCH/Align/src/MCHAlignLinkDef.h @@ -0,0 +1,19 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#ifdef __CLING__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; +#pragma link C++ class o2::mch::Aligner + ; + +#endif diff --git a/Detectors/MUON/MCH/Align/src/alignment-workflow.cxx b/Detectors/MUON/MCH/Align/src/alignment-workflow.cxx new file mode 100644 index 0000000000000..d74cb5fbc6442 --- /dev/null +++ b/Detectors/MUON/MCH/Align/src/alignment-workflow.cxx @@ -0,0 +1,90 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file alignment-workflow.cxx +/// \brief Implementation of a DPL device to perform alignment for muon spectrometer +/// +/// \author Chi ZHANG, CEA-Saclay + +#include "MCHAlign/AlignmentSpec.h" + +#include "CCDB/BasicCCDBManager.h" +#include "CCDB/CCDBTimeStampUtils.h" +#include "CommonUtils/ConfigurableParam.h" +#include "DetectorsRaw/HBFUtils.h" +#include "Headers/STFHeader.h" +#include "Framework/CallbackService.h" +#include "Framework/ConcreteDataMatcher.h" +#include "Framework/ConfigParamRegistry.h" +#include "Framework/ControlService.h" +#include "Framework/DataProcessorSpec.h" +#include "Framework/Lifetime.h" +#include "Framework/Output.h" +#include "Framework/Task.h" +#include "Framework/Logger.h" + +using namespace o2::framework; +using namespace std; + +namespace o2::mch +{ + +class SeederTask : public Task +{ + public: + void run(ProcessingContext& pc) final + { + const auto& hbfu = o2::raw::HBFUtils::Instance(); + auto& tinfo = pc.services().get(); + if (hbfu.startTime != 0) { + tinfo.creation = hbfu.startTime; + } else { + tinfo.creation = std::chrono::time_point_cast(std::chrono::system_clock::now()).time_since_epoch().count(); + } + if (hbfu.orbitFirstSampled != 0) { + tinfo.firstTForbit = hbfu.orbitFirstSampled; + } else { + tinfo.firstTForbit = 0; + } + auto& stfDist = pc.outputs().make(Output{"FLP", "DISTSUBTIMEFRAME", 0}); + pc.services().get().endOfStream(); + pc.services().get().readyToQuit(QuitRequest::Me); + } +}; + +} // namespace o2::mch + +o2::framework::DataProcessorSpec getSeederSpec() +{ + return DataProcessorSpec{ + "seeder", + Inputs{}, + Outputs{{"FLP", "DISTSUBTIMEFRAME", 0}}, + AlgorithmSpec{o2::framework::adaptFromTask()}, + Options{}}; +} + +// we need to add workflow options before including Framework/runDataProcessing +void customize(vector& workflowOptions) +{ + workflowOptions.emplace_back("configKeyValues", VariantType::String, "", + ConfigParamSpec::HelpString{"Semicolon separated key=value strings"}); + workflowOptions.emplace_back("disable-input-from-ccdb", VariantType::Bool, false, + ConfigParamSpec::HelpString{"Do not read magnetic field and geometry from CCDB"}); +} + +#include "Framework/runDataProcessing.h" +WorkflowSpec defineDataProcessing(const ConfigContext& configcontext) +{ + o2::conf::ConfigurableParam::updateFromString(configcontext.options().get("configKeyValues")); + bool disableCCDB = configcontext.options().get("disable-input-from-ccdb"); + return WorkflowSpec{o2::mch::getAlignmentSpec(disableCCDB), getSeederSpec()}; +} diff --git a/Detectors/MUON/MCH/CMakeLists.txt b/Detectors/MUON/MCH/CMakeLists.txt index f34cb0fe437a8..21fe99b5e3e21 100644 --- a/Detectors/MUON/MCH/CMakeLists.txt +++ b/Detectors/MUON/MCH/CMakeLists.txt @@ -8,7 +8,7 @@ # In applying this license CERN does not waive the privileges and immunities # granted to it by virtue of its status as an Intergovernmental Organization # or submit itself to any jurisdiction. - +add_subdirectory(Align) add_subdirectory(Base) add_subdirectory(Constants) add_subdirectory(Contour)