00001 #include "EMFluctuation.h"
00002 #include "MessageService/MsgService.h"
00003 #include <fstream>
00004
00005 using namespace std;
00006
00007 CVSID("$Id: EMFluctuation.cxx,v 1.1 2005/04/14 18:10:31 zheng Exp $");
00008
00009 EMFluctuation::EMFluctuation(Double_t inputEnergy){
00010
00011 fPol4a = new TF1("Pol4a","pol4",0,20);
00012 fPol4b = new TF1("Pol4b","pol4",0,20);
00013 fUppPars_a = new double[5];
00014 fLowPars_a = new double[5];
00015 fUppPars_b = new double[5];
00016 fLowPars_b = new double[5];
00017 fUppEn = 0;
00018 fLowEn = 0;
00019 fInputEnergy = inputEnergy;
00020
00021 if(fInputEnergy>0) Init();
00022
00023 }
00024
00025 EMFluctuation::~EMFluctuation(){
00026 delete fPol4a;
00027 delete fPol4b;
00028 delete [] fUppPars_a;
00029 delete [] fLowPars_a;
00030 delete [] fUppPars_b;
00031 delete [] fLowPars_b;
00032 }
00033
00034 void EMFluctuation::ReInit(Double_t inputEnergy){
00035
00036 if(inputEnergy>0){
00037 fInputEnergy = inputEnergy;
00038 Init();
00039 }
00040 else MSG("EMFluctuation",Msg::kWarning)
00041 <<"Trying to ReInit with energy <=0 ... Sticking with previous energy "
00042 << fInputEnergy << endl;
00043
00044 }
00045
00046 void EMFluctuation::Init(){
00047
00049 int max = 12;
00050 double EnergyArray[12] = {1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0,15.0,20.0};
00051
00052 double LongPars[12][5] = {{1.06013,-0.468047,0.137782,-0.01263,0.000568776},
00053 {1.12785,-0.559151,0.145845,-0.0135959,0.00054578},
00054 {1.0812,-0.501449,0.117231,-0.0100283,0.000380176},
00055 {1.06712,-0.485045,0.106381,-0.00867457,0.000317992},
00056 {1.06773,-0.488333,0.104753,-0.00839656,0.000290656},
00057 {1.03906,-0.446111,0.0896104,-0.00680519,0.000234797},
00058 {1.04251,-0.456328,0.0911958,-0.0068878,0.000227099},
00059 {1.01184,-0.419094,0.0785208,-0.0054554,0.000173014},
00060 {1.08737,-0.488646,0.0980623,-0.00765436,0.000253581},
00061 {1.04148,-0.447376,0.0838367,-0.00591113,0.000179403},
00062 {1.05618,-0.460514,0.083876,-0.00590791,0.000178138},
00063 {1.06594,-0.496382,0.0932775,-0.00666602,0.000185565}};
00064
00065 double TranPars[12][5] = {{0.603255,-0.0112602,0.329161,0.00178921,-0.00318977},
00066 {0.43331,-0.0564151,0.239078,0.00330727,-0.00221715},
00067 {0.369569,0.00351878,0.186755,-0.000443012,-0.00154753},
00068 {0.322359,-0.0117416,0.160638,0.000834619,-0.00126155},
00069 {0.288158,-0.00238299,0.147826,0.000676277,-0.0012804},
00070 {0.263413,-0.00610981,0.132551,0.000452005,-0.00112661},
00071 {0.246907,-0.00220951,0.120681,0.000213735,-0.000922181},
00072 {0.231862,-0.00447519,0.113009,0.000142793,-0.000912504},
00073 {0.210691,0.00563696,0.111799,-1.47068e-05,-0.000990548},
00074 {0.205146,0.00240673,0.101715,-5.32772e-05,-0.000820729},
00075 {0.163379,-0.00469016,0.0837834,0.000147515,-0.00065458},
00076 {0.12919,-5.3113e-05,0.0770989,0.000163551,-0.000762262}};
00077
00079
00080
00081
00082 int LowVal = -1;
00083 for(int i=0;i<max;i++){
00084 if(fInputEnergy>EnergyArray[i]) LowVal = i;
00085 else break;
00086 }
00087
00088 for(int i=0;i<5;i++){
00089 if(LowVal == -1) {
00090 fLowPars_a[i] = LongPars[0][i];
00091 fUppPars_a[i] = LongPars[0][i];
00092 fLowPars_b[i] = TranPars[0][i];
00093 fUppPars_b[i] = TranPars[0][i];
00094 fLowEn = EnergyArray[0];
00095 fUppEn = EnergyArray[0];
00096 }
00097 else if(LowVal == max-1) {
00098 fLowPars_a[i] = LongPars[LowVal][i];
00099 fUppPars_a[i] = LongPars[LowVal][i];
00100 fLowPars_b[i] = TranPars[LowVal][i];
00101 fUppPars_b[i] = TranPars[LowVal][i];
00102 fLowEn = EnergyArray[LowVal];
00103 fUppEn = EnergyArray[LowVal];
00104 }
00105 else {
00106 fLowPars_a[i] = LongPars[LowVal][i];
00107 fUppPars_a[i] = LongPars[LowVal+1][i];
00108 fLowPars_b[i] = TranPars[LowVal][i];
00109 fUppPars_b[i] = TranPars[LowVal+1][i];
00110 fLowEn = EnergyArray[LowVal];
00111 fUppEn = EnergyArray[LowVal+1];
00112 }
00113 }
00114 }
00115
00116 Double_t EMFluctuation::CalcLongFluc(Double_t t){
00117
00118
00119 t/=0.06;
00120
00121 for(Int_t i=0;i<5;i++) fPol4a->SetParameter(i,fLowPars_a[i]);
00122 Double_t fracErr_low = fPol4a->Eval(t);
00123 if(fUppEn==fLowEn) return fracErr_low;
00124
00125 for(Int_t i=0;i<5;i++) fPol4a->SetParameter(i,fUppPars_a[i]);
00126 Double_t fracErr_upp = fPol4a->Eval(t);
00127
00128 Double_t fracErr = fracErr_low +
00129 (fInputEnergy - fLowEn)*(fracErr_upp-fracErr_low)/(fUppEn-fLowEn);
00130 return fracErr;
00131
00132 }
00133
00134 Double_t EMFluctuation::CalcTranFluc(Double_t r){
00135
00136
00137 r/=0.041;
00138
00139
00140
00141 if(r>5) r = 5;
00142
00143 for(Int_t i=0;i<5;i++) fPol4b->SetParameter(i,fLowPars_b[i]);
00144 Double_t fracErr_low = fPol4b->Eval(r);
00145 if(fUppEn==fLowEn) return fracErr_low;
00146
00147 for(Int_t i=0;i<5;i++) fPol4b->SetParameter(i,fUppPars_b[i]);
00148 Double_t fracErr_upp = fPol4b->Eval(r);
00149
00150 Double_t fracErr = fracErr_low +
00151 (fInputEnergy - fLowEn)*(fracErr_upp-fracErr_low)/(fUppEn-fLowEn);
00152 return fracErr;
00153
00154 }