class TGraphErrors ;

const int mxT=720;
class Tower;
class Tower{
  public:
  char  *name, *pname;
  float phix;
  float sl,ersl;
};

Tower towerA[mxT];


slope2gain(int etaBin=4){
  TGraphErrors *grI=new TGraphErrors; // halds all, for presentation
  TGraphErrors *grIR=new TGraphErrors; // reversed X& Y
  TGraphErrors *grC=new TGraph; // list of corrected gains

  grI->SetMarkerStyle(24);
  grI->SetLineColor(kGreen);
  grC->SetMarkerStyle(3);
  grC->SetLineColor(kRed);
  grIR->SetMarkerStyle(3);
  grIR->SetLineColor(kGreen);

  TString tit="Slopes for eta bin="; tit+=etaBin; tit+="; phi index  ~ TPC sector";
  TH1F* hSl=new TH1F("sl",tit,10,0,14);
  hSl->SetMinimum(-0.2);  hSl->SetMaximum(-0.1);

  tit="weighted, truncated slope distribution; slope";
  TH1F* hAvr=new TH1F("avr",tit,30,-0.21,-0.05);
  
  // grI->SetMarkerColor(kGreen);
  //  grI->SetMarkerSize(1);
  
  int np= readSlopes(etaBin,grI,grIR);
  hSl->SetEntries(np);

  float avrSl=findAverage(grIR,hAvr);
  printf("avrSl=%f\n",avrSl);
  float relGainMarg=0.05;
  doGainCorrection(np,avrSl,relGainMarg,grC,etaBin);


  //
  c=new TCanvas();
  c->Divide(2,1);
  c->cd(1);  hSl->Draw(); grI->Draw("P");  grC->Draw("P");
  c->cd(2);  hAvr->Draw();

  char fname[100];
  sprintf(fname,"etaBin%02d.ps",etaBin);
  //  c->Print(fname);

  //  grI->Print();
  
}

//================================================
//================================================
int
doGainCorrection(int np, float avrSl, float relGainMarg, TGraphErrors *grC,int etaBin) {
  int i;
  char fname[100];
  sprintf(fname,"gainCorr.eta%02d",etaBin);
  FILE *fd=fopen(fname,"w"); assert(fd);
  fprintf(fd,"# HV corrections for etaBin=%d  goalSlope=%.4f\n",etaBin,avrSl);
  fprintf(fd,"#name   pname  gainCorr hvCorr  currSlope  phiIndex\n");
  for(i=0;i<np;i++) {
    Tower *tw=towerA+i; 
    float gainCorr= tw->sl/avrSl;
    int x=0;
    if(fabs(1.-gainCorr)<relGainMarg) { 
      gainCorr= 1.;
    } else {
      int n=grC->GetN();
      grC->SetPoint(n,tw->phix,tw->sl);
    }
    float hvCorr=pow(gainCorr,1./8.3);
    fprintf(fd,"%s %s %.3f %.3f %.4f %f\n",tw->name,tw->pname,gainCorr, hvCorr,tw->sl,tw->phix);
  }
  
  fclose(fd);
}

//================================================
//================================================
float
findAverage( TGraphErrors *grIR, TH1F* hAvr){
  //grIR->Draw("AP");
  int nSkip=3;
   grIR->Sort();
  int i;
  int n=grIR->GetN();
  double *er=grIR->GetEX();
  double *val=grIR->GetX();
  
  for(i=nSkip; i<n-nSkip; i++){
    //  printf("i=%d sl=%f er=%f\n",i,val[i],er[i]);
    float w=1./er[i]/er[i];
    hAvr->Fill(val[i],w);
  }
  hAvr->Fit("gaus");
  TF1 *ff=hAvr->GetFunction("gaus"); assert(ff);
  return ff->GetParameter(1);
}


//================================================
//================================================
int 
readSlopes(int etaBin, TGraphErrors *grI, TGraphErrors *grIR) {

  char *fname="slopes.dat";
  FILE *fd=fopen(fname,"r");
  assert(fd);

  int nt=0; 
  while(1) {
    const int mx=1000;
    char buf[mx];
    char  *pname=new char [mx];
    char  *name1=new char [mx];
    float sl, ersl;
    char *ret=fgets(buf,mx,fd);
    if(ret==0) break;
    if(buf[0]=='#') continue;
    //printf("=%s=\n",buf);
    int n=sscanf(buf,"%s %s %f %f",name1, pname,&sl,&ersl);
    assert(n==4);
    assert(name1[2]=='T');
    int sec=atoi(name1);
    int eta=atoi(name1+4);
    int isub=name1[3]-'A';
    if(eta!=etaBin) continue;
    
    Tower *tw=towerA+nt;
    nt++;
    tw->name=name1;
    tw->pname=pname;
    tw->sl=sl;
    tw->ersl=ersl;
    float phix=sec+isub/5.+(eta-1)/5./12.;
    tw->phix=phix;

    //  printf("name=%s=%s= eta=%d  isub=%d sec=%d sl=%f\n",name1,tw->name, eta,isub,sec, sl);
    //  break;
    //add point to tgraph
    int n=grI->GetN();
    grI->SetPoint(n,phix,sl);
    grI->SetPointError(n,0.,ersl);

    //add point to reversed tgraph
    int n2=grIR->GetN();
    grIR->SetPoint(n2,sl,phix);
    grIR->SetPointError(n2,ersl,0.);
    // if(nt>10) break;
     
    // break;
  }
  printf("found %d towers for eta bin=%d\n",nt,etaBin);
  fclose(fd);

  return nt;
}