#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 <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
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

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 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 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
///////////////////////////////////////////////////////////////////////////////////
};
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 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 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 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 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;
}  


  
};


char *module_name[7] = { 
  "HEC", 
  "EMEC", 
  "FCal025", 
  "FCal010", 
  "Beam",
  "Monitor",
  "Integrator"
};
std::ifstream indata;
int myLocation=0;

///=============================================================================
void HiLumRead2014_04_29() {
///=============================================================================
indata.open("location.dat", ios::in);
if(!indata.is_open()) {
  printf("[HiLumRead2014] location.dat, bye\n");
}
indata>>myLocation;
if (myLocation==0) {
  printf("Running on the Sasha's laptop\n");
} else {
  printf("Running at location#%d\n", myLocation);
}

TH1F   *dummy     = new TH1F("dummy",     "", 1025, -12.5, 512*12.5);
TH1F   *dummyRate = new TH1F("dummyRate", "", 1025, -12.5, 512*12.5);

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

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

TH2F XdiffHE[10];
TH2F XdiffFE[10];
TH2F YdiffHE[10];
TH2F YdiffFE[10];

TH1F CherRateHE[10];
TH1F S123RateHE[10];

TH1F dX_HE[10];
TH1F dY_HE[10];

TH2F XdiffHEevt[10];
TH2F XdiffFEevt[10];
TH2F YdiffHEevt[10];
TH2F YdiffFEevt[10];

TH1F CherRateHEevt[10];
TH1F S123RateHEevt[10];

TH1F dX_HEevt[10];
TH1F dY_HEevt[10];

//double AminRange[10]={ 100.,  200.,  300.,  500.,  700.,
//                      1000., 2000., 3000., 5000., 7000.} ;
double AminRange[5]={  200.,  500.,  1000., 3000., 7000.};
int nRanges=5;

for (int i=0; i<nRanges; i++) {
  sprintf(myTex1, "%.0f<Amin", AminRange[i]);
  if (i<nRanges-1) {
    sprintf(myTex2, "<%.0f", AminRange[i+1]);
  } else {
    sprintf(myTex2, " ");
  }
  sprintf(myText, "Bunch X HEC-EMEC, %s%s", myTex1, myTex2);
  XdiffHE[i].SetNameTitle(myText, myText);
  XdiffHE[i].GetXaxis()->SetTitle("EMEC X, a.u.");
  XdiffHE[i].GetYaxis()->SetTitle("HEC X, a.u.");
  sprintf(myText, "Bunch Y HEC-EMEC, %s%s", myTex1, myTex2);
  YdiffHE[i].SetNameTitle(myText, myText);
  YdiffHE[i].GetXaxis()->SetTitle("EMEC Y, a.u.");
  YdiffHE[i].GetYaxis()->SetTitle("HEC Y, a.u.");
  sprintf(myText, "Bunch X FCal-EMEC, %s%s", myTex1, myTex2);
  XdiffFE[i].SetNameTitle(myText, myText);
  XdiffFE[i].GetXaxis()->SetTitle("EMEC X, a.u.");
  XdiffFE[i].GetYaxis()->SetTitle("FCal X, a.u.");
  sprintf(myText, "Bunch Y FCal-EMEC, %s%s", myTex1, myTex2);
  YdiffFE[i].SetNameTitle(myText, myText);
  YdiffFE[i].GetXaxis()->SetTitle("EMEC Y, a.u.");
  YdiffFE[i].GetYaxis()->SetTitle("FCal Y, a.u.");
  XdiffHE[i].SetBins(51, -1.02, 1.02, 51, -1.02, 1.02);
  XdiffFE[i].SetBins(51, -1.02, 1.02, 51, -1.02, 1.02);
  YdiffHE[i].SetBins(51, -1.02, 1.02, 51, -1.02, 1.02);
  YdiffFE[i].SetBins(51, -1.02, 1.02, 51, -1.02, 1.02);

  sprintf(myText, "Event X HEC-EMEC, %s%s", myTex1, myTex2);
  XdiffHEevt[i].SetNameTitle(myText, myText);
  XdiffHEevt[i].GetXaxis()->SetTitle("EMEC X, a.u.");
  XdiffHEevt[i].GetYaxis()->SetTitle("HEC X, a.u.");
  sprintf(myText, "Event Y HEC-EMEC, %s%s", myTex1, myTex2);
  YdiffHEevt[i].SetNameTitle(myText, myText);
  YdiffHEevt[i].GetXaxis()->SetTitle("EMEC Y, a.u.");
  YdiffHEevt[i].GetYaxis()->SetTitle("HEC Y, a.u.");
  sprintf(myText, "Event X FCal-EMEC, %s%s", myTex1, myTex2);
  XdiffFEevt[i].SetNameTitle(myText, myText);
  XdiffFEevt[i].GetXaxis()->SetTitle("EMEC X, a.u.");
  XdiffFEevt[i].GetYaxis()->SetTitle("FCal X, a.u.");
  sprintf(myText, "Event Y FCal-EMEC, %s%s", myTex1, myTex2);
  YdiffFEevt[i].SetNameTitle(myText, myText);
  YdiffFEevt[i].GetXaxis()->SetTitle("EMEC Y, a.u.");
  YdiffFEevt[i].GetYaxis()->SetTitle("FCal Y, a.u.");
  XdiffHEevt[i].SetBins(51, -1.02, 1.02, 51, -1.02, 1.02);
  XdiffFEevt[i].SetBins(51, -1.02, 1.02, 51, -1.02, 1.02);
  YdiffHEevt[i].SetBins(51, -1.02, 1.02, 51, -1.02, 1.02);
  YdiffFEevt[i].SetBins(51, -1.02, 1.02, 51, -1.02, 1.02);
  XdiffHEevt[i].SetMarkerColor(kRed);
  YdiffHEevt[i].SetMarkerColor(kRed);
  XdiffFEevt[i].SetMarkerColor(kRed);
  YdiffFEevt[i].SetMarkerColor(kRed);

  sprintf(myText, "Bunch dX HEC-EMEC, %s%s", myTex1, myTex2);
  dX_HE[i].SetNameTitle(myText, myText);
  dX_HE[i].GetXaxis()->SetTitle("X_{HEC}-X_{EMEC}, a.u.");
  dX_HE[i].SetBins(101, -2.02, 2.02);
  sprintf(myText, "Bunch dY HEC-EMEC, %s%s", myTex1, myTex2);
  dY_HE[i].SetNameTitle(myText, myText);
  dY_HE[i].GetXaxis()->SetTitle("Y_{HEC}-Y_{EMEC}, a.u.");
  dY_HE[i].SetBins(101, -2.02, 2.02);

  sprintf(myText, "Event dX HEC-EMEC, %s%s", myTex1, myTex2);
  dX_HEevt[i].SetNameTitle(myText, myText);
  dX_HEevt[i].GetXaxis()->SetTitle("X_{HEC}-X_{EMEC}, a.u.");
  dX_HEevt[i].SetBins(101, -2.02, 2.02);
  sprintf(myText, "Event dY HEC-EMEC, %s%s", myTex1, myTex2);
  dY_HEevt[i].SetNameTitle(myText, myText);
  dY_HEevt[i].GetXaxis()->SetTitle("Y_{HEC}-Y_{EMEC}, a.u.");
  dY_HEevt[i].SetBins(101, -2.02, 2.02);
  dX_HEevt[i].SetLineColor(kRed);
  dY_HEevt[i].SetLineColor(kRed);

  sprintf(myText, "Cherenkov rate HEC-EMEC, %s%s", myTex1, myTex2);
  CherRateHE[i].SetNameTitle(myText, myText);
  CherRateHE[i].GetXaxis()->SetTitle("Cherenkov tracks per bunch");
  CherRateHE[i].SetBins(400, -100., 3900.);
    
  sprintf(myText, "Cherenkov rate HEC-EMEC,sum, %s%s", myTex1, myTex2);
  CherRateHEevt[i].SetNameTitle(myText, myText);
  CherRateHEevt[i].GetXaxis()->SetTitle("Cherenkov tracks per event");
  CherRateHEevt[i].SetBins(400, -100., 3900.);
  CherRateHEevt[i].SetLineColor(kRed);

  sprintf(myText, "S123 rate HEC-EMEC, %s%s", myTex1, myTex2);
  S123RateHE[i].SetNameTitle(myText, myText);
  S123RateHE[i].GetXaxis()->SetTitle("S123 tracks per bunch/cluster");
  S123RateHE[i].SetBins(400, -100., 3900.);

  sprintf(myText, "S123 rate HEC-EMEC,sum, %s%s", myTex1, myTex2);
  S123RateHEevt[i].SetNameTitle(myText, myText);
  S123RateHEevt[i].GetXaxis()->SetTitle("S123 tracks per event");
  S123RateHEevt[i].SetBins(400, -100., 3900.);
  S123RateHEevt[i].SetLineColor(kRed); 
}

//---------------------------------------------------------------
miniEventHeader     MH;
miniBurstHeader     BH;
eventEssentials     evtEssentials;
cherenkovEssentials cherEssentials;

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

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

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

double time[515];
//-------------------------------------------------------------------------
// histograms and covariances for S123<->GADC vs rate prehistory
double historyTime[10]={0.01, 0.03, 0.1, 0.3, 1.0, 3.0, 10., 30., 100., 0.}; // prehistory integration, msec
double historyRate[10]={10*0.0};
int    nLaps[9]={9*0};
TH1F prehistorySpectrum[10];
for (int i=0; i<9; i++) {
  sprintf(myText, "Average rate over %.2fms", historyTime[i]);
  prehistorySpectrum[i].SetNameTitle(myText, myText);
}
for (int i=0; i<10; i++) {
  prehistorySpectrum[i].SetBins(30000, 0., 300.);
  prehistorySpectrum[i].SetStats(111110);
  prehistorySpectrum[i].GetYaxis()->SetTitle("Bunches");
  prehistorySpectrum[i].GetXaxis()->SetTitle("Rate, 10^{6} p/sec");
}
sprintf(myText, "Instant rate over 0ms");
prehistorySpectrum[9].SetNameTitle(myText, myText);

double RateBins[1000]; // rate bin upper edge
double RateBinsSum0[1000];
double RateBinsSum1[1000];
int    nRateBins[10];
TH2F         cherVStracks2D[1000];
myCovariance cherVStracksCov[1000]; // covariance structure 
for (int iRate=0; iRate<10; iRate++) {
  for (int ib=0; ib<100; ib++) {
    sprintf(myText, "Rate#%d_Bin#%d", iRate, ib);
    cherVStracks2D[100*iRate+ib].SetNameTitle(myText, myText);
    cherVStracks2D[100*iRate+ib].SetBins(100, 0., 300., 100, -100., 3500.);
    cherVStracks2D[100*iRate+ib].GetXaxis()->SetTitle("S123 tracks");
    cherVStracks2D[100*iRate+ib].GetYaxis()->SetTitle("Cherenkov GADC, ADC counts");
    cherVStracks2D[100*iRate+ib].Reset();
    cherVStracksCov[100*iRate+ib].Clear();
  }
}

//-------------------------------------------------------------------------
int iqq=1;
int iqqOld=1;

iqq=3;
iqqOld=3;

for (int i=0; i<512; i++) {time[i]=12.5*i;}


for (int iG=0; iG<2; iG++) { // set X and Y for module channels //
  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, 900, 900);
TCanvas *moduleCanvas1  = new TCanvas("moduleCanvas1",    "Module Currents",  100, 0, 900, 900);
TCanvas *spillCanvas    = new TCanvas("spillCanvas",      "Spill Structure",  200, 0, 900, 600);
TCanvas *positionCanvas0 = new TCanvas("positionCanvas0", "HEC-EMEC, X",        0, 0, 500, 900);
TCanvas *positionCanvas1 = new TCanvas("positionCanvas1", "HEC-EMEC, Y",      500, 0, 500, 900);
TCanvas *positionCanvas2 = new TCanvas("positionCanvas2", "FCal-EMEC, X",    1600, 0, 500, 900);
TCanvas *positionCanvas3 = new TCanvas("positionCanvas3", "FCal-EMEC, Y",    2100, 0, 500, 900);
TCanvas *rateCanvas0     = new TCanvas("rateCanvas0",     "FCal-EMEC rates", 1800, 0, 1200, 1000);
TCanvas *rateCanvas1     = new TCanvas("rateCanvas1",     "Intensiy Integrals",    0, 0, 1200, 1000);

TCanvas *slopeCanvas00   = new TCanvas("slopeCanvas00", "slopeCanvas00", 1200, 0, 1000, 1200);
TCanvas *slopeCanvas01   = new TCanvas("slopeCanvas01", "slopeCanvas01", 3000, 0, 1000, 1200);


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

int myRun0=850;
int myRun1=1200;

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

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

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

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

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

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

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

myRun0=1050;
myRun1=1500;

myRun0=850;
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==0) {
  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;

for (int i=0; i<nRanges; i++) {
  XdiffHE[i].Reset();
  XdiffFE[i].Reset();
  YdiffHE[i].Reset();
  YdiffFE[i].Reset();

  XdiffHEevt[i].Reset();
  XdiffFEevt[i].Reset();
  YdiffHEevt[i].Reset();
  YdiffFEevt[i].Reset();

  dX_HE[i].Reset();
  dY_HE[i].Reset();

  dX_HEevt[i].Reset();
  dY_HEevt[i].Reset();

  CherRateHE[i].Reset();
  S123RateHE[i].Reset();

  CherRateHEevt[i].Reset();
  S123RateHEevt[i].Reset();
}

for (int iRate=0; iRate<10; iRate++) {
  prehistorySpectrum[iRate].Reset();
}

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",         &s123mips.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("mipCount",        &evtEssentials.mipCount[0]);
eventTree -> SetBranchAddress("mipCount2",       &evtEssentials.mipCount2[0]);
eventTree -> SetBranchAddress("mipCount3",       &evtEssentials.mipCount3[0]);

eventTree -> SetBranchAddress("isCluster0",      &evtEssentials.isCluster0[0]);
eventTree -> SetBranchAddress("isCluster1",      &evtEssentials.isCluster1[0]);
eventTree -> SetBranchAddress("isIsolated",      &evtEssentials.isIsolated[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);
    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);
      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, s123mips.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("HV EMEC is Zero, skipping this run\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 ///
  evtEssentials.Clear();
  for (int i=0; i<50; i++) Channel[i].Clear();
  eventTree->GetEntry(iEnt);

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

  if (iqq<3) { // print event RF bunch structure (filled, empty, cluster)
    printf("run#%d  burst=%4d\t event=%5d(%3d) empty=%d notEmpty=%d spillCore=%d signal=%.2f  \n", 
            myRun, MH.burst, MH.event, MH.burstEvent, 
            MH.isEmpty, MH.notEmpty, MH.filledSegment, MH.SignalFactor);
  }

  if (MH.burst != burstSet) { // first or new burst, load burst entry
    if (burstEnt<burstEntries) {
      printf("Loading burst 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->0 for exit\n");
        printf("        1 - event-by-event\n");
        printf("        2 - non-stop run show\n");
        printf("        3 - for quiet run-by-run\n");
        printf("        4 - for quiet all runs\n --> ");
        cin >> iqq;
        if (iqq==0) break;
      }

      if (iqq<3) {      
        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    nLaps=500; // how many laps to average for display
        for (int i=0; i<cherEssentials.ns30; i++) {
          rate+=cherEssentials.sum30[i];
          if (i>0 && i%nLaps==0) {
            rate/=nLaps;
            if (myLapMax<rate) myLapMax=rate;
            rate=0.;
          }
        }      
      
//      dummy->GetYaxis()->SetRangeUser(0., BH.ChGADC*cherEssentials.lapMax); // maximum protons/lap
        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);
        spillCanvas->Clear();
        spillCanvas->Divide(1,2);
        spillCanvas->cd(1);
        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%nLaps==0) {
            tLap/=nLaps;
            rate/=nLaps;
            myLine->DrawLine(tLap-0.5*nLaps*lapTime, BH.ChGADC*rate,
                             tLap+0.5*nLaps*lapTime, BH.ChGADC*rate);
            tLap=0.;
            rate=0.;
          }
        }      
        spillCanvas->Update();
      }
    } else {
      printf("Burst Entry #%d is out of range: nBursts=%d\n", burstEnt, burstEntries);
      continue;
    }
  }

  if (!MH.filledSegment) 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
  int cluster=0;
  int nClusters=0;

  for (int iCl=0; iCl<5; iCl++) {
    if (evtEssentials.RFclust0[iCl]==-1) break;
    nClusters++;
  }

  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();
  }
// compute 9 prehistory rates
  for (int iRate=0; iRate<9; iRate++) {
    historyRate[iRate]=0.;
    nLaps[iRate]=0;
  }
  bool added=true;
  for (int jLap=0; jLap<myLap && added; jLap++) {
    added=false;
    for (int iRate=8; iRate>=0; iRate--) {
      if (jLap*lapTime>historyTime[iRate]) break; // do not integrate if time is longer than needed
      added=true;
      nLaps[iRate]++;
      historyRate[iRate]+=cherEssentials.sum30[myLap-jLap];
    }
  }



  if (iqq<3) {
    for (int iRate=0; iRate<9; iRate++) {
      printf("  rate bin#%d T=%6.2fms nLaps=%5d Integral=%7.0fp Rate=%5.0fp/lap = %fp/ms \n", 
             iRate, historyTime[iRate], nLaps[iRate],   BH.ChGADC*historyRate[iRate], 
             BH.ChGADC*historyRate[iRate]/nLaps[iRate], BH.ChGADC*historyRate[iRate]/nLaps[iRate]/lapTime);
    }    
  }  
  
// 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++) {

    //if (evtEssentials.RFclust0[iClust]!=evtEssentials.RFclust1[iClust]) continue; // use singles only!

    for (int j=evtEssentials.RFclust0[iClust]; j<=evtEssentials.RFclust1[iClust]; j++) {
      if (MaxTracks<evtEssentials.mipCount2[j]) MaxTracks=evtEssentials.mipCount2[j];
      evtMIPs+=evtEssentials.mipCount2[j];
      evtGADC+=evtEssentials.GADCevent[j];
      prehistorySpectrum[9].Fill(BH.ChGADC*evtEssentials.GADCevent[j]); // fill GADC filled bunches spectrum
      for (int iRate=0; iRate<9; iRate++) {
        prehistorySpectrum[iRate].Fill(BH.ChGADC*historyRate[iRate]/nLaps[iRate]/lapTime/1000.); // /1000. is conversion from p/ms to 10^6p/sec
      }
    }
  }
  evtGADC*=BH.ChGADC;
  //if (evtMIPs>100) iqq=1;

  for (int i=nRanges-1; i>=0; i--) {
    if (qHEC>AminRange[i] && qEMEC>AminRange[i]) {
      XdiffHEevt[i].Fill(xEMEC, xHEC);
      YdiffHEevt[i].Fill(yEMEC, yHEC);

      dX_HEevt[i].Fill(xHEC-xEMEC);
      dY_HEevt[i].Fill(yHEC-yEMEC);

      CherRateHEevt[i].Fill(evtGADC);
      S123RateHEevt[i].Fill(evtMIPs);
      break;
    }
  }
  for (int i=nRanges-1; i>=0; i--) {
    if (qFCal>AminRange[i] && qEMEC>AminRange[i]) {
      XdiffFEevt[i].Fill(xEMEC, xFCal);
      YdiffFEevt[i].Fill(yEMEC, yFCal);
      break;
    }
  }

  
// cluster-by-cluster
  for (int iClust=0; iClust<nClusters; iClust++) {
    int iHEC=0;
    int iEMEC=1;
    int iFCal=2;
    if (iFCal==blankModule) iFCal=3;

    double qHEC =myModule[iHEC] .Qcl[iClust];
    double qEMEC=myModule[iEMEC].Qcl[iClust];
    double qFCal=myModule[iFCal].Qcl[iClust];

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

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

    double bunchMIPs=0.;
    double bunchGADC=0.;
    for (int j=evtEssentials.RFclust0[iClust]; j<=evtEssentials.RFclust1[iClust]; j++) {
      bunchMIPs+=evtEssentials.mipCount3[j];
      bunchGADC+=evtEssentials.GADCevent[j];
    }
    bunchGADC*=BH.ChGADC; // calibration

    for (int i=nRanges-1; i>=0; i--) {
      if (qHEC>AminRange[i] && qEMEC>AminRange[i]) {
        XdiffHE[i].Fill(xEMEC, xHEC);
        YdiffHE[i].Fill(yEMEC, yHEC);

        dX_HE[i].Fill(xHEC-xEMEC);
        dY_HE[i].Fill(yHEC-yEMEC);

        CherRateHE[i].Fill(bunchGADC);
        S123RateHE[i].Fill(bunchMIPs);        
        break;
      }
    }
    for (int i=nRanges-1; i>=0; i--) {
      if (qFCal>AminRange[i] && qEMEC>AminRange[i]) {
        XdiffFE[i].Fill(xEMEC, xFCal);
        YdiffFE[i].Fill(yEMEC, yFCal);
        break;
      }
    }


  } // clusters  

  if (MH.isEmpty || iqq>2) continue; // skip the event display

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

  for (int iClust=0; iClust<nClusters; iClust++) {
    printf("Cluster#%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.mipCount2[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("Sum 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);
    }
  }

  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.;

//    for (int ich=0; ich<5; ich++) {
//      pedWindows[ich+5*iMod].Print();
//    }
//    continue;

    int id1=iChD+5*iMod;
    int nPedWindows=pedWindows[id1].nPedWindows;
    int is0=pedWindows[id1].firstSample[0];
    int is1=pedWindows[id1].lastSample[nPedWindows-1];
//    printf(" module#%d nPedWindows=%d samples range=%d...%d\n", iMod, nPedWindows, is0, is1);
//    pedWindows[id1].Print();

    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->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++) {
      if (myModule[iMod].Qcl[iClust]>0.) {
         sprintf(myText, "Q=%.0ffC X=%.2f Y=%.2f \n", 
                 myModule[iMod].Qcl[iClust], myModule[iMod].Xcl[iClust], 
                                             myModule[iMod].Ycl[iClust]);
      } else {
         sprintf(myText, "Q=%.0ffC \n", myModule[iMod].Qcl[iClust]);
      }
      myLatex->DrawLatex(T0+165.*evtEssentials.RFclust0[iClust], amaxF+0.06*df, myText);
      
    }
    

    iZone++;

  }

// GADC response //
  double GADCmin=0.;
  double GADCmax=0.;
//    for (int ich=0; ich<5; ich++) {
//      pedWindows[ich+5*iMod].Print();
//    }
//    continue;
  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);
  for (int i=0; i<40 && 165.*(i+1)<6400.; i++) {
    myLine->DrawLine(T0+165.*i, evtEssentials.GADCevent[i], 
                     T0+165.*(i+1), evtEssentials.GADCevent[i]);
    if (i>0) {
      myLine->DrawLine(T0+165.*i, evtEssentials.GADCevent[i-1], 
                       T0+165.*i, evtEssentials.GADCevent[i]);
    }
  }
  myLine->SetLineColor(kBlue);
  myLine->SetLineWidth(2);
  for (int i=0; i<40 && 165.*(i+1)<6400.; i++) {
    myLine->DrawLine(T0+165.*i,     evtEssentials.GADCevent1[i], 
                     T0+165.*(i+1), evtEssentials.GADCevent1[i]);
    if (i>0) {
      myLine->DrawLine(T0+165.*i, evtEssentials.GADCevent1[i-1], 
                       T0+165.*i, evtEssentials.GADCevent1[i]);
    }
  }
  myLine->SetLineColor(kRed);
  for (int i=0; i<40 && 165.*(i+1)<6400.; i++) {
    myLine->DrawLine(T0+165.*i, evtEssentials.MPICevent[i], 
                     T0+165.*(i+1), evtEssentials.MPICevent[i]);
    if (i>0) {
      myLine->DrawLine(T0+165.*i, evtEssentials.MPICevent[i-1], 
                       T0+165.*i, evtEssentials.MPICevent[i]);
    }
  }
  
  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 (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 tracs");
  dummy->DrawCopy();
  myLine->SetLineColor(kBlack);
  for (int i=0; i<40 && 165.*(i+1)<6400.; i++) {
    myLine->DrawLine(T0+165.*i, evtEssentials.mipCount[i], 
                     T0+165.*(i+1), evtEssentials.mipCount[i]);
    if (i>0) {
      myLine->DrawLine(T0+165.*i, evtEssentials.mipCount[i-1], 
                       T0+165.*i, evtEssentials.mipCount[i]);
    }
  }
  myLine->SetLineColor(kRed);
  for (int i=0; i<40 && 165.*(i+1)<6400.; i++) {
    myLine->DrawLine(T0+165.*i-10.,     evtEssentials.mipCount2[i], 
                     T0+165.*(i+1)-10., evtEssentials.mipCount2[i]);
    if (i>0) {
      myLine->DrawLine(T0+165.*i-10., evtEssentials.mipCount2[i-1], 
                       T0+165.*i-10., evtEssentials.mipCount2[i]);
    }
  }
  myLine->SetLineColor(kBlue);
  for (int i=0; i<40 && 165.*(i+1)<6400.; i++) {
    myLine->DrawLine(T0+165.*i+10.,     evtEssentials.mipCount3[i], 
                     T0+165.*(i+1)+10., evtEssentials.mipCount3[i]);
    if (i>0) {
      myLine->DrawLine(T0+165.*i+10., evtEssentials.mipCount3[i-1], 
                       T0+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];
//    for (int ir=ir0; ir<ir1+1; ir++) sumS123+=evtEssentials.mipCount3[ir];
//    for (int ir=ir0; ir<ir1+1; ir++) sumS123+=evtEssentials.mipCount[ir];
    sprintf(myText, "%d", sumS123);
    myLatex->DrawLatex(T0+82.5*(ir0+ir1), 1.07*mipCountMAX, myText);
  }



  myLine->SetLineColor(kBlack);
  
  moduleCanvas ->Update();
  moduleCanvas1->Update();
  if (iqq < 2) {
    printf("Next Event!\n");
    printf("  push->0 for exit\n");
    printf("        1 - event-by-event\n");
    printf("        2 - non-stop run show\n");
    printf("        3 - for quiet run-by-run\n");
    printf("        4 - for quiet all runs\n --> ");
    cin >> iqq;
    if (iqq==0) break;
  }
}
///if (moduleHV[1]>0.) {
positionCanvas0->Clear();
positionCanvas0->Divide(3, 5);
positionCanvas1->Clear();
positionCanvas1->Divide(3, 5);
positionCanvas2->Clear();
positionCanvas2->Divide(2, 5);
positionCanvas3->Clear();
positionCanvas3->Divide(2, 5);
sprintf(myText, "Run#%d HEC-EMEC, X", myRun);
positionCanvas0->SetTitle(myText);
sprintf(myText, "Run#%d HEC-EMEC, Y", myRun);
positionCanvas1->SetTitle(myText);
sprintf(myText, "Run#%d FCal-EMEC, X", myRun);
positionCanvas2->SetTitle(myText);
sprintf(myText, "Run#%d FCal-EMEC, Y", myRun);
positionCanvas3->SetTitle(myText);

sprintf(myText, "Run#%d HEC-EMEC, rates", myRun);
rateCanvas0->Clear();
rateCanvas0->Divide(4,5);
rateCanvas0->SetTitle(myText);


for (int i=0; i<nRanges; i++) {
  positionCanvas0->cd(3*i+1)->SetGrid();
  XdiffHE[i].   DrawCopy("box");
  XdiffHEevt[i].DrawCopy("same");
  positionCanvas1->cd(3*i+1)->SetGrid();
  YdiffHE[i]   .DrawCopy("box");
  YdiffHEevt[i].DrawCopy("same");
  positionCanvas2->cd(2*i+1)->SetGrid();
  XdiffFE[i].   DrawCopy("box");
  XdiffFEevt[i].DrawCopy("same");
  positionCanvas3->cd(2*i+1)->SetGrid();
  YdiffFE[i].   DrawCopy("box");
  YdiffFEevt[i].DrawCopy("same");

  positionCanvas0->cd(3*i+2)->SetGrid();
  dX_HE[i]   .DrawCopy();  
  positionCanvas0->cd(3*i+3)->SetGrid();
  dX_HEevt[i].DrawCopy();
 
  positionCanvas1->cd(3*i+2)->SetGrid();
  dY_HE[i]   .DrawCopy();  
  positionCanvas1->cd(3*i+3)->SetGrid();
  dY_HEevt[i].DrawCopy(); 

  rateCanvas0->cd(4*i+1)->SetGrid();
  S123RateHE[i]   .DrawCopy();  
  rateCanvas0->cd(4*i+2)->SetGrid();
  S123RateHEevt[i].DrawCopy(); 
  rateCanvas0->cd(4*i+3)->SetGrid();
  CherRateHE[i].   DrawCopy();  
  rateCanvas0->cd(4*i+4)->SetGrid();
  CherRateHEevt[i].DrawCopy();  
}

rateCanvas1->Clear();
rateCanvas1->Divide(3,3);
for (int iRate=0; iRate<9; iRate++) {
  rateCanvas1->cd(iRate+1)->SetGrid(1);
  prehistorySpectrum[iRate].DrawCopy();
}
printf("maximum S123 mipCount2=%d\n", MaxTracks);
/// create rate/prehistory bins ////

for (int iRate=0; iRate<10; iRate++) {
  nRateBins[iRate]=0;
  for (int ib=0; ib<100; ib++) {
    RateBins[100*iRate+ib]=1.0e10; // set upper bin to infinity
  }

  printf(" %s\n", prehistorySpectrum[iRate].GetTitle());
  double sBin=0.;
  int ibin0=0;
  int nBins=0; 
  for (int i=1; i<prehistorySpectrum[iRate].GetNbinsX(); i++) {
    sBin+=prehistorySpectrum[iRate].GetBinContent(i);
    if (sBin>1000.) {
      printf("  nBins=%3d    bin0#%4d=(%7.3f) bin1#%4d(%7.3f) Sum=%.1f\n", nBins, ibin0, prehistorySpectrum[iRate].GetBinCenter(ibin0), 
                                                                  i,     prehistorySpectrum[iRate].GetBinCenter(i), sBin);
      RateBins[100*iRate+nBins]=prehistorySpectrum[iRate].GetBinCenter(i) + 
                            0.5*prehistorySpectrum[iRate].GetBinWidth(i); // upper edge
      ibin0=i+1;
      sBin=0.;
      nBins++;
      if (nBins>99) nBins=99;
    }
  }
  //RateBins[100*iRate+nBins]=prehistorySpectrum[iRate].GetXaxis()->GetXmax();
  nRateBins[iRate]=nBins+1;
  printf("======================== rate#%d will use %d bins\n\n", iRate, nBins+1);  
}

for (int iRate=0; iRate<10; iRate++) {
  for (int ib=0; ib<nRateBins[iRate]+1; ib++) {
    int id=100*iRate+ib;
    RateBinsSum0[id]=0.;
    RateBinsSum1[id]=0.;
    cherVStracks2D[id].Reset();
    int nxb=MaxTracks+1;
    if (nxb>300) {
      nxb=300;
    }
    cherVStracks2D[id].SetBins(nxb, -0.5, MaxTracks+0.5, 1000, -100., 3500.);
    cherVStracksCov[id].Clear();
    printf("arrays and hists for bin#%3d of rate%2d is ready, id=%d\n", ib, iRate, id); 
  }
}

printf("Updating plots...\n");

positionCanvas0->Update();
positionCanvas1->Update();
positionCanvas2->Update();
positionCanvas3->Update();

rateCanvas0->Update();
rateCanvas1->Update();

printf("Saving plots...\n");

sprintf(myText, "reports/BeamPosition/ModuleCorrelationRun%04devtSum0429pass1.ps(", myRun);
positionCanvas0->Print(myText);
sprintf(myText, "reports/BeamPosition/ModuleCorrelationRun%04devtSum0429pass1.ps", myRun);
positionCanvas1->Print(myText);
positionCanvas2->Print(myText);
positionCanvas3->Print(myText);
rateCanvas0->Print(myText);
sprintf(myText, "reports/BeamPosition/ModuleCorrelationRun%04devtSum0429pass1.ps)", myRun);
rateCanvas1->Print(myText);

printf("\n\n\n Pass#2 \n\n\n");
//iqq=1;
//for (int iRate=0; iRate<9; iRate++) {
//  printf("Integration time[%d]=%fms\n", iRate, historyTime[iRate]);
//}

burstEnt=0;
for (int iEnt=0; iEnt<eventEntries; iEnt++) {  //  second event loop ///
  evtEssentials.Clear();
  for (int i=0; i<50; i++) Channel[i].Clear();
//  printf("trying to load event entry#%d\n", iEnt); 
  eventTree->GetEntry(iEnt);
//  printf("   event entry#%d loaded!\n", iEnt); 

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

  if (MH.burst != burstSet) { // first or new burst, load burst entry
    if (burstEnt<burstEntries) {
      printf("Loading burst 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->0 for exit\n");
        printf("        1 - event-by-event\n");
        printf("        2 - non-stop run show\n");
        printf("        3 - for quiet run-by-run\n");
        printf("        4 - for quiet all runs\n --> ");
        cin >> iqq;
        if (iqq==0) break;
      }

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

  if (!MH.filledSegment) 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];
  }
  if (HGoverflow && iqq<3) {iqqOld=iqq; iqq=1;} 
//  else if (iqq!=3)  iqq=iqqOld;
//

// compute centers of gravity and some sums
  int cluster=0;
  int nClusters=0;

  for (int iCl=0; iCl<5; iCl++) {
    if (evtEssentials.RFclust0[iCl]==-1) break;
    nClusters++;
  }

  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
  

// compute 9 prehistory rates
  for (int iRate=0; iRate<9; iRate++) {
    historyRate[iRate]=0.;
    nLaps[iRate]=0;
  }
  bool added=true; // prepare history rates //
  int myLap=MH.firstGated/30;
  for (int jLap=0; jLap<myLap && added; jLap++) {
    added=false;
    for (int iRate=8; iRate>=0; iRate--) {
      if (jLap*lapTime>historyTime[iRate]) break; // do not integrate if time is longer than needed
      added=true;
      nLaps[iRate]++;
      historyRate[iRate]+=cherEssentials.sum30[myLap-jLap];
    }
  }

  for (int iRate=0; iRate<9; iRate++) {
    historyRate[iRate]*=BH.ChGADC/nLaps[iRate]/lapTime/1000.; // conversion from GADC to 10^6p/sec;
  }

  for (int iRate=0; iRate<9 && iqq==1; iRate++) {
    printf("Integration time[%d]=%fms\n", iRate, historyTime[iRate]);
  }

  for (int iClust=0; iClust<nClusters; iClust++) {

    for (int j=evtEssentials.RFclust0[iClust]; j<=evtEssentials.RFclust1[iClust]; j++) {
      
      historyRate[9]=BH.ChGADC*evtEssentials.GADCevent1[j]; // bunch rate in GADC tracks
      
      for (int iRate=0; iRate<10; iRate++) {
        for (int ib=0; ib<nRateBins[iRate]; ib++) {
          int id=100*iRate+ib;
          if (historyRate[iRate]<RateBins[id]){
            cherVStracks2D[id] .Fill((double)evtEssentials.mipCount2[j], evtEssentials.GADCevent1[j]);
            cherVStracksCov[id].Fill((double)evtEssentials.mipCount2[j], evtEssentials.GADCevent1[j]);
            RateBinsSum0[id]+=1.0;
            RateBinsSum1[id]+=historyRate[iRate];
            break;
          }
        }
      }  // rates loop
    }    // samples-in-cluster loop
    
/*
    double sumMIPS=0.;
    double sumGADC=0.;
    double clustMax=0.;
    for (int j=evtEssentials.RFclust0[iClust]; j<=evtEssentials.RFclust1[iClust]; j++) {
      sumMIPS+=(double)evtEssentials.mipCount2[j];
      sumGADC+=evtEssentials.GADCevent1[j];
      if (clustMax<evtEssentials.GADCevent1[j])  clustMax=evtEssentials.GADCevent1[j];
    }  
    historyRate[9]=BH.ChGADC*clustMax; // bunch rate in GADC tracks
      
    for (int iRate=0; iRate<10; iRate++) {
      for (int ib=0; ib<nRateBins[iRate]; ib++) {
        int id=100*iRate+ib;
        if (historyRate[iRate]<RateBins[id]){
          cherVStracks2D[id] .Fill(sumMIPS, sumGADC);
          cherVStracksCov[id].Fill(sumMIPS, sumGADC);
          RateBinsSum0[id]+=1.0;
          RateBinsSum1[id]+=historyRate[iRate];
          break;
        }
      }
    } // rates loop
*/
  } // clusters loop
  

 

} // end of second loop over events
//

slopeCanvas01->Clear();
slopeCanvas01->Divide(4,5);
int iZone=0;
int iqq1=1;
double maxRate=0.;
bool firstPlot=true;

for (int iRate=0; iRate<9; iRate++) {
//  printf("Integration time[%d]=%fms\n", iRate, historyTime[iRate]);
  for (int ib=0; ib<nRateBins[iRate]; ib++) {
    int id=100*iRate+ib;
    if (RateBinsSum0[id]<10.) continue;
    RateBinsSum1[id]/=RateBinsSum0[id];
    if (maxRate<RateBinsSum1[id]) maxRate=RateBinsSum1[id];
  }
  printf("Integration time[%d]=%fms \n", iRate, historyTime[iRate]);
}
printf(" maximum rate=%.2fx10^{6}p/sec\n", maxRate);
dummyRate->SetBins(200, 0., 1.05*maxRate);

for (int iRate=0; iRate<9; iRate++) {
  slopeCanvas00->Clear();
  slopeCanvas00->Divide(3,3);
  slopeCanvas01->cd(2*iRate+1)->SetGrid();
  sprintf(myText, "Integration time=%.2fms \n", historyTime[iRate]);
  dummyRate->SetTitle(myText);
  dummyRate->GetYaxis()->SetTitle("Correlation");
  dummyRate->GetYaxis()->SetRangeUser(-1., 1.);
  dummyRate->SetStats(0);
  dummyRate->DrawCopy();
  slopeCanvas01->cd(2*iRate+2)->SetGrid();
  dummyRate->GetYaxis()->SetTitle("GADC/track");
  dummyRate->GetYaxis()->SetRangeUser(0., 3./BH.ChGADC);
  dummyRate->DrawCopy();
  myLine->SetLineWidth(4);
  myLine->SetLineColor(kGray);
  myLine->DrawLine(0., 1./BH.ChGADC, 1.05*maxRate, 1./BH.ChGADC);
  myLatex->SetTextSize(0.06);
  myLatex->SetTextAlign(12);
  myLatex->SetTextColor(kRed);
  myLatex->DrawLatex(0.05*maxRate, 2.75/BH.ChGADC, "Red:  (0,0) fit");
  myLatex->SetTextColor(kRed);
  myLatex->DrawLatex(0.05*maxRate, 2.50/BH.ChGADC, "Black: free fit");

  iZone=0;
  for (int ib=0; ib<nRateBins[iRate]; ib++) {
    int id=100*iRate+ib;
    if (RateBinsSum0[id]<10.) continue;
    sprintf(myText, "Rate(%.2fms)=%.2fx10^{6}p/sec", historyTime[iRate], RateBinsSum1[id]);
    cherVStracks2D[id].SetTitle(myText);
    iZone++;
    slopeCanvas00->cd(iZone)->SetGrid();
    gStyle->SetOptStat(111110);
    int lastX=cherVStracks2D[id].GetNbinsX();
    int lastY=cherVStracks2D[id].GetNbinsX();
    double sY=0.;
    double sX=0.;
    for (int ix=1; ix<cherVStracks2D[id].GetNbinsX(); ix++) {
      for (int iy=1; iy<cherVStracks2D[id].GetNbinsY(); iy++) {
        sY+=cherVStracks2D[id].GetBinContent(ix, iy);
      }
      if (sY>0.99*cherVStracks2D[id].GetEntries()) {
        lastX=ix;
        break;
      }
    }
    for (int iy=1; iy<cherVStracks2D[id].GetNbinsY(); iy++) {
      for (int ix=1; ix<cherVStracks2D[id].GetNbinsX(); ix++) {
        sX+=cherVStracks2D[id].GetBinContent(ix, iy);
      }
      if (sX>0.99*cherVStracks2D[id].GetEntries()) {
        lastY=iy;
        break;
      }
    }

    double x0=cherVStracks2D[id].GetXaxis()->GetXmin();
    double x1=cherVStracks2D[id].GetXaxis()->GetBinCenter(lastX);
    double y0=cherVStracks2D[id].GetYaxis()->GetXmin();
    double y1=cherVStracks2D[id].GetYaxis()->GetBinCenter(lastY);
    y0=-10.;
    cherVStracks2D[id].GetXaxis()->SetRangeUser(x0, x1);
    cherVStracks2D[id].GetYaxis()->SetRangeUser(y0, y1);
    cherVStracks2D[id].DrawCopy("box");

    myLatex->SetTextSize(0.06);
    sprintf(myText, "corr=%.2f", cherVStracksCov[id].GetCorrelation(10));
    myLatex->SetTextAlign(12);
    myLatex->DrawLatex(0.95*x0+0.05*x1, 0.95*y1+0.05*y0, myText); 
    sprintf(myText, " slopeX=%.2f",    cherVStracksCov[id].GetSlopeX(10));
    myLatex->DrawLatex(0.95*x0+0.05*x1, 0.85*y1+0.15*y0, myText); 
    sprintf(myText, "slopeX0=%.2f",    cherVStracksCov[id].GetSlopeX0(10));
    myLatex->DrawLatex(0.95*x0+0.05*x1, 0.75*y1+0.25*y0, myText); 

    myLatex->SetTextSize(0.05);
    myLatex->SetTextAlign(22);
    myLatex->SetTextColor(kBlack);
    slopeCanvas01->cd(2*iRate+1);
    myLatex->DrawLatex(RateBinsSum1[id],     cherVStracksCov[id].GetCorrelation(10), "x");
    slopeCanvas01->cd(2*iRate+2);
    if (cherVStracksCov[id].GetSlopeX(10)>3./BH.ChGADC) {
      myLatex->DrawLatex(RateBinsSum1[id],    2.9/BH.ChGADC ,      "#uparrow");
    } else if (cherVStracksCov[id].GetSlopeX(10)<0.) {
      myLatex->DrawLatex(RateBinsSum1[id],    0.1/BH.ChGADC ,      "#downarrow");
    } else {
      myLatex->DrawLatex(RateBinsSum1[id],     cherVStracksCov[id].GetSlopeX(10),      "x");
    }
    myLatex->SetTextColor(kRed);
    if (cherVStracksCov[id].GetSlopeX0(10)>3./BH.ChGADC) {
      myLatex->DrawLatex(RateBinsSum1[id],    2.9/BH.ChGADC ,      "#uparrow");
    } else if (cherVStracksCov[id].GetSlopeX0(10)<0.) {
      myLatex->DrawLatex(RateBinsSum1[id],    0.2/BH.ChGADC ,      "#downarrow");
    } else {
      myLatex->DrawLatex(RateBinsSum1[id],     cherVStracksCov[id].GetSlopeX0(10),     "x");
    }

    
    myLatex->SetTextSize(0.06);
    slopeCanvas00->Update();
    slopeCanvas01->Update();
//    printf("push for next set ->");
//    cin>>iqq1;
    if (iZone==9 || ib==nRateBins[iRate]-1) {
      iZone=0;
//      printf("push for next set ->");
//      cin>>iqq1;
//      if (iqq1==-1) exit(-1);
//      if (iqq1==0)  break;
      if (firstPlot) {
        firstPlot=false;
        sprintf(myText, "reports/BeamPosition/ModuleCorrelationRun%04devtSum0429pass2.ps(", myRun);
        slopeCanvas00->Print(myText);
      } else {
        sprintf(myText, "reports/BeamPosition/ModuleCorrelationRun%04devtSum0429pass2.ps", myRun);
        slopeCanvas00->Print(myText);
      }
      slopeCanvas00->Clear();
      slopeCanvas00->Divide(3,3);
    }
  }
  if (iqq1==0) break;
}
sprintf(myText, "reports/BeamPosition/ModuleCorrelationRun%04devtSum0429pass2.ps)", myRun);
slopeCanvas01->Print(myText);

//
if (iqq==0) break; // exit run loop
if (iqq < 4) {
  printf("new run!\n");
  printf("  push->0 for exit\n");
  printf("        1 - event-by-event\n");
  printf("        2 - non-stop run show\n");
  printf("        3 - for quiet run-by-run\n");
  printf("        4 - for quiet all runs\n --> ");
  cin >> iqq;
  if (iqq==0) break;  // exit run loop
}
f->Close();


} // end of one run loop

}


/*


        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 {
              }
            }
          }
        }
*/