tpc/generator/generateTruthIbfGainMap.C

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