RooDataHist * genHistFromModelPdf(const char * name, RooAbsPdf *model,  RooRealVar *var,    double ScaleLumi,  int range, int rebin, int seed ) {
  double genEvents =  model->expectedEvents(*var);
  TRandom3 *rndm = new TRandom3();
  rndm->SetSeed(seed);
  double nEvt = rndm->PoissonD( genEvents) ;
  int intEvt = ( (nEvt- (int)nEvt) >= 0.5) ? (int)nEvt +1 : int(nEvt);
  RooDataSet * data = model->generate(*var ,   intEvt   );
  cout<< " expected events for " << name << " = "<< genEvents << endl; 
  cout<< " data->numEntries() for name " << name << " == " << data->numEntries()<< endl;
  // cout<< " nEvt from PoissonD for" << name << " == " << nEvt<< endl;
  //cout<< " cast of nEvt  for" << name << " == " << intEvt<< endl; 
  RooAbsData *binned_data = data->binnedClone();
  TH1 * toy_hist = binned_data->createHistogram( name, *var, Binning(range/rebin )  );
  for(int i = 1; i <= toy_hist->GetNbinsX(); ++i) {
    toy_hist->SetBinError( i,  sqrt( toy_hist->GetBinContent(i)) );
    if(toy_hist->GetBinContent(i) == 0.00) {
      cout<< " WARNING: histo " << name << " has 0 enter in bin number " << i << endl;   
    }
    if(toy_hist->GetBinContent(i) < 0.1) {
      toy_hist->SetBinContent(i, 0.0);
      toy_hist->SetBinError(i, 0.0);
      cout<< " WARNING: setting value 0.0 to histo " << name << " for bin number " << i << endl;   
    }  
  }
  RooDataHist * toy_rooHist = new RooDataHist(name, name , RooArgList(*var), toy_hist );
  return toy_rooHist; 
}

void gtime2(Int_t nsteps = 200, Int_t np=5000) {
   if (np > 5000) np = 5000;
   Int_t color[5000];
   Double_t cosphi[5000], sinphi[5000], speed[5000];
   TRandom3 r;
   Double_t xmin = 0, xmax = 10, ymin = -10, ymax = 10;
   TGraphTime *g = new TGraphTime(nsteps,xmin,ymin,xmax,ymax);
   g->SetTitle("TGraphTime demo 2;X;Y");
   Int_t i,s;
   Double_t phi,fact = xmax/Double_t(nsteps);
   for (i=0;i<np;i++) { //calculate some object parameters
      speed[i]  = r.Uniform(0.5,1);
      phi       = r.Gaus(0,TMath::Pi()/6.);
      cosphi[i] = fact*speed[i]*TMath::Cos(phi);
      sinphi[i] = fact*speed[i]*TMath::Sin(phi);
      Double_t rc = r.Rndm();
      color[i] = kRed;
      if (rc > 0.3) color[i] = kBlue;
      if (rc > 0.7) color[i] = kYellow;
   }
   for (s=0;s<nsteps;s++) { //fill the TGraphTime step by step
      for (i=0;i<np;i++) {
         Double_t xx = s*cosphi[i];
         if (xx < xmin) continue;
         Double_t yy = s*sinphi[i];
         TMarker *m = new TMarker(xx,yy,25);
         m->SetMarkerColor(color[i]);
         m->SetMarkerSize(1.5 -s/(speed[i]*nsteps));
         g->Add(m,s);
      }
      g->Add(new TPaveLabel(.70,.92,.98,.99,Form("shower at %5.3f nsec",3.*s/nsteps),"brNDC"),s);
   }
   g->Draw();
}

int decroissance_pi(float _p_pi = 1.0){

  //	random.seed(60934386)
  TRandom3 rand;
  rand.SetSeed();
  //float _p_pi;
  //cout << "Entrez l'impulsion des pions (en GeV) : p = ";
  //cin >> _p_pi;

  double x[1000] = {0.0};
  double y[1000] = {0.0};

  for (int i=0;i<1000;i++){
    //double _theta_cm_mu = rand.Uniform(TMath::Pi());
    double _theta_cm_mu = TMath::ACos(gRandom->Uniform(2.0)-1.0);
    double _phi_cm_mu = rand.Uniform(2*TMath::Pi());
    if (i < 10){
      cout << "i: " << i <<  " th_cm_mu = " << _theta_cm_mu << endl;
    }
    double _theta_lab_mu = 0.0, _p_lab_mu = 0.0;
    ThetaLab_mu(_p_pi, _theta_cm_mu, _theta_lab_mu, _p_lab_mu);
    cout << "th_lab_mu= " << _theta_lab_mu << " _p_lab_mu= " << _p_lab_mu << endl;
    x[i] = 1000*_theta_lab_mu;
    y[i] = _p_lab_mu;
  }

  TGraph *graph = new TGraph(1000, x, y);
  graph->SetTitle("Theta vs p_{LAB}");

  graph->Draw("A*");

  return 0;

}

int main(int argc, char *argv[]) {

    TRandom3 *ran = new TRandom3(0);

    vektor b = new double[N];
    vektor x = new double[N];
//     vektor b = malloc_vektor(N);
//     vektor x = malloc_vektor(N);

//     matrix A = new double[N][N];
    matrix A = malloc_matrix(N);
    double dummy = 0.;

    for (int i=0; i<N; ++i) {
        for (int j=0; j<N; ++j) {
            dummy = ran->Uniform();
            A[i][j] = dummy;
            A[j][i] = dummy;
        }
        b[i] = ran->Uniform();
        x[i] = ran->Uniform();
    }

//     print_matrix(A);
    vektor x2 = malloc_vektor(N);
    x2 = vec_copy(cg(N, A, x, b, 1000, 1e-10, 0));
//     print_vektor(x2);

    delete[] b;
    delete[] x;
//     delete[] A;
    return 0;
}

void plot_distortion(){
    
    TRandom3* rand = new TRandom3();

//    TF1 *func = new TF1("func","TMath::Abs([0]+[1]*x)",0,12);
//
//    rand->SetSeed(0);
//    Double_t slope=0.2*rand->Gaus(0,1);
//    Double_t anchor_point=3.5;
//    //want offset to be set by requiring func at anchor point to be 1
//    Double_t offset=(1-slope*anchor_point);
//    func->SetParameter(0,offset);
//    func->SetParameter(1,slope);
    
    //FN distortion
    TF1 *func = new TF1("func","[0]/(0.2*pow(x,0.1))+[1]",0,12);
    rand->SetSeed(0);
    Double_t scaling =0.2*rand->Gaus(0,1);
    func->SetParameter(0,scaling);
    func->SetParameter(1,0);
    //set offset so that func(FinalEnergy)=1;
    func->SetParameter(1,1-1*func->Eval(12));
    
    func->Draw();
    
    
   }

//____________________________________________________________________
void fitCircle(Int_t n=10000) {
   //generates n points around a circle and fit them
   TCanvas *c1 = new TCanvas("c1","c1",600,600);
   c1->SetGrid();
   gr = new TGraph(n);
   if (n> 999) gr->SetMarkerStyle(1);
   else        gr->SetMarkerStyle(3);
   TRandom3 r;
   Double_t x,y;
   for (Int_t i=0;i<n;i++) {
      r.Circle(x,y,r.Gaus(4,0.3));
      gr->SetPoint(i,x,y);
   }
   c1->DrawFrame(-5,-5,5,5);
   gr->Draw("p");

   //Fit a circle to the graph points
   TVirtualFitter::SetDefaultFitter("Minuit");  //default is Minuit
   TVirtualFitter *fitter = TVirtualFitter::Fitter(0, 3);
   fitter->SetFCN(myfcn);

   fitter->SetParameter(0, "x0",   0, 0.1, 0,0);
   fitter->SetParameter(1, "y0",   0, 0.1, 0,0);
   fitter->SetParameter(2, "R",    1, 0.1, 0,0);

   Double_t arglist[1] = {0};
   fitter->ExecuteCommand("MIGRAD", arglist, 0);

   //Draw the circle on top of the points
   TArc *arc = new TArc(fitter->GetParameter(0),
      fitter->GetParameter(1),fitter->GetParameter(2));
   arc->SetLineColor(kRed);
   arc->SetLineWidth(4);
   arc->Draw();
}

//____________________________________________________________________
void fitCircle(Int_t n=10000) {
   //generates n points around a circle and fit them
   TCanvas *c1 = new TCanvas("c1","c1",600,600);
   c1->SetGrid();
   gr = new TGraph(n);
   if (n> 999) gr->SetMarkerStyle(1);
   else        gr->SetMarkerStyle(3);
   TRandom3 r;
   Double_t x,y;
   for (Int_t i=0;i<n;i++) {
      r.Circle(x,y,r.Gaus(4,0.3));
      gr->SetPoint(i,x,y);
   }
   c1->DrawFrame(-5,-5,5,5);
   gr->Draw("p");


   auto chi2Function = [&](const Double_t *par) {
      //minimisation function computing the sum of squares of residuals
      // looping at the graph points
      Int_t np = gr->GetN();
      Double_t f = 0;
      Double_t *x = gr->GetX();
      Double_t *y = gr->GetY();
      for (Int_t i=0;i<np;i++) {
         Double_t u = x[i] - par[0];
         Double_t v = y[i] - par[1];
         Double_t dr = par[2] - std::sqrt(u*u+v*v);
         f += dr*dr;
      }
      return f;
   };

   // wrap chi2 funciton in a function object for the fit
   // 3 is the number of fit parameters (size of array par)
   ROOT::Math::Functor fcn(chi2Function,3);
   ROOT::Fit::Fitter  fitter;


   double pStart[3] = {0,0,1};
   fitter.SetFCN(fcn, pStart);
   fitter.Config().ParSettings(0).SetName("x0");
   fitter.Config().ParSettings(1).SetName("y0");
   fitter.Config().ParSettings(2).SetName("R");

   // do the fit 
   bool ok = fitter.FitFCN();
   if (!ok) {
      Error("line3Dfit","Line3D Fit failed");
   }   

   const ROOT::Fit::FitResult & result = fitter.Result();
   result.Print(std::cout);

   //Draw the circle on top of the points
   TArc *arc = new TArc(result.Parameter(0),result.Parameter(1),result.Parameter(2));
   arc->SetLineColor(kRed);
   arc->SetLineWidth(4);
   arc->Draw();
}

// Code from http://www.hongliangjie.com/2012/12/19/how-to-generate-gamma-random-variables/
// Parameter b could be theta...
double gsl_ran_gamma(const double a, const double b, TRandom3 &rand)
{
  
  if (a < 1)
  {
    double u = rand.Uniform(1);
    return gsl_ran_gamma(1.0 + a, b, rand) * pow (u, 1.0 / a);
  }

  double x, v, u;
  double d = a - 1.0 / 3.0;
  double c = (1.0 / 3.0) / sqrt (d);
  
  while (1) 
  {
    do 
    {
      x = rand.Gaus(0, 1.0);
      v = 1.0 + c * x;
    }
    while (v <= 0);
      
    v = v * v * v;
    u = rand.Uniform(1);

    if (u < 1 - 0.0331 * x * x * x * x) 
      break;

    if (log (u) < 0.5 * x * x + d * (1 - v + log (v)))
      break;
  }
    
  return b * d * v;
}

// ------------------------------------------------------------------------------------------------------//
// ------------------------------------------------------------------------------------------------------//
void frac_error(double a, double b, double *ret){

	double c = (b-a);
	gStyle->SetOptFit(0011);
	gStyle->SetOptStat(0); 
	TCanvas *can = new TCanvas("can","can",900,600);

	double frac = a/b;
	TH1F *h_    = new TH1F("h_","h_",100,frac*0.98,frac*1.02);
	TRandom3 *r  = new TRandom3();
	
	for(int i=0;i<50000;i++){
	 
	  double ap = r->Poisson(a);
	  double cp = r->Poisson(c);
	  h_->Fill((double)ap/(ap+cp));
	  
	}
	
	h_->Fit("gaus","Q","Q");
		
	h_->Draw();
	can->SaveAs("statEP.png");
	double val = h_->GetFunction("gaus")->GetParameter(2);
	delete h_;
	
	ret[0] = val;
}

void eta()
{
  // declare and define two histograms (their number of bins and boundaries)
  TH1F* h_eta = new TH1F("h_eta","",100,-10,10);
  TH1F* h_theta = new TH1F("h_theta","",100,0,TMath::Pi());
  
  // start producing distributions
  // set the number of events to be 10000
  const int neve=10000;
  TRandom3* myrndm = new TRandom3();

  for(unsigned int i=0; i<neve; i++){
    
    // produce a flat theta distributions
    float theta = myrndm->Rndm()*TMath::Pi();
    h_theta->Fill(theta);

    float eta = -TMath::Log(TMath::Tan(0.5*theta));
    h_eta->Fill(eta);
  }
  
  h_theta->Draw();
  h_eta->Draw();

  // save the histograms in a root file
  TFile* outFile = new TFile("eta_theta.root","recreate");
  
  h_theta->Write();
  h_eta->Write();

  outFile->Close();


}

// ------------------------------------------------------------------------------------------------------//
// ------------------------------------------------------------------------------------------------------//
void bias_corr_stat_error(double Sa, double Sb, double Sc, double Sd,
		          double Qa, double Qb, double Qc, double Qd, 
		          double Scale, double *ret){
 
      	double ea = (Sa)/(Sa+Sb);
      	double ed = (Sd)/(Sd+Sc);
      	double ew = (Sa+Sb)/(Sa+Sb+Sd+Sc);
 	double st = Sa+Sb+Sc+Sd;
 
 	double na = (Qa + Sa);
      	double nb = (Qb + Sb);
      	double nc = (Qc + Sc);     		
      	double nd = (Qd + Sd);
 	
 	double S[4] = {0,0,0,0};
 	signal(na,nb,nc,nd,ea,ed,ew,S);
 	
 	double s1 = S[0];
 	double s2 = S[2];
 	
 	double cen1 = st/s1;
 	double cen2 = st/s2;
 	
      	TH1F *tmp_u = new TH1F("tmp_u","tmp_u",100,cen1*(1-0.01),cen1*(1+0.01)); 
      	TH1F *tmp_d = new TH1F("tmp_d","tmp_d",100,cen2*(1-0.01),cen2*(1+0.01));
      	
      	TRandom3 *rnd = new TRandom3();
      	      	
      	for(int i=0;i<50000; i++){
      	
      		double qa = Scale*rnd->Poisson((1./Scale)*Qa);
      		double qb = Scale*rnd->Poisson((1./Scale)*Qb);
      		double qc = Scale*rnd->Poisson((1./Scale)*Qc);
      		double qd = Scale*rnd->Poisson((1./Scale)*Qd);
      		
      		na = (qa + Sa);
      		nb = (qb + Sb);
      		nc = (qc + Sc);     		
      		nd = (qd + Sd);
      		
      		signal(na,nb,nc,nd,ea,ed,ew,S);

      		tmp_u->Fill(st/S[0]);      		     
      		tmp_d->Fill(st/S[2]);      		               	
          	
        }
        
	tmp_u->Fit("gaus","Q","Q");
	tmp_d->Fit("gaus","Q","Q");

	double val1 = tmp_u->GetFunction("gaus")->GetParameter(2);
	double val2 = tmp_d->GetFunction("gaus")->GetParameter(2);
	delete tmp_u;
	delete tmp_d;
	
	ret[0] = val1;
	ret[1] = val2;
}

void macro6(){

    TH1F* sig_h=new TH1F("sig_h","Signal Histo",50,0,10);
    TH1F* gaus_h1=new TH1F("gaus_h1","Gauss Histo 1",30,0,10);
    TH1F* gaus_h2=new TH1F("gaus_h2","Gauss Histo 2",30,0,10);
    TH1F* bkg_h=new TH1F("exp_h","Exponential Histo",50,0,10);

    // simulate the measurements
    TRandom3 rndgen;
    for (int imeas=0;imeas<4000;imeas++){
        bkg_h->Fill(rndgen.Exp(4));
        if (imeas%4==0) gaus_h1->Fill(rndgen.Gaus(5,2));
        if (imeas%4==0) gaus_h2->Fill(rndgen.Gaus(5,2));
        if (imeas%10==0)sig_h->Fill(rndgen.Gaus(5,.5));}

    // Format Histograms
    TH1F* histos[4]={sig_h,bkg_h,gaus_h1,gaus_h2};
    for (int i=0;i<4;++i){
        histos[i]->Sumw2(); // *Very* Important
        format_h(histos[i],i+1);
        }

    // Sum
    TH1F* sum_h= new TH1F(*bkg_h);
    sum_h->Add(sig_h,1.);
    sum_h->SetTitle("Exponential + Gaussian");
    format_h(sum_h,kBlue);

    TCanvas* c_sum= new TCanvas();
    sum_h->Draw("hist");
    bkg_h->Draw("SameHist");
    sig_h->Draw("SameHist");

    // Divide
    TH1F* dividend=new TH1F(*gaus_h1);
    dividend->Divide(gaus_h2);

    // Graphical Maquillage
    dividend->SetTitle(";X axis;Gaus Histo 1 / Gaus Histo 2");
    format_h(dividend,kOrange);
    gaus_h1->SetTitle(";;Gaus Histo 1 and Gaus Histo 2");
    gStyle->SetOptStat(0);

    TCanvas* c_divide= new TCanvas();
    c_divide->Divide(1,2,0,0);
    c_divide->cd(1);
    c_divide->GetPad(1)->SetRightMargin(.01);
    gaus_h1->DrawNormalized("Hist");
    gaus_h2->DrawNormalized("HistSame");

    c_divide->cd(2);
    dividend->GetYaxis()->SetRangeUser(0,2.49);
    c_divide->GetPad(2)->SetGridy();
    c_divide->GetPad(2)->SetRightMargin(.01);
    dividend->Draw();
}

void histo_graph(TH1D *histo,TGraph *graph,Int_t N_to_fill) {
  Int_t N_Points=graph->GetN();
  Double_t x,y;

  Double_t y_min=0.0;
  Double_t y_max=0.0;
  Double_t x_min=0.0;
  Double_t x_max=0.0;

  for (Int_t i=0;i<N_Points;i++) {
    graph->GetPoint(i,x,y);
    if (y>y_max) y_max=y;
    if (x>x_max) x_max=x;
    if (y<y_min) y_min=y;
    if (x<x_min) x_min=x;
  }
  Double_t histo_x_min=histo->GetXaxis()->GetXmin();
  Double_t histo_x_max=histo->GetXaxis()->GetXmax();
  if (histo_x_min>x_min) x_min=histo_x_min;
  if (histo_x_max<x_max) x_max=histo_x_max;

  //reverse graph order if necessary
  TGraph graph2(N_Points);
  Double_t x0,xNm1;
  graph->GetPoint(0,x0,y);
  graph->GetPoint(N_Points-1,xNm1,y);
  if (xNm1<x0) {
    for (Int_t iPoint=0;iPoint<N_Points;iPoint++) {
      Double_t x,y;
      graph->GetPoint(N_Points-iPoint,x,y);
      graph2.SetPoint(iPoint,x,y);
    }
    graph=&graph2;
  }

  TRandom3 r;
  r.SetSeed(0);

  Int_t N_tried=0;
  Int_t N_filled=0;
  while (N_filled<N_to_fill) {
    N_tried++;
    //rejection method
    x=r.Uniform(x_min,x_max);
    y=r.Uniform(0.0,y_max);
    
    if (y<graph->Eval(x)) {
      histo->Fill(x);
      N_filled++;
    }
  }
  printf("x_max=%8.6g; x_min=%8.6g; y_max=%8.6g; N_filled=%d; N_tried=%d\n",x_max,x_min,y_max,N_filled,N_tried);
  Double_t norm=(x_max-x_min)*y_max*N_filled/N_tried;
  histo->Scale(norm/histo->Integral("width"));

}

//ks-test for Gaussian Random samples
void test(Int_t num)
{
  Double_t *d=new Double_t[num];
  Double_t p,t;
  TRandom3 r;
  for (int i=0;i<num;i++)
    d[i]=r.Gaus();
  ROOT::Math::GoFTest *tt=new ROOT::Math::GoFTest(num,d,ROOT::Math::GoFTest::kGaussian);
  tt->KolmogorovSmirnovTest(p,t);
  cout<<"p:"<<p<<",t:"<<t<<endl;
}

void generateGaussRandom(){

  TRandom3 *gr = new TRandom3(12345);
  cout<<"generating gaus rand for use "<< Nrand<<endl; 
  for(int j=0;j<Nrand;j++){
    if(j%100000000==0) cout<<" j " << j <<endl; 
    float tmp = gr->Gaus();
    vrandgaus.push_back(tmp);
  }
  
}

//______________________________________________________________________________
// Fill y and error vectors with random points from TRandom3
void FillRandVectors(int nPoints, vector<double> &xVector, vector< double > &yVector, vector< double > &yErrorVector, int seed = 250, double lowerBound= 10, double upperBound= 20, double lowerErrorBound = 1, double upperErrorBound= 2)
{
	//Call TRandom3
	TRandom3 *jrand = new TRandom3(seed);
	for (int i = 0; i < nPoints; i++)
	{
		xVector.push_back(i);
		yVector.push_back(jrand->Uniform(lowerBound, upperBound));
		yErrorVector.push_back(jrand->Uniform(lowerErrorBound, upperErrorBound));
	}
	delete jrand;
}

Double_t Reconstruct( Double_t xt, TRandom3& R )
{
   // apply some Gaussian smearing + bias and efficiency corrections to fake reconstruction
   const Double_t cutdummy = -99999.0;
   Double_t xeff = 0.3 + (1.0 - 0.3)/20.0*(xt + 10.0);  // efficiency
   Double_t x    = R.Rndm();
   if (x > xeff) return cutdummy;
   else {
     Double_t xsmear= R.Gaus(-2.5,0.2); // bias and smear
     return xt+xsmear;
   }
}

// Read in an event spectrum from a text file containing
// spectral data, and then generate some random events.
// This assumes that the input data file contains two columns,
// the first with (evenly spaced!) energy bins in units of eV, and the second
// with a number of events per kg per day per keV
void FillRecoilSpectrumFromFile(TH1F* recoilHisto, Double_t time, Double_t mass, const char* filename, Int_t seed)
{
	FILE* dataFile;
	dataFile = fopen(filename,"r");
	if (dataFile == NULL) {
		cout << "The file " << filename << " is not found!!! Quitting now." << endl;
		exit(1);
	}
	
	// Count the number of lines to set the array sizes
	fstream lineCountStream;
	lineCountStream.open(filename, fstream::in);
	Int_t lineCount = 0;
	while (lineCountStream.peek() != EOF) {
		lineCountStream.ignore(128, '\n');
		lineCount++;
	}
	lineCountStream.close();
	cout << "The file has " << lineCount << " lines." << endl;
	
	// Should add a section here to double check the formatting of the file
	
	// Warning message about units
	cout << "Reading data and filling histogram from " << filename << "..." << endl;
	cout << "I am assuming that the first column in the file is in units of MeV "
			 << "and that the second column is in units of 1/(kg day keV)!" << endl << endl;
	
	// Read the data
	const Int_t nLines = lineCount;
	Float_t energies[nLines];
	Float_t rates[nLines];
	for (Int_t thisLine = 0; thisLine < nLines; thisLine++)
		fscanf(dataFile, "%E %E", &energies[thisLine], &rates[thisLine]);
	fclose(dataFile);
	
	// Set up the source histogram
	TH1F SpectralDataHistogram("SpectralDataHistogram", "Spectrum of Events from File", nLines-1, energies);
	for (Int_t thisBin = 0; thisBin < (nLines - 1); thisBin++) {
		SpectralDataHistogram.SetBinContent(thisBin,rates[thisBin]);
	}
	
	// Fill the recoil histogram with the correct number
	// of events, randomly chosen from the source histogram.
	TRandom3* randGen = new TRandom3(seed);
	Float_t meanNEvents = time * mass * SpectralDataHistogram.ComputeIntegral();
	Int_t nEvents = randGen->Poisson(meanNEvents);
	cout << "nEvents = " << nEvents << endl;
	for (Int_t thisEvent; thisEvent<nEvents; thisEvent++)
		recoilHisto->Fill(SpectralDataHistogram.GetRandom());
}

void TestBkg(){

  
  TFile * f = new TFile("./TestBkg.root","recreate");
  TTree * tree = new TTree("flt","data");
			
  LendaEvent *e = new LendaEvent();
  tree->Branch("Event",&e);
			
  TRandom3 * r = new TRandom3();
  for (int i=0;i<100000;i++){

    Double_t v =    r->Uniform()*30.0;

    e->pushTime(30);
    e->pushTime(30);
    
    e->pushTime(v);
    
    e->pushEnergy(1);
    e->pushEnergy(1);
    e->pushEnergy(1);

    e->pushLongGate(2);
    e->pushLongGate(2);
    e->pushLongGate(2);


    e->pushShortGate(2);
    e->pushShortGate(2);
    e->pushShortGate(2);

    e->pushCubicTime(e->times[0]);
    e->pushCubicTime(e->times[1]);
    e->pushCubicTime(e->times[2]);

    
    e->pushSoftwareCFD(1);    e->pushSoftwareCFD(1);    e->pushSoftwareCFD(1);
    
    
    e->Finalize();

    tree->Fill();
    e->Clear();
  }

  tree->Write();

  f->Close();
}

void write_ntuple_to_file(){

    TFile ofile("conductivity_experiment.root","RECREATE");

    // Initialise the TNtuple
    TNtuple cond_data("cond_data",
                      "Example N-Tuple",
                      "Potential:Current:Temperature:Pressure");

    // Fill it randomly to fake the acquired data
    TRandom3 rndm;
    float pot,cur,temp,pres;
    for (int i=0;i<10000;++i){
        pot=rndm.Uniform(0.,10.);      // get voltage
        temp=rndm.Uniform(250.,350.);  // get temperature
        pres=rndm.Uniform(0.5,1.5);    // get pressure
        cur=pot/(10.+0.05*(temp-300.)-0.2*(pres-1.)); // current
// add some random smearing (measurement errors)
        pot*=rndm.Gaus(1.,0.01); // 1% error on voltage
        temp+=rndm.Gaus(0.,0.3); // 0.3 abs. error on temp.
        pres*=rndm.Gaus(1.,0.02);// 1% error on pressure
        cur*=rndm.Gaus(1.,0.01); // 1% error on current
// write to ntuple
        cond_data.Fill(pot,cur,temp,pres);
        }

    // Save the ntuple and close the file
    cond_data.Write();
    ofile.Close();
}