packages/tpc/TpcLoadDistortionCorrection.cc

0001 /*!
0002  * \file TpcLoadDistortionCorrection.cc
0003  * \brief loads distortion correction histogram from file to DistortionCorrectionObject and stores on node tree
0004  * \author Hugo Pereira Da Costa <hugo.pereira-da-costa@cea.fr>
0005  */
0006
0007 #include "TpcLoadDistortionCorrection.h"
0008 #include "TpcDistortionCorrectionContainer.h"
0009
0010 #include <fun4all/Fun4AllReturnCodes.h>
0011 #include <phool/PHCompositeNode.h>
0012 #include <phool/PHDataNode.h>
0013 #include <phool/getClass.h>
0014
0015 #include <TFile.h>
0016 #include <TH1.h>
0017
0018 namespace
0019 {
0020
0021   // print histogram
0022   void print_histogram(TH1* h)
0023   {
0024     std::cout << "TpcLoadDistortionCorrection::InitRun - name: " << h->GetName() << std::endl;
0025     for (const auto& axis : {h->GetXaxis(), h->GetYaxis(), h->GetZaxis()})
0026     {
0027       if (axis)
0028       {
0029         std::cout
0030             << "  " << axis->GetName()
0031             << " bins: " << axis->GetNbins()
0032             << " min: " << axis->GetXmin()
0033             << " max: " << axis->GetXmax()
0034             << std::endl;
0035       }
0036     }
0037     std::cout << std::endl;
0038   }
0039
0040 }  // namespace
0041
0042 //_____________________________________________________________________
0043 TpcLoadDistortionCorrection::TpcLoadDistortionCorrection(const std::string& name)
0044   : SubsysReco(name)
0045 {
0046 }
0047
0048 //_____________________________________________________________________
0049 int TpcLoadDistortionCorrection::InitRun(PHCompositeNode* topNode)
0050 {
0051   // look for distortion calibration object
0052   PHNodeIterator iter(topNode);
0053
0054   std::cout << "TpcLoadDistortionCorrection::InitRun - m_flags: (i,in_use,radians,interpolate_z):";
0055   for (int i=0; i<4; i++)
0056   {
0057     std::cout << "("<< i <<", "<<m_correction_in_use[i] << ", " << m_phi_hist_in_radians[i] << ", " << m_interpolate_z[i] << ")"<< std::endl;
0058   }
0059
0060   /// Get the RUN node and check
0061   auto runNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "RUN"));
0062   if (!runNode)
0063   {
0064     std::cout << "TpcLoadDistortionCorrection::InitRun - RUN Node missing, quitting" << std::endl;
0065     return Fun4AllReturnCodes::ABORTRUN;
0066   }
0067
0068   // create and populate the nodes for each distortion, if present:
0069   for (int i = 0; i < 4; i++)
0070   {
0071     if (!m_correction_in_use[i])
0072     {
0073       continue;
0074     }
0075
0076     // get distortion correction object and create if not found
0077     auto distortion_correction_object = findNode::getClass<TpcDistortionCorrectionContainer>(topNode, m_node_name[i]);
0078     if (!distortion_correction_object)
0079     {
0080       std::cout << "TpcLoadDistortionCorrection::InitRun - creating TpcDistortionCorrectionContainer in node " << m_node_name[i] << std::endl;
0081       distortion_correction_object = new TpcDistortionCorrectionContainer;
0082       auto node = new PHDataNode<TpcDistortionCorrectionContainer>(distortion_correction_object, m_node_name[i]);
0083       runNode->addNode(node);
0084     }
0085
0086     std::cout << "TpcLoadDistortionCorrection::InitRun - reading corrections from " << m_correction_filename[i] << std::endl;
0087     auto distortion_tfile = TFile::Open(m_correction_filename[i].c_str());
0088     if (!distortion_tfile)
0089     {
0090       std::cout << "TpcLoadDistortionCorrection::InitRun - cannot open " << m_correction_filename[i] << std::endl;
0091       exit(1);
0092     }
0093
0094     const std::array<const std::string, 2> extension = {{"_negz", "_posz"}};
0095     for (int j = 0; j < 2; ++j)
0096     {
0097       distortion_correction_object->m_hDPint[j] = dynamic_cast<TH1*>(distortion_tfile->Get((std::string("hIntDistortionP")+extension[j]).c_str()));
0098       assert(distortion_correction_object->m_hDPint[j]);
0099       distortion_correction_object->m_hDRint[j] = dynamic_cast<TH1*>(distortion_tfile->Get((std::string("hIntDistortionR")+extension[j]).c_str()));
0100       assert(distortion_correction_object->m_hDRint[j]);
0101       distortion_correction_object->m_hDZint[j] = dynamic_cast<TH1*>(distortion_tfile->Get((std::string("hIntDistortionZ")+extension[j]).c_str()));
0102       assert(distortion_correction_object->m_hDZint[j]);
0103     }
0104
0105     // assign correction object dimension from histograms dimention, assuming all histograms have the same
0106     distortion_correction_object->m_dimensions = distortion_correction_object->m_hDPint[0]->GetDimension();
0107
0108     // only dimensions 2 or 3 are supported
0109     assert(distortion_correction_object->m_dimensions == 2 || distortion_correction_object->m_dimensions == 3);
0110
0111     // assign whether phi corrections (DP) should be read as radians or mm
0112     distortion_correction_object->m_phi_hist_in_radians = m_phi_hist_in_radians[i];
0113
0114     // assign whether 2D corrections should be interpolated to zero at readout or not (has no effect on 3D corrections)
0115     distortion_correction_object->m_interpolate_z = m_interpolate_z[i];
0116
0117     //assign whether the correction should be scaled, and if so, by how much
0118     distortion_correction_object->m_use_scalefactor = m_use_scalefactor[i];
0119     distortion_correction_object->m_scalefactor = m_scalefactor[i];
0120
0121
0122     if (Verbosity())
0123     {
0124       for (const auto& h : {
0125                distortion_correction_object->m_hDPint[0], distortion_correction_object->m_hDPint[1],
0126                distortion_correction_object->m_hDRint[0], distortion_correction_object->m_hDRint[1],
0127                distortion_correction_object->m_hDZint[0], distortion_correction_object->m_hDZint[1]})
0128       {
0129         print_histogram(h);
0130       }
0131     }
0132   }
0133   return Fun4AllReturnCodes::EVENT_OK;
0134 }
0135
0136 //_____________________________________________________________________
0137 int TpcLoadDistortionCorrection::process_event(PHCompositeNode* /*unused*/)
0138 {
0139   return Fun4AllReturnCodes::EVENT_OK;
0140 }