int const MAXA=5;
int const MAXB=5;
int const MAX=MAXA+MAXB;

int list[MAX] ={4139,4140,4141,4142,4143,4145,4146,4147,4148,4149};//omit 4144
int listA[MAXA] ={4139,4140,4141,4142,4143};
int listB[MAXB] = {4145,4146,4147,4148,4149};

int const MAX=10;


gStyle->SetPalette(1);


void replay(int runnumber=1184,int option = 1){

  // output variable define in output.def
  char datafile[300] = Form("/cache/mss/halla/e07006/raw/e07006_%d.dat.0",runnumber);
  if(option==1) char rootfile[300] = Form("pedestal_%d.root",runnumber); //use this to run pedestal
  if(option==2)  char rootfile[300] = Form("bpm_%d.root",runnumber); //use this to run bpm
  if(option==3)  char rootfile[300] = Form("raster_%d.root",runnumber); //raster beam

  if(option==3){
    THaApparatus* BEAM2 = new THaRasteredBeam("rb","Beamline");
    gHaApps->Add( BEAM2 );
  }
  
  THaApparatus* BEAM1 = new THaUnRasteredBeam("urb","Beamline");
  gHaApps->Add( BEAM1 );
 
 
  THaEvent* event = new THaEvent;
  
  THaAnalyzer* analyzer = new THaAnalyzer;

  THaRun* run = new THaRun(datafile);

  analyzer->SetEvent( event );
  analyzer->SetOutFile(rootfile);

  run->SetLastEvent( -1 );
  
  analyzer->Process(*run);

}

////////////////

void draw_bpm_pedestals(int run=1184){
  TFile *file = new TFile(Form("pedestal_%d.root",run));
  TTree *T = (TTree*)file->Get("T");

  //histogram & Fit prepare//
  TF1 *fit[2][4];
  double mean[2][4];
  TH1F *hist[2][4];
  char hist_name[2][4][300];
  char plot_var[2][4][300];
  char plot_hist[2][4][300];

  char BPM[2][300] = {"BPMA","BPMB"};

  for(int i=0;i<2;i++){
    for(int j=0;j<4;j++){
      sprintf(plot_var[i][j],"urb.%s.rawcur.%d",BPM[i],j+1);
      sprintf(hist_name[i][j],"%s_rawcur_%d",BPM[i],j+1);
      hist[i][j] = new TH1F(hist_name[i][j],hist_name[i][j],100,0,1000);
      sprintf(plot_hist[i][j],"%s >> %s",plot_var[i][j],hist_name[i][j]);
    }
  }


  //plot//  
  TCanvas *c = new TCanvas("BPM_ped",Form("BPM pedestals run %d",run), 1000,800);
  c->Divide(2,4);
  int count=0;
  for(int i=0;i<2;i++){
    for(int j=0;j<4;j++){
      count++;
      c->cd(count);
      T->Draw(plot_hist[i][j]);
      hist[i][j]->Fit("gaus");
      fit[i][j] = hist[i][j]->GetFunction("gaus");
      mean[i][j] = fit[i][j]->GetParameter("Mean");
      c->Update();
    }
  }

  c->Print(Form("BPM_pedestals_run_%d.png",run));

  //yeild//
  for(int i=0;i<2;i++){
    cout<< "change the pedestal for"<<  BPM[i]<< " to :"<<endl;
    for(int j=0;j<4;j++){
      cout<<(int)mean[i][j]<<" ";
    }
    cout<<endl;
  }

  //record//
  if(1){
    char outfile[300] =Form("BPM_pedestals_%d.dat",run);
    FILE *fw = fopen(outfile,"wb");
    for(int i=0;i<2;i++){
      fprintf(fw, "change the pedestal for %s to :\n",BPM[i]);
      for(int j=0;j<4;j++){
	fprintf(fw,"%d   ",(int)mean[i][j]);
      }
      fprintf(fw,"\n");
    }
  }
  
  fclose(fw);

}



void plot(int option=1){
  gStyle->SetPalette(1);
  ///(1)all file uses (2)A (3)B /////////////////////////////////

  //2d hist  
  int const hnum = 4;

  TH2F *hist[hnum];
  char hist_name[hnum][300]={"BPMA","BPMB","epicA","epicB"};
  char plot_hist[hnum][300];
  char plot_2d[hnum][300];
  char plot_var[2*hnum][300];
  char BPM[2*hnum][300] = {"BPMA","BPMA","BPMB","BPMB","IPM1H04A","IPM1H04A","IPM1H04B","IPM1H04B"};
  char id[2*hnum][300] = {"rotpos1","rotpos2","rotpos1","rotpos2","XPOS","YPOS","XPOS","YPOS"};
  int bin[2*hnum]={1000,1000,1000,1000,100,100,100,100};
  int min[2*hnum] = {-4,-4,-3, 5,-3,-3,-5,-1};
  int max[2*hnum] = { 4,5, 5,13, 3, 3, 2, 5};
  
  for(int i=0;i<2*hnum;i++){ 
    if(i<hnum)    sprintf(plot_var[i],"urb.%s.%s*1e3",BPM[i],id[i]);
    else    sprintf(plot_var[i],"%s.%s",BPM[i],id[i]);
  }

  for(int i=0;i<hnum;i++){
    hist[i] = new TH2F(hist_name[i],hist_name[i],bin[2*i],min[2*i],max[2*i],bin[2*i+1],min[2*i+1],max[2*i+1]);
    sprintf(plot_2d[i],"%s:%s",plot_var[2*i+1],plot_var[2*i]); 
    sprintf(plot_hist[i],"%s>> %s",plot_2d[i],hist_name[i]); 

    //cout<<plot_hist[i]<<endl;
    hist[i]->GetXaxis()->SetTitle(Form("%s [mm]",id[2*i]));
    hist[i]->GetYaxis()->SetTitle(Form("%s [mm]",id[2*i+1]));
  }
  //

    TCanvas *c3 = new TCanvas("2d","2d",1000,800);
    c3->Divide(2,2);

    ///////////////////////////////
  if(option==1){
    TChain *T = new TChain("T");
    TChain *E = new TChain("E");
  
    for(int k=0;k<MAX;k++){
      T->Add(Form("bpm_%d.root",list[k]));
      E->Add(Form("bpm_%d.root",list[k]));
    } 
   
  ///plots///
  
    for(int i=0;i<hnum/2;i++){
	c3->cd(i+1); T->Draw(plot_hist[i],"","colz");  
	c3->cd(i+3);E->Draw(plot_hist[i+2],"","");
	hist[i+2]->SetMarkerStyle(3);
	hist[i+2]->Draw();

    }
    c3->Print("all_runs.png");
  }
  ////////////////////////////////
  else if(option==2 || option==3){

    TH1F *h1d[2*hnum];
    char h1d_name[2*hnum][300];
    char plot_h1d[2*hnum][300];

    for(int i=0;i<2*hnum;i++){
      sprintf(h1d_name[i],"%s_%s",hist_name[i/2],id[i]);
      h1d[i] = new TH1F(h1d_name[i],h1d_name[i],bin[i],min[i],max[i]);
      sprintf(plot_h1d[i],"%s>>%s",plot_var[i],h1d_name[i]);
      h1d[i]->GetXaxis()->SetTitle(Form("%s [mm]",id[i]));
    }
    int max_run;
    if(option==2) max_run=MAXA; 
    if(option==3) max_run=MAXB;

    TChain *T[MAX];
    TChain *E[MAX];

    TCanvas *c1 = new TCanvas("1d","1d",1000,800);
    c1->Divide(4,2);

    for(int k=0;k<max_run;k++){

      T[k] = new TChain("T");
      if(option==2)  T[k]->Add(Form("bpm_%d.root",listA[k]));
      if(option==3) T[k]->Add(Form("bpm_%d.root",listB[k]));
      T[k]->SetMarkerColor(k);
      T[k]->SetLineColor(k);
      
      E[k] = new TChain("E");
      if(option==2) E[k]->Add(Form("bpm_%d.root",listA[k]));
      if(option==3) E[k]->Add(Form("bpm_%d.root",listB[k]));
      E[k]->SetMarkerColor(k);
      E[k]->SetLineColor(k);
      E[k]->SetMarkerStyle(3);

      ///////plot 2d
      for(int i=0;i<hnum/2;i++){
	if(k==0){
	  c3->cd(i+1);T[k]->Draw(plot_hist[i],"","colz");
	  c3->cd(i+2+1);E[k]->Draw(plot_hist[i+2],"","");
	  hist[i+2]->SetMarkerStyle(3);
	  hist[i+2]->Draw();
	  c3->Update();
	}
	else {
	  c3->cd(i+1);T[k]->Draw(plot_2d[i],"","same");
	  c3->cd(i+2+1);E[k]->Draw(plot_2d[i+2],"","same");
	  c3->Update();
	}
      }//plot

      for(int i=0;i<hnum;i++){
	if(k==0){
	  c1->cd(i+1);T[k]->Draw(plot_h1d[i],"","");
	  c1->cd(i+4+1);E[k]->Draw(plot_h1d[i+4],"","");
	  h1d[i+4]->SetMarkerStyle(3);
	  h1d[i+4]->Draw();
	  c1->Update();
	}
	else{
	  c1->cd(i+1);T[k]->Draw(plot_var[i],"","same");
	  c1->cd(i+4+1);E[k]->Draw(plot_var[i+4],"","same");
	  c1->Update();
	}
      }


    }
    if(option==2){
     c1->Print("1d_rotepic_BPMA.png");
     c3->Print("2d_rotepic_BPMA.png");
   }
   if(option==3){
     c1->Print("1d_rotepic_BPMB.png");
     c3->Print("2d_rotepic_BPMB.png");
   }
  }
  
} 
/////////////////////////////////

void plot_all(){
  plot(1);
}

void plot_AB(int choice=1){
  plot(choice+1);
}


////////////////////////////////

void get_rot_epic(int runnumber){
  //get mean /////////////////////////////////
  FILE *fw = fopen("epic_rot.dat","a");
  TChain *T;
  TChain *E;
  
  
  TCanvas *c1 = new TCanvas("1d","1d",1000,800);
  c1->Divide(4,2);
  
  TF1 *fit;
  double mrax,mray,mrbx,mrby,meax,meay,mebx,meby;
  
  
  T = new TChain("T");
  T->Add(Form("bpm_%d.root",runnumber));
  
  E = new TChain("E");
  E->Add(Form("bpm_%d.root",runnumber));
  E->SetMarkerStyle(3);
  
  c1->cd(1);T->Draw("urb.BPMA.rotpos1*1e3>>rax(1000,-4,4)");
  rax->Fit("gaus"); fit = rax->GetFunction("gaus"); mrax = fit->GetParameter("Mean");
  
  c1->cd(2);T->Draw("urb.BPMA.rotpos2*1e3>>ray(1000,-4,5)");
  
  ray->Fit("gaus"); fit = ray->GetFunction("gaus"); mray = fit->GetParameter("Mean");
  
  c1->cd(3);T->Draw("urb.BPMB.rotpos1*1e3>>rbx(1000,-3,5)");
  rbx->Fit("gaus"); fit = rbx->GetFunction("gaus"); mrbx = fit->GetParameter("Mean");
  
  c1->cd(4);T->Draw("urb.BPMB.rotpos2*1e3>>rby(1000,5,13)");
  rby->Fit("gaus"); fit = rby->GetFunction("gaus"); mrby = fit->GetParameter("Mean");
  
  c1->cd(5); E ->Draw("IPM1H04A.XPOS>>eax(100,-3,3)");
  meax = eax->GetMean();
  
  c1->cd(6); E ->Draw("IPM1H04A.YPOS>>eay(100,-3,3)");
  meay = eay->GetMean();
  c1->cd(7); E ->Draw("IPM1H04B.XPOS>>ebx(100,-5,1)");
  mebx = ebx->GetMean();
  
  c1->cd(8); E ->Draw("IPM1H04B.YPOS>>eby(100,-1,5)");
  meby = eby->GetMean();
  
  c1->Print(Form("rotepic_mean_%d.png",runnumber));
  
  cout<<runnumber<<"\t"
      <<mrax<<"\t"<<mray<<"\t"<<mrbx<<"\t"<<mrby<<"\t"
      <<meax<<"\t"<<meay<<"\t"<<mebx<<"\t"<<meby<<endl;
  
  fprintf(fw,"%d\t %lf\t %lf\t %lf\t %lf\t %lf\t %lf\t %lf\t %lf\n",runnumber,mrax,mray,mrbx,mrby,meax,meay,mebx,meby);
  fclose(fw);
}
////////////////////////////////


void repeat(){
  for(int i=0;i<MAXA;i++)  get_rot_epic(listA[i]);
  for(int i=0;i<MAXB;i++)  get_rot_epic(listB[i]);

}

void calculate(){

    FILE *fo = fopen("epic_rot.dat","r");
    int runnumber[10];
    double rax[10],ray[10],rbx[10],rby[10],eax[10],eay[10],ebx[10],eby[10];
    char buff[300];
    char *filestatus;
    int count=0;
    while(fgets(buff,255,fo)!=NULL){
    
      sscanf(buff,"%d\t %lf\t %lf\t %lf\t %lf\t %lf\t %lf\t %lf\t %lf\n",
	     &runnumber[count],
	     &rax[count],&ray[count],&rbx[count],&rby[count],
	     &eax[count],&eay[count],&ebx[count],&eby[count]);
      count++;

    }
    fclose(fo);

    double srx=0,srxx=0,sry=0,sryy=0,srxy=0;
    double sxe=0,sxerx=0,sxery=0;
    double sye=0,syerx=0,syery=0;

    TMatrix MA(3,3);
    TVector ex(3),ey(3);
    for(int i=0;i<MAXA;i++){
      //calA
      srx +=rax[i]; srxx+=rax[i]*rax[i];
      sry +=ray[i]; sryy+=ray[i]*ray[i];
      srxy+=rax[i]*ray[i];

      sxe +=eax[i]; sxerx+=eax[i]*rax[i]; sxery+=eax[i]*ray[i];
      sye +=eay[i]; syerx+=eay[i]*rax[i]; syery+=eay[i]*ray[i];
    }

    MA(0,0) = srx;  MA(0,1) = sry;  MA(0,2) = MAXA;
    MA(1,0) = srxx; MA(1,1) = srxy; MA(1,2) = srx;
    MA(2,0) = srxy; MA(2,1) = sryy; MA(2,2) = sry;

    ex(0) = sxe; ex(1) = sxerx; ex(2) = sxery;
    ey(0) = sye; ey(1) = syerx; ey(2) = syery;

    TMatrix iMA(MA); iMA.Invert();
    TVector solA1(ex); solA1*=iMA;
    TVector solA2(ey); solA2*=iMA;

    /////////////////////////////
    //reset
    srx=0;srxx=0;sry=0;sryy=0;srxy=0;
    sxe=0;sxerx=0;sxery=0;
    sye=0;syerx=0;syery=0;
    

    TMatrix MB(3,3);
    TVector ex(3),ey(3);
    int runB = 5;
    for(int i=MAXA;i<MAXA+MAXB;i++){
      //calB
      srx +=rbx[i]; srxx+=rbx[i]*rbx[i];
      sry +=rby[i]; sryy+=rby[i]*rby[i];
      srxy+=rbx[i]*rby[i];

      sxe +=ebx[i]; sxerx+=ebx[i]*rbx[i]; sxery+=ebx[i]*rby[i];
      sye +=eby[i]; syerx+=eby[i]*rbx[i]; syery+=eby[i]*rby[i];
    }


    MB(0,0) = srx;  MB(0,1) = sry;  MB(0,2) = MAXB;
    MB(1,0) = srxx; MB(1,1) = srxy; MB(1,2) = srx;
    MB(2,0) = srxy; MB(2,1) = sryy; MB(2,2) = sry;

    ex(0) = sxe; ex(1) = sxerx; ex(2) = sxery;
    ey(0) = sye; ey(1) = syerx; ey(2) = syery;

    TMatrix iMB(MB); iMB.Invert();
    TVector solB1(ex); solB1*=iMB;
    TVector solB2(ey); solB2*=iMB;

    double A[6], B[6];



    A[0] = solA1(0);
    A[1] = solA1(1);
    A[2] = solA2(0);
    A[3] = solA2(1);
    A[4] = solA1(2);
    A[5] = solA2(2);

    B[0] = solB1(0);
    B[1] = solB1(1);
    B[2] = solB2(0);
    B[3] = solB2(1);
    B[4] = solB1(2);
    B[5] = solB2(2);

    cout<<A[0]<<"\t"<<A[1]<<"\t"<<A[2]<<"\t"<<A[3]<<"\t"<<A[4]<<"\t"<<A[5]<<"\n";
    cout<<B[0]<<"\t"<<B[1]<<"\t"<<B[2]<<"\t"<<B[3]<<"\t"<<B[4]<<"\t"<<B[5]<<"\n";

    FILE *fw = fopen("BPM_constant.dat","wb");
    fprintf(fw,"%lf\t %lf\t %lf\t %lf\t %lf\t %lf\n",A[0],A[1],A[2],A[3],A[4],A[5]);
    fprintf(fw,"%lf\t %lf\t %lf\t %lf\t %lf\t %lf\n",B[0],B[1],B[2],B[3],B[4],B[5]);
    fclose(fw);
}


void test(){
  FILE *fo = fopen("epic_rot.dat","r");
  FILE *f2 = fopen("BPM_constant.dat","r");    
  int runnumber[10];
  double rax[10],ray[10],rbx[10],rby[10],eax[10],eay[10],ebx[10],eby[10];
  char buff[300];
  char *filestatus;
  int count=0;
  while(fgets(buff,255,fo)!=NULL){
    
    sscanf(buff,"%d\t %lf\t %lf\t %lf\t %lf\t %lf\t %lf\t %lf\t %lf\n",
	   &runnumber[count],
	   &rax[count],&ray[count],&rbx[count],&rby[count],
	   &eax[count],&eay[count],&ebx[count],&eby[count]);
    count++;
    
  }
  fclose(fo);
  
  double A[6], B[6];
  
  fscanf(f2,"%lf\t %lf\t %lf\t %lf\t %lf\t %lf\n",&A[0],&A[1],&A[2],&A[3],&A[4],&A[5]);
  fscanf(f2,"%lf\t %lf\t %lf\t %lf\t %lf\t %lf\n",&B[0],&B[1],&B[2],&B[3],&B[4],&B[5]);
  fclose(f2);

  for(int i=0;i<6;i++){
    cout<<A[i]<<"\t"<<B[i]<<endl;
  }

  TMatrix MA(2,2),MB(2,2);
  TVector offA(2),offB(2);

  MA(0,0) = A[0]; MA(0,1) = A[1]; MA(1,0) = A[2]; MA(1,1) = A[3]; offA(0) = A[4]; offA(1) = A[5];
  MB(0,0) = B[0]; MB(0,1) = B[1]; MB(1,0) = B[2]; MB(1,1) = B[3]; offB(0) = B[4]; offB(1) = B[5];


  TVector rotA(2),rotB(2);
  double calAx[10],calAy[10],calBx[10],calBy[10];

  for(int i=0;i<count;i++){
    rotA(0) = rax[i];
    rotA(1) = ray[i];

    rotB(0) = rbx[i];
    rotB(1) = rby[i];

    TVector eA(rotA); eA*=MA;
    TVector eB(rotB); eB*=MB;
    calAx[i] = eA(0)+offA(0);
    calAy[i] = eA(1)+offA(1);
    calBx[i] = eB(0)+offB(0);
    calBy[i] = eB(1)+offB(1);

    cout<<runnumber[i]<<"\t";
    cout<<calAx[i]<<":"<<eax[i]<<"\t";
    cout<<calAy[i]<<":"<<eay[i]<<"\t";
    cout<<calBx[i]<<":"<<ebx[i]<<"\t";
    cout<<calBy[i]<<":"<<eby[i]<<"\n";

 }

  TGraph *Ax = new TGraph(count,calAx,eax);Ax->SetTitle("x_calculated vs x_epic:BPMA"); Ax->GetYaxis()->SetTitle("epic x (mm)");Ax->GetXaxis()->SetTitle("calculate x (mm)");
  TGraph *Ay = new TGraph(count,calAy,eay);Ay->SetTitle("y_calculated vs y_epic:BPMA"); Ay->GetYaxis()->SetTitle("epic y (mm)");Ay->GetXaxis()->SetTitle("calculate y (mm)");
  TGraph *Bx = new TGraph(count,calBx,ebx);Bx->SetTitle("x_calculated vs x_epic:BPMB");Bx->GetYaxis()->SetTitle("epic x (mm)");Bx->GetXaxis()->SetTitle("calculate x (mm)");
  TGraph *By = new TGraph(count,calBy,eby);By->SetTitle("y_calculated vs y_epic:BPMB");By->GetYaxis()->SetTitle("epic y (mm)");By->GetXaxis()->SetTitle("calculate y (mm)");

  gStyle->SetOptFit(0111);
  TCanvas *c = new TCanvas("c","c",1000,800);
  c->Divide(2,2);
  c->cd(1); Ax->Draw("A*");Ax->Fit("pol1");
  c->cd(2); Ay->Draw("A*");Ay->Fit("pol1");
  c->cd(3); Bx->Draw("A*");Bx->Fit("pol1");
  c->cd(4); By->Draw("A*");By->Fit("pol1");

  c->Print("constant_test.png");
}



void cal_raster(int runnumber){

  TFile *file = new TFile(Form("raster_%d.root",runnumber));
  TTree *T = (TTree*)file->Get("T");

  const int id1=5;
  const int id2=3;
  //histogram & Fit prepare//
  TF1 *fit[id1][id2];
  double mean[id1][id2];
  double sigma[id1][id2];
  TH1F *hist[id1][id2];
  char hist_name[id1][id2][300];
  char plot_var[id1][id2][300];
  char plot_hist[id1][id2][300];

  char BPM[id1][300] = {"BPMA","BPMB","target","rawcur","rawslope"};
  char cord[id2][100]={"x","y","z"};

  for(int i=0;i<id1;i++){
    for(int j=0;j<id2;j++){

      if(i<2) {
	sprintf(plot_var[i][j],"urb.%s.%s",BPM[i],cord[j]);
	sprintf(hist_name[i][j],"urb_%s_raw_cur_%d",BPM[i],j+1);
	hist[i][j] = new TH1F(hist_name[i][j],hist_name[i][j],4000,-0.005,0.005);
      }
      else if(i==2){
	sprintf(plot_var[i][j],"urb.%s",cord[j]);
	sprintf(hist_name[i][j],"urb_raw_cur_%d",j+1);
      hist[i][j] = new TH1F(hist_name[i][j],hist_name[i][j],4000,-0.005,0.005);
      }
      else{
	sprintf(plot_var[i][j],"rb.Raster.%s.%s",BPM[i],cord[j]);
	sprintf(hist_name[i][j],"raster_%s_%d",BPM[i],j+1);
	hist[i][j] = new TH1F(hist_name[i][j],hist_name[i][j],4000,1000,5000);
      }
      sprintf(plot_hist[i][j],"%s >> %s",plot_var[i][j],hist_name[i][j]);
    }
  }


  //plot//  
  TCanvas *c = new TCanvas("Raster",Form("Raster Constant for run %d",runnumber), 1000,800);
  c->Divide(2,5);
  int count=0;
  for(int i=0;i<id1;i++){
    for(int j=0;j<id2-1;j++){
      count++;
      c->cd(count);
      T->Draw(plot_hist[i][j]);
      mean[i][j] = hist[i][j]->GetMean();
      sigma[i][j] = hist[i][j]->GetRMS();
      c->Update();
    }
  }

  ////////////////////////

  double  rax = mean[3][0];
  double  ray = mean[3][1];
  double  bax = mean[0][0];
  double  bay = mean[0][1];
  double  bbx = mean[1][0];
  double  bby = mean[1][1];
  double  tax = mean[2][0];
  double  tay = mean[2][1];
  
  double drax = sigma[3][0];
  double dray = sigma[3][1];
  double  dbax = sigma[0][0];
  double  dbay = sigma[0][1];
  double  dbbx = sigma[1][0];
  double  dbby = sigma[1][1];
  double  dtax = sigma[2][0];
  double  dtay = sigma[2][1];
 
  double kappax = 1.0;
  double kappay = -1.0;
 
  Double_t bpma_ofx= bax-rax*dbax/drax*kappax ;  
  Double_t bpma_slx= dbax/drax*kappax ; 
  Double_t bpma_ofy= bay-ray*dbay/dray*kappay; 
  Double_t bpma_sly= dbay/dray*kappay; 
  
  Double_t bpmb_ofx= bbx-rax*dbbx/drax*kappax ; 
  Double_t bpmb_slx= dbbx/drax*kappax ; 
  Double_t bpmb_ofy= bby-ray*dbby/dray*kappay; 
  Double_t bpmb_sly= dbby/dray*kappay; 
  
  Double_t tar_ofx= tax-rax*dtax/drax*kappax ; 
  Double_t tar_slx= dtax/drax*kappax ; 
  Double_t tar_ofy= tay-ray*dtay/dray*kappay; 
  Double_t tar_sly= dtay/dray*kappay; 
  
  
  
  cout<<" Raster Constants: (please modify database accordingly)"<<endl;
  cout<<bpma_ofx<<" "<<bpma_ofy<<" "<<bpma_slx<<" "<<bpma_sly<<" 0. 0."<<endl;
  cout<<bpmb_ofx<<" "<<bpmb_ofy<<" "<<bpmb_slx<<" "<<bpmb_sly<<" 0. 0."<<endl;
  cout<<tar_ofx<<" "<<tar_ofy<<" "<<tar_slx<<" "<<tar_sly<<" 0. 0."<<endl;

    
}

void singlestep_raster(int runnumber){
  replay(runnumber,3);
  cal_raster(runnumber);
}