/* This file is part of MAUS: http://micewww.pp.rl.ac.uk:8080/projects/maus

 *

 * MAUS is free software: you can redistribute it and/or modify

 * it under the terms of the GNU General Public License as published by

 * the Free Software Foundation, either version 3 of the License, or

 * (at your option) any later version.

 *

 * MAUS is distributed in the hope that it will be useful,

 * but WITHOUT ANY WARRANTY; without even the implied warranty of

 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

 * GNU General Public License for more details.

 *

 * You should have received a copy of the GNU General Public License

 * along with MAUS.  If not, see <http://www.gnu.org/licenses/>.

 *

 */

#include <iostream>

#include <fstream>

#include <cmath>

#include "TFile.h"

#include "TH1F.h"

#include "TH2F.h"

#include "TTree.h"

#include "G4TransportationManager.hh"

#include "G4FieldManager.hh"

#include "G4Field.hh"

#include "G4HCofThisEvent.hh"

#include "G4Step.hh"

#include "G4ThreeVector.hh"

#include "G4SDManager.hh"

#include "G4ios.hh"

#include "Interface/MICEEvent.hh"

#include "Config/MiceModule.hh"

#include "SciFiSD.hh"

SciFiSD::SciFiSD(MiceModule* mod) : MAUSSD(mod) {

}

G4bool SciFiSD::ProcessHits(G4Step* aStep, G4TouchableHistory* ROhist) {

/* edep reads the total energy deposited by an event.  If the total energy is zero

   then this function returns false, no hit has occured.  If there was some energy

   deposited in the tracker then space is made in the vector and the event information

   is stored in the ith position   */

  G4double edep = aStep->GetTotalEnergyDeposit();

  double pid = aStep->GetTrack()->GetDefinition()->GetPDGEncoding();

  if ( edep == 0. ) return false;

  int old_chanNo = legacy_chanNo(aStep);

  //Every time the function is called hit_i is assigned to the final spot

  //in the vector and then the size of the vector is increased by one.

  int hit_i = _hits.size();

  _hits.push_back(Json::Value());

  //All the information for the tracker event is stored in _hits

  G4TouchableHandle theTouchable = aStep->GetPreStepPoint()->GetTouchableHandle();

  G4int fiberNumber = theTouchable->GetCopyNumber();           // get the fiber copy number

  G4ThreeVector Pos = aStep->GetPreStepPoint()->GetPosition(); // true MC position

  G4ThreeVector Mom = aStep->GetTrack()->GetMomentum();        // true momentum

  Json::Value channel_id;

  channel_id["type"]             = "Tracker";

  channel_id["fiber_number"]     = fiberNumber;

  channel_id["tracker_number"] = _module->propertyInt("Tracker");

  channel_id["station_number"] = _module->propertyInt("Station");

  channel_id["plane_number"]   = _module->propertyInt("Plane");

  _hits[hit_i]["channel_id"] = channel_id;

  _hits[hit_i]["track_id"] = aStep->GetTrack()->GetTrackID();

  _hits[hit_i]["energy"]  = aStep->GetTrack()->GetTotalEnergy();

  _hits[hit_i]["charge"]  =  aStep->GetTrack()->GetDefinition()->GetPDGCharge();

  _hits[hit_i]["pid"] = pid;

  _hits[hit_i]["mass"]= aStep->GetTrack()->GetDefinition()->GetPDGMass();

  _hits[hit_i]["time"]= aStep->GetPreStepPoint()->GetGlobalTime();

  _hits[hit_i]["energy_deposited"] = edep;

  _hits[hit_i]["momentum"]["x"]=Mom.x();

  _hits[hit_i]["momentum"]["y"]=Mom.y();

  _hits[hit_i]["momentum"]["z"]=Mom.z();

  _hits[hit_i]["hit_position"]["x"]=Pos.x();

  _hits[hit_i]["hit_position"]["y"]=Pos.y();

  _hits[hit_i]["hit_position"]["z"]=Pos.z();

  /*

  int parent = aStep->GetTrack()->GetParentID();

  int track  = aStep->GetTrack()->GetTrackID();

  double alarm = 0;

  if (_hits[hit_i]["momentum"]["z"].asDouble() < 0.0 ) {

    alarm=-999;

  }

    std::ofstream out1("negative_mom.txt", std::ios::out | std::ios::app);

    out1 << aStep->GetTrack()->GetCurrentStepNumber() << " " << _hits[hit_i]["momentum"]["x"].asDouble() << " " << _hits[hit_i]["momentum"]["y"].asDouble() << " " << _hits[hit_i]["momentum"]["z"].asDouble() << " "

         << _hits[hit_i]["hit_position"]["x"].asDouble() << " " << _hits[hit_i]["hit_position"]["y"].asDouble() << " " << _hits[hit_i]["hit_position"]["z"].asDouble()

         << " " << edep << " " << pid << " " << parent << " " << track << " " <<  alarm << "\n";

    out1.close();

  */  

  //Determines the channel number the fiber is located in.

  //It does not look like this is saved to a vector or used above, is this data used anywhere?

  int chanNo;

  int numbFibers = 7*2*(_module->propertyDouble("CentralFibre")+0.5);

  if ( _module->propertyInt("Tracker") == 0 ) {

    chanNo = floor((numbFibers-fiberNumber)/7.0);

  } else {

    chanNo = floor(fiberNumber/7.0); // -1 because we start numbering from 0

  }

  return true;

}

// Shouldn't this be used to clear out variables before calling the above again?

void SciFiSD::EndOfEvent(G4HCofThisEvent* HCE) {

  // do nothing

}

// This probably should not be here in the final version.

// If we need a check on the above code a unit test should be written and placed there.

int SciFiSD::legacy_chanNo(G4Step* aStep) {

  Hep3Vector pos = _module->globalPosition();

  Hep3Vector perp(-1., 0., 0.);

  double chanWidth = 1.4945; // effective channel width without overlap

  perp *= _module->globalRotation();

  Hep3Vector delta = aStep->GetPreStepPoint()->GetPosition() - pos;

  double dist = delta.x() * perp.x() + delta.y() * perp.y() + delta.z() * perp.z();

  double fibre = _module->propertyDouble("CentralFibre") + dist * 2.0 /

                ( _module->propertyDouble("Pitch") * 7.0 );

  double centFibre = _module->propertyDouble("CentralFibre");

  int  firstChan = (int)(fibre + 0.5);

  int secondChan(-1);

  double overlap = chanWidth - (0.1365/2.0);

  double underlap = (0.1365/2.0);

  // int nChans(1);

  nChans = 1;

  double distInChan = ((firstChan - centFibre)*chanWidth) - dist;

  // distance in channel greater than section which does not overlap

  if ((sqrt(distInChan*distInChan) > overlap) || (sqrt(distInChan*distInChan) < underlap)) {

    nChans = 2;

    if (distInChan < 0) {

      secondChan = firstChan - 1;

    } else {

      secondChan = firstChan + 1;

    }

  }

  return firstChan;

}