#include "GenNoiseElectrons.h"
#include "FairLogger.h"
#include "TClonesArray.h"
#include "CbmMuchPoint.h"
#include "TH1F.h"
#include "TRandom3.h"
#include "CbmMuchGeoScheme.h"
#include "TGeoManager.h"
#include "CbmMCTrack.h"
#include "CbmTrdPoint.h"

ClassImp(LxGenNoiseElectrons)

LxGenNoiseElectrons::LxGenNoiseElectrons() : fNofNoiseE(1), fNofStations(4), fMCTracks(0), fMuchPoints(0), fTrdPoints(0),
   fOutMCTracks(0), fOutMuchPoints(0), fOutTrdPoints(0)
{
}

LxGenNoiseElectrons::~LxGenNoiseElectrons()
{
   delete fOutMCTracks;
   delete fOutMuchPoints;
   delete fOutTrdPoints;
}

static Double_t x_rmss[4][3];
static Double_t y_rmss[4][3];

static Double_t GetRMS(const char* name)
{
   Double_t result = -1.;
   char fileName[128];
   sprintf(fileName, "%s.root", name);
   TFile* curFile = TFile::CurrentFile();
   TFile* f = new TFile(fileName);
   
   if (!f->IsZombie())
   {
      TH1F* h = static_cast<TH1F*> (f->Get(name));
   
      if (0 != h)
         result = h->GetRMS();
   }
   
   delete f;
   TFile::CurrentFile() = curFile;
   return result;
}

InitStatus LxGenNoiseElectrons::Init()
{
   FairRootManager* ioman = FairRootManager::Instance();
    
   if (0 == ioman)
      LOG(fatal) <<  "No FairRootManager";
   
   fMCTracks = static_cast<TClonesArray*> (ioman->GetObject("MCTrack"));
   fMuchPoints = static_cast<TClonesArray*> (ioman->GetObject("MuchPoint"));
   fTrdPoints = static_cast<TClonesArray*> (ioman->GetObject("TrdPoint"));
   
   if (0 == fMCTracks || 0 == fMuchPoints || 0 == fTrdPoints)
      LOG(fatal) <<  "No MC tracks or points";
   
   char name[128];
   
   for (int i = 0; i < 4; ++i)
   {
      for (int j = 0; j < 3; ++j)
      {
         sprintf(name, "noise_e_x_%d_%d", i, j);
         x_rmss[i][j] = GetRMS(name);
         
         if (x_rmss[i][j] < 0)
         {
            sprintf(name, "Couldn't read noise_e_x_%d_%d", i, j);
            LOG(fatal) <<  name;
         }
         
         sprintf(name, "noise_e_y_%d_%d", i, j);
         y_rmss[i][j] = GetRMS(name);
         
         if (y_rmss[i][j] < 0)
         {
            sprintf(name, "Couldn't read noise_e_y_%d_%d", i, j);
            LOG(fatal) <<  name;
         }
      }
   }
   
   time_t initTime;
   gRandom->SetSeed(time(&initTime));
   
   fOutMCTracks = new TClonesArray("CbmMCTrack", 1000);
   ioman->Register("MCTrack", "MC", fOutMCTracks, kTRUE);
   fOutMuchPoints = new TClonesArray("CbmMuchPoint", 1000);
   ioman->Register("MuchPoint", "MUCH", fOutMuchPoints, kTRUE);
   fOutTrdPoints = new TClonesArray("CbmTrdPoint", 1000);
   ioman->Register("TrdPoint", "TRD", fOutTrdPoints, kTRUE);
   
   return kSUCCESS;
}

void LxGenNoiseElectrons::Exec(Option_t* /*opt*/)
{
   fOutMCTracks->Delete();
   fOutMuchPoints->Delete();
   fOutTrdPoints->Delete();
   
   Int_t nofMCTracks = fMCTracks->GetEntriesFast();
   TClonesArray& mctref = *fOutMCTracks;
   
   for (Int_t i = 0; i < nofMCTracks; ++i)
   {
      const CbmMCTrack* mcTrack = static_cast<const CbmMCTrack*> (fMCTracks->At(i));
      new(mctref[i]) CbmMCTrack(*mcTrack);
   }
   
   Int_t nofMuchPoints = fMuchPoints->GetEntriesFast();
   TClonesArray& mpref = *fOutMuchPoints;
   
   for (Int_t i = 0; i < nofMuchPoints; ++i)
   {
      const CbmMuchPoint* muchPoint = static_cast<const CbmMuchPoint*> (fMuchPoints->At(i));
      new(mpref[i]) CbmMuchPoint(*muchPoint);
   }
   
   Int_t nofTrdPoints = fTrdPoints->GetEntriesFast();
   TClonesArray& tpref = *fOutTrdPoints;
   
   for (Int_t i = 0; i < nofTrdPoints; ++i)
   {
      const CbmTrdPoint* trdPoint = static_cast<const CbmTrdPoint*> (fTrdPoints->At(i));
      new(tpref[i]) CbmTrdPoint(*trdPoint);
   }
   
   Int_t ind = nofMuchPoints;
   
   for (Int_t i = 0; i < nofMuchPoints; ++i)
   {
      const CbmMuchPoint* pMuchPt = static_cast<const CbmMuchPoint*> (fMuchPoints->At(i));
      Int_t stationNumber = CbmMuchGeoScheme::GetStationIndex(pMuchPt->GetDetectorId());
      Int_t layerNumber = CbmMuchGeoScheme::GetLayerIndex(pMuchPt->GetDetectorId());
      
      for (int j = 0; j < fNofNoiseE; ++j)
      {
         TVector3 posIn;
         pMuchPt->PositionIn(posIn);
         TVector3 posOut;
         pMuchPt->PositionOut(posOut);
         TVector3 posDelta = posOut - posIn;
         
         for (int k = 0; k < 20; ++k)
         {
            TVector3 noiseDelta(gRandom->Gaus(0, x_rmss[stationNumber][layerNumber]), gRandom->Gaus(0, y_rmss[stationNumber][layerNumber]), 0);
            TVector3 newIn = posIn + noiseDelta;
            TGeoNode* node = gGeoManager->FindNode(newIn.X(), newIn.Y(), newIn.Z());
         
            if (0 == node)
               continue;
         
            TString nodeName(node->GetName());
         
            if (!nodeName.Contains("active", TString::kIgnoreCase))
               continue;
            
            TVector3 newOut = newIn + 0.5 * posDelta;
            node = gGeoManager->FindNode(newOut.X(), newOut.Y(), newOut.Z());
         
            if (0 == node)
               continue;
         
            TString nodeName2(node->GetName());
            
            if (!nodeName2.Contains("active", TString::kIgnoreCase))
               continue;
            
            CbmMuchPoint noiseMuchPt(*pMuchPt);
            noiseMuchPt.SetPosition(newIn);
            noiseMuchPt.SetPositionOut(newOut);
            noiseMuchPt.SetTrackID(-1);
            new(mpref[ind++]) CbmMuchPoint(noiseMuchPt);
            break;
         }
      }
   }
}

void LxGenNoiseElectrons::Finish()
{
}