void heat_transfer(void)
{
// 
// plot the number of fits to a cdc segment that result in prob > prob_cut (currently set to 0.01).
//

// #include <TMath.h>
#include <TRandom.h>

gROOT->Reset();
//TTree *Bfield = (TTree *) gROOT->FindObject("Bfield");
gStyle->SetPalette(1,0);
gStyle->SetOptStat(kFALSE);
gStyle->SetOptFit(kFALSE);
// gStyle->SetOptFit(1111);
gStyle->SetPadRightMargin(0.15);
gStyle->SetPadLeftMargin(0.15);
gStyle->SetPadBottomMargin(0.15);
gStyle->SetFillColor(0);
//

   char string[256];
   char filename[80];
   Int_t j,jj;
   
   #define npts 3;
   Int_t ngen=10000;

   // input data set first point (actually -100) to 100 so that it is not plotted.

   // use offset on x-axis to make data visible and give negative y-values for reference (not plotted)

   Double_t dummyx[npts]={-100,0,2000};
   Double_t dummyy[npts]={-1,-1,-1};

   //
   TCanvas *c1 = new TCanvas("c1","c1 heat transfer",200,10,700,700);
   c1->SetBorderMode(0);
   c1->SetFillColor(0);

   c1->SetGridx();
   c1->SetGridy();
   c1->SetBorderMode(0);
   c1->SetFillColor(0);
   c1->SetLogx();
   c1->SetLogy();


   Double_t tmin=0;
   Double_t tmax=10000;
   Double_t ymin=0.01;
   Double_t ymax=100;

   // dummy to draw axes

   TGraph *dummy = new TGraph (npts,dummyx,dummyy);
   TLegend *leg = new TLegend(0.35,0.74,0.85,0.9);
   dummy->SetMarkerColor(1);
   dummy->SetMarkerStyle(21);
   t1 = new TLatex(0.20,0.91,"Shapes");
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw(); 
   dummy->SetTitle("");
   dummy->GetXaxis()->SetRangeUser(tmin,tmax);
   dummy->GetYaxis()->SetRangeUser(ymin,ymax);
   dummy->GetXaxis()->SetTitleSize(0.04);
   dummy->GetYaxis()->SetTitleSize(0.04);
   dummy->GetYaxis()->SetTitleOffset(1.5);
   dummy->GetXaxis()->SetTitle("Time (s)");
   dummy->GetYaxis()->SetTitle("Temperature (^{o}C)");
   dummy->Draw("Ap");
   dummy->GetXaxis()->SetNdivisions(505);

   // plot temperature function. Consider x as a parameter and plot time dependence.

   TF1 *temp = new TF1("temp_func",temp_func,tmin,tmax,2);
   Double_t alpha=1.12e-4;    // unis of alpha is  m2/s,    1.12e-4 at room temperature
   // Double_t alpha=0.01;
   Double_t x=0;         // x in cm
   temp->SetParameter(0,alpha);
   temp->SetParameter(1,x);

   // Double_t sum = temp->Integral(tmin,tmax)*1e-3*1e-9/termination;
   sprintf (string,"x=%.2f m, #alpha=%g  m^{2}/s \n",x,alpha);
   TF1 *temp1 = temp->DrawCopy("sameC");
   leg->AddEntry(temp1,string,"l");
   temp1->SetLineColor(2);

   Double_t x=0.10;   
   temp->SetParameter(0,alpha);    
   temp->SetParameter(1,x); 
   sprintf (string,"x=%.2f m, #alpha=%g  m^{2}/s \n",x,alpha);
   TF1 *temp2 = temp->DrawCopy("sameC");
   leg->AddEntry(temp2,string,"l");
   temp2->SetLineColor(4);

   Double_t x=0.5;   
   temp->SetParameter(0,alpha);    
   temp->SetParameter(1,x); 
   sprintf (string,"x=%.2f m, #alpha=%g  m^{2}/s \n",x,alpha);
   TF1 *temp3 = temp->DrawCopy("sameC");
   leg->AddEntry(temp3,string,"l");
   temp3->SetLineColor(1);

   Double_t x=1.0;   
   temp->SetParameter(0,alpha);    
   temp->SetParameter(1,x); 
   sprintf (string,"x=%.2f m, #alpha=%g  m^{2}/s \n",x,alpha);
   TF1 *temp3 = temp->DrawCopy("sameC");
   leg->AddEntry(temp3,string,"l");
   temp3->SetLineColor(3);

   leg->Draw();

   sprintf(string,"#alpha =0.3f m^{2}/s\n",alpha);
   t1 = new TLatex(0.50,0.65,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.025);
   // t1->Draw(); 
  
   sprintf(filename,"heat_transfer_c1.pdf");
   c1->SaveAs(filename);
   sprintf(filename,"heat_transfer_c1.png");
   c1->SaveAs(filename);

}
Double_t temp_func (Double_t *x, Double_t *par) {


   // assume that initial conditions consist of u(x,0) = T0 for -a < x < a

   Double_t alpha=par[0];
   Double_t x1=par[1];
   Double_t t1=x[0];
   Double_t pi=3.14159;
   Double_t a=0.01;    // half-width of initial constant Temperature region, a is in meters
   Double_t T0= 100;

   char string[256];

   Double_t arg1 = (x1+a)/sqrt(2*alpha*t1);
   Double_t arg2 = (x1-a)/sqrt(2*alpha*t1);
   Double_t u1 = 0.5*T0*TMath::Erf(arg1);
   Double_t u2 = 0.5*T0*TMath::Erf(arg2);
   Double_t u;
   if (arg1 > 0 && arg2 >= 0) {
       u = u1 - u2;
       }
   else if (arg1<= 0 && arg2 <0) {
       u = -(u2 + u1);
       }
   else if (arg1 > =0 && arg2 <0) {
       u = u1 - u2;
       }
   else {
       u = 0;
   }
   
   
   /*sprintf (string,"t1=%f x=%f, arg1=%f, arg2=%f, u1=%f, u2=%f, u=%f\n",t1,x1,arg1,arg2,u1,u2,u);
   printf ("string=%s",string);*/

   return u;
}