TGeoVolume* KVSpectroDetector::GetGeoVolume(const Char_t* name, const Char_t* material, const Char_t* shape_name, const Char_t* params){
	// Construct a TGeoVolume shape which can be used to represent
	// a detector in the current geometry managed by gGeoManager.
	// If the argument material is empty, the name of the detector is used.
	// Input: name - name given to the volume.
	//        material - material of the volume. The list of available
	//                   materials can be found with 
	//                   det->GetRangeTable()->GetListOfMaterials()->ls()
	//                   where det is a KVSpectroDetector or another
	//                   object inheriting from KVMaterial.
	//        shape_name - name of the shape associated to the volum
	//                     (Box, Arb8, Cone, Sphere, ...),  given
	//                     by the short name of the shape used in
	//                     the methods XXX:
	//                     TGeoManger::MakeXXX(...)


	TGeoMedium *med = GetGeoMedium(material);
	if(!med) return NULL;
	TString method = Form("Make%s",shape_name);
	TString  parameters = Form("%p,%p,%s",name,med,params);

	Info("GetGeoVolume","Trying to run the command gGeoManager->%s(%s)",method.Data(),parameters.Data());

	 gGeoManager->Execute(method.Data(),parameters.Data());
	 TGeoVolume* vol = (TGeoVolume*)gGeoManager->GetListOfVolumes()->Last(); 
	if(vol) vol->SetLineColor(med->GetMaterial()->GetDefaultColor());
	return vol;
}

void KVSpectroDetector::UpdateVolumeAndNodeNames(){
	// Update the names of Volumes and the names of the Nodes of this
	// detector.
	// The name of the volume representing the detector (returned
	// by GetAbsGeoVolume()) is DET_<detector name>. 
	// The name of the active volumes is ACTIVE_<detector name>_<material name>.
	// The name of the other volumes is <detector name>_<material name>.

	GetAbsGeoVolume()->SetName( Form("DET_%s", GetName() ) );
	TObjArray  *nodes = GetAbsGeoVolume()->GetNodes();
	TGeoNode   *node  = NULL;
	TGeoVolume *vol   = NULL;

	TIter next( nodes );
	while( (node = (TGeoNode *)next()) ){
		TString name, nname;
		vol = node->GetVolume();
		name.Form("%s_%s", GetName(), vol->GetMaterial()->GetName());
		if( GetActiveVolumes()->Contains( vol ) )
			name.Prepend("ACTIVE_");
		vol->SetName( name.Data() );
		nname = name;
		Int_t i=0;
		while( nodes->FindObject( nname.Data() ) )
			nname.Form("%s_%d",name.Data(), ++i);
		node->SetName( nname.Data() );
		node->SetNumber( i );
	}
}

UChar_t KVSpectroDetector::GetPosition( Double_t *XYZf, Int_t idx ){
	// Returns in the array 'XYZf', the coordinates (in cm) of a point randomly drawn in the 
	// active volume with index 'idx'. We assume that the shape of this
	// volume is a TGeoBBox. These coordinates are given in the focal-plane
	// of reference.
	// The function returns a binary code where:
	//   - bit 0 = 1 if X is OK  (001);
	//   - bit 1 = 1 if Y is OK  (010);
	//   - bit 2 = 1 if Z is OK  (100);
	//
	//  If no coordinates are OK, the returned value is null (000) and if 
	// X, Y and Z are OK then the returned value is equal 7 (111).

	TGeoVolume *vol = GetActiveVolume(idx);
	if( !vol ||  !PositionIsOK() ) return 0;

   	const TGeoBBox *box = (TGeoBBox *)vol->GetShape();
   	Double_t dx = box->GetDX();
   	Double_t dy = box->GetDY();
   	Double_t dz = box->GetDZ();
   	Double_t ox = (box->GetOrigin())[0];
   	Double_t oy = (box->GetOrigin())[1];
   	Double_t oz = (box->GetOrigin())[2];
   	Double_t xyz[3];
  	xyz[0] = ox+( TestBit(kRdmPos) ? dx*(2*gRandom->Rndm()-1) : 0.);
   	xyz[1] = oy+( TestBit(kRdmPos) ? dy*(2*gRandom->Rndm()-1) : 0.);
   	xyz[2] = oz+( TestBit(kRdmPos) ? dz*(2*gRandom->Rndm()-1) : 0.);
	if( ActiveVolumeToFocal( xyz , XYZf, idx ) ) return 7;
	return 0;
}

int main() {

  TGeoNode *node =  NULL;
  TGeoVolume *vol = NULL;
  
  LMCgeomN *g = new LMCgeomN("Telescope");
  TObjArray *oa = g->GetGeoManager()->GetListOfNodes();
  for (int i=0; i<oa->GetEntries(); i++) {
    node =  (TGeoNode*)oa->At(i);
    vol = node->GetVolume();
    cout << "= " << node->GetName() << "  " << vol->GetName()  <<endl;
    TObjArray *vnodes = vol->GetNodes();
    cout << vnodes->GetEntries() << endl;
    for (int j=0; j<vnodes->GetEntries(); j++) {
      node = (TGeoNode*)vnodes->At(j);
      cout << "== " << node->GetName() << endl;
      vol = node->GetVolume();
      TObjArray *vnodes1 = vol->GetNodes();
      for (int k=0; k<vnodes1->GetEntries(); k++) {
	node = (TGeoNode*)vnodes1->At(k);
	cout << "=== " << node->GetName() << endl;
	vol = node->GetVolume();
	TObjArray *vnodes2 = vol->GetNodes();
	if(!vnodes2) continue;
	for (int q=0; q<vnodes2->GetEntries(); q++) {
	   node = (TGeoNode*)vnodes2->At(q);
	   cout << "==== " << node->GetName() << endl;
	 }
	
      }
    }
  }

  return 0;
}

void fw_simGeo_set_volume_color_by_material(const char* material_re, Bool_t use_names, Color_t color, Char_t transparency=-1)
{
   // Note: material_re is a perl regexp!
   // If you want exact match, enclose in begin / end meta characters (^ / $):
   //   set_volume_color_by_material("^materials:Silicon$", kRed);

   TPMERegexp re(material_re, "o");
   TGeoMaterial *m;
   TIter it(FWGeometryTableViewManager_GetGeoManager()->GetListOfMaterials());
   while ((m = (TGeoMaterial*) it()) != 0)
   {
      if (re.Match(use_names ? m->GetName() : m->GetTitle()))
      {
         if (transparency != -1)
         {
            m->SetTransparency(transparency);
         }
         TGeoVolume *v;
         TIter it2(FWGeometryTableViewManager_GetGeoManager()->GetListOfVolumes());
         while ((v = (TGeoVolume*) it2()) != 0)
         {
            if (v->GetMaterial() == m)
            {
               v->SetLineColor(color);
            }
         }
      }
   }
}

void Visualizer::Show(){

  std::cout<<"=============================================="<<std::endl;
  std::cout<<"========= Inside Expected SHOW() ============="<<std::endl;
  std::cout<<"=============================================="<<std::endl;

  TGeoVolume *top = gGeoManager->MakeBox("Top", NULL, kInfinity, kInfinity, kInfinity);
  gGeoManager->SetTopVolume(top);
  for(int i = 0 ; i < fVolumes.size() ; i++){
  top->AddNode(std::get<0>(fVolumes[i]), 1 , std::get<1>(fVolumes[i]));

  }

  top->SetLineColor(kGreen);
  gGeoManager->CloseGeometry();
  #ifndef USE_OGL
  top->Draw();
  #else
  top->Draw("ogl"); //to display the geometry using openGL
  #endif
  //
  //TPad::x3d("OPENGL");
  gGeoManager->Export("plane.root");
  //top->Export("planeTop.root");
  //fApp->Run();
}

TGeoVolume* create_tof_pole()
{
  // needed materials
  TGeoMedium* boxVolMed   = gGeoMan->GetMedium(BoxVolumeMedium);
  TGeoMedium* airVolMed   = gGeoMan->GetMedium(KeepingVolumeMedium);
   
  Float_t dx=Pole_Size_X;
  Float_t dy=Pole_Size_Y;
  Float_t dz=Pole_Size_Z;
  Float_t width_alux=Pole_Thick_X;
  Float_t width_aluy=Pole_Thick_Y;
  Float_t width_aluz=Pole_Thick_Z;
  
  TGeoVolume* pole = new TGeoVolumeAssembly("Pole");
  TGeoBBox*   pole_alu_box = new TGeoBBox("", dx/2., dy/2., dz/2.);
  TGeoVolume* pole_alu_vol = 
    new TGeoVolume("pole_alu", pole_alu_box, boxVolMed);
  pole_alu_vol->SetLineColor(kGreen); // set line color for the alu box
  pole_alu_vol->SetTransparency(20); // set transparency for the TOF
  TGeoTranslation* pole_alu_trans 
    = new TGeoTranslation("", 0., 0., 0.);
  pole->AddNode(pole_alu_vol, 0, pole_alu_trans);

  TGeoBBox* pole_air_box = new TGeoBBox("", dx/2.-width_alux, dy/2.-width_aluy, dz/2.-width_aluz);
  TGeoVolume* pole_air_vol = 
    new TGeoVolume("pole_air", pole_air_box, airVolMed);
  pole_air_vol->SetLineColor(kYellow); // set line color for the alu box
  pole_air_vol->SetTransparency(70); // set transparency for the TOF
  TGeoTranslation* pole_air_trans 
    = new TGeoTranslation("", 0., 0., 0.);
  pole_alu_vol->AddNode(pole_air_vol, 0, pole_air_trans);

  return pole;
}

TGeoVolume* create_tof_bar(Float_t dx, Float_t dy, Float_t dz)
{
  // needed materials
  TGeoMedium* boxVolMed = gGeoMan->GetMedium(BoxVolumeMedium);
  TGeoMedium* airVolMed = gGeoMan->GetMedium(KeepingVolumeMedium);
   
  Float_t width_alux=Pole_Thick_X;
  Float_t width_aluy=Pole_Thick_Y;
  Float_t width_aluz=Pole_Thick_Z;
  
  TGeoVolume* bar = new TGeoVolumeAssembly("Bar");
  TGeoBBox*   bar_alu_box = new TGeoBBox("", dx/2., dy/2., dz/2.);
  TGeoVolume* bar_alu_vol = 
    new TGeoVolume("bar_alu", bar_alu_box, boxVolMed);
  bar_alu_vol->SetLineColor(kGreen); // set line color for the alu box
  bar_alu_vol->SetTransparency(20); // set transparency for the TOF
  TGeoTranslation* bar_alu_trans 
    = new TGeoTranslation("", 0., 0., 0.);
  bar->AddNode(bar_alu_vol, 0, bar_alu_trans);

  TGeoBBox* bar_air_box = new TGeoBBox("", dx/2.-width_alux, dy/2.-width_aluy, dz/2.-width_aluz);
  TGeoVolume* bar_air_vol = 
    new TGeoVolume("bar_air", bar_air_box, airVolMed);
  bar_air_vol->SetLineColor(kYellow); // set line color for the alu box
  bar_air_vol->SetTransparency(70);   // set transparency for the TOF
  TGeoTranslation* bar_air_trans 
    = new TGeoTranslation("", 0., 0., 0.);
  bar_alu_vol->AddNode(bar_air_vol, 0, bar_air_trans);

  return bar;
}

TGeoVolume* create_tof_module(Int_t modType)
{
  Int_t cType = CounterTypeInModule[modType];
  Float_t dx=Module_Size_X[modType];
  Float_t dy=Module_Size_Y[modType];
  Float_t dz=Module_Size_Z[modType];
  Float_t width_aluxl=Module_Thick_Alu_X_left;
  Float_t width_aluxr=Module_Thick_Alu_X_right;
  Float_t width_aluy=Module_Thick_Alu_Y;
  Float_t width_aluz=Module_Thick_Alu_Z;

  Float_t shift_gas_box = (Module_Thick_Alu_X_right - Module_Thick_Alu_X_left)/2;

  Float_t dxpos=CounterXDistance[modType];
  Float_t startxpos=CounterXStartPosition[modType];
  Float_t dzoff=CounterZDistance[modType];
  Float_t rotangle=CounterRotationAngle[modType];

  TGeoMedium* boxVolMed          = gGeoMan->GetMedium(BoxVolumeMedium);
  TGeoMedium* noActiveGasVolMed  = gGeoMan->GetMedium(NoActivGasMedium);

  TString moduleName = Form("module_%d", modType);
  TGeoVolume* module = new TGeoVolumeAssembly(moduleName);

  TGeoBBox* alu_box = new TGeoBBox("", dx/2., dy/2., dz/2.);
  TGeoVolume* alu_box_vol = 
    new TGeoVolume("alu_box", alu_box, boxVolMed);
  alu_box_vol->SetLineColor(kGreen); // set line color for the alu box
  alu_box_vol->SetTransparency(20); // set transparency for the TOF
  TGeoTranslation* alu_box_trans 
    = new TGeoTranslation("", 0., 0., 0.);
  module->AddNode(alu_box_vol, 0, alu_box_trans);

  TGeoBBox* gas_box = new TGeoBBox("", (dx-(width_aluxl+width_aluxr))/2., (dy-2*width_aluy)/2., (dz-2*width_aluz)/2.);
  TGeoVolume* gas_box_vol = 
    new TGeoVolume("gas_box", gas_box, noActiveGasVolMed);
  gas_box_vol->SetLineColor(kYellow); // set line color for the gas box
  gas_box_vol->SetTransparency(70); // set transparency for the TOF
  TGeoTranslation* gas_box_trans 
    = new TGeoTranslation("", shift_gas_box, 0., 0.);
  alu_box_vol->AddNode(gas_box_vol, 0, gas_box_trans);
  
  for (Int_t j=0; j<5; j++){ //loop over counters (modules)
    Float_t zpos;
    if (0 == modType) {
      zpos = dzoff *=-1;
    } else {
      zpos = 0.;
    }
    TGeoTranslation* counter_trans 
      = new TGeoTranslation("", startxpos+ j*dxpos , 0.0 , zpos);

    TGeoRotation* counter_rot = new TGeoRotation();
    counter_rot->RotateY(rotangle);
    TGeoCombiTrans* counter_combi_trans = new TGeoCombiTrans(*counter_trans, *counter_rot);
    gas_box_vol->AddNode(gCounter[cType], j, counter_combi_trans);
  }

  return module;
}

void set_volume_color_by_material(const char* material_re, Color_t color, Char_t transparency=-1)
{
   // Note: material_re is a perl regexp!
   // If you want exact match, enclose in begin / end meta characters (^ / $):
   //   set_volume_color_by_material("^materials:Silicon$", kRed);

   TPMERegexp re(material_re, "o");
   TGeoMaterial *m;
   TIter it(gGeoManager->GetListOfMaterials());
   while ((m = (TGeoMaterial*) it()) != 0)
   {
      if (re.Match(m->GetName()))
      {
         if (transparency != -1)
         {
            m->SetTransparency(transparency);
         }
         TGeoVolume *v;
         TIter it2(gGeoManager->GetListOfVolumes());
         while ((v = (TGeoVolume*) it2()) != 0)
         {
            if (v->GetMaterial() == m)
            {
               v->SetLineColor(color);
            }
         }
      }
   }
   full_update();
}

void geom_cms()
{
   TEveManager::Create();

   TFile::SetCacheFileDir(".");
   gGeoManager = gEve->GetGeometry("http://root.cern.ch/files/cms.root");
   gGeoManager->DefaultColors();

   TGeoVolume* top = gGeoManager->GetTopVolume()->FindNode("CMSE_1")->GetVolume();

   TEveGeoTopNode* trk = new TEveGeoTopNode(gGeoManager, top->FindNode("TRAK_1"));
   trk->SetVisLevel(6);
   gEve->AddGlobalElement(trk);

   TEveGeoTopNode* calo = new TEveGeoTopNode(gGeoManager, top->FindNode("CALO_1"));
   calo->SetVisLevel(3);
   gEve->AddGlobalElement(calo);

   TEveGeoTopNode* muon = new TEveGeoTopNode(gGeoManager, top->FindNode("MUON_1"));
   muon->SetVisLevel(4);
   gEve->AddGlobalElement(muon);

   gEve->FullRedraw3D(kTRUE);

   // EClipType not exported to CINT (see TGLUtil.h):
   // 0 - no clip, 1 - clip plane, 2 - clip box
   TGLViewer *v = gEve->GetDefaultGLViewer();
   v->GetClipSet()->SetClipType(TGLClip::EType(1));
   v->ColorSet().Background().SetColor(kMagenta+4);
   v->SetGuideState(TGLUtil::kAxesEdge, kTRUE, kFALSE, 0);
   v->RefreshPadEditor(v);

   v->CurrentCamera().RotateRad(-1.2, 0.5);
   v->DoDraw();
}

TGeoHMatrix &KVSpectroDetector::GetActiveVolToFocalMatrix(Int_t i ) const{ 
	// Returns the matrix which transforms coordinates form the reference
	// frame of the active volume 'i' to the reference frame of the focal
	// plan.
	
	static TGeoHMatrix mm;
	mm.Clear();

	TGeoVolume *vol = (TGeoVolume *)fActiveVolumes->At(i);

	if( vol ){

		Int_t prev_id = gGeoManager->GetCurrentNodeId();
		// To be sure that the matrices will be calculated again
		gGeoManager->CdTop();

		gGeoManager->CdNode( vol->GetUniqueID() );

		// matrix 'volume to target'
		TGeoMatrix *vol_to_tgt = gGeoManager->GetCurrentMatrix();

		if( fFocalToTarget ) mm = fFocalToTarget->Inverse()*(*vol_to_tgt);
		else                 mm = *vol_to_tgt;

		mm.SetName( Form("%s_to_focal",vol->GetName()) );

		// just to not create problems if this method is called during a tracking
		gGeoManager->CdNode( prev_id );
	}
	else mm.SetName("Identity_matrix");

		return mm;
}

GeoHandler& GeoHandler::i_collect(const TGeoNode* current, int level, Region rg, LimitSet ls) {
  TGeoVolume* volume = current->GetVolume();
  TObjArray* nodes = volume->GetNodes();
  int num_children = nodes ? nodes->GetEntriesFast() : 0;
  Volume vol(volume);
  Region   region = vol.region();
  LimitSet limits = vol.limitSet();

  if ( m_propagateRegions )  {
    if ( !region.isValid() && rg.isValid() )   {
      region = rg;
      vol.setRegion(region);
    }
    if ( !limits.isValid() && ls.isValid() )  {
      limits = ls;
      vol.setLimitSet(limits);
    }
  }
  (*m_data)[level].insert(current);
  //printf("GeoHandler: collect level:%d %s\n",level,current->GetName());
  if (num_children > 0) {
    for (int i = 0; i < num_children; ++i) {
      TGeoNode* node = (TGeoNode*) nodes->At(i);
      i_collect(node, level + 1, region, limits);
    }
  }
  return *this;
}

void Visualizer::Show(TGeoVolume *vol){
  TGeoVolume *top = gGeoManager->MakeBox("Top", NULL, kInfinity, kInfinity, kInfinity);
  gGeoManager->SetTopVolume(top);
  //TGeoVolume *vol = fGeoManager->MakeSphere("SPHERE", NULL, 30, 40, 0, 180, 0, 360);
  top->AddNode(vol, 1);
  gGeoManager->CloseGeometry();
  top->Draw();
  //fApp->Run();
}

void s_intersection()
{
   gROOT->GetListOfCanvases()->Delete();
   TCanvas *c = new TCanvas("composite shape", "Intersection boolean operation", 700, 1000);

   c->Divide(1,2,0,0);
   c->cd(2);
   gPad->SetPad(0,0,1,0.4);
   c->cd(1);
   gPad->SetPad(0,0.4,1,1);
   
   if (gGeoManager) delete gGeoManager;
   
   new TGeoManager("xtru", "poza12");
   TGeoMaterial *mat = new TGeoMaterial("Al", 26.98,13,2.7);
   TGeoMedium *med = new TGeoMedium("MED",1,mat);
   TGeoVolume *top = gGeoManager->MakeBox("TOP",med,100,100,100);
   gGeoManager->SetTopVolume(top);

   // define shape components with names
   TGeoBBox *box = new TGeoBBox("bx", 40., 40., 40.); 
   TGeoSphere *sph = new TGeoSphere("sph", 40., 45.);
   // define named geometrical transformations with names
   TGeoTranslation *tr = new TGeoTranslation(0., 0., 45.);
   tr->SetName("tr");
   // register all used transformations
   tr->RegisterYourself();
   // create the composite shape based on a Boolean expression
   TGeoCompositeShape *cs = new TGeoCompositeShape("mir", "sph:tr * bx");

   TGeoVolume *vol = new TGeoVolume("COMP2",cs);
   top->AddNode(vol,1);
   gGeoManager->CloseGeometry();
   gGeoManager->SetNsegments(100);
   top->Draw();
   MakePicture();

   c->cd(2);

   TPaveText *pt = new TPaveText(0.01,0.01,0.99,0.99);

   pt->SetLineColor(1);

   TText *text = pt->AddText("TGeoCompositeShape - composite shape class");

   text->SetTextColor(2);
   pt->AddText("----- Here is an example of boolean intersection operation : A * B");
   pt->AddText("----- A == sphere (with inner radius non-zero), B == box");
   pt->AddText(" ");
   pt->SetAllWith("-----","color",4);
   pt->SetAllWith("-----","font",72);
   pt->SetAllWith("-----","size",0.04);
   pt->SetTextAlign(12);
   pt->SetTextSize(0.044);
   pt->Draw();
   c->cd(1);
}

void Create_TOF_Geometry_v12a() {
  // Load the necessary FairRoot libraries 
  gROOT->LoadMacro("$VMCWORKDIR/gconfig/basiclibs.C");
  basiclibs();
  gSystem->Load("libGeoBase");
  gSystem->Load("libParBase");
  gSystem->Load("libBase");

  // Load needed material definition from media.geo file
  create_materials_from_media_file();

  // Get the GeoManager for later usage
  gGeoMan = (TGeoManager*) gROOT->FindObject("FAIRGeom");
  gGeoMan->SetVisLevel(7);  

  // Create the top volume 
  /*
  TGeoBBox* topbox= new TGeoBBox("", 1000., 1000., 1000.);
  TGeoVolume* top = new TGeoVolume("top", topbox, gGeoMan->GetMedium("air"));
  gGeoMan->SetTopVolume(top);
  */

  TGeoVolume* top = new TGeoVolumeAssembly("TOP");
  gGeoMan->SetTopVolume(top);
 
  TGeoVolume* tof = new TGeoVolumeAssembly(geoVersion);
  top->AddNode(tof, 1);
  
  for(Int_t counterType = 0; counterType < NumberOfDifferentCounterTypes; counterType++) { 
    gCounter = create_new_counter(counterType);
  }

  for(Int_t moduleType = 0; moduleType < NofModuleTypes; moduleType++) { 
    gModules[moduleType] = create_new_tof_module(moduleType);
  }

  gPole = create_tof_pole();

  position_tof_poles();
  position_inner_tof_modules();
  position_outer_tof_modules();
  

  gGeoMan->CloseGeometry();
  gGeoMan->CheckOverlaps(0.001);
  gGeoMan->PrintOverlaps();
  gGeoMan->Test();

  TFile* outfile = new TFile(FileName,"RECREATE");
  top->Write();
  //gGeoMan->Write();
  outfile->Close();
  //  top->Draw("ogl");
  //top->Raytrace();

}

void visibility_all_volumes(Bool_t vis_state)
{
   TGeoVolume *v;
   TIter it(gGeoManager->GetListOfVolumes());
   while ((v = (TGeoVolume*) it()) != 0)
   {
      v->SetVisibility(vis_state);
   }
   full_update();
}

void KVFAZIALNS2016::BuildFAZIA()
{
   //Build geometry of FAZIASYM
   //All telescopes are : Si(300µm)-Si(500µm)-CsI(10cm)
   //No attempt has been made to implement real thicknesses
   //
   Info("BuildFAZIA", "Compact geometry, %f cm from target",
        fFDist);

   TGeoVolume* top = gGeoManager->GetTopVolume();

   Double_t distance_block_cible = fFDist * KVUnits::cm;
   Double_t thick_si1 = 300 * KVUnits::um;
   TGeoTranslation trans;
   trans.SetDz(distance_block_cible + thick_si1 / 2.);

   KVFAZIABlock* block = new KVFAZIABlock;

   TGeoRotation rot1, rot2;
   TGeoHMatrix h;
   TGeoHMatrix* ph = 0;
   Double_t theta = 0;
   Double_t phi = 0;

   Double_t theta_min = fFThetaMin;//smallest lab polar angle in degrees
   Double_t centre_hole = 2.*tan(theta_min * TMath::DegToRad()) * distance_block_cible;
   Double_t dx = (block->GetTotalSideWithBlindage()) / 2.;

   TVector3 centre;
   for (Int_t bb = 0; bb < fNblocks; bb += 1) {

      if (bb == 1)        centre.SetXYZ(-1 * (dx - centre_hole / 2), -dx - centre_hole / 2, distance_block_cible);
      else if (bb == 2)   centre.SetXYZ(-1 * (dx + centre_hole / 2), dx - centre_hole / 2, distance_block_cible);
      else if (bb == 3)   centre.SetXYZ(-1 * (-dx + centre_hole / 2), dx + centre_hole / 2, distance_block_cible);
      else if (bb == 0)   centre.SetXYZ(-1 * (-dx - centre_hole / 2), -dx + centre_hole / 2, distance_block_cible);
      else if (bb == 4)   centre.SetXYZ(-1 * (-dx - centre_hole / 2), -3 * dx + centre_hole / 2, distance_block_cible); //centre.SetXYZ(-1 * (dx - centre_hole / 2), -3 * dx - centre_hole / 2, distance_block_cible);
      else {
         Warning("BuildFAZIA", "Block position definition is done only for %d blocks", fNblocks);
      }
      theta = centre.Theta() * TMath::RadToDeg();
      phi = centre.Phi() * TMath::RadToDeg();
      printf("BLK #%d => theta=%1.2lf - phi=%1.2lf\n", bb, theta, phi);

      rot2.SetAngles(phi + 90., theta, 0.);
      rot1.SetAngles(-1.*phi, 0., 0.);
      h = rot2 * trans * rot1;
      ph = new TGeoHMatrix(h);
      top->AddNode(block, bb, ph);
   }

   // add telescope for elastic scattering monitoring
//   RutherfordTelescope();
   // Change default geometry import angular range for rutherford telescope
   SetGeometryImportParameters(.25, 1., 1.84);
}

void geom_cms_stereo()
{
   TEveManager::Create();

   TFile::SetCacheFileDir(".");
   gGeoManager = gEve->GetGeometry("http://root.cern.ch/files/cms.root");
   gGeoManager->DefaultColors();

   TGeoVolume* top = gGeoManager->GetTopVolume()->FindNode("CMSE_1")->GetVolume();

   TEveGeoTopNode* trk = new TEveGeoTopNode(gGeoManager, top->FindNode("TRAK_1"));
   trk->SetVisLevel(6);
   gEve->AddGlobalElement(trk);

   TEveGeoTopNode* calo = new TEveGeoTopNode(gGeoManager, top->FindNode("CALO_1"));
   calo->SetVisLevel(3);
   gEve->AddGlobalElement(calo);

   TEveGeoTopNode* muon = new TEveGeoTopNode(gGeoManager, top->FindNode("MUON_1"));
   muon->SetVisLevel(4);
   gEve->AddGlobalElement(muon);

   // --- Stereo ---

   TEveWindowSlot* slot = 0;
   slot = TEveWindow::CreateWindowInTab(gEve->GetBrowser()->GetTabRight());

   TEveViewer* sv = new TEveViewer("Stereo GL", "Stereoscopic view");
   sv->SpawnGLViewer(gEve->GetEditor(), kTRUE);
   sv->AddScene(gEve->GetGlobalScene());

   slot->ReplaceWindow(sv);

   gEve->GetViewers()->AddElement(sv);

   gEve->GetBrowser()->GetTabRight()->SetTab(1);

   // --- Redraw ---

   gEve->FullRedraw3D(kTRUE);
   gEve->EditElement(sv);

   // --- Fix editor ---

   // EClipType not exported to CINT (see TGLUtil.h):
   // 0 - no clip, 1 - clip plane, 2 - clip box
   TGLViewer *v = gEve->GetDefaultGLViewer();
   v->GetClipSet()->SetClipType(TGLClip::EType(1));
   v->ColorSet().Background().SetColor(kMagenta+4);
   v->SetGuideState(TGLUtil::kAxesEdge, kTRUE, kFALSE, 0);
   v->RefreshPadEditor(v);

   v->CurrentCamera().RotateRad(-1.2, 0.5);
   v->DoDraw();
}

void s_difference()
{
   gROOT->GetListOfCanvases()->Delete();
   TCanvas *c = new TCanvas("composite shape", "Difference boolean operation", 700, 1000);

   c->Divide(1,2,0,0);
   c->cd(2);
   gPad->SetPad(0,0,1,0.4);
   c->cd(1);
   gPad->SetPad(0,0.4,1,1);
   
   if (gGeoManager) delete gGeoManager;
   
   new TGeoManager("xtru", "poza12");
   TGeoMaterial *mat = new TGeoMaterial("Al", 26.98,13,2.7);
   TGeoMedium *med = new TGeoMedium("MED",1,mat);
   TGeoVolume *top = gGeoManager->MakeBox("TOP",med,100,100,100);
   gGeoManager->SetTopVolume(top);

   // define shape components with names
   TGeoTorus *tor = new TGeoTorus("tor", 45., 15., 20., 45., 145.); 
   TGeoSphere *sph = new TGeoSphere("sph", 20., 45., 0., 180., 0., 270.);
   // create the composite shape based on a Boolean expression
   TGeoCompositeShape *cs = new TGeoCompositeShape("mir", "sph - tor");

   TGeoVolume *vol = new TGeoVolume("COMP3",cs);
   top->AddNode(vol,1);
   gGeoManager->CloseGeometry();
   gGeoManager->SetNsegments(60);
   top->Draw();
   MakePicture();

   c->cd(2);

   TPaveText *pt = new TPaveText(.01, .01, .99, .99);

   pt->SetLineColor(1);

   TText *text = pt->AddText("TGeoCompositeShape - composite shape class");

   text->SetTextColor(2);

   pt->AddText("----- It's an example of boolean difference: A - B");
   pt->AddText("----- A == part of sphere (0-180, 0-270), B == partial torus (45-145)");
   pt->AddText(" ");
   pt->SetAllWith("-----","color",4);
   pt->SetAllWith("-----","font",72);
   pt->SetAllWith("-----","size",0.04);
   pt->SetTextAlign(12);
   pt->SetTextSize(0.044);
   pt->Draw();
   c->cd(1);
}