From 66a072005a48afd81384d8e9a3fc371e1ba5a97e Mon Sep 17 00:00:00 2001
From: agheata <Andrei.Gheata@cern.ch>
Date: Fri, 4 Jul 2014 15:24:09 +0200
Subject: [PATCH] Fixes for the new parallel world navigation feature with a
 tutorial macro

(cherry picked from commit 2c7ca5526c5bbf40951dd40da6dc88af118dfbd6)
---
 geom/geom/inc/TGeoParallelWorld.h   |  6 +-
 geom/geom/inc/TGeoPhysicalNode.h    |  2 +
 geom/geom/inc/TGeoVolume.h          |  9 ++-
 geom/geom/src/TGeoManager.cxx       |  1 +
 geom/geom/src/TGeoParallelWorld.cxx | 45 +++++++++----
 geom/geom/src/TGeoPhysicalNode.cxx  | 15 +++++
 geom/geom/src/TGeoVolume.cxx        |  2 +
 tutorials/geom/parallel_world.C     | 97 +++++++++++++++++++++++++++++
 8 files changed, 161 insertions(+), 16 deletions(-)
 create mode 100644 tutorials/geom/parallel_world.C

diff --git a/geom/geom/inc/TGeoParallelWorld.h b/geom/geom/inc/TGeoParallelWorld.h
index 712c2450400..2ea9e68609c 100644
--- a/geom/geom/inc/TGeoParallelWorld.h
+++ b/geom/geom/inc/TGeoParallelWorld.h
@@ -35,13 +35,14 @@ protected :
    TObjArray         *fPhysical;       // array of physical nodes
    TGeoVolume        *fVolume;         // helper volume
    Bool_t             fIsClosed;       // Closed flag
+   Bool_t             fUseOverlaps;    // Activated if user defined overlapping candidates
 
    TGeoParallelWorld(const TGeoParallelWorld&); 
    TGeoParallelWorld& operator=(const TGeoParallelWorld&);
 
 public:
    // constructors
-   TGeoParallelWorld() : TNamed(),fGeoManager(0),fPhysical(0),fVolume(0),fIsClosed(kFALSE) {}
+   TGeoParallelWorld() : TNamed(),fGeoManager(0),fPhysical(0),fVolume(0),fIsClosed(kFALSE),fUseOverlaps(kFALSE) {}
    TGeoParallelWorld(const char *name, TGeoManager *mgr);
 
    // destructor
@@ -49,6 +50,9 @@ public:
    // API for adding components nodes
 //   void              AddNode(TGeoVolume *vol, TGeoMatrix *matrix=0);
    void              AddNode(TGeoPhysicalNode *pnode);
+   // Adding overlap candidates can improve performance
+   void              AddOverlap(TGeoVolume *vol);
+   
    
    // Closing a parallel geometry is mandatory
    Bool_t            CloseGeometry();
diff --git a/geom/geom/inc/TGeoPhysicalNode.h b/geom/geom/inc/TGeoPhysicalNode.h
index 5a3bb2794e5..06f93abedb0 100644
--- a/geom/geom/inc/TGeoPhysicalNode.h
+++ b/geom/geom/inc/TGeoPhysicalNode.h
@@ -30,6 +30,7 @@ class TGeoMatrix;
 class TGeoVolume;
 class TGeoNode;
 class TGeoShape;
+class TGeoNavigator;
 
 //////////////////////////////////////////////////////////////////////////////
 //                                                                          //
@@ -81,6 +82,7 @@ public:
    
  
    Bool_t            IsAligned() const {return TObject::TestBit(kGeoPNodeAligned);}
+   Bool_t            IsMatchingState(TGeoNavigator *nav) const;
    Bool_t            IsVolAttributes() const {return TObject::TestBit(kGeoPNodeVolAtt);}
    Bool_t            IsVisible() const {return TObject::TestBit(kGeoPNodeVisible);}
    Bool_t            IsVisibleFull() const {return TObject::TestBit(kGeoPNodeFull);}
diff --git a/geom/geom/inc/TGeoVolume.h b/geom/geom/inc/TGeoVolume.h
index b6c3ac027c6..6daa8fe1018 100644
--- a/geom/geom/inc/TGeoVolume.h
+++ b/geom/geom/inc/TGeoVolume.h
@@ -96,10 +96,11 @@ public:
       kVolumeOverlap =     BIT(17),
       kVolumeImportNodes = BIT(18),
       kVolumeMulti   =     BIT(19),
-      kVoxelsXYZ     =     BIT(20),
-      kVoxelsCyl     =     BIT(21),
+      kVoxelsXYZ     =     BIT(20), // not used
+      kVoxelsCyl     =     BIT(21), // not used
       kVolumeClone   =     BIT(22),
-      kVolumeAdded   =     BIT(23)
+      kVolumeAdded   =     BIT(23),
+      kVolumeOC      =     BIT(21)  // overlapping candidates
    };
    // constructors
    TGeoVolume();
@@ -158,6 +159,7 @@ public:
    Bool_t          IsActive() const {return TGeoAtt::IsActive();}
    Bool_t          IsActiveDaughters() const {return TGeoAtt::IsActiveDaughters();}
    Bool_t          IsAdded()     const {return TObject::TestBit(kVolumeAdded);}
+   Bool_t          IsOverlappingCandidate() const {return TObject::TestBit(kVolumeOC);}
    Bool_t          IsReplicated() const {return TObject::TestBit(kVolumeReplicated);}
    Bool_t          IsSelected() const  {return TObject::TestBit(kVolumeSelected);}
    Bool_t          IsCylVoxels() const {return TObject::TestBit(kVoxelsCyl);}
@@ -226,6 +228,7 @@ public:
    void            SetCurrentPoint(Double_t x, Double_t y, Double_t z);
    void            SetCylVoxels(Bool_t flag=kTRUE) {TObject::SetBit(kVoxelsCyl, flag); TObject::SetBit(kVoxelsXYZ, !flag);}
    void            SetNodes(TObjArray *nodes) {fNodes = nodes; TObject::SetBit(kVolumeImportNodes);}
+   void            SetOverlappingCandidate(Bool_t flag) {TObject::SetBit(kVolumeOC,flag);}
    void            SetShape(const TGeoShape *shape);
    void            SetTransparency(Char_t transparency=0) {if (fMedium) fMedium->GetMaterial()->SetTransparency(transparency);} // *MENU*
    void            SetField(TObject *field)          {fField = field;}
diff --git a/geom/geom/src/TGeoManager.cxx b/geom/geom/src/TGeoManager.cxx
index 3e0f7772e05..1bdda0e8486 100644
--- a/geom/geom/src/TGeoManager.cxx
+++ b/geom/geom/src/TGeoManager.cxx
@@ -3123,6 +3123,7 @@ void TGeoManager::RefreshPhysicalNodes(Bool_t lock)
    TIter next(gGeoManager->GetListOfPhysicalNodes());
    TGeoPhysicalNode *pn;
    while ((pn=(TGeoPhysicalNode*)next())) pn->Refresh();
+   if (fParallelWorld && fParallelWorld->IsClosed()) fParallelWorld->RefreshPhysicalNodes();
    if (lock) LockGeometry();
 }
 
diff --git a/geom/geom/src/TGeoParallelWorld.cxx b/geom/geom/src/TGeoParallelWorld.cxx
index b3315a4d3fc..9dec8878f0e 100644
--- a/geom/geom/src/TGeoParallelWorld.cxx
+++ b/geom/geom/src/TGeoParallelWorld.cxx
@@ -35,8 +35,9 @@ TGeoParallelWorld::TGeoParallelWorld(const char *name, TGeoManager *mgr)
                   : TNamed(name,""),
                     fGeoManager(mgr),
                     fPhysical(0),
-                    fVolume(new TGeoVolume()),
-                    fIsClosed(kFALSE)
+                    fVolume(new TGeoVolumeAssembly(name)),
+                    fIsClosed(kFALSE),
+                    fUseOverlaps(kFALSE)
 {
 // Default constructor
 }
@@ -46,7 +47,6 @@ TGeoParallelWorld::~TGeoParallelWorld()
 {
 // Destructor
    delete fPhysical;
-   delete fVolume;
 }
 
 //_____________________________________________________________________________
@@ -54,15 +54,20 @@ void TGeoParallelWorld::AddNode(TGeoPhysicalNode *pnode)
 {
 // Add a node normally to this world. Overlapping nodes not allowed
    if (fIsClosed) Fatal("AddNode", "Cannot add nodes to a closed parallel geometry");
-   if (pnode->GetVolume()->IsAssembly()) {
-      Error("AddNode", "Assemblies not allowed in parallel world, Physical node %s not added",
-             pnode->GetName());
-      return;
-   }   
    if (!fPhysical) fPhysical = new TObjArray(256);
    fPhysical->Add(pnode);
 }
 
+//_____________________________________________________________________________
+void TGeoParallelWorld::AddOverlap(TGeoVolume *vol)
+{
+// To use this optimization, the user should declare the full list of volumes
+// which may overlap with any of the physical nodes of the parallel world. Better
+// be done before misalignment
+   fUseOverlaps = kTRUE;
+   vol->SetOverlappingCandidate(kTRUE);
+}
+
 //_____________________________________________________________________________
 Bool_t TGeoParallelWorld::CloseGeometry()
 {
@@ -92,9 +97,10 @@ void TGeoParallelWorld::RefreshPhysicalNodes()
    TIter next(fPhysical);
    Int_t copy = 0;
    while ((pnode = (TGeoPhysicalNode*)next())) {
-      fVolume->AddNode(pnode->GetVolume(), copy++, pnode->GetMatrix());
+      fVolume->AddNode(pnode->GetVolume(), copy++, new TGeoHMatrix(*pnode->GetMatrix()));
    }
    // Voxelize the volume
+   fVolume->GetShape()->ComputeBBox();
    fVolume->Voxelize("ALL");
 }   
 
@@ -103,11 +109,14 @@ TGeoPhysicalNode *TGeoParallelWorld::FindNode(Double_t point[3])
 {
 // Finds physical node containing the point
    if (!fIsClosed) Fatal("FindNode", "Parallel geometry must be closed first");
+   TGeoNavigator *nav = fGeoManager->GetCurrentNavigator();
+   // Fast return if not in an overlapping candidate
+   if (fUseOverlaps && !nav->GetCurrentVolume()->IsOverlappingCandidate()) return 0;
    TGeoVoxelFinder *voxels = fVolume->GetVoxels();
    Int_t id;
    Int_t ncheck = 0;
    // get the list of nodes passing thorough the current voxel
-   TGeoNodeCache *cache = fGeoManager->GetCurrentNavigator()->GetCache();
+   TGeoNodeCache *cache = nav->GetCache();
    TGeoStateInfo &info = *cache->GetInfo();
    Int_t *check_list = voxels->GetCheckList(point, ncheck, info);
    cache->ReleaseInfo(); // no hierarchical use
@@ -135,6 +144,16 @@ TGeoPhysicalNode *TGeoParallelWorld::FindNextBoundary(Double_t point[3], Double_
 // Same functionality as TGeoNavigator::FindNextDaughterBoundary for the
 // parallel world
    if (!fIsClosed) Fatal("FindNode", "Parallel geometry must be closed first");
+   TGeoPhysicalNode *pnode = 0;
+   TGeoNavigator *nav = fGeoManager->GetCurrentNavigator();
+   // Fast return if not in an overlapping candidate
+   if (fUseOverlaps && !nav->GetCurrentVolume()->IsOverlappingCandidate()) return 0;
+   TIter next(fPhysical);
+   // Ignore the request if the current state in the main geometry matches one
+   // of the physical nodes in the parallel geometry
+   while ((pnode = (TGeoPhysicalNode*)next())) {
+      if (pnode->IsMatchingState(nav)) return 0;
+   }   
    Double_t snext = TGeoShape::Big();
    step = stepmax;
    TGeoVoxelFinder *voxels = fVolume->GetVoxels();
@@ -142,7 +161,6 @@ TGeoPhysicalNode *TGeoParallelWorld::FindNextBoundary(Double_t point[3], Double_
    Int_t nd = fVolume->GetNdaughters();
    Int_t i;
    TGeoNode *current;
-   TGeoPhysicalNode *pnode = 0;
    Double_t lpoint[3], ldir[3];
    const Double_t tolerance = TGeoShape::Tolerance();
    if (nd<5) {
@@ -170,7 +188,7 @@ TGeoPhysicalNode *TGeoParallelWorld::FindNextBoundary(Double_t point[3], Double_
    Int_t ncheck = 0;
    Int_t sumchecked = 0;
    Int_t *vlist = 0;
-   TGeoNodeCache *cache = fGeoManager->GetCurrentNavigator()->GetCache();
+   TGeoNodeCache *cache = nav->GetCache();
    TGeoStateInfo &info = *cache->GetInfo();
    cache->ReleaseInfo(); // no hierarchical use
    voxels->SortCrossedVoxels(point, dir, info);
@@ -198,6 +216,9 @@ TGeoPhysicalNode *TGeoParallelWorld::FindNextBoundary(Double_t point[3], Double_
 Double_t TGeoParallelWorld::Safety(Double_t point[3], Double_t safmax)
 {
 // Compute safety for the parallel world
+   TGeoNavigator *nav = fGeoManager->GetCurrentNavigator();
+   // Fast return if not in an overlapping candidate
+   if (fUseOverlaps && !nav->GetCurrentVolume()->IsOverlappingCandidate()) return TGeoShape::Big();
    Double_t local[3];
    Double_t safe = safmax;
    Double_t safnext;
diff --git a/geom/geom/src/TGeoPhysicalNode.cxx b/geom/geom/src/TGeoPhysicalNode.cxx
index 5e8dae01dbb..1f61c51a536 100644
--- a/geom/geom/src/TGeoPhysicalNode.cxx
+++ b/geom/geom/src/TGeoPhysicalNode.cxx
@@ -501,6 +501,21 @@ Bool_t TGeoPhysicalNode::SetPath(const char *path)
    return kTRUE;
 }
 
+//_____________________________________________________________________________
+Bool_t TGeoPhysicalNode::IsMatchingState(TGeoNavigator *nav) const
+{
+// Checks if a given navigator state matches this physical node
+   TGeoNodeCache *cache = nav->GetCache();
+   if (!cache) {
+      Fatal("SetBranchAsState","no state available");
+      return kFALSE;
+   }
+   TGeoNode **branch = (TGeoNode **) cache->GetBranch();
+   for (Int_t i=0; i<=fLevel; i++)
+      if (fNodes->At(i) != branch[i]) return kFALSE;
+   return kTRUE;   
+}   
+
 ClassImp(TGeoPNEntry)
 
 //_____________________________________________________________________________
diff --git a/geom/geom/src/TGeoVolume.cxx b/geom/geom/src/TGeoVolume.cxx
index 3030b9b6f7e..0d5fab64159 100644
--- a/geom/geom/src/TGeoVolume.cxx
+++ b/geom/geom/src/TGeoVolume.cxx
@@ -1698,6 +1698,7 @@ TGeoVolume *TGeoVolume::CloneVolume() const
    // copy extensions
    vol->SetUserExtension(fUserExtension);
    vol->SetFWExtension(fFWExtension);
+   vol->SetOverlappingCandidate(IsOverlappingCandidate());
    return vol;
 }
 
@@ -1767,6 +1768,7 @@ TGeoVolume *TGeoVolume::MakeCopyVolume(TGeoShape *newshape)
    CloneNodesAndConnect(vol);
 //   ((TObject*)vol)->SetBit(kVolumeImportNodes);
    ((TObject*)vol)->SetBit(kVolumeClone);
+   vol->SetOverlappingCandidate(IsOverlappingCandidate());
    return vol;       
 }    
 
diff --git a/tutorials/geom/parallel_world.C b/tutorials/geom/parallel_world.C
new file mode 100644
index 00000000000..62720411ee8
--- /dev/null
+++ b/tutorials/geom/parallel_world.C
@@ -0,0 +1,97 @@
+//______________________________________________________________________________
+void parallel_world(Bool_t usepw=kTRUE, Bool_t useovlp=kTRUE)
+{
+// Misaligning geometry generate in many cases overlaps, due to the idealization
+// of the design and the fact that in real life movements of the geometry volumes
+// have constraints and are correlated. This typically generates inconsistent
+// response of the navigation methods, leading to inefficiencies during tracking,
+// errors in the material budget calculations, and so on. Among those, there are
+// dangerous cases when the hidden volumes are sensitive.
+// This macro demonstrates how to use the "parallel world" feature to assign
+// highest navigation priority to some physical paths in geometry.
+// 
+
+   TGeoManager *geom = new TGeoManager("parallel_world", "Showcase for prioritized physical paths");
+   TGeoMaterial *matV = new TGeoMaterial("Vac", 0,0,0);
+   TGeoMedium *medV = new TGeoMedium("MEDVAC",1,matV);
+   TGeoMaterial *matAl = new TGeoMaterial("Al", 26.98,13,2.7);
+   TGeoMedium *medAl = new TGeoMedium("MEDAL",2,matAl);
+   TGeoMaterial *matSi = new TGeoMaterial("Si", 28.085,14,2.329);
+   TGeoMedium *medSi = new TGeoMedium("MEDSI",3,matSi);
+   TGeoVolume *top = gGeoManager->MakeBox("TOP",medV,100,400,1000);
+   gGeoManager->SetTopVolume(top);
+
+   // Shape for the support block
+   TGeoBBox *sblock = new TGeoBBox("sblock", 20,10,2);
+   // The volume for the support
+   TGeoVolume *support = new TGeoVolume("block",sblock, medAl);
+   support->SetLineColor(kGreen);
+   
+   // Shape for the sensor to be prioritized in case of overlap
+   TGeoBBox *ssensor = new TGeoBBox("sensor", 19,9,0.2);
+   // The volume for the sensor
+   TGeoVolume *sensor = new TGeoVolume("sensor",ssensor, medSi);
+   sensor->SetLineColor(kRed);
+   
+   // Chip assembly of support+sensor
+   TGeoVolumeAssembly *chip = new TGeoVolumeAssembly("chip");
+   chip->AddNode(support, 1);
+   chip->AddNode(sensor,1, new TGeoTranslation(0,0,-2.1));
+   
+   // A ladder that normally sags 
+   TGeoBBox *sladder = new TGeoBBox("sladder", 20,300,5);
+   // The volume for the ladder
+   TGeoVolume *ladder = new TGeoVolume("ladder",sladder, medAl);   
+   ladder->SetLineColor(kBlue);
+   
+   // Add nodes
+   top->AddNode(ladder,1);
+   for (Int_t i=0; i<10; i++)
+      top->AddNode(chip, i+1, new TGeoTranslation(0, -225.+50.*i, 10));
+
+   gGeoManager->CloseGeometry();
+   TGeoParallelWorld *pw = 0;
+   if (usepw) pw = gGeoManager->CreateParallelWorld("priority_sensors");
+// Align chips
+   align();
+   if (usepw) {
+      if (useovlp) pw->AddOverlap(ladder);
+      pw->CloseGeometry();
+      gGeoManager->SetUseParallelWorldNav(kTRUE);
+   } 
+   TString cname;
+   cname = usepw ? "cpw" : "cnopw";
+   TCanvas *c = (TCanvas*)gROOT->GetListOfCanvases()->FindObject(cname);
+   if (c) c->cd();
+   else   c = new TCanvas(cname, "",800,600);
+   top->Draw();
+//   top->RandomRays(0,0,0,0,sensor->GetName());
+   // Track random "particles" coming from the block side and draw only the tracklets
+   // actually crossing one of the sensors. Note that some of the tracks coming
+   // from the outer side may see the full sensor, while the others only part of it.
+   TStopwatch timer;
+   timer.Start();
+   top->RandomRays(100000,0,0,-30,sensor->GetName());
+   timer.Stop();
+   timer.Print();
+   TView3D *view = (TView3D*)gPad->GetView();
+   view->SetParallel();
+   view->Side();
+}
+
+//______________________________________________________________________________
+void align()
+{
+// Aligning 2 sensors so they will overlap with the support. One sensor is positioned
+// normally while the other using the shared matrix
+   TGeoPhysicalNode *node;
+   TGeoParallelWorld *pw = gGeoManager->GetParallelWorld();
+   Double_t sag;
+   for (Int_t i=0; i<10; i++) {
+      node = gGeoManager->MakePhysicalNode(TString::Format("/TOP_1/chip_%d",i+1));
+      sag = 8.-0.494*(i-4.5)*(i-4.5);
+      TGeoTranslation *tr = new TGeoTranslation(0., -225.+50.*i, 10-sag);
+      node->Align(tr);
+      if (pw) pw->AddNode(node);
+   }   
+}
-- 
GitLab