doxy: all TPC root converted
[u/mrichter/AliRoot.git] / TPC / LandauTest.C
CommitLineData
1c53abe2 1#include "TH1.h"
2#include "TH2.h"
3#include "TFile.h"
4#include "TTree.h"
5
6#include "TRandom.h"
7#include "TPad.h"
8#include "TCanvas.h"
9
3a4edebe 10/// \file LandauTest.C
1c53abe2 11
12class TLandauMean: public TObject {
13public:
14 void Init(Int_t n, Float_t mean, Float_t sigma); // initial parameters
15 void Gener(); // gener sample
a7a1dd76 16
3a4edebe 17 //void Anal();
a7a1dd76 18
1c53abe2 19 Int_t fNSample; // number of samples
20 Float_t fLMean; // landau mean
21 Float_t fLSigma; // landau sigma
22 //
23 Float_t fTM_0_6[3]; // truncated method - first 3 momenta
24 Float_t fTM_0_7[3]; // truncated method - first 3 momenta
25 Float_t fTM_0_8[3]; // truncated method - first 3 momenta
26 Float_t fTM_0_10[3]; // truncated method - first 3 momenta
27 //
28 Float_t fLM_0_6[3]; // truncated log. method - first 3 momenta
29 Float_t fLM_0_7[3]; // truncated log. method - first 3 momenta
30 Float_t fLM_0_8[3]; // truncated log. method - first 3 momenta
31 Float_t fLM_0_10[3]; // truncated log. method - first 3 momenta
32
33 Float_t fMedian3; // median 3 value
34private:
35 Float_t Moment3(Float_t sum1, Float_t sum2, Float_t sum3, Int_t n, Float_t m[3]);
36 ClassDef(TLandauMean,1)
37};
38
3a4edebe 39/// \cond CLASSIMP
1c53abe2 40ClassImp(TLandauMean)
3a4edebe 41/// \endcond
1c53abe2 42
43void TLandauMean::Init(Int_t n, Float_t mean, Float_t sigma)
44{
a7a1dd76 45 /// init parameters
46
1c53abe2 47 fNSample = n;
48 fLMean = mean;
49 fLSigma = sigma;
50}
51
52Float_t TLandauMean::Moment3(Float_t sumi1, Float_t sumi2, Float_t sumi3, Int_t sum, Float_t m[3])
53{
54 Float_t m3=0;
55
56 // m3 = (sumi3-3*pos*sumi2+3*pos*pos*sumi-pos*pos*pos*sum)/sum;
57 Float_t pos = sumi1/sum;
58 m[0] = pos;
59 m[1] = sumi2/sum-pos*pos;
60 if (m[1]<=0){
61 printf("pici pici\n");
62 }
63 else
64 m[1] = TMath::Sqrt(m[1]);
65 m3 = (sumi3-3*pos*sumi2+3*pos*pos*sumi1-pos*pos*pos*sum)/sum;
66 Float_t sign = m3/TMath::Abs(m3);
67 m3 = TMath::Power(sign*m3,1/3.);
68 m3*=sign;
69
70 m[2] = m3;
71 return m3;
72}
73
74void TLandauMean::Gener()
75{
a7a1dd76 76 /// generate sample
77
1c53abe2 78 Float_t * buffer = new Float_t[fNSample];
79
80 for (Int_t i=0;i<fNSample;i++) {
81 buffer[i] = gRandom->Landau(fLMean,fLSigma);
82 if (buffer[i]>1000) buffer[i]=1000;
83 }
84
85 Int_t *index = new Int_t[fNSample];
86 TMath::Sort(fNSample,buffer,index,kFALSE);
87
88 //
89 Float_t median = buffer[index[fNSample/3]];
90 //
91 Float_t sum06[4] = {0.,0.,0.,0.};
92 Float_t sum07[4] = {0.,0.,0.,0.};
93 Float_t sum08[4] = {0.,0.,0.,0.};
94 Float_t sum010[4] = {0.,0.,0.,0.};
95 //
96 Float_t suml06[4] = {0.,0.,0.,0.};
97 Float_t suml07[4] = {0.,0.,0.,0.};
98 Float_t suml08[4] = {0.,0.,0.,0.};
99 Float_t suml010[4] = {0.,0.,0.,0.};
100 //
101
102 for (Int_t i =0; i<fNSample; i++){
103 Float_t amp = buffer[index[i]];
104 Float_t lamp = median*TMath::Log(1.+amp/median);
105 if (i<0.6*fNSample){
106 sum06[0]+= amp;
107 sum06[1]+= amp*amp;
108 sum06[2]+= amp*amp*amp;
109 sum06[3]++;
110 suml06[0]+= lamp;
111 suml06[1]+= lamp*lamp;
112 suml06[2]+= lamp*lamp*lamp;
113 suml06[3]++;
114 }
115
116 if (i<0.7*fNSample){
117 sum07[0]+= amp;
118 sum07[1]+= amp*amp;
119 sum07[2]+= amp*amp*amp;
120 sum07[3]++;
121 suml07[0]+= lamp;
122 suml07[1]+= lamp*lamp;
123 suml07[2]+= lamp*lamp*lamp;
124 suml07[3]++;
125 }
126 if (i<0.8*fNSample){
127 sum08[0]+= amp;
128 sum08[1]+= amp*amp;
129 sum08[2]+= amp*amp*amp;
130 sum08[3]++;
131 suml08[0]+= lamp;
132 suml08[1]+= lamp*lamp;
133 suml08[2]+= lamp*lamp*lamp;
134 suml08[3]++;
135 }
136 if (i<1*fNSample){
137 sum010[0]+= amp;
138 sum010[1]+= amp*amp;
139 sum010[2]+= amp*amp*amp;
140 sum010[3]++;
141 suml010[0]+= lamp;
142 suml010[1]+= lamp*lamp;
143 suml010[2]+= lamp*lamp*lamp;
144 suml010[3]++;
145
146 }
147 }
148 //
149 fMedian3 = median;
150 //
151 Moment3(sum06[0],sum06[1],sum06[2],sum06[3],fTM_0_6);
152 Moment3(sum07[0],sum07[1],sum07[2],sum07[3],fTM_0_7);
153 Moment3(sum08[0],sum08[1],sum08[2],sum08[3],fTM_0_8);
154 Moment3(sum010[0],sum010[1],sum010[2],sum010[3],fTM_0_10);
155 //
156
157 Moment3(suml06[0],suml06[1],suml06[2],suml06[3],fLM_0_6);
158 Moment3(suml07[0],suml07[1],suml07[2],suml07[3],fLM_0_7);
159 Moment3(suml08[0],suml08[1],suml08[2],suml08[3],fLM_0_8);
160 Moment3(suml010[0],suml010[1],suml010[2],suml010[3],fLM_0_10);
161 //
162 fLM_0_6[0] = (TMath::Exp(fLM_0_6[0]/median)-1.)*median;
163 fLM_0_7[0] = (TMath::Exp(fLM_0_7[0]/median)-1.)*median;
164 fLM_0_8[0] = (TMath::Exp(fLM_0_8[0]/median)-1.)*median;
165 fLM_0_10[0] = (TMath::Exp(fLM_0_10[0]/median)-1.)*median;
166 //
167 delete [] buffer;
168}
169
170
171void GenerLandau(Int_t nsamples)
172{
173 TLandauMean * landau = new TLandauMean;
174 TFile f("Landau.root","recreate");
175 TTree * tree = new TTree("Landau","Landau");
176 tree->Branch("Landau","TLandauMean",&landau);
177
178 for (Int_t i=0;i<nsamples;i++){
179 Int_t n = 20 + Int_t(gRandom->Rndm()*150);
180 Float_t mean = 40. +gRandom->Rndm()*50.;
181 Float_t sigma = 5. +gRandom->Rndm()*15.;
182 landau->Init(n, mean, sigma);
183 landau->Gener();
184 tree->Fill();
185 }
186 tree->Write();
187 f.Close();
188
189}
190
191
192
193
194
195TH1F * LandauTest(Float_t meano, Float_t sigma, Float_t meanlog0, Int_t n,Float_t ratio)
196{
3a4edebe 197 /// test for different approach of de dx resolution
198 /// meano, sigma - mean value of Landau distribution and sigma
199 /// meanlog0 - scaling factor for logarithmic mean value
200 /// n - number of used layers
201 /// ratio - ratio of used amplitudes for truncated mean
202
1c53abe2 203
204 TCanvas * pad = new TCanvas("Landau test");
205 pad->Divide(2,2);
206 TH1F * h1 = new TH1F("h1","Logarithmic mean",300,0,4*meano);
207 TH1F * h2 = new TH1F("h2","Logarithmic amplitudes",300,0,8*meano);
208 TH1F * h3 = new TH1F("h3","Mean",300,0,4*meano);
209 TH1F * h4 = new TH1F("h4","Amplitudes",300,0,8*meano);
210
211 for(Int_t j=0;j<10000;j++){
212 //generate sample and sort it
213 Float_t * buffer = new Float_t[n];
214 Float_t * buffer2= new Float_t[n];
215
216 for (Int_t i=0;i<n;i++) {
217 buffer[i] = gRandom->Landau(meano,sigma);
218 buffer2[i] = buffer[i];
219 }
220 //add crosstalk
221 for (Int_t i=1;i<n-1;i++) {
222 buffer[i] = buffer2[i]*1.0+ buffer2[i-1]*0.0+ buffer2[i+1]*0.0;
223 buffer[i] = TMath::Min(buffer[i],1000.);
224 }
225 Int_t *index = new Int_t[n];
226 TMath::Sort(n,buffer,index,kFALSE);
227
228 //calculate mean
229 Float_t sum;
230 sum=0;
231 Float_t mean;
232 Float_t used = 0;
233 for (Int_t i=0;i<n*ratio;i++) {
234 if (buffer[index[i]]<1000.){
235 Float_t amp = meanlog0*TMath::Log(1+buffer[index[i]]/meanlog0);
236 sum += amp;
237 used++;
238 }
239 }
240 mean = sum/used;
241 //
242 sum=0;
243 used=0;
244 Float_t sum2=0;
245 Float_t meanlog =meanlog0;
246 for (Int_t i=0;i<n*ratio;i++) {
247 if (buffer[index[i]]<1000.){
248 Float_t amp = meanlog*TMath::Log(1.+buffer[index[i]]/(meanlog));
249 sum +=amp;
250 sum2+=buffer[index[i]];
251 used++;
252 h2->Fill(amp);
253 h4->Fill(buffer[index[i]]);
254 }
255 }
256 mean = sum/used;
257 mean = (TMath::Exp(mean/meanlog)-1)*meanlog;
258 Float_t mean2 = sum2/used;
259 //mean2 = (mean+mean2)/2.;
260 h1->Fill(mean);
261 h3->Fill(mean2);
262 }
263
264 pad->cd(1);
265 h1->Draw();
266 pad->cd(2);
267 h2->Draw();
268 pad->cd(3);
269 h3->Draw();
270 pad->cd(4);
271 h4->Draw();
272
273
274 return h1;
275
276}
277