00001
00002
00003
00004 #include "Stiostream.h"
00005 #include <assert.h>
00006 #include "PAI.h"
00007 #include "StarCallf77.h"
00008 #include "TMath.h"
00009 #include "TRandom.h"
00010 #include "TError.h"
00011 #include "TH2.h"
00012 #include "TString.h"
00013
00014 #include "TGraph.h"
00015 #include "StMessMgr.h"
00016 #define xgasini F77_NAME(xgasini,XGASINI)
00017 #define xgastab F77_NAME(xgastab,XGASTAB)
00018 #define xfintera F77_NAME(xfintera,XFINTERA)
00019 PAI* PAI::fgInstance = 0;
00020
00021
00022 extern "C"
00023 {
00024 void type_of_call xgasini(Int_t &lun);
00025 void type_of_call xgastab(Float_t &g, Float_t &mdNdX, Int_t &mNoInTable, Float_t *mEnergy, Float_t *mdNdE);
00026 Float_t type_of_call xfintera(Float_t &x, Float_t *A, Float_t *F, Int_t &N);
00027 }
00028
00029 Float_t type_of_call xfintera(Float_t &x, Float_t *A, Float_t *F, Int_t &N) {
00030 Int_t bin = TMath::BinarySearch(N,A,x);
00031 if (bin < 0) bin = 0;
00032 if (bin > N - 3) bin = N-3;
00033 Double_t xlow = A[bin];
00034 Double_t xup = A[bin+1];
00035 Double_t dx = xup - xlow;
00036 if (dx < 1.e-5) {
00037 LOG_WARN << "bin\t" << bin << "\tN\t" << N << "\tx\t" << x << "\txlow\t" << xlow << "\txup\t" << xup << endm;
00038 }
00039 Double_t ylow = F[bin];
00040 Double_t yup = F[bin+1];
00041 Float_t y = ((xup*ylow-xlow*yup) + x*(yup-ylow))/dx;
00042 #ifdef DEBUG
00043 if (y > 20) {
00044 LOG_DEBUG << "bin\t" << bin << "\tN\t" << N << "\tx\t" << x
00045 << "\txlow\t" << xlow << "\txup\t" << xup << "\ty\t" << y << endm;
00046 }
00047 #endif
00048 return y;
00049 }
00050
00051 PAI *PAI::Instance(Int_t NoBetaGammas, Int_t NoEntries, Double_t BetaGammaLog10Min, Double_t BetaGammaLog10Max) {
00052 if (fgInstance) return fgInstance;
00053 fgInstance = new PAI(NoBetaGammas, NoEntries, BetaGammaLog10Min, BetaGammaLog10Max);
00054 return fgInstance;
00055 }
00056
00057 PAI::PAI(Int_t NoBetaGammas, Int_t NoEntries, Double_t BetaGammaLog10Min, Double_t BetaGammaLog10Max) :
00058 mNoBetaGammas(NoBetaGammas), mNoEntries(NoEntries), mNoInTable(0),
00059 mBetaGammaLog10Min(BetaGammaLog10Min), mBetaGammaLog10Max(BetaGammaLog10Max),
00060 mDeltaBetaGammaLog10(0),mEnergy(0), mdNdE(0), mdNdX(0), mGraphs(0) {
00061 Int_t lun = 6;
00062 xGasIni(lun);
00063 mEnergy = new Float_t [mNoEntries];
00064 mdNdE = new Float_t [mNoBetaGammas*mNoEntries];
00065 mdNdX = new Float_t [mNoBetaGammas];
00066 mDeltaBetaGammaLog10 = ( mBetaGammaLog10Max - mBetaGammaLog10Min )/ (mNoBetaGammas - 1);
00067 for (Int_t jg = 0; jg < mNoBetaGammas; jg++) {
00068 Float_t BetaGammaLog10 = mBetaGammaLog10Min + mDeltaBetaGammaLog10 * jg;
00069 Float_t BetaGamma = TMath::Power(10.,BetaGammaLog10);
00070 Float_t gamma = TMath::Sqrt(BetaGamma*BetaGamma + 1);
00071 xGasTab(gamma,mdNdX[jg],mNoInTable,&mEnergy[0],&mdNdE[jg*mNoEntries]);
00072 #ifdef DEBUG
00073 if (! mGraphs) {mGraphs = new TGraph*[ mNoBetaGammas ]; memset (mGraphs, 0, mNoBetaGammas*sizeof(TGraph *));}
00074 if (mGraphs[jg]) delete mGraphs[jg];
00075 mGraphs[jg] = new TGraph(mNoInTable);
00076 Int_t N = 0;
00077 LOG_DEBUG << "jg\t" << jg << "\tBetaGammaLog10\t" << BetaGammaLog10
00078 << "\tBetaGamma\t" << BetaGamma << "\tgamma\t" << gamma
00079 << "\tmdNdX\t" << mdNdX[jg]
00080 << "\tmNoInTable\t" << mNoInTable << endm;
00081 for (int i = 0; i < mNoInTable; i++) {
00082 LOG_DEBUG << "\t" << i << "\tmEnergy\t" << mEnergy[i] << "\tmdNdE\t" << mdNdE[jg*mNoEntries +i] << endm;
00083 }
00084 for (int i = 0; i < mNoInTable; i++) {
00085
00086 Double_t y = TMath::Exp(mEnergy[i]);
00087 Double_t z = mdNdE[jg*mNoEntries+i]/y;
00088 mGraphs[jg]->SetPoint(N++,y,z);
00089 }
00090 #endif
00091 }
00092 }
00093
00094 PAI::~PAI() {
00095 fgInstance = 0;
00096 delete [] mEnergy;
00097 delete [] mdNdE;
00098 delete [] mdNdX;
00099 #ifdef DEBUG
00100 for (Int_t i = 0; i < mNoBetaGammas; i++) SafeDelete(mGraphs[i]);
00101 delete [] mGraphs;
00102 #endif
00103 }
00104
00105 void PAI::xGasIni(Int_t &lun) {xgasini(lun);}
00106
00107 void PAI::xGasTab(Float_t &g, Float_t &mdNdX, Int_t &mNoInTable, Float_t *mEnergy, Float_t *mdNdE) {
00108 xgastab(g,mdNdX, mNoInTable, mEnergy, mdNdE);
00109 }
00110
00111 Float_t PAI::xFintera(Float_t &x, Float_t *A, Float_t *F, Int_t &N) {
00112 return xfintera(x, A, F, N);
00113 }
00114
00115 void PAI::xNext(Float_t BetaGamma, Float_t &dX, Float_t &dE) {
00116
00117 Float_t BetaGammaL = TMath::Log10(BetaGamma);
00118 Int_t jg = (Int_t) ((BetaGammaL - mBetaGammaLog10Min)/mDeltaBetaGammaLog10);
00119 if (jg < 0) jg = 0;
00120 if (jg >= mNoBetaGammas) jg = mNoBetaGammas-1;
00121 dX = - TMath::Log(gRandom->Rndm())/mdNdX[jg];
00122 Float_t random = gRandom->Rndm()*mdNdX[jg];
00123 dE = TMath::Exp(xFintera(random,&mdNdE[mNoEntries*jg],mEnergy,mNoInTable));
00124 #ifdef DEBUG
00125 LOG_DEBUG << "BetaGamma\t" << BetaGamma << "\trandom\t" << random << "\tdX\t" << dX << "\tdE\t" << dE << endm;
00126 #endif
00127 }
00128
00129 void PAI::xGenerate(Int_t NsSteps, Int_t Nevents) {
00130 Double_t *step = new Double_t [NsSteps];
00131 for (int i = 0; i < NsSteps; i++) step[i] = TMath::Power(2.,(i - 2.)/3.);
00132 for (int jg = 0; jg < mNoBetaGammas; jg++) {
00133 Double_t betaGammaLog10 = mBetaGammaLog10Min + mDeltaBetaGammaLog10*jg;
00134 Double_t betaGamma = TMath::Power(10., betaGammaLog10);
00135 TH1D **hist = new TH1D * [NsSteps];
00136 for (int i = 0; i < NsSteps; i++) {
00137 hist[i] = new TH1D(Form("dNdE_g%i_s%i",jg,i),
00138 Form("log10(dE/dx (keV/cm)) distribution for beta*gamma=%f and dX=%f",
00139 betaGamma,step[i]),500, -0.5, 4.5);
00140 LOG_INFO << "Generate Histograms\t" << hist[i]->GetName() << "\t" << hist[i]->GetTitle() << endm;
00141 }
00142 Float_t dX, dE;
00143 for (int iev = 0; iev < Nevents; iev++) {
00144 Double_t x = 0;
00145 Double_t E = 0;
00146 Int_t js = 0;
00147 while (js < NsSteps) {
00148 xNext(betaGamma, dX, dE);
00149
00150 x += dX;
00151 E += dE;
00152 for (int j = js; j < NsSteps; j++) {
00153 if (x > step[j]) {
00154 js++;
00155 Double_t dEdx = 1.e-3*E/step[j];
00156 if (dEdx > 0) hist[j]->Fill(TMath::Log10(dEdx));
00157 }
00158 }
00159 }
00160 }
00161 }
00162 }
00163
00164
00165
00166
00167
00168
00169
