]>
Commit | Line | Data |
---|---|---|
2803ac99 | 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 | } |