Adding simple example to load default debug streamer
[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
10
11class TLandauMean: public TObject {
12public:
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
32private:
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
37ClassImp(TLandauMean)
38
39void 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
48Float_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
70void 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
167void 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
191TH1F * 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