C++ TRandom::Rndm方法代码示例

TCanvas* graph2dfit_test()
{
	gStyle->SetOptStat(0);
	gStyle->SetOptFit();

	TCanvas *c = new TCanvas("c", "Graph2D example", 0, 0, 600, 800);
	
	Double_t rnd, x, y, z;
	Double_t e = 0.3;
	Int_t nd = 400;
	Int_t np = 10000;

	TRandom r;
	Double_t fl = 6;
	TF2  *f2 = new TF2("f2", "1000*(([0]*sin(x)/x)*([1]*sin(y)/y))+200", -fl, fl, -fl, fl);
	f2->SetParameters(1, 1);
	TGraph2D *dt = new TGraph2D();

	// Fill the 2D graph
	Double_t zmax = 0;
	for (Int_t N = 0; N<nd; N++) {
		f2->GetRandom2(x, y);
		// Generate a random number in [-e,e]
		rnd = 2 * r.Rndm()*e - e;
		z = f2->Eval(x, y)*(1 + rnd);
		if (z>zmax) zmax = z;
		dt->SetPoint(N, x, y, z);
	}


	f2->SetParameters(0.5, 1.5);
	dt->Fit(f2);
	TF2 *fit2 = (TF2*)dt->FindObject("f2");

	f2->SetParameters(1, 1);

	for (Int_t N = 0; N<np; N++) {
		f2->GetRandom2(x, y);
		// Generate a random number in [-e,e]
		rnd = 2 * r.Rndm()*e - e;
		z = f2->Eval(x, y)*(1 + rnd);
		h1->Fill(f2->Eval(x, y) - z);
		z = dt->Interpolate(x, y);
		h2->Fill(f2->Eval(x, y) - z);
		z = fit2->Eval(x, y);
		h3->Fill(f2->Eval(x, y) - z);
	}

	gStyle->SetPalette(1);
	
	f2->SetTitle("Original function with Graph2D points on top");
	f2->SetMaximum(zmax);
	gStyle->SetHistTopMargin(0);
	f2->Draw("surf1");
	dt->Draw("same p0");



	return c;
}

void GenerateNonStd(unsigned int Npart)
{
  frodo *fr = frodo::instance();

  static TRandom randy(0);

  fr->HBDParticles.clear();
  
  double s = randy.Rndm();
  //cout<<"s= "<<s<<endl;
  double s1 = 0;
  if(s < 0.5) s1 = 1;
  else s1 = -1;
  
  double a = randy.Rndm();
  //cout<<"a= "<<a <<endl;
  double Cx = 2*fr->Ring[0].C()*a - fr->Ring[0].C();
  double Cy = s1*sqrt(fr->Ring[0].C()*fr->Ring[0].C()-Cx*Cx);
  fr->CX.push_back(Coordinates(Cx));
  fr->CY.push_back(Coordinates(Cy));
  //cout<<"ring center is: ( "<< Cx <<" , " << Cy << " ) " <<endl;

  for (unsigned int i=0; i<Npart; i++)
    {
      //cout<<"Npart= " <<i<<endl;
      double S = randy.Rndm();
      //cout<<"S= "<<S<<endl;
      double s2 = 0;
      if(S < 0.5) s2 = 1;
      else s2 = -1;

      double b = randy.Rndm();
      //cout<<"b= "<<b<<endl;
      
      double x = Cx + 2*fr->Ring[0].R()*b - fr->Ring[0].R();      
      double y = Cy + s2*sqrt(fr->Ring[0].R()*fr->Ring[0].R() - (x-Cx)*(x-Cx));
      
      if(fabs(x) <= 1000 || fabs(y) <= 1125)
	{
	  i--;
	  //cout<<"reject!"<<endl;
	  continue;
	}
      double q = -1;
      while (q<=0) q=randy.Exp(80000);
      //cout<<"q= " <<q<<endl;
      
      fr->HBDParticles.push_back( AParticle(x,y,q) );
    }
  //cout<<"# of particles: " <<fr->HBDParticles.size()<<endl;
  
  
}

//______________________________________________________________________
void generate_random(Int_t i)
{
   const Double_t dr = 3.5;

   r.Rannor(r1, r4);
   r.Rannor(r7, r9);

   r2 = (2 * dr * r.Rndm(i)) - dr;
   r3 = (2 * dr * r.Rndm(i)) - dr;
   r5 = (2 * dr * r.Rndm(i)) - dr;
   r6 = (2 * dr * r.Rndm(i)) - dr;
   r8 = (2 * dr * r.Rndm(i)) - dr;
}

void binomialFancy() {
  Double_t x;
  Double_t y;
  Double_t res1;
  Double_t res2;
  Double_t err;
  Double_t serr=0;
  const Int_t nmax=10000;
  printf("\nTMath::Binomial fancy test\n");
  printf("Verify Newton formula for (x+y)^n\n");
  printf("x,y in [-2,2] and n from 0 to 9  \n");
  printf("=================================\n");
  TRandom r;
  for(Int_t i=0; i<nmax; i++) {
    do {
        x=2*(1-2*r.Rndm());
        y=2*(1-2*r.Rndm());
    } while (TMath::Abs(x+y)<0.75); //Avoid large cancellations
    for(Int_t j=0; j<10; j++) {
       res1=TMath::Power(x+y,j);
       res2=0;
       for(Int_t k=0; k<=j; k++)
          res2+=TMath::Power(x,k)*TMath::Power(y,j-k)*TMath::Binomial(j,k);
       if((err=TMath::Abs(res1-res2)/TMath::Abs(res1))>1e-10)
 	      printf("res1=%e res2=%e x=%e y=%e err=%e j=%d\n",res1,res2,x,y,err,j);
       serr +=err;
     }
  }
  printf("Average Error = %e\n",serr/nmax);
}

void multicolor(Int_t isStack=0) {
   TCanvas *c1 = new TCanvas;
   Int_t nbins = 20;
   TH2F *h1 = new TH2F("h1","h1",nbins,-4,4,nbins,-4,4);
   h1->SetFillColor(kBlue);
   TH2F *h2 = new TH2F("h2","h2",nbins,-4,4,nbins,-4,4);
   h2->SetFillColor(kRed);
   TH2F *h3 = new TH2F("h3","h3",nbins,-4,4,nbins,-4,4);
   h3->SetFillColor(kYellow);
   THStack *hs = new THStack("hs","three plots");
   hs->Add(h1);
   hs->Add(h2);
   hs->Add(h3);
   TRandom r;
   Int_t i;
   for (i=0;i<20000;i++) h1->Fill(r.Gaus(),r.Gaus());
   for (i=0;i<200;i++) {
      Int_t ix = (Int_t)r.Uniform(0,nbins);
      Int_t iy = (Int_t)r.Uniform(0,nbins);
      Int_t bin = h1->GetBin(ix,iy);
      Double_t val = h1->GetBinContent(bin);
      if (val <= 0) continue;
      if (!isStack) h1->SetBinContent(bin,0);
      if (r.Rndm() > 0.5) {
         if (!isStack) h2->SetBinContent(bin,0);
         h3->SetBinContent(bin,val);
      } else {
         if (!isStack) h3->SetBinContent(bin,0);
         h2->SetBinContent(bin,val);
      }
   }
   hs->Draw("lego1");
}

void graph2derrorsfit()
{
   TCanvas *c1 = new TCanvas("c1");

   Double_t rnd, x, y, z, ex, ey, ez;
   Double_t e = 0.3;
   Int_t nd = 500;

   TRandom r;
   TF2  *f2 = new TF2("f2","1000*(([0]*sin(x)/x)*([1]*sin(y)/y))+200",-6,6,-6,6);
   f2->SetParameters(1,1);
   TGraph2DErrors *dte = new TGraph2DErrors(nd);

   // Fill the 2D graph
   for (Int_t i=0; i<nd; i++) {
      f2->GetRandom2(x,y);      
      rnd = r.Uniform(-e,e); // Generate a random number in [-e,e]
      z = f2->Eval(x,y)*(1+rnd);
      dte->SetPoint(i,x,y,z);
      ex = 0.05*r.Rndm();
      ey = 0.05*r.Rndm();
      ez = TMath::Abs(z*rnd);
      dte->SetPointError(i,ex,ey,ez);
   }

   f2->SetParameters(0.5,1.5);
   dte->Fit(f2);
   TF2 *fit2 = (TF2*)dte->FindObject("f2");
   fit2->SetTitle("Minuit fit result on the Graph2DErrors points");
   fit2->Draw("surf1");
   dte->Draw("axis p0");
}

void triangles(Int_t ntriangles=50) {
  TCanvas *c1 = new TCanvas("c1","triangles",10,10,700,700);
  TRandom r;
  Double_t dx = 0.2; Double_t dy = 0.2;
  Int_t ncolors = gStyle->GetNumberOfColors();
  Double_t x[4],y[4];
  TColor *c;
  Int_t ci;
  for (Int_t i=0;i<ntriangles;i++) {
     x[0] = r.Uniform(.05,.95); y[0] = r.Uniform(.05,.95);
     x[1] = x[0] + dx*r.Rndm(); y[1] = y[0] + dy*r.Rndm();
     x[2] = x[1] - dx*r.Rndm(); y[2] = y[1] - dy*r.Rndm();
     x[3] = x[0];               y[3] = y[0];
     TPolyLine *pl = new TPolyLine(4,x,y);
     pl->SetUniqueID(i);
     ci = ncolors*r.Rndm();
     c  = gROOT->GetColor(TColor::GetColorPalette(ci));
     c->SetAlpha(r.Rndm());
     pl->SetFillColor(ci);
     pl->Draw("f");
  }
  c1->AddExec("ex","TriangleClicked()");
}

void simpleTree()
{
    using std::cout;
    using std::cerr;
    using std::endl;

    int event;
    double px;
    double py;
    double pz;
    double random;

    TFile output("simpleTree.root", "RECREATE");
    if (!output.IsOpen())
    {
        cerr << "failed to open output file." << endl;

        return;
    }
    
    // Create Tree
    TTree tree("T", "simple tree");
    tree.Branch("event", &event, "event/I");
    tree.Branch("px", &px, "px/F");
    tree.Branch("py", &py, "py/F");
    tree.Branch("pz", &pz, "pz/F");
    tree.Branch("random", &random, "random/F");

    // Initialize Random generator
    TRandom rnd;

    // Fill Tree
    for(event = 0; 1000 > event; ++event)
    {
        rnd.Rannor(px, py);
        pz = px * px + py * py;
        random = rnd.Rndm();

        tree.Fill();
    }

    tree.Write();
 }

void circular() {
   auto T = new TTree("T","test circular buffers");
   TRandom r;
   Float_t px,py,pz;
   Double_t randomNum;
   UShort_t i;
   T->Branch("px",&px,"px/F");
   T->Branch("py",&py,"px/F");
   T->Branch("pz",&pz,"px/F");
   T->Branch("random",&randomNum,"random/D");
   T->Branch("i",&i,"i/s");
   T->SetCircular(20000); //keep a maximum of 20000 entries in memory
   for (i = 0; i < 65000; i++) {
      r.Rannor(px,py);
      pz = px*px + py*py;
      randomNum = r.Rndm();
      T->Fill();
   }
   T->Print();
}

double sampleXgenFromFunction1d(TRandom& rnd, TF1* function, double xMin, double xMax, double pMax)
{
  //std::cout << "<sampleXgenFromFunction1d>:" << std::endl;
  //std::cout << " function = " << function->GetTitle() << std::endl;
  //int numParameter = function->GetNpar();
  //for ( int iParameter = 0; iParameter < numParameter; ++iParameter ) {
  //  std::cout << " parameter #" << iParameter << " = " << function->GetParameter(iParameter) << std::endl;
  //}
  //assert(0);
  double x;
  bool isDone = false;
  while ( !isDone ) {
    x = rnd.Uniform(xMin, xMax);
    double u = rnd.Rndm();
    double f_x = function->Eval(x);
    //std::cout << "x = " << x << ": f(x) = " << f_x << " (pMax = " << pMax << ")" << std::endl;
    isDone = ((u*pMax) < f_x);
  }
  return x;
}

TCanvas* graph2dfit()
{
   gStyle->SetOptStat(0);
   gStyle->SetOptFit();

   TCanvas *c = new TCanvas("c","Graph2D example",0,0,600,800);
   c->Divide(2,3);

   Double_t rnd, x, y, z;
   Double_t e = 0.3;
   Int_t nd = 400;
   Int_t np = 10000;

   TRandom r;
   Double_t fl = 6;
   TF2  *f2 = new TF2("f2","1000*(([0]*sin(x)/x)*([1]*sin(y)/y))+200",
      -fl,fl,-fl,fl);
   f2->SetParameters(1,1);
   TGraph2D *dt = new TGraph2D();

   // Fill the 2D graph
   Double_t zmax = 0;
   for (Int_t N=0; N<nd; N++) {
      f2->GetRandom2(x,y);
      // Generate a random number in [-e,e]
      rnd = 2*r.Rndm()*e-e;
      z = f2->Eval(x,y)*(1+rnd);
      if (z>zmax) zmax = z;
      dt->SetPoint(N,x,y,z);
   }

   Double_t hr = 350;
   TH1D *h1 = new TH1D("h1",
   "#splitline{Difference between Original}{#splitline{function and Function}{with noise}}",
   100, -hr, hr);
   TH1D *h2 = new TH1D("h2",
   "#splitline{Difference between Original}{#splitline{function and Delaunay triangles}{interpolation}}",
   100, -hr, hr);
   TH1D *h3 = new TH1D("h3",
   "#splitline{Difference between Original}{function and Minuit fit}",
   500, -hr, hr);

   f2->SetParameters(0.5,1.5);
   dt->Fit(f2);
   TF2 *fit2 = (TF2*)dt->FindObject("f2");

   f2->SetParameters(1,1);

   for (Int_t N=0; N<np; N++) {
      f2->GetRandom2(x,y);
      // Generate a random number in [-e,e]
      rnd = 2*r.Rndm()*e-e;
      z = f2->Eval(x,y)*(1+rnd);
      h1->Fill(f2->Eval(x,y)-z);
      z = dt->Interpolate(x,y);
      h2->Fill(f2->Eval(x,y)-z);
      z = fit2->Eval(x,y);
      h3->Fill(f2->Eval(x,y)-z);
   }

   c->cd(1);
   f2->SetTitle("Original function with Graph2D points on top");
   f2->SetMaximum(zmax);
   gStyle->SetHistTopMargin(0);
   f2->Draw("surf1");
   dt->Draw("same p0");

   c->cd(3);
   dt->SetMargin(0.1);
   dt->SetFillColor(36);
   dt->SetTitle("Histogram produced with Delaunay interpolation");
   dt->Draw("surf4");

   c->cd(5);
   fit2->SetTitle("Minuit fit result on the Graph2D points");
   fit2->Draw("surf1");

   h1->SetFillColor(47);
   h2->SetFillColor(38);
   h3->SetFillColor(29);

   c->cd(2); h1->Fit("gaus","Q") ; h1->Draw();
   c->cd(4); h2->Fit("gaus","Q") ; h2->Draw();
   c->cd(6); h3->Fit("gaus","Q") ; h3->Draw();
   c->cd();
   return c;
}

void Trasporta(Int_t s,Int_t Rhum, TRandom *GeneratoreEsterno=0,Int_t Noise_Medio = 20) {

	TStopwatch tempo;

	tempo.Start(kTRUE);

	cout<<endl<<"Sto trasportando attraverso i rivelatori: attendere... "<<endl;
	TRandom *smear;
	if(GeneratoreEsterno == 0){
		smear = new TRandom3();
		cout<<"Generatore Interno";
	}else{
		smear = GeneratoreEsterno;
		cout<<"Generatore Esterno";
	}
	
	cout<<" FirstRNDM: "<<smear->Rndm()<<endl;
	
	//////////////////////////////////////////////////////
	//Creo un nuovo file e
	//Definisco Struct per salvare i nuovi dati x y z 
	//////////////////////////////////////////////////////

	//Definisco il nuovo albero per salvare i punti di hit	
	TFile sfile("trasporto_tree.root","RECREATE");
  
	TTree *trasporto = new TTree("Ttrasporto","TTree con 3 branches");
	 //Punti sul layer
	TTree *Rel_Lay1 = new TTree("Layer1","TTree con 1 branch");
	TTree *Rel_Lay2 = new TTree("Layer2","TTree con 1 branch");
	 //rumore
	TTree *Noise = new TTree("Rumore","TTree con 1 branch");

	typedef struct {
		Double_t X,Y,Z;
		Int_t Flag;		
	} HIT; 
	static HIT beam;  	
	static HIT lay1;  
	static HIT lay2;

	typedef struct {
		Int_t event;
		Int_t tipo;
		Int_t Noiselay1;
		Int_t Noiselay2;
	} infoRumore;
	static infoRumore InfoR;


	//Dichiaro i rami dei tree
	trasporto->Branch("BeamPipe",&beam.X,"X/D:Y:Z:Flag/I");  
	trasporto->Branch("Layer1",&lay1.X,"X/D:Y:Z:Flag/I"); 
	trasporto->Branch("Layer2",&lay2.X,"X/D:Y:Z:Flag/I");  

	Rel_Lay1->Branch("RealLayer1",&lay1.X,"X/D:Y:Z:Flag/I"); 
	Rel_Lay2->Branch("RealLayer2",&lay2.X,"X/D:Y:Z:Flag/I"); 
	
	Noise->Branch("Rumore",&InfoR,"event/I:tipo:Noiselay1:Noiselay2"); 
	Double_t temp_phi = 0;
	Int_t Nnoise=0;	

  ////////////////////////////////
  //Acquisizione Vertici
  ///////////////////////////////
  TClonesArray *dir = new TClonesArray("Direction",100);	
  typedef struct {
    Double_t X,Y,Z;
    Int_t N;
  }SINGLE_EVENT;
  static SINGLE_EVENT event;    //struct con molteplicita' e vertice di un singolo evento
	
  TFile hfile("event_tree.root");


  TTree *Born = (TTree*)hfile.Get("T");       
  TBranch *b1=Born->GetBranch("Event"); 
  TBranch *b2=Born->GetBranch("Direzioni");  //acquisisco i due branches


  b1->SetAddress(&event.X); //passo l'indirizzo del primo oggetto della struct e assegno tutto a b1
  b2->SetAddress(&dir); // lo stesso per il vettore 




  /////////////////////////
  //Geometria del rivelatore
  /////////////////////////
  Double_t R1=3;	//raggio 3 cm beam pipe
  Double_t R2=4;	//raggio 4 cm primo layer
  Double_t R3=7;	//raggio 7 cm secondo layer

  Double_t limit = 8.; //lunghezza layer su z-> z in [-8,8]

  //Variabili Varie
	Double_t Xo=0.;Double_t Yo=0.;Double_t Zo=0.;
  	Double_t X1=0.;Double_t Y1=0.;Double_t Z1=0.;
	Double_t X2=0.;Double_t Y2=0.;Double_t Z2=0.;

//.........这里部分代码省略.........

Int_t runPrefetchReading(bool caching = false)
{
   //const char *options = 0;
   Int_t freq = 1000; 
   Int_t cachesize = -1;
   Float_t percententries = 1.00;
   Float_t percentbranches = 1.00;
   TStopwatch sw;
   
   //set the reading mode to async prefetching
   gEnv->SetValue("TFile.AsyncPrefetching", 1);

   //enable the local caching of blocks
   TString cachedir="file:/tmp/xcache/";
   // or using xrootd on port 2000
   // TString cachedir="root://localhost:2000//tmp/xrdcache1/";
   if (caching) gEnv->SetValue("Cache.Directory", cachedir.Data());
   
   // open the local if any
   TString filename("atlasFlushed.root");
   if (gSystem->AccessPathName(filename,kReadPermission) && filename.Index(":") == kNPOS) {
      // otherwise open the http file
      filename.Prepend("https://root.cern.ch/files/");
      //filename.Prepend("root://cache01.usatlas.bnl.gov//data/test1/");
      //filename.Prepend( "root://pcitdss1401//tmp/" );
      //filename.Prepend("http://www-root.fnal.gov/files/");
      //filename.Prepend("http://oink.fnal.gov/distro/roottest/");
   }
   
   TString library("atlasFlushed/atlasFlushed");
   fprintf(stderr,"Starting to load the library\n");
   gSystem->Load(library);

   fprintf(stderr,"Starting to open the file\n");
   TFile *file = TFile::Open( filename, "TIMEOUT=30" );
   if (!file || file->IsZombie()) {
      Error("runPrefetchReading","Could not open the file %s within 30s",filename.Data());
      return 1;
   }
   fprintf(stderr,"The file has been opened, setting up the TTree\n");

   // file->MakeProject("atlasFlushed","*","RECREATE+");

   // Try the known names :)
   const char *names [] = { "E","Events","CollectionTree","ntuple","T" };
   TTree *T = NULL;
   for (unsigned int i = 0; i < sizeof(names)/sizeof(names[0]); ++i) {
      file->GetObject(names[i], T);
      if (T) break;
   }
   if (T==0) {
     Error("runPrefetchReading","Could not find a tree which the conventional names in %s.",filename.Data());
     return 2;
   }
   TFile::SetReadaheadSize(0);  // (256*1024);
   Long64_t nentries = T->GetEntries();

   int efirst = 0;
   int elast  = efirst+nentries;
   if (cachesize == -2) {
      gEnv->SetValue("TFile.AsyncReading", 0);
      cachesize = -1;
   }
   T->SetCacheSize(cachesize);

   if (cachesize != 0) {
      T->SetCacheEntryRange(efirst,elast);
      if (percentbranches < 1.00) {
         int nb = T->GetListOfBranches()->GetEntries();
         int incr = nb * percentbranches;
         for(int b=0;b < nb; b += incr) T->AddBranchToCache(((TBranch*)T->GetListOfBranches()->At(b)),kTRUE);
      } else {
         T->AddBranchToCache("*");
      }
      T->StopCacheLearningPhase();
   }

   //...........................................................................
   // First read, with saving the info in cache
   //...........................................................................
   fprintf(stderr,"Setup done. Starting to read the entries\n");
   TRandom r;
   for (Long64_t i = efirst; i < elast; i++) {
     //if (i%100 == 0 || i>2000) fprintf(stderr,"i.debug = %lld\n",i);
     // if (i==2000) gDebug = 7;
     if (i % freq == 0){
       // for (Long64_t i=elast-1;i>=efirst;i--) {
       if (i%freq == 0 || i==(elast-1)) fprintf(stderr,"i = %lld\n",i);
       if (r.Rndm() > percententries) continue; 
       T->LoadTree(i);
       if (percentbranches < 1.00) {
         int nb = T->GetListOfBranches()->GetEntries();
         int incr = nb * percentbranches;
         for(int b=0;b<nb; b += incr) ((TBranch*)T->GetListOfBranches()->At(b))->GetEntry(i);   
         int count = 0;
         int maxcount = 1000 + 100 ;
         for(int x = 0; x < maxcount; ++x ) { /* waste cpu */ count = sin(cos((double)count)); }
       } else {
         T->GetEntry(i);
       }
//.........这里部分代码省略.........

void test()
{
  TRandom *rnd = new TRandom(time(0));
  Double_t x = rnd->Rndm();
  cout << "x = " << x << endl;
}

int main(int argc, char* argv[]) {
  TRandom* r = new TRandom1();

  // Test 1: test multiple files
  // Produces a flat background 0-10 MeV and 0-6000mm radius
  TFile* f1 = new TFile("data1_1.root", "RECREATE");
  TNtuple* n1 = new TNtuple("ntuple", "data from ascii file", "energy:radius");
  for (int i=0; i<10000; i++) {
    n1->Fill(r->Rndm() * 10, pow(r->Rndm(), 0.3333) * 6000);
  }
  f1->Write();
  f1->Close();

  TFile* f2 = new TFile("data1_2.root", "RECREATE");
  TNtuple* n2 = new TNtuple("ntuple", "data from ascii file", "energy:radius");
  for (int i=0; i<10000; i++){
    n2->Fill(r->Rndm() * 10, pow(r->Rndm(), 0.3333) * 6000);
  }
  f2->Write();
  f2->Close();

  // Test 2: test multiple contributions
  // Produces a flat background 0-5 MeV and a spike at 8 MeV near the edge
  TFile* f3 = new TFile("data2_1.root", "RECREATE");
  TNtuple* n3 = new TNtuple("ntuple", "data from ascii file", "energy:radius");
  for (int i=0; i<10000; i++){
    n3->Fill(r->Rndm()*5, pow(r->Rndm(), 0.3333) * 6000);
  }
  f3->Write();
  f3->Close();

  TFile* f4 = new TFile("data2_2.root", "RECREATE");
  TNtuple* n4 = new TNtuple("ntuple", "data from ascii file", "energy:radius");
  for (int i=0; i<10000; i++){
    n4->Fill(8, pow(r->Rndm(), 0.1) * 6000);
  }
  f4->Write();
  f4->Close();

  // Test 3: test signal with zero expected
  // Produces a flat background from 0-1 MeV
  TFile* f5 = new TFile("data3_1.root", "RECREATE");
  TNtuple* n5 = new TNtuple("ntuple", "data from ascii file", "energy:radius");
  for (int i=0; i<10000; i++){
    n5->Fill(r->Rndm(), pow(r->Rndm(), 0.33333) * 6000);
  }
  f5->Write();
  f5->Close();

  // Test 4: test signal with expected but all cut
  // Produces a flat background 0-5 MeV all at the edge
  TFile* f6 = new TFile("data4_1.root", "RECREATE");
  TNtuple* n6 = new TNtuple("ntuple", "data from ascii file", "energy:radius");
  for (int i=0; i<10000; i++){
    n6->Fill(r->Rndm() * 5, 6000);
  }
  f6->Write();
  f6->Close();

  return 0;
}