#include <iostream>
#include <fstream>
#include <sstream>
#include <string>
#include <iomanip>
#include <limits>
#include <math.h>
#include <TArc.h>
#include <TArrow.h>
#include <TBox.h>
#include <TBrowser.h>
#include <TCanvas.h>
#include <TCondition.h>
#include <TDatime.h>
#include <TDirectory.h>
#include <TEllipse.h>
#include <TGraphErrors.h>
#include <TGButton.h>
#include <TGButtonGroup.h>
#include <TGClient.h>
#include <TGComboBox.h>
#include <TGProgressBar.h>
#include <TGLabel.h>
#include <TGTextEntry.h>
#include <TGFrame.h>
#include <TGraph.h>
#include <TGraphErrors.h>
#include <TF1.h>
#include <TFile.h>
#include <TH1D.h>
#include <TH2D.h>
#include <TH3I.h>
#include <TLatex.h>
#include <TLeaf.h>
#include <TLine.h>
#include "TMath.h"
#include <TObjectTable.h>
#include <TPaveText.h>
#include <TPolyLine.h>
#include <TProfile.h>
#include <TRandom.h>
#include <TRootEmbeddedCanvas.h>
#include <TStyle.h>
#include "TSystem.h"
#include <TTree.h>
#include <TVectorD.h>
#include <TDatime.h>
#include <sys/types.h>
#include <time.h>
#include <TROOT.h>
#include <TError.h>
#include "TVirtualFFT.h"
#define SizeOfArray(a) (sizeof(a)/sizeof(a[0]))
/*
#include <classes/miniBurstHeader.C>
#include <classes/miniEventHeader.C>
#include <classes/eventEssentials.C>
#include <classes/cherenkovEssentials.C>
*/
const double lapTime=165.*30.*1.0e-6; // accelerator lap in msec
int iCol=1;

struct miniEventHeader {
  Int_t burst;          // burst number in run (some may be missing due to 
                        // missing or rejected GADC data) 
  Int_t event;          // event# in run
  Int_t burstEvent;     // event# in burst
  Int_t scaler10mhz;    // 10MHz counter
  Int_t scaler6mhz;     // 6MHz (165ns, accelerator RF) counter
  Int_t firstGated;     // same as firstGADCsample above...???

  Bool_t   isEmpty;     // no activity in S123 and Cherenkov for low rate, 
                        // no activity in (some of) modules and Cherenkov for high rate
  Bool_t   notEmpty;    // some activity detected!
  Bool_t   S123inBeam;  // S123 are placed in the beam
  Bool_t   filledSegment;  // event is in 'active' part of GatedADC record

  Double_t RFperiod;    // RF A*sin(omega*(t+t0)) fit
  Double_t RFfreq;      // RF A*sin(omega*(t+t0)) fit
  Double_t RFamp;       // RF A*sin(omega*(t+t0)) fit
  Double_t TimeGlobal;  // RF A*sin(omega*(t+t0)) fit
  Double_t SignalFactor;// effective signal/noise ratio in GlobalPhase finder

  void Print() {
    printf("[miniEventHeader->Print] burst#%03d event#%05d(%03d) scaler10=%07d scaler6=%07d firstGADCsample=%07d\n",
            burst, event, burstEvent, scaler10mhz, scaler6mhz, firstGated);
    printf("isEmpty=%1d notEmpty=%1d S123beam=%1d filledSegmen=%1d\n ", isEmpty, notEmpty, S123inBeam, filledSegment);
    printf("RFperiod=%.1lfns RFfreq=%6.1e RFamp=%5.1e Tglobal=%5.1e SignalFactor=%5.1e\n", 
            (Double_t)RFperiod, (Double_t)RFfreq, (Double_t)RFamp, (Double_t)TimeGlobal, (Double_t)SignalFactor);
  }
///////////////////////////////////////////////////////////////////////////////////
};


struct eventEssentials {

int firstGADCsample; // dead! supposed to be same as miniEventHeader.firstGated!

float GADCevent[40];  // 40 Gated ADC samples (RF cycles), starting with firstGADCsample
float GADCevent1[40]; // 40 Gated ADC samples (RF cycles), starting with firstGADCsample, tails subtracted
float MPICevent[40]; // 40 MPI ADC samples for Cherenkov. NOT ALL FILLED! Integral after FFT
float tRFb[40];      // time fof RF bunches in the "MPI time frame"

int bunchAction[40]; // 0 - useless bunch 1 - signal 2 - pedestal
                     // ATTN~! for true electronc pedestals use events with isEmpty=true!

int mipCount[40];    // likelihood track counting; wokrs for Ntracks<100 actually <~80
int mipCount2[40];   // sqrt(Sx*Sy), where Sx and Sy are two smallest scintillator signals (in MIPs)
int mipCount3[40];   // (S1*S2*S3)^(1/3), Si in MIPs
//int mipCountX[40];   // hybrid of mipCount[] at low rate and  mipCount2[] at high rate

bool isIsolated[40]; // indicates isolated filled bunch (isCluster0=isCluster1=true);
bool isCluster0[40]; // first in cluster flag
bool isCluster1[40]; // last in cluster flag
bool inScope[40];

float scintFFT[120]; // scintillators in MPI.FFT units. id=3*iRF+iScint iRF RF cycle in event

bool  HGoverflow[4]; // module has HighGain overflow
bool  LGoverflow[4]; // module has LowGain overflow
int   RFclust0[5];   // RF# for clusters, as indicated by isCluster0 and isCluster1 
int   RFclust1[5];
float miniChannel[200]; //fC, 4mod*5chan*2gain*5clust

/*
void computeMIPcountX() {
  for (int i=0; i<40; i++) {
    if (mipCount[i]<40) {
      mipCountX[i]=mipCount[i];
    } else if (mipCount[i]<60) {
      mipCountX[i]=int(0.5+((60-i)*mipCount[i]+(i-40)*mipCount2[i])/20.);
    } else {
      mipCountX[i]=mipCount2[i];
    }
  }
}
*/
int mipCountX(int i) {
  if (i<0 || i>39) {
    printf("[eventEssentials->mipCountX] bad index#%d\n", i);
    exit(-1);
  }
  if (mipCount[i]<60) {
    return mipCount[i];
  } else if (mipCount[i]<80) {
    return int(0.5+((80-i)*mipCount[i]+(i-60)*mipCount2[i])/20.);
  }
  return mipCount2[i];
}

void CorrectGADCtail() {

  for (int i=0; i<40; i++) {
    GADCevent1[i]=GADCevent[i];
  } 
  for (int i=0; i<38; i++) { 
    GADCevent1[i+1]-=0.16*GADCevent[i];
    GADCevent1[i+2]-=0.04*GADCevent[i];
  }
}

void Clear () {
  firstGADCsample=-1000;
  for (int i=0; i<120; i++) {
    scintFFT[i]=0.;
  }
  for (int i=0; i<40; i++) {
    GADCevent[i]=0.;
    MPICevent[i]=0.;

    tRFb[i]=-10000.;

    mipCount[i]=0;
    mipCount2[i]=0;
    mipCount3[i]=0;

    isIsolated[i]=false;
    isCluster0[i]=false;
    isCluster1[i]=false;
    inScope[i]=false;
  }
  for (int i=0; i<4; i++) {
    HGoverflow[i]=false;
    LGoverflow[i]=false;
  }
  for (int i=0; i<5; i++) {
    RFclust0[i]=-1;
    RFclust1[i]=-1;
  }
  for (int i=0; i<200; i++) {
    miniChannel[i]=0.;
  }
}
};


//////////////////////////////////////////////////////////////////////////
struct miniBurstHeader {
  Int_t   burst;  // burst number, must be same as in Event Header. Some tricks used in code
  Float_t ChMPIF; // Cherenkov calibration, MPI FFT
  Float_t ChGADC; // Cherenkov calibration, GADC
};

///////////////////////////////////////////////////////////////////////////////////
struct cherenkovEssentials {
//   Int_t buffSize=1024*1024*16/30+1;
   Int_t     ns30;          // number of accelerator beam 'laps', 30 RFcycles/lap
   Float_t   sum30[559241]; // reserved buffer for prehistory, sum of filled bunches in each "lap"  
   Bool_t    lapBool[30];   // filled lap flag
   Double_t  lapFill[30];   // sum over all laps with same N%30
   Int_t     lapActn[30];   // what to do with that lap (see bunchAction[40] above)
   
   float lapMax;            // maximum lap sum, GADC 

   void ComputeLapMax() {
     for (int i=0; i<ns30; i++) {if (lapMax<sum30[i]) lapMax=sum30[i];}
   }

   void Compute() {
     ComputeLapMax();
   }

   void Clear() {
     ns30=0;
     lapMax=0.;
     for (int i=0; i<559241; i++) {
       sum30[i]=0.;
     }
     for (int i=0; i<30; i++) {
       lapBool[i]=false;
       lapFill[i]=0.0;
       lapActn[i]=0;
     }
   }
}; 
/////////////////////////////////////////////////////////////////////////////////
struct scintMIPs {
  double scintMIP[3]; //scintillator MIPs for the Run in MPI FFT units (integral)
};
////////////////////////////////////////////////////////////////////////////////
struct pedestalWindows {
  int nPedWindows;      // number of valid Pedestal windows in the event
                        // those fancy yellow bars on the event display
  int firstSample[25];  // first sample for pedestal window
  int lastSample[25];   // last  sample for pedestal window

  void Print() {
    printf("nPed=%d\n", nPedWindows);
    for (int i=0; i<nPedWindows; i++) {
      printf(   "#%d %3d...%3d\n", i, firstSample[i], lastSample[i]);
    }
  }

};

struct moduleSignalsF {
  float ampI[512]; // signal shape in module channel, raw-ped
  float ampF[512]; // FTT-reconstructed curremt, microAmps
  double X;        // X center of gravity,  defined by geometry
  double Y;        // Y center of gravity,  defined by geometry
  double eventSum;
  void Clear() {eventSum=0.;}
};

struct Module { 
  double X;        // X center of gravity,  entire event
  double Y;        // Y center of gravity,  entire event
  double Q;        // module Sum, entire event;
  double Xcl[5];   // X center of gravity, per cluster
  double Ycl[5];   // Y center of gravity, per cluster
  double Qcl[5];

  void Add(double A, double x, double y, int iCl) {
    X+=A*x;
    Y+=A*y;
    Q+=A;
    Xcl[iCl]+=A*x;
    Ycl[iCl]+=A*y;
    Qcl[iCl]+=A;    
  }

  void Compute() {
    if (Q>0.) {
      X/=Q;
      Y/=Q;
    } else {
      X=-1.e10;
      Y=-1.e10;
    }
    for (int i=0; i<5; i++) {
      if (Qcl[i]>0.) {
        Xcl[i]/=Qcl[i];
        Ycl[i]/=Qcl[i];
      } else {
        Xcl[i]=-1.e10;
        Ycl[i]=-1.e10;
      }      
    }
  }

  void Clear() {
    X=0.;
    Y=0.;
    Q=0;
    for (int i=0; i<5; i++) {Xcl[i]=0.; Ycl[i]=0.; Qcl[i]=0.;}
  }
};

// computes covariance, corrleation coefficient, slope and variance of all for 2 variables
struct myCovariance {
  double sumW[11]; //0...9 subsets, 10 - entire dataset
  double sumX[11];
  double sumY[11];
  double sumXX[11];
  double sumYY[11];
  double sumXY[11];
  int    indx;

void Clear() {
  for (int i=0; i<11; i++) {
    sumW[i]=0.;
    sumX[i]=0.;
    sumY[i]=0.;
    sumXX[i]=0.;
    sumYY[i]=0.;
    sumXY[i]=0.;
    indx=0;
  }
}

void Fill(double x, double y, double w) {
  sumW[10]+=w;
  sumX[10]+=x*w;
  sumY[10]+=y*w;
  sumXX[10]+=x*x*w;
  sumYY[10]+=y*y*w;
  sumXY[10]+=x*y*w;
  sumW[indx]+=w;
  sumX[indx]+=x*w;
  sumY[indx]+=y*w;
  sumXX[indx]+=x*x*w;
  sumYY[indx]+=y*y*w;
  sumXY[indx]+=x*y*w;
//  printf("[myCovariance->Fill] filling indx=%d w=%.1f x=%1.f y=%.1f\n", indx, w, x, y);
  indx++;
  if (indx==10) indx=0;
}

void Fill(double x, double y) {
  Fill(x, y, 1.0);
}

double GetWeight(int i) {
  if (i<0 || i>10) {
    printf("[myCovariance->GetWeight] wrong index %d\n", i);
    exit(-1);
  }
  return sumW[i];
}
double GetAverageX(int i) {
  if (i<0 || i>10 || sumW[i]==0.) {
    printf("[myCovariance->GetAverageX] wrong weight=%f or index %d\n", sumW[i], i);
    exit(-1);
  }
  return sumX[i]/sumW[i];
}
double GetAverageY(int i) {
  if (i<0 || i>10 || sumW[i]==0.) {
    printf("[myCovariance->GetAverageY] wrong weight=%f or index %d\n", sumW[i], i);
    exit(-1);
  }
  return sumY[i]/sumW[i];
}
double GetAverageX2(int i) {
  if (i<0 || i>10 || sumW[i]==0.) {
    printf("[myCovariance->GetAverageX2] wrong weight=%f or index %d\n", sumW[i], i);
    exit(-1);
  }
  return sumXX[i]/sumW[i];
}
double GetAverageY2(int i) {
  if (i<0 || i>10 || sumW[i]==0.) {
    printf("[myCovariance->GetAverageY2] wrong weight=%f or index %d\n", sumW[i], i);
    exit(-1);
  }
  return sumYY[i]/sumW[i];
}
double GetSigmaX(int i) {
  if (i<0 || i>10 || sumW[i]==0.) {
    printf("[myCovariance->GetSigmaX] wrong weight=%f or index %d\n", sumW[i], i);
    exit(-1);
  }
  double result=sumXX[i]/sumW[i]-GetAverageX(i)*GetAverageX(i);
  if (result<0.) result=0.0;
  return sqrt(result);
}
double GetSigmaX2(int i) {
  if (i<0 || i>10 || sumW[i]==0.) {
    printf("[myCovariance->GetSigmaX2] wrong weight=%f or index %d\n", sumW[i], i);
    exit(-1);
  }
  double result=sumXX[i]/sumW[i]-GetAverageX(i)*GetAverageX(i);
  if (result<0.) result=0.0;
  return result;
}
double GetSigmaY(int i) {
  if (i<0 || i>10 || sumW[i]==0.) {
    printf("[myCovariance->GetSigmaY] wrong weight=%f or index %d\n", sumW[i], i);
    exit(-1);
  }
  double result=sumYY[i]/sumW[i]-GetAverageY(i)*GetAverageY(i);
  if (result<0.) result=0.0;
  return sqrt(result);
}
double GetSigmaY2(int i) {
  if (i<0 || i>10 || sumW[i]==0.) {
    printf("[myCovariance->GetSigmaY] wrong weight=%f or index %d\n", sumW[i], i);
    exit(-1);
  }
  double result=sumYY[i]/sumW[i]-GetAverageY(i)*GetAverageY(i);
  if (result<0.) result=0.0;
  return result;
}
double GetCorrelation(int i) {
  if (i<0 || i>10 || sumW[i]==0.) {
    printf("[myCovariance->GetCorrelation] wrong weight=%f or index %d\n", sumW[i], i);
    exit(-1);
  }
//  double result=sumW[i]/(sumW[i]-1.)*(sumXY[i]/sumW[i]-GetAverageY(i)*GetAverageX(i))/GetSigmaX(i)/GetSigmaY(i);
  double result=(sumXY[i]/sumW[i]-GetAverageY(i)*GetAverageX(i))/GetSigmaX(i)/GetSigmaY(i);
  return result;
}  
double GetCorrelationVariance() {
  double c;
  double c0=0.;
  double c1=0.;
  double c2=0.;
  for (int i=0; i<10; i++) {
    c=GetCorrelation(i);
    c0+=sumW[i];
    c1+=c*sumW[i];
    c2+=c*c*sumW[i];
  }
  c1/=c0;
  c2/=c0;
  c2-=c1*c1;
  if (c2>0.) c2=sqrt(c2);
  else       c2=0.;
  return c2;
}

double GetSlope(int i) {
  if (i<0 || i>10 || sumW[i]==0.) {
    printf("[myCovariance->GetSlope] wrong weight=%f or index %d\n", sumW[i], i);
    exit(-1);
  }
  double x=GetSigmaX(i)*GetSigmaX(i)-GetSigmaY(i)*GetSigmaY(i);
  double y=2.0*(sumXY[i]/sumW[i]-GetAverageY(i)*GetAverageX(i));
  double alpha=0.5*atan2(y, x);
  double slope=tan(alpha);
  return slope;
}  
double GetAngle(int i) {
  if (i<0 || i>10 || sumW[i]==0.) {
    printf("[myCovariance->GetAngle] wrong weight=%f or index %d\n", sumW[i], i);
    exit(-1);
  }
  double x=GetSigmaX(i)*GetSigmaX(i)-GetSigmaY(i)*GetSigmaY(i);
  double y=2.0*(sumXY[i]/sumW[i]-GetAverageY(i)*GetAverageX(i));
  double alpha=0.5*atan2(y, x);
  return alpha;
}
  
double GetSlope0(int i) { // special case with line, going through (0,0)
  if (i<0 || i>10 || sumW[i]==0.) {
    printf("[myCovariance->GetSlope0] wrong weight=%f or index %d\n", sumW[i], i);
    exit(-1);
  }
  double x=sumXX[i]/sumW[i]-sumYY[i]/sumW[i];
  double y=2.0*sumXY[i]/sumW[i];
  double alpha=0.5*atan2(y, x);
  double slope=tan(alpha);
  return slope;
}  

double GetAngle0(int i) {// special case with line, going through (0,0)
  if (i<0 || i>10 || sumW[i]==0.) {
    printf("[myCovariance->GetAngle0] wrong weight=%f or index %d\n", sumW[i], i);
    exit(-1);
  }
  double x=GetSigmaX(i)*GetSigmaX(i)-GetSigmaY(i)*GetSigmaY(i);
  double y=2.0*sumXY[i]/sumW[i];
  double alpha=0.5*atan2(y, x);
  return alpha;
}  

double GetSlopeX(int i) { // classic linear fit y=b+a*(x-<x>);
  if (i<0 || i>10 || sumW[i]==0.) {
    printf("[myCovariance->GetSlopeX] wrong weight=%f or index %d\n", sumW[i], i);
    exit(-1);
  }
  double slope=(sumXY[i]/sumW[i]-GetAverageY(i)*GetAverageX(i))/GetSigmaX2(i);
  return slope;
}
double GetSlopeXVariance() {
  double c;
  double c0=0.;
  double c1=0.;
  double c2=0.;
  for (int i=0; i<10; i++) {
    c=GetSlopeX(i);
    c0+=sumW[i];
    c1+=c*sumW[i];
    c2+=c*c*sumW[i];
  }
  c1/=c0;
  c2/=c0;
  c2-=c1*c1;
  if (c2>0.) c2=sqrt(c2);
  else       c2=0.;
  return c2;
}
  
double GetSlopeY(int i) { // classic linear fit y=b+a*(x-<x>);
  if (i<0 || i>10 || sumW[i]==0.) {
    printf("[myCovariance->GetSlopeX] wrong weight=%f or index %d\n", sumW[i], i);
    exit(-1);
  }
  double slope=(sumXY[i]/sumW[i]-GetAverageY(i)*GetAverageX(i))/GetSigmaY2(i);
  return slope;
}  
double GetSlopeYVariance() {
  double c;
  double c0=0.;
  double c1=0.;
  double c2=0.;
  for (int i=0; i<10; i++) {
    c=GetSlopeY(i);
    c0+=sumW[i];
    c1+=c*sumW[i];
    c2+=c*c*sumW[i];
  }
  c1/=c0;
  c2/=c0;
  c2-=c1*c1;
  if (c2>0.) c2=sqrt(c2);
  else       c2=0.;
  return c2;
}

double GetSlopeX0(int i) { // classic linear fit y=b+a*(x-<x>);
  if (i<0 || i>10 || sumW[i]==0.) {
    printf("[myCovariance->GetSlopeX0] wrong weight=%f or index %d\n", sumW[i], i);
    exit(-1);
  }
  double slope=sumXY[i]/sumW[i]/GetAverageX2(i);
  return slope;
}  
double GetSlopeX0Variance() {
  double c;
  double c0=0.;
  double c1=0.;
  double c2=0.;
  for (int i=0; i<10; i++) {
    c=GetSlopeX0(i);
    c0+=sumW[i];
    c1+=c*sumW[i];
    c2+=c*c*sumW[i];
  }
  c1/=c0;
  c2/=c0;
  c2-=c1*c1;
  if (c2>0.) c2=sqrt(c2);
  else       c2=0.;
  return c2;
}

double GetSlopeY0(int i) { // classic linear fit y=b+a*(x-<x>);
  if (i<0 || i>10 || sumW[i]==0.) {
    printf("[myCovariance->GetSlopeX0] wrong weight=%f or index %d\n", sumW[i], i);
    exit(-1);
  }
  double slope=sumXY[i]/sumW[i]/GetAverageY2(i);
  return slope;
}  
double GetSlopeY0Variance() {
  double c;
  double c0=0.;
  double c1=0.;
  double c2=0.;
  for (int i=0; i<10; i++) {
    c=GetSlopeY0(i);
    c0+=sumW[i];
    c1+=c*sumW[i];
    c2+=c*c*sumW[i];
  }
  c1/=c0;
  c2/=c0;
  c2-=c1*c1;
  if (c2>0.) c2=sqrt(c2);
  else       c2=0.;
  return c2;
}


  
};


char *module_name[7] = { 
  "HEC", 
  "EMEC", 
  "FCal025", 
  "FCal010", 
  "Beam",
  "Monitor",
  "Integrator"
};
std::ifstream indata;
int myLocation=0;
int displayWidth;
int displayHeight;
int canvasWidth=600;
int canvasHeight=1000;
char postscriptDirectorry[200];
int eventToCatch=-1;

///=============================================================================
void HiLumReadBareBones2014_05_23() {
///=============================================================================
miniEventHeader     MH;

//printf("point#0\n");
//MH.Print();

TGClient *TC = new TGClient();
displayWidth  = TC->GetDisplayWidth();
displayHeight = TC->GetDisplayHeight();
if (canvasWidth>displayWidth/4)    canvasWidth=displayWidth/4;
if (canvasHeight>displayHeight-50) canvasHeight=displayHeight-50;

printf("w=%d h=%d\n", displayWidth, displayHeight);

indata.open("location.dat", ios::in);
if(!indata.is_open()) {
  printf("[HiLumRead2014] location.dat, bye\n");
}
indata>>myLocation;
indata.close();
indata.clear();

if (myLocation==1) {
  printf("Running on the Sasha's laptop\n");
  sprintf(postscriptDirectorry, "plots2014_05/");
} else if (myLocation==0) {
  sprintf(postscriptDirectorry, "/raid04/HiLum/analysis2012/reports/MIP/spectra/");
  printf("Running on the Sasha's desktop\n");
} else {
  printf("Running at location#%d\n", myLocation);
  sprintf(postscriptDirectorry, " ");
}

TH1F   *dummy     = new TH1F("dummy",     "", 1025, -12.5, 512*12.5);
TH1F   *dummyRate = new TH1F("dummyRate", "", 1025, -12.5, 512*12.5);
TH1F   *summaryGainDummy     = new TH1F("summaryGainDummy",     "", 2000, 0.5, 2000.5);
TH1F   *summaryPositionDummy = new TH1F("summaryPositionDummy", "", 2000, 0.5, 2000.5);
summaryGainDummy->GetXaxis()->SetTitle("Run Index");
summaryGainDummy->GetYaxis()->SetTitle("GADC/Calibration");
summaryGainDummy->GetYaxis()->SetRangeUser(0., 1.25);
summaryGainDummy->GetXaxis()->SetRangeUser(-0.5, 0.5);
summaryPositionDummy->GetXaxis()->SetTitle("Run Index");
summaryPositionDummy->GetYaxis()->SetTitle("Beam position, EMEC balance");
summaryPositionDummy->GetYaxis()->SetRangeUser(-0.5, 0.5);
summaryPositionDummy->GetXaxis()->SetRangeUser(-0.5, 0.5);


dummy    ->SetLineStyle(2);
dummyRate->SetLineStyle(2);
dummyRate->GetXaxis()->SetTitle("Prehistory Rate, 10^{6}p/sec");

TLine  *myLine  = new TLine();
TLatex *myLatex = new TLatex();
myLatex->SetTextSize(0.06);
char   myText[500];
char   myTex1[500];
char   myTex2[500];


//--------------------------------------------------------------

miniBurstHeader     BH;
eventEssentials     evtEssentials;
cherenkovEssentials cherEssentials;

pedestalWindows  pedWindows[25];
moduleSignalsF   Channel[50];
Module           myModule[4];
scintMIPs        mip123;

double moduleHV[4];
int    blankModule;
bool   S123inBeam;

char   rootName[200];
char   treeName[200];
char   branchName[200];
char   titleH[200];

double time[515];
///-------------------------------------------------------------------------
int iqq=1;
int iqqOld=1;

for (int i=0; i<512; i++) {time[i]=12.5*i;} // set 512*12.5ns time grid for signals 
for (int iG=0; iG<2; iG++) { // set X and Y for module channels, based on HiLum LogBook records //
  int iCh, iMod;
  iMod=0; // HEC
  iCh=0; Channel[iG+2*iCh+10*iMod].X=-1.; Channel[iG+2*iCh+10*iMod].Y= 1.;
  iCh=1; Channel[iG+2*iCh+10*iMod].X= 1.; Channel[iG+2*iCh+10*iMod].Y= 1.;
  iCh=2; Channel[iG+2*iCh+10*iMod].X=-1.; Channel[iG+2*iCh+10*iMod].Y=-1.;
  iCh=3; Channel[iG+2*iCh+10*iMod].X= 1.; Channel[iG+2*iCh+10*iMod].Y=-1.;
  iMod=1; // EMEC
  iCh=0; Channel[iG+2*iCh+10*iMod].X=-1.; Channel[iG+2*iCh+10*iMod].Y= 1.;
  iCh=1; Channel[iG+2*iCh+10*iMod].X=-1.; Channel[iG+2*iCh+10*iMod].Y=-1.;
  iCh=2; Channel[iG+2*iCh+10*iMod].X= 1.; Channel[iG+2*iCh+10*iMod].Y= 1.;
  iCh=3; Channel[iG+2*iCh+10*iMod].X= 1.; Channel[iG+2*iCh+10*iMod].Y=-1.;
  iMod=2; // FCal-250
  iCh=0; Channel[iG+2*iCh+10*iMod].X= 1.; Channel[iG+2*iCh+10*iMod].Y= 1.;
  iCh=1; Channel[iG+2*iCh+10*iMod].X=-1.; Channel[iG+2*iCh+10*iMod].Y= 1.;
  iCh=2; Channel[iG+2*iCh+10*iMod].X= 1.; Channel[iG+2*iCh+10*iMod].Y=-1.;
  iCh=3; Channel[iG+2*iCh+10*iMod].X=-1.; Channel[iG+2*iCh+10*iMod].Y=-1.;
  iMod=3; // FCal-100
  iCh=0; Channel[iG+2*iCh+10*iMod].X= 1.; Channel[iG+2*iCh+10*iMod].Y= 1.;
  iCh=1; Channel[iG+2*iCh+10*iMod].X=-1.; Channel[iG+2*iCh+10*iMod].Y= 1.;
  iCh=2; Channel[iG+2*iCh+10*iMod].X= 1.; Channel[iG+2*iCh+10*iMod].Y=-1.;
  iCh=3; Channel[iG+2*iCh+10*iMod].X=-1.; Channel[iG+2*iCh+10*iMod].Y=-1.;
}

TCanvas *moduleCanvas   = new TCanvas("moduleCanvas",     "Module Signals",             0, 0, canvasWidth, canvasHeight);
TCanvas *moduleCanvas1  = new TCanvas("moduleCanvas1",    "Module Currents",  canvasWidth, 0, canvasWidth, canvasHeight);
moduleCanvas ->Connect("ProcessedEvent(Int_t,Int_t,Int_t,TObject*)", 0, 0, "MyExecuteEvent(Int_t,Int_t,Int_t,TObject*)");
moduleCanvas1->Connect("ProcessedEvent(Int_t,Int_t,Int_t,TObject*)", 0, 0, "MyExecuteEvent(Int_t,Int_t,Int_t,TObject*)");

TCanvas *spillCanvas    = new TCanvas("spillCanvas",      "Spill Structure",          2*canvasWidth, 0, canvasWidth, canvasHeight);
TCanvas *spillCanvas1   = new TCanvas("spillCanvas1",     "Average Spill Structure",  3*canvasWidth, 0, canvasWidth, canvasHeight);
spillCanvas ->Connect("ProcessedEvent(Int_t,Int_t,Int_t,TObject*)", 0, 0, "MyExecuteEvent(Int_t,Int_t,Int_t,TObject*)");
spillCanvas1->Connect("ProcessedEvent(Int_t,Int_t,Int_t,TObject*)", 0, 0, "MyExecuteEvent(Int_t,Int_t,Int_t,TObject*)");

moduleCanvas->Iconify();
moduleCanvas1->Iconify();

int MaxTracks=0; // maximum tracks seen by S123 using count2 method;
int runIndex=0;

double HGmax=0.;
double HGmin=0.;

int myRun0=850;
int myRun1=1200;

myRun0=1098;
myRun1=1098;

//myRun0=1099;
//myRun1=1099;

//myRun0=1024;
//myRun1=1024;

//myRun0=1093;
//myRun1=1098;

//myRun0=1041;
//myRun1=1041;

//myRun0=1019;
//myRun1=1019;

//myRun0=1117;
//myRun1=1117;

//myRun0=1098;
//myRun1=1099;

//myRun0=860;
//myRun1=860;

//myRun0=1050;
//myRun1=1500;

for (int myRun=myRun0; myRun<=myRun1; myRun++) { // runs looper
//*********************************************************************
int burstSet=-1;
MaxTracks=0;
// trying to open ROOT file with data trees
if (myLocation==1) {
  sprintf(rootName, "reports/ntuple/run%04dv20140404tree.root", myRun);
} else {
  sprintf(rootName, "/raid04/HiLum/analysis2012/reports/ntuples/run%04dv20140404tree.root", myRun);
}

TFile *f = new TFile(rootName,"READONLY");
if (!f->IsOpen()) continue;

///================================================================
/// good place to reset per-run histograms
///================================================================
//
printf("using file %s\n", rootName);
// extract an event TTree from file
sprintf(treeName, "run%04devents", myRun);
TTree *eventTree = (TTree*)f->Get(treeName);

// extract a spill TTree from file
sprintf(treeName, "run%04dspills", myRun);
TTree *burstTree = (TTree*)f->Get(treeName);

// extract a run TTree from file
sprintf(treeName, "run%04d", myRun);
TTree *runTree = (TTree*)f->Get(treeName);

Int_t eventEntries = (Int_t)eventTree->GetEntries();
Int_t burstEntries = (Int_t)burstTree->GetEntries();
Int_t runEntries   = (Int_t)runTree->GetEntries();
printf("\n\nEntries: Events=%d Bursts=%d Runs=%d\n", eventEntries, burstEntries, runEntries);
printf("\n\nIgnore that croaking about Bool_t, it works anyway!\n");
//  connect to Run branches

runTree -> SetBranchAddress("ScintMIPs",         &mip123.scintMIP[0]);
runTree -> SetBranchAddress("HVmod",             &moduleHV[0]);
runTree -> SetBranchAddress("AntennaModule",     &blankModule);
runTree -> SetBranchAddress("S123inBeam",        &S123inBeam);
runTree -> GetEntry(0);

//------------------------------------------------------------------------
//
// connect to Burst tree branches

burstTree -> SetBranchAddress("miniBirstHeader", &BH.burst);
burstTree -> SetBranchAddress("nLaps",           &cherEssentials.ns30);
burstTree -> SetBranchAddress("LapSum",          &cherEssentials.sum30[0]);
burstTree -> SetBranchAddress("LapFill",         &cherEssentials.lapFill[0]);
burstTree -> SetBranchAddress("LapBool",         &cherEssentials.lapBool[0]);
burstTree -> SetBranchAddress("LapActn",         &cherEssentials.lapActn[0]);


// connect to Event tree Branches
eventTree  -> SetBranchAddress("EventHeaderBranch", &MH.burst);

eventTree  -> SetBranchAddress("EventGADC",      &evtEssentials.GADCevent[0]);
eventTree  -> SetBranchAddress("MPICevent",      &evtEssentials.MPICevent[0]);

eventTree -> SetBranchAddress("firstGADCsample", &evtEssentials.firstGADCsample);
eventTree -> SetBranchAddress("bunchAction",     &evtEssentials.bunchAction[0]);

eventTree -> SetBranchAddress("scintFFT",        &evtEssentials.scintFFT[0]);

eventTree -> SetBranchAddress("mipCount",        &evtEssentials.mipCount[0]);
eventTree -> SetBranchAddress("mipCount2",       &evtEssentials.mipCount2[0]);
eventTree -> SetBranchAddress("mipCount3",       &evtEssentials.mipCount3[0]);

eventTree -> SetBranchAddress("isIsolated",      &evtEssentials.isIsolated[0]);
eventTree -> SetBranchAddress("isCluster0",      &evtEssentials.isCluster0[0]);
eventTree -> SetBranchAddress("isCluster1",      &evtEssentials.isCluster1[0]);
eventTree -> SetBranchAddress("inScope",         &evtEssentials.inScope[0]);


eventTree -> SetBranchAddress("HGoverflow",      &evtEssentials.HGoverflow[0]);
eventTree -> SetBranchAddress("LGoverflow",      &evtEssentials.LGoverflow[0]);

eventTree -> SetBranchAddress("RFclust0",        &evtEssentials.RFclust0[0]);
eventTree -> SetBranchAddress("RFclust1",        &evtEssentials.RFclust1[0]);

eventTree -> SetBranchAddress("miniModule",      &evtEssentials.miniChannel[0]);

for (int iMod=0; iMod<5; iMod++) {
  int iCh1=5;
  if (iMod==4) iCh1=4;
  for (int iCh=0; iCh<iCh1; iCh++) {
    int id=iCh+5*iMod; // pedestal windows data goes to special structure
    sprintf(branchName, "nPedWindows_mod%dchan%d", iMod, iCh);
    if (iqq<3) printf("setting branch address for pedWindows[%d], branchName=%s\n", id, branchName);

    eventTree -> SetBranchAddress(branchName, &pedWindows[id].nPedWindows);
    sprintf(branchName, "firstPedWindowSample_mod%dchan%d", iMod, iCh);
    eventTree -> SetBranchAddress(branchName, &pedWindows[id].firstSample[0]);
    sprintf(branchName, "lastPedWindowSample_mod%dchan%d",  iMod, iCh);
    eventTree -> SetBranchAddress(branchName, &pedWindows[id].lastSample[0]);

    for (int iG=0; iG<2; iG++) {
      id=iG+2*iCh+10*iMod;
      sprintf(branchName, "signal_mod%dchan%dgain%d", iMod, iCh, iG);
      if (iqq<3) printf("setting branch address for signals[%d], branchName=%s \n", id, branchName);
      eventTree -> SetBranchAddress(branchName, &Channel[id].ampI[0]);
      if (iMod != 4) { // for scintillators - different story!!!!
        sprintf(branchName, "current_mod%dchan%dgain%d", iMod, iCh, iG);
        eventTree -> SetBranchAddress(branchName, &Channel[id].ampF[0]);
      } else {         // FFT signal for PMTs not implemented yet!
      }
    }
  }
}

printf("\n\n");
printf(" Scintillator MIPs (FFT) ->");
for (int i=0; i<3; i++) {printf(" S%d=%.1f", i+1, mip123.scintMIP[i]);} printf("\n");
printf(" HV(V)->");
for (int iMod=0; iMod<4; iMod++) printf("%s=%.1f ", module_name[iMod], moduleHV[iMod]); printf("\n"); 

//if (moduleHV[1]==0.) { // HV EMEC not set
//  printf("\n\nHV EMEC is Zero, skipping this run\n\n");
//  continue;
//} 


printf(" blank module#%d\n", blankModule);
printf(" S123 in beam =%d\n\n", S123inBeam);
//==============================================================
int    burstEnt=0;
double myLapMax=0;


for (int iEnt=0; iEnt<eventEntries; iEnt++) {  //  event loop ///

//  printf("point#2\n");
//  MH.Print();

  if (iEnt==eventToCatch) iqq=1;
  evtEssentials.Clear();
  for (int i=0; i<50; i++) Channel[i].Clear();

//  printf("point#3\n");
//  MH.Print();


//**************************
  eventTree->GetEntry(iEnt);
//**************************

//   printf("\npoint#4\n");
//   MH.Print();



  double T0=MH.TimeGlobal;
  evtEssentials.CorrectGADCtail();

  if (iqq<3) { // print event RF bunch structure (filled, empty, cluster)
    spillCanvas  ->Show();
    spillCanvas1 ->Show();
    moduleCanvas ->Show();
    moduleCanvas1->Show();    
  }

  if (MH.burst != burstSet) { // first or new burst, load burst entry
    if (burstEnt<burstEntries) {
      printf("Loading data for burst#%d\n", MH.burst);
      do burstTree->GetEntry(burstEnt++); // load burst data until burst# is right
      while (BH.burst != MH.burst);

      if (BH.burst != MH.burst) {
        printf("Wrong burst number in Burst Tree BH.burst=%d\n", BH.burst);
        exit(-1);
      }
      burstSet=BH.burst;
      cherEssentials.Compute(); // compute what you need on the GADC record (laps!!! not bunches);

      if (iqq < 3) {
        printf("new spill!\n");
        printf("  push-> -1 for quick exit\n");
        printf("          0 - exit, but keep plots\n");
        printf("          1 - event-by-event\n");
        printf("         -N - run until bumping over evt#N\n");
        printf("          2 - non-stop spill show\n");
        printf("          3 - for quiet run-by-run\n");
        printf("          4 - for quiet all runs\n --> ");
        cin >> iqq;

        if (iqq==-1) exit(-1);
        if (iqq==0) break;
        if (iqq<0) {
          eventToCatch=-iqq;
          iqq=3;
        }
        if (iqq<3) {
          moduleCanvas ->Show();
          moduleCanvas1->Show();
          spillCanvas  ->Show();
          spillCanvas1 ->Show();
        } else {
          moduleCanvas ->Iconify();
          moduleCanvas1->Iconify();
          spillCanvas  ->Iconify();
          spillCanvas1 ->Iconify();
        }

        dummy->SetStats(0);
        dummy->SetTitle("Spill profile");
        dummy->GetXaxis()->SetRangeUser(0., lapTime*cherEssentials.ns30);
        dummy->GetXaxis()->SetTitle("Spill Time, ms");
        double rate=0.;
        double tLap=0.;
        int    nLapsShow=500; // how many laps to average for display
        for (int i=0; i<cherEssentials.ns30; i++) {
          rate+=cherEssentials.sum30[i];
          if (i>0 && i%nLapsShow==0) {
            rate/=nLapsShow;
            if (myLapMax<rate) myLapMax=rate;
            rate=0.;
          }
        }      
      
//      dummy->GetYaxis()->SetRangeUser(0., BH.ChGADC*cherEssentials.lapMax); // maximum protons/lap
        spillCanvas->Clear();
        spillCanvas->Divide(1,3);
        spillCanvas->cd(1);
        dummy->GetYaxis()->SetRangeUser(0., 1.1*BH.ChGADC*myLapMax); // maximum protons/lap
        dummy->GetYaxis()->SetTitle("<Rate> (500 laps) , protons/lap");
        sprintf(myText, "Spill profile for Run#%d spill#%d", myRun, burstSet);
        dummy->SetTitle(myText);
        dummy->DrawCopy();
        myLine->SetLineWidth(2);
        myLine->SetLineColor(kBlack);
        for (int i=0; i<cherEssentials.ns30; i++) {
          tLap+=i*lapTime;
          rate+=cherEssentials.sum30[i];
          if (i>0 && i%nLapsShow==0) {
            tLap/=nLapsShow;
            rate/=nLapsShow;
            myLine->DrawLine(tLap-nLapsShow*lapTime, BH.ChGADC*rate,
                             tLap+nLapsShow*lapTime, BH.ChGADC*rate);
            tLap=0.;
            rate=0.;
          }
        }

        spillCanvas->cd(2);
        dummy->GetYaxis()->SetRangeUser(-1.1, 1.1); // maximum protons/lap
        dummy->GetYaxis()->SetTitle("Balance, X");
        sprintf(myText, "Position for Run#%d spill#%d", myRun, burstSet);
        dummy->SetTitle(myText);
        dummy->DrawCopy();

        spillCanvas->cd(3);
        dummy->GetYaxis()->SetRangeUser(0., 2.); // maximum protons/lap
        dummy->GetYaxis()->SetTitle("GADC/S_{123}");
        sprintf(myText, "Run#%d spill#%d", myRun, burstSet);
        dummy->SetTitle(myText);
        dummy->DrawCopy();               
        spillCanvas->Update();

      }
    } else {
      printf("Burst Entry #%d is out of range: nBursts=%d\n", burstEnt, burstEntries);
      continue;
    }
  }

  int cluster=0;
  int nClusters=0;

  for (int iCl=0; iCl<5; iCl++) {
    if (evtEssentials.RFclust0[iCl]==-1) break;
    nClusters++;
  }
  int i1=0;
  int i2=0.;
  for (int i=0; i<40; i++) {
    if (evtEssentials.MPICevent[i]!=0.0) {
      i1=i;
      break;
    }
  }
  for (int i=39; i>-1; i--) {
    if (evtEssentials.MPICevent[i]!=0.0) {
      i2=i;
      break;
    }
  }

  if (!MH.filledSegment) {
    if (MH.isEmpty) { // fill nTrack=0 GADC spectra as noise

// good place to fill the electroinic n0ise spectra out-of-beam---

      if (iqq==1) {
        for (int i=0; i<40; i++) {
          printf("empty event, RFbunch#%2d inScope=%1d MPI.Cher=%.1f\n", i, evtEssentials.inScope[i], evtEssentials.MPICevent[i]);
        }        
      }
    }
    continue;
  }
  bool HGoverflow=false;
  bool LGoverflow=false; 
//---> check for miniModules overflow in High and Low gain MPI data
  for (int iMod=0; iMod<4; iMod++) {
    HGoverflow = HGoverflow || evtEssentials.HGoverflow[iMod];
    LGoverflow = LGoverflow || evtEssentials.LGoverflow[iMod];
  }
//----> Set conditions for signal shapes display. Here - HG overflow in ANY of the miniModules
  if (HGoverflow && iqq<3) {iqqOld=iqq; iqq=1;} 
//  else if (iqq!=3)  iqq=iqqOld;
//
  
// compute centers of gravity and some sums

  for (int iMod=0; iMod<4; iMod++) {
    if (iMod==blankModule) continue; // skip module at HV=0;
    int iG=1;
    if (evtEssentials.HGoverflow[iMod]) iG=0;
    cluster=0;
    myModule[iMod].Clear();
    
    for (int iClust=0; iClust<nClusters; iClust++) {
      if (evtEssentials.RFclust0[iClust]==-1) break; // no more clusters //
      for (int iCh=0; iCh<4; iCh++) {
        if (HGmax<evtEssentials.miniChannel[1+2*iCh+10*iMod+40*iClust]) 
            HGmax=evtEssentials.miniChannel[1+2*iCh+10*iMod+40*iClust];
        if (HGmin>evtEssentials.miniChannel[1+2*iCh+10*iMod+40*iClust]) 
            HGmin=evtEssentials.miniChannel[1+2*iCh+10*iMod+40*iClust];
          //---------------------------------------------------------------
        double a=evtEssentials.miniChannel[iG+2*iCh+10*iMod+40*iClust];
        double x=Channel[iG+2*iCh+10*iMod].X;
        double y=Channel[iG+2*iCh+10*iMod].Y;
        myModule[iMod].Add(a, x, y, iClust);
          //---------------------------------------------------------------
      } // channel loop
    } // clusters
    myModule[iMod].Compute(); 
  } // module loop
// now myModule[iMod].X contains <X> over all VALID events signals;
//     myModule[iMod].X[iCl] has <X> over each VALID event signal (cluster)


//  int myLap=evtEssentials.firstGADCsample/30;
  int myLap=MH.firstGated/30;

  if (iqq<3) {
//  printf("myLap=%d\n", myLap);

    spillCanvas->cd(1);
    myLatex->SetTextColor(kRed);
    myLatex->SetTextSize(0.04);
    myLatex->SetTextAlign(22);
    myLatex->DrawLatex(myLap*lapTime, BH.ChGADC*cherEssentials.sum30[myLap], "o");
    myLatex->SetTextAlign(11);
    myLatex->SetTextSize(0.06);
    spillCanvas->Update();
  }



  
  for (int iClust=0; iClust<nClusters; iClust++) { // remedy problematic spill run#1117spill#61 and such
    for (int j=evtEssentials.RFclust0[iClust]; j<=evtEssentials.RFclust1[iClust]; j++) {
       evtEssentials.isIsolated[j]=(evtEssentials.RFclust0[iClust]==evtEssentials.RFclust1[iClust]);
    }
  }
// for entire event!
  int iHEC=0;
  int iEMEC=1;
  int iFCal=2;
  if (iFCal==blankModule) iFCal=3;
  double qHEC =myModule[iHEC] .Q;
  double qEMEC=myModule[iEMEC].Q;
  double qFCal=myModule[iFCal].Q;

  double xHEC =myModule[iHEC] .X;
  double xEMEC=myModule[iEMEC].X;
  double xFCal=myModule[iFCal].X;

  double yHEC =myModule[iHEC] .Y;
  double yEMEC=myModule[iEMEC].Y;
  double yFCal=myModule[iFCal].Y;

  double evtMIPs=0.;
  double evtGADC=0.;


  for (int iClust=0; iClust<nClusters; iClust++) {
    for (int j=evtEssentials.RFclust0[iClust]; j<=evtEssentials.RFclust1[iClust]; j++) {
      if (MaxTracks<evtEssentials.mipCountX(j)) MaxTracks=evtEssentials.mipCountX(j);
      evtMIPs+=evtEssentials.mipCountX(j); // sum signal over filled bunches 
      evtGADC+=evtEssentials.GADCevent[j];

      if (evtEssentials.MPICevent[j]<-2500.) iqq=1; // good place to set up a trap for special 
                                                    // condition(s), requiring a clean event display

      double a1=evtEssentials.scintFFT[3*j]    /mip123.scintMIP[0];
      double a2=evtEssentials.scintFFT[3*j+1]  /mip123.scintMIP[1];
      double a3=evtEssentials.scintFFT[3*j+2]  /mip123.scintMIP[2];

      if (evtEssentials.isIsolated[j]) { // special treatment for isolated bunched
      }

    }
  }
  evtGADC*=BH.ChGADC;

//
  if (iqq<3 && evtMIPs>0) {

    spillCanvas  ->Show();
    spillCanvas1 ->Show();
    moduleCanvas ->Show();
    moduleCanvas1->Show();
 
    spillCanvas->cd(2);
    myLatex->SetTextColor(kBlack);
    myLatex->SetTextSize(0.04);
    myLatex->SetTextAlign(22);
    myLatex->DrawLatex(myLap*lapTime, xEMEC, "o");
    myLatex->SetTextColor(kRed);
    myLatex->DrawLatex(myLap*lapTime, xHEC,  "o");

    spillCanvas->cd(3);
    myLatex->SetTextColor(kBlack);
    myLatex->SetTextSize(0.04);
    myLatex->SetTextAlign(22);
    myLatex->DrawLatex(myLap*lapTime, evtGADC/evtMIPs, "o");
    myLatex->SetTextAlign(11);
    myLatex->SetTextSize(0.06);
    spillCanvas->Update();
  }

  
// cluster-by-cluster
  for (int iClust=0; iClust<nClusters; iClust++) {
    if (iFCal==blankModule) iFCal=3;
    double bunchMIPs=0.;
    double bunchGADC=0.;
    for (int j=evtEssentials.RFclust0[iClust]; j<=evtEssentials.RFclust1[iClust]; j++) {
      bunchMIPs+=evtEssentials.mipCount2[j];
      bunchGADC+=evtEssentials.GADCevent[j];
    }
    bunchGADC*=BH.ChGADC; // calibration
  } // clusters  

  if (iqq>2) continue; // skip the event display

  printf("Event Display---->\n");

  for (int iClust=0; iClust<nClusters; iClust++) {
    printf("\nCluster#%d %2d...%2d\n", iClust, evtEssentials.RFclust0[iClust], evtEssentials.RFclust1[iClust]);
    
    double bMIPs=0.;
    double bGADC=0.;
    for (int j=evtEssentials.RFclust0[iClust]; j<=evtEssentials.RFclust1[iClust]; j++) {
      printf(" j=%2d S123=(%5.0f %5.0f %5.0f)  GADC=%5.0f\n", j, (float)evtEssentials.mipCount[j], 
                                                                 (float)evtEssentials.mipCount2[j], 
                                                                 (float)evtEssentials.mipCount3[j], 
                                                                 (float)evtEssentials.GADCevent[j]);
      bMIPs+=evtEssentials.mipCountX(j);
      bGADC+=evtEssentials.GADCevent[j];
    }
    printf(" GADC=%6.0f S132=%.0f ", bGADC, bMIPs);
    for (int iMod=0; iMod<4; iMod++) {
      if (iMod==blankModule) continue; // skip module at HV=0;
      if (myModule[iMod].Qcl[iClust]>0.) {
         printf("%7s sum=%6.0f X=%5.2f Y=%5.2f \n", module_name[iMod], 
                 myModule[iMod].Qcl[iClust], myModule[iMod].Xcl[iClust], 
                                             myModule[iMod].Ycl[iClust]);
      } else {
         printf("%7s sum=%6.0f \n", module_name[iMod], myModule[iMod].Qcl[iClust]);
      }
    }
  } // clusters

  printf("\nSum over Event\n");
  for (int iMod=0; iMod<4; iMod++) {
    if (iMod==blankModule) continue; // skip module at HV=0;
    if (myModule[iMod].Q>0.) { 
      printf("  %7s sum=%6.0f X=%5.2f Y=%5.2f\n", module_name[iMod], 
                 myModule[iMod].Q, myModule[iMod].X, myModule[iMod].Y);
    } else {  
      printf("  %7s sum=%6.0f\n", module_name[iMod], myModule[iMod].Q);
    }
  }
  printf("\n");

  moduleCanvas->Clear();
  moduleCanvas->Divide(1,4);
  moduleCanvas1->Clear();
  sprintf(myText, "Reconstructed Modules and S123 for Run#%d spill#%d", myRun, burstSet);
  moduleCanvas1->SetTitle(myText);
  moduleCanvas1->Divide(1,4);
  int iZone=1;
  myLatex->SetTextColor(kBlack);

  for (int iMod=0; iMod<4; iMod++) { // do not go into scintillators iMod=4;
    
    if (iMod==blankModule) continue; // skip module at HV=0;

    int iG=1;
    int iChD=4;
    if (evtEssentials.HGoverflow[iMod]) iG=0;

    int id=iG+2*iChD+10*iMod;
    double amin=0.;
    double amax=0.;
    double aminF=0.;
    double amaxF=0.;
    for (int i=0; i<512; i++) {
//      printf("id=%2d is=%3d ampI=%5.1f ampF=%5.1f\n", id, i, Channel[id].ampI[i], Channel[id].ampF[i]);
      if (amin >Channel[id].ampI[i]) amin =Channel[id].ampI[i]; 
      if (amax <Channel[id].ampI[i]) amax =Channel[id].ampI[i];
      if (aminF>Channel[id].ampF[i]) aminF=Channel[id].ampF[i]; 
      if (amaxF<Channel[id].ampF[i]) amaxF=Channel[id].ampF[i];
    }
    double dd=amax -amin;
    double df=amaxF-aminF;

    double gain=1.;
    if (iG==0) gain=10.;

    int id1=iChD+5*iMod;
    int nPedWindows=pedWindows[id1].nPedWindows;
    int is0=pedWindows[id1].firstSample[0];
    int is1=pedWindows[id1].lastSample[nPedWindows-1];

    moduleCanvas->cd(iZone);
    dummy->GetYaxis()->SetRangeUser(amin-0.05*dd, amax+0.05*dd);
    dummy->GetYaxis()->SetTitle("MPI ADC, mV");
    dummy->GetXaxis()->SetRangeUser(-12.5, 512*12.5);
    dummy->GetXaxis()->SetTitle("Time, ns");
    sprintf(titleH, "%s gain#%d", module_name[iMod], iG);
    dummy->SetTitle(titleH);
    dummy->SetStats(0);
    dummy->DrawCopy();

    myLine->SetLineWidth(13);
    myLine->SetLineColor(kYellow);
    for (int iWin=0; iWin<nPedWindows; iWin++) {
      myLine->DrawLine(time[pedWindows[id1].firstSample[iWin]], 0., 
                       time[pedWindows[id1].lastSample[iWin]],  0.);
    }
    myLine->SetLineColor(kGray);
    for (int i=0; i<40; i++) {
      int is=MH.firstGated+i;
      if (cherEssentials.lapBool[is%30]) {
        myLine->DrawLine(165*i, amin-0.05*dd, 165*i, amax+0.05*dd);
      }
    }

    myLine->SetLineWidth(2);
    myLine->SetLineColor(kBlack);
    for (int i=1; i<512; i++) {
      myLine->DrawLine(time[i-1], Channel[id].ampI[i-1], time[i], Channel[id].ampI[i] ); 
    }

    
    moduleCanvas1->cd(iZone);
    dummy->GetYaxis()->SetRangeUser(aminF-0.05*df, amaxF+0.05*df);
    sprintf(titleH, "%s gain#%d", module_name[iMod], iG);
    dummy->GetYaxis()->SetTitle("Current, #muA");
    dummy->GetXaxis()->SetTitle("Time, ns");
    dummy->SetTitle(titleH);
    dummy->DrawCopy();

    myLine->SetLineWidth(13);
    myLine->SetLineColor(kYellow);
    for (int iWin=0; iWin<nPedWindows; iWin++) {
      myLine->DrawLine(time[pedWindows[id1].firstSample[iWin]], 0., 
                       time[pedWindows[id1].lastSample[iWin]],  0.);
    }

    myLine->SetLineWidth(2);
    myLine->SetLineColor(kBlack);
    for (int i=1; i<is1-is0; i++) { // FFT reconstracted signal has different sample indexing!
      myLine->DrawLine(time[is0+i-1], Channel[id].ampF[i-1], time[is0+i], Channel[id].ampF[i] ); 
    }
    

    myLine->SetLineWidth(1);
    for (int iCh1=0; iCh1<4; iCh1++) {
      myLine->SetLineColor(iCh1+1);
      int id2=iG+2*iCh1+10*iMod;
      for (int i=1; i<is1-is0; i++) {
        myLine->DrawLine(time[is0+i-1], Channel[id2].ampF[i-1], time[is0+i], Channel[id2].ampF[i] ); 
      }
    }

    for (int iClust=0; iClust<nClusters; iClust++) {
      sprintf(myText, "Q=%.0ffC", myModule[iMod].Qcl[iClust]);
      myLatex->DrawLatex(165.*evtEssentials.RFclust0[iClust], amaxF+0.04*df, myText);
      if (myModule[iMod].Qcl[iClust]>0.) {
         sprintf(myText, "X=%.2f", myModule[iMod].Xcl[iClust]);
         myLatex->DrawLatex(165.*evtEssentials.RFclust0[iClust], amaxF-0.03*df, myText);
         sprintf(myText, "Y=%.2f", myModule[iMod].Ycl[iClust]);
         myLatex->DrawLatex(165.*evtEssentials.RFclust0[iClust], amaxF-0.10*df, myText);
      }
      
    }
    

    iZone++;

  }


// GADC/MPI-C response //
  double GADCmin=0.;
  double GADCmax=0.;
  double ddGADC=0.;
  for (int i=0; i<40 && 165.*(i+1)<6400.; i++) {    
    if (GADCmin>evtEssentials.GADCevent[i]) GADCmin=evtEssentials.GADCevent[i];
    if (GADCmax<evtEssentials.GADCevent[i]) GADCmax=evtEssentials.GADCevent[i];
  }
  ddGADC=GADCmax-GADCmin;
  dummy->GetYaxis()->SetRangeUser(GADCmin-0.05*ddGADC, GADCmax+0.05*ddGADC);
  moduleCanvas ->cd(iZone);
  dummy->SetTitle("Gated Cherenkov");
  dummy->GetYaxis()->SetTitle("Gated ADC, counts");
  dummy->GetXaxis()->SetTitle("Time, ns");
  dummy->DrawCopy();
  myLine->SetLineColor(kBlue-10);
  myLine->SetLineWidth(2);
  
  printf("********************************\n");
  printf("*  Cherenkov and S123 samples  *\n");
  printf("********************************\n");
  printf("* 'in'  = sample within the FFT reconstruction range *\n");
  printf("* 'act' = 0 for useless sample *\n");
  printf("* 'act' = 1 for signal sample  *\n");
  printf("* 'act' = 2 for noise sample   *\n");
  printf("********************************\n");
   
  for (int i=0; i<40 && 165.*(i+1)<6300.; i++) {
    printf("#%2d T=%6.0fns GADC=%7.1f in=%d", 
            i, 165.*i, evtEssentials.GADCevent[i], evtEssentials.inScope[i]); 
    if (evtEssentials.inScope[i]) {
      double a1=evtEssentials.scintFFT[3*i]    /mip123.scintMIP[0];
      double a2=evtEssentials.scintFFT[3*i+1]  /mip123.scintMIP[1];
      double a3=evtEssentials.scintFFT[3*i+2]  /mip123.scintMIP[2];
      printf("  MPIC=%8.1f ", evtEssentials.MPICevent[i]/33.);
      printf("  tracks: S123=(%5.1f %5.1f %5.1f) N=%3d  ", a1, a2, a3, evtEssentials.mipCount2[i]);
      printf(" act=%d", evtEssentials.bunchAction[i]);
      if (a3<0.25 && a1>0.25 && a2>0.25 && sqrt(a1*a2)>1.75) printf("odd");
      if (evtEssentials.MPICevent[i]<-1500.) printf(" NEGATIVE");
    }
    printf("\n");
    myLine->DrawLine(165.*i,     evtEssentials.GADCevent[i], 
                     165.*(i+1), evtEssentials.GADCevent[i]);
    if (i>0) {
      myLine->DrawLine(165.*i, evtEssentials.GADCevent[i-1], 
                       165.*i, evtEssentials.GADCevent[i]);
    }
  }
  myLine->SetLineColor(kBlue);
  myLine->SetLineWidth(2);
  for (int i=0; i<40 && 165.*(i+1)<6400.; i++) {
    myLine->DrawLine(165.*i,     evtEssentials.GADCevent1[i], 
                     165.*(i+1), evtEssentials.GADCevent1[i]);
    if (i>0) {
      myLine->DrawLine(165.*i, evtEssentials.GADCevent1[i-1], 
                       165.*i, evtEssentials.GADCevent1[i]);
    }
  }
  myLine->SetLineColor(kRed);
  for (int i=0; i<39 && 165.*(i+1)<6400.; i++) {
    if (!evtEssentials.inScope[i]) continue;
    myLine->DrawLine(165.*i,     evtEssentials.MPICevent[i]/33., 
                     165.*(i+1), evtEssentials.MPICevent[i]/33.);
    if (i>0) {
     if (!evtEssentials.inScope[i-1]) continue;
      myLine->DrawLine(165.*i, evtEssentials.MPICevent[i-1]/33., 
                       165.*i, evtEssentials.MPICevent[i]/33.);
    }
  }
  
  for (int iCl=0; iCl<nClusters; iCl++) { // show GADC track "counting" on the plot
    int ir0=evtEssentials.RFclust0[iCl];
    int ir1=evtEssentials.RFclust1[iCl];
    double sumGADC=0.;
    for (int ir=ir0; ir<ir1+1; ir++) sumGADC+=evtEssentials.GADCevent[ir];
    double tracksGADC=sumGADC*BH.ChGADC;
    sprintf(myText, "%.0f", tracksGADC);
    myLatex->DrawLatex(82.5*(ir0+ir1), GADCmax+0.07*ddGADC, myText);
  }

// S123 track counting
  double mipCountMAX=3.;
  for (int i=0; i<40 && 165.*(i+1)<6400.; i++) { 
    if (!evtEssentials.inScope[i]) continue;
    if (mipCountMAX<evtEssentials.mipCount[i])  mipCountMAX=evtEssentials.mipCount[i];
    if (mipCountMAX<evtEssentials.mipCount2[i]) mipCountMAX=evtEssentials.mipCount2[i];
    if (mipCountMAX<evtEssentials.mipCount3[i]) mipCountMAX=evtEssentials.mipCount3[i];
  }    
  dummy->GetYaxis()->SetRangeUser(0., 1.05*mipCountMAX);
  moduleCanvas1 ->cd(iZone);
  dummy->SetTitle("S123 track counting");
  dummy->GetYaxis()->SetTitle("Tracks");
  dummy->DrawCopy();

  for (int i=0; i<39 && 165.*(i+1)<6400.; i++) {
    if (!evtEssentials.inScope[i]) continue;
    myLine->SetLineColor(kBlack);
    myLine->SetLineWidth(2);
    myLine->DrawLine(165.*i,     evtEssentials.mipCount[i], 
                     165.*(i+1), evtEssentials.mipCount[i]);
    if (i>0) {
      if (!evtEssentials.inScope[i-1]) continue;
      myLine->DrawLine(165.*i, evtEssentials.mipCount[i-1], 
                       165.*i, evtEssentials.mipCount[i]);
    }

    double a1=evtEssentials.scintFFT[3*i]    /mip123.scintMIP[0];
    double a2=evtEssentials.scintFFT[3*i+1]  /mip123.scintMIP[1];
    double a3=evtEssentials.scintFFT[3*i+2]  /mip123.scintMIP[2];
    myLine->SetLineColor(kGreen);
    myLine->SetLineWidth(5);
    myLine->DrawLine(165.*i,     a1,   165.*(i+1), a1);
    myLine->DrawLine(165.*i,     a2,   165.*(i+1), a2);
    myLine->DrawLine(165.*i,     a3,   165.*(i+1), a3);

    if (a3<0.25 && a1>0.25 && a2>0.25 && sqrt(a1*a2)>1.75) {
//      printf("*\n");
      myLatex->SetTextSize(0.08);
      myLatex->DrawLatex(165.*i, 0.0, "#Uparrow"); 
      myLatex->SetTextSize(0.06);
    } else {
//      printf(" \n");
    }

  }

  myLine->SetLineWidth(2);
  myLine->SetLineColor(kRed);
  for (int i=0; i<39 && 165.*(i+1)<6400.; i++) {
    if (!evtEssentials.inScope[i]) continue;
    myLine->DrawLine(165.*i-10.,     evtEssentials.mipCount2[i], 
                     165.*(i+1)-10., evtEssentials.mipCount2[i]);
    if (i>0) {
      if (!evtEssentials.inScope[i-1]) continue;
      myLine->DrawLine(165.*i-10., evtEssentials.mipCount2[i-1], 
                       165.*i-10., evtEssentials.mipCount2[i]);
    }
  }

  for (int i=0; i<39 && 165.*(i+1)<6400.; i++) {
    if (!evtEssentials.inScope[i]) continue;
    myLine->SetLineColor(kBlue);
    myLine->DrawLine(165.*i+10.,     evtEssentials.mipCount3[i], 
                     165.*(i+1)+10., evtEssentials.mipCount3[i]);
    if (i>0) {
      if (!evtEssentials.inScope[i-1]) continue;
      myLine->DrawLine(165.*i+10., evtEssentials.mipCount3[i-1], 
                       165.*i+10., evtEssentials.mipCount3[i]);
    }
  }
  
  for (int iCl=0; iCl<nClusters; iCl++) {
    int ir0=evtEssentials.RFclust0[iCl];
    int ir1=evtEssentials.RFclust1[iCl];
    int sumS123=0;
    for (int ir=ir0; ir<ir1+1; ir++) sumS123+=evtEssentials.mipCount2[ir];
    sprintf(myText, "%d", sumS123);
    myLatex->DrawLatex(82.5*(ir0+ir1), 1.07*mipCountMAX, myText);
  }


  moduleCanvas ->Update();
  moduleCanvas1->Update();

  myLine->SetLineColor(kBlack);
  
  if (iqq < 2) {
    printf("Next Event!\n");
    printf("  push-> -1 for quick exit\n");
    printf("          0 - exit, but keep plots\n");
    printf("          1 - event-by-event\n");
    printf("         -N - run until bumping over evt#N\n");
    printf("          2 - non-stop spill show\n");
    printf("          3 - for quiet run-by-run\n");
    printf("          4 - for quiet all runs\n --> ");
    cin >> iqq;
    if (iqq==-1) exit(-1);
    if (iqq==0) break;
    if (iqq<0) {
      eventToCatch=-iqq;
      iqq=3;
    }
  }
} // end of run events looper
//========================================================================
// do some histogram output, save run plots and data
//========================================================================
if (iqq==0) break; // exit run loop
if (iqq < 4) {
  printf("new run!\n");
  printf("  push-> -1 for quick exit\n");
  printf("          0 - exit, but keep plots\n");
  printf("          1 - event-by-event\n");
  printf("         -N - run until bumping over evt#N\n");
  printf("          2 - non-stop spill show\n");
  printf("          3 - for quiet run-by-run\n");
  printf("          4 - for quiet all runs\n --> ");
  cin >> iqq;
  if (iqq==-1) exit(-1);
  if (iqq==0) break;  // exit run loop
  if (iqq<0) {
    eventToCatch=-iqq;
    iqq=3;
  }
}
f->Close();

runIndex++;

} // end of one run loop
//========================================================================
// do some histogram "global" output, save run plots and data
//========================================================================

}

// function to "pop-up" a pad from canvas
void MyExecuteEvent(Int_t event, Int_t px, Int_t py, TObject *sel) {
  if(sel) {
    //
  }
  TCanvas *c = (TCanvas *) gTQSender;

  if ( (event == kButton1Double) ) {
    char cname[256];
    TPad *pad = c->Pick(px, py, 0);
    //printf("d-clicked on pad %s #%d\n", pad->GetName(), pad->GetNumber());
    if(pad == 0) {
      std::cout << "Can;t get pad at " << px << " " << py << std::endl;
      return;
    }
    if(pad == c) {
      std::cout << "You have selected parent object" << std::endl;
      return;
    }
    float xyratio =  float(pad->GetAbsWNDC()*pad->GetCanvas()->GetWindowWidth()) / float(pad->GetAbsHNDC()*pad->GetCanvas()->GetWindowHeight());
    //std::cout << " xyratio: " << xyratio << std::endl;
    TPad *pad_new = (TPad *) pad->Clone();
    pad_new->SetPad(0.0, 0.0, 1.0, 1.0);
    //pad_new -> SetEditable(true);

    sprintf(cname,"%s_%d",pad->GetTitle(), (int)time(0));
    int iY=600;
    int iX=(int)(600*xyratio);
    if (iX<600) {
      iY=iY*600/iX;
      iX=600;
    }
    TCanvas *c1 = new TCanvas(cname, cname, iX, iY);
    c1->cd();
    if (cname[0]=='S' && cname[1]=='p' && cname[2]=='e' && cname[3]=='c') {
      gStyle->SetOptStat(111110);
    } 
    pad_new->SetEditable(kTRUE);
    pad_new->Draw();    
    pad_new->Modified();
    pad_new->Update();
//    pad_new->SetEditable(kFALSE);
//     c1->RaiseWindow();
//     c1->Modified();
//     c1->Update();
    c1->SetWindowSize(iX+20, iY+40);
    c1->Update();
    gSystem->ProcessEvents();
  }
}
/*


        sprintf(treeName, "run%04d", myRun);
        runTree = new TTree(treeName, treeName);
        runTree -> Branch("S123inBeam", RH->GetS123pointer(),                 "S123inBeam/B");
        runTree -> Branch("HVmod",      RH->GetHVsetPointer(),                "HVmod[4]/D");
        runTree -> Branch("AntennaModule",  RH->GetAntennaPointer(),          "AntennaModule/I");
        runTree -> Branch("CryoX",      RH->GetCryoXPointer(),                "CryoX[3]/D");
        runTree -> Branch("CryoPos",    RH->GetCryoPositionPointer(),         "CryoPos[3]/D");
        runTree -> Branch("ScintMIPs",  myMIPcounter->GetMIPaveragePointer(), "ScintMIPs[3]/D");
       

        sprintf(treeName, "run%04dspills", myRun);
        spillsTree = new TTree(treeName, treeName);
        spillsTree -> Branch("miniBirstHeader", &miniBH.burst,                "burst/I:ChMPIF/F:ChGADC/F");
        spillsTree -> Branch("nLaps",   CherBuffer->GetNlapsPointer(),        "nLaps/I");
        spillsTree -> Branch("LapSum",  CherBuffer->GetLapsSumPointer(),      "lapSum[nLaps]/F"); // all filled bunches, summed over laps
        spillsTree -> Branch("LapFill", CherBuffer->GetFilledSumPointer(),    "lapFill[30]/D"); // accumulated amplitudes for 30 bunches
        spillsTree -> Branch("LapBool", CherBuffer->GetFilledSamplePointer(), "lapBool[30]/B"); // true for filled, false for empty
        spillsTree -> Branch("LapActn", ProcStatus->GetBunchActionPointer(),  "lapActn[30]/I"); // 0 - no action 1 - signal 2 - ped
        spillsTree -> Branch("HVbeam",  RH->GetHVbeamPointer(),               "HVbeam[4]/D");


        sprintf(treeName, "run%04devents", myRun);
        eventsTree = new TTree(treeName, treeName);
        eventsTree -> Branch("EventHeaderBranch", &miniEH.burst, 
     "burst/I:event:burstEvent:scaler10mhz:scaler6mhz:firstGated:isEmpty/O:notEmpty:S123inBeam:filledSegment:RFperiod/D:RFfreq:RFamp:TimeGlobal:SignalFactor");// branch for event header;
        eventsTree -> Branch("EventGADC",       &evtEssentials.GADCevent[0],    "GADCevent[40]/F");
        eventsTree -> Branch("MPICevent",       &evtEssentials.MPICevent[0],    "MPICevent[40]/F");
        eventsTree -> Branch("firstGADCsample", &evtEssentials.firstGADCsample, "firstGADCsample/I");
        eventsTree -> Branch("bunchAction",     &evtEssentials.bunchAction[0],  "bunchAction[40]/I");

        eventsTree -> Branch("scintFFT",        &evtEssentials.scintFFT[0],     "scintFFT[120]/F");

        eventsTree -> Branch("mipCount",        &evtEssentials.mipCount[0],     "mipCount[40]/I");
        eventsTree -> Branch("mipCount2",       &evtEssentials.mipCount2[0],    "mipCount2[40]/I");
        eventsTree -> Branch("mipCount3",       &evtEssentials.mipCount3[0],    "mipCount3[40]/I");

        eventsTree -> Branch("isIsolated",      &evtEssentials.isIsolated[0],   "isIsolated[40]/B");
        eventsTree -> Branch("isCluster0",      &evtEssentials.isCluster0[0],   "isCluster0[40]/B");
        eventsTree -> Branch("isCluster1",      &evtEssentials.isCluster1[0],   "isCluster1[40]/B");
        eventsTree -> Branch("inScope",         &evtEssentials.inScope[0],      "inScope[40]/B");
       

        eventsTree -> Branch("HGoverflow",      &evtEssentials.HGoverflow[0],   "HGoverflow[4]/B");
        eventsTree -> Branch("LGoverflow",      &evtEssentials.LGoverflow[0],   "LGoverflow[4]/B");

        eventsTree -> Branch("RFclust0",        &evtEssentials.RFclust0[0],     "RFclust0[5]/I");
        eventsTree -> Branch("RFclust1",        &evtEssentials.RFclust1[0],     "RFclust1[5]/I");
  
        eventsTree -> Branch("miniModule",      &evtEssentials.miniChannel[0],   "miniChannel[200]/F");


        for (int iMod=0; iMod<5; iMod++) {
          TModule *Module = &Setup[iMod];
          int iChan1=5;
          if (iMod==4) iChan1=4; // no sum for counters
          for (int iChan=0; iChan<iChan1; iChan++) {
            TChannel *Channel = Module->GetChannel(iChan);
            for (int iGain=0; iGain<2; iGain++) {
              TSignal *Signal = Channel->GetSignalLong(iGain);
              if (iGain==0) { // Pedestal Windows are same for both gains, do it once
                sprintf(branchName, "nPedWindows_mod%dchan%d", iMod, iChan);
                eventsTree -> Branch(branchName, Signal->GetNPedWindowsPointer(), "nPedWindows/I");

                sprintf(branchName, "firstPedWindowSample_mod%dchan%d", iMod, iChan);
                eventsTree -> Branch(branchName, Signal->GetFirstSamplePointer(), "firstSample[25]/I");

                sprintf(branchName, "lastPedWindowSample_mod%dchan%d",  iMod, iChan);
                eventsTree -> Branch(branchName, Signal->GetLastSamplePointer(),  "lastSample[25]/I");
              }

              sprintf(branchName, "signal_mod%dchan%dgain%d", iMod, iChan, iGain); // MPI ADC signals, mV
              eventsTree -> Branch(branchName, Signal->GetAmpIFpoint(), "ampIF[512]/F"); 
              if (iMod != 4) {
                sprintf(branchName, "current_mod%dchan%dgain%d", iMod, iChan, iGain); // Module currents, mA
                eventsTree -> Branch(branchName, Signal->GetAmpFinalFpoint(), "signalFinalF[512]/F");
              } else {
              }
            }
          }
        }
*/