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 "KVNucleus.h"
#include "Binary_t.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 KVDetector: public KVMaterial, public KVPosition {
 
private:
    KVPosition fEWPosition;
    KVUniqueNameList fParentStrucList;
    KVGeoDetectorNode fNode;
   static Int_t fDetCounter;
    KVMaterial* fActiveLayer;
    KVList* fIDTelescopes;
    KVList* fIDTelAlign;
    TList* fIDTele4Ident;
    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)
   {
      KVPosition::SetMatrix(m);
   }
   void SetShape(TGeoBBox* s)
   {
      KVPosition::SetShape(s);
   }
   TGeoHMatrix* GetMatrix() const
   {
      return KVPosition::GetMatrix();
   }
   TGeoBBox* GetShape() const
   {
      return KVPosition::GetShape();
   }
   TVector3 GetRandomPointOnSurface() const
   {
      return KVPosition::GetRandomPointOnSurface();
   }
   TVector3 GetSurfaceCentre() const
   {
      return KVPosition::GetSurfaceCentre();
   }
   TVector3 GetVolumeCentre() const
   {
      return KVPosition::GetVolumeCentre();
   }
   TVector3 GetSurfaceNormal() const
   {
      return KVPosition::GetSurfaceNormal();
   }
   Double_t GetSurfaceArea(int npoints = 100000) const
   {
      return KVPosition::GetSurfaceArea(npoints);
   }
   Double_t GetMisalignmentAngle() const
   {
      return KVPosition::GetMisalignmentAngle();
   }
 
    KVUniqueNameList fDetSignals;
 
   void remove_signal_for_calibrator(KVCalibrator* K);
 
protected:
 
    TString fFName;
    KVList* fCalibrators;
    KVList* fParticles;
    KVList* fAbsorbers;
    UShort_t fSegment;
    Double_t fGain;
    Int_t fCalWarning;
 
    Double_t fTotThickness;
    Double_t fDepthInTelescope;
 
   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;
    TList* fAlignedDetectors[2];
 
    Bool_t fSimMode;
    Bool_t fPresent;
    Bool_t fDetecting;
 
    Bool_t fSingleLayer;
 
   void AddDetectorSignal(KVDetectorSignal* ds)
   {
      ds->SetDetector(this);
      fDetSignals.Add(ds);
   }
 
public:
   KVDetector();
   KVDetector(const Char_t* type, const Float_t thick = 0.0);
   KVDetector(const KVDetector&);
   void init();
   virtual ~ KVDetector();
 
#if ROOT_VERSION_CODE >= ROOT_VERSION(3,4,0)
   virtual void Copy(TObject& obj) const;
#else
   virtual void Copy(TObject& obj);
#endif
 
   virtual void SetMaterial(const Char_t* type);
   void AddAbsorber(KVMaterial*);
   void SetActiveLayer(KVMaterial* actif)
   {
      fActiveLayer = actif;
      SetBit(kActiveSet);
   }
   KVMaterial* GetActiveLayer() const
   {
      return fActiveLayer;
   }
   KVMaterial* GetAbsorber(Int_t i) const;
   KVMaterial* GetAbsorber(const Char_t* name) const
   {
      return (KVMaterial*)(fAbsorbers ? fAbsorbers->FindObject(name) : 0);
   }
   KVList* GetListOfAbsorbers() const
   {
      return fAbsorbers;
   }
   Int_t GetNumberOfAbsorberLayers() const
   {
      return fAbsorbers->GetEntries();
   }
   virtual const Char_t* GetArrayName();
   virtual Double_t GetDepthInTelescope() const
   {
      return fDepthInTelescope;
   }
 
   Double_t GetTotalThicknessInCM()
   {
 
      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();
   };
   virtual void DetectParticle(KVNucleus*, TVector3* norm = 0);
   virtual Double_t GetELostByParticle(KVNucleus*, TVector3* norm = 0);
   virtual Double_t GetParticleEIncFromERes(KVNucleus*, TVector3* norm = 0);
 
   void SetGain(Double_t gain);
   Double_t GetGain() const;
 
   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);
   }
   virtual Double_t GetEnergyLoss() const
   {
      return GetEnergy();
   }
   virtual void SetEnergyLoss(Double_t e) const
   {
      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
   {
      return (HasDetectorSignal("Energy") && IsOK());
   }
   Bool_t IsCalibrated(const KVNameValueList& params) const;
 
   virtual void Clear(Option_t* opt = "");
   virtual void Reset(Option_t* opt = "")
   {
      Clear(opt);
   }
   virtual void Print(Option_t* option = "") const;
 
   void AddHit(KVNucleus* part)
   {
 
      if (!fParticles) {
         fParticles = new KVList(kFALSE);
         fParticles->SetCleanup();
      }
      fParticles->Add(part);
      SetAnalysed();
   }
 
   void RemoveHit(KVNucleus* part)
   {
 
      fParticles->Remove(part);
      if (!fParticles->GetEntries()) SetAnalysed(kFALSE);
   }
 
   KVList* GetHits() const
   {
      return fParticles;
   }
   void ClearHits()
   {
      if (fParticles) fParticles->Clear();
   }
   Int_t GetNHits() const
   {
      return (fParticles ? fParticles->GetEntries() : 0);
   }
 
   inline UShort_t GetSegment() const;
   inline virtual void SetSegment(UShort_t s);
   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()->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);
   }
 
 
   virtual void AddIDTelescope(TObject* idt);
   KVList* GetIDTelescopes()
   {
      return fIDTelescopes;
   }
   KVList* GetAlignedIDTelescopes();
   TList* GetTelescopesForIdentification();
 
   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 TGeoVolume* GetGeoVolume();
   virtual void AddToGeometry();
   virtual void GetVerticesInOwnFrame(TVector3* /*corners[8]*/, Double_t /*depth*/, Double_t /*layer_thickness*/);
   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();
   }
 
   virtual Double_t GetMaxDeltaE(Int_t Z, Int_t A);
   virtual Double_t GetEIncOfMaxDeltaE(Int_t Z, Int_t A);
   virtual Double_t GetDeltaE(Int_t Z, Int_t A, Double_t Einc);
   virtual Double_t GetTotalDeltaE(Int_t Z, Int_t A, Double_t Einc);
   virtual Double_t GetERes(Int_t Z, Int_t A, Double_t Einc);
   virtual Double_t GetIncidentEnergy(Int_t Z, Int_t A, Double_t delta_e =
                                         -1.0, enum SolType type = kEmax);
   /*virtual Double_t GetEResFromDeltaE(...)  - DON'T IMPLEMENT, CALLS GETINCIDENTENERGY*/
   virtual Double_t GetDeltaEFromERes(Int_t Z, Int_t A, Double_t Eres);
   virtual Double_t GetIncidentEnergyFromERes(Int_t Z, Int_t A, Double_t Eres);
   virtual Double_t GetRange(Int_t Z, Int_t A, Double_t Einc);
   virtual Double_t GetLinearRange(Int_t Z, Int_t A, Double_t Einc);
   virtual Double_t GetPunchThroughEnergy(Int_t Z, Int_t A);
   virtual TGraph* DrawPunchThroughEnergyVsZ(Int_t massform = KVNucleus::kBetaMass);
   virtual 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);
 
   virtual void SetEResAfterDetector(Double_t e)
   {
      fEResforEinc = e;
   }
   virtual Double_t GetEResAfterDetector() const
   {
      return fEResforEinc;
   }
 
   virtual void ReadDefinitionFromFile(const Char_t*);
 
   virtual TList* GetAlignedDetectors(UInt_t direction = /*KVGroup::kBackwards*/ 1);
   void ResetAlignedDetectors(UInt_t direction = /*KVGroup::kBackwards*/ 1);
 
   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);
   }
 
   virtual void DeduceACQParameters(KVEvent*, KVNumberList&) {}
 
   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
   {
      if (ROOTGeo()) return fEWPosition.GetSolidAngle();
      return KVPosition::GetSolidAngle();
   }
   TVector3 GetRandomDirection(Option_t* t = "isotropic")
   {
      if (ROOTGeo()) return fEWPosition.GetRandomDirection(t);
      return KVPosition::GetRandomDirection(t);
   }
   void GetRandomAngles(Double_t& th, Double_t& ph, Option_t* t = "isotropic")
   {
      if (ROOTGeo()) fEWPosition.GetRandomAngles(th, ph, t);
      else KVPosition::GetRandomAngles(th, ph, t);
   }
   TVector3 GetDirection()
   {
      if (ROOTGeo()) return fEWPosition.GetDirection();
      return KVPosition::GetDirection();
   }
   Double_t GetDistance() const
   {
      if (ROOTGeo()) return fEWPosition.GetDistance();
      return KVPosition::GetDistance();
   }
   Double_t GetTheta() const
   {
      if (ROOTGeo()) return fEWPosition.GetTheta();
      return KVPosition::GetTheta();
   }
   Double_t GetSinTheta() const
   {
      if (ROOTGeo()) return fEWPosition.GetSinTheta();
      return KVPosition::GetSinTheta();
   }
   Double_t GetCosTheta() const
   {
      if (ROOTGeo()) return fEWPosition.GetCosTheta();
      return KVPosition::GetCosTheta();
   }
   Double_t GetPhi() const
   {
      if (ROOTGeo()) return fEWPosition.GetPhi();
      return KVPosition::GetPhi();
   }
 
   void SetThickness(Double_t thick);
   Bool_t IsSingleLayer() const
   {
      return fSingleLayer;
   }
   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;
   }
   KVDetectorSignal* AddDetectorSignal(const KVString& type)
   {
 
      auto signal = new KVDetectorSignal(type, this);
      fDetSignals.Add(signal);
      return signal;
   }
   Bool_t AddDetectorSignalExpression(const KVString& type, const KVString& _expr);
 
   virtual Int_t GetIndex() const
   {
 
      return 0;
   }
 
   ClassDef(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;
}
 
inline UShort_t KVDetector::GetSegment() const
{
   return fSegment;
}
 
inline void KVDetector::SetSegment(UShort_t s)
{
   fSegment = s;
}
 
inline void KVDetector::SetGain(Double_t gain)
{
   fGain = gain;
}
 
inline Double_t KVDetector::GetGain() const
{
   return fGain;
}
#endif