AliRoot Core  3dc7879 (3dc7879)
LandauTest.C
Go to the documentation of this file.
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 
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
16 
17  //void Anal();
18 
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 
40 ClassImp(TLandauMean)
42 
43 void TLandauMean::Init(Int_t n, Float_t mean, Float_t sigma)
44 {
46 
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 
75 {
77 
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 {
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 {
202 
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 
printf("Chi2/npoints = %f\n", TMath::Sqrt(chi2/npoints))
void Init(Int_t n, Float_t mean, Float_t sigma)
Definition: LandauTest.C:43
Float_t fLM_0_7[3]
Definition: LandauTest.C:29
Float_t fTM_0_7[3]
Definition: LandauTest.C:24
Double_t landau(Double_t *xp, Double_t *pp)
Definition: DrawESD.C:38
Float_t fLM_0_8[3]
Definition: LandauTest.C:30
Float_t fLSigma
Definition: LandauTest.C:21
Float_t fMedian3
Definition: LandauTest.C:33
void Gener()
Definition: LandauTest.C:74
Float_t fLM_0_6[3]
Definition: LandauTest.C:28
TH1F * LandauTest(Float_t meano, Float_t sigma, Float_t meanlog0, Int_t n, Float_t ratio)
Definition: LandauTest.C:195
void sum()
TTree * tree
Float_t fLMean
Definition: LandauTest.C:20
TF1 * f
Definition: interpolTest.C:21
Float_t fTM_0_10[3]
Definition: LandauTest.C:26
void GenerLandau(Int_t nsamples)
Definition: LandauTest.C:171
Float_t fTM_0_6[3]
Definition: LandauTest.C:23
Int_t fNSample
Definition: LandauTest.C:19
Float_t fLM_0_10[3]
Definition: LandauTest.C:31
Float_t fTM_0_8[3]
Definition: LandauTest.C:25
Float_t Moment3(Float_t sum1, Float_t sum2, Float_t sum3, Int_t n, Float_t m[3])
Definition: LandauTest.C:52