void thermal_properties(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 jpts 19;
   #define mpts 25;
   #define npts 3;
   #define max 201;
   # define npts_cp 10;
   Int_t ngen=10000;

   // temperature is in K, thermal conductivity in W/m-K, specific heat in J/kg-K

   // Double_t temperature[mpts]={1,2,3,4,5,6,7,8,9,10,15,20,30,40,50,60,70,80,90,100,150,200,250,300,350,400,500,600,800,1000,1200};
   // Double_t conductivity[mpts]={4220,8400,12500,16200,19500,22200,23900,24800,24900,24300,17100,10800,4450,2170,1250,829,647,557,508,482,429,413,406,401,396,393,386,379,366,352,339};
   Double_t temperature[mpts]={1,4,6,8,10,15,20,25,30,35,40,50,60,70,76,80,90,100,120,140,160,180,200,250,300};
   Double_t conductivity[mpts]={-1,70,96,120,134,120,88,60,40,28,20,12,8,6.2,5.7,5.2,4.7,4.5,4.3,4.2,4.1,4.0,4.0,4.0,4.0};
   Double_t conductivity_ofhc[mpts]={-1,2.4,3.7,4.7,6.0,8.5,11,12,12,11,10,7.7,6.2,5.5,5.2,4.9,4.7,4.5,4.3,4.2,4.1,4.0,4.0,4.0,4.0};
   Double_t specific_heat[mpts];
   Double_t thermal_diffusivity[mpts];
   Double_t thermal_diffusivity_ofhc[mpts];
   Double_t Cu_rho=8.96e3 ; // units are kg/m3
   Double_t cu_cpt[npts_cp]={1,2,4,6,8,10,20,30,100,300};
   Double_t cu_cp[npts_cp]={0.012,0.028,0.091,0.23,0.47,0.86,7.7,27,254,386};   // BNL Cryo p. VIII-B-2 J/kg-K

   Double_t etemp[jpts]={1,10,20,40,60,80,100,150,200,273,293,298,300,400,500,600,700,800,900};
   Double_t electrical_resistivity[jpts]={0.002,0.00202,0.0028,0.0239,0.0971,0.215,0.348,0.699,1.046,1.543,1.678,1.712,1.725,2.402,3.09,3.792,4.514,5.262,6.041,};

   Double_t L_WK = 2.44e-8;  //  Lorentz number, proportionality constant in the Wiedemann-Franz Law.
   Double_t electrical_resistivity_WK[mpts];    // calculation of resistivity from thermal condctivity and Wiedemann-Franz law.
   Double_t electrical_resistivity_WK_ofhc[mpts];    // calculation of resistivity from thermal condctivity and Wiedemann-Franz law.
   Double_t resistivity_formula[mpts];  // use formula from George's reference on MIITS calculation. (4.3-13)

   //
   TCanvas *c1 = new TCanvas("c1","c1 thermal_properties",200,10,700,700);

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

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


   Double_t TD=343.5*0.85;    // Debye temperature in Kelvin, emperical factor to reproduce data at about 4 K.
   Double_t xmin=0;
   Double_t xmax=500;
   Double_t ymin=0;
   Double_t ymax=10;

   // 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 RRR=75;
   for (j=0;j<mpts;j++) {
       Double_t T = temperature[j];
       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;
       Cv = (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);
       specific_heat[j] = Cv;
       conductivity[j]=100*conductivity[j]; // convert from W/cm-K to W/m-K
       conductivity_ofhc[j]=100*conductivity_ofhc[j]; // convert from W/cm-K to W/m-K
       // conductivity_ofhc[j]=50*conductivity_ofhc[j]; // convert from W/cm-K to W/m-K and /2 for magneto resistive effect
       thermal_diffusivity[j] = conductivity[j]/(specific_heat[j]*Cu_rho);  // use pure copper numbers
       thermal_diffusivity_ofhc[j] = conductivity_ofhc[j]/(specific_heat[j]*Cu_rho);  // OFHC copper numbers
       electrical_resistivity_WK[j] = 0.943* L_WK*T*1e8/conductivity[j];   // units of 10^-8 scale by 0.943 so at 300 it is 1.725
       electrical_resistivity_WK_ofhc[j] = 0.943*L_WK*T*1e8/conductivity_ofhc[j];   // units of 10^-8  scale by 0.943 so at 300 it is 1.725
       resistivity_formula[j] = 1.545/RRR*(1+ 1/(2.32547e17*pow(T,-5) + 9.57137e13*pow(T,-3) + 1.62735e10*T));
       printf (" T=%g, heat=%g, resistivity=%g, diffusivity=%g formula=%g\n",T,conductivity[j],electrical_resistivity_WK[j],thermal_diffusivity[j],resistivity_formula[j]);
    }
  
   // fudge on last point to adjust y scale
   // conductivity[mpts-1] = 100;

   Double_t ymin=100;
   Double_t ymax=20000;
   TGraph *conduct = new TGraph (mpts,temperature,conductivity);
   TGraph *conduct_ofhc = new TGraph (mpts,temperature,conductivity_ofhc);
   conduct->SetTitle("");
   conduct->GetXaxis()->SetRangeUser(xmin,xmax);
   conduct->GetYaxis()->SetRangeUser(ymin,ymax);
   conduct->GetXaxis()->SetTitleSize(0.04);
   conduct->GetYaxis()->SetTitleSize(0.04);
   conduct->GetYaxis()->SetTitleOffset(1.5);
   // conduct->GetYaxis()->SetTitleOffset(2.0);
   conduct->GetXaxis()->SetTitle("Temperature (K)");
   conduct->GetYaxis()->SetTitle("Thermal Conductivity k (W/m-K)");
   conduct->GetXaxis()->SetNdivisions(505);
   conduct->SetMarkerColor(1);
   conduct->SetMarkerStyle(20);
   conduct->SetMarkerSize(0.8);
   conduct->Draw("Ap");

   conduct_ofhc->SetMarkerColor(2);
   conduct_ofhc->SetMarkerStyle(20);
   conduct_ofhc->SetMarkerSize(0.8);
   conduct_ofhc->Draw("Samep");

   // plot k from empirical resistivity and the Wiedemann-Franz law

   Double_t thermal_empirical[1];     
   Double_t temp_empirical[1]={4};
   Double_t T = temp_empirical[0];
   Double_t res = 1.7*(2/75.);     // units if 10-8, room temperature value x magneto resistive effect / RRR.
   thermal_empirical[0] = L_WK*T*1e8/res;
   TGraph *thermal_emp = new TGraph(1,temp_empirical,thermal_empirical);
   thermal_emp->SetMarkerColor(2);
   thermal_emp->SetMarkerStyle(21);
   thermal_emp->Draw("Samep");
   printf ("T=%g, L_WK=%g, res=%g, thermal_emp=%g\n",T,L_WK,res,thermal_empirical[0]);

   //
   // TCanvas *c2 = new TCanvas("c2","c2 thermal_properties",200,10,700,700);

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

   TLegend *legcp = new TLegend(0.30,0.15,0.85,0.25);
   TGraph *specific = new TGraph (mpts,temperature,specific_heat);
   specific->SetMarkerColor(2);
   specific->SetMarkerStyle(20);
   specific->SetMarkerSize(0.8);
   specific->SetTitle("");
   specific->GetXaxis()->SetRangeUser(xmin,xmax);
   // specific->GetYaxis()->SetRangeUser(ymin,ymax);
   specific->GetXaxis()->SetTitleSize(0.04);
   specific->GetYaxis()->SetTitleSize(0.04);
   specific->GetYaxis()->SetTitleOffset(1.5);
   specific->GetXaxis()->SetTitle("Temperature (K)");
   specific->GetYaxis()->SetTitle("Specific Heat c_{p} (J/kg-K)");
   specific->GetXaxis()->SetNdivisions(505);
   legcp->AddEntry(specific,"Debye Formula Cu","p");
   specific->Draw("Ap");

   TGraph *bnl_cp = new TGraph(npts_cp,cu_cpt,cu_cp);
   bnl_cp->SetMarkerColor(4);
   bnl_cp->SetMarkerStyle(20);
   bnl_cp->SetMarkerSize(0.8);
   bnl_cp->Draw("Samep");
   legcp->AddEntry(bnl_cp,"BNL Data Cu","p");
   legcp->Draw("");
   
   


   //
   // TCanvas *c3 = new TCanvas("c3","c3 thermal_properties",200,10,700,700);

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

   Double_t ymin=0.0001;
   Double_t ymax=20;
   TGraph *diffusivity = new TGraph (mpts,temperature,thermal_diffusivity);
   TGraph *diffusivity_ofhc = new TGraph (mpts,temperature,thermal_diffusivity_ofhc);
   diffusivity->SetTitle("");
   diffusivity->GetXaxis()->SetRangeUser(xmin,xmax);
   diffusivity->GetYaxis()->SetRangeUser(ymin,ymax);
   diffusivity->GetXaxis()->SetTitleSize(0.04);
   diffusivity->GetYaxis()->SetTitleSize(0.04);
   diffusivity->GetYaxis()->SetTitleOffset(1.5);
   diffusivity->GetXaxis()->SetTitle("Temperature (K)");
   diffusivity->GetYaxis()->SetTitle("Thermal Diffusivity #alpha (m^{2}/s)");
   diffusivity->GetXaxis()->SetNdivisions(505);
   diffusivity->SetMarkerColor(1);
   diffusivity->SetMarkerStyle(20);
   diffusivity->SetMarkerSize(0.8);
   diffusivity->Draw("Ap");

   diffusivity_ofhc->SetMarkerColor(2);
   diffusivity_ofhc->SetMarkerStyle(20);
   diffusivity_ofhc->SetMarkerSize(0.8);
   diffusivity_ofhc->Draw("Samep");

   //
   // TCanvas *c4 = new TCanvas("c4","c4 thermal_properties",200,10,700,700);

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

   TLegend *leg = new TLegend(0.15,0.85,0.85,1.0);
   TGraph *resistivity = new TGraph (jpts,etemp,electrical_resistivity);
   TGraph *resistivity_WK = new TGraph (mpts,temperature,electrical_resistivity_WK);
   TGraph *resistivity_WK_ofhc = new TGraph (mpts,temperature,electrical_resistivity_WK_ofhc);
   Double_t resistivity_empirical[1]={1.7*2/75};     // units if 10-8, room temperature value x magneto resistive effect / RRR.
   Double_t temp_empirical[1]={4};
   TGraph *resistivity_emp = new TGraph(1,temp_empirical,resistivity_empirical);
   TGraph *resistivity_for = new TGraph(mpts,temperature,resistivity_formula);

   Double_t ymin=0.001;
   Double_t ymax=5;
   resistivity_WK->SetTitle("");
   resistivity_WK->GetXaxis()->SetRangeUser(xmin,xmax);
   resistivity_WK->GetYaxis()->SetRangeUser(ymin,ymax);
   resistivity_WK->GetXaxis()->SetTitleSize(0.04);
   resistivity_WK->GetYaxis()->SetTitleSize(0.04);
   resistivity_WK->GetYaxis()->SetTitleOffset(1.5);
   resistivity_WK->GetXaxis()->SetTitle("Temperature (K)");
   resistivity_WK->GetYaxis()->SetTitle("Electrical Resistivity #rho (10^{-8}#Omega-m)");
   resistivity_WK->GetXaxis()->SetNdivisions(505);
   resistivity_WK->SetMarkerColor(4);
   resistivity_WK ->SetMarkerStyle(20);
   resistivity_WK->SetMarkerSize(0.8);
   resistivity_WK->Draw("Ap");
   sprintf (string,"k & Wiedemann-Franz Law Pure Cu\n");
   leg->AddEntry(resistivity_WK,string,"p");

   resistivity_WK_ofhc->SetMarkerColor(2);
   resistivity_WK_ofhc ->SetMarkerStyle(20);
   resistivity_WK_ofhc->SetMarkerSize(0.8);
   resistivity_WK_ofhc->Draw("Samep");
   sprintf (string,"k & Wiedemann-Franz Law OFHC\n");
   leg->AddEntry(resistivity_WK_ofhc,string,"p");

   /*resistivity_for->SetMarkerColor(1);
   resistivity_for->SetMarkerStyle(20);
   resistivity_for->Draw("Samep");
   sprintf (string,"Formula 4.3-13\n");
   leg->AddEntry(resistivity_for,string,"p");*/

   resistivity->SetMarkerColor(1);
   resistivity->SetMarkerStyle(20);
   resistivity->SetMarkerSize(0.8);
   resistivity->Draw("Samep");
   sprintf (string,"Data on pure Cu\n");
   leg->AddEntry(resistivity,string,"p");

   resistivity_emp->SetMarkerColor(2);
   resistivity_emp->SetMarkerStyle(21);
   // resistivity_emp->SetMarkerSize(0.8);
   resistivity_emp->Draw("Samep");
   sprintf (string,"Scaled from Room T by 2/75\n");
   leg->AddEntry(resistivity_emp,string,"p");

  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(); 


   TCanvas *c4 = new TCanvas("c4","c4 thermal_properties",200,10,700,350);

   c4->SetGridx();
   c4->SetGridy();
   c4->SetBorderMode(0);
   c4->SetFillColor(0);
   c4->SetLogx();
   c4->SetLogy();
   c4->Divide(2,1);

   c4->cd(1);
   c4_1->SetGridx();
   c4_1->SetGridy();
   c4_1->SetBorderMode(0);
   c4_1->SetFillColor(0);
   c4_1->SetLogx();
   c4_1->SetLogy();


   Double_t xmin=1;
   Double_t xmax=500;
   Double_t ymin=100;
   Double_t ymax=20000;
   conduct_ofhc->SetTitle("");
   conduct_ofhc->GetXaxis()->SetRangeUser(xmin,xmax);
   conduct_ofhc->GetYaxis()->SetRangeUser(ymin,ymax);
   conduct_ofhc->GetXaxis()->SetTitleSize(0.04);
   conduct_ofhc->GetYaxis()->SetTitleSize(0.04);
   conduct_ofhc->GetYaxis()->SetTitleOffset(1.5);
   // conduct_ofhc->GetYaxis()->SetTitleOffset(2.0);
   conduct_ofhc->GetXaxis()->SetTitle("Temperature (K)");
   conduct_ofhc->GetYaxis()->SetTitle("Thermal Conductivity k (W/m-K)");
   conduct_ofhc->GetXaxis()->SetNdivisions(505);
   conduct_ofhc->SetMarkerColor(2);
   conduct_ofhc->SetMarkerStyle(21);

   Double_t p0=log(1000);
   Double_t p1=-1.5;
   Double_t p2=1;
   Double_t p3=-1;
   Double_t p4=0.;
   Double_t fitmin=4;
   Double_t fitmax=70;
   TF1 *fit = new TF1("fit_func",fit_func,fitmin,fitmax,5);
   fit->SetParameters(p0,p1,p2,p3,p4);
   fit->FixParameter(4,p4);

   conduct_ofhc->Fit("fit_func","","",fitmin,fitmax);
   conduct_ofhc->Draw("Ap");
   p0 = fit->GetParameter(0);
   p1 = fit->GetParameter(1);
   p2 = fit->GetParameter(2);
   p3 = fit->GetParameter(3);
   p4 = fit->GetParameter(4);
   printf ("Fit parameters: p0 %f p1 %f  p2 %f p3 %f p4 %f\n",p0,p1,p2,p3,p4);

   // TCanvas *c5 = new TCanvas("c5","c5 thermal_properties",200,10,700,700);


   c4->cd(2);
   c4_2->SetGridx();
   c4_2->SetGridy();
   c4_2->SetBorderMode(0);
   c4_2->SetFillColor(0);
   c4_2->SetLogx();
   c4_2->SetLogy();


   /*c5->SetGridx();
   c5->SetGridy();
   c5->SetBorderMode(0);
   c5->SetFillColor(0);
   c5->SetLogx();
   c5->SetLogy();*/

   Double_t xmin=0.01;
   Double_t xmax=100;
   Double_t ymin=0.001;
   Double_t ymax=10;

   Double_t ACmag=0.5;   // Units are W/cm/K  
   Double_t ACslope=0.09;    // units are per deg K
   Double_t ACoption=1;   // use to choose either uppper or lower half of heat transfer curve

   TF1 *transfer = new TF1("transfer_func",transfer_func,xmin,xmax,3);
   transfer->SetParameters(ACmag,ACslope,ACoption);
   transfer->SetTitle("");
   transfer->GetXaxis()->SetRangeUser(xmin,xmax);
   transfer->GetYaxis()->SetRangeUser(ymin,ymax);
   transfer->GetXaxis()->SetTitleSize(0.04);
   transfer->GetYaxis()->SetTitleSize(0.04);
   transfer->GetYaxis()->SetTitleOffset(1.5);
   transfer->GetXaxis()->SetTitle("#Delta T (K)");
   transfer->GetYaxis()->SetTitle("Heat Transfer q (W/cm^{2})");
   transfer->GetXaxis()->SetNdivisions(505);
   transfer->SetMarkerColor(2);
   transfer->SetMarkerStyle(21);

   transfer->Draw();
   ACoption = 0;
   TF1 *transfer1 = transfer->DrawCopy("samel");
   transfer1->SetParameter(2,ACoption);
   transfer1->SetLineStyle(2);

   TCanvas *c5 = new TCanvas("c5","c5 thermal_properties",200,10,700,700);

   c5->SetGridx();
   c5->SetGridy();
   c5->SetBorderMode(0);
   c5->SetFillColor(0);
   c5->SetLogx();
   c5->SetLogy();
   // c5->Divide(2,1);

   /*c5->cd(1);
   c5_1->SetGridx();
   c5_1->SetGridy();
   c5_1->SetBorderMode(0);
   c5_1->SetFillColor(0);
   c5_1->SetLogx();
   c5_1->SetLogy();*/

   // resistivity_WK->Draw("Ap");

   Double_t ymin=0.001;
   Double_t ymax=5;
   Double_t xmin=0.001;
   Double_t xmax=300;
   resistivity->SetTitle("");
   resistivity->GetXaxis()->SetRangeUser(xmin,xmax);
   resistivity->GetYaxis()->SetRangeUser(ymin,ymax);
   resistivity->GetXaxis()->SetTitleSize(0.04);
   resistivity->GetYaxis()->SetTitleSize(0.04);
   resistivity->GetYaxis()->SetTitleOffset(1.5);
   resistivity->GetXaxis()->SetTitle("Temperature (K)");
   resistivity->GetYaxis()->SetTitle("Electrical Resistivity #rho (10^{-8}#Omega-m)");
   resistivity->GetXaxis()->SetNdivisions(505);
   resistivity->SetMarkerColor(4);
   resistivity ->SetMarkerStyle(20);
   resistivity->SetMarkerSize(0.8);
   resistivity->Draw("Ap");

   TLegend *leg1 = new TLegend(0.17,0.58,0.45,0.88);

   /*sprintf (string,"k & WF Law for Pure Cu\n");
   leg1->AddEntry(resistivity_WK,string,"p");*/

   sprintf (string,"WF Law & k for OFHC Cu\n");
   leg1->AddEntry(resistivity_WK_ofhc,string,"p");

   sprintf (string,"Data on pure Cu\n");
   leg1->AddEntry(resistivity,string,"p");

  leg1->Draw();


   Double_t rho=4.244;
   Double_t slope=1.425;
   Double_t cc2=-0.3722;
   Double_t cc3=0.27;
   Double_t cc4=0;
   Double_t scale=1;
   Double_t mean=log(75);
   Double_t RRR=42;
   Double_t fitmin=4;
   Double_t fitmax=300;
   TF1 *fit_rho = new TF1("fit_rho",fit_rho_func,fitmin,fitmax,7);
   fit_rho->SetParameters(rho,slope,cc2,cc3,scale,mean,RRR);
   // fit_rho->FixParameter(0,rho);
   // fit_rho->FixParameter(1,slope);
   // fit_rho->FixParameter(2,cc2);
   // fit_rho->FixParameter(3,cc3);
   fit_rho->FixParameter(4,scale);
   fit_rho->FixParameter(5,mean);
   fit_rho->FixParameter(6,RRR);
   // resistivity_WK_ofhc->Fit("fit_rho","WW","",fitmin,fitmax);
   resistivity_WK_ofhc->Draw("samep");

   rho = fit_rho->GetParameter(0);
   slope = fit_rho->GetParameter(1);
   cc2 = fit_rho->GetParameter(2);
   cc3 = fit_rho->GetParameter(3);
   scale = fit_rho->GetParameter(4);
   mean = fit_rho->GetParameter(5);
   RRR = fit_rho->GetParameter(6);
   printf ("Fit parameters: rho %f slope %f  cc2 %f cc3 %f cc4 %f mean %f RRR %f\n",rho,slope,cc2,cc3,cc4,mean,RRR);

 
   /*Double_t rho=8.893;
   Double_t slope=1.349;
   Double_t cc2=-0.7459;
   Double_t cc3=0.27;
   Double_t cc4=0;
   Double_t scale=1;
   Double_t mean=log(75);
   Double_t RRR=850;
   Double_t fitmin=4;
   Double_t fitmax=300;
   TF1 *fit_rho2 =  new TF1("fit_rho2",fit_rho_func,fitmin,fitmax,7);
   fit_rho2->SetParameters(rho,slope,cc2,cc3,scale,mean,RRR);
   fit_rho2->FixParameter(0,rho);   
   fit_rho2->FixParameter(1,slope);
   fit_rho2->FixParameter(2,cc2);
   fit_rho2->FixParameter(3,cc3);
   fit_rho2->FixParameter(4,scale);
   fit_rho2->FixParameter(5,mean);
   fit_rho2->FixParameter(6,RRR);
   resistivity->Fit("fit_rho2","WW","",fitmin,fitmax);
   resistivity->Draw("samep");*/

   Double_t rho=4.202;
   Double_t slope=1.61;
   Double_t cc2=-0.6352;
   Double_t cc3=0.4155;
   Double_t scale=1;
   Double_t mean=log(75);

   Double_t Bfield=0;

   TF1 *Rcu=  new TF1("Rcu",cu_resistance,fitmin,fitmax,2);
   /*TF1 *RRR10 = Rcu->DrawCopy("sameC");
   RRR = 10;
   RRR10->SetParameters(RRR,Bfield);
   RRR10->SetLineColor(2);*/

   TF1 *RRR50 = Rcu->DrawCopy("sameC");
   RRR = 50;
   Bfield = 0;
   RRR50->SetParameters(RRR,Bfield);
   RRR50->SetLineColor(4);
   sprintf (string,"RRR=%.0f B=%0.f T\n",RRR,Bfield);
   leg1->AddEntry(RRR50,string,"l");

   TF1 *RRR50_B2 = Rcu->DrawCopy("sameC");
   RRR = 50;
   Bfield = 2;
   RRR50_B2->SetParameters(RRR,Bfield);
   RRR50_B2->SetLineColor(4);
   RRR50_B2->SetLineStyle(2);
   sprintf (string,"RRR=%.0f B=%0.f T\n",RRR,Bfield);
   leg1->AddEntry(RRR50_B2,string,"l");


   /*TF1 *RRR42 = Rcu->DrawCopy("sameC");
   RRR = 42;
   Bfield = 0;
   RRR42->SetParameters(RRR,Bfield);
   RRR42->SetLineColor(2);
   sprintf (string,"RRR=%.0f B=%0.f T\n",RRR,Bfield);
   leg1->AddEntry(RRR42,string,"l");*/

   TF1 *RRR100 = Rcu->DrawCopy("sameC");
   RRR = 100;
   Bfield = 0;
   RRR100->SetParameters(RRR,Bfield);
   RRR100->SetLineColor(1);
   sprintf (string,"RRR=%.0f B=%0.f T\n",RRR,Bfield);
   leg1->AddEntry(RRR100,string,"l");

   TF1 *RRR100_B2 = Rcu->DrawCopy("sameC");
   RRR = 100;
   Bfield = 2;
   RRR100_B2->SetParameters(RRR,Bfield);
   RRR100_B2->SetLineColor(1);
   RRR100_B2->SetLineStyle(2);
   sprintf (string,"RRR=%.0f B=%0.f T\n",RRR,Bfield);
   leg1->AddEntry(RRR100_B2,string,"l");

   TF1 *RRR200 = Rcu->DrawCopy("sameC");
   RRR = 200;
   Bfield = 0;
   RRR200->SetParameters(RRR,Bfield);
   RRR200->SetLineColor(2);
   sprintf (string,"RRR=%.0f B=%0.f T\n",RRR,Bfield);
   leg1->AddEntry(RRR200,string,"l");

   TF1 *RRR200_B2 = Rcu->DrawCopy("sameC");
   RRR = 200;
   Bfield=2;
   RRR200_B2->SetParameters(RRR,Bfield);
   RRR200_B2->SetLineColor(2);
   RRR200_B2->SetLineStyle(2);
   sprintf (string,"RRR=%.0f B=%0.f T\n",RRR,Bfield);
   leg1->AddEntry(RRR200_B2,string,"l");

   /*TF1 *RRR1000 = Rcu->DrawCopy("sameC");
   RRR = 1000;
   RRR1000->SetParameters(RRR,Bfield);
   RRR1000->SetLineColor(2);
   sprintf (string,"RRR=%.0f B=%0.f T\n",RRR,Bfield);
   leg1->AddEntry(RRR1000,string,"l");*/
   

   sprintf(filename,"thermal_properties_c1.eps");
   c1->SaveAs(filename);
   sprintf(filename,"thermal_properties_c1.png");
   c1->SaveAs(filename);
   /*sprintf(filename,"thermal_properties_c2.pdf");
   c2->SaveAs(filename);
   sprintf(filename,"thermal_properties_c2.png");
   c2->SaveAs(filename);
   sprintf(filename,"thermal_properties_c3.pdf");
   c3->SaveAs(filename);
   sprintf(filename,"thermal_properties_c3.png");
   c3->SaveAs(filename);*/
   sprintf(filename,"thermal_properties_c4.eps");
   c4->SaveAs(filename);
   sprintf(filename,"thermal_properties_c4.png");
   c4->SaveAs(filename);
   sprintf(filename,"thermal_properties_c5.eps");
   c5->SaveAs(filename);
   sprintf(filename,"thermal_properties_c5.png");
   c5->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;
}
Double_t fit_func (Double_t *x, Double_t *par) {

   // Fit ln(thermal conductivity) as a funciton of ln(Temperature)

   Double_t p0=par[0];
   Double_t p1=par[1];
   Double_t p2=par[2];
   Double_t p3=par[3];
   Double_t p4=par[4];

   Double_t x1=log(x[0]);
   Double_t pi=3.14159;

   char string[256];

   Double_t func;

   func = p0 + p1*x1 + p2*x1*x1 + p3*x1*x1*x1 + p4*x1*x1*x1*x1;

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

   if (func < 50) {
       return exp(func);
      }
   else {
      return exp(50);
     }
}
Double_t fit_rho_func (Double_t *x, Double_t *par) {

   // Fit ln(electrical conductivity) as a funciton of ln(Temperature)
   // Use Fermi function to fit resistivity vs temperature

   // Note that function fits data best for RRR ~ 20 - 200, as the scaled function is scaled from a RRR=42, but the data dependence
  //           change slightly as a function of RRR.


   Double_t rho=4.202;
   Double_t slope=1.61;
   Double_t cc2=-0.6352;
   Double_t cc3=0.4155;
   Double_t scale=1;
   Double_t mean=log(75);
   Double_t RRR=par[6];

   /*Double_t rho=par[0];
   Double_t slope=par[1];
   Double_t cc2=par[2];
   Double_t cc3=par[3];
   Double_t scale=par[4];
   Double_t mean=par[5];
   Double_t RRR=par[6];*/

   Double_t Rcu_300=1.545;    // in units of 10^-8 Ohm-m  // Note that RRR should be referenced to T=273K. (not 300 K).

   Double_t x1=log(x[0]);
   Double_t pi=3.14159;
   /*Double_t RRR850= 850;  // use this RRR value as reference.
   Double_t logRRR850 = log(Rcu_300/RRR850);*/
   Double_t RRR42= 42;  // use this RRR value as reference.
   Double_t logRRR42 = log(Rcu_300/RRR42);

   char string[256];

   Double_t func;
   Double_t arg;

   // arg = -(slope*(x1-mean)+cc2_scale*(x1-mean)*(x1-mean)+cc3*(x1-mean)*(x1-mean)*(x1-mean)+cc4*(x1-mean)*(x1-mean)*(x1-mean)*(x1-mean) );
   arg = -(slope*(x1-mean)+cc2*(x1-mean)*(x1-mean)+cc3*(x1-mean)*(x1-mean)*(x1-mean) );

    if (arg < -100) {
       func = rho;
       }
   else if (arg > 100) {
      func = logRRR42;
      }
   else  {
      func = logRRR42 + rho/(1+exp(arg));
      }

   // reference to room temperature

   scale = log(RRR)/log(RRR42);
   // printf ("x[0]=%f RRR %f RRR42 %f scale=%f \n",x[0],RRR,RRR42,scale);
   Double_t scaled = exp(log(Rcu_300)-scale*(log(Rcu_300)-func));
   return scaled;
}
Double_t cu_resistance (Double_t *x, Double_t *par) {

   // x[0] is the temperature, cu_resistance is in 10^-8 Ohm-m.
   // par[0] is the RRR
   // par[1] is the B field in Tesla
   // Fit ln(electrical resistivity) as a funciton of ln(Temperature) for copper
   // Use Fermi function to fit resistivity vs temperature

   // Note that function fits data best for RRR ~ 20 - 200, as the scaled function is scaled from a RRR=42, but the data dependence
  //           change slightly as a function of RRR.

   Double_t rho=4.202;
   Double_t slope=1.61;
   Double_t cc2=-0.6352;
   Double_t cc3=0.4155;
   Double_t scale=1;
   Double_t mean=log(75);
   Double_t RRR=par[0]*1.72/1.545;         //Note that RRR should be referenced to T=273K. (not 300 K). Adjust for difference
   Double_t B[1];
   B[0] = par[1];   // pass on the value of the B-field.

   // scale RRR by the B field factor.

  Double_t deltaB = cu_magnetoresistivity (B,  par);
  RRR = RRR/(1 + deltaB);

   Double_t Rcu_300=1.72;    // in units of 10^-8 Ohm-m
   // Double_t Rcu_300=1.545;    // in units of 10^-8 Ohm-m.  Note that RRR should be referenced to T=273K. (not 300 K).

   Double_t x1=log(x[0]);
   Double_t pi=3.14159;
   /*Double_t RRR850= 850;  // use this RRR value as reference.
   Double_t logRRR850 = log(Rcu_300/RRR850);*/
   Double_t RRR42= 42;  // use this RRR value as reference.
   Double_t logRRR42 = log(Rcu_300/RRR42);

   char string[256];

   Double_t func;
   Double_t arg;

   arg = -(slope*(x1-mean)+cc2*(x1-mean)*(x1-mean)+cc3*(x1-mean)*(x1-mean)*(x1-mean) );

    if (arg < -100) {
       func = rho;
       }
   else if (arg > 100) {
      func = logRRR42;
      }
   else  {
      func = logRRR42 + rho/(1+exp(arg));
      }

   // reference to room temperature

   scale = log(RRR)/log(RRR42);
   // printf ("x[0]=%f RRR %f RRR42 %f scale=%f \n",x[0],RRR,RRR42,scale);
   Double_t scaled = exp(log(Rcu_300)-scale*(log(Rcu_300)-func));

   // printf("x[0] %f B[0] %f par[0] %f delta B %f\n",x[0],B[0],par[0],deltaB);

   return scaled;
}
Double_t transfer_func (Double_t *x, Double_t *par) {

   // parameterization of cooling over a large temperature range due to nucleate and film boiling.
   // note that ACconst = heat transfer q (W/cm2)*deltax/(Delta T) 

   Double_t ACmag=par[0];
   Double_t ACslope=par[1];
   Double_t option=par[2];   // option=1 upper curve, option=0 lower curve

   Double_t x1=x[0];    // Note that input temperature should be given relative to the helium bath 
   Double_t pi=3.14159;
   Double_t Tbath=4.2;
   Double_t T1=1;   // degrees relative to bath temperature
   Double_t T2=10;   // degrees relative to bath temperature
   Double_t deltax=0.533; // width of conductor (cm)

   char string[256];

   if (x1 < 0.001) {
      return 0;
      }

   Double_t ACconst;
   Double_t factor;
   if (option) {
      factor = 1;
      }
   else if {
     factor=10;
     }

     if (x1 <  T1/factor) {
       ACconst = ACmag*(1-ACslope*T1);
       }
   else if (x1 < T2/factor) {
       ACconst = ACmag*(1-ACslope*T1)/(x1*factor);
       }
   else {
       ACconst = ACmag*(1-ACslope*T2);
       }

   Double_t q = ACconst*x1/deltax;    // note that ACconst * Delta T = heat transfer q (W/cm)/deltax = q(W/cm2), p. 202 of Helium Cryogenics by Van Sciver

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

   return q;
}
Double_t cu_magnetoresistivity (Double_t *x, Double_t *par) {

   // Parameterization of the magneto resistance of Cu from Benz Journal of Applied Physics vol 40, number 1 April 1969, p. 2003

   Double_t S0=par[0];  // Resistivity ratio at H=0;

   // Double_t x1=x[0]*1e-6;  // independent variable is H*S0 x 10^-6 G Eq. (3).
   Double_t x1=x[0]*S0*1e-2;  // independent variable is H*S0 x 10^-6 G Eq. (3), Use H in Tesla (10^4 G)
   Double_t pi=3.14159;
   Double_t R0=1.72/S0;     // Units are 10^-8 Ohm-m; Resistivity at 4.2 deg K.

   Double_t k1=2.368e-4;    // Constants in polynomial fit to data X<2
   Double_t k2=7.389e-2;
   Double_t k3=6.505e-2;
   Double_t k4= -7.121e-3;

   Double_t c1=-0.1649;    // Constants in polynomial fit to data X>2
   Double_t c2=0.2603;
   Double_t c3=-1.129e-3;
   Double_t c4= 4.011e-6;

   char string[256];

   Double_t func=0;

   if (x1 < 2) {
       func = k1 + k2*x1 + k3*x1*x1 + k4*x1*x1*x1;
      }
  else if (x1 < 130) {
       func = c1 + c2*x1 + c3*x1*x1 + c4*x1*x1*x1;
       }

   // Return value of resistivity instead of fractional change from H=0;

  // func = R0*(1+func);  // units are 10^-8 Ohm-m;

   /*sprintf (string,"x[0]=%f, x1=%f S0=%f, R0=%f, func=%f\n",x[0],x1,S0, R0,func);
   printf ("string=%s",string);*/

   return func;
}