//*****************************************************************************
//       FitBeam <file name> <runnumber> <history> <dbid>
//  this executable opens a ROOT file with the summary histograms and 
//  fits the beamline for the COT and the silicon detectors. It then
//  puts the fit results into the online production database for the
//  specified runnumber. 
//
//      05/31/01  Hartmut Stadie  IEKP Karlsruhe
//*****************************************************************************
//=============================================================================
// RCS Current Revision Record
//-----------------------------------------------------------------------------
// $Source: /cdf/code/cdfcvs/run2/Alignment/test/FitBeam.cc,v $
// $Revision: 1.26 $
// $Date: 2003/05/15 04:46:22 $
// $Author: rlc $
// $State: Exp $
// $Locker:  $
//*****************************************************************************

#include <cmath>
#include <iostream>
#include <sstream>
#include <map>
#include <string>

#include "TROOT.h"
#include "TFile.h"
#include "TH1.h"
#include "TH2.h"
#include "TProfile.h"
#include "TF1.h"
#include "TGraphErrors.h"

#include "CalibDB/RunListKey.hh"
#include "AlignmentDB/CotBeamPosition.Defs.hh"        
#include "AlignmentDB/SvxBeamPosition.Defs.hh"
#include "AlignmentDB/SiAlignFrame.Defs.hh"
#include "CLHEP/Matrix/Vector.h"
#include "CLHEP/Matrix/SymMatrix.h"
#include "CalibDB/Helpers.hh"
#include "CalibDB/RunListDefs.hh"


struct beamline
{
  double x0,y0,slopex,slopey;
  double sigmax0,sigmay0,sigmaslopex,sigmaslopey;
  double sxy, sxsx, sxsy,sysx,sysy,ssxsy;
  double widthx,widthy;
  double emittxfit, bstarxfit, zzeroxfit, emittyfit, bstaryfit, zzeroyfit;
  double beamz, widthz;   
  double minimzfit, widthzfit, bstarzfit, zzerozfit;
  double statistics0;
  int algo;
  int history;
  int widthquality1, widthquality2;
  void reset()
  {
    x0          = 0;
    y0          = 0;
    slopex      = 0;
    slopey      = 0;
    sigmax0     = 0;
    sigmay0     = 0;
    sigmaslopex = 0;
    sigmaslopey = 0;
    sxy         = 0;
    sxsx        = 0;
    sxsy        = 0;
    sysx        = 0;
    sysy        = 0;
    ssxsy       = 0;
    widthx      = 0;
    widthy      = 0;
    beamz       = 0;
    widthz      = 0;   
    minimzfit   = 0;
    widthzfit   = 0;
    bstarzfit   = 0;
    zzerozfit   = 0;
    emittxfit = 0;
    bstarxfit = 0;
    zzeroxfit = 0;
    emittyfit = 0;
    bstaryfit = 0;
    zzeroyfit = 0;
    statistics0 = 0;
    algo        = 0;
    widthquality1 = 1;
    widthquality2 = 1;
  }
}; // end struct beamline


const int ZNUM = 33;
TH1F *_svxxwidth[ZNUM], *_svxywidth[ZNUM];



struct globalalignment
{
  double tx,ty,slopex,slopey;
};
 

// names of the alignment tables
const int NUM_ALIGNMENT_NAMES = 5;
const int MAX_DIV = 15;

std::string alignmentNames[NUM_ALIGNMENT_NAMES] = { "SiAlignFrame","SiAlignBarrel","SiAlignLadder",
					       "SiAlignWafer","SiAlignWarp" };

std::vector<CID> getAlignmentCIDs(std::string dbid, int run, int version, std::string status)
{ 
  std::vector<CID> cids;
  RunList_mgr runmanager(dbid,"RunList"); 
  for(int i = 0 ; i < NUM_ALIGNMENT_NAMES ; ++i)
    {
      // get the key    
      RunListKey rlk(status,alignmentNames[i],run,version);
      RunList_var row;  
      if( (runmanager.get(rlk,row) != Result::success) ||
	  (row->size() != 1) ) continue;
      RunList::iterator first = row->begin();
      cids.push_back(first->id());
    }
  return cids;
}

bool readGlobalAlignment(std::string dbid, int run, int version, std::string status,globalalignment* align)
{ 
  SiAlignFrame_mgr alignmanager(dbid,"SiAlignFrame");
  RunListKey rlk(status,"SiAlignFrame",run,version);
  SiAlignFrameContainer_ptr aligncont;
  if( alignmanager.get(rlk,aligncont) != Result::success ) return false;
  SiAlignFrame ali = *( aligncont->begin());
  align->tx = ali.Tx();
  align->ty = ali.Ty();
  align->slopex =  ali.alpha_y();
  align->slopey = -ali.alpha_x();
  return true;
}

bool fitDPhi(beamline *beam, TH1F* hdphi)
{
   HepSymMatrix V(4,0);
   HepVector b(4,0);
   float ntracks = 0;
   
   // fill vector and matrix
   for(int i = 0, k = 1; i < 4 ; ++i)
    {
      b[i] = hdphi->GetBinContent(hdphi->FindBin(12+i));
      if(b[i] == 0) return false;
      for(int j = i; j < 4 ; ++j)
	{
	  V[i][j] = hdphi->GetBinContent(hdphi->FindBin(k));
	  if(V[i][j] == 0) return false;
	  ++k;
	}
    }
   //    for(int i = 0; i < 4 ; ++i)
   //      {
   //        b[i] = hdphi->GetBinContent(hdphi->FindBin(12+i));	
   //        for(int j = i; j < 4 ; ++j)
   // 	 {
   // 	   V[i][j] = hdphi->GetBinContent(hdphi->FindBin(i * 4 + j));
   // 	 }
   //      }
   ntracks = hdphi->GetBinContent(hdphi->FindBin(16));
  
   //fit the beamline
   int ierr;
   HepSymMatrix G = V.inverse(ierr);
   if( ierr != 0 ) 
     {
       return false;
     }
   HepVector x = G * b;
   //std::cout << x << G; 
   beam->x0          = x[0];
   beam->slopex      = x[2];
   beam->sigmax0     = G[0][0];
   beam->sigmaslopex = G[2][2];
   beam->y0          = x[1];
   beam->slopey      = x[3];
   beam->sigmay0     = G[1][1];
   beam->sigmaslopey = G[3][3];
   beam->statistics0 = ntracks;
   beam->sxy         = G[0][1];
   beam->sxsx        = G[0][2];
   beam->sxsy        = G[0][3];
   beam->sysx        = G[1][2];
   beam->sysy        = G[1][3];
   beam->ssxsy       = G[2][3];
   beam->algo        = 10; //d-phi
   return true;
}


bool accSvxBeamWidth(beamline *beam,int div)
{
  
  if( ! gROOT->FindObject("svxxminus42")) {
    beam->widthquality1 = false;
    beam->widthquality2 = false;
    return false;
  }  
  
  
  float z[ZNUM] = {-42,-40,-38,-36,-34,-32,-30,-28,-26,-24,-22,-10,-8,-6,-4,-2,
                 0,2,4,6,8,10,22,24,26,28,30,32,34,36,38,40,42}; 
  float xpos[ZNUM], ypos[ZNUM];
  
  for(int i=0; i<ZNUM; ++i)
  {
    std::ostringstream sx,sy;
    sx << "svxx";
    sy << "svxy";
    if(z[i] < 0 )
      {
    	sx << "minus";
    	sy << "minus";
      }
    sx << abs(z[i]) << "Div" << div << std::ends;
    sy << abs(z[i]) << "Div" << div << std::ends;
    
    
    xpos[i] = (beam->slopex)*z[i] + beam->x0;
    ypos[i] = (beam->slopey)*z[i] + beam->y0;
    
  
    for(int bin=0; bin < 400; ++bin)
    {
      int xentries, yentries;
      float bincentrex, bincentrey;
      
      bincentrex = ((TH1F*)gROOT->FindObject(sx.str().c_str()))->GetBinCenter(bin);
      bincentrey = ((TH1F*)gROOT->FindObject(sy.str().c_str()))->GetBinCenter(bin);
      xentries = ((TH1F*)gROOT->FindObject(sx.str().c_str()))->GetBinContent(bin);
      yentries = ((TH1F*)gROOT->FindObject(sy.str().c_str()))->GetBinContent(bin);
      
      for (int j=0; j<xentries; ++j) _svxxwidth[i]->Fill(bincentrex - xpos[i]);
      
      for (int j=0; j<yentries; ++j) _svxywidth[i]->Fill(bincentrey - ypos[i]);
    
    
    }// end for bin<400
  
  
  }//end for i<ZNUM
  

  return true;  
    
}// end accSvxBeamWidth


bool fitSvxBeamWidth(beamline *beam, float scalex, float scaley, float scalexerr, float scaleyerr)
{
  float z[ZNUM] = {-42,-40,-38,-36,-34,-32,-30,-28,-26,-24,-22,-10,-8,-6,-4,-2,
                 0,2,4,6,8,10,22,24,26,28,30,32,34,36,38,40,42};
  float zerr[ZNUM] = {1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1};
  float wx[ZNUM],wxerr[ZNUM];
  float wy[ZNUM],wyerr[ZNUM];
  float wxcor[ZNUM], wxcorh[ZNUM], wxcorl[ZNUM], wxcormean[ZNUM], wxcorsd[ZNUM], wxcorerr[ZNUM];
  float wycor[ZNUM], wycorh[ZNUM], wycorl[ZNUM], wycormean[ZNUM], wycorsd[ZNUM], wycorerr[ZNUM];
  float sigmasqx[ZNUM],sigmasqy[ZNUM];
  
  
  float zfitchisquarex = 0.;
  float zfitchisquarey = 0.;
  
  
  if( ! gROOT->FindObject("svxsigmasqxminus42")) {
    
    beam->widthquality1 = 0;
    beam->widthquality2 = 0;
    return false;
  }  
  
  
  
      
  for(int i = 0 ; i < ZNUM ; ++i)
  {
        
    if(_svxxwidth[i]->GetEntries() > 160 && _svxywidth[i]->GetEntries() > 160)
    {
      _svxxwidth[i]->Fit("gaus","IQ");
      _svxywidth[i]->Fit("gaus","IQ");
    
      if(_svxxwidth[i]->GetFunction("gaus")->GetChisquare() != 0 && 
         _svxywidth[i]->GetFunction("gaus")->GetChisquare() != 0)
      {
        wx[i] = _svxxwidth[i]->GetFunction("gaus")->GetParameter(2);
        wxerr[i] = _svxxwidth[i]->GetFunction("gaus")->GetParError(2);
	
	wy[i] = _svxywidth[i]->GetFunction("gaus")->GetParameter(2);
        wyerr[i] = _svxywidth[i]->GetFunction("gaus")->GetParError(2);
        
      }
      else 
      {
        beam->widthquality1 = 0;
        std::cout << "Chi square exactly zero, no width fit input" << std::endl;       
      }
    }
    else 
    {
      beam->widthquality1 = 0;	
      std::cout << "not enough entries in histogram, no width fit input" << std::endl;
    }  

  } // for i<ZNUM     
  
  
  for(int i=0; i<ZNUM; ++i)
  { 
    std::ostringstream sx,sy;
    sx << "svxsigmasqx";
    sy << "svxsigmasqy";
    if(z[i] < 0)
    {
      sx << "minus";
      sy << "minus";
    }
    sx << abs(z[i]) << std::ends;
    sy << abs(z[i]) << std::ends;
    sigmasqx[i] = ((TH1F*)gROOT->FindObject(sx.str().c_str()))->GetMean();
    sigmasqy[i] = ((TH1F*)gROOT->FindObject(sy.str().c_str()))->GetMean();
  } // for i<ZNUM     
      
    
  for(int i = 0; i < ZNUM; i++) 
  {
    if (wx[i] > (scalex+scalexerr)*sqrt(sigmasqx[i])) 
    {
      wxcor[i] = sqrt(wx[i]*wx[i] - scalex*scalex*sigmasqx[i]);
      wxcorh[i] = sqrt(wx[i]*wx[i] - (scalex+scalexerr)*(scalex+scalexerr)*sigmasqx[i]);
      wxcorl[i] = sqrt(wx[i]*wx[i] - (scalex-scalexerr)*(scalex-scalexerr)*sigmasqx[i]);
      wxcormean[i] = (wxcor[i]+wxcorh[i]+wxcorl[i])/3;
      wxcorsd[i] = sqrt(((wxcor[i]-wxcormean[i])*(wxcor[i]-wxcormean[i])+
    		      (wxcorh[i]-wxcormean[i])*(wxcorh[i]-wxcormean[i])+
    		      (wxcorl[i]-wxcormean[i])*(wxcorl[i]-wxcormean[i]))/2);
      wxcorerr[i] = sqrt(wxcorsd[i]*wxcorsd[i]+wxerr[i]*wxerr[i]);
      
    }  
    else 
    {
      if (wx[i] > scalex*sqrt(sigmasqx[i]))
      {
        wxcor[i] = sqrt(wx[i]*wx[i] - scalex*scalex*sigmasqx[i]);
        wxcorerr[i] = wxerr[i]*2;
      }
      else
      {
        wxcor[i] = 0.;
        wxcorerr[i] = wxerr[i]*2;
      }
      	
    }

    if (wy[i] > (scaley+scaleyerr)*sqrt(sigmasqy[i])) 
    {
      wycor[i] = sqrt(wy[i]*wy[i] - scaley*scaley*sigmasqy[i]);
      wycorh[i] = sqrt(wy[i]*wy[i] - (scaley+scaleyerr)*(scaley+scaleyerr)*sigmasqy[i]);
      wycorl[i] = sqrt(wy[i]*wy[i] - (scaley-scaleyerr)*(scaley-scaleyerr)*sigmasqy[i]);
      wycormean[i] = (wycor[i]+wycorh[i]+wycorl[i])/3;
      wycorsd[i] = sqrt(((wycor[i]-wycormean[i])*(wycor[i]-wycormean[i])+
    		      (wycorh[i]-wycormean[i])*(wycorh[i]-wycormean[i])+
    		      (wycorl[i]-wycormean[i])*(wycorl[i]-wycormean[i]))/2);
      wycorerr[i] = sqrt(wycorsd[i]*wycorsd[i]+wyerr[i]*wyerr[i]);
    }  
    else 
    {
      if (wy[i] > scaley*sqrt(sigmasqy[i]))
      {
        wycor[i] = sqrt(wy[i]*wy[i] - scaley*scaley*sigmasqy[i]);
        wycorerr[i] = wyerr[i]*2;
      }
      else
      {
        wycor[i] = 0.;
        wycorerr[i] = wyerr[i]*2;
      }	
     
    }
  }// end for i < ZNUM
  
  
  //computing weighted average width in x
  float beamwidthx = 0;
  float numbwx = 0; 
  float denbwx = 0;
  for (int i = 0; i < ZNUM; i++)
  {
    numbwx = numbwx + wxcor[i]/(wxcorerr[i]*wxcorerr[i]);
    denbwx = denbwx + 1/(wxcorerr[i]*wxcorerr[i]);
  }
  beamwidthx = numbwx/denbwx;
  if(beamwidthx > 0.005 || beamwidthx < 0.0010)
    {
      beam->widthquality2 = 0;
      std::cout << "bad range for avg beam width x, no avg width input" << std::endl;
      return true;
    }
  beam->widthx = beamwidthx;
  
  //computing weighted average width in y
  float beamwidthy = 0;
  float numbwy = 0; 
  float denbwy = 0;
  for (int i = 0; i < ZNUM; i++)
  {
    numbwy = numbwy + wycor[i]/(wycorerr[i]*wycorerr[i]);
    denbwy = denbwy + 1/(wycorerr[i]*wycorerr[i]);
  }
  beamwidthy = numbwy/denbwy;  
  if(beamwidthy > 0.005 || beamwidthy < 0.0010)
    {
      beam->widthquality2 = 0;
      std::cout << "bad range for avg beam width y, no avg width input" << std::endl;
      return true;
    }
  beam->widthy = beamwidthy;
  
  
  //getting beam width parameters in x
  TGraphErrors* gx = new TGraphErrors(ZNUM,z,wxcor,zerr,wxcorerr);
  TF1 *fx = new TF1("fx","sqrt([0]*([1]+(x-[2])*(x-[2])/[1]))",-50.,50.);
  fx->SetParameters(0.0000002,45.,0.);
  gx->Fit("fx","IQ");
  beam->emittxfit = gx->GetFunction("fx")->GetParameter(0);
  beam->bstarxfit = gx->GetFunction("fx")->GetParameter(1);
  beam->zzeroxfit = gx->GetFunction("fx")->GetParameter(2);
  zfitchisquarex = gx->GetFunction("fx")->GetChisquare();
  
  //getting beam width parameters in y
  TGraphErrors* gy = new TGraphErrors(ZNUM,z,wycor,zerr,wycorerr);
  TF1 *fy = new TF1("fy","sqrt([0]*([1]+(x-[2])*(x-[2])/[1]))",-50.,50.);
  fy->SetParameters(0.0000002,45.,0.);
  gy->Fit("fy","IQ");
  beam->emittyfit = gy->GetFunction("fy")->GetParameter(0);
  beam->bstaryfit = gy->GetFunction("fy")->GetParameter(1);
  beam->zzeroyfit = gy->GetFunction("fy")->GetParameter(2);
  zfitchisquarey = gy->GetFunction("fy")->GetChisquare();
  
  if (zfitchisquarex > 160 || zfitchisquarey > 160) 
  {
    beam->widthquality1 = 0;
    std::cout << "large chi square, no width fit input" << std::endl;
    std::cout << "chi square x is " << zfitchisquarex << std::endl;
    std::cout << "chi square y is " << zfitchisquarey << std::endl;
  }  
  if (beam->emittxfit > 0.0000003 || beam->emittxfit < 0.00000005 || beam->emittyfit > 0.0000003 || beam->emittyfit < 0.00000005) 
  {
    beam->widthquality1 = 0;
    std::cout << "emittance out of range, no width fit input" << std::endl;
    std::cout << "emit x is " << beam->emittxfit << std::endl;
    std::cout << "emit y is " << beam->emittyfit << std::endl;
  }
  if (beam->bstarxfit > 50 || beam->bstarxfit < 20 || beam->bstaryfit > 50 || beam->bstaryfit < 20) 
  {
    beam->widthquality1 = 0;
    std::cout << "amplitude function out of range, no width fit input" << std::endl;
    std::cout << "beta* x is " << beam->bstarxfit << std::endl;
    std::cout << "beta* y is " << beam->bstaryfit << std::endl;
  }  
  if (beam->zzeroxfit > 30 || beam->zzeroxfit < -5 || beam->zzeroyfit > 10 || beam->zzeroyfit < -25) 
  {
    beam->widthquality1 = 0;
    std::cout << "beam foci out of range, no width fit input" << std::endl;
    std::cout << "z0 x is " << beam->zzeroxfit << std::endl;
    std::cout << "z0 y is " << beam->zzeroyfit << std::endl;
  }  
 
  
  if (beam->widthquality1 != 0) std::cout << "criteria passed" << std::endl;
  return true;     
 
}//end fitSvxBeamWidth


bool fitBeam(beamline *beam, TProfile* hxz, TProfile* hyz, TH1F* hz)
{ 
  //std::cout << "fitBeam(beamline *beam, TProfile* hxz, TProfile* hyz, TH1F* hz)" << std::endl;
  TF1 beamfit("beamfit","[0] + x * [1]");
  TF1 zvtxdist("zvtxdist","[0] * exp(- (x- [1]) * (x - [1]) / 2 / [2] /[2]) / ( 1 + (x - [4]) * (x - [4])/ [3] / [3])");                     
  
  if(hxz->GetEntries() <  1000) return false;
 

  hxz->Fit("beamfit","IQ","",-40,40);
  hxz->Fit("beamfit","EQ","",-40,40);
  beam->x0          = beamfit.GetParameter(0);
  beam->slopex      = beamfit.GetParameter(1);
  beam->sigmax0     = beamfit.GetParError(0)*beamfit.GetParError(0);
  beam->sigmaslopex = beamfit.GetParError(1)*beamfit.GetParError(1);
      
  hyz->Fit("beamfit","IQ","",-40,40);
  hyz->Fit("beamfit","EQ","",-40,40);
  beam->y0          = beamfit.GetParameter(0);
  beam->slopey      = beamfit.GetParameter(1);
  beam->sigmay0     = beamfit.GetParError(0)*beamfit.GetParError(0);
  beam->sigmaslopey = beamfit.GetParError(1)*beamfit.GetParError(1);
 
  
  hz->Fit("gaus","IQ","",-40,40);
  beam->beamz  = hz->GetFunction("gaus")->GetParameter(1);
  beam->widthz = hz->GetFunction("gaus")->GetParameter(2);

  zvtxdist.SetParameters(hz->GetMaximum(),beam->beamz,37,35,beam->beamz);
  hz->Fit("zvtxdist","IQ","",-50,50);
  //get fit results minimzfit, widthzfit, bstarzFit, zzerozfit;
  beam->widthzfit = zvtxdist.GetParameter(2);
  beam->bstarzfit = zvtxdist.GetParameter(3);
  beam->zzerozfit = zvtxdist.GetParameter(4);
  beam->minimzfit = zvtxdist.GetParameter(1);
  
  beam->statistics0 = hz->GetEntries();
  beam->algo = 20;//vertices

  return true;
}

bool fitBeam(beamline *beam, TH2F* hxz, TH2F* hyz, TH1F* hz)
{ 
  //std::cout << "fitbeam(beamline *beam, TH2F* hxz, TH2F* hyz, TH1F* hz)" << std::endl;
  TProfile *htempx  = hxz->ProfileX("htempx",0,hxz->GetNbinsY(),"e");
  TProfile *htempy  = hyz->ProfileX("htempy",0,hyz->GetNbinsY(),"e");
      
  bool result = fitBeam(beam,htempx,htempy,hz);
  
  delete htempx;
  delete htempy;
  
  return result;
}

//stores the beam lines into the database
bool fillSvxDB(beamline beamlines[MAX_DIV],int entries, int runnumber, globalalignment *align, std::string dbid,
	       const std::vector<CID>& cids, std::string productionversion)
{ 
  //get the run list key
  SvxBeamPositionContainer contDB;
  SvxBeamPosition_mgr iomSvxBeam(dbid,"SvxBeamPosition");
  if(iomSvxBeam.isValid()==false)
    {
      std::cout <<"FitBeam: SvxBeamPosition_mgr initialization failure" << std::endl;
      return  false;
    }                                                                

  RunListKey qk(runnumber,Version::generate); //just add a new entry
  qk.setDataStatus("GOOD");
  qk.newVersion();
  qk.setAlgorithm(productionversion);
  contDB.clear();
  //fill the container
  for(int i = 0 ; i < entries ; ++i)
    {
      beamline* beam = &beamlines[i];

      if( beam->statistics0 == 0 ) continue;
      std::ostringstream os;
      os << (long) i + 1 << " ";
      os << beam->x0 << " " << beam->y0 << " " << beam->slopex << " " 
	 << beam->slopey << " ";
      if(beam->widthquality2==1 && i==0) {
        os << beam->widthx << " " << beam->widthy << " " << (double)0 << " ";
      }	
      else os << 0 << " " << 0 << " " << (double)0 << " ";
      
      os << beam->beamz << " " << beam->widthz << " ";
      os << 0 << " " << 0 << " " << 0 << " " << 0 << " ";
      os << 0 << " " << 0 << " " << 0 << " ";
      os << 0 << " " << 0 << " "; 
      if((std::fabs(beam->minimzfit-beam->beamz) > 15) || (beam->bstarzfit>100) 
	 || (beam->bstarzfit<0) || (beam->widthzfit < 0))
	{
	  // ignore fit results
	  os << 0 << " " << 0 << " " << 0 << " " << 0 << " ";
	  os << 0 << " " << 0 << " " << 0 << " " << 0 << " ";
	}
      else
	{
	  os << beam->minimzfit << " " << beam->widthzfit << " " 
	     << beam->bstarzfit << " " << beam->zzerozfit << " ";
	  os << 0 << " " << 0 << " " << 0 << " "  << 0 << " ";
	}
      // the transverse emittance and beta* from the fit
      // to the beam width as a function of z, as in CDF 4066
      if(beam->widthquality1==1 && i==0) {
        std::cout << "writing x and y width parameters" << std::endl;
        os << beam->emittxfit << " " << beam->bstarxfit << " " << beam->zzeroxfit << " ";
        os << beam->emittyfit << " " << beam->bstaryfit << " " << beam->zzeroyfit << " ";
      }
      else {
        os << 0 << " " << 0 << " " << 0 << " ";
        os << 0 << " " << 0 << " " << 0 << " ";
      }	
      os << 0 << " " << 0 << " " << 0 << " ";
      os << 0 << " " << 0 << " " << 0 << " ";               
      //cov
      os << beam->sigmax0 << " " << beam->sxy << " " << beam->sxsx << " " << beam->sxsy << " ";
      os << beam->sigmay0 << " " << beam->sysx << " " << beam->sysy << " ";
      os << beam->sigmaslopex << " " << beam->ssxsy << " ";
      os << beam->sigmaslopey << " ";
      


      //statistics
      os << beam->statistics0 << " " << 0 << " ";
      //flags
      long flag0,flag1;
      long hist = beam->history;
      long algo = beam->algo;
      long srce = 10;// 5 online, 10 offline
      flag0 = (algo << 16) + hist;
      flag1 =  srce;
      os << flag0 << " " << flag1 << " ";
      //alignment
      os << align->tx << " " << align->ty << " " << align->slopex << " "<< align->slopey << " ";
      
      os << std::ends;
      std::cout << "SVX entry:" << i << std::endl;
      std::cout << os.str() << std::endl;
  
      SvxBeamPosition svxbeampos;
      std::istringstream is(os.str());    
      is >> svxbeampos;
  
      contDB.push_back(svxbeampos);
    }
  if( contDB.empty() ) return false;
  //store result
  for(std::vector<CID>::const_iterator iter = cids.begin() ; iter != cids.end(); ++iter)
    contDB.addParent(RunListKey(*iter));
  if(iomSvxBeam.put(qk,contDB) != Result::success)
    {
      std::cout << "FitBeam: Put failed for SvxBeamPosition Run = " << runnumber
		<< " in database " << dbid << "." << std::endl; 
      return false;
    }
  return true;
}

bool fillCotDB(beamline beamlines[MAX_DIV],int entries, int runnumber, std::string dbid, std::string productionversion)
{ 
  CotBeamPositionContainer contDB;
  CotBeamPosition_mgr iomCotBeam(dbid,"CotBeamPosition");
  if(iomCotBeam.isValid()==false)
    {
      std::cout <<"FitBeam: CotBeamPosition_mgr initialization failure" << std::endl;
      return  false;
    }                                                                
  RunListKey qk(runnumber,Version::generate); //just add a new entry
  qk.setDataStatus("GOOD");
  qk.newVersion(); 
  qk.setAlgorithm(productionversion);
  contDB.clear();

  for( int i = 0 ; i < entries ; ++i )
    {
      beamline* beam = &beamlines[i];
      
      if( beam->statistics0 == 0 ) continue;
      
      std::ostringstream os;
      
      os << (long) i + 1 << " ";
      os << beam->x0 << " " << beam->y0 << " " << beam->slopex << " " 
	 << beam->slopey << " ";
      os << (double)0 << " " << (double)0 << " " << (double)0 << " ";
      os << beam->beamz << " " << beam->widthz << " ";
      os << 0 << " " << 0 << " " << 0 << " " << 0 << " ";
      os << 0 << " " << 0 << " " << 0 << " ";
      os << 0 << " " << 0 << " ";
      if((std::fabs(beam->minimzfit-beam->beamz) > 15) || (beam->bstarzfit>100) 
	 || (beam->bstarzfit<0) || (beam->widthzfit < 0))
	{
	  // ignore fit results
	  os << 0 << " " << 0 << " " << 0 << " " << 0 << " ";
	  os << 0 << " " << 0 << " " << 0 << " " << 0 << " ";
	}
      else
	{
	  os << beam->minimzfit << " " << beam->widthzfit << " " 
	     << beam->bstarzfit << " " << beam->zzerozfit << " ";
	  os << 0 << " " << 0 << " " << 0 << " "  << 0 << " ";
	}
      os << beam->sigmax0 << " " << beam->sxy << " " << beam->sxsx << " " << beam->sxsy << " ";
      os << beam->sigmay0 << " " << beam->sysx << " " << beam->sysy << " ";
      os << beam->sigmaslopex << " " << beam->ssxsy << " ";
      os << beam->sigmaslopey << " ";
      //statistics
      os << beam->statistics0 << " " << 0 << " ";
      //flags
      long flag0,flag1;
      long hist = beam->history;
      long algo = beam->algo;
      long srce = 10;// 5 online, 10 offline
      flag0 = (algo << 16) + hist;
      flag1 =  srce;
      os << flag0 << " " << flag1 << " ";
      //alignment
      os << 0 << " " << 0 << " " << 0 << " "<< 0 << " ";
      
      os << std::ends;

      std::cout << "Cot entry:" << i << std::endl;
      std::cout << os.str() << std::endl;
      
      std::istringstream is(os.str());    
      CotBeamPosition cotbeampos;
      is >> cotbeampos;
      contDB.push_back(cotbeampos);
    } 
  if( contDB.empty() ) return false;
  if(iomCotBeam.put(qk,contDB) != Result::success)
    {
      std::cout << "FitBeam: Put failed for CotBeamPosition Run = " << runnumber
		<< " in database " << dbid << "." << std::endl; 
      return false;
    }

  return true;
}

int main(int argc, char* argv[])
{
  if((argc != 16) && (argc != 17)) 
    {
      std::cout << "USAGE: FitBeam <root hist file> <runnumber> <history> <dbid> <aligndbid> <alignmentrun> <alignment version> <alignment status> <production version> <vertices/dphi/both> <width/nowidth> <scale factor x> <scale factor y> <scale factor error x> <scale factor error y> [nocot/nosvx]" << std::endl;
      return 1;
    }
  TROOT root("Root","ROOT for FitBeam");
  gROOT->SetBatch(true);
  TFile* file = new TFile(argv[1]);
  beamline beam[MAX_DIV];
  int entries = 0;
  int history = atoi(argv[3]);
  int runnumber = atoi(argv[2]);
  std::string dbid = argv[4];
  globalalignment align;
  std::string alignmentdbid = argv[5];
  int alignmentrun = atoi(argv[6]);
  int alignmentversion = atoi(argv[7]);
  std::string alignmentstatus = argv[8];
  std::string productionversion = argv[9];
  std::string option = argv[10];
  bool makedphi = (option == "dphi") || (option == "both");
  bool makevertices =  (option == "vertices") || (option == "both"); 
  std::string option2 = argv[11];
  bool inputwidth = (option2 == "width");
  float scalex = atof(argv[12]);
  float scaley = atof(argv[13]);
  float scalexerr = atof(argv[14]);
  float scaleyerr = atof(argv[15]);
  
  bool makesvx = true;
  bool makecot = true;
  if(argc == 17)
    {
      std::string option = argv[16];
      if(option == "nosvx") makesvx = false;
      else if(option == "nocot") makecot = false;
    }

 
  for(int i = 0 ; i < MAX_DIV ; ++i) beam[i].history = history;
  
  
  std::vector<CID> alignmentCIDs = getAlignmentCIDs(alignmentdbid,alignmentrun,alignmentversion,
						    alignmentstatus);
  		    
  if( alignmentCIDs.size() != NUM_ALIGNMENT_NAMES)
    {
      std::cout << "FitBeam: Failed to get all the CIDs for the alignment." << std::endl;
      return 1;
    }
   
//   std::cout << "alignment CIDs: ";
//   for(std::vector<CID>::const_iterator iter = alignmentCIDs.begin() ; 
//       iter != alignmentCIDs.end() ; ++iter)
//     std::cout << *iter << " , ";
//   std::cout << std::endl;
  if(! readGlobalAlignment(alignmentdbid,alignmentrun,alignmentversion,alignmentstatus,&align))
    {
      std::cout << "FitBeam: Failed to read the global alignment." << std::endl;
      return 1;
    }
//   std::cout << "global alignment:" << align.tx << " , " << align.ty << " , " << align.slopex 
// 	       << " , " << align.slopey << " , " << std::endl;
  
  if(! file)
    {
      std::cout << "FitBeam: could not open ROOT file." << std::endl;
      return 1;
    }
  file->cd("FitBeamModule");

  bool result = false;
  if(makevertices)
    {
      if(makesvx)
	{
	  std::cout << "FitBeam: fitting SVX beam line...   ";
	  
	  for(int i=0; i<ZNUM; ++i)
	  {
	    _svxxwidth[i] = new TH1F("svxxwidth", "beam width in x", 400, -0.2, 0.2);
	    _svxywidth[i] = new TH1F("svxywidth", "beam width in y", 400, -0.2, 0.2);
	    
	  }//end for i<ZNUM
	  
	  
	  for(int i = 0; i < MAX_DIV ; ++i) beam[i].reset();
	  entries = 0;
	  if( gROOT->FindObject("svxxz")->InheritsFrom("TProfile") )
	    {
	      result = fitBeam(&beam[entries],(TProfile*)gROOT->FindObject("svxxz"),
			       (TProfile*)gROOT->FindObject("svxyz"),
			       (TH1F*)gROOT->FindObject("svxz"));
	      
	      ++entries;
	      
	      //fit divisions
	      while( entries < MAX_DIV )
		{
		  std::ostringstream namexz,nameyz,namez;
		  namexz << "svxxzDiv" << entries << std::ends;
		  nameyz << "svxyzDiv" << entries << std::ends;
		  namez << "svxzDiv" << entries << std::ends;
		  TH1F* hsvxz = (TH1F*)gROOT->FindObject(namez.str().c_str());
		  if(! hsvxz) break;
		  if(hsvxz->GetEntries() == 0) break;
		  bool res = fitBeam(&beam[entries],(TProfile*)gROOT->FindObject(namexz.str().c_str()),
				     (TProfile*)gROOT->FindObject(nameyz.str().c_str()),hsvxz);
		  if( ! res ) beam[entries].reset();
		  else if (inputwidth)
		  {
		    file->cd("BeamWidthModule");
		    accSvxBeamWidth(&beam[entries],entries);
		    file->cd("FitBeamModule");
		  }  
		  ++entries;
		}
	    }
	  else
	    {
	      result = fitBeam(&beam[entries],(TH2F*)gROOT->FindObject("svxxz"),
			       (TH2F*)gROOT->FindObject("svxyz"),
			       (TH1F*)gROOT->FindObject("svxz"));
	      
	      ++entries;
	      
	    }
	  if( result ) 
	    {
	      
	      if (inputwidth)
	      {
	        file->cd("BeamWidthModule");
		fitSvxBeamWidth(&beam[0], scalex, scaley, scalexerr, scaleyerr); 
		file->cd("FitBeamModule");
	      }
	      
	      std::cout << "done." << std::endl;
	      std::cout << "FitBeam: making SVX DB entry...   " << std::endl;
	      if( fillSvxDB(beam,entries,runnumber,&align,dbid,alignmentCIDs,productionversion) )
		std::cout << "done." << std::endl;
	      else 
		std::cout << "failed." << std::endl;
	    }
	  else std::cout << "failed." <<  std::endl;
	}
      if(makecot)
	{
	  std::cout << "FitBeam: fitting COT beam line...   "; 
	  for(int i = 0; i < MAX_DIV ; ++i) beam[i].reset();
	  entries = 0;
	  if( gROOT->FindObject("cotxz")->InheritsFrom("TProfile") )
	    {
	      result = fitBeam(&beam[entries],(TProfile*)gROOT->FindObject("cotxz"),
			       (TProfile*)gROOT->FindObject("cotyz"),
			       (TH1F*)gROOT->FindObject("cotz"));
	      ++entries; 
	      //fit divisions
	      while( entries < MAX_DIV )
		{
		  std::ostringstream namexz,nameyz,namez;
		  namexz << "cotxzDiv" << entries << std::ends;
		  nameyz << "cotyzDiv" << entries << std::ends;
		  namez << "cotzDiv" << entries << std::ends;
		  TH1F* hcotz = (TH1F*)gROOT->FindObject(namez.str().c_str());
		  if(! hcotz) break;
		  if(hcotz->GetEntries() == 0) break;
		  bool res = fitBeam(&beam[entries],(TProfile*)gROOT->FindObject(namexz.str().c_str()),
				     (TProfile*)gROOT->FindObject(nameyz.str().c_str()),hcotz);
		  if(( ! res  ) || (hcotz->GetEntries() < 5000) )  
		    beam[entries].reset();
		  ++entries;
		}
	    }
	  else
	    {
	      result = fitBeam(&beam[entries],(TH2F*)gROOT->FindObject("cotxz"),
			       (TH2F*)gROOT->FindObject("cotyz"),
			       (TH1F*)gROOT->FindObject("cotz"));
	      ++entries; 
	    }
	  if( (((TH1F*)gROOT->FindObject("cotz"))->GetEntries() > 5000) && result )
	    {
	      std::cout << "done." << std::endl;
	      std::cout << "FitBeam: making COT DB entry...   " << std::endl;
	      if( fillCotDB(beam,entries,runnumber,dbid,productionversion) )
		std::cout << "done" << std::endl;	
	      else 
		std::cout << "failed." << std::endl;
	    }
	  else std::cout << "failed." <<  std::endl;
	}
    }
  if(makedphi)
    {
      if(makesvx)
	{
	  std::cout << "FitBeam: fitting SVX beam line using d-phi...   ";
	  for(int i = 0; i < MAX_DIV ; ++i) beam[i].reset();
	  entries = 0;
	  if( gROOT->FindObject("svxdphi") && fitDPhi(&beam[entries],(TH1F*)gROOT->FindObject("svxdphi")) )
	    {
	      //make divisions
	      ++entries;
	      while( entries < MAX_DIV )
		{
		  std::ostringstream namedphi;
		  namedphi << "svxdphiDiv" << entries << std::ends;
		  TH1F* hdphi = (TH1F*)gROOT->FindObject(namedphi.str().c_str());
		  if(!  hdphi) break;
		  if( hdphi->GetBinContent(hdphi->FindBin(16)) == 0 ) break;
		  if( ! fitDPhi(&beam[entries],hdphi) ) beam[entries].reset();
		  ++entries;
		}
	      std::cout << "done." << std::endl;
	      std::cout << "FitBeam: making SVX DB entry for d-phi result...   " << std::endl;
	      if( fillSvxDB(beam,entries,runnumber,&align,dbid,alignmentCIDs,productionversion) )
		std::cout << "done." << std::endl;
	      else 
		std::cout << "failed." << std::endl;
	    }
	  else 
	    std::cout << "failed." <<  std::endl;
	}
      if(makecot)
	{
	  std::cout << "FitBeam: fitting COT beam line using d-phi...   ";
	  for(int i = 0; i < MAX_DIV ; ++i) beam[i].reset();
	  entries = 0;
	  if( gROOT->FindObject("cotdphi") && fitDPhi(&beam[entries],(TH1F*)gROOT->FindObject("cotdphi")) )
	    {
	      //make divisions
	      ++entries;
	      while( entries < MAX_DIV )
		{
		  std::ostringstream namedphi;
		  namedphi << "cotdphiDiv" << entries << std::ends;
		  TH1F* hdphi = (TH1F*)gROOT->FindObject(namedphi.str().c_str());
		  if(!  hdphi) break;
		  if( hdphi->GetBinContent(hdphi->FindBin(16)) == 0 ) break;
		  if( ! fitDPhi(&beam[entries],hdphi) ) beam[entries].reset();
		  ++entries;
		}
	      std::cout << "done." << std::endl;
	      std::cout << "FitBeam: making COT DB entry for d-phi result...   " << std::endl;
	      if( fillCotDB(beam,entries,runnumber,dbid,productionversion) )
		std::cout << "done" << std::endl;	
	      else 
		std::cout << "failed." << std::endl;
	    }
	  else 
	    std::cout << "failed." <<  std::endl;
	}
    }
  file->Close();
  std::cout << "FitBeam: finished." << std::endl;
}