#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
//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 (mipCount[i]<40) {
  return mipCount[i];
} else if (mipCount[i]<60) {
  return int(0.5+((60-i)*mipCount[i]+(i-40)*mipCount2[i])/20.);
}
return mipCount2[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 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 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 HiLumRead2014_05_12() {
///=============================================================================

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[200];
char   myTex1[200];
char   myTex2[200];


//---------------------------------------------------------------
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("Prehistory 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];
double rateCut[10]={10*0.0};
double slopePass0[10]={10*0.0};
double slopePass1[10]={10*0.0};

TH2F         cherVStracks2D[1000];
myCovariance cherVStracksCov[1000]; // covariance structure 
for (int iRate=0; iRate<10; iRate++) {
  for (int ib=0; ib<100; ib++) {
    sprintf(myText, "Prehistory_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();
  }
}
//--------- same, but for Min(HEC, EMEC) for position resolution study
TH1F *minHECEMEC = new TH1F("HEC_EMEC minumum", "HEC_EMEC minumum, fC", 10000, -200., 39800.);
TH1F *HECspec    = new TH1F("HEC spectrum",     "HEC spectrum, fC",     10000, -200., 39800.);
TH1F *EMECspec   = new TH1F("EMEC spectrum",    "EMEC spectrum, fC",    10000, -200., 39800.);
minHECEMEC->GetXaxis()->SetTitle("min(HEC, EMEC) event integral, fC");
minHECEMEC->GetYaxis()->SetTitle("Events");
HECspec ->GetXaxis()->SetTitle("HEC event integral, fC");
HECspec ->GetYaxis()->SetTitle("Events");
EMECspec->GetXaxis()->SetTitle("EMEC event integral, fC");
EMECspec->GetYaxis()->SetTitle("Events");
TH2F HECvsEMECx[100];
TH2F HECvsEMECy[100];
TH1F HECdEMECx[100];
TH1F HECdEMECy[100];
myCovariance hecVSemecX[100];
myCovariance hecVSemecY[100];
double HECvsEMECbins[100]; // upper boundary
double HECvsEMECsum0[100];
double HECvsEMECsum1[100]; // bin 'centers'
int nHECEMECbins=0;
double HECbins[100]; // upper boundary
double HECsum0[100];
double HECsum1[100]; // bin 'centers'
int nHECbins=0;
double EMECbins[100]; // upper boundary
double EMECsum0[100];
double EMECsum1[100]; // bin 'centers'
int nEMECbins=0;

for (int i=0; i<100; i++) {
  sprintf(myText, "HEC-EMEC X position correlation bin#%d", i);
  HECvsEMECx[i].SetNameTitle(myText, myText);
  HECvsEMECx[i].SetBins(120, -1.2, 1.2, 120, -1.2, 1.2);
  HECvsEMECx[i].GetXaxis()->SetTitle("EMEC position, X");
  HECvsEMECx[i].GetYaxis()->SetTitle("HEC position, X");
  sprintf(myText, "HEC-EMEC Y position correlation bin#%d", i);
  HECvsEMECy[i].SetNameTitle(myText, myText);
  HECvsEMECy[i].SetBins(120, -1.2, 1.2, 120, -1.2, 1.2);
  HECvsEMECy[i].GetXaxis()->SetTitle("EMEC position, Y");
  HECvsEMECy[i].GetYaxis()->SetTitle("HEC position, Y");
  sprintf(myText, "HEC-EMEC dX bin#%d", i);
  HECdEMECx[i].SetNameTitle(myText, myText);
  HECdEMECx[i].SetBins(101, -2.02, 2.02);
  HECdEMECx[i].GetXaxis()->SetTitle("HEC-EMEC dX, a.u.");
  sprintf(myText, "HEC-EMEC dY bin#%d", i);
  HECdEMECy[i].SetNameTitle(myText, myText);
  HECdEMECy[i].SetBins(101, -2.02, 2.02);
  HECdEMECy[i].GetXaxis()->SetTitle("HEC-EMEC dY, a.u.");
  
}

//-------------------------------------------------------------------------
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, 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*)");
TCanvas *positionCanvas0 = new TCanvas("positionCanvas0", "HEC-EMEC, X",               0, 0, canvasWidth, canvasHeight);
TCanvas *positionCanvas1 = new TCanvas("positionCanvas1", "HEC-EMEC, Y",     canvasWidth, 0, canvasWidth, canvasHeight);
TCanvas *positionCanvas2 = new TCanvas("positionCanvas2", "FCal-EMEC, X",    canvasWidth, 0, canvasWidth, canvasHeight);
TCanvas *positionCanvas3 = new TCanvas("positionCanvas3", "FCal-EMEC, Y",    canvasWidth, 0, canvasWidth, canvasHeight);
positionCanvas0->Connect("ProcessedEvent(Int_t,Int_t,Int_t,TObject*)", 0, 0, "MyExecuteEvent(Int_t,Int_t,Int_t,TObject*)");
positionCanvas1->Connect("ProcessedEvent(Int_t,Int_t,Int_t,TObject*)", 0, 0, "MyExecuteEvent(Int_t,Int_t,Int_t,TObject*)");
positionCanvas2->Connect("ProcessedEvent(Int_t,Int_t,Int_t,TObject*)", 0, 0, "MyExecuteEvent(Int_t,Int_t,Int_t,TObject*)");
positionCanvas3->Connect("ProcessedEvent(Int_t,Int_t,Int_t,TObject*)", 0, 0, "MyExecuteEvent(Int_t,Int_t,Int_t,TObject*)");

TCanvas *rateCanvas0     = new TCanvas("rateCanvas0",     "FCal-EMEC rates",    canvasWidth, 0, canvasWidth,  canvasHeight);
TCanvas *rateCanvas1     = new TCanvas("rateCanvas1",     "Intensiy Integrals",           0, 0, canvasWidth,  canvasHeight);
rateCanvas0->Connect("ProcessedEvent(Int_t,Int_t,Int_t,TObject*)", 0, 0, "MyExecuteEvent(Int_t,Int_t,Int_t,TObject*)");
rateCanvas1->Connect("ProcessedEvent(Int_t,Int_t,Int_t,TObject*)", 0, 0, "MyExecuteEvent(Int_t,Int_t,Int_t,TObject*)");


TCanvas *slopeCanvas00   = new TCanvas("slopeCanvas00", "slopeCanvas00", (int)2.0*canvasWidth, 0, (int)1.5*canvasWidth, canvasHeight);
TCanvas *slopeCanvas01   = new TCanvas("slopeCanvas01", "slopeCanvas01", (int)3.5*canvasWidth, 0, (int)1.5*canvasWidth, canvasHeight);
slopeCanvas00->Connect("ProcessedEvent(Int_t,Int_t,Int_t,TObject*)", 0, 0, "MyExecuteEvent(Int_t,Int_t,Int_t,TObject*)");
slopeCanvas01->Connect("ProcessedEvent(Int_t,Int_t,Int_t,TObject*)", 0, 0, "MyExecuteEvent(Int_t,Int_t,Int_t,TObject*)");

TCanvas *slopeCanvasXX = new TCanvas("slopeCanvasXX", "slopeCanvasXX", 2*canvasWidth, 0, 900, 900); 

TCanvas *GADCspectra00   = new TCanvas("GADCspectra00", "GADCspectra00", 1600, 0, canvasWidth, canvasHeight);
GADCspectra00->Connect("ProcessedEvent(Int_t,Int_t,Int_t,TObject*)", 0, 0, "MyExecuteEvent(Int_t,Int_t,Int_t,TObject*)");

TCanvas *globalSummary = new TCanvas("globalSummary", "globalSummary", 1600, 0, canvasWidth, canvasHeight);
globalSummary->Connect("ProcessedEvent(Int_t,Int_t,Int_t,TObject*)", 0, 0, "MyExecuteEvent(Int_t,Int_t,Int_t,TObject*)");

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

slopeCanvas00->Iconify();
slopeCanvas01->Iconify();
GADCspectra00->Iconify();

positionCanvas2->Iconify();
positionCanvas3->Iconify();

rateCanvas0->Iconify();


TH1F trackSpectrum[36]; // 0...10 tracks spectra for GADC, 1..3 samples summed
for (int nSum=0; nSum<3; nSum++) {
  for (int nTrack=0; nTrack<12; nTrack++) {
    int id=nTrack*3+nSum;
    if (nTrack==11) {
      sprintf(myText, "GADC noise (empty events), %d samples", nSum+1);
      trackSpectrum[id].SetLineColor(kRed);
      trackSpectrum[id].SetFillColor(kRed);
    } else {
      sprintf(myText, "GADC %d tracks, %d samples", nTrack, nSum+1);
      trackSpectrum[id].SetFillColor(kGray);
    }
    trackSpectrum[id].SetNameTitle(myText, myText);
    trackSpectrum[id].SetBins(5000, -999.5, 3999.5);
  }
}
//=============================================================================

double slopeVSrun[2000];
double xVSrun[2000];
double yVSrun[2000];
int    myRunX[2000];
TH1F *xProfile = new TH1F("xProfile", "xProfile", 250, -2.5, 2.5);
TH1F *yProfile = new TH1F("yProfile", "yProfile", 250, -2.5, 2.5);
//=============================================================================
int MaxTracks=0; // maximum tracks seen by S123 using count2 method;
int runIndex=0;

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=1098;
//myRun1=1099;

//myRun0=850;
//myRun1=1500;

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

for (int id=0; id<36; id++) {
  trackSpectrum[id].Reset();
}

xProfile->Reset();
yProfile->Reset();
minHECEMEC->Reset();
HECspec   ->Reset();
EMECspec  ->Reset();
nHECEMECbins=0;
nHECbins=0;
nEMECbins=0;

for (int i=0; i<100; i++) {
  HECdEMECx[i].Reset();
  HECdEMECy[i].Reset();
  HECvsEMECx[i].Reset();
  HECvsEMECy[i].Reset();
  HECvsEMECbins[i]=0.;
  HECvsEMECsum0[i]=0.;
  HECvsEMECsum1[i]=0.;
  hecVSemecX[i].Clear();
  hecVSemecY[i].Clear();
}


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("\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 ///
  evtEssentials.Clear();
  for (int i=0; i<50; i++) Channel[i].Clear();
  eventTree->GetEntry(iEnt);

  if (iEnt==eventToCatch) iqq=1;

  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);
    spillCanvas  ->Show();
    spillCanvas1 ->Show();
    moduleCanvas ->Show();
    moduleCanvas1->Show();
    
    positionCanvas0->Iconify();
    positionCanvas1->Iconify();
    rateCanvas0    ->Iconify();
    rateCanvas1    ->Iconify();

    GADCspectra00->Iconify();    
  }

  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->0 for exit\n");
        printf("        1 - event-by-event\n");
        printf("        2 - spill-by-spill\n");
        printf("        3 - for quiet run-by-run\n");
        printf("        4 - for quiet all runs\n --> ");
        cin >> iqq;
        if (iqq==0) break;
        if (iqq<0) {
          eventToCatch=-iqq;
          iqq=3;
        }
        if (iqq<3) {
          moduleCanvas ->Show();
          moduleCanvas1->Show();
          spillCanvas  ->Show();
          spillCanvas1 ->Show();
          positionCanvas0 ->Iconify();
          positionCanvas1 ->Iconify();
          slopeCanvas00    ->Iconify();
          slopeCanvas01    ->Iconify();
          slopeCanvasXX    ->Iconify();
        } else {
          moduleCanvas ->Iconify();
          moduleCanvas1->Iconify();
          spillCanvas  ->Iconify();
          spillCanvas1 ->Iconify();
          positionCanvas0 ->Show();
          positionCanvas1 ->Show();
          slopeCanvas00    ->Show();
          slopeCanvas01    ->Show();
          slopeCanvasXX    ->Show();
        }

        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
        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%nLaps==0) {
            tLap/=nLaps;
            rate/=nLaps;
            myLine->DrawLine(tLap-nLaps*lapTime, BH.ChGADC*rate,
                             tLap+nLaps*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;
    }
  }

  if (!MH.filledSegment) {
    if (MH.isEmpty) { // fill nTrack=0 GADC spectra as noise
      for (int i=0; i<36; i+=3) {
        double s=0.;
        for (int nSum=0; nSum<3; nSum++) {
          s+=evtEssentials.GADCevent[i+nSum];
          trackSpectrum[33+nSum].Fill(s);
        }        
      }
    }
    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++) {
    for (int j=evtEssentials.RFclust0[iClust]; j<=evtEssentials.RFclust1[iClust]; j++) {
      if (MaxTracks<evtEssentials.mipCountX(j)) MaxTracks=evtEssentials.mipCountX(j);
      evtMIPs+=evtEssentials.mipCountX(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
      }
      if (evtEssentials.isIsolated[j] && evtEssentials.mipCount[j]<11) {
        double s=0.;
        for (int nSum=0; nSum<3; nSum++) {
          s+=evtEssentials.GADCevent[j+nSum];
          trackSpectrum[evtEssentials.mipCount[j]*3+nSum].Fill(s);
        }              
      }
    }
  }
  evtGADC*=BH.ChGADC;

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

    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");
//    myLatex->SetTextColor(kRed);
//    myLatex->DrawLatex(myLap*lapTime, 0.5*(yHEC+yEMEC), "o");
//    myLatex->DrawLatex(myLap*lapTime, yHEC-yEMEC, "-");

    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();
  }

  //if (evtMIPs>100) iqq=1;
  double qMIN=qHEC;
  if (qEMEC<qMIN) qMIN=qEMEC;
  minHECEMEC->Fill(qMIN);  
  HECspec   ->Fill(qHEC);  
  EMECspec  ->Fill(qEMEC);  



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


  } // 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.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("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);
  }


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

  myLine->SetLineColor(kBlack);
  
  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 (iqq<0) {
      eventToCatch=-iqq;
      iqq=3;
    }
  }
}
///if (moduleHV[1]>0.) {

sprintf(myText, "Run#%d HEC & EMEC, rates", myRun);
rateCanvas0->Clear();
rateCanvas0->Divide(1,3);
rateCanvas0->SetTitle(myText);
rateCanvas0->Show();
rateCanvas0->cd(1);
minHECEMEC->DrawCopy();
rateCanvas0->cd(2);
HECspec->DrawCopy();
rateCanvas0->cd(3);
EMECspec->DrawCopy();
rateCanvas0->Update();
rateCanvas0->Show();

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

  sampleSize=prehistorySpectrum[iRate].GetEntries()/20;
  if (sampleSize<2000) sampleSize=2000;
  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>sampleSize) {
      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((int)nxb, -0.5, MaxTracks+0.5, 500, -100., 3500.);
    cherVStracksCov[id].Clear();
    //printf("arrays and hists for bin#%3d of rate%2d is ready, id=%d\n", ib, iRate, id); 
  }
}
// compute response bins for HEC-EMEC position study
for (int ib=0; ib<100; ib++) {
  HECvsEMECbins[ib]=1.0e10;
  HECvsEMECsum0[ib]=0.;
  HECvsEMECsum1[ib]=0.;
  HECbins[ib]=1.0e10;
  HECsum0[ib]=0.;
  HECsum1[ib]=0.;
  EMECbins[ib]=1.0e10;
  EMECsum0[ib]=0.;
  EMECsum1[ib]=0.;
  nHECEMECbins=0;
  nHECbins=0;
  nEMECbins=0;
}
int ibin0=minHECEMEC->FindBin(0.);
for (int i=ibin0; i<minHECEMEC->GetNbinsX(); i++) {
  HECvsEMECsum0[nHECEMECbins]+=minHECEMEC->GetBinContent(i); 
  HECvsEMECsum1[nHECEMECbins]+=minHECEMEC->GetBinContent(i)*minHECEMEC->GetBinCenter(i);
  if (HECvsEMECsum0[nHECEMECbins]>1000. || i==minHECEMEC->GetNbinsX()-1) {
    if (HECvsEMECsum0[nHECEMECbins]>0.) { // may be Zero for last bin!
      HECvsEMECsum1[nHECEMECbins]/=HECvsEMECsum0[nHECEMECbins];
      HECvsEMECbins[nHECEMECbins]=minHECEMEC->GetBinCenter(i)+0.5*minHECEMEC->GetBinWidth(i);
      printf("  HEC_EMEC bin#%d, N=%.0f <A>=%.1f, \tAmax=%.1f\n", 
              nHECEMECbins, HECvsEMECsum0[nHECEMECbins], 
              HECvsEMECsum1[nHECEMECbins], HECvsEMECbins[nHECEMECbins]);
      nHECEMECbins++;
    }
  } 
}
ibin0=HECspec->FindBin(0.);
for (int i=ibin0; i<HECspec->GetNbinsX(); i++) {
  HECsum0[nHECbins]+=HECspec->GetBinContent(i); 
  HECsum1[nHECbins]+=HECspec->GetBinContent(i)*HECspec->GetBinCenter(i);
  if (HECsum0[nHECbins]>1000. || i==HECspec->GetNbinsX()-1) {
    if (HECsum0[nHECbins]>0.) { // may be Zero for last bin!
      HECsum1[nHECbins]/=HECsum0[nHECbins];
      HECbins[nHECbins]=HECspec->GetBinCenter(i)+0.5*HECspec->GetBinWidth(i);
      printf("  HEC bin#%d, N=%.0f <A>=%.1f, \tAmax=%.1f\n", 
              nHECbins, HECsum0[nHECbins], 
              HECsum1[nHECbins], HECvsEMECbins[nHECbins]);
      nHECbins++;
    }
  } 
}
ibin0=EMECspec->FindBin(0.);
for (int i=1; i<EMECspec->GetNbinsX(); i++) {
  EMECsum0[nEMECbins]+=EMECspec->GetBinContent(i); 
  EMECsum1[nEMECbins]+=EMECspec->GetBinContent(i)*EMECspec->GetBinCenter(i);
  if (EMECsum0[nEMECbins]>1000. || i==EMECspec->GetNbinsX()-1) {
    if (EMECsum0[nEMECbins]>0.) { // may be Zero for last bin!
      EMECsum1[nEMECbins]/=EMECsum0[nEMECbins];
      EMECbins[nEMECbins]=EMECspec->GetBinCenter(i)+0.5*EMECspec->GetBinWidth(i);
      printf("  EMEC bin#%d, N=%.0f <A>=%.1f, \tAmax=%.1f\n", 
              nEMECbins, EMECsum0[nEMECbins], 
              EMECsum1[nEMECbins], EMECbins[nEMECbins]);
      nEMECbins++;
    }
  } 
}

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

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

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

// save GADC spectra
sprintf(myText, "%srun%04dv20140514spectraPass1.root", postscriptDirectorry, myRun);
TFile fHist(myText, "RECREATE");
printf("[RunReport] root file: %s\n", myText);
for (int id=0; id<33; id++) {
  trackSpectrum[id].Write();
}
minHECEMEC->Write();
HECspec   ->Write();
EMECspec  ->Write();
fHist.Close();

// show GADC spectra
for (int nSum=0; nSum<3; nSum++) {
  GADCspectra00->Clear();
  GADCspectra00->Divide(2,5);
  
  for (int nTrack=0; nTrack<11; nTrack++) {
    GADCspectra00->cd(nTrack+1)->SetGrid();
    int id=nTrack*3+nSum;
    trackSpectrum[id].GetXaxis()->SetRangeUser(-200.,1000.);
    trackSpectrum[id].DrawCopy();
    double thisMax=trackSpectrum[id].GetMaximum();
    double zeroMax=trackSpectrum[33+nSum].GetMaximum();
    trackSpectrum[33+nSum].Scale(thisMax/zeroMax);
    trackSpectrum[33+nSum].DrawCopy("same");
  }


  GADCspectra00->Update();
  if (iqq>0) {
    if (nSum==0) {
      sprintf(myText, "%sModuleCorrelationRun%04d_GADC_0514.ps(", postscriptDirectorry, myRun);
    } else if (nSum==1) {
      sprintf(myText, "%sModuleCorrelationRun%04d_GADC_0514.ps",  postscriptDirectorry, myRun);
    } else {
      sprintf(myText, "%sModuleCorrelationRun%04d_GADC_0514.ps)", postscriptDirectorry, myRun);
    }
    GADCspectra00->Print(myText);
  }
}
//
//
//
//
//
//
printf("\n\n\n Pass#2 \n\n\n");
//iqq=1;
//
//
//
//
//
//
//
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;
        if (iqq<0) {
          eventToCatch=-iqq;
          iqq=3;
        }
      }

    } 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

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




//=============================================== 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++) {
      
      evtMIPs+=evtEssentials.mipCountX(j);
      evtGADC+=evtEssentials.GADCevent[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.mipCountX(j), evtEssentials.GADCevent1[j]);
            cherVStracksCov[id].Fill((double)evtEssentials.mipCountX(j), evtEssentials.GADCevent1[j]);
            RateBinsSum0[id]+=1.0;
            RateBinsSum1[id]+=historyRate[iRate];
            break;
          }
        }
      }  // rates loop
    }    // samples-in-cluster loop
    
  } // clusters loop
  
  evtGADC*=BH.ChGADC;
//=============================================== fill HEC-EMEC position correlations  
  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 qMIN=qHEC;
  if (qEMEC<qMIN) qMIN=qEMEC;
  if (qMIN<=0.) continue;
  int ibin=0;
  for (int i=0; i<nHECEMECbins; i++) { 
    ibin=i;
    if (qMIN<HECvsEMECbins[i]) break;
  }
  if (iqq==1) {
    printf(" HEC=%.1f\t EMEC=%.1f bin#%d\n", qHEC, qEMEC, ibin); 
  }
  HECdEMECx[ibin] .Fill(xHEC-xEMEC);
  HECdEMECy[ibin] .Fill(yHEC-yEMEC);
  HECvsEMECx[ibin].Fill(xEMEC, xHEC);
  HECvsEMECy[ibin].Fill(yEMEC, yHEC);
  
  hecVSemecX[ibin].Fill(xEMEC, xHEC);
  hecVSemecY[ibin].Fill(yEMEC, yHEC);
  if (ibin>1) { // select large signals!
    xProfile->Fill(xEMEC);
    yProfile->Fill(yEMEC);
  }

} // end of second loop over events
//

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]<100.) 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);

slopeCanvasXX->Clear();
slopeCanvasXX->Divide(2,2);
int iRate0=1;
int iRate1=3;
int iRate2=7;


for (int iRate=0; iRate<9; iRate++) {  
  slopeCanvas00->Clear();
  slopeCanvas00->Divide(2,2);
  slopeCanvas01->Clear();
  slopeCanvas01->Divide(1,2);
  slopeCanvas01->cd(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)->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(kBlack);
  myLatex->DrawLatex(0.05*maxRate, 2.50/BH.ChGADC, "Black: free fit");

  TH2F *temp = (TH2F*) cherVStracks2D[100*iRate].Clone("temp");
  temp->Reset();
  int iXmax=0; 
  int iYmax=0;
  temp->SetStats(0);

  iZone=0;
  int iqq2=1;
  rateCut[iRate]=1.0e10;
  slopePass0[iRate]=0.;
  slopePass1[iRate]=0.;

  for (int ib=0; ib<nRateBins[iRate]; ib++) {
    int id=100*iRate+ib;
    if (RateBinsSum0[id]<100.) continue;
    if (ib>3) {
      double slopeProbe=slopePass1[iRate]/slopePass0[iRate];
      if (cherVStracksCov[id].GetSlopeX0(10)>2.0*slopeProbe) {
        rateCut[iRate]=RateBinsSum1[id];
        break;  
      }
    }
    slopePass0[iRate]+=RateBinsSum0[id];
    slopePass1[iRate]+=cherVStracksCov[id].GetSlopeX0(10)*RateBinsSum0[id];
  }  
  printf("Prehistory#%d rateCut=%.2fx10^{6}p/sec slope=%.1f\n", iRate,  rateCut[iRate], slopePass1[iRate]/slopePass0[iRate]);
  if (iRate==0) {
    slopeVSrun[runIndex]=BH.ChGADC*slopePass1[iRate]/slopePass0[iRate];
  }

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

    myLine->SetLineWidth(3);
    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);
    myLatex->SetTextColor(kBlack);
    myLine ->SetLineColor(kBlack);
    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); 
    double dx=0.05*(x1-x0);
    for (double xx0=2.*x0-x1; xx0<2*x1-x0; xx0+=dx) {
      double xx1=xx0+dx;
      double yy0= cherVStracksCov[id].GetAverageY(10)+cherVStracksCov[id].GetSlopeX(10)*(xx0-cherVStracksCov[id].GetAverageX(10));   
      double yy1= cherVStracksCov[id].GetAverageY(10)+cherVStracksCov[id].GetSlopeX(10)*(xx1-cherVStracksCov[id].GetAverageX(10));
      if (xx0>x0 && yy0>y0 && xx1<x1 && yy1<y1) myLine->DrawLine(xx0, yy0, xx1, yy1);   
    }

    myLatex->SetTextColor(kRed);
    myLine ->SetLineColor(kRed);
    sprintf(myText, "slopeX0=%.2f",    cherVStracksCov[id].GetSlopeX0(10));
    myLatex->DrawLatex(0.95*x0+0.05*x1, 0.75*y1+0.25*y0, myText);
    dx=0.05*(x1-x0);
    for (double xx0=2.*x0-x1; xx0<2*x1-x0; xx0+=dx) {
      double xx1=xx0+dx;
      double yy0= cherVStracksCov[id].GetSlopeX0(10)*xx0;   
      double yy1= cherVStracksCov[id].GetSlopeX0(10)*xx1;
      if (xx0>x0 && yy0>y0 && xx1<x1 && yy1<y1) myLine->DrawLine(xx0, yy0, xx1, yy1);   
    }

    myLine->SetLineWidth(1);
    myLatex->SetTextSize(0.05);
    myLatex->SetTextAlign(22);
    myLatex->SetTextColor(kBlack);
    myLine ->SetLineColor(kBlack);
    slopeCanvas01->cd(1);
    myLatex->DrawLatex(RateBinsSum1[id],     cherVStracksCov[id].GetCorrelation(10), "x");
    myLine ->DrawLine(RateBinsSum1[id], cherVStracksCov[id].GetCorrelation(10)+cherVStracksCov[id].GetCorrelationVariance(),
                      RateBinsSum1[id], cherVStracksCov[id].GetCorrelation(10)-cherVStracksCov[id].GetCorrelationVariance());
    slopeCanvas01->cd(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");
      myLine ->DrawLine(RateBinsSum1[id],      cherVStracksCov[id].GetSlopeX(10)-cherVStracksCov[id].GetSlopeXVariance(),
                        RateBinsSum1[id],      cherVStracksCov[id].GetSlopeX(10)+cherVStracksCov[id].GetSlopeXVariance());
    }
    myLatex->SetTextColor(kRed);
    myLine ->SetLineColor(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.1/BH.ChGADC ,      "#Downarrow");
    } else {
      myLatex->DrawLatex(RateBinsSum1[id],     cherVStracksCov[id].GetSlopeX0(10),     "x");
      myLine ->DrawLine(RateBinsSum1[id],      cherVStracksCov[id].GetSlopeX0(10)-cherVStracksCov[id].GetSlopeX0Variance(),
                        RateBinsSum1[id],      cherVStracksCov[id].GetSlopeX0(10)+cherVStracksCov[id].GetSlopeX0Variance());
    }

    
    myLatex->SetTextSize(0.06);
    slopeCanvas00->Update();
    slopeCanvas01->Update();
    if (iZone==4 || ib==nRateBins[iRate]-1) {
      iZone=0;
      if (iqq>0 && iqq1>0) {
        if (firstPlot) {
          firstPlot=false;
          sprintf(myText, "%sModuleCorrelationRun%04devtSum0514pass2.ps(", postscriptDirectorry, myRun);
        } else {
          sprintf(myText, "%sModuleCorrelationRun%04devtSum0514pass2.ps", postscriptDirectorry, myRun);
        }
        slopeCanvas00->Print(myText);
        slopeCanvas00->Clear();
        slopeCanvas00->Divide(2,2);
      }
    }
  } // rate bins loop
  slopeCanvas01->Print(myText);

  double x0t=temp->GetXaxis()->GetXmin();
  double x1t=temp->GetXaxis()->GetBinCenter(iXmax);
  double y0t=temp->GetYaxis()->GetXmin();
  double y1t=temp->GetYaxis()->GetBinCenter(iYmax);
  y0t=-10.;
  temp->GetXaxis()->SetRangeUser(x0t, x1t);
  temp->GetYaxis()->SetRangeUser(y0t, y1t);

  bool plotItXX=false;

  if (iRate==iRate0) {
     slopeCanvasXX->cd(1);
     gPad->SetLogz(1);
     temp->SetTitle("All rates");
     temp->DrawCopy("");
     temp->Reset();
     slopeCanvasXX->cd(2);
     plotItXX=true;
  } else if (iRate==iRate1) {
     slopeCanvasXX->cd(3);
     plotItXX=true;
  } else if (iRate==iRate2) {
     slopeCanvasXX->cd(4);
     plotItXX=true;
  }
  if (plotItXX) {
     sprintf(myText, "Integration Time = %.2fms", historyTime[iRate]); 
     temp->SetTitle(myText); 
     temp->DrawCopy("box");
     int iCol=1;
     int nPlots=0;
     for (int ib=nRateBins[iRate]-1; ib>0 && nPlots<8; ib-=2) {
       int id=100*iRate+ib;
       if (RateBinsSum0[id]<100.) continue;
       cherVStracks2D[id].SetMarkerColor(iCol);
       cherVStracks2D[id].SetFillColor(iCol);
       cherVStracks2D[id].GetXaxis()->SetRangeUser(x0t, x1t);
       cherVStracks2D[id].GetYaxis()->SetRangeUser(y0t, y1t);
       cherVStracks2D[id].DrawCopy("same");
       nPlots++;
       slopeCanvasXX->Update();
       if (iqq2==1) {
//         printf("Color=%d push->", iCol);
//         cin>>iqq2;
         if (iqq2==0) {
           iqq=0;
           break;
         } else if (iqq2==-1) {
           exit(-1);
         }
       }
       if (nPlots<6) iCol++;
       if (iCol==5) iCol++;
       if (iCol==8) iCol=1;
       slopeCanvasXX->Update();
     }
  }

  temp->Delete();
   
  if (iqq1==0 || iqq==0) break;
} //  integration time loop
sprintf(myText, "%sModuleCorrelationRun%04devtSum0514pass2.ps)", postscriptDirectorry, myRun);
slopeCanvasXX->Print(myText);
if (iqq==0) break;


positionCanvas0->Clear();
positionCanvas0->Divide(2, 5);
positionCanvas1->Clear();
positionCanvas1->Divide(2, 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);
positionCanvas0->Show();
positionCanvas1->Show();

iZone=1;
iqq1=2;
bool firsttime=true;
myLatex->SetTextSize(0.07);
for (int ib=0; ib<nHECEMECbins; ib++) {
  positionCanvas0->cd(2*iZone-1)->SetGrid();
  gStyle->SetOptStat(111110);
  sprintf(myText, "min(HEC, EMEC)=%.1ffC", HECvsEMECsum1[ib]);
  HECvsEMECx[ib].SetTitle(myText);
  HECvsEMECx[ib].DrawCopy("box");
  sprintf(myText, "corr=%.2f", hecVSemecX[ib].GetCorrelation(10));
  myLatex->SetTextAlign(12);
  myLatex->DrawLatex(-1., 1.0, myText); 
  sprintf(myText, "slope=%.2f#pm%.2f", hecVSemecX[ib].GetSlopeX(10), hecVSemecX[ib].GetSlopeXVariance());
  myLatex->SetTextAlign(12);
  myLatex->DrawLatex(-1., 0.7, myText); 
  
  positionCanvas0->cd(2*iZone)->SetGrid();
  gStyle->SetOptStat(111110);
  HECdEMECx[ib].SetTitle("");
  HECdEMECx[ib].DrawCopy();

  positionCanvas1->cd(2*iZone-1)->SetGrid();
  gStyle->SetOptStat(111110);
  sprintf(myText, "min(HEC, EMEC)=%.1ffC", HECvsEMECsum1[ib]);
  HECvsEMECx[ib].SetTitle(myText);
  HECvsEMECy[ib].DrawCopy("box");
  sprintf(myText, "corr=%.2f", hecVSemecY[ib].GetCorrelation(10));
  myLatex->SetTextAlign(12);
  myLatex->DrawLatex(-1., 1.0, myText); 
  sprintf(myText, "slope=%.2f#pm%.2f", hecVSemecY[ib].GetSlopeX(10), hecVSemecY[ib].GetSlopeXVariance());
  myLatex->SetTextAlign(12);
  myLatex->DrawLatex(-1., 0.7, myText); 

  positionCanvas1->cd(2*iZone)->SetGrid();
  gStyle->SetOptStat(111110);
  HECdEMECy[ib].DrawCopy();
  positionCanvas0->Update();
  positionCanvas1->Update();

  if (iZone==5 || ib==nHECEMECbins-1) {
    iZone=0;
    
    sprintf(myTex1, "%sModuleCorrelationRun%04devtSum0514_xy.ps", postscriptDirectorry, myRun);
    sprintf(myTex2, "%sModuleCorrelationRun%04devtSum0514_xy.ps", postscriptDirectorry, myRun);
    if (firsttime) {
      firsttime=false;
      sprintf(myTex1, "%sModuleCorrelationRun%04devtSum0514_xy.ps(", postscriptDirectorry, myRun);
    }
    if (ib==nHECEMECbins-1) {
      sprintf(myTex2, "%sModuleCorrelationRun%04devtSum0514_xy.ps)", postscriptDirectorry, myRun);
    }
    positionCanvas0->Print(myTex1);
    positionCanvas1->Print(myTex2);
    
    if (iqq1==1) {
      printf(" HEC-EMEC position correlations show, push->");
      cin>>iqq1;
    }
    if (iqq1==0) {
      iqq=0;
      break;
    }
    positionCanvas0->Clear();
    positionCanvas1->Clear();
    positionCanvas0->Divide(2, 5);
    positionCanvas1->Divide(2, 5);
  } else {
    if (iqq1==1) {
      printf("push->");
      cin>>iqq1;
    }
    if (iqq1==0) {
      iqq=0;
      break;
    }
  }
  iZone++;
}

//========================================================================
myRunX[runIndex]=myRun;
xVSrun[runIndex]=xProfile->GetMean();
yVSrun[runIndex]=yProfile->GetMean();
globalSummary->Clear();
globalSummary->Divide(1,2);
globalSummary->cd(1)->SetGrid(1);
summaryGainDummy->SetBins(runIndex+1, -0.5, runIndex+0.5);
summaryGainDummy->SetStats(0);
summaryGainDummy->SetLineStyle(2);
summaryGainDummy->DrawCopy();
myLatex->SetTextAlign(22);

for (int i=0; i<=runIndex; i++) {
  myLatex->SetTextColor(kBlack);
  if (myRunX[i]<1050) myLatex->SetTextColor(kBlue);
  myLatex->DrawLatex((double)i, slopeVSrun[i], "x");
}

globalSummary->cd(2)->SetGrid(1);
summaryPositionDummy->SetBins(runIndex+1, -0.5, runIndex+0.5);
summaryPositionDummy->SetStats(0);
summaryPositionDummy->DrawCopy();

myLatex->SetTextColor(kBlack);
if (moduleHV[1]==0.) { // HV EMEC not set?
  myLatex->SetTextColor(kGray);
} 
for (int i=0; i<=runIndex; i++) {
  myLatex->DrawLatex((double)i, xVSrun[i], "x");
}

myLatex->SetTextColor(kRed);
if (moduleHV[1]==0.) { // HV EMEC not set?
  myLatex->SetTextColor(kGray);
} 
for (int i=0; i<=runIndex; i++) {
  myLatex->DrawLatex((double)i, yVSrun[i], "y");
}
globalSummary->Show();
globalSummary->Update();
sprintf(myText, "%sglobalSummary.ps", postscriptDirectorry);
globalSummary->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
  if (iqq<0) {
    eventToCatch=-iqq;
    iqq=3;
  }
}
f->Close();

runIndex++;

} // end of one run loop

}


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