source: Idarraga/mafalda/MAFTools/TimePixCalibrationHandler.cpp @ 294

/*
 * TimePixCalibrationHandler.cpp
 *
 *  Created on: March, 2012
 *      Author: John Idarraga <idarraga@cern.ch>
 */


#include "TimePixCalibrationHandler.h"

#include <iostream>
#include <fstream>
#include <math.h>
#include <stdlib.h>

using namespace std;

namespace MAFTools {

TimePixCalibrationHandler::TimePixCalibrationHandler(const char * af,
                const char * bf,
                const char * cf,
                const char * tf,
                int width, int height){

        m_width = width;
        m_height = height;

        ifstream afs(af, fstream::in);
        ifstream bfs(bf, fstream::in);
        ifstream cfs(cf, fstream::in);
        ifstream tfs(tf, fstream::in);

        ////////////////////////////////////
        // a
        m_a = new double[m_width*m_height];
        double temp = 0.;
        int i = 0;
        while ( afs.good() ) {
                afs >> temp;
                if(afs.eof()) break;
                m_a[i++] = temp;
        }
        afs.close();
        if(i != m_width*m_height){
                cout << "Calibration does not seem to be complete for paremeter a." << endl;
                cout << "Got " << i << " items.  " << m_width << "*" << m_height
                                << " were requested.  Giving up." << endl;
                exit(1);
        }

        ////////////////////////////////////
        // b
        m_b = new double[m_width*m_height];
        i = 0;
        while ( bfs.good() ) {
                bfs >> temp;
                if(bfs.eof()) break;
                m_b[i++] = temp;
        }
        bfs.close();
        if(i != m_width*m_height){
                cout << "Calibration does not seem to be complete for paremeter b." << endl;
                cout << "Got " << i << " items.  " << m_width << "*" << m_height
                                << " were requested.  Giving up." << endl;
                exit(1);
        }

        ////////////////////////////////////
        // c
        m_c = new double[m_width*m_height];
        i = 0;
        while ( cfs.good() ) {
                cfs >> temp;
                if(cfs.eof()) break;
                m_c[i++] = temp;
        }
        cfs.close();
        if(i != m_width*m_height){
                cout << "Calibration does not seem to be complete for paremeter c." << endl;
                cout << "Got " << i << " items.  " << m_width << "*" << m_height
                                << " were requested.  Giving up." << endl;
                exit(1);
        }

        ////////////////////////////////////
        // t
        m_t = new double[m_width*m_height];
        i = 0;
        while ( tfs.good() ) {
                tfs >> temp;
                if(tfs.eof()) break;
                m_t[i++] = temp;
        }
        tfs.close();
        if(i != m_width*m_height){
                cout << "Calibration does not seem to be complete for paremeter t." << endl;
                cout << "Got " << i << " items.  " << m_width << "*" << m_height
                                << " were requested.  Giving up." << endl;
                exit(1);
        }

}

TimePixCalibrationHandler::~TimePixCalibrationHandler(){
}

double TimePixCalibrationHandler::GetE(pair<int,int> pix, int tot){

        //TF1 * surr = GetSurrogateTF1(pix);

        int index = XYtoC(pix, m_width);

        //double par[] = { m_a[index], m_b[index], m_c[index], m_t[index] };

        // cuadratic solution
        double a = m_a[index];
        double b = m_b[index] - m_a[index]*m_t[index] - (double)tot;
        double c = -m_c[index] - m_t[index]*m_b[index];

        double sol1 = -b + sqrt(b*b - 4*a*c);
        sol1 /= 2*a;
        double sol2 = -b - sqrt(b*b - 4*a*c);
        sol2 /= 2*a;

        cout << "sol1 = " << sol1 << " , sol2 = " << sol2 << endl;

        if(sol2 <= 0) return sol1;
        else return sol2;

        return 0.;

        /*
        if(pix.first == 111 && pix.second == 70) {
                std::cout << "a = " << par[0] << std::endl;
                std::cout << "b = " << par[1] << std::endl;
                std::cout << "c = " << par[2] << std::endl;
                std::cout << "t = " << par[3] << std::endl;
        }*/

/*
        double evalTOT = 0.;
        double prev_e = 0.;
        double dist = 0.;
        double prev_dist = 10000;

        // search
        for(double e = 3.0; e <= 10000. ; e+=0.05) { // steps of 0.05 keV, start search

                //evalTOT = surr->Eval(e);
                evalTOT = SurrogateCalibFunction(e, par);
                if(evalTOT < 0) continue; // don't consider those values. Not physical.
                dist = TMath::Abs(evalTOT - tot);
                //if(pix.first == 111 && pix.second == 70) cout << "dist = " << dist << " | e = " << e << " | evalTOT = " << evalTOT << " | tot = " << tot << endl;
                if( dist > prev_dist ) break; // if I am going far, break.

                prev_dist = dist;
                prev_e = e; // keep the last energy
        }

        return prev_e;
        */
}

int TimePixCalibrationHandler::GetTOT (pair<int, int> pix, double E) {

        int index = XYtoC(pix, m_width);
        double par[] = { m_a[index], m_b[index], m_c[index], m_t[index] };

        return SurrogateCalibFunction(E, par);
        //return GetSurrogateTF1(pix)->Eval(E);
}

double TimePixCalibrationHandler::SurrogateCalibFunction(double x, double * par){
        return par[0]*x + par[1] - (par[2]/(x-par[3]));
}

#ifndef __MAFALDA_EMBEDDED
// Only if needed one can fetch a few functions as TF1 objects
// Not included for embedded applications
TF1 * TimePixCalibrationHandler::GetSurrogateTF1(pair<int, int> pix){

        // check first if the function has been defined.  Use the Set.
        m_surrogateItr = m_surrogateSet.find(pix);

        // If surrogate TF1 does not exist yet, build the function first
        if( m_surrogateItr == m_surrogateSet.end() ) {
                TString ftitle = "surrogate_";
                ftitle += pix.first; ftitle += "_"; ftitle += pix.second;
                // Usually the parameters come given in keV !!! WARNING !!!
                m_surrogateMap[pix] = new TF1(ftitle,"[0]*x + [1] - ([2]/(x-[3]))",0,1000); // 1 MeV
                // Populate parameters.  Look for right index first
                int index = XYtoC(pix, m_width);
                m_surrogateMap[pix]->SetParameters(m_a[index], m_b[index], m_c[index], m_t[index]);
        }

        return m_surrogateMap[pix];
}
#endif

// These two functions are already defined in MAFTools but for certain tests
// I want this class to be self-contained.
int TimePixCalibrationHandler::XYtoC(pair<int, int> position, int dimX) {
        return XYtoC(position.first, position.second, dimX);
}
int TimePixCalibrationHandler::XYtoC(int x, int y, int dimX){
        return y * dimX + x;
}

} // end namespace MAFTools