---

CalProfile.cxx

Go to the documentation of this file.

// 
//  Original author: Regis Terrier
//                   terrier@cdf.in2p3.fr
//                                       December 1999


#include "reconstruction/CalProfile.h"

#include "reconstruction/GlastRecon.h"
#include "reconstruction/LbldData.h"
#include "reconstruction/CalRecon.h"
#include "reconstruction/TrackerRecon.h"
#include "reconstruction/gamma.h"
#include "reconstruction/ReconVisitor.h"

#include "data/GlastData.h"

#include "analysis/Midnight.h" 

#include "DetectorData.h" // for Tower, Calorimeter data

// The incomplete gamma function
//extern "C"{
//      extern double Gamma_(double* a,double* x);
//}

// The fitting function
// Declared statis below
//extern void fcn(int &npar, double *gin, double &f, double *u, int flag);

int nbins;
std::vector<double> layer;
std::vector<double> delta;
double slope,size,radlen;

double fit_energy;
double ki2;
double fit_alpha;
double fit_lambda;
double fit_start;
double energy_err;
double alpha_err;
double lambda_err;
double start_err;


static double gam_prof(double *par, int i)
{
        double result =0; 

        double length = ((2.3*i+par[2])/1.85)/slope;

        // Evaluation of the parameters of CsI at this energy
        // douple alpha = par[1];
        
        double alpha = 2.65 * exp(0.15*log(par[0]));
        double lambda = 2.29 * exp(-0.031*log(par[0]));

        double x=length/lambda;
        double dx = 2.3 / (1.85 *lambda)/slope;
        
        double gamma1 =0;
        double gamma2 = 0;

        // Now we will calculate the gamma incomplete function

        gamma1 = Gamma(alpha,x);
        x += dx;

        gamma2 = Gamma(alpha,x);        
        result = par[0]*(gamma2 - gamma1);

//      std::cout << " relstat :" << result << "\n";
        return result;
}



static void fcn(int &npar, double *gin, double &f, double *par, int iflag)
{
     int i;
//calculate chisquare
   double chisq = 0;
   double dlt;
   for (i=0;i<nbins; i++) {
     dlt  = (layer[i]-gam_prof(par,i));
     if(layer[i]>0.0001) chisq += dlt*dlt/layer[i];
   }
   f = chisq;
  // std::cout << " chisq=" << f << "\n";
}

void CalProfile::reconstruct(const GlastData* gdata)
{
        if( eTotal==0 ) return; //protect against numerical errors.

        nbins = static_cast<int>( num_layers );
        std::vector<double> E_layer = cal->Get_Elayer();
        layer.resize(num_layers);

        // layer goes from bottom to top

        for(int jlayer = 0; jlayer < num_layers; jlayer++) { // From top to bottom
              int ilayer = static_cast<int>( num_layers-1 - jlayer );
              layer[jlayer] = E_layer[ilayer];  // Energy per layer
            }

        delta.resize(num_layers);

        slope = fabs(dirZ);
        radlen = 1.85;
        size = 1.96;  

        // Vector of step, initial min and max value
        float vstrt[5];
        float stp[5];
        float bmin[5];
        float bmax[5];

        // par[0] is energy
        // par[1] is alpha
        // par[2] is the starting point
        
        // Init start values and bounds
        vstrt[0] = eTotal;
        if(eTotal>0.0001)       vstrt[1] = 2.65 * exp(0.15*log(eTotal));
        else vstrt[1]=0;
        vstrt[2] = 1.8f;                        // eq to 1X0 in CsI

        stp[0]=0.1f;
        stp[1]= 0.01f;
        stp[2] = 0.2f;

        bmin[0] = eTotal;
        bmin[1] = 1;
        bmin[2] = 0.5;

        bmax[0] = 5000;
        bmax[1] = 15;
        bmax[2] = 5;

        // Defines a object to do the minimization
        minuit = new Midnight(5);
        //Sets the function to be minimized
        minuit->SetFCN(fcn); 

    double arglist[10];
    int ierflg = 0;

    arglist[0] = 1;
    minuit->mnexcm("SET ERR", arglist ,1,ierflg);


        // Defines parameters
        
   minuit->mnparm(0, "Energy",    vstrt[0], stp[0], 0,0,ierflg);
   minuit->mnparm(1, "Alpha",    vstrt[1], stp[1], 0,0,ierflg);
   minuit->mnparm(2, "Starting point",  vstrt[2], stp[2], 0,0,ierflg);

   // Calls Migrad
   arglist[0]=500;
   arglist[1]=1;
   minuit->mnexcm("MIGRAD",arglist,2,ierflg);

   // Retrieve parameter information 

   minuit->GetParameter( 0, fit_energy,energy_err ); 
   minuit->GetParameter( 1, fit_alpha,alpha_err ); 
   minuit->GetParameter( 2, fit_start,start_err ); 
   minuit->GetParameter( 3, fit_lambda,lambda_err ); 

   arglist[0] = 3;
   minuit->mnexcm("CALL FCN", arglist ,1,ierflg);

   // Get chi-square
   double edm,errdef;
   int nvpar,nparx,icstat;
   minuit->mnstat(ki2,edm,errdef,nvpar,nparx,icstat);

   // Fills data
        data[i_CsI_FitEnergy] = fit_energy;
        data[i_Prof_Chisq] = ki2;
        data[i_CsI_Alpha] = fit_alpha;
        data[i_CsI_Lambda] = fit_lambda;



}


// implementation of the CalProfile class

CalProfile::CalProfile(CalRecon *c, TrackerRecon* t)
:Recon(),cal(c),trk(t)
{
    getParameters();
        // Energy from profile fitting
        i_CsI_FitEnergy= data.addElement("CsI_Fit_Energy");
        // Chi_square from profile fitting
        i_Prof_Chisq = data.addElement("Prof_Chisq");
        // Alpha profile parameter
        i_CsI_Alpha = data.addElement("CsI_Alpha");
        //Lambda profile parameter
        i_CsI_Lambda = data.addElement("CsI_Lambda");
        //Starting point profile parameter
        i_CsI_Start = data.addElement("CsI_Start");
        
        Init();

}

void CalProfile::getParameters ()
{

    num_layers = DetectorData::Calorimeter::numberOfLayers();
    s_rad_len = DetectorData::Calorimeter::radLen();
    mod_width = DetectorData::Tower::mod_width;
    s_wall_thickness = DetectorData::Tower::wall_thickness;
    s_wall_gap = DetectorData::Tower::wall_gap;
    s_xNum = DetectorData::Glast::xNum;

        // space between CALs parameters
        s_inner_thickness = DetectorData::Calorimeter::gridInnerThickness();
        s_side_thickness = DetectorData::Calorimeter::sideThickness();
        s_total_thickness = s_inner_thickness + s_side_thickness;
}


void CalProfile::Init()
{
        // Defines a object to do the minimization
        minuit = new Midnight(5);

        //Sets the function to be minimized
        minuit->SetFCN(fcn); 

        minuit->mninit(16,17,18);


        // Retrieve information from tracker and calorimeter reconstruction

        // First the reconstructed energy
        eTotal = getCal()->data.getData("CsI_Energy_Sum");    

        // Then the information on the shower direction and position
        dirX = getTrk()->data.getData("Gamma_Xdir");
        dirY = getTrk()->data.getData("Gamma_Ydir");
        dirZ = getTrk()->data.getData("Gamma_Zdir");

        centerX = getCal()->data.getData("CsI_X");
        centerY = getCal()->data.getData("CsI_Y");
        centerZ = getCal()->data.getData("CsI_Z");

}

void CalProfile::accept (ReconVisitor& rv)
{
    rv.visit(this);
}

void CalProfile::clear ()
{

    data.clear();
    eTotal = 0.0;
    nbins=0;
        layer;
        delta;
        slope=0;
        size=0;
        radlen=0;
        fit_energy=0;
        ki2=0;
        fit_alpha=0;
        fit_lambda=0;
}

---