BinFluctuationEM.cxx

Go to the documentation of this file.

#include "BinFluctuationEM.h"
#include "MessageService/MsgService.h"
#include <fstream>

using namespace std;

CVSID("$Id: BinFluctuationEM.cxx,v 1.2 2004/09/29 00:07:29 cbs Exp $");

BinFluctuationEM::BinFluctuationEM(Double_t inputEnergy){

  fPol4a = new TF1("Pol4a","pol4",0,20);
  fPol4b = new TF1("Pol4b","pol4",0,20);
  fUppPars_a = new double[5];
  fLowPars_a = new double[5];
  fUppPars_b = new double[5];
  fLowPars_b = new double[5];
  fUppEn = 0;
  fLowEn = 0;
  fInputEnergy = inputEnergy;

  if(fInputEnergy>0) Init();

}

BinFluctuationEM::~BinFluctuationEM(){
  delete fPol4a;
  delete fPol4b;
  delete [] fUppPars_a;
  delete [] fLowPars_a;
  delete [] fUppPars_b;
  delete [] fLowPars_b;
}

void BinFluctuationEM::ReInit(Double_t inputEnergy){
  
  if(inputEnergy>0){
    fInputEnergy = inputEnergy;
    Init();
  }
  else MSG("FitShowerEM",Msg::kWarning)
    <<"Trying to ReInit with energy <=0 ... Sticking with previous energy "
    << fInputEnergy << endl;
  
}

void BinFluctuationEM::Init(){
  
  int max = 9;
  double EnergyArray[9] = {1.0,2.0,4.0,5.0,6.0,7.0,8.0,10.0,15.0};
  
  double a[9][5] = {{1.815,-0.8182,0.259,-0.02858,0.001549},     //1GeV
                    {1.971,-1.047,0.2918,-0.03065,0.001351},     //2GeV
                    {1.736,-0.8027,0.1856,-0.01704,0.0007054},   //4 GeV
                    {1.775,-0.8145,0.182,-0.01619,0.0006334},    //5GeV
                    {1.82,-0.8121,0.1768,-0.01586,0.0006187},    //6GeV
                    {1.642,-0.7629,0.1722,-0.01585,0.0006154},   //7 GeV
                    {1.888,-0.889,0.1883,-0.01597,0.0005586},    //8GeV
                    {1.889,-0.9553,0.2084,-0.01823,0.0006304},   //10GeV
                    {1.808,-0.8572,0.1682,-0.01332,0.0004303}};  //15GeV
  
  double b[9][5] = {{1.276,-0.5111,0.1116,-0.01024,0.0003537},   //1GeV
                    {1.112,-0.4453,0.08835,-0.007366,0.0002318}, //2GeV
                    {0.8736,-0.3052,0.04869,-0.003169,8.005e-5}, //4 GeV  
                    {0.8617,-0.3493,0.06396,-0.004886,0.0001392},//5GeV
                    {0.7246,-0.2662,0.04372,-0.002875,6.989e-5}, //6GeV
                    {0.824,-0.3194,0.05317,-0.003619,9.042e-5},  //7 GeV
                    {0.6148,-0.1982,0.02825,-0.001624,3.785e-5}, //8GeV
                    {0.6669,-0.2314,0.03431,-0.002067,4.71e-5},  //10GeV
                    {0.5692,-0.2021,0.03125,-0.002036,5.037e-5}}; //15GeV
  
  
  //choose Pol4's  
  int LowVal = -1;
  for(int i=0;i<max;i++){
    if(fInputEnergy>EnergyArray[i]) LowVal = i;
    else break;
  }
  
  for(int i=0;i<5;i++){
    if(LowVal == -1) {
      fLowPars_a[i] = a[0][i];
      fUppPars_a[i] = a[0][i];
      fLowPars_b[i] = b[0][i];
      fUppPars_b[i] = b[0][i];
      fLowEn = EnergyArray[0];
      fUppEn = EnergyArray[0];
    }
    else if(LowVal == max-1) {
      fLowPars_a[i] = a[LowVal][i];
      fUppPars_a[i] = a[LowVal][i];
      fLowPars_b[i] = b[LowVal][i];
      fUppPars_b[i] = b[LowVal][i];
      fLowEn = EnergyArray[LowVal];
      fUppEn = EnergyArray[LowVal];
    }
    else {
      fLowPars_a[i] = a[LowVal][i];
      fUppPars_a[i] = a[LowVal+1][i];
      fLowPars_b[i] = b[LowVal][i];
      fUppPars_b[i] = b[LowVal+1][i];
      fLowEn = EnergyArray[LowVal];
      fUppEn = EnergyArray[LowVal+1];
    }
  }
}

Double_t BinFluctuationEM::CalcFluctuation(Double_t t, Double_t r){
  
  for(int i=0;i<5;i++){
    fPol4a->SetParameter(i,fLowPars_a[i]);
    fPol4b->SetParameter(i,fLowPars_b[i]);
  }
  double p0_low = fPol4a->Eval(t);
  double p2_low = fPol4b->Eval(t);
  double fracErr_low = p0_low + r*r*p2_low;  

  if(fUppEn==fLowEn) return fracErr_low;
  
  for(int i=0;i<5;i++){
    fPol4a->SetParameter(i,fUppPars_a[i]);
    fPol4b->SetParameter(i,fUppPars_b[i]);
  }
  double p0_upp = fPol4a->Eval(t);
  double p2_upp = fPol4b->Eval(t);
  double fracErr_upp = p0_upp + r*r*p2_upp;
  
  double fracErr = fracErr_low
    + (fInputEnergy - fLowEn)*(fracErr_upp-fracErr_low)/(fUppEn-fLowEn);
  
  return fracErr;

}

---