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