void debye_temp(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;
   #define max 201;
   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]={0,5,20};
   Double_t dummyy[npts]={-1,-1,-1};

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

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


   Double_t TD=343.5;    // Debye temperature in Kelvin
   Double_t xmin=0;
   Double_t xmax=1;
   Double_t ymin=0;
   Double_t ymax=10;

   // 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(xmin,xmax);
   dummy->GetYaxis()->SetRangeUser(ymin,ymax);
   dummy->GetXaxis()->SetTitleSize(0.04);
   dummy->GetYaxis()->SetTitleSize(0.04);
   dummy->GetYaxis()->SetTitleOffset(1.5);
   dummy->GetXaxis()->SetTitle("T/T_{D}");
   dummy->GetYaxis()->SetTitle("C_{V} ()");
   dummy->Draw("Ap");
   dummy->GetXaxis()->SetNdivisions(505);

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

   TF1 *debye = new TF1("debye_func",debye_func,xmin,xmax,1);
   Double_t x=0;         // x is T/TD
   debye->SetParameter(0,TD);

   sprintf (string,"T_{D}=%.2f ^{o}K",TD);
   TF1 *debye1 = debye->DrawCopy("sameC");
   leg->AddEntry(debye1,string,"l");
   debye1->SetLineColor(2);

   Double_t Cv;
   Double_t R=8.31; // Gas constant in J/Kmole
   Double_t Acu=63.546e-3;    // Atomic weight of Cu in kg
   Double_t fudge=1.05;
   for (j=0;j<max;j++) {
       Double_t T = j*xmax*TD/(max-1);
       if (T < TD/500) T = TD/500;
       Double_t TDmax = TD/T;
       Double_t sum=debye->Integral(0,TDmax);
       Cv = fudge*(R/Acu)*9*(T/TD)*(T/TD)*(T/TD)*sum;
       printf ("j=%d T=%f TD=%f sum=%f, Cv=%f\n",j,T,TD,sum,Cv);
    }
    Double_t T=298;   // room temperature
    Double_t TDmax = TD/T;
    Double_t sum=debye->Integral(0,TDmax);
    Cv = fudge*(R/Acu)*9*(T/TD)*(T/TD)*(T/TD)*sum;
    printf ("\nCv at room temperature for Copper = %f J/kg/K \n\n",Cv);

   Double_t x=2.;   
   debye->SetParameter(0,TD); 
   sprintf (string,"T_{D}=%.2f ^{o}K",TD);
   // TF1 *debye2 = new TF1("debye2",debye1->Integral(0,x),xmin,xmax,1);
   /*TF1 *debye2 = debye->DrawCopy("sameC");
   leg->AddEntry(debye2,string,"l");
   debye2->SetLineColor(4);*/

   leg->Draw();

   Double_t alpha=1;
   sprintf(string,"#alpha/k =%.3f Kcm3/J\n",alpha);
   t1 = new TLatex(0.50,0.65,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.025);
   // t1->Draw(); 
  
   sprintf(filename,"debye_temp_c1.pdf");
   c1->SaveAs(filename);
   sprintf(filename,"debye_temp_c1.png");
   c1->SaveAs(filename);

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

   // Compute function inside integral in Debye function 

   Double_t TD=par[0];
   Double_t x1=x[0];
   Double_t pi=3.14159;

   char string[256];

   Double_t func;
   if (x1 < 0.0001) {
       func = x1*x1;
       }
   else if {
       func = (x1*x1*x1*x1*exp(x1))/((exp(x1)-1)*(exp(x1)-1));
       }

   /*sprintf (string,"x1=%f func=%f\n",x1,func);
   printf ("string=%s",string);*/

   return func;
}