// Include files
#include <iostream>
#include <string>
#include <vector>
#include "time.h"
#include "boost/lexical_cast.hpp"
#include <fstream> 
#include <stdlib.h>
#include "dim/dis.hxx"
#include "TSystem.h"
#include <TROOT.h>
#include <TH1D.h>
#include <TH2D.h>
#include <TProfile.h>
#include <TFile.h>
#include <TBenchmark.h>
#include <TROOT.h>
#include <TRandom.h>
#include <TObject.h>
#include <TGraph.h>
#include <TMath.h>
#include "TImage.h"
#include "TCanvas.h"
#include "TArrayD.h"
#include "TColor.h"
#include "TAttImage.h"
#include "TEnv.h"
#include "TASImage.h"
#include "TSpectrum.h"
#include "TString.h"
#include "TF1.h"
#include "TF2.h"

TArrayD xVal;
TArrayD yVal;
TArrayD Intensity;
TString date;

ofstream intensity;
ofstream writefile;
ofstream LogFile;

struct LAMSCommand{
  int debugmode;
  char filename[50];
};

struct ImgAnalysisResults{
  double refXCoord;
  double refYCoord;
  double refLumi;
  double reflXCoord;
  double reflYCoord;
  double reflLumi;
  //  double radSep;
  double imgLumi;
  int quality;
};

ImgAnalysisResults Results;
ImgAnalysisResults compute(struct LAMSCommand command);

int main()
{ 
  gROOT->Reset();
  //create DIM command from PVSS
  DimCommand PVSSCommand("RichLAMS/ImageFilename", "I:1;C");

  // create DIM services, one per camera
  ImgAnalysisResults result= {0,0,0,0,0,0,0,999};
  DimService LAMSResults( "RichLAMS/ImgAnalResults", "D:7;I:1", (void*)&result, sizeof(result) );  
  
  // start the DIM server
  DimServer::start("RichLAMS");

  void* voidCommandPointer;
  LAMSCommand command;
  // enter infinite loop waiting for commands
  while(1) {
  
    while (PVSSCommand.getNext() ) 
    {
      // command arrived
      voidCommandPointer = PVSSCommand.getData();
      LAMSCommand* comPointer = static_cast<LAMSCommand*>(voidCommandPointer);
      command = *comPointer;
      //      std::cout << "Filename " << command.filename << std::endl;

      result = compute(command);
      //   ImgAnalysisResults Results;
      // update the service realted to this camera
      if(Results.quality<10)LAMSResults.updateService();
    } 
    sleep(1);
  }
  return 0;
}


ImgAnalysisResults compute(struct LAMSCommand command){
    
  TString camName="";
  Int_t h = 0; 
  Int_t w = 0 ; 
  Int_t bufferSize= 0;
  bool toDraw = true;
  
  // get all the info from the filename
  TString filename= command.filename;
  TString copyfilename= filename;
  TString rich= filename;
  TString time= filename;
  TString exp= filename;
  TString gain= filename;
  TString date= filename;
  TString camera= filename;
  
  //get rich number
  rich.Remove(2,40);
  TString Rich = rich;
  Rich.Remove(0,1);
  //rich1
  if(Rich == "1"){
    camera.Remove(7,27);
    camera.Remove(0,2);
    date.Remove(16,27);
    date.Remove(0,8);
    time.Remove(23,27);
    time.Remove(0,17);
    // copyfilename.Remove(23,4); 
  } 
  //rich2
  if(Rich == "2"){
    camera.Remove(7,40);
    camera.Remove(0,2);
    date.Remove(16,40);
    date.Remove(0,8);
    time.Remove(23,40);
    time.Remove(0,17);
    exp.Remove(30,40);
    exp.Remove(0,24);
    //    gain.Remove(36,40);
    //   if(camera.Contains("5")) gain.Remove(36,40);
    if(!(camera.Contains("5"))) gain.Remove(36,40);
    gain.Remove(0,30);
    //   if(camera.Contains("5")) copyfilename.Remove(36,40);
    // else copyfilename.Remove(36,4); 
  } 
  //get camera value without "cam"
  TString cam = camera;
  int quality = -1;
  cam.Remove(0,3); 
  if(cam.Contains("0")) cam.Remove(0,1);
  TString camID = cam;

  TString path="/group/rich/LAMS/RICH";
  path.Append(Rich);
  path.Append("/");
  path.Append(date);
  path.Append("/");

  //Output data files
  TString tobewritten= path;
  tobewritten.Append("data");
  if(command.debugmode == 0){tobewritten.Append("-R");}
  if(command.debugmode == 1){tobewritten.Append("-DebugR");}
  tobewritten.Append(Rich); 
  tobewritten.Append("-cam"); 
  tobewritten.Append(cam); 
  tobewritten.Append(".txt");
  writefile.open(tobewritten, std::ios::app); 
  writefile<<date<<" "<<time<<" ";

  TString LogFileName=path;
  LogFileName.Append("RICH");
  LogFileName.Append(Rich);
  if(command.debugmode == 0){LogFileName.Append("_AnalysisLogfile_");}
  if(command.debugmode == 1){LogFileName.Append("_AnalysisDebugLogfile_");} // this is the debug file
  LogFileName.Append(date); 
  LogFileName.Append(".txt");
  LogFile.open(LogFileName,std::ios::app);
  gSystem->RedirectOutput(LogFileName); 
  
  //convert
  //  TString toConvert ="convert ";
  //   toConvert.Append(path);
  //   toConvert.Append(filename);
  //   toConvert.Append(" ");
  //   toConvert.Append(path);
  //   toConvert.Append(copyfilename);
  //   toConvert.Append(".bmp");
  //   std::cout<<"the convert command "<< toConvert<<std::endl; 
  //  gSystem->RedirectOutput(ConvertOutput);
  // Int_t execCode = gSystem->Exec(toConvert);
  // std::cout<<"gSystem->execute returns  "<<execCode<<std::endl;
  //  gSystem->RedirectOutput(LogFileName); 
  //  ifstream test(ConvertOutput, std::ios::binary);
  //  test.seekg(0, std::ios::end);
  //  if (int(test.tellg()) !=0)   quality=4;
  //compose .bmp Filename
  TString Filename = path;
  Filename.Append(copyfilename);
  std::cout<<"************ Start of log from this analysis ************"<<std::endl;
  std::cout <<"filename from PVSS:   "<<copyfilename<<std::endl; 
  if(Rich == "2") copyfilename.Remove(36,4);
  if(Rich == "1") copyfilename.Remove(23,4);
  std::cout<<"filename after remove  "<<copyfilename<<std::endl; 
  
  Filename.Append(".bmp");
  std::cout<<"Looking at the file ---> "<<Filename<<std::endl;
  if(command.debugmode == 1){
    std::cout  << "Data extracted from filename: rich "<< rich << " Camera " << camera << " Date " << date << " Time " << time << std::endl; 
  }
  if(Rich == "1")
    {
      bufferSize = 921600;
      h=480; w=640;
    } 
  else 
    if(Rich =="2")
      {
	bufferSize=3932160;
	h=1024; w=1280;           
	if(command.debugmode == 1){ std::cout  << " Exp " << exp << " Gain " << gain << std::endl; }
      }
  
  if(command.debugmode == 1){std::cout<<"Filename "<<Filename<<std::endl;}       
 
  //now analyze the image
  //TCanvas* c1 = new TCanvas("image","image",1200, 700);
  //if(toDraw) c1->Divide(4,3);
  
  TString oneOrtwo = Rich; 
  TString camStr = cam;
  Int_t xbin=0;
  Int_t ybin=0;
  Int_t zbin=0;
  Int_t bin=0;
  TArrayD xVal=TArrayD(10);
  TArrayD yVal=TArrayD(10);
  TArrayD Intensity=TArrayD(10);
  Int_t counter=0;
  Double_t xbuff=0;
  Double_t ybuff=0;
  Double_t cbuff=0;
  Double_t xvalue=0;
  Double_t yvalue=0;
  Int_t point1=0;
  Int_t point2=0;
  Double_t buffermax=0;
  Double_t threshold=0.5;
  Double_t thresholdII= 0.5;
  int radius1= 0;
  int radius2= 0;
  int deltaRad= 0;
  Double_t overallcourse=0;
  Double_t overallfine=0;
  int a=0;
  Double_t x1coarse;
  Double_t x1fine;
  Double_t y1coarse;
  Double_t y1fine;
  Double_t intensityC;
  Double_t intensityF;
  Int_t myno=0;
  Int_t max=0;
  Int_t n=0;
  Int_t nII=0;
  Int_t nIII=0;
  TString toDel="rm "; 
  
  TCanvas *c1 = new TCanvas("c1","c1",1200, 700);
  TH2D* total =new TH2D ("total", "total",w-2, 0.5, w-1.5, h-2 , 0.5, h-1.5);
  TH2D* ImageHistogram =new TH2D("ImageHistogram", "ImageHistogram", 
                                 w-2, 1-w/2 , w/2-1, h-2, 1-h/2, h/2-1);
  TH2D* BufferHistogram=new TH2D("BufferHistogram","BufferHistogram",w-2, 1-w/2 , w/2-1, h-2, 1-h/2, h/2-1); 
  TH2D* ImgHistRebinned = new TH2D ("ImgHistRebinned","ImgHistRebinned", w/10 -2, 10-w/2, w/2-10, h/10 -2 , 10- h/2, h/2 -10);
  TH2D* BuffHist = new TH2D ("BuffHist", "BuffHist",w/10-2, 10-w/2 , w/2-10, h/10-2, 10-h/2, h/2-10); 
  
  ifstream *myfile = new ifstream (Filename, std::ios::binary);
  //  if(quality==4) goto BadImage;
  char buff[bufferSize];//1310720-54]; 
  for(Int_t check=0; check<3; check++)
    {
      if(!myfile){ std::cout<<"trying to open "<<Filename<<"  for the "<<check+1<<"  time "<<std::endl; myfile->open(Filename,std::ios::binary);}
      else check=3;
      if (!myfile && check==2){ std::cout<<"error opening your file"<<std::endl;  quality=4;  goto BadImage;}
    }
  if(myfile)
    {
      myfile->seekg(54, std::ios::beg);
      myfile->read(buff, sizeof(buff));
      myfile->close();
    }
  
  if(toDraw) c1->Divide(4,3);
  
  for(Int_t g=0; g<bufferSize/3; g++)
    total->SetBinContent(g,0);
  
  for(Int_t f=1; f<=bufferSize; f++) // if(buff[f] & 0)
    {
      myno=0; // cout<<"f is "<<f<<" ------ "<<buff[f-1] <<endl;
      if(buff[f-1] & 1 ) {myno+=1;}// cout<<"1 ";}  else cout<<"0 ";
      if(buff[f-1] & 2 ) {myno+=2;}// cout<<"1 ";}  else cout<<"0 ";
      if(buff[f-1] & 4 ) { myno+=4;}// cout<<"1 ";}  else cout<<"0 ";
      if(buff[f-1] & 8 )  {myno+=8; }//cout<<"1 ";}  else cout<<"0 ";
      if(buff[f-1] & 16 ) {myno+=16;}// cout<<"1 ";} else cout<<"0 ";
      if(buff[f-1] & 32 ) { myno+=32;}// cout<<"1 ";} else cout<<"0 ";
      if(buff[f-1] & 64 ) {myno+=64;}// cout<<"1 ";} else cout<<"0 ";
      if(buff[f-1] & 128 ) {myno+= 128;}// cout<<"1 ";} else cout<<"0 "<<endl;
      if(myno > 255)if(command.debugmode == 1){std::cout<<"greater "<<myno<<std::endl;}
      //  if(myno <50) hist->Fill(myno);  if(myno >=50 && myno<240) histII->Fill(myno);  if(myno >=240) histIII->Fill(myno);
      if((f-1)%3==0) { /*hist->Fill(myno);*/ total->AddBinContent(n, myno);
      /* TwoD->SetBinContent(n,myno);*/ n++;} 
      if(f%3 ==0) {/* histIII->Fill(myno)*/ ; total->AddBinContent(nIII, myno);
      /*TwoDIII->SetBinContent(nIII,myno);*/ nIII++;}  
      if((f-2)%3==0){/* histII->Fill(myno);*/ total->AddBinContent(nII, myno);
      /* TwoDII->SetBinContent(nII,myno);*/ nII++;}
      if(myno>max) max=myno; 
    }
  
  if(toDraw){c1->cd(1);total->DrawCopy("colz");} //Draw original image
  
  for(Int_t i=0; i<bufferSize/3; i++)
    {
      Int_t x = i%w;
      Int_t y = (int) (i/w);
      ImageHistogram->SetBinContent(ImageHistogram->GetBin(x, y) , total->GetBinContent(i)); 
    }
  
  if(toDraw){c1->cd(2); ImageHistogram->DrawCopy("colz");}
  
  
  for(Int_t i=0; i<int(bufferSize/300); i++)
    ImgHistRebinned->SetBinContent(i,0);
  
  for(Int_t i=0; i<int (bufferSize/3); i++)
    {
      ImageHistogram->GetBinXYZ(i, xbin, ybin, zbin);
      bin= ImgHistRebinned->GetBin((int)(xbin/10), (int)(ybin/10));
      ImgHistRebinned->AddBinContent(bin, ImageHistogram->GetBinContent(i));   //filling ImgHistRebinned
    }
  
  if(toDraw)  {    c1->cd(3); ImgHistRebinned->DrawCopy("colz");    }
  
  if(oneOrtwo == '1') 
    {
      radius1=60; radius2=40; deltaRad=150; thresholdII=0.8;
      overallcourse=(ImgHistRebinned->Integral())/10000000;
      overallfine=   (ImageHistogram->Integral())/10000000; 
      writefile<<overallcourse<<" "<<overallfine<<" "; //--->write 
    }
  if(oneOrtwo == "2") 
    {
      if(camStr=="14") { radius1=300; radius2=280;}
      if(camStr=="15") { radius1=320; radius2=300;}
      if(camStr=="0") { radius1=270; radius2=250;} 
      if(camStr=="2") { radius1=270; radius2=250;} 
      if(camStr=="9") { radius1=300; radius2=280;} 
      if(camStr=="8") { radius1=280; radius2=260;} 
      if(camStr=="3") { radius1=320; radius2=300;} 
      if(camStr=="4") { radius1=260; radius2=240;} 
      if(camStr=="5") { radius1=320; radius2=300;} 
      if(camStr=="6") { radius1=350; radius2=320;} 
      if(camStr=="7") { radius1=300; radius2=280;} 
      thresholdII=0.5;
      overallcourse=(ImgHistRebinned->Integral())/1000000000;
      overallfine=   (ImageHistogram->Integral())/1000000000; 
      writefile<<overallcourse<<" "<<overallfine<<" "; 
    }
  std::cout<<"Info from the filename: RICH"<<oneOrtwo<<"  camera "<<camStr<<" radii to use in analysis "<<radius1<<"   "<<radius2<<std::endl;
  
  //Spot1   
  if(ImgHistRebinned->GetMaximum()>0)
    {
      point1= ImgHistRebinned->GetMaximumBin();
      //Spot 1 from ImgHistRebinned
      Double_t content= ImgHistRebinned->GetBinContent(point1);
      ImgHistRebinned->GetMaximumBin(xbin, ybin, zbin);
      yvalue=ImgHistRebinned->GetYaxis()->GetBinCenter(ybin);
      xvalue=ImgHistRebinned->GetXaxis()->GetBinCenter(xbin);       // F.Soomro 6 April 09
      //  if(oneOrtwo=='1' && yvalue>0 && camStr=="4")         
      //    xvalue=ImgHistRebinned->GetXaxis()->GetBinCenter(xbin)-(content-39000)/325;// ------> for R1 box 4
      //    else xvalue=ImgHistRebinned->GetXaxis()->GetBinCenter(xbin);
      // if(oneOrtwo=='1' && camStr =='4'){                                       // commented out: F.Soomro 6 April 09  
      //      if(command.debugmode == 1){
      //        std::cout<<"minus "<<content<<" "<<(content-39000)/325<<std::endl;
      //        std::cout<<"For R1 cam4: content is "<<content<<"  so subtracting  "<<(content-39000)/325<<std::endl;
      //      }
      //    }
      //   std::cout<<"For Spot 1 from the rebinned histogram --- max bin no: "<<point1<<", content of that bin: "<<
      // ImgHistRebinned->GetMaximum()<<",  x,y centers of the bin: "<<xvalue<<" "<<yvalue<<std::endl;
      
      for(Int_t i=0; i< w/10; i++)
	for(Int_t j=0; j<h/10; j++)
	  {
	    Double_t xx=ImgHistRebinned->GetXaxis()->GetBinCenter(i);
	    Double_t yy=ImgHistRebinned->GetYaxis()->GetBinCenter(j);      
	    Double_t cc=ImgHistRebinned->GetBinContent(i,j);
	    if( sqrt(pow(xvalue-xx,2)+pow(yvalue-yy,2))< radius1 )  
	      {
		xbuff+=xx*cc;
		ybuff+=yy*cc;
		cbuff+=cc;
		BuffHist->SetBinContent(i, j, cc);
		ImgHistRebinned->SetBinContent(i,j,0);
	      }
	  }
  
      xbuff=xbuff/cbuff;
      ybuff=ybuff/cbuff;
      std::cout<<"Preliminary calc of Spot 1 from the rebinned histogram --- ("<<xbuff<<","<<ybuff<<")"<<std::endl; 
      
      xvalue=xbuff;
      yvalue=ybuff;
      xbuff=0;
      ybuff=0;
      cbuff=0;
      buffermax=BuffHist->GetMaximum()*threshold;
      
      for(Int_t i=0; i<w/10; i++) {
	for(Int_t j=0; j<h/10; j++){
	  Double_t cc=BuffHist->GetBinContent(i,j); 
	  if(cc < buffermax) { BuffHist->SetBinContent(i,j,0);}
	  else
	    {
	      Double_t xx=ImgHistRebinned->GetXaxis()->GetBinCenter(i);
	      Double_t yy=ImgHistRebinned->GetYaxis()->GetBinCenter(j);
	      if( sqrt(pow(xvalue-xx,2)+pow(yvalue-yy,2))<radius2)
		{
		  xbuff+=xx*cc;
		  ybuff+=yy*cc;
		  cbuff+=cc;
		  a++;
		}
	    }//else
	}//int j
      } //int i
      
      xbuff=xbuff/cbuff;
      ybuff=ybuff/cbuff;
      std::cout<<"Content of the max bin of BuffHist "<<BuffHist->GetMaximum();
      std::cout<<" (should be equal to the content printed two lines before)!!"<<std::endl;
      std::cout<<"Final calc of Spot 1: from the rebinned histogram --- ("<<xbuff<<","<<ybuff<<")";
      std::cout<<",  number of bins used for it: "<<a<<std::endl; 
      // Donot remove writefile<<BuffHist->GetMaximum()<<" ";
      Intensity[counter]=cbuff/a;
      Intensity[counter+1]=cbuff;
      xVal[counter]=xbuff;
      yVal[counter]=ybuff;
      counter++;
  
      
      if(toDraw) 
	{
	  c1->cd(4); ImgHistRebinned->DrawCopy("colz");
	  c1->cd(5); BuffHist->DrawCopy("colz");
	}
      
      
      //Spot 1 from ImageHistogram 
      if(oneOrtwo=='1' && camStr=="3" && yvalue<0) threshold=0.8;        //----------------->box 4       
      xvalue=xbuff;
      yvalue=ybuff;
      xbuff=0;
      ybuff=0;
      cbuff=0;
      a = 0;
      
      for(Int_t i=0; i<w; i++){
	for(Int_t j=0; j<h; j++)
	  {
	    Double_t xx=ImageHistogram->GetXaxis()->GetBinCenter(i);
	    Double_t yy=ImageHistogram->GetYaxis()->GetBinCenter(j);    
	    Double_t cc=ImageHistogram->GetBinContent(i,j);
	    if(sqrt(pow(xvalue-xx,2)+pow(yvalue-yy,2))<radius1)
	      {
		BufferHistogram->SetBinContent(i,j, cc);
		ImageHistogram->SetBinContent(i,j,0);
	      }
	  }//int j
      }//int i
      buffermax=BufferHistogram->GetMaximum()*threshold;
      
      for(Int_t i=0; i<w; i++) {
	for(Int_t j=0; j<h; j++){
	  Double_t xx=ImageHistogram->GetXaxis()->GetBinCenter(i);
	  Double_t yy=ImageHistogram->GetYaxis()->GetBinCenter(j);
	  Double_t cc=BufferHistogram->GetBinContent(i,j);
	  if( cc < buffermax) { BufferHistogram->SetBinContent(i,j,0); }
	  else
	    if( sqrt(pow(xvalue-xx,2)+pow(yvalue-yy,2))<radius2) 
	      {
		xbuff+=xx*cc;
		ybuff+=yy*cc;               
		cbuff+=cc;
		a++;
	      }
	}
      }
      
      xbuff=xbuff/cbuff;
      ybuff=ybuff/cbuff;
      std::cout<<"Final calc of Spot 1: from the original histogram --- ("<<xbuff<<","<<ybuff<<")";
      std::cout<<",  number of bins used for it: "<<a<<std::endl;
      
      // Donot remove  writefile<<BufferHistogram->GetMaximum()<<" ";
      // Donot remove Intensity[counter]=cbuff/a;// BufferHistogram->GetMaximum(); 
      xVal[counter]=xbuff;
      yVal[counter]=ybuff;
      counter++;
      
      if(toDraw)
	{
	  c1->cd(6); ImageHistogram->DrawCopy("colz");
	  c1->cd(7); BufferHistogram->DrawCopy("colz"); 
	}
      
    }          // matches if max>0
  if(oneOrtwo=='1')  point1=int(ybuff);   
  if(oneOrtwo=='2')  point1=int(xbuff);         //-------------------------->addition
  
  //Spot2
  BuffHist->Reset(); 
  if(ImgHistRebinned->GetMaximum()>0)
    {
      if(oneOrtwo=='1' )
	{
	  BuffHist->Reset();
	  Int_t xbinII=0; Int_t ybinII=0;
	  point2= ImgHistRebinned->GetMaximumBin(xbinII, ybinII, zbin);
	  Double_t deltaX=ImgHistRebinned->GetXaxis()->GetBinCenter(xbin)- ImgHistRebinned->GetXaxis()->GetBinCenter(xbinII);
	  Double_t deltaY=ImgHistRebinned->GetYaxis()->GetBinCenter(ybin)- ImgHistRebinned->GetYaxis()->GetBinCenter(ybinII);
	  Double_t deltaR=sqrt(pow(deltaX,2)+pow(deltaY,2));
	  while (deltaR<=deltaRad)
	    {
	      ImgHistRebinned->SetBinContent(point2,0);
	      point2= ImgHistRebinned->GetMaximumBin(xbinII, ybinII, zbin); 
	      deltaX=ImgHistRebinned->GetXaxis()->GetBinCenter(xbin)-ImgHistRebinned->GetXaxis()->GetBinCenter(xbinII);
	      deltaY=ImgHistRebinned->GetYaxis()->GetBinCenter(ybin)-ImgHistRebinned->GetYaxis()->GetBinCenter(ybinII);
	      deltaR=sqrt(pow(deltaX,2)+pow(deltaY,2));
	    }
	}
      
      //Spot 2 from ImgHistRebinned
      xbuff=0;
      ybuff=0;
      cbuff=0;
      ImgHistRebinned->GetMaximumBin(xbin, ybin, zbin);
      xvalue=ImgHistRebinned->GetXaxis()->GetBinCenter(xbin);
      yvalue=ImgHistRebinned->GetYaxis()->GetBinCenter(ybin);
      std::cout<<"For Spot 2 from the rebinned histogram --- max bin no: "<<point2<<", content of that bin: ";
      std::cout<<ImgHistRebinned->GetMaximum()<<",  x,y centers of the bin: "<<xvalue<<" "<<yvalue<<std::endl;
      
      for(Int_t i=0; i< w/10; i++) {
	for(Int_t j=0; j<h/10; j++)
	  {    
	    Double_t xx=ImgHistRebinned->GetXaxis()->GetBinCenter(i);    
	    Double_t yy=ImgHistRebinned->GetYaxis()->GetBinCenter(j);
	    Double_t cc=ImgHistRebinned->GetBinContent(i,j);
	    //  if(cc>=ImgHistRebinned->GetMaximum()*threshold) 
	    if( sqrt(pow(xvalue-xx,2)+pow(yvalue-yy,2))<radius1)  
	      {
		xbuff+=xx*cc;
		ybuff+=yy*cc;
		cbuff+=cc;
		BuffHist->SetBinContent(i,j,cc);
		ImgHistRebinned->SetBinContent(i,j,0);
	      }
	  }//int j
      }//int i
      
      xbuff=xbuff/cbuff;
      ybuff=ybuff/cbuff;
      
      std::cout<<"Preliminary calc of Spot 2 from the rebinned histogram --- ("<<xbuff<<","<<ybuff<<")"<<std::endl; 
      
      if(oneOrtwo=='1' && camStr=="3"&& yvalue>0) thresholdII=0.5;                   // ----------> R1cam03 - previosuly box 4
      xvalue=xbuff;
      yvalue=ybuff;
      xbuff=0;
      ybuff=0;
      cbuff=0;
      a = 0;
      buffermax=BuffHist->GetMaximum()*threshold; 
      
      for(Int_t i=0; i< w/10; i++){
	for(Int_t j=0; j<h/10; j++)
	  {
	    Double_t xx=ImgHistRebinned->GetXaxis()->GetBinCenter(i);    
	    Double_t yy=ImgHistRebinned->GetYaxis()->GetBinCenter(j);
	    Double_t cc=BuffHist->GetBinContent(i,j); 
	    if(cc < buffermax) { BuffHist->SetBinContent(i,j,0);}
	    else if( sqrt(pow(xvalue-xx,2)+pow(yvalue-yy,2))< radius2) 
	      {
		xbuff+=xx*cc;
		ybuff+=yy*cc;
		cbuff+=cc;
		a++;
	      }
	  }//int j
      } //int i
      
      xbuff=xbuff/cbuff;
      ybuff=ybuff/cbuff;
      std::cout<<"Content of the max bin of BuffHist "<<BuffHist->GetMaximum();
      std::cout<<"  (should be equal to the content printed two lines before)!!"<<std::endl;
      std::cout<<"Final calc of Spot 2: from the rebinned histogram --- ("<<xbuff<<","<<ybuff<<")";
      std::cout<<",  number of bins used for it: "<<a<<std::endl; 

      Intensity[counter]=cbuff/a;//BuffHist->GetMaximum();
      Intensity[counter+1]=cbuff;
      xVal[counter]=xbuff;
      yVal[counter]=ybuff;
      counter++;
      if(oneOrtwo=='1')   point2 = int(ybuff); 
      if(oneOrtwo=='2')   point2 = int(xbuff);            //-----------------.addition

      if(toDraw)
	{
	  c1->cd(8); ImgHistRebinned->DrawCopy("colz");
	  c1->cd(9); BuffHist->DrawCopy("colz");
	}

      //Spot 2 from ImageHistogram
      
      xvalue=xbuff;
      yvalue=ybuff;
      xbuff=0;
      ybuff=0;
      cbuff=0;a= 0;
      BufferHistogram->Reset();
      
      for(Int_t i=0; i<w; i++){
	for(Int_t j=0; j<h; j++)
	  {
	    Double_t xx=ImageHistogram->GetXaxis()->GetBinCenter(i);
	    Double_t yy=ImageHistogram->GetYaxis()->GetBinCenter(j);    
	    Double_t cc=ImageHistogram->GetBinContent(i,j); 
	    if( sqrt(pow(xvalue-xx,2)+pow(yvalue-yy,2))< radius1)
	      {
		BufferHistogram->SetBinContent(i,j,cc);
		ImageHistogram->SetBinContent(i,j,0);
	      }
	  }//int j
      }//int i
      
      buffermax=BufferHistogram->GetMaximum()*thresholdII;

      for(Int_t i=0; i<w; i++){
	for(Int_t j=0; j<h; j++)
	  {
	    Double_t xx=ImageHistogram->GetXaxis()->GetBinCenter(i);
	    Double_t yy=ImageHistogram->GetYaxis()->GetBinCenter(j);    
	    Double_t cc=BufferHistogram->GetBinContent(i,j);
	    if( cc< buffermax)  {  BufferHistogram->SetBinContent(i,j,0);}
	    else
	      if( sqrt(pow(xvalue-xx,2)+pow(yvalue-yy,2))< radius2) 
		{
		  xbuff+=xx*cc;
		  ybuff+=yy*cc;           
		  cbuff+=cc;
		  a++;
		}
	  }//int j
      }//int i
      
      xbuff=xbuff/cbuff;
      ybuff=ybuff/cbuff;
      std::cout<<"Final calc of Spot 2: from the original histogram --- ("<<xbuff<<","<<ybuff<<")"<<",  number of bins used for it: "<<a<<std::endl;
      
      //  writefile<<BufferHistogram->GetMaximum()<<" ";
      //  Intensity[counter]=BufferHistogram->GetMaximum();
      xVal[counter]=xbuff;
      yVal[counter]=ybuff;
      counter++;
      
      if(toDraw)
	{
	  c1->cd(10); ImageHistogram->DrawCopy("colz");
	  c1->cd(11); BufferHistogram->DrawCopy("colz");
	}
      std::cout<<"Thresholds used on spot1 and spot2: "<<threshold<<" "<<thresholdII<<std::endl; 
    }// matches if max>0
  
  //**************ordering the values*********************//
  if(counter>4) std::cout<<"too many spots"<<std::endl;
  if(counter<4) std::cout<<"less than two spots"<<std::endl;
  if(counter==4)
    {
      std::cout<<"exactly two spots found"<<std::endl;
      //---------------------RICH 1--------------------//
      if(oneOrtwo=='1')
	{
	  if( (point1>point2 && camStr!='3')    || (point2>point1 && camStr=='3') )
	    {
	      std::cout<<"Rich"<<oneOrtwo<<", camera: "<<camStr<<" point1 is "<<point1<<" and point2 is "<<point2<<". Because ";
	      std::cout<<(point1>point2 ? point1 : point2) <<" is greater, changing the order now"<<std::endl;
	      x1coarse=xVal[0];
	      x1fine=xVal[1];
	      xVal[0]=xVal[2];
	      xVal[1]=xVal[3];
	      xVal[2]=x1coarse;
	      xVal[3]=x1fine;
	      y1coarse=yVal[0];
	      y1fine=yVal[1];
	      yVal[0]=yVal[2];
	      yVal[1]=yVal[3];
	      yVal[2]=y1coarse;
	      yVal[3]=y1fine;
	      intensityC=Intensity[0];
	      intensityF=Intensity[1];
	      Intensity[0]=Intensity[2];         // reference goes first,
	      Intensity[1]=Intensity[3]; 
	      Intensity[2]=intensityC;
	      Intensity[3]=intensityF;
	    } 
	  
	  //at this point, we are ready with our complete data point
	  int checked[4];
	  for(int i=0; i<4; i++) checked[i]=0;
	  if(camStr=='0') goto R1Cam1;
	  if(camStr=='1') goto R1Cam2;
	  if(camStr=='2') goto R1Cam3;
	  if(camStr=='3') goto R1Cam4;
	  
	checking:  //  ------------------ check if cam ID is correct --------------------//  
	  {
	  R1Cam1:
	    if(checked[0]!=1)
	      {
		TArrayD cam1=TArrayD(10);
		cam1[0]=131;
		cam1[1]=30;
		cam1[2]=5;
		cam1[3]=30;
		cam1[4]=-122;
		cam1[5]=30;
		cam1[6]=81;
		cam1[7]=20;
		cam1[8]=0.77;
		cam1[9]=0.3;
		checked[0]=1;
		std::cout<<"Checking for it to match R1cam00 (prev:cam1) "<<std::endl;
		std::cout<<xVal[1]-cam1[0]<<" " << xVal[3]-cam1[2] <<" " << yVal[1]-cam1[4]<<" " << yVal[3]-cam1[6]<<std::endl;
		if( TMath::Abs(xVal[1]-cam1[0])<cam1[1] && TMath::Abs(xVal[3]-cam1[2])<cam1[3] && 
		    TMath::Abs(yVal[1]-cam1[4])<cam1[5] && TMath::Abs(yVal[3]-cam1[6])<cam1[7] )
		  {
		    std::cout<<"Identified as R1cam00 (prev:cam1)"<<std::endl;
		    quality=0; camID='0'; goto PrintOutValues;    
		  }
		else
		  {
		    if( ( (TMath::Abs(xVal[1]-cam1[2])<cam1[3] && TMath::Abs(yVal[1]-cam1[6])<cam1[7]) ||  
			  (TMath::Abs(xVal[3]-cam1[0])<cam1[1] && TMath::Abs(yVal[3]-cam1[4])<cam1[5]) )
			|| ( (TMath::Abs(xVal[1]-cam1[0])<cam1[1] && TMath::Abs(yVal[1]-cam1[4])<cam1[5]) ||  
			     (TMath::Abs(xVal[3]-cam1[2])<cam1[3] && TMath::Abs(yVal[3]-cam1[6])<cam1[7]) ) )
		      {
			std::cout<<"Identified as R1cam00 (prev:cam1) but one spot only"<<std::endl;
			quality=2; camID='0';
			if(camStr=='2')  {  if (TMath::Abs(overallfine-cam1[8])<cam1[9]) goto PrintOutValues; }
			else { goto checking;}
		      }
		    else { goto checking; }
		  }
	      } // if checked[0]!=1

	  R1Cam2:
	    if(checked[1]!=1)
	      {
		TArrayD cam2=TArrayD(10);
		cam2[0]=-224;
		cam2[1]=30;
		cam2[2]=69;
		cam2[3]=30;
		cam2[4]=-162;
		cam2[5]=30;
		cam2[6]=93;
		cam2[7]=25;
		cam2[8]=1.11;
		cam2[9]=0.35;
		checked[1]=1;
		std::cout<<"Checking for it to match R1cam01 (prev:cam2) "<<std::endl;
		std::cout<<xVal[1]-cam2[0]<<" " << xVal[3]-cam2[2] <<" " << yVal[1]-cam2[4]<<" " << yVal[3]-cam2[6]<<std::endl;
		if( TMath::Abs(xVal[1]-cam2[0])<cam2[1] && TMath::Abs(xVal[3]-cam2[2])<cam2[3] && 
		    TMath::Abs(yVal[1]-cam2[4])<cam2[5] && TMath::Abs(yVal[3]-cam2[6])<cam2[7] )
		  {
		    std::cout<<"Identified as R1cam01 (prev:cam2)"<<std::endl;
		    quality=0; camID='1';   goto PrintOutValues;  
		  }
		else
		  {
		    if( ( (TMath::Abs(xVal[1]-cam2[2])<cam2[3] && TMath::Abs(yVal[1]-cam2[6])<cam2[7]) || 
			  (TMath::Abs(xVal[3]-cam2[0])<cam2[1] && TMath::Abs(yVal[3]-cam2[4])<cam2[5]) )
			||( (TMath::Abs(xVal[1]-cam2[0])<cam2[1] && TMath::Abs(yVal[1]-cam2[4])<cam2[5]) || 
			    (TMath::Abs(xVal[3]-cam2[2])<cam2[3] && TMath::Abs(yVal[3]-cam2[6])<cam2[7]) ) )
		      {
			std::cout<<"Identified as R1cam01 (prev:cam2) but one spot only"<<std::endl;
			quality=2; camID='1';  goto PrintOutValues;  
		      }
		    else { goto checking;}
		  }
	      }//if checked[1]!=1

	  R1Cam3:
	    if(checked[2]!=1)
	      {
		TArrayD cam3=TArrayD(10);
		cam3[0]=88;
		cam3[1]=25;
		cam3[2]=-32;
		cam3[3]=30;
		cam3[4]=-140;
		cam3[5]=25;
		cam3[6]=73;
		cam3[7]=25;
		cam3[8]=1.07;
		cam3[9]=0.22;
		checked[2]=1;
		std::cout<<"Checking for it to match R1cam02 (prev cam3) "<<std::endl;
		std::cout<<xVal[1]-cam3[0]<<" " << xVal[3]-cam3[2] <<" " << yVal[1]-cam3[4]<<" " << yVal[3]-cam3[6]<<std::endl;
		if( TMath::Abs(xVal[1]-cam3[0])<cam3[1] && TMath::Abs(xVal[3]-cam3[2])<cam3[3] && 
		    TMath::Abs(yVal[1]-cam3[4])<cam3[5] && TMath::Abs(yVal[3]-cam3[6])<cam3[7] )
		  {
		    std::cout<<"Identified as R1cam02 (prev cam3)"<<std::endl;
		    quality=0; camID='2';  goto PrintOutValues;  
		  } 
		else 
		  {
		    if( ( (TMath::Abs(xVal[1]-cam3[2])<cam3[3] && TMath::Abs(yVal[1]-cam3[6])<cam3[7]) || 
			  (TMath::Abs(xVal[3]-cam3[0])<cam3[1] && TMath::Abs(yVal[3]-cam3[4])<cam3[5]) )
			|| ( (TMath::Abs(xVal[1]-cam3[0])<cam3[1] && TMath::Abs(yVal[1]-cam3[4])<cam3[5]) || 
			     (TMath::Abs(xVal[3]-cam3[2])<cam3[3] && TMath::Abs(yVal[3]-cam3[6])<cam3[7]) ) )
		      {
			std::cout<<"Identified as R1cam02 (prev cam3) but one spot only"<<std::endl;
			quality=2; camID='2'; 
			if (camStr=='1') { if (TMath::Abs(overallfine-cam3[8]<cam3[9]))  goto PrintOutValues; }
			else { goto checking; }
		      }
		    else { goto checking; }
		  }
	      } //if checked[2]!=1
	    
	  R1Cam4:
	    if(checked[3]!=1)
	      {
		TArrayD cam4=TArrayD(10);
		cam4[0]=-165; 
		cam4[1]=30;
		cam4[2]=136;     
		cam4[3]=30;
		cam4[4]=103; 
		cam4[5]=20;
		cam4[6]=-116;   
		cam4[7]=30;
		cam4[8]=3.5;
		cam4[9]=0.33;
		checked[3]=1;
		std::cout<<"Checking for it to be R1cam03 (prev cam4) "<<std::endl;
		std::cout<<xVal[3]-cam4[0]<<" " << xVal[1]-cam4[2] <<" " << yVal[3]-cam4[4]<<" " << yVal[1]-cam4[6]<<std::endl;
		if( TMath::Abs(xVal[1]-cam4[0])<cam4[1] && TMath::Abs(xVal[3]-cam4[2])<cam4[3] && 
		    TMath::Abs(yVal[1]-cam4[4])<cam4[5] && TMath::Abs(yVal[3]-cam4[6])<cam4[7] )
		  {
		    std::cout<<"Identified as R1cam03 (prev cam4)"<<std::endl;
		    quality=0; camID='3';  goto PrintOutValues;  
		  } 
		else 
		  {
		    if( ( (TMath::Abs(xVal[1]-cam4[2])<cam4[3] && TMath::Abs(yVal[1]-cam4[6])<cam4[7]) || 
			  (TMath::Abs(xVal[3]-cam4[0])<cam4[1] && TMath::Abs(yVal[3]-cam4[4])<cam4[5]) )
			||( (TMath::Abs(xVal[1]-cam4[0])<cam4[1] && TMath::Abs(yVal[1]-cam4[4])<cam4[5]) || 
			    (TMath::Abs(xVal[3]-cam4[2])<cam4[3] && TMath::Abs(yVal[3]-cam4[6])<cam4[7]) ) )
		      {
			std::cout<<"Identified as R1cam03 (prev cam4) but one spot only"<<std::endl;
			quality=2; camID='3';  goto PrintOutValues;  
		      }
		    else  { goto checking; }
		  }
	      }//if checked[3]!=1
	    
	  } //matches checking  //------------------ check if cam ID is correct ends here--------------------// 
     } //matches if oneOrTwo==1
      
      
      //---------------------RICH 2--------------------//  
      if(oneOrtwo=='2')
	{
	  if(camStr=='5' || camStr=='3' || camStr=='4'|| camStr=='2' || 
	     camStr=="14"  || camStr=="15" || camStr=='6'  || camStr=='8' || camStr=='7')    
	    if(point1<point2)
              {
                std::cout<<"camera is "<<oneOrtwo<<", point1 is "<<point1<<" and point2 is "<<point2<<". Because ";
		std::cout<<(point1>point2 ? point1 : point2) <<"is greater, and camStr is "<<camStr<<", changing the order now"<<std::endl;
		x1coarse=xVal[0];
                x1fine=xVal[1];
                xVal[0]=xVal[2];
                xVal[1]=xVal[3];
                xVal[2]=x1coarse;
                xVal[3]=x1fine;
                y1coarse=yVal[0];
                y1fine=yVal[1];
                yVal[0]=yVal[2];
                yVal[1]=yVal[3];
                yVal[2]=y1coarse;
                yVal[3]=y1fine;
                intensityC=Intensity[0];
                intensityF=Intensity[1];
                Intensity[0]=Intensity[2];         
                Intensity[1]=Intensity[3]; 
                Intensity[2]=intensityC;
                Intensity[3]=intensityF;
              } 
 
	  if(camStr=='0'|| camStr=='9')    
	    if(point1>point2)
              {
                std::cout<<"camera is "<<oneOrtwo<<", point1 is "<<point1<<" and point2 is "<<point2<<". Because ";
		std::cout<<(point1>point2 ? point1 : point2) <<"is greater, and camStr is "<<camStr<<", changing the order now"<<std::endl;
		x1coarse=xVal[0];
                x1fine=xVal[1];
                xVal[0]=xVal[2];
                xVal[1]=xVal[3];
                xVal[2]=x1coarse;
                xVal[3]=x1fine;
                y1coarse=yVal[0];
                y1fine=yVal[1];
                yVal[0]=yVal[2];
                yVal[1]=yVal[3];
                yVal[2]=y1coarse;
                yVal[3]=y1fine;
                intensityC=Intensity[0];
                intensityF=Intensity[1];
                Intensity[0]=Intensity[2];         //in x and y, reference does first, but in Intensity, reflected goes first
                Intensity[1]=Intensity[3]; 
                Intensity[2]=intensityC;
                Intensity[3]=intensityF;
              } 
	  
	}// if oneOrTwo==2
      
    }  //counter ==4
  
  //**************finished ordering the values***************//
  
  //***********output values on screen and file*************//
 PrintOutValues:
  if(oneOrtwo=='1')
    {
      if(camID!=camStr)
	{
	  std::cout<<"camIDendtified as "<<camID<<"  but camStr was "<<camStr<<std::endl;
	  if(camID.Length()==0) {  quality=1; std::cout<<"so setting quality to 1"<<std::endl; }
	  else { quality=3; std::cout<<"camID.Length() was not 0, so setting quality to 3"<<"  camID was "<<camID<<" and camStr was "<<camStr<<std::endl;}
	}
      std::cout<<"camera# from file name:  "<<camStr<<", camera# identified from this code:  "<<camID<<
	" ,so setting the quality to "<<quality<<std::endl;
    } // matches oneOrtwo==1
  
  std::cout<<"The numbers as written into the data file:"<<std::endl;
  std::cout<<"Intensity Reference and Reflected"<<std::endl;
  
  for(Int_t i =0; i<counter; i++)
    {
      std::cout<<Intensity[i]<<" ";
      writefile<<Intensity[i]<<" ";
    }
  
  std::cout<<std::endl<<"The X and Y coordinates "<<std::endl;
  
  for(Int_t i =0; i<counter; i++)
    {
      std::cout<<xVal[i]<<" "<<yVal[i]<<std::endl;
      writefile<<xVal[i]<<" "<<yVal[i]<<" ";
    }
  
  if(oneOrtwo=='2')
    { 
      if(quality==-1) quality=0;
      if(overallfine<0.03){ quality = -1; std::cout<<"R2 message: the image was probably dark"<<std::endl;}
      if(overallfine>0.5) { quality = -1; std::cout<<"R2 message: the image was probably too bright"<<std::endl;}
    }
  //*****************************************************// 
  writefile<<quality<<std::endl;
  toDel.Append(Filename);
  std::cout<<"deleting the file  "<<Filename<<std::endl;
  gSystem->Exec(toDel);
  ImgAnalysisResults Results;
  Results.refXCoord = xVal[1];
  Results.refYCoord = yVal[1];
  Results.refLumi = Intensity[1];
  Results.reflXCoord = xVal[3];
  Results.reflYCoord = yVal[3];
  Results.reflLumi = Intensity[3];
  //   Results.radSep = sqrt(pow(xVal[3]-xVal[1],2)+pow(yVal[3]-yVal[1],2));      
  Results.imgLumi = overallfine;
  Results.quality = quality;
  std::cout<<"deleting the ifstream "<<std::endl;
  std::cout  << "******* The Output to PVSS ********* " << std::endl;
  std::cout << "Overall Image Lumi " << Results.imgLumi << std::endl;
  std::cout << "Reference Lumi " << Results.refLumi << std::endl;
  std::cout << "Reflected lumi " <<  Results.reflLumi << std::endl;
  std::cout << "Reference X " << Results.refXCoord << " and Y " <<  Results.refYCoord << std::endl;
  std::cout << "Reflected X " << Results.reflXCoord << " and Y " <<  Results.reflYCoord << std::endl;
  std::cout << "Image Quality " << Results.quality <<". (Remember that only 0 is good, the rest isn't )"<<std::endl;
  // std::cout << "Radial Separation " << Results.radSep <<std::endl;
  //  std::cout<<"************ End of log from this analysis ************ \n\n"<<std::endl;
  // return Results;
  goto theEnd;
  
 BadImage:
  {
    toDel.Append(Filename);
    std::cout<<"deleting the file  "<<Filename<<std::endl;
    gSystem->Exec(toDel);
    // delete c1;
    if(quality==4)
      {
	// std::cout<<"convert command not happy with the image\n\n"<<std::endl;
	// std::cout<<"the corrupt bmp file: "<<Filename<<std::endl;
	std::cout<<"Couldnt open the bmp file and deleting the ifstream now"<<std::endl;
      }
    if(quality != 4)  
      std::cout<<"convert command was happy but still a file opening error!!\n\n"<<std::endl;
    
    ImgAnalysisResults Results;
    Results.refXCoord = xVal[1];
    Results.refYCoord = yVal[1];
    Results.refLumi = Intensity[1];
    Results.reflXCoord = xVal[3];
    Results.reflYCoord = yVal[3];
    Results.reflLumi = Intensity[3];
    // Results.radSep = sqrt(pow(xVal[3]-xVal[1],2)+pow(yVal[3]-yVal[1],2));
    Results.imgLumi = overallfine;
    Results.quality = quality;
    std::cout  << "******* The Output to PVSS ********* " << std::endl;
    std::cout << "Overall Image Lumi " << Results.imgLumi << std::endl;
    std::cout << "Reference Lumi " << Results.refLumi << std::endl;
    std::cout << "Reflected lumi " <<  Results.reflLumi << std::endl;
    std::cout << "Reference X " << Results.refXCoord << " and Y " <<  Results.refYCoord << std::endl;
    std::cout << "Reflected X " << Results.reflXCoord << " and Y " <<  Results.reflYCoord << std::endl;
    std::cout << "Image Quality (0 is good, the rest isn't ) " << Results.quality <<std::endl;
    //std::cout << "Radial Separation " << Results.radSep <<std::endl;
    for(int i=0; i<14; i++)   writefile<<100.01<<"  ";
    writefile<<quality<<std::endl;
    // return Results;   
    goto theEnd;
  }//bad image
 theEnd:
  {
    std::cout<<"************ End of log from this analysis ************ \n\n"<<std::endl;
    gSystem->RedirectOutput(0);
    delete ImageHistogram;
    delete ImgHistRebinned;
    delete BufferHistogram;
    delete BuffHist;
    delete total;
    delete c1;
    delete myfile;
    writefile.close();
    LogFile.close();
    return Results;  
   }
  
  
}