void dpwGds() {
  gROOT->SetStyle("Plain");

  TCanvas *c1=new TCanvas("c1","Canvas",800,600);
  gStyle->SetOptStat(0);
  gStyle->SetPalette(1);
  //gStyle->SetPadLeftMargin(0.12);
  //gStyle->SetPadRightMargin(0.12);
  //c1->UseCurrentStyle();

  TH2F *h=new TH2F("label",";MIP x position (#mum);MIP y position (#mum);Fraction of charge",
		   31,-77.5,77.5,31,-77.5,77.5);

  h->GetXaxis()->SetTitleOffset(1.3);
  h->GetYaxis()->SetTitleOffset(1.5);
  h->GetZaxis()->SetTitleOffset(1.3);

  ifstream fin("GDS_DPW.dat");

  double p,s;
  fin >> p >> s;
  assert(p==0.050);
  assert(s==0.012);

  unsigned np(0);

  for(unsigned i(0);i<21;i++) {
    for(unsigned j(0);j<25;j++) {
      unsigned jx(j/5);
      unsigned jy(j%5);
      
      fin >> p;


      double cx(0.0);double cy(0.0); unsigned nk(1);
      if(i== 1) nk=4;
      if(i== 2) nk=4;
      if(i== 3) nk=4;
      if(i== 4) nk=8;
      if(i== 5) nk=4;
      if(i== 6) nk=4;
      if(i== 7) nk=8;
      if(i== 8) nk=8;
      if(i== 9) nk=4;
      if(i==10) nk=4;
      if(i==11) nk=8;
      if(i==12) nk=8;
      if(i==13) nk=8;
      if(i==14) nk=4;
      if(i==15) nk=2;
      if(i==16) nk=4;
      if(i==17) nk=4;
      if(i==18) nk=4;
      if(i==19) nk=4;
      if(i==20) nk=1;

      for(unsigned k(0);k<nk;k++) {
	if(i== 1) {
	  if(k==0) {cx=  5.0;cy=  0.0;}
	  if(k==1) {cx= -5.0;cy=  0.0;}
	  if(k==2) {cx=  0.0;cy=  5.0;}
	  if(k==3) {cx=  0.0;cy= -5.0;}
	}
	if(i== 2) {
	  if(k==0) {cx=  5.0;cy=  5.0;}
	  if(k==1) {cx= -5.0;cy=  5.0;}
	  if(k==2) {cx=  5.0;cy= -5.0;}
	  if(k==3) {cx= -5.0;cy= -5.0;}
	}
	if(i== 3) {
	  if(k==0) {cx= 10.0;cy=  0.0;}
	  if(k==1) {cx=-10.0;cy=  0.0;}
	  if(k==2) {cx=  0.0;cy= 10.0;}
	  if(k==3) {cx=  0.0;cy=-10.0;}
	}
	if(i== 5) {
	  if(k==0) {cx= 10.0;cy= 10.0;}
	  if(k==1) {cx=-10.0;cy= 10.0;}
	  if(k==2) {cx= 10.0;cy=-10.0;}
	  if(k==3) {cx=-10.0;cy=-10.0;}
	}
	if(i== 6) {
	  if(k==0) {cx= 15.0;cy=  0.0;}
	  if(k==1) {cx=-15.0;cy=  0.0;}
	  if(k==2) {cx=  0.0;cy= 15.0;}
	  if(k==3) {cx=  0.0;cy=-15.0;}
	}
	if(i== 9) {
	  if(k==0) {cx= 15.0;cy= 15.0;}
	  if(k==1) {cx=-15.0;cy= 15.0;}
	  if(k==2) {cx= 15.0;cy=-15.0;}
	  if(k==3) {cx=-15.0;cy=-15.0;}
	}
	if(i==10) {
	  if(k==0) {cx= 20.0;cy=  0.0;}
	  if(k==1) {cx=-20.0;cy=  0.0;}
	  if(k==2) {cx=  0.0;cy= 20.0;}
	  if(k==3) {cx=  0.0;cy=-20.0;}
	}
	if(i==14) {
	  if(k==0) {cx= 20.0;cy= 20.0;}
	  if(k==1) {cx=-20.0;cy= 20.0;}
	  if(k==2) {cx= 20.0;cy=-20.0;}
	  if(k==3) {cx=-20.0;cy=-20.0;}
	}
	if(i==15) {
	  if(k==0) {cx=  0.0;cy=-25.0;}
	  if(k==1) {cx=-25.0;cy=  0.0;}
	}
	if(i==20) {
	  if(k==0) {cx=-25.0;cy=-25.0;}
	}
	
	if(i== 2) {cx= 5.0;cy= 5.0;}
	if(i== 4) {cx=10.0;cy= 5.0;}
	if(i== 5) {cx=10.0;cy=10.0;}
	if(i== 7) {cx=15.0;cy= 5.0;}
	if(i== 8) {cx=15.0;cy=10.0;}
	if(i== 9) {cx=15.0;cy=15.0;}
	if(i==11) {cx=20.0;cy= 5.0;}
	if(i==12) {cx=20.0;cy=10.0;}
	if(i==13) {cx=20.0;cy=15.0;}
	if(i==14) {cx=20.0;cy=20.0;}
	if(i==16) {cx=25.0;cy= 5.0;}
	if(i==17) {cx=25.0;cy=10.0;}
	if(i==18) {cx=25.0;cy=15.0;}
	if(i==19) {cx=25.0;cy=20.0;}
	if(i==20) {cx=25.0;cy=25.0;}

	if(jx>=1 && jx<=3 && jy>=1 && jy<=3) {
	  np++;
	  h->Fill(cx+50.0*(jx-2.0),cy+50.0*(2.0-jy),p);
	}
      }
    }
  }

  cout << "Number of points = " << np << endl;

  h->GetZaxis()->SetRangeUser(0,0.58);
  h->Draw("lego2");
  //h->Draw("colz");

  c1->Update();
  c1->Print("../dpsGds.png");
}