void basic() {
//  Read data from an ascii file and create a root file with an histogram and an ntuple.
//   see a variant of this macro in basic2.C
//Author: Rene Brun


// read file $ROOTSYS/tutorials/tree/basic.dat
// this file has 3 columns of float data
   TString dir = gSystem->UnixPathName(__FILE__);
   dir.ReplaceAll("basic.C","");
   dir.ReplaceAll("/./","/");
   ifstream in;
   in.open(Form("%sbasic.dat",dir.Data()));

   Float_t x,y,z;
   Int_t nlines = 0;
   TFile *f = new TFile("basic.root","RECREATE");
   TH1F *h1 = new TH1F("h1","x distribution",100,-4,4);
   TNtuple *ntuple = new TNtuple("ntuple","data from ascii file","x:y:z");

   while (1) {
      in >> x >> y >> z;
      if (!in.good()) break;
      if (nlines < 5) printf("x=%8f, y=%8f, z=%8f\n",x,y,z);
      h1->Fill(x);
      ntuple->Fill(x,y,z);
      nlines++;
   }
   printf(" found %d points\n",nlines);

   in.close();

   f->Write();
}

void MakeTTreeFromBinfit() {
//  Read data from an ascii file and create a root file with an histogram and an ntuple.
//   see a variant of this macro in basic2.C
      

// read file $ROOTSYS/tutorials/tree/basic.dat
// this file has 3 columns of float data
   TString dir = gSystem->UnixPathName(gInterpreter->GetCurrentMacroName());
   dir.ReplaceAll("MakeTTreeFromBinfit.C","");
   dir.ReplaceAll("/./","/");
   ifstream in;
   in.open(Form("%sBinfit_results_reco_ptGT8.txt",dir.Data()));

   Float_t fracd,fracs, fracn, chi;
   Int_t nlines;
   string result;
   TFile *f = new TFile("Binfit_results_spsdpsnlo_ptYdy_reco8gev.root","RECREATE");
   TNtuple *ntuple = new TNtuple("ntuple","data from ascii file","fracd:fracs:fracn:chi");

   while (1) {
      in >> fracd >> fracs >> fracn >> chi;
      if (!in.good()) break;
      if (nlines < 5) printf("fracd=%8f, chi=%8f \n",fracd,chi);
      ntuple->Fill(fracd,fracs,fracn,chi);
      nlines++;
   }
   printf(" found %d points\n",nlines);

   in.close();

   f->Write();
}

double  Getmean(TString sys1, TString sys2, float cent1, float cent2){

  gROOT->SetStyle("Plain");
  TFile *infile = TFile::Open(Form("glau_%s%s_ntuple_1M.root",sys1.Data(),sys2.Data()));
  infile->cd();
  TH2F *pEcc2B = new TH2F("pEcc2B","#epsilon_{2} vs impact paramter in 200 GeV Collisions from Glauber MC;B;#epsilon_{n}",25200,-0.5,251.5,100,0,1);
  if(sys2.Contains("smeared")){
  TNtuple *nt = (TNtuple*)infile->Get("nt");
  nt->Project("pEcc2B","Ecc3G:B");
  }
  else{
  TNtuple *nt = (TNtuple*)infile->Get(Form("nt_%s_%s",sys1.Data(),sys2.Data()));
  nt->Project("pEcc2B","Ecc2:B");
  }
  TH1D* hB = (TH1D*)pEcc2B->ProjectionX("hB",0,-1);
  for(int ibin=0;ibin<hB->GetNbinsX();ibin++){
      if(hB->Integral(0,ibin) >= cent1 /100. * hB->Integral()){ int bin1 = ibin; break;}
  }
  for(int ibin=0;ibin<hB->GetNbinsX();ibin++){
      if(hB->Integral(0,ibin) >= cent2 /100. * hB->Integral()){ int bin2 = ibin; break;}
  }
  cout<< bin1 << "\t" << bin2 << endl;
  TH1D* hEcc2 = (TH1D*)pEcc2B->ProjectionY("hEcc2",bin1,bin2);
//  hEcc2->Draw();
  cout << hEcc2 -> GetEntries() << endl;
  double mean = hEcc2->GetMean();
  return mean;
}

		double pValueFromPoint2PointDissimilarity(IDataSet * data, IDataSet * mcData)
		{
            std::cout << "GC: calculating T stat for data" << std::endl;
            double T = calculateTstatistic( data, mcData );
			char buffer[20];
			sprintf( buffer, "Tdata = %f", T );
			cout << buffer << endl;

			int nPerm = 25;
			vector<double> Tvalues = permutation( data, mcData, nPerm );

            int count = 0;
			vector<double>::iterator Titer;
			for ( Titer = Tvalues.begin(); Titer != Tvalues.end(); ++Titer ){
				if ( T < *Titer ) count++;
			}

			double pvalue = count/double(nPerm);

	    string fileName = ResultFormatter::GetOutputFolder();	
	    fileName.append("/tvalues.root");
	    	
            TFile * outputFile = new TFile(fileName.c_str(), "RECREATE");
            TNtuple * ntuple = new TNtuple("tvalues", "tvalues", "T:Tdata:pvalue");
            for ( int i = 0; i < nPerm; i++ ) ntuple->Fill(Tvalues[i], T, pvalue);
            ntuple->Write();
            outputFile->Close();
            delete outputFile;

       		return pvalue;
        }

void read_ntuple_from_file(){

  int i,j,k,n;

  TFile *in = new TFile("ntupleoutputsample.root");
  TNtuple *data = (TNtuple*) in->GetObjectChecked("data","TNtuple");

  double pot,cur,temp,pres;
  float *row_content; //Must necessarily be float and not a double... WHY?

  TH1D *histo = new TH1D("histo","HISTO",100,0,10);

  cout << "Potential\tCurrent\tTemperature\tPressure" << endl;
  for(i=0;i<data->GetEntries();++i){
      data->GetEntry(i);
      row_content = data->GetArgs();
      pot = row_content[0];
      cur = row_content[1];
      temp = row_content[2];
      pres = row_content[3];
      cout << pot << "\t" << cur << "\t" << temp << "\t" << pres << endl;
      histo->Fill(pot);
  }

  histo->Draw();
}

void mkTree() {

   std::string file("lumis.dat");

   cout << file << endl;

   ifstream in;
   in.open("lumis.dat"); 

   Float_t run,ls,lumiDelivered, lumiReported;
   Int_t nlines = 0;
   TFile *f = new TFile("lumis.root","RECREATE");
   TNtuple *ntuple = new TNtuple("ntuple","data from ascii file","run:ls:lumiDelivered:lumiReported");

   while (1) {
      in >> run >> ls >> lumiDelivered >> lumiReported;
      if (!in.good()) break;
      if (nlines < 5) printf("run=%8f, ls=%8f, lumiDelivered=%8f, lumiReported=%8f\n",run, ls, lumiDelivered, lumiReported);
      ntuple->Fill(run, ls, lumiDelivered, lumiReported);
      nlines++;
   }
   printf(" found %d points\n",nlines);

   in.close();

   f->Write();

   gROOT->ProcessLine(".q");
}

void makeTTree(){

	//	Declaring 
	Double_t num;
	Double_t xNumBins, yNumBins;

	//	Setting up Canvas to put plot the 2D histogram (from Tamii) and the 2 Bracnh Tree (created by this program)
	TCanvas *canvas = new TCanvas("canvas","Canvas");
	canvas->Divide(2,1);
	canvas->cd(1);


	//	Opening up histogram file from tamii's analyzer
	//	Remember to do h2root before this program runs.
	TFile *f1 = new TFile("run6106_test44.root");
	TH2F *hist2D = (TH2F*)f1->Get("h159");
	//	Plotting the 2D histogram on the first part of the Canvas
	hist2D->Draw("COLZ");
	
	//	Going through 2D hist and getting the total number of x and y bins.
	xNumBins = hist2D->GetNbinsX();
	yNumBins = hist2D->GetNbinsY();
	
	//	Creating a new file to save the 2 branch tree to.
	TFile *t1 = new TFile("~/e329/root/13C/runs/run3014.root","recreate");

	//	This is the actual tree being created.
	//	!!! For you YingYing, instead of x vs theta, you would have theta(fp or target) vs Yfp !!!
	TNtuple *DATA = new TNtuple("DATA","Xpos vs Theta","Xpos:Theta");
	

	//	This is the main nested for loop that goes through each bin in the 2D hist 
	//	takes the number of counts in that (x,y)bin and creates that many events
	//	in the 2 branch tree  
	for (int binx = 0; binx < xNumBins+2; binx++){

		for (int biny = 0; biny < yNumBins+2; biny++){

			//getting the number of events in a particular (x,y)bin
			num = hist2D->GetBinContent(binx,biny);
			double x = hist2D->GetXaxis()->GetBinCenter(binx); 
			double y = hist2D->GetYaxis()->GetBinCenter(biny); 
			//cout << num << " at x= " << x << ", y= " << y << endl;

			//Filling the Tree 'num' times with xbin and ybin events
			for (int i = 0; i < num; i++){
				DATA->Fill(x,y);
			}
		}
	}
	
	//Saving the Tree that was just created and filled
	t1->Write();

	//Ploting the Tree in the second part of the canvas.
	canvas->cd(2);
	DATA->Draw("Theta:Xpos>>(2000,-600,600,1000,-3,3)","","COLZ");
}

	void jackknife( I_XMLConfigReader * xmlFile, MinimiserConfiguration * theMinimiser, 
		FitFunctionConfiguration * theFunction, ParameterSet* argumentParameterSet, vector<string> CommandLineParam, int start, int stop )
	{
		cout << "Starting JackKnife" << endl;
		PDFWithData  * pdfAndData  = xmlFile->GetPDFsAndData()[0];
		ParameterSet * parSet      = xmlFile->GetFitParameters( CommandLineParam );
		PhaseSpaceBoundary * phase = xmlFile->GetPhaseSpaceBoundaries()[0];
		pdfAndData->SetPhysicsParameters( parSet );
			
		MemoryDataSet * dataset = (MemoryDataSet*) pdfAndData->GetDataSet();
		vector<IDataSet*> data;
		data.push_back(dataset);
	
		IPDF * pdf = pdfAndData->GetPDF();
		vector<IPDF*> vectorPDFs;
		vectorPDFs.push_back(pdf);
	
		// Repeat nominal fit to get the nominal value of the physics parameter	
		FitResult * nominal = FitAssembler::DoFit( theMinimiser, theFunction, argumentParameterSet, vectorPDFs, data, xmlFile->GetConstraints() );
		double nominal_value = nominal->GetResultParameterSet()->GetResultParameter("tau")->GetValue();
		cout << nominal_value << endl;
		double jackknifed_value = 0.;

		MemoryDataSet * subset = new MemoryDataSet( phase );
		int nData = dataset->Yield();	
                string fileName = ResultFormatter::GetOutputFolder();
                fileName.append("/jackknife.root");         
		TFile * outputFile = new TFile("jackknife.root", "RECREATE");	
                TNtuple * jack = new TNtuple("jack", "jackknifed - nominal", "diff_jackknifed_nominal:reco_time:true_time");

		double reco_time = 0.;
		double true_time = 0.;

		for ( int i = start; i < stop; i++ )
		{
			for ( int j = 0; j < nData; j++ )
			{
				if ( j == i ) {
					reco_time = dataset->GetDataPoint( j )->GetObservable( "time" )->GetValue();
					true_time = dataset->GetDataPoint( j )->GetObservable( "truetime" )->GetValue();
					continue;
				}
				subset->AddDataPoint( dataset->GetDataPoint( j ) );
			}
			cout << "Creating subset " << i << " containing " << subset->Yield() << " candidates" << endl;
			vector<IDataSet*> vectorData;
			vectorData.push_back(subset);
			FitResult * result = FitAssembler::DoFit( theMinimiser, theFunction, argumentParameterSet, vectorPDFs, vectorData, xmlFile->GetConstraints() );
			if ( result->GetFitStatus() == 3 ) {
				jackknifed_value = result->GetResultParameterSet()->GetResultParameter("tau")->GetValue();
				jack->Fill( (Float_t)(jackknifed_value - nominal_value), (Float_t)reco_time, (Float_t)true_time);
			}
			subset->Clear();
			delete result;
		}
		outputFile->Write();
		outputFile->Close();
	}

int execQuery() 
{
   TNtuple *ed = new TNtuple("ed","ed","pt:val:stat");
   ed->SetScanField(0);
   ed->Fill(1.0, 5.0, 4.0);
   TSQLRow *row = ed->Query("pt:val:stat")->Next();
   fprintf(stdout,"1: %s 2: %s 3: %s\n",row->GetField(0),row->GetField(1),row->GetField(2));
   return 0;
}

void Example2(){
    gROOT->Reset();
    gROOT->Clear();
    gROOT->SetStyle("Plain");

    gStyle->SetTextSize(0.01908148);
    gStyle->SetTitleFontSize(0.07);
    gStyle->SetOptTitle(1);
    gStyle->SetOptStat(1110);
    gStyle->SetOptFit(1111);
    gStyle->SetTitleXOffset(1.1);
    gStyle->SetTitleYOffset(1.55);
    gStyle->SetPadTopMargin(0.15);
    gStyle->SetPadBottomMargin(0.15);
    gStyle->SetPadLeftMargin(0.15);

    // select the one of the versus : Hf energy, Centrality, # Charged , # Tracks
    int Vselect = 1;

    // number of Trigger you will use , 
    const int Ntr = 2; 
    //if you want to use more than two trigger please specify the Ntr2 for dividivg the canvas
    const int Ntr2 = 1;
    TFile* inf = new TFile("openhlt2.root");

    string triggers[6] = {"L1Tech_BSC_minBias_threshold1.v0","L1Tech_BSC_minBias_threshold2.v0","L1_SingleJet30","L1_TripleJet30","L1_QuadJet15","L1_DoubleJet70"};
    string vers[4] = {"hiHF","hiBin","Ncharged","hiNtracks"};

    double bins [4] = {170,40,240,140};
    double limits[4] = {170000,40,24000,1400};

    TCanvas* c1 = new TCanvas("c1","c1",800,400);
    c1->Divide(Ntr,Ntr2);

    TNtuple* nt = (TNtuple*)inf->Get("HltTree");

    TProfile* p1[Ntr];

    for(unsigned int i =0; i<Ntr ; i++){
        c1->cd(i+1);

        p1[i] = new TProfile("p1",Form(";%s;%s",vers[Vselect].data(),triggers[i].data()),bins[1],0,limits[1]);
        nt->SetAlias("trigger",triggers[i].data());

        nt->SetAlias("versus",vers[Vselect].data());
        nt->Draw("trigger:versus>>p1","","prof");

    }

    c1->Print(Form("%s_vsCent.gif",triggers[Vselect].data()));

}

	void writeVolumeWeightedProfile() {
		std::cout << "** Write Volume Weighted Profile" << std::endl;

		const size_t steps = 50;
		const float rMin = 5;
		const float rMax = ComaRadiusKpc;
		const size_t nPoints = 10000;

#if GADGET_ROOT_ENABLED
		TFile *file = new TFile("coma_profile_volume_weighted.root", "RECREATE",
				"Rho Weighted");
		if (file->IsZombie())
		throw std::runtime_error(
				"Root output file cannot be properly opened");
		TNtuple *ntuple = new TNtuple("events", "Rho", "r:B");
#endif
		std::ofstream out("coma_profile_volume_weighted.csv");
		out << "r B" << std::endl;

		float stepLog = (log10(rMax) - log10(rMin)) / (steps - 1);
		for (size_t i = 0; i < steps; i++) {
			std::cout << ".";
			std::cout.flush();

			float r = pow(10, log10(rMin) + stepLog * i);

			float avgB = 0;
			for (size_t i = 0; i < nPoints; i++) {
				Vector3f p = randomUnitVector();
				p *= r;
				p += ComaPositionKpc;
				Vector3f b;
				bool isGood = dmf->getField(p, b);
				avgB += b.length(); // / nPoints;
			}
			avgB /= nPoints;

#if GADGET_ROOT_ENABLED
			ntuple->Fill(r, avgB);
#endif
			out << r << " " << avgB << std::endl;

		}

#if GADGET_ROOT_ENABLED
		file->Write();
		file->Close();
#endif
		std::cout << std::endl;
	}

double getUncertaintyValue(double m0,double m12, int channel, string valuename, TFile * input)
{
  bool DEBUGuncertaintyValue=false;
/*
  string buffer = doubletostr(m0)+doubletostr(m12)+inttostr(channel)+valuename;
  map <std::string, double>::iterator iter;
  iter=UncertaintyValueBuffer.find(buffer);
  if (bufferValue && iter!=UncertaintyValueBuffer.end()) return iter->second;
*/
  double returnvalue2 ;
  double min=-10;
  double max=+10;

  //cout << "v" ;
  string leaveName="";
  if (valuename=="K") leaveName="K";
  if (valuename=="crossSection") {leaveName="crossSection"; min=0; max=100000;};
  if (valuename=="PDFUncertainty") leaveName="relUncPDF";
  if (valuename=="ScaleUncertainty2Q") leaveName="relScaleUnc2Q";
  if (valuename=="ScaleUncertainty12Q") leaveName="relScaleUncHalfQ";
  if (leaveName=="") throw "no leave found for that value";

  string commandstring=leaveName+">>h1";
  string conditionstring= "finalState>" + doubletostr(channel-0.1)  +"&& finalState<"+ doubletostr(channel+0.1) +"&& m0>" + doubletostr(m0-0.5) + " && m0<"+doubletostr(m0+0.5) +" && m12>" + doubletostr(m12-0.5) + " && m12<"+doubletostr(m12+0.5);

  if (DEBUGuncertaintyValue){
  	  cout << "	 commandstring " << commandstring << endl;
  	  cout << "	 conditionstring " << conditionstring << endl;
  	  cout << "	 leaveName " << leaveName << endl;
  	  cout << "	 valuename :" << valuename << " channel: " << channel << " m0: " << m0 << " m12: " << m12   << endl;
  	  }

  TNtuple *ntuple = (TNtuple*)file_robin->Get("SignalUncertainties");
  if (ntuple==0) throw "Ntuple not found";

  if (DEBUGuncertaintyValue) cout << "(min,max) "<<min << " "<<max<<endl;
  TH1F * h1 = new TH1F("h1","h1",100000,min,max);
  ntuple->Draw(commandstring.c_str(),conditionstring.c_str());

  if (DEBUGuncertaintyValue) cout << "mean " << h1->GetMean()<<endl;
  
  returnvalue2 = fabs (h1->GetMean());
//  UncertaintyValueBuffer[buffer]=returnvalue2;
  
  delete h1;

  return returnvalue2;
}

void Hijing()
{
  TFile *fin = TFile::Open("/sphenix/user/pinkenbu/jade/hijingdat2.root");
  gROOT->cd();
  TH1 *hjbkg = new TH1F("hjbkg","Hijing Background",100,0.,0.5);
  TNtuple *de = (TNtuple *) fin->Get("de");
  de->Project("hjbkg","dtotal","(ID<=2)*edep/250.");

       char fname [100];
       sprintf(fname, "Gamma_Neutron_Hijing_Energy_Graphs.root");
       TFile *fout = TFile::Open(fname,"UPDATE");
       hjbkg->Write();
       fout->Write();
       fout->Close();
  fin->Close();
}

void parallelcoord() {

   TNtuple *nt = NULL;

   Double_t s1x, s1y, s1z;
   Double_t s2x, s2y, s2z;
   Double_t s3x, s3y, s3z;
   r = new TRandom();;

   new TCanvas("c1", "c1",0,0,800,700);
   gStyle->SetPalette(1);

   nt = new TNtuple("nt","Demo ntuple","x:y:z:u:v:w");

   for (Int_t i=0; i<20000; i++) {
      r->Sphere(s1x, s1y, s1z, 0.1);
      r->Sphere(s2x, s2y, s2z, 0.2);
      r->Sphere(s3x, s3y, s3z, 0.05);

      generate_random(i);
      nt->Fill(r1, r2, r3, r4, r5, r6);

      generate_random(i);
      nt->Fill(s1x, s1y, s1z, s2x, s2y, s2z);

      generate_random(i);
      nt->Fill(r1, r2, r3, r4, r5, r6);

      generate_random(i);
      nt->Fill(s2x-1, s2y-1, s2z, s1x+.5, s1y+.5, s1z+.5);

      generate_random(i);
      nt->Fill(r1, r2, r3, r4, r5, r6);

      generate_random(i);
      nt->Fill(s1x+1, s1y+1, s1z+1, s3x-2, s3y-2, s3z-2);

      generate_random(i);
      nt->Fill(r1, r2, r3, r4, r5, r6);
   }
   nt->Draw("x:y:z:u:v:w","","para",5000);
   TParallelCoord* para = (TParallelCoord*)gPad->GetListOfPrimitives()->FindObject("ParaCoord");
   para->SetDotsSpacing(5);
   TParallelCoordVar* firstaxis = (TParallelCoordVar*)para->GetVarList()->FindObject("x");
   firstaxis->AddRange(new TParallelCoordRange(firstaxis,0.846018,1.158469));
   para->AddSelection("violet");
   para->GetCurrentSelection()->SetLineColor(kViolet);
   firstaxis->AddRange(new TParallelCoordRange(firstaxis,-0.169447,0.169042));
   para->AddSelection("Orange");
   para->GetCurrentSelection()->SetLineColor(kOrange+9);
   firstaxis->AddRange(new TParallelCoordRange(firstaxis,-1.263024,-0.755292));
}

int main() {
  std::vector<std::pair<std::string, std::string>> sames = {
    std::make_pair("/tmp/rcas/tmp/f.root", "n"),
    std::make_pair("/tmp/rcas/tmp/nt.root", "nt")
  };
  std::map<std::string, std::future<FiledNtuple>> nt_map;
  std::vector<std::future<FiledNtuple>> nts;

  for (int n=0; n<2000; ++n) {
    for (unsigned i=0; i<2; ++i) {
      nts.push_back(
        std::async([i, sames]() {
          TFile* f(TFile::Open(sames[i].first.c_str()));
          TNtuple* nt;
          f->GetObject(sames[i].second.c_str(), nt);

          return FiledNtuple {f, nt};
        }));
      nt_map[sames[i].first] =
        std::async([i, sames]() {
          TFile* f(TFile::Open(sames[i].first.c_str()));
          TNtuple* nt;
          f->GetObject(sames[i].second.c_str(), nt);

          nt->SetDirectory(nullptr);
          f->Close();
          f = nullptr;

          return FiledNtuple {f, nt};
        }
      );
    }
  }

  for (auto& nti : nt_map) {
    std::cout << nti.first <<"\t" << nti.second.get().nt->GetEntries() << std::endl;
  }

  std::cout << "##############################\n";

  for (auto it=nts.rbegin(); it!=nts.rend(); ++it) {
    auto nt = (*it).get();
    std::cout << nt.nt->GetName() <<" "<< nt.nt->GetEntries() << std::endl;
  }
}

void decayang1()
{
   gROOT->Reset();

  float E, p, px, py, pz;
  float phi,cos_theta;

  const float PI=3.14159;
  const float mpi=0.13957018;
  const float mrho=0.77549; 
  TLorentzVector P1, P2;

  TFile *f = new TFile("decayang.root","RECREATE");

  TNtuple *cmspi = new TNtuple("cmspi","Pi decay to two photons in CMS","px1:py1:pz1:E1:px2:py2:pz2:E2");
  TNtuple *cmsrho = new TNtuple("cmsrho","Rho decay to two pions in CMS","px1:py1:pz1:E1:px2:py2:pz2:E2");
  TNtuple *boostpi = new TNtuple("boostpi","Pi decay to two photons in lab frames","b1px1:b1py1:b1pz1:b1e1:b1px2:b1py2:b1pz2:b1e2");


  for (int i=0; i<10000; i++) {
    phi = 2.0 * PI *(gRandom->Rndm());                // generate phi
    cos_theta = 2.0 * (gRandom->Rndm())-1;            //generate cos(theta)

     //Start with the pi to photon decays
    E=mpi/2.0;
    p=E;
    px=p*cos(phi)*sqrt(1.0-pow(cos_theta,2));
    py=p*sin(phi)*sqrt(1.0-pow(cos_theta,2));
    pz=p*cos_theta;
    cmspi->Fill(px, py, pz, E, -px, -py, -pz, E);
    
    

     //Now do the rho to pi decays
    E=mrho/2.0;
    p=sqrt(E*E-mpi*mpi);
    px = p * cos(phi) * sqrt( 1.0 - pow(cos_theta,2));
    py = p * sin(phi) * sqrt( 1.0 - pow(cos_theta,2));
    pz = p * cos_theta;
    cmsrho->Fill(px, py, pz, E, -px, -py, -pz, E);
  }

  f->Write(); 
}

void CheckDALIDet(char * ridffile="test.ridf"){

  gSystem->Load("libanacore.so");

  TArtEventStore *estore = new TArtEventStore();
  estore->Open(ridffile);
  TArtRawEventObject *rawevent = estore->GetRawEventObject();

  TFile *fout = new TFile("dataid.root","RECREATE");
  TNtuple *ntp = new TNtuple("ntp","ntp","geo:ch:val");

  int neve = 0;
  while(estore->GetNextEvent() && neve<100){

    int eve_number=rawevent->GetEventNumber();
    for(int i=0;i<rawevent->GetNumSeg();i++){
      TArtRawSegmentObject *seg = rawevent->GetSegment(i);
      int device = seg->GetDevice();
      int fp = seg->GetFP();
      int detector = seg->GetDetector();
      int module = seg->GetModule();
      if(DALI==device&&DALIA==detector){
	cout << "    seg:"<< i <<" dev:"<< device 
	     << " fp:"<<fp<< " det:"<<detector<< " mod:"<<module
	     << " #data=" << seg->GetNumData() << endl;
	for(int j=0;j<seg->GetNumData();j++){
	  TArtRawDataObject *d = seg->GetData(j);
	  int geo = d->GetGeo();
	  int ch = d->GetCh();
	  unsigned int val = d->GetVal();
	  cout << "       geo:" << geo 
	       << " ch:" << ch << " val:" << val << endl;
	  ntp->Fill((float)geo,(float)ch,(float)val);
	}
      }
    }

    rawevent->Clear();
    neve ++;
  }
  fout->Write();
  fout->Close();

}

void RIDF2DataId(){

  gSystem->Load("libanacore.so");

  TArtEventStore *estore = new TArtEventStore();
  estore->Open("test.ridf");
  TArtRawEventObject *rawevent = estore->GetRawEventObject();

  TFile *fout = new TFile("dataid.root","RECREATE");
  TNtuple *ntp = new TNtuple("ntp","ntp","sys:cat:det:dtype:val");

  int neve = 0;
  while(estore->GetNextEvent() && neve < 10000){

    for(int i=0;i<rawevent->GetNumSeg();i++){
      TArtRawSegmentObject *seg = rawevent->GetSegment(i);
      int device = seg->GetDevice();
      int fp = seg->GetFP();
      int detector = seg->GetDetector();
      int module = seg->GetModule();

      for(int j=0;j<seg->GetNumData();j++){
	TArtRawDataObject *d = seg->GetData(j);
	int geo = d->GetGeo(); 
	int ch = d->GetCh();
	int val = d->GetVal(); 
	int cat = d->GetCategoryID(); 
	int det = d->GetDetectorID(); 
	int id = d->GetDatatypeID(); 

	ntp->Fill((float)detector,(float)cat,(float)det,(float)id,(float)val);
      }
    }

    estore->ClearData();
    neve ++;
  }

  fout->Write();
  fout->Close();

}

int main(){

  std::vector<double> initparams;
  initparams.push_back(5.0);

  FFF::FitSimple fit(&myNLL,&myMC);
  FFF::DataWrapper* data = fit.GenerateData(initparams);
  std::cout << "The correct answer is 5.0, with a likelihood of 0.0" << std::endl;


  std::vector<double> results = fit.MigradFit(initparams,data);
  std::cout << "Migrad: " << results[0] << " L=" << results[1] << std::endl;

  TNtuple *lspace = fit.MCMCMapLikelihood(initparams,data);
  TFile f("lspace.root","RECREATE");
  lspace->Write();
  f.Close();
  results = fit.MCMCFit(initparams,data);
  std::cout << "MCMC: " << results[0] << " L=" << results[1] << std::endl;

  results = fit.SAnnealFit(initparams,data);
  std::cout << "SAnneal: " << results[0] << " L=" << results[1] << std::endl;

  // You can also do normal minuit2 stuff
  ROOT::Minuit2::FCNBase* migradfunc = fit.GetMinuit2FCNBase(data);
  std::vector<double> verrors;
  for (size_t i=0;i<initparams.size();i++)
    verrors.push_back(0.0); //FIXME
  ROOT::Minuit2::MnUserParameters mnParams(initparams,verrors);
  ROOT::Minuit2::MnMigrad migrad(*migradfunc,mnParams);  
  ROOT::Minuit2::FunctionMinimum theMin = migrad(); //FIXME maxfcn, tol
  ROOT::Minuit2::MnUserParameters migradresult = theMin.UserParameters();
  ROOT::Minuit2::MnMinos minos(*migradfunc,theMin);
  std::pair<double, double> errors = minos(0);

  std::cout << "Minos errors:" << std::endl;
  std::cout << migradresult.Params()[0] << " +" << errors.second << " -" << -1*errors.first << std::endl;

  return 0;

}

TNtuple* LikelihoodSpace::get_contour(float delta) {
  TNtuple* contour = (TNtuple*) this->samples->Clone("lscontour");
  contour->Reset();

  // Get list of branch names
  std::vector<std::string> names;
  for (int i=0; i<this->samples->GetListOfBranches()->GetEntries(); i++) {
    std::string name = this->samples->GetListOfBranches()->At(i)->GetName();
    if (name == "likelihood") {
      continue;
    }
    names.push_back(name);
  }

  float* params_branch = new float[names.size()];
  for (size_t i=0; i<names.size(); i++) {
    this->samples->SetBranchAddress(names[i].c_str(), &params_branch[i]);
  }

  float ml_branch;
  this->samples->SetBranchAddress("likelihood", &ml_branch);

  // Build a new TNtuple with samples inside the contour
  float* v = new float[names.size() + 1];
  for (int i=0; i<this->samples->GetEntries(); i++) {
    this->samples->GetEntry(i);
    if (ml_branch < this->ml + delta) {
      for (size_t j=0; j<names.size(); j++) {
        v[j] = params_branch[j];
      }
      v[names.size()] = ml_branch;
      contour->Fill(v);
    }
  }

  this->samples->ResetBranchAddresses();
  delete[] v;

  return contour;
}