//void cere_calib()
{ 
// gas.C   R.Michaels  Dec 2003
// Plot the gas Cerenkov.  Must run analyzer first.
// Do this: ".x gas.C" in root.  Assumes the input 
// is Afile.root (usually this is a link).

// cere_calib.C  X.Zheng  Apr 2004
// This program does the calibration for R-arm gas cherenkov
// and L-arm A1, A2 aerogel cherenkov detectors.
//
// calib.C Y.Qiang  May 2004
// Modification of peak searching algorithm in Gas Cherenkov detector 
// calibration has been made.  
// Some other minor changes.
//
// The user can choose whether the pedastals should be checked, then
//   the program fits to single ph.e. peaks of each channal. 
//   For aerogel detectors a TDC cut is used to clean up the ADC signal;
//   For gas cherenkov, sometimes the fit is not good due to low statistics, 
//       the user can choose either to accept or to decline a fit during 
//       the procedure.
//
// The inputs are 
//   (1) root files containing blocks L.a1.*, L.a2.* and R.gas.*;
//   (2) database files, you need to define the date in 
//       this script to select proper database.
//
// The outputs are
//   (1) new pedastals (set to old values is pedastals were NOT checked);
//   (2) new amplifications for each channel, (set to old values if the 
//       fit failed);
//
//   The outputs are printed to the screen as well as saved in a new database
//   file. The histograms are saved to postscript files.
//
// Make sure the data you are using with correct date

  gROOT->Reset();
  gStyle->SetOptStat(0);
  gStyle->SetLabelSize(0.1,"x");
  gStyle->SetLabelSize(0.1,"y");
  gStyle->SetPadBottomMargin(0.16);
  gStyle->SetPadTopMargin(0.1);
  gStyle->SetPadLeftMargin(0.16);
  gStyle->SetPadRightMargin(0.1);
  gStyle->SetTitleH(0.10);
  gStyle->SetTitleW(0.40);
  
  //  TPostScript ps("./ps/cere_calib.ps",112);  
  
  char filename[100],dummy[200];
  TDatime date("2004-06-01 14:15:00"); //date for data we are using
 
  Int_t nrun, nrunped, ped;

  Int_t nchevt = 500000;
  Int_t naevt = 50000;

  cout << "Enter the run number for pedastal run, or -1 if do not check pedastal" << endl;
  cin >> nrunped;
  if (nrunped<0) ped=0;
  else ped=1;

  TTree *t1=0, *t2=0;
  if(ped){
    sprintf(filename,"./rootfiles/ped_%d.root",nrunped);  // pedistal data file
    printf("%s openned\n",filename);
    f1 = new TFile(filename);
    t1 =  (TTree *)f1->Get("T"); 
  }

  cout << "Enter run number to use" << endl;
  cin >> nrun;
  if (nrun<0) break;
  sprintf(filename,"./rootfiles/cere_%d.root",nrun);
  printf("%s openned\n",filename);
  f2 = new TFile(filename);
  t2 =  (TTree *)f2->Get("T");   

  Int_t detid, i, j, chtot;
  Float_t ped_old[50], ped_new[50];
  Float_t amp_old[50], amp_new[50];
  Float_t peak_new[50], fit_okay[50];
  Float_t peak_mult[50], peak_mult_err[50];
  Float_t hv_corr[50];
  Float_t tdcp0,tdcp1,tdcp2;

  // align the single ph.e. peak to this channel 100 for A1, 60 for A2, 100 for GasC
  Float_t peak_pos[3]={100.,60.,100.}; 

  Int_t hxmin,hxmax, imin, imax;
  FILE *dbfile;
  char varname[100],cutname[100];
  cout << "Which detector do you want to check? (1=Aero1, 2=Aero2, 3=Gas, 0=Quit)" <<endl;
  cin >> detid;

  if ( detid<=0 || detid>3 ) break;  


  //////////////////////////////////////////////////////////////
  /// Gas Cherenkov Calibration
  //////////////////////////////////////////////////////////////

  if (detid==3){   
    cout << "Check Gas Cerenkov Calibration:\n" << endl;
    // 111=portrait, 112=landscape, 113=eps

    FILE* dbfile = THaAnalysisObject::OpenFile("R.cer",date,"","r",1);
    char buf[256];
    for (j=0;j<2;j++) fgets(buf,256,dbfile); // gets one line from file
    fscanf(dbfile,"%s",dummy);
    chtot=atoi(dummy);
    printf("total %d channels\n\n",chtot);
    for (j=0;j<14;j++) {
      fgets(buf,256,dbfile); 
      //printf("r0 %s\n",buf);
    }
    for (j=0;j<chtot;j++){
      fscanf(dbfile,"%s",dummy);
      ped_old[j]=atof(dummy);
    }
    for (j=0;j<2;j++) fgets(buf,256,dbfile);
    for (j=0;j<chtot;j++){
      fscanf(dbfile,"%s",dummy);
      amp_old[j]=atof(dummy);
      //    printf("%d %d %.3f\n",j,ped_old[j], amp_old[j]);
    }

    int i=0;
    Double_t p0,p1,p2;
    TCanvas *c1 =  new TCanvas("c1","Gas Cerenkov ADC",200,10,800,600);
    c1->SetFillColor(20);
    if (ped) c1->Divide(2,1);
    TPad *c1_1 = (TPad*)c1->GetPad(1);
    TH1F* h1 = new TH1F("h1","",100,400,600);
    TH1F* h2 = new TH1F("h2","",50,500,1000);

    while (i<chtot){
      h1->Reset();
      h1->SetBins(100,ped_old[i]-50,ped_old[i]+50);
      sprintf(cutname,"ADC %d, pedastal",i+1);
      h1->SetTitle(cutname);
      TAxis *axisx=h1->GetXaxis();
      axisx->SetLabelSize(0.05);
      TAxis *axisy=h1->GetYaxis();
      axisy->SetLabelSize(0.05);
      
      h2->Reset();
      h2->SetBins(50,ped_old[i]+20,ped_old[i]+420);
      sprintf(cutname,"ADC %d, single ph.e. peak",i+1);
      h2->SetTitle(cutname);
      axisx=h2->GetXaxis();
      axisx->SetLabelSize(0.05);
      axisy=h2->GetYaxis();
      axisy->SetLabelSize(0.05);  

      ped_new[i]=ped_old[i];
      // pedastals

      if (ped) {
	c1->cd(1);
	printf("now look for pedastal for ADC ch %d, from %d to %d\n",i+1, ped_old[i]-50, ped_old[i]+50);
	sprintf(varname,"R.cer.a[%d]>>h1",i);
	sprintf(cutname,"R.cer.a[%d]>%d && R.cer.a[%d]<%d",i,ped_old[i]-50,i,ped_old[i]+50);
	t1->Draw(varname,cutname,"",5000);
      
	imin=ped_old[i]-10; 
	imax=ped_old[i]+10;
	h1->Fit("gaus","","",imin,imax);
    
	TF1 *fit1 = h1->GetFunction("gaus");
	p0 = fit1->GetParameter(0);
	p1 = fit1->GetParameter(1);
	ped_new[i]=p1;

	TLine *a0 = new TLine(ped_new[i],0.,ped_new[i],1.0*(p0));
	a0->SetLineColor(3);
	a0->Draw();

	TLine *a1 = new TLine(ped_old[i],0.,ped_old[i],0.9*(p0));
	a1->SetLineColor(2);
	a1->Draw();
      }

      // end of pedestals
 
      c1->cd(1+ped);

      sprintf(varname,"R.cer.a[%d]>>h2",i);
      sprintf(cutname,"R.cer.t[%d]>1400&&R.cer.t[%d]<1500",i,i);
      //      sprintf(cutname,"(R.ps.asum_c+0.8*R.sh.asum_c)>1500 && R.ps.asum_c>400",i,i);
      t2->Draw(varname,cutname,"",nchevt);
      imin=ped_old[i]+(100./amp_old[i]*(1.-0.6));
      if (imin<ped_old[i]+20) imin=ped_old[i]+20;
      imax=ped_old[i]+(100./amp_old[i]*(1.+0.6));
      if (imax>ped_old[i]+520) imax=ped_old[i]+520;
      printf("first fit ch %d from %d to %d\n", i+1, imin, imax);
      h2->Fit("gaus","","",imin,imax);
      TF1 *fit2 = h2->GetFunction("gaus");
      Double_t p0 = fit2->GetParameter(0);
      Double_t p1 = fit2->GetParameter(1);
      Double_t p2 = fit2->GetParameter(2);
      printf("Coarse fit result: peak: %.2f, sigma: %.2f\n",p1, p2);
      
      for (int ii = 0; ii <= 10; ii++){
	imin = p1-p2*1.5;
	if (imin<ped_old[i]+10) imin=ped_old[i]+10;
	imax = p1+p2*0.8;
	if (imax>ped_old[i]+510) imax=ped_old[i]+510;
	printf("No.%d fine fit ch %d from %d to %d\n", ii+1, i+1, imin, imax);
	h2->Fit("gaus","","",imin,imax);
	fit2 = h2->GetFunction("gaus");
	p0 = fit2->GetParameter(0);
	p1 = fit2->GetParameter(1);
	p2 = fit2->GetParameter(2);
	printf("No.%d fit result: peak: %.2f, sigma: %.2f\n", ii+1, p1, p2);
      }
      
      // plot new peak position
      amp_new[i]=100./(p1-ped_new[i]);
      peak_new[i]=p1;
      
      TLine *a0 = new TLine(p1,0.,p1,1.0*(p0));
      a0->SetLineColor(3);
      a0->Draw();

      TLine *a1 = new TLine(ped_old[i]+(100./amp_old[i]),0.,ped_old[i]+(100./amp_old[i]),0.9*(p0));
      a1->SetLineColor(2);
      a1->Draw();
      c1->Update();
      
      if (ped) sprintf(varname,"./ps/R_gas_wped_%d.ps",i+1);
      else sprintf(varname,"./ps/R_gas_%d.ps",i+1);
      c1->Draw();
      c1->Print(varname);
      printf("%s created\n", varname);
      
      cout << "Enter -1 to exit, 0 to deny the fit, or any character to continue" <<endl;
      cin >> dummy;
      if (atoi(dummy)==-1) break;
      else if (atoi(dummy)==0){fit_okay[i]=0; amp_new[i]=amp_old[i];}
      else fit_okay[i]=1;
      i++;
    }
    //    delete h1,h2;
  
    for (i=0;i<chtot;i++){
      if (fit_okay[i]==1){
	if (amp_new[i]) 
	  printf("channel %d: ped_o=%d ped_n=%d amp_o=%.3f amp_n=%.3f\n",
		 i+1, ped_old[i], ped_new[i], amp_old[i], amp_new[i]);
      }
      else printf("channel %d: ped_o=%d ped_n=%d amp_o=%.3f amp_new failed\n",
		  i+1, ped_old[i], ped_new[i], amp_old[i]);
    }


    sprintf(filename,"db_R.cer.dat.new");
    FILE* newdbfile = fopen(filename,"w+");
    {
      fprintf(newdbfile, "Pedestal values of Gas Cherenkov channels --------------------------------------\n");
      for(i=0;i<10;i++)    fprintf(newdbfile, "%d   ",ped_new[i]);
      fprintf(newdbfile, "\nAmplitude transformation coefficients of Gas Cherenkov channels ------------------\n");
      for(i=0;i<10;i++)     fprintf(newdbfile, "%.3f  ",amp_new[i]);
      fprintf(newdbfile,"\n################################################################################\n");
    }
    fclose(newdbfile);
    

    printf("Replace the constants in file <db_R.cer.dat> by the following:\n\n");
    printf("Pedestal values of Gas Cherenkov channels --------------------------------------\n");
    for(i=0;i<10;i++)  printf("%d   ",ped_new[i]);
    cout << endl;
    printf("\nAmplitude transformation coefficients of Gas Cherenkov channels ------------------\n");
    for(i=0;i<10;i++) printf("%.3f  ",amp_new[i]);
    cout << endl;
    printf("################################################################################\n");
    cout << "also saved in db_R.cer.dat.new"<<endl;
    
  }

  ///////////////////////////////////////////////////////////////////////
  /// AeroGel Detectors
  ///////////////////////////////////////////////////////////////////////

if (detid==1 || detid==2){
  gStyle->SetOptLogy(0);

  printf("Check Aero %d Calibration:\n\n", detid);

  sprintf(varname,"L.a%d",detid);
  FILE* dbfile = THaAnalysisObject::OpenFile(varname,date,"","r",1);
  //FILE* dbfile = fopen(Form("DB/%s/db_L.a%d.dat",DBdate,detid),"r");
  //cout<<"DB/"<<DBdate<<"/db_L.a"<<detid<<".dat opened"<<endl;
  // now read the file
  char buf[256];
  for (j=0;j<2;j++) fgets(buf,256,dbfile); // gets one line from file
  //  printf("%s",buf);
    // keep reading until you get the pedestals
    //...
  // dbfile = fopen(Form("DB/20040413/db_L.%d.dat",detid),"r");
    //  if (dbfile != NULL){
    //  for (j=0;j<13;j++) {fscanf(dbfile,"%s",dummy);}
  fscanf(dbfile,"%s",dummy);
  chtot=atoi(dummy);
  printf("%d channels\n\n",chtot);
  for (j=0;j<17+detid;j++) fgets(buf,256,dbfile);
  for (j=0;j<chtot;j++){
    fscanf(dbfile,"%s",dummy);
    ped_old[j]=atof(dummy);
  }
  for (j=0;j<2;j++) fgets(buf,256,dbfile);
  for (j=0;j<chtot;j++){
    fscanf(dbfile,"%s",dummy);
    amp_old[j]=atof(dummy);
    //    printf("%d %d %.3f\n",j,ped_old[j], amp_old[j]);
  }

  i=0;

  if (ped){
    c1 = new TCanvas("c1","",200,10,900,600);
    c1->Divide(3,4);
  }
  else{
    c1 = new TCanvas("c1","",200,10,800,600);
    c1->Divide(2,4);
  }

  while (i<chtot){    
    
    if ((int(i/4)*4+4)>chtot)
      sprintf(varname,"Aero%d Cerenkov ADC %d to %d",
	      detid, int(i/4)*4+1,chtot);
    else sprintf(varname,"Aero%d Cerenkov ADC %d to %d",
		 detid, int(i/4)*4+1, int(i/4)*4+4);
    c1->SetTitle(varname);
    c1->SetFillColor(20);

    TH1F* hist[12];
    TH1F* hist_plot[12];

    int ih = 0;
    for (j=0;j<4;j++){
      
      hist[ih] = new TH1F(Form("h%d",ih), cutname, 600, 1200,1800);
      TAxis *axisx=hist[ih]->GetXaxis();
      axisx->SetLabelSize(0.1);
      TAxis *axisy=hist[ih]->GetYaxis();
      axisy->SetLabelSize(0.1);
      hist[ih]->SetMinimum(1.);
      hist[ih]->SetFillColor(0);
      ih++;
      
      hist[ih] = new TH1F(Form("h%d",ih), cutname, 250, ped_old[i]-100,ped_old[i]+400);
      TAxis *axisx=hist[ih]->GetXaxis();
      axisx->SetLabelSize(0.1);
      TAxis *axisy=hist[ih]->GetYaxis();
      axisy->SetLabelSize(0.1);
      hist[ih]->SetMinimum(1.);
      hist[ih]->SetFillColor(0);
      ih++;

      hist[ih] = new TH1F(Form("h%d",ih), cutname, 40, ped_old[i]-20, ped_old[i]+20);
      TAxis *axisx=hist[ih]->GetXaxis();
      axisx->SetLabelSize(0.1);
      TAxis *axisy=hist[ih]->GetYaxis();
      axisy->SetLabelSize(0.1);
      hist[ih]->SetMinimum(1.);
      ih++;
      
      if (i<chtot){
	
	// TDC spectrum
	
	c1->cd(j*(ped+2)+1);
	
	//	sprintf(varname, "h%d", 3*j );
	sprintf(cutname, "Aero %d TDC %d, raw", detid, i+1);
	hist[3*j]->SetTitle(cutname);
	hist[3*j]->SetLineColor(1);
	
	sprintf(varname,"L.a%d.t[%d]>>h%d",detid,i,3*j);
	sprintf(cutname,"L.a%d.t[%d]>1200 && L.a%d.t[%d]<1800",detid,i,detid,i);
	
	t2->Draw(varname,cutname,"",naevt);
	
	TH1F *hist_plot[3*j+0] = (TH1F*)hist[3*j+0]->Clone("h0new");
	hist_plot[3*j+0]->Draw();

	imin=1400; imax=1600;
	hist[3*j]->SetLineColor(1);
	hist[3*j]->Fit("gaus","QR0","",imin,imax);
	
	// plot the fit result in green
	TF1 *fit1=hist[3*j]->GetFunction("gaus");
	Double_t p0 = fit1->GetParameter(0);
	Double_t p1 = fit1->GetParameter(1);
	Double_t p2 = fit1->GetParameter(2);
	
	tdcp0=p0; tdcp1=p1; tdcp2=p2;
	
	imin=tdcp1-tdcp2*5.; imax=tdcp1+tdcp2*5.;
	printf("TDC clean region: ch %d to %d\n", imin, imax);
	TF1 *h1fit=new TF1("h1fit","gaus",imin,imax);
	h1fit->SetParameters(p0,p1,p2);
	h1fit->SetLineColor(3);
	h1fit->SetLineWidth(2);
	h1fit->Draw("same");
	
	// ADC spectrum with TDC cut
	
	c1->cd(j*(ped+2)+2);
	
	//	sprintf(varname, "h%d", 3*j+1);
	sprintf(cutname, "Aero %d ADC %d, w/ TDC cut", detid, i+1);
	hist[3*j+1]->SetTitle(cutname);
	hist[3*j+1]->SetLineColor(1);
	
	sprintf(varname,"L.a%d.a[%d]>>h%d",detid,i,3*j+1);
	sprintf(cutname,"L.a%d.a[%d]>10 && L.a%d.a[%d]<1200 && L.a%d.t[%d]>%d && L.a%d.t[%d]<%d",detid,i,detid,i,detid,i,imin,detid,i,imax);
	
	t2->Draw(varname,cutname,"",naevt);
	
	TH1F *hist_plot[3*j+1] = (TH1F*)hist[3*j+1]->Clone("h0new");
	hist_plot[3*j+1]->Draw();

	if (detid==1){
	  imin=ped_old[i]+(peak_pos[detid-1]/amp_old[i]*(1.-0.5));
	  imax=ped_old[i]+(peak_pos[detid-1]/amp_old[i]*(1.+0.5));
	}
	else{
	  // for a2, fit wider range
	  imin=ped_old[i]+(peak_pos[detid-1]/amp_old[i]*(1.-0.6));
	  imax=ped_old[i]+(peak_pos[detid-1]/amp_old[i]*(1.+0.6));
	}
	printf("fit ch %d from %d to %d\n",i, imin, imax);
	hist[3*j+1]->SetLineColor(1);
	hist[3*j+1]->Fit("gaus","QR0","",imin,imax);
	
	// plot the fit result in green
	TF1 *fit2=hist[3*j+1]->GetFunction("gaus");
	Double_t p0 = fit2->GetParameter(0);
	Double_t err0 = fit2->GetParError(0);
	Double_t p1 = fit2->GetParameter(1);
	Double_t err1 = fit2->GetParError(1);
	Double_t p2 = fit2->GetParameter(2);
	Double_t err2 = fit2->GetParError(2);
	//     break;
	
	TF1 *h2fit=new TF1("h2fit","gaus",imin,imax);
	h2fit->SetParameters(p0,p1,p2);
	h2fit->SetLineColor(3);
	h2fit->SetLineWidth(2);
	h2fit->Draw("same");
	
	
	// pedastals
	
	if (ped){
	  
	  c1->cd(j*(ped+2)+3);
	  sprintf(cutname,"ADC %d, pedastal",i+1);
	  hist[3*j+2]->SetTitle(cutname);
	  
	  sprintf(varname,"L.a%d.a[%d]>>h%d",detid,i,3*j+2);
	  sprintf(cutname,"L.a%d.a[%d]>10 && L.a%d.a[%d]<1000",detid,i,detid,i);
	  t1->Draw(varname,cutname,"",naevt);
	  
	  TH1F *hist_plot[3*j+2] = (TH1F*)hist[3*j+2]->Clone("h0new");
	  hist_plot[3*j+2]->Draw();

	  imin=ped_old[i]-10; imax=ped_old[i]+10;
	  hist[3*j+2]->Fit("gaus","RQ0","",imin,imax);
	  
	  TF1 *fit3 = hist[3*j+2]->GetFunction("gaus");
	  Double_t p0_single=p0;
	  Double_t p0 = fit3->GetParameter(0);
	  Double_t err0 = fit3->GetParError(0);
	  Double_t p1 = fit3->GetParameter(1);
	  Double_t err1 = fit3->GetParError(1);
	  Double_t p2 = fit3->GetParameter(2);
	  Double_t err2 = fit3->GetParError(2);
	  
	  // plot the fit result in green
	  TF1 *h3fit=new TF1("h3fit","gaus",imin,imax);
	  h3fit->SetParameters(p0,p1,p2);
	  h3fit->SetLineColor(3);
	  h3fit->SetLineWidth(1);
	  h3fit->Draw("same");
	  
	  
	  ped_new[i]=p1;
	  
	  hist[3*j+2]->SetLineColor(1);
	  hist[3*j+2]->SetAxisRange(p1-20,p1+20);
	  
	  // new pedastal in GREEN
	  TLine *a0 = new TLine(ped_new[i],0.,
				ped_new[i],0.8*p0);
	  a0->SetLineColor(3);
	  a0->SetLineWidth(4);
	  a0->Draw();
	  
	  // old pedastal in RED
	  TLine *a1 = new TLine(ped_old[i],0.,
				ped_old[i],1.1*p0);
	  a1->SetLineColor(2);
	  a1->SetLineWidth(2);
	  a1->Draw();
	}
	else{
	  ped_new[i]=ped_old[i];
	}

	// end of pedestals
	
	c1->cd(j*(ped+2)+1);
	Double_t p0 = fit2->GetParameter(0);
	Double_t err0 = fit2->GetParError(0);
	Double_t p1 = fit2->GetParameter(1);
	Double_t err1 = fit2->GetParError(1);
	Double_t p2 = fit2->GetParameter(2);
	Double_t err2 = fit2->GetParError(2);
	hist[3*j+1]->SetMaximum(p0*1.4);
	
	// new pedastal in GREEN
	TLine *a0 = new TLine(ped_new[i],1.,ped_new[i],p0*0.8);
	a0->SetLineColor(4);
	a0->Draw("same");
	
	// old pedastal in RED
	TLine *a1 = new TLine(ped_old[i],1.,ped_old[i],p0);
	a1->SetLineColor(2);
	a1->Draw("same");
	
	// old peak position in RED
	TLine *a2 = new TLine(ped_old[i]+peak_pos[detid-1]/amp_old[i],1.,
			      ped_old[i]+peak_pos[detid-1]/amp_old[i],p0);
	a2->SetLineColor(2);
	a2->Draw("same");
	
	// new peak position in GREEN
	amp_new[i]=peak_pos[detid-1]/(p1-ped_new[i]);
	Float_t vx[4]={p1-err1/2.,p1-err1/2.,p1+err1/2.,p1+err1/2.};
	Float_t vy[4]={1.,p0/2.,p0/2.,1.};
	TPolyLine *a3 = new TPolyLine(4,vx,vy);
	a3->SetFillColor(3);
	a3->SetLineColor(3);
	a3->SetLineWidth(1);
	a3->Draw("same");
	
	i++;
      }
      
      c1->Update();
    }
    //    c1->Draw();
    sprintf(varname,"./ps/L_a%d_%d.ps",detid,int(i/4));
    c1->Print(varname);
    
    if (i<chtot){ 
      cout << "red: old calibration, green: new fit" <<endl;
      cout << "Enter 0 to exit, or any char to continue" <<endl;
      cin >> dummy;
      if (atoi(dummy)==0){
	for(j=0; j<12; j++){
	  delete hist[j];
	}
	break;
      }
    }
    for(j=0; j<12; j++){
      delete hist[j];
    }
    
  }//end while(i<chtot)
  
  cout << "Found >=3sigma shifts in calibration constants for following channels:\n"<<endl;
  
  for(i=0;i<chtot;i++){
      printf("channel \#%d: ped_o = %.0f ped_n = %.0f amp_o=%.3f amp_n=%.3f\n",
	     i, ped_old[i], ped_new[i],
	     amp_old[i],amp_new[i]);
  } 

  sprintf(filename,"db_L.a%d.dat.new",detid);
  FILE* newdbfile = fopen(filename,"w+");
  {

    fprintf(newdbfile, "New pedastal values from run %d and new gain from run %d\n\n", nrunped, nrun);
    fprintf(newdbfile, "Pedestal values of Aerogel Cherenkov channels -----------------------------------\n");
    for(i=0;i<10;i++)    fprintf(newdbfile, "%d   ",ped_new[i]);
    fprintf(newdbfile,"\n");
    for(i=10;i<20;i++)    fprintf(newdbfile, "%d   ",ped_new[i]);
    fprintf(newdbfile,"\n");
    for(i=20;i<chtot;i++)    fprintf(newdbfile, "%d   ",ped_new[i]);
    fprintf(newdbfile,"\n");
    
    fprintf(newdbfile, "Amplitude transformation coefficients of Aerogel ADC channels ------------------\n");
    for(i=0;i<10;i++)     fprintf(newdbfile, "%.3f  ",amp_new[i]);
    fprintf(newdbfile,"\n");
    for(i=10;i<20;i++)     fprintf(newdbfile, "%.3f  ",amp_new[i]);
    fprintf(newdbfile,"\n");
    for(i=20;i<chtot;i++)     fprintf(newdbfile, "%.3f  ",amp_new[i]);
    fprintf(newdbfile,"\n");
    fprintf(newdbfile,"################################################################################\n");
  }
  fclose(newdbfile);

  printf("Replace the constants in file <db_L.a%d.dat> by the following:\n\n", detid);
  printf("Pedestal values of Aerogel Cherenkov channels -----------------------------------\n");
  for(i=0;i<10;i++){
    printf("%d   ",ped_old[i]);
  }
  cout << endl;
  for(i=10;i<20;i++){
    printf("%d   ",ped_old[i]);
  }
  cout << endl;
  for(i=20;i<chtot;i++){
    printf("%d   ",ped_old[i]);
  }
  printf("\nAmplitude transformation coefficients of Aerogel ADC channels ------------------\n");
  for(i=0;i<10;i++)     printf("%.3f  ",amp_new[i]);
  cout << endl;
  for(i=10;i<20;i++)     printf("%.3f  ",amp_new[i]);
  cout << endl;
  for(i=20;i<chtot;i++)     printf("%.3f  ",amp_new[i]);
  cout << endl;
  printf("################################################################################\n");
  cout << "also saved in "<<filename<<endl;
   
  printf(" \n\nCorrection coeff for HV :\n");
  printf("\nNote: this is a suggestion only,  approval from Hall A exp coordinator and spokesperson needed for HV adjustment\n\n");

  for(int i=0; i<chtot ;i++) {
    hv_corr[i] = 100./exp(log(1./amp_new[i])/15.);
  }
  for(i=0;i<10;i++)     printf("%.3f  ",hv_corr[i]);
  cout << endl;
  for(i=10;i<20;i++)     printf("%.3f  ",hv_corr[i]);
  cout << endl;
  for(i=20;i<chtot;i++)     printf("%.3f  ",hv_corr[i]);
  cout << endl;


  sprintf(filename,"db_L.a%d.hv.new",detid);
  FILE* hvfile = fopen(filename,"w+");
  {
    fprintf(hvfile, "Suggested HV adjustment (%) from pedastal run %d and run %d\n\n", nrunped, nrun);
    fprintf(hvfile, "Note: approval from Hall A exp coordinator and spokesperson needed for HV adjustment\n\n", nrunped, nrun);
    fprintf(hvfile, " Aerogel Cherenkov channels -----------------------------------\n");
    for(i=0;i<10;i++)    fprintf(hvfile, "%8.2f   ",hv_corr[i]);
    fprintf(hvfile,"\n");
    for(i=10;i<20;i++)    fprintf(hvfile, "%8.2f   ",hv_corr[i]);
    fprintf(hvfile,"\n");
    for(i=20;i<chtot;i++)    fprintf(hvfile, "%8.2f   ",hv_corr[i]);
    fprintf(hvfile,"\n");
    
    fprintf(hvfile,"################################################################################\n");
  }
  fclose(hvfile);

  cout << "also saved in "<<filename<<endl;




}
break;

}