Updates to Trains. create a job-script to help
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / IdentifiedHighPt / macros / my_functions.C
CommitLineData
4ebdd20e 1#include <TMath.h>
2#include <TH2.h>
3#include <TProfile.h>
4
5#include <iostream>
6
7using namespace std;
8
9const Double_t eMass = 0.000511; //electron mass
10const Double_t piMass = 0.13957; //pion mass
11const Double_t kMass = 0.493676999999999977; //kaon mass
12const Double_t pMass = 0.938271999999999995; //proton mass
13
14
15Double_t fitf3G(Double_t* xx, Double_t* par);
16Double_t FitFunc(Double_t* x, Double_t *par);
17Double_t SigmaFunc(Double_t* x, Double_t *par);
18
19TH2D* hDeDxVsP = 0;
20TProfile* hMeanP = 0;
21Double_t fixMIP = 50.0;
22Double_t fixPlateau = 75.0;
23
24//______________________________________________________________________________
25Double_t fitf3G(Double_t* xx, Double_t* par)
26{
27 //
28 // Could speed up fit by forcing it to use <p>. In that way the parameters
29 // could be amde statis cand only changed when going to a new p bin
30 //
31 Double_t p = xx[0];
32 Double_t dedx = xx[1];
33
34 const Int_t bin = hDeDxVsP->GetXaxis()->FindBin(p);
35
36 if(hMeanP) {
37 // cout << "p before: " << p;
38 p = hMeanP->GetBinContent(bin);
39 // cout << ", p after: " << p << endl;
40 }
41 const Int_t binStart = Int_t(par[0]);
42 const Int_t binStop = Int_t(par[1]);
43
44 if(bin<binStart || bin>binStop) {
45
46 cout << "Error: bin " << bin << " not inside inteval [" << binStart
47 << "; " << binStop << "]" << endl;
48 return 0;
49 }
50
51 const Int_t nParDeDx = Int_t(par[2]);
52
53 Double_t* parDeDx = &par[3];
54
55 Int_t offset = 4 + nParDeDx; // binStart + binStop + nParDeDx + optionDedx + nParDeDx parameters
56
57 const Int_t nParSigma = Int_t(par[offset]);
58 offset += 1; // nParSigma
59
60 Double_t* parSigma = &par[offset];
61 offset += 1 + nParSigma; // optionSigma + nParSigma parameters
62
63 Double_t piMean = FitFunc(&p, parDeDx);
64 Double_t pKeff = p*piMass/kMass; // corresponding p of a pion with same dE/dx
65 Double_t kMean = FitFunc(&pKeff, parDeDx);
66 Double_t pPeff = p*piMass/pMass; // corresponding p of a pion with same dE/dx
67 Double_t pMean = FitFunc(&pPeff, parDeDx);
68
69 const Double_t piSigma = SigmaFunc(&piMean, parSigma);
70 const Double_t kSigma = SigmaFunc(&kMean, parSigma);
71 const Double_t pSigma = SigmaFunc(&pMean, parSigma);
72
73
74 const Int_t j = bin - binStart;
75 const Double_t piYield = par[j * 3 + offset + 0];
76 const Double_t kYield = par[j * 3 + offset + 1];
77 const Double_t pYield = par[j * 3 + offset + 2];
78
79 return piYield* TMath::Gaus(dedx, piMean, piSigma, kTRUE)
80 + kYield * TMath::Gaus(dedx, kMean, kSigma, kTRUE)
81 + pYield * TMath::Gaus(dedx, pMean, pSigma, kTRUE);
82}
83
84
85//______________________________________________________________________________
86Double_t FitFunc(Double_t* x, Double_t *par)
87{
88 static const Double_t bgMIP = 0.5/piMass;
89 static const Double_t beta2MIP = bgMIP*bgMIP / (1.0+bgMIP*bgMIP);
90 // static const Double_t betapowMIP = TMath;
91 static const Double_t logMIP = TMath::Log(1+bgMIP);
92
93 Int_t option = Int_t(par[0]);
94 Int_t specie = option;
95 option = option%10;
96 specie -= option;
97 specie /= 10;
98
99
100 Double_t bg = 0;
101 switch (specie) {
102
103 case 0: // pion
104 bg = x[0]/piMass;
105 break;
106 case 1: // kaon
107 bg = x[0]/kMass;
108 break;
109 case 2: // proton
110 bg = x[0]/pMass;
111 break;
112 case 3: // electron
113 bg = x[0]/eMass;
114 break;
115 default:
116 cout << "Error in FitFunc: specie " << specie << " not supported!!!!!" << endl;
117 return 0;
118 break;
119 }
120
121 if(bg > 10000.0)
122 bg = 10000.0;
123
124 const Double_t beta2 = bg*bg / (1.0+bg*bg);
125
126 switch (option) {
127
128 case 1: // standard parametrisation
129 {
130 /*
131 c0/beta^2 + c1 * log (1+x)
132 */
133 const Double_t c0 = par[1];
134 const Double_t c1 = par[2];
135
136 const Double_t value = c0/beta2 + c1*TMath::Log(1+bg);
137 return value;
138 }
139 break;
140 case 2: // fix the dE/dx to 50 at 0.5 GeV/c
141 {
142 const Double_t c1 = par[1];
143 const Double_t c0 = (fixMIP-par[1]*logMIP) * beta2MIP;
144
145 const Double_t value = c0/beta2 + c1*TMath::Log(1+bg);
146 return value;
147 }
148 break;
149 case 3: // fix the dE/dx to 50 at 0.5 GeV/c and the plateau to 75
150 {
151 /*
152 a/beta^2 + b/c*log( (1+x)^c / (1 + d*(1+x)^c) )
153
154 Assymptotic behavior:
155
156 1) Small bg (and d small so that d*(1+x)^c << 1)
157
158 a/beta^2 + b * log (1+x)
159
160 So this is the same beavior as the standard expression.
161
162 2) Large bg where d*(1+x)^c >> 1
163 a - b/c*log(d) = plateau
164 -> d = exp(c*(a-plateau)/b)
165
166 */
167 const Double_t b = par[1];
168 const Double_t a = (fixMIP-par[1]*logMIP) * beta2MIP;
169 const Double_t c = par[2];
170 const Double_t d = TMath::Exp(c*(a-fixPlateau)/b);
171
172 // cout << bg << ": " << a << ", " << b << ", " << c << ", " << d << endl;
173
174 const Double_t powbg = TMath::Power(1.0+bg, c);
175
176 const Double_t value = a/beta2 + b/c*TMath::Log(powbg/(1.0 + d*powbg));
177 return value;
178 }
179 break;
180 case 4: // fix the dE/dx to 50 at 0.5 GeV/c and the plateau to 75
181 {
182 /*
183 a/beta^2 + b/c*log( (1+x)^c / (1 + d*(1+x)^c) )
184
185 Assymptotic behavior:
186
187 1) Small bg (and d small so that d*(1+x)^c << 1)
188
189 a/beta^2 + b * log (1+x)
190
191 So this is the same beavior as the standard expression.
192
193 2) Large bg where d*(1+x)^c >> 1
194 a - b/c*log(d) = plateau
195 -> d = exp(c*(a-plateau)/b)
196
197 */
198 const Double_t a = par[1];
199 const Double_t b = par[2];
200 const Double_t c = par[3];
201 const Double_t d = TMath::Exp(c*(a-fixPlateau)/b);
202
203 // cout << bg << ": " << a << ", " << b << ", " << c << ", " << d << endl;
204
205 const Double_t powbg = TMath::Power(1.0+bg, c);
206
207 const Double_t value = a/beta2 + b/c*TMath::Log(powbg/(1.0 + d*powbg));
208 return value;
209 }
210 break;
211 // case 3: // fix the dE/dx to 50 at 0.5 GeV/c + powerlaw
212
213 // static const bgMIP = 0.5/piMass;
214 // static const beta2MIP = bgMIP*bgMIP / (1.0+bgMIP*bgMIP);
215 // static const logMIP = TMath::Log(1+bgMIP);
216
217 // const Double_t c1 = par[1];
218 // const Double_t c2 = par[2]; // expect it to be 0.75 from Bichsel (beta**-1.5 instead of -2)
219 // const Double_t c0 = (50.0-par[1]*logMIP) * beta2MIP;
220
221 // const Double_t value = TMathh::Power(c0/beta2, c2) + c1*TMath::Log(1+bg);
222 // return value;
223
224 default:
225 break;
226 }
227
228 cout << "Error in FitFunc: option " << option << " not supported!!!!!" << endl;
229 return 0;
230}
231
232//______________________________________________________________________________
233Double_t SigmaFunc(Double_t* x, Double_t *par)
234{
235 Int_t option = Int_t(par[0]);
236
237 switch (option) {
238 case 1: // fixed sigma
239 return par[1];
240 break;
241 case 2: // relative sigma
242 return par[1]*x[0];
243 break;
244 case 3: // relative sigma + extrapolation
245 return (par[1] + (x[0]-fixMIP)*par[2])*x[0];
246 break;
247 case 4: // relative sigma with dE/dx to some power close to 1
248 return par[1]*TMath::Power(x[0], par[2]);
249 break;
250 case 5: // relative sigma with dE/dx to some power close to 1
251 return TMath::Sqrt(par[1]*par[1]*x[0]*x[0] + par[2]*par[2]);
252 break;
253 default:
254 break;
255 }
256
257 cout << "Error in SigmaFunc: option " << option << " not supported!!!!!" << endl;
258 return 0;
259}