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