ATO-98 - connect distortion trees - with custom description ()
[u/mrichter/AliRoot.git] / TPC / LandauTest.C
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
16   // void Anal();
17
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 {
42   /// init parameters
43
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 {
73   /// generate sample
74
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