// 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 <fstream>  
#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 compute(struct LAMSCommand command);

int main()
{ 
  //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,0};
  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);
      // update the service realted to this camera
      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);
    gain.Remove(0,30);
    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("output");
  gSystem->Exec(toConvert);
  gSystem->RedirectOutput(LogFileName); 
  ifstream test("output", std::ios::binary);
  test.seekg(0, std::ios::end);
  if (int(test.tellg()) !=0)   quality=4;
  gSystem->Exec("rm output");
  //compose .bmp Filename
  TString Filename = path;
  Filename.Append(copyfilename);
  Filename.Append(".bmp");
  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 "; 
  
  ifstream myfile(Filename, std::ios::binary);
  //  if(quality==4) goto BadImage;
  char buff[bufferSize];//1310720-54]; 
  if(!myfile) { std::cout<<"error opening your file"<<std::endl; goto BadImage;}
  if(myfile)
    {
      myfile.seekg(54, std::ios::beg);
      myfile.read(buff, sizeof(buff));
      myfile.close();
    }

   TCanvas* c1 = new TCanvas("image","image",1200, 700);
  if(toDraw) c1->Divide(4,3);
  
  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); 
  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");}
  

 
  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); 
  
  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++;
       }
      }
    }
     
  
  xbuff=xbuff/cbuff;
  ybuff=ybuff/cbuff;
  std::cout<<"Content of the max bin of BuffHist "<<BuffHist->GetMaximum()<<
    "  (should be equal to the content printed two lines before)!!"<<std::endl<<
    "Final calc of Spot 1: from the rebinned histogram --- ("<<xbuff<<","<<ybuff<<")"<<
    ",  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);
      }
    }

  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<<")"<<
    ",  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: "<<
     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);
      }
    }
   
  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++;
      }
    }
  
  xbuff=xbuff/cbuff;
  ybuff=ybuff/cbuff;
  std::cout<<"Content of the max bin of BuffHist "<<BuffHist->GetMaximum()<<
    "  (should be equal to the content printed two lines before)!!"<<std::endl<<
    "Final calc of Spot 2: from the rebinned histogram --- ("<<xbuff<<","<<ybuff<<")"<<
    ",  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);
      }
    }
   

 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++;
      }
    }
   
  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)
  {
    //---------------------RICH 1--------------------//
    if(oneOrtwo=='1')
      {
        if( (point1>point2 && camStr!='3')    || (point2>point1 && camStr=='3') )
           {
             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]=116;
           cam1[1]=20;
           cam1[2]=24;
           cam1[3]=25;
           cam1[4]=-125;
           cam1[5]=20;
           cam1[6]=78;
           cam1[7]=15;
           cam1[8]=2.67;
           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;
         }
       R1Cam2:
         if(checked[1]!=1)
         {
           TArrayD cam2=TArrayD(10);
           cam2[0]=-202;
           cam2[1]=20;
           cam2[2]=50;
           cam2[3]=35; //25; // F.Soomro 6 April 09
           cam2[4]=-161;
           cam2[5]=30;
           cam2[6]=90;
           cam2[7]=25;
           cam2[8]=2.7;
           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;
         }
       R1Cam3:
         if(checked[2]!=1)
         {
           TArrayD cam3=TArrayD(10);
           cam3[0]=74;
           cam3[1]=15;
           cam3[2]=-18;
           cam3[3]=50;//25; // F.Soomro 6 April 09
           cam3[4]=-140;
           cam3[5]=20;
           cam3[6]=64;
           cam3[7]=20;
           cam3[8]=1.04;
           cam3[9]=0.2;
           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;
         }
       R1Cam4:
         if(checked[3]!=1)
         {
           TArrayD cam4=TArrayD(10);
           cam4[0]=-148; //-124;   F.Soomro 6 apr09   
           cam4[1]=30;
           cam4[2]=115;//65;       F.Soomro 6 apr09   
           cam4[3]=20;//25;        F.Soomro 6 apr09   
           cam4[4]=103;//37;     F.Soomro 6 apr09   
           cam4[5]=10;
           cam4[6]=-112;//-159;    F.Soomro 6 apr09   
           cam4[7]=30;//15;        F.Soomro 6 apr09  
           cam4[8]=3.44;
           cam4[9]=0.3;
           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;
         }
       } //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)
              {
                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)
              {
                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)
       if(camID.Length()==0)  quality=1; 
       else quality=3;
     std::cout<<"camera# from file name:  "<<camStr<<", camera# identified from this code:  "<<camID<<
       " ,so setting the quality to "<<quality<<std::endl;
   }
   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;
   std::cout<<"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;}
     }
   //*****************************************************// 
   writefile<<quality<<std::endl;
   writefile.close();
   LogFile.close(); 
   delete ImageHistogram;
   delete ImgHistRebinned;
   delete BufferHistogram;
   delete BuffHist;
   delete total;
   delete c1;
   toDel.Append(Filename);
   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  << "******* 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<<"\n\n"<<std::endl;
   return Results;
   goto theEnd;
  
 BadImage:
   {
     // delete c1;
     toDel.Append(Filename);
     gSystem->Exec(toDel);
     if(quality==4) std::cout<<"convert command not happy with the image\n\n"<<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<<"\n\n"<<std::endl;
     for(int i=0; i<14; i++)
       writefile<<100.01<<"  ";
     writefile<<quality<<std::endl;
     writefile.close();
     LogFile.close();
     return Results;   
     goto theEnd;
   }
 theEnd:
   {
     gSystem->RedirectOutput(0);
   }


}