#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 <TPolyMarker3D.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();
    TPolyMarker3D *marker;
    

    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)
    {
        marker = new TPolyMarker3D(TotalEvents);
        marker->SetMarkerColor(kRed);
        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;
                marker->SetNextPoint(hitPos[t][0],hitPos[t][1], hitPos[t][2]);
                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");
        TCanvas *geomCanvas=new TCanvas();
        TString VolList[4]={"Ch4","Ch5","Ch8","Ch9"};
        gGeoManager->GetVolume(VolList[XCh])->Draw();
        marker->Draw("SAME");
        
    }