source: JEM-EUSO/esaf_lal/tags/v1_r0/esaf/packages/simulation/detector/G4Detector/G4fresnellens/src/FresnelSurfaceHelper.cc @ 117

#include "FresnelSurfaceHelper.h"
#include "ArraySurface.h"
#include "SecondOrderSurface.h"
#include "ParametricSurface.h"
#include "FresnelSurface.h"

#include <iostream>
#include <fstream>
#include <string.h>
#include <float.h>

using namespace std;
using namespace G4FresnelLens;

//______________________________________________________________________________
FresnelSurfaceHelper::FresnelSurfaceHelper(const char* filename){
    Init();
    if (filename) Read(filename);
}

//______________________________________________________________________________
FresnelSurfaceHelper::FresnelSurfaceHelper(int n, double *rho, double *z){
    Init();
    Read(n, rho, z);
}

//______________________________________________________________________________
FresnelSurfaceHelper::~FresnelSurfaceHelper(){
    Clear();
}

//______________________________________________________________________________
void FresnelSurfaceHelper::Clear(){
    if ( fZ ) delete [] fZ;
    if ( fRho ) delete [] fRho;
    fZ=fRho=0;
    fN=0;

    if (fOwnSurfaces) {
        for (unsigned i(0); i<fSurfaces.size(); i++){
            delete fSurfaces[i];
        }
    }
    fSurfaces.clear();
    fOwnSurfaces=true;
}

//______________________________________________________________________________
void FresnelSurfaceHelper::Init(){
    fN=0;
    fZ=fRho=0;
    fZfilter=DBL_MAX;
    fShift=0.;
    fScale=1.;
    fOwnSurfaces=true;
    fFitAccuracy=-1.;
}

//______________________________________________________________________________
void FresnelSurfaceHelper::Read(const char* file){
    Clear();
    //printf("reading %s\n",file);

    int nread=10000;
    fRho=new double[nread];
    fZ=new double[nread];

    ifstream infile(file);
    if( !infile.good() ) {
        printf("Bad file %s\n", file);
        exit(1);
    }
    string line;
    fN=0;
    while( getline(infile,line,'\n') ){
        double rho,z;
        if( 2!=sscanf(line.data(), "%lg %lg", &rho, &z) ) continue;
        if ( fN>=nread ){
            double *newRho=new double[nread*2];
            double *newZ=new double[nread*2];
            memcpy(newRho, fRho, nread*sizeof(double));
            memcpy(newZ, fZ, nread*sizeof(double));
            nread*=2;
            delete [] fRho;
            delete [] fZ;
            fRho=newRho;
            fZ=newZ;
        }
        AddPoint(rho,z);
    }
    infile.close();
    Convert();
}

//______________________________________________________________________________
void FresnelSurfaceHelper::Read(int n, double *rho, double *z){
    Clear();
    fRho=new double[n];
    fZ=new double[n];
    fN=0;
    for(int i(0); i<n; i++) AddPoint(rho[i],z[i]);
    Convert();
}

//______________________________________________________________________________
void FresnelSurfaceHelper::Convert(){
    int previ=0;
    //printf("converting %i\n", fN);
    for(int i(0); i<fN; i++){
        if ( i<(fN-1) && fRho[i]!=fRho[i+1] ) continue;
        SSurface* surf=0;
        if ( fFitAccuracy>0. ){
            SecondOrderSurface* fit=new SecondOrderSurface();
            double dev=fit->GuessParameters(i-previ+1, fRho+previ, fZ+previ);
            if ( dev<fFitAccuracy ) surf=(SSurface*)fit;
        }
        if (!surf) surf=new ArraySurface(i-previ+1, fRho+previ, fZ+previ);
        //printf("%i -> %i\n", previ, i-previ+1);
        fSurfaces.push_back(surf);
        previ=i+1;
    }
    fOwnSurfaces=true;
}

//______________________________________________________________________________
FresnelSurface* FresnelSurfaceHelper::CreateFresnelLens(ParametricSurface* downsurf, ParametricSurface* upsurf){
    fOwnSurfaces=false;
    return new FresnelSurface(fSurfaces, downsurf, upsurf);
}