analysis/ZDCNeutronLocPol/ZDCNeutronLocPol.cc

0001 #include "ZDCNeutronLocPol.h"
0002
0003 #include <calobase/TowerInfoDefs.h>
0004 #include <caloreco/CaloWaveformFitting.h>
0005
0006 #include <calobase/TowerInfoDefs.h>
0007 #include <caloreco/CaloWaveformFitting.h>
0008 #include <ffarawobjects/CaloPacketContainerv1.h>
0009 #include <ffarawobjects/CaloPacketv1.h>
0010 #include <ffarawobjects/Gl1Packetv1.h>
0011 #include <fun4all/Fun4AllReturnCodes.h>
0012
0013 #include <Event/packet.h>
0014
0015 #include <fun4all/Fun4AllReturnCodes.h>
0016 #include <phool/PHCompositeNode.h>
0017 #include <phool/PHIODataNode.h>    // for PHIODataNode
0018 #include <phool/PHNodeIterator.h>  // for PHNodeIterator
0019 #include <phool/PHObject.h>        // for PHObject
0020 #include <phool/getClass.h>
0021 #include <phool/phool.h>
0022 #include <phool/recoConsts.h>
0023
0024 #include <TFile.h>
0025 #include <TH1.h>
0026 #include <TH2.h>
0027 #include <TTree.h>
0028 #include <phool/PHCompositeNode.h>
0029 #include <cmath>
0030 #include <fstream>
0031
0032 R__LOAD_LIBRARY(libuspin.so)
0033
0034 //____________________________________________________________________________..
0035 ZDCNeutronLocPol::ZDCNeutronLocPol(const std::string &name)
0036   : SubsysReco(name)
0037 {
0038   std::cout << "ZDCNeutronLocPol::ZDCNeutronLocPol(const std::string &name) Calling ctor" << std::endl;
0039 }
0040
0041 //____________________________________________________________________________..
0042 ZDCNeutronLocPol::~ZDCNeutronLocPol()
0043 {
0044   std::cout << "ZDCNeutronLocPol::~ZDCNeutronLocPol() Calling dtor" << std::endl;
0045 }
0046
0047 //____________________________________________________________________________..
0048 int ZDCNeutronLocPol::Init(PHCompositeNode * /*topNode*/)
0049 {
0050   std::cout << "ZDCNeutronLocPol::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0051
0052   smdHits = new TTree();
0053   smdHits = new TTree("smdHits", "smdHits");
0054   smdHits->SetDirectory(0);
0055   smdHits->Branch("bunchnumber", &bunchnumber, "bunchnumber/I");
0056   // smdHits -> Branch("evtseq_gl1",  &evtseq_gl1, "evtseq_gl1/I");
0057   // smdHits -> Branch("evtseq_zdc",  &evtseq_zdc, "evtseq_zdc/I");
0058   // smdHits -> Branch("BCO_gl1",  &BCO_gl1, "BCO_gl1/l");
0059   // smdHits -> Branch("BCO_zdc",  &BCO_zdc, "BCO_zdc/l");
0060   smdHits->Branch("showerCutN", &showerCutN, "showerCutN/I");
0061   smdHits->Branch("showerCutS", &showerCutS, "showerCutS/I");
0062   smdHits->Branch("n_x", &n_x, "n_x/F");
0063   smdHits->Branch("n_y", &n_y, "n_y/F");
0064   smdHits->Branch("s_x", &s_x, "s_x/F");
0065   smdHits->Branch("s_y", &s_y, "s_y/F");
0066   smdHits->Branch("zdcN1_adc", &zdcN1_adc, "zdcN1_adc/F");
0067   smdHits->Branch("zdcN2_adc", &zdcN2_adc, "zdcN2_adc/F");
0068   smdHits->Branch("zdcN3_adc", &zdcN3_adc, "zdcN3_adc/F");
0069   smdHits->Branch("zdcS1_adc", &zdcS1_adc, "zdcS1_adc/F");
0070   smdHits->Branch("zdcS2_adc", &zdcS2_adc, "zdcS2_adc/F");
0071   smdHits->Branch("zdcS3_adc", &zdcS3_adc, "zdcS3_adc/F");
0072   smdHits->Branch("veto_NF", &veto_NF, "veto_NF/F");
0073   smdHits->Branch("veto_NB", &veto_NB, "veto_NB/F");
0074   smdHits->Branch("veto_SF", &veto_SF, "veto_SF/F");
0075   smdHits->Branch("veto_SB", &veto_SB, "veto_SB/F");
0076   smdHits->Branch("smdS1_adc", &smdS1_adc, "smdS1_adc/F");
0077   smdHits->Branch("smdS2_adc", &smdS2_adc, "smdS2_adc/F");
0078   smdHits->Branch("smdS6_adc", &smdS6_adc, "smdS6_adc/F");
0079   smdHits->Branch("smdS7_adc", &smdS7_adc, "smdS7_adc/F");
0080
0081   smdHits->Branch("smdN1_adc", &smdN1_adc, "smdN1_adc/F");
0082   smdHits->Branch("smdN2_adc", &smdN2_adc, "smdN2_adc/F");
0083   smdHits->Branch("smdN6_adc", &smdN6_adc, "smdN6_adc/F");
0084   smdHits->Branch("smdN7_adc", &smdN7_adc, "smdN7_adc/F");
0085
0086   smdHits->Branch("smdS1_v_adc", &smdS1_v_adc, "smdS1_v_adc/F");
0087   smdHits->Branch("smdS2_v_adc", &smdS2_v_adc, "smdS2_v_adc/F");
0088   smdHits->Branch("smdS7_v_adc", &smdS7_v_adc, "smdS7_v_adc/F");
0089   smdHits->Branch("smdS8_v_adc", &smdS8_v_adc, "smdS8_v_adc/F");
0090
0091   smdHits->Branch("smdN1_v_adc", &smdN1_v_adc, "smdN1_v_adc/F");
0092   smdHits->Branch("smdN2_v_adc", &smdN2_v_adc, "smdN2_v_adc/F");
0093   smdHits->Branch("smdN7_v_adc", &smdN7_v_adc, "smdN7_v_adc/F");
0094   smdHits->Branch("smdN8_v_adc", &smdN8_v_adc, "smdN8_v_adc/F");
0095
0096   // waveform
0097   h_waveformZDC = new TH2F("h_waveformZDC", "h_waveformZDC", 17, -0.5, 16.5, 512, -500, 20000);
0098   h_waveformSMD_North = new TH2F("h_waveformSMD_North", "h_waveformSMD_North", 17, -0.5, 16.5, 512, -500, 20000);
0099   h_waveformSMD_South = new TH2F("h_waveformSMD_South", "h_waveformSMD_South", 17, -0.5, 16.5, 512, -500, 20000);
0100   h_waveformVeto_North = new TH2F("h_waveformVeto_North", "h_waveformVeto_North", 17, -0.5, 16.5, 512, -500, 20000);
0101   h_waveformVeto_South = new TH2F("h_waveformVeto_South", "h_waveformVeto_South", 17, -0.5, 16.5, 512, -500, 20000);
0102   h_waveformAll = new TH2F("h_waveformAll", "h_waveformAll", 17, -0.5, 16.5, 512, -500, 20000);
0103
0104   for (int i = 0; i < 32; i++)
0105   {
0106     h_rawADC[i] = new TH1F(Form("rawADC%d", i), Form("rawADC%d", i), 128, 0, 4000);
0107     h_pedADC[i] = new TH1F(Form("pedADC%d", i), Form("pedADC%d", i), 128, 0, 4000);
0108     h_signalADC[i] = new TH1F(Form("signalDC%d", i), Form("signalADC%d", i), 128, 0, 4000);
0109   }
0110
0111   return Fun4AllReturnCodes::EVENT_OK;
0112 }
0113
0114 //____________________________________________________________________________..
0115 int ZDCNeutronLocPol::InitRun(PHCompositeNode * /*topNode*/)
0116 {
0117   std::cout << "ZDCNeutronLocPol::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
0118
0119   evtcnt = 0;
0120
0121   return Fun4AllReturnCodes::EVENT_OK;
0122 }
0123
0124 //____________________________________________________________________________..
0125 int ZDCNeutronLocPol::process_event(PHCompositeNode *topNode)
0126 {
0127   // std::cout << "ZDCNeutronLocPol::process_event(PHCompositeNode *topNode) Processing Event" << std::endl;
0128   if (evtcnt % 10000 == 0)
0129   {
0130     std::cout << "Event: " << evtcnt << std::endl;
0131   }
0132
0133   p_gl1 = findNode::getClass<Gl1Packetv1>(topNode, "GL1Packet");
0134   zdc_cont = findNode::getClass<CaloPacketContainerv1>(topNode, "ZDCPackets");
0135
0136   showerCutN = 0;
0137   showerCutS = 0;
0138
0139   if (p_gl1)
0140   {
0141     bunchnumber = p_gl1->getBunchNumber();
0142     if (evtcnt % 1000 == 0)
0143     {
0144       std::cout << bunchnumber << std::endl;
0145     }
0146     if (zdc_cont->get_npackets() != 1)
0147     {
0148       return Fun4AllReturnCodes::EVENT_OK;
0149     }
0150
0151     CaloPacket *p_zdc = zdc_cont->getPacket(0);
0152
0153     if (p_zdc)
0154     {
0155       int nThresholdNVer = 0;
0156       int nThresholdNHor = 0;
0157       int nThresholdSVer = 0;
0158       int nThresholdSHor = 0;
0159
0160       // evtseq_gl1 = p_gl1->getEvtSequence();
0161       // evtseq_zdc = p_zdc->getEvtSequence();
0162       // BCO_gl1 = p_gl1->getBCO();
0163       // BCO_zdc = p_zdc->getBCO();
0164
0165       // in this for loop we get: zdc_adc and smd_adc
0166       for (int c = 0; c < p_zdc->iValue(0, "CHANNELS"); c++)
0167       {
0168         /*
0169         double baseline = 0.;
0170         double baseline_low = 0.;
0171         double baseline_high = 0.;
0172
0173         for(int s = 0; s < 3; s++)
0174         {
0175           baseline_low += p_zdc->iValue(s, c);
0176         }
0177
0178         baseline_low /= 3.;
0179
0180         for (int s = p_zdc->iValue(0, "SAMPLES")-3; s < p_zdc->iValue(0, "SAMPLES"); s++)
0181         {
0182           baseline_high += p_zdc->iValue(s,c);
0183         }
0184
0185         baseline_high /=3.;
0186
0187         baseline = baseline_low;
0188
0189         if(baseline_high < baseline_low) {baseline = baseline_high;}
0190
0191         for (int s = 0; s < p_zdc->iValue(0, "SAMPLES"); s++)
0192         {
0193           h_waveformAll->Fill(s, p_zdc->iValue(s, c) - baseline);
0194
0195           if (c < 16) //-->[0,15]
0196           {
0197             h_waveformZDC->Fill(s, p_zdc->iValue(s, c) - baseline);
0198           }
0199
0200           if ((c > 47) && (c < 64)) //-->[48,63]
0201           {
0202             h_waveformSMD_North->Fill(s, p_zdc->iValue(s, c) - baseline);
0203           }
0204
0205           if ((c > 111) && (c < 128)) //-->[112,127]
0206           {
0207             h_waveformSMD_South->Fill(s, p_zdc->iValue(s, c) - baseline);
0208           }
0209
0210           if ((c > 15) && (c < 18)) //-->[16,17]
0211           {
0212             h_waveformVeto_North->Fill(s, p_zdc->iValue(s, c) - baseline);
0213           }
0214
0215           if ((c > 79) && (c < 82)) //-->[80,81]
0216           {
0217             h_waveformVeto_South->Fill(s, p_zdc->iValue(s, c) - baseline);
0218           }
0219
0220         }
0221
0222         */
0223
0224         std::vector<float> resultFast = anaWaveformFast(p_zdc, c);  // fast waveform fitting
0225         float signalFast = resultFast.at(0);
0226         float timingFast = resultFast.at(1);
0227         float pedFast = resultFast.at(2);
0228
0229         ////////// ZDC channels ///////
0230         if (c >= 0 && c < 16)
0231         {
0232           signalFast = (4 < timingFast && timingFast < 10) ? signalFast : 0;
0233           unsigned int towerkey = TowerInfoDefs::decode_zdc(c);
0234           int zdc_side = TowerInfoDefs::get_zdc_side(towerkey);
0235           if (zdc_side == 0)
0236           {
0237             if (c == 0)
0238             {
0239               zdcS1_adc = signalFast;
0240               continue;
0241             }
0242             if (c == 2)
0243             {
0244               zdcS2_adc = signalFast;
0245               continue;
0246             }
0247             if (c == 4)
0248             {
0249               zdcS3_adc = signalFast;
0250               continue;
0251             }
0252           }
0253           else if (zdc_side == 1)
0254           {
0255             if (c == 8)
0256             {
0257               zdcN1_adc = signalFast;
0258               continue;
0259             }
0260             if (c == 10)
0261             {
0262               zdcN2_adc = signalFast;
0263               continue;
0264             }
0265             if (c == 12)
0266             {
0267               zdcN3_adc = signalFast;
0268               continue;
0269             }
0270           }
0271         }
0272         /////////////end ZDC channels//////////////////
0273
0274         //////////// North Veto counters //////////////////
0275         if (c == 16)
0276         {
0277           signalFast = (6 < timingFast && timingFast < 11) ? signalFast : 0;
0278           veto_NF = signalFast;
0279           continue;
0280         }
0281         if (c == 17)
0282         {
0283           signalFast = (6 < timingFast && timingFast < 11) ? signalFast : 0;
0284           veto_NB = signalFast;
0285           continue;
0286         }
0287         //////////// end North veto counters /////////
0288
0289         ///////// North SMD ////////
0290         if (c >= 48 && c < 63)
0291         {
0292           signalFast = (9 < timingFast) ? signalFast : 0;
0293
0294           rawADC[c - 48] = signalFast + pedFast;
0295           pedADC[c - 48] = pedFast;
0296           signalADC[c - 48] = signalFast;
0297
0298           h_rawADC[c - 48]->Fill(signalFast + pedFast);
0299           h_pedADC[c - 48]->Fill(pedFast);
0300           h_signalADC[c - 48]->Fill(signalFast);
0301           if (c == 48)
0302           {
0303             smdN1_v_adc = signalFast;
0304           }
0305           if (c == 49)
0306           {
0307             smdN2_v_adc = signalFast;
0308           }
0309           if (c == 54)
0310           {
0311             smdN7_v_adc = signalFast;
0312           }
0313           if (c == 55)
0314           {
0315             smdN8_v_adc = signalFast;
0316           }
0317           if (c == 56)
0318           {
0319             smdN1_adc = signalFast;
0320           }
0321           if (c == 57)
0322           {
0323             smdN2_adc = signalFast;
0324           }
0325           if (c == 61)
0326           {
0327             smdN6_adc = signalFast;
0328           }
0329           if (c == 62)
0330           {
0331             smdN7_adc = signalFast;
0332           }
0333           if (c < 56 && signalFast > 5.)
0334           {
0335             nThresholdNHor++;
0336           }
0337           if (c >= 56 && signalFast > 5.)
0338           {
0339             nThresholdNVer++;
0340           }
0341           continue;
0342         }
0343         //////// end North SMD///////////////
0344
0345         //////////// South Veto counters //////////////////
0346         if (c == 80)
0347         {
0348           signalFast = (6 < timingFast && timingFast < 11) ? signalFast : 0;
0349           veto_SF = signalFast;
0350           continue;
0351         }
0352         if (c == 81)
0353         {
0354           signalFast = (6 < timingFast && timingFast < 11) ? signalFast : 0;
0355           veto_SB = signalFast;
0356           continue;
0357         }
0358         //////////// end South veto counters /////////
0359         //// South SMD ////////
0360         if (c >= 112 && c < 127)
0361         {
0362           signalFast = (7 < timingFast && timingFast < 12) ? signalFast : 0;
0363           rawADC[c - 96] = signalFast + pedFast;
0364           pedADC[c - 96] = pedFast;
0365           signalADC[c - 96] = signalFast;
0366
0367           h_rawADC[c - 96]->Fill(signalFast + pedFast);
0368           h_pedADC[c - 96]->Fill(pedFast);
0369           h_signalADC[c - 96]->Fill(signalFast);
0370           if (c == 112)
0371           {
0372             smdS1_v_adc = signalFast;
0373           }
0374           if (c == 113)
0375           {
0376             smdS2_v_adc = signalFast;
0377           }
0378           if (c == 118)
0379           {
0380             smdS7_v_adc = signalFast;
0381           }
0382           if (c == 119)
0383           {
0384             smdS8_v_adc = signalFast;
0385           }
0386           if (c == 120)
0387           {
0388             smdS1_adc = signalFast;
0389           }
0390           if (c == 121)
0391           {
0392             smdS2_adc = signalFast;
0393           }
0394           if (c == 125)
0395           {
0396             smdS6_adc = signalFast;
0397           }
0398           if (c == 126)
0399           {
0400             smdS7_adc = signalFast;
0401           }
0402           if (c < 120 && signalFast > 5.)
0403           {
0404             nThresholdSHor++;
0405           }
0406           if (c >= 120 && signalFast > 5.)
0407           {
0408             nThresholdSVer++;
0409           }
0410           continue;
0411         }
0412         /////// end South SMD //////////////
0413
0414       }  // end channel loop
0415
0416       // shower cut (SMD noise threshold cut)
0417       if (nThresholdNHor > 1 && nThresholdNVer > 1)
0418       {
0419         showerCutN = 1;
0420       }
0421
0422       if (nThresholdSHor > 1 && nThresholdSVer > 1)
0423       {
0424         showerCutS = 1;
0425       }
0426
0427       ////////// Compute SMD position and fill hists
0428       CompSmdAdc();
0429       CompSmdPos();
0430       n_y = smd_pos[0];
0431       n_x = smd_pos[1];
0432       s_y = smd_pos[2];
0433
0434       s_x = smd_pos[3];
0435       smdHits->Fill();
0436     }  // end if p_zdc good
0437   }  // end if p_gl1 good
0438
0439   evtcnt++;
0440
0441   return Fun4AllReturnCodes::EVENT_OK;
0442 }
0443
0444 //____________________________________________________________________________..
0445 int ZDCNeutronLocPol::End(PHCompositeNode * /*topNode*/)
0446 {
0447   std::cout << "ZDCNeutronLocPol::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0448   TFile *ofile = new TFile(outfile.c_str(), "RECREATE");
0449
0450   smdHits->Write();
0451
0452   // h_waveformZDC->Write();
0453   // h_waveformSMD_North->Write();
0454   // h_waveformSMD_South->Write();
0455   // h_waveformVeto_North->Write();
0456   // h_waveformVeto_South->Write();
0457   // h_waveformAll->Write();
0458   /*
0459   for (int i = 0; i < 32; i++)
0460   {
0461     h_rawADC[i]->Write();
0462     h_pedADC[i]->Write();
0463     h_signalADC[i]->Write();
0464   }
0465   */
0466   ofile->Write();
0467   ofile->Close();
0468   delete (smdHits);
0469
0470   return Fun4AllReturnCodes::EVENT_OK;
0471 }
0472
0473 std::vector<float> ZDCNeutronLocPol::anaWaveformFast(CaloPacket *p, const int channel)
0474 {
0475   std::vector<float> waveform;
0476   // Chris: preallocation = speed improvement
0477   waveform.reserve(p->iValue(0, "SAMPLES"));
0478   for (int s = 0; s < p->iValue(0, "SAMPLES"); s++)
0479   {
0480     waveform.push_back(p->iValue(s, channel));
0481   }
0482   std::vector<std::vector<float>> multiple_wfs;
0483   multiple_wfs.push_back(waveform);
0484
0485   std::vector<std::vector<float>> fitresults_zdc;
0486   fitresults_zdc = WaveformProcessingFast->calo_processing_fast(multiple_wfs);
0487
0488   std::vector<float> result;
0489   result = fitresults_zdc.at(0);
0490   return result;
0491 }
0492
0493 void ZDCNeutronLocPol::CompSmdAdc()  // mulitplying by relative gains
0494 {
0495   for (int i = 0; i < 15; i++)  // last one is reserved for analog sum
0496   {
0497     // multiply SMD channels with their gain factors
0498     // to get the absolute ADC values in the same units
0499     // rgains come from CompSmdAdc()
0500     signalADC[i] = signalADC[i] * smd_north_rgain[i];
0501     signalADC[i + 16] = signalADC[i + 16] * smd_south_rgain[i];
0502   }
0503 }
0504
0505 void ZDCNeutronLocPol::CompSmdPos()  // computing position with weighted averages
0506 {
0507   float w_ave[4];  // 0 -> north hor (y); 1 -> noth vert (x); 2 -> south hor (y); 3 -> south vert. (x)
0508   float weights[4] = {0};
0509   memset(weights, 0, sizeof(weights));  // memset float works only for 0
0510   float w_sum[4];
0511   memset(w_sum, 0, sizeof(w_sum));
0512
0513   float smd_threshold = 5.;
0514
0515   // these constants convert the SMD channel number into real dimensions (cm's)
0516   const float hor_scale = 2.0 * 11.0 / 10.5 * sin(3.1415 / 4);  // from gsl_math.h
0517   const float ver_scale = 1.5 * 11.0 / 10.5;
0518   float hor_offset = (hor_scale * 8 / 2.0) * (7.0 / 8.0);
0519   float ver_offset = (ver_scale * 7 / 2.0) * (6.0 / 7.0);
0520
0521   for (int i = 0; i < 8; i++)
0522   {
0523     if (signalADC[i] > smd_threshold)
0524     {
0525       weights[0] += signalADC[i];            // summing weights
0526       w_sum[0] += (float) i * signalADC[i];  // summing for the average
0527     }
0528
0529     if (signalADC[i + 16] > smd_threshold)
0530     {
0531       weights[2] += signalADC[i + 16];
0532       w_sum[2] += ((float) i + 16.) * signalADC[i + 16];
0533     }
0534   }
0535   for (int i = 0; i < 7; i++)
0536   {
0537     if (signalADC[i + 8] > smd_threshold)
0538     {
0539       weights[1] += signalADC[i + 8];
0540       w_sum[1] += ((float) i + 8.) * signalADC[i + 8];
0541     }
0542
0543     if (signalADC[i + 24] > smd_threshold)
0544     {
0545       weights[3] += signalADC[i + 24];
0546       w_sum[3] += ((float) i + 24.) * signalADC[i + 24];
0547     }
0548   }
0549
0550   if (weights[0] > 0.0)
0551   {
0552     w_ave[0] = w_sum[0] / weights[0];  // average = sum / sumn of weights...
0553     smd_pos[0] = hor_scale * w_ave[0] - hor_offset;
0554   }
0555   else
0556   {
0557     smd_pos[0] = 0;
0558   }
0559
0560   if (weights[1] > 0.0)
0561   {
0562     w_ave[1] = w_sum[1] / weights[1];
0563     smd_pos[1] = ver_scale * (w_ave[1] - 8.0) - ver_offset;
0564   }
0565   else
0566   {
0567     smd_pos[1] = 0;
0568   }
0569
0570   if (weights[2] > 0.0)
0571   {
0572     w_ave[2] = w_sum[2] / weights[2];
0573     smd_pos[2] = hor_scale * (w_ave[2] - 16.0) - hor_offset;
0574   }
0575   else
0576   {
0577     smd_pos[2] = 0;
0578   }
0579
0580   if (weights[3] > 0.0)
0581   {
0582     w_ave[3] = w_sum[3] / weights[3];
0583     smd_pos[3] = ver_scale * (w_ave[3] - 24.0) - ver_offset;
0584   }
0585   else
0586   {
0587     smd_pos[3] = 0;
0588   }
0589 }
0590
0591 void ZDCNeutronLocPol::setFileName(const std::string &fname)
0592 {
0593   outfile = fname;
0594 }
0595
0596 void ZDCNeutronLocPol::setGainMatch(const std::string &gainfile)
0597 {
0598   // getting gains
0599   float col1;
0600   std::string col2;
0601    //std::string gainfile = gainname;
0602    std::ifstream gain_infile(gainfile);
0603
0604   if (!gain_infile)
0605   {
0606     std::cout << "No relative gain file: All relative gains set to 1." << std::endl;
0607     return;
0608   }
0609
0610   std::cout << "SMD Relative Gains: " << std::endl;
0611
0612   for (float &i : gain)
0613   {
0614     gain_infile >> col1 >> col2;
0615     i = col1;
0616   }
0617
0618   for (int i = 0; i < 16; i++)  // relative gains of SMD north channels
0619   {
0620     smd_north_rgain[i] = gain[i];  // 0-7: y channels, 8-14: x channels, 15: analog sum
0621     std::cout << "N" << i << " " << gain[i] << std::endl;
0622   }
0623
0624   for (int i = 0; i < 16; i++)  // relative gains of SMD south channels
0625   {
0626     smd_south_rgain[i] = gain[i + 16];  // 0-7: y channels, 8-14: x channels, 15: analog sum
0627     std::cout << "S" << i << " " << gain[i] << std::endl;
0628   }
0629
0630   gain_infile.close();
0631 }