Added photon count approximation based on gain curves

Analysis/include/Containers.h

@@ -60,6 +60,8 @@ class Channel {
     std::string PMTcode;
     /** Polarity of the channel */
     bool pos_polarity;
+    /** Gain of the PMT for this channel */
+    float gain;
     /*
      *  Waveform storage
      */
@@ -121,6 +123,8 @@ class Channel {
     std::vector< double > Diff_max;
     /** Vector of integral of PWF_histo in each hit window in pC */
     std::vector< double > Charge;
+    /** Approximate number of photons detected by this channel */
+    std::vector< float > nPhotons;
     /** Pedestal of raw waveform */
     double PedMean;
     /** Pedestal RMS of raw waveform */

Analysis/src/DataReader2021.cpp

@@ -204,6 +204,8 @@ void DataReader2021::LoadConfigurationFile(std::string _inFile ){
     std::cout << "Loading .xml Configuration File..." << std::endl;
     std::cout << "Found " << m_XMLparser->getBaseNodeCount("channel") << " channel entries " << std::endl;
 
+    TFile gainFile( (std::getenv("JZCaPA") + std::string("/Utils/PMTgain.root")).c_str(), "read");
+
     std::vector < Channel* > channelEntries;
     int first_run, last_run;
 
@@ -228,6 +230,13 @@ void DataReader2021::LoadConfigurationFile(std::string _inFile ){
         m_XMLparser->getChildValue("channel",i,"is_on",buffer_ch->is_on);
         m_XMLparser->getChildValue("channel",i,"Vop",buffer_ch->Vop);
 
+        //Get the voltage vs gain plot for this PMT. If it exists, set the gain value for the channel
+        TGraph *gainPlot = (TGraph*)gainFile.Get(buffer_ch->PMTcode.c_str());
+        if(!gainPlot){
+          std::cout << Form("No gain plot found for %s row %d col %d", buffer_ch->detector.c_str(), buffer_ch->mapping_row, buffer_ch->mapping_column) << std::endl;
+        }
+        buffer_ch->gain = (gainPlot) ? gainPlot->Eval(buffer_ch->HV) : 0;
+
         bool isNew(true);
         for( int k = 0; k < channelEntries.size(); k++){
             if(buffer_ch->name == channelEntries.at(k)->name){

Analysis/src/RPDAnalysis2021.cpp

@@ -129,6 +129,7 @@ void RPDAnalysis2021::SetBranches( TTree* _tree ){
     for(int row = 0; row < 4; row++){
         for(int col = 0; col < 4; col++){
             m_AnalysisTree->Branch( Form("rpd%d_%d_Charge", row, col),           "std::vector<double>", &rpd[row][col]->Charge            );
+            m_AnalysisTree->Branch( Form("rpd%d_%d_nPhotons", row, col),         "std::vector<float>",  &rpd[row][col]->nPhotons          );
             m_AnalysisTree->Branch( Form("rpd%d_%d_Peak_max", row, col),         "std::vector<double>", &rpd[row][col]->Peak_max          );
             m_AnalysisTree->Branch( Form("rpd%d_%d_Diff_max", row, col),         "std::vector<double>", &rpd[row][col]->Diff_max          );
             m_AnalysisTree->Branch( Form("rpd%d_%d_Peak_center", row, col),      "std::vector<int>",    &rpd[row][col]->Peak_center       );

Analysis/src/WFAnalysis.cpp

@@ -145,6 +145,7 @@ void WFAnalysis::AnalyzeEvent( const std::vector< Channel* > vCh ){
     Ch->Peak_max.clear();
     Ch->Diff_max.clear();
     Ch->Charge.clear();
+    Ch->nPhotons.clear();
 
     //retrieving information for each channel as a histogram
     TH1D* h = Ch->WF_histo;
@@ -206,6 +207,7 @@ void WFAnalysis::AnalyzeEvent( const std::vector< Channel* > vCh ){
 
       for(int hit = 0; hit < Ch->nHits; hit++){
         Ch->Charge.push_back( GetCharge( Ch, Ch->hit_window.first[hit], Ch->hit_window.second[hit] ) );
+        Ch->nPhotons.push_back( -Ch->Charge.back()*6.241e6/Ch->gain ); //6.241e6 is electrons per pC
 
         //Get peak information from the waveform using the hit window
         Ch->Peak_center.push_back( GetMaximumBin( hProcessed, Ch->hit_window.first[hit], Ch->hit_window.second[hit] ) );