This is the code for chi-squared program. There are definitely ways to make this more efficient.

// New Chi-squared Program


// This is a C++ program to be used with ROOT that subracts pedestals (run1060) from two quartz fiber
// runs, one with one fiber illuminated and the others with all 8 fibers illuminated and finds their 
// ratio. The average ratio (over all eta) for a particular phi bin is then calculated and the ratio of 
//one fiber to all fibers illuminated is divided by  the average ratio. The graphs show the chi-squared
// distribution for a fiber.
// chi-squared = (average of how much the ratio/avg ratio deviates from one).


//====================================
//====================================
//====================================

float meancalc(int cent,float yval[],int offset){
  float meanval=0;
  float N=0;
  int low=50;
  
  if (cent=0.02 || deviation[jj]>1)
	  printf("chi-squared= %f, phi bin %d%c\n",deviation[jj],kk,cch);
	jj++;
      }
}
	
//====================================
//====================================
//====================================

newchiSq() {
  int runa=1071;
  int runb=1078;
  int peds=1060;

  gStyle->SetOptStat(1111111);
  TCanvas *c1=new TCanvas("layer","Laser Data: Single quartz-fiber tests",800,600);  
  c1->SetBorderMode(0);

 
  TString hFilea="/N/u1/bullinger/outLaser/pi";
  hFilea+=runa;  hFilea+=".hist.root";
 
  TString hFileb="/N/u1/bullinger/outLaser/pi";
  hFileb+=runb; hFileb+=".hist.root";
  
   TString hFilec="/N/u1/bullinger/outLaser/pi";
   hFilec+=peds; hFilec+=".hist.root";

  fda=new TFile(hFilea);
  assert(fda->IsOpen());

  fdb=new TFile(hFileb);
  assert(fdb->IsOpen());
  
  fdc=new TFile(hFilec);
  assert(fdc->IsOpen());
  
  int j=1;
  int i=1;
  int k,section;
  int count=1,ind=1;
  float N;
  float centa=0;
  float centb=0;
  float centc=0;
  float meana[800];
  float meanb[800];
  float meanc[800];
  float neta[65];
  float dev[65];
  int div = 12;
  char ch;
  float ratio[800];
  float ratior[800];
  float average[65];
  int index=1;

  for (int d=0; d<12;d++){
    neta[d]=d+1;
  }
   //Loop through all sec values 1 to 12
  for (section=1;section<=12;section++)
      // Loop through all sub values 'A' to 'E'
    for (ch='A';ch<='E';ch++)
	// Loop through all eta values 1 to 12
	for(index=1;index<=12;index++)
	  { 
 
	    int sec=section,eta=index;
	    char sub = ch;
	    char namea[100];
            char nameb[100];
	    
	    sprintf(namea,"%02dT%c%02d",sec,sub,index);
	    sprintf(nameb,"%02dT%c%02d",sec,sub,index);
	    
            
	    TH1F *ha=fda->Get(namea);
	    assert(ha);
  
            TH1F *hb=fdb->Get(nameb);
            assert(hb);
	    
	    TH1F *hc=fdc->Get(namea);
	    assert(hc);

	    centa=ha->GetMaximumBin();
            centb=hb->GetMaximumBin();
	    centc=hc->GetMaximumBin();

	    float *ya=ha->GetArray();
            float *yb=hb->GetArray();
	    float *yc=hc->GetArray();
	  

		
	    meana[j]=meancalc(centa,ya,50);	 
            meanb[j]=meancalc(centb,yb,50);

	    meanc[j]=meancalc(centc,yc,5);

	    // Subtract pedestals
	    meana[j]=meana[j]-meanc[j]; 
	    meanb[j]=meanb[j]-meanc[j];

	    // Flagging bad towers (with just pedestals). If mean a or b is around +/-10, set ratio to 0
	      if (((meana[j]>-10)&&(meana[j]<10))||((meanb[j]>-10)&&(meanb[j]<10)))
		{
		  ratio[j]=0;
		}
	     else
	       ratio[j]=meanb[j]/meana[j];
	    
	    // If ratio is negative or greater than 2 or equal to 1, output tower info and make zero
	    // cutoff = 1.6 for fiber 2
	    int cutoff = 2;
	    //if (runb==1073)
	    //  cutoff=1.6;
	    if (ratio[j] < 0 || ratio[j] > cutoff || ratio[j]==1)
	      {     
		ratio[j]=0;
		//	printf("r=%d tower %s\n",ratio[j],namea);

	      }

	    j++;
	  }

  // Calculate the average ratio for one phi bin (the average ratio for all 12 eta)
  avgcalc(average,ratio);

  // Calculate the "ratio of ratios":ratior=ratio/[avg ratio * average]
  ratiorcalc(average,ratio,ratior);
  
  // Needed to change name in graph according to which run is used
  char graphnm[100];
  int fib=0;

  if (runb==1072)
    fib=1;
  if (runb==1073)
    fib=2;
  if (runb==1074)
    fib=3;
  if (runb==1075)
    fib=4;
  if (runb==1076)
    fib=5;
  if (runb==1078)
    fib=6;
  if (runb==1079)
    fib=7;
  if (runb==1080)
    fib=8;
  
  sprintf(graphnm,"Chi-squared Distribution: Fiber %d",fib);

  TH1F *hx = new TH1F("dev",graphnm,500,0,0.5);
  hx->SetXTitle("Chi-squared");
  hx->SetYTitle("Frequency");

  // Calculate chi squared and store in dev[j]
  div=12;
  j=1;
  count=1;
  float temp[15];
  temp[0]=0;
  temp[13]=0;
  float tp;
  
  for (int i=1; i<=60; i++)
    { 

      dev[i]=0;
      for (int k=1;k<=12;k++)
	{
	  if (ratior[j]!=0)
	    temp[k] = 1.0-ratior[j];
	  else 
	    temp[k]=0;
          
	  
	  tp=temp[k]*temp[k];     

   
	  dev[i]+=tp;  

	  if (ratio[j]==0)
	    {
	      div=div-1; 
	    }
	
	  if (k==12)
	    {
	      dev[i]=dev[i]/div;
	      div=12;
              count=1;
	    }
	 
	  j++;
	}  	
 
      for (int m=0;m<=11;m++)
	{	  
          float avgval=(temp[m]+temp[m+2])/2;
	  if (temp[m]==0)
	    avgval=temp[m+2];
	  if (temp[m+2]==0)
	    avgval=temp[m];
     
	  if (m==0) avgval=temp[2];
	  if (m==11)avgval=temp[12];
	  float diffsq = (avgval-temp[m+1]);
          if (diffsq <0)
	    diffsq=-diffsq;
	  if (diffsq>=0.2||diffsq<=-0.2)
	    //printf("i is %d m is %d diff is %f\n",i,m+1,diffsq);
	}


      // Attempt to pick out non-smooth curves (broken tiles)	    

      // printf("dev %f\n",dev[i]);   
      hx->Fill(dev[i]);
    }

  // Print the phi bins with chi-squared above 0.02
  phiPrint(dev);            
  hx->Draw();
  

   return 0;
}