]>
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 | ||
3a4edebe | 10 | /// \file LandauTest.C |
1c53abe2 | 11 | |
12 | class TLandauMean: public TObject { | |
13 | public: | |
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 | |
34 | private: | |
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 | 40 | ClassImp(TLandauMean) |
3a4edebe | 41 | /// \endcond |
1c53abe2 | 42 | |
43 | void 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 | ||
52 | Float_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 | ||
74 | void 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 | ||
171 | void 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 | ||
195 | TH1F * 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 |