DrumCalib.cc

Go to the documentation of this file.

#include <TFile.h>
#include <TChain.h>
#include <TTree.h>
#include <TGraph.h>

#include <cmath>
#include <iomanip>
#include <cstdlib>
#include <iostream>
#include <string>
#include <map>
#include <vector>
#include <sstream>
#include <fstream>
using namespace std;

struct Parameters {
  string fileNames;
} gParameters;

struct RunInfo{
  RunInfo() : nTotalPhotons(0) {}
  string offlineRev;
  string configList;
  string md5;
  unsigned int nTotalPhotons;
};

struct CalibResult {
  double calib;
  double calibErr;
  double nPhotGen;
  double nPhotPixel;
  double electronicsFactor;
};

map<int, RunInfo> ReadRunInfo(const string& fileNames);
void CalibrateTel(ostream& calibXml, const int index, const RunInfo runInfo);
map<int, CalibResult> CalibrateTelescope(TChain* data);



// ------------------------------------------------------------
int main(int argc, char** argv) {

  if (argc<2) {
    cout << " Please run: DrumCalib <input-files>\n" << endl;
    return 5;
  }
  gParameters.fileNames = string(argv[1]);
  cout << "Reading from: \"" << gParameters.fileNames << "\"" <<endl;
  map<int, RunInfo> runs = ReadRunInfo(gParameters.fileNames);


  ofstream calibXml("FSimulationCalibConfig.generated.xml.in");
  calibXml << "<?xml version=\"1.0\" encoding=\"iso-8859-1\"?>\n\n"
           << "<FSimulationCalibConfig xmlns:xsi=\"http://www.w3.org/2001/XMLSchema-instance\"\n"
           << "                        xsi:noNamespaceSchemaLocation='@XMLSCHEMALOCATION@/FSimulationCalibConfig.xsd'>\n\n"
           << "  <!-- these are used for the flatfielding -->\n\n";
  
  
  for (map<int, RunInfo>::const_iterator iRun = runs.begin();
       iRun != runs.end(); ++iRun) {
    
    const int index = iRun->first;
    const RunInfo run = iRun->second;


    calibXml << "  <!-- \n"
             << "    This set of calibration constants was generated using Offline revision "
             << run.offlineRev << "\n"
             << "    with the following set of configuration parameters:\n";
    calibXml << "\"" << run.configList << "\"\n"
             << "    The following Md5-hex-digest encodes this list:\n"
             << "    -->\n";
    calibXml << "  <EndToEnd signature=\"" << run.md5 << "\"> \n"
             << "    <calibration>" << endl;
    

    CalibrateTel(calibXml, index, run);    
  }

  calibXml << "</FSimulationCalibConfig>\n";
  calibXml.close();
  
  cout << "-----------------------------------------------------------------------------" << endl;
  cout << "Telescope(s) calibrated! Output in file \"FSimulationCalibConfig.generated.xml.in\"." << endl;
  
  return 0;
}


// ------------------------------------------------------------
map<int, RunInfo> ReadRunInfo(const string& fileNames) {
  
  TChain* runInfo = new TChain("runInfo");
  runInfo->Add(fileNames.c_str());
  
  int index;
  unsigned int nPhotTotal;
  char md5[1000];
  char offline_rev[100];
  char config[1000000];
  runInfo->SetBranchAddress("index", &index);
  runInfo->SetBranchAddress("nPhotonsTot", &nPhotTotal);
  runInfo->SetBranchAddress("md5", md5);
  runInfo->SetBranchAddress("offline_rev", offline_rev);
  runInfo->SetBranchAddress("config", config);
  
  map<int, RunInfo> run;
  for (int i=0; i<runInfo->GetEntries(); ++i) {
    
    runInfo->GetEntry(i);
    
    if (run.count(index)) {
      run[index].nTotalPhotons += nPhotTotal;
      if (run[index].offlineRev != string(offline_rev) ||
          run[index].configList != config || 
          run[index].md5 != md5) {
        cout << " ERROR: Input data mismatch ! " << endl;
        exit(1);
      }
      
    } else {
      
      run[index].nTotalPhotons = nPhotTotal;
      run[index].offlineRev = offline_rev;
      run[index].configList = config;
      run[index].md5 = md5;
      
    }
        
  }
  
  return run;
}



// -----------------------------------------------------------------------
void CalibrateTel(ostream& calibXml, const int index, const RunInfo runInfo) {

  const bool writeErrorsToXml = false;
  
  cout << "Calibrate telescope, index = " << index << endl;
  
  const int eyeId = index / 100;
  const int telId = index % 100;
  
  ostringstream dataName;
  dataName << "drumData_" << index;
  
  TChain* drumData = new TChain(dataName.str().c_str());
  drumData->Add(gParameters.fileNames.c_str());
  
  map<int, CalibResult> constants = CalibrateTelescope(drumData);

  /*
  // find largst pixel id
  int maxPixelId = 0;
  for (map<int, CalibResult>::const_iterator iCal = constants.begin();
       iCal != constants.end(); ++iCal) {
    if (iCal->first>maxPixelId) 
      maxPixelId = iCal->first;
  }
  vector<double> 
  */

  double meanError = 0;
  for (map<int, CalibResult>::const_iterator iCal = constants.begin();
       iCal != constants.end(); ++iCal) {

    const int pixelId = iCal->first;
    const CalibResult calData = iCal->second;
    
    const double calib = calData.calib;
    const double calibErr = calData.calibErr;
            

    cout << " pixel id=" << setw(4) << pixelId
         << " s=" << setw(10) << fixed << calib 
         << " +- " << setw(10) << calibErr  << " [ph/ADC]"
         << "\n";
    
    
    calibXml << "      " << setprecision(16) << setw(18) << calib;
    if (writeErrorsToXml)
      calibXml << "  <!-- +- " << setw(18) << calibErr << " -->";
    calibXml << "\n";
    
    meanError += calibErr;
    
  } // end loop pixels

  if (constants.size())
    meanError /= constants.size();
  
  calibXml << "  <!-- mean calib. uncertainty: " << setw(18) << meanError << " -->";
  
  
  calibXml << "    </calibration>\n"
              "  </EndToEnd>\n\n";
  
  
} // end of CalculateCalibrationConstants



// ---------------------------------------------------------------------------------
map<int, CalibResult> CalibrateTelescope(TChain* drumData) { 

  int pixelId;
  double nGen, nPix, calib, calibErr, electronicsFactor;
  drumData->SetBranchAddress("id", &pixelId);
  drumData->SetBranchAddress("nGen", &nGen);
  drumData->SetBranchAddress("nPix", &nPix);
  drumData->SetBranchAddress("calib", &calib);
  drumData->SetBranchAddress("calibErr", &calibErr);
  drumData->SetBranchAddress("electronicsFactor", &electronicsFactor);
  
  
  map<int, int> counter;
  map<int, CalibResult> constants;
  
  // count photons
  for (int i=0; i<drumData->GetEntries(); ++i) {

    drumData->GetEntry(i);
    
    if (counter.count(pixelId)) {
      counter[pixelId]++;
      constants[pixelId].nPhotGen += nGen;
      constants[pixelId].nPhotPixel += nPix;
      if (constants[pixelId].electronicsFactor != electronicsFactor) {
        cout << "ERROR: incompatible pixel data found!" << endl;
        exit(2);
      }
    } else {
      counter[pixelId] = 1;
      constants[pixelId].nPhotGen = nGen;
      constants[pixelId].nPhotPixel = nPix;
      constants[pixelId].electronicsFactor = electronicsFactor;
      constants[pixelId].calib = 0;
      constants[pixelId].calibErr = 0;
    }
  }

  // check pixel data
  bool first = true;
  double fact = 0;
  for (map<int, int>::const_iterator iPix = counter.begin(); 
       iPix != counter.end(); ++iPix ) {
    if (first) {
      fact = iPix->second;
    } else {
      if (fact != iPix->second) {
        cout << "ERROR: incompatible number of pixels found in input data " << endl;
        exit(3);
      }
    }
  }
  
  // compute calibration constants
  for (map<int, CalibResult>::iterator iPix = constants.begin(); 
       iPix != constants.end(); ++iPix ) {
    
    CalibResult& dat = iPix->second;
    
    const double epsilon = double(dat.nPhotPixel) / dat.nPhotGen;
    const double epsilonErr = sqrt((epsilon - epsilon*epsilon) / dat.nPhotGen);
    
    // const double pixelSignal = pixelPhotons * qEff * gain; // add electronics [ph -> ADC]
    // const double calib = flux_drum_perp * diaArea * pixelSolid / pixelSignal;
    
    const double calib = 1 / (epsilon * dat.electronicsFactor);
    const double calibErr = epsilonErr / (epsilon * epsilon * dat.electronicsFactor);
    
    dat.calib = calib;
    dat.calibErr = calibErr;

    cout << "gen=" << dat.nPhotGen << " pix=" << dat.nPhotPixel << " eps=" << epsilon << endl;
  }    
  
  return constants;
}