Seeking guidance on implementing a user defined equation in ROOT

From: Russell Leslie <u4504546_at_anu.edu.au>
Date: Thu, 19 Feb 2009 12:40:26 +1100

Dear all,  

NOTE - message resent because I sent this originally from another account and it did not turn up on the roottalk list.  

I am still feeling my way forward on ROOT, but I have found it to be a very useful tool-set. One problem I am currently working on is to try to implement the following equations in ROOT (LaTex form of the equations below).

$$
G_{k}(t) = e^{-\Lambda^{*} t}G_{k}^{(fluct.)}(t)+ \int_{0}^{t} \Lambda^{*} G_{k}^{(fluct.)}(t')\cdotp e^{-\Lambda^{*} t'}G_{k}^{(stat.)} (t-t')dt' $$  

Where Gk(stat) is give by the following equation $$
G_{k}^{(stat.)} (t) = \left\langle \alpha_{k} \right\rangle + ( 1 - \left\langle \alpha_{k}\right\rangle )e^{- \Gamma t} $$  

and Gk(fluct) is given by the this equation $$
G_{k}^{(fluct.)}(t) = ( 1 - \lambda_{k}\tau_{c} )e^{ - \lambda _{k} t } $$  

The equations are from the manual for the Coulex code GOSIA and I want to use them as the fitting function for the output of a Monte Carlo simulation that we use.  

I have had no problems producing the curves for Gk(stat) and Gk(fluct) - but I have not been able to implement the integral term in the Gk equation (the first equation above). The code I have so far is given below  

{

//Calculate Gks
//parameter[0] = alpha_k
//parameter[1] = gamma

    TF1 *Gks = new TF1("Gks","[0]+(1-[0])*exp(-[1]*x)",0,15);     Gks->SetParameters(0.221243,0.4);
//Calculate Gkf
//parameter[0] = lambda_k
//parameter[1] = tau_c

    TF1 *Gkf = new TF1("Gkf","(1-[0]*[1])*exp(-[0]*x)",0,15);     Gkf->SetParameters(0.172,0.5);
//Calculate expLamba*t
//parameter[0] = 1

  //parameter[1] = Lambda*
    TF1 *eLst = new TF1("eLst","[0]*exp(-[1]*x)",0,15);     eLst->SetParameters(1,0.0345);     

//Combine Gks and Gkf to for Gkcalc

    TF1 *Gkcalc = new TF1("Gkcalc","Gkf*eLst*[0]*[1]",0,15); // This is the bit I can't get!!!!!!!!     

//Create a canvas to draw Gkcalc

    TCanvas *Gkcalc1 = new TCanvas("Gkcalc1","Gkcalc1",10,10,1400,1200);
//Draw and format Gkcalc

    Gkcalc->Draw();
     Gkcalc->SetLineWidth(2.0);

    Gkcalc->SetLineColor(kBlue);
    Gkcalc->GetXaxis()->SetTitle("time");
    Gkcalc->GetYaxis()->SetTitle("G_{k}(t)");
    Gkcalc->SetTitle("G_{k}(t)  = e^{-#Lambda^{*} t}G_{k}^{(fluct.)}(t)+ #int_{0}^{t} #Lambda^{*} G_{k}^{(fluct.)}(t').e^{-#Lambda^{*} t'}G_{k}^{(stat.)} (t-t')dt' ");
    Gkcalc->GetYaxis()->SetLabelFont(22);
    Gkcalc->GetYaxis()->SetTitleFont(22);
    Gkcalc->GetXaxis()->SetLabelFont(22);
    Gkcalc->GetXaxis()->SetTitleFont(22);

//Add other lines

    Gks->Draw("same");
     Gks->SetLineWidth(2.0);
    Gks->SetLineColor(kRed);
    Gkf->Draw("same");
     Gkf->SetLineWidth(2.0);
    Gkf->SetLineColor(kGreen);
    eLst->Draw("same");
     eLst->SetLineWidth(2.0);
    eLst->SetLineColor(kMagenta);
    TLegend *legend = new TLegend(.75,.80,.95,.95);

    legend->AddEntry(Gkcalc,"G_{k}(t)");
    legend->AddEntry(Gks,"G_{k}^{(stat)}(t)");
    legend->AddEntry(Gkf,"G_{k}^{(fluct)}(t)");
    legend->AddEntry(Gks,"e^{-#Lambda^{*} t}");
    legend->Draw();
   

}
My root details as as follows ROOT 5.18/00b (branches/v5-18-00-patches_at_22563, Oct 19 2008, 22:04:00 on linuxx8664gcc4.3)  

Any help or guidance you can give on this would be greatly appreciated.  

Russell Leslie
Nuclear Physics
R.S.Phys.S.E.
ANU Received on Thu Feb 19 2009 - 02:40:37 CET