INDRAGeometryBuilder.cpp

//Created by KVClassFactory on Wed Feb  1 11:09:41 2012
//Author: John Frankland,,,
 
#include "INDRAGeometryBuilder.h"
#include "KVTarget.h"
#include "KVUnits.h"
#include "TEnv.h"
#include "TMath.h"
#include "TString.h"
#include "TGeoMedium.h"
#include "TRotation.h"
#include "TGeoMatrix.h"
#include "TGeoManager.h"
#include "TGeoCompositeShape.h"
#include "TGeoEltu.h"
#include "KVINDRA.h"
#include "KVDataSet.h"
 
using namespace std;
 
//#define DEBUG 1
 
ClassImp(INDRAGeometryBuilder)
 
 
 
void INDRAGeometryBuilder::read_indra_struct_file()
{
   // Read infos in dataset file
   //
   //   indra-struct.env
   //
   // or alternatively, by defining the variable
   //
   //   [dataset].INDRA.StructureFile: [name of dataset]
   //
   // use the file from another dataset.
   //
   // These files were previously used for all geometry specifications;
   // now they mainly serve to define which detectors/rings are present,
   // and change the default nominal thickness of silicon detectors or
   // define the default gas pressures in the ionization chambers
 
   TString alt_dataset = GetDataSetEnv(fDataSet, "INDRA.StructureFile", "");
   if(alt_dataset=="") alt_dataset=fDataSet;
 
   auto path = KVDataSet::GetFullPathToDataSetFile(alt_dataset,"indra-struct.env");
 
   Info("read_indra_struct_file", "Reading INDRA structure from file : %s", path.Data());
   fStrucInfos.ReadFile(path, kEnvAll);
 
   auto run_index_multiplier = GetDataSetEnv(alt_dataset, "DataSet.RunFileIndexMultiplier.raw", 1);
 
   //test if additional geometrical specifications for certain runs exists
   KVString lruns = fStrucInfos.GetValue("AddOnForRuns", "");
   if (lruns != "") {
      lruns.Begin(",");
      while (!lruns.End()) {
         KVString sruns = lruns.Next();
         KVNumberList nlr(sruns.Data());
         //the current run needs specific geometry
         if (nlr.Contains(fCurrentRun*run_index_multiplier))
         {
            path = fStrucInfos.GetValue(sruns.Data(), "");
            Info("BuildGeometry", "Additional geometry for run=%d in file #%s#", fCurrentRun, path.Data());
            path = KVDataSet::GetFullPathToDataSetFile(alt_dataset,path.Data());
            fStrucInfos.ReadFile(path, kEnvChange);
         }
      }
   }
 
   // to be used as name of KVINDRA object
   SetName(fStrucInfos.GetValue("INDRA.Name", ""));
   SetTitle(fStrucInfos.GetValue("INDRA.Title", ""));
 
   KVString layers = fStrucInfos.GetValue("INDRA.Layers", "");
   layers.Begin(" ");
   while (!layers.End()) read_layer_infos(layers.Next());
 
}
 
 
 
 
void INDRAGeometryBuilder::read_layer_infos(const Char_t* name)
{
   // read infos on layer 'name' in file "$KVROOT/KVFiles/data/indra-struct.[dataset].env"
 
   Int_t number = fStrucInfos.GetValue(Form("INDRA.Layer.%s.Number", name), 0);
   Info("read_layer_infos", "Building layer %d:%s", number, name);
   KVNumberList rings = fStrucInfos.GetValue(Form("INDRA.Layer.%s.Rings", name), "");
   rings.Begin();
   while (!rings.End()) {
      Int_t ring_num = rings.Next();
      TString prefix = Form("INDRA.Layer.%s.Ring.%d", name, ring_num);
      read_ring_infos(ring_num, prefix);
   }
}
 
 
 
 
void INDRAGeometryBuilder::read_ring_infos(Int_t ring_number, const Char_t* prefix)
{
   // read infos on ring in file "$KVROOT/KVFiles/data/indra-struct.[dataset].env"
 
   Info("read_ring_infos", "Building ring %d (%s)", ring_number, prefix);
 
   Int_t ntel = fStrucInfos.GetValue(Form("%s.Ntel", prefix), 0);
   Int_t step = fStrucInfos.GetValue(Form("%s.Step", prefix), 1);
   Int_t num_first = fStrucInfos.GetValue(Form("%s.Nfirst", prefix), 1);
   Double_t thmin = fStrucInfos.GetValue(Form("%s.ThMin", prefix), 0.0);
   Double_t thmax = fStrucInfos.GetValue(Form("%s.ThMax", prefix), 0.0);
   Double_t phi_0 = fStrucInfos.GetValue(Form("%s.Phi0", prefix), 0.0);
   Double_t dphi = fStrucInfos.GetValue(Form("%s.Dphi", prefix), -1.0);
   dphi = (dphi < 0. ? 360. / ntel : dphi);
   KVNumberList absent = fStrucInfos.GetValue(Form("%s.Absent", prefix), "");
 
   Double_t phi = phi_0;
   UInt_t counter = num_first;
   auto det_props = read_telescope_infos(prefix, ring_number, counter, ntel, thmin, thmax, phi);
 
   for (int i = 0; i < ntel; i++) {
      if (!absent.Contains(counter)) {
         for (auto& d : det_props) detector_map[d.name] = d;
      }
      phi = (phi + dphi > 360. ? (phi + dphi - 360.) : (phi + dphi));
      counter += step;
      if (i + 1 < ntel) det_props = read_telescope_infos(prefix, ring_number, counter, ntel, thmin, thmax, phi);
   }
}
 
 
 
 
std::vector<INDRAGeometryBuilder::detector_properties>
INDRAGeometryBuilder::read_telescope_infos(
        const Char_t* prefix, Int_t ring, Int_t module, Int_t ntel,
        double thmin, double thmax, double phi)
{
   // read telescope infos in file "$KVROOT/KVFiles/data/indra-struct.[dataset].env"
 
   KVString telescopes = fStrucInfos.GetValue(Form("%s.Telescope", prefix), "");
   telescopes.Begin(" ");
   KVString name;
   while (!telescopes.End()) {
      name = telescopes.Next();
      KVNumberList mods = fStrucInfos.GetValue(Form("%s.Telescope.%s.Modules", prefix, name.Data()), "1-100");
      if (mods.Contains(module)) break;
   }
 
   KVString detectors = fStrucInfos.GetValue(Form("INDRA.Telescope.%s.Detectors", name.Data()), "");
   double aziwidth = fStrucInfos.GetValue(Form("INDRA.Telescope.%s.AziWidth", name.Data()), 0.0);
   double phimin = phi-0.5*aziwidth;
   double phimax = phi+0.5*aziwidth;
 
   detectors.Begin(" ");
   Int_t idet = 1;
   std::vector<detector_properties> det_props;
   while (!detectors.End()) {
      KVString det = detectors.Next();
      det.Begin("(,)");
      KVString dettype = det.Next();
 
      // set flag if PHOSwich ring is used
      if (dettype == "PHOS") fWithPhoswich = kTRUE;
 
      Double_t detthick = det.Next().Atof();
      if (dettype == "SILI" || dettype == "SI75")
         det_props.push_back(detector_properties{Form("%s_%02d", dettype.Data(), ring), detthick, ring, 0, {thmin, thmax}, {phimin, phimax}, true});
      else if (dettype == "PHOS") {
         // for phoswich thicknesses are given as PHOS(0.05,25)
         // for geometry definition, '25cm' layer is set as active layer - so take second thickness value
         detthick = det.Next().Atof();
         det_props.push_back(detector_properties{Form("%s_%02d", dettype.Data(), module), detthick, 1, module, {thmin, thmax}, {phimin, phimax}, true});
      }
      else if (dettype == "SI")
         // silicon detector thicknesses are given in microns
         det_props.push_back(detector_properties{Form("%s_%02d%02d", dettype.Data(), ring, module), detthick * KVUnits::um, ring, module, {thmin, thmax}, {phimin, phimax}, true});
      else
         det_props.push_back(detector_properties{Form("%s_%02d%02d", dettype.Data(), ring, module), detthick, ring, module, {thmin, thmax}, {phimin, phimax}, true});
 
      idet++;
   }
   return det_props;
}
 
 
 
 
void INDRAGeometryBuilder:: CorrectCoordinates(Double_t* coo, Double_t& offX, Double_t& offY)
{
   // calculate offset in X and Y
   // correct coordinates for offset
   // return offset TGeoTranslation
 
   offX = offY = 0;
   for (int i = 0; i < 8; i++) {
      offX += coo[2 * i];
      offY += coo[2 * i + 1];
   }
   offX /= 8.;
   offY /= 8.;
   for (int i = 0; i < 8; i++) {
      coo[2 * i] -= offX;
      coo[2 * i + 1] -= offY;
   }
}
 
 
 
 
INDRAGeometryBuilder::INDRAGeometryBuilder(const TString& _data_set, Int_t current_run)
   : fDataSet(_data_set), fCurrentRun(current_run)
{
   // Default constructor
   fEtalonVol = 0;
}
 
 
 
 
void INDRAGeometryBuilder::CalculateBackPlaneCoordinates(TVector3* frontcoords, TVector3 normal, Double_t depth, TVector3* backcoords)
{
   // Given an array of ncoords points frontcoords[ncoords] which lie in a plane
   // calculate the corresponding points on a parallel plane which is
   // at a distance from the origin which is greater by depth
   // (distance plane-origin = length of normal to plane passing through origin)
   // 'normal' is vector point of intersection plane-normal to origin
 
   Double_t planeDist = normal.Mag();//distance original plane-origin
   for (int i = 0; i < 4; i++) {
      Double_t rap = (planeDist + depth) / planeDist;
      backcoords[i] = rap * frontcoords[i];
   }
}
 
 
 
 
void INDRAGeometryBuilder::ReadDetCAO(const Char_t* detname, Int_t ring)
{
   // Read infos in CAO file for this detector.
   // We look for a TEnv file with name "detname.cao"
 
   TEnv infos;
   TString path = Form("%s.cao", detname);
   SearchKVFile(path.Data(), path, "data");
   infos.ReadFile(path, kEnvAll);
 
   Double_t unit_converter = 0.1;// distances in file given in mm: convert to cm
 
   fInnerFront[0].SetX(infos.GetValue("INNER.A.X", 0.0)*unit_converter);
   fInnerFront[0].SetY(infos.GetValue("INNER.A.Y", 0.0)*unit_converter);
   fInnerFront[0].SetZ(infos.GetValue("INNER.A.Z", 0.0)*unit_converter);
   fInnerFront[3].SetX(infos.GetValue("INNER.B.X", 0.0)*unit_converter);
   fInnerFront[3].SetY(infos.GetValue("INNER.B.Y", 0.0)*unit_converter);
   fInnerFront[3].SetZ(infos.GetValue("INNER.B.Z", 0.0)*unit_converter);
   fInnerFront[2].SetX(infos.GetValue("INNER.C.X", 0.0)*unit_converter);
   fInnerFront[2].SetY(infos.GetValue("INNER.C.Y", 0.0)*unit_converter);
   fInnerFront[2].SetZ(infos.GetValue("INNER.C.Z", 0.0)*unit_converter);
   fInnerFront[1].SetX(infos.GetValue("INNER.D.X", 0.0)*unit_converter);
   fInnerFront[1].SetY(infos.GetValue("INNER.D.Y", 0.0)*unit_converter);
   fInnerFront[1].SetZ(infos.GetValue("INNER.D.Z", 0.0)*unit_converter);
 
   // To get back to INDRA phi-coordinates, all points need to be rotated
   // by -150 degrees around the Z-axis
   for (int i = 0; i <= 3; i++) {
      fInnerFront[i].RotateZ(-150.*TMath::DegToRad());
   }
   CalculateCentre(fInnerFront, fInnerCentre);
 
   fOuterFront[0].SetX(infos.GetValue("OUTER.A.X", 0.0)*unit_converter);
   fOuterFront[0].SetY(infos.GetValue("OUTER.A.Y", 0.0)*unit_converter);
   fOuterFront[0].SetZ(infos.GetValue("OUTER.A.Z", 0.0)*unit_converter);
   fOuterFront[3].SetX(infos.GetValue("OUTER.B.X", 0.0)*unit_converter);
   fOuterFront[3].SetY(infos.GetValue("OUTER.B.Y", 0.0)*unit_converter);
   fOuterFront[3].SetZ(infos.GetValue("OUTER.B.Z", 0.0)*unit_converter);
   fOuterFront[2].SetX(infos.GetValue("OUTER.C.X", 0.0)*unit_converter);
   fOuterFront[2].SetY(infos.GetValue("OUTER.C.Y", 0.0)*unit_converter);
   fOuterFront[2].SetZ(infos.GetValue("OUTER.C.Z", 0.0)*unit_converter);
   fOuterFront[1].SetX(infos.GetValue("OUTER.D.X", 0.0)*unit_converter);
   fOuterFront[1].SetY(infos.GetValue("OUTER.D.Y", 0.0)*unit_converter);
   fOuterFront[1].SetZ(infos.GetValue("OUTER.D.Z", 0.0)*unit_converter);
 
   // To get back to INDRA phi-coordinates, all points need to be rotated
   // by -150 degrees around the Z-axis
   for (int i = 0; i <= 3; i++) {
      fOuterFront[i].RotateZ(-150.*TMath::DegToRad());
   }
   CalculateCentre(fOuterFront, fOuterCentre);
 
   phi0 = infos.GetValue("PHIZERO", 30.0);
   thetaMin = infos.GetValue("THETAMIN", 0.0);
   thetaMax = infos.GetValue("THETAMAX", 0.0);
   Ndets = infos.GetValue("NDETS", 1);
   deltaPhi = 360. / Ndets;
 
   // calculate theoretical geometry:
   // points of intersection with the plane defined by real Huguet geometry
   // and the 4 vectors to theoretical thetamin/max phimin/max
   Double_t Ph0 = fInnerCentre.Phi() * TMath::RadToDeg();
   CalculateCornersInPlane(fInnerFront, thetaMin, thetaMax,
                           Ph0 - deltaPhi / 2.,
                           Ph0 + deltaPhi / 2.,
                           fFrameFront);
   // calculate centre of frame = direction of detector
   CalculateCentre(fFrameFront, fFrameCentre);
   fFrameMat.SetMaterial(infos.GetValue("FRAME.MATERIAL", "Al"));
 
   fInnerPads = infos.GetValue("INNER.PADS", 1);
   fOuterPads = infos.GetValue("OUTER.PADS", 0);
   fInnerRing = infos.GetValue("INNER.RING", ring);
   fOuterRing = infos.GetValue("OUTER.RING", 0);
   fInnerDmod = infos.GetValue("INNER.DMOD", 1);
   fOuterDmod = infos.GetValue("OUTER.DMOD", 1);
   fInnerModmin = infos.GetValue("INNER.MODMIN", 1);
   fOuterModmin = infos.GetValue("OUTER.MODMIN", 1);
   fModmax = infos.GetValue("MODMAX", 24);
 
   // read infos on internal structure of detector
   auto nlay = infos.GetValue("LAYERS", 1);
 
   fActiveLayer = infos.GetValue("ACTIVELAYER", 1);
   for (int lay = 1; lay <= nlay; ++lay) {
      auto mat = new KVMaterial;
      mat->SetMaterial(infos.GetValue(Form("LAYER.%d.MATERIAL", lay), ""));
      auto gas_pressure = infos.GetValue(Form("LAYER.%d.PRESSURE", lay), -1.0);
      if (gas_pressure > 0) mat->SetPressure(gas_pressure * KVUnits::torr); // gas pressure in Torr
      auto thick = infos.GetValue(Form("LAYER.%d.THICK", lay), -1.0);
      if (thick > 0) mat->SetThickness(thick * KVUnits::cm); // linear thickness in cm
      auto dens = infos.GetValue(Form("LAYER.%d.DENSITY", lay), -1.0);
      if (dens > 0) mat->SetAreaDensity(dens * KVUnits::ug / KVUnits::cm / KVUnits::cm); // area density in ug/cm2
 
      fCurrentDetector.AddAbsorber(mat);
   }
   if (fActiveLayer > 1) fCurrentDetector.SetActiveLayer(fActiveLayer - 1);
}
 
 
 
 
 
void INDRAGeometryBuilder::Print(Option_t*) const
{
   printf("Detector name : %s\n", fDetName.Data());
   printf("Total thickness : %f cm\n", fCurrentDetector.GetTotalThicknessInCM());
   printf("Frame centre:\n");
   fFrameCentre.Print();
   printf("Inner centre:\n");
   fInnerCentre.Print();
   printf("Outer centre: \n");
   fOuterCentre.Print();
   printf("Coordinates of outer pads :\n");
   for (int i = 0; i < 4; i++) fOuterFront[i].Print();
   printf("Coordinates of inner pads :\n");
   for (int i = 0; i < 4; i++) fInnerFront[i].Print();
   printf("Material of outer frame : %s\n", fFrameMat.GetTitle());
}
 
 
 
 
 
void INDRAGeometryBuilder::MakeFrame(TString det_type, Int_t ring_num)
{
   // Create the TGeoVolume corresponding to the outer casing
   // (dead zone) of the detector
 
   TVector3 fFrameBack[4];
   TVector3 fFrameBackCentre;
   //CalculateBackPlaneCoordinates(fFrameFront, fFrameCentre, fTotalThickness, fFrameBack);
   // make deadzone be 1um thick
   Double_t dz = 5.e-5;
   CalculateBackPlaneCoordinates(fFrameFront, fFrameCentre, 2 * dz, fFrameBack);
   CalculateCentre(fFrameBack, fFrameBackCentre);
 
   TVector3 corners[8]; // 8 vertices of the volume
   TransformToOwnFrame(fFrameFront, fFrameCentre, corners);
   TransformToOwnFrame(fFrameBack, fFrameBackCentre, &corners[4]);
 
   Double_t vertices[16];
   for (int i = 0; i < 8; i++) {
      vertices[2 * i] = corners[i].X();
      vertices[2 * i + 1] = corners[i].Y();
   }
 
   //Double_t dz = 0.99 * fTotalThickness / 2.;
   TString vol_name;
   vol_name.Form("STRUCT_%s_%02d", det_type.Data(), ring_num);
   fFrameVolume = gGeoManager->MakeVolumeAssembly(vol_name);
   fFrameVolume->SetMedium(gGeoManager->GetMedium("Vacuum"));
   vol_name.Form("DEADZONE_%s_%02d", det_type.Data(), ring_num);
   TGeoMedium* med = fFrameMat.GetGeoMedium();
   Double_t offX, offY;
   CorrectCoordinates(vertices, offX, offY);
   // place the deadzone 1um (=1.e-4 cm) behind the detector entrance window
   TGeoTranslation* offset  = new TGeoTranslation(offX, offY, -(0.5 * fCurrentDetector.GetTotalThicknessInCM() - (1.e-4 + dz)));
   TGeoVolume* frame = gGeoManager->MakeArb8(vol_name.Data(), med, dz, vertices);
   frame->SetLineColor(med->GetMaterial()->GetDefaultColor());
   fFrameVolume->AddNode(frame, 1, offset);
}
 
 
 
 
 
void INDRAGeometryBuilder::TransformToOwnFrame(TVector3* orig, TVector3& centre, TVector3* ownframe)
{
   // takes the four points (TVector3) in the array and transforms them
   // to a reference frame in which the +ve z-axis goes through their centre
 
   TRotation rot_to_frame;
   rot_to_frame.SetYEulerAngles(-centre.Phi(), -centre.Theta(), 0.);
   TVector3 displZ(0, 0, centre.Mag());
   for (int i = 0; i < 4; i++) {
      ownframe[i] = rot_to_frame * orig[i] - displZ;
   }
}
 
 
 
 
Bool_t INDRAGeometryBuilder::SetDetectorThickness(const TString& detname, Double_t thickness)
{
   // Called by KVINDRA::set_detector_thicknesses().
   //
   // Is 'detname' is present in the current geometry, we override its nominal thickness
   // (read from the indra-struct* file) with the given value and return kTRUE.
   //
   // If 'detname' is not present in the current geometry, we print a warning message and return
   // kFALSE (incoherency between the indra-struct* present detectors
   // & the list of detector thicknesses).
 
   auto det = detector_map.find(detname.Data());
   if (det != detector_map.end()) {
      det->second.thick = thickness;
      return kTRUE;
   }
   Warning("SetDetectorThickness", "New thickness value %f given for %s,"
           " but detector is apparently absent from geometry?",
           thickness, detname.Data());
   return kFALSE;
}
 
 
 
 
 
void INDRAGeometryBuilder::PlaceFrame(Double_t phi, Int_t copy_no)
{
   // position frame (dead zone) volume in geometry
 
   Double_t theta = fFrameCentre.Theta() * TMath::RadToDeg();
   TGeoRotation rot1, rot2;
   rot2.SetAngles(phi + 90., theta, 0.);
   rot1.SetAngles(-90., 0., 0.);
   // distance to frame centre
   Double_t trans_z = fFrameCentre.Mag() + fCurrentDetector.GetTotalThicknessInCM() / 2.;
 
   TGeoTranslation tran(0, 0, trans_z);
   TGeoHMatrix h = rot2 * tran * rot1;
   TGeoHMatrix* ph = new TGeoHMatrix(h);
 
   // add detector volume to geometry
   gGeoManager->GetTopVolume()->AddNode(fFrameVolume, copy_no, ph);
}
 
 
 
 
 
void INDRAGeometryBuilder::PlaceDetector()
{
   // position detector inside frame
 
   fFrameVolume->AddNode(fDetVolume, 1, fDetectorPosition);
   fFrameVolume->GetNode(Form("%s_1", fDetVolume->GetName()))->SetName(fDetVolume->GetName());
}
 
 
 
 
Bool_t INDRAGeometryBuilder::CheckDetectorPresent(TString detname)
{
#ifdef DEBUG
   Info("CheckDetectorPresent", "%s - present:%d thick:%f",
        detname.Data(), detector_map[detname.Data()].present, detector_map[detname.Data()].thick);
#endif
   auto it = detector_map.find(detname.Data());
   return (it !=  detector_map.end()) && it->second.present;
}
 
 
 
 
void INDRAGeometryBuilder::MakeRing(const Char_t* det, int ring)
{
 
   fCurrentDetectorType = det;
   fCurrentDetector.RemoveAllAbsorbers();
 
   if (!strcmp(det, "PHOS"))
      fDetName.Form("%s_01", det);
   else
      fDetName.Form("%s_%02d01", det, ring);
   ReadDetCAO(fDetName, ring);
   //Print();
   Double_t phi = phi0;
   Int_t innerMod = fInnerModmin;
   Int_t outerMod = fOuterModmin;
   int i = 1;
 
   cout << fDetName.Data()  << "  " << fInnerCentre.Mag() << "  " <<
        0.5 * (fInnerFront[0].Theta() + fInnerFront[3].Theta())*TMath::RadToDeg() << "  " <<
        0.5 * (fInnerFront[1].Theta() + fInnerFront[2].Theta())*TMath::RadToDeg() << "  " <<
        0.5 * (fInnerFront[2].Phi() + fInnerFront[3].Phi())*TMath::RadToDeg() << "  " <<
        0.5 * (fInnerFront[0].Phi() + fInnerFront[1].Phi())*TMath::RadToDeg() << endl;
 
// check presence of all detectors of ring
   Int_t npresent = 0;
   Int_t ntotal = 0;
 
   for (i = 1; i <= Ndets; i++) {
      if (!strcmp(det, "PHOS")) fDetName.Form("%s_%02d", det, innerMod);
      else fDetName.Form("%s_%02d%02d", det, fInnerRing, innerMod);
      ntotal++;
      npresent += (Int_t)CheckDetectorPresent(fDetName);
 
      innerMod += fInnerDmod;
 
      if (fInnerPads == 3) {
         if (!strcmp(det, "PHOS")) fDetName.Form("%s_%02d", det, innerMod);
         else fDetName.Form("%s_%02d%02d", det, fInnerRing, innerMod);
         ntotal++;
         npresent += (Int_t)CheckDetectorPresent(fDetName);
         innerMod += fInnerDmod;
         if (!strcmp(det, "PHOS")) fDetName.Form("%s_%02d", det, innerMod);
         else fDetName.Form("%s_%02d%02d", det, fInnerRing, innerMod);
         ntotal++;
         npresent += (Int_t)CheckDetectorPresent(fDetName);
         innerMod += fInnerDmod;
      }
      else {
         if (fInnerPads == 2) {
            if (!strcmp(det, "PHOS")) fDetName.Form("%s_%02d", det, innerMod);
            else fDetName.Form("%s_%02d%02d", det, fInnerRing, innerMod);
            ntotal++;
            npresent += (Int_t)CheckDetectorPresent(fDetName);
            innerMod += fInnerDmod;
         }
         if (fOuterPads) {
            if (!strcmp(det, "PHOS")) fDetName.Form("%s_%02d", det, outerMod);
            else fDetName.Form("%s_%02d%02d", det, fOuterRing, outerMod);
            ntotal++;
            npresent += (Int_t)CheckDetectorPresent(fDetName);
            outerMod += fOuterDmod;
            if (fOuterPads == 2) {
               if (!strcmp(det, "PHOS")) fDetName.Form("%s_%02d", det, outerMod);
               else fDetName.Form("%s_%02d%02d", det, fOuterRing, outerMod);
               ntotal++;
               npresent += (Int_t)CheckDetectorPresent(fDetName);
               outerMod += fOuterDmod;
            }
         }
      }
   }
   if (!npresent) {
      Info("MakeRing", "Detector type %s  ring number %d : ABSENT", det, ring);
      return;
   }
   Info("MakeRing", "Detector type %s  ring number %d : %d/%d detectors present", det, ring, npresent, ntotal);
 
   innerMod = fInnerModmin;
   outerMod = fOuterModmin;
 
   bool made_prototype_block{false};
 
   for (i = 1; i <= Ndets; i++) {
 
      // check absence of whole block
      ntotal = 0;
      npresent = 0;
      if (!strcmp(det, "PHOS")) fDetName.Form("%s_%02d", det, innerMod);
      else fDetName.Form("%s_%02d%02d", det, fInnerRing, innerMod);
      ntotal++;
      auto first_is_present = CheckDetectorPresent(fDetName);
      npresent += (Int_t)first_is_present;
      // check if thickness/pressure of first detector in block has changed compared to current prototype
      if (first_is_present) {
         if (fCurrentDetector.IsGasDetector()) {
//            if(detector_map[fDetName.Data()].thick != fCurrentDetector.GetPressure()) {
//               // change pressure of gas
//               fCurrentDetector.SetPressure(detector_map[fDetName.Data()].thick);
//               made_prototype_block = false; // set flag to generate new prototype
//            }
         }
         else {
            if (detector_map[fDetName.Data()].thick != fCurrentDetector.GetThickness()) {
               // change thickness of current active layer
               fCurrentDetector.SetThickness(detector_map[fDetName.Data()].thick);
               made_prototype_block = false; // set flag to generate new prototype
            }
         }
      }
 
      innerMod += fInnerDmod;
 
      if (fInnerPads == 3) {
         if (!strcmp(det, "PHOS")) fDetName.Form("%s_%02d", det, innerMod);
         else fDetName.Form("%s_%02d%02d", det, fInnerRing, innerMod);
         ntotal++;
         npresent += (Int_t)CheckDetectorPresent(fDetName);
         innerMod += fInnerDmod;
         if (!strcmp(det, "PHOS")) fDetName.Form("%s_%02d", det, innerMod);
         else fDetName.Form("%s_%02d%02d", det, fInnerRing, innerMod);
         ntotal++;
         npresent += (Int_t)CheckDetectorPresent(fDetName);
         innerMod += fInnerDmod;
      }
      else {
         if (fInnerPads == 2) {
            if (!strcmp(det, "PHOS")) fDetName.Form("%s_%02d", det, innerMod);
            else fDetName.Form("%s_%02d%02d", det, fInnerRing, innerMod);
            ntotal++;
            npresent += (Int_t)CheckDetectorPresent(fDetName);
            innerMod += fInnerDmod;
         }
         if (fOuterPads) {
            if (!strcmp(det, "PHOS")) fDetName.Form("%s_%02d", det, outerMod);
            else fDetName.Form("%s_%02d%02d", det, fOuterRing, outerMod);
            ntotal++;
            npresent += (Int_t)CheckDetectorPresent(fDetName);
            outerMod += fOuterDmod;
            if (fOuterPads == 2) {
               if (!strcmp(det, "PHOS")) fDetName.Form("%s_%02d", det, outerMod);
               else fDetName.Form("%s_%02d%02d", det, fOuterRing, outerMod);
               ntotal++;
               npresent += (Int_t)CheckDetectorPresent(fDetName);
               outerMod += fOuterDmod;
            }
         }
      }
      if (!npresent) Info("MakeRing", "\tBlock %d: ABSENT", i);
      else if (npresent < ntotal) {
         Info("MakeRing", "\tBlock %d: %d/%d detectors present", i, npresent, ntotal);
         Info("MakeRing", "\tCASE NOT TREATED. ALL DETECTORS PRESENT.");
      }
      if (npresent) {
 
         if (!made_prototype_block) {
            // Create prototype for ring
            made_prototype_block = true;
 
            Info("MakeRing", "DET:%s RING:%d  Making new prototype group for i=%d/%d",
                 det, ring, i, Ndets);
 
            /* build and add frame */
            MakeFrame(det, ring);
 
            // make pads in inner ring
            MakeDetector("A1", fInnerFront, fInnerCentre);
            PlaceDetector();
            if (fInnerPads == 3) {
               MakeDetector("A2", fOuterFront, fOuterCentre);
               PlaceDetector();
               TVector3 reflex[4];
               ReflectPad(fInnerFront, fFrameCentre.Phi(), reflex);
               TVector3 refcent;
               CalculateCentre(reflex, refcent);
               MakeDetector("A3", reflex, refcent);
               PlaceDetector();
            }
            else {
               if (fInnerPads == 2) {
                  TVector3 reflex[4];
                  ReflectPad(fInnerFront, fFrameCentre.Phi(), reflex);
                  TVector3 refcent;
                  CalculateCentre(reflex, refcent);
                  MakeDetector("A2", reflex, refcent);
                  PlaceDetector();
               }
               if (fOuterPads) {
                  // make pads in outer ring
                  MakeDetector("B1", fOuterFront, fOuterCentre);
                  PlaceDetector();
                  if (fOuterPads == 2) {
                     TVector3 reflex[4];
                     ReflectPad(fOuterFront, fFrameCentre.Phi(), reflex);
                     TVector3 refcent;
                     CalculateCentre(reflex, refcent);
                     MakeDetector("B2", reflex, refcent);
                     PlaceDetector();
                  }
               }
            }
         }
 
         PlaceFrame(phi, i);
      }
 
      if (((ring > 9 && ring < 13) && i == 1) || (ring > 12 && i == 2)) {
         fEtalonTheta[ring - 10] = fFrameCentre.Theta() * TMath::RadToDeg();
         fEtalonPhi[ring - 10] = phi;
         fEtalonTheta[ring - 9] = fFrameCentre.Theta() * TMath::RadToDeg();
         fEtalonPhi[ring - 9] = phi;
      }
      phi += deltaPhi;
   }
 
}
 
 
 
 
 
void INDRAGeometryBuilder::MakeDetector(const Char_t* det, TVector3* som, TVector3 cen)
{
 
   // make volume corresponding to the actual detector
 
   TVector3 corners[8]; // 8 vertices of the volume
   Double_t vertices[16];
 
   // Check whether this is a real multi-layer detector
   // i.e. make sure all "layers" have non-zero thickness
   Bool_t multi_layer = fCurrentDetector.IsMultiLayer();
   if (multi_layer) {
      TIter next(fCurrentDetector.GetListOfAbsorbers());
      KVMaterial* abs;
      Int_t no_abs = 0;
      while ((abs = (KVMaterial*)next())) {
         Double_t thick = abs->GetThickness();
         if (thick > 0.0) ++no_abs;
      }
      multi_layer = (no_abs > 1);
   }
   if (multi_layer) {
      fDetVolume = gGeoManager->MakeVolumeAssembly(Form("DET_%s", det));
      fDetVolume->SetMedium(gGeoManager->GetMedium("Vacuum"));
   }
   TVector3 frontPlane[4], backPlane[4], frontCentre, backCentre;
   // front plane of first absorber is front plane of detector
   for (int i = 0; i < 4; i++) frontPlane[i] = som[i];
   frontCentre = cen;
 
   /**** BUILD & ADD LAYERS ****/
   TIter next(fCurrentDetector.GetListOfAbsorbers());
   KVMaterial* abs;
   Double_t depth_in_det = 0.;
   Int_t no_abs = 1;
 
   Double_t offX, offY; //offset of each layer
   fDetectorPosition = 0;
 
   auto fTotalThickness = fCurrentDetector.GetTotalThicknessInCM();
 
   while ((abs = (KVMaterial*)next())) {
      // get medium for absorber
      TGeoMedium* med = abs->GetGeoMedium();
      Double_t thick = abs->GetThickness();
      if (thick == 0.0) {
         no_abs++;   // ignore zero thickness layers
         continue;
      }
 
      // calculate coordinates of back plane
      CalculateBackPlaneCoordinates(frontPlane, frontCentre, thick, backPlane);
      // calculate centre of back plane
      CalculateCentre(backPlane, backCentre);
      // calculate corners of volume
      TransformToOwnFrame(frontPlane, fFrameCentre, corners);
      TransformToOwnFrame(backPlane, fFrameCentre, &corners[4]);
      for (int i = 0; i < 8; i++) {
         vertices[2 * i] = corners[i].X();
         vertices[2 * i + 1] = corners[i].Y();
      }
 
      Double_t dz = thick / 2.;
      TString vol_name;
      if (multi_layer) {
         if (no_abs == fActiveLayer) vol_name = Form("ACTIVE_%s", det);
         else vol_name = Form("%s_%d_%s", det, no_abs, abs->GetName());
      }
      else
         vol_name = Form("DET_%s", det);
 
      CorrectCoordinates(vertices, offX, offY);
      TGeoVolume* vol =
         gGeoManager->MakeArb8(vol_name.Data(), med, dz, vertices);
      vol->SetLineColor(med->GetMaterial()->GetDefaultColor());
      if (multi_layer) {
         /*** position absorber in mother ***/
         Double_t trans_z = -fTotalThickness / 2. + depth_in_det + dz; // (note: reference is CENTRE of absorber)
         TGeoTranslation* tr = new TGeoTranslation(offX, offY, trans_z);
         fDetVolume->AddNode(vol, no_abs, tr);
      }
      else {
         // single absorber: mother is absorber is detector is mother is ...
         fDetVolume = vol;
         fDetectorPosition = new TGeoTranslation(offX, offY, 0);
      }
 
      depth_in_det += thick;
      no_abs++;
 
      // front of next absorber is back of current absorber
      for (int i = 0; i < 4; i++) frontPlane[i] = backPlane[i];
      frontCentre = backCentre;
   }
}
 
 
 
 
void INDRAGeometryBuilder::CalculateCentre(TVector3* corners, TVector3& inter)
{
   // TVector3 corners[4] array of corner coordinates in usual order:
   //       corners[3] : theta-min, phi-min
   //       corners[2] : theta-max, phi-min
   //       corners[1] : theta-max, phi-max
   //       corners[0] : theta-min, phi-max
   //
   // Calculate the intersection of the normal to the plane passing through the origin.
 
   TVector3 a10 = corners[1] - corners[0];
   TVector3 a21 = corners[2] - corners[1];
   TVector3 normal = a10.Cross(a21);
   Double_t d = (corners[3] * normal) / normal.Mag2();
   inter = normal * d;
}
 
 
 
 
 
void INDRAGeometryBuilder::CalculateCornersInPlane(TVector3* plane, Double_t thetamin,
      Double_t thetamax, Double_t phimin, Double_t phimax, TVector3* corners)
{
   // fill TVector3 corners[4] array in usual order:
   //       corners[3] : theta-min, phi-min
   //       corners[2] : theta-max, phi-min
   //       corners[1] : theta-max, phi-max
   //       corners[0] : theta-min, phi-max
   // theta, phi in degrees
   // These points lie in the plane defined by the equivalent TVEctor3 plane[4] corners
 
   TVector3 a10 = plane[1] - plane[0];
   TVector3 a21 = plane[2] - plane[1];
   TVector3 normal = a10.Cross(a21);
   Double_t d0 = plane[3] * normal;
 
   TVector3 l;
   l.SetMagThetaPhi(1., thetamin * TMath::DegToRad(), phimin * TMath::DegToRad());
   corners[3] = l * (d0 / (l * normal));
   l.SetMagThetaPhi(1., thetamax * TMath::DegToRad(), phimin * TMath::DegToRad());
   corners[2] = l * (d0 / (l * normal));
   l.SetMagThetaPhi(1., thetamax * TMath::DegToRad(), phimax * TMath::DegToRad());
   corners[1] = l * (d0 / (l * normal));
   l.SetMagThetaPhi(1., thetamin * TMath::DegToRad(), phimax * TMath::DegToRad());
   corners[0] = l * (d0 / (l * normal));
}
 
 
 
 
std::vector<INDRAGeometryBuilder::theta_phi> INDRAGeometryBuilder::GetDirectionsToTestForImport() const
{
   std::vector<std::vector<theta_phi>> CSI_COORDS = {
      // Ring 1
      {  {2.40628, 30}, {2.40628, 60}, {2.40628, 90}, {2.40628, 120}, {2.40628, 150}, {2.40628, 180}, {2.40628, 210}, {2.40628, 240},
         {2.40628, 270}, {2.40628, 300}, {2.40628, 330}, {2.40628, 0}
      },
      // Ring 2
      {  {3.59843, 30}, {3.59843, 60}, {3.59843, 90}, {3.59843, 120},
         {3.59843, 150}, {3.59843, 180}, {3.59843, 210}, {3.59843, 240}, {3.59843, 270}, {3.59843, 300}, {3.59843, 330}, {3.59843, 0}
      },
      // Ring 3
      {  {5.6605, 20.427}, {5.6605, 39.573}, {5.6605, 50.427}, {5.6605, 69.573}, {5.6605, 80.427}, {5.6605, 99.573}, {5.6605, 110.427},
         {5.6605, 129.573}, {5.6605, 140.427}, {5.6605, 159.573}, {5.6605, 170.427}, {5.6605, 189.573}, {5.6605, 200.427}, {5.6605, 219.573},
         {5.6605, 230.427}, {5.6605, 249.573}, {5.6605, 260.427}, {5.6605, 279.573}, {5.6605, 290.427}, {5.6605, 309.573}, {5.6605, 320.427},
         {5.6605, 339.573}, {5.6605, 350.427}, {5.6605, 9.57299}
      },
      // Ring 4
      {  {8.23882, 20.2681}, {8.23882, 39.7319}, {8.23882, 50.2681}, {8.23882, 69.7319},
         {8.23882, 80.2681}, {8.23882, 99.7319}, {8.23882, 110.268}, {8.23882, 129.732}, {8.23882, 140.268}, {8.23882, 159.732},
         {8.23882, 170.268}, {8.23882, 189.732}, {8.23882, 200.268}, {8.23882, 219.732}, {8.23882, 230.268}, {8.23882, 249.732},
         {8.23882, 260.268}, {8.23882, 279.732}, {8.23882, 290.268}, {8.23882, 309.732}, {8.23882, 320.268}, {8.23882, 339.732},
         {8.23882, 350.268}, {8.23882, 9.73186}
      },
      // Ring 5
      {  {11.7348, 20.452}, {11.7348, 39.548}, {11.7348, 50.452}, {11.7348, 69.548}, {11.7348, 80.452},
         {11.7348, 99.548}, {11.7348, 110.452}, {11.7348, 129.548}, {11.7348, 140.452}, {11.7348, 159.548}, {11.7348, 170.452},
         {11.7348, 189.548}, {11.7348, 200.452}, {11.7348, 219.548}, {11.7348, 230.452}, {11.7348, 249.548}, {11.7348, 260.452},
         {11.7348, 279.548}, {11.7348, 290.452}, {11.7348, 309.548}, {11.7348, 320.452}, {11.7348, 339.548}, {11.7348, 350.452},
         {11.7348, 9.54796}
      },
      // Ring 6
      {  {16.4298, 20.1517}, {16.4298, 39.8483}, {16.4298, 50.1517}, {16.4298, 69.8483}, {16.4298, 80.1517},
         {16.4298, 99.8483}, {16.4298, 110.152}, {16.4298, 129.848}, {16.4298, 140.152}, {16.4298, 159.848}, {16.4298, 170.152},
         {16.4298, 189.848}, {16.4298, 200.152}, {16.4298, 219.848}, {16.4298, 230.152}, {16.4298, 249.848}, {16.4298, 260.152},
         {16.4298, 279.848}, {16.4298, 290.152}, {16.4298, 309.848}, {16.4298, 320.152}, {16.4298, 339.848}, {16.4298, 350.152},
         {16.4298, 9.84832}
      },
      // Ring 7
      {  {23.1125, 20.5305}, {23.1125, 39.4695}, {23.1125, 50.5305}, {23.1125, 69.4695}, {23.1125, 80.5305},
         {23.1125, 99.4695}, {23.1125, 110.531}, {23.1125, 129.469}, {23.1125, 140.531}, {23.1125, 159.469}, {23.1125, 170.531},
         {23.1125, 189.469}, {23.1125, 200.531}, {23.1125, 219.469}, {23.1125, 230.531}, {23.1125, 249.469}, {23.1125, 260.531},
         {23.1125, 279.469}, {23.1125, 290.531}, {23.1125, 309.469}, {23.1125, 320.531}, {23.1125, 339.469}, {23.1125, 350.531},
         {23.1125, 9.46947}
      },
      // Ring 8
      {  {30.1022, 20.2707}, {30.1022, 39.7293}, {30.1022, 50.2707}, {30.1022, 69.7293}, {30.1022, 80.2707},
         {30.1022, 99.7293}, {30.1022, 110.271}, {30.1022, 129.729}, {30.1022, 140.271}, {30.1022, 159.729}, {30.1022, 170.271},
         {30.1022, 189.729}, {30.1022, 200.271}, {30.1022, 219.729}, {30.1022, 230.271}, {30.1022, 249.729}, {30.1022, 260.271},
         {30.1022, 279.729}, {30.1022, 290.271}, {30.1022, 309.729}, {30.1022, 320.271}, {30.1022, 339.729}, {30.1022, 350.271},
         {30.1022, 9.72928}
      },
      // Ring 9
      {  {39.6964, 20.3843}, {39.6964, 39.6157}, {39.6964, 50.3843}, {39.6964, 69.6157}, {39.6964, 80.3843},
         {39.6964, 99.6157}, {39.6964, 110.384}, {39.6964, 129.616}, {39.6964, 140.384}, {39.6964, 159.616}, {39.6964, 170.384},
         {39.6964, 189.616}, {39.6964, 200.384}, {39.6964, 219.616}, {39.6964, 230.384}, {39.6964, 249.616}, {39.6964, 260.384},
         {39.6964, 279.616}, {39.6964, 290.384}, {39.6964, 309.616}, {39.6964, 320.384}, {39.6964, 339.616}, {39.6964, 350.384},
         {39.6964, 9.61571}
      },
      // Ring 10
      {  {49.4614, 20.6178}, {46., 39.3822}, // trajectory missing etalon
         {49.4614, 50.6178}, {49.4614, 69.3822}, {49.4614, 80.6178},
         {49.4614, 99.3822}, {49.4614, 110.618}, {49.4614, 129.382}, {49.4614, 140.618}, {49.4614, 159.382}, {49.4614, 170.618},
         {49.4614, 189.382}, {49.4614, 200.618}, {49.4614, 219.382}, {49.4614, 230.618}, {49.4614, 249.382}, {49.4614, 260.618},
         {49.4614, 279.382}, {49.4614, 290.618}, {49.4614, 309.382}, {49.4614, 320.618}, {49.4614, 339.382}, {49.4614, 350.618},
         {49.4614, 9.38224}
      },
      // Ring 11
      {  {63.6629, 20.7968}, {60., 39.2032}, // trajectory missing etalon
         {63.6629, 50.7968}, {63.6629, 69.2032}, {63.6629, 80.7968},
         {63.6629, 99.2032}, {63.6629, 110.797}, {63.6629, 129.203}, {63.6629, 140.797}, {63.6629, 159.203}, {63.6629, 170.797},
         {63.6629, 189.203}, {63.6629, 200.797}, {63.6629, 219.203}, {63.6629, 230.797}, {63.6629, 249.203}, {63.6629, 260.797},
         {63.6629, 279.203}, {63.6629, 290.797}, {63.6629, 309.203}, {63.6629, 320.797}, {63.6629, 339.203}, {63.6629, 350.797},
         {63.6629, 9.20323}
      },
      // Ring 12
      {  {79.188, 21.4427}, {75., 38.5573}, // trajectory missing etalon
         {79.188, 51.4427}, {79.188, 68.5573}, {79.188, 81.4427}, {79.188, 98.5573},
         {79.188, 111.443}, {79.188, 128.557}, {79.188, 141.443}, {79.188, 158.557}, {79.188, 171.443}, {79.188, 188.557}, {79.188, 201.443},
         {79.188, 218.557}, {79.188, 231.443}, {79.188, 248.557}, {79.188, 261.443}, {79.188, 278.557}, {79.188, 291.443}, {79.188, 308.557},
         {79.188, 321.443}, {79.188, 338.557}, {79.188, 351.443}, {79.188, 8.55734}
      },
      // Ring 13
      {  {101.234, 28.5099}, {101.703, 45}, {101.234, 61.4901}, {95., 73.5099}, // trajectory missing etalon
         {101.703, 90}, {101.234, 106.49}, {101.234, 118.51}, {101.703, 135}, {101.234, 151.49}, {101.234, 163.51},
         {101.703, 180}, {101.234, 196.49}, {101.234, 208.51}, {101.703, 225}, {101.234, 241.49}, {101.234, 253.51}, {101.703, 270},
         {101.234, 286.49}, {101.234, 298.51}, {101.703, 315}, {101.234, 331.49}, {101.234, 343.51}, {101.703, 0}, {101.234, 16.4901}
      },
      // Ring 14
      {  {118.722, 30.9923}, {118.722, 59.0077}, {118.722, 75.9923}, // trajectory missing etalon
         {118.722, 104.008}, {118.722, 120.992}, {118.722, 149.008},
         {118.722, 165.992}, {118.722, 194.008}, {118.722, 210.992}, {118.722, 239.008}, {118.722, 255.992}, {118.722, 284.008},
         {118.722, 300.992}, {118.722, 329.008}, {118.722, 345.992}, {118.722, 14.0077}
      },
      // Ring 15
      {  {136.289, 30.5072}, {136.289, 59.4928}, {136.289, 75.5072}, // trajectory missing etalon
         {136.289, 104.493}, {136.289, 120.507}, {136.289, 149.493}, {136.289, 165.507}, {136.289, 194.493},
         {136.289, 210.507}, {136.289, 239.493}, {136.289, 255.507}, {136.289, 284.493}, {136.289, 300.507}, {136.289, 329.493},
         {136.289, 345.507}, {136.289, 14.4928}
      },
      // Ring 16
      {  {150.75, 45}, {150.75, 80},  // trajectory missing etalon
         {150.75, 135}, {150.75, 180}, {150.75, 225}, {150.75, 270}, {150.75, 315}, {150.75, 0}
      },
      // Ring 17
      {  {168.561, 45}, {174., 90}, // trajectory missing etalon
         {168.561, 135}, {168.561, 180}, {168.561, 225}, {168.561, 270}, {168.561, 315}, {168.561, 0.}
      }
   };
   std::vector<theta_phi> ETALON_COORDS = {
      {50.389351, 36.335958}, // ring 10
      {63.6629, 39.2032}, // ring 11
      {79.188, 38.5573}, // ring 12
      {101.234, 73.5099}, // ring 13
      {117.825, 83.5543}, // ring 14
      {133.642, 82.1937}, // ring 15
      {149.375, 90}, // ring 16
      {164.703, 90} // ring 17
   };
   // small utility functions to access arrays, taking into account the fact that ring 2
   // modules are numbered 1, 3, 5, ...
   auto get_direction = [&](int ring, int module) {
      if (ring == 2) return CSI_COORDS[1][module / 2];
      return CSI_COORDS[ring - 1][module - 1];
   };
   auto get_etalon_direction = [&](int ring) {
      return ETALON_COORDS[ring - 10];
   };
   // now scan map of present detectors and add corresponding directions to vector for KVGeoImport
   std::vector<theta_phi> directions_for_import;
 
   using ring_t = std::unordered_map<int, int>;
   std::unordered_map<int, ring_t> ring_module;
   for (auto& p : detector_map) {
      if (p.second.module > 0) { //avoid etalon telescopes with module=0
         auto ring = p.second.ring;
         auto module = p.second.module;
         if (!ring_module[ring][module]) { // only add one direction for each ring/module pair
            directions_for_import.push_back(get_direction(ring, module));
            ++ring_module[ring][module];
         }
      }
      else { // etalon telescopes
         auto ring = p.second.ring;
         directions_for_import.push_back(get_etalon_direction(ring));
      }
   }
   return directions_for_import;
}
 
 
 
 
Bool_t INDRAGeometryBuilder::CheckImportResults(const KVSeqCollection* detlist)
{
   // Make sure that all detectors in detector map are in the list of detectors
   // imported into the multidet array
    //
    // We also set the (simplified geometry) thetamin/max & phimin/max for each detector.
    // This is used e.g. in the GUI for visualising scalers during experiments.
 
   Info("CheckImportResults", "Checking results of import of %lu expected detectors",
        detector_map.size());
   bool ok = true;
   unsigned long found = 0;
   for (auto& det : detector_map) {
       auto _det = detlist->FindObject(det.first.c_str());
      if (!_det) {
         Error("CheckImportResults", "Detector %s was not imported",
               det.first.c_str());
         ok = false;
      }
      else {
         ++found;
          dynamic_cast<KVDetector*>(_det)->SetPhiMinMax(det.second.phiminmax.first,det.second.phiminmax.second);
          dynamic_cast<KVDetector*>(_det)->SetPolarMinMax(det.second.thetaminmax.first,det.second.thetaminmax.second);
      }
   }
   if (found < detector_map.size())
      Warning("CheckImportResults", "Found %lu imported detectors (%lu missing)",
              found, detector_map.size() - found);
   else
      Info("CheckImportResults", "All %lu expected detectors were imported", found);
   return ok;
}
 
 
 
 
void INDRAGeometryBuilder::CloseAndDraw()
{
   gGeoManager->CloseGeometry();
   //gGeoManager->GetTopVolume()->SetVisContainers();
   gGeoManager->DefaultColors();
   gGeoManager->GetTopVolume()->Draw("ogl");
}
 
 
 
 
void INDRAGeometryBuilder::ReflectPad(TVector3* orig, Double_t phicentre, TVector3* newpad)
{
   for (int i = 0; i < 4; i++) {
      newpad[3 - i].SetMagThetaPhi(orig[i].Mag(), orig[i].Theta(), 2.*phicentre - orig[i].Phi());
   }
}
 
 
 
 
void INDRAGeometryBuilder::BuildTarget()
{
   //TO DO: multi-layer targets + target with angle to beam
 
   // target holder
   KVMaterial target_holder_mat("Al");
   new TGeoBBox("TARGET_FRAME", 3., 3., 0.1 / 2.);
   new TGeoEltu("TARGET_HOLE", 2., 2., 0.1 / 2.);
   TGeoCompositeShape* cs = new TGeoCompositeShape("TARGET_FRAME", "TARGET_FRAME - TARGET_HOLE");
   TGeoVolume* target_frame = new TGeoVolume("TARGET_FRAME", cs, target_holder_mat.GetGeoMedium());
   gGeoManager->GetTopVolume()->AddNode(target_frame, 1);
 
   // add target
   // only single-layer target treated correctly
   KVTarget* T = gIndra->GetTarget();
   if (T) {
      KVMaterial* targMat = (KVMaterial*)T->GetLayers()->First();
      TGeoVolume* target = gGeoManager->MakeEltu("TARGET", targMat->GetGeoMedium(), 2., 2., targMat->GetThickness() / 2.);
      gGeoManager->GetTopVolume()->AddNode(target, 1);
   }
}
 
 
 
 
void INDRAGeometryBuilder::Build(Bool_t withTarget, Bool_t closeGeometry)
{
   if (!gGeoManager) {
      Error("Build", "You must define gGeoManager with KVMultiDetArray::CreateGeoManager before calling this method");
      return;
   }
 
   // read infos on present/absent detectors & nominal thicknesses/pressures
   read_indra_struct_file();
 
   // read specific dataset detector thicknesses
   gIndra->SetDetectorThicknesses();
 
   if (fWithPhoswich) MakeRing("PHOS", 1);
   else {
      MakeRing("SI", 1);
      MakeRing("CSI", 1);
   }
   for (int ring = 2; ring <= 16; ring += 2) {
      MakeRing("CI", ring);
      if (ring > 1 && ring < 10)MakeRing("SI", ring);
      MakeRing("CSI", ring);
      if (ring == 12) {
         MakeRing("CI", ring + 1);
         MakeRing("CSI", ring + 1);
      }
   }
   for (int ring = 10; ring <= 17; ring++) {
      if (CheckDetectorPresent(Form("SI75_%d", ring))) MakeEtalon(ring);
   }
   if (withTarget) BuildTarget();
   if (closeGeometry) {
      gGeoManager->DefaultColors();
      gGeoManager->CloseGeometry();
   }
 
}
 
 
 
void INDRAGeometryBuilder::Build(const KVNumberList& rings, const KVNameValueList& detectors)
{
   // Build geometry, using only the ring numbers in the KVNumberList, and
   // only the detectors which are in the KVNameValueList:
   //
   // KVNameValueList dets;
   // dets.SetValue("CI",1);
   // dets.SetValue("SI",1);
   // igb.Build("2-9", dets); // use INDRAGeometryBuilder object 'igb' to build
   //                         // ChIo and Silicon detectors of rings 2-9
   //
   // Possible detector types are "CI", "SI", "CSI", "PHOS"
   // If you want to see the target in the previous example, use dets.SetValue("TARGET",1)
   //
   // Include "ring+100" in number list to build etalon telescope for "ring"
   //
   // This method will not destroy any pre-existing geometry (gGeoManager),
   // but will create one if none exists.
   // Also it will not close the geometry, leaving the possibility to add ancillary
   // detectors.
 
   if (!gIndra) {
      Error("Build", "You must build the geometry with gDataSet->BuildMultiDetector() before calling this method");
      return;
   }
 
   if (!gGeoManager) {
      Error("Build", "You must defined gGeoManager with KVMultiDetArray::CreateGeoManager before calling this method");
      return;
   }
 
   if (rings.Contains(1)) {
      if (detectors.HasParameter("PHOS")) MakeRing("PHOS", 1);
      if (detectors.HasParameter("SI"))MakeRing("SI", 1);
      if (detectors.HasParameter("CSI"))MakeRing("CSI", 1);
   }
   for (int ring = 2; ring <= 16; ring += 2) {
      if (rings.Contains(ring) && detectors.HasParameter("CI")) MakeRing("CI", ring);
      if (ring > 1 && ring < 10) {
         if (rings.Contains(ring) && detectors.HasParameter("SI"))MakeRing("SI", ring);
      }
      if (rings.Contains(ring) && detectors.HasParameter("CSI"))MakeRing("CSI", ring);
      if (ring == 12) {
         if (rings.Contains(ring + 1) && detectors.HasParameter("CI"))MakeRing("CI", ring + 1);
         if (rings.Contains(ring + 1) && detectors.HasParameter("CSI"))MakeRing("CSI", ring + 1);
      }
   }
   for (int ring = 10; ring <= 17; ring++) {
      if (rings.Contains(ring + 100)) MakeEtalon(ring);
   }
 
   if (detectors.HasParameter("TARGET"))
      BuildTarget();
}
 
 
 
 
void INDRAGeometryBuilder::BuildEtalonVolumes()
{
   KVMaterial mat_si("Si");
   KVMaterial mat_li("Li");
   KVMaterial mat_al("Al");
   TGeoMedium* Silicon = mat_si.GetGeoMedium();
   TGeoMedium* Lithium = mat_li.GetGeoMedium();
   TGeoMedium* Alu = mat_al.GetGeoMedium();
 
   Double_t sili_diameter_total = 2.54;
   Double_t holder_thickness = 0.05;
   Double_t holder_length = 1.0;
   Double_t sili_diameter_active = 2.36;
   Double_t sili_silicon_thickness = 0.22;
   Double_t sili_lithium_thickness = 0.0044;
   Double_t si75_thickness = 0.008;
   Double_t si75_diameter_active = 2.2;
 
   fEtalonVol = new TGeoVolumeAssembly("STRUCT_ETALON");
   fEtalonVol->SetMedium(gGeoManager->GetMedium("Vacuum"));
 
   TGeoVolume* holder = gGeoManager->MakeTube("DEADZONE_ETALON_HOLDER", Alu,
                        (sili_diameter_total) / 2., (sili_diameter_total + 2 * holder_thickness) / 2., holder_length / 2.);
   holder->SetLineColor(holder->GetMaterial()->GetDefaultColor());
   fEtalonVol->AddNode(holder, 1);
 
   Double_t w = sili_silicon_thickness + sili_lithium_thickness;
   TGeoVolume* sili_dz = gGeoManager->MakeTube("DEADZONE_SILI", Alu,
                         sili_diameter_active / 2., (sili_diameter_total) / 2., w / 2.);
   sili_dz->SetLineColor(sili_dz->GetMaterial()->GetDefaultColor());
 
   fEtalonVol->AddNode(sili_dz, 1, new TGeoTranslation(0, 0, holder_length / 4.));
 
   TGeoVolumeAssembly* sili = new TGeoVolumeAssembly("DET_SILI");
   sili->SetMedium(gGeoManager->GetMedium("Vacuum"));
   TGeoVolume* sili_si = gGeoManager->MakeTube("ACTIVE_SILI", Silicon,
                         0., (sili_diameter_active) / 2., (sili_silicon_thickness) / 2.);
   TGeoVolume* sili_li = gGeoManager->MakeTube("SILI_LI", Lithium,
                         0., (sili_diameter_active) / 2., (sili_lithium_thickness) / 2.);
   sili_si->SetLineColor(sili_si->GetMaterial()->GetDefaultColor());
   sili_li->SetLineColor(sili_li->GetMaterial()->GetDefaultColor());
   sili->AddNode(sili_si, 1, new TGeoTranslation(0, 0, -(w / 2. - sili_silicon_thickness / 2.)));
   sili->AddNode(sili_li, 1, new TGeoTranslation(0, 0, w / 2. - sili_lithium_thickness / 2.));
   fEtalonVol->AddNode(sili, 1, new TGeoTranslation(0, 0, holder_length / 4.));
 
   TGeoVolume* si75 = gGeoManager->MakeTube("DET_SI75", Silicon,
                      0., (si75_diameter_active) / 2., (si75_thickness) / 2.);
   TGeoVolume* si75_dz = gGeoManager->MakeTube("DEADZONE_SI75", Alu,
                         (si75_diameter_active) / 2., (sili_diameter_total) / 2., (si75_thickness) / 2.);
 
   fEtalonVol->AddNode(si75_dz, 1);
   fEtalonVol->AddNode(si75, 1);
 
}
 
 
 
 
void INDRAGeometryBuilder::MakeEtalon(int RING)
{
   // Build and add etalon telescope for ring
   //Double_t theta[] = {51.075000,63.340000,79.435000,100.685000,118.235000,133.905000,149.790000,166.435000};
   //Double_t phi[] = {37.500000,37.500000,37.500000,90.000000,78.750000,78.750000,90.000000,90.000000};
   Double_t dist[] = {17.4005, 17.4005, 17.4005, 16.7005, 16.7005, 16.7005, 17.4005, 17.4005,};
 
   if (!fEtalonVol) BuildEtalonVolumes();
 
   cout << "Etalons " << RING << " : " << fEtalonTheta[RING - 10] << "  " << fEtalonPhi[RING - 10] << endl;
   TGeoTranslation trans;
   trans.SetDz(dist[RING - 10]);
   TGeoRotation rot1, rot2;
   TGeoHMatrix h;
   TGeoHMatrix* ph = 0;
 
//    rot2.SetAngles(phi[RING-10]+90., theta[RING-10], 0.);
//    rot1.SetAngles(-1.*phi[RING-10], 0., 0.);
   rot2.SetAngles(fEtalonPhi[RING - 10] + 90., fEtalonTheta[RING - 10], 0.);
   rot1.SetAngles(-1.*fEtalonPhi[RING - 10], 0., 0.);
//    if(RING==13){
//       TGeoTranslation p(4.5,0,0);
//       h = p * rot2 * trans * rot1;
//    }
   if (RING == 10) {
      TGeoTranslation p(1.5, 2.5, 0);
      h = rot2 * trans * p * rot1;
   }
   else if (RING == 11) {
      TGeoTranslation p(1.5, -1.5, 0);
      h = rot2 * trans * p * rot1;
   }
   else if (RING == 12) {
      TGeoTranslation p(1.5, 0, 0);
      h = rot2 * trans * p * rot1;
   }
   else if (RING == 13) {
      TGeoTranslation p(-4, -0.5, 0);
      h = rot2 * trans * p * rot1;
   }
   else if (RING == 14) {
      TGeoTranslation p(-1.7, 1.75, 0);
      h = rot2 * trans * p * rot1;
   }
   else if (RING == 15) {
      TGeoTranslation p(-1.7, -3, 0);
      h = rot2 * trans * p * rot1;
   }
   else if (RING == 16) {
      TGeoTranslation p(0, 2.5, 0);
      h = rot2 * trans * p * rot1;
   }
   else {
      TGeoTranslation p(0, -2.25, 0);
      h = rot2 * trans * p * rot1;
   }
   ph = new TGeoHMatrix(h);
   gGeoManager->GetTopVolume()->AddNode(fEtalonVol, RING, ph);
}