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; }