KVDetector.h

/***************************************************************************
                          kvdetector.h  -  description
                             -------------------
    begin                : Thu May 16 2002
    copyright            : (C) 2002 by J.D. Frankland
    email                : frankland@ganil.fr
 
$Id: KVDetector.h,v 1.71 2009/05/22 14:45:40 ebonnet Exp $
 ***************************************************************************/
 
/***************************************************************************
 *                                                                         *
 *   This program 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 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/
 
#ifndef KVDETECTOR_H
#define KVDETECTOR_H
 
#ifndef KVD_RECPRC_CNXN
#define KVD_RECPRC_CNXN 1
#endif
#ifndef KVD_NORECPRC_CNXN
#define KVD_NORECPRC_CNXN 0
#endif
 
#include "KVMaterial.h"
#include "KVPosition.h"
#include "KVList.h"
#include "KVUnownedList.h"
#include "KVNucleus.h"
#include "KVGeoDetectorNode.h"
#include "KVUniqueNameList.h"
#include "KVDetectorSignal.h"
#include "KVEvent.h"
 
class KVGeoStrucElement;
class KVGroup;
class KVCalibrator;
class TGeoVolume;
class TTree;
class TGraph;
class KVDetectorSignalTrace;
 
class KVDetector: public KVMaterial, public KVPosition {
 
private:
    KVPosition fEWPosition;
    KVUniqueNameList fParentStrucList;
    KVGeoDetectorNode fNode;
   static Int_t fDetCounter;
    KVMaterial* fActiveLayer;
    TString fNameOfArray;
 
   enum {
      kIsAnalysed = BIT(14),    //for reconstruction of particles
      kActiveSet = BIT(15),     //internal - flag set true when SetActiveLayer called
      kUnidentifiedParticle = BIT(16),  //set if detector is in an unidentified particle's list
      kIdentifiedParticle = BIT(17),    //set if detector is in an identified particle's list
   };
 
    Int_t fIdentP;
    Int_t fUnidentP;
 
   void SetMatrix(const TGeoHMatrix* m) override
   {
      KVPosition::SetMatrix(m);
   }
   void SetShape(TGeoBBox* s) override
   {
      KVPosition::SetShape(s);
   }
   TGeoHMatrix* GetMatrix() const override
   {
      return KVPosition::GetMatrix();
   }
   TGeoBBox* GetShape() const override
   {
      return KVPosition::GetShape();
   }
   TVector3 GetRandomPointOnSurface() const override
   {
      return KVPosition::GetRandomPointOnSurface();
   }
   TVector3 GetSurfaceCentre() const override
   {
      return KVPosition::GetSurfaceCentre();
   }
   TVector3 GetVolumeCentre() const override
   {
      return KVPosition::GetVolumeCentre();
   }
   TVector3 GetSurfaceNormal() const override
   {
      return KVPosition::GetSurfaceNormal();
   }
   Double_t GetSurfaceArea(int npoints = 100000) const override
   {
      return KVPosition::GetSurfaceArea(npoints);
   }
   Double_t GetMisalignmentAngle() const override
   {
      return KVPosition::GetMisalignmentAngle();
   }
 
    KVUniqueNameList fDetSignals;
    mutable KVUniqueNameList fDetSignalsForRawTree;
 
   void remove_signal_for_calibrator(KVCalibrator* K);
 
   template<typename AbsorberStack>
   Double_t get_corrected_energy(AbsorberStack* stack, KVNucleus* nuc, Double_t e, Bool_t transmission)
   {
      auto z = nuc->GetZ();
      auto a = nuc->GetA();
      enum SolType solution = kEmax;
      if (!transmission) solution = kEmin;
      Double_t EINC, ERES = GetEResAfterDetector();
      if (transmission && ERES > 0.) { // if residual energy is known we use it to calculate EINC.
 
         EINC = stack->GetIncidentEnergyFromERes(z, a, ERES);
         if (EINC < stack->GetEIncOfMaxDeltaE(z, a)) {// if EINC < max of dE curve, we change solution
            solution = kEmin;
         } // we could keep the EINC value calculated using ERES, but then the corrected dE of this detector would not depend on the measured dE !
      }
      EINC = stack->GetIncidentEnergy(z, a, e, solution);
      if (EINC < 0) {
         SetEResAfterDetector(-1.);
         return EINC;
      }
      ERES = stack->GetERes(z, a, EINC);
 
      SetEResAfterDetector(-1.);
      return (EINC - ERES);
   }
 
protected:
 
    TString fFName;
    KVList* fCalibrators;
    KVUnownedList fParticles;
    KVList fAbsorbers;
 
   Double_t ELossActive(Double_t* x, Double_t* par);
   Double_t EResDet(Double_t* x, Double_t* par);
   Double_t RangeDet(Double_t* x, Double_t* par);
 
    TF1* fELossF;
    TF1* fEResF;
    TF1* fRangeF;
 
    Double_t fEResforEinc;
 
    Bool_t fSimMode;
    Bool_t fPresent;
    Bool_t fDetecting;
 
    Bool_t fSingleLayer;
 
   void AddDetectorSignal(KVDetectorSignal* ds)
   {
      ds->SetDetector(this);
      fDetSignals.Add(ds);
   }
   virtual Bool_t use_signal_for_raw_data_tree(const TString&) const
   {
      return true;
   }
 
public:
   KVDetector();
   KVDetector(const Char_t* type, const Float_t thick = 0.0);
   KVDetector(const Char_t* gas, const Double_t thick, const Double_t pressure, const Double_t temperature = 19.0);
   KVDetector(const KVDetector&);
   void init();
   virtual ~ KVDetector();
 
   void Copy(TObject& obj) const override;
 
   void SetMaterial(const Char_t* type) override;
   void AddAbsorber(KVMaterial*);
   void SetActiveLayer(KVMaterial* actif)
   {
      fActiveLayer = actif;
      SetBit(kActiveSet);
   }
   void SetActiveLayer(Int_t i)
   {
      if (auto mat = GetAbsorber(i)) SetActiveLayer(mat);
   }
   KVMaterial* GetActiveLayer() const override
   {
      return fActiveLayer;
   }
   KVMaterial* GetAbsorber(Int_t i) const;
   KVMaterial* GetAbsorber(const Char_t* name) const
   {
      return (KVMaterial*)fAbsorbers.FindObject(name);
   }
   const KVList* GetListOfAbsorbers() const
   {
      return &fAbsorbers;
   }
   Int_t GetNumberOfAbsorberLayers() const
   {
      return fAbsorbers.GetEntries();
   }
   void RemoveAllAbsorbers();
 
   Double_t GetTotalThicknessInCM() const
   {
 
      double fTotThickness = 0;
      TIter next(&fAbsorbers);
      KVMaterial* mat;
      while ((mat = (KVMaterial*)next())) fTotThickness += mat->GetThickness();
      return fTotThickness;
   }
   KVGeoDetectorNode* GetNode()
   {
      return &fNode;
   }
 
   const Char_t* GetMaterialName() const
   {
      if (GetActiveLayer())
         return GetActiveLayer()->GetName();
      return KVMaterial::GetName();
   }
   void DetectParticle(KVNucleus*, TVector3* norm = 0) override;
   Double_t GetELostByParticle(KVNucleus*, TVector3* norm = 0) override;
   Double_t GetParticleEIncFromERes(KVNucleus*, TVector3* norm = 0) override;
 
   virtual Double_t GetCalibratedEnergy() const
   {
      return GetDetectorSignalValue("Energy");
   }
   virtual Double_t GetEnergy() const
   {
      Double_t ELoss = GetActiveLayer() ? GetActiveLayer()->GetEnergyLoss() : KVMaterial::GetEnergyLoss();
      if (IsSimMode()) return ELoss; // in simulation mode, return calculated energy loss in active layer
      if (ELoss > 0) return ELoss;
      ELoss = GetCalibratedEnergy();
      if (ELoss < 0) ELoss = 0;
      SetEnergy(ELoss);
      return ELoss;
   }
   virtual void SetEnergy(Double_t e) const
   {
      if (GetActiveLayer()) GetActiveLayer()->SetEnergyLoss(e);
      else KVMaterial::SetEnergyLoss(e);
   }
   Double_t GetEnergyLoss() const override
   {
      return GetEnergy();
   }
   void SetEnergyLoss(Double_t e) const override
   {
      SetEnergy(e);
   }
   virtual Double_t GetCorrectedEnergy(KVNucleus*, Double_t e =
                                          -1., Bool_t transmission = kTRUE);
   virtual Int_t FindZmin(Double_t ELOSS = -1., Char_t mass_formula = -1);
 
   Bool_t AddCalibrator(KVCalibrator* cal, const KVNameValueList& opts = "");
   Bool_t ReplaceCalibrator(const Char_t* type, KVCalibrator* cal, const KVNameValueList& opts = "");
   KVCalibrator* GetCalibrator(const Char_t* name,
                               const Char_t* type) const;
   KVCalibrator* GetCalibrator(const Char_t* type) const;
   KVList* GetListOfCalibrators() const
   {
      return fCalibrators;
   }
   Bool_t IsCalibrated(const KVNameValueList& params = {}) const;
 
   void Clear(Option_t* opt = "") override;
   virtual void Reset(Option_t* opt = "")
   {
      Clear(opt);
   }
   void Print(Option_t* option = "") const override;
 
   void AddHit(KVNucleus* part)
   {
 
      fParticles.Add(part);
   }
 
   auto RemoveHit(KVNucleus* part)
   {
 
      auto p = fParticles.Remove(part);
      if (!fParticles.GetEntries()) SetAnalysed(kFALSE);
      return p;
   }
 
   void ClearHits()
   {
      fParticles.Clear();
   }
   Int_t GetNHits() const
   {
      return fParticles.GetEntries();
   }
 
   Bool_t IsAnalysed()
   {
      return TestBit(kIsAnalysed);
   }
   void SetAnalysed(Bool_t b = kTRUE)
   {
      SetBit(kIsAnalysed, b);
   }
   virtual Bool_t Fired(Option_t* opt = "any") const
   {
 
      if (!IsDetecting()) return kFALSE; //detector not working, no answer at all
      if (IsSimMode()) return (GetActiveLayer() ? GetActiveLayer()->GetEnergyLoss() > 0. : KVMaterial::GetEnergyLoss() > 0.); // simulation mode: detector fired if energy lost in active layer
 
      TString OPT(opt);
      OPT.ToLower();
      Bool_t all = (OPT == "all");
 
      TIter raw_it(&GetListOfDetectorSignals());
      KVDetectorSignal* ds;
      int count_raw = 0;
      while ((ds = (KVDetectorSignal*)raw_it())) {
         if (ds->IsRaw()) {
            ++count_raw;
            if (ds->IsFired()) {
               if (!all) return kTRUE;
            }
            else {
               if (all) return kFALSE;
            }
         }
      }
      return all && count_raw;
   }
   virtual void RemoveCalibrators();
 
   Double_t GetDetectorSignalValue(const KVString& type, const KVNameValueList& params = "") const
   {
 
      KVDetectorSignal* s = GetDetectorSignal(type);
      return (s ? s->GetValue(params) : 0);
   }
   void SetDetectorSignalValue(const KVString& type, Double_t val) const
   {
 
      KVDetectorSignal* s = GetDetectorSignal(type);
      if (s) s->SetValue(val);
   }
   Double_t GetInverseDetectorSignalValue(const KVString& output, Double_t value, const KVString& input, const KVNameValueList& params = "") const
   {
 
      KVDetectorSignal* s = GetDetectorSignal(output);
      return (s ? s->GetInverseValue(value, input, params) : 0);
   }
   virtual KVDetectorSignal* GetDetectorSignal(const KVString& type) const
   {
 
      return fDetSignals.get_object<KVDetectorSignal>(type);
   }
   Bool_t HasDetectorSignal(const KVString& type) const
   {
      return (GetDetectorSignal(type) != nullptr);
   }
 
   inline void IncrementUnidentifiedParticles(Int_t n = 1)
   {
      fUnidentP += n;
      fUnidentP = (fUnidentP > 0) * fUnidentP;
      SetBit(kUnidentifiedParticle, (Bool_t)(fUnidentP > 0));
   }
   inline void IncrementIdentifiedParticles(Int_t n = 1)
   {
      fIdentP += n;
      fIdentP = (fIdentP > 0) * fIdentP;
      SetBit(kIdentifiedParticle, (Bool_t)(fIdentP > 0));
   }
   Bool_t BelongsToUnidentifiedParticle() const
   {
      return TestBit(kUnidentifiedParticle);
   }
   Bool_t BelongsToIdentifiedParticle() const
   {
      return TestBit(kIdentifiedParticle);
   }
 
   static KVDetector* MakeDetector(const Char_t* name, Float_t thick);
 
   virtual Double_t GetEntranceWindowSurfaceArea();
   TVector3 GetActiveLayerSurfaceCentre() const
   {
      return GetSurfaceCentre();
   }
   TVector3 GetActiveLayerVolumeCentre() const
   {
      return KVPosition::GetVolumeCentre();
   }
   TGeoBBox* GetActiveLayerShape() const
   {
      return GetShape();
   }
   TGeoHMatrix* GetActiveLayerMatrix() const
   {
      return GetMatrix();
   }
 
   Double_t GetMaxDeltaE(Int_t Z, Int_t A) override;
   Double_t GetEIncOfMaxDeltaE(Int_t Z, Int_t A) override;
   Double_t GetDeltaE(Int_t Z, Int_t A, Double_t Einc, Double_t = 0.) override;
   virtual Double_t GetTotalDeltaE(Int_t Z, Int_t A, Double_t Einc);
   Double_t GetERes(Int_t Z, Int_t A, Double_t Einc, Double_t = 0.) override;
   Double_t GetIncidentEnergy(Int_t Z, Int_t A, Double_t delta_e =
                                 -1.0, enum SolType type = kEmax) override;
   /*virtual Double_t GetEResFromDeltaE(...)  - DON'T IMPLEMENT, CALLS GETINCIDENTENERGY*/
   Double_t GetDeltaEFromERes(Int_t Z, Int_t A, Double_t Eres) override;
   Double_t GetIncidentEnergyFromERes(Int_t Z, Int_t A, Double_t Eres) override;
   Double_t GetRange(Int_t Z, Int_t A, Double_t Einc) override;
   Double_t GetLinearRange(Int_t Z, Int_t A, Double_t Einc) override;
   Double_t GetPunchThroughEnergy(Int_t Z, Int_t A) override;
   virtual KVDrawable<TGraph> DrawPunchThroughEnergyVsZ(Int_t massform = KVNucleus::kBetaMass);
   virtual KVDrawable<TGraph> DrawPunchThroughEsurAVsZ(Int_t massform = KVNucleus::kBetaMass);
 
   virtual TF1* GetEResFunction(Int_t Z, Int_t A);
   virtual TF1* GetELossFunction(Int_t Z, Int_t A);
   virtual TF1* GetRangeFunction(Int_t Z, Int_t A);
 
   virtual Double_t GetSmallestEmaxValid(Int_t Z, Int_t A) const;
 
   virtual void SetEResAfterDetector(Double_t e)
   {
      fEResforEinc = e;
   }
   virtual Double_t GetEResAfterDetector() const
   {
      return fEResforEinc;
   }
 
   virtual void ReadDefinitionFromFile(const Char_t*);
 
   virtual void SetSimMode(Bool_t on = kTRUE)
   {
      fSimMode = on;
   }
   virtual Bool_t IsSimMode() const
   {
      return fSimMode;
   }
 
   virtual Bool_t IsPresent() const
   {
      return fPresent;
   }
   void SetPresent(Bool_t yes = kTRUE)
   {
      fPresent = yes;
   }
   virtual Bool_t IsDetecting() const
   {
      return fDetecting;
   }
   void SetDetecting(Bool_t yes = kTRUE)
   {
      fDetecting = yes;
   }
 
   virtual Bool_t IsOK() const
   {
      return (fPresent && fDetecting);
   }
 
   KVGroup* GetGroup() const;
   UInt_t GetGroupNumber();
 
   void AddParentStructure(KVGeoStrucElement* elem);
   void RemoveParentStructure(KVGeoStrucElement* elem);
   KVGeoStrucElement* GetParentStructure(const Char_t* type, const Char_t* name = "") const;
 
   void SetActiveLayerMatrix(const TGeoHMatrix*);
   void SetActiveLayerShape(TGeoBBox*);
   void SetEntranceWindowMatrix(const TGeoHMatrix*);
   void SetEntranceWindowShape(TGeoBBox*);
   const KVPosition& GetEntranceWindow() const
   {
      return fEWPosition;
   }
   const TVector3 GetCentreOfEntranceWindow() const
   {
      return GetEntranceWindow().GetSurfaceCentre();
   }
   Double_t GetSolidAngle() const override
   {
      if (ROOTGeo()) return fEWPosition.GetSolidAngle();
      return KVPosition::GetSolidAngle();
   }
   TVector3 GetRandomDirection(Option_t* t = "isotropic") override
   {
      if (ROOTGeo()) return fEWPosition.GetRandomDirection(t);
      return KVPosition::GetRandomDirection(t);
   }
   void GetRandomAngles(Double_t& th, Double_t& ph, Option_t* t = "isotropic") override
   {
      if (ROOTGeo()) fEWPosition.GetRandomAngles(th, ph, t);
      else KVPosition::GetRandomAngles(th, ph, t);
   }
   TVector3 GetDirection() override
   {
      if (ROOTGeo()) return fEWPosition.GetDirection();
      return KVPosition::GetDirection();
   }
   Double_t GetDistance() const override
   {
      if (ROOTGeo()) return fEWPosition.GetDistance();
      return KVPosition::GetDistance();
   }
   Double_t GetTheta() const override
   {
      if (ROOTGeo()) return fEWPosition.GetTheta();
      return KVPosition::GetTheta();
   }
   Double_t GetSinTheta() const override
   {
      if (ROOTGeo()) return fEWPosition.GetSinTheta();
      return KVPosition::GetSinTheta();
   }
   Double_t GetCosTheta() const override
   {
      if (ROOTGeo()) return fEWPosition.GetCosTheta();
      return KVPosition::GetCosTheta();
   }
   Double_t GetPhi() const override
   {
      if (ROOTGeo()) return fEWPosition.GetPhi();
      return KVPosition::GetPhi();
   }
 
   void SetThickness(Double_t thick) override;
   Bool_t IsSingleLayer() const
   {
      return fSingleLayer;
   }
   Bool_t IsMultiLayer() const
   {
      return !IsSingleLayer();
   }
   Bool_t IsGasDetector() const
   {
      return IsGas();
   }
   Bool_t HasSameStructureAs(const KVDetector*) const;
   void SetNameOfArray(const TString& n)
   {
      fNameOfArray = n;
   }
   const Char_t* GetNameOfArray() const
   {
      return fNameOfArray;
   }
 
   const KVSeqCollection& GetListOfDetectorSignals() const
   {
      return fDetSignals;
   }
   const KVSeqCollection& GetListOfDetectorSignalsForRawTree() const
   {
 
      if(fDetSignalsForRawTree.IsEmpty())
      {
         for(auto s : fDetSignals)
         {
            auto sig = dynamic_cast<KVDetectorSignal*>(s);
            if(use_signal_for_raw_data_tree(sig->GetName()))
               fDetSignalsForRawTree.Add(sig);
         }
      }
      return fDetSignalsForRawTree;
   }
   KVDetectorSignal* AddDetectorSignal(const KVString& type);
   KVDetectorSignalTrace* AddDetectorSignalTrace();
   Bool_t AddDetectorSignalExpression(const KVString& type, const KVString& _expr);
 
   virtual Int_t GetIndex() const
   {
 
      return 0;
   }
 
   void SetPressure(Double_t P) override;
   void SetTemperature(Double_t T) override;
 
   virtual Bool_t IsSegmented() const
   {
      return kFALSE;
   }
   virtual void AddEnergyLossInSubDetector(int, double) {}
 
   ClassDefOverride(KVDetector, 10)      //Base class for the description of detectors in multidetector arrays
};
 
inline KVCalibrator* KVDetector::GetCalibrator(const Char_t* name,
      const Char_t* type) const
{
   if (fCalibrators)
      return (KVCalibrator*) fCalibrators->FindObjectWithNameAndType(name, type);
   return 0;
}
 
inline KVCalibrator* KVDetector::GetCalibrator(const Char_t* type) const
{
   if (fCalibrators)
      return (KVCalibrator*) fCalibrators->FindObjectByType(type);
   return 0;
}
#endif