CALIBRATION/hms_cer_calib/calibration.h

+//////////////////////////////////////////////////////////
+// This class has been automatically generated on
+// Wed Aug 15 17:10:26 2018 by ROOT version 6.10/02
+// from TTree T/Hall A Analyzer Output DST
+// found on file: ../../ROOTfiles/coin_replay_production_3424_-1.root
+//////////////////////////////////////////////////////////
+
+#ifndef calibration_h
+#define calibration_h
+
+#include <TROOT.h>
+#include <TChain.h>
+#include <TFile.h>
+#include <TSelector.h>
+#include <TTreeReader.h>
+#include <TTreeReaderValue.h>
+#include <TTreeReaderArray.h>
+
+// Headers needed by this particular selector
+
+
+class calibration : public TSelector {
+public :
+   TTreeReader     fReader;  //!the tree reader
+   TTree          *fChain = 0;   //!pointer to the analyzed TTree or TChain
+   Bool_t          fFullShow;
+   
+
+   //Declare histograms
+   TH1F          **fPulseInt; 
+   TH1F         ***fPulseInt_quad;
+   TH1F           *fCut_everything;
+   TH1F           *fCut_electron;
+   TH1F           *fBeta_Cut;
+   TH1F           *fBeta_Full;
+   TH1F          **fTiming_Cut;
+   TH1F          **fTiming_Full;
+
+   // Readers to access the data (delete the ones you do not need).
+   TTreeReaderArray<Double_t> H_cer_goodAdcPulseInt    = {fReader, "H.cer.goodAdcPulseInt"};
+   TTreeReaderArray<Double_t> H_cer_goodAdcTdcDiffTime = {fReader, "H.cer.goodAdcTdcDiffTime"};
+   TTreeReaderValue<Double_t> H_cer_xAtCer             = {fReader, "H.cer.xAtCer"};
+   TTreeReaderValue<Double_t> H_cer_yAtCer             = {fReader, "H.cer.yAtCer"};
+
+   TTreeReaderArray<Double_t> H_tr_beta                = {fReader, "H.tr.beta"};
+   TTreeReaderValue<Int_t> Ndata_H_tr_beta             = {fReader, "Ndata.H.tr.beta"};
+
+   TTreeReaderValue<Double_t> H_cal_etotnorm           = {fReader, "H.cal.etotnorm"};
+
+
+   calibration(TTree * /*tree*/ =0) {fPulseInt=0,fPulseInt_quad=0,fCut_everything=0,fCut_electron=0,fBeta_Cut=0,fBeta_Full=0,fTiming_Cut=0,fTiming_Full=0,fFullShow=kFALSE;}
+   virtual ~calibration() { }
+   virtual Int_t   Version() const { return 2; }
+   virtual void    Begin(TTree *tree);
+   virtual void    SlaveBegin(TTree *tree);
+   virtual void    Init(TTree *tree);
+   virtual Bool_t  Notify();
+   virtual Bool_t  Process(Long64_t entry);
+   virtual Int_t   GetEntry(Long64_t entry, Int_t getall = 0) { return fChain ? fChain->GetTree()->GetEntry(entry, getall) : 0; }
+   virtual void    SetOption(const char *option) { fOption = option; }
+   virtual void    SetObject(TObject *obj) { fObject = obj; }
+   virtual void    SetInputList(TList *input) { fInput = input; }
+   virtual TList  *GetOutputList() const { return fOutput; }
+   virtual void    SlaveTerminate();
+   virtual void    Terminate();
+
+   ClassDef(calibration,0);
+
+};
+
+#endif
+
+#ifdef calibration_cxx
+void calibration::Init(TTree *tree)
+{
+   // The Init() function is called when the selector needs to initialize
+   // a new tree or chain. Typically here the reader is initialized.
+   // It is normally not necessary to make changes to the generated
+   // code, but the routine can be extended by the user if needed.
+   // Init() will be called many times when running on PROOF
+   // (once per file to be processed).
+
+   fReader.SetTree(tree);
+}
+
+Bool_t calibration::Notify()
+{
+   // The Notify() function is called when a new file is opened. This
+   // can be either for a new TTree in a TChain or when when a new TTree
+   // is started when using PROOF. It is normally not necessary to make changes
+   // to the generated code, but the routine can be extended by the
+   // user if needed. The return value is currently not used.
+
+   return kTRUE;
+}
+
+
+#endif // #ifdef calibration_cxx

CALIBRATION/hms_cer_calib/run_cal.C

+#include <TProof.h>
+#include <iostream>
+#include <fstream>
+#include <string>
+#include <stdio.h>
+
+void run_cal(Int_t RunNumber = 0, Int_t NumEvents = 0, Int_t coin = 0)
+{
+  if (RunNumber == 0) {
+    cout << "Enter a Run Number (-1 to exit): ";
+    cin >> RunNumber;
+    if (RunNumber <= 0) return;
+  }
+  if (NumEvents == 0) {
+    cout << "\nNumber of Events to analyze: ";
+    cin >> NumEvents;
+  }
+  if (coin == 0) {
+    cout << "\nIf this is a coincident run enter 1: ";
+    cin >> coin;
+  }
+
+  cin.ignore(numeric_limits<streamsize>::max(), '\n');
+
+  string calib_raw;
+  cout << "\nEnter options for calibration  (enter NA to skip): ";
+  getline(std::cin, calib_raw);
+  TString calib_option = calib_raw;
+
+  cout << "\n\n";
+
+  TChain ch("T");
+  if (coin == 1) ch.Add(Form("../../ROOTfiles/coin_replay_production_%d_%d.root", RunNumber, NumEvents));
+  else ch.Add(Form("../../ROOTfiles/hms_replay_production_all_%d_%d.root", RunNumber, NumEvents));
+  TProof *proof = TProof::Open("");
+  //proof->SetProgressDialog(0);  
+  ch.SetProof();
+
+  if (calib_option != "NA") {
+    //Start calibration process
+    ch.Process("calibration.C+",calib_option);
+
+    cout << "\n\nUpdate calibration constants with the better estimate (y/n)? ";
+      
+    TString user_input;
+    cin >> user_input;
+    if (user_input == "y") {
+      ifstream temp;
+      temp.open("calibration_temp.txt", ios::in);
+      if (temp.is_open()) {
+	rename("calibration_temp.txt", Form("../../PARAM/HMS/CER/hcer_calib_%d.param", RunNumber));
+      }
+
+      else cout << "Error opening calibration constants, may have to update constants manually!" << endl;
+	   
+    }
+
+    else {
+      remove("calibration_temp.txt"); 
+    }
+  }
+}

CALIBRATION/hms_hodo_calib/.gitignore

+*.pdf
+*.root
+*.so
+*.pcm

CALIBRATION/hms_hodo_calib/.rootrc

+Rint.History: .root_history

CALIBRATION/hms_hodo_calib/README.md

+# Hall C HMS Hodo Calibration Codes
+
+A set of ROOT scripts to determine the time offset parameters for the HMS hodoscope. To use these parameters a  parameter flag ptofusinginvadc is set to 0. The first script determines the time walk correction factor. The next script determines the effective propagation speeed in the paddle, the time difference between the positive and negative PMTs and then the relative time difference of all paddles compared to paddle 7 in plane S1X.
+
+
+This supercedes an older calibration code. The parameters from that code are
+used if a parameter flag ptofusinginvadc is set to 1. 
+
+The two codes have different parameters and it is possible to switch between the parameters of the two calibration codes using the  parameter flag ptofusinginvadc
+
+
+## Instructions
+
+1. Replay the data with ptofusinginvadc=0 and need to have T.hms.* and H.hod.* in the tree. 
+   Go to: PARAM/HMS/GEN/h_fadc_debug.param, and set hhodo_debug_adc = 1 --->  Use the correct hodo raw leafs variables 
+
+(In case you are calibrating HMS hodo using a coincidence run, then make sure to include T.coin.*)
+
+
+2. Determine the time walk correction parameters
+
+     a. Start "root -l" and then  .x timeWalkHistos.C+("current_dir/to/ROOT_filename.root", Run_Number, "hms") ---> If doing coincidence, then "hms"->"coin"
+
+     b. This creats the file: timeWalkHistos.root
+
+     c. Start "root -l" and then .x timeWalkCalib.C+
+
+     d. This creates the parameter file "../../PARAM/HMS/HODO/hhodo_TWcalib_runnumber.param"
+
+3.  Replay the data with ptofusinginvadc=0 and the new parameter files (the simplest is to copy phodo_TWcalib_runnumber.param to phodo_TWcalib.param).
+
+4. Determine the effective propagation speed in the paddle, the time difference between the positive and negative PMTs and then the relative time difference of all paddles compared to paddle 7 in plane S1X. The script
+puts cuts on H.cal.etracknorm, H.hgcer.npeSum and H.hod.betanotrack to select electrons. These cuts are hard coded as  etrknrm_low_cut = 0.7, npngcer_npeSum_low_cut = 0.7 , betanotrack_low_cut = 0.5 and betanotrack_hi_cut = 1.5. These may need to be modified. The event must have a track. 
+
+*** Make sure to have included H.hgcer.* and H.cal.etracknorm in the ROOT tree
+
+     a.  Start "root -l" and then  .x  fitHodoCalib.C+("current_dir/to/ROOT_filename.root",Run_Number)     
+
+     b.  This creates the parameter file "../../PARAM/HMS/HODO/hhodo_Vpcalib_runnumber.param"
+
+     c.  It also creates the root file HodoCalibPlots_runnumber.root
+
+     d.  To analyze cosmic data :  .x  fitHodoCalib.C+("current_dir/to/ROOT_filename.root",Run_Number,kTRUE) 
+
+     e. For cosmic data the speed of light is set to -30 cm/ns and the PID cut is just on P.hod.betanotrack with the default of betanotrack_low_cut = -1.2 and betanotrack_hi_cut = -.7

CALIBRATION/set_peddefault/README.md

+1) First need to prepare histograms of the goodAdcPed for each detector.
+2) the 2d histograms of goodADCPed verus the paddle number are contained
+   in HMS_PedDefault.def and SHMS_PedDefault.def
+   Move files to DEF-files/SHMS/PRODUCTION/SHMS_PedDefault.def
+   and DEF-files/HMS/PRODUCTION/HMS_PedDefault.def
+3) Include the files in your  tree/histogram def file set in replay script.
+   #include "DEF-files/SHMS/PRODUCTION/SHMS_PedDefault.def"
+   #include "DEF-files/HMS/PRODUCTION/HMS_PedDefault.def"
+4) Replay the data
+5) Start root
+   a) .L run_ped_default.C
+   b) run_shms_ped_default("entirepath/DirName/filename.root")
+   c) run_hms_ped_default("entirepath/DirName/filename.root")
+6) The SHMS does HGCER, NGCER, AERO, Preshower and Shower
+7) HSM does CER and CAL.
+8) By hand copy each set of Pedestal defaults into the detector "cuts" file
+  For example:
+phgcer_PedDefault= 2086, 2153, 2320, 1987
+
+  goes into phgcer_cuts.param or the equivelent
+

CALIBRATION/set_peddefault/run_ped_default.C

+#include "TFile.h"
+#include "TH1D.h"
+#include <iostream>
+void calc_ped_default(TString , TString ,TString, Double_t);
+void run_shms_ped_default(TString);
+void run_hms_ped_default(TString);
+
+//
+void run_shms_ped_default(TString file_name = "") {
+  calc_ped_default(file_name,"hgcer","p",0);
+  calc_ped_default(file_name,"ngcer","p",0);
+  calc_ped_default(file_name,"aero","p",1);
+  calc_ped_default(file_name,"aero","p",2);
+  calc_ped_default(file_name,"cal_prshwr","p",1);
+  calc_ped_default(file_name,"cal_prshwr","p",2);
+  calc_ped_default(file_name,"cal_shwr","p",0);
+}
+//
+void run_hms_ped_default(TString file_name = "") {
+  calc_ped_default(file_name,"cer","h",0);
+  calc_ped_default(file_name,"cal_hA","h",1);
+  calc_ped_default(file_name,"cal_hB","h",1);
+  calc_ped_default(file_name,"cal_hC","h",1);
+  calc_ped_default(file_name,"cal_hD","h",1);
+  calc_ped_default(file_name,"cal_hA","h",2);
+  calc_ped_default(file_name,"cal_hB","h",2);
+
+}
+
+
+void calc_ped_default(TString golden_file = "", TString detector = "",
+                           TString spect = "", Double_t polarity = -1) {
+
+  if (golden_file == "") {
+    cout << "Enter golden run root file name: " << endl;
+    cin >> golden_file;
+  }
+  if (detector == "") {
+    cout << "Enter detector prefix (hodo_1x (etc.), hgcer, aero, cal_prshwr, "
+            "cal_shwr, ngcer, cal_hA (etc.)): "
+         << endl;
+    cin >> detector;
+  }
+  if (spect == "") {
+    cout << "Enter a spectrometer (p or h): " << endl;
+    cin >> spect;
+  }
+  if (polarity == -1) {
+    cout << "Enter a detector polarity (pos = 1, neg = 2): " << endl;
+    cin >> polarity;
+  }
+
+  TString histname = Form("%s%s", spect.Data(), detector.Data());
+  if (histname.Contains("hodo") && histname.Contains("1x") && polarity == 1)
+    histname = Form("%s%s", histname.Data(), "_good_pped_vs_pmt_pos");
+  if (histname.Contains("hodo") && histname.Contains("1x") && polarity == 2)
+    histname = Form("%s%s", histname.Data(), "_good_pped_vs_pmt_neg");
+  if (histname.Contains("hodo") && histname.Contains("1y") && polarity == 1)
+    histname = Form("%s%s", histname.Data(), "_good_pped_vs_pmt_pos");
+  if (histname.Contains("hodo") && histname.Contains("1y") && polarity == 2)
+    histname = Form("%s%s", histname.Data(), "_good_pped_vs_pmt_neg");
+  if (histname.Contains("hodo") && histname.Contains("2x") && polarity == 1)
+    histname = Form("%s%s", histname.Data(), "_good_pped_vs_pmt_pos");
+  if (histname.Contains("hodo") && histname.Contains("2x") && polarity == 2)
+    histname = Form("%s%s", histname.Data(), "_good_pped_vs_pmt_neg");
+  if (histname.Contains("hodo") && histname.Contains("2y") && polarity == 1)
+    histname = Form("%s%s", histname.Data(), "_good_pped_vs_pmt_pos");
+  if (histname.Contains("hodo") && histname.Contains("2y") && polarity == 2)
+    histname = Form("%s%s", histname.Data(), "_good_pped_vs_pmt_neg");
+  if (histname.Contains("aero") && polarity == 1)
+    histname = Form("%s%s", histname.Data(), "_good_pped_vs_pmt_pos");
+  if (histname.Contains("aero") && polarity == 2)
+    histname = Form("%s%s", histname.Data(), "_good_pped_vs_pmt_neg");
+  if (histname.Contains("hcer"))
+    histname = Form("%s%s", histname.Data(), "_good_pped_vs_pmt");
+  if (histname.Contains("hgcer"))
+    histname = Form("%s%s", histname.Data(), "_good_pped_vs_pmt");
+  if (histname.Contains("ngcer"))
+    histname = Form("%s%s", histname.Data(), "_good_pped_vs_pmt");
+  if (histname.Contains("_prshwr") && polarity == 1)
+    histname = Form("%s%s", histname.Data(), "_good_pped_vs_pmt_pos");
+  if (histname.Contains("_prshwr") && polarity == 2)
+    histname = Form("%s%s", histname.Data(), "_good_pped_vs_pmt_neg");
+  if (histname.Contains("_shwr"))
+    histname = Form("%s%s", histname.Data(), "_good_pped_vs_pmt");
+  if (histname.Contains("hcal_h") && polarity == 1)
+    histname = Form("%s%s", histname.Data(), "_good_pped_vs_pmt_pos");
+  if (histname.Contains("hcal_h") && polarity == 2)
+    histname = Form("%s%s", histname.Data(), "_good_pped_vs_pmt_neg");
+
+  TH2F* H1_ped_vs_pmt;
+
+  TFile* f1 = new TFile(golden_file, "READ");
+  if (f1->IsZombie()) {
+    cout << "Cannot find : " << golden_file << endl;
+    return;
+  }
+
+  f1->GetObject(histname.Data(), H1_ped_vs_pmt);
+
+  TH1D* H1_pmt;
+
+  H1_pmt = (TH1D*)H1_ped_vs_pmt->ProjectionX("H1_pmt", 1,
+                                             H1_ped_vs_pmt->GetSize() - 2);
+
+  TH1D* H1_ped[H1_pmt->GetSize() - 2];
+
+  for (Int_t ipmt = 0; ipmt < (H1_pmt->GetSize() - 2); ipmt++) {
+    H1_ped[ipmt] = (TH1D*)H1_ped_vs_pmt->ProjectionY(
+        Form("H1_ped_pmt%d", ipmt + 1), ipmt + 1, ipmt + 1);
+  }
+
+  Double_t H1_ped_peak[H1_pmt->GetSize() - 2];
+  Double_t H1_ped_fitpeak[H1_pmt->GetSize() - 2];
+  for (Int_t ipmt = 0; ipmt < (H1_pmt->GetSize() - 2); ipmt++) {
+
+    if (H1_ped[ipmt]->GetEntries() > 25) {
+      TSpectrum* s = new TSpectrum(1);
+      gSystem->RedirectOutput("/dev/null", "a");
+      s->Search(H1_ped[ipmt], 1.0, "nobackground&&nodraw", 0.001);
+      gSystem->RedirectOutput(0);
+      TList* functions = H1_ped[ipmt]->GetListOfFunctions();
+      TPolyMarker* pm = (TPolyMarker*)functions->FindObject("TPolyMarker");
+      if (pm) {
+        H1_ped_peak[ipmt] = *pm->GetX();
+      } else {
+        H1_ped_peak[ipmt] = 0;
+      }
+    } else {
+      H1_ped_peak[ipmt] = 0;
+    }
+
+  }
+
+  TF1* Gaussian =
+      new TF1("Gaussian", "[0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2]))");
+  Gaussian->SetParLimits(0, 0, 1000);
+  Gaussian->SetParLimits(2, 0, 2);
+
+  for (Int_t ipmt = 0; ipmt < (H1_pmt->GetSize() - 2); ipmt++) {
+    if (H1_ped_peak[ipmt] != 0) {
+    Gaussian->SetParameter(1, H1_ped_peak[ipmt]);
+    gSystem->RedirectOutput("/dev/null", "a");
+    H1_ped[ipmt]->Fit(Gaussian, "QMN");
+    gSystem->RedirectOutput(0);
+    H1_ped_peak[ipmt]=int(Gaussian->GetParameter(1)*4*4.096);
+    }
+  }
+  TString side="";
+  if (polarity ==1)  side="Pos";
+  if (polarity ==2)  side="Neg";
+  TString PedDefName = spect+detector+"_Ped"+side+"Default= ";
+  if (detector == "cal_prshwr") {
+    PedDefName = spect+"cal_Ped"+side+"Default= ";
+  }
+  if (detector == "cal_shwr") {
+    PedDefName = spect+"cal_arr_Ped"+side+"Default= ";
+  }
+  if (detector == "cal_hA") {
+    PedDefName = spect+"cal_Ped"+side+"Default= ";
+  }
+  if (detector == "cal_hB" || detector == "cal_hC" || detector == "cal_hD") {
+  } else {
+  cout << PedDefName ;
+  }
+   for (Int_t ipmt = 0; ipmt < (H1_pmt->GetSize() - 2); ipmt++) {
+     cout << H1_ped_peak[ipmt] ;
+     if (ipmt <(H1_pmt->GetSize() - 3))  cout << ", " ;
+       if (ipmt!=0 && ipmt%16 == 0 )  cout << endl;
+   } 
+   cout << endl;
+   if (detector == "cal_hB" && polarity==2) {
+    for (Int_t ipmt = 0; ipmt < (H1_pmt->GetSize() - 2); ipmt++) {
+     cout << 0 ;
+     if (ipmt <(H1_pmt->GetSize() - 3))  cout << ", " ;
+       if (ipmt!=0 && ipmt%16 == 0 )  cout << endl;
+     } 
+   cout << endl;
+     for (Int_t ipmt = 0; ipmt < (H1_pmt->GetSize() - 2); ipmt++) {
+     cout << 0 ;
+     if (ipmt <(H1_pmt->GetSize() - 3))  cout << ", " ;
+       if (ipmt!=0 && ipmt%16 == 0 )  cout << endl;
+     } 
+    cout << endl;
+    
+   }
+
+}

CALIBRATION/shms_aero_calib/gain.r

+paero_neg_gain = 1./13.0732, 1./11.1519, 1./12.7171, 1./11.7946, 1./10.7654, 1./8.72566, 1./12.9509, 
+paero_pos_gain = 1./7.45058, 1./10.1029, 1./10.7762, 1./13.0162, 1./10.5675, 1./10.9487, 1./10.8312, 

CALIBRATION/shms_aero_calib/gain_4985.r

+paero_neg_gain = 1./13.5843, 1./11.1229, 1./12.5078, 1./11.5904, 1./11.019, 1./8.54602, 1./12.5383, 
+paero_pos_gain = 1./7.3637, 1./9.83655, 1./10.4292, 1./12.8879, 1./10.45, 1./10.5741, 1./10.6287, 

CALIBRATION/shms_cal_calib/THcPShTrack.h

@@ -166,6 +166,7 @@ void THcPShTrack::SetEs(Double_t* alpha) {
       UInt_t ncol = 1;
       if (nblk > fNrows_pr) ncol = 2;
       (*iter)->SetEdep(adc*Ycor(Y,ncol)*alpha[nblk-1]);
+      //(*iter)->SetEdep(adc*alpha[nblk-1]);
     }
     else
       //Shower block, no coordinate correction.

CALIBRATION/shms_cal_calib/input.dat

+-15. 27.	Delta range, %
+0.5 1.5	Beta range
+0.  	Heavy Gas Cherenkov, threshold on signals in p.e.
+10.  	Noble Gas Cherenkov, threshold on signals in p.e.
+5.	Minimum number of hits per channel required to be calibrated
+0. 3.	Range of uncalibrated energy deposition histogram
+500	Binning of uncalibrated energy deposition histogram
+-.1 1.25	Gaussian fit range around e- peak in the uncalibrated Edep histogram
+2.5	Gaussian fit range around mean e- peak in the uncalibrated Edep histogram 
+
+; Calibration constants for file shms_defocused.root (SHMS 9643, 9644, 9645, 9646, 9648, 9649, 9650, 9652 and 9654) 306065 events processed
+
+pcal_neg_gain_cor = 23.98, 26.50, 27.96, 17.52, 22.82, 16.46, 22.16, 29.95, 30.88, 43.38, 21.45, 22.89, 21.67, 55.90,
+pcal_pos_gain_cor = 11.72, 53.90, 23.71, 24.34, 28.44, 30.40, 23.79, 26.81, 28.26, 18.71, 26.45, 29.51, 22.84,  0.00,
+pcal_arr_gain_cor =  0.00,  0.00,  0.00,  8.94, 26.83, 22.93, 56.80, 38.60, 28.83, 62.52, 44.36,  0.00,  0.00,  0.00,  0.00,  0.00,
+                     0.00, 27.21, 29.68, 35.31, 27.88, 36.87, 33.00, 30.82, 35.13, 30.02, 34.52, 26.17,  0.00,  0.00,  0.00,  0.00,
+                     0.00, 45.97, 27.32, 37.10, 29.30, 31.68, 28.79, 31.26, 30.81, 26.81, 30.30, 33.84, 27.26,  0.00,  0.00,  0.00,
+                     0.00, 37.45, 27.89, 37.59, 25.33, 27.25, 34.88, 29.72, 32.25, 29.53, 31.04, 26.59, 25.45, 37.25, 18.66,  0.00,
+                     0.00, 39.42, 32.60, 29.71, 30.02, 30.66, 27.36, 31.91, 31.44, 27.97, 29.05, 32.15, 31.40, 39.06, 39.77,  0.00,
+                     0.00, 40.89, 28.40, 30.11, 31.55, 33.34, 26.43, 28.06, 32.35, 28.88, 32.93, 33.22, 36.02, 25.97, 33.14,  0.00,
+                     0.00, 33.96, 28.59, 29.46, 29.55, 34.78, 30.58, 39.43, 41.19, 40.14, 35.08, 34.17, 28.21, 28.97, 28.91,  0.00,
+                     18.83, 29.93, 28.87, 28.48, 27.35, 29.80, 32.57, 35.65, 43.36, 44.82, 49.44, 30.55, 33.73, 32.49, 37.63,  8.54,
+                     24.72, 29.79, 28.03, 28.83, 30.09, 30.23, 30.97, 32.67, 31.17, 34.81, 27.39, 30.92, 32.68, 32.03, 36.40,  0.00,
+                     18.14, 26.43, 28.73, 29.75, 27.81, 28.17, 28.83, 32.40, 30.13, 34.18, 31.61, 33.15, 32.58, 32.02, 29.99,  0.00,
+                     9.24, 25.62, 26.91, 24.97, 26.43, 31.29, 30.08, 27.68, 31.54, 33.71, 31.36, 32.03, 29.98, 33.83, 34.00,  0.00,
+                     0.00, 26.08, 23.04, 26.06, 29.22, 30.54, 28.60, 27.62, 26.95, 27.16, 27.36, 30.72, 28.00, 20.02,  0.00,  0.00,
+                     0.00, 31.38, 26.62, 32.41, 57.55, 28.71, 23.28, 26.16, 27.63, 22.19, 22.46, 36.94,  0.00,  0.00,  0.00,  0.00,
+                     0.00,  0.00,  0.00,  0.00, 42.28,  5.81, 45.36, 25.30, 21.53, 46.54, 16.55,  0.00,  0.00,  0.00,  0.00,  0.00,
+