From ade17cd641f677475ad66099dad6280b077aa0b4 Mon Sep 17 00:00:00 2001 From: Matthias Kleiner Date: Fri, 22 Jan 2021 17:14:21 +0100 Subject: [PATCH 1/2] adding macro to create space charge density maps from simulations For each created space charge density map the IDCs (3D, 1D, 0D) and the distortions/corrections are calculated and stored in root files. This macro can for example be used to create the input data for the ML algorithms which will be used to correct the space charge distortions. The input to create the density maps are hits from an O2 simulation ("o2sim_HitsTPC.root"). --- Detectors/TPC/spacecharge/CMakeLists.txt | 13 + .../macro/createSCHistosFromHits.C | 815 ++++++++++++++++++ 2 files changed, 828 insertions(+) create mode 100644 Detectors/TPC/spacecharge/macro/createSCHistosFromHits.C diff --git a/Detectors/TPC/spacecharge/CMakeLists.txt b/Detectors/TPC/spacecharge/CMakeLists.txt index a1eb18a71242d..df3b60aa764ff 100644 --- a/Detectors/TPC/spacecharge/CMakeLists.txt +++ b/Detectors/TPC/spacecharge/CMakeLists.txt @@ -37,6 +37,19 @@ o2_add_test_root_macro(macro/createResidualDistortionObject.C O2::CommonUtils LABELS tpc) + +install(FILES macro/createSCHistosFromHits.C + DESTINATION share/macro/) + +o2_add_test_root_macro(macro/createSCHistosFromHits.C + PUBLIC_LINK_LIBRARIES O2::TPCSpaceCharge + O2::CommonUtils + O2::TPCBase + O2::TPCSimulation + O2::DataFormatsTPC + O2::MathUtils + LABELS tpc) + o2_add_test(PoissonSolver COMPONENT_NAME spacecharge PUBLIC_LINK_LIBRARIES O2::TPCSpaceCharge diff --git a/Detectors/TPC/spacecharge/macro/createSCHistosFromHits.C b/Detectors/TPC/spacecharge/macro/createSCHistosFromHits.C new file mode 100644 index 0000000000000..33d549d6eed15 --- /dev/null +++ b/Detectors/TPC/spacecharge/macro/createSCHistosFromHits.C @@ -0,0 +1,815 @@ +// Copyright CERN and copyright holders of ALICE O2. This software is +// distributed under the terms of the GNU General Public License v3 (GPL +// Version 3), copied verbatim in the file "COPYING". +// +// See http://alice-o2.web.cern.ch/license for full licensing information. +// +// 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 createSCHistosFromHits +/// \brief This macro implements the creation of space charge density maps from hit files. + +/* + Macro to create random space-charge density distributions from simulated TPC hits and calculated the corresponding integrated digital currents. The density distributions are used to calculated the space-charge distortions. + + Parallization (omp) can only be used if compiled with g++ and not with ROOT! Uncommented for now + Example: +g++ -o createSCHistosFromHits createSCHistosFromHits.C -I ~/alice/sw/osx_x86-64/FairLogger/latest/include -L ~/alice/sw/osx_x86-64/FairLogger/latest/lib -I /Users/matthias/alice/sw/osx_x86-64/FairRoot/latest/include/ -I$ROOTSYS/include -L$ROOTSYS/lib -lRIO -lCore -std=c++17 -I$O2_ROOT/include -L$O2_ROOT/lib -I$O2_ROOT/include/GPU -I /Users/matthias/alice/sw/osx_x86-64/Vc/latest/include/ -lO2TPCBase -lO2TPCSimulation -lO2CommonUtils -lHist -lMathCore -O3 -Xpreprocessor -fopenmp -I/usr/local/include -L/usr/local/lib -lomp -lfairlogger -lboost_system -lboost_filesystem + + Input: + - One o2sim_HitsTPC.root file for first studies + - List of o2sim_HitsTPC.root files + - Number of ion pile-up events to be used + - Ion drift velocity + + Outputs: + - CalDet Objects with 3D-IDCs, std::vector with 1D-IDCs, std::vector with 0D-IDCs (the vector has size 1), + - Granularity? + - Rows, pads at first + - Histograms with space-charge density + - Granularity + - Different for IDCs + - Same as distortions + - Objects of type SpaceCharge<> containing the density, global distortions, global corrections + + Algorithm: + 1) Choose n events randomly from list of o2sim_HitsTPC.root files + - n = number of ion pile-up events + - Assign global indices to single files + + 2) Place r-phi projection of hits randomly in z + - Take into account A, C side + - Electron transport + - Added: diffusion, average space-charge distortions (if provided) + - Could be added: static distortions, electron drift-time + - Ion transport + - Along straight lines + - Could be added: average space-charge distortions, static distortions + - Space charge: + - Apply epsilon + - Convert to charge (C) + - epsilon variations from IBF-map (has to be provided as input) + - IDC: + - Apply effective gain and digitize + - Convert z position to time of IDC measurement + - Gain variations from gain-map (has to be provided as input) + - Grouping of IDCs has not yet be determined + - Same or different for IDC and SC density? + - Integration over 1 ms steps in time + - r, phi: rows, pads + - Average over several rows, pads + + 3) Primary ionization + - Transport of ions + - Along straight lines + - Could be added: average space-charge distortions, static distortions + - 3D hit positions as starting points + - Calculate end points using z position of IBF slice + - Could be added: electron drift-time + + 4) Calculate distortions and corrections + - Corrections to be used for calibration studies + - Active area / volume only? + - Dead regions due to radial distortions + - Do we care about |eta|>0.9 ? + + TODO: + - replace hard coded constants (named constXXX and global constants) by O2 constants + - CalDet objects for IDCs? + - Array of nZ objects, one for each 1 ms + - Variation of ion drift time in IDCs and SC density? + - Make many smaller (~1000 event) CalDets / histos / maps and combine later? DONE + - For the final data for ML, first use only A side values for IDCs, density, corrections + - Later provide data for both sides? Fixed boundary at z = 0 required. + + Usage: + .L $SC/createSCHistosFromHits.C+ + + root -b -l -q $SC/createSCHistosFromHits.C+ +*/ + +#if !defined(__CLING__) || defined(__ROOTCLING__) +// root includes +#include "TFile.h" +#include "TH3F.h" +#include "TMath.h" +#include "TObjArray.h" +#include "TRandom.h" +#include "TString.h" +#include "TSystem.h" +#include "TTree.h" + +// O2 includes +#include "CommonUtils/TreeStreamRedirector.h" +#include "MathUtils/Utils.h" +#include "TPCBase/CalDet.h" +#include "TPCBase/CRU.h" +#include "TPCBase/Mapper.h" +#include "TPCBase/PadPos.h" +#include "TPCBase/ParameterDetector.h" +#include "TPCBase/ParameterGEM.h" +#include "TPCBase/ParameterElectronics.h" +#include "TPCBase/CDBInterface.h" +#include "TPCSimulation/ElectronTransport.h" +#include "TPCSimulation/GEMAmplification.h" +#include "TPCSimulation/SAMPAProcessing.h" +#include "TPCSimulation/Point.h" +#include "TPCSpaceCharge/SpaceCharge.h" +#include "TPCSpaceCharge/PoissonSolverHelpers.h" +#include "DataFormatsTPC/Defs.h" +#include "TPCSpaceCharge/SpaceChargeHelpers.h" + +#include +// #include // comment in if using multithreading +#endif + +using namespace o2::tpc; +using CalPad = CalDet; +using CalPadArr = CalArray; + +// Physics parameters +const float mEpsilon0 = o2::tpc::TPCParameters::E0 * 0.01; //8.854187817e-14; // vacuum permittivity [A·s/(V·cm)] + +// Gas parameters +// const float mK0 = 2.92f; // in cm^2 / (Vs), reduced ion mobility K0 for Ne-CO2-N2 (90-10-5) with H2O content 130 ppm. Deisting thesis page 88 +// const float mTNull = 273.15f; // in K, zero temperature +// const float mTGas = 294.15f; // in K, TPC gas temperature corresponding to 21 degrees Celsius +// const float mVDriftIon = 1.2577968; // in cm / ms, ion drift velocity = K0 * TGas / TNull * 1 atm / Pmeasured * Ez0 * 1e-3. For Ne-CO2-N2 (90-10-5) with H2O content 130 ppm at Pmeasured = 1 atm. + +// TPC parameters +// const float mEz0 = 400.f; // in V / cm, nominal drift field in z direction +const float mZROC = o2::tpc::TPCParameters::TPCZ0; // absolute - position of G1T +const float mRIFC = o2::tpc::TPCParameters::IFCRADIUS; // inner field cage radius in cm +const float mROFC = o2::tpc::TPCParameters::OFCRADIUS; // outer field cage radius in cm + +const char* outfnameHists = "spaceChargeDensityHist"; // name of the output file for the histograms +const char* outfnameIDC = "spaceChargeDensityIDC"; // name of the output file for the IDCs +const char* hisSCRandomName = "hisSCRandom"; + +CalPad loadMap(std::string mapfile, std::string mapName); +void normalizeHistoQVEps0(TH3& histoIonsPhiRZ); + +template +void setOmegaTauT1T2(o2::tpc::SpaceCharge& sc); + +/// Create SC density histograms and IDC containers from simulated TPC hits +/// An interaction rate of 50 kHz is assumed. Therefore, the ion drift time determines the number of ion pile-up events. +/// \param ionDriftTime ion drift time in ms. The value determines the number of bins in z/time direction of the histograms and IDCs ( 1 bin / ms / side). +/// \param nEvIon number of ion pile-up events +/// \param debug debug info streaming level +/// \param sides set which sides will be simulated. sides=0: A- and C-Side, sides=1: A-Side only, sides=2: C-Side only +/// \param inputfolder folder to the directory where the gain and epsilon map is stored. If an average distortion map is used, the distortion.root file should also be located there +/// \param distortionType sets the type of the electron distortions: 0->no distortions of electrons are applied, 1->average distortion of electrons. Distortions can be created by the makeDistortionsCorrections() function. +void createSCHistosFromHits(const int ionDriftTime = 200, const int nEvIon = 1, const int debug = 1, const int sides = 0, const char* inputfolder = "", const int distortionType = 0 /*, const int nThreads = 1*/) +{ + // load average distortions of electrons + SpaceCharge spacecharge; + if (distortionType == 1) { + const std::string inpFileDistortions = Form("%sdistortions.root", inputfolder); + TFile fInp(inpFileDistortions.data(), "READ"); + if (sides == 1) { // A-side + spacecharge.setGlobalDistortionsFromFile(fInp, Side::A); + } else if (sides == 2) { // C-side + spacecharge.setGlobalDistortionsFromFile(fInp, Side::C); + } else { + spacecharge.setGlobalDistortionsFromFile(fInp, Side::A); + spacecharge.setGlobalDistortionsFromFile(fInp, Side::C); + } + } + + auto startTotal = std::chrono::high_resolution_clock::now(); + gRandom->SetSeed(0); + std::cout << "Seed is: " << gRandom->GetSeed() << std::endl; + + auto& cdb = CDBInterface::instance(); + cdb.setUseDefaults(); + + const static Mapper& mapper = Mapper::instance(); + + static GEMAmplification& gemAmplification = GEMAmplification::instance(); + gemAmplification.updateParameters(); + + static ElectronTransport& electronTransport = ElectronTransport::instance(); + electronTransport.updateParameters(); + + static SAMPAProcessing& sampaProcessing = SAMPAProcessing::instance(); + sampaProcessing.updateParameters(); + + // load gain map + const std::string gainMapFile = Form("%sGainMap.root", inputfolder); + const std::string gainMapName = "GainMap"; + const CalPad gainMap = loadMap(gainMapFile, gainMapName); + + // load ibf map + const std::string ibfMapFile = Form("%sIBFMap.root", inputfolder); + const std::string ibfMapName = "Gain"; + const CalPad ibfMap = loadMap(ibfMapFile, ibfMapName); + + const std::string hitFileList = Form("%so2sim_HitsTPC.list", inputfolder); + TObjArray* arrHitFiles = (TObjArray*)gSystem->GetFromPipe(TString::Format("cat %s", hitFileList.data()).Data()).Tokenize("\n"); + const int nHitFiles = arrHitFiles->GetEntries(); // number of files with TPC hits + std::cout << "Number of Hit Files: " << nHitFiles << std::endl; + + // histograms for space-charge density + const int nPhiBins = 360; + const int nRBins = 257; + const int nZBins = 257 * 2; + + TH3F hisSCRandom(hisSCRandomName, hisSCRandomName, nPhiBins, 0, ::TWOPI, nRBins, mRIFC, mROFC, nZBins, -mZROC, mZROC); + TH3F hisIBF("hisIBF", "hisIBF", nPhiBins, 0, ::TWOPI, nRBins, mRIFC, mROFC, nZBins, -mZROC, mZROC); + TH3F hisPI("hisPI", "hisPI", nPhiBins, 0, ::TWOPI, nRBins, mRIFC, mROFC, nZBins, -mZROC, mZROC); + TH1F hisEpsilon("hisEpsilon", "hisEpsilon", 50, 0, 50); + + const int nZBinsSide = ionDriftTime; + // vector with CalDet objects for IDCs (1 / ms) + std::vector vecIDC(nZBinsSide); + for (auto& calpadIDC : vecIDC) { + calpadIDC = CalPad("IDC", PadSubset::ROC); + } + + for (int iev = 0; iev < nEvIon; ++iev) { + std::cout << " " << std::endl; + std::cout << "event: " << iev + 1 << " from " << nEvIon << " events" << std::endl; + const int indexHitFile = gRandom->Uniform(0, nHitFiles - 1); + std::cout << "indexHitFile: " << indexHitFile << std::endl; + const std::string hitFileName = arrHitFiles->At(indexHitFile)->GetName(); + std::cout << "hitFileName: " << hitFileName.data() << std::endl; + + std::unique_ptr hitFile = std::unique_ptr(TFile::Open(hitFileName.data())); + std::unique_ptr hitTree = std::unique_ptr((TTree*)hitFile->Get("o2sim")); + const int nEvents = hitTree->GetEntries(); // number of simulated events per hit file + + // vector of HitGroups per sector + std::vector<::HitGroup>* arrSectors[::Sector::MAXSECTOR]; + for (int isec = 0; isec < ::Sector::MAXSECTOR; ++isec) { + arrSectors[isec] = nullptr; + std::stringstream sectornamestr; + sectornamestr << "TPCHitsShiftedSector" << isec; + hitTree->SetBranchAddress(sectornamestr.str().c_str(), &arrSectors[isec]); + } + + // 1) Choose n events randomly from list of o2sim_HitsTPC.root files + const int indexEv = gRandom->Uniform(0, nEvents - 1); + std::cout << "index of event which is chosen: " << indexEv << std::endl; + hitTree->GetEntry(indexEv); + + // randomize ion drift length for position of ibf disk and dirft of primary ionization + const float driftLIons = gRandom->Uniform(0., mZROC); + + // rotate event by random phi angle + const float phiRot = gRandom->Uniform(0, ::TWOPI); + + // z position of ions + float zIonsIBF = mZROC - driftLIons; + const int zbinIDC = static_cast(zIonsIBF / mZROC * ionDriftTime); + + // in case driftLIons = 0 avoid seg fault + if (zbinIDC == nZBinsSide) { + continue; + } + + int startSec = 0; + int endSec = ::Sector::MAXSECTOR; + + // set sector loop depending on the side which was set + if (sides == 1) { // A-side + endSec *= 0.5; + } else if (sides == 2) { // C-side + startSec = endSec * 0.5; + } + + // #pragma omp parallel for num_threads(nThreads) // comment in for using multi threading + for (int isec = startSec; isec < endSec; ++isec) { // loop over sectors + auto vecTracks = arrSectors[isec]; + for (auto& hitsTrack : *vecTracks) { + for (size_t ihit = 0; ihit < hitsTrack.getSize(); ++ihit) { + const auto& elHit = hitsTrack.getHit(ihit); + GlobalPosition3D posHit(elHit.GetX(), elHit.GetY(), elHit.GetZ()); + const int nEle = static_cast(elHit.GetEnergyLoss()); + if (nEle <= 0) { + continue; + } + + // phi rotation of hit + const float rHit = posHit.rho(); + float phiHit = posHit.phi() + phiRot; + o2::math_utils::bringTo02PiGen(phiHit); + posHit.SetX(rHit * std::cos(phiHit)); + posHit.SetY(rHit * std::sin(phiHit)); + + // z position of ions + float zIonsPI = std::abs(posHit.Z()) - driftLIons; + if (posHit.Z() < 0) { + zIonsIBF *= -1; + zIonsPI *= -1; + } + + // Primary ionization + if (std::signbit(zIonsPI) == std::signbit(posHit.Z())) { + hisSCRandom.Fill(phiHit, rHit, zIonsPI, nEle); + hisPI.Fill(phiHit, rHit, zIonsPI, nEle); + } + + // apply distortion of electron if specified + if (distortionType == 1) { + spacecharge.distortElectron(posHit); + } + + // IBF: Place r-phi projection of hits randomly in z + float driftTimeEle = 0.f; + int epsilonSum = 0; + for (int iele = 0; iele < nEle; iele++) { + const GlobalPosition3D posHitDiff = electronTransport.getElectronDrift(posHit, driftTimeEle); + float phiHitDiff = posHitDiff.phi(); + o2::math_utils::bringTo02PiGen(phiHitDiff); + + const DigitPos digiPadPos = mapper.findDigitPosFromGlobalPosition(posHitDiff); + if (!digiPadPos.isValid()) { + continue; + } + + // Attachment + if (electronTransport.isElectronAttachment(driftTimeEle)) { + continue; + } + + const auto padPos = digiPadPos.getPadPos(); + const auto row = static_cast(padPos.getRow()); + const auto pad = static_cast(padPos.getPad()); + const CRU cru = digiPadPos.getCRU(); + + const int gain = static_cast(gemAmplification.getEffectiveStackAmplification() * gainMap.getValue(cru, row, pad)); + if (gain == 0) { + continue; + } + const int epsilon = static_cast(gain * ibfMap.getValue(cru, row, pad) * 0.01); // IBF value is in % -> convert to absolute value + epsilonSum += epsilon; + + hisSCRandom.Fill(phiHitDiff, posHitDiff.rho(), zIonsIBF, epsilon); + hisIBF.Fill(phiHitDiff, posHitDiff.rho(), zIonsIBF, epsilon); + + // fill pads with adc value + auto padPosGlobal = digiPadPos.getGlobalPadPos(); + auto rowRoc = static_cast(padPosGlobal.getRow()); + if (cru.roc().isOROC()) { + rowRoc -= mapper.getNumberOfRowsROC(ROC(0)); + } + + const float charge = vecIDC[zbinIDC].getValue(cru, row, pad) + gain; + ((CalPadArr&)(vecIDC[zbinIDC].getCalArray(static_cast(cru.roc().getRoc())))).setValue(rowRoc, pad, charge); + } // electron loop + if (nEle > 0) { + hisEpsilon.Fill(epsilonSum / nEle); + } + } // hit loop + } // track loop + } // sector loop + + for (int isec = 0; isec < ::Sector::MAXSECTOR; ++isec) { + delete arrSectors[isec]; + } + } // event loop + + // normalize histograms to Q / cm^3 / epsilon0 and IDC to Q + normalizeHistoQVEps0(hisSCRandom); + normalizeHistoQVEps0(hisIBF); + normalizeHistoQVEps0(hisPI); + for (auto& calpadIDC : vecIDC) { + calpadIDC *= TMath::Qe(); + } + + auto stopTotal = std::chrono::high_resolution_clock::now(); + std::chrono::duration elapsedTotal = stopTotal - startTotal; + printf("Total time: %f sec for %d events\n", elapsedTotal.count(), nEvIon); + + o2::utils::TreeStreamRedirector pcstream(Form("%s.root", outfnameHists), "recreate"); + pcstream.GetFile()->cd(); + if (debug > 0) { + printf("Dumping space-charge density to debug tree...\n"); + for (int iphi = 1; iphi <= hisSCRandom.GetNbinsX(); ++iphi) { + float phi = hisSCRandom.GetXaxis()->GetBinCenter(iphi); + + for (int ir = 1; ir <= hisSCRandom.GetNbinsY(); ++ir) { + float r = hisSCRandom.GetYaxis()->GetBinCenter(ir); + + for (int iz = 1; iz <= hisSCRandom.GetNbinsZ(); ++iz) { + float z = hisSCRandom.GetZaxis()->GetBinCenter(iz); + + float density = hisSCRandom.GetBinContent(iphi, ir, iz); + float densityIBF = hisIBF.GetBinContent(iphi, ir, iz); + float densityPI = hisPI.GetBinContent(iphi, ir, iz); + + pcstream << "density" + << "iphi=" << iphi + << "ir=" << ir + << "iz=" << iz + << "phi=" << phi + << "r=" << r + << "z=" << z + << "scdensity=" << density + << "scibf=" << densityIBF + << "scpi=" << densityPI + << "\n"; + } + } + } + } + hisSCRandom.Write(); + hisIBF.Write(); + hisPI.Write(); + hisEpsilon.Write(); + pcstream.Close(); + + // write idcs in different files than the histogram to be able to use "hadd" for merging + // make 1D-IDCs and 0D-IDCs + std::vector idc1DASide(nZBinsSide); // 1D-merged idc + std::vector idc1DCSide(nZBinsSide); // 1D-merged idc + + std::vector idc0DASide(1); // 0D-merged idc (A single float value cannot be written to file) TODO find better solution + std::vector idc0DCSide(1); // 0D-merged idc (A single float value cannot be written to file) TODO find better solution + + // calculate 1D IDC + for (unsigned long iSlice = 0; iSlice < vecIDC.size(); ++iSlice) { + const auto vecCalArr = vecIDC[iSlice].getData(); + const int maxrocs = ROC::MaxROC; + for (int iROC = 0; iROC < maxrocs; ++iROC) { + ROC roc(iROC); + const Side side = roc.side(); + if (side == Side::A) { + // 1D IDC for A side + idc1DASide[iSlice] += vecCalArr[iROC].getSum(); + } else { + // 1D IDC for C side + idc1DCSide[iSlice] += vecCalArr[iROC].getSum(); + } + } + } + + // calculate 0D IDC + idc0DASide[0] = std::accumulate(idc1DASide.begin(), idc1DASide.end(), (float)0); + idc0DCSide[0] = std::accumulate(idc1DCSide.begin(), idc1DCSide.end(), (float)0); + + std::cout << "output path is: " << outfnameIDC << std::endl; + TFile fIDC(Form("%s.root", outfnameIDC), "RECREATE"); + fIDC.WriteObject(&vecIDC, "IDC"); + fIDC.WriteObject(&idc1DASide, "IDC_1D_A_Side"); + fIDC.WriteObject(&idc1DCSide, "IDC_1D_C_Side"); + fIDC.WriteObject(&idc0DASide, "IDC_0D_A_Side"); + fIDC.WriteObject(&idc0DCSide, "IDC_0D_C_Side"); +} + +/// Normalize histogram with number of ions per bin to charge per bin volume +/// \param histoIonsPhiRZ TH3 histogram. Content: number of ions. Axes: phi, r, z +void normalizeHistoQVEps0(TH3& histoIonsPhiRZ) +{ + for (int iphi = 1; iphi <= histoIonsPhiRZ.GetNbinsX(); ++iphi) { + const float deltaPhi = histoIonsPhiRZ.GetXaxis()->GetBinWidth(iphi); + + for (int ir = 1; ir <= histoIonsPhiRZ.GetNbinsY(); ++ir) { + const float r0 = histoIonsPhiRZ.GetYaxis()->GetBinLowEdge(ir); + const float r1 = histoIonsPhiRZ.GetYaxis()->GetBinUpEdge(ir); + + for (int iz = 1; iz <= histoIonsPhiRZ.GetNbinsZ(); ++iz) { + const float deltaZ = histoIonsPhiRZ.GetZaxis()->GetBinWidth(iz); + const float volume = deltaPhi * 0.5 * (r1 * r1 - r0 * r0) * deltaZ; + + const float charge = histoIonsPhiRZ.GetBinContent(iphi, ir, iz) * TMath::Qe(); + histoIonsPhiRZ.SetBinContent(iphi, ir, iz, charge / (volume * mEpsilon0)); + } + } + } +} + +/// load gain or ibf map +/// \param mapfile file to the map +/// \param mapName name of the object +CalPad loadMap(std::string mapfile, std::string mapName) +{ + TFile f(mapfile.data(), "READ"); + CalPad* map = nullptr; + f.GetObject(mapName.data(), map); + + if (!map) { + std::cout << mapfile.data() << " NOT FOUND! RETURNING! " << std::endl; + } + + f.Close(); + return *map; +} + +/// merge the idc values from different root files +/// \param files vetor of paths to files containing the IDCs which will be averaged +/// \param outFile output filename +void makeAverageIDCs(const std::vector& files, const char* outFile = outfnameIDC) +{ + // vector containing the path of the relevant files which will be merged + std::cout << "merge IDCs for average map" << std::endl; + + // merged idc CalPads + std::vector idc3D; // 3D-merged idc + std::vector idc1DASide; // 1D-merged idc + std::vector idc1DCSide; // 1D-merged idc + + std::vector idc0DASide(1); // 0D-merged idc (A single float value cannot be written to file) TODO find better solution + std::vector idc0DCSide(1); // 0D-merged idc (A single float value cannot be written to file) TODO find better solution + + // merge the 3D IDC values + const int nMaps = files.size(); + for (int iFile = 0; iFile < nMaps; ++iFile) { + const auto str = files[iFile]; + std::cout << "merging file: " << str.data() << std::endl; + + TFile finp(str.data(), "READ"); + std::vector* idcTmp = nullptr; + finp.GetObject("IDC", idcTmp); + + // number of z-bins + const auto nSlices = idcTmp->size(); + + // init vector for first file + if (iFile == 0) { + idc3D.resize(nSlices); + for (auto& calpadIDC : idc3D) { + calpadIDC = CalPad("IDC", PadSubset::ROC); + } + idc1DASide.resize(nSlices); + idc1DCSide.resize(nSlices); + } + + // merge the idc values + for (unsigned long iSlice = 0; iSlice < nSlices; ++iSlice) { + idc3D[iSlice] += (*idcTmp)[iSlice]; + } + delete idcTmp; + } + + // if the average IDCs are calculated the IDCs have to be normalized to the number of maps + const int nSlices = idc3D.size(); + + // sum up all z-slices for A- and C-side + CalPad idcZSlicesSummedTmp("IDC", PadSubset::ROC); + for (int iSlice = 0; iSlice < nSlices; ++iSlice) { + idc3D[iSlice] /= nMaps; + idcZSlicesSummedTmp += idc3D[iSlice]; + } + // normalize to number of z slices + idcZSlicesSummedTmp /= nSlices; + + // setting each z-slice to average IDC + for (int iSlice = 0; iSlice < nSlices; ++iSlice) { + idc3D[iSlice] = idcZSlicesSummedTmp; + } + + // calculate 1D IDC + for (unsigned long iSlice = 0; iSlice < idc3D.size(); ++iSlice) { + const auto vecCalArr = idc3D[iSlice].getData(); + const int maxrocs = ROC::MaxROC; + for (int iROC = 0; iROC < maxrocs; ++iROC) { + ROC roc(iROC); + const Side side = roc.side(); + if (side == Side::A) { + // 1D IDC for A side + idc1DASide[iSlice] += vecCalArr[iROC].getSum(); + } else { + // 1D IDC for C side + idc1DCSide[iSlice] += vecCalArr[iROC].getSum(); + } + } + } + + // calculate 0D IDC + idc0DASide[0] = std::accumulate(idc1DASide.begin(), idc1DASide.end(), (float)0); + idc0DCSide[0] = std::accumulate(idc1DCSide.begin(), idc1DCSide.end(), (float)0); + + std::cout << "output path is: " << outFile << std::endl; + TFile fMergedIDC(outFile, "RECREATE"); + fMergedIDC.WriteObject(&idc3D, "IDC"); + fMergedIDC.WriteObject(&idc1DASide, "IDC_1D_A_Side"); + fMergedIDC.WriteObject(&idc1DCSide, "IDC_1D_C_Side"); + fMergedIDC.WriteObject(&idc0DASide, "IDC_0D_A_Side"); + fMergedIDC.WriteObject(&idc0DCSide, "IDC_0D_C_Side"); +} + +/// \param outFileDistortions path to the file for output distortions (which can be read in from a SpaceCharge object) +/// \param sides set for which sides the distortions/corrections will be calculated. sides=0: A- and C-Side, sides=1: A-Side only, sides=2: C-Side only +/// \param inpFile name of the root file containing the space charge density histogram +/// \param histName name of the space charge density histogram in the root file +template +void makeDistortionsCorrections(const char* outFileDistortions = "distortions.root", const int sides = 0, const char* inpFile = "", const char* histName = hisSCRandomName) +{ + TFile fSCDensity(inpFile, "READ"); + std::cout << "input file: " << inpFile << std::endl; + std::cout << "output file: " << outFileDistortions << std::endl; + + using SC = o2::tpc::SpaceCharge; + SC spacecharge; + setOmegaTauT1T2(spacecharge); + spacecharge.fillChargeDensityFromFile(fSCDensity, histName); + + if (sides != 2) { + spacecharge.calculateDistortionsCorrections(Side::A); + } + if (sides != 1) { + spacecharge.calculateDistortionsCorrections(Side::C); + } + + int sideStart = 0; + int sideEnd = 2; + if (sides == 1) { + sideEnd = 1; + } else if (sides == 2) { + sideStart = 1; + } + + // dump distortion object to file if output file is specified + TFile fOut(outFileDistortions, "RECREATE"); + for (int iSide = sideStart; iSide < sideEnd; ++iSide) { + const Side side = iSide == 0 ? Side::A : Side::C; + spacecharge.dumpGlobalCorrections(fOut, side); + spacecharge.dumpGlobalDistortions(fOut, side); + spacecharge.dumpLocalCorrections(fOut, side); + spacecharge.dumpDensity(fOut, side); + } +} + +/// make average distortion map from random maps +/// \param files vector with files which contain the random space charge maps +/// \param histName name of the space charge histogram in the root files +void makeAverageDensityMap(const std::vector files, const char* histName = hisSCRandomName) +{ + // 1. loop over the maps and create the average map (still z dependent) + TH3F averageMap; + const int nMaps = files.size(); + for (int iFile = 0; iFile < nMaps; ++iFile) { + const auto str = files[iFile]; + std::cout << "using density map: " << str.data() << std::endl; + TFile fInp(str.data(), "READ"); + TH3F* densMapTmp = (TH3F*)fInp.Get(histName); + if (iFile == 0) { + averageMap = *densMapTmp; + } else { + averageMap.Add(densMapTmp); + } + delete densMapTmp; + } + averageMap.Scale(1. / nMaps); + + // 2.a sum up all z-slices to remove z dependence + const int nBinsPhi = averageMap.GetNbinsX(); + const int nBinsR = averageMap.GetNbinsY(); + const int nBinsZ = averageMap.GetNbinsZ(); + + for (int iPhi = 1; iPhi <= nBinsPhi; ++iPhi) { + for (int iR = 1; iR <= nBinsR; ++iR) { + const int nBinsHalfZ = 0.5 * nBinsZ; + + // C-Side (?) + const int startBinCSide = 1; + const int endBinCSide = nBinsHalfZ; + const float meanDensCSide = averageMap.Integral(iPhi, iPhi, iR, iR, startBinCSide, endBinCSide) / nBinsHalfZ; // integral over all z bins for each r and phi bin normalized to number of z slices + for (int iZ = startBinCSide; iZ <= endBinCSide; ++iZ) { + averageMap.SetBinContent(iPhi, iR, iZ, meanDensCSide); + } + + // A-Side (?) + const int startBinASide = 0.5 * nBinsZ + 1; + const int endBinASide = nBinsZ; + const float meanDensASide = averageMap.Integral(iPhi, iPhi, iR, iR, startBinASide, endBinASide) / nBinsHalfZ; // integral over all z bins for each r and phi bin normalized to number of z slices + for (int iZ = startBinASide; iZ <= endBinASide; ++iZ) { + averageMap.SetBinContent(iPhi, iR, iZ, meanDensASide); + } + } + } + + TFile fOut("spaceChargeDensityHist_average.root", "RECREATE"); + averageMap.Write(); + fOut.Close(); +} + +/// use the average space charge density map, scale it and calculat the corrections +/// \param inpFile input density file +/// \param outFile output scaled density file +/// \param sides set for which sides will be processed. sides=0: A- and C-Side, sides=1: A-Side only, sides=2: C-Side only +/// \param scaleFac multiply sigma by this value. The resulting scaling is "1 + scaleFac * sigmaScale" +/// \param sigmaScale sigma of the scaling +template +void createScaledMeanMap(const std::string inpFile, const std::string outFile, const int sides, const int scaleFac = 1, const float sigmaScale = 0.03f) +{ + // load the mean histo + using SC = o2::tpc::SpaceCharge; + SC scOriginal; + SC scScaled; + + TFile fInp(inpFile.data(), "READ"); + if (sides != 2) { + scOriginal.setDensityFromFile(fInp, Side::A); + } + if (sides != 1) { + scOriginal.setDensityFromFile(fInp, Side::C); + } + + setOmegaTauT1T2(scScaled); + int sideStart = 0; + int sideEnd = 2; + if (sides == 1) { + sideEnd = 1; + } else if (sides == 2) { + sideStart = 1; + } + + for (int iSide = sideStart; iSide < sideEnd; ++iSide) { + const Side side = iSide == 0 ? Side::A : Side::C; + for (size_t iZ = 0; iZ < nZ; ++iZ) { + for (size_t iR = 0; iR < nR; ++iR) { + for (size_t iPhi = 0; iPhi < nPhi; ++iPhi) { + const DataT density = scOriginal.getDensity(iZ, iR, iPhi, side); + const float scaleVal = 1 + scaleFac * sigmaScale; + scScaled.fillDensity(density * scaleVal, iZ, iR, iPhi, side); + } + } + } + scScaled.setDensityFilled(side); + } + + if (sides != 2) { + scScaled.calculateDistortionsCorrections(Side::A); + } + if (sides != 1) { + scScaled.calculateDistortionsCorrections(Side::C); + } + + // dump distortion object to file if output file is specified + for (int iSide = sideStart; iSide < sideEnd; ++iSide) { + const Side side = iSide == 0 ? Side::A : Side::C; + TFile fOut(outFile.data(), "RECREATE"); + scScaled.dumpGlobalCorrections(fOut, side); + scScaled.dumpGlobalDistortions(fOut, side); + scScaled.dumpLocalCorrections(fOut, side); + scScaled.dumpDensity(fOut, side); + } +} + +/// scale the IDCs from the average (input) map +/// \param inpIDCs input IDC File +/// \param outFile output file name +/// \param scaleFac multiply sigma by this value. The resulting scaling is "1 + scaleFac * sigmaScale" +/// \param sigmaScale sigma of the scaling +void scaleIDCs(const char* inpIDCs, const char* outFile, const int scaleFac = 1, const float sigmaScale = 0.03f) +{ + const float scaleVal = 1 + scaleFac * sigmaScale; + std::cout << "scaling IDC map: " << inpIDCs << std::endl; + + TFile finp(inpIDCs, "READ"); + std::vector* idc3D = nullptr; + finp.GetObject("IDC", idc3D); + + std::vector* idc1DASide = nullptr; + std::vector* idc1DCSide = nullptr; + finp.GetObject("IDC_1D_A_Side", idc1DASide); + finp.GetObject("IDC_1D_C_Side", idc1DCSide); + + std::vector* idc0DASide = nullptr; + std::vector* idc0DCSide = nullptr; + finp.GetObject("IDC_0D_A_Side", idc0DASide); + finp.GetObject("IDC_0D_C_Side", idc0DCSide); + + // scale the 3d idcs + for (unsigned long iSlice = 0; iSlice < idc3D->size(); ++iSlice) { + (*idc3D)[iSlice] *= scaleVal; + } + + // scale the 1d idcs + for (unsigned long iSlice = 0; iSlice < idc1DASide->size(); ++iSlice) { + (*idc1DASide)[iSlice] *= scaleVal; + (*idc1DCSide)[iSlice] *= scaleVal; + } + + // scale the 0d idcs + for (unsigned long iSlice = 0; iSlice < idc0DASide->size(); ++iSlice) { + (*idc0DASide)[iSlice] *= scaleVal; + (*idc0DCSide)[iSlice] *= scaleVal; + } + + std::cout << "output path is: " << outFile << std::endl; + TFile fMergedIDC(outFile, "RECREATE"); + fMergedIDC.WriteObject(idc3D, "IDC"); + fMergedIDC.WriteObject(idc1DASide, "IDC_1D_A_Side"); + fMergedIDC.WriteObject(idc1DCSide, "IDC_1D_C_Side"); + fMergedIDC.WriteObject(idc0DASide, "IDC_0D_A_Side"); + fMergedIDC.WriteObject(idc0DCSide, "IDC_0D_C_Side"); + + delete idc3D; + delete idc1DASide; + delete idc1DCSide; + delete idc0DASide; + delete idc0DCSide; +} + +/// helper function to set omegatau for the space charge class +template +void setOmegaTauT1T2(o2::tpc::SpaceCharge& sc) +{ + sc.setOmegaTauT1T2(0.32f, 1, 1); +} From 77096b0c31e1340be7a1165c01b19dbe6fdeebc1 Mon Sep 17 00:00:00 2001 From: Matthias Kleiner Date: Sat, 23 Jan 2021 10:01:42 +0100 Subject: [PATCH 2/2] setting macro to compile only --- Detectors/TPC/spacecharge/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Detectors/TPC/spacecharge/CMakeLists.txt b/Detectors/TPC/spacecharge/CMakeLists.txt index df3b60aa764ff..3f529cdb82a3c 100644 --- a/Detectors/TPC/spacecharge/CMakeLists.txt +++ b/Detectors/TPC/spacecharge/CMakeLists.txt @@ -48,7 +48,7 @@ o2_add_test_root_macro(macro/createSCHistosFromHits.C O2::TPCSimulation O2::DataFormatsTPC O2::MathUtils - LABELS tpc) + LABELS tpc COMPILE_ONLY) o2_add_test(PoissonSolver COMPONENT_NAME spacecharge