//Newest version of the toy tests, some changes in the files and stuff, so
//replace old stuff rather than adding in here. This comes from unftoys_basic
//originally (oct/12/2006).
//last modified oct/12/2006
  void unftoys_new(int tau, int lastmombin, int toynum, int mb, int a, int mes, int ipcut, double gwidth, int nevents, int sqrtNerrs=1, int extradata){

  int bias=5; 
  char base[200];
  char mcbase[200];
  char VVFbase[200];
  char outbase[200];

  char smes[10], sipcut[5];

  if(mes==1)
    sprintf(smes,"sig");
  else if(mes==2)
    sprintf(smes,"sigpeak");
  else 
    sprintf(smes,"bq");

  sprintf(sipcut,"0%d",ipcut);

  sprintf(base,"_%d_%d_%s_%s",mb,a,smes,sipcut);
  sprintf(mcbase,"_%d_%d_%s_%s",466,131,smes,sipcut);
  if(ipcut)
    sprintf(VVFbase,"_%d_ip%s",mb,sipcut);
  else 
    sprintf(VVFbase,"_%d",mb);
  sprintf(outbase,"%s_%d_%d_%d_%d_%.1f_%d_%d_%d",base,tau,bias,lastmombin,toynum,gwidth,nevents,sqrtNerrs,extradata);

  cout << base << endl;

  cout << "mb         " << mb << endl;
  cout << "a          " << a << endl;
  cout << "mes        " << mes << endl;
  cout << "ipcut      " << ipcut << endl;
  cout << "tau        " << tau << endl;
  cout << "lastmombin " << lastmombin << endl;
  cout << "toynum     " << toynum << endl;
  cout << "gwidth     " << gwidth << endl;

  int numbins=18;
  double max =14.4;
  const int maxdepth=1;
  int mdepth=1;
  int oritau=tau;
  int gauss=0;
  int trutoys = 1000000;
  int mctoys=nevents*10;
  int dattoys=nevents;
  double mx, mxrec, hit, hit2;
  int offset=0;
  char fname[200];
  int cov1[5], cov2[5], cov3[5];

  for(int i=0; i<5; i++){
    cov1[i] = cov2[i] = cov3[i] = 0;
  }


  //Histos for the toys - MC and data
  TH1D *mcgen = new TH1D("mcgen_20","",numbins,0.,max);
  TH1D *mcreco = new TH1D("mcreco_11","",numbins,0.,max);
  TH2D *detmat_mc = new TH2D("detmat_mc_4000","",numbins,0.,max,numbins,0.,max);
  TH2D *probmat_mc = new TH2D("probmat_mc_4000","",numbins,0.,max,numbins,0.,max);

  TH1D *mcgen_ori = new TH1D("mcgen_20_ori","",numbins,0.,max);
  TH2D *detmat_mc_ori = new TH2D("detmat_mc_4000_ori","",numbins,0.,max,numbins,0.,max);

  TH1D *datgen = new TH1D("datgen","",numbins,0.,max);
  TH1D *datreco = new TH1D("datreco_30","",numbins,0.,max);
  TH1D *dattru = new TH1D("dattru_31","",numbins,0.,max);
  TH1D *daterr = new TH1D("daterr_130","",numbins,0.,max);
  TH1D *biasreco = new TH1D("biasreco","",numbins,0.,max);
  TH1D *dathighstat = new TH1D("dathighstat","",numbins,0.,max);
  TH1D *biastru = new TH1D("biastru","",numbins,0.,max);

  TH1D *chi2_andreas = new TH1D("chi2_andreas","",100,0.,100.);
  TH1D *chi2_andreas_large = new TH1D("chi2_andreas_large","",100,0.,2000.);
  TH1D *condition = new TH1D("condition","",100,0.,100);

  mcgen->Sumw2();
  mcgen_ori->Sumw2();
  datgen->Sumw2();
  dattru->Sumw2();
  datreco->Sumw2();

  TH1D *moment1 = new TH1D("moment1","",20,2.,4.5);
  TH1D *moment2 = new TH1D("moment2","",20,1.,3.5);
  TH1D *moment3 = new TH1D("moment3","",20,0.,3.);

  TH1D *mom1pull = new TH1D("mom1pull","",20,-5.,5.);
  TH1D *mom2pull = new TH1D("mom2pull","",20,-5.,5.);
  // TH1D *mom2pull = new TH1D("mom2pull","",40,-10.,10.);
  TH1D *mom3pull = new TH1D("mom3pull","",20,-5.,5.);
  //  TH1D *mom3pull = new TH1D("mom3pull","",80,-20.,20.);

  TH1D *mom1pull1 = new TH1D("mom1pull1","",20,-5.,5.);
  TH1D *mom2pull1 = new TH1D("mom2pull1","",20,-5.,5.);
  TH1D *mom3pull1 = new TH1D("mom3pull1","",20,-5.,5.);
  TH1D *mom1pull2 = new TH1D("mom1pull2","",20,-5.,5.);
  TH1D *mom2pull2 = new TH1D("mom2pull2","",20,-5.,5.);
  TH1D *mom3pull2 = new TH1D("mom3pull2","",20,-5.,5.);
  TH1D *mom1pull3 = new TH1D("mom1pull3","",20,-5.,5.);
  TH1D *mom2pull3 = new TH1D("mom2pull3","",20,-5.,5.);
  TH1D *mom3pull3 = new TH1D("mom3pull3","",20,-5.,5.);

  TH1D *mom1diff = new TH1D("mom1diff","",30,-1.5,1.5);
  TH1D *mom2diff = new TH1D("mom2diff","",30,-1.5,1.5);
  TH1D *mom3diff = new TH1D("mom3diff","",30,-2.5,2.5);

  TH1D *mom1diff1 = new TH1D("mom1diff1","",30,-1.5,1.5);
  TH1D *mom2diff1 = new TH1D("mom2diff1","",30,-1.5,1.5);
  TH1D *mom3diff1 = new TH1D("mom3diff1","",30,-2.5,2.5);
  TH1D *mom1diff2 = new TH1D("mom1diff2","",30,-1.5,1.5);
  TH1D *mom2diff2 = new TH1D("mom2diff2","",30,-1.5,1.5);
  TH1D *mom3diff2 = new TH1D("mom3diff2","",30,-2.5,2.5);
  TH1D *mom1diff3 = new TH1D("mom1diff3","",30,-1.5,1.5);
  TH1D *mom2diff3 = new TH1D("mom2diff3","",30,-1.5,1.5);
  TH1D *mom3diff3 = new TH1D("mom3diff3","",30,-2.5,2.5);

  TH1D *meanres = new TH1D("meanres","",numbins,0.,max);
  TH1D *meanerr = new TH1D("meanerr","",numbins,0.,max);
  TH2D *meancov = new TH2D("meancov","",numbins,0.,max,numbins,0.,max);
  TH2D *meancor = new TH2D("meancor","",numbins,0.,max,numbins,0.,max);


  TH1D *meanres1 = new TH1D("meanres1","",numbins,0.,max);
  TH1D *meanres2 = new TH1D("meanres2","",numbins,0.,max);
  TH1D *meanres3 = new TH1D("meanres3","",numbins,0.,max);
  TH1D *meanres4 = new TH1D("meanres4","",numbins,0.,max);
  TH1D *meanerr1 = new TH1D("meanerr1","",numbins,0.,max);
  TH1D *meanerr2 = new TH1D("meanerr2","",numbins,0.,max);
  TH1D *meanerr3 = new TH1D("meanerr3","",numbins,0.,max);
  TH1D *meanerr4 = new TH1D("meanerr4","",numbins,0.,max);


  TH1D *specs2 = new TH1D("specs2","",maxdepth,1.,(double)maxdepth+1.); 
  TH1D *mom1s2 = new TH1D("mom1s2","",maxdepth,1.,(double)maxdepth+1.); 
  TH1D *mom2s2 = new TH1D("mom2s2","",maxdepth,1.,(double)maxdepth+1.); 
  TH1D *mom3s2 = new TH1D("mom3s2","",maxdepth,1.,(double)maxdepth+1.); 
  TH1D *sums2 = new TH1D("sums2","",maxdepth,1.,(double)maxdepth+1.); 

  TLegend *leg = new TLegend(0.5,0.5,0.85,0.85);
  leg->SetBorderSize(0);
  leg->SetFillColor(10);

  //Get the distributions from files
  sprintf(fname,"toyfiles/toyinput%s.root",mcbase);
  cout << "Opening " << fname << " for MC input" << endl;
  TFile *fmchisto = new TFile(fname,"r");
  mchisto = (TH1D*)fmchisto->Get("mxtrue");
  mchisto->Sumw2();
  mchisto->Scale(1/mchisto->Integral());

  sprintf(fname,"toyfiles/toyinput%s.root",base);
  cout << "Opening " << fname << " for data truth input" << endl;
  TH1D *datatruehisto;
  TFile *fdatatruehisto = new TFile(fname,"r");
  datatruehisto = (TH1D*)fdatatruehisto->Get("mxtrue");
  datatruehisto->Sumw2();
  datatruehisto->Scale(1/datatruehisto->Integral());

  sprintf(fname,"toyfiles/mestoys%sFitres.root",VVFbase);
  cout << "Opening " << fname << " for data input" << endl;
  TH1D *datahisto;
  TFile *fdatahisto = new TFile(fname,"r");
  datahisto = (TH1D*)fdatahisto->Get("vubchop");
  datahisto->Sumw2();
  datahisto->Scale(1/datahisto->Integral());

  //Get the resolution histos and normalize to unit area
  TH1D *mcreso[13];
  char hname[20];
  Double_t mchistomax=0;
  sprintf(fname,"toyfiles/detresponse%s.root",mcbase);
  cout << "Opening " << fname << " for MC resolution input" << endl;
  TFile *mcfres = new TFile(fname,"r");
  for(int i=0; i<13; i++){
    sprintf(hname, "%s%d%s","bin",i+1,"histo");
    mcreso[i] = (TH1D*)mcfres->Get(hname);
    if(mchistomax<mcreso[i]->GetMaximum()){
      mchistomax=mcreso[i]->GetMaximum();
    }
  }

  cout << mchistomax << endl;

  TH1D *datareso[13];
  Double_t datahistomax=0;
  sprintf(fname,"toyfiles/detresponse%s.root",base);
  cout << "Opening " << fname << " for data resolution input" << endl;
  TFile *datafres = new TFile(fname,"r");
  for(int i=0; i<13; i++){
    sprintf(hname, "%s%d%s","bin",i+1,"histo");
    datareso[i] = (TH1D*)datafres->Get(hname);
    if(datahistomax<datareso[i]->GetMaximum()){
      datahistomax=datareso[i]->GetMaximum();
    }
  }

  cout << datahistomax << endl;


  //Other variables
  //Some matrices, vectors and histos for the unfolding
  TMatrixD datacov(numbins,numbins), udatacov(numbins,numbins), udetcov(numbins,numbins);
  TVectorD unfres(numbins), unfbias(numbins);
  TH1D* hunfres = new TH1D("result","result",numbins,0.,max);
  TH1D* hunfraw = new TH1D("raw result","raw result",numbins,0.,max);

  Double_t chi2, cchi2, tchi2, achi2; 
  Int_t mnum=numbins;
  Int_t trubin, missed;
  missed=0;
  TMatrixD tcov(mnum,mnum), tcov2(mnum,mnum), tmat(mnum,mnum) ;
  const int nb = numbins;
  double meanchi2[20], weights[nb];
  double werr, wtmp;
  for(int i=0; i<20; i++){
    meanchi2[i]=0.;
  }
  for(int i=0; i<numbins; i++){
    weights[i]=0.;
  }
  double specchi2, mom1chi2, mom2chi2, mom3chi2;
  const int maxtoynum=toynum;
  double specchi2vec[maxdepth][maxtoynum], mom1chi2vec[maxdepth][maxtoynum],  mom2chi2vec[maxdepth][maxtoynum], mom3chi2vec[maxdepth][maxtoynum];

  TUnfHisto *unfhist = new TUnfHisto(numbins,numbins);

  Double_t datamom1, datamom2, datamom3, dataerrmom1, dataerrmom2, dataerrmom3;
  Double_t datamomcor12, datamomcor13, datamomcor23;
  datamom1 = datamom2 = datamom3 = dataerrmom1 = dataerrmom2 = dataerrmom3 = datamomcor12 = datamomcor13 = datamomcor23 = 0;
  Double_t tmpmom1, tmpmom2, tmpmom3, tmperrmom1, tmperrmom2, tmperrmom3;
  Double_t tmpmomcor12, tmpmomcor13, tmpmomcor23;
  tmpmom1 = tmpmom2 = tmpmom3 = tmperrmom1 = tmperrmom2 = tmperrmom3 = tmpmomcor12 = tmpmomcor13 = tmpmomcor23 = 0;

  double dmom1[maxdepth], dmom2[maxdepth], dmom3[maxdepth], derrmom1[maxdepth], derrmom2[maxdepth], derrmom3[maxdepth];
  double dmomcor12[maxdepth], dmomcor13[maxdepth], dmomcor23[maxdepth];

  for(int i=0; i<mdepth; i++){
    dmom1[i] = dmom2[i] = dmom3[i] = derrmom1[i] = derrmom2[i] = derrmom3[i] = dmomcor12[i] = dmomcor13[i] = dmomcor23[i] = 0.;
  }

  TMatrixD datatucov(numbins+2, numbins+2);
  TMatrixD dataderiv(3,numbins+2), dataderivT(numbins+2,3), datatempmat(3,numbins+2), datamomcov(3,3);

  ISLQuantity datamom;
  TH1MomentsPlain *dataMoments0 = new TH1MomentsPlain(hunfres,1,lastmombin);

  //TRandom for the smearing of the weights
  TRandom grandom(4702475);
  //TRandom for the toys, use TRandom2 later for higher periodicity
  TRandom2 random(4796140);
  int trubin;
  int i=0;
  while(i<mctoys){
    mx=mchisto->GetRandom();
    mcgen->Fill(mx);
    mcgen_ori->Fill(mx);
    if(gwidth!=0.){
      mxrec = random.Gaus(mx, gwidth);
    }else{
      trubin=mchisto->FindBin(mx);
      if(trubin>13){
	missed++;
	continue;
      }
      mxrec=random.Uniform(0.,max);
      hit2=random.Uniform(0.,mchistomax);
      if(hit2 < mcreso[trubin-1]->GetBinContent(mcreso[trubin-1]->FindBin(mxrec))){
      } else {
	continue;
      }
    }
    detmat_mc->Fill(mxrec,mx);
    detmat_mc_ori->Fill(mxrec,mx);
    mcreco->Fill(mxrec);
    i++;
  }

  cout << "Entries in MC reco spectrum " << mcreco->Integral() << endl;

  for(i=0; i<trutoys; i++){
    mx=datatruehisto->GetRandom();
    dattru->Fill(mx);
  }

  cout << mcgen->Integral() << endl;
  cout << mcreco->Integral() << endl;
  cout << missed << endl;

  //The moments for the MC distribution
  //Now get the moments using Thorstens class
  Double_t tmom1, tmom2, tmom3;
  Double_t tmom1err, tmom2err, tmom3err;
  tmom1 = tmom2 = tmom3 = 0;

  //Fill our covariance matrix into a matrix like it is needed for Thorsten's
  //class, i.e. add zeroes.
  TMatrixD mctucov(numbins+2, numbins+2);
  
  for(int i=1; i<=numbins; i++){
    mctucov(i,i) = dattru->GetBinError(i)*dattru->GetBinError(i); 
  }
  
  //The derivative matrices for the covariance matrix on the moments and the
  //covariance matrix itself (see Thorsten's BAD 636, p.72)
  TMatrixD mcderiv(3,numbins+2), mcderivT(numbins+2,3), mctempmat(3,numbins+2), mcmomcov(3,3);

  //Needed to get the moments
  ISLQuantity tmom;
  ISLQuantity tmpmom;

  //Pass the histo and the first and last bin that should be used for the moments
  TH1MomentsPlain *tMoments0 = new TH1MomentsPlain(dattru,1,lastmombin);
  tMoments0->Initialise(dattru);

  tmom = tMoments0->M1();  //The first moment
  tmom1 = tmom.GetValue(); 
  tmom = tMoments0->VAR(); //The second central moment
  tmom2 = tmom.GetValue();
  tmom = tMoments0->V3();
  tmom3 = tmom.GetValue();

  //Compute the covariance matrix: momcov = deriv * tucov * derivT
  for(int i=0; i<numbins+2; i++){
    mcderiv(0,i) = tMoments0->M1Deriv(i);
    mcderiv(1,i) = tMoments0->VARDeriv(i);
    mcderiv(2,i) = tMoments0->V3Deriv(i);
    mcderivT(i,0) = tMoments0->M1Deriv(i);
    mcderivT(i,1) = tMoments0->VARDeriv(i); 
    mcderivT(i,2) = tMoments0->V3Deriv(i);
  }
  
  mctempmat.Mult(mcderiv,mctucov);
  mcmomcov.Mult(mctempmat,mcderivT);
  tmom1err = sqrt(mcmomcov(0,0));
  tmom2err = sqrt(mcmomcov(1,1));
  tmom3err = sqrt(mcmomcov(2,2));

  cout << "Mom1 " << tmom1 << "+-" << tmom1err << endl;
  cout << "Mom2 " << tmom2 << "+-" << tmom2err << endl;
  cout << "Mom3 " << tmom3 << "+-" << tmom3err << endl;
  
    for(int tnum=0; tnum<toynum; tnum++){
      cout << tnum << endl;
      //  Reset everything
      for(int i=0; i<numbins; i++){
	datgen->SetBinContent(i+1,0.);
	datgen->SetBinError(i+1,0.);
	datreco->SetBinContent(i+1,0.);
	datreco->SetBinError(i+1,0.);
	biasreco->SetBinContent(i+1,0.);
	biasreco->SetBinError(i+1,0.);
	hunfres->SetBinContent(i+1,0.);
	hunfres->SetBinError(i+1,0.);
	unfres(i)=unfbias(i)=unftheo(i)=0.;
	for(int j=0; j<numbins; j++){
	  datacov(i,j) = udatacov(i,j) = udetcov(i,j) = utheocov(i,j) = 0.;
	}
      }
      chi2=cchi2=tchi2=achi2=0.;
      tmpmom1 = tmpmom2 = tmpmom3 = tmperrmom1 = tmperrmom2 = tmperrmom3 = 0.; 
      tmpmomcor12 = tmpmomcor13 = tmpmomcor23 = 0;

      for(int i=0; i<mnum; i++)
	for(int j=0; j<mnum; j++)
	  tcov(i,j)=tcov2(i,j)=tmat(i,j)=0.;

      //Reset the truth histos
      for(int b=0; b<numbins; b++){
	mcgen->SetBinContent(b+1, mcgen_ori->GetBinContent(b+1));
	mcgen->SetBinError(b+1, mcgen_ori->GetBinError(b+1));
	for(int bb=0; bb<numbins; bb++){
	  detmat_mc->SetBinContent(b+1,bb+1, detmat_mc_ori->GetBinContent(b+1,bb+1));
	  detmat_mc->SetBinError(b+1,bb+1, detmat_mc_ori->GetBinError(b+1,bb+1));
	}
      }

      missed=0;
      i=0;
      while(i<dattoys){
	mxrec = datahisto->GetRandom();
	mx=datatruehisto->GetRandom();
	datgen->Fill(mx);
	if(!extradata){
	  if(gwidth!=0.){
	    mxrec = random.Gaus(mx, gwidth);
	  }else{
	    trubin=datatruehisto->FindBin(mx);
	    if(trubin>13){
	      missed++;
	      continue;
	    }
	    mxrec=random.Uniform(0.,max);
	    hit2=random.Uniform(0.,datahistomax);
	    if(hit2 < datareso[trubin-1]->GetBinContent(datareso[trubin-1]->FindBin(mxrec))){
	    } else {
	      continue;
	    }
	  }
	}
	datreco->Fill(mxrec);
	i++;
      }

      cout << datgen->Integral() << endl;
      cout << datreco->Integral() << endl;
      cout << missed << endl;

      for(int i=0; i<numbins; i++){
	if(sqrtNerrs){
	  datreco->SetBinError(i+1, sqrt(datreco->GetBinContent(i+1)));
	  daterr->SetBinContent(i+1, sqrt(datreco->GetBinContent(i+1)));
	}else{
	  datreco->SetBinError(i+1, datahisto->GetBinError(i+1));
	  daterr->SetBinContent(i+1, datahisto->GetBinError(i+1));
	}
	mcreco->SetBinError(i+1, sqrt(mcreco->GetBinContent(i+1)));
      }

      //      for(int i=4; i<14; i++){
      for(int i=0; i<numbins; i++){
	datacov(i,i) = datreco->GetBinContent(i+1);
      }

//       //A test plot
//       TCanvas *canv2 = new TCanvas("canv2","canv2",500,400);
//       datgen->Draw();
//       datreco->SetLineColor(2);
//       datreco->Draw("same");
 
      //These are the depth of iteration - do this before the init, because we need to feed in the reweighted distributions and detmat
      for(int depth=1; depth<=mdepth; depth++){

	if(depth==1){
	  tau=oritau;
	}
 
	//Now go to the unfolding
	unfhist->init(datreco, mcreco, mcgen, detmat_mc);
	
	datgen->Scale(1/datgen->Integral());
	mcgen->Scale(datgen->Integral()/mcgen->Integral());
	mcgen_ori->Scale(datgen->Integral()/mcgen_ori->Integral());
	dattru->Scale(1/dattru->Integral());
	
	unfres = unfhist->Unfold(tau);

	//Get the covariance matrix for statistical uncertainties on measured spectrum
	udatacov = unfhist->GetCov(datacov, datreco, 1000, tau);

	//Get the covariance matrix for statistical uncertainties on signal MC
	udetcov = unfhist->GetMatStatCov(1000,tau);

	for(int i=0; i<numbins; i++){
	  for(int j=0; j<numbins; j++){
	    //	    if(i==0 && j==0) cout << "stat " << udatacov(i,j) << endl;
	    udatacov(i,j) = udatacov(i,j) + udetcov(i,j);	
	    //	    if(i==0 && j==0) cout << "det " << udetcov(i,j) << endl;

	  }
	}


	for(Int_t i=0; i<numbins; i++){
	  hunfres->SetBinContent(i+1,unfres(i));
	  hunfres->SetBinError(i+1, sqrt(udatacov(i,i)));
	}
      
	for(int i=0; i<numbins; i++){
	  hunfraw->SetBinContent(i+1, hunfres->GetBinContent(i+1));
	  meanres->SetBinContent(i+1, meanres->GetBinContent(i+1) + hunfres->GetBinContent(i+1)/toynum/mdepth);
	  if(depth==1){
	    meanres1->SetBinContent(i+1, meanres1->GetBinContent(i+1) + hunfres->GetBinContent(i+1)/toynum);
	    meanerr1->SetBinContent(i+1, meanerr1->GetBinContent(i+1) + hunfres->GetBinError(i+1)/toynum);
	  }
	  if(depth==2){
	    meanres2->SetBinContent(i+1, meanres2->GetBinContent(i+1) + hunfres->GetBinContent(i+1)/toynum);
	    meanerr2->SetBinContent(i+1, meanerr2->GetBinContent(i+1) + hunfres->GetBinError(i+1)/toynum);
	  }
	  if(depth==3){
	    meanres3->SetBinContent(i+1, meanres3->GetBinContent(i+1) + hunfres->GetBinContent(i+1)/toynum);
	    meanerr3->SetBinContent(i+1, meanerr3->GetBinContent(i+1) + hunfres->GetBinError(i+1)/toynum);
	  }
	  if(depth==4){
	    meanres4->SetBinContent(i+1, meanres4->GetBinContent(i+1) + hunfres->GetBinContent(i+1)/toynum);
	    meanerr4->SetBinContent(i+1, meanerr4->GetBinContent(i+1) + hunfres->GetBinError(i+1)/toynum);
	  }
	  meanerr->SetBinContent(i+1, meanerr->GetBinContent(i+1) + hunfres->GetBinError(i+1)/toynum);
	  for(int j=0; j<numbins; j++){
	    meancov->SetBinContent(i+1,j+1, meancov->GetBinContent(i+1,j+1) + udatacov(i,j)/toynum);
	  }
	}
      

	//Get the chi2 taking into account bin-by-bin correlations
	for(int i=0; i<mnum; i++)
	  for(int j=0; j<mnum; j++){
	    tcov(i,j)=udatacov(i+offset,j+offset);
	    tcov2(i,j)=udatacov(i+offset,j+offset);
	    tmat(i,j)=udatacov(i+offset,j+offset);
	  }
	const Double_t det = tcov2.Determinant();
	cout << "det: " << det <<endl;
  
	//In case of a singular/near singular matrix SVD is the 
	//  right decomposition choice
	TDecompSVD svd(tcov2);
	const TVectorD singularVal = svd.GetSig();
      
	//We can set the the threshold for the singular values with SetTol
	//				       Those smaller than tol*max(singularVal) are set to 0 for
	//					 Inversion, Determinant, equation solving
	//						    svd.SetTol(1.0e-12);
	//      singularVal.Print();

	//	Let's calculate the matrix condition number, large numbers
	//number indicate  trouble. In case of SVD it is simply max(singular)/min(singular)
	cout << "condition: " <<  svd.Condition() << endl;
	if(svd.Condition()>0)
	  condition->Fill(log10(svd.Condition()));
	svd.Invert(tcov2);
	//    tcov2.Print();
	//    tcov.Print();
	//      (tcov2*tcov).Print();
	//    (tcov*tcov2).Print();
	//cout << "Now get the chi2" << endl;
	chi2 = 0.0;
	cchi2 = 0.;
	for(int l=0; l<mnum; l++){
	  for(int m=0; m<mnum; m++){
	    chi2 += (hunfres->GetBinContent(l+1) - dattru->GetBinContent(l+1)) * tcov2(l,m) * (hunfres->GetBinContent(m+1) - dattru->GetBinContent(m+1));
	  }
	}
// 	cout << "chi2 for raw spectrum " << chi2 << endl;

	for(int l=0; l<mnum-1; l++){
	  for(int m=0; m<mnum-1; m++){
	    cchi2 += (hunfres->GetBinContent(l+1) - dattru->GetBinContent(l+1)) * tcov2(l,m) * (hunfres->GetBinContent(m+1) - dattru->GetBinContent(m+1));
	  }
	}
	//cout << "check chi2 for raw spectrum " << cchi2 << endl;
  
	//Test chi2
	for(int m=0; m<numbins; m++){
	  if(udatacov(m,m)!=0.){
	    tchi2+= (hunfres->GetBinContent(m+1) - dattru->GetBinContent(m+1)) * (hunfres->GetBinContent(m+1) - dattru->GetBinContent(m+1))/udatacov(m,m);
	  }
	}

	//Andreas' chi2
	for(int m=0; m<numbins; m++){
	  if(dattru->GetBinError(m+1)){
	    achi2+= (hunfres->GetBinContent(m+1) - dattru->GetBinContent(m+1)) * (hunfres->GetBinContent(m+1) - dattru->GetBinContent(m+1)) / dattru->GetBinError(m+1) / dattru->GetBinError(m+1);
	  }
	  //      dattru->SetBinError(m+1, datgen->GetBinError(m+1));
	}


// 	cout << "Old chi2 " << tchi2 << endl;
// 	cout << "Old chi2 - corrected " << tbchi2 << endl;  
// 	cout << "Andreas chi2 " << achi2 << endl;
// 	cout << "Area of unfolded spectrum " << hunfres->Integral() << endl;
// 	cout << "Area of true spectrum " << datgen->Integral() << endl;
// 	cout << "Area of ini spectrum " << mcgen->Integral() << endl;
    
	chi2_andreas->Fill(achi2);
	chi2_andreas_large->Fill(achi2);

	//The moments for the data distribution

	//Fill our covariance matrix into a matrix like it is needed for Thorsten's
	//class, i.e. add zeroes.
	for(int i=1; i<=numbins; i++){
	  for(int j=1; j<=numbins; j++){
	    datatucov(i,j) = udatacov(i-1,j-1); 
	  }
	}
    
	//Pass the histo and the first and last bin that should be used for the moments
	dataMoments0->Initialise(hunfres);
    
	tmpmom = dataMoments0->M1();  //The first moment
	tmpmom1 = tmpmom.GetValue(); 
	tmpmom = dataMoments0->VAR(); //The second central moment
	tmpmom2 = tmpmom.GetValue();
	tmpmom = dataMoments0->V3();
	tmpmom3 = tmpmom.GetValue();
    
	//Compute the covariance matrix: momcov = deriv * tucov * derivT
	for(int i=0; i<numbins+2; i++){
	  dataderiv(0,i) = dataMoments0->M1Deriv(i);
	  dataderiv(1,i) = dataMoments0->VARDeriv(i);
	  dataderiv(2,i) = dataMoments0->V3Deriv(i);
	  dataderivT(i,0) = dataMoments0->M1Deriv(i);
	  dataderivT(i,1) = dataMoments0->VARDeriv(i); 
	  dataderivT(i,2) = dataMoments0->V3Deriv(i);
	}
    
	datatempmat.Mult(dataderiv,datatucov);
	datamomcov.Mult(datatempmat,dataderivT);
	tmperrmom1 = sqrt(datamomcov(0,0));
	tmperrmom2 = sqrt(datamomcov(1,1));
	tmperrmom3 = sqrt(datamomcov(2,2));
	tmpmomcor12 = datamomcov(1,0)/tmperrmom1/tmperrmom2;
	tmpmomcor13 = datamomcov(2,0)/tmperrmom1/tmperrmom3;
	tmpmomcor23 = datamomcov(1,2)/tmperrmom2/tmperrmom3;

	//Compute the effective chi2 to decide about iteration
	specchi2=0.;
	//The spectrum part:
	for(int m=0; m<numbins; m++){
	  if(dattru->GetBinError(m+1)){
	    specchi2+= (hunfres->GetBinContent(m+1) - dattru->GetBinContent(m+1)) * (hunfres->GetBinContent(m+1) - dattru->GetBinContent(m+1)) / dattru->GetBinError(m+1) / dattru->GetBinError(m+1);
	  }
	}
	meanchi2[depth-1] += specchi2/toynum;
	specchi2vec[depth-1][tnum] = specchi2;

	//The moments part:
	mom1chi2 = (tmpmom1-tmom1)*(tmpmom1-tmom1)/tmom1err/tmom1err;
	mom2chi2 = (tmpmom2-tmom2)*(tmpmom2-tmom2)/tmom2err/tmom2err;
	mom3chi2 = (tmpmom3-tmom3)*(tmpmom3-tmom3)/tmom3err/tmom3err;

	mom1chi2vec[depth-1][tnum] = mom1chi2;
	mom2chi2vec[depth-1][tnum] = mom2chi2;
	mom3chi2vec[depth-1][tnum] = mom3chi2;

	meanchi2[depth-1]+=(tmpmom1-tmom1)*(tmpmom1-tmom1)/tmom1err/tmom1err/toynum;
	meanchi2[depth-1]+=(tmpmom2-tmom2)*(tmpmom2-tmom2)/tmom2err/tmom2err/toynum;
	meanchi2[depth-1]+=(tmpmom3-tmom3)*(tmpmom3-tmom3)/tmom3err/tmom3err/toynum;
       
	//	cout << specs2->GetBinContent(depth) << endl;
	specs2->SetBinContent(depth, specs2->GetBinContent(depth) + specchi2/toynum);
	//	cout << specs2->GetBinContent(depth) << endl;
	mom1s2->SetBinContent(depth, mom1s2->GetBinContent(depth) + mom1chi2/toynum);
	mom2s2->SetBinContent(depth, mom2s2->GetBinContent(depth) + mom2chi2/toynum);
	mom3s2->SetBinContent(depth, mom3s2->GetBinContent(depth) + mom3chi2/toynum);
	sums2->SetBinContent(depth, meanchi2[depth-1]);

	cout << endl;
	cout << "The mean chi2 - spectrum and moments combined - " << depth << ". iteration - " << specchi2+mom1chi2+mom2chi2+mom3chi2 << endl;

	//Reweight according to the least result
	for(int iw=0; iw<numbins; iw++){
	  if(mcgen->GetBinContent(iw+1)>0. && hunfres->GetBinContent(iw+1)>=0.){
	    //	    cout << iw+1 << "  " << hunfres->GetBinContent(iw+1) << "  " << mcgen->GetBinContent(iw+1) << endl;
	    weights[iw] = (double)hunfres->GetBinContent(iw+1)/mcgen->GetBinContent(iw+1);
	    werr = (double)hunfres->GetBinError(iw+1)/mcgen->GetBinContent(iw+1);
	    //	    cout << iw+1 << "  " << weights[iw] << "+-" << werr << endl;
	    //	    wtmp = grandom.Gaus(weights[iw],werr);
	    //	    weights[iw] = wtmp;
	    //	    cout << iw+1 << "  " << weights[iw] << endl;
	  }else{
	    weights[iw]=1.;
	  }
	}
	for(int iw=0; iw<numbins; iw++){
	  mcgen->SetBinContent(iw+1, mcgen->GetBinContent(iw+1)*weights[iw]);
// 	  if(depth==2){
//	  cout << mcgen->GetBinContent(iw+1) << endl;
// 	  }
	}
	for(int iw=0; iw<numbins; iw++){
	  for(int jw=0; jw<numbins; jw++){
	    //	    cout << iw+1 << "  " << jw+1 << "  " << detmat_mc->GetBinContent(iw+1, jw+1) << endl;
	    detmat_mc->SetBinContent(iw+1,jw+1, detmat_mc->GetBinContent(iw+1,jw+1)*weights[jw]);
	    //	    cout << iw+1 << "  " << jw+1 << "  " << detmat_mc->GetBinContent(iw+1, jw+1) << endl;
// 	    if(depth==2){
// 	      cout << iw+1 << "  " << jw+1 << "  " << detmat_mc->GetBinContent(iw+1,jw+1) << endl;
// 	    }
	  }
	}

	dmom1[depth-1] = dmom1[depth-1] + tmpmom1;  
	dmom2[depth-1] = dmom2[depth-1] + tmpmom2;  
	dmom3[depth-1] = dmom3[depth-1] + tmpmom3; 
	derrmom1[depth-1] = derrmom1[depth-1] + tmperrmom1;
	derrmom2[depth-1] = derrmom2[depth-1] + tmperrmom2;
	derrmom3[depth-1] = derrmom3[depth-1] + tmperrmom3;
	dmomcor12[depth-1] = dmomcor12[depth-1] + tmpmomcor12;
	dmomcor13[depth-1] = dmomcor13[depth-1] + tmpmomcor13;
	dmomcor23[depth-1] = dmomcor23[depth-1] + tmpmomcor23;

	if(depth==1){
	  mom1pull1->Fill((tmom1-tmpmom1)/tmperrmom1);
	  mom2pull1->Fill((tmom2-tmpmom2)/tmperrmom2);
	  mom3pull1->Fill((tmom3-tmpmom3)/tmperrmom3);
	  mom1diff1->Fill((tmom1-tmpmom1));
	  mom2diff1->Fill((tmom2-tmpmom2));
	  mom3diff1->Fill((tmom3-tmpmom3));
	}
	if(depth==2){
	  mom1pull2->Fill((tmom1-tmpmom1)/tmperrmom1);
	  mom2pull2->Fill((tmom2-tmpmom2)/tmperrmom2);
	  mom3pull2->Fill((tmom3-tmpmom3)/tmperrmom3);
	  mom1diff2->Fill((tmom1-tmpmom1));
	  mom2diff2->Fill((tmom2-tmpmom2));
	  mom3diff2->Fill((tmom3-tmpmom3));
	}
	if(depth==3){
	  mom1pull3->Fill((tmom1-tmpmom1)/tmperrmom1);
	  mom2pull3->Fill((tmom2-tmpmom2)/tmperrmom2);
	  mom3pull3->Fill((tmom3-tmpmom3)/tmperrmom3);
	  mom1diff3->Fill((tmom1-tmpmom1));
	  mom2diff3->Fill((tmom2-tmpmom2));
	  mom3diff3->Fill((tmom3-tmpmom3));
	}

      }

	datamom1 = datamom1 + tmpmom1;  
	datamom2 = datamom2 + tmpmom2;  
	datamom3 = datamom3 + tmpmom3; 
	dataerrmom1 = dataerrmom1 + tmperrmom1;
	dataerrmom2 = dataerrmom2 + tmperrmom2;
	dataerrmom3 = dataerrmom3 + tmperrmom3;
	datamomcor12 = datamomcor12 + tmpmomcor12;
	datamomcor13 = datamomcor13 + tmpmomcor13;
	datamomcor23 = datamomcor23 + tmpmomcor23;
   
	moment1->Fill(tmpmom1);
	moment2->Fill(tmpmom2);
	moment3->Fill(tmpmom3);
	mom1pull->Fill((tmom1-tmpmom1)/tmperrmom1);
	mom2pull->Fill((tmom2-tmpmom2)/tmperrmom2);
	mom3pull->Fill((tmom3-tmpmom3)/tmperrmom3);
	mom1diff->Fill((tmom1-tmpmom1));
	mom2diff->Fill((tmom2-tmpmom2));
	mom3diff->Fill((tmom3-tmpmom3));

// 	cout << "for the diff " << tmom1 << "  " << tmpmom1 << "  " << tmperrmom1 << endl;
// 	cout << (tmom1-tmpmom1) << endl;
// 	cout << mom1chi2 << endl;
        
	if(TMath::Abs(int((tmom1-tmpmom1)/tmperrmom1))<5)
	  cov1[TMath::Abs(int((tmom1-tmpmom1)/tmperrmom1))]++;
	if(TMath::Abs(int((tmom2-tmpmom2)/tmperrmom2))<5)
	  cov2[TMath::Abs(int((tmom2-tmpmom2)/tmperrmom2))]++;
	if(TMath::Abs(int((tmom3-tmpmom3)/tmperrmom3))<5)
	  cov3[TMath::Abs(int((tmom3-tmpmom3)/tmperrmom3))]++;

    }


    double meanspec[maxdepth], meanmom1[maxdepth], meanmom2[maxdepth], meanmom3[maxdepth], meansum[maxdepth];
    double varspec[maxdepth], varmom1[maxdepth], varmom2[maxdepth], varmom3[maxdepth], varsum[maxdepth];
    for(int d=0; d<maxdepth; d++){
      meanspec[d] = meanmom1[d] = meanmom2[d] = meanmom3[d] = meansum[d] = 0.;
      varspec[d] = varmom1[d] = varmom2[d] = varmom3[d] = varsum[d] = 0.;
    }
    for(int d=0; d<maxdepth; d++){
      for(int t=0; t<toynum; t++){
	meanspec[d] += specchi2vec[d][t]/toynum;
	meanmom1[d] += mom1chi2vec[d][t]/toynum;
	//	cout << d << " " << t << "  " << mom1chi2vec[d][t] << endl;
	meanmom2[d] += mom2chi2vec[d][t]/toynum;
	meanmom3[d] += mom3chi2vec[d][t]/toynum;
      }
      meansum[d] = meanspec[d] + meanmom1[d] + meanmom2[d] + meanmom3[d];
    }
    for(int d=0; d<maxdepth; d++){
      for(int t=0; t<toynum; t++){
	varspec[d] += (specchi2vec[d][t] - meanspec[d])*(specchi2vec[d][t] - meanspec[d])/(toynum-1);
	varmom1[d] += (mom1chi2vec[d][t] - meanmom1[d])*(mom1chi2vec[d][t] - meanmom1[d])/(toynum-1);
	varmom2[d] += (mom2chi2vec[d][t] - meanmom2[d])*(mom2chi2vec[d][t] - meanmom2[d])/(toynum-1);
	varmom3[d] += (mom3chi2vec[d][t] - meanmom3[d])*(mom3chi2vec[d][t] - meanmom3[d])/(toynum-1);
      }
      varsum[d] = varspec[d] + varmom1[d] + varmom2[d] + varmom3[d];
    }

    for(int d=0; d<maxdepth; d++){
      specs2->SetBinError(d+1, sqrt(varspec[d]));
      mom1s2->SetBinError(d+1, sqrt(varmom1[d]));
      mom2s2->SetBinError(d+1, sqrt(varmom2[d]));
      mom3s2->SetBinError(d+1, sqrt(varmom3[d]));
      sums2->SetBinError(d+1,  sqrt(varsum[d]));
    }

    char outname[200];
    sprintf(outname, "toyres/qualstoys_res%s.root",outbase);
    TFile *fout = new TFile(outname, "RECREATE");
  
    datreco->SetLineColor(52);
    datreco->SetLineWidth(2);
    datreco->SetLineStyle(3);
    datreco->SetXTitle("mX^2 / GeV^2");
    leg->AddEntry(dattru,"True","l");
    leg->AddEntry(datreco, "Smeared", "l");
    leg->AddEntry(hunfres, "Unfolded", "p");
    datreco->Scale(1/datreco->Integral());
    datreco->Write();
    dattru->Write();
    hunfres->Write();
    leg->Write();

    chi2_andreas->Write();
    chi2_andreas_large->Write();

    condition->Write();
    moment1->Write();
    moment2->Write();
    moment3->Write();

    mom1pull->Write();
    mom2pull->Write();
    mom3pull->Write();
	
    mom1diff->Write();
    mom2diff->Write();
    mom3diff->Write();

    mom1pull1->Write();
    mom2pull1->Write();
    mom3pull1->Write();
    mom1diff1->Write();
    mom2diff1->Write();
    mom3diff1->Write();
    mom1pull2->Write();
    mom2pull2->Write();
    mom3pull2->Write();
    mom1diff2->Write();
    mom2diff2->Write();
    mom3diff2->Write();
    mom1pull3->Write();
    mom2pull3->Write();
    mom3pull3->Write();
    mom1diff3->Write();
    mom2diff3->Write();
    mom3diff3->Write();

    specs2->Write();
    mom1s2->Write();
    mom2s2->Write();
    mom3s2->Write();
    sums2->Write();
	  
    for(int i=0; i<numbins; i++){
      meanres->SetBinError(i+1, meanerr->GetBinContent(i+1));
      meanres1->SetBinError(i+1, meanerr1->GetBinContent(i+1));
      meanres2->SetBinError(i+1, meanerr2->GetBinContent(i+1));
      meanres3->SetBinError(i+1, meanerr3->GetBinContent(i+1));
      meanres4->SetBinError(i+1, meanerr4->GetBinContent(i+1));
      for(int j=0; j<numbins; j++){
	if(meancov->GetBinContent(i+1,i+1)*meancov->GetBinContent(j+1,j+1)){
	  meancor->SetBinContent(i+1,j+1, meancov->GetBinContent(i+1,j+1)/sqrt(meancov->GetBinContent(i+1,i+1))/sqrt(meancov->GetBinContent(j+1,j+1)));
	}
      }
    } 

    mcgen_ori->SetLineColor(52);
    mcgen_ori->SetTitle("");
    mcgen_ori->SetXTitle("mX^2 / GeV^2");
    mcgen_ori->SetLineStyle(3);
    dattru->SetLineColor(8);
    mcgen_ori->SetLineWidth(2);
    dattru->SetLineWidth(2);
    dattru->SetTitle("");
    dattru->SetXTitle("mX^2 / GeV^2");
    meanres->SetLineColor(94);
    meanres->SetLineWidth(2);
    meanres->SetMarkerStyle(8);
    meanres->SetMarkerColor(94);
    meanres->SetMarkerSize(0.8);

    meanres1->SetLineColor(94);
    meanres1->SetLineWidth(2);
    meanres1->SetMarkerStyle(8);
    meanres1->SetMarkerColor(94);
    meanres1->SetMarkerSize(0.8);
    meanres2->SetLineColor(94);
    meanres2->SetLineWidth(2);
    meanres2->SetMarkerStyle(8);
    meanres2->SetMarkerColor(94);
    meanres2->SetMarkerSize(0.8);
    meanres3->SetLineColor(94);
    meanres3->SetLineWidth(2);
    meanres3->SetMarkerStyle(8);
    meanres3->SetMarkerColor(94);
    meanres3->SetMarkerSize(0.8);
    meanres4->SetLineColor(94);
    meanres4->SetLineWidth(2);
    meanres4->SetMarkerStyle(8);
    meanres4->SetMarkerColor(94);
    meanres4->SetMarkerSize(0.8);

    mcgen_ori->Write();
    meanres->Write();
    meanres1->Write();
    meanres2->Write();
    meanres3->Write();
    meanres4->Write();
    meancov->Write();
    meancor->Write();

    fout->Close();

    TCanvas *canv1 = new TCanvas("canv1","canv1",1000,1000);
    canv1->Divide(2,2);

    canv1->cd(1);
    meanres->Draw();
    mcgen_ori->Draw("samehist");
    dattru->Draw("samehist");
    meanres->Draw("same");

    canv1->cd(2);
    datgen->Draw("hist");
    datreco->Draw("same");

    canv1->cd(3);
    detmat_mc_ori->Draw();

    specs2->SetXTitle("Depth of iteration");
    specs2->SetYTitle("S^{2}_{spectrum}");
    specs2->SetTitleOffset(1.1,"Y");
    sums2->SetXTitle("Depth of iteration");
    sums2->SetYTitle("S^{2}_{comb}");
    sums2->SetTitleOffset(1.1,"Y");
    mom1s2->SetXTitle("Depth of iteration");
    mom1s2->SetYTitle("S^{2}_{M1}");
    mom1s2->SetTitleOffset(1.1,"Y");
    mom2s2->SetXTitle("Depth of iteration");
    mom2s2->SetYTitle("S^{2}_{M'2}");
    mom2s2->SetTitleOffset(1.1,"Y");
    mom3s2->SetXTitle("Depth of iteration");
    mom3s2->SetYTitle("S^{2}_{M'3}");
    mom3s2->SetTitleOffset(1.1,"Y");

    TCanvas *canv2 = new TCanvas("canv2","canv2",1000,1000);
    canv2->Divide(2,2);
    canv2->cd(1);
    specs2->Draw();
    canv2->cd(2);
    mom1s2->Draw();
    canv2->cd(3);
    mom2s2->Draw();
    canv2->cd(4); 
    mom3s2->Draw();

    TCanvas *canv3 = new TCanvas("canv3","canv3",500,400);
    sums2->Draw();

    //The moments for the MC distribution
    //Now get the moments using Thorstens class
    Double_t mcmom1, mcmom2, mcmom3, mcerrmom1, mcerrmom2, mcerrmom3;
    Double_t mcmomcor12, mcmomcor13, mcmomcor23;
    mcmom1 = mcmom2 = mcmom3 = mcerrmom1 = mcerrmom2 = mcerrmom3 = mcmomcor12 = mcmomcor13 = mcmomcor23 = 0;

    for(int i=1; i<=numbins; i++){
      mctucov(i,i) = mcgen_ori->GetBinError(i)*mcgen_ori->GetBinError(i); 
      //      cout << i << "  " << mcgen->GetBinError(i) << "  " << mcgen_ori->GetBinError(i) << "  " << dattru->GetBinError(i) << "  " << hunfres->GetBinError(i) << endl;
    }


    //Needed to get the moments
    ISLQuantity mcmom;

    cout << "Moments for mcgen" << endl;

    //Pass the histo and the first and last bin that should be used for the moments
    TH1MomentsPlain *mcMoments0 = new TH1MomentsPlain(mcgen_ori,1,lastmombin);
    mcMoments0->Initialise(mcgen_ori);

    mcmom = mcMoments0->M1();  //The first moment
    mcmom1 = mcmom.GetValue();
    mcmom = mcMoments0->VAR(); //The second central moment
    mcmom2 = mcmom.GetValue();
    mcmom = mcMoments0->V3();
    mcmom3 = mcmom.GetValue();

    //Compute the covariance matrix: momcov = deriv * tucov * derivT
    for(int i=0; i<numbins+2; i++){
      mcderiv(0,i) = mcMoments0->M1Deriv(i);
      mcderiv(1,i) = mcMoments0->VARDeriv(i);
      mcderiv(2,i) = mcMoments0->V3Deriv(i);
      mcderivT(i,0) = mcMoments0->M1Deriv(i);
      mcderivT(i,1) = mcMoments0->VARDeriv(i); 
      mcderivT(i,2) = mcMoments0->V3Deriv(i);
    }

    mctempmat.Mult(mcderiv,mctucov);
    mcmomcov.Mult(mctempmat,mcderivT);
    mcerrmom1 = sqrt(mcmomcov(0,0));
    mcerrmom2 = sqrt(mcmomcov(1,1));
    mcerrmom3 = sqrt(mcmomcov(2,2));
    mcmomcor12 = mcmomcov(1,0)/mcerrmom1/mcerrmom2;
    mcmomcor13 = mcmomcov(2,0)/mcerrmom1/mcerrmom3;
    mcmomcor23 = mcmomcov(1,2)/mcerrmom2/mcerrmom3;

    cout << "1. Moment       " << mcmom1 << " +- " << mcerrmom1 << endl;
    cout << "2. Moment       " << mcmom2 << " +- " << mcerrmom2 << endl;
    cout << "3. Moment       " << mcmom3 << " +- " << mcerrmom3 << endl;
    cout << "Correlation 1 2 " << mcmomcor12 << endl;
    cout << "Correlation 1 3 " << mcmomcor13 << endl;
    cout << "Correlation 2 3 " << mcmomcor23 << endl;

    printf("&%.4f&%.4f&%.4f&%.4f&%.4f&%.4f\n&%.3f&%.3f&%.3f\n",mcmom1,mcerrmom1,mcmom2,mcerrmom2,mcmom3,mcerrmom3,mcmomcor12,mcmomcor13,mcmomcor23);
    printf("MOM %.6f %.6f %.6f %.1f %.1f %.1f\n", mcmom1, mcmom2, mcmom3, 0., 0., 0.);
    printf("ERRMOM %.6f %.6f %.6f %.1f %.1f %.1f\n", mcerrmom1, mcerrmom2, mcerrmom3, 0., 0., 0.);
    printf("CORRMX %.6f %.6f %.6f\n", mcmomcor12, mcmomcor13, mcmomcor23);

    //The true data moments
    mcmom1 = mcmom2 = mcmom3 = mcerrmom1 = mcerrmom2 = mcerrmom3 = mcmomcor12 = mcmomcor13 = mcmomcor23 = 0;
    for(int i=1; i<=numbins; i++){
      mctucov(i,i) = dattru->GetBinError(i)*dattru->GetBinError(i); 
    }

    cout << endl;
    cout << "Moments for dattru" << endl;
    mcMoments0->Initialise(dattru);

    mcmom = mcMoments0->M1();  //The first moment
    mcmom1 = mcmom.GetValue(); 
    mcmom = mcMoments0->VAR(); //The second central moment
    mcmom2 = mcmom.GetValue();
    mcmom = mcMoments0->V3();
    mcmom3 = mcmom.GetValue();

    //Compute the covariance matrix: momcov = deriv * tucov * derivT
    for(int i=0; i<numbins+2; i++){
      mcderiv(0,i) = mcMoments0->M1Deriv(i);
      mcderiv(1,i) = mcMoments0->VARDeriv(i);
      mcderiv(2,i) = mcMoments0->V3Deriv(i);
      mcderivT(i,0) = mcMoments0->M1Deriv(i);
      mcderivT(i,1) = mcMoments0->VARDeriv(i); 
      mcderivT(i,2) = mcMoments0->V3Deriv(i);
    }

    mctempmat.Mult(mcderiv,mctucov);
    mcmomcov.Mult(mctempmat,mcderivT);
    mcerrmom1 = sqrt(mcmomcov(0,0));
    mcerrmom2 = sqrt(mcmomcov(1,1));
    mcerrmom3 = sqrt(mcmomcov(2,2));
    mcmomcor12 = mcmomcov(1,0)/mcerrmom1/mcerrmom2;
    mcmomcor13 = mcmomcov(2,0)/mcerrmom1/mcerrmom3;
    mcmomcor23 = mcmomcov(1,2)/mcerrmom2/mcerrmom3;

    cout << "1. Moment       " << mcmom1 << " +- " << mcerrmom1 << endl;
    cout << "2. Moment       " << mcmom2 << " +- " << mcerrmom2 << endl;
    cout << "3. Moment       " << mcmom3 << " +- " << mcerrmom3 << endl;
    cout << "Correlation 1 2 " << mcmomcor12 << endl;
    cout << "Correlation 1 3 " << mcmomcor13 << endl;
    cout << "Correlation 2 3 " << mcmomcor23 << endl;

    printf("&%.4f&%.4f&%.4f&%.4f&%.4f&%.4f\n",mcmom1,mcerrmom1,mcmom2,mcerrmom2,mcmom3,mcerrmom3);
    printf("&%.3f&%.3f&%.3f\n",mcmomcor12,mcmomcor13,mcmomcor23);
    printf("MOM %.6f %.6f %.6f %.1f %.1f %.1f\n", mcmom1, mcmom2, mcmom3, 0., 0., 0.);
    printf("ERRMOM %.6f %.6f %.6f %.1f %.1f %.1f\n", mcerrmom1, mcerrmom2, mcerrmom3, 0., 0., 0.);
    printf("CORRMX %.6f %.6f %.6f\n", mcmomcor12, mcmomcor13, mcmomcor23);


    //The data moments averages
    if(toynum){
      datamom1 = datamom1/toynum;
      datamom2 = datamom2/toynum;
      datamom3 = datamom3/toynum;
      dataerrmom1 = dataerrmom1/toynum;
      dataerrmom2 = dataerrmom2/toynum;
      dataerrmom3 = dataerrmom3/toynum;
      datamomcor12 = datamomcor12/toynum;
      datamomcor13 = datamomcor13/toynum;
      datamomcor23 = datamomcor23/toynum;
    }

    cout << endl;
    cout << "The data moments" << endl;
    cout << "1. Moment       " << datamom1 << " +- " << dataerrmom1 << endl;
    cout << "2. Moment       " << datamom2 << " +- " << dataerrmom2 << endl;
    cout << "3. Moment       " << datamom3 << " +- " << dataerrmom3 << endl;
    cout << "Correlation 1 2 " << datamomcor12 << endl;
    cout << "Correlation 1 3 " << datamomcor13 << endl;
    cout << "Correlation 2 3 " << datamomcor23 << endl;

    printf("&%.4f&%.4f&%.4f&%.4f&%.4f&%.4f\n",datamom1,dataerrmom1,datamom2,dataerrmom2,datamom3,dataerrmom3);
    printf("&%.3f&%.3f&%.3f\n",datamomcor12,datamomcor13,datamomcor23);
    printf("MOM %.6f %.6f %.6f %.1f %.1f %.1f\n", datamom1, datamom2, datamom3, 0., 0., 0.);
    printf("ERRMOM %.6f %.6f %.6f %.1f %.1f %.1f\n", dataerrmom1, dataerrmom2, dataerrmom3, 0., 0., 0.);
    printf("CORRMX %.6f %.6f %.6f\n", datamomcor12, datamomcor13, datamomcor23);

    printf("%.4f +- %.4f <br>\n", datamom1, dataerrmom1);
    printf("(%.4f, %.4f) <br>\n", sqrt( (datamom1-mcmom1)*(datamom1-mcmom1) + dataerrmom1*dataerrmom1 ), mcmom1-datamom1 );
    printf("%.4f +- %.4f <br>\n", datamom2, dataerrmom2);
    printf("(%.4f, %.4f) <br>\n", sqrt( (datamom2-mcmom2)*(datamom2-mcmom2) + dataerrmom2*dataerrmom2 ), mcmom2-datamom2 );
    printf("%.4f +- %.4f <br>\n", datamom3, dataerrmom3);
    printf("(%.4f, %.4f) <br>\n", sqrt( (datamom3-mcmom3)*(datamom3-mcmom3) + dataerrmom3*dataerrmom3 ), mcmom3-datamom3 );
    printf("%.4f <br>\n", datamomcor12);
    printf("%.4f <br>\n", datamomcor13);
    printf("%.4f <br>\n", datamomcor23);
    printf("%.0f <br>\n", specs2->GetBinContent(1));
    printf("%.0f <br>\n", mom1s2->GetBinContent(1));
    printf("%.0f <br>\n", mom2s2->GetBinContent(1));
    printf("%.0f <br>\n", mom3s2->GetBinContent(1));

    cout << endl;
    cout << "Coverage tests" << endl;
    cout << cov1[0] << "  " << cov1[1] << "  " << cov1[2] << "  " << cov1[3] << "  " << cov1[4] << endl;
    cout << cov2[0] << "  " << cov2[1] << "  " << cov2[2] << "  " << cov2[3] << "  " << cov2[4] << endl;
    cout << cov3[0] << "  " << cov3[1] << "  " << cov3[2] << "  " << cov3[3] << "  " << cov3[4] << endl;

    printf("&%d&%d&%d&%d&%d\n",cov1[0],cov1[1],cov1[2],cov1[3],cov1[4]);
    printf("&%d&%d&%d&%d&%d\n",cov2[0],cov2[1],cov2[2],cov2[3],cov2[4]);
    printf("&%d&%d&%d&%d&%d\n",cov3[0],cov3[1],cov3[2],cov3[3],cov3[4]);

    //For the different depth of iteration
    for(int i=0; i<mdepth; i++){
      if(toynum){
	dmom1[i] = dmom1[i]/toynum;
	dmom2[i] = dmom2[i]/toynum;
	dmom3[i] = dmom3[i]/toynum;
	derrmom1[i] = derrmom1[i]/toynum;
	derrmom2[i] = derrmom2[i]/toynum;
	derrmom3[i] = derrmom3[i]/toynum;
	dmomcor12[i] = dmomcor12[i]/toynum;
	dmomcor13[i] = dmomcor13[i]/toynum;
	dmomcor23[i] = dmomcor23[i]/toynum;
      }

      cout << endl;
      cout << "The data moments" << endl;
      cout << "1. Moment       " << dmom1[i] << " +- " << derrmom1[i] << endl;
      cout << "2. Moment       " << dmom2[i] << " +- " << derrmom2[i] << endl;
      cout << "3. Moment       " << dmom3[i] << " +- " << derrmom3[i] << endl;
      cout << "Correlation 1 2 " << dmomcor12[i] << endl;
      cout << "Correlation 1 3 " << dmomcor13[i] << endl;
      cout << "Correlation 2 3 " << dmomcor23[i] << endl;
      
      printf("&%.4f&%.4f&%.4f&%.4f&%.4f&%.4f\n",dmom1[i],derrmom1[i],dmom2[i],derrmom2[i],dmom3[i],derrmom3[i]);
      printf("&%.3f&%.3f&%.3f\n",dmomcor12[i],dmomcor13[i],dmomcor23[i]);
      printf("MOM %.6f %.6f %.6f %.1f %.1f %.1f\n", dmom1[i], dmom2[i], dmom3[i], 0., 0., 0.);
      printf("ERRMOM %.6f %.6f %.6f %.1f %.1f %.1f\n", derrmom1[i], derrmom2[i], derrmom3[i], 0., 0., 0.);
      printf("CORRMX %.6f %.6f %.6f\n", dmomcor12[i], dmomcor13[i], dmomcor23[i]);
    }
 
}