EXPERT ROOT

// -------------------------------------------------------------------------

// -----                        ERND source file                   -----

// -----                  Created data  by developerName               -----

// -------------------------------------------------------------------------

#include "ERND.h"

#include "TClonesArray.h"

#include "TParticle.h"

#include "TVirtualMC.h"

#include "TString.h"

#include "FairRootManager.h"

#include "FairLogger.h"

// -----   Default constructor   -------------------------------------------

ERND::ERND() : 

  FairDetector("ERND", kTRUE),

   fNDPoints(NULL),

   fStep(0)

{

  LOG(INFO) << "  ERND::ERND()" << FairLogger::endl;

  ResetParameters();

  fNDPoints = new TClonesArray("ERNDPoint");

  flGeoPar = new TList();

  flGeoPar->SetName( GetName());

  fVerboseLevel = 1;

  fVersion = 1;

}

// -------------------------------------------------------------------------

// -----   Standard constructor   ------------------------------------------

ERND::ERND(const char* name, Bool_t active, Int_t verbose) 

  : FairDetector(name, active,verbose),

   fNDPoints(NULL)

  {

  ResetParameters();

  fNDPoints = new TClonesArray("ERNDPoint");

  flGeoPar = new TList();

  flGeoPar->SetName( GetName());

  fVersion = 1;

}

ERND::~ERND() {

  if (fNDPoints) {

    fNDPoints->Delete();

    delete fNDPoints;

  }

}

void ERND::Initialize()

{

  FairDetector::Initialize();

}

Bool_t ERND::ProcessHits(FairVolume* vol) {

  // Set constants for Birk's Law implentation

  /*

  static const Double_t dP = 1.032 ;

  static const Double_t BirkC1 =  0.013/dP;

  static const Double_t BirkC2 =  9.6e-6/(dP * dP);

  */

  //Birks constants from Craun, R. L.; Smith, D. L. NIM 80,2, p. 239, 1970

  static const Double_t dP = 0.97;

 // static const Double_t BirkC1 =  0.00856/dP;

  static const Double_t BirkC1 =  0.001/dP;

  static const Double_t BirkC2 =  4.99e-6/(dP * dP);

  static Int_t          eventID;           //!  event index

  static Int_t          trackID;           //!  track index

  static Int_t          mot0TrackID;       //!  0th mother track index

  static Int_t       pdg;              //!  pdg

  static TLorentzVector posIn, posOut;    //!  position

  static TLorentzVector momIn, momOut;    //!  momentum

  static Double32_t     time;              //!  time

  static Double32_t     length;            //!  length

  static Double32_t     eLoss;             //!  energy loss

  static Int_t          stilbenNr;

  static Double_t       lightYield;

  gMC->SetMaxStep(0.01);

  if ( gMC->IsTrackEntering() ) { // Return true if this is the first step of the track in the current volume

    eLoss  = 0.;

    lightYield = 0.;

    eventID = gMC->CurrentEvent();

    gMC->TrackPosition(posIn);

    gMC->TrackMomentum(momIn);

    trackID  = gMC->GetStack()->GetCurrentTrackNumber();

    time   = gMC->TrackTime() * 1.0e09;  // Return the current time of flight of the track being transported

    length = gMC->TrackLength(); // Return the length of the current track from its origin (in cm)

    mot0TrackID  = gMC->GetStack()->GetCurrentTrack()->GetMother(0);

    pdg = gMC->TrackPid(); // GeV/c2

    gMC->CurrentVolOffID(3,stilbenNr); 

  }

  eLoss += gMC->Edep(); // GeV //Return the energy lost in the current step

  Double_t curLightYield = 0.;

  // Apply Birk's law ( Adapted from G3BIRK/Geant3)

  // Correction for all charge states

  if (gMC->TrackCharge()!=0) { // Return the charge of the track currently transported

    Double_t BirkC1Mod = 0;

    // Apply correction for higher charge states

      if (TMath::Abs(gMC->TrackCharge())>=2)

        BirkC1Mod=BirkC1*7.2/12.6;

      else

        BirkC1Mod=BirkC1;

    if (gMC->TrackStep()>0)

    {

      Double_t dedxcm=gMC->Edep()*1000./gMC->TrackStep(); //[MeV/cm]

      curLightYield=gMC->Edep()*1000./(1.+BirkC1Mod*dedxcm+BirkC2*dedxcm*dedxcm); //[MeV]

      curLightYield /= 1000.; //[GeV]

      lightYield+=curLightYield;

    }

  }

        if (gMC->IsTrackExiting()    || //Return true if this is the last step of the track in the current volume 

            gMC->IsTrackStop()       || //Return true if the track energy has fallen below the threshold

            gMC->IsTrackDisappeared()) 

        { 

    gMC->TrackPosition(posOut);

    gMC->TrackMomentum(momOut);

          if (eLoss > 0. && gMC->TrackCharge()!=0){

      AddPoint( eventID, trackID, mot0TrackID, pdg,

                TVector3(posIn.X(),   posIn.Y(),   posIn.Z()),

                TVector3(posOut.X(),  posOut.Y(),  posOut.Z()),

                TVector3(momIn.Px(),  momIn.Py(),  momIn.Pz()),

                TVector3(momOut.Px(), momOut.Py(), momOut.Pz()),

                time, length, eLoss, stilbenNr, lightYield);

    }

  }

  return kTRUE;

}

// -----   Public method EndOfEvent   -----------------------------------------

void ERND::BeginEvent() {

}

void ERND::EndOfEvent() {

  LOG(DEBUG) << "ND Points Count: " << fNDPoints->GetEntriesFast() << FairLogger::endl;

  Print();

  Reset();

}

// -----   Public method Register   -------------------------------------------

void ERND::Register() {

  FairRootManager* ioman = FairRootManager::Instance();

  if (!ioman)

        Fatal("Init", "IO manager is not set");        

  ioman->Register("NDPoint","ND", fNDPoints, kTRUE);

}

// ----------------------------------------------------------------------------

// -----   Public method GetCollection   --------------------------------------

TClonesArray* ERND::GetCollection(Int_t iColl) const {

  if (iColl == 0) 

    return fNDPoints;

  else 

    return NULL;

}

// ----------------------------------------------------------------------------

// -----   Public method Print   ----------------------------------------------

void ERND::Print(Option_t *option) const

{

  for (Int_t i_point = 0; i_point < fNDPoints->GetEntriesFast(); i_point++){

    ERNDPoint* point = (ERNDPoint*)fNDPoints->At(i_point);

    point->Print();

  }

}

// ----------------------------------------------------------------------------

// -----   Public method Reset   ----------------------------------------------

void ERND::Reset() {

  fNDPoints->Clear();

  ResetParameters();

}

// ----------------------------------------------------------------------------

// -----   Public method CopyClones   -----------------------------------------

void ERND::CopyClones(TClonesArray* cl1, TClonesArray* cl2, Int_t offset) {

  LOG(DEBUG) << "   ERND::CopyClones(TClonesArray* cl1, TClonesArray* cl2, Int_t offset)" 

            << FairLogger::endl;

  Int_t nEntries = cl1->GetEntriesFast();

  LOG(DEBUG) << "decector: " << nEntries << " entries to add" << FairLogger::endl;

  TClonesArray& clref = *cl2;

  ERNDPoint* oldpoint = NULL;

  for (Int_t i=0; i<nEntries; i++) {

  oldpoint = (ERNDPoint*) cl1->At(i);

   Int_t index = oldpoint->GetTrackID() + offset;

   oldpoint->SetTrackID(index);

   new (clref[cl2->GetEntriesFast()]) ERNDPoint(*oldpoint);

  }

  LOG(DEBUG) << "decector: " << cl2->GetEntriesFast() << " merged entries" << FairLogger::endl;

}

// ----------------------------------------------------------------------------

// -----   Private method AddPoint   --------------------------------------------

ERNDPoint* ERND::AddPoint(Int_t eventID, Int_t trackID,

                                    Int_t mot0trackID,

                                    Int_t pdg,

                                    TVector3 posIn,

                                    TVector3 posOut, TVector3 momIn,

                                    TVector3 momOut, Double_t time,

                                    Double_t length, Double_t eLoss, Int_t stilbenNr, Float_t lightYield) {

  TClonesArray& clref = *fNDPoints;

  Int_t size = clref.GetEntriesFast();

  return new(clref[size]) ERNDPoint(eventID, trackID, mot0trackID, pdg,

                                          posIn, posOut, momIn, momOut, time, length, eLoss,stilbenNr,lightYield);

}

// ----------------------------------------------------------------------------

// -----   Public method ConstructGeometry   ----------------------------------

void ERND::ConstructGeometry() {

  TString fileName = GetGeometryFileName();

  if(fileName.EndsWith(".root")) {

    LOG(DEBUG) << "Constructing ERND geometry from ROOT file " << fileName.Data() << FairLogger::endl;

    ConstructRootGeometry();

  } else {

    LOG(DEBUG) << "Geometry file name is not set" << FairLogger::endl;

  }

}

// ----------------------------------------------------------------------------

// ----------------------------------------------------------------------------

Bool_t ERND::CheckIfSensitive(std::string name)

{

  TString volName = name;

  if(volName.Contains("crystal")) {

    return kTRUE;

  }

  return kFALSE;

}

// ----------------------------------------------------------------------------

// ----------------------------------------------------------------------------

void ERND::ResetParameters() {

};

// ----------------------------------------------------------------------------

ClassImp(ERND)