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