// - can be given a set of run numbers (which need not necessarily exist)
// - Determines [and plots] tdc offsets (after subtraction of walk effect) for each of the L20 PMTs
// - [ TDC_raw - ( TDC_param(p0,p1) + Offset ) = 1000,  where  
//   TDC_param(p0,p1) = p0/sqrt(ADC-ADC_pedestal-p1) ]
// - ADC pedestals and TDC walk effect parameterization are read from datafiles
// - Reads in Walk/L20-walk.dat for walkparameters and Pedestals/L20-adc-ped-####.dat for pedestals
// - Writes to L20_cal-####.dat (format like blast.sc_cal)
//   and to L20-tdc-after-walk-####.dat (walk-corrected) and L20-tdc-####.dat (uncorrected)
// - Uses gRaw->GetL20 in TBLRaw.cc with cut on flasher events

// Revisions:
// - Original code updated (4/9/04)

{
  gROOT->Reset();

  // read in the data file numbers, put them in list[]
  char tmp[250];
  int i, j, k, n1, n2, list[2000], nl, il;
  int arc=gApplication->Argc(); 
  char** arv=gApplication->Argv(); 

  for (i=0;i<arc;i++){ // scan command line

    TString arg(arv[i]);
    k = sscanf(arv[i],"%d-%d",&n1,&n2);

    if (k == 1 && n1 > 0) list[nl++] = n1;
    if (k == 2 && n1 > 0 && n2>n1)
      for (int j=n1;j<n2+1;j++) {
	if (j<301) { cout<<"Do not use this script for run numbers < 301"<<endl; exit(0);}
	list[nl++] = j;
      }
  }

  // ********** read in walk effect parameters **********
  
  char data[100];

  float walk_par0[2][2][14];  //  left=0,right=1 / top=0,bottom=1 / bars 0-13
  float walk_par1[2][2][14];  //  left=0,right=1 / top=0,bottom=1 / bars 0-13

  char walkfilename[100];
  sprintf(walkfilename,"../Walk/L20-walk.dat");

  FILE *test;
  test = fopen(walkfilename,"rb");
  if (test==NULL) 
  {
    cout<<endl;
    cout<<"Error: Walk effect parameter file "<<walkfilename<<" not found."<<endl;
    cout<<endl;
//    exit(0);
  }
  else
  {
    fclose(test);
      
    cout<<"Processing walk effect parameter file "<<walkfilename<<" ..."<<endl;
    ifstream walkfile(walkfilename, ios::in);
      
    int lr,tb,bar;
    float p0,p1,p2;
    while(!walkfile.eof())
    {
      walkfile.getline(data,100);
      if (data[0]!='#' && data[0]!=NULL)  // filter comment lines (starting with "#")
      {
        // cout<<"'"<<data<<"'"<<endl;
        sscanf(data,"%d\t%d\t%d\t%f\t%f\t%f",&lr,&tb,&bar,&p0,&p1,&p2);
        walk_par0[lr][tb][bar]=p0;
        walk_par1[lr][tb][bar]=p1;
        // cout<<"lr: "<<lr<<"  / tb: "<<tb<<"  / bar: "<<bar<<"  / p0: "<<p0<<"  / p1: "<<p1<<"  / p2: "<<p2<<endl;
      }  //  end of comment line filter
    }  //  end of "eof check" while loop
    walkfile.close();

    gSystem.Load("./libBlast.so");
    Setup();

    for (int il=0;il<gOpt->nList();il++) {

      // ********** read in ADC pedestals **********

      int ped[2][2][16];  //  left=0,right=1 / top=0,bottom=1 / bars 0-15
      int ped_pmt[28][2];
      double par_pmt[28][2][3];

      char pedfilename[100];
      sprintf(pedfilename,"../Pedestals/data/L20-adc-ped-%d.dat",list[il]);

      FILE *test;
      test = fopen(pedfilename,"rb");
      if (test==NULL) 
      {
        cout<<endl;
        cout<<"Error: Pedestal file "<<pedfilename<<" not found."<<endl;
        cout<<endl;
//        exit(0);
      }
      else
      {
        fclose(test);
      
        cout<<"Processing ADC pedestal file "<<pedfilename<<" ..."<<endl;
        ifstream pedfile(pedfilename, ios::in);
      
        int lr,tb,bar,p;
        while(!pedfile.eof())
        { 
          pedfile.getline(data,100);
          if (data[0]!=33 && data[0]!=NULL)  // filter comment lines (starting with "!")
          {
            sscanf(data,"%d\t%d\t%d\t%d",&lr,&tb,&bar,&p);
            ped[lr][tb][bar]=p;
            // cout<<list[i]<<": "<<lr<<" / "<<tb<<" / "<<bar<<" = "<<ped[i][lr][tb][bar]<<endl;
          }  //  end of comment line filter
        }  //  end of "eof check" while loop
        pedfile.close();
 
        char title[125];
        int k;
        int sect, pad, tb;

        // make arrays of pointers for all the histograms.
        TH1F *tdcLEFT_top[14], *tdcLEFT_bot[14] ;
        TH1F *tdcRIGHT_top[14], *tdcRIGHT_bot[14] ;

        TH1F *tdcrawLEFT_top[14], *tdcrawLEFT_bot[14] ;
        TH1F *tdcrawRIGHT_top[14], *tdcrawRIGHT_bot[14] ;

        // book the histograms
        for(k=0;k<14;k++) {
          sprintf(title,"L20[%d]_LEFT_top_tdc",k);
          TH1F *tdcLEFT_top[k]  = new TH1F(title,title,4097,-0.5,4096.5);
          sprintf(title,"L20[%d]_LEFT_bot_tdc",k);
          TH1F *tdcLEFT_bot[k]  = new TH1F(title,title,4097,-0.5,4096.5);
          sprintf(title,"L20[%d]_RIGHT_top_tdc",k);
          TH1F *tdcRIGHT_top[k] = new TH1F(title,title,4097,-0.5,4096.5);
          sprintf(title,"L20[%d]_RIGHT_bot_tdc",k);
          TH1F *tdcRIGHT_bot[k] = new TH1F(title,title,4097,-0.5,4096.5);

          sprintf(title,"L20[%d]_LEFT_top_tdc_raw",k);
          TH1F *tdcrawLEFT_top[k]  = new TH1F(title,title,4097,-0.5,4096.5);
          sprintf(title,"L20[%d]_LEFT_bot_tdc_raw",k);
          TH1F *tdcrawLEFT_bot[k]  = new TH1F(title,title,4097,-0.5,4096.5);
          sprintf(title,"L20[%d]_RIGHT_top_tdc_raw",k);
          TH1F *tdcrawRIGHT_top[k] = new TH1F(title,title,4097,-0.5,4096.5);
          sprintf(title,"L20[%d]_RIGHT_bot_tdc_raw",k);
          TH1F *tdcrawRIGHT_bot[k] = new TH1F(title,title,4097,-0.5,4096.5);
        } 

	FILE *test2;
	test2 = fopen(gOpt->DataFile(il),"rb");
	if (test2==NULL) {
	  cout<<endl;
	  cout<<"Error: Data file "<<gOpt->DataFile(il)<<" not found."<<endl;
	  cout<<endl;
	  //    exit(0);
	}
	else {
	  fclose(test2);
	  
	  cout<<"\nProcessing "<<gOpt->DataFile(il)<<" ...\n"<<endl;
	  
	  gRaw->Open(gOpt->DataFile(il));
	  gRaw->ReadNext(1);
	  
	  int nL20;
	  for (int ctr=0 ; !gRaw->ReadNext() && ctr < gOpt->StopAtEvent(); ctr++) {
	    if (ctr%50000==0) printf("\n%9d ", ctr);
	    else if (ctr%10000==0) cout<<":"<<flush;
	    else if (ctr%1000==0) cout<<"."<<flush;
	    
	    nL20 = gRaw->GetL20nhits();
	    //      cout<<" ************* L20 hits "<<nL20<<" ******** L15 hits "<<gRaw->GetL15nhits()<<endl;
	    if (gRaw.GetEventClass() == 8){ //cut on flasher events
	      for (k=0;k<nL20;k++){
		sect = gRaw->GetL20sector(k);
		pad  = gRaw->GetL20paddle(k);
		tb   = gRaw->GetL20tb(k);
		
		//patch for new LADS numbering
		if (pad < 14) {
		  sect = 0;
		}
		else if (pad > 13 && pad < 28) {
		  pad -= 14;
		}
		
		// **** raw TDC flasher offset ****
		if (sect == 0 && tb == 0 ) {
		  //      && gRaw->GetL20tdc(k) > 0 && gRaw->GetL20tdc(k) < 4095) {
		  tdcrawLEFT_top[pad]->Fill(gRaw->GetL20tdc(k));
		}
		else if (sect == 0 && tb == 1 ) {
		  //      && gRaw->GetL20tdc(k) > 0 && gRaw->GetL20tdc(k) < 4095) {
		  tdcrawLEFT_bot[pad]->Fill(gRaw->GetL20tdc(k));
		}
		else if (sect == 1 && tb == 0 ) {
		  //      && gRaw->GetL20tdc(k) > 0 && gRaw->GetL20tdc(k) < 4095) {
		  tdcrawRIGHT_top[pad]->Fill(gRaw->GetL20tdc(k));
		}
		else if (sect == 1 && tb == 1 ) {
		  //      && gRaw->GetL20tdc(k) > 0 && gRaw->GetL20tdc(k) < 4095) {
		  tdcrawRIGHT_bot[pad]->Fill(gRaw->GetL20tdc(k));
		}
		
		// **** walk corrected TDC flasher offset ****
		double tdc_sqrt2 = gRaw->GetL20adc(k)-ped[sect][tb][pad]-walk_par1[sect][tb][pad];
		if (tdc_sqrt2>0)  {
		  if (sect == 0 && tb == 0 ) {
		    //     && gRaw->GetL20tdc(k) > 0 && gRaw->GetL20tdc(k) < 4095) {
		    tdcLEFT_top[pad]->Fill(gRaw->GetL20tdc(k)-walk_par0[sect][tb][pad]/sqrt(tdc_sqrt2)-1000);
		  }
		  else if (sect == 0 && tb == 1 ) {
		    //     && gRaw->GetL20tdc(k) > 0 && gRaw->GetL20tdc(k) < 4095) {
		    tdcLEFT_bot[pad]->Fill(gRaw->GetL20tdc(k)-walk_par0[sect][tb][pad]/sqrt(tdc_sqrt2)-1000);
		  }
		  else if (sect == 1 && tb == 0 ) {
		    //     && gRaw->GetL20tdc(k) > 0 && gRaw->GetL20tdc(k) < 4095) {
		    tdcRIGHT_top[pad]->Fill(gRaw->GetL20tdc(k)-walk_par0[sect][tb][pad]/sqrt(tdc_sqrt2)-1000);
		  }
		  else if (sect == 1 && tb == 1 ) {
		    //   && gRaw->GetL20tdc(k) > 0 && gRaw->GetL20tdc(k) < 4095) {
		    tdcRIGHT_bot[pad]->Fill(gRaw->GetL20tdc(k)-walk_par0[sect][tb][pad]/sqrt(tdc_sqrt2)-1000);
		  }
		  // cout<<sect<<" "<<pad<<" "<<tb<<endl;
		} // end of tdc_sqrt if
		
	      } // end of flasher if
	    } // end loop over L20 hits

	  } // loop on events
	  gRaw->CodaRead.closedatafile();
	}  //  end of check if data file exists

        // **** data file output ****

        TString outfile("data/L20-tdc-after-walk-");
        outfile += list[il];
        outfile += ".dat";

        ofstream out(outfile, ios::out);

        TString outfile3("data/L20-tdc-");
        outfile3 += list[il];
        outfile3 += ".dat";
 
        ofstream outraw(outfile3, ios::out);

        int mean = 0;
        side_t side = 0;
        pmt_t pmt = 0;
        char txt[200];
        char txt2[200];
   
//        TCanvas *c6 = new TCanvas("l20.LEFT.top.tdc.offset","l20.LEFT_top_tdc", 10, 10, 960, 720);
//        c6->Divide(4,4);
//        c6->Iconify();

        char head[100];
        sprintf(head,"lr\ttb\tbar\tgausmax\tmax-err\tgaus-sigma");
        out<<"! L20 left top  :  TDC offset (after walk effect correction)"<<endl;
        out<<"! "<<head<<endl;
        outraw<<"! L20 left top  :  TDC offset (raw)"<<endl;
        outraw<<"! "<<head<<endl;

        for (k=0; k<14; k++) {
//          c6->cd(k+1);

          mean = tdcLEFT_top[k]->GetMaximumBin();// GetMean();
          gaus1 = new TF1("gaus1","gaus",mean-60,mean+60);    
//          gaus1->SetLineColor(kBlue);

          tdcLEFT_top[k]->GetXaxis()->SetRangeUser(mean-50,mean+50);
//          gPad->Modified();

//          tdcLEFT_top[k]->Draw();

          tdcLEFT_top[k]->Fit("gaus1","WRN");
    
          Double_t par[5]; 
          gaus1->GetParameters(par); 
//          int posmax = par[1];
//          double posmaxerr = gaus1->GetParError(1);
//          double sigma = par[2];

//          char postext[30];
//          sprintf(postext,"%d +/- %.2f",posmax,posmaxerr);

//          int maxi = tdcLEFT_top[k]->GetMaximum();
//          pave = new TPaveText(mean-48,maxi*.9,mean-15,maxi); 
//          pave->SetFillColor(22);
//          TText *t1=pave->AddText(postext);
//          t1->SetTextSize(.05);  
    
//          pave->Draw();

          sprintf(txt,"0\t0\t%d\t%d\t%.3f\t%.3f",k,par[1],gaus1->GetParError(1),par[2]);
          out<<txt<<endl;

          mean = tdcrawLEFT_top[k]->GetMaximumBin();
          gaus1 = new TF1("gaus1","gaus",mean-60,mean+60);    
          tdcrawLEFT_top[k]->GetXaxis()->SetRangeUser(mean-50,mean+50);
          tdcrawLEFT_top[k]->Fit("gaus1","WRN");
          Double_t par[5]; 
          gaus1->GetParameters(par); 
          sprintf(txt,"0\t0\t%d\t%d\t%.3f\t%.3f",k,par[1],gaus1->GetParError(1),par[2]);
          outraw<<txt<<endl;

          //left-top
          side = 0;
          pmt = 0;
          ped_pmt[k+side*14][pmt] = ped[side][pmt][k];
          par_pmt[k+side*14][pmt][0] = walk_par0[side][pmt][k];
          par_pmt[k+side*14][pmt][1] = walk_par1[side][pmt][k];
          par_pmt[k+side*14][pmt][2] = par[1];
        }

//        TCanvas *c7 = new TCanvas("l20.LEFT.bot.tdc.offset","l20.LEFT_bot_tdc", 10, 10, 960, 720);
//        c7->Divide(4,4);
//        c7->Iconify();

        out<<"!"<<endl;
        out<<"! L20 left bottom  :  TDC offset (after walk effect correction)"<<endl;
        out<<"! "<<head<<endl;

        outraw<<"!"<<endl;
        outraw<<"! L20 left bottom  :  TDC offset (raw)"<<endl;
        outraw<<"! "<<head<<endl;

        for (k=0; k<14; k++) {
//          c7->cd(k+1);

          mean = tdcLEFT_bot[k]->GetMaximumBin();// GetMean();
          gaus1 = new TF1("gaus1","gaus",mean-60,mean+60);    
//          gaus1->SetLineColor(kBlue);

//          tdcLEFT_bot[k]->Draw();

          tdcLEFT_bot[k]->GetXaxis()->SetRangeUser(mean-50,mean+50);
//          gPad->Modified();

          tdcLEFT_bot[k]->Fit("gaus1","WRN");
      
          Double_t par[5]; 
          gaus1->GetParameters(par); 
//          int posmax = par[1];
//          double posmaxerr = gaus1->GetParError(1);
//          double sigma = par[2];

//          sprintf(postext,"%d +/- %.2f",posmax,posmaxerr);

//          int maxi = tdcLEFT_bot[k]->GetMaximum();
//          pave = new TPaveText(mean-48,maxi*.9,mean-15,maxi); 
//          pave->SetFillColor(22);
//          TText *t1=pave->AddText(postext);
//          t1->SetTextSize(.05);  
    
//          pave->Draw();

          sprintf(txt,"0\t1\t%d\t%d\t%.3f\t%.3f",k,par[1],gaus1->GetParError(1),par[2]);
          out<<txt<<endl;

          mean = tdcrawLEFT_bot[k]->GetMaximumBin();
          gaus1 = new TF1("gaus1","gaus",mean-60,mean+60);    
          tdcrawLEFT_bot[k]->GetXaxis()->SetRangeUser(mean-50,mean+50);
          tdcrawLEFT_bot[k]->Fit("gaus1","WRN");
          Double_t par[5]; 
          gaus1->GetParameters(par); 
          sprintf(txt,"0\t1\t%d\t%d\t%.3f\t%.3f",k,par[1],gaus1->GetParError(1),par[2]);
          outraw<<txt<<endl;

          //left-bottom
          side = 0;
          pmt = 1;
          ped_pmt[k+side*14][pmt] = ped[side][pmt][k];
          par_pmt[k+side*14][pmt][0] = walk_par0[side][pmt][k];
          par_pmt[k+side*14][pmt][1] = walk_par1[side][pmt][k];
          par_pmt[k+side*14][pmt][2] = par[1];
        }

//        TCanvas *c8 = new TCanvas("l20.RIGHT.top.tdc.offset","l20.RIGHT_top_tdc", 30, 30, 960, 720);
//        c8->Divide(4,4);
//        c8->Iconify();
  
        out<<"!"<<endl;
        out<<"! L20 right top  :  TDC offset (after walk effect correction)"<<endl;
        out<<"! "<<head<<endl;

        outraw<<"!"<<endl;
        outraw<<"! L20 right top  :  TDC offset (raw)"<<endl;
        outraw<<"! "<<head<<endl;

        for (k=0; k<14; k++) {
//          c8->cd(k+1);

          mean = tdcRIGHT_top[k]->GetMaximumBin();// GetMean();
          gaus1 = new TF1("gaus1","gaus",mean-60,mean+60);    
//          gaus1->SetLineColor(kBlue);

//          tdcRIGHT_top[k]->Draw();

          tdcRIGHT_top[k]->GetXaxis()->SetRangeUser(mean-50,mean+50);
//          gPad->Modified();

          tdcRIGHT_top[k]->Fit("gaus1","WRN");
   
          Double_t par[5]; 
          gaus1->GetParameters(par); 
//          int posmax = par[1];
//          double posmaxerr = gaus1->GetParError(1);
//          double sigma = par[2];

//          sprintf(postext,"%d +/- %.2f",posmax,posmaxerr);

//          int maxi = tdcRIGHT_top[k]->GetMaximum();
//          pave = new TPaveText(mean-48,maxi*.9,mean-15,maxi); 
//          pave->SetFillColor(22);
//          TText *t1=pave->AddText(postext);
//          t1->SetTextSize(.05);  
    
//          pave->Draw();

          sprintf(txt,"1\t0\t%d\t%d\t%.3f\t%.3f",k,par[1],gaus1->GetParError(1),par[2]);
          out<<txt<<endl;

          mean = tdcrawRIGHT_top[k]->GetMaximumBin();// GetMean();
          gaus1 = new TF1("gaus1","gaus",mean-60,mean+60);    
          tdcrawRIGHT_top[k]->GetXaxis()->SetRangeUser(mean-50,mean+50);
          tdcrawRIGHT_top[k]->Fit("gaus1","WRN");
          Double_t par[5]; 
          gaus1->GetParameters(par); 
          sprintf(txt,"1\t0\t%d\t%d\t%.3f\t%.3f",k,par[1],gaus1->GetParError(1),par[2]);
          outraw<<txt<<endl;

          //right-top
          side = 1;
          pmt = 0;
          ped_pmt[k+side*14][pmt] = ped[side][pmt][k];
          par_pmt[k+side*14][pmt][0] = walk_par0[side][pmt][k];
          par_pmt[k+side*14][pmt][1] = walk_par1[side][pmt][k];
          par_pmt[k+side*14][pmt][2] = par[1];
        }

//        TCanvas *c9 = new TCanvas("l20.RIGHT.bot.tdc.offset","l20.RIGHT_bot_tdc", 30, 30, 960, 720);
//        c9->Divide(4,4);
//        c9->Iconify();

        out<<"!"<<endl;
        out<<"! L20 right bottom  :  TDC offset (after walk effect correction)"<<endl;
        out<<"! "<<head<<endl;

        outraw<<"!"<<endl;
        outraw<<"! L20 right bottom  :  TDC offset (raw)"<<endl;
        outraw<<"! "<<head<<endl;

        for (k=0; k<14; k++) {
//          c9->cd(k+1);

          mean = tdcRIGHT_bot[k]->GetMaximumBin();// GetMean();
          gaus1 = new TF1("gaus1","gaus",mean-60,mean+60);    
//          gaus1->SetLineColor(kBlue);

//          tdcRIGHT_bot[k]->Draw();

          tdcRIGHT_bot[k]->GetXaxis()->SetRangeUser(mean-50,mean+50);
//          gPad->Modified();

          tdcRIGHT_bot[k]->Fit("gaus1","WRN");
    
          Double_t par[5]; 
          gaus1->GetParameters(par); 
//          int posmax = par[1];
//          double posmaxerr = gaus1->GetParError(1);
//          double sigma = par[2];

//          sprintf(postext,"%d +/- %.2f",posmax,posmaxerr);

//          int maxi = tdcRIGHT_bot[k]->GetMaximum();
//          pave = new TPaveText(mean-48,maxi*.9,mean-15,maxi); 
//          pave->SetFillColor(22);
//          TText *t1=pave->AddText(postext);
//          t1->SetTextSize(.05);  
    
//          pave->Draw();

          sprintf(txt,"1\t1\t%d\t%d\t%.3f\t%.3f",k,par[1],gaus1->GetParError(1),par[2]);
          out<<txt<<endl;

          mean = tdcrawRIGHT_bot[k]->GetMaximumBin();
          gaus1 = new TF1("gaus1","gaus",mean-60,mean+60); 
          tdcrawRIGHT_bot[k]->GetXaxis()->SetRangeUser(mean-50,mean+50);
          tdcrawRIGHT_bot[k]->Fit("gaus1","WRN"); 
          Double_t par[5]; 
          gaus1->GetParameters(par); 
          sprintf(txt,"1\t1\t%d\t%d\t%.3f\t%.3f",k,par[1],gaus1->GetParError(1),par[2]);
          outraw<<txt<<endl;
   
          //right-bottom
          side = 1;
          pmt = 1;
          ped_pmt[k+side*14][pmt] = ped[side][pmt][k];
          par_pmt[k+side*14][pmt][0] = walk_par0[side][pmt][k];
          par_pmt[k+side*14][pmt][1] = walk_par1[side][pmt][k];
          par_pmt[k+side*14][pmt][2] = par[1];
        }

        out.close();
        outraw.close();

        // **** cal-file output ****

        ofstream fp_out;
        TString outfile2("data/L20_cal-");
        outfile2 += list[il];
        outfile2 += ".dat";
        fp_out.open(outfile2, ios::out);
        fp_out << "# ****************************L20**************************" << endl;
        fp_out << "# ****                                                 ****" << endl;
        fp_out << "# ****  THIS is a valid replacement for the entire L20 ****" << endl;
        fp_out << "# ****  within the blast.sc_cal file                   ****" << endl;
        fp_out << "# ****  produced by l20_tdc_offsets_walk.C             ****" << endl;
        fp_out << "# ****                                                 ****" << endl;
        fp_out << "# *********************************************************" << endl;
        fp_out << "#" << endl;
        fp_out << "#   sect     pad     t1off    t2off     t1sca   t2sca   a1off"
               << "   a2off   a1sca   a2sca   t1par0    t1par1    t2par0    t2par1" << endl;

        for (int yy = 0; yy < 28; yy++){      // loop over the pads in l20 canvas
          sprintf(txt2,"       1%8d%10d%10d      50      50%8d%8d      50      50%10.2lf%10d%10.2lf%10d",
                  yy, par_pmt[yy][0][2], par_pmt[yy][1][2], ped_pmt[yy][0], ped_pmt[yy][1], par_pmt[yy][0][0],
                  par_pmt[yy][0][1], par_pmt[yy][1][0], par_pmt[yy][1][1]);
          fp_out << txt2 << endl;
        }
        fp_out.close();

      }  //  end of check if ADC pedestal file exists

    } // loop on files

  }  //  end of check if walkparameter file exists

}