Hi Kevin,
As I wrote you last time, computing distances based on TGeoShape
methods makes sense only if you do not have rotations/translations - in
fact a geometry. These methods (Contains(), DistToIn/Out()) are
implemented at shape level, but globally handled and optimized by the
TGeoManager class who takes care of geometrical transformations and
hierarchy.
If you will want to compute some distance in your example geometry for
instance, you will deal now with nodes and volumes rather than shapes. A
volume is a shape associated with a medium and possibly containing other
positioned volumes (nodes) inside. These objects are linked togeather
making a hierarchy : the same object in your geometry might be
positioned several times,so the only way to identify an unique object is
by full path : as in your example /TOP_1/aper_1 is different from
/TOP_1/aper_2 - they are the same object but positioned with different
translations.
Now here is how you can compute distances in your geometry : you have
to define a point and a direction in the master reference system, which
is the one attached to the TOP volume; each volume have a LOCAL frame
which is the one attached to its shape; positioned volumes inside have
local transformations defined with respect to this LOCAL frame; the
global transformation matrix associated to a given path (object)in your
geometry is obtained by multiplying all local transformations from TOP
to bottom in the corresponding branch.
Double_t point[3];
Double_t dir[3];
// define them here
...
// first we have to initialize the modeler with this point and
// direction and to find out where we are (which is the deepest path)
// containing the point
TGeoNode *current = gGeoManager->InitTrack(point,dir);
// this method not only initializes the state, but also performs a
// search (TGeoManager::FindNode())in the tree converting always the
// point to the local frame of a given volume (in fact shape) and
// calling TGeoShape::Contains(localpoint).
// print this out :
if (gGeoManager->IsOutside()) printf("Point is outside geometry\n");
else printf("Current path: %s\n", gGeoManager->GetPath());
// compute the distance to the next boundary
gGeoManager->FindNextBoundary();
Double_t dist = gGeoManager->GetStep();
// you might also want to compute the safe distance (maximum distance
// in any direction not crossing any boundary)
gGeoManager->Safety();
Double_t safety = gGeoManager->GetSafeDistance();
// or you might want to make the step=dist and see where you end up
TGeoNode *next = gGeoManager->Step();
if (gGeoManager->IsEntering())
printf("we have crossed a boundary at d=%f\n",dist);
if (gGeoManager->IsOutside()) printf("we have exited geometry\n");
printf("Current path: %s\n", gGeoManager->GetPath());
I hope this helps.
N.B. I cannot check the geometry in your example right now, but I have
noticed that you used AddNodeOverlap for positioning some volumes. Note
that you should use this ONLY if these objects really overlap each
other, otherwise nothing bad will happen, but you really decrease the
performance of the search algorithms, so I advice you to us AddNode()
instead.
Regards,
Andrei
Kevin Reil wrote:
> Hi all,
>
> If I create a TShape *shape= new TGeoBBox(dx,dy,dz,origin) I can find the
> distance to that shape using shape->DistToIn(point,dir); However, I cannot
> rotate the TShape at its location.
>
> Instead I use the TGeoManager and.... now I can position and rotate but
> cannot get shape->DistToIn(point,dir) since the position and rotation are
> not part of the TShape (they describe the node inside a volume).
>
> ...
>
> I create a TopVolume (top) and inside it place several volumes at
> different positions and rotations. I wish to know the distance from a
> point to the surface of
>
> if (pos not contained) {
> dist=top->GetNode(i)->GetVolume()->GetShape()->DistToIn(pos,dir);
> }
>
> But the TGeoMatrix used to position the node inside of TopVolume does not
> do anything to the shape.
>
> For nonrotated objects I can adjust the pos by
> pos[i]-=top->GetNode(i)->GetGetMatrix()->GetTranslation()[i];
>
> But when the object has been rotated the distance to the surface isn't
> correct.
>
> Is there a way to use DistToIn(pos,dir) in a TGeoVolume->AddNode() manner?
>
> An example script is attached.
>
> Cheers,
> Kevin
>
> |------------------------------------|---------------------------------|
> | Kevin Reil | 2575 Sand Hill Road, MS 26 |
> | X2447, 103D A&E Bldg. 041 | Menlo Park, CA 94025 |
> |------------------------------------|---------------------------------|
> | http://www.slac.stanford.edu/~reil | Office (650) 926-2447 |
> | reil@slac.stanford.edu | Home (650) 938-1767 |
> | http://reil.no-ip.org | Fax (650) 926-5368 |
> |----------------------------------------------------------------------|
> | And my father dwelt in a tent. |
> |----------------------------------------------------------------------|
>
> On Fri, 30 May 2003, Axel Naumann wrote:
>
>
>>Hi,
>>
>>Int_t GetNtuple(Int_t stage, TFile* &fin, TNtuple* &nin)
>>
>>You'll need to pass the arg by ref (better but uncommonway of stating
>>that: You'll need to define the method such that the arg gets passed by
>>ref), to allow GetNtuple to change the arg.
>>
>>Cheers, Axel.
>>
>>
>>>Hi all,
>>>
>>>I should know this but
>>>
>>>I have files
>>>photon_1.root with ntuple phot1
>>>photon_2.root with ntuple phot2
>>>photon_3.root with ntuple phot3
>>>
>>>If the file is there I want to return a pointer to it and a pointer to the
>>>ntuple it contains.
>>>
>>>If it is not there, I want to return a pointer to a file ready to write to
>>>and a pointer to an ntuple ready to fill.
>>>
>>>the new TFile and new TNtuple are lost when I go out of scope (I believe)
>>>so when I do nin->GetEntries() before leaving I get the anwer when I do it
>>>after leaving the function, I get a seg fault. I now its a scope thing but
>>>I cannot recall the solution.
>>>
>>>Core example it is below.
>>>
>>>Thanks,
>>>Kevin
>>>
>>>
>>>
>>>Int_t GetNtuple(Int_t stage, TFile *fin, TNtuple *nin)
>>>{
>>> char fname[256];
>>> sprintf(fname,"photons_%d.root",stage);
>>> char ntname[256];
>>> sprintf(ntname,"phot%d",stage);
>>> cout << ntname << " " << fname << endl;
>>> TFile *tfin = new TFile("photons_1.root","READ");
>>> if (!fin->IsOpen()) {
>>> fin = new TFile(fname,"RECREATE");
>>> nin=new TNtuple(ntname,"photons","x:y:z:dx:dy:dz:lambda");
>>> } else {
>>> tnin = (TNtuple*)gROOT->FindObject(ntname);
>>> }
>>> fin=tfin->Copy((TObject*)tfin);
>>> nin=tnin->Copy((TObject*)tnin);
>>> cout << nin->GetEntries() << endl;
>>> //works here (there are 10,000 if file is there else 0)
>>> return 0;
>>>}
>>>
>>>int main()
>>>{
>>>TNtuple *phot1;
>>>TFile *f1;
>>>
>>>Int_t nocode=GetNtuple(1,f1,phot1);
>>>cout << phot1->GetEntries() << endl; // segfault
>>>}
>>>
>>>|------------------------------------|---------------------------------|
>>>| Kevin Reil | 2575 Sand Hill Road, MS 26 |
>>>| X2447, 103D A&E Bldg. 041 | Menlo Park, CA 94025 |
>>>|------------------------------------|---------------------------------|
>>>| http://www.slac.stanford.edu/~reil | Office (650) 926-2447 |
>>>| reil@slac.stanford.edu | Home (650) 938-1767 |
>>>| http://reil.no-ip.org | Fax (650) 926-5368 |
>>>|----------------------------------------------------------------------|
>>>| And my father dwelt in a tent. |
>>>|----------------------------------------------------------------------|
>>>
>>>On Thu, 29 May 2003, Mayly Sanchez wrote:
>>>
>>>
>>>>Hi Rene,
>>>>the root file is here:
>>>>http://minos.phy.tufts.edu/msanchez/dchisq90_c.root
>>>>Mayly
>>>>
>>>>On Thursday, May 29, 2003, at 01:47 PM, Rene Brun wrote:
>>>>
>>>>
>>>>>Hi Mayly,
>>>>>
>>>>>Could you send the canvas.root file instead of canvas.gif?
>>>>>
>>>>>Rene Brun
>>>>>
>>>>>On Thu, 29 May
>>>>>2003, Mayly Sanchez wrote:
>>>>>
>>>>>
>>>>>>Thanks Rene, that fixed the contours but now it has a funny problem in
>>>>>>the lego plots. For some cases (I still have to determine the
>>>>>>pattern),
>>>>>>it draws lego plots where part of it shows as a wireframe while the
>>>>>>rest seems ok. The best way to explain is a picture:
>>>>>>http://minos.phy.tufts.edu/msanchez/dchisq90_lego.gif
>>>>>>
>>>>>>Any ideas?
>>>>>>Mayly
>>>>>>
>>>>>>On Wednesday, May 21, 2003, at 06:42 PM, Rene Brun wrote:
>>>>>>
>>>>>>
>>>>>>>Hi,
>>>>>>>
>>>>>>>Some changes have been made in the past few weeks by Olivier Couet
>>>>>>>to solve the kind of problems you are reporting.
>>>>>>>Take the version from CVS head and install from source.
>>>>>>>
>>>>>>>Rene Brun
>>>>>>>
>>>>>>>On Wed,
>>>>>>>21 May 2003, Mayly Sanchez wrote:
>>>>>>>
>>>>>>>
>>>>>>>>Hi,
>>>>>>>>I have a macro for drawing non-equidistant 2d contours that used to
>>>>>>>>work before but fails in most recent versions of root. The latest
>>>>>>>>test
>>>>>>>>has been done with 3.05/05.
>>>>>>>>
>>>>>>>>Here are the steps I follow:
>>>>>>>>{TH2F *ch = new
>>>>>>>>TH2F("ch","",nbinx,nbinxmin,nbinxmax,nbiny,nbinymin,nbinymax);
>>>>>>>>ch->Fill(x,y,weight);
>>>>>>>>ch->SetContour(2);
>>>>>>>>ch->SetContourLevel(0,0.0)
>>>>>>>>ch->SetContourLevel(1,2.0)
>>>>>>>>
>>>>>>>>ch->Draw("cont2");
>>>>>>>>}
>>>>>>>>
>>>>>>>>Now what happens is: that if I do a lego2 plot before drawing the
>>>>>>>>contours, the levels are set at the proper heights. Once I have run
>>>>>>>>Draw("cont2") or any of the other cont options it recalculates the
>>>>>>>>levels in a weird way. It seems to put the first level at the
>>>>>>>>GetMinimum value and the second half way between the maximum and the
>>>>>>>>minimum, as if it wanted to do equidistant contours.
>>>>>>>>
>>>>>>>>Did something change? Can I force it to do non-equidistant contours
>>>>>>>>again? I need this urgently so any workarounds are welcome, thanks,
>>>>>>>>
>>>>>>>>Mayly
>>>>>>>>
>>>>>>>>ps. the method with SetContourLevels(2,vector) was also tested and
>>>>>>>>gave
>>>>>>>>the same results
>>>>>>>>
>>>>>>>>
>>>>>>>
>>>>
>>
>>
>>------------------------------------------------------------------------
>>
>>//void rootgeom()
>>{
>> // DISTANCE UNITS ARE mm
>> //--- Definition of a simple geometry
>> gSystem->Load("libGeom");
>> TGeoManager *geom = new TGeoManager("simple1", "Simple geometry");
>>
>> //--- define some materials
>> TGeoMaterial *matVacuum = new TGeoMaterial("Vacuum", 0,0,0);
>>
>> //--- define some media
>> TGeoMedium *med;
>> TGeoMedium *Vacuum = new TGeoMedium("Vac",1, matVacuum);
>>
>> //--- make the top container volume
>> TGeoVolume *top = geom->MakeBox("TOP", Vacuum, 270., 270., 120.);
>> geom->SetTopVolume(top); // mandatory !
>>
>> // A 45 degree about y and 90 degree about y rotation
>> TGeoRotation *rot1 = new TGeoRotation("rot1", 0., -45, 0.);
>> TGeoRotation *rot2 = new TGeoRotation("rot2", 0., 90., 0.);
>>
>> // Location of apertures and mirrors
>>
>> // 3 apertures in a row in x
>> TGeoTranslation *al1 = new TGeoTranslation(100.0, 0.0, 50.0);
>> TGeoTranslation *al2 = new TGeoTranslation(200.0, 0.0, 50.0);
>> TGeoTranslation *al3 = new TGeoTranslation(300.0, 0.0, 50.0);
>> // 45 degree mirror to reflect into z
>> TGeoCombiTrans *ml1 = new TGeoCombiTrans(350.0, 0.0, 50.0,rot1);
>> // 3 apertures in z
>> TGeoCombiTrans *al4 = new TGeoCombiTrans(350.0, 0.0, 0.0,rot2);
>> TGeoCombiTrans *al5 = new TGeoCombiTrans(350.0, 0.0, -100.0,rot2);
>> TGeoCombiTrans *al6 = new TGeoCombiTrans(350.0, 0.0, -200.0,rot2);
>> // 45 degree mirror to reflect back to z
>> TGeoCombiTrans *ml2 = new TGeoCombiTrans(350.0, 0.0, -250.0,rot1);
>> // 2 more apertures in x
>> TGeoTranslation *al7 = new TGeoTranslation(400.0, 0.0, -250.0);
>> TGeoTranslation *al8 = new TGeoTranslation(500.0, 0.0, -250.0);
>>
>> // Full aperture dimensions. Will be /2 later to get span from center
>> Double_t aperdim[3]={ 0.001, 10.0, 100.0};
>> Double_t mir1dim[3] ={0.001, 100.0*sqrt(2), 100.0*sqrt(2)};
>>
>> aper = new TGeoBBox(aperdim[0]/2,aperdim[1]/2,aperdim[2]/2);
>> mir1 = new TGeoBBox(mir1dim[0]/2,mir1dim[1]/2,mir1dim[2]/2);
>>
>> TGeoVolume *v1=new TGeoVolume("aper",aper);
>> TGeoVolume *v2=new TGeoVolume("mir1",mir1);
>>
>> top->AddNodeOverlap(v1,1,al1);
>> top->AddNodeOverlap(v1,2,al2);
>> top->AddNodeOverlap(v1,3,al3);
>> top->AddNodeOverlap(v2,1,ml1);
>> top->AddNodeOverlap(v1,4,al4);
>> top->AddNodeOverlap(v1,5,al5);
>> top->AddNodeOverlap(v1,6,al6);
>> top->AddNodeOverlap(v2,2,ml2);
>> top->AddNodeOverlap(v1,7,al7);
>> top->AddNodeOverlap(v1,8,al8);
>>
>> Int_t nnodes=top->GetNdaughters();
>> for (Int_t i=0; i<nnodes; ++i) {
>> Double_t pos[3]={0,0,50};
>> Double_t dir[3]={1,0,0};
>> TGeoShape *cshape=top->GetNode(i)->GetVolume()->GetShape();
>> if(cshape->Contains(pos)) {
>> cout << "X=0 inside node " << i << endl;
>> } else {
>> cout << "This shape is ";
>> cout << cshape->DistToIn(pos,dir) << " mm away from x=0" << endl;
>> }
>> }
>>
>> //--- close the geometry
>> geom->CloseGeometry();
>>
>> //--- draw the ROOT box
>> geom->SetVisLevel(4);
>> top->Draw();
>> //if (gPad) gPad->x3d();
>>}
>>
>>
>>
>
This archive was generated by hypermail 2b29 : Thu Jan 01 2004 - 17:50:12 MET