TH1F *hx[4];

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

// Creates histogram of centroid frequency with 10 values in a bin
int initHisto(){
  int nblaser = 400; // Number bins for laser run
  int nbpeds = 4000; //Number bins for pedestal run

  hx[0]=new TH1F("centr","Mean Distribution-peds subtracted, run=1066",nblaser,0,4000);
    hx[1]=new TH1F("mn","Mean Distribution-peds subtracted, run 1060",nbpeds,0,4000);
}

//====================================
//====================================
//====================================
int finish(){
  printf("Finish, nCentroids=%d\n",hx[0]->GetEntries());
  hx[0]->SetXTitle("Mean");
  hx[0]->Draw(); 
  }
//====================================
//====================================
//====================================
int finhxone(){
  printf("Finish hx 1  nCentroids=%d\n",hx[1]->GetEntries());
  hx[1]->SetXTitle("Mean");
  hx[1]->Draw();
}

//====================================
//====================================
//====================================
#if 1
float meancalc(int cent,float yval[],int offset){
  float meanval=0;
  float N=0;
  int low=120;
  
  if (centSetOptStat(1111111);
  initHisto();

  TString hFilea="/N/u1/bullinger/outLaser/pi";
  TString hFileb="/N/u1/bullinger/outLaser/pi";
  hFilea+=runa;  hFilea+=".hist.root";
  hFileb+=runb;  hFileb+=".hist.root";
  
  fda=new TFile(hFilea);
  assert(fda->IsOpen());

  fdb=new TFile(hFileb);
  assert(fdb->IsOpen());

 

  int j,i,k;
  char ch;
  int section;
  int offlaser = 50;  // For run 1066: range around MaxBin to include in mean calculation 
  int high = 3000; // High limit for sending error message 
  int low = 60;   // Low limit for sending error message
  float N;  
  
  int cenlaser[800];
  int cenpeds[800];
  float meanlaser[800];
  float meanpeds[800];
  float sigma[800];
  float subpeds[800];
  
  
  int index;   
  
  j=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 name[100];

	    sprintf(name,"%02dT%c%02d",sec,sub,index);
	    //	    printf("searching for histo '%s'\n",name);

	    TH1F *ha=fda->Get(name);
	    assert(ha);
  
            TH1F *hb=fdb->Get(name);
    
	    float *y=ha->GetArray();
            float *ya=ha->GetArray();
  
	    cenlaser[j]=ha->GetMaximumBin();
            cenpeds[j]=hb->GetMaximumBin();

	    meanlaser[j]=0;
	    meanpeds[j]=0;
	    
	    sigma[j]=0;

            // Calculates mean
	    meanlaser[j]=meancalc(cenlaser[j],y,offlaser);
            meanpeds[j]=meancalc(cenpeds[j],ya,5);
	    // Calculates RMS or sigma
	    sigma[j]=sigmacalc(cenlaser[j],y,offlaser,meanlaser[j]);

	    float centr=meanlaser[j]-meanpeds[j]+50;
            float mn=meanpeds[j];
	    hx[0]->Fill(centr);
            hx[1]->Fill(mn);


	    // if (meanpeds[j]>high)
	    //  printf("High Mean '%s', MaxBin=%d, mean=%f\n",name,centroid[j],mean[j]);
	    //            if (meanpeds[j]< 5)
	    //	      printf("Low Mean '%s', MaxBin=%d, mean=%f\n",name,cenpeds[j],meanpeds[j]);
#if 1
            if ((meanlaser[j]-meanpeds[j]+50)>high)
	      printf("High Mean '%s', MaxBin=%d, mean=%f\n",name,cenlaser[j],meanlaser[j]-meanpeds[j]+50);
            if ((meanlaser[j]-meanpeds[j]+50)< low)
	      printf("Low Mean '%s', MaxBin=%d, mean=%f\n",name,cenlaser[j],meanlaser[j]-meanpeds[j]+50);
#endif
	    //printf("centroid= %d,mean= %f,sigma= %f\n",centroid[j],mean[j],sigma[j]);

	    j++;
	  }
    
  
  finish();
 
  return;


}