#include <iostream>

#include <TROOT.h>
#include <TFile.h>
#include <TTree.h>
#include <TStyle.h>
#include <TGeoManager.h>
#include <TSpline.h>
#include <TH1D.h>
#include <TH2D.h>
#include <TCanvas.h>
#include <TPostScript.h>
#include <TRandom.h>

#include <ICOMETInput.hxx>
#include <ICOMETOutput.hxx>
#include <ICOMETEvent.hxx>
#include <IOADatabase.hxx>
#include <IHitSelection.hxx>
#include <IHit.hxx>

#include "TElecSim.hxx"
#include "TAMPDSim.hxx"

TROOT root("ROOT","Root of Everything");


void testOptics(){
  
  Bool_t verbose = kTRUE;
  /// write hitograms output in ps file
  TPostScript ampdData("OpticalResponse.ps",111);
  const Int_t Ncosmics = 1000;
  TH1D* hdir = new TH1D("hdir","hdir",100,0.,1.);
  TH1D* hedep = new TH1D("hedep","hedep",100,0.,10.);
  TH1D* hedep1 = new TH1D("hedep1","hedep1",100,0.,10.);
  TH1D* hnphot = new TH1D("hnphot","hnphot",100,0.,100000.);
  TH1D* hnphet = new TH1D("hnphet","hnphet",100,0.,100.);
  TAMPDSim ampd;
  // setting ampd basic parameters   
  ampd.SetNPixels(100);
  ampd.SetPDE(0.88, 0.95, 0.3);
  ampd.SetCrosstalk(0.2);
  ampd.SetPeWidth(0.02);
  if (verbose){
    std::cout<<"******** AMPDSim parameters : "<<std::endl; 
    std::cout<<"Integration time : "<<ampd.GetIntegrationTime()
	     <<" ns / Threshold : "<<ampd.GetThreshold()<<" pe"<<std::endl; 
  }
  Double_t orientation = 0;
  Double_t E0 = 0;
  Double_t nPhotons = 0;
  Double_t nPheTheo = 0;
  TArrayD ae(Ncosmics);
  TArrayD aphot(Ncosmics);
  TArrayD aphe(Ncosmics);
  for (Int_t j = 0; j <  Ncosmics; ++j){ 
    // one end readout  !!! 
    orientation = TMath::Cos(gRandom->Uniform(TMath::Pi()/4));
    orientation *=orientation;
    hdir->Fill(orientation);
    // Landau for MIP MPV = 2.01 MeV and FWHM = 0.285
    E0 = gRandom->Landau(2.01, 0.285);
    hedep->Fill(E0);
    E0 *= orientation; 
    hedep1->Fill(E0);
    nPhotons = ampd.GetNScintPhoton(E0)/2;
    hnphot->Fill(nPhotons);
    nPheTheo = ampd.GetNpheIdeal(nPhotons);    
    hnphet->Fill(nPheTheo);
    ae.AddAt(E0, j);
    aphot.AddAt(nPhotons, j);
    aphe.AddAt(nPheTheo, j);
  }
  TGraph* gbirks= new TGraph(Ncosmics, ae.GetArray(), aphot.GetArray());
  TGraph* gphe = new TGraph(Ncosmics, aphot.GetArray(), aphe.GetArray());
  TCanvas* cdat = new TCanvas("cdata","cdata",10,10,450,600);
  ampdData.NewPage();  
  cdat->Divide(2,3);  
  cdat->cd(1);
  hedep->Draw("");
  hedep->GetXaxis()->SetTitle("dEdx/cm MIP [MeV]");
  cdat->cd(2);
  hedep1->Draw("");
  hedep1->GetXaxis()->SetTitle("dEdx/cm/d#Omega_{hod} MIP [MeV]");
  cdat->cd(3);
  hnphot->Draw("");
  hnphot->GetXaxis()->SetTitle("#gamma_{Scint}(MIP)");
  cdat->cd(4);
  hnphet->Draw("");
  hnphet->GetXaxis()->SetTitle("Nphe ideal");
  cdat->Update();  
  cdat->cd(5);
  gbirks->SetMarkerStyle(4);
  gbirks->SetMarkerColor(2);
  gbirks->Draw("AP");
  gbirks->GetYaxis()->SetTitle("Number of #gamma_{scint}");
  gbirks->GetXaxis()->SetTitle("dEdx/cm/d#Omega_{hod} MIP [MeV]");
  cdat->cd(6);
  gphe->SetMarkerStyle(4);
  gphe->SetMarkerColor(4);
  gphe->Draw("AP");
  gphe->GetYaxis()->SetTitle("Nphe ideal");
  gphe->GetXaxis()->SetTitle("Number of #gamma_{scint}");
  ampdData.Close();  
  cdat->Close();
}

void testLEDResponse( Int_t pixels, Double_t pewidth, Double_t crosstalk,
		      Int_t intensity, Int_t npulse, const Char_t* filename){
  Bool_t verbose = kFALSE;

  /// write hitograms output in ps file
  TString outname = TString(filename) + TString(".ps");  
  TPostScript ampdData(outname,111);

  std::string outputName = std::string(filename) +  ".root";
  TFile* output = new TFile(outputName.c_str(), "NEW");
  if(!output->IsOpen()) {
    std::cout<<"Unable to open output file "<<outputName.c_str()<<std::endl;
    exit(1);
  }

  const Int_t NLEDPulse = npulse;
 
  TH1D* hnphot = new TH1D("hnphot","hnphot",100,0.,50000.);
  TH1D* hnphet = new TH1D("hnphet","hnphet",100,0.,10.);
  TH1D* hnpix = new TH1D("hnpix","hnpix",10,0.,10.);
  TH1D* hnpix_ct = new TH1D("hnpix_ct","hnpix_ct",100,0.,20.);
  TH1D* hnphereal_wo_ct = new TH1D("hnphereal_wo_ct","hnphereal_wo_ct",100,0.,10.);
  TH1D* hnphereal = new TH1D("hnphereal","hnphereal",100,0.,10.);

  TH2D* h2pix = new TH2D("h2pix","h2pix",100,0.,100.,100,0.,100.);
  TH2D* h2pixct = new TH2D("h2pixct","h2pixct",100,0.,100.,100,0.,100.);
  TH2D* h2phe = new TH2D("h2phe","h2phe",100,0.,100.,500,0.,100.);
  TH2D* h2phect = new TH2D("h2phect","h2phect",100,0.,100.,500,0.,100.);

  TAMPDSim ampd;
  // setting ampd basic parameters   
  ampd.SetNPixels(pixels);
  ampd.SetPDE(0.88, 0.95, 0.3);
  ampd.SetCrosstalk(crosstalk);
  ampd.SetPeWidth(pewidth);
  
  if (verbose){
    std::cout<<"******** AMPDSim parameters : "<<std::endl; 
    std::cout<<"Integration time : "<<ampd.GetIntegrationTime()
	     <<" ns / Threshold : "<<ampd.GetThreshold()<<" pe"<<std::endl; 
    std::cout<<"AMPD Npixels : "<<ampd.GetNPixels()<<std::endl; 
    std::cout<<"AMPD PDE "<<ampd.GetPDE()<<std::endl; 
    std::cout<<"AMPD 1 Phe resol "<<ampd.GetPeWidth()<<std::endl; 
    std::cout<<"AMPD Crosstalk level "<<ampd.GetCrosstalk()<<std::endl; 
  }

    Double_t nPhotons = 0;
    Double_t nPheTheo = 0;
    Double_t nPixFired = 0;
    Double_t nPixCrossTalk = 0;
    Double_t nPhe_w_ct = 0;
    Double_t nPheReal = 0;
    Double_t nPheReal_wo_ct = 0;

   TArrayD aphot(NLEDPulse);
   TArrayD aphe(NLEDPulse);
   TArrayD apix(NLEDPulse);
   TArrayD apixct(NLEDPulse);
   TArrayD aphereal(NLEDPulse);
   TArrayD apherealct(NLEDPulse);
   
   for (Int_t j = 0; j <  NLEDPulse; ++j){ 
     // LED pulse is considered as a gaussian signal with a mean value 
     // arbitrary set to give 4 phe in average 
     // 10% signal width
     nPhotons = gRandom->Gaus(70000.,30000.)*ampd.GetPDE()/2;
     hnphot->Fill(nPhotons); 
//      nPheTheo = ampd.GetNpheIdeal(nPhotons);
     nPheTheo = Double_t(gRandom->Poisson(3));
     hnphet->Fill(nPheTheo);
     aphot.AddAt(nPhotons, j);
     aphe.AddAt(nPheTheo, j);
     nPixFired = ampd.GetNpixelsFired(nPheTheo);
     hnpix->Fill(nPixFired);
     h2pix->Fill(nPheTheo, nPixFired);
     apix.AddAt(nPixFired, j);
     nPixCrossTalk = ampd.GetNpheCrosstalk(nPixFired);
     if (nPixCrossTalk > 0){
       hnpix_ct->Fill(nPixCrossTalk);
       h2pixct->Fill(nPheTheo, nPixCrossTalk);
       apixct.AddAt(nPixCrossTalk, j);
     }
     nPhe_w_ct =  nPixFired + nPixCrossTalk;
     nPheReal_wo_ct =  ampd.ResolutionSmearing(nPixFired);
     if (nPheReal_wo_ct > 0) {
       hnphereal_wo_ct->Fill(nPheReal_wo_ct); 
       aphereal.AddAt(nPheReal_wo_ct, j);      
       h2phe->Fill(nPheTheo, nPheReal_wo_ct);
     }
     nPheReal = ampd.ResolutionSmearing(nPhe_w_ct);
     if (nPheReal > 0) {
       hnphereal->Fill(nPheReal); 
       apherealct.AddAt(nPheReal, j);
       h2phe->Fill(nPheTheo, nPheReal);
     }     
     if (verbose){
       std::cout<<"-------------------------------------"<<std::endl; 
       std::cout<<"nPhotons 1 end RO "<<nPhotons<<std::endl; 
       std::cout<<"nPheTheo "<<nPheTheo<<std::endl; 
       std::cout<<"nPixFired "<<nPixFired<<std::endl; 
       std::cout<<"nPixCrossTalk "<<nPixCrossTalk<<std::endl; 
       std::cout<<"nPhe_w_ct "<<nPhe_w_ct<<std::endl; 
       std::cout<<"nPheReal_wo_ct "<<nPheReal_wo_ct<<std::endl; 
       std::cout<<"nPheReal "<<nPheReal<<std::endl; 
     }
   }
   
   TGraph* gphe = new TGraph(NLEDPulse, aphot.GetArray(), aphe.GetArray());
   TGraph* gpix= new TGraph(NLEDPulse, aphe.GetArray(), apix.GetArray());
   TGraph* gpixct = new TGraph(NLEDPulse, aphe.GetArray(), apixct.GetArray());
   TGraph* gphereal= new TGraph(NLEDPulse, aphe.GetArray(), aphereal.GetArray());
   TGraph* gpherealct = new TGraph(NLEDPulse, aphe.GetArray(), apherealct.GetArray());
   

   TCanvas* cdat = new TCanvas("cdata","cdata",10,10,550,400);
   ampdData.NewPage();  
   cdat->cd();
   hnphot->Draw("");
   hnphot->GetXaxis()->SetTitle("#gamma_{LED}()");
   cdat->Update();  

   ampdData.NewPage();  
   hnphet->Draw("");
   hnphet->GetXaxis()->SetTitle("Nphe ideal");
   cdat->Update();
  
   ampdData.NewPage();  
   hnpix->Draw("");
   hnpix->GetXaxis()->SetTitle("Nphe ideal");
   cdat->Print("npix_led.gif");
   cdat->Update();
   
   ampdData.NewPage();  
   hnpix_ct->Draw("");
   hnpix_ct->GetXaxis()->SetTitle("Induced crosstalk");
   cdat->Print("ct_led.gif");
   cdat->Update();

   ampdData.NewPage();  
   hnphereal_wo_ct->Draw("");
   hnphereal_wo_ct->GetXaxis()->SetTitle("Nphe wo CT ");
   cdat->Print("npherealwoct.gif");
   cdat->Update();

   ampdData.NewPage();  
   hnphereal->Draw("");
   hnphereal->GetXaxis()->SetTitle("Nphe w CT");
   cdat->Print("nphereal.gif");
   cdat->Update();
   
   ampdData.NewPage();  
   cdat->Clear();
   cdat->Divide(2,3);  
   cdat->cd(1);
   gphe->SetMarkerStyle(7);
   gphe->SetMarkerColor(2);
   gphe->Draw("AP");
   gphe->GetYaxis()->SetTitle("Nphe ideal");
   gphe->GetXaxis()->SetTitle("Number of #gamma_{scint}");
   cdat->cd(2);
   gpix->SetMarkerStyle(7);
   gpix->SetMarkerColor(2);
   gpix->Draw("AP");
   gpix->GetYaxis()->SetTitle("Pixels fired");
   gpix->GetXaxis()->SetTitle("Nphe ideal");
   cdat->cd(3);
   gphereal->SetMarkerStyle(7);
   gphereal->SetMarkerColor(6);
   gphereal->Draw("AP");
   gpix->GetYaxis()->SetTitle("NPhe real");
   gphereal->GetXaxis()->SetTitle("Nphe ideal");
   cdat->cd(4);
   gpherealct->SetMarkerStyle(7);
   gpherealct->SetMarkerColor(4);
   gpherealct->Draw("AP");
   gpherealct->GetYaxis()->SetTitle("NPhe real w CT");
   gpherealct->GetXaxis()->SetTitle("Nphe ideal");
   cdat->cd(5);
   gpixct->SetMarkerStyle(7); 
   gpixct->SetMarkerColor(2);
   gpixct->Draw("AP");
   cdat->Update();
   ampdData.Close();  
   if (output) {
     output->cd();
     hnphot->Write();
     hnphet->Write();
     hnpix->Write();
     hnpix_ct->Write();
     hnphereal->Write();
     hnphereal_wo_ct->Write();
     h2pix->Write();
     h2pixct->Write();
     h2phe->Write();
     h2phect->Write();
     output->Close();
   }
   cdat->Close();
}

void testMIPResponse(){

  Bool_t verbose = kTRUE;

  /// write hitograms output in ps file
  TPostScript ampdData("MIPresponse.ps",111);

  std::string outputName = "MIPresponse.root";
  TFile* output = new TFile(outputName.c_str(), "NEW");
  if(!output->IsOpen()) {
    std::cout<<"Unable to open output file "<<outputName.c_str()<<std::endl;
    exit(1);
  }

  const Int_t Ncosmics = 5000;
 

  TH1D* hdir = new TH1D("hdir","hdir",100,0.,1.);
  TH1D* hedep = new TH1D("hedep","hedep",100,0.,10.);
  TH1D* hedep1 = new TH1D("hedep1","hedep1",100,0.,10.);
  TH1D* hnphot = new TH1D("hnphot","hnphot",100,0.,100000.);
  TH1D* hnphet = new TH1D("hnphet","hnphet",100,0.,100.);
  TH1D* hnpix = new TH1D("hnpix","hnpix",40,0.,40.);
  TH1D* hnpix_ct = new TH1D("hnpix_ct","hnpix_ct",100,0.,20.);
  TH1D* hnphereal_wo_ct = new TH1D("hnphereal_wo_ct","hnphereal_wo_ct",600,0.,60.);
  TH1D* hnphereal = new TH1D("hnphereal","hnphereal",600,0.,60.);

  TH2D* h2pix = new TH2D("h2pix","h2pix",100,0.,100.,100,0.,100.);
  TH2D* h2pixct = new TH2D("h2pixct","h2pixct",100,0.,100.,100,0.,100.);
  TH2D* h2phe = new TH2D("h2phe","h2phe",100,0.,100.,500,0.,100.);
  TH2D* h2phect = new TH2D("h2phect","h2phect",100,0.,100.,500,0.,100.);

  TAMPDSim ampd;
  // setting ampd basic parameters   
  ampd.SetNPixels(400);
  ampd.SetPDE(0.88, 0.95, 0.3);
  ampd.SetCrosstalk(0.2);
  ampd.SetPeWidth(0.1);
  
  if (verbose){
    std::cout<<"******** AMPDSim parameters : "<<std::endl; 
    std::cout<<"Integration time : "<<ampd.GetIntegrationTime()
	     <<" ns / Threshold : "<<ampd.GetThreshold()<<" pe"<<std::endl; 
    std::cout<<"AMPD Npixels : "<<ampd.GetNPixels()<<std::endl; 
    std::cout<<"AMPD PDE "<<ampd.GetPDE()<<std::endl; 
    std::cout<<"AMPD Crosstalk level "<<ampd.GetCrosstalk()<<std::endl; 
  }

    Double_t orientation = 0;
    Double_t E0 = 0;
    Double_t nPhotons = 0;
    Double_t nPheTheo = 0;
    Double_t nPixFired = 0;
    Double_t nPixCrossTalk = 0;
    Double_t nPhe_w_ct = 0;
    Double_t nPheReal = 0;
    Double_t nPheReal_wo_ct = 0;

   TArrayD ae(Ncosmics);
   TArrayD aphot(Ncosmics);
   TArrayD aphe(Ncosmics);
   TArrayD apix(Ncosmics);
   TArrayD apixct(Ncosmics);
   TArrayD aphereal(Ncosmics);
   TArrayD apherealct(Ncosmics);
   
   for (Int_t j = 0; j <  Ncosmics; ++j){ 
      // Energy deposited for a MIP in 1cm of scintillator plastic
      // 1.68 MeV.g-1.cm2 and density of scintillator 1.032
      // add a cosine smearing between 0 and 45 deg
      orientation = TMath::Cos(gRandom->Uniform(TMath::Pi()/4));
      orientation *=orientation;
      hdir->Fill(orientation);
      // Landau for MIP MPV = 2.01 MeV and FWHM = 0.285
      E0 = gRandom->Landau(2.01, 0.285);
      hedep->Fill(E0);
      E0 *= orientation; 
      hedep1->Fill(E0);
      // one end readout  !!! 
      nPhotons = ampd.GetNScintPhoton(E0)/2;
      hnphot->Fill(nPhotons);
      
      nPheTheo = ampd.GetNpheIdeal(nPhotons);
      hnphet->Fill(nPheTheo);
      ae.AddAt(E0, j);
      aphot.AddAt(nPhotons, j);
      aphe.AddAt(nPheTheo, j);
    
      nPixFired = ampd.GetNpixelsFired(nPheTheo);
      hnpix->Fill(nPixFired);
      h2pix->Fill(nPheTheo, nPixFired);
      
      apix.AddAt(nPixFired, j);
      
      nPixCrossTalk = ampd.GetNpheCrosstalk(nPixFired);
      if (nPixCrossTalk > 0){
	hnpix_ct->Fill(nPixCrossTalk);
	h2pixct->Fill(nPheTheo, nPixCrossTalk);
	apixct.AddAt(nPixCrossTalk, j);
      }
      nPhe_w_ct =  nPixFired + nPixCrossTalk;
      
      nPheReal_wo_ct =  ampd.ResolutionSmearing(nPixFired);
      if (nPheReal_wo_ct > 0) {
	hnphereal_wo_ct->Fill(nPheReal_wo_ct); 
	aphereal.AddAt(nPheReal_wo_ct, j);      
	h2phe->Fill(nPheTheo, nPheReal_wo_ct);

      }

      
      nPheReal = ampd.ResolutionSmearing(nPhe_w_ct);
      if (nPheReal > 0) {
	hnphereal->Fill(nPheReal); 
	apherealct.AddAt(nPheReal, j);
	h2phe->Fill(nPheTheo, nPheReal);
      }
      
      if (verbose){
	std::cout<<"-------------------------------------"<<std::endl; 
	std::cout<<"orientation "<<orientation<<std::endl; 
	std::cout<<"E0 "<<E0<<std::endl; 
	std::cout<<"nPhotons 1 end RO "<<nPhotons<<std::endl; 
	std::cout<<"nPheTheo "<<nPheTheo<<std::endl; 
	std::cout<<"nPixFired "<<nPixFired<<std::endl; 
	std::cout<<"nPixCrossTalk "<<nPixCrossTalk<<std::endl; 
	std::cout<<"nPhe_w_ct "<<nPhe_w_ct<<std::endl; 
	std::cout<<"nPheReal_wo_ct "<<nPheReal_wo_ct<<std::endl; 
	std::cout<<"nPheReal "<<nPheReal<<std::endl; 
      }
    }
    
    TGraph* gbirk= new TGraph(Ncosmics, ae.GetArray(), aphot.GetArray());
    TGraph* gphe = new TGraph(Ncosmics, aphot.GetArray(), aphe.GetArray());
    TGraph* gpix= new TGraph(Ncosmics, aphe.GetArray(), apix.GetArray());
    TGraph* gpixct = new TGraph(Ncosmics, aphe.GetArray(), apixct.GetArray());
    TGraph* gphereal= new TGraph(Ncosmics, aphe.GetArray(), aphereal.GetArray());
    TGraph* gpherealct = new TGraph(Ncosmics, aphe.GetArray(), apherealct.GetArray());
 

    TCanvas* cdat = new TCanvas("cdata","cdata",10,10,450,600);
    ampdData.NewPage();  
    cdat->Divide(2,2);  
    cdat->cd(1);
    hdir->Draw("");
    cdat->cd(2);
    hedep1->Draw("");
    hedep1->GetXaxis()->SetTitle("dEdx/cm/d#Omega_{hod} MIP [MeV]");
    cdat->cd(3);
    hnphot->Draw("");
    hnphot->GetXaxis()->SetTitle("#gamma_{Scint}(MIP)");
    cdat->cd(4);
    hnphet->Draw("");
    hnphet->GetXaxis()->SetTitle("Nphe ideal");
    cdat->Update();  
      
    ampdData.NewPage();  
    cdat->cd(1);
    hnpix->Draw("");
    hnpix->GetXaxis()->SetTitle("Nphe ideal");
    cdat->cd(2);
    hnpix_ct->Draw("");
    hnpix_ct->GetXaxis()->SetTitle("Induced crosstalk");
    cdat->Update();
    cdat->cd(3);
    hnphereal_wo_ct->Draw("");
    hnphereal_wo_ct->GetXaxis()->SetTitle("Nphe wo CT ");
    cdat->cd(4);
    hnphereal->Draw("");
    hnphereal->GetXaxis()->SetTitle("Nphe w CT");
    cdat->Update();
  
   ampdData.NewPage();  
   cdat->Clear();
   cdat->cd();
   gbirk->SetMarkerStyle(7);
   gbirk->SetMarkerColor(2);
   gbirk->Draw("AP");
   gbirk->GetYaxis()->SetTitle("Number of #gamma_{scint}");
   gbirk->GetXaxis()->SetTitle("dEdx/cm/d#Omega_{hod} MIP [MeV]");
   cdat->Update();
  
   ampdData.NewPage();  
   cdat->Clear();
   gphe->SetMarkerStyle(7);
   gphe->SetMarkerColor(2);
   gphe->Draw("AP");
   gphe->GetYaxis()->SetTitle("Nphe ideal");
   gphe->GetXaxis()->SetTitle("Number of #gamma_{scint}");
   cdat->Update();

   ampdData.NewPage();  
   cdat->Clear();
   gpix->SetMarkerStyle(7);
   gpix->SetMarkerColor(2);
   gpix->Draw("AP");
   gpix->GetXaxis()->SetTitle("Nphe ideal");
   gpix->GetYaxis()->SetTitle("Pixels fired");
   cdat->Update();

   
   ampdData.NewPage();  
   cdat->Clear();
   gphereal->SetMarkerStyle(7);
   gphereal->SetMarkerColor(6);
   gphereal->Draw("AP");
   gphereal->GetYaxis()->SetTitle("NPhe real");
   gphereal->GetXaxis()->SetTitle("Nphe ideal");
   cdat->Update();

   ampdData.NewPage();  
   cdat->Clear();
   gpherealct->SetMarkerStyle(7);
   gpherealct->SetMarkerColor(4);
   gpherealct->Draw("AP");
   gpherealct->GetYaxis()->SetTitle("NPhe real w CT");
   gpherealct->GetXaxis()->SetTitle("Nphe ideal");
   cdat->Update();

   ampdData.NewPage();  
    cdat->Clear();
  gpixct->SetMarkerStyle(7); 
   gpixct->SetMarkerColor(2);
   gpixct->Draw("AP");
   gpixct->GetYaxis()->SetTitle("CT level (phe unit)");
   gpixct->GetXaxis()->SetTitle("Nphe ideal");
   cdat->Update();
   ampdData.Close();  
   if (output) {
     output->cd();
     hdir->Write();
     hedep->Write();
     hedep1->Write();
     hnphot->Write();
     hnphet->Write();
     hnpix->Write();
     hnpix_ct->Write();
     hnphereal->Write();
     hnphereal_wo_ct->Write();
     h2pix->Write();
     h2pixct->Write();
     h2phe->Write();
     h2phect->Write();
     output->Close();
   }
   cdat->Close();
}
  
int main(int argc, char **argv) {
  
 testOptics();
 testLEDResponse(500,0.08,0.1,3,10000,"LED_AMPD_500pix");
 testMIPResponse();
 return 0;
}