doxy: all TPC root converted
[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 /// \file LandauTest.C
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
39 /// \cond CLASSIMP
40 ClassImp(TLandauMean)
41 /// \endcond
42
43 void TLandauMean::Init(Int_t n, Float_t mean, Float_t sigma)
44 {
45   /// init parameters
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
74 void TLandauMean::Gener()
75 {
76   /// generate sample
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 {
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
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
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