Fix in AliTrackerBase::PropagateTo... routines: length integral was not correctly...
[u/mrichter/AliRoot.git] / PWGUD / 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 }