Updates to Trains. create a job-script to help
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / IdentifiedHighPt / macros / my_functions.C
1 #include <TMath.h>
2 #include <TH2.h>
3 #include <TProfile.h>
4
5 #include <iostream>
6
7 using namespace std;
8
9 const Double_t eMass  = 0.000511;               //electron mass
10 const Double_t piMass = 0.13957;               //pion mass
11 const Double_t kMass  = 0.493676999999999977;  //kaon mass
12 const Double_t pMass  = 0.938271999999999995;  //proton mass
13
14
15 Double_t fitf3G(Double_t* xx, Double_t* par);
16 Double_t FitFunc(Double_t* x, Double_t *par);
17 Double_t SigmaFunc(Double_t* x, Double_t *par);
18
19 TH2D* hDeDxVsP = 0;
20 TProfile* hMeanP = 0;
21 Double_t fixMIP     = 50.0;
22 Double_t fixPlateau = 75.0;
23
24 //______________________________________________________________________________
25 Double_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 //______________________________________________________________________________
86 Double_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 //______________________________________________________________________________
233 Double_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 }