finished lecture 4

Browse code finished lecture 4 master
1 parent 366cd20 commit 8f246add76a4f85b5e84aebfba002ae008fd5258 mchrzasz authored on 8 Oct 2016

+Double_t sqr(Double_t x){return x*x;}
+Double_t Camel2(Int_t nDim, Double_t *Xarg){
+  // 2-dimensional distribution for FOAM, normalized to one (within 1e-5)
+  Double_t x=Xarg[0];
+  Double_t y=Xarg[1];
+  Double_t GamSq= sqr(0.100e0);
+  Double_t Dist=exp(-(sqr(x-1./3) +sqr(y-1./3))/GamSq)/GamSq/TMath::Pi();
+  Dist        +=exp(-(sqr(x-2./3) +sqr(y-2./3))/GamSq)/GamSq/TMath::Pi();
+  return 0.5*Dist;
+}// Camel2
+Int_t kanwa(){
+  // TH2D  *hst_xy = new TH2D("hst_xy" ,  "x-y plot", 50,0,1.0, 50,0,1.0);
+  Double_t *MCvect =new Double_t[2]; // 2-dim vector generated in the MC run
+  TRandom3  *PseRan   = new TRandom3();  // Create random number generator
+  PseRan->SetSeed(4357);                // Set seed
+  TFoam   *FoamX    = new TFoam("FoamX");   // Create Simulator
+  FoamX->SetkDim(2);          // No. of dimensions, obligatory!
+  FoamX->SetnCells(500);      // No. of cells, can be omitted, default=2000
+  FoamX->SetRhoInt(Camel2);   // Set 2-dim distribution, included below
+  FoamX->SetPseRan(PseRan);   // Set random number generator
+  FoamX->Initialize();        // Initialize simulator, takes a few seconds...
+  // From now on FoamX is ready to generate events according to Camel2(x,y)
+  for(Long_t loop=0; loop<100000; loop++){
+    FoamX->MakeEvent();          // generate MC event
+    FoamX->GetMCvect( MCvect);   // get generated vector (x,y)
+    Double_t x=MCvect[0];
+    Double_t y=MCvect[1];
+    if(loop<10) cout<<"(x,y) =  ( "<< x <<", "<< y <<" )"<<endl;
+    //hst_xy->Fill(x,y);           // fill scattergram
+  }// loop
+  Double_t mcResult, mcError;
+  FoamX->GetIntegMC( mcResult, mcError);  // get MC integral, should be one
+  cout << " mcResult= " << mcResult << " +- " << mcError <<endl;
+  // now hst_xy will be plotted visualizing generated distribution
+  // TCanvas *cKanwa = new TCanvas("cKanwa","Canvas for plotting",600,600);
+  //cKanwa->cd();
+  //   hst_xy->Draw("lego2");
+  return 1;
+}//kanwa
+double phi(int n, double *Xarg)
+{
+  double x=Xarg[0];
+  return exp(-x*x/2.)/sqrt(TMath::Pi()*2);
+}
+void PhimFoam()
+{
+  TRandom3  *PseRan   = new TRandom3();  // Create random number generator
+  TFoam   *FoamX    = new TFoam("FoamX");   // Create Simulator
+  FoamX->SetkDim(1);          // No. of dimensions, obligatory!
+  FoamX->SetRhoInt(phi);   // Set 1-dim distribution, included below
+  FoamX->SetPseRan(PseRan);   // Set random number generator
+  FoamX->Initialize();        // Initialize simulator, takes a few seconds...
+  for(Long_t loop=0; loop<100000; loop++){
+    FoamX->MakeEvent();          // generate MC event
+    //FoamX->GetMCvect( MCvect);   // get generated vector (x,y)
+    //Double_t x=MCvect[0];
+    //Double_t y=MCvect[1];
+    //if(loop<10) cout<<"(x,y) =  ( "<< x <<", "<< y <<" )"<<endl;
+    // hst_xy->Fill(x,y);           // fill scattergram
+  }
+  Double_t mcResult, mcError;
+  FoamX->GetIntegMC(mcResult,mcError);  // get MC integral, should be one
+  cout<<"wartosc= "<<mcResult<<" ("<<mcError<<")"<<endl;
+}
+double fsubg(double x)
+{
+  double a=2.*(exp(-0.5)/sqrt(2*TMath::Pi())-1.);
+  double b=2.-exp(-0.5)/sqrt(2*TMath::Pi());
+  // double a=-1.516;
+  //double b=1.758;
+  double f=(exp(-x*x/2.))/sqrt(2.*TMath::Pi());
+  double g=a*x+b;
+  return (f-g);
+}
+double f(double x)
+{
+  return (exp(-x*x/2.))/sqrt(2.*TMath::Pi());
+}
+void ZA(int n, double * wyniki)
+{
+  if(wyniki==NULL){return;}
+  // double a=2.*(exp(-0.5)/sqrt(2*TMath::Pi())-1.);
+  //double b=2.-exp(-0.5)/sqrt(2*TMath::Pi());
+  double wartosc=0.;
+  double sigma=0.;
+  double *wartf=new double[n];
+  double x;
+  for(int i=0;i<n;++i)
+    {
+      x=gRandom->Rndm();
+      wartf[i]=(f(x)+f(1.-x))/2.;
+    }
+  for(int i=0;i<n;++i)
+    {
+      wartosc+=wartf[i];
+    }
+  wartosc/=(double)n;
+  for(int i=0;i<n;++i)
+    {
+      sigma+=pow((wartf[i]-wartosc),2.);
+    }
+  sigma/=(double)(n-1.);
+  sigma=sqrt(sigma/(double)n);
+  wyniki[0]=wartosc;
+  //  cout<<2*(double)n/trafienia<<endl;
+  wyniki[2]=sigma;
+  wyniki[1]=wartosc-0.34134474606854294859;
+}
+void MZK(int n, double * wyniki)
+{
+  if(wyniki==NULL){return;}
+  double a=2.*(exp(-0.5)/sqrt(2*TMath::Pi())-1.);
+  double b=2.-exp(-0.5)/sqrt(2*TMath::Pi());
+  double wartosc=0.;
+  double sigma=0.;
+  double *wartf=new double[n];
+  double x;
+  for(int i=0;i<n;++i)
+    {
+      x=gRandom->Rndm();
+      wartf[i]=fsubg(x);
+    }
+  for(int i=0;i<n;++i)
+    {
+      wartosc+=wartf[i];
+    }
+  wartosc/=(double)n;
+  for(int i=0;i<n;++i)
+    {
+      sigma+=pow((wartf[i]-wartosc),2.);
+    }
+  sigma/=(double)(n-1.);
+  sigma=sqrt(sigma/(double)n);
+  wyniki[0]=wartosc+0.5*a+b;
+  //  cout<<2*(double)n/trafienia<<endl;
+  wyniki[2]=sigma;
+  wyniki[1]=wartosc-0.34134474606854294859;
+}
+double generatorG()
+{
+  double q;//q=(gRandom->Rndm());
+  double P;//P=(exp(-0.5)-1)*q*q/(2*sqrt(2*TMath::Pi()))+q/(2*sqrt(2*TMath::Pi()));
+  double a=2.*(exp(-0.5)/sqrt(2*TMath::Pi())-1.);
+  double b=2.-exp(-0.5)/sqrt(2*TMath::Pi());
+  while(1)
+    {
+      q=gRandom->Rndm();
+      //    P=(exp(-0.5)-1)*q*q/(2*sqrt(2*TMath::Pi()))+q/(2*sqrt(2*TMath::Pi()));
+      P=a*q*q/2.+b*q;
+      //   if((gRandom->Rndm()*(1+exp(-0.5))/(2*sqrt(2*TMath::Pi())))<P)
+      if(gRandom->Rndm()<P)
+	{
+	  return q;
+	}
+    }
+}