tpc/generator/generateTruthIbfGainMap.C

0001 #include <TFile.h>
0002 #include <TH2.h>
0003 #include <TH3.h>
0004 #include <TString.h>
0005
0006 #include <format>
0007 #include <iostream>
0008
0009 bool generateTruthIbfGainMap(const char* adcFile, const char *adcName, const char* ionFile, const char* ibfName, const char* primName, const char*  /*outputFile*/){
0010
0011  //load the adc-per-bin data from the specified file.
0012   TFile* adcInputFile = TFile::Open(adcFile, "READ");
0013   TH3* hAdc=nullptr;
0014   adcInputFile->GetObject(adcName, hAdc);
0015   if (hAdc == nullptr)  {
0016     std::cout << "ADC hist " <<  adcName << " or file "
0017           <<  adcFile << " not found!" << std::endl;
0018     return false;
0019   }
0020
0021   //load the ions-per-bin data from the specified file
0022   TFile* ibfInputFile = TFile::Open(ionFile, "READ");
0023   TH3* hIbf=nullptr;
0024   TH3* hPrimaries=nullptr;
0025   ibfInputFile->GetObject(ibfName,hIbf);
0026   ibfInputFile->GetObject(primName,hPrimaries);
0027   if (hIbf == nullptr) {
0028     std::cout << "IBF hist " << ibfName << " or file "
0029           << ionFile << " not found!" << std::endl;
0030     return false;
0031   }
0032   if (hPrimaries == nullptr) {
0033     std::cout << "IBF hist " <<  primName << " IN file "
0034           << ionFile << " not found!" << std::endl;
0035     return false;
0036   }
0037
0038 //  TFile *output=TFile::Open(outputFile,"RECREATE");
0039   //generate 2D histograms with the IBF, and IBF + Primaries, per side:
0040
0041   const int nSections=3;
0042   TH2 *hFlatTotal[nSections];
0043   TH2 *hFlatIbf[nSections]; //flat charge sums in north and south
0044   int zbins;
0045   int neg;
0046   int cm;
0047   int pos;
0048   zbins=hAdc->GetZaxis()->GetNbins();
0049   neg=1;
0050   cm=zbins/2;
0051   pos=zbins;
0052
0053   TString suffix[]={"","_negz","_posz"};
0054   int low[]={neg,neg,cm+1};
0055   int high[]={pos,cm,pos};
0056
0057   for (int i=0;i<nSections;i++){
0058     hPrimaries->GetZaxis()->SetRange(low[i],high[i]);
0059     hFlatTotal[i]=(TH2*)hPrimaries->Project3D(std::format("xy{}",suffix[i].Data()).c_str());
0060     hIbf->GetZaxis()->SetRange(low[i],high[i]);
0061     hFlatIbf[i]=(TH2*)hIbf->Project3D(std::format("xy{}",suffix[i].Data()).c_str());
0062     hFlatTotal[i]->Add(hFlatIbf[i]);
0063   }
0064
0065
0066   //sanity check that the bounds are the same
0067
0068   TAxis* ax[3] = {nullptr, nullptr, nullptr};
0069   ax[0] = hPrimaries->GetXaxis();
0070   ax[1] = hPrimaries->GetYaxis();
0071   ax[2] = hPrimaries->GetZaxis();
0072
0073   TAxis* ax2[3] = {nullptr, nullptr, nullptr};
0074   ax2[0] = hAdc->GetXaxis();
0075   ax2[1] = hAdc->GetYaxis();
0076   ax2[2] = hAdc->GetYaxis();
0077
0078
0079   int nbins[3];
0080   for (int i = 0; i < 3; i++)//note these are dimensions, not sections.
0081   {
0082     nbins[i] = ax[i]->GetNbins();  //number of bins, not counting under and overflow.
0083     if (nbins[i]!=ax2[i]->GetNbins()){
0084       std::cout << "Primaries and Adc bins are different in axis " << i
0085         << ". (" << nbins[i] << " vs "<< ax2[i]->GetNbins() << ").  Aborting unless in z." << std::endl;
0086       if (i!=2)
0087       {
0088     return false;
0089       }
0090       std::cout << "this is a z difference, which we can recover." << std::endl;
0091     }
0092   }
0093
0094
0095   //project into 2D per rphi, and divide Ions by Adc to get the Ion/Adc gain.
0096   TH2 * hIonGain[nSections];
0097   TH2 *hIbfGain[nSections];
0098   TH2 *hFlatAdc[nSections];
0099
0100   for (int i=0;i<nSections;i++){
0101     hIonGain[i]=static_cast<TH2*>(hFlatTotal[i]->Clone(std::format("hIonGain{}",suffix[i].Data()).c_str()));//with primaries
0102     hIbfGain[i]=static_cast<TH2*>(hFlatIbf[i]->Clone(std::format("hIbfGain{}",suffix[i].Data()).c_str()));//with only ibf
0103     hAdc->GetZaxis()->SetRange(low[i],high[i]);
0104     hFlatAdc[i]=static_cast<TH2*>(hAdc->Project3D(std::format("xy{}",suffix[i].Data()).c_str()));
0105
0106     hIonGain[i]->Divide(hFlatAdc[i]);
0107     hIbfGain[i]->Divide(hFlatAdc[i]);
0108
0109     hIonGain[i]->Write();
0110     hIbfGain[i]->Write();
0111   }
0112
0113   //
0114   return true;
0115 }