Tracking Problem with TGeoPgon

Hello,

I’m using the geometry package to take a data set containing tracks to determine if and where a detector was hit. The detectors are small hexagonal cells that I’m representing with TGeoPgon objects. I have drawn the volume and ray traced ( with top->Raytrace(true) ) and I’m getting hexagonal cells arranged the way I want.

When I wrote a function to tell me if and where a cell was hit, it would return hits on cells that are not within the cell area. I didn’t know what I did wrong, so I shrunk my cells so they fit entirely within the region with the most tracks. What I saw was that it appeared that my hexagon was built out of squares with the 4 diagonal sides made up from 3 small squares. I tried increasing the size of my cells relative to the base coordinates (my cells are 9 cm^2) by changing the base unit from meters to cm to no effect. Is there a way to fix this? Am I using a function for tracking that I shouldn’t?

Code building the geometry:

void Scope::BuildDetector()
{
    //---Units: using centimeters as base, i.e. 1 position unit=1 centimeter
    Double_t cm=1;//,micron=1e-4;
    Double_t AdjustX=0*cm,AdjustY=0*cm;

    geom = new TGeoManager("Assemblies", "Geometry using assemblies");
    //--- define some materials
    matVacuum = new TGeoMaterial("Vacuum", 0,0,0);
    pvt=new TGeoMaterial("polyvinyltoluene");
    //--- define some media
    Vacuum = new TGeoMedium("Vacuum",1, matVacuum);
    PVT=new TGeoMedium("PVT",1,pvt);
    top = geom->MakeBox("TOP", Vacuum, 100., 100., 100.);
    geom->SetTopVolume(top);

    fixer=geom->MakeBox("fixer",Vacuum,99.9,99.9,99.9);
    top->AddNode(fixer,1);

    TGeoPgon *hex=new TGeoPgon(0,360,6,2);
    hex->DefineSection(0,-.25*cm,0,.5*cm);//1.611*cm);
    hex->DefineSection(1,.25*cm,0,.5*cm);//1.611*cm);

    Ch4=new TGeoVolume("Ch4",hex,PVT);
    fixer->AddNode(Ch4,1,new TGeoTranslation(-0.3*cm+AdjustX,-0.4*cm+AdjustY,-11.8*cm));

    Ch5=new TGeoVolume("Ch5",hex,PVT);
    fixer->AddNode(Ch5,1,new TGeoTranslation(-0.3*cm+AdjustX,-0.4*cm+AdjustY,-3.0*cm));

    Ch8=new TGeoVolume("Ch8",hex,PVT);
    fixer->AddNode(Ch8,1,new TGeoTranslation(-0.3*cm+AdjustX,-0.4*cm+AdjustY,0.2*cm));

    Ch9=new TGeoVolume("Ch9",hex,PVT);
    fixer->AddNode(Ch9,1,new TGeoTranslation(-0.3*cm+AdjustX,-0.4*cm+AdjustY,8.2*cm));

    geom->CloseGeometry();

}

Here’s the function I’m using for tracking. The tree it’s pulling from contains a vector of the track direction and the z=0 intercept of the track:

vector<vector<Float_t> > Scope::CellHits(UInt_t evt, UInt_t ch)
{
    Double_t cm=1, micron=1e-4;
    TGeoNode *next;
    Double_t startX,startY,startZ,norm,xmag,ymag,zmag,lpoint[3];
    TString cName,DectList[4]={"Ch4_1","Ch5_1","Ch8_1","Ch9_1"};
    vector<UInt_t> hits,trignum;
    vector<Float_t> HitPos(3,0);
    vector<vector<Float_t> > rtnVector;
    CosmicData->GetEntry(evt);
    for(Int_t n=0;n<CurrentEvent->GetNTracks();n++)
    {
        startX=CurrentEvent->GetXPos(n)*micron;
        startY=CurrentEvent->GetYPos(n)*micron;
        startZ=CurrentEvent->GetZPos(n)*micron;
        xmag=CurrentEvent->GetXSlope(n);
        ymag=CurrentEvent->GetYSlope(n);
        zmag=CurrentEvent->GetZSlope(n);
        norm=Sqrt(xmag*xmag+ymag*ymag+zmag*zmag);
        xmag=xmag/norm;
        ymag=ymag/norm;
        zmag=zmag/norm;

        gGeoManager->SetCurrentPoint(-50*xmag*cm+startX,-50*ymag*cm+startY,-50*zmag*cm+startZ);
        gGeoManager->SetCurrentDirection(xmag,ymag,zmag);
        gGeoManager->FindNode();
        next=gGeoManager->FindNextBoundaryAndStep(TGeoShape::Big());
        while(next!=NULL)
        {
            cName=next->GetName();
            if(cName==DectList[ch])
            {                
                gGeoManager->MasterToLocal(gGeoManager->GetCurrentPoint(),lpoint);
                HitPos[0]=lpoint[0];
                HitPos[1]=lpoint[1];
                HitPos[2]=lpoint[2];
                rtnVector.push_back(HitPos);
            }
            next=gGeoManager->FindNextBoundaryAndStep(TGeoShape::Big());
        }
        if(rtnVector.size()==0)
        {
            HitPos[0]=999;
            rtnVector.push_back(HitPos);
        }
    }  
    return rtnVector;

}


Hi,

Your tracking algorithm is correct, this is more or less what raytracing does, but there you do not see such artifact. Likely to be a rounding issue since you use a float array to store the hits positions. Try using double precision everywhere.

Cheers,

Hello again,

Thanks for the reply. I tried using Double_ts instead, unfortunately it didn’t work. It might not be too much of a problem since it appears that the “corners” have a lower hit density than the rest of the shape. If I must, I can probably write another function to cut them out. Any idea why they don’t show up in ray tracing?

Thanks,
Stephen

EDIT: After a bit more searching, I think this has to do with how my geometry is being voxelized. Is there a way to tell the TGeoManager to use smaller voxels? I can see that there are different options available for the CloseGeometry() function, but I’m having a hard time finding what they all are.

Hi,
The fact that the problem does not appear in raytracing is the proof that there is an artefact somehow in your procedure - that I cannot reproduce since I do not have all your code. This is not due to voxelisation, which is only supposed to optimize the navigation time without affecting numerically the results.

Maybe you can post a simplified version that could reproduce your plot.

Cheers,

Hello again,
Here’s a simplified version of my code that produces the graph I posted. I used TRandom3->Gaus() to create the track z=0 intercepts and just set all the slopes to dx=dy=0, dz=1. The blue line it draws is the size of the cell face. I tested this script by compiling with ACLiC, then typing ScatterPlot(2,100000).

Thanks,
Stephen

#include <TH2.h>
#include <TMath.h>
#include <TClonesArray.h>
#include <TGraph.h>
#include <TLeaf.h>
#include <TCanvas.h>
#include <TFormula.h>
#include <Riostream.h>
#include <TGeoVolume.h>
#include <TGeoPgon.h>
#include <TGeoMatrix.h>
#include <TGeoBBox.h>
#include <TGeoManager.h>
#include <TGeoNavigator.h>
#include <TGeoNode.h>
#include <TVirtualGeoTrack.h>
#include <TRandom3.h>

#include <stdio.h>
#include <stdlib.h>
#include <vector>
#include <exception>

using namespace std;


TGeoManager *geom;
TGeoMaterial *matVacuum,*pvt;
TGeoMedium *Vacuum,*PVT;
TGeoVolume *top,*fixer,*Ch4,*Ch5,*Ch8,*Ch9;
vector<Double_t> xHit,yHit;
TRandom3 *Rand=new TRandom3();

void BuildDetector()
{
    //---Units: using centimeters as base, i.e. 1 position unit=1 centimeter
    Double_t cm=1;//,micron=1e-4;
    Double_t AdjustX=.3*cm,AdjustY=.4*cm;

    geom = new TGeoManager("Assemblies", "Geometry using assemblies");
    //--- define some materials
    matVacuum = new TGeoMaterial("Vacuum", 0,0,0);
    pvt=new TGeoMaterial("polyvinyltoluene");
    //--- define some media
    Vacuum = new TGeoMedium("Vacuum",1, matVacuum);
    PVT=new TGeoMedium("PVT",1,pvt);
    top = geom->MakeBox("TOP", Vacuum, 100., 100., 100.);
    geom->SetTopVolume(top);

    fixer=geom->MakeBox("fixer",Vacuum,99.9,99.9,99.9);
    top->AddNode(fixer,1);

    TGeoPgon *hex=new TGeoPgon(0,360,6,2);
    hex->DefineSection(0,-.25*cm,0,1.611*cm);
    hex->DefineSection(1,.25*cm,0,1.611*cm);

    Ch4=new TGeoVolume("Ch4",hex,PVT);
    fixer->AddNode(Ch4,1,new TGeoTranslation(-0.3*cm+AdjustX,-0.4*cm+AdjustY,-11.8*cm));

    Ch5=new TGeoVolume("Ch5",hex,PVT);
    fixer->AddNode(Ch5,1,new TGeoTranslation(-0.3*cm+AdjustX,-0.4*cm+AdjustY,-3.0*cm));

    Ch8=new TGeoVolume("Ch8",hex,PVT);
    fixer->AddNode(Ch8,1,new TGeoTranslation(-0.3*cm+AdjustX,-0.4*cm+AdjustY,0.2*cm));

    Ch9=new TGeoVolume("Ch9",hex,PVT);
    fixer->AddNode(Ch9,1,new TGeoTranslation(-0.3*cm+AdjustX,-0.4*cm+AdjustY,8.2*cm));

    geom->CloseGeometry();

}

vector<vector<Double_t> > CellHits(UInt_t evt,UInt_t ch)
{
    Double_t cm=1, micron=1e-4;
    TGeoNode *next;
    Double_t startX,startY,startZ,norm,xmag,ymag,zmag,lpoint[3];
    TString cName,DectList[4]={"Ch4_1","Ch5_1","Ch8_1","Ch9_1"};
    vector<Double_t> HitPos(3,0);
    vector<vector<Double_t> > rtnVector;
    Int_t nTracks=1;
    for(Int_t n=0;n<nTracks;n++)
    {
        startX=xHit[evt];
        startY=yHit[evt];
        startZ=0;
        xmag=0;
        ymag=0;
        zmag=1;
        /*norm=TMath::Sqrt(xmag*xmag+ymag*ymag+zmag*zmag);
        xmag=xmag/norm;
        ymag=ymag/norm;
        zmag=zmag/norm;*/

        gGeoManager->SetCurrentPoint(-50*xmag*cm+startX,-50*ymag*cm+startY,-50*zmag*cm+startZ);
        gGeoManager->SetCurrentDirection(xmag,ymag,zmag);
        gGeoManager->FindNode();
        next=gGeoManager->FindNextBoundaryAndStep(TGeoShape::Big());
        while(next!=NULL)
        {
            cName=next->GetName();
            if(cName==DectList[ch])
            {
                gGeoManager->MasterToLocal(gGeoManager->GetCurrentPoint(),lpoint);
                HitPos[0]=lpoint[0];
                HitPos[1]=lpoint[1];
                HitPos[2]=lpoint[2];
                rtnVector.push_back(HitPos);
            }
            next=gGeoManager->FindNextBoundaryAndStep(TGeoShape::Big());
        }
        if(rtnVector.size()==0)
        {
            HitPos[0]=999;
            rtnVector.push_back(HitPos);
        }
    }
    return rtnVector;

}

void ScatterPlot(Int_t XCh,UInt_t TotalEvents)
{
    for(UInt_t m=0;m<TotalEvents;m++)
    {
        xHit.push_back(Rand->Gaus());
        yHit.push_back(Rand->Gaus());
    }
    BuildDetector();
    Double_t cm=1;
    TH2D *Scatter=new TH2D("ScatterPlot","Cell hit positions",40,-2,2,40,-2,2);
    vector<vector<Double_t> > hitPos;
    Int_t ChNumSwitch[4]={4,5,8,9};

    for(UInt_t m=0;m<TotalEvents;m++)
    {
        hitPos=CellHits(m,XCh);
        for(UInt_t t=0;t<hitPos.size();t++)
        {
            if(hitPos[t][0]==999)continue;
            Scatter->Fill(hitPos[t][0],hitPos[t][1]);
        }
    }
    Scatter->SetXTitle(Form("Channel %i X position (cm)",ChNumSwitch[XCh]));
    Scatter->SetYTitle(Form("Channel %i Y position (cm)",ChNumSwitch[XCh]));
    TCanvas *ScatterCanvas=new TCanvas();
    ScatterCanvas->SetWindowSize(900+(900-ScatterCanvas->GetWw()),900+(900-ScatterCanvas->GetWh()));
    Scatter->Draw();
    TGraph *outline=new TGraph();
    outline->SetPoint(0,-1.86*cm,0);
    outline->SetPoint(1,-.93*cm,1.611*cm);
    outline->SetPoint(2,.93*cm,1.611*cm);
    outline->SetPoint(3,1.86*cm,0);
    outline->SetPoint(4,.93*cm,-1.611*cm);
    outline->SetPoint(5,-.93*cm,-1.611*cm);
    outline->SetPoint(6,-1.86*cm,0);
    outline->SetLineColor(4);
    outline->Draw("L");
}


Hi again Stephen,

As I suspected, there is an artifact, this is not from your code but rather coming from improper usage of TH2, which is not supposed to represent 2D points with precision and makes a binning. I just filled a 3D polymarker with your collected data and plotted it on top of the chamber geometry to make this clear. See attached plots.

Cheers,
Andrei





test1.C (5.65 KB)

Hello Andrei,

Thanks a lot for the help. I had just assumed that since the TH2 scatter option plotted points, that binning was ignored. I’ll try using a TGraph instead.

Thanks again,
Stephen