#include <TH1D.h>
#include <TH2D.h>
#include <TPad.h>
#include <TDirectory.h>
#include <TEntryList.h>
#include <TFile.h>
#include <TMath.h>
#include <TTree.h>
#include <TChain.h>
#include <THStack.h>
#include <TCanvas.h>
#include <TF1.h>
#include <TLegend.h>
#include <stdio.h>
#include <iostream>
#include <deque>
#include <list>
#include <queue>
#include <fstream>

#include "EtaVEL.h"
#include "EtaVEL.cpp"
/*
Zum erstellen der correction map ist createGainAndEtaFile(...) in EVELAlg.C der entry point.
Zum erstellen des HR images ist createImage(...) der entry point.
*/
int etabins = 25;
int nEtas = 25;
Double_t dum[3][3];
Int_t x,y,f,q;

int counter[5];
int remoteCounter[5];

//TH2D *sum = new TH2D("sum","sum",3,-0.1,2.1,3,-0.1,2.1);
//TH2F *subPos = new TH2F("subPos","subPos", 100, -1.,1. ,100, -1.,1.);
TH2D *subPosAEta = new TH2D("subPosAEta","subPosAEta", 50, -.5,1.5 ,50, -.5,1.5);
TH2D *subPosBEta = new TH2D("subPosBEta","subPosBEta", 50, -.5,1.5 ,50, -.5,1.5);



TH1D *cE = new TH1D("clusterEnergy","clusterEnergy",400, 0.,4000.);
//TH1D *cES = new TH1D("clusterEnergyS","clusterEnergyS",400, 0.,4000.);


TH2D *cES3vs2 = new TH2D("clusterEnergy3vs2","clusterEnergy3vs2",800, 0.,8000.,600,0.,6000.);
TH2D *cES3vs2S = new TH2D("clusterEnergy3vs2S","clusterEnergy3vs2S",800, 0.,8000.,600,0.,6000.);

double th = 0.99;
double sigmas = 1.0;

TH2D *imgRLR = new TH2D("imgRLR","imgRLR",160,0.0,160.0 ,160 ,0.0,160.0);
TH2D *imgLR = new TH2D("imgLR","imgLR",160*2,0.0,160.0 ,160*2 ,0.0,160.0);

TH2D *clusHist= new TH2D("clusHist","clusHist",3,-0.5,2.5,3,-0.5,2.5);
TH2D *clusHistC= new TH2D("clusHistC","clusHistC",3,-0.5,2.5,3,-0.5,2.5);

int **imgArray;

int findShape(Double_t cluster[3][3], double sDum[2][2]){
  int corner = -1;

  double sum = 	cluster[0][0] + cluster[1][0] + cluster[2][0] + cluster[0][1] + cluster[1][1] + cluster[2][1] + cluster[0][2] + cluster[1][2] + cluster[2][2];
	
  double sumTL = cluster[0][0] + cluster[1][0] + cluster[0][1] + cluster[1][1]; //2 ->BL
  double sumTR = cluster[1][0] + cluster[2][0] + cluster[2][1] + cluster[1][1]; //0 ->TL
  double sumBL = cluster[0][1] + cluster[0][2] + cluster[1][2] + cluster[1][1]; //3 ->BR
  double sumBR = cluster[1][2] + cluster[2][1] + cluster[2][2] + cluster[1][1]; //1 ->TR
  double sumMax = 0;


  //double **sDum = subCluster;
  Double_t ssDum[2][2];

  // if(sumTL  >= sumMax){
  sDum[0][0] = cluster[0][0]; sDum[1][0] = cluster[1][0];
    sDum[0][1] = cluster[0][1]; sDum[1][1] = cluster[1][1];

    ssDum[0][0] = cluster[0][0]; ssDum[1][0] = cluster[0][1];
    ssDum[0][1] = cluster[1][0]; ssDum[1][1] = cluster[1][1];

    corner = 2;
    sumMax=sumTL;
    // } 

  if(sumTR  >= sumMax){
    sDum[0][0] = cluster[1][0]; sDum[1][0] = cluster[2][0];
    sDum[0][1] = cluster[1][1]; sDum[1][1] = cluster[2][1];

    ssDum[0][0] = cluster[2][0]; ssDum[1][0] = cluster[2][1];
    ssDum[0][1] = cluster[1][0]; ssDum[1][1] = cluster[1][1];

    corner = 0;
    sumMax=sumTR;
  } 

  if(sumBL  >= sumMax){
    sDum[0][0] = cluster[0][1]; sDum[1][0] = cluster[1][1];
    sDum[0][1] = cluster[0][2]; sDum[1][1] = cluster[1][2];

    ssDum[0][0] = cluster[0][2]; ssDum[1][0] = cluster[0][1];
    ssDum[0][1] = cluster[1][2]; ssDum[1][1] = cluster[1][1];

    corner = 3;
    sumMax=sumBL;
  }

  if(sumBR  >= sumMax){
    sDum[0][0] = cluster[1][1]; sDum[1][0] = cluster[2][1];
    sDum[0][1] = cluster[1][2]; sDum[1][1] = cluster[2][2];

    ssDum[0][0] = cluster[2][2]; ssDum[1][0] = cluster[2][1];
    ssDum[0][1] = cluster[1][2]; ssDum[1][1] = cluster[1][1];

    corner = 1;
    sumMax=sumBR;
  }

  switch(corner){
  case 0:
    cES3vs2->Fill(sum,sumTR); break;
  case 1:
    cES3vs2->Fill(sum,sumBR); break;
  case 2:
    cES3vs2->Fill(sum,sumTL); break;
  case 3:
    cES3vs2->Fill(sum,sumBL); break;
  }
  
  counter[corner]++;
  remoteCounter[q]++;

  //  cout << "local corner is: " << corner << " remote corner is: " << q << endl;

  return corner;
}




int placePhoton( TH2D *img, double subCluster[2][2], int cX, int cY, int corner, double *sX, double *sY, double *scX, double *scY){
  double tot = subCluster[0][0] + subCluster[0][1] + subCluster[1][0] + subCluster[1][1];
  double t = subCluster[1][0] + subCluster[1][1];
  double r = subCluster[0][1] + subCluster[1][1];
  
  double xHitC = r/tot;
  double yHitC = t/tot;

  imgRLR->Fill(cX,cY);

  cE->Fill(tot);

  double dX, dY;

  //before looking at annas code
  /*  if(corner == 0){ dX=-1.; dY=-1.;   }
  if(corner == 1){ dX=-1.; dY=+1.;  }
  if(corner == 2){ dX=+1.; dY=-1.;  }
  if(corner == 3){ dX=+1.; dY=+1.;    }*/

  if(corner == 0){ dX=-1.; dY=+1.;   } //top left
  if(corner == 1){ dX=+1.; dY=+1.;  } //top right
  if(corner == 2){ dX=-1.; dY=-1.;  } //bottom left
  if(corner == 3){ dX=+1.; dY=-1.;    } //bottom right

  imgLR->Fill(cX+0.25*dX,cY+0.25*dY);

  double posX = ((double)cX) + 0.5*dX + xHitC;
  double posY = ((double)cY) + 0.5*dY + yHitC;

  subPosBEta->Fill(xHitC ,yHitC);
  if(img){
    img->Fill(posX,posY);
  }

  if(xHitC < 0.02&& yHitC < 0.02){

    cES3vs2S->Fill(dum[0][0]+dum[0][1]+dum[0][2]+dum[1][0]+dum[1][1]+dum[1][2]+dum[2][0]+dum[2][1]+dum[2][2],subCluster[0][0]+subCluster[0][1]+subCluster[1][0]+subCluster[1][1]);
  }
  

  if(sX && sY && scX && scY){
    *sX = xHitC; //0.5 + 0.5*dX + xHitC;
    *sY = yHitC; //0.5 + 0.5*dY + yHitC;
    *scX = ((double)cX) + 0.5*dX;
    *scY = ((double)cY) + 0.5*dY;
  }
  return 1;
}



void placePhotonCorr(TH2D *img, EtaVEL *e,double sX, double sY, double scX, double scY){
  int bin = e->findBin(sX,sY);
  if(bin <= 0) return;
  double subX = ((double)(e->getXBin(bin))+.5)/((double)e->getNPixels());
  double subY = ((double)(e->getYBin(bin))+.5)/((double)e->getNPixels());

  if(img!=NULL){
    img->Fill(scX+ subX , scY+ subY);
  }
  subPosAEta->Fill(subX,subY);

  int iscx = scX;
  int iscy = scY;
  if(iscx >=nx || iscx<0 || iscy >=ny || iscy<0) return;
  //cout << iscx*e->getNPixels()+e->getXBin(bin) << " " << iscy*e->getNPixels()+e->getXBin(bin) << endl;
  if(img==NULL) return;
  imgArray[iscx*e->getNPixels()+e->getXBin(bin)][iscy*e->getNPixels()+e->getYBin(bin)]++;
}

void gainCorrection(Double_t corrected[3][3], TH2D *gainMap){

  for(int xx = 0; xx < 3; xx++)
    for(int yy = 0; yy < 3; yy++){
      if(gainMap && gainMap->GetBinContent(x+xx+2,y+yy+2) != 0){
	corrected[xx][yy] = dum[xx][yy] / gainMap->GetBinContent(x+xx+2,y+yy+2);
	clusHistC->Fill(xx,yy,corrected[xx][yy]);
      }
      else
	corrected[xx][yy] = dum[xx][yy];

      clusHist->Fill(xx,yy,dum[xx][yy]);
    }
}


EtaVEL *plotEtaDensity(TChain* tree2, TEntryList *el, EtaVEL *oldEta = NULL, TH2D **img = NULL, TH2D *gainMap=NULL, int nPixels=25) {


  
  EtaVEL *newEta = new EtaVEL(25,-0.02,1.02);

  Long64_t listEntries=el->GetN();
  Long64_t treeEntry;
  Long64_t chainEntry;
  
  Int_t treenum=0;
  tree2->SetEntryList(el);   

  double gainCorrC[3][3];
  double subCluster[2][2];
  double sX, sY, scX, scY;

  cout << "Events: " << listEntries << endl;
  if(oldEta == NULL){  cout << "Old Eta is NULL " <<  endl; }
  for(int i = 0; i<4; i++){ counter[i] = 0; remoteCounter[i] = 0; }

   for (Long64_t il =0; il<listEntries;il++) {
     treeEntry = el->GetEntryAndTree(il,treenum);
      
      chainEntry = treeEntry+tree2->GetTreeOffset()[treenum];
      if (tree2->GetEntry(chainEntry)) {

	gainCorrection(gainCorrC,gainMap);
	//cout << gainCorrC[1][1] << endl;

	//finds corner
	int corner = findShape(gainCorrC,subCluster);

	int validEvent;
	
	
	if(img){
	  validEvent = placePhoton(img[0],subCluster,x,y, corner, &sX, &sY, &scX, &scY);
	}else{
	  //calc etaX, etaY
	  validEvent = placePhoton(NULL,subCluster,x,y, corner, &sX, &sY, &scX, &scY);
	}

	//fill etavel
	newEta->fill(sX,sY);
	



	if(oldEta && img && img[1]){
	  placePhotonCorr(img[1],oldEta, sX,sY, scX, scY);
	}else{
	  placePhotonCorr(NULL,newEta,sX,sY,scX,scY);
	}
	
	
      }
      //cout << il << endl;
      int ssize = 500000;
      if(il % ssize == 0 && il != 0 && oldEta==NULL){

	cout << " -------------- "<< endl;
	newEta->updatePixelPos();
	

	//newEta->resolveSelfIntersect();
	char tit[1000];
	/*	TFile *ff = new TFile("/scratch/Spider.root","UPDATE");
	sprintf(tit,"subPosAEta%i",newEta->getIt()); subPosAEta->SetName(tit);
	subPosAEta->Write(); subPosAEta->Reset();
	sprintf(tit,"subPosBEta%i",newEta->getIt()); subPosBEta->SetName(tit);
	subPosBEta->Write(); subPosBEta->Reset();
	sprintf(tit,"Eta%i",newEta->getIt()); newEta->Write(tit);
	ff->Close(); */
	//il = 0;
      }
      
      if(il % ssize == ssize-1){
	double prog = (double)il/(double)listEntries*100.;
	cout << prog << "%" << endl;
	//if(prog > 19.) return newEta;
	if(newEta->converged == 1){ cout << "converged ... " << endl; return newEta; }
      }
      
   }

   cout << "local corners: " ;
   for(int i = 0; i<4; i++) cout << i << ": " << counter[i] << " || " ;
   cout << endl;

   //cout << "remote corners: " ;
   //for(int i = 0; i<4; i++) cout << i << ": " << remoteCounter[i] << " || " ;
   //cout << endl;
   
   return newEta;
}




TChain *openTree(char *tname, char *fname,double lEc, double hEc, double rms=5., char *chainName=">>thischan"){
  TChain *tree2;
  //  TH1D **etaDI;
  char cut[1000];

  tree2=new TChain(tname);
  tree2->Add(fname);
  tree2->Print();

  //sprintf(cut,"(x<=40) && (data[%d][%d]>%f*rms) && Sum$(data)<%f && Sum$(data)>%f",1,1,rms, hEc, lEc);
  //  sprintf(cut,"(x<=40) && (data[%d][%d]>%f*rms)",1,1,rms);// && Sum$(data)<%f && Sum$(data)>%f",1,1,rms, hEc, lEc);
  sprintf(cut,"(x<=40) && Sum$(data)<%f && Sum$(data)>%f", hEc, lEc);
  //    sprintf(cut,"");
  cout << cut << endl;

  tree2->Draw(chainName, cut, "entrylist");


  tree2->SetBranchAddress("iFrame",&f);
  tree2->SetBranchAddress("x",&x);
  tree2->SetBranchAddress("y",&y);
  tree2->SetBranchAddress("data",dum);
  //tree2->SetBranchAddress("q",&q);
  
  cout << "openTree : end" << endl;
  return tree2;
}

EtaVEL *etaDensity(char *tname, char *fname, double lEc = 1000, double hEc=3000, TH2D *gainMap=NULL, int nPixels=25) {
  /** open tree and make selection */
  TChain *tree2 = openTree(tname,fname,lEc,hEc);
  TEntryList *elist = (TEntryList*)gDirectory->Get("thischan");
  if(elist == NULL) { cout << "could not open tree " << endl; return NULL; }
  
  EtaVEL *etaDen = plotEtaDensity(tree2,elist,NULL,NULL,gainMap,nPixels);

  
  //etaDen->Draw("colz");
  cout << "done" << endl;

  return etaDen;
}

void interpolate(char *tname, char *fname, EtaVEL *etaDI, double lEc = 1000, double hEc=3000, TH2D *gainMap=NULL) {

  TChain *tree2 = openTree(tname,fname,lEc,hEc,5.,">>intChain");
  TEntryList *elist = (TEntryList*)gDirectory->Get("intChain");
  if(elist == NULL) { cout << "could not open tree " << endl; return; }

  double nPixels = (double)etaDI->getNPixels();

  TH2D **img = new TH2D*[3];
  img[0] = new TH2D("img","img",nPixels*160,0.0,160.0 ,nPixels*160 ,0.0,160.0);
  img[1] = new TH2D("imgE","imgE",nPixels*160,0.0,160.0 ,nPixels*160 ,0.0,160.0);

  int inPixels = etaDI->getNPixels();

  imgArray = new int*[inPixels*160];
  for(int i = 0; i < inPixels*160; i++){
    imgArray[i] = new int[inPixels*160];
    for(int j = 0; j < inPixels*160; j++){
      imgArray[i][j] = 0;
    }
  }

  cout << "starting" << endl;
  plotEtaDensity(tree2,elist, etaDI,img,gainMap);

  //img->Draw("colz");
}


TH2D *createGainMap(char *tname, char *fname, double lEc = 0,double hEc=10000){
  char name[100];
  TH1D *avgSpec3 = new TH1D("avgSpec3", "avgSpec3",hEc/20,0,hEc); 
  TH1D ***specs3 = new TH1D**[160];
  TH1D ***specs1 = new TH1D**[160];
  for(int xx = 0; xx < 160; xx++){
    specs3[xx] = new TH1D*[160];
    specs1[xx] = new TH1D*[160];
    for(int yy = 0; yy < 160; yy++){
      sprintf(name,"S3x%iy%i",xx,yy);
      specs3[xx][yy] = new TH1D(name,name,hEc/20,0,hEc);
      sprintf(name,"S1x%iy%i",xx,yy);
      specs1[xx][yy] = new TH1D(name,name,hEc/20,0,hEc);
    }
  }
  

  TChain *tree2 = openTree(tname,fname,0,hEc,5.,">>gainChan");
  TEntryList *elist = (TEntryList*)gDirectory->Get("gainChan");
  if(elist == NULL) { cout << "could not open tree " << endl; return NULL; }

  Long64_t listEntries=elist->GetN();
  Long64_t treeEntry;
  Long64_t chainEntry;
  
  Int_t treenum=0;
  tree2->SetEntryList(elist);   

  cout << "Events: " << listEntries << endl;
  for(int i = 0; i<4; i++) counter[i] = 0;
   for (Long64_t il =0; il<listEntries;il++) {
     treeEntry = elist->GetEntryAndTree(il,treenum);
     chainEntry = treeEntry+tree2->GetTreeOffset()[treenum];
     
     if (tree2->GetEntry(chainEntry)) {
       double sum = 0;
       for(int xx = 0; xx < 3; xx++)
	 for(int yy = 0; yy < 3; yy++)
	   sum += dum[xx][yy];
       specs3[x][y]->Fill(sum);
       specs1[x][y]->Fill(dum[1][1]);
       avgSpec3->Fill(sum);
     }
   }  

   TH2D *gainMap3 = new TH2D("gainMap3","gainMap3",160,-0.5,160.-0.5,160,-.5,160.-.5);
   TH2D *gainMap1 = new TH2D("gainMap1","gainMap1",160,-0.5,160.-0.5,160,-.5,160.-.5);
   for(int xx = 0; xx < 160; xx++){
     for(int yy = 0; yy < 160; yy++){
        TF1 *gf3 = new TF1("gf3","gaus", lEc, hEc);
	specs3[xx][yy]->Fit(gf3,"Q");
	double e3 = gf3->GetParameter(1);
	gainMap3->Fill(xx,yy,e3);

        TF1 *gf1 = new TF1("gf1","gaus", lEc, hEc);
	specs1[xx][yy]->Fit(gf1,"Q");
	double e1 = gf1->GetParameter(1);
	gainMap1->Fill(xx,yy,e1);

     }
   }

   return gainMap3;
}

void writeMatlab2DHisto(int xx, int yy,char *outFileName){
  ofstream outFile;
  outFile.open (outFileName);

  cout << "create matlab file with " << xx << " xbins and " << yy << " ybins" << endl;
  
  for(int y = 0; y < yy; y++){
      for(int x = 0; x < xx; x++){
	outFile << imgArray[x][y] << "\t";
      }
      outFile << endl;
  }
  
  outFile.close();
}

//COMPLETE STUFF

void createImage(char *tdir, char *tname, char *ftname, char *ifname = NULL, int useGM=0, double lEth=-1., double hEth=-1.){
  imgRLR->Reset();
  imgLR->Reset();

  char fname[1000];
  char inFName[1000];
  char outFName[1000];
  char moutFName[1000];
  if(ifname == NULL){
    sprintf(fname,"%s/%s_*.root",tdir,tname);
  }else{
    sprintf(fname,"%s",ifname);
  }

  if(useGM)  sprintf(inFName,"%s/%s-PlotsWGMVEL.root",tdir,ftname);
  else  sprintf(inFName,"%s/%s-PlotsVEL.root",tdir,ftname);

  sprintf(outFName,"%s/%s-ImgVEL.root",tdir,tname);
  sprintf(moutFName,"%s/%s-ImgVEL.mf",tdir,tname);

  TFile *inFile = new TFile(inFName,"READ");
  
  cout << "Image Tree File Name: " << fname << endl;
  cout << "Eta File Name: " << inFName << endl;
  cout << "Out File Name: " << outFName << endl;
  cout << "Matlab Out File Name: " << moutFName << endl;

  TH2D *gm = NULL;
  if(useGM){
    cout << "Load gain map" << endl;
    gm = (TH2D *)gDirectory->Get("gainMap");
    if(gm == NULL){ cout << "can not find gainMap in file" << endl; return; }
  }

  cout << "Load eta" << endl;
  EtaVEL *ee = (EtaVEL *)gDirectory->Get("etaDist");

  cout << "Select Energy BW" << endl;
  TH1D *spec = (TH1D *)gDirectory->Get("avgSpec3");
  if(spec == NULL){ cout << "can not find avgSpec3" << endl; return; }

  TF1 *gf3 = new TF1("gf3","gaus", 0, 10000);
  spec->Fit(gf3,"Q");
  double avgE = gf3->GetParameter(1);
  double sigE = gf3->GetParameter(2);
  cout << "avgE: " << avgE << " sigE: " << sigE << endl;
  cout << endl;

  if(lEth == -1.) lEth = avgE-5.*sigE;
  if(hEth == -1.) hEth = avgE+5.*sigE;   
  cout << lEth << " < E < " << hEth << " (eV)" << endl;

  cout << "start with interpolation" << endl;
  interpolate( tname, fname, ee,lEth,hEth ,gm);  

    
  TH2D *img = (TH2D *)gDirectory->Get("img");
  if(img == NULL){ cout << "could not find 2d-histogram: img " << endl; return; }


  TH2D *imgE = (TH2D *)gDirectory->Get("imgE");
  if(imgE == NULL){ cout << "could not find 2d-histogram: imgE " << endl; return; }


  //TH2D *imgEOM = (TH2D *)gDirectory->Get("imgEOM");
  //if(imgEOM == NULL){ cout << "could not find 2d-histogram: imgEOM " << endl; return; }

  TFile *outFile = new TFile(outFName,"UPDATE");
  imgLR->Write();
  imgRLR->Write();
  imgE->Write();
  //imgEOM->Write();
  img->Write();
  outFile->Close();
  inFile->Close(); 
  cout << "writing matlab file: " << moutFName << endl;
  writeMatlab2DHisto(160*ee->getNPixels(),160*ee->getNPixels(),moutFName);
  cout << "Done : " << outFName << endl;
  
}

/**
   \par tdir input tree directory
   \par tname input tree name
   \par ifname input file name if different than tdir/tname_*.root
   \par useGM use gain map
   \par maxExpEinEv spectrum maximum
   \par nPixels sub-pixels bins
   \par lEth low threshold
   \par hEth high threshold

 */


EtaVEL *createGainAndEtaFile(char *tdir, char *tname, char *ifname=NULL, int useGM=0, double maxExpEinEv=25000., int nPixels =25, double lEth=-1., double hEth=-1.){
  char fname[1000];
  char outFName[1000];


  if(ifname == NULL){
    sprintf(fname,"%s/%s_*.root",tdir,tname);
  }else{
    sprintf(fname,"%s",ifname);
  }

  if(useGM)  sprintf(outFName,"%s/%s-PlotsWGVEL.root",tdir,tname);
  else  sprintf(outFName,"%s/%s-PlotsVEL.root",tdir,tname);

  
  cout << "Tree File Name: " << fname << endl;
  cout << "Output File Name: " << outFName << endl;

  /** creates gain map and 3x3 spectrum */
  cout << "Creating gain map: " << endl;
  TH2D *gm = createGainMap(tname,fname,0,maxExpEinEv/10.);
  gm->SetName("gainMap");


  /** gets average 3x3 spectrum and fits it with a gaus */
  TH1D *spec = (TH1D *)gDirectory->Get("avgSpec3");
  if(spec == NULL){ cout << "can not find avgSpec3" << endl; return NULL; }
  TF1 *gf3 = new TF1("gf3","gaus", 0, maxExpEinEv/10.);
  spec->Fit(gf3,"Q");
  double avgE = gf3->GetParameter(1);
  double sigE = gf3->GetParameter(2);  
  cout << "avgE: " << avgE << " sigE: " << sigE << endl;
  cout << endl;


  /** sets high and low threshold if not given by the user */
  if(lEth == -1.) lEth = avgE-5.*sigE;
  if(hEth == -1.) hEth = avgE+5.*sigE;
  cout << lEth << " < E < " << hEth << " (eV)" << endl;




  cout << "calculating eta stuff" << endl;

  EtaVEL *newEta;
  if(useGM) newEta = etaDensity(tname,fname,lEth,hEth,gm,nPixels);
  else newEta = etaDensity(tname,fname,lEth,hEth,NULL,nPixels);

  cout << "writing to file " << outFName << endl;

  TFile *outFile = new TFile(outFName,"UPDATE");

  newEta->Write("etaDist");
  
  gm->Write();
  spec->Write();
  subPosAEta->Write();
  cES3vs2->Write();

  outFile->Close();
  cout << "Done : " << outFName << endl;
  return newEta;
}

void exportSpec(char *tdir, char *tname){
  char tfname[1000];
  char ofname[1000];
  char cleanName[1000];

  for(int p = 0; p < strlen(tname);p++){
    cleanName[p+1] = '\0';
    cleanName[p] = tname[p];

    if(tname[p] == '-') cleanName[p] = '_';
  }
  
  sprintf(tfname,"%s/%s-PlotsVEL.root",tdir,tname);
  sprintf(ofname,"%s/%s_SpecVEL.m",tdir,cleanName);
  TFile *tf = new TFile(tfname);
  TH1D *spec = (TH1D *)gDirectory->Get("avgSpec3");
  
  ofstream outFile;
  outFile.open (ofname);

  if(outFile.fail()){
    cout << "Could not open file : " << ofname << endl;
    return;
  }

  cout << "create matlab file with with spec " << ofname << endl;
  

  outFile << cleanName << " = [ " << endl;
  for(int i = 0; i < spec->GetNbinsX(); i++){ 
    outFile << i << " " << spec->GetBinCenter(i) << " " << spec->GetBinContent(i) << " ; " << endl; 
  }
  
  outFile << " ] ; " << endl;

  outFile.close();
}