PWGLF/SPECTRA/ChargedHadrons/dNdPt/macros/plots/divide.C

   1 //------------------------------------------------------------------------------
   2 // divide.C
   3 //
   4 // macro to divide two TGraphs using Exponential interpolation
   5 //------------------------------------------------------------------------------
   6
   7 using namespace std;
   8
   9 Double_t ExtrapolateExp(Double_t x0, Double_t y0, Double_t x1, Double_t y1, Double_t x)
  10 {
  11         // extra/interpolates using an exponential defined by (x0,y0) and (x1,y1) to x
  12
  13         Double_t b = (TMath::Log(y0) - TMath::Log(y1)) / (x0 - x1);
  14         Double_t a = TMath::Log(y0) - b * x0;
  15         Double_t a2 = TMath::Log(y1) - b * x1;
  16
  17         //Printf("%f %f %f", b, a, a2);
  18
  19         return TMath::Exp(a + b * x);
  20 }
  21
  22
  23 TGraphErrors* divide(TGraph* dividend, TGraph* divisor, Bool_t exp = kTRUE)
  24 {
  25
  26 cout << "---------------------------------------------------------" << endl;
  27 cout << "starting macro: << divide.C >> " << endl;
  28 cout << "---------------------------------------------------------" << endl;
  29
  30 Double_t x=0;
  31 Double_t y=0;
  32 Double_t y1;
  33 Double_t y2;
  34 Double_t ex;
  35 Double_t ey;
  36 Double_t e1;
  37 Double_t e2;
  38
  39
  40 divisor->GetPoint(0,x,y);
  41 Double_t minx = x;
  42 cout << "divisor: minx = " << minx << endl;
  43
  44 divisor->GetPoint(divisor->GetN()-1,x,y);
  45 Double_t maxx = x;
  46 cout << "divisor: maxx = " << maxx << endl;
  47
  48 dividend->GetPoint(0,x,y);
  49 if (x > minx) { minx = x; }
  50
  51 dividend->GetPoint(dividend->GetN()-1,x,y);
  52 if (x < maxx) { maxx = x; }
  53
  54 cout << "minx = " << minx << endl;
  55 cout << "maxx = " << maxx << endl;
  56
  57 Int_t bins = 0;
  58
  59 for (int i=0; i < dividend->GetN(); i++) {
  60     dividend->GetPoint(i,x,y);
  61     if ((x <= maxx) && (x >= minx)) {bins++;}
  62 }
  63
  64 TGraphErrors* result = new TGraphErrors(bins);
  65 TGraph* errors = new TGraph(divisor->GetN(),divisor->GetX(),divisor->GetEY());
  66
  67 Int_t j=0;
  68
  69 for (int i=0; i < dividend->GetN(); i++) {
  70     dividend->GetPoint(i,x,y1);
  71     if (!((x <= maxx) && (x >= minx))) { continue; }
  72     if (exp) {
  73         Double_t x1temp = 0;
  74         Double_t y1temp = 0;
  75         for (int n=0; n < divisor->GetN(); n++) {
  76             divisor->GetPoint(n,x1temp,y1temp);
  77
  78             if (x1temp >= x)  break;
  79         }
  80         cout << "x=" << x;
  81         cout << "x1temp=" << x1temp;
  82         if (x1temp == x) { y2 = y1temp; }
  83         if (x1temp > x) {
  84             Double_t x0temp = 0;
  85             Double_t y0temp = 0;
  86             cout << "nfortemp=" << n << endl;
  87             divisor->GetPoint(n-1,x0temp,y0temp);
  88             y2 = ExtrapolateExp(x0temp,y0temp,x1temp,y1temp,x);
  89             cout << "exp extrapolation used" << endl;
  90         }
  91     } else {
  92         y2 = divisor->Eval(x);
  93     }
  94
  95     y = y1/y2;
  96     e1 = dividend->GetErrorY(i);
  97     //e2 = errors->Eval(x);
  98     e2 = 0;
  99     ex = 0;
 100     ey = y*sqrt((e1/y1)*(e1/y1) + (e2/y2)*(e2/y2));
 101     result->SetPoint(j,x,y);
 102     result->SetPointError(j,ex,ey);
 103     j++;
 104     cout << "data point # " << j << endl;
 105     cout << "   x= " << x << " +- " << ex <<endl;
 106     cout << "   y= " << y << " +- " << ey <<endl;
 107 }
 108 delete errors;
 109 errors = 0;
 110
 111
 112
 113 cout << "min:" << minx <<endl;
 114 cout << "max:" << maxx <<endl;
 115
 116 cout << "---------------------------------------------------------" << endl;
 117 cout << "finished macro: << divide.C >> " << endl;
 118 cout << "---------------------------------------------------------" << endl;
 119 return result;
 120 }