UserModule.cc

Go to the documentation of this file.

#include <sstream>

#include "UserModule.h"

#include <det/Detector.h>

#include <fdet/FDetector.h>
#include <fdet/Telescope.h>
#include <fdet/Eye.h>

#include <utl/ErrorLogger.h>
#include <utl/TimeStamp.h>
#include <utl/TimeInterval.h>
#include <utl/TabulatedFunction.h>
#include <utl/MultiTabulatedFunction.h>
#include <utl/Photon.h>
#include <utl/AugerUnits.h>
#include <utl/MathConstants.h>
#include <utl/PhysicalConstants.h>
#include <utl/RandomEngine.h>
#include <utl/Trace.h>
#include <utl/TraceAlgorithm.h>
#include <utl/Reader.h>

#include <evt/Event.h>
#include <evt/ShowerSimData.h>

#include <fevt/FEvent.h>
#include <fevt/Eye.h>
#include <fevt/TelescopeSimData.h>
#include <fevt/Telescope.h>
#include <fevt/PixelSimData.h>
#include <fevt/Pixel.h>
#include <fevt/ChannelSimData.h>
#include <fevt/Channel.h>

#include <fwk/CentralConfig.h>

#include <TFile.h>
#include <TObjArray.h>
#include <TGraph.h>
#include <TGraphErrors.h>
#include <TMultiGraph.h>
#include <TF1.h>
#include <TH2D.h>
#include <TVector3.h>
#include <TMarker.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TTree.h>
#include <TFolder.h>
#include <TString.h>
#include <TPaveText.h>
#include <TColor.h>
#include <TROOT.h>
#include <TPolyLine.h>

#include <CLHEP/Random/Randomize.h>

using namespace std;
using namespace utl;
using namespace evt;
using namespace fevt;
using namespace fwk;
using namespace det;


UserModule::UserModule() :
  fOutputFile(0),
  fTree(0),
  fLongitudinalProfileCanvas(0),
  fEnergyDepositCanvas(0),
  fFluorLightCanvas(0),
  fCherBeamCanvas(0),
  fLightFluxCanvas(0),
  fFluorFluxCanvas(0),
  fCherMieFluxCanvas(0),
  fCherRaylFluxCanvas(0),
  fCherDirFluxCanvas(0),
  fLightOnCameraCanvas(0),
  fLightAtDiaCanvas(0),
  fCameraCanvas(0),
  fPhotonGraphs(0),
  fFADCGraphs(0),
  fEventId(0)
{ }


UserModule::~UserModule()
{
  // delete fOutputFile;
  // delete fTree;
  // delete fLongitudinalProfileCanvas;
  // delete fEnergyDepositCanvas;
  // delete fFluorLightCanvas;
  // delete fCherBeamCanvas;
  // delete fLightFluxCanvas;
  // delete fFluorFluxCanvas;
  // delete fCherMieFluxCanvas;
  // delete fCherRaylFluxCanvas;
  // delete fCameraCanvas;
  // delete fLightOnCameraCanvas;
  // delete fLightAtDiaCanvas;
  // delete fEventId;
  // delete fPhotonGraphs;
  // delete fFADCGraphs;
}


VModule::ResultFlag
UserModule::Init()
{
  //
  // Initialize your module here. This method
  // is called once at the beginning of the run.
  // The eSuccess flag indicates the method ended
  // successfully.  For other possible return types,
  // see the VModule documentation.
  //
  INFO("UserModule::Init()");

  CentralConfig* const cc = CentralConfig::GetInstance();

  Branch topB = cc->GetTopBranch("UserModule");

  Branch eyeB = topB.GetChild("eye");
  if (!eyeB) {
    ERROR("Could not find branch eye");
    throw utl::NonExistentComponentException();
  }
  eyeB.GetData(fEye);

  Branch telB = topB.GetChild("telescope");
  if (!telB) {
    ERROR("Could not find branch telescope");
    throw utl::NonExistentComponentException();
  }
  telB.GetData(fTelescope);

  fOutputFile = new TFile("OfflineOutput.root","recreate");
  fLongitudinalProfileCanvas = new TCanvas("LongitudinalProfileC",
                                           "LongitudinalProfileC");
  fEnergyDepositCanvas = new TCanvas("EnergyDepositC",
                                     "EnergyDepositC");
  fFluorLightCanvas = new TCanvas("FluorLightC", "FluorLightC");
  fCherBeamCanvas = new TCanvas("CherBeamC", "CherBeamC");

  fLightFluxCanvas = new TCanvas("LightFlux", "LightFlux");

  fFluorFluxCanvas = new TCanvas("FluorFlux", "FluorFlux");

  fCherMieFluxCanvas = new TCanvas("CherMieFlux", "CherMieFlux");

  fCherRaylFluxCanvas = new TCanvas("CherRaylFlux", "CherRaylFlux");

  fCherDirFluxCanvas = new TCanvas("CherDirFlux", "CherDirFlux");

  fLightOnCameraCanvas = new TCanvas("LightCamera", "LightCamera");

  fLightAtDiaCanvas = new TCanvas("LightAtDiaphragm", "LightAtDiaphragm");

  fCameraCanvas = new TCanvas("Camera", "Camera",60,60,600,600);

  fPhotonGraphs = new TObjArray;

  fFADCGraphs = new TObjArray;

  fEventId = new TString;

  return eSuccess;
}


VModule::ResultFlag
UserModule::Run(Event& event)
{
  INFO("UserModule::Run()");

  if (event.HasSimShower()) {
    // AnalyzeShower(event.GetSimShower());
    AnalyzeLightAtDiaphragm(event.GetFEvent());
    AnalyzeCamera(event.GetFEvent());
    return eSuccess;
  }

  ERROR("Event does not contain a simulated shower.\n"
        "Going to the next event.");
  return eContinueLoop;
}


VModule::ResultFlag
UserModule::Finish()
{
  INFO("UserModule::Finish()");

  fOutputFile->cd();

  fOutputFile->Close();

  return eSuccess;
}


void
UserModule::ClearAllPlots()
{
  fLongitudinalProfileCanvas->Clear();
  fEnergyDepositCanvas->Clear();
  fCherBeamCanvas->Clear();
  fFluorLightCanvas->Clear();
  fLightFluxCanvas->Clear();
  fFluorFluxCanvas->Clear();
  fCherMieFluxCanvas->Clear();
  fCherRaylFluxCanvas->Clear();
  fCherDirFluxCanvas->Clear();
  fCameraCanvas->Clear();
  fLightOnCameraCanvas->Clear();
  fLightAtDiaCanvas->Clear();
}


void
UserModule::AnalyzeShower(const evt::ShowerSimData& shower)
{
  //utl::TabulatedFunction fluor = shower.GetFluorescencePhotons(0);
  //utl::TabulatedFunction cher = shower.GetCherenkovPhotons(0);

  // GH data of the shower
  /*
    cout << "---------------------------" << endl;
    cout << "GH Parameters of the Shower" << endl;
    cout << "---------------------------" << endl;
    cout << "Nmax     = " << shower.GetGHParameters().GetNMax() << endl;
    cout << "Xmax     = " << shower.GetGHParameters().GetXMax()/(g/(cm*cm))
    << " g/cm²" << endl;
    //cout << "X0       = " << gh.GetXZero()/(g/(cm*cm)) << " g/cm²" << endl;
    cout << "chi²     = " << shower.GetGHParameters().GetChiSquare() << endl
    << endl;
  */

  PlotLongitudinalProfile(shower);
  PlotEnergyDeposit(shower);
  PlotTotalFluorLight(shower);
  PlotTotalCherBeam(shower);
}


void
UserModule::PlotLongitudinalProfile(const evt::ShowerSimData& shower)
{
  utl::TabulatedFunction profile = shower.GetLongitudinalProfile();
  unsigned int size = profile.GetNPoints();
  double x[size];
  double y[size];
  for (unsigned int i = 0; i < size; ++i) {
    x[i] = profile[i].X()/(g/cm2);
    y[i] = profile[i].Y();
  }

  fLongitudinalProfileCanvas->cd();

  TGraph* const graph = new TGraph(size, x, y);
  graph->SetMarkerStyle(8);
  graph->Draw("APL");
  graph->SetTitle("Longitudinal Charge Profile");
  //BRD 18/2/05
  //  graph->GetXaxis()->SetTitle("Vertical depth (g/cm²)");
  graph->GetXaxis()->SetTitle("slant depth (g/cm²)");
  graph->GetYaxis()->SetTitle("#Particles");

  TString s;
  fOutputFile->cd();
  s = "LongitudinalProfileC_";
  s += *fEventId;
  fLongitudinalProfileCanvas->SetName(s.Data());
  fLongitudinalProfileCanvas->Write();
}


void
UserModule::PlotEnergyDeposit(const evt::ShowerSimData& shower)
{
  if (!shower.HasdEdX())
    return;
  utl::TabulatedFunction edep = shower.GetdEdX();
  unsigned int size = edep.GetNPoints();
  double x[size];
  double y[size];
  //double delta = abs(edep[0].X() - edep[1].X());
  for (unsigned int i = 0; i < size; ++i) {
    x[i] = edep[i].X()/(g/(cm*cm));
    //BRD 18/2/05
    //      y[i] = edep[i].Y()/(eV/(g/cm/cm))/(delta/(g/(cm*cm)));
    y[i] = edep[i].Y()/(eV/(g/cm2));
  }

  fEnergyDepositCanvas->cd();

  TGraph* const graph = new TGraph(size, x, y);
  graph->SetMarkerStyle(8);
  graph->Draw("APL");
  graph->SetTitle("EnergyDeposit");
  //BRD 18/2/05
  //  graph->GetXaxis()->SetTitle("Vertical depth (g/cm²)");
  graph->GetXaxis()->SetTitle("slant depth (g/cm^{2})");
  graph->GetYaxis()->SetTitle("dE/dX (eV/(g/cm^{2}))");

  TString s;
  fOutputFile->cd();
  s = "EnergyDepositC_";
  s += *fEventId;
  fEnergyDepositCanvas->SetName(s.Data());
  fEnergyDepositCanvas->Write();
}


void
UserModule::PlotTotalFluorLight(const evt::ShowerSimData& shower)
{
  Detector& theDetector = Detector::GetInstance();

  const std::vector<double>& waves =
    theDetector.GetAtmosphere().GetWavelengths();

  const unsigned int nWaves = waves.size();
  utl::TabulatedFunction part1 = shower.GetFluorescencePhotons(0);
  unsigned int size = part1.GetNPoints();
  double x[size];
  double y[size];
  for (unsigned int i = 0; i < size; ++i) {
    x[i] = 0;
    y[i] = 0;
  }
  for (unsigned int iwl = 0; iwl < nWaves; ++iwl) {
    utl::TabulatedFunction part = shower.GetFluorescencePhotons(iwl);
    for (unsigned int i = 0; i < size; ++i) {
      if (!iwl)
        x[i] = part[i].X()/(g/cm2);
      y[i] += part[i].Y();
    }
  }

  fFluorLightCanvas->cd();

  TGraph* const graph = new TGraph(size, x, y);
  graph->SetMarkerStyle(8);
  graph->Draw("APL");
  graph->SetTitle("Fluorescence light");
  //BRD 18/2/05
  //  graph->GetXaxis()->SetTitle("Vertical depth (g/cm²)");
  graph->GetXaxis()->SetTitle("slant depth (g/cm^{2})");
  graph->GetYaxis()->SetTitle("#Photons");

  TString s;
  fOutputFile->cd();
  s = "FluorLightC_";
  s += *fEventId;
  fFluorLightCanvas->SetName(s.Data());
  fFluorLightCanvas->Write();
}


void
UserModule::PlotTotalCherBeam(const evt::ShowerSimData& shower)
{
  Detector& theDetector = Detector::GetInstance();

  const std::vector<double>& waves =
    theDetector.GetAtmosphere().GetWavelengths();

  const unsigned int nWaves = waves.size();

  utl::TabulatedFunction ph = shower.GetCherenkovPhotons(0);
  unsigned int size = ph.GetNPoints();
  double x[size];
  double y[size];
  for (unsigned int i = 0; i < size; ++i) {
    //x[i] = theDetector.GetAtmosphere().EvaluateHeight(ph[i].X());
    x[i] = ph[i].X()/(g/cm2);
    y[i] = 0;
  }

  for (unsigned int iwl = 0; iwl < nWaves; ++iwl) {

    utl::TabulatedFunction ph = shower.GetCherenkovPhotons(iwl);

    for (unsigned int i = 0; i < size; ++i) {
      y[i] += ph[i].Y();
    }
  }

  fCherBeamCanvas->cd();

  TGraph* const graph = new TGraph(size, x, y);
  graph->SetMarkerStyle(8);
  graph->Draw("APL");
  graph->SetTitle("Cherenkov beam");
  //BRD 18/2/05
  //  graph->GetXaxis()->SetTitle("Vertical depth (g/cm²)");
  graph->GetXaxis()->SetTitle("slant depth (g/cm^{2})");
  graph->GetYaxis()->SetTitle("#Photons");

  TString s;
  fOutputFile->cd();
  s = "CherenkovBeamC_";
  s += *fEventId;
  fCherBeamCanvas->SetName(s.Data());
  fCherBeamCanvas->Write();
}


void
UserModule::AnalyzeLightAtDiaphragm(const fevt::FEvent& fevent)
{
  // Get the simulated longitudinal shower profile
  if (fevent.HasEye(fEye))
    if (fevent.GetEye(fEye).HasTelescope(fTelescope))
      if (fevent.GetEye(fEye).GetTelescope(fTelescope).HasSimData()) {
        const fevt::TelescopeSimData& sim =
          fevent.GetEye(fEye).GetTelescope(fTelescope).GetSimData();
        PlotLightFlux(sim);
        PlotFluorFlux(sim);
        PlotCherMieFlux(sim);
        PlotCherRaylFlux(sim);
        PlotCherDirFlux(sim);
        ClearAllPlots();
      }
}


void
UserModule::PlotLightFlux(const fevt::TelescopeSimData& sim)
{
  const int iwl = 0; // for realistic answer you have to loop over all wavelength and sum
                     // them up taking into account the wavelength dependence of the overall
                     // telescope efficiency

  if (!sim.HasPhotonTrace(fevt::FdConstants::eFluorDirect, iwl))
    return;

  const utl::TraceD& fluorTrace    = sim.GetPhotonTrace(fevt::FdConstants::eFluorDirect, iwl);
  const utl::TraceD& cherDirTrace  = sim.GetPhotonTrace(fevt::FdConstants::eCherDirect, iwl);
  const utl::TraceD& cherRaylTrace =
    sim.GetPhotonTrace(fevt::FdConstants::eCherRayleighScattered, iwl);
  const utl::TraceD& cherMieTrace  =
    sim.GetPhotonTrace(fevt::FdConstants::eCherMieScattered, iwl);
  double bin = fluorTrace.GetBinning();
  const int size = int(bin*fluorTrace.GetSize()/100.0);
  double x[size];
  double yt[size];
  double yf[size];
  double ycd[size];
  double ycr[size];
  double ycm[size];

  for (int i = 0; i < size; ++i) {
    x[i]   = i;
    yt[i]  = 0;
    yf[i]  = 0;
    ycd[i] = 0;
    ycr[i] = 0;
    ycm[i] = 0;
  }

  int j = 0;
  for (unsigned int i = fluorTrace.GetStart(); i < fluorTrace.GetStop(); ++i) {

    yf[j] += fluorTrace[i];
    ycd[j] += cherDirTrace[i];
    ycr[j] += cherRaylTrace[i];
    ycm[j] += cherMieTrace[i];

    if (!(i % int(100/bin)))
      ++j;

    if (j > size-1)
      break;

  }

  // Total flux

  for (int i = 0; i < size; ++i) {
    yt[i] = yf[i] + ycm[i] + ycr[i] + ycd[i];
  }

  fLightFluxCanvas->cd();

  TGraph* const gt = new TGraph(size, x, yt);
  gt->SetLineColor(1);
  gt->Draw("AL");
  gt->SetTitle("Light flux");
  gt->GetXaxis()->SetTitle("bin [100 ns]");
  gt->GetXaxis()->SetRange(0,size);
  gt->GetYaxis()->SetTitle("#Photons");

  TGraph* const gf = new TGraph(size, x, yf);
  gf->SetLineColor(3);
  gf->Draw("L");

  TGraph* const gcm = new TGraph(size, x, ycm);
  gcm->SetLineColor(6);
  gcm->Draw("L");

  TGraph* const gcr = new TGraph(size, x, ycr);
  gcr->SetLineColor(4);
  gcr->Draw("L");

  TGraph* const gcd = new TGraph(size, x, ycd);
  gcd->SetLineColor(2);
  gcd->Draw("L");

  TString s;
  fOutputFile->cd();
  s = "LightFluxC_";
  s += *fEventId;
  fLightFluxCanvas->SetName(s.Data());
  fLightFluxCanvas->Write();
}


void
UserModule::PlotFluorFlux(const fevt::TelescopeSimData& sim)
{
  const int iwl = 0; // for realistic answer you have to loop over all wavelength and sum
                     // them up taking into account the wavelength dependence of the overall
                     // telescope efficiency

  if (!sim.HasPhotonTrace(fevt::FdConstants::eFluorDirect, iwl))
    return;

  const utl::TraceD& fluorTrace = sim.GetPhotonTrace(fevt::FdConstants::eFluorDirect, iwl);
  const double bin = fluorTrace.GetBinning();
  const int size = int(bin*fluorTrace.GetSize()/100.0);
  double x[size];
  double y[size];

  for (int i = 0; i < size; ++i) {
    x[i] = i;
    y[i] = 0;
  }

  int j = 0;
  for (unsigned int i = fluorTrace.GetStart(); i < fluorTrace.GetStop(); ++i) {

    y[j] += fluorTrace[i];
    if (i%(int(100/bin)) == 0)
      ++j;

    if (j > size-1)
      break;
  }

  fFluorFluxCanvas->cd();

  TGraph* const gt = new TGraph(size, x, y);
  gt->SetLineColor(1);
  gt->Draw("AL");
  gt->SetTitle("Fluorescence light flux");
  gt->GetXaxis()->SetTitle("bin [100 ns]");
  gt->GetXaxis()->SetRange(0,size);
  gt->GetYaxis()->SetTitle("#Photons");

  TString s;
  fOutputFile->cd();
  s = "FluorFluxC_";
  s += *fEventId;
  fFluorFluxCanvas->SetName(s.Data());
  fFluorFluxCanvas->Write();
}


void
UserModule::PlotCherMieFlux(const fevt::TelescopeSimData& sim)
{
  const int iwl = 0; // for realistic answer you have to loop over all wavelength and sum
                     // them up taking into account the wavelength dependence of the overall
                     // telescope efficiency

  if (!sim.HasPhotonTrace(fevt::FdConstants::eCherMieScattered, iwl))
    return;

  const utl::TraceD& trace = sim.GetPhotonTrace(fevt::FdConstants::eCherMieScattered, iwl);
  const double bin = trace.GetBinning();
  const int size = int(bin*trace.GetSize()/100.0);
  double x[size];
  double y[size];

  for (int i = 0; i < size; ++i) {
    x[i] = i;
    y[i] = 0;
  }

  int j = 0;
  for (unsigned int i = trace.GetStart(); i < trace.GetStop(); ++i) {

    y[j] += trace[i];
    if (i%(int(100/bin)) == 0)
      ++j;

    if (j > size-1)
      break;
  }

  fCherMieFluxCanvas->cd();

  TGraph* const gt = new TGraph(size, x, y);
  gt->SetLineColor(1);
  gt->Draw("AL");
  gt->SetTitle("Cherenkov Mie-scattered light flux");
  gt->GetXaxis()->SetTitle("bin [100 ns]");
  gt->GetXaxis()->SetRange(0,size);
  gt->GetYaxis()->SetTitle("#Photons");

  TString s;
  fOutputFile->cd();
  s = "CherMieFluxC_";
  s += *fEventId;
  fCherMieFluxCanvas->SetName(s.Data());
  fCherMieFluxCanvas->Write();
}


void
UserModule::PlotCherRaylFlux(const fevt::TelescopeSimData& sim)
{
  const int iwl = 0; // for realistic answer you have to loop over all wavelength and sum
                     // them up taking into account the wavelength dependence of the overall
                     // telescope efficiency

  if (!sim.HasPhotonTrace(fevt::FdConstants::eCherRayleighScattered, iwl))
    return;

  const utl::TraceD& trace = sim.GetPhotonTrace(fevt::FdConstants::eCherRayleighScattered, iwl);
  const double bin = trace.GetBinning();
  const int size = int(bin*trace.GetSize()/100.0);
  double x[size];
  double y[size];

  for (int i = 0; i < size; ++i) {
    x[i] = i;
    y[i] = 0;
  }

  int j = 0;
  for (unsigned int i = trace.GetStart(); i < trace.GetStop(); ++i) {

    y[j] += trace[i];
    if (i%(int(100/bin)) == 0)
      ++j;

    if (j > size-1)
      break;
  }

  fCherRaylFluxCanvas->cd();

  TGraph* const gt = new TGraph(size, x, y);
  gt->SetLineColor(1);
  gt->Draw("AL");
  gt->SetTitle("Cherenkov Rayleigh-scattered light flux");
  gt->GetXaxis()->SetTitle("bin [100 ns]");
  gt->GetXaxis()->SetRange(0,size);
  gt->GetYaxis()->SetTitle("#Photons");

  TString s;
  fOutputFile->cd();
  s = "CherRayleighFluxC_";
  s += *fEventId;
  fCherRaylFluxCanvas->SetName(s.Data());
  fCherRaylFluxCanvas->Write();
}


void
UserModule::PlotCherDirFlux(const fevt::TelescopeSimData& sim)
{
  const int iwl = 0; // for realistic answer you have to loop over all wavelength and sum
                     // them up taking into account the wavelength dependence of the overall
                     // telescope efficiency

  if (!sim.HasPhotonTrace(fevt::FdConstants::eCherDirect, iwl))
    return;

  const utl::TraceD& trace = sim.GetPhotonTrace(fevt::FdConstants::eCherDirect, iwl);
  const double bin = trace.GetBinning();
  const int size = int(bin*trace.GetSize()/100.0);
  double x[size];
  double y[size];
  for (int i = 0; i < size; ++i) {
    x[i] = i;
    y[i] = 0;
  }
  int j = 0;
  for (unsigned int i = trace.GetStart(); i < trace.GetStop(); ++i) {

    y[j] += trace[i];
    if (i%(int(100/bin)) == 0)
      ++j;

    if (j > size-1)
      break;
  }

  fCherDirFluxCanvas->cd();

  TGraph* const gt = new TGraph(size, x, y);
  gt->SetLineColor(1);
  gt->Draw("AL");
  gt->SetTitle("Cherenkov direct light flux");
  gt->GetXaxis()->SetTitle("bin [100 ns]");
  gt->GetXaxis()->SetRange(0,size);
  gt->GetYaxis()->SetTitle("#Photons");

  TString s;
  fOutputFile->cd();
  s = "CherDirectFluxC_";
  s += *fEventId;
  fCherDirFluxCanvas->SetName(s.Data());
  fCherDirFluxCanvas->Write();
}


void
UserModule::AnalyzeCamera(const fevt::FEvent& fevent)
{
  // Get the simulated longitudinal shower profile
  if (fevent.HasEye(fEye))
    if (fevent.GetEye(fEye).HasTelescope(fTelescope)) {
      const fevt::Telescope& tel = fevent.GetEye(fEye).GetTelescope(fTelescope);
      PlotLightOnCamera(tel);
      PlotCamera(tel);
      PlotPhotonTraces(tel);
      PlotFADCTraces(tel);
      //AnalyzeLightAtDiaphragm(fevent);
      //  PlotLightAtDia(tel);
      ClearAllPlots();
    }
}


void
UserModule::PlotLightAtDia(const fevt::Telescope& tel)
{
  const int iwl = 0; // for realistic answer you have to loop over all wavelength and sum
                     // them up taking into account the wavelength dependence of the overall
                     // telescope efficiency

  cout << " graph plot at dia 0 " << endl;
  const fevt::TelescopeSimData& telSimData = tel.GetSimData();
  cout << " graph plot at dia 1 " << endl;
  if (!telSimData.HasPhotonTrace(fevt::FdConstants::eTotal, iwl))
    cout << " problem!" << endl;
  //    telSimData.MakePhotonTrace(fevt::FdConstants::eTotal);
  const TraceD& totalTrace = telSimData.GetPhotonTrace(fevt::FdConstants::eTotal, iwl);
  cout << " graph plot at dia 2 " << endl;

  fLightAtDiaCanvas->cd();
  cout << " graph plot at dia 3 " << endl;

  const int size2 = int(totalTrace.GetSize());

  double x[size2];
  double y[size2];

  for (int i = 0; i < size2; ++i) {
    x[i] = i;
    y[i] = totalTrace[i];
  }

  //  unsigned int size = totalTrace.GetSize();
  /* double tracebin   =  totalTrace.GetBinning();

  TraceD* auxTrace = new TraceD(2000,100);
  for (unsigned int i=0; i<size2; i++) {
    unsigned int bin = int(i*tracebin/100.0);
    (*auxTrace)[bin] =
      photonTrace[i];
  }
  for (int i=0; i<2000; ++i) {
    x[i] = i;
    y[i] = auxTrace[i];
  }
  */

  /*
    int k=0;
    int j=0;
    while (j<1000) {
    j = int(k*tracebin/100.0);
    x[j] = j;
    y[j] += totalTrace[k];
    k++;
  }
  */
  cout << " graph plot at dia 4 " << endl;

  //Detector& theDetector = Detector::GetInstance();
  const double diaphragmArea =
    Detector::GetInstance().GetFDetector().GetEye(tel.GetEyeId()).GetTelescope(tel.GetId()).GetDiaphragmArea();

  cout << " graph plot at dia 5 " << endl;
  TGraph* const graph = new TGraph(size2, x, y);

  graph->SetMarkerStyle(8);
  graph->Draw("APL");
  TString s1;
  s1 = "Light At Diaphragm - area";
  s1 += diaphragmArea;
  s1 += "Eye ";
  s1 += tel.GetEyeId();
  s1 += "Tel ";
  s1 += tel.GetId();

  //  graph->SetTitle("Light At Diaphragm");
  graph->SetTitle(s1.Data());
  //BRD 18/2/05
  //  graph->GetXaxis()->SetTitle("Vertical depth (g/cm²)");
  graph->GetXaxis()->SetTitle("trace bin");
  graph->GetYaxis()->SetTitle("photons/m2");

  TString s;
  fOutputFile->cd();
  s = "LightAtDia_";
  s += *fEventId;
  fLightAtDiaCanvas->SetName(s.Data());
  fLightAtDiaCanvas->Write();
}


void
UserModule::PlotLightOnCamera(const fevt::Telescope& tel)
{
  fLightOnCameraCanvas->cd();
  const int size = 440;

  bool first = true;
  TraceD trace;
  for (int i = 0; i < size; ++i) {
    if (tel.HasPixel(i+1)) {
      const fevt::Pixel& pix = tel.GetPixel(i+1);
      if (pix.HasSimData()) {
        const fevt::PixelSimData& pxsimdata = pix.GetSimData();
        if (pxsimdata.HasPhotonTrace(fevt::FdConstants::eTotal)){
          if (first) {
            trace = pxsimdata.GetPhotonTrace(fevt::FdConstants::eTotal);
            first = false;
          } else {
            trace += pxsimdata.GetPhotonTrace(fevt::FdConstants::eTotal);
          }
        }

      }
    }
  }
  const int size2 = int(trace.GetSize());
  double x[size2];
  double y[size2];
  for (int i = 0; i < size2; ++i) {
    x[i] = i;
    y[i] = 0;
  }

  for (unsigned int i = trace.GetStart(); i < trace.GetStop(); ++i) {
    y[i] = trace[i];
  }

  unsigned int j = 0;
  for (unsigned int i = trace.GetStop()-1; i >= trace.GetStart(); i--) {
    if (y[i]) {
      j = i;
      break;
    }
  }

  TGraph* const gt = new TGraph(j, x, y);
  gt->SetLineColor(1);
  gt->Draw("AL");
  gt->SetTitle("Photons(370nm) on camera");
  gt->GetXaxis()->SetTitle("bin [100 ns]");
  gt->GetXaxis()->SetRange(0,j);
  gt->GetYaxis()->SetTitle("#Photons");

  TString s;
  fOutputFile->cd();
  s = "LightOnCameraC_";
  s += *fEventId;
  s += "Eye ";
  s += tel.GetEyeId();
  s += "Tel ";
  s += tel.GetId();
  fLightOnCameraCanvas->SetName(s.Data());
  fLightOnCameraCanvas->Write();
}


void
UserModule::PlotCamera(const fevt::Telescope& tel)
{
  fCameraCanvas->cd();
  const int size = 440;
  double xs[size];
  double ys[size];
  //  double xfluor[size];
  //  double yfluor[size];

  double integral[size];
  bool signal[size];

  for (int i = 0; i < size; ++i){
    integral[i] = 0;
    signal[i] = false;
  }

  double maxint = 0;
  double minint = 0;
  bool first = true;
  //int sizeSignal = 0;
  //int j = 0;

  for (int i=0; i<size; ++i) {

    if (tel.HasPixel(i+1)) {

      const fevt::Pixel& pix = tel.GetPixel(i+1);

      if (pix.HasSimData()) {
        const fevt::PixelSimData& pxsimdata = pix.GetSimData();


        bool firstTrace=true;
        TraceD trace;

        if (pxsimdata.HasPhotonTrace(fevt::FdConstants::eTotal)) {
          if (firstTrace) {
            firstTrace=false;
            trace = pxsimdata.GetPhotonTrace(fevt::FdConstants::eTotal);
          }
          signal[i] = true;

        } else {

          if (pxsimdata.HasPhotonTrace(fevt::FdConstants::eFluorDirect)){
            if (firstTrace) {
              firstTrace=false;
              trace = pxsimdata.GetPhotonTrace(fevt::FdConstants::eFluorDirect);
            }
            signal[i] = true;
          }

          if (pxsimdata.HasPhotonTrace(fevt::FdConstants::eCherDirect)) {
            if (firstTrace) {
              firstTrace=false;
              trace = pxsimdata.GetPhotonTrace(fevt::FdConstants::eCherDirect);
            } else
              trace += pxsimdata.GetPhotonTrace(fevt::FdConstants::eCherDirect);
            signal[i] = true;
          }

          if (pxsimdata.HasPhotonTrace(fevt::FdConstants::eCherMieScattered)) {
            if (firstTrace) {
              firstTrace=false;
              trace = pxsimdata.GetPhotonTrace(fevt::FdConstants::eCherMieScattered);
            } else
              trace += pxsimdata.GetPhotonTrace(fevt::FdConstants::eCherMieScattered);
            signal[i] = true;
          }

          if (pxsimdata.HasPhotonTrace(fevt::FdConstants::eCherRayleighScattered)){
            if (firstTrace) {
              firstTrace=false;
              trace = pxsimdata.GetPhotonTrace(fevt::FdConstants::eCherRayleighScattered);
            } else
              trace += pxsimdata.GetPhotonTrace(fevt::FdConstants::eCherRayleighScattered);
            signal[i] = true;
          }
        }

        /*
        if (pxsimdata.HasPhotonTrace(fevt::FdConstants::eBackground)) {
          if (firstTrace) {
            firstTrace=false;
            trace = pxsimdata.GetPhotonTrace(fevt::FdConstants::eBackground);
          } else
            trace += pxsimdata.GetPhotonTrace(fevt::FdConstants::eBackground);
        }
        */

        if (!firstTrace) {

          const TraceD::SizeType size = trace.GetSize();
          integral[i] = TraceAlgorithm::Sum(trace, 0, size-1);

          if (first){
            minint = integral[i];
            first = false;
          }
          //cout << i << "->" << integral[i] << endl;
          if (integral[i] > maxint)
            maxint = integral[i];

          if (integral[i] < minint)
            minint = integral[i];
        }

        //cout << "pixel with signal: " << i << endl;
        int nx = int(i/22.0);
        int ny = i%22;

        if (ny%2)
          xs[i] = nx;
        else
          xs[i] = nx + 0.5;

        ys[i] = 22 - ny;
      }
    }
  }

  double deltacolor = (maxint - minint)/50;
  double x[size];
  double y[size];

  for (int i = 0; i < size; ++i) {

    int nx = int(i/22.0);
    int ny = i%22;

    if(ny%2)
      x[i] = 1 + nx;
    else
      x[i] = 1 + nx + 0.5;

    y[i] = 1 + 22 - ny;

  }

  TGraph* const graph = new TGraph(size, x, y);
  graph->SetMarkerStyle(4);
  graph->Draw("AP");
  graph->SetTitle("Pixels with photons");
  graph->GetXaxis()->SetNdivisions(0);
  graph->GetYaxis()->SetNdivisions(0);
  //double x2[1];
  //double y2[1];
  //x2[0] = x[0];
  //y2[0] = y[0];

  //Rainbow colors
  gStyle->SetPalette(1);
  // Create the hexagon points
  float yy[] = { 0.5, 0.25,  -0.25 , -0.5,  -0.25, 0.25, 0.5 };
  float xx[] = { 0, 0.433, 0.433, 0,  -0.433, -0.433, 0 };

  for (int i = 0; i < size; ++i) {
    float yyaux[7];
    float xxaux[7];
    for (int l = 0; l < 7; ++l) {
      yyaux[l] = 1 + yy[l] + ys[i];
      xxaux[l] = 1 + xx[l] + xs[i];
    }

    TPolyLine* const pline = new TPolyLine(7, xxaux, yyaux);
    int cl = 51 + int(ceil((integral[i]-minint)/deltacolor));
    if (cl > 100)
      cl = 100;
    pline->SetFillColor(cl);
    pline->SetLineColor(1);
    pline->SetLineWidth(1);
    if (signal[i]){
      pline->SetLineColor(2);
      pline->SetLineWidth(2);
    }
    pline->Draw("f");
    pline->Draw();
  }

  TString s;
  fOutputFile->cd();
  s = "CameraC_";
  s += *fEventId;
  fCameraCanvas->SetName(s.Data());
  fCameraCanvas->Write();
}


void
UserModule::PlotPhotonTraces(const fevt::Telescope& tel)
{
  for (fevt::Telescope::ConstPixelIterator pixiter = tel.PixelsBegin(ComponentSelector::eHasData);
       pixiter!= tel.PixelsEnd(ComponentSelector::eHasData); ++pixiter) {

    const fevt::Pixel& pixel = *pixiter;
    if (pixel.HasSimData()) {
      const fevt::PixelSimData& pixelsimdata = pixel.GetSimData();
      if (pixelsimdata.HasPhotonTrace(fevt::FdConstants::eTotal)) {
        const TraceD& trace = pixelsimdata.GetPhotonTrace(fevt::FdConstants::eTotal);

        // Copy the trace
        const int size = trace.GetSize();
        double x[size];
        for (int i = 0; i < size; ++i)
          x[i] = i;
        double y[size];

        std::copy(trace.Begin(), trace.End(),y);

        // Graph for the photon trace
        TGraph* const graph = new TGraph(500,x,y);
        ostringstream title;
        title << "Tel " << pixel.GetTelescopeId()
              << " Pix "<< pixel.GetId();

        graph->SetName(title.str().c_str());
        graph->SetTitle(title.str().c_str());

        fPhotonGraphs->Add(graph);
      }
    }
  }

  fOutputFile->cd();

  TFolder folder;

  string s = "Pixels_";
  s += *fEventId;
  folder.SetName(s.c_str());
  folder.Add(fPhotonGraphs);
  folder.AddFolder("Total","Total",fPhotonGraphs);

  folder.Write();
}


void
UserModule::PlotFADCTraces(const fevt::Telescope& tel)
{
  /*double x[3000];
  for (int i=0; i<3000; ++i)
    x[i]=i;*/

  for (fevt::Telescope::ConstChannelIterator chiter = tel.ChannelsBegin();
       chiter != tel.ChannelsEnd(); ++chiter){

    const fevt::Channel& channel = *chiter;
    if (channel.HasSimData()) {
      //cout << "Channel " << channel.GetId() << " has sim data" << endl;
      const fevt::ChannelSimData& channelsimdata = channel.GetSimData();
      if (channelsimdata.HasFADCTrace(fevt::FdConstants::eTotal)) {
        //cout << "Channel " << channel.GetId() << " has FADC trace" << endl;
        const TraceI& trace = channelsimdata.GetFADCTrace(fevt::FdConstants::eTotal);
        //cout << "Channel " << channel.GetId() << " has data" << endl;
        // Copy the trace
        const int size = trace.GetSize();
        double x[size];
        for (int i=0; i<size; ++i)
          x[i]=i;
        double y[size];

        std::copy(trace.Begin(), trace.End(),y);

        // Graph for the photon trace
        TGraph* graph = new TGraph(500,x,y);
        ostringstream title;
        title << "Tel "<<  channel.GetTelescopeId()
              << " Channel "<< channel.GetId();

        graph->SetName(title.str().c_str());
        graph->SetTitle(title.str().c_str());

        fFADCGraphs->Add(graph);
      }
    }
  }

  fOutputFile->cd();

  TFolder folder;

  string s = "Channels_";
  s += *fEventId;
  folder.SetName(s.c_str());
  folder.Add(fFADCGraphs);
  folder.AddFolder("FADC", "FADC", fFADCGraphs);
  folder.Write();
}