From e059886f37eb875b9d9e597709f999547ba9f919 Mon Sep 17 00:00:00 2001 From: sevdokim <30727271+sevdokim@users.noreply.github.com> Date: Sun, 11 Apr 2021 15:49:36 +0200 Subject: [PATCH] CPV fix simulation workflow (#5887) * Fixed CPV simulation workflow * applied clang-format * fix for ~statement should be inside braces~ mumbling Co-authored-by: sevdokim --- .../CPV/base/include/CPVBase/CPVSimParams.h | 40 +++--- .../CPV/calib/include/CPVCalib/Pedestals.h | 13 +- Detectors/CPV/calib/src/CalibParams.cxx | 2 +- Detectors/CPV/calib/src/Pedestals.cxx | 38 +++++- .../include/CPVSimulation/Digitizer.h | 13 +- .../include/CPVSimulation/RawWriter.h | 2 + Detectors/CPV/simulation/src/Detector.cxx | 1 - Detectors/CPV/simulation/src/Digitizer.cxx | 80 ++++++++---- Detectors/CPV/simulation/src/RawWriter.cxx | 24 +++- Detectors/CPV/testsimulation/CMakeLists.txt | 5 + Detectors/CPV/testsimulation/plot_clu_cpv.C | 118 ++++++++++++++++++ Detectors/CPV/testsimulation/plot_dig_cpv.C | 55 +++++--- Detectors/CPV/testsimulation/plot_hit_cpv.C | 2 +- .../CPVWorkflow/RawToDigitConverterSpec.h | 3 + .../workflow/src/RawToDigitConverterSpec.cxx | 47 ++++++- 15 files changed, 359 insertions(+), 84 deletions(-) create mode 100644 Detectors/CPV/testsimulation/plot_clu_cpv.C diff --git a/Detectors/CPV/base/include/CPVBase/CPVSimParams.h b/Detectors/CPV/base/include/CPVBase/CPVSimParams.h index 818b498662c5e..eb83ad2431ce3 100644 --- a/Detectors/CPV/base/include/CPVBase/CPVSimParams.h +++ b/Detectors/CPV/base/include/CPVBase/CPVSimParams.h @@ -30,7 +30,7 @@ struct CPVSimParams : public o2::conf::ConfigurableParamHelper { float mPadSizeX = 1.13; ///< overall size of CPV active size float mPadSizeZ = 2.1093; ///< in phi and z directions float mDetR = 0.1; ///< Relative energy fluctuation in track for 100 e- - float mdEdx = 4.0; ///< Average energy loss in CPV; + float mdEdx = 400.0; ///< Average energy loss in CPV (arbitrary units); int mNgamz = 5; ///< Ionization size in Z int mNgamx = 9; ///< Ionization size in Phi float mCPVGasThickness = 1.3; ///< width of ArC02 gas gap @@ -38,26 +38,28 @@ struct CPVSimParams : public o2::conf::ConfigurableParamHelper { float mB = 0.7; ///< Parameter to model CPV response //Parameters used in electronic noise calculation and thresholds (Digitizer) - float mReadoutTime = 5.; ///< Read-out time in ns for default simulaionts - float mDeadTime = 20.; ///< PHOS dead time (includes Read-out time i.e. mDeadTime>=mReadoutTime) - float mReadoutTimePU = 2000.; ///< Read-out time in ns if pileup simulation on in DigitizerSpec - float mDeadTimePU = 30000.; ///< PHOS dead time if pileup simulation on in DigitizerSpec - bool mApplyDigitization = true; ///< if energy digitization should be applied - float mZSthreshold = 0.01; ///< Zero Suppression threshold - float mADCWidth = 0.005; ///< Widht of ADC channel used for energy digitization - float mNoise = 0.01; ///< charge noise in one pad - float mCoeffToNanoSecond = 1.e+9; ///< Conversion for time units - float mSortingDelta = 0.1; ///< used in sorting clusters inverse sorting band in cm + float mReadoutTime = 5.; ///< Read-out time in ns for default simulaionts + float mDeadTime = 20.; ///< PHOS dead time (includes Read-out time i.e. mDeadTime>=mReadoutTime) + float mReadoutTimePU = 2000.; ///< Read-out time in ns if pileup simulation on in DigitizerSpec + float mDeadTimePU = 30000.; ///< PHOS dead time if pileup simulation on in DigitizerSpec + bool mApplyDigitization = true; ///< if energy digitization should be applied + //float mZSthreshold = 0.01; ///< Zero Suppression threshold + float mZSnSigmas = 3.; ///< Zero Suppression threshold + //float mADCWidth = 0.005; ///< Widht of ADC channel used for energy digitization + //float mNoise = 0.01; ///< charge noise in one pad + //float mCoeffToNanoSecond = 1.e+9; ///< Conversion for time units + float mSortingDelta = 0.1; ///< used in sorting clusters inverse sorting band in cm //Parameters used in clusterization - float mDigitMinEnergy = 0.01; ///< Minimal amplitude of a digit to be used in cluster - float mClusteringThreshold = 0.050; ///< Seed digit minimal amplitude - float mUnfogingEAccuracy = 1.e-3; ///< Accuracy of energy calculation in unfoding prosedure (GeV) - float mUnfogingXZAccuracy = 1.e-1; ///< Accuracy of position calculation in unfolding procedure (cm) - float mLocalMaximumCut = 0.030; ///< Threshold to separate local maxima - float mLogWeight = 4.5; ///< weight in cluster center of gravity calculation - int mNMaxIterations = 10; ///< Maximal number of iterations in unfolding procedure - bool mUnfoldClusters = false; ///< Perform cluster unfolding? + //float mDigitMinEnergy = 0.01; ///< Minimal amplitude of a digit to be used in cluster + float mDigitMinEnergy = 5.; ///< Minimal amplitude of a digit to be used in cluster + float mClusteringThreshold = 10.; ///< Seed digit minimal amplitude + float mUnfogingEAccuracy = 1.e-3; ///< Accuracy of energy calculation in unfoding prosedure (GeV) + float mUnfogingXZAccuracy = 1.e-1; ///< Accuracy of position calculation in unfolding procedure (cm) + float mLocalMaximumCut = 0.030; ///< Threshold to separate local maxima + float mLogWeight = 4.5; ///< weight in cluster center of gravity calculation + int mNMaxIterations = 10; ///< Maximal number of iterations in unfolding procedure + bool mUnfoldClusters = false; ///< Perform cluster unfolding? inline float CellWr() const { return 0.5 * mPadSizeX; } ///< Distance between wires (2 wires above 1 pad) diff --git a/Detectors/CPV/calib/include/CPVCalib/Pedestals.h b/Detectors/CPV/calib/include/CPVCalib/Pedestals.h index a23f85536c901..911f03054c2db 100644 --- a/Detectors/CPV/calib/include/CPVCalib/Pedestals.h +++ b/Detectors/CPV/calib/include/CPVCalib/Pedestals.h @@ -22,6 +22,7 @@ #include "TObject.h" class TH1; +class TH1F; namespace o2 { @@ -45,22 +46,26 @@ class Pedestals /// \param cellID Absolute ID of cell /// \return pedestal for the cell short getPedestal(short cellID) const { return short(mPedestals.at(cellID)); } + float getPedSigma(short cellID) const { return mPedSigmas.at(cellID); } /// \brief Set pedestal /// \param cellID Absolute ID of cell /// \param c is the pedestal (expected to be in range <254) - void setPedestal(short cellID, short c) { mPedestals[cellID] = static_cast(c); } + void setPedestal(short cellID, short c) { mPedestals[cellID] = (c > 0 && c < 511) ? c : mPedestals[cellID]; } + void setPedSigma(short cellID, float c) { mPedSigmas[cellID] = (c > 0) ? c : mPedSigmas[cellID]; } /// \brief Set pedestals from 1D histogram with cell absId in x axis /// \param 1D(NCHANNELS) histogram with calibration coefficients /// \return Is successful bool setPedestals(TH1* h); + bool setPedSigmas(TH1F* h); private: - static constexpr short NCHANNELS = 23040; ///< Number of channels in 3 modules starting from 0 - std::array mPedestals; ///< Container for pedestals + static constexpr short NCHANNELS = 23040; ///< Number of channels in 3 modules starting from 0 + std::array mPedestals; ///< Container for pedestals + std::array mPedSigmas; ///< Container for pedestal sigmas - ClassDefNV(Pedestals, 1); + ClassDefNV(Pedestals, 2); }; } // namespace cpv diff --git a/Detectors/CPV/calib/src/CalibParams.cxx b/Detectors/CPV/calib/src/CalibParams.cxx index 1be71c1b260b2..4276a5d62b15f 100644 --- a/Detectors/CPV/calib/src/CalibParams.cxx +++ b/Detectors/CPV/calib/src/CalibParams.cxx @@ -22,7 +22,7 @@ using namespace o2::cpv; CalibParams::CalibParams(short /*dummy*/) { //produce reasonable objest for test purposes - mGainCalib.fill(0.01); + mGainCalib.fill(1.); } bool CalibParams::setGain(TH2* h, short module) diff --git a/Detectors/CPV/calib/src/Pedestals.cxx b/Detectors/CPV/calib/src/Pedestals.cxx index eb04ce5ff594e..3f4f6ee39e811 100644 --- a/Detectors/CPV/calib/src/Pedestals.cxx +++ b/Detectors/CPV/calib/src/Pedestals.cxx @@ -11,15 +11,17 @@ #include "CPVCalib/Pedestals.h" #include "FairLogger.h" #include +#include using namespace o2::cpv; Pedestals::Pedestals(int /*dummy*/) { //produce reasonable objest for test purposes - mPedestals.fill(40); + mPedestals.fill(200); //typical pedestal value + mPedSigmas.fill(1.5); //typical pedestal sigma } - +//______________________________________________________________________________ bool Pedestals::setPedestals(TH1* h) { if (!h) { @@ -33,11 +35,37 @@ bool Pedestals::setPedestals(TH1* h) } for (short i = 1; i <= NCHANNELS; i++) { - if (h->GetBinContent(i) > 255) { - LOG(ERROR) << "pedestal value too large:" << h->GetBinContent(i) << "can not be stored in char"; + if (h->GetBinContent(i) > 511) { + LOG(ERROR) << "setPedestals : pedestal value = " << h->GetBinContent(i) + << " in channel " << i + << " exceeds max possible value 511 (limited by CPV electronics)"; + continue; + } + mPedestals[i] = short(h->GetBinContent(i)); + } + return true; +} +//_______________________________________________________________________________ +bool Pedestals::setPedSigmas(TH1F* h) +{ + if (!h) { + LOG(ERROR) << "no input histogam"; + return false; + } + + if (h->GetNbinsX() != NCHANNELS) { + LOG(ERROR) << "Wrong dimentions of input histogram:" << h->GetNbinsX() << " instead of " << NCHANNELS; + return false; + } + + for (short i = 1; i <= NCHANNELS; i++) { + if (h->GetBinContent(i) < 0) { + LOG(ERROR) << "pedestal sigma = " << h->GetBinContent(i) + << " in channel " << i + << " cannot be less than 0"; continue; } - mPedestals[i] = char(h->GetBinContent(i)); + mPedSigmas[i] = float(h->GetBinContent(i)); } return true; } diff --git a/Detectors/CPV/simulation/include/CPVSimulation/Digitizer.h b/Detectors/CPV/simulation/include/CPVSimulation/Digitizer.h index 10dfc44b893b3..adb199b353a27 100644 --- a/Detectors/CPV/simulation/include/CPVSimulation/Digitizer.h +++ b/Detectors/CPV/simulation/include/CPVSimulation/Digitizer.h @@ -14,6 +14,8 @@ #include "DataFormatsCPV/Digit.h" #include "CPVBase/Geometry.h" #include "CPVCalib/CalibParams.h" +#include "CPVCalib/Pedestals.h" +#include "CPVCalib/BadChannelMap.h" #include "CPVBase/Hit.h" #include "SimulationDataFormat/MCCompLabel.h" #include "SimulationDataFormat/MCTruthContainer.h" @@ -39,15 +41,16 @@ class Digitizer : public TObject int source, int entry, double dt); protected: - float uncalibrate(float e, int absId); - float simulateNoise(); + float simulatePedestalNoise(int absId); private: static constexpr short NCHANNELS = 23040; //128*60*3: toatl number of CPV channels std::unique_ptr mCalibParams; /// Calibration coefficients - std::array mArrayD; - - ClassDefOverride(Digitizer, 2); + std::unique_ptr mPedestals; /// Pedestals + std::unique_ptr mBadMap; /// Bad channel map + std::array mArrayD; ///array of digits (for inner use) + std::array mDigitThresholds; + ClassDefOverride(Digitizer, 3); }; } // namespace cpv } // namespace o2 diff --git a/Detectors/CPV/simulation/include/CPVSimulation/RawWriter.h b/Detectors/CPV/simulation/include/CPVSimulation/RawWriter.h index 0083572838c84..ec407f91bf0dd 100644 --- a/Detectors/CPV/simulation/include/CPVSimulation/RawWriter.h +++ b/Detectors/CPV/simulation/include/CPVSimulation/RawWriter.h @@ -26,6 +26,7 @@ #include "DataFormatsCPV/Digit.h" #include "DataFormatsCPV/TriggerRecord.h" #include "CPVCalib/CalibParams.h" +#include "CPVCalib/Pedestals.h" namespace o2 { @@ -78,6 +79,7 @@ class RawWriter FileFor_t mFileFor = FileFor_t::kFullDet; ///< Granularity of the output files std::string mOutputLocation = "./"; ///< Rawfile name std::unique_ptr mCalibParams; ///< CPV calibration + std::unique_ptr mPedestals; ///< CPV pedestals std::vector mPayload; ///< Payload to be written gsl::span mDigits; ///< Digits input vector - must be in digitized format including the time response std::unique_ptr mRawWriter; ///< Raw writer diff --git a/Detectors/CPV/simulation/src/Detector.cxx b/Detectors/CPV/simulation/src/Detector.cxx index 0bc36e4b0e9ac..98ce6f20f3f38 100644 --- a/Detectors/CPV/simulation/src/Detector.cxx +++ b/Detectors/CPV/simulation/src/Detector.cxx @@ -261,7 +261,6 @@ Bool_t Detector::ProcessHits(FairVolume* v) // Now calculate pad response double qpad = padResponseFunction(qhit, zg, xg); - qpad += cpvparam.mNoise * rnor2; if (qpad < 0) { continue; } diff --git a/Detectors/CPV/simulation/src/Digitizer.cxx b/Detectors/CPV/simulation/src/Digitizer.cxx index 2d0bc19d1227a..a17a9ddf62391 100644 --- a/Detectors/CPV/simulation/src/Digitizer.cxx +++ b/Detectors/CPV/simulation/src/Digitizer.cxx @@ -31,7 +31,8 @@ void Digitizer::init() mCalibParams.reset(new CalibParams(1)); // test default calibration LOG(INFO) << "[CPVDigitizer] No reading calibration from ccdb requested, set default"; } else { - LOG(INFO) << "[CPVDigitizer] can not gey calibration object from ccdb yet"; + LOG(INFO) << "[CPVDigitizer] can not get calibration object from ccdb yet. Using default"; + mCalibParams.reset(new CalibParams(1)); // test default calibration // o2::ccdb::CcdbApi ccdb; // std::map metadata; // do we want to store any meta data? // ccdb.init("http://ccdb-test.cern.ch:8080"); // or http://localhost:8080 for a local installation @@ -41,6 +42,45 @@ void Digitizer::init() // } } } + if (!mPedestals) { + if (o2::cpv::CPVSimParams::Instance().mCCDBPath.compare("localtest") == 0) { + mPedestals.reset(new Pedestals(1)); // test default calibration + LOG(INFO) << "[CPVDigitizer] No reading calibration from ccdb requested, set default"; + } else { + LOG(INFO) << "[CPVDigitizer] can not get pedestal object from ccdb yet. Using default"; + mPedestals.reset(new Pedestals(1)); // test default calibration + // o2::ccdb::CcdbApi ccdb; + // std::map metadata; // do we want to store any meta data? + // ccdb.init("http://ccdb-test.cern.ch:8080"); // or http://localhost:8080 for a local installation + // mPedestals = ccdb.retrieveFromTFileAny("CPV/Calib", metadata, mEventTime); + // if (!mPedestals) { + // LOG(FATAL) << "[CPVDigitizer] can not get calibration object from ccdb"; + // } + } + } + if (!mBadMap) { + if (o2::cpv::CPVSimParams::Instance().mCCDBPath.compare("localtest") == 0) { + mBadMap.reset(new BadChannelMap(1)); // test default calibration + LOG(INFO) << "[CPVDigitizer] No reading calibration from ccdb requested, set default"; + } else { + LOG(INFO) << "[CPVDigitizer] can not get bad channel map object from ccdb yet. Using default"; + mBadMap.reset(new BadChannelMap(1)); // test default calibration + // o2::ccdb::CcdbApi ccdb; + // std::map metadata; // do we want to store any meta data? + // ccdb.init("http://ccdb-test.cern.ch:8080"); // or http://localhost:8080 for a local installation + // mBadMap = ccdb.retrieveFromTFileAny("CPV/Calib", metadata, mEventTime); + // if (!mBadMap) { + // LOG(FATAL) << "[CPVDigitizer] can not get calibration object from ccdb"; + // } + } + } + + //signal thresolds for digits + //note that digits are calibrated objects + for (int i = 0; i < NCHANNELS; i++) { + mDigitThresholds[i] = o2::cpv::CPVSimParams::Instance().mZSnSigmas * + mPedestals->getPedSigma(i) * mCalibParams->getGain(i); + } } //_______________________________________________________________________ @@ -60,14 +100,11 @@ void Digitizer::processHits(const std::vector* hits, const std::vector o2::cpv::CPVSimParams::Instance().mZSthreshold) { - mArrayD[i].setAmplitude(energy); - mArrayD[i].setAbsId(i); - } + float amplitude = simulatePedestalNoise(i); + mArrayD[i].setAmplitude(amplitude); + mArrayD[i].setAbsId(i); } } else { //if digits exist, no noise should be added for (auto& dBg : digitsBg) { //digits are sorted and unique @@ -79,12 +116,12 @@ void Digitizer::processHits(const std::vector* hits, const std::vector 0) { - mArrayD[i].setAmplitude(mArrayD[i].getAmplitude() + uncalibrate(h.GetEnergyLoss(), i)); + mArrayD[i].setAmplitude(mArrayD[i].getAmplitude() + h.GetEnergyLoss()); } else { - mArrayD[i].setAmplitude(uncalibrate(h.GetEnergyLoss(), i)); + mArrayD[i].setAmplitude(h.GetEnergyLoss()); mArrayD[i].setAbsId(i); } - if (mArrayD[i].getAmplitude() > o2::cpv::CPVSimParams::Instance().mZSthreshold) { + if (mArrayD[i].getAmplitude() > mDigitThresholds[i]) { int labelIndex = mArrayD[i].getLabel(); if (labelIndex == -1) { //no digit or noisy labelIndex = labels.getIndexedSize(); @@ -109,23 +146,22 @@ void Digitizer::processHits(const std::vector* hits, const std::vector o2::cpv::CPVSimParams::Instance().mZSthreshold) { + if (!mBadMap->isChannelGood(i)) { + continue; //bad channel -> skip this digit + } + if (mArrayD[i].getAmplitude() > mDigitThresholds[i]) { digitsOut.push_back(mArrayD[i]); } } } -float Digitizer::simulateNoise() { return gRandom->Gaus(0., o2::cpv::CPVSimParams::Instance().mNoise); } - -//_______________________________________________________________________ -float Digitizer::uncalibrate(const float e, const int absId) +float Digitizer::simulatePedestalNoise(int absId) { - // Decalibrate CPV digit, i.e. transform from amplitude to ADC counts a factor read from CDB - float calib = mCalibParams->getGain(absId); - if (calib > 0) { - return e / calib; - } else { - return 0; // TODO apply de-calibration from OCDB + //this function is to simulate pedestal and its noise (ADC counts) + if (absId < 0 || absId >= NCHANNELS) { + return 0.; } + return gRandom->Gaus(0, mPedestals->getPedSigma(absId) * mCalibParams->getGain(absId)); } diff --git a/Detectors/CPV/simulation/src/RawWriter.cxx b/Detectors/CPV/simulation/src/RawWriter.cxx index 533e48d2a6079..ce77135949549 100644 --- a/Detectors/CPV/simulation/src/RawWriter.cxx +++ b/Detectors/CPV/simulation/src/RawWriter.cxx @@ -60,6 +60,28 @@ void RawWriter::digitsToRaw(gsl::span digitsbranch, gsl::span(1); // test default calibration + LOG(INFO) << "[RawWriter] No reading calibration from ccdb requested, set default"; + } else { + LOG(INFO) << "[RawWriter] getting calibration object from ccdb"; + o2::ccdb::CcdbApi ccdb; + std::map metadata; + ccdb.init("http://ccdb-test.cern.ch:8080"); // or http://localhost:8080 for a local installation + auto tr = triggerbranch.begin(); + double eventTime = -1; + // if(tr!=triggerbranch.end()){ + // eventTime = (*tr).getBCData().getTimeNS() ; + // } + //add copy constructor if necessary + // mPedestals = std::make_unique(ccdb.retrieveFromTFileAny("CPV/Calib", metadata, eventTime)); + if (!mPedestals) { + LOG(FATAL) << "[RawWriter] can not get calibration object from ccdb"; + } + } + } + for (auto trg : triggerbranch) { processTrigger(digitsbranch, trg); } @@ -84,7 +106,7 @@ bool RawWriter::processTrigger(const gsl::span digitsbranch, con o2::cpv::Geometry::absIdToHWaddress(absId, ccId, dil, gas, pad); //Convert Amp to ADC counts - short charge = dig.getAmplitude() / mCalibParams->getGain(absId); + short charge = dig.getAmplitude() / mCalibParams->getGain(absId) + mPedestals->getPedestal(absId); if (charge > 4095) { charge = 4095; } diff --git a/Detectors/CPV/testsimulation/CMakeLists.txt b/Detectors/CPV/testsimulation/CMakeLists.txt index c79c68e23f613..8a8dd8f9da37c 100644 --- a/Detectors/CPV/testsimulation/CMakeLists.txt +++ b/Detectors/CPV/testsimulation/CMakeLists.txt @@ -21,3 +21,8 @@ o2_add_test_root_macro(plot_hit_cpv.C o2_add_test_root_macro(plot_dig_cpv.C PUBLIC_LINK_LIBRARIES FairRoot::Base O2::CPVSimulation LABELS cpv) + +o2_add_test_root_macro(plot_clu_cpv.C + PUBLIC_LINK_LIBRARIES FairRoot::Base O2::CPVSimulation + LABELS cpv) + diff --git a/Detectors/CPV/testsimulation/plot_clu_cpv.C b/Detectors/CPV/testsimulation/plot_clu_cpv.C new file mode 100644 index 0000000000000..a3d899640c1bc --- /dev/null +++ b/Detectors/CPV/testsimulation/plot_clu_cpv.C @@ -0,0 +1,118 @@ +#if !defined(__CLING__) || defined(__ROOTCLING__) +#include +#include + +#include "TROOT.h" +#include +#include "TCanvas.h" +#include "TH2.h" +//#include "DataFormatsParameters/GRPObject.h" +#include "FairFileSource.h" +#include "FairLogger.h" +#include "FairRunAna.h" +//#include "FairRuntimeDb.h" +#include "FairParRootFileIo.h" +#include "FairSystemInfo.h" +#include "SimulationDataFormat/MCCompLabel.h" + +#include "DataFormatsCPV/Cluster.h" +#include "CPVBase/Geometry.h" +#endif + +void plot_clu_cpv(std::string inputfile = "cpvclusters.root", int ifirst = 0, int ilast = -1) +{ + // macros to plot CPV clusters + + // Clusters + TFile* file0 = TFile::Open(inputfile.data()); + std::cout << " Open clusters file " << inputfile << std::endl; + TTree* cluTree = (TTree*)gFile->Get("o2sim"); + std::vector* mClustersArray = nullptr; + cluTree->SetBranchAddress("CPVCluster", &mClustersArray); + + if (!mClustersArray) { + std::cout << "CPV clusters not found in the file. Exiting ..." << std::endl; + return; + } + + TH1F* hClusterTotEnergy[5]; + TH1F* hClusterMaxEnergy[5]; + TH1F* hClusterSize[5]; + TH1F* hClusterSizeX[5]; + TH1F* hClusterSizeZ[5]; + + TH2D* vMod[5][1000] = {0}; + int primLabels[5][1000]; + for (int mod = 2; mod < 5; mod++) { + hClusterTotEnergy[mod] = new TH1F(Form("hClusterTotEnergy%d", mod), + Form("Cluster Total Energy mod %d", mod), 10000, 0, 10000); + hClusterMaxEnergy[mod] = new TH1F(Form("hClusterMaxEnergy%d", mod), + Form("Cluster Max Energy mod %d", mod), 10000, 0, 10000); + hClusterSize[mod] = new TH1F(Form("hClusterSize%d", mod), + Form("Cluster Size mod %d", mod), 100, 0, 100); + hClusterSizeX[mod] = new TH1F(Form("hClusterSizeX%d", mod), + Form("Cluster SizeX mod %d", mod), 100, 0, 100); + hClusterSizeZ[mod] = new TH1F(Form("hClusterSizeZ%d", mod), + Form("Cluster SizeZ mod %d", mod), 100, 0, 100); + + for (int j = 0; j < 100; j++) + primLabels[mod][j] = -1; + } + + int nEntries = cluTree->GetEntriesFast(); + if (ilast < 0) + ilast = nEntries; + if (ilast > nEntries) + ilast = nEntries; + + for (int ievent = ifirst; ievent < ilast; ievent++) { + cluTree->GetEvent(ievent); + + std::vector::const_iterator it; + std::cout << "I start cluster cycling (record #" << ievent << " in o2sim tree)" << std::endl; + + for (it = mClustersArray->begin(); it != mClustersArray->end(); it++) { + float en = (*it).getEnergy(); + float posX, posZ; + (*it).getLocalPosition(posX, posZ); + int cluSize = (*it).getMultiplicity(); + int mod = (*it).getModule(); + if (!vMod[mod][0]) { + gROOT->cd(); + vMod[mod][0] = + new TH2D(Form("hMod%d_prim%d", mod, 0), Form("hMod%d_prim%d", mod, 0), + 100, -100., 100., 100, -100., 100.); + } + vMod[mod][0]->Fill(posX, posZ, en); + hClusterTotEnergy[mod]->Fill(en); + hClusterSize[mod]->Fill(cluSize); + } + + std::cout << "I finish cycling clusters" << std::endl; + } + TCanvas* c[5]; + TH2D* box = new TH2D("box", "CPV module", 100, -100., 100., 100, -100., 100.); + TCanvas* cTotEn = new TCanvas(); + cTotEn->Divide(3, 1); + TCanvas* cSize = new TCanvas(); + cSize->Divide(3, 1); + + for (int mod = 2; mod < 5; mod++) { + c[mod] = + new TCanvas(Form("ClusterInMod%d", mod), Form("CPV clusters in module %d", mod), 10 * mod, 0, 600 + 10 * mod, 400); + box->Draw(); + int j = 0; + while (vMod[mod][j]) { + vMod[mod][j]->SetLineColor(j + 1); + if (j == 0) + vMod[mod][j]->Draw("box"); + else + vMod[mod][j]->Draw("boxsame"); + j++; + } + cTotEn->cd(mod - 1); + hClusterTotEnergy[mod]->Draw(); + cSize->cd(mod - 1); + hClusterSize[mod]->Draw(); + } +} diff --git a/Detectors/CPV/testsimulation/plot_dig_cpv.C b/Detectors/CPV/testsimulation/plot_dig_cpv.C index 23e4b98963db0..e66b4ef438221 100644 --- a/Detectors/CPV/testsimulation/plot_dig_cpv.C +++ b/Detectors/CPV/testsimulation/plot_dig_cpv.C @@ -6,6 +6,7 @@ #include #include "TCanvas.h" #include "TH2.h" +#include "TH1.h" //#include "DataFormatsParameters/GRPObject.h" #include "FairFileSource.h" #include "FairLogger.h" @@ -21,11 +22,11 @@ void plot_dig_cpv(int ievent = 0, std::string inputfile = "o2dig.root") { - // macros to plot CPV hits + // macros to plot CPV digits - // Hits - TFile* file0 = TFile::Open("o2dig.root"); - std::cout << " Open hits file " << inputfile << std::endl; + // Digits + TFile* file0 = TFile::Open(inputfile.data()); + std::cout << " Open digits file " << inputfile << std::endl; TTree* hitTree = (TTree*)gFile->Get("o2sim"); std::vector* mDigitsArray = nullptr; hitTree->SetBranchAddress("CPVDigit", &mDigitsArray); @@ -36,47 +37,58 @@ void plot_dig_cpv(int ievent = 0, std::string inputfile = "o2dig.root") } hitTree->GetEvent(ievent); - TH2D* vMod[5][100] = {0}; - int primLabels[5][100]; - for (int mod = 1; mod < 5; mod++) + TH1F* hDigitAmplitude[5]; + TH2D* vMod[5][1000] = {0}; + int primLabels[5][1000]; + for (int mod = 1; mod < 5; mod++) { + hDigitAmplitude[mod] = new TH1F(Form("hDigitAmplitude_mod%d", mod), + Form("Digit amplitudes in module %d", mod), + 4096, 0., 4096.); for (int j = 0; j < 100; j++) primLabels[mod][j] = -1; + } std::vector::const_iterator it; short relId[3]; + std::cout << "I start digit cycling" << std::endl; for (it = mDigitsArray->begin(); it != mDigitsArray->end(); it++) { short absId = (*it).getAbsId(); float en = (*it).getAmplitude(); int lab = (*it).getLabel(); //TODO + //std::cout << "label = " << lab << std::endl; o2::cpv::Geometry::absToRelNumbering(absId, relId); + hDigitAmplitude[relId[0]]->Fill(en); // check, if this label already exist int j = 0; bool found = false; - while (primLabels[relId[0]][j] >= -2) { - if (primLabels[relId[0]][j] == lab) { - found = true; - break; - } else { - j++; - } - } + //no labels for the time being + // while (primLabels[relId[0]][j] >= -2) { + // if (primLabels[relId[0]][j] == lab) { + // found = true; + // break; + // } else { + // j++; + // } + // } if (!found) { primLabels[relId[0]][j] = lab; } if (!vMod[relId[0]][j]) { gROOT->cd(); vMod[relId[0]][j] = - new TH2D(Form("hMod%d_prim%d", relId[0], j), Form("hMod%d_prim%d", relId[0], j), 60, 0., 60., 128, 0., 128.); + new TH2D(Form("hMod%d_prim%d", relId[0], j), Form("hMod%d_prim%d", relId[0], j), 128, 0., 128., 60, 0., 60.); } vMod[relId[0]][j]->Fill(relId[1] - 0.5, relId[2] - 0.5, en); } - TCanvas* c[5]; - TH2D* box = new TH2D("box", "CPV module", 60, 0., 60., 128, 0., 128.); - for (int mod = 1; mod < 5; mod++) { + std::cout << "I finish cycling digits" << std::endl; + + TCanvas *c[5], *c_ampl[5]; + TH2D* box = new TH2D("box", "CPV module", 128, 0., 128., 60, 0., 60.); + for (int mod = 2; mod < 5; mod++) { c[mod] = - new TCanvas(Form("DigitInMod%d", mod), Form("CPV hits in module %d", mod), 10 * mod, 0, 600 + 10 * mod, 400); + new TCanvas(Form("DigitInMod%d", mod), Form("CPV digits in module %d", mod), 10 * mod, 0, 600 + 10 * mod, 400); box->Draw(); int j = 0; while (vMod[mod][j]) { @@ -87,5 +99,8 @@ void plot_dig_cpv(int ievent = 0, std::string inputfile = "o2dig.root") vMod[mod][j]->Draw("boxsame"); j++; } + c_ampl[mod] = new TCanvas(Form("DigitAmplitudes_%d", mod), Form("DigitAmplitudes_%d", mod), + 10 * mod + 800, 0, 600 + 10 * mod + 200, 400); + hDigitAmplitude[mod]->Draw(); } } diff --git a/Detectors/CPV/testsimulation/plot_hit_cpv.C b/Detectors/CPV/testsimulation/plot_hit_cpv.C index cb883f2fa5f91..5be311e7eb953 100644 --- a/Detectors/CPV/testsimulation/plot_hit_cpv.C +++ b/Detectors/CPV/testsimulation/plot_hit_cpv.C @@ -72,7 +72,7 @@ void plot_hit_cpv(int ievent = 0, std::string inputprefix = "o2sim") if (!vMod[relId[0]][j]) { gROOT->cd(); vMod[relId[0]][j] = - new TH2D(Form("hMod%d_prim%d", relId[0], j), Form("hMod%d_prim%d", relId[0], j), 60, 0., 60., 128, 0., 128.); + new TH2D(Form("hMod%d_prim%d", relId[0], j), Form("hMod%d_prim%d", relId[0], j), 128, 0., 128., 60, 0., 60.); } vMod[relId[0]][j]->Fill(relId[1] - 0.5, relId[2] - 0.5, en); } diff --git a/Detectors/CPV/workflow/include/CPVWorkflow/RawToDigitConverterSpec.h b/Detectors/CPV/workflow/include/CPVWorkflow/RawToDigitConverterSpec.h index 78764982c9ebe..d44ce0aa8bc6b 100644 --- a/Detectors/CPV/workflow/include/CPVWorkflow/RawToDigitConverterSpec.h +++ b/Detectors/CPV/workflow/include/CPVWorkflow/RawToDigitConverterSpec.h @@ -16,6 +16,7 @@ #include "DataFormatsCPV/TriggerRecord.h" #include "CPVCalib/CalibParams.h" #include "CPVCalib/BadChannelMap.h" +#include "CPVCalib/Pedestals.h" #include "CPVReconstruction/RawDecoder.h" namespace o2 @@ -62,7 +63,9 @@ class RawToDigitConverterSpec : public framework::Task private: int mDDL = 15; + bool mIsPedestalData; ///< No subtract pedestals if true std::unique_ptr mCalibParams; ///< CPV calibration + std::unique_ptr mPedestals; ///< CPV pedestals std::unique_ptr mBadMap; ///< BadMap std::vector mOutputDigits; ///< Container with output cells std::vector mOutputTriggerRecords; ///< Container with output cells diff --git a/Detectors/CPV/workflow/src/RawToDigitConverterSpec.cxx b/Detectors/CPV/workflow/src/RawToDigitConverterSpec.cxx index 3c18ba36b3a83..365421ed4c066 100644 --- a/Detectors/CPV/workflow/src/RawToDigitConverterSpec.cxx +++ b/Detectors/CPV/workflow/src/RawToDigitConverterSpec.cxx @@ -29,6 +29,15 @@ void RawToDigitConverterSpec::init(framework::InitContext& ctx) { mDDL = ctx.options().get("DDL"); LOG(DEBUG) << "Initialize converter "; + + //Read command-line options + //Pedestal flag (on/off) + std::string optPedestal(""); + if (ctx.options().isSet("pedestal")) { + optPedestal = ctx.options().get("pedestal"); + } + LOG(INFO) << "Pedestal data: " << optPedestal; + mIsPedestalData = optPedestal == "on" ? true : false; } void RawToDigitConverterSpec::run(framework::ProcessingContext& ctx) @@ -81,6 +90,28 @@ void RawToDigitConverterSpec::run(framework::ProcessingContext& ctx) } } + if (!mPedestals && !mIsPedestalData) { + if (o2::cpv::CPVSimParams::Instance().mCCDBPath.compare("localtest") == 0) { + mPedestals = std::make_unique(1); // test default calibration + LOG(INFO) << "No reading calibration from ccdb requested, set default"; + } else { + LOG(INFO) << "Getting calibration object from ccdb"; + //TODO: configuring ccdb address from config file, readign proper calibration/BadMap and updateing if necessary + o2::ccdb::CcdbApi ccdb; + std::map metadata; + ccdb.init("http://ccdb-test.cern.ch:8080"); // or http://localhost:8080 for a local installation + // auto tr = triggerbranch.begin(); + // double eventTime = -1; + // if(tr!=triggerbranch.end()){ + // eventTime = (*tr).getBCData().getTimeNS() ; + // } + // mPedestals = ccdb.retrieveFromTFileAny("CPV/Calib", metadata, eventTime); + // if (!mPedestals) { + // LOG(FATAL) << "Can not get calibration object from ccdb"; + // } + } + } + for (const auto& rawData : framework::InputRecordWalker(ctx.inputs())) { // enum RawErrorType_t { @@ -155,12 +186,18 @@ void RawToDigitConverterSpec::run(framework::ProcessingContext& ctx) for (uint32_t adch : decoder.getDigits()) { AddressCharge ac = {adch}; unsigned short absId = ac.Address; - //test bad map - if (mBadMap->isChannelGood(absId)) { - if (ac.Charge > o2::cpv::CPVSimParams::Instance().mZSthreshold) { - float amp = mCalibParams->getGain(absId) * ac.Charge; - currentDigitContainer->emplace_back(absId, amp, -1); + //if we deal with non-pedestal data? + if (!mIsPedestalData) { //not a pedestal data + //test bad map + if (mBadMap->isChannelGood(absId)) { + //we need to subtract pedestal from amplidute and calibrate it + float amp = mCalibParams->getGain(absId) * (ac.Charge - mPedestals->getPedestal(absId)); + if (amp > 0) { + currentDigitContainer->emplace_back(absId, amp, -1); + } } + } else { //pedestal data, no calibration needed. + currentDigitContainer->emplace_back(absId, (float)ac.Charge, -1); } } //Check and send list of hwErrors