[u/mrichter/AliRoot.git] / TPC / LandauTest.C

#include "TH1.h"
#include "TH2.h"
#include "TFile.h"
#include "TTree.h"

#include "TRandom.h"
#include "TPad.h"
#include "TCanvas.h"

/// \file LandauTest.C

class TLandauMean: public TObject {
public:
  void Init(Int_t n, Float_t mean, Float_t sigma);  // initial parameters
  void Gener();          // gener sample 

  //void Anal();

  Int_t fNSample;      // number of samples
  Float_t fLMean;        // landau mean
  Float_t fLSigma;       // landau sigma
  //
  Float_t fTM_0_6[3];    // truncated method  - first 3 momenta
  Float_t fTM_0_7[3];    // truncated method  - first 3 momenta
  Float_t fTM_0_8[3];    // truncated method  - first 3 momenta
  Float_t fTM_0_10[3];   // truncated method  - first 3 momenta
  //
  Float_t fLM_0_6[3];    // truncated log.  method  - first 3 momenta
  Float_t fLM_0_7[3];    // truncated log.  method  - first 3 momenta
  Float_t fLM_0_8[3];    // truncated log.  method  - first 3 momenta
  Float_t fLM_0_10[3];   // truncated log.  method  - first 3 momenta

  Float_t fMedian3;      // median 3 value
private:
  Float_t Moment3(Float_t sum1, Float_t sum2, Float_t sum3, Int_t n, Float_t m[3]);
  ClassDef(TLandauMean,1)
};

/// \cond CLASSIMP
ClassImp(TLandauMean)
/// \endcond

void TLandauMean::Init(Int_t n, Float_t mean, Float_t sigma)
{
  /// init parameters

  fNSample = n;
  fLMean   = mean;
  fLSigma  = sigma;
}

Float_t TLandauMean::Moment3(Float_t sumi1, Float_t sumi2, Float_t sumi3, Int_t sum, Float_t m[3])
{
  Float_t m3=0;

  //  m3 = (sumi3-3*pos*sumi2+3*pos*pos*sumi-pos*pos*pos*sum)/sum; 
  Float_t pos = sumi1/sum;
  m[0] = pos;
  m[1] = sumi2/sum-pos*pos;
  if (m[1]<=0){
    printf("pici pici\n");
  }
  else
    m[1] = TMath::Sqrt(m[1]);
  m3 = (sumi3-3*pos*sumi2+3*pos*pos*sumi1-pos*pos*pos*sum)/sum; 
  Float_t sign = m3/TMath::Abs(m3);
  m3 = TMath::Power(sign*m3,1/3.);
  m3*=sign;

  m[2] = m3;  
  return m3;
}

void TLandauMean::Gener()
{
  /// generate sample

  Float_t * buffer = new Float_t[fNSample];
  
  for (Int_t i=0;i<fNSample;i++) {
    buffer[i] = gRandom->Landau(fLMean,fLSigma); 
    if (buffer[i]>1000) buffer[i]=1000;
  }

  Int_t *index = new Int_t[fNSample];
  TMath::Sort(fNSample,buffer,index,kFALSE);

  //
  Float_t median = buffer[index[fNSample/3]];
  //
  Float_t sum06[4]  = {0.,0.,0.,0.};
  Float_t sum07[4]  = {0.,0.,0.,0.};
  Float_t sum08[4]  = {0.,0.,0.,0.};
  Float_t sum010[4]  = {0.,0.,0.,0.};
  //
  Float_t suml06[4] = {0.,0.,0.,0.};
  Float_t suml07[4] = {0.,0.,0.,0.};
  Float_t suml08[4] = {0.,0.,0.,0.};
  Float_t suml010[4] = {0.,0.,0.,0.};
  //

  for (Int_t i =0; i<fNSample; i++){
    Float_t amp  = buffer[index[i]];
    Float_t lamp = median*TMath::Log(1.+amp/median);
    if (i<0.6*fNSample){
      sum06[0]+= amp;
      sum06[1]+= amp*amp;
      sum06[2]+= amp*amp*amp;
      sum06[3]++;
      suml06[0]+= lamp;
      suml06[1]+= lamp*lamp;
      suml06[2]+= lamp*lamp*lamp;
      suml06[3]++;
    }

    if (i<0.7*fNSample){
      sum07[0]+= amp;
      sum07[1]+= amp*amp;
      sum07[2]+= amp*amp*amp;
      sum07[3]++;
      suml07[0]+= lamp;
      suml07[1]+= lamp*lamp;
      suml07[2]+= lamp*lamp*lamp;
      suml07[3]++;
    }
    if (i<0.8*fNSample){
      sum08[0]+= amp;
      sum08[1]+= amp*amp;
      sum08[2]+= amp*amp*amp;
      sum08[3]++;
      suml08[0]+= lamp;
      suml08[1]+= lamp*lamp;
      suml08[2]+= lamp*lamp*lamp;
      suml08[3]++;
    }
    if (i<1*fNSample){
      sum010[0]+= amp;
      sum010[1]+= amp*amp;
      sum010[2]+= amp*amp*amp;
      sum010[3]++;
      suml010[0]+= lamp;
      suml010[1]+= lamp*lamp;
      suml010[2]+= lamp*lamp*lamp;
      suml010[3]++;

    }
  }
  //  
  fMedian3 = median;
  //
  Moment3(sum06[0],sum06[1],sum06[2],sum06[3],fTM_0_6);  
  Moment3(sum07[0],sum07[1],sum07[2],sum07[3],fTM_0_7);  
  Moment3(sum08[0],sum08[1],sum08[2],sum08[3],fTM_0_8);  
  Moment3(sum010[0],sum010[1],sum010[2],sum010[3],fTM_0_10);  
  //

  Moment3(suml06[0],suml06[1],suml06[2],suml06[3],fLM_0_6);  
  Moment3(suml07[0],suml07[1],suml07[2],suml07[3],fLM_0_7);  
  Moment3(suml08[0],suml08[1],suml08[2],suml08[3],fLM_0_8);  
  Moment3(suml010[0],suml010[1],suml010[2],suml010[3],fLM_0_10);  
  //
  fLM_0_6[0] = (TMath::Exp(fLM_0_6[0]/median)-1.)*median;
  fLM_0_7[0] = (TMath::Exp(fLM_0_7[0]/median)-1.)*median;
  fLM_0_8[0] = (TMath::Exp(fLM_0_8[0]/median)-1.)*median;
  fLM_0_10[0] = (TMath::Exp(fLM_0_10[0]/median)-1.)*median;
  //
  delete [] buffer;
}   


void GenerLandau(Int_t nsamples)
{
  TLandauMean * landau = new TLandauMean;
  TFile f("Landau.root","recreate");
  TTree * tree = new TTree("Landau","Landau");
  tree->Branch("Landau","TLandauMean",&landau); 
  
  for (Int_t i=0;i<nsamples;i++){
    Int_t   n     = 20 + Int_t(gRandom->Rndm()*150);
    Float_t mean  = 40. +gRandom->Rndm()*50.;
    Float_t sigma = 5.  +gRandom->Rndm()*15.;
    landau->Init(n, mean, sigma);
    landau->Gener();
    tree->Fill();
  }
  tree->Write();
  f.Close();

}


TH1F *  LandauTest(Float_t meano,  Float_t sigma, Float_t meanlog0, Int_t n,Float_t ratio)
{ 
  /// test for different approach of de dx resolution
  /// meano, sigma  - mean value of Landau distribution and sigma
  /// meanlog0      - scaling factor for logarithmic mean value
  /// n             - number of used layers
  /// ratio         - ratio of used amplitudes for truncated mean


  TCanvas * pad = new TCanvas("Landau test");
  pad->Divide(2,2);
  TH1F * h1 = new TH1F("h1","Logarithmic mean",300,0,4*meano);
  TH1F * h2 = new TH1F("h2","Logarithmic amplitudes",300,0,8*meano);
  TH1F * h3 = new TH1F("h3","Mean",300,0,4*meano);
  TH1F * h4 = new TH1F("h4","Amplitudes",300,0,8*meano);

  for(Int_t j=0;j<10000;j++){
    //generate sample and sort it
    Float_t * buffer = new Float_t[n];
    Float_t * buffer2= new Float_t[n];
    
    for (Int_t i=0;i<n;i++) {
      buffer[i] = gRandom->Landau(meano,sigma); 
      buffer2[i] = buffer[i];    
    }
    //add crosstalk
    for (Int_t i=1;i<n-1;i++) {
      buffer[i] =    buffer2[i]*1.0+ buffer2[i-1]*0.0+ buffer2[i+1]*0.0;
      buffer[i] = TMath::Min(buffer[i],1000.); 
    }
    Int_t *index = new Int_t[n];
    TMath::Sort(n,buffer,index,kFALSE);

    //calculate mean
    Float_t sum;
    sum=0;
    Float_t mean;
    Float_t used = 0;
    for (Int_t i=0;i<n*ratio;i++) {
      if (buffer[index[i]]<1000.){
	Float_t amp = meanlog0*TMath::Log(1+buffer[index[i]]/meanlog0);
	sum += amp;            
	used++;
      }
    }
    mean = sum/used;
    //
    sum=0;
    used=0;
    Float_t sum2=0;
    Float_t meanlog =meanlog0;
    for (Int_t i=0;i<n*ratio;i++) {
      if (buffer[index[i]]<1000.){
	Float_t amp = meanlog*TMath::Log(1.+buffer[index[i]]/(meanlog));
	sum +=amp;
	sum2+=buffer[index[i]];
	used++;
	h2->Fill(amp);
	h4->Fill(buffer[index[i]]);
      }
    }
    mean = sum/used;
    mean = (TMath::Exp(mean/meanlog)-1)*meanlog;
    Float_t mean2 = sum2/used;
    //mean2 = (mean+mean2)/2.;
    h1->Fill(mean);    
    h3->Fill(mean2);
  }

  pad->cd(1);
  h1->Draw();
  pad->cd(2);
  h2->Draw();
  pad->cd(3);
  h3->Draw();
  pad->cd(4);
  h4->Draw();


  return h1;

}
Commit	Line	Data
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
3a4edebe	10	/// \file LandauTest.C
1c53abe2	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
a7a1dd76	16
3a4edebe	17	//void Anal();
a7a1dd76	18
1c53abe2	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
3a4edebe	39	/// \cond CLASSIMP
1c53abe2	40	ClassImp(TLandauMean)
3a4edebe	41	/// \endcond
1c53abe2	42
	43	void TLandauMean::Init(Int_t n, Float_t mean, Float_t sigma)
	44	{
a7a1dd76	45	/// init parameters
a7a1dd76	46
1c53abe2	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-3possumi2+3pospossumi-pospospossum)/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-3possumi2+3pospossumi1-pospospossum)/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	{
a7a1dd76	76	/// generate sample
a7a1dd76	77
1c53abe2	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]+= ampampamp;
	109	sum06[3]++;
	110	suml06[0]+= lamp;
	111	suml06[1]+= lamp*lamp;
	112	suml06[2]+= lamplamplamp;
	113	suml06[3]++;
	114	}
	115
	116	if (i<0.7*fNSample){
	117	sum07[0]+= amp;
	118	sum07[1]+= amp*amp;
	119	sum07[2]+= ampampamp;
	120	sum07[3]++;
	121	suml07[0]+= lamp;
	122	suml07[1]+= lamp*lamp;
	123	suml07[2]+= lamplamplamp;
	124	suml07[3]++;
	125	}
	126	if (i<0.8*fNSample){
	127	sum08[0]+= amp;
	128	sum08[1]+= amp*amp;
	129	sum08[2]+= ampampamp;
	130	sum08[3]++;
	131	suml08[0]+= lamp;
	132	suml08[1]+= lamp*lamp;
	133	suml08[2]+= lamplamplamp;
	134	suml08[3]++;
	135	}
	136	if (i<1*fNSample){
	137	sum010[0]+= amp;
	138	sum010[1]+= amp*amp;
	139	sum010[2]+= ampampamp;
	140	sum010[3]++;
	141	suml010[0]+= lamp;
142	suml010[1]+= lamp*lamp;
143	suml010[2]+= lamplamplamp;
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	{
3a4edebe	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
1c53abe2	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