KVFAZIA3D.cpp

//Author: gruyer,,,
 
#include "KVFAZIA3D.h"
#include "KVUnits.h"
#include "KVFAZIA3DBlock.h"
#include <map>
 
#include <KVDataSet.h>
 
ClassImp(KVFAZIA3D)
 
// BEGIN_HTML <!--
/* -->
<h2>KVFAZIA3D</h2>
<h4>description of the FAZIA-12B demonstrator</h4>
<!-- */
// --> END_HTML
 
 
 
KVFAZIA3D::KVFAZIA3D()
   : KVFAZIA()
{
   // Default constructor
   SetTitle(ClassName());
}
 
 
 
 
void KVFAZIA3D::MakeConfigFile(const char* filename)
{
// Read a csv file with the quartet coordinates to generate the FAZIA3D.Quartet.Position.env file
// used to build the FAZIA detector.
//
// the csv file structure should be (no header):
// quartet_index,corner_index,quartet_index*100+corner_index,X,Y(Z),Z
//
// example for the first quartet of Block 0:
//    1,11,111,-66.8,887.2,114.3
//    1,12,112,-46.8,887.9,114.2
//    1,13,113,-47.0,889.9,94.3
//    1,14,114,-67.0,889.2,94.4
//    1,21,121,-44.6,888.0,114.2
//    1,22,122,-24.6,888.7,114.0
//    1,23,123,-24.8,890.7,94.2
//    1,24,124,-44.8,890.0,94.3
//    1,31,131,-44.8,890.2,92.1
//    1,32,132,-24.8,890.9,92.0
//    1,33,133,-25.0,893.0,72.1
//    1,34,134,-45.0,892.3,72.2
//    1,41,141,-67.0,889.5,92.2
//    1,42,142,-47.0,890.2,92.1
//    1,43,143,-47.2,892.2,72.2
//    1,44,144,-67.2,891.5,72.3
// example for the last quartet of Block 11:
//    114,11,11411,-207.0,876.4,-19.4
//    114,12,11412,-187.3,880.0,-19.6
//    114,13,11413,-187.4,879.3,-39.6
//    114,14,11414,-207.0,875.7,-39.4
//    114,21,11421,-185.2,880.4,-19.6
//    114,22,11422,-165.5,884.1,-19.8
//    114,23,11423,-165.5,883.3,-39.8
//    114,24,11424,-185.2,879.7,-39.6
//    114,31,11431,-185.2,879.6,-41.8
//    114,32,11432,-165.5,883.2,-42.0
//    114,33,11433,-165.6,882.5,-62.0
//    114,34,11434,-185.3,878.8,-61.8
//    114,41,11441,-207.0,875.6,-41.6
//    114,42,11442,-187.4,879.2,-41.8
//    114,43,11443,-187.4,878.4,-61.8
//    114,44,11444,-207.1,874.8,-61.6
//
// The output file contains, for each quartet, the center coordinates (X,Y,Z)
// and the three rotation angles (Psi,Theta,Phi)
 
   KVFileReader rd;
   rd.OpenFileToRead(filename);
 
   KVEnv fout;
   fout.AddCommentLine("File containing the position and orientation of the FAZIA quartets");
   fout.AddCommentLine("as measured by the GANIL geometers");
   fout.AddCommentLine(Form("Created from '%s' using KVFAZIA3D::MakeConfigFile()", filename));
 
   std::map<int, TVector3> pos;
 
   while (rd.IsOK()) {
      rd.ReadLine(",");
      if (!rd.GetNparRead()) continue;
      int idp = rd.GetReadPar(2).Atoi();
      int idt = idp / 10;
      int idq = idt / 10;
      pos[idp] = {rd.GetReadPar(5).Atof() / 10, rd.GetReadPar(3).Atof() / 10, rd.GetReadPar(4).Atof() / 10}; //-89.74
   }
 
   for (int bb = 0; bb < 12; bb++) {
      for (auto qq : {
               1, 2, 3, 4
            }) {
         int id = 10 * bb + qq;
 
         // vector to the external corners of the silicon detectors
         TVector3 v1 = pos[id * 100 + 11];
         TVector3 v2 = pos[id * 100 + 22];
         TVector3 v3 = pos[id * 100 + 33];
         TVector3 v4 = pos[id * 100 + 44];
 
         // vector to the center of the quartet for translation
         TVector3 center = v1 + v2 + v3 + v4;
         center *= .25;
         fout.SetValue(Form("B%03d-Q%d.X", bb, qq), center.X());
         fout.SetValue(Form("B%03d-Q%d.Y", bb, qq), center.Y());
         fout.SetValue(Form("B%03d-Q%d.Z", bb, qq), center.Z());
 
         // bring the quartet center at (0,0,0)
         v1 -= center;
         v2 -= center;
         v3 -= center;
         v4 -= center;
 
         TVector3 e1 = v2 - v1;
         TVector3 e2 = v3 - v2;
 
         // vector normal to the quartet surface for rotation
         TVector3 normal = e1.Cross(e2).Unit();
 
         double theta = normal.Theta() * TMath::RadToDeg();
         double phi = normal.Phi() * TMath::RadToDeg();
 
         // bring the quartet in the XY plane
         TRotation rr;
         rr.RotateZ(-normal.Phi());
         TRotation rr2;
         rr2.RotateY(-normal.Theta());
 
         v1 = rr * v1;
         v2 = rr * v2;
         v1 = rr2 * v1;
         v2 = rr2 * v2;
 
         // recompute the upper edge of the quartet to find the final rotation
         e1 = v2 - v1;
         TVector3 oy(0, 1, 0);
 
         // compute the final rotation angle
         double psi = TMath::Sign(oy.Angle(e1), e1.X()) * TMath::RadToDeg();
         fout.SetValue(Form("B%03d-Q%d.Psi", bb, qq), -psi);
         fout.SetValue(Form("B%03d-Q%d.Theta", bb, qq), theta);
         fout.SetValue(Form("B%03d-Q%d.Phi", bb, qq), phi);
      }
   }
   fout.WriteFile("FAZIA3D.Quartet.Position.env");
 
}
 
 
 
 
void KVFAZIA3D::GetGeometryParameters()
{
   //defined configuration of the blocks
   fNblocks = KVBase::GetDataSetEnv(fDataSet, "FAZIA.NBlocks", 12.0);
   fFDist = KVBase::GetDataSetEnv(fDataSet, "FAZIA.DistanceZ", 100.0);
   SetGeometryImportParameters(.25, 1, 1, 0, 14, 360,
                               KVBase::GetDataSetEnv(fDataSet, "FAZIA.DistanceX", 0.0),
                               KVBase::GetDataSetEnv(fDataSet, "FAZIA.DistanceY", 0.0),
                               fFDist - 100.0);
}
 
 
 
 
 
 
void KVFAZIA3D::BuildFAZIA()
{
   // Build the FAZIA demonstrator geometry.
   //
   // Distance from target and minimum polar angle between innermost blocks
   // can be changed for different datasets using:
   //
   // [dataset].FAZIA.DistanceZ:  95.4
   // [dataset].FAZIA.DistanceX:  1.
   // [dataset].FAZIA.DistanceY:  0.
   //
   // Number of blocks in demonstrator, their numbers, and positions can be modified for
   // different datasets using the following environment variables:
   //
   // [dataset].FAZIA.NBlocks: 12
   // [dataset].FAZIA.BlockNumbers: 0,1,2,3,4,5,6,7,8,9,10,11
   // [dataset].FAZIA.BlockPositions: 0,1,2,3,4,5,6,7,8,9,10,11
 
   Info("BuildFAZIA", "Exact geometry using 'FAZIA3D.Quartet.Position.env', %f cm from target.", fFDist);
 
   TGeoVolume* top = gGeoManager->GetTopVolume();
 
   // block numbering and positions
   auto block_numbers = KVBase::GetDataSetEnv<KVString>(fDataSet, "FAZIA.BlockNumbers", "0,1,2,3,4,5,6,7,8,9,10,11");
   block_numbers.Begin(",");
   auto block_positions = KVBase::GetDataSetEnv<KVString>(fDataSet, "FAZIA.BlockPositions", "0,1,2,3,4,5,6,7,8,9,10,11");
   block_positions.Begin(",");
 
   for (Int_t i = 0; i < fNblocks; i += 1) {
 
      Int_t block_number = block_numbers.Next().Atoi();
      Int_t bb = block_positions.Next().Atoi();
 
      KVFAZIA3DBlock* block = new KVFAZIA3DBlock(block_number);
      block->Build();
 
      TGeoTranslation final_trans;
      final_trans.SetDz(fFDist * KVUnits::cm);
      if (KVBase::GetDataSetEnv(fDataSet, "FAZIETO.DistanceX", 0.0) != 0.) final_trans.SetDx(KVBase::GetDataSetEnv(fDataSet, "FAZIETO.DistanceX", 0.0));
      if (KVBase::GetDataSetEnv(fDataSet, "FAZIETO.DistanceY", 0.0) != 0.) final_trans.SetDy(KVBase::GetDataSetEnv(fDataSet, "FAZIETO.DistanceY", 0.0));
 
      TGeoHMatrix* ph = new TGeoHMatrix(final_trans);
      top->AddNode(block, block_number, ph);
   }
 
}