{

// Dimensions of cylinder
Double_t a=0.8.,b=2.5.,h=2.5; 
Int_t nDiscrCyl = 10;//3;

//binwidth
Double_t dr = (b-a)/(nDiscrCyl-1);
Double_t dphi = TMath::Pi()*2/nDiscrCyl;
Double_t dz = h/(nDiscrCyl-1);

TH3F *h3cyl = new TH3F("h3","test", nDiscrCyl,a-dr/2,b+dr/2, nDiscrCyl+1,0-dphi/2,TMath::Pi()*2+dphi/2, nDiscrCyl,0-dz/2,h+dz/2);
TRandom3 rand;
for (int i=0;i<100000;i++) {
  double r = rand.Gaus(1.5,0.2);
  double phi = rand.Gaus(3.,1.5);
  double z = rand.Gaus(1.4,0.2);
  h3cyl->Fill(r,phi,z);
} 

h3cyl->Draw("glbox");
cout<<h3cyl->Interpolate(1.5,2,0.1)<<endl;



// REORDERING INTO CARTESIAN COORDINATES
Int_t nDiscr = 3;//6;
Double_t dx = 2*b/(nDiscr-1)/2; 
Double_t dz = h/(nDiscr-1)/2; 

TH3F *h3cart = new TH3F("h3cart","test1",nDiscr,-b-dx,b+dx,nDiscr,-b-dx,b+dx,nDiscr,-dz,h+dz);
for (Int_t zi = 0; zi<nDiscr; zi++) {
  for (Int_t yi = 0; yi<nDiscr; yi++) {
    for (Int_t xi = 0; xi<nDiscr; xi++) {
      
      Double_t xp = (xi)*2*b/(nDiscr-1)-b;
      Double_t yp = (yi)*2*b/(nDiscr-1)-b;
      Double_t zp = (zi)*h/(nDiscr-1);
      
      Double_t r = TMath::Sqrt(xp*xp+yp*yp);
      Double_t phi = TMath::ATan2(yp,xp);
      if (phi<0) phi += TMath::Pi()*2;
      if (phi>2*TMath::Pi()) phi -= TMath::Pi()*2;


      //if ( r>=b || r<=a || zp<=0 || zp>=h ) { // within domain
      if ( r>b || r<a || zp<0 || zp>h ) { // within domain + border

	h3cart->SetBinContent(xi+1,yi+1,zi+1,0);
      } else {
	// cout<<r<<" "<<phi<<" "<<zp<<" | "<<xp<<" "<<yp<<" "<<zp<<endl;
	h3cart->SetBinContent(xi+1,yi+1,zi+1,h3cyl->Interpolate(r,phi,zp));
      }
    }
  }
}
h3cart->Draw("glbox");
}