// $Id$
//
//    File: DEventProcessor_photons.cc
// Created: Tue May 24 10:30:32 EDT 2011
// Creator: zihlmann (on Linux chaos 2.6.35-28-generic unknown)
//

#include "DEventProcessor_photons.h"
using namespace jana;

#include "TRACKING/DMCThrown.h"
#include "BCAL/DBCALShower.h"
#include "BCAL/DBCALCluster.h"
#include "FCAL/DFCALShower.h"
#include "FCAL/DFCALCluster.h"



// Routine used to create our DEventProcessor
#include <JANA/JApplication.h>
extern "C"{
void InitPlugin(JApplication *app){
	InitJANAPlugin(app);
	app->AddProcessor(new DEventProcessor_photons());
}
} // "C"


//------------------
// DEventProcessor_photons (Constructor)
//------------------
DEventProcessor_photons::DEventProcessor_photons()
{

}

//------------------
// ~DEventProcessor_photons (Destructor)
//------------------
DEventProcessor_photons::~DEventProcessor_photons()
{

}

//------------------
// init
//------------------
jerror_t DEventProcessor_photons::init(void)
{
	// Create histograms here
 
  ROOTFile = new TFile("photonhists.root", "RECREATE");
  ROOTFile->cd();

  hist2d[0] = new TH2F("hist2d_0","Energy vs Angle gen", 30, 6.5, 21.5, 30, 0.1, 3.1);
  hist2d[0]->GetXaxis()->SetTitle("Angle generated [deg]");
  hist2d[0]->GetYaxis()->SetTitle("Energy generated [GeV]");
  
  hist2d[1] = new TH2F("hist2d_1","Energy vs Angle rec", 30, 6.5, 21.5, 30, 0.1, 3.1);
  hist2d[1]->GetXaxis()->SetTitle("Angle reconstructed [deg]");
  hist2d[1]->GetYaxis()->SetTitle("Energy reconstructed [GeV]");

  hist2d[2] = new TH2F("hist2d_2","dE/E vs Angle rec", 30, 6.5, 21.5, 50, -1., 1.);
  hist2d[2]->GetXaxis()->SetTitle("Angle reconstructed [deg]");
  hist2d[2]->GetYaxis()->SetTitle("(Egen-Erec)/Egen");

  hist2d[3] = new TH2F("hist2d_3","dE/E vs Energy rec", 30, 0.1, 3.1, 50, -1., 1.);
  hist2d[3]->GetXaxis()->SetTitle("Energy reconstructed [GeV]");
  hist2d[3]->GetYaxis()->SetTitle("(Egen-Erec)/Egen");

  hist2d[4] =  new TH2F("hist2d_4","Generated Energy vs Reconstr. Energy BCAL", 
			64, 0., 3.2, 64, 0., 3.2);
  hist2d[4]->GetYaxis()->SetTitle("Energy generated [GeV]");
  hist2d[4]->GetXaxis()->SetTitle("Energy reconstr. [GeV]");

  hist2d[5] =  new TH2F("hist2d_5","Generated Energy vs Reconstr. Energy FCAL", 
			64, 0., 3.2, 64, 0., 3.2 );
  hist2d[5]->GetYaxis()->SetTitle("Energy generated [GeV]");
  hist2d[5]->GetXaxis()->SetTitle("Energy reconstr. [GeV]");

  hist2d[6] =  new TH2F("hist2d_6","Generated Angle vs Reconstr. Angle BCAL", 
			240, 5., 25., 240, 5., 25.);
  hist2d[6]->GetYaxis()->SetTitle("Angle generated [deg.]");
  hist2d[6]->GetXaxis()->SetTitle("Angle reconstr. [deg.]");

  hist2d[7] =  new TH2F("hist2d_7","Generated Angle vs Reconstr. Angle FCAL", 
			240, 5., 25., 240, 5., 25.);
  hist2d[7]->GetYaxis()->SetTitle("Angle generated [GeV]");
  hist2d[7]->GetXaxis()->SetTitle("Angle reconstr. [GeV]");

  hist2d[8] = new TH2F("hist2d_8","dE/E vs Angle gen", 30, 6.5, 21.5, 50, -1., 1.);
  hist2d[8]->GetXaxis()->SetTitle("Angle generated [deg]");
  hist2d[8]->GetYaxis()->SetTitle("(Egen-Erec)/Egen");

  hist2d[9] = new TH2F("hist2d_9","dE/E vs Energy gen", 30, 0.1, 3.1, 50, -1., 1.);
  hist2d[9]->GetXaxis()->SetTitle("Energy generated [GeV]");
  hist2d[9]->GetYaxis()->SetTitle("(Egen-Erec)/Egen");

  EventTree = new TTree("EventTree", "Raw data");

  EventTree->Branch("EventNum",    &EventNum,    "EventNum/I"); 
  EventTree->Branch("NGen",    &NGen,    "NGen/I"); 
  EventTree->Branch("GenE",    GenE,    "GenE[NGen]/F"); 
  EventTree->Branch("GenA",    GenA,    "GenA[NGen]/F"); 

  EventTree->Branch("NRecB",    &NRecB,    "NRecB/I"); 
  EventTree->Branch("BcalE",    BcalE,    "BcalE[NRecB]/F"); 
  EventTree->Branch("BcalA",    BcalA,    "BcalA[NRecB]/F"); 

  EventTree->Branch("NRecF",    &NRecF,    "NRecF/I"); 
  EventTree->Branch("FcalE",    FcalE,    "FcalE[NRecF]/F"); 
  EventTree->Branch("FcalA",    FcalA,    "FcalA[NRecF]/F"); 

  pthread_mutex_init(&mutex, NULL);

  return NOERROR;
}

//------------------
// brun
//------------------
jerror_t DEventProcessor_photons::brun(JEventLoop *eventLoop, int runnumber)
{
	return NOERROR;
}

//------------------
// evnt
//------------------
jerror_t DEventProcessor_photons::evnt(JEventLoop *loop, int eventnumber)
{
	// Fill histograms here

  vector<const DMCThrown*> ThrownPhotons;
  loop->Get(ThrownPhotons);

  vector <const DBCALShower*> BcalPhotons1;
  loop->Get(BcalPhotons1);

  vector <const DFCALShower*> FcalPhotons;
  loop->Get(FcalPhotons);

  Float_t Egen = ThrownPhotons[0]->energy();
  Float_t Agen = ThrownPhotons[0]->momentum().Theta()*180.0/M_PI;

  pthread_mutex_lock(&mutex);

  GenE[0] = Egen;
  GenA[0] = Agen;
  NGen = 1;

  double zVertex = ThrownPhotons[0]->z(); 

  hist2d[0]->Fill(Agen, Egen);

  EventNum = eventnumber;

  Int_t nBcal = BcalPhotons1.size();

  Int_t nFcal = FcalPhotons.size();

  Float_t Erec = 0;
  Float_t Arec = 0;

  NRecB = nBcal; 
  for (Int_t j=0;j<nBcal;j++){

    Double_t l = sqrt(BcalPhotons1[j]->x*BcalPhotons1[j]->x + 
		      BcalPhotons1[j]->y*BcalPhotons1[j]->y);
    
    Erec = BcalPhotons1[j]->E;
    Double_t d = BcalPhotons1[j]->z-zVertex;
    if (d!=0.0){
      Arec =  atan(l/d)*180.0/M_PI;
    } else {
      Arec = 90.;
    }
    
    BcalE[j] = Erec;
    BcalA[j] = Arec;
    
    hist2d[4]->Fill(Erec, Egen);
    hist2d[6]->Fill(Arec, Agen);
  }
  
  NRecF = nFcal;
  for (Int_t j=0;j<nFcal;j++){
    
    Erec = FcalPhotons[j]->getEnergy();
    Double_t l = sqrt(FcalPhotons[j]->getPosition().X()*FcalPhotons[j]->getPosition().X() + 
		      FcalPhotons[j]->getPosition().Y()*FcalPhotons[j]->getPosition().Y());

    Double_t d = FcalPhotons[j]->getPosition().Z()-65.;
    if (d!=0.0){
      Arec =  atan(l/d)*180.0/M_PI;
    } else {
      Arec = 90.;
    }

    FcalE[j] = Erec;
    FcalA[j] = Arec;
    hist2d[5]->Fill(Erec, Egen);
    hist2d[7]->Fill(Arec, Agen);
  }
  
  if ( ((nBcal==0)&&(nFcal==1)) ||
       ((nBcal==1)&&(nFcal==0))) {
    hist2d[1]->Fill(Arec, Erec);
    Float_t dE = (Egen-Erec)/Egen;
    hist2d[2]->Fill(Arec, dE);
    hist2d[3]->Fill(Erec, dE);
    hist2d[8]->Fill(Agen, dE);
    hist2d[9]->Fill(Egen, dE);
  }
  
  EventTree->Fill();

  pthread_mutex_unlock(&mutex);

  return NOERROR;
}

//------------------
// erun
//------------------
jerror_t DEventProcessor_photons::erun(void)
{
	// Any final calculations on histograms (like dividing them)
	// should be done here. This may get called more than once.
	return NOERROR;
}

//------------------
// fini
//------------------
jerror_t DEventProcessor_photons::fini(void)
{
	// Called at very end. This will be called only once

 
  ROOTFile->cd();
  for (Int_t i=0;i<10;i++)
    hist2d[i]->Write();

  EventTree->Write();

  ROOTFile->Close();
   
  cout<<endl<<"Close ROOT file"<<endl;
  
  return NOERROR;
}