#include "Stntuple/photon/TDiPhoMat.hh"

using std::cout;
using std::endl;


TDiPhoMat::TDiPhoMat():m4(TVector(4)),ebew(TVector(4)),ebeww(TVector(4)) {
  Clear();
}

TDiPhoMat::~TDiPhoMat() {}

void TDiPhoMat::Clear() {
  m4.Zero();
  ebew.Zero();
  ebeww.Zero();
  nonfiducial = 0;
}


void TDiPhoMat::Print() {

  double total = m4(0)+m4(1)+m4(2)+m4(3);
  double scale = double(total+nonfiducial)/total;

  cout << "total: " << total << "   nonfid: " << nonfiducial 
       << "  scale: " << scale << endl;
  cout << " meas    ans    sig    frac     sig     scale     sig  " << endl;

  for(int i=0; i<4; i++) {
    cout << m4(i) << "     " << ebew(i) << "     " << sqrt(ebeww(i)) 
	 << "     " << ebew(i)/total << "     " << sqrt(ebeww(i))/total 
	 << "     " << scale*ebew(i) << "     " << scale*sqrt(ebeww(i)) 
	 << endl;
  }

}

double TDiPhoMat::AddEvent(double es1, double eb1, double es2, double eb2,
			bool pass1, bool pass2) {
  if(pass1) {
    if(pass2) {
      m4(0)++;
    } else {
      m4(1)++;
    }
  } else {
    if(pass2) {
      m4(2)++;
    } else {
      m4(3)++;
    }
  }

  TMatrix temp4(4,4);

  // first index pp,pf,fp,ff, second index: gg,gj,jg,jj
  temp4(0,0) = es1*es2;              // gg
  temp4(1,0) = es1*(1.0-es2);      
  temp4(2,0) = (1.0-es1)*es2;      
  temp4(3,0) = (1.0-es1)*(1.0-es2);
  temp4(0,1) = es1*eb2;            
  temp4(1,1) = es1*(1.0-eb2);      
  temp4(2,1) = (1.0-es1)*eb2;      
  temp4(3,1) = (1.0-es1)*(1.0-eb2);
  temp4(0,2) = eb1*es2;            
  temp4(1,2) = eb1*(1.0-es2);      
  temp4(2,2) = (1.0-eb1)*es2;      
  temp4(3,2) = (1.0-eb1)*(1.0-es2);
  temp4(0,3) = eb1*eb2;            
  temp4(1,3) = eb1*(1.0-eb2);      
  temp4(2,3) = (1.0-eb1)*eb2;      
  temp4(3,3) = (1.0-eb1)*(1.0-eb2);


  TMatrix temp4i(TMatrix::kInverted,temp4);

  int ind;
  if(pass1) {
    if(pass2) {
      ind = 0;
    } else {
      ind = 1;
    }
  } else {
    if(pass2) {
      ind = 2;
    } else {
      ind = 3;
    }
  }
  for(int i=0; i<4; i++) {
    ebew(i) += temp4i(i,ind);
    ebeww(i) += temp4i(i,ind)*temp4i(i,ind);
  }
  return temp4i(0,ind);
}