/*
  g++ -Wall -std=c++0x -I. code/CalibrationTest.cc -o CalibrationTest.exe `root-config --libs --cflags`
*/

/*
 0 = SCs
 1 = SCs with perfect ADC

 0 = All TCs
 0 = True energy TCs
 1 = All TCs with cell centring
 2 = Threshold TCs
 3 = Equal share (0 bits)
 4 = Equal share if > 0 (4 bits)
 5 = SuperTCs (2 bits)

 10-17 = Fractions 1-8 bits/TC (5-26 bits)
 18 = HGCROC-like implementation of fraction 1 bit
 19 = HGCROC-like implementation of fraction 2 bits

 20-25 = Compression 8-3 bits/TC, drop LSBs (16-32 bits)
 30-35 = Compression 8-3 bits/TC, saturate (16-32 bits)
 
 40 = Fractions 1 bit/TC + 4E/5M
 41 = Fractions 2 bits/TC + 3E/3M
 42 = Fractions 2 bits/TC + 4E/2M
 43 = Fractions 2 bits/TC + 5E/1M


 Work to limits when compressing:
 SC 25 x 2^12 < 2^5 x 2^12 = 2^17
 NTC lo = 4 x 25 x 2^12 = 25 x 2^14 < 2^5 x 2^15 = 2^19
 NTC hi = 9 x 25 x 2^12 = 225 x 2^12 < 2^8 x 2^12 = 2^20
 CTC = 4 x 225 x 2^12 = 225 x 2^14 < 2^8 x 2^14 = 2^22
 DONE ! Do two TOT ranges properly
 DONE ! Fix Gaussian to integrate properly
 Make PU from other Gaussians plus model as density not per sensor cell
 Set up various positions lo/hi radius for both densities
 DONE! Change numbers of cells to work for both 2.2 and 3.3 cm sizes
 Make sure we catch equal edges of fractions i.e. 1/4, 1/2, etc. correctly.
 Connect up fSinTheta with fEta

*/


#include "Hack.hh"

#include <signal.h>

#include <cassert>
#include <iomanip>
#include <fstream>
#include <sstream>

#include "TFile.h"
#include "TF1.h"
#include "TH1D.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TProfile.h"
//#include "TMath.h"
#include "TE1F.hh"
#include "TH1A.hh"
#include "Random.hh"
#include "PackerFloat22.hh"
#include "TFileHandler.hh"

const unsigned kNumberOfRadii(24);

double cutRadius(unsigned r) {
  assert(r<kNumberOfRadii);
  //return 0.5*(r+6);
  //return 0.0025*(r+4);
  return 0.06*(r+0.5)/kNumberOfRadii;
}

class PhotonModel {
public:
  PhotonModel() {
  }
  
  PhotonModel(unsigned m, double a) : fModel(m), fA(a), fNorm(1.0/(2.0*fA*fA*acos(-1.0))) {
  }

  double density(double r) const {
    if(fModel==2) return fNorm*exp(-0.5*r*r/(fA*fA));
    if(fModel==1) return fNorm*exp(-sqrt(r/fA))/12.0;
    return fNorm*exp(-r/fA);
  }

private:
  unsigned fModel;
  double fA;
  double fNorm;
};


class Array {
public:
  enum {
    kScale=4,
    kHiDensityCtcs=3*kScale,
    kLoDensityCtcs=2*kScale,
    kHiDensityNtcs=2*kHiDensityCtcs,
    kLoDensityNtcs=2*kLoDensityCtcs,
    kHiDensityScs=3*kHiDensityNtcs,
    kLoDensityScs=3*kLoDensityNtcs,


    //kCoarseCells=10,
    //kTriggerCells=2*kCoarseCells,
    //kHiDensityCells=3*kTriggerCells,
    //kLoDensityCells=2*kTriggerCells,
  };
    
  enum ScMethods {
    kScPerfect,
    kScReal,

    kScTrue,
    kNumberOfScMethods
  };

  enum NtcMethods {
    kNtcReal,
    kNtcEqualShare,
    kNtcThreshold,
    kNtcThresholdEqualShare,
    kNtcSuperTcs,

    kNtcFractions, // 5-12
    kNtcCompressions=kNtcFractions+8, // 13-20
    kNtcElinks=kNtcCompressions+8, // 21-24
    kNtcDropBits=kNtcElinks+4, // 25-32
    
    kNtcTrue=kNtcDropBits+8, // 33
    kNumberOfNtcMethods

    //kNumberOfCtcMethods=4,
    //kNumberOfHtcMethods=1
  };

  
  Array(unsigned p=0, bool e=true, bool pu=true) : fHiDensity(p==0), fEnergy(e), fPileup(pu), fEta(fHiDensity?2.9:1.6),
						   fNumberOfScs(fHiDensity?kHiDensityScs:kLoDensityScs),
						   fNumberOfNtcs(fHiDensity?kHiDensityNtcs:kLoDensityNtcs),
						   fNumberOfCtcs(fHiDensity?kHiDensityCtcs:kLoDensityCtcs),
						   
						   //fTriggerCellSide(sqrt(4.908)),
						   fTriggerCellSide(fHiDensity?2.2:3.3),
						   fCoarseCellSide(2.0*fTriggerCellSide),
						   fSensorCellSide(fTriggerCellSide/3.0),
						   fArraySide(fTriggerCellSide*fNumberOfNtcs),
						   fArrayHalfSide(0.5*fArraySide),
						   
						   fNoise(fHiDensity?1/2.0:1/6.0), fMip(10.0), fSinTheta(1.0/cosh(fEta)),
						   fPacker5E1M(5,1),fPacker3E3M(3,3),fPacker4E6M(4,6),fPacker4E5M(4,5),
						   fPacker4E4M(4,4),fPacker4E3M(4,3),fPacker4E2M(4,2),fPacker4E1M(4,1) {
    std::cout << "Array::ctor() place = " << p
	      << ", high density = " << (fHiDensity?"true":"false")
	      << ", energy = " << (fEnergy?"true":"false")
	      << ", pileup = " << (fPileup?"true":"false")
	      << ", eta = " << fEta
	      << ", sin(theta) = " << fSinTheta << std::endl;

    std::cout << " kNtcEqualShare = " << kNtcEqualShare << std::endl;
    std::cout << " kNtcThreshold = " << kNtcThreshold << std::endl;
    std::cout << " kNtcThresholdEqualShare = " << kNtcThresholdEqualShare << std::endl;
    std::cout << " kNtcSuperTcs = " << kNtcSuperTcs << std::endl;
    std::cout << " kNtcFractions = " << kNtcFractions << std::endl;
    std::cout << " kNtcCompressions = " << kNtcCompressions << std::endl;
    std::cout << " kNtcElinks = " << kNtcElinks << std::endl;
    std::cout << " kNtcDropBits = " << kNtcDropBits << std::endl;
    std::cout << " kNtcTrue = " << kNtcTrue << std::endl;


    
    hScAdcConversion=new TH2F("ArrayScAdcConversion",";Original ADC value;Final ADC value",1000,0.0,100000.0,1000,0.0,100000.0);
    hScAdcDeltaConversion=new TH2F("ArrayScAdcDeltaConversion",";Original ADC value;Final-Original ADC value",1000,0.0,100000.0,1000,-500.0,500.0);

    for(unsigned m(0);m<kNumberOfNtcMethods;m++) {
      std::ostringstream sout;
      sout << "ArrayNtcMethod" << std::setw(2) << std::setfill('0') << m;
      
      hDeltaE1[m]=new TH2F((sout.str()+"DeltaE1").c_str(),
			   ";Total ADC sum;Delta E/E",100,1000.0,1100.0,100,-0.25,0.25);
      hDeltaE2[m]=new TH2F((sout.str()+"DeltaE2").c_str(),
			   ";Total ADC sum;Delta E/E",100,0.0,20000.0,100,-0.25,0.25);

      hNtcAdcConversion0[m]=new TH2F((sout.str()+"AdcConversion0").c_str(),
				     ";Original ADC value;Final ADC value",1000,0.0,1000.0,1000,0.0,1000.0);
      hNtcAdcConversion1[m]=new TH2F((sout.str()+"AdcConversion1").c_str(),
				    ";Original ADC value;Final ADC value",1000,0.0,50000.0,1000,0.0,50000.0);
      hNtcAdcConversion2[m]=new TH2F((sout.str()+"AdcConversion2").c_str(),
				    ";Original ADC value;Final ADC value",1000,0.0,500000.0,1000,0.0,500000.0);
      hNtcAdcDeltaConversion0[m]=new TH2F((sout.str()+"AdcDeltaConversion0").c_str(),
					 ";Original ADC value;Final-Original ADC value",1000,0.0,1000.0,1000,-5000.0,5000.0);
      hNtcAdcDeltaConversion1[m]=new TH2F((sout.str()+"AdcDeltaConversion1").c_str(),
					 ";Original ADC value;Final-Original ADC value",1000,0.0,50000.0,1000,-5000.0,5000.0);
      hNtcAdcDeltaConversion2[m]=new TH2F((sout.str()+"AdcDeltaConversion2").c_str(),
					 ";Original ADC value;Final-Original ADC value",1000,0.0,500000.0,1000,-5000.0,5000.0);
      hNtcAdcFracConversion0[m]=new TH2F((sout.str()+"AdcFracConversion0").c_str(),
					";Original ADC value;Final-Original/Original ADC value",1000,0.0,1000.0,1000,-1.0,1.0);
      hNtcAdcFracConversion1[m]=new TH2F((sout.str()+"AdcFracConversion1").c_str(),
					";Original ADC value;Final-Original/Original ADC value",1000,0.0,50000.0,1000,-1.0,1.0);
      hNtcAdcFracConversion2[m]=new TH2F((sout.str()+"AdcFracConversion2").c_str(),
					";Original ADC value;Final-Original/Original ADC value",1000,0.0,500000.0,1000,-1.0,1.0);

      hNtcAdcChiSq0[m]=new TH1F((sout.str()+"AdcChiSq0").c_str(),
				";ChiSq0",100,0.0,1000000.0);
      hNtcAdcChiSq1[m]=new TH1F((sout.str()+"AdcChiSq1").c_str(),
				";ChiSq1",100,0.0,1000000.0);


    numChiSq0[m]=0.0;
    sumChiSq0[m]=0.0;
    numChiSq1[m]=0.0;
    sumChiSq1[m]=0.0;

    }

    hNtcAdcChiSqAvg0=new TH1F("ArrayAdcChiSqAvg0",
			      ";ChiSqAvg0",kNumberOfNtcMethods,0,kNumberOfNtcMethods);
    hNtcAdcChiSqAvg1=new TH1F("ArrayAdcChiSqAvg1",
			      ";ChiSqAvg1",kNumberOfNtcMethods,0,kNumberOfNtcMethods);

  }

  unsigned numberOfSensorCells() const {
    return fNumberOfScs;
  }

  unsigned numberOfTriggerCells() const {
    return fNumberOfNtcs;
  }

  unsigned numberOfCoarseCells() const {
    return fNumberOfCtcs;
  }
  /*
  unsigned numberOfHgcrocXCells() const {
    return (fHiDensity?kCoarseCells:kCoarseCells/2);
  }

  unsigned numberOfHgcrocYCells() const {
    return kCoarseCells/2;
  }
  */
  double sensorCellPosition(unsigned i) const {
    assert(i<fNumberOfScs);
    return fSensorCellSide*(i+0.5-0.5*fNumberOfScs);
  }

  unsigned sensorCellNumber(double p) const {
    assert(fabs(p)<fArrayHalfSide);
    unsigned i((p+fArrayHalfSide)/fSensorCellSide);
    assert(i<fNumberOfScs);
    return i;
  }

  double triggerCellPosition(unsigned i) const {
    assert(i<fNumberOfNtcs);
    return fTriggerCellSide*(i+0.5-0.5*fNumberOfNtcs);
  }

  unsigned triggerCellNumber(double p) const {
    assert(fabs(p)<fArrayHalfSide);
    unsigned i((p+fArrayHalfSide)/fTriggerCellSide);
    assert(i<fNumberOfNtcs);
    return i;
  }
  /*
  double coarseCellPosition(unsigned i) const {
    assert(i<numberOfCoarseCells());
    return fCoarseCellSide*(i+0.5-0.5*numberOfCoarseCells());
  }

  double hgcrocXCellPosition(unsigned i) const {
    assert(i<numberOfHgcrocXCells());
    if(fHiDensity) return fCoarseCellSide*(i+0.5-0.5*numberOfHgcrocXCells());
    return 2.0*fCoarseCellSide*(i+0.5-0.5*numberOfHgcrocXCells());
  }

  double hgcrocYCellPosition(unsigned i) const {
    assert(i<numberOfHgcrocYCells());
    return 2.0*fCoarseCellSide*(i+0.5-0.5*numberOfHgcrocYCells());
  }
  */
  void initialise() {
    for(unsigned ix(0);ix<numberOfSensorCells();ix++) {
      for(unsigned iy(0);iy<numberOfSensorCells();iy++) {
	fNoiseArray[ix][iy]=Random::random().Gaus(0.0,fNoise);

	if(!fPileup || Random::random().Uniform()>(fHiDensity?0.5:0.05)) fPileupArray[ix][iy]=0.0;
	else fPileupArray[ix][iy]=0.5+Random::random().Exp(5.0);

	fEnergyArray[ix][iy]=0.0;
      }
    }
  }
    
  void increment(unsigned ix, unsigned iy, double e) {
    assert(ix<numberOfSensorCells());
    assert(iy<numberOfSensorCells());
    if(fEnergy) fEnergyArray[ix][iy]+=e;
  }

  void setFractions(unsigned ix, unsigned iy, unsigned m, unsigned a=3) {
    assert(ix<numberOfCoarseCells());
    assert(iy<numberOfCoarseCells());

    unsigned eMax(0);
    unsigned sum(0);
    unsigned mx,my;
    
    for(unsigned jx(2*ix);jx<2*(ix+1);jx++) {
      for(unsigned jy(2*iy);jy<2*(iy+1);jy++) {
	sum+=fNtcAdc[0][jx][jy];
	if(eMax<=fNtcAdc[kNtcReal][jx][jy]) {
	  eMax=fNtcAdc[kNtcReal][jx][jy];
	  mx=jx;
	  my=jy;
	}
      }
    }

    unsigned nx((mx%2)==0?mx+1:mx-1);
    unsigned ny((my%2)==0?my+1:my-1);

    unsigned fraction;

    if(m<8) {
      unsigned bins(1<<(m+1));
      if(m>=4) bins=(1<<(m-3))-1;
      
      for(unsigned jx(2*ix);jx<2*(ix+1);jx++) {
	for(unsigned jy(2*iy);jy<2*(iy+1);jy++) {
	  if(!(jx==mx && jy==my)) {
	    if(fNtcAdc[kNtcReal][jx][jy]==0) fraction=0;
	    else fraction=(2*bins*fNtcAdc[kNtcReal][jx][jy]-1)/sum;
	    
	    assert(fraction<bins);
	    //if(fraction==bins) {
	    //assert(2*fNtcAdc[kNtcReal][jx][jy]==sum);
	    //fraction=bins-1;
	    //}

	    fNtcAdc[a][jx][jy]=((1+2*fraction)*sum+2*bins)/(4*bins);
	    if(m>=4 && fraction==0) fNtcAdc[a][jx][jy]=0;
	    
	    hDeltaE1[kNtcFractions+m]->Fill(sum,(double(fNtcAdc[kNtcFractions+m][jx][jy])-double(fNtcAdc[kNtcReal][jx][jy]))/sum);
	    hDeltaE2[kNtcFractions+m]->Fill(sum,(double(fNtcAdc[kNtcFractions+m][jx][jy])-double(fNtcAdc[kNtcReal][jx][jy]))/sum);
	  }
	}
      }
      
    } else if(m==8) {
      unsigned sum2=sum;

      //std::cout << "Calling 4E5M with data = " << (sum+1)/2 << std::endl;
      fPacker4E5M.pack((sum+1)/2);
      sum=2*fPacker4E5M.unpack();

      // Fractions 1 bit/TC + 4E/5M
      unsigned bins(2);      
      for(unsigned jx(2*ix);jx<2*(ix+1);jx++) {
	for(unsigned jy(2*iy);jy<2*(iy+1);jy++) {
	  if(!(jx==mx && jy==my)) {
	    if(fNtcAdc[kNtcReal][jx][jy]==0) fraction=0;
	    else fraction=(2*bins*fNtcAdc[kNtcReal][jx][jy]-1)/sum2;	    
	    assert(fraction<bins);

	    fNtcAdc[a][jx][jy]=((1+2*fraction)*sum+2*bins)/(4*bins);
	  }
	}
      }

    } else if(m==9) {
      unsigned sum2=sum;

      //std::cout << "Calling 3E3M with sum " << sum << ", data = " << (sum+1024)/2048 << std::endl;
      fPacker3E3M.pack((sum+1024)/2048);
      sum=2048*fPacker3E3M.unpack();
      //std::cout << "Called 3E3M with output sum = " << sums << std::endl;

      // Fractions 2 bits/TC + 3E/3M
      unsigned bins(4);      
      for(unsigned jx(2*ix);jx<2*(ix+1);jx++) {
	for(unsigned jy(2*iy);jy<2*(iy+1);jy++) {
	  if(!(jx==mx && jy==my)) {
	    if(fNtcAdc[kNtcReal][jx][jy]==0) fraction=0;
	    else fraction=(2*bins*fNtcAdc[kNtcReal][jx][jy]-1)/sum2;	    
	    assert(fraction<bins);

	    fNtcAdc[a][jx][jy]=((1+2*fraction)*sum+2*bins)/(4*bins);
	  }
	}
      }

    } else if(m==10) {
      unsigned sum2=sum;
      
      //std::cout << "Calling 4E2M with data = " << (sum+8)/16 << std::endl;
      fPacker4E2M.pack((sum+8)/16);
      sum=16*fPacker4E2M.unpack();

      // Fractions 2 bits/TC + 4E/2M
      unsigned bins(4);      
      for(unsigned jx(2*ix);jx<2*(ix+1);jx++) {
	for(unsigned jy(2*iy);jy<2*(iy+1);jy++) {
	  if(!(jx==mx && jy==my)) {
	    if(fNtcAdc[kNtcReal][jx][jy]==0) fraction=0;
	    else fraction=(2*bins*fNtcAdc[kNtcReal][jx][jy]-1)/sum2;	    
	    assert(fraction<bins);

	    fNtcAdc[a][jx][jy]=((1+2*fraction)*sum+2*bins)/(4*bins);
	  }
	}
      }

    } else if(m==11) {
      unsigned sum2=sum;

      //std::cout << "Calling 5E1M with data = " << sum << std::endl;
      fPacker5E1M.pack(sum);
      sum=fPacker5E1M.unpack();

      // Fractions 2 bits/TC + 5E/1M
      unsigned bins(4);      
      for(unsigned jx(2*ix);jx<2*(ix+1);jx++) {
	for(unsigned jy(2*iy);jy<2*(iy+1);jy++) {
	  if(!(jx==mx && jy==my)) {
	    if(fNtcAdc[kNtcReal][jx][jy]==0) fraction=0;
	    else fraction=(2*bins*fNtcAdc[kNtcReal][jx][jy]-1)/sum2;	    



	    //std::cout << "Orig " << fNtcAdc[0][jx][jy] << " sum " << sum << " bins " << bins << " fraction " << fraction << std::endl;
	   


	    assert(fraction<bins);
	    fNtcAdc[a][jx][jy]=((1+2*fraction)*sum+2*bins)/(4*bins);
	  }
	}
      }

    } else {
      assert(false);
    }

    fNtcAdc[a][mx][my]=sum-fNtcAdc[a][nx][my]-fNtcAdc[a][mx][ny]-fNtcAdc[a][nx][ny];

    if(m<8) assert((fNtcAdc[kNtcReal][2*ix  ][2*iy  ]+fNtcAdc[kNtcReal][2*ix+1][2*iy  ]+fNtcAdc[kNtcReal][2*ix  ][2*iy+1]+fNtcAdc[kNtcReal][2*ix+1][2*iy+1])==
		   (fNtcAdc[a][2*ix  ][2*iy  ]+fNtcAdc[a][2*ix+1][2*iy  ]+fNtcAdc[a][2*ix  ][2*iy+1]+fNtcAdc[a][2*ix+1][2*iy+1]));
    
    if(false || fNtcAdc[a][mx][my]>100000) {
      std::cout << "a = " << a << ", E1 = " << fNtcAdc[kNtcReal][mx][my] << " ~ " << fNtcAdc[a][mx][my] << " = " << fNtcAdc[a][2*ix  ][2*iy  ]
		<< ", E2 = " << fNtcAdc[kNtcReal][mx][ny] << " ~ " << fNtcAdc[a][mx][ny] << " = " << fNtcAdc[a][2*ix+1][2*iy  ]
		<< ", E3 = " << fNtcAdc[kNtcReal][nx][my] << " ~ " << fNtcAdc[a][nx][my] << " = " << fNtcAdc[a][2*ix  ][2*iy+1]
		<< ", E4 = " << fNtcAdc[kNtcReal][nx][ny] << " ~ " << fNtcAdc[a][nx][ny] << " = " << fNtcAdc[a][2*ix+1][2*iy+1]
		<< ", sum = " << fNtcAdc[kNtcReal][2*ix  ][2*iy  ]+fNtcAdc[kNtcReal][2*ix+1][2*iy  ]+fNtcAdc[kNtcReal][2*ix  ][2*iy+1]+fNtcAdc[kNtcReal][2*ix+1][2*iy+1]
		<< " ~ " << fNtcAdc[a][2*ix  ][2*iy  ]+fNtcAdc[a][2*ix+1][2*iy  ]+fNtcAdc[a][2*ix  ][2*iy+1]+fNtcAdc[a][2*ix+1][2*iy+1]
		<< " ~ " << sum
		<< std::endl;
      assert(false);
    }
  
    if(m<8) assert((fNtcAdc[kNtcReal][2*ix][2*iy]+fNtcAdc[kNtcReal][2*ix+1][2*iy]+fNtcAdc[kNtcReal][2*ix][2*iy+1]+fNtcAdc[kNtcReal][2*ix+1][2*iy+1])==sum);
    assert((fNtcAdc[a][2*ix][2*iy]+fNtcAdc[a][2*ix+1][2*iy]+fNtcAdc[a][2*ix][2*iy+1]+fNtcAdc[a][2*ix+1][2*iy+1])==sum);
  }

  
  void digitise(unsigned mxx) {
    unsigned sigma3(unsigned(3.0*fNoise*fMip+0.5));
    //std::cout << "sigma3 = " << sigma3 << std::endl;
    unsigned halfMip(unsigned(0.5*fMip+0.5));
    //std::cout << "halfMip = " << halfMip << std::endl;
    unsigned nGap(300);

    // Sensor cells

    for(unsigned ix(0);ix<fNumberOfScs;ix++) {
      for(unsigned iy(0);iy<fNumberOfScs;iy++) {
	fScAdc[kScTrue][ix][iy]=0;
	
	double e(fNoiseArray[ix][iy]+fPileupArray[ix][iy]+fEnergyArray[ix][iy]);

	if(e<0.0) {
	  for(unsigned m(0);m<kNumberOfScMethods;m++) fScAdc[m][ix][iy]=0;

	} else {
	  fScAdc[kScPerfect][ix][iy]=unsigned(fMip*e+0.5);
	  fScAdc[kScReal][ix][iy]=fScAdc[kScPerfect][ix][iy];

	  // ADC-TOT transition
	  if(fScAdc[kScReal][ix][iy]<1000+2*nGap) {

	    // Catch transition region
	    if(fScAdc[kScReal][ix][iy]>1000) fScAdc[kScReal][ix][iy]=1000+nGap;
	    
	  } else {
	    
	    // Do two-range 9-bit TOT
	    fScAdc[kScReal][ix][iy]=(fScAdc[kScReal][ix][iy]+13)/25;
	    if(fScAdc[kScReal][ix][iy]<0x200) {
	      fScAdc[kScReal][ix][iy]*=25;
	    } else {
	      fScAdc[kScReal][ix][iy]=(fScAdc[kScReal][ix][iy]+2)/4;
	      if(fScAdc[kScReal][ix][iy]>0x1ff) fScAdc[kScReal][ix][iy]=0x1ff;
	      fScAdc[kScReal][ix][iy]*=100;

	      //std::cout << "ix,iy = " << ix << ", " << iy << " ADC = " << fScAdc[kScReal][ix][iy] << std::endl;
	    }
	  }

	  hScAdcConversion->Fill(fScAdc[kScPerfect][ix][iy],fScAdc[kScReal][ix][iy]);
	  hScAdcDeltaConversion->Fill(fScAdc[kScPerfect][ix][iy],int(fScAdc[kScReal][ix][iy])-int(fScAdc[kScPerfect][ix][iy]));

	  // Thresholds
	  for(unsigned m(0);m<kNumberOfScMethods;m++) {
	    if(fScAdc[m][ix][iy]<halfMip) fScAdc[m][ix][iy]=0;
	  }
	}
      }
    }
    
    // Trigger cells
    
    unsigned nScSum(0);

    unsigned mipT2(unsigned(0.5+2.0*fMip/fSinTheta));
    //std::cout << "MipT2 = " << mipT2 << std::endl;
    
    for(unsigned ix(0);ix<fNumberOfNtcs;ix++) {
      for(unsigned iy(0);iy<fNumberOfNtcs;iy++) {
	fNtcAdc[kNtcTrue][ix][iy]=0;
	
	// Do perfect first
	fNtcAdc[kNtcReal][ix][iy]=0;
	
	for(unsigned jx(3*ix);jx<3*(ix+1);jx++) {
	  for(unsigned jy(3*iy);jy<3*(iy+1);jy++) {
	    if(fScAdc[kScReal][jx][jy]>=sigma3) {
	      nScSum+=fScAdc[kScReal][jx][jy];

	      //if(fScAdc[kScReal][jx][jy]>25*0x200) std::cout << "Summing ix,iy = " << jx << ", " << jy << " ADC = " << fScAdc[kScReal][jx][jy] << std::endl;


	      fNtcAdc[kNtcReal][ix][iy]+=fScAdc[kScReal][jx][jy];
	    }
	  }
	}
	
	for(unsigned m(0);m<kNumberOfNtcMethods;m++) {
	  if(m!=kNtcReal) {

	    // Do thresholds
	    if(fNtcAdc[kNtcReal][ix][iy]>=mipT2) {
	      fPacker4E3M.pack((fNtcAdc[kNtcReal][ix][iy]+1)/2);
	      fNtcAdc[kNtcThreshold][ix][iy]=2*fPacker4E3M.unpack();
	    } else {
	      fNtcAdc[kNtcThreshold][ix][iy]=0;
	    }

	    // Do packing
	    fPacker4E4M.pack(fNtcAdc[kNtcReal][ix][iy]);
	    fNtcAdc[kNtcCompressions  ][ix][iy]=fPacker4E4M.unpack();
	    fPacker4E3M.pack(fNtcAdc[kNtcReal][ix][iy]/2);
	    fNtcAdc[kNtcCompressions+1][ix][iy]=2*fPacker4E3M.unpack();
	    fPacker4E2M.pack((fNtcAdc[kNtcReal][ix][iy]+2)/4);
	    fNtcAdc[kNtcCompressions+2][ix][iy]=4*fPacker4E2M.unpack();
	    fPacker4E1M.pack((fNtcAdc[kNtcReal][ix][iy]+4)/8);
	    fNtcAdc[kNtcCompressions+3][ix][iy]=8*fPacker4E1M.unpack();

	    fPacker4E4M.pack(fNtcAdc[kNtcReal][ix][iy]&fPacker4E4M.dataMax());
	    fNtcAdc[kNtcCompressions+4][ix][iy]=fPacker4E4M.unpack();
	    fPacker4E3M.pack(fNtcAdc[kNtcReal][ix][iy]&fPacker4E3M.dataMax());
	    fNtcAdc[kNtcCompressions+5][ix][iy]=fPacker4E3M.unpack();
	    fPacker4E2M.pack(fNtcAdc[kNtcReal][ix][iy]&fPacker4E2M.dataMax());
	    fNtcAdc[kNtcCompressions+6][ix][iy]=fPacker4E2M.unpack();
	    fPacker4E1M.pack(fNtcAdc[kNtcReal][ix][iy]&fPacker4E1M.dataMax());
	    fNtcAdc[kNtcCompressions+7][ix][iy]=fPacker4E1M.unpack();
	    /*
	    std::cout << "Orig = " << fNtcAdc[kNtcReal][ix][iy]
		      << ", packers = " << fNtcAdc[kNtcCompressions  ][ix][iy]
		      << ", "  << fNtcAdc[kNtcCompressions+1][ix][iy]
		      << ", "  << fNtcAdc[kNtcCompressions+2][ix][iy]
		      << ", "  << fNtcAdc[kNtcCompressions+3][ix][iy]
		      << ", "  << fNtcAdc[kNtcCompressions+4][ix][iy]
		      << ", "  << fNtcAdc[kNtcCompressions+5][ix][iy]
		      << ", "  << fNtcAdc[kNtcCompressions+6][ix][iy]
		      << ", "  << fNtcAdc[kNtcCompressions+7][ix][iy]
		      << std::endl;
	    */

	    // Do bit dropping
	    for(unsigned b(0);b<8;b++) {
	      unsigned half(1<<b);
	      fNtcAdc[kNtcDropBits+b][ix][iy]=2*half*((fNtcAdc[kNtcReal][ix][iy]+half)/(2*half));
	    }
	  }
	}
      }
    }

    //std::cout << "Done normal trigger cells" << std::endl;

    for(unsigned ix(0);ix<fNumberOfNtcs;ix+=2) {
      for(unsigned iy(0);iy<fNumberOfNtcs;iy+=2) {

    	unsigned eMax(0);
	unsigned sum(0);
	unsigned nNtcTh(0);
	unsigned mx,my;

	for(unsigned jx(ix);jx<ix+2;jx++) {
	  for(unsigned jy(iy);jy<iy+2;jy++) {
	    if(fNtcAdc[kNtcReal][jx][jy]>0) nNtcTh++;

	    sum+=fNtcAdc[kNtcReal][jx][jy];
	    if(eMax<=fNtcAdc[kNtcReal][jx][jy]) {
	      eMax=fNtcAdc[kNtcReal][jx][jy];
	      mx=jx;
	      my=jy;
	    }
	  }
	}

	unsigned nx((mx%2)==0?mx+1:mx-1);
	unsigned ny((my%2)==0?my+1:my-1);

	//std::cout << "Doing normal trigger cells " << mx << ", " << my << std::endl;
	
	fPacker4E3M.pack((sum+4)/8);
	fPacker4E5M.pack((sum+1)/2);
	
	unsigned iNtcTh(0);
	for(unsigned jx(ix);jx<ix+2;jx++) {
	  for(unsigned jy(iy);jy<iy+2;jy++) {
	    if(jx==mx && jy==my) fNtcAdc[kNtcSuperTcs][jx][jy]=8*fPacker4E3M.unpack();
	    else fNtcAdc[kNtcSuperTcs][jx][jy]=0;

	    fNtcAdc[kNtcEqualShare][jx][jy]=(2*fPacker4E5M.unpack()+2*(jy-iy)+(jx-ix))/4;
	    
	    if(fNtcAdc[kNtcReal][jx][jy]>0) {
	      fNtcAdc[kNtcThresholdEqualShare][jx][jy]=(2*fPacker4E5M.unpack()+iNtcTh)/nNtcTh;
	      iNtcTh++;
	    } else fNtcAdc[kNtcThresholdEqualShare][jx][jy]=0;
	  }
	}
	assert(iNtcTh==nNtcTh);

      }
    }

    //std::cout << "Done coarse trigger cells" << std::endl;
    
    for(unsigned ix(0);ix<numberOfCoarseCells();ix++) {
      for(unsigned iy(0);iy<numberOfCoarseCells();iy++) {
	for(unsigned f(0);f<(kNtcCompressions-kNtcFractions);f++) setFractions(ix,iy,f,kNtcFractions+f);
	for(unsigned f(0);f<(kNtcDropBits-kNtcElinks);f++) setFractions(ix,iy,f+8,kNtcElinks+f);
      }
    }



    
    // Check
    unsigned nTcSum[kNumberOfNtcMethods];
    for(unsigned m(0);m<kNumberOfNtcMethods;m++) {
      //for(unsigned m(0);m<kNtcCompressions;m++) {
      nTcSum[m]=0;
      for(unsigned ix(0);ix<fNumberOfNtcs;ix++) {
	for(unsigned iy(0);iy<fNumberOfNtcs;iy++) {
	  nTcSum[m]+=fNtcAdc[m][ix][iy];
	  
	  hNtcAdcConversion0[m]->Fill(fNtcAdc[kNtcReal][ix][iy],fNtcAdc[m][ix][iy]);
	  hNtcAdcConversion1[m]->Fill(fNtcAdc[kNtcReal][ix][iy],fNtcAdc[m][ix][iy]);
	  hNtcAdcConversion2[m]->Fill(fNtcAdc[kNtcReal][ix][iy],fNtcAdc[m][ix][iy]);
	  hNtcAdcDeltaConversion0[m]->Fill(fNtcAdc[kNtcReal][ix][iy],int(fNtcAdc[m][ix][iy])-int(fNtcAdc[kNtcReal][ix][iy]));
	  hNtcAdcDeltaConversion1[m]->Fill(fNtcAdc[kNtcReal][ix][iy],int(fNtcAdc[m][ix][iy])-int(fNtcAdc[kNtcReal][ix][iy]));
	  hNtcAdcDeltaConversion2[m]->Fill(fNtcAdc[kNtcReal][ix][iy],int(fNtcAdc[m][ix][iy])-int(fNtcAdc[kNtcReal][ix][iy]));
	  if(fNtcAdc[kNtcReal][ix][iy]>0) {
	    hNtcAdcFracConversion0[m]->Fill(fNtcAdc[kNtcReal][ix][iy],(int(fNtcAdc[m][ix][iy])-int(fNtcAdc[kNtcReal][ix][iy]))/(1.0*fNtcAdc[kNtcReal][ix][iy]));
	    hNtcAdcFracConversion1[m]->Fill(fNtcAdc[kNtcReal][ix][iy],(int(fNtcAdc[m][ix][iy])-int(fNtcAdc[kNtcReal][ix][iy]))/(1.0*fNtcAdc[kNtcReal][ix][iy]));
	    hNtcAdcFracConversion2[m]->Fill(fNtcAdc[kNtcReal][ix][iy],(int(fNtcAdc[m][ix][iy])-int(fNtcAdc[kNtcReal][ix][iy]))/(1.0*fNtcAdc[kNtcReal][ix][iy]));
	  }
	}
      }
      //std::cout << "Ntc method " << m << " sum = " << nTcSum[m] << ", Sc sum = " << nScSum << std::endl;
      if(m!=kNtcEqualShare && m!=kNtcThreshold && m!=kNtcThresholdEqualShare && m!=kNtcSuperTcs && m<kNtcCompressions) assert(nScSum==nTcSum[m]);
    }

    for(unsigned m(0);m<kNumberOfNtcMethods;m++) {
      for(unsigned ix(0);ix<fNumberOfNtcs;ix+=2) {
	for(unsigned iy(0);iy<fNumberOfNtcs;iy+=2) {
	  double chiSq[2]={0.0,0.0};
	  for(unsigned jx(ix);jx<ix+2;jx++) {
	    for(unsigned jy(iy);jy<iy+2;jy++) {
	      chiSq[0]+=pow(int(fNtcAdc[m][jx][jy])-int(fNtcAdc[kNtcReal][jx][jy]),2);
	      if(fNtcAdc[kNtcReal][jx][jy]>0) chiSq[1]+=pow((int(fNtcAdc[m][jx][jy])-int(fNtcAdc[kNtcReal][jx][jy]))/(1.0*fNtcAdc[kNtcReal][jx][jy]),2);
	    }
	  }
	  hNtcAdcChiSq0[m]->Fill(chiSq[0]);
	  hNtcAdcChiSq1[m]->Fill(chiSq[1]);

	  numChiSq0[m]++;
	  sumChiSq0[m]+=chiSq[0];
	  numChiSq1[m]++;
	  sumChiSq1[m]+=chiSq[1];
	}
      }

      hNtcAdcChiSqAvg0->SetBinContent(m+1,sumChiSq0[m]/numChiSq0[m]);
      hNtcAdcChiSqAvg1->SetBinContent(m+1,sumChiSq1[m]/numChiSq1[m]);
    }


#ifdef CRAP


    
    unsigned nSum[5]={0,0,0,0,0};
    unsigned nSum2[kNumberOfNtcMethods]={0,0,0,0,0,0,0};
    unsigned nSum3[kNumberOfCtcMethods]={0,0,0,0};



    for(unsigned ix(0);ix<numberOfTriggerCells();ix++) {
      for(unsigned iy(0);iy<numberOfTriggerCells();iy++) {
	nSum[0]+=fTcAdc[ix][iy];
	
	//fStcAdc[ix][iy]=fTcAdc[ix][iy];
	//if(fTcAdc[ix][iy]<mipT2) fTcAdc[ix][iy]=0; // Do later

	if(fNtcAdc[0][ix][iy]>=mipT2) fNtcAdc[1][ix][iy]=fNtcAdc[0][ix][iy];
      }
    }
    
    for(unsigned ix(0);ix<numberOfTriggerCells();ix+=2) {
      for(unsigned iy(0);iy<numberOfTriggerCells();iy+=2) {
	
	unsigned eMax(0);
	unsigned sum(0);
	unsigned mx,my;
	for(unsigned jx(ix);jx<ix+2;jx++) {
	  for(unsigned jy(iy);jy<iy+2;jy++) {
	    sum+=fTcAdc[jx][jy];
	    if(eMax<=fTcAdc[jx][jy]) {
	      eMax=fTcAdc[jx][jy];
	      mx=jx;
	      my=jy;
	    }
	    //fStcAdc[jx][jy]=0;
	  }
	}
	fStcAdc[mx][my]=sum;
	fNtcAdc[4][mx][my]=sum;

	//std::cout << "TCs = " << fTcAdc[ix][iy] << ", " << fTcAdc[ix+1][iy] << ", " << fTcAdc[ix][iy+1] << ", " << fTcAdc[ix+1][iy+1] << ", sum = " << sum
	//	  << ", STCs = " << fStcAdc[ix][iy] << ", " << fStcAdc[ix+1][iy] << ", " << fStcAdc[ix][iy+1] << ", " << fStcAdc[ix+1][iy+1] << std::endl;
	assert((fTcAdc[ix][iy]+fTcAdc[ix+1][iy]+fTcAdc[ix][iy+1]+fTcAdc[ix+1][iy+1])==sum);

	for(unsigned jx(ix);jx<ix+2;jx++) {
	  for(unsigned jy(iy);jy<iy+2;jy++) {
	    fNtcAdc[2][jx][jy]=(sum+2*(jy-iy)+(jx-ix))/4;
	  }
	}

	unsigned nx((mx%2)==0?mx+1:mx-1);
	unsigned ny((my%2)==0?my+1:my-1);

	unsigned e32(32*fTcAdc[mx][my]);
	unsigned e64(64*fTcAdc[mx][my]);
	
	unsigned ratio(0);
	for(unsigned r(1);r<8;r++) {
	  if(e32<((23-2*r)*sum)) ratio=r;
	  //if(e64<((61-6*r)*sum)) ratio=r;
	}

	/for(unsigned r(1);r<16;r++) {
	//if(e64<((23-2*r)*sum)) ratio=r;
	//}

	//std::cout << "ETot = " << sum << ", E1 = " << fTcAdc[mx][my] << ", ratio = " << ratio
	//	  << " reco E1 = " << ((23-2*ratio)*sum)/32 << std::endl;	
	//	  << " reco E1 = " << ((61-6*ratio)*sum)/64 << std::endl;	

	unsigned sum4(sum/4);
	unsigned sum8(sum/8);
	//unsigned sum9(sum/9);
	unsigned sum16(sum/16);

	//if(sum==fTcAdc[mx][my]) {
	//fFtcAdc[mx][my]=sum;
	//fFtcAdc[nx][my]=0;
	//fFtcAdc[mx][ny]=0;
	//fFtcAdc[nx][ny]=0;
	  
	//fNtcAdc[3][mx][my]=sum;
	//} else {
	/*
	  fNtcAdc[3][mx][my]=sum-2*sum4-sum8;
	  fNtcAdc[3][nx][my]=sum4;
	  fNtcAdc[3][mx][ny]=sum4;
	  fNtcAdc[3][nx][ny]=sum8;	  
	*/
	/*	
	  fNtcAdc[3][mx][my]=sum-3*sum4;
	  fNtcAdc[3][nx][my]=sum4;
	  fNtcAdc[3][mx][ny]=sum4;
	  fNtcAdc[3][nx][ny]=sum4;	  
	*/

	  // f = 1/3
	  //fNtcAdc[3][nx][my]=3*sum16;
	  //fNtcAdc[3][mx][ny]=3*sum16;
	  //fNtcAdc[3][nx][ny]=sum16;

	/*
	// BY HAND
	fNtcAdc[3][nx][my]=sum4;
	fNtcAdc[3][mx][ny]=sum4;
	if(fTcAdc[nx][ny]<sum8) fNtcAdc[3][nx][ny]=  sum16;
	else                    fNtcAdc[3][nx][ny]=3*sum16;
	*/

	// Three bits
	// E2,E3 = 1/8 -> 1/4, E4 = 0 -> 1/4
	fNtcAdc[3][nx][my]=((17+2*ratio)*sum)/128;
	fNtcAdc[3][mx][ny]=((17+2*ratio)*sum)/128;
	fNtcAdc[3][nx][ny]=(( 2+4*ratio)*sum)/128;

	// E2,E3,E4 = 0 -> 1/4
	//fNtcAdc[3][nx][my]=((1+2*ratio)*sum)/64;
	//fNtcAdc[3][mx][ny]=((1+2*ratio)*sum)/64;
	//fNtcAdc[3][nx][ny]=((1+2*ratio)*sum)/64;

	// Check
	//fNtcAdc[3][nx][my]=fTcAdc[nx][my];
	//fNtcAdc[3][mx][ny]=fTcAdc[mx][ny];
	//fNtcAdc[3][nx][ny]=fTcAdc[nx][ny];
	
	fNtcAdc[3][mx][my]=sum-fNtcAdc[3][nx][my]-fNtcAdc[3][mx][ny]-fNtcAdc[3][nx][ny];
	//assert(fNtcAdc[3][mx][my]==fNtcAdc[0][mx][my]);
	
	//std::cout << "E1 = " << fTcAdc[mx][my] << " ~ " << fNtcAdc[3][mx][my]
	//	  << ", E2 = " << fTcAdc[mx][ny] << " ~ " << fNtcAdc[3][mx][ny]
	//	  << ", E3 = " << fTcAdc[nx][my] << " ~ " << fNtcAdc[3][nx][my]
	//	  << ", E4 = " << fTcAdc[nx][ny] << " ~ " << fNtcAdc[3][nx][ny] << std::endl;
	
	  // f = 1/2
	  //fFtcAdc[mx][my]=sum-5*sum9;
	  //fFtcAdc[nx][my]=2*sum9;
	  //fFtcAdc[mx][ny]=2*sum9;
	  //fFtcAdc[nx][ny]=sum9;
	  //}

	  fFtcAdc[mx][my]=fNtcAdc[3][mx][my];
	  fFtcAdc[nx][my]=fNtcAdc[3][nx][my];
	  fFtcAdc[mx][ny]=fNtcAdc[3][mx][ny];
	  fFtcAdc[nx][ny]=fNtcAdc[3][nx][ny];

	  
	fCtcAdc[0][ix/2][iy/2]=sum;
      }
    }


    for(unsigned ix(0);ix<numberOfCoarseCells();ix++) {
      for(unsigned iy(0);iy<numberOfCoarseCells();iy++) {
	setFractions(ix,iy,7,3);
      }
    }

    

    for(unsigned ix(0);ix<numberOfCoarseCells();ix+=2) {
      for(unsigned iy(0);iy<numberOfCoarseCells();iy+=2) {
	
	unsigned eMax(0);
	unsigned sum(0);
	unsigned mx,my;
	for(unsigned jx(ix);jx<ix+2;jx++) {
	  for(unsigned jy(iy);jy<iy+2;jy++) {
	    sum+=fCtcAdc[0][jx][jy];
	    if(eMax<=fCtcAdc[0][jx][jy]) {
	      eMax=fCtcAdc[0][jx][jy];
	      mx=jx;
	      my=jy;
	    }
	  }
	}
	fCtcAdc[3][mx][my]=sum;

	for(unsigned jx(ix);jx<ix+2;jx++) {
	  for(unsigned jy(iy);jy<iy+2;jy++) {
	    fCtcAdc[1][jx][jy]=(sum+2*(jy-iy)+(jx-ix))/4;
	  }
	}

	unsigned nx((mx%2)==0?mx+1:mx-1);
	unsigned ny((my%2)==0?my+1:my-1);

	unsigned sum4(sum/4);
	unsigned sum8(sum/8);
	//unsigned sum9(sum/9);
	unsigned sum16(sum/16);

	fCtcAdc[2][nx][my]=sum4;
	fCtcAdc[2][mx][ny]=sum4;
	if(fCtcAdc[0][nx][ny]<sum8) fCtcAdc[2][nx][ny]=  sum16;
	else                        fCtcAdc[2][nx][ny]=3*sum16;
	
	fCtcAdc[2][mx][my]=sum-fCtcAdc[2][nx][my]-fCtcAdc[2][mx][ny]-fCtcAdc[2][nx][ny];
      }
    }

    for(unsigned ix(0);ix<numberOfCoarseCells();ix++) {
      for(unsigned iy(0);iy<numberOfCoarseCells();iy++) {
	for(unsigned m(0);m<kNumberOfCtcMethods;m++) nSum3[m]+=fCtcAdc[m][ix][iy];
      }
    }

    if(fHiDensity) {
      for(unsigned ix(0);ix<numberOfTriggerCells();ix+=2) {
	for(unsigned iy(0);iy<numberOfTriggerCells();iy+=4) {
	  unsigned sum(0);
	  unsigned eMax(0);
	  unsigned mx,my;
	  for(unsigned jx(ix);jx<ix+2;jx++) {
	    for(unsigned jy(iy);jy<iy+4;jy++) {
	      sum+=fTcAdc[jx][jy];
	      if(eMax<=fTcAdc[jx][jy]) {
		eMax=fTcAdc[jx][jy];
		mx=jx;
		my=jy;
	      }
	    }
	  }
	  fNtcAdc[6][mx][my]=sum;
	  
	  for(unsigned jx(ix);jx<ix+2;jx++) {
	    for(unsigned jy(iy);jy<iy+4;jy++) {
	      fNtcAdc[5][jx][jy]=(sum+2*(jy-iy)+(jx-ix))/8;
	    }
	  }
	}
      }
      
    } else {
      for(unsigned ix(0);ix<numberOfTriggerCells();ix+=4) {
	for(unsigned iy(0);iy<numberOfTriggerCells();iy+=4) {
	  unsigned sum(0);
	  unsigned eMax(0);
	  unsigned mx,my;
	  for(unsigned jx(ix);jx<ix+4;jx++) {
	    for(unsigned jy(iy);jy<iy+4;jy++) {
	      sum+=fTcAdc[jx][jy];
	      if(eMax<=fTcAdc[jx][jy]) {
		eMax=fTcAdc[jx][jy];
		mx=jx;
		my=jy;
	      }
	    }
	  }
	  fNtcAdc[6][mx][my]=sum;

	  for(unsigned jx(ix);jx<ix+4;jx++) {
	    for(unsigned jy(iy);jy<iy+4;jy++) {
	      fNtcAdc[5][jx][jy]=(sum+4*(jy-iy)+(jx-ix))/16;
	    }
	  }
	}
      }
    }

    if(fHiDensity) {
      for(unsigned ix(0);ix<numberOfTriggerCells();ix+=2) {
	for(unsigned iy(0);iy<numberOfTriggerCells();iy+=4) {
	  fHtcAdc[ix/2][iy/4]=0;
	}
      }
      for(unsigned ix(0);ix<numberOfTriggerCells();ix+=2) {
	for(unsigned iy(0);iy<numberOfTriggerCells();iy+=2) {
	  fHtcAdc[ix/2][iy/4]+=fCtcAdc[0][ix/2][iy/2];
	}
      }
    } else {
      for(unsigned ix(0);ix<numberOfTriggerCells();ix+=4) {
	for(unsigned iy(0);iy<numberOfTriggerCells();iy+=4) {
	  fHtcAdc[ix/4][iy/4]=0;
	}
      }
      for(unsigned ix(0);ix<numberOfTriggerCells();ix+=2) {
	for(unsigned iy(0);iy<numberOfTriggerCells();iy+=2) {
	  fHtcAdc[ix/4][iy/4]+=fCtcAdc[0][ix/2][iy/2];
	}
      }
    }

    // Cross check

    for(unsigned ix(0);ix<numberOfTriggerCells();ix++) {
      for(unsigned iy(0);iy<numberOfTriggerCells();iy++) {
	nSum[1]+=fStcAdc[ix][iy];
	nSum[2]+=fNtcAdc[3][ix][iy];

	assert(fNtcAdc[0][ix][iy]==fTcAdc[ix][iy]);
	for(unsigned m(0);m<kNumberOfNtcMethods;m++) nSum2[m]+=fNtcAdc[m][ix][iy];
      }
    }
    for(unsigned ix(0);ix<numberOfCoarseCells();ix++) {
      for(unsigned iy(0);iy<numberOfCoarseCells();iy++) {
	nSum[3]+=fCtcAdc[0][ix][iy];
      }
    }
    for(unsigned ix(0);ix<numberOfHgcrocXCells();ix++) {
      for(unsigned iy(0);iy<numberOfHgcrocYCells();iy++) {
	nSum[4]+=fHtcAdc[ix][iy];
      }
    }

    //std::cout << "Check sums = " << nSum[0] << ", " << nSum[1] << ", " << nSum[2] << ", " << nSum[3] << ", " << nSum[4] << std::endl;
    //std::cout << "Check sums2 = " << nSum2[0] << ", " << nSum2[1] << ", " << nSum2[2] << ", " << nSum2[3] << ", " << nSum2[4]  << ", " << nSum2[5] << ", " << nSum2[6] << std::endl;
    assert(nSum[0]==nSum[1]);
    assert(nSum[0]==nSum[2]);
    assert(nSum[0]==nSum[3]);
    assert(nSum[0]==nSum[4]);

    assert(nSum[0]==nSum2[0]);
    for(unsigned m(2);m<kNumberOfNtcMethods;m++) assert(nSum[0]==nSum2[m]);

    for(unsigned m(0);m<kNumberOfCtcMethods;m++) assert(nSum[0]==nSum3[m]);
#endif

    }
  
  double scNoiseEnergy(unsigned ix, unsigned iy) const {
    assert(ix<numberOfSensorCells());
    assert(iy<numberOfSensorCells());
    return fNoiseArray[ix][iy];
  }
  
  double scPhysicalEnergy(unsigned ix, unsigned iy) const {
    assert(ix<numberOfSensorCells());
    assert(iy<numberOfSensorCells());
    return fPileupArray[ix][iy]+fEnergyArray[ix][iy];
  }
  
  double scTotalEnergy(unsigned ix, unsigned iy) const {
    assert(ix<numberOfSensorCells());
    assert(iy<numberOfSensorCells());
    return fNoiseArray[ix][iy]+fPileupArray[ix][iy]+fEnergyArray[ix][iy];
  }

  double scTrueEnergy(unsigned ix, unsigned iy) const {
    assert(ix<numberOfSensorCells());
    assert(iy<numberOfSensorCells());
    return fEnergyArray[ix][iy];
  }

  double scAdcEnergy(unsigned ix, unsigned iy, unsigned m) const {
    assert(ix<numberOfSensorCells());
    assert(iy<numberOfSensorCells());
    assert(m<kNumberOfScMethods);

    if(m<(kNumberOfScMethods-1)) return fScAdc[m][ix][iy]/fMip;
    return scTotalEnergy(ix,iy);
  }
  
  unsigned scAdc(unsigned ix, unsigned iy, unsigned m) const {
    assert(ix<numberOfSensorCells());
    assert(iy<numberOfSensorCells());
    assert(m<kNumberOfScMethods);
    return fScAdc[m][ix][iy];
  }
  
  double tcPhysicalEnergy(unsigned ix, unsigned iy) const {
    assert(ix<numberOfSensorCells());
    assert(iy<numberOfSensorCells());

    unsigned nSc(fHiDensity?3:2);
    double e(0.0);
    for(unsigned iScX(nSc*ix);iScX<nSc*(ix+1);iScX++) {
      for(unsigned iScY(nSc*iy);iScY<nSc*(iy+1);iScY++) {
	e+=scPhysicalEnergy(iScX,iScY);
      }
    }
    return e;
  }

  double tcTotalEnergy(unsigned ix, unsigned iy) const {
    assert(ix<numberOfTriggerCells());
    assert(iy<numberOfTriggerCells());

    unsigned nSc(3);
    double e(0.0);
    for(unsigned iScX(nSc*ix);iScX<nSc*(ix+1);iScX++) {
      for(unsigned iScY(nSc*iy);iScY<nSc*(iy+1);iScY++) {
	e+=scTotalEnergy(iScX,iScY);
      }
    }
    return e;
  }

  double ctcTotalEnergy(unsigned ix, unsigned iy) const {
    assert(ix<fNumberOfCtcs);
    assert(iy<fNumberOfCtcs);

    unsigned nSc(fHiDensity?6:4);
    double e(0.0);
    for(unsigned iScX(nSc*ix);iScX<nSc*(ix+1);iScX++) {
      for(unsigned iScY(nSc*iy);iScY<nSc*(iy+1);iScY++) {
	e+=scTotalEnergy(iScX,iScY);
      }
    }
    return e;
  }
  /*
  double htcTotalEnergy(unsigned ix, unsigned iy) const {
    assert(ix<numberOfHgcrocXCells());
    assert(iy<numberOfHgcrocYCells());

    unsigned nScX(fHiDensity?6:8);
    unsigned nScY(fHiDensity?12:8);
    double e(0.0);
    for(unsigned iScX(nScX*ix);iScX<nScX*(ix+1);iScX++) {
      for(unsigned iScY(nScY*iy);iScY<nScY*(iy+1);iScY++) {
	e+=scTotalEnergy(iScX,iScY);
      }
    }
    return e;
  }
  */
  double tcAdcEnergy(unsigned ix, unsigned iy, unsigned m) const {
    assert(ix<numberOfTriggerCells());
    assert(iy<numberOfTriggerCells());
    assert(m<Array::kNumberOfNtcMethods);

    if(m<(kNumberOfNtcMethods-1)) return fNtcAdc[m][ix][iy]/fMip;
    return tcTotalEnergy(ix,iy);
  }

  double ntcAdcEnergy(unsigned ix, unsigned iy, unsigned m) const {
    assert(ix<numberOfTriggerCells());
    assert(iy<numberOfTriggerCells());
    assert(m<Array::kNumberOfNtcMethods);

    if(m<(kNumberOfNtcMethods-1)) return fNtcAdc[m][ix][iy]/fMip;
    return tcTotalEnergy(ix,iy);
  }
  /*  
  double stcAdcEnergy(unsigned ix, unsigned iy) const {
    assert(ix<numberOfTriggerCells());
    assert(iy<numberOfTriggerCells());
    return fStcAdc[ix][iy]/fMip;
  }

  double ftcAdcEnergy(unsigned ix, unsigned iy) const {
    assert(ix<numberOfTriggerCells());
    assert(iy<numberOfTriggerCells());
    return fFtcAdc[ix][iy]/fMip;
  }

  double ctcAdcEnergy(unsigned ix, unsigned iy, unsigned m) const {
    assert(ix<numberOfCoarseCells());
    assert(iy<numberOfCoarseCells());
    assert(m<Array::kNumberOfCtcMethods);

    //if(m==1) {
    //  if(fHiDensity) return fHtcAdc[ix][iy/2]/(2.0*fMip);
    //  return fHtcAdc[ix/2][iy/2]/(4.0*fMip);
    //}

    return fCtcAdc[m][ix][iy]/fMip;
  }

  double htcAdcEnergy(unsigned ix, unsigned iy) const {
    assert(ix<numberOfHgcrocXCells());
    assert(iy<numberOfHgcrocYCells());
    return fHtcAdc[ix][iy]/fMip;
  }
*/

  const bool fHiDensity;
  const bool fEnergy;
  const bool fPileup;
  const double fEta;

  const unsigned fNumberOfScs;
  const unsigned fNumberOfNtcs;
  const unsigned fNumberOfCtcs;
  
  const double fTriggerCellSide;
  const double fCoarseCellSide;
  const double fSensorCellSide;
  const double fArraySide;
  const double fArrayHalfSide;
  const double fNoise;
  const double fMip;
  const double fSinTheta;

  TH2F *hDeltaE1[kNumberOfNtcMethods];
  TH2F *hDeltaE2[kNumberOfNtcMethods];
  TH2F *hScAdcConversion;
  TH2F *hScAdcDeltaConversion;
  TH2F *hNtcAdcConversion0[kNumberOfNtcMethods];
  TH2F *hNtcAdcConversion1[kNumberOfNtcMethods];
  TH2F *hNtcAdcConversion2[kNumberOfNtcMethods];
  TH2F *hNtcAdcDeltaConversion0[kNumberOfNtcMethods];
  TH2F *hNtcAdcDeltaConversion1[kNumberOfNtcMethods];
  TH2F *hNtcAdcDeltaConversion2[kNumberOfNtcMethods];
  TH2F *hNtcAdcFracConversion0[kNumberOfNtcMethods];
  TH2F *hNtcAdcFracConversion1[kNumberOfNtcMethods];
  TH2F *hNtcAdcFracConversion2[kNumberOfNtcMethods];
  TH1F *hNtcAdcChiSq0[kNumberOfNtcMethods];
  TH1F *hNtcAdcChiSq1[kNumberOfNtcMethods];
  TH1F *hNtcAdcChiSqAvg0;
  TH1F *hNtcAdcChiSqAvg1;
  
  double numChiSq0[kNumberOfNtcMethods];
  double sumChiSq0[kNumberOfNtcMethods];
  double numChiSq1[kNumberOfNtcMethods];
  double sumChiSq1[kNumberOfNtcMethods];
  
  PackerFloat22 fPacker5E1M;
  PackerFloat22 fPacker3E3M;
  PackerFloat22 fPacker4E6M;
  PackerFloat22 fPacker4E5M;
  PackerFloat22 fPacker4E4M;
  PackerFloat22 fPacker4E3M;
  PackerFloat22 fPacker4E2M;
  PackerFloat22 fPacker4E1M;

  double fNoiseArray[kHiDensityScs][kHiDensityScs];
  double fPileupArray[kHiDensityScs][kHiDensityScs];
  double fEnergyArray[kHiDensityScs][kHiDensityScs];

  unsigned fScAdc[kNumberOfScMethods][kHiDensityScs][kHiDensityScs];
		  		  
  //unsigned fTcAdc[kTriggerCells][kTriggerCells];
  unsigned fNtcAdc[kNumberOfNtcMethods][kHiDensityNtcs][kHiDensityNtcs];
  //unsigned fStcAdc[kTriggerCells][kTriggerCells];
  //unsigned fFtcAdc[kTriggerCells][kTriggerCells];
  //unsigned fCtcAdc[kNumberOfCtcMethods][kCoarseCells][kCoarseCells];
  //unsigned fHtcAdc[kTriggerCells/2][kTriggerCells/4];
};

void histEvent(unsigned n, const Array &a) {
  std::ostringstream sout;
  sout << "Evt" << n;
  
  TH2F *h2[2][2];
  h2[0][0]=new TH2F((sout.str()+"ScEnergy").c_str(),
		    ";Sensor cell x;Sensor cell y",
		    a.numberOfSensorCells(),0,a.numberOfSensorCells(),
		    a.numberOfSensorCells(),0,a.numberOfSensorCells());
  h2[0][1]=new TH2F((sout.str()+"ScAdcEnergy").c_str(),
		    ";Sensor cell x;Sensor cell y",
		    a.numberOfSensorCells(),0,a.numberOfSensorCells(),
		    a.numberOfSensorCells(),0,a.numberOfSensorCells());
  h2[1][1]=new TH2F((sout.str()+"TcAdcEnergy").c_str(),
		    ";Trigger cell x;Trigger cell y",
		    a.numberOfTriggerCells(),0,a.numberOfTriggerCells(),
		    a.numberOfTriggerCells(),0,a.numberOfTriggerCells());

  TH2F *hTrueEnergyTcCm=new TH2F((sout.str()+"TrueEnergyTcCm").c_str(),
				 ";x (cm);y (cm)",
				 a.numberOfTriggerCells(),-0.5*a.fArraySide,0.5*a.fArraySide,
				 a.numberOfTriggerCells(),-0.5*a.fArraySide,0.5*a.fArraySide);
  
  TH2F *hTcAdcEnergy[Array::kNumberOfNtcMethods],*hTcAdcEnergyCm[Array::kNumberOfNtcMethods];
  for(unsigned m(0);m<Array::kNumberOfNtcMethods;m++) {
    std::ostringstream sout2;
    sout2 << "Method" << std::setw(2) << std::setfill('0') << m;
    
    hTcAdcEnergy[m]=new TH2F((sout.str()+sout2.str()+"TcAdcEnergy").c_str(),
			     ";Trigger cell x;Trigger cell y",
			     a.numberOfTriggerCells(),0,a.numberOfTriggerCells(),
			     a.numberOfTriggerCells(),0,a.numberOfTriggerCells());
    hTcAdcEnergyCm[m]=new TH2F((sout.str()+sout2.str()+"TcAdcEnergyCm").c_str(),
			       ";x (cm);y (cm)",
			       a.numberOfTriggerCells(),-0.5*a.fArraySide,0.5*a.fArraySide,
			       a.numberOfTriggerCells(),-0.5*a.fArraySide,0.5*a.fArraySide);
  }

  /*
  TH2F *hCtcAdcEnergy[Array::kNumberOfCtcMethods];
  for(unsigned m(0);m<Array::kNumberOfCtcMethods;m++) {
    std::ostringstream sout2;
    sout2 << "Method" << m;
    
    hCtcAdcEnergy[m]=new TH2F((sout.str()+sout2.str()+"CtcAdcEnergy").c_str(),
			      ";Coarse cell x;Coarse cell y",
			      a.numberOfCoarseCells(),0,a.numberOfCoarseCells(),
			      a.numberOfCoarseCells(),0,a.numberOfCoarseCells());
  }
  
  TH2F *hHtcAdcEnergy[1];
  for(unsigned m(0);m<1;m++) {
    std::ostringstream sout2;
    sout2 << "Method" << m;
    
    hHtcAdcEnergy[m]=new TH2F((sout.str()+sout2.str()+"HtcAdcEnergy").c_str(),
			      ";Hgcroc cell x;Hgcroc cell y",
			      a.numberOfHgcrocXCells(),0,a.numberOfHgcrocXCells(),
			      a.numberOfHgcrocYCells(),0,a.numberOfHgcrocYCells());
  }
  */
  /*
    TH2F *hFtcAdcEnergy=new TH2F((sout.str()+"FtcAdcEnergy").c_str(),
    ";Trigger cell x;Trigger cell y",
    a.numberOfTriggerCells(),0,a.numberOfTriggerCells(),
    a.numberOfTriggerCells(),0,a.numberOfTriggerCells());
    TH2F *hCtcAdcEnergy=new TH2F((sout.str()+"CtcAdcEnergy").c_str(),
    ";Trigger cell x;Trigger cell y",
    a.numberOfTriggerCells()/2,0,a.numberOfTriggerCells()/2,
    a.numberOfTriggerCells()/2,0,a.numberOfTriggerCells()/2);
  */

  double tcArea(a.fHiDensity?9.0:4.0);
  
  for(unsigned ix(0);ix<a.numberOfSensorCells();ix++) {
    for(unsigned iy(0);iy<a.numberOfSensorCells();iy++) {
      h2[0][0]->Fill(ix,iy,a.scTotalEnergy(ix,iy));
      hTrueEnergyTcCm->Fill(a.sensorCellPosition(ix),a.sensorCellPosition(iy),a.scTrueEnergy(ix,iy)/tcArea);

      h2[0][1]->Fill(ix,iy,a.scAdcEnergy(ix,iy,0));
    }
  }
  
  for(unsigned ix(0);ix<a.numberOfTriggerCells();ix++) {
    for(unsigned iy(0);iy<a.numberOfTriggerCells();iy++) {
      for(unsigned m(0);m<Array::kNumberOfNtcMethods;m++) {
	hTcAdcEnergy[m]->Fill(ix,iy,a.tcAdcEnergy(ix,iy,m)/tcArea);
	hTcAdcEnergyCm[m]->Fill(a.triggerCellPosition(ix),a.triggerCellPosition(iy),a.tcAdcEnergy(ix,iy,m)/tcArea);
      }
    }
  }

  tcArea*=4.0;
  /*  
  for(unsigned ix(0);ix<a.numberOfCoarseCells();ix++) {
    for(unsigned iy(0);iy<a.numberOfCoarseCells();iy++) {
      for(unsigned m(0);m<Array::kNumberOfCtcMethods;m++) {
	hCtcAdcEnergy[m]->Fill(ix,iy,a.ctcAdcEnergy(ix,iy,m)/tcArea);
      }
    }
  }

  tcArea*=(a.fHiDensity?2.0:4.0);
  
  for(unsigned ix(0);ix<a.numberOfHgcrocXCells();ix++) {
    for(unsigned iy(0);iy<a.numberOfHgcrocYCells();iy++) {
      for(unsigned m(0);m<1;m++) {
	hHtcAdcEnergy[m]->Fill(ix,iy,a.htcAdcEnergy(ix,iy)/tcArea);
      }
    }
  }
  */
}

int main(int argc, const char **argv) {
  unsigned place(0);
  bool ez(true);
  bool pu(true);
  double fLayerZ(330.0);

  double fPt=40.0;
  //fPt=1000.0;
  double fPtSpread=0.5;
  //double fEta=3.0;
  double fMipsPerGeV(10.0);
  
  if(argc>4) {
    std::cerr << "Usage: " << argv[0] << " h or l" << std::endl;
    return 1;
  }

  if(argc==2) {
    std::istringstream sin(argv[1]);
    sin >> place;
  }

  if(argc==3) {
    if(argv[2][0]=='e') ez=true;
    if(argv[2][0]=='z') ez=false;
  }

  if(argc==4) {
    if(argv[3][0]=='p') pu=true;
    if(argv[3][0]=='u') pu=false;
  }

  TFileHandler tfh(std::string("SuperTCs")+(place==0?"0":"1"));

  PackerFloat22 aPacker(4,7);

  unsigned bins(1<<22);
  TH1F *hPackerA=new TH1F("ScPackerA",";Input;Output",bins,0,bins);
  TH1F *hPackerD=new TH1F("ScPackerD",";Input;Output-Input",bins,0,bins);
  TH1F *hPackerR=new TH1F("ScPackerR",";Input;(Output-Input)/Input",bins,0,bins);
  
  for(unsigned i(0);i<bins;i++) {
    aPacker.pack(i);
    unsigned op(aPacker.unpack());
    hPackerA->SetBinContent(i+1,op);
    hPackerD->SetBinContent(i+1,int(op)-int(i));
    if(i>0) hPackerR->SetBinContent(i+1,(int(op)-int(i))/(1.0*i));
  }

  PackerFloat22 bPacker(4,3);

  bins=bPacker.dataMax()+1;
  TH1F *hBPackerA=new TH1F("ScBPackerA",";Input;Output",bins,-0.5,bins-0.5);
  TH1F *hBPackerD=new TH1F("ScBPackerD",";Input;Output-Input",bins,-0.5,bins-0.5);
  TH1F *hBPackerR=new TH1F("ScBPackerR",";Input;(Output-Input)/Input",bins,-0.5,bins-0.5);
  TH1F *hCPackerA=new TH1F("ScCPackerA",";Input;Output",bins,-0.5,bins-0.5);
  TH1F *hCPackerD=new TH1F("ScCPackerD",";Input;Output-Input",bins,-0.5,bins-0.5);
  TH1F *hCPackerR=new TH1F("ScCPackerR",";Input;(Output-Input)/Input",bins,-0.5,bins-0.5);
  
  for(unsigned i(0);i<bins;i++) {
    bPacker.pack(i);
    unsigned op(bPacker.unpack());
    hBPackerA->SetBinContent(i+1,op);
    hBPackerD->SetBinContent(i+1,int(op)-int(i));
    if(i>0) hBPackerR->SetBinContent(i+1,(int(op)-int(i))/(1.0*i));

    op=bPacker.unpackTimes2();
    hCPackerA->SetBinContent(i+1,op/2.0);
    hCPackerD->SetBinContent(i+1,int(op)/2.0-int(i));
    if(i>0) hCPackerR->SetBinContent(i+1,(int(op)/2.0-int(i))/(1.0*i));
  }
  
  //const double dTriggerCells(sqrt(4.908));
  //const unsigned nSensorCells[2]={36,24};
  //const double dSensorCells[2]={dTriggerCells/3.0,dTriggerCells/2.0};

  //const unsigned nDensity(0);

  Array theArray(place,ez,pu);

  PhotonModel thePhotonModel[3];
  thePhotonModel[0]=PhotonModel(0,0.61);
  //thePhotonModel[1]=PhotonModel(1,3.0);
  thePhotonModel[1]=PhotonModel(1,0.059);
  thePhotonModel[2]=PhotonModel(2,1.05);
  unsigned usePhotonModel(1);
  
  double nrg(200.0);
  
  TH2F *hCentre=new TH2F("Centre",";x (cm);y(cm)",
			 10*theArray.numberOfTriggerCells(),
			 -theArray.fArrayHalfSide,theArray.fArrayHalfSide,
			 10*theArray.numberOfTriggerCells(),
			 -theArray.fArrayHalfSide,theArray.fArrayHalfSide);
  
  unsigned nScEventTotal(theArray.numberOfSensorCells()*theArray.numberOfSensorCells());
  unsigned nTcEventTotal(theArray.numberOfTriggerCells()*theArray.numberOfTriggerCells());

  TH2F *hNtcCountVsR[2];
  hNtcCountVsR[0]=new TH2F("NtcCountVsR0",
			   ";Cut radius #rho;Number of trigger cells",
			   kNumberOfRadii,0.0,0.06,300,0,300);
  hNtcCountVsR[1]=new TH2F("NtcCountVsR1",
			   ";Cut radius #rho;Number of trigger cells",
			   kNumberOfRadii,0.0,0.06,300,0,300);

  for(unsigned r(0);r<kNumberOfRadii;r++) {
    unsigned n[2];
    double dx,dy,dr;

    n[0]=0;
    for(unsigned ix(0);ix<theArray.numberOfTriggerCells();ix++) {
      for(unsigned iy(0);iy<theArray.numberOfTriggerCells();iy++) {
	dx=theArray.triggerCellPosition(ix)-theArray.triggerCellPosition(theArray.numberOfTriggerCells()/2);
	dy=theArray.triggerCellPosition(iy)-theArray.triggerCellPosition(theArray.numberOfTriggerCells()/2);
	dr=sqrt(dx*dx+dy*dy);

	if(dr<fLayerZ*cutRadius(r)) n[0]++;
      }
    }
    hNtcCountVsR[0]->Fill(cutRadius(r),n[0]);

    for(unsigned jx(0);jx<100;jx++) {
      for(unsigned jy(0);jy<100;jy++) {
	
	n[1]=0;
	for(unsigned ix(0);ix<theArray.numberOfTriggerCells();ix++) {
	  for(unsigned iy(0);iy<theArray.numberOfTriggerCells();iy++) {
	    dx=theArray.triggerCellPosition(ix)-theArray.fCoarseCellSide*(jx-49.5)/100.0;
	    dy=theArray.triggerCellPosition(iy)-theArray.fCoarseCellSide*(jy-49.5)/100.0;
	    dr=sqrt(dx*dx+dy*dy);
	    
	    if(dr<fLayerZ*cutRadius(r)) n[1]++;
	  }
	}
	hNtcCountVsR[1]->Fill(cutRadius(r),n[1],0.0001);
      }
    }
  }
  
  
  TH1F *hScReadout=new TH1F("ScReadout",";Number of sensor cells",nScEventTotal,0,nScEventTotal);

  TH1F *hScEnergyTrue=new TH1F("ScEnergyTrue",";Energy (MIPs)",100,0.0,50*nrg);
  TE1F *hScEnergyIntegral=new TE1F("ScEnergyIntegral",";Radius (cm);Energy (MIPs)",100,0.0,5.0,false);
  TH1F *hScEnergySpectrum=new TH1F("ScEnergySpectrum",";Energy (MIPs)",120,-2.0,10.0);

  //  unsigned nr(10);
  TH2F *hScResponse[Array::kNumberOfScMethods][kNumberOfRadii];
  TH1A *hAScResponse[Array::kNumberOfScMethods][kNumberOfRadii];
  TH2F *hNtcResponse[2][Array::kNumberOfNtcMethods][kNumberOfRadii];
  TH1A *hANtcResponse[2][Array::kNumberOfNtcMethods][kNumberOfRadii];

  /*
  TH2F *hStcResponse[kNumberOfRadii],*hFtcResponse[kNumberOfRadii];
  TH2F *hTrueCtcResponse[kNumberOfRadii],*hCtcResponse[Array::kNumberOfCtcMethods][kNumberOfRadii];
  TH2F *hTrueHtcResponse[kNumberOfRadii],*hHtcResponse[kNumberOfRadii];
  
  TH1A *hATrueScResponse[kNumberOfRadii],*hAScResponse[kNumberOfRadii];
  TH1A *hATrueTcResponse[kNumberOfRadii],*hATcResponse[Array::kNumberOfNtcMethods][kNumberOfRadii];
  TH1A *hAStcResponse[kNumberOfRadii],*hAFtcResponse[kNumberOfRadii];
  TH1A *hATrueCtcResponse[kNumberOfRadii],*hACtcResponse[Array::kNumberOfCtcMethods][kNumberOfRadii];
  TH1A *hATrueHtcResponse[kNumberOfRadii],*hAHtcResponse[kNumberOfRadii];
  */
  
  for(unsigned m(0);m<Array::kNumberOfScMethods;m++) {
    for(unsigned r(0);r<kNumberOfRadii;r++) {
      std::ostringstream sout;
      sout << "Method" << std::setw(2) << std::setfill('0') << m << "Radius" << std::setw(2) << std::setfill('0') << r;
      
      hScResponse[m][r]=new TH2F((std::string("ScResponse")+sout.str()).c_str(),
				 ";Energy (MIPs);Energy (MIPs)",
				 100,0.0,50*nrg,140,0.0,70*nrg);
      hAScResponse[m][r]=new TH1A((std::string("ScResponse")+sout.str()).c_str(),
				  ";Energy (MIPs);Energy (MIPs)",
				  100,0.0,50*nrg);
    }
  }

  for(unsigned m(0);m<Array::kNumberOfNtcMethods;m++) {
    for(unsigned r(0);r<kNumberOfRadii;r++) {
      std::ostringstream sout;
      sout << "Method" << std::setw(2) << std::setfill('0') << m << "Radius" << std::setw(2) << std::setfill('0') << r;
      
      hNtcResponse[0][m][r]=new TH2F((std::string("NtcResponse0")+sout.str()).c_str(),
				  ";Energy (MIPs);Energy (MIPs)",
				  100,0.0,50*nrg,140,0.0,70*nrg);
      hANtcResponse[0][m][r]=new TH1A((std::string("NtcResponse0")+sout.str()).c_str(),
				  ";Energy (MIPs);Energy (MIPs)",
				   100,0.0,50*nrg);
      hNtcResponse[1][m][r]=new TH2F((std::string("NtcResponse1")+sout.str()).c_str(),
				  ";Energy (MIPs);Energy (MIPs)",
				  100,0.0,50*nrg,140,0.0,70*nrg);
      hANtcResponse[1][m][r]=new TH1A((std::string("NtcResponse1")+sout.str()).c_str(),
				  ";Energy (MIPs);Energy (MIPs)",
				   100,0.0,50*nrg);
    }
  }

    /*
    hStcResponse[r]=new TH2F((std::string("StcResponse")+sout.str()).c_str(),
			     ";Energy (MIPs);Energy (MIPs)",
			     100,0.0,50*nrg,140,0.0,70*nrg);
    hAStcResponse[r]=new TH1A((std::string("StcResponse")+sout.str()).c_str(),
			      ";Energy (MIPs);Energy (MIPs)",
			      100,0.0,50*nrg);
    hFtcResponse[r]=new TH2F((std::string("FtcResponse")+sout.str()).c_str(),
			     ";Energy (MIPs);Energy (MIPs)",
			     100,0.0,50*nrg,140,0.0,70*nrg);
    hAFtcResponse[r]=new TH1A((std::string("FtcResponse")+sout.str()).c_str(),
			      ";Energy (MIPs);Energy (MIPs)",
			      100,0.0,50*nrg);

    hTrueCtcResponse[r]=new TH2F((std::string("TrueCtcResponse")+sout.str()).c_str(),
				";Energy (MIPs);Energy (MIPs)",
				100,0.0,50*nrg,140,0.0,70*nrg);
    hATrueCtcResponse[r]=new TH1A((std::string("TrueCtcResponse")+sout.str()).c_str(),
				 ";Energy (MIPs);Energy (MIPs)",
				 100,0.0,50*nrg);

    for(unsigned m(0);m<Array::kNumberOfCtcMethods;m++) {
      std::ostringstream sout2;
      sout2 << "Method" << m;
      
      hCtcResponse[m][r]=new TH2F((std::string("CtcResponse")+sout2.str()+sout.str()).c_str(),
				  ";Energy (MIPs);Energy (MIPs)",
				  100,0.0,50*nrg,140,0.0,70*nrg);
      hACtcResponse[m][r]=new TH1A((std::string("CtcResponse")+sout2.str()+sout.str()).c_str(),
				   ";Energy (MIPs);Energy (MIPs)",
				   100,0.0,50*nrg);
    }

    hTrueHtcResponse[r]=new TH2F((std::string("TrueHtcResponse")+sout.str()).c_str(),
				";Energy (MIPs);Energy (MIPs)",
				100,0.0,50*nrg,140,0.0,70*nrg);
    hATrueHtcResponse[r]=new TH1A((std::string("TrueHtcResponse")+sout.str()).c_str(),
				 ";Energy (MIPs);Energy (MIPs)",
				 100,0.0,50*nrg);

    hHtcResponse[r]=new TH2F((std::string("HtcResponse")+sout.str()).c_str(),
			    ";Energy (MIPs);Energy (MIPs)",
			    100,0.0,50*nrg,140,0.0,70*nrg);
    hAHtcResponse[r]=new TH1A((std::string("HtcResponse")+sout.str()).c_str(),
			     ";Energy (MIPs);Energy (MIPs)",
			     100,0.0,50*nrg);
    */


  
  TH1F *hScAdcSpectrum=new TH1F("ScAdcSpectrum",";ADC units",100,0.0,100.0);
  TH1F *hScAdcEnergySpectrum=new TH1F("ScAdcEnergySpectrum",";Energy (MIPs)",100,0.0,10.0);

  TH2F *hScAdcVsEnergy0=new TH2F("ScAdcVsEnergy0",";Energy (MIPs);Energy (MIPs)",
				 1100,-2.0,20.0,200,0.0,20.0);
  TH2F *hScAdcVsEnergy1=new TH2F("ScAdcVsEnergy1",";Energy (MIPs);Energy (MIPs)",
				 1000,0.0,1000.0,1000,0.0,1000.0);
  TH2F *hScAdcVsEnergy2=new TH2F("ScAdcVsEnergy2",";Energy (MIPs);Energy (MIPs)",
				 1000,0.0,100000.0,1000,0.0,100000.0);
  TH2F *hScAdcVsEnergy3=new TH2F("ScAdcVsEnergy3",";Energy (MIPs);Energy (MIPs)",
				 1000,0.0,1000000.0,1000,0.0,1000000.0);


  TH1F *hTcReadout=new TH1F("TcReadout",";Number of trigger cells",nScEventTotal,0,nScEventTotal);

  TH1F *hTcAdcRatio0=new TH1F("TcAdcRatio0",
			      ";Side ratio",
			      110,0.0,1.1);
  TH2F *hTcAdcRatio0VsStcAdc=new TH2F("TcAdcRatio0VsStcAdc",
				      ";Energy (ADC);Side ratio",
				      200,0.0,20000.0,110,0.0,1.1);
  TH1F *hTcAdcRatio1=new TH1F("TcAdcRatio1",
			      ";Opposite ratio",
			      110,0.0,1.1);
  TH2F *hTcAdcRatio1VsStcAdc=new TH2F("TcAdcRatio1VsStcAdc",
				      ";Energy (ADC);Opposite ratio",
				      200,0.0,20000.0,110,0.0,1.1);
  TH1F *hTcAdcRatio2=new TH1F("TcAdcRatio2",
			      ";Opposite ratio",
			      110,0.0,1.1);
  TH2F *hTcAdcRatio2VsStcAdc=new TH2F("TcAdcRatio2VsStcAdc",
				      ";Energy (ADC);Opposite ratio",
				      200,0.0,20000.0,110,0.0,1.1);
  
  TH1F *hTcAdcEnergySpectrum[Array::kNumberOfNtcMethods];
  TH2F *hTcAdcVsEnergy0[Array::kNumberOfNtcMethods];
  TH2F *hTcAdcVsEnergy1[Array::kNumberOfNtcMethods];
  TH2F *hTcAdcVsEnergy2[Array::kNumberOfNtcMethods];
  TH2F *hTcAdcVsEnergy3[Array::kNumberOfNtcMethods];

  for(unsigned m(0);m<Array::kNumberOfNtcMethods;m++) {
    std::ostringstream sout2;
    sout2 << "Method" << std::setw(2) << std::setfill('0') << m;
    
    hTcAdcEnergySpectrum[m]=new TH1F((sout2.str()+"TcAdcEnergySpectrum").c_str(),
				     ";Energy (MIPs)",100,0.0,100.0);
    
    hTcAdcVsEnergy0[m]=new TH2F((sout2.str()+"TcAdcVsEnergy0").c_str(),
				";Energy (MIPs);Energy (MIPs)",
				1000,0.0,20.0,200,0.0,20.0);
    hTcAdcVsEnergy1[m]=new TH2F((sout2.str()+"TcAdcVsEnergy1").c_str(),
				";Energy (MIPs);Energy (MIPs)",
				1000,0.0,100.0,1000,0.0,100.0);
    hTcAdcVsEnergy2[m]=new TH2F((sout2.str()+"TcAdcVsEnergy2").c_str(),
				";Energy (MIPs);Energy (MIPs)",
				1000,0.0,1000.0,1000,0.0,1000.0);
    hTcAdcVsEnergy3[m]=new TH2F((sout2.str()+"TcAdcVsEnergy3").c_str(),
				";Energy (MIPs);Energy (MIPs)",
				1000,0.0,10000.0,1000,0.0,10000.0);
  }
  /*  
  TH1F *hStcAdcEnergySpectrum=new TH1F("StcAdcEnergySpectrum",";Energy (MIPs)",100,0.0,100.0);
  TH2F *hStcAdcVsEnergy0=new TH2F("StcAdcVsEnergy0",";Energy (MIPs);Energy (MIPs)",
				  1000,0.0,20.0,200,0.0,20.0);
  TH2F *hStcAdcVsEnergy1=new TH2F("StcAdcVsEnergy1",";Energy (MIPs);Energy (MIPs)",
				  1000,0.0,100.0,1000,0.0,100.0);
  TH2F *hStcAdcVsEnergy2=new TH2F("StcAdcVsEnergy2",";Energy (MIPs);Energy (MIPs)",
				  1000,0.0,1000.0,1000,0.0,1000.0);
  TH2F *hStcAdcVsEnergy3=new TH2F("StcAdcVsEnergy3",";Energy (MIPs);Energy (MIPs)",
				  1000,0.0,10000.0,1000,0.0,10000.0);

  TH1F *hFtcAdcEnergySpectrum=new TH1F("FtcAdcEnergySpectrum",";Energy (MIPs)",100,0.0,100.0);
  TH2F *hFtcAdcVsEnergy0=new TH2F("FtcAdcVsEnergy0",";Energy (MIPs);Energy (MIPs)",
				  1000,0.0,20.0,200,0.0,20.0);
  TH2F *hFtcAdcVsEnergy1=new TH2F("FtcAdcVsEnergy1",";Energy (MIPs);Energy (MIPs)",
				  1000,0.0,100.0,1000,0.0,100.0);
  TH2F *hFtcAdcVsEnergy2=new TH2F("FtcAdcVsEnergy2",";Energy (MIPs);Energy (MIPs)",
				  1000,0.0,1000.0,1000,0.0,1000.0);
  TH2F *hFtcAdcVsEnergy3=new TH2F("FtcAdcVsEnergy3",";Energy (MIPs);Energy (MIPs)",
				  1000,0.0,10000.0,1000,0.0,10000.0);

  TH1F *hCtcAdcEnergySpectrum=new TH1F("CtcAdcEnergySpectrum",";Energy (MIPs)",100,0.0,100.0);
  TH2F *hCtcAdcVsEnergy0=new TH2F("CtcAdcVsEnergy0",";Energy (MIPs);Energy (MIPs)",
				  1000,0.0,20.0,200,0.0,20.0);
  TH2F *hCtcAdcVsEnergy1=new TH2F("CtcAdcVsEnergy1",";Energy (MIPs);Energy (MIPs)",
				  1000,0.0,100.0,1000,0.0,100.0);
  TH2F *hCtcAdcVsEnergy2=new TH2F("CtcAdcVsEnergy2",";Energy (MIPs);Energy (MIPs)",
				  1000,0.0,1000.0,1000,0.0,1000.0);
  TH2F *hCtcAdcVsEnergy3=new TH2F("CtcAdcVsEnergy3",";Energy (MIPs);Energy (MIPs)",
				  1000,0.0,10000.0,1000,0.0,10000.0);

  */
  
  unsigned nScTotal(0);
  unsigned nScReadout(0);
  unsigned nTcTotal(0);
  unsigned nTcReadout(0);

  double centre[2];

  double nScMs[2]={0.0,0.0};
  double xScMs[2]={0.0,0.0}; 

  double nTcMs[Array::kNumberOfNtcMethods]={0.0,0.0,0.0,0.0,0.0};
  double xTcMs[Array::kNumberOfNtcMethods]={0.0,0.0,0.0,0.0,0.0}; 
  /*
  double nCtcMs[Array::kNumberOfCtcMethods]={0.0,0.0,0.0,0.0};
  double xCtcMs[Array::kNumberOfCtcMethods]={0.0,0.0,0.0,0.0};

  double nHtcMs[2]={0.0,0.0};
  double xHtcMs[2]={0.0,0.0}; 
  */
  
  // KEEP!
  double scResponse[Array::kNumberOfScMethods][kNumberOfRadii];
  double ntcResponse[2][Array::kNumberOfNtcMethods][kNumberOfRadii];

  
  for(unsigned nEvt(0);nEvt<10;nEvt++) {
    centre[0]=theArray.fCoarseCellSide*(Random::random().Uniform()-0.5);
    centre[1]=theArray.fCoarseCellSide*(Random::random().Uniform()-0.5);
    hCentre->Fill(centre[0],centre[1]);

    unsigned nScEventReadout(0);
    unsigned nTcEventReadout(0);

    // Process the event
    theArray.initialise();

    double sigma(2.0);
    double sigma2(2.0*sigma*sigma);
    
    double trueEnergy(0.0);
    double thisEnergy(fMipsPerGeV*fPt*cosh(theArray.fEta));
    //double thisEnergy((0.1+0.9*Random::random().Uniform())*nrg);

    unsigned modelBins(4);
    double modelBinsSize(theArray.fSensorCellSide/modelBins);
    double modelBinsArea(modelBinsSize*modelBinsSize);
    
    for(unsigned ix(0);ix<theArray.numberOfSensorCells();ix++) {
      for(unsigned iy(0);iy<theArray.numberOfSensorCells();iy++) {
	
	double e(0.0);
	for(unsigned nx(0);nx<modelBins;nx++) {
	  double dx(theArray.sensorCellPosition(ix)+(nx+0.5-(modelBins/2))*modelBinsSize-centre[0]);
	  
	  for(unsigned ny(0);ny<modelBins;ny++) {
	    double dy(theArray.sensorCellPosition(iy)+(ny+0.5-(modelBins/2))*modelBinsSize-centre[1]);
	    double dr(sqrt(dx*dx+dy*dy));
	  
	    //e+=thisEnergy*exp(-(dx*dx+dy*dy)/sigma2)); // Gaussian
	    //e+=thisEnergy*exp(-sqrt(dx*dx+dy*dy)/4.0)); // Exponential
	    double de(thisEnergy*thePhotonModel[usePhotonModel].density(dr)*modelBinsArea);
	    hScEnergyIntegral->Fill(dr,de);
	    
	    e+=de;
	  }
	}

	e*=std::max(Random::random().Gaus(1.0,fPtSpread),0.0);
	trueEnergy+=e;
	theArray.increment(ix,iy,e);
      }
    }

    if(nEvt<10) std::cout << "This, true energy = " << thisEnergy << ", " << trueEnergy << " MIPs" << std::endl;
    hScEnergyTrue->Fill(trueEnergy);    
    
    theArray.digitise(0);
    if(nEvt<10) std::cout << "Completed digitisation" << std::endl;
    
    // Analyse the event
    if(nEvt<5) histEvent(nEvt,theArray);

    for(unsigned r(0);r<kNumberOfRadii;r++) {
      for(unsigned m(0);m<Array::kNumberOfScMethods;m++) scResponse[m][r]=0.0;
      for(unsigned m(0);m<Array::kNumberOfNtcMethods;m++) {
	ntcResponse[0][m][r]=0.0;
	ntcResponse[1][m][r]=0.0;
      }
    }
    
    double trueTcR[kNumberOfRadii],tcR[Array::kNumberOfNtcMethods][kNumberOfRadii];
    double stcR[kNumberOfRadii],ftcR[kNumberOfRadii];
    //double trueCtcR[kNumberOfRadii],ctcR[Array::kNumberOfCtcMethods][kNumberOfRadii];
    //double trueHtcR[kNumberOfRadii],htcR[kNumberOfRadii];
    /*
    for(unsigned r(0);r<kNumberOfRadii;r++) {
      trueScR[r]=0.0;
      scR[r]=0.0;

      trueTcR[r]=0.0;
      for(unsigned m(0);m<Array::kNumberOfNtcMethods;m++) tcR[m][r]=0.0;

      stcR[r]=0.0;
      ftcR[r]=0.0;

      trueCtcR[r]=0.0;
      for(unsigned m(0);m<Array::kNumberOfCtcMethods;m++) ctcR[m][r]=0.0;

      trueHtcR[r]=0.0;
      htcR[r]=0.0;
    }
    */
    
    for(unsigned ix(0);ix<theArray.numberOfSensorCells();ix++) {
      for(unsigned iy(0);iy<theArray.numberOfSensorCells();iy++) {
	if(theArray.fScAdc[Array::kScReal][ix][iy]>0) nScEventReadout++;

	unsigned m=Array::kScReal;
	  
	hScEnergySpectrum->Fill(theArray.scTotalEnergy(ix,iy));
	hScAdcSpectrum->Fill(theArray.fScAdc[m][ix][iy]);
	hScAdcEnergySpectrum->Fill(theArray.scAdcEnergy(ix,iy,m));
	
	hScAdcVsEnergy0->Fill(theArray.scTotalEnergy(ix,iy),theArray.scAdcEnergy(ix,iy,m));
	hScAdcVsEnergy1->Fill(theArray.scTotalEnergy(ix,iy),theArray.scAdcEnergy(ix,iy,m));
	hScAdcVsEnergy2->Fill(theArray.scTotalEnergy(ix,iy),theArray.scAdcEnergy(ix,iy,m));
	hScAdcVsEnergy3->Fill(theArray.scTotalEnergy(ix,iy),theArray.scAdcEnergy(ix,iy,m));

	double dx(theArray.sensorCellPosition(ix)-centre[0]);
	double dy(theArray.sensorCellPosition(iy)-centre[1]);
	double dr(sqrt(dx*dx+dy*dy));

	for(unsigned m(0);m<Array::kNumberOfScMethods;m++) {
	  for(unsigned r(0);r<kNumberOfRadii;r++) {
	    if(dr<fLayerZ*cutRadius(r)) {
	      scResponse[m][r]+=theArray.scAdcEnergy(ix,iy,m);
	    }
	  }
	}

	/*
	if(dr<7.5) {
	  nScMs[0]+=2.0*theArray.scTrueEnergy(ix,iy);
	  xScMs[0]+=(dx*dx+dy*dy)*theArray.scTrueEnergy(ix,iy);
	  nScMs[1]+=2.0*theArray.scAdcEnergy(ix,iy);
	  xScMs[1]+=(dx*dx+dy*dy)*theArray.scAdcEnergy(ix,iy);
	}
	*/
      }
    }

    if(nEvt<10) std::cout << "Completed sensor cells" << std::endl;

    hScReadout->Fill(nScEventReadout);

    nScTotal+=nScEventTotal;
    nScReadout+=nScEventReadout;

    for(unsigned m(0);m<Array::kNumberOfScMethods;m++) {
      for(unsigned r(0);r<kNumberOfRadii;r++) {
	hScResponse[m][r]->Fill(trueEnergy,scResponse[m][r]);
	hAScResponse[m][r]->Fill(trueEnergy,scResponse[m][r]);
      }
    }
    
    ////////////////////////////////////////////////////////////////

    for(unsigned ix(0);ix<theArray.numberOfTriggerCells();ix++) {
      for(unsigned iy(0);iy<theArray.numberOfTriggerCells();iy++) {
	if(theArray.ntcAdcEnergy(ix,iy,Array::kNtcReal)>0.0) nTcEventReadout++;
	/*
	if(theArray.fStcAdc[ix][iy]>0) {
	  unsigned nx((ix%2)==0?ix+1:ix-1);
	  unsigned ny((iy%2)==0?iy+1:iy-1);
	  assert(theArray.fStcAdc[nx][iy]==0);
	  assert(theArray.fStcAdc[ix][ny]==0);
	  assert(theArray.fStcAdc[nx][ny]==0);
	  
	  hTcAdcRatio0->Fill(1.0*theArray.fTcAdc[ix][iy]/theArray.fStcAdc[ix][iy]);
	  hTcAdcRatio0VsStcAdc->Fill(theArray.fStcAdc[ix][iy],1.0*theArray.fTcAdc[ix][iy]/theArray.fStcAdc[ix][iy]);

	  hTcAdcRatio1->Fill(1.0*theArray.fTcAdc[nx][iy]/theArray.fTcAdc[ix][iy]);
	  hTcAdcRatio1->Fill(1.0*theArray.fTcAdc[ix][ny]/theArray.fTcAdc[ix][iy]);
	  hTcAdcRatio1VsStcAdc->Fill(theArray.fStcAdc[ix][iy],1.0*theArray.fTcAdc[nx][iy]/theArray.fTcAdc[ix][iy]);
	  hTcAdcRatio1VsStcAdc->Fill(theArray.fStcAdc[ix][iy],1.0*theArray.fTcAdc[ix][ny]/theArray.fTcAdc[ix][iy]);

	  hTcAdcRatio2->Fill(1.0*theArray.fTcAdc[nx][ny]/theArray.fTcAdc[ix][iy]);
	  hTcAdcRatio2VsStcAdc->Fill(theArray.fStcAdc[ix][iy],1.0*theArray.fTcAdc[nx][ny]/theArray.fTcAdc[ix][iy]);
	}
	*/
	for(unsigned m(0);m<Array::kNumberOfNtcMethods;m++) {
	  hTcAdcEnergySpectrum[m]->Fill(theArray.tcAdcEnergy(ix,iy,m));
	  hTcAdcVsEnergy0[m]->Fill(theArray.tcTotalEnergy(ix,iy),theArray.tcAdcEnergy(ix,iy,m));
	  hTcAdcVsEnergy1[m]->Fill(theArray.tcTotalEnergy(ix,iy),theArray.tcAdcEnergy(ix,iy,m));
	  hTcAdcVsEnergy2[m]->Fill(theArray.tcTotalEnergy(ix,iy),theArray.tcAdcEnergy(ix,iy,m));
	  hTcAdcVsEnergy3[m]->Fill(theArray.tcTotalEnergy(ix,iy),theArray.tcAdcEnergy(ix,iy,m));
	}
	/*	
	hStcAdcEnergySpectrum->Fill(theArray.stcAdcEnergy(ix,iy));
	hStcAdcVsEnergy0->Fill(theArray.tcTotalEnergy(ix,iy),theArray.stcAdcEnergy(ix,iy));
	hStcAdcVsEnergy1->Fill(theArray.tcTotalEnergy(ix,iy),theArray.stcAdcEnergy(ix,iy));
	hStcAdcVsEnergy2->Fill(theArray.tcTotalEnergy(ix,iy),theArray.stcAdcEnergy(ix,iy));
	hStcAdcVsEnergy3->Fill(theArray.tcTotalEnergy(ix,iy),theArray.stcAdcEnergy(ix,iy));

	hFtcAdcEnergySpectrum->Fill(theArray.ftcAdcEnergy(ix,iy));
	hFtcAdcVsEnergy0->Fill(theArray.tcTotalEnergy(ix,iy),theArray.ftcAdcEnergy(ix,iy));
	hFtcAdcVsEnergy1->Fill(theArray.tcTotalEnergy(ix,iy),theArray.ftcAdcEnergy(ix,iy));
	hFtcAdcVsEnergy2->Fill(theArray.tcTotalEnergy(ix,iy),theArray.ftcAdcEnergy(ix,iy));
	hFtcAdcVsEnergy3->Fill(theArray.tcTotalEnergy(ix,iy),theArray.ftcAdcEnergy(ix,iy));
	*/

	double dx;
	double dy;

	for(unsigned ic(0);ic<2;ic++) {
	  if(ic==0) {
	    dx=theArray.triggerCellPosition(ix)-centre[0];
	    dy=theArray.triggerCellPosition(iy)-centre[1];
	  } else {
	    unsigned icx(theArray.triggerCellNumber(centre[0]));
	    unsigned icy(theArray.triggerCellNumber(centre[1]));
	    dx=theArray.triggerCellPosition(ix)-theArray.triggerCellPosition(icx);
	    dy=theArray.triggerCellPosition(iy)-theArray.triggerCellPosition(icy);
	  }
	  double dr(sqrt(dx*dx+dy*dy));
	  
	  for(unsigned m(0);m<Array::kNumberOfNtcMethods;m++) {
	    for(unsigned r(0);r<kNumberOfRadii;r++) {
	      if(dr<fLayerZ*cutRadius(r)) {
		ntcResponse[ic][m][r]+=theArray.ntcAdcEnergy(ix,iy,m);
	      }
	    }
	  }
	}
	/*
	if(dr<8.0) {
	  for(unsigned m(0);m<Array::kNumberOfNtcMethods;m++) {
	    nTcMs[m]+=2.0*theArray.tcAdcEnergy(ix,iy,m);
	    xTcMs[m]+=(dx*dx+dy*dy)*theArray.tcAdcEnergy(ix,iy,m);
	  }
	}
	*/
      }
    }
    
    hTcReadout->Fill(nTcEventReadout);

    nTcTotal+=nTcEventTotal;
    nTcReadout+=nTcEventReadout;

    for(unsigned m(0);m<Array::kNumberOfNtcMethods;m++) {
      for(unsigned r(0);r<kNumberOfRadii;r++) {
	hNtcResponse[0][m][r]->Fill(trueEnergy,ntcResponse[0][m][r]);
	hANtcResponse[0][m][r]->Fill(trueEnergy,ntcResponse[0][m][r]);
	hNtcResponse[1][m][r]->Fill(trueEnergy,ntcResponse[1][m][r]);
	hANtcResponse[1][m][r]->Fill(trueEnergy,ntcResponse[1][m][r]);
      }
    }

    ////////////////////////////////////////////////////////////////
    
#ifdef CRAP
    for(unsigned ix(0);ix<theArray.numberOfCoarseCells();ix++) {
      for(unsigned iy(0);iy<theArray.numberOfCoarseCells();iy++) {

	for(unsigned m(0);m<Array::kNumberOfCtcMethods;m++) {
	  hCtcAdcEnergySpectrum->Fill(theArray.ctcAdcEnergy(ix,iy,m));
	  hCtcAdcVsEnergy0->Fill(theArray.tcTotalEnergy(ix,iy),theArray.ctcAdcEnergy(ix,iy,m));
	  hCtcAdcVsEnergy1->Fill(theArray.tcTotalEnergy(ix,iy),theArray.ctcAdcEnergy(ix,iy,m));
	  hCtcAdcVsEnergy2->Fill(theArray.tcTotalEnergy(ix,iy),theArray.ctcAdcEnergy(ix,iy,m));
	  hCtcAdcVsEnergy3->Fill(theArray.tcTotalEnergy(ix,iy),theArray.ctcAdcEnergy(ix,iy,m));
	}
	
	double dx(theArray.coarseCellPosition(ix)-centre[0]);
	double dy(theArray.coarseCellPosition(iy)-centre[1]);
	double dr(sqrt(dx*dx+dy*dy));

	for(unsigned r(0);r<kNumberOfRadii;r++) {
	  if(dr<fLayerZ*cutRadius(r)) {
	    trueCtcR[r]+=theArray.ctcTotalEnergy(ix,iy);
	    for(unsigned m(0);m<Array::kNumberOfCtcMethods;m++) ctcR[m][r]+=theArray.ctcAdcEnergy(ix,iy,m);
	  }
	}

	if(dr<8.0) {
	  for(unsigned m(0);m<Array::kNumberOfCtcMethods;m++) {
	    nCtcMs[m]+=2.0*theArray.ctcAdcEnergy(ix,iy,m);
	    xCtcMs[m]+=(dx*dx+dy*dy)*theArray.ctcAdcEnergy(ix,iy,m);
	  }
	}
      }
    }

    for(unsigned r(0);r<kNumberOfRadii;r++) {
      hTrueCtcResponse[r]->Fill(trueEnergy,trueCtcR[r]);
      hATrueCtcResponse[r]->Fill(trueEnergy,trueCtcR[r]);

      for(unsigned m(0);m<Array::kNumberOfCtcMethods;m++) {
	hCtcResponse[m][r]->Fill(trueEnergy,ctcR[m][r]);
	hACtcResponse[m][r]->Fill(trueEnergy,ctcR[m][r]);
      }
    }

    ////////////////////////////////////////////////////////////////
    
    for(unsigned ix(0);ix<theArray.numberOfHgcrocXCells();ix++) {
      for(unsigned iy(0);iy<theArray.numberOfHgcrocYCells();iy++) {
	double dx(theArray.hgcrocXCellPosition(ix)-centre[0]);
	double dy(theArray.hgcrocYCellPosition(iy)-centre[1]);
	double dr(sqrt(dx*dx+dy*dy));

	for(unsigned r(0);r<kNumberOfRadii;r++) {
	  if(dr<fLayerZ*cutRadius(r)) {
	    trueHtcR[r]+=theArray.htcTotalEnergy(ix,iy);
	    htcR[r]+=theArray.htcAdcEnergy(ix,iy);
	  }
	}

	if(dr<8.0) {
	  for(unsigned m(0);m<1;m++) {
	    nHtcMs[m]+=2.0*theArray.htcAdcEnergy(ix,iy);
	    xHtcMs[m]+=(dx*dx+dy*dy)*theArray.htcAdcEnergy(ix,iy);
	  }
	}
      }
    }
    
    for(unsigned r(0);r<kNumberOfRadii;r++) {
      hTrueHtcResponse[r]->Fill(trueEnergy,trueHtcR[r]);
      hATrueHtcResponse[r]->Fill(trueEnergy,trueHtcR[r]);

      hHtcResponse[r]->Fill(trueEnergy,htcR[r]);
      hAHtcResponse[r]->Fill(trueEnergy,htcR[r]);
    }
#endif

  } // End of event loop

  ////////////////////////////////////////////////////////////////
    
  std::cout << "Total, read out sensor cells = " << nScTotal
	    << ", " << nScReadout << ", occupancy = "
	    << 1.0*nScReadout/nScTotal << std::endl;
  std::cout << "Total, read out trigger cells = " << nTcTotal
	    << ", " << nTcReadout << ", occupancy = "
	    << 1.0*nTcReadout/nTcTotal << std::endl;
  std::cout << std::endl;

  double trueRms(sqrt(xScMs[0]/nScMs[0]));

  std::cout << "Number, sumSq and RMS for true energy sensor cells = " << nScMs[0]
	    << ", " << xScMs[0] << ", " << sqrt(xScMs[0]/nScMs[0])
	    << ", " << sqrt(xScMs[0]/nScMs[0])/trueRms << std::endl;
  std::cout << "Number, sumSq and RMS for ADC energy sensor cells = " << nScMs[1]
	    << ", " << xScMs[1] << ", " << sqrt(xScMs[1]/nScMs[1])
	    << ", " << sqrt(xScMs[1]/nScMs[1])/trueRms << std::endl;
  std::cout << std::endl;

  for(unsigned m(0);m<Array::kNumberOfNtcMethods;m++) {
    std::cout << "Number, sumSq and RMS for method " << m << " TC energy sensor cells = " << nTcMs[m]
	      << ", " << xTcMs[m] << ", " << sqrt(xTcMs[m]/nTcMs[m])
	      << ", " << sqrt(xTcMs[m]/nTcMs[m])/trueRms << std::endl;
  }
  std::cout << std::endl;

  /*
  for(unsigned m(0);m<Array::kNumberOfCtcMethods;m++) {
    std::cout << "Number, sumSq and RMS for method " << m << " CTC energy sensor cells = " << nCtcMs[m]
	      << ", " << xCtcMs[m] << ", " << sqrt(xCtcMs[m]/nCtcMs[m])
	      << ", " << sqrt(xCtcMs[m]/nCtcMs[m])/trueRms << std::endl;
  }
  std::cout << std::endl;

  for(unsigned m(0);m<1;m++) {
    std::cout << "Number, sumSq and RMS for method " << m << " HTC energy sensor cells = " << nHtcMs[m]
	      << ", " << xHtcMs[m] << ", " << sqrt(xHtcMs[m]/nHtcMs[m])
	      << ", " << sqrt(xHtcMs[m]/nHtcMs[m])/trueRms << std::endl;
  }
  */
  
  if(false) {
    double parAvg[10],parRms[10];
    std::string sFit("pol1");

    double minRms(1.0e12);
    double minR(-1.0);

    TH1F *hScEsVsR[Array::kNumberOfScMethods];
    TH1F *hScClVsR[Array::kNumberOfScMethods];
    TH1F *hScPuVsR[Array::kNumberOfScMethods];

    for(unsigned m(0);m<Array::kNumberOfScMethods;m++) {
      std::ostringstream sout;
      sout << "ScMethod" << m;

      hScEsVsR[m]=new TH1F((sout.str()+"EsVsR").c_str(),
			   ";Cut radius #rho;Energy resolution (MIPs)",
			    kNumberOfRadii,0.00875,0.05875);
      hScClVsR[m]=new TH1F((sout.str()+"ClVsR").c_str(),
			   ";Cut radius #rho;Calibration (ADC/MIP)",
			    kNumberOfRadii,0.00875,0.05875);
      hScPuVsR[m]=new TH1F((sout.str()+"PuVsR").c_str(),
			   ";Cut radius #rho;Pileup (MIPs)",
			   kNumberOfRadii,0.00875,0.05875);
      
      minRms=1.0e12;
      minR=-1.0;

      for(unsigned r(0);r<kNumberOfRadii;r++) {
	TH1F *hAvg=hAScResponse[m][r]->Avg();
	TH1F *hRms=hAScResponse[m][r]->Rms();
	
	std::cout << "Sc method " << m << ", radius " << fLayerZ*cutRadius(r) << ", entries and integral = "
		  << hAvg->GetEntries() << ", " << hAvg->Integral() << std::endl;
								       
	if(hAvg->Integral()>100) {
	hAvg->Fit("pol1");
	hAvg->GetFunction("pol1")->GetParameters(parAvg);
	
	hScClVsR[m]->SetBinContent(r+1,parAvg[1]);
	hScPuVsR[m]->SetBinContent(r+1,parAvg[0]/parAvg[1]);
	
	hRms->Fit(sFit.c_str());
	hRms->GetFunction(sFit.c_str())->GetParameters(parRms);
	
	double typical=(*(hRms->GetFunction(sFit.c_str())))(4000)/parAvg[1];
	std::cout << "Sc method " << m << ", radius " << fLayerZ*cutRadius(r) << ", typical RMS = " << typical;
	
	hScEsVsR[m]->SetBinContent(r+1,typical);

	hAvg->RecursiveRemove(hAvg->GetFunction("pol1"));
	hRms->RecursiveRemove(hRms->GetFunction(sFit.c_str()));
    
	if(minRms>typical) {
	  minRms=typical;
	  minR=fLayerZ*cutRadius(r);
	  std::cout << " new min" << std::endl;
	} else {
	  std::cout << std::endl;
	}
	}
      }
      std::cout << "Sc, method " << m << ", min radius " << minR << ", typical RMS = " << minRms << " <<<=======" << std::endl << std::endl;
    }


    TH1F *hNtcEsVsR[2][Array::kNumberOfNtcMethods];
    TH1F *hNtcClVsR[2][Array::kNumberOfNtcMethods];
    TH1F *hNtcPuVsR[2][Array::kNumberOfNtcMethods];

    for(unsigned m(0);m<Array::kNumberOfNtcMethods;m++) {
      std::ostringstream sout;
      sout << "NtcMethod" << std::setw(2) << std::setfill('0') << m;

      hNtcEsVsR[0][m]=new TH1F((sout.str()+"EsVsR0").c_str(),
			       ";Cut radius #rho;Energy resolution (MIPs)",
			       //kNumberOfRadii,0.00875,0.05875);
			       kNumberOfRadii,0.0,0.06);
      hNtcClVsR[0][m]=new TH1F((sout.str()+"ClVsR0").c_str(),
			       ";Cut radius #rho;Calibration (ADC/MIP)",
			       kNumberOfRadii,0.0,0.06);
      hNtcPuVsR[0][m]=new TH1F((sout.str()+"PuVsR0").c_str(),
			       ";Cut radius #rho;Pileup (MIPs)",
			       kNumberOfRadii,0.0,0.06);
      hNtcEsVsR[1][m]=new TH1F((sout.str()+"EsVsR1").c_str(),
			       ";Cut radius #rho;Energy resolution (MIPs)",
			       kNumberOfRadii,0.0,0.06);
      hNtcClVsR[1][m]=new TH1F((sout.str()+"ClVsR1").c_str(),
			       ";Cut radius #rho;Calibration (ADC/MIP)",
			       kNumberOfRadii,0.0,0.06);
      hNtcPuVsR[1][m]=new TH1F((sout.str()+"PuVsR1").c_str(),
			       ";Cut radius #rho;Pileup (MIPs)",
			       kNumberOfRadii,0.0,0.06);

      for(unsigned ic(0);ic<2;ic++) {
	minRms=1.0e12;
	minR=-1.0;
	
	for(unsigned r(0);r<kNumberOfRadii;r++) {
	  TH1F *hAvg=hANtcResponse[ic][m][r]->Avg();
	  TH1F *hRms=hANtcResponse[ic][m][r]->Rms();
	  
	  std::cout << "Ntc method " << m << ", radius " << fLayerZ*cutRadius(r) << ", entries and integral = "
		    << hAvg->GetEntries() << ", " << hAvg->Integral() << std::endl;
	  
	  if(hAvg->Integral()>100) {
	    hAvg->Fit("pol1");
	    hAvg->GetFunction("pol1")->GetParameters(parAvg);
	    
	    hNtcClVsR[ic][m]->SetBinContent(r+1,parAvg[1]);
	    hNtcPuVsR[ic][m]->SetBinContent(r+1,parAvg[0]/parAvg[1]);
	    
	    hRms->Fit(sFit.c_str());
	    hRms->GetFunction(sFit.c_str())->GetParameters(parRms);
	    
	    double typical=(*(hRms->GetFunction(sFit.c_str())))(4000)/parAvg[1];
	    std::cout << "Ntc method " << m << ", radius " << fLayerZ*cutRadius(r) << ", typical RMS = " << typical;
	    
	    hNtcEsVsR[ic][m]->SetBinContent(r+1,typical);
	    
	    hAvg->RecursiveRemove(hAvg->GetFunction("pol1"));
	    hRms->RecursiveRemove(hRms->GetFunction(sFit.c_str()));
	    
	    if(minRms>typical) {
	      minRms=typical;
	      minR=fLayerZ*cutRadius(r);
	      std::cout << " new min" << std::endl;
	    } else {
	      std::cout << std::endl;
	    }
	  }
	}
	
	std::cout << "Centre method " << ic << ", Ntc method " << m << ", min radius " << minR << ", typical RMS = " << minRms << " <<<=======" << std::endl << std::endl;
      }
    }

    /*
    for(unsigned m(0);m<Array::kNumberOfCtcMethods;m++) {
      minRms=1.0e12;
      minR=-1.0;
    
      std::ostringstream sout2;
      sout2 << "Method" << m;
    
      TH1F *hCtcEsVsR=new TH1F((std::string("Ctc")+sout2.str()+"EsVsR").c_str(),
			       ";Cut radius #rho;Energy resolution (MIPs)",
			       kNumberOfRadii,0.0075,0.0575);
      TH1F *hCtcPuVsR=new TH1F((std::string("Ctc")+sout2.str()+"PuVsR").c_str(),
			       ";Cut radius #rho;Pileup (MIPs)",
			       kNumberOfRadii,0.0075,0.0575);

      for(unsigned r(0);r<kNumberOfRadii;r++) {
	TH1F *hAvg=hACtcResponse[m][r]->Avg();
	TH1F *hRms=hACtcResponse[m][r]->Rms();

	hAvg->Fit("pol1");
	hAvg->GetFunction("pol1")->GetParameters(parAvg);
	hRms->Fit(sFit.c_str());
	hRms->GetFunction(sFit.c_str())->GetParameters(parRms);

	hCtcClVsR->SetBinContent(r+1,parAvg[0]);
	hCtcPuVsR->SetBinContent(r+1,parAvg[0]);
    
	double typical=(*(hRms->GetFunction(sFit.c_str())))(4000)/parAvg[1];
	std::cout << "Sc, radius " << fLayerZ*cutRadius(r) << ", typical RMS = " << typical;

	hCtcEsVsR->SetBinContent(r+1,typical);
	hAvg->RecursiveRemove(hAvg->GetFunction("pol1"));
	hRms->RecursiveRemove(hRms->GetFunction(sFit.c_str()));
    
	if(minRms>typical) {
	  minRms=typical;
	  minR=fLayerZ*cutRadius(r);
	  std::cout << " new min" << std::endl;
	} else {
	  std::cout << std::endl;
	}
      }
      std::cout << "Ctc method " << m << ", min radius " << minR << ", typical RMS = " << minRms << " <<<=======" << std::endl << std::endl;
    }


    minRms=1.0e12;
    minR=-1.0;
  
    TH1F *hHtcEsVsR=new TH1F("HtcEsVsR",";Cut radius #rho;Energy resolution (MIPs)",
			    kNumberOfRadii,0.0075,0.0575);
    TH1F *hHtcPuVsR=new TH1F("HtcPuVsR",";Cut radius #rho;Pileup (MIPs)",
			    kNumberOfRadii,0.0075,0.0575);

    for(unsigned r(0);r<kNumberOfRadii;r++) {
      TH1F *hAvg=hAHtcResponse[r]->Avg();
      TH1F *hRms=hAHtcResponse[r]->Rms();

      hAvg->Fit("pol1");
      hAvg->GetFunction("pol1")->GetParameters(parAvg);
      hRms->Fit(sFit.c_str());
      hRms->GetFunction(sFit.c_str())->GetParameters(parRms);

      hHtcPuVsR->SetBinContent(r+1,parAvg[0]);
    
      double typical=(*(hRms->GetFunction(sFit.c_str())))(4000)/parAvg[1];
      std::cout << "Htc, radius " << fLayerZ*cutRadius(r) << ", typical RMS = " << typical;

      hHtcEsVsR->SetBinContent(r+1,typical);

      hAvg->RecursiveRemove(hAvg->GetFunction("pol1"));
      hRms->RecursiveRemove(hRms->GetFunction(sFit.c_str()));
    
      if(minRms>typical) {
	minRms=typical;
	minR=fLayerZ*cutRadius(r);
	std::cout << " new min" << std::endl;
      } else {
	std::cout << std::endl;
      }
    }

    std::cout << "Htc, min radius " << minR << ", typical RMS = " << minRms << " <<<=======" << std::endl << std::endl;
  }
    */  
  }
  return 0;
}