added charge scale factor in channel status (#206)

* added charge scale factor in channel status

* moved definition of chargeScale per James's comment

* addressing Tanner's formatting comments

* Move chargeScale into AddMCPhoton

* Switched from PMT ID to LCN, also added charge calibration constants into the output ntuple.

* Switched from PMT ID to LCN, also added charge calibration constants into the output ntuple.

* removed unused line

ratdb/CHANNEL_STATUS_DEFAULTS.ratdb

@@ -12,3 +12,10 @@
   "run_range": [0, 0],
   "default_value": 1,
 }
+
+{
+  "name": "charge_scale",
+  "index": "",
+  "run_range": [0, 0],
+  "default_value": 1.0,
+}

src/core/include/RAT/Gsim.hh

@@ -74,15 +74,13 @@ class Gsim : public Producer, G4UserRunAction, G4UserEventAction, G4UserTracking
   void Init();  // the real constructor
   void AddMCPhoton(DS::MCPMT *rat_mcpmt, const GLG4HitPhoton *photon, EventInfo *exinfo = NULL,
                    std::string process = "unknown");
-
   /* Storing optical creation track ID and step */
   void PhotonRecurse(std::vector<int> &PhotonIDs, int trackID, int &parentID, int &firstCreatedID);
   void SetOpticalPhotonIDs(std::string particle_type, int trackID, int parentID);
   std::vector<int> OpticalPhotonIDs;
 
   G4RunManager *theRunManager;
   GLG4DebugMessenger *theDebugMessenger;
-
   RAT::DS::PMTInfo *fPMTInfo;
   std::vector<RAT::PMTTime *> fPMTTime;      //< PMT transit time/delay calculator (indexed by modeltype)
   std::vector<RAT::PMTCharge *> fPMTCharge;  //< PMT single-pe charge calculator (indexed by modeltype)

src/core/src/Gsim.cc

@@ -438,7 +438,6 @@ void Gsim::MakeEvent(const G4Event *g4ev, DS::Root *ds) {
   ds->SetRunID(theRunManager->GetCurrentRun()->GetRunID());
   mc->SetID(g4ev->GetEventID());
   mc->SetUTC(exinfo->utc);
-
   // Vertex Info
   for (int ivert = 0; ivert < g4ev->GetNumberOfPrimaryVertex(); ivert++) {
     G4PrimaryVertex *pv = g4ev->GetPrimaryVertex(ivert);
@@ -555,7 +554,6 @@ void Gsim::MakeEvent(const G4Event *g4ev, DS::Root *ds) {
     rat_mcpmt->SetType(fPMTInfo->GetType(a_pmt->GetID()));
 
     numPE += a_pmt->GetEntries();
-
     /** Add "real" hits from actual simulated photons */
     for (int i = 0; i < a_pmt->GetEntries(); i++) {
       // Find the optical process responsible
@@ -573,10 +571,10 @@ void Gsim::MakeEvent(const G4Event *g4ev, DS::Root *ds) {
 
 void Gsim::AddMCPhoton(DS::MCPMT *rat_mcpmt, const GLG4HitPhoton *photon, EventInfo * /*exinfo*/, std::string process) {
   DS::MCPhoton *rat_mcphoton = rat_mcpmt->AddNewMCPhoton();
+  double chargeScale = DS::RunStore::GetCurrentRun()->GetChannelStatus()->GetChargeScaleByPMTID(rat_mcpmt->GetID());
   // Only real photons are added in Gsim, noise and afterpulsing handled in processors
   rat_mcphoton->SetDarkHit(false);
   rat_mcphoton->SetAfterPulse(false);
-
   rat_mcphoton->SetLambda(photon->GetWavelength());
 
   double x, y, z;
@@ -592,7 +590,9 @@ void Gsim::AddMCPhoton(DS::MCPMT *rat_mcpmt, const GLG4HitPhoton *photon, EventI
   rat_mcphoton->SetCreationTime(photon->GetCreationTime());
 
   rat_mcphoton->SetFrontEndTime(fPMTTime[fPMTInfo->GetModel(rat_mcpmt->GetID())]->PickTime(photon->GetTime()));
-  rat_mcphoton->SetCharge(fPMTCharge[fPMTInfo->GetModel(rat_mcpmt->GetID())]->PickCharge());
+  // Set the charge for the photoelectron, scaled by an optional calibration parameter chargeScale with a default value
+  // of one
+  rat_mcphoton->SetCharge(chargeScale * fPMTCharge[fPMTInfo->GetModel(rat_mcpmt->GetID())]->PickCharge());
   rat_mcphoton->SetCreatorProcess(process);
 }
 

src/ds/include/RAT/DS/ChannelStatus.hh

@@ -28,10 +28,11 @@ class ChannelStatus : public TObject {
   ChannelStatus() : TObject() {}
   virtual ~ChannelStatus() {}
 
-  virtual void AddChannel(int lcn, int is_online, double offset) {
+  virtual void AddChannel(int lcn, int is_online, double offset, double chargescale) {
     lcns.push_back(lcn);
     online.push_back(is_online);
     cable_offset.push_back(offset);
+    charge_scale.push_back(chargescale);
     lcn_to_index[lcn] = lcns.size() - 1;
   }
 
@@ -41,10 +42,13 @@ class ChannelStatus : public TObject {
   virtual double GetCableOffsetByChannel(int lcn) const { return cable_offset.at(lcn_to_index.at(lcn)); }
   virtual double GetCableOffsetByPMTID(int pmtid) const { return cable_offset.at(pmtid_to_index.at(pmtid)); }
 
+  virtual double GetChargeScaleByChannel(int lcn) const { return charge_scale.at(lcn_to_index.at(lcn)); }
+  virtual double GetChargeScaleByPMTID(int pmtid) const { return charge_scale.at(pmtid_to_index.at(pmtid)); }
+
   virtual void LinkPMT(int pmtid, int lcn) {
     // create entry with default values if none are specified
     if (lcn_to_index.find(lcn) == lcn_to_index.end()) {
-      AddChannel(lcn, default_is_online, default_offset);
+      AddChannel(lcn, default_is_online, default_offset, default_charge_scale);
     }
     pmtid_to_index[pmtid] = lcn_to_index[lcn];
   }
@@ -54,6 +58,8 @@ class ChannelStatus : public TObject {
     default_offset = lCableOffset->GetD("default_value");
     DBLinkPtr lChannelOnline = DB::Get()->GetLink("channel_online", index);
     default_is_online = lChannelOnline->GetD("default_value");
+    DBLinkPtr lChargeScale = DB::Get()->GetLink("charge_scale", index);
+    default_charge_scale = lChannelOnline->GetD("default_value");
     for (int pmtid = 0; pmtid < pmtinfo->GetPMTCount(); pmtid++) {
       int lcn = pmtinfo->GetChannelNumber(pmtid);
       LinkPMT(pmtid, lcn);
@@ -66,6 +72,14 @@ class ChannelStatus : public TObject {
     } catch (DBNotFoundError& e) {
       warn << "Cable offset table not found! Looking for table cable_offset[" << index << "]\n";
     }
+    // charge scale
+    try {
+      std::vector<int> lcns = get_lcns(lChargeScale);
+      std::vector<double> values = lChargeScale->GetDArray("value");
+      insert_values(lcns, values, &charge_scale);
+    } catch (DBNotFoundError& e) {
+      warn << "Charge Scale table not found! Looking for table charge_scale[" << index << "]\n";
+    }
     // dead channels
     try {
       std::vector<int> lcns = get_lcns(lChannelOnline);
@@ -119,8 +133,10 @@ class ChannelStatus : public TObject {
   std::vector<int> lcns;
   std::vector<int> online;
   std::vector<double> cable_offset;
+  std::vector<double> charge_scale;
   double default_offset;
   int default_is_online;
+  double default_charge_scale;
 };
 
 }  // namespace DS

src/io/include/RAT/OutNtupleProc.hh

@@ -81,6 +81,7 @@ class OutNtupleProc : public Processor {
   std::vector<int> pmtChannel;
   std::vector<bool> pmtIsOnline;
   std::vector<double> pmtCableOffset;
+  std::vector<double> pmtChargeScale;
   std::vector<double> pmtX;
   std::vector<double> pmtY;
   std::vector<double> pmtZ;

src/io/src/OutNtupleProc.cc

@@ -81,6 +81,7 @@ bool OutNtupleProc::OpenFile(std::string filename) {
   metaTree->Branch("pmtChannel", &pmtChannel);
   metaTree->Branch("pmtIsOnline", &pmtIsOnline);
   metaTree->Branch("pmtCableOffset", &pmtCableOffset);
+  metaTree->Branch("pmtChargeScale", &pmtChargeScale);
   metaTree->Branch("pmtX", &pmtX);
   metaTree->Branch("pmtY", &pmtY);
   metaTree->Branch("pmtZ", &pmtZ);
@@ -620,6 +621,7 @@ OutNtupleProc::~OutNtupleProc() {
       pmtChannel.push_back(channel);
       pmtIsOnline.push_back(ch_status->GetOnlineByPMTID(id));
       pmtCableOffset.push_back(ch_status->GetCableOffsetByPMTID(id));
+      pmtChargeScale.push_back(ch_status->GetChargeScaleByPMTID(id));
       pmtX.push_back(position.X());
       pmtY.push_back(position.Y());
       pmtZ.push_back(position.Z());