#include <stdlib.h>
 
#include <stdio.h>
 
#include "TH1F.h"
 
#include "TH2F.h"
 
#include "TH1.h"
 
#include "TStyle.h"
 
#include "TString.h"
 
#include "TGraph2D.h"
 
#include "TGraph.h"
 
#include "TCanvas.h"
 
#include "TROOT.h"
 
#include "TFile.h"
 
#include "TCanvas.h"
 
 
 
int qeplot(char *fname="2010_05_19_mcp9000596_d.dat",int nx=9, int ny=9,int wl=4000, const char *snc="9000596" ){
 
/*
 
hdr[0]=1;//recid
 
hdr[1]=(nref*5+5)*sizeof(int);
 
hdr[2]=j;
 
hdr[3]=xpos;
 
hdr[4]=ypos;
 
 
 
hdr[5]=xpos;//recid
 
hdr[6]=ypos;
 
fhdr[7]=bgrtokref;
 
fhdr[8]=bgrtok;
 
hdr[9]=time(NULL);
 
 
 
 
 
sdata[0]=wl[0][i];
 
sdata[1]=tok[0][i];
 
sdata[2]=tok[j][i];
 
sdata[3]=qe[0][i];
 
sdata[4]=qe[j][i];
 
   */
 
gROOT->Reset();
 
char rname[256];
 
FILE 
*fp
=fopen(fname
,"r"); 
if (!fp) return -1;
 
TFile *f=new TFile(rname,"RECREATE");
 
int hdr[5];
 
int hdrp[5];
 
float sdata[5];
 
float *fhdr=(float *) hdrp;
 
TGraph2D * gr2d=0;
 
TGraph   * noised=0;
 
TH2F * h2d=0;
 
TH1F * h1d=0;
 
TGraph *noise= new TGraph();
 
 
 
const float fac=0.3595e-3;
 
 
 
TString sn(snc);
 
int nwl=46;
 
float minwl=200;
 
float maxwl=660;
 
TString swl;
 
swl.Form(" %d nm",int(wl/10.));
 
 
 
if (nx==1) {
 
   noised=new TGraph();
 
   h1d=new TH1F("h1d",(sn+swl+TString(";y [steps]")).Data(),ny,-0.5,ny-0.5);
 
   h2d=new TH2F("h2d",(sn+swl+TString(";y [steps];wavelength [nm]")).Data(),ny,-0.5,ny-0.5,nwl,minwl-1,maxwl-1);
 
}
 
if (ny==1) {
 
   noised=new TGraph();
 
   h1d=new TH1F("h1d",(sn+swl+TString(";y [steps]")).Data(),nx,-0.5,nx-0.5);
 
   h2d=new TH2F("h2d",(sn+swl+TString(";y [steps];wavelength [nm]")).Data(),nx,-0.5,nx-0.5,nwl,minwl-1,maxwl-1);
 
}
 
if (ny>1 && nx>1) {
 
  gr2d=new TGraph2D(ny*nx);
 
 
 
  h2d=new TH2F("h2d",(sn+swl+TString(";x [steps];y [steps]")).Data(),nx,-0.5,nx-0.5,ny,-0.5,ny-0.5);
 
} 
 
int posk0=0;
 
 
 
int icount=0;
 
  int nb
=fread(hdr
,sizeof(int),5,fp
);  
//  printf("[%d] %d %d %d %d %d\n",nb, hdr[0], hdr[1], hdr[2], hdr[3],hdr[4]);
 
  if (nb!=5) break;
 
  switch (hdr[0]){
 
     case 3:{
 
      int nb
=fread(hdrp
,sizeof(int),5,fp
);  
      printf("x=%u y=%u  | i0=%g i1=%g t=%u\t",hdrp
[0],hdrp
[1],fhdr
[2],fhdr
[3],hdrp
[4]);   
      hdr[3]=hdrp[0];
 
      hdr[4]=hdrp[1];
 
      noise->SetPoint(icount++,hdrp[4],-fhdr[3]);
 
     }
 
     case 2:
 
     case 1:{
 
        int nref=hdr[1]/sizeof(int)/5-1;
 
        if (hdr[0]==3) nref--;
 
        printf("%d [%d] x=%d y=%d\n",hdr
[2],nref
, hdr
[3],hdr
[4]);   
        int posk=hdr[2]-1;
 
        for (int i=0;i<nref;i++) {
 
           nb
=fread(sdata
,sizeof(float),5,fp
); 
           
 
           if (ny>1 && nx>1) { 
 
             if (sdata[0]==wl){
 
               // printf("%d x=%d y=%d z=%f\n",posk, hdr[3],hdr[4],sdata[4]); 
 
                if (hdr[3]==93933) hdr[3]=225001;
 
                if (hdr[3]==154024) hdr[3]=285000;
 
                if (hdr[3]==15350)  hdr[3]=277501;
 
                if (hdr[3]==282289) hdr[3]=270001;
 
                
 
                h2d->Fill(float(posk/ny),float(posk%ny),sdata[4]);
 
                gr2d->SetPoint(posk0,hdr[3]*fac,hdr[4]*fac,sdata[4]);
 
                posk0++;
 
                
 
             }
 
           } else  {
 
             if (sdata[0]==wl) {
 
               if (hdr[4]==81679)  hdr[4]= 126000;
 
               h1d->Fill(hdr[2],sdata[4]);
 
               if (nx==1) noised->SetPoint(posk0,hdr[4]*fac,sdata[4]);
 
               if (ny==1) noised->SetPoint(posk0,hdr[3]*fac,sdata[4]);
 
               posk0++;
 
             }
 
             h2d->Fill(hdr[2],sdata[0]/10.,sdata[4]);
 
             
 
             //printf("%5.0f %g %g %g %g\n",sdata[0],sdata[1],sdata[2],sdata[3],sdata[4]);     
 
           }
 
           
 
           
 
        }
 
        break;
 
}
 
  }
 
}
 
 
 
 
 
 gStyle->SetStatColor(0);
 
  //gStyle->SetTitleColor(0);
 
  gStyle->SetCanvasColor(0);
 
  gStyle->SetPadColor(0);
 
  gStyle->SetPadBorderMode(0);
 
  gStyle->SetCanvasBorderMode(0);
 
  gStyle->SetFrameBorderMode(0);
 
  
 
gStyle->SetPalette(1);
 
gStyle->SetOptStat(0);
 
 
 
 
 
if (h2d) h2d->Draw("colz");
 
if (h1d) h1d->Draw();
 
if (noise) {
 
TCanvas *c0 = new TCanvas("c0","dark noise vs time");
 
 
 
noise->SetNameTitle("noise",sn+TString(";time(s);dark current(A)"));
 
noise->Draw("AWL");
 
noise->GetXaxis()->SetTimeDisplay(kTRUE); 
 
 
 
noise->SetLineWidth(1);
 
noise->SetLineColor(kBlue);
 
noise->SetMarkerStyle(21);
 
noise->GetXaxis()->SetTimeDisplay(1);  // The X axis is a time axis
 
noise->GetYaxis()->SetLabelSize(0.02);
 
          
 
noise->GetXaxis()->SetLabelSize(0.02);
 
noise->GetXaxis()->SetTimeFormat("%d.%m.%H:%M");
 
noise->GetXaxis()->SetTitle("Time");
 
 
 
noise->Draw("AWL");
 
 TString str;
 
 str.Form("_darknoise");
 
 if (nx==1) str+=TString("_y"); 
 
 if (ny==1) str+=TString("_x");
 
 if (ny>1 && nx>1) str+=TString("_xy"); 
 
 c0->SaveAs((sn+str+TString(".png")).Data());
 
 c0->SaveAs((sn+str+TString(".pdf")).Data());
 
 
 
}
 
if (noised && ( nx==1 || ny==1)) {
 
 TCanvas *c = new TCanvas("c","QE1d");
 
 noised->Draw("AWL");
 
 
 
 
 
 TString str;
 
 str.Form("_%d_nm",int(wl/10.));
 
 
 
 
 
 if (nx==1) {
 
   str+=TString("_y"); 
 
   noised->SetNameTitle("noised",sn+swl+TString(";y [mm]"));
 
 } else {
 
   str+=TString("_x");
 
   noised->SetNameTitle("noised",sn+swl+TString(";x [mm]"));
 
 } 
 
 c->SaveAs((sn+str+TString(".pdf")).Data());
 
 c->SaveAs((sn+str+TString(".png")).Data());
 
}
 
 
 
if (gr2d && nx>1 && ny>1) {
 
 TCanvas *c = new TCanvas("c","QE2d",0,0,600,600);
 
 gr2d->Draw("colz");
 
 gr2d->SetNameTitle("gr2d",sn+swl+TString(";x [mm];y [mm]"));
 
;
 
 TString str;
 
 str.Form("_%d_nm",int(wl/10.));
 
 c->SaveAs((sn+str+TString("_xy.pdf")).Data());
 
 c->SaveAs((sn+str+TString("_xy.png")).Data());
 
}
 
f->Write();
 
//f->Close();
 
//delete f;
 
return 0; 
 
}
 
 
 
TH2F *wlnormalize(TH2F *h,int n){
 
 
 
TH2F *hn=(TH2F *)h->Clone();
 
int nx=h->GetNbinsX();
 
int ny=h->GetNbinsY();
 
float max=0;
 
int nymax=1;
 
for (int i=1;i<ny+1;i++){
 
  int nxy=h->GetBin(n,i);
 
  float val=h->GetBinContent(nxy);
 
  if (val>max){
 
     nymax=i;
 
     max=val;
 
  }
 
}
 
printf("Maximum at %d %f\n", nymax
,max
);  
for (int j=1;j<nx+1;j++){
 
int nxy0=h->GetBin(j,nymax);
 
float norm=h->GetBinContent(nxy0);
 
 
 
for (int i=1;i<ny+1;i++){
 
  int nxy=h->GetBin(j,i);
 
  float val=hn->GetBinContent(nxy);
 
  if (norm) hn->SetBinContent(nxy,val/norm);
 
  else  hn->SetBinContent(nxy,val/norm);
 
}
 
}
 
return hn;
 
}