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

// #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.1,2,5};
   Double_t dummyy[npts]={0.0001,0.0001,200};

   //
   TCanvas *c1 = new TCanvas("c1","c1 magnetoresistivity",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 xmin=0;
   Double_t xmax=2;
   Double_t ymin=0.001;
   Double_t ymax=10;

   // dummy to draw axes

   TGraph *dummy = new TGraph (npts,dummyx,dummyy);
   TLegend *leg = new TLegend(0.2,0.7,0.6,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("B (T)");
   dummy->GetYaxis()->SetTitle("#rho (10^{-8} #Omega-m)");
   dummy->Draw("Ap");
   dummy->GetXaxis()->SetNdivisions(505);

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

   TF1 *magnetoR = new TF1("magnetoresistivity_func",magnetoresistivity_func,xmin,xmax,1);
   sprintf(string,"Benz J App Phys 40 (1969) 2003\n");
   t1 = new TLatex(0.20,0.65,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.025);
   t1->Draw(); 

   Double_t S0=40;
   magnetoR->SetParameter(0,S0);
   TF1 *magnetoR1 = magnetoR->DrawCopy("sameC");
   sprintf(string,"Resistivity ratio=%f\n",S0);
   leg->AddEntry(magnetoR1,string,"l");
   magnetoR1->SetLineColor(2);

   Double_t S0=75;
   magnetoR->SetParameter(0,S0);
   TF1 *magnetoR2 = magnetoR->DrawCopy("sameC");
   sprintf(string,"Resistivity ratio=%f\n",S0);
   leg->AddEntry(magnetoR2,string,"l");
   magnetoR2->SetLineColor(4);

   Double_t S0=200;
   magnetoR->SetParameter(0,S0);
   TF1 *magnetoR3 = magnetoR->DrawCopy("sameC");
   sprintf(string,"Resistivity ratio=%f\n",S0);
   leg->AddEntry(magnetoR3,string,"l");
   magnetoR3->SetLineColor(1);

   Double_t S0=1000;
   magnetoR->SetParameter(0,S0);
   TF1 *magnetoR4 = magnetoR->DrawCopy("sameC");
   sprintf(string,"Resistivity ratio=%f\n",S0);
   leg->AddEntry(magnetoR4,string,"l");
   magnetoR4->SetLineColor(3);


   leg->Draw();

  
   sprintf(filename,"magnetoresistivity_c1.pdf");
   c1->SaveAs(filename);
   sprintf(filename,"magnetoresistivity_c1.png");
   c1->SaveAs(filename);

}
Double_t magnetoresistivity_func (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 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;
}