EXPERT ROOT

/*

Neurad prototype for Wuppertal tests. With one 64-pixel MAPMT H12700  

Todo SB: create or find material for black paper, currently polypropilen is used!

*/

#include <iomanip>

#include <iostream>

#include "TGeoManager.h"

#include "TMath.h"

void create_NeuRad_Wupper_Proto_geo()

{

  // ---------------  INIT ----------------------------------------------------

  // Create a global translation

  Float_t trans_X = 0,trans_Y = 0, trans_Z = 50.; 

  // Create a zero rotation

  TGeoRotation *fZeroRotation = new TGeoRotation();

  fZeroRotation->RotateX(0.);

  fZeroRotation->RotateY(0.);

  fZeroRotation->RotateZ(0.);

  Double_t fiber_X = 0.3; //cm

  Double_t fiber_Y = 0.3; //cm

  Double_t fiber_Z = 25.; //cm

  Double_t dead_cladding_thick = 0.005; //cm dead layer between fibers (cladidng and painting)

  Double_t submodule_wrapping = 0.01; //cm Teflon wrapping around 8x8 fibers assembly

  Double_t housing_Thick=0.01; //cm Paper wrapping (housing) around module

  Double_t ersatz_pmt_Z = 0.24; //cm simulation of the overall pmt material

  Double_t cover_Z = 0.1; //cm simulation of the plastic cover at the opposite (w.r.t PMT) end 

  Int_t fibers_in_pixel_X = 2;

  Int_t fibers_in_pixel_Y = 2;

  Int_t pixels_in_submodule_X = 4;

  Int_t pixels_in_submodule_Y = 4;

  Int_t submodules_in_module_X = 2;

  Int_t submodules_in_module_Y = 2;

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

  TGeoManager*   gGeoMan = NULL;

  // -------   Load media from media file   -----------------------------------

  FairGeoLoader*    geoLoad = new FairGeoLoader("TGeo","FairGeoLoader");

  FairGeoInterface* geoFace = geoLoad->getGeoInterface();

  TString geoPath = gSystem->Getenv("VMCWORKDIR");

  TString medFile = geoPath + "/geometry/media.geo";

  geoFace->setMediaFile(medFile);

  geoFace->readMedia();

  gGeoMan = gGeoManager;

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

  // -------   Geometry file name (output)   ----------------------------------

  TString geoFileName = geoPath + "/geometry/NeuRad.v4.geo.root";

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

  // -----------------   Get and create the required media    -----------------

  FairGeoMedia*   geoMedia = geoFace->getMedia();

  FairGeoBuilder* geoBuild = geoLoad->getGeoBuilder();

  FairGeoMedium* mBC408      = geoMedia->getMedium("BC408");

  if ( ! mBC408 ) Fatal("Main", "FairMedium BC408 not found");

  geoBuild->createMedium(mBC408);

  TGeoMedium* pBC408 = gGeoMan->GetMedium("BC408");

  if ( ! pBC408 ) Fatal("Main", "Medium BC408 not found");

   FairGeoMedium* mTefl      = geoMedia->getMedium("Tefflon");

  if ( ! mTefl ) Fatal("Main", "FairMedium Tefflon not found");

  geoBuild->createMedium(mTefl);

  TGeoMedium* pTefl = gGeoMan->GetMedium("Tefflon");

  if ( ! pTefl ) Fatal("Main", "Medium Tefflon not found");

   FairGeoMedium* mPolypr      = geoMedia->getMedium("polypropylene");

  if ( ! mPolypr ) Fatal("Main", "FairMedium polypropylene not found");

  geoBuild->createMedium(mPolypr);

  TGeoMedium* pPolypr = gGeoMan->GetMedium("polypropylene");

  if ( ! pPolypr ) Fatal("Main", "Medium polypropylene not found");

  FairGeoMedium* mAl      = geoMedia->getMedium("aluminium");

  if ( ! mAl ) Fatal("Main", "FairMedium aluminium not found");

  geoBuild->createMedium(mAl);

  TGeoMedium* pAl = gGeoMan->GetMedium("aluminium");

  if ( ! pAl ) Fatal("Main", "Medium aluminium not found");

  FairGeoMedium* vacuum      = geoMedia->getMedium("vacuum");

  if ( ! vacuum ) Fatal("Main", "FairMedium vacuum not found");

  geoBuild->createMedium(vacuum);

  TGeoMedium* pMed0 = gGeoMan->GetMedium("vacuum");

  if ( ! pMed0 ) Fatal("Main", "Medium vacuum not found");

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

  //------------------------- VOLUMES -----------------------------------------

  // --------------   Create geometry and top volume  -------------------------

  gGeoMan = (TGeoManager*)gROOT->FindObject("FAIRGeom");

  gGeoMan->SetName("DETgeom");

  TGeoVolume* top = new TGeoVolumeAssembly("TOP");

  gGeoMan->SetTopVolume(top);

  TGeoVolume* NeuRad = new TGeoVolumeAssembly("NeuRad");

  //------------------ BC408  fiber  -----------------------------------------

  TGeoVolume* fiber = gGeoManager->MakeBox("fiber", pBC408, fiber_X/2.-dead_cladding_thick, fiber_Y/2.-dead_cladding_thick, fiber_Z/2.);

  TGeoVolume* dead_cladding = gGeoManager->MakeBox("dead_cladding", pBC408, fiber_X/2., fiber_Y/2., fiber_Z/2.);  

  //------------------  pixel  vacuum box -----------------------------------------

  Double_t pixel_X = fiber_X *fibers_in_pixel_X;

  Double_t pixel_Y = fiber_Y *fibers_in_pixel_Y;

  TGeoVolume* pixel = gGeoManager->MakeBox("pixel", pMed0, pixel_X/2., pixel_Y/2., fiber_Z/2.);

  //-------------------------submodule-------------------

  Double_t submodule_X = pixel_X*pixels_in_submodule_X + 2.*submodule_wrapping;

  Double_t submodule_Y = pixel_Y*pixels_in_submodule_Y + 2.*submodule_wrapping;

  TGeoVolume* submodule = gGeoManager->MakeBox("submodule", pTefl, submodule_X/2., submodule_Y/2., fiber_Z/2.);

  //------------------ volume assmbly module  -----------------------------------------

  Double_t module_X = submodule_X*submodules_in_module_X+2.*housing_Thick;

  Double_t module_Y = submodule_Y*submodules_in_module_Y+2.*housing_Thick;

  // Double_t module_Z = fiber_Z+2.*(ersatz_pmt_Z+cover_Z);

 // TGeoVolume* module =  gGeoManager->MakeBox("module", pTefl, module_X/2., module_Y/2., module_Z/2.);

  TGeoVolume* module = new TGeoVolumeAssembly("module");

   //------------------  ersatz PMT, paper wrapping (module housing) and cover -----------------------------------------

  TGeoVolume* ersatz_pmt = gGeoManager->MakeBox("ersatz_pmt", pAl, module_X/2., module_Y/2., ersatz_pmt_Z/2.);

  TGeoVolume* cover = gGeoManager->MakeBox("cover", pPolypr, module_X/2., module_Y/2., cover_Z/2.);

  TGeoBBox* sh_module_housing_out = new TGeoBBox("sh_module_housing_out", module_X/2.,  module_Y/2., fiber_Z/2.);

  TGeoBBox* sh_module_housing_in = new TGeoBBox("sh_module_housing_in", module_X/2.-housing_Thick, module_Y/2.-housing_Thick, fiber_Z/2.);

  TGeoCompositeShape* sh_module_housing = new TGeoCompositeShape("sh_module_housing","sh_module_housing_out-sh_module_housing_in");

  TGeoVolume* module_housing = new TGeoVolume("module_housing",sh_module_housing,pPolypr);

  //------------------ STRUCTURE  -----------------------------------------

  dead_cladding->AddNode(fiber, 0, new TGeoCombiTrans(.0,.0,.0, 

                                                  fZeroRotation));

  //------------------ Add fibers to pixel  -----------------------------

  Int_t i_fiber = 0;

  for (Int_t i_Y_fiber = 0; i_Y_fiber < fibers_in_pixel_Y; i_Y_fiber++){

    for (Int_t i_X_fiber = 0; i_X_fiber < fibers_in_pixel_X; i_X_fiber++){

      Double_t fiber_in_pixel_X_trans = -pixel_X/2. + fiber_X/2. + fiber_X * i_X_fiber;

      Double_t fiber_in_pixel_Y_trans = -pixel_Y/2. + fiber_Y/2. + fiber_Y * i_Y_fiber; 

      Double_t fiber_in_pixel_Z_trans = 0.;

      pixel->AddNode( dead_cladding, i_fiber, new TGeoCombiTrans(fiber_in_pixel_X_trans, 

                                                            fiber_in_pixel_Y_trans,

                                                            fiber_in_pixel_Z_trans, 

                                                            fZeroRotation));

      i_fiber++;

    }

  }

  //----------- Add pixels to submodule

  Int_t i_pixel = 0;

  for (Int_t i_Y_pixel = 0; i_Y_pixel < pixels_in_submodule_Y; i_Y_pixel++){

    for (Int_t i_X_pixel = 0; i_X_pixel < pixels_in_submodule_X; i_X_pixel++){

      //Double_t pixel_in_submodule_X_trans = pixel_X*(pixels_in_submodule_X/2.-i_X_pixel-0.5);

      //Double_t pixel_in_submodule_Y_trans = pixel_Y*(pixels_in_submodule_Y/2.-i_Y_pixel-0.5);

      Double_t pixel_in_submodule_X_trans = - submodule_X/2 + submodule_wrapping + pixel_X/2 + pixel_X * i_X_pixel;

      Double_t pixel_in_submodule_Y_trans = - submodule_Y/2 + submodule_wrapping + pixel_Y/2 + pixel_Y * i_Y_pixel;

      Double_t pixel_in_submodule_Z_trans = 0.;

      submodule->AddNode( pixel, i_pixel, new TGeoCombiTrans(pixel_in_submodule_X_trans, 

                                                            pixel_in_submodule_Y_trans,

                                                            pixel_in_submodule_Z_trans, 

                                                            fZeroRotation));

      i_pixel++;

    }

  }

// ------------Add submodules to module

Int_t i_submodule = 0;

  for (Int_t i_Y_submodule = 0; i_Y_submodule < submodules_in_module_Y; i_Y_submodule++){

    for (Int_t i_X_submodule = 0; i_X_submodule <submodules_in_module_X; i_X_submodule++){

      //Double_t submodule_in_module_X_trans = submodule_X*(submodules_in_module_X/2.-i_X_submodule-0.5);

      //Double_t submodule_in_module_Y_trans = submodule_Y*(submodules_in_module_Y/2.-i_Y_submodule-0.5);

      Double_t submodule_in_module_X_trans = - module_X/2 + housing_Thick + submodule_X/2 + submodule_X * i_X_submodule;

      Double_t submodule_in_module_Y_trans = - module_Y/2 + housing_Thick + submodule_Y/2 + submodule_Y * i_Y_submodule;

      Double_t submodule_in_module_Z_trans = 0.;

      module->AddNode( submodule, i_submodule, new TGeoCombiTrans(submodule_in_module_X_trans, 

                                                            submodule_in_module_Y_trans,

                                                            submodule_in_module_Z_trans, 

                                                            fZeroRotation));

      i_submodule++;

    }

  }

//------------------ Add ersatz_pmt, paper wrapping (module housing)and cover to module -----------

  Double_t ersatz_pmt_in_module_X_trans = 0.;

  Double_t ersatz_pmt_in_module_Y_trans = 0.;

  Double_t ersatz_pmt_in_module_Z_trans = fiber_Z/2. + ersatz_pmt_Z/2.;

  module->AddNode(ersatz_pmt, 0, new TGeoCombiTrans(ersatz_pmt_in_module_X_trans, 

                                                  ersatz_pmt_in_module_Y_trans, 

                                                  ersatz_pmt_in_module_Z_trans, 

                                                  fZeroRotation));

  cover_in_module_Z_trans = - fiber_Z/2. - cover_Z/2.;

  cover_in_module_X_trans = ersatz_pmt_in_module_X_trans;

  cover_in_module_Y_trans = ersatz_pmt_in_module_Y_trans;

  module->AddNode(cover, 0, new TGeoCombiTrans( cover_in_module_X_trans, 

                                                   cover_in_module_Y_trans, 

                                                   cover_in_module_Z_trans, 

                                                  fZeroRotation));

   module->AddNode(module_housing, 0, new TGeoCombiTrans(.0,.0,.0, fZeroRotation));

  //------------------ top structure -----------------------------

  NeuRad->AddNode(module, 0, new TGeoCombiTrans(.0,.0,.0, fZeroRotation));

  top->AddNode(NeuRad, 0, new TGeoCombiTrans(trans_X,trans_Y,trans_Z, fZeroRotation));

  // ---------------   Finish   -----------------------------------------------

  gGeoMan->CloseGeometry();

  gGeoMan->CheckOverlaps(0.001);

  gGeoMan->PrintOverlaps();

  gGeoMan->Test();

  TFile* geoFile = new TFile(geoFileName, "RECREATE");

  top->Write();

  geoFile->Close();

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

}