void candlehisto()
{
   TCanvas *c1 = new TCanvas("c1","Candle Presets",800,600);
   c1->Divide(3,2);

   TRandom *rand = new TRandom();
   TH2I *h1 = new TH2I("h1","Sin",18,0,360,100,-1.5,1.5);
   h1->GetXaxis()->SetTitle("Deg");

   float myRand;
   for (int i = 0; i < 360; i+=10) {
      for (int j = 0; j < 100; j++) {
         myRand = rand->Gaus(sin(i*3.14/180),0.2);
         h1->Fill(i,myRand);
      }
   }

   for (int i = 1; i < 7; i++) {
      c1->cd(i);
      char str[16];
      sprintf(str,"CANDLEX%d",i);
      TH2I * myhist = (TH2I*)h1->DrawCopy(str);
      myhist->SetTitle(str);
   }

   TCanvas *c2 = new TCanvas("c2","Violin Presets",800,300);
   c2->Divide(2,1);

   for (int i = 1; i < 3; i++) {
      c2->cd(i);
      char str[16];
      sprintf(str,"VIOLINX%d",i);
      TH2I * myhist = (TH2I*)h1->DrawCopy(str);
      myhist->SetFillColor(kGray+2);
   }

   TCanvas *c3 = new TCanvas("c3","Playing with candle and violin-options",800,600);
   c3->Divide(3,2);
   char myopt[6][16] = {"1000000","2000000","3000000","1112111","112111","112111"};
   for (int i = 0; i < 6; i++) {
      c3->cd(i+1);
      char str[16];
      sprintf(str, "candlex(%s)",myopt[i]);
      TH2I * myhist = (TH2I*)h1->DrawCopy(str);
      myhist->SetFillColor(kYellow);
      if (i == 4) {
         TH2I * myhist2 = (TH2I*)h1->DrawCopy("candlex(1000000) same");
         myhist2->SetFillColor(kRed);
      }
      if (i == 5) {
         myhist->SetBarWidth(0.2);
         myhist->SetBarOffset(0.25);
         TH2I * myhist2 = (TH2I*)h1->DrawCopy("candlex(2000000) same");
         myhist2->SetFillColor(kRed);
         myhist2->SetBarWidth(0.6);
         myhist2->SetBarOffset(-0.5);
      }
      myhist->SetTitle(str);
   }
}

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);
}

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 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 write(int n) {

   TRandom R;
   TStopwatch timer;

   TFile f1("mathcoreVectorIO_F.root","RECREATE");

   // create tree
   TTree t1("t1","Tree with new Float LorentzVector");

   XYZTVectorF *v1 = new XYZTVectorF();
   t1.Branch("LV branch","ROOT::Math::XYZTVectorF",&v1);

   timer.Start();
   for (int i = 0; i < n; ++i) {
      double Px = R.Gaus(0,10);
      double Py = R.Gaus(0,10);
      double Pz = R.Gaus(0,10);
      double E  = R.Gaus(100,10);
      v1->SetCoordinates(Px,Py,Pz,E);
      t1.Fill();
   }

   f1.Write();
   timer.Stop();
   std::cout << " Time for new Float Vector " << timer.RealTime() << "  " << timer.CpuTime() << std::endl;
   t1.Print();
}

void candleplotoption()
{
   TCanvas *c1 = new TCanvas("c1","Candle Presets",1200,800);
   c1->Divide(3,2);

   TRandom *randnum = new TRandom();
   TH2I *h1 = new TH2I("h1","Sin",18,0,360,300,-1.5,1.5);
   h1->GetXaxis()->SetTitle("Deg");
   float myRand;
   for (int i = 0; i < 360; i+=10) {
      for (int j = 0; j < 100; j++) {
         myRand = randnum->Gaus(sin(i*3.14/180),0.2);
         h1->Fill(i,myRand);
      }
   }
   for (int i = 1; i < 7; i++) {
      c1->cd(i);
      char str[16];
      sprintf(str,"candlex%d",i);
      TH2I * myhist = (TH2I*)h1->DrawCopy(str);
      myhist->SetTitle(str);
   }

   TCanvas *c2 = new TCanvas("c2","Candle Individual",1200,800);
   c2->Divide(4,4);
   char myopt[16][8] = {"0","1","11","21","31","30","111","311","301","1111","2321","12111","112111","212111","312111"};
   for (int i = 0; i < 15; i++) {
      c2->cd(i+1);
      char str[16];
      sprintf(str, "candlex(%s)",myopt[i]);
      TH2I * myhist = (TH2I*)h1->DrawCopy(str);
      myhist->SetTitle(str);
   }
}

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 Interpolation()
{
   ROOT::R::TRInterface &r=ROOT::R::TRInterface::Instance();
   //Creating points
   TRandom rg;
   std::vector<Double_t> x(10),y(10);
   for(int i=0;i<10;i++)
   {
      x[i]=i;
      y[i]=rg.Gaus();
   }

   r["x"]=x;
   r["y"]=y;


   // do plotting only in non-batch mode
   if (!gROOT->IsBatch() )  {

      r<<"dev.new()";//Required to activate new window for plot
      //Plot parameter. Plotting using two rows and one column
      r<<"par(mfrow = c(2,1))";

      //plotting the points
      r<<"plot(x, y, main = 'approx(.) and approxfun(.)')";

      //The function "approx" returns a list with components x and y
      //containing n coordinates which interpolate the given data points according to the method (and rule) desired.
      r<<"points(approx(x, y), col = 2, pch = '*')";
      r<<"points(approx(x, y, method = 'constant'), col = 4, pch = '*')";
   }
   else {
      r << "print('Interpolated points')";
      r << "print(approx(x,y,n=20))";
   }

   //The function "approxfun" returns a function performing (linear or constant)
   //interpolation of the given data.
   //For a given set of x values, this function will return the corresponding interpolated values.
   r<<"f <- approxfun(x, y)";
   //using approxfun with const method
   r<<"fc <- approxfun(x, y, method = 'const')";

   if (!gROOT->IsBatch() ) {
      r<<"curve(f(x), 0, 11, col = 'green2')";
      r<<"points(x, y)";

      r<<"curve(fc(x), 0, 10, col = 'darkblue', add = TRUE)";
      // different interpolation on left and right side :
      r<<"plot(approxfun(x, y, rule = 2:1), 0, 11,col = 'tomato', add = TRUE, lty = 3, lwd = 2)";
      r<<"dev.off()";//Required to close new window for plot
   }
   else {
      r << "x2=x+0.5";
      r << "print('Result of approxfun with default method')";
      r << "print(paste('x = ',x,'  f(x) = ',f(x2)))";
      r << "print('Result of approxfun with const method')";
      r << "print(paste('x = ',x,'  f(x) = ',fc(x2)))";
   }
}

void MixingResult::randomise () {//randomizer
  static TRandom donram(42); 
  std::cout << "Randomising " << getName() << "; old value "
	    << measurement;
  double delta = donram.Gaus()*error;
  measurement += delta;
  std::cout << " delta " << delta << " new value " << measurement << std::endl; 
}

void th2polyHoneycomb(){
   gStyle->SetCanvasPreferGL(true);
   TH2Poly *hc = new TH2Poly();
   hc->Honeycomb(0,0,.1,25,25);

   TRandom ran;
   for (int i = 0; i<30000; i++) {
      hc->Fill(ran.Gaus(2.,1), ran.Gaus(2.,1));
   }

   hc->Draw("gllego");
}

void makeDecalibCDB(Int_t firstRun, Int_t lastRun, Float_t decalib = 0.065)
{
  //Generates a random decalibration factors O(1)
  //to be applied in the anchor run simulations with raw:// .
  //If decalib<0, no decalibration generated, all factors=1.
  
  //Run range is [firstRun,lastRun] and gaussian sigma = decalib.
  //Author: Boris Polishchuk.
  
  AliCDBManager::Instance()->SetDefaultStorage("raw://");
  AliCDBManager::Instance()->SetRun(firstRun);

  TString emcPath("PHOS/Calib/EmcGainPedestals");
  AliCDBEntry* entryEmc = AliCDBManager::Instance()->Get(emcPath.Data(),-1);

  AliPHOSEmcCalibData* clb=0;

  if(entryEmc) clb = (AliPHOSEmcCalibData*)entryEmc->GetObject();
  else { printf("CDB entry not found. Exit.\n"); return; }
  
  if(!clb) { printf("Calibration parameters for PHOS EMC not found.\n"); return; }
  
  printf("\t\tEMC calibration object found: FirstRun=%d LastRun=%d Version=%d SubVersion=%d\n",
         entryEmc->GetId().GetFirstRun(), entryEmc->GetId().GetLastRun(),
         entryEmc->GetId().GetVersion(),entryEmc->GetId().GetSubVersion());
  
  
  TRandom rn;
  rn.SetSeed(0); //the seed is set to the current  machine clock

  Float_t adcChannelEmc;

  for(Int_t module=1; module<6; module++) {
    for(Int_t column=1; column<57; column++) {
      for(Int_t row=1; row<65; row++) {

	if(decalib<0.) adcChannelEmc = 1.;
	else
	  adcChannelEmc =rn.Gaus(1.,decalib);

        clb->SetADCchannelEmcDecalib(module,column,row,adcChannelEmc);
      }
    }
  }

  AliCDBManager::Instance()->SetDefaultStorage("local://./");
  AliCDBStorage* storage = AliCDBManager::Instance()->GetDefaultStorage();
  
  AliCDBMetaData *md = new AliCDBMetaData();
  AliCDBId id(emcPath.Data(),firstRun,lastRun);
  storage->Put(clb,id, md);

}

void JEC_fit_Uncertainty(int N)
{
  TF1 *func[1000];
  TFile *inf = new TFile("L3Graphs_test_Icone5.root");
  TGraphErrors *g = (TGraphErrors*)inf->Get("Correction_vs_CaloPt");
  TGraphErrors *vg[1000];
  int i,k;
  double x[20],y[20],ex[20],ey[20];
  double vx[20],vy[20],vex[20],vey[20];
  for(i=0;i<g->GetN();i++)
    {
      g->GetPoint(i,x[i],y[i]);
      ex[i]=g->GetErrorX(i);
      ey[i]=g->GetErrorY(i); 
    }  
  TRandom *rnd = new TRandom();
  rnd->SetSeed(0);
  for(k=0;k<N;k++)
    {
      for(i=0;i<g->GetN();i++)
        {	
          vx[i] = rnd->Gaus(x[i],ex[i]);
          //vx[i] = x[i];
          vy[i] = rnd->Gaus(y[i],ey[i]);
          vex[i] = ex[i];
          vey[i] = ey[i];
        }
      vg[k] = new TGraphErrors(g->GetN(),vx,vy,vex,vey);
      func[k] = new TF1("func","[0]+[1]/(pow(log10(x),[2])+[3])",1,2000);
      func[k]->SetParameters(1,3,6,5);
      vg[k]->Fit(func[k],"RQ");     	
    }
  
  TCanvas *c = new TCanvas("c","c");
  gPad->SetLogx();
  g->SetMarkerStyle(20);
  g->SetMaximum(3.5);
  g->Draw("AP");
  for(k=0;k<N;k++)
    {
      func[k]->SetLineColor(5);
      func[k]->SetLineWidth(1);
      cout<<func[k]->GetChisquare()<<endl;
      vg[k]->SetMarkerColor(2);
      vg[k]->SetLineColor(2);
      vg[k]->SetMarkerStyle(21);
      //if (func[k]->GetChisquare()<0.1)
        //vg[k]->Draw("sameP");
      func[k]->Draw("same");  
    }  	 
}  

double write(int n) {



  TRandom R;
  TStopwatch timer;


  TFile f1("mathcoreLV.root","RECREATE");

  // create tree
  TTree t1("t1","Tree with new LorentzVector");

  std::vector<ROOT::Math::XYZTVector>  tracks;
  std::vector<ROOT::Math::XYZTVector> * pTracks = &tracks;
  t1.Branch("tracks","std::vector<ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > >",&pTracks);

  double M = 0.13957;  // set pi+ mass

  timer.Start();
  double sum = 0;
  for (int i = 0; i < n; ++i) {
    int nPart = R.Poisson(5);
    pTracks->clear();
    pTracks->reserve(nPart);
    for (int j = 0; j < nPart; ++j) {
      double px = R.Gaus(0,10);
      double py = R.Gaus(0,10);
      double pt = sqrt(px*px +py*py);
      double eta = R.Uniform(-3,3);
      double phi = R.Uniform(0.0 , 2*TMath::Pi() );
      RhoEtaPhiVector vcyl( pt, eta, phi);
      // set energy
      double E = sqrt( vcyl.R()*vcyl.R() + M*M);
      XYZTVector q( vcyl.X(), vcyl.Y(), vcyl.Z(), E);
      // fill track vector
      pTracks->push_back(q);
      // evaluate sum of components to check
      sum += q.x()+q.y()+q.z()+q.t();
    }
    t1.Fill();
  }

  f1.Write();
  timer.Stop();
  std::cout << " Time for new Vector " << timer.RealTime() << "  " << timer.CpuTime() << std::endl;

  t1.Print();
  return sum;
}

TMap* CreateDCSAliasMap()
{
  // Creates a DCS structure
  // The structure is the following:
  //   TMap (key --> value)
  //     <DCSAlias> --> <valueList>
  //     <DCSAlias> is a string
  //     <valueList> is a TObjArray of AliDCSValue
  //     An AliDCSValue consists of timestamp and a value in form of a AliSimpleValue

  // In this example 6 aliases exists: DCSAlias1 ... DCSAlias6
  // Each contains 1000 values randomly generated by TRandom::Gaus + 5*nAlias

  TMap* aliasMap = new TMap;
  aliasMap->SetOwner(1);
  TRandom random;

	FILE *fp = fopen("./DCSValues.txt","r");

	char name[50];
	Float_t val;
	while(!(EOF == fscanf(fp,"%s %f",name,&val))){
		TObjArray* valueSet = new TObjArray;
		valueSet->SetOwner(1);

		TString aliasName=name;
		
		//printf("alias: %s\t\t",aliasName.Data());

		int timeStamp=10;
		
		
		if(aliasName.Contains("HV")) {
			for(int i=0;i<200;i++){
				dcsVal = new AliDCSValue((Float_t) (val+random.Gaus(0,val*0.1)), timeStamp+10*i);
				valueSet->Add(dcsVal);
			}
		} else {
			for(int i=0;i<2;i++){
				AliDCSValue* dcsVal = new AliDCSValue((UInt_t) (val), timeStamp+10*i);
				valueSet->Add(dcsVal);
			}
		}
		
		aliasMap->Add(new TObjString(aliasName), valueSet);

	}
	fclose(fp);
  return aliasMap;
}

void drawing_pion_decay(){

	double Pi_lifetime = 0.26033e-6; //s
	double c = 3e8;//m/s
	TCanvas *c1 = new TCanvas("test", "test");
	TRandom *r = new TRandom();
	

	TView *view = TView::CreateView(1, 0, 0);
	view->ShowAxis();
	view->SetRange(-5, -5, 0, 5, 5, 10);
	
	

	for (int i = 0; i < 100; ++i)
	{
		double px = 0.3;
		double py = 0.3;
		double pz = -1; // GeV
		double energy = sqrt(px*px+py*py+pz*pz+M_pion*M_pion);
		
		//life time
		double t = r->Exp(Pi_lifetime);
		//decay length
		double gamma = energy/M_pion;
		double length = c*t*gamma/1e3;
		// decay position
		double vx = 0;
		double vy = 0;
		double vz = 10 - length;

		TLorentzVector pion(px, py, pz, energy);

		double decay_particle_mass[2] = {M_muon, M_neu};
		TGenPhaseSpace event;
		event.SetDecay(pion, 2, decay_particle_mass);
		event.Generate();

		TLorentzVector Muon = *(event.GetDecay(0));
		TLorentzVector Neu = *(event.GetDecay(1));
		
		plot_particle(Muon, vx, vy, vz, 2);
		plot_particle(Neu, vx, vy, vz, 4);


	}


}

double error(double a, double b) {

 TRandom rndm;
 rndm.SetSeed(12345);
 
 if (a<=0.000000001) return 0.;

 double eff = a/b;
 int n0 = rndm.Poisson((a+b));
 // int n1 = rndm.Binomial(n0,eff);

 double err=sqrt(eff*(1-eff)/n0);

 return err;
}

void fmpmt(void){

  TTree *t = new TTree("tree","fine mesh");
  double adc;
  t->Branch("adc",&adc,"adc/D");

  double Mult(double seed){
    do{
      double _fac=2.79360;
      TRandom _rand;
      _rand.SetSeed(0);
      double _judge = _rand.Uniform(0,1);
      double _x = _rand.Uniform(0.001,10);
    } while (TMath::Gaus(_x, _fac, _fac/20) < _judge);
    return seed*_x;
  }


  double prob = 0.7;
  double npe = 0.1;
  double judge, tmp;
  TRandom rand;
  rand.SetSeed(0);

  for(int i=0; i<4000; ++i){
    if(i%100==0){
      cout<<"=====  "<<i<<"  ====="<<endl;
    }
    do{
      tmp = rand.Uniform(0,4);
      judge = rand.Uniform(0,1);
    } while (TMath::Poisson(tmp, npe) < judge);
    for(int j=0; j<19; ++j){
      if(j==0){
        judge = rand.Uniform(0,1);
        if(judge>prob){
          continue;
        }
      }
      tmp=Mult(tmp);
    }
    adc=tmp;
    t->Fill();
  }
  t->Draw("adc");


}

void tv3Write() {
   //creates the Tree
   TVector3 *v = new TVector3();
   TVector3::Class()->IgnoreTObjectStreamer();
   TFile *f = new TFile("v3.root","recreate");
   TTree *T = new TTree("T","v3 Tree");
   T->Branch("v3","TVector3",&v,32000,1);
   TRandom r;
   for (Int_t i=0;i<10000;i++) {
      v->SetXYZ(r.Gaus(0,1),r.Landau(0,1),r.Gaus(100,10));
      T->Fill();
   }
   T->Write();
   T->Print();
   delete f;
}

void TestDistribution(TRandom& aRandom, TInt aSamples)
	{
	TUint32* data = new TUint32[aSamples];
	test_NotNull(data);

	TInt i;
	TReal mean = 0.0;
	for (i = 0 ; i < aSamples ; ++i)
		{
		data[i] = aRandom.Next();
		mean += (TReal)data[i] / aSamples;
		}

	TReal sum2 = 0.0;
	for (i = 0 ; i < aSamples ; ++i)
		{
		TReal d = (TReal)data[i] - mean;
		sum2 += d * d;
		}
	TReal variance = sum2 / (aSamples - 1);

	test.Printf(_L("  mean == %f\n"), mean);
	test.Printf(_L("  variance == %f\n"), variance);

	delete [] data;
	}

void generate_imt_tree()
{
  // Create the file and the tree
  TFile hfile("ttree_read_imt.root", "RECREATE", "File for IMT test");
  TTree tree("TreeIMT", "TTree for IMT test");

  int nvbranches = 50, nabranches = 50;
  int nentries = 1000, nvelems = 100;
  std::vector<std::vector<Double_t>> vbranches(nvbranches);
  std::vector<A> abranches(nabranches);

  // Create the tree branches
  for (int i = 0; i < nvbranches; ++i) {
    vbranches[i] = std::vector<Double_t>(nvelems);

    std::string branchname("Vbranch");
    branchname += std::to_string(i);
    branchname += std::string("."); // make the top-level branch name be included in the sub-branch names

    tree.Branch(branchname.c_str(), &vbranches[i]);
  }
  for (int i = 0; i < nabranches; ++i) {
    std::string branchname("Abranch");
    branchname += std::to_string(i);
    branchname += std::string("."); // make the top-level branch name be included in the sub-branch names

    tree.Branch(branchname.c_str(), &abranches[i]);
  }

  // Fill the tree
  TRandom rand;
  for (int i = 0; i < nentries; i++) {
    for (int i = 0; i < nvbranches; ++i) {
      for (int j = 0; j < nvelems; ++j) {
        vbranches[i][j] = rand.Uniform();
      }
    }
    for (int i = 0; i < nabranches; ++i) {
      abranches[i].Build();
    }
    Int_t nb = tree.Fill();
  }

  // Write the file 
  hfile.Write();
  hfile.Close();
}