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

struct eventEssentials {

int firstGADCsample; // first Gated ADC sample for current event

float GADCevent[40]; // 40 Gated ADC samples (RF cycles), starting with firstGADCsample
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 miniModule[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++) {
    miniModule[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)

   void Clear() {
     ns30=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, entire event
  double Y;        // Y center of gravity, entire event
  double Xrf[5];   // X center of gravity, per cluster
  double Yrf[5];   // Y center of gravity, per cluster

};


char *module_name[7] = { 
  "HEC", 
  "EMEC", 
  "FCal025", 
  "FCal010", 
  "Beam",
  "Monitor",
  "Integrator"
};

///=============================================================================
void HiLumRead2014_04_14() {


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

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

pedestalWindows  pedWindows[25];
moduleSignalsF   Module[50];

scintMIPs        s123mips;

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

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

double time[512];
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; Module[iG+2*iCh+10*iMod].X=-1.; Module[iG+2*iCh+10*iMod].Y= 1.;
  iCh=1; Module[iG+2*iCh+10*iMod].X= 1.; Module[iG+2*iCh+10*iMod].Y= 1.;
  iCh=2; Module[iG+2*iCh+10*iMod].X=-1.; Module[iG+2*iCh+10*iMod].Y=-1.;
  iCh=3; Module[iG+2*iCh+10*iMod].X= 1.; Module[iG+2*iCh+10*iMod].Y=-1.;
  iMod=1; // EMEC
  iCh=0; Module[iG+2*iCh+10*iMod].X=-1.; Module[iG+2*iCh+10*iMod].Y= 1.;
  iCh=1; Module[iG+2*iCh+10*iMod].X=-1.; Module[iG+2*iCh+10*iMod].Y=-1.;
  iCh=2; Module[iG+2*iCh+10*iMod].X= 1.; Module[iG+2*iCh+10*iMod].Y= 1.;
  iCh=3; Module[iG+2*iCh+10*iMod].X= 1.; Module[iG+2*iCh+10*iMod].Y=-1.;
  iMod=2; // FCal-250
  iCh=0; Module[iG+2*iCh+10*iMod].X= 1.; Module[iG+2*iCh+10*iMod].Y= 1.;
  iCh=1; Module[iG+2*iCh+10*iMod].X=-1.; Module[iG+2*iCh+10*iMod].Y= 1.;
  iCh=2; Module[iG+2*iCh+10*iMod].X= 1.; Module[iG+2*iCh+10*iMod].Y=-1.;
  iCh=3; Module[iG+2*iCh+10*iMod].X=-1.; Module[iG+2*iCh+10*iMod].Y=-1.;
  iMod=3; // FCal-100
  iCh=0; Module[iG+2*iCh+10*iMod].X= 1.; Module[iG+2*iCh+10*iMod].Y= 1.;
  iCh=1; Module[iG+2*iCh+10*iMod].X=-1.; Module[iG+2*iCh+10*iMod].Y= 1.;
  iCh=2; Module[iG+2*iCh+10*iMod].X= 1.; Module[iG+2*iCh+10*iMod].Y=-1.;
  iCh=3; Module[iG+2*iCh+10*iMod].X=-1.; Module[iG+2*iCh+10*iMod].Y=-1.;
}

TCanvas *moduleCanvas  = new TCanvas("moduleCanvas",  "Module Signals",    0, 0, 800, 1100);
TCanvas *moduleCanvas1 = new TCanvas("moduleCanvas1", "Module Currents", 800, 0, 800, 1100);


int myRun0=850;
int myRun1=1200;

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

myRun0=1098;
myRun1=1099;

for (int myRun=myRun0; myRun<myRun1; myRun++) { // runs looper
int burstSet=-1;

// create a new ROOT file
sprintf(rootName, "reports/ntuple/run%04dv20140404tree.root", myRun);


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

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.miniModule[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, &Module[id].ampI[0]);
      if (iMod != 4) { // for scintillators - different story!!!!
        sprintf(branchName, "current_mod%dchan%dgain%d", iMod, iCh, iG);
        eventTree -> SetBranchAddress(branchName, &Module[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"); 
printf(" blank module#%d\n", blankModule);
printf(" S123 in beam =%d\n\n", S123inBeam);
//==============================================================
int burstEnt=0;

for (int iEnt=0; iEnt<eventEntries; iEnt++) {
  evtEssentials.Clear();
  eventTree->GetEntry(iEnt);
  if (iqq<3) {
    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);
    for (int i=0; i<40; i++) printf("%3d ", evtEssentials.isCluster0[i]);
    printf("\n");
    for (int i=0; i<40; i++) printf("%3d ", evtEssentials.isCluster1[i]);
    printf("\n");
    for (int i=0; i<40; i++) printf("%3d ", evtEssentials.isIsolated[i]);
    printf("\n");
  }

  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;
    } else {
      printf("Burst Entry #%d is out of range: nBursts=%d\n", burstEnt, burstEntries);
      continue;
    }
  }

  bool HGoverflow=false;
  for (int iMod=0; iMod<4; iMod++) {
    HGoverflow = HGoverflow || evtEssentials.HGoverflow[iMod];
  }
  if (HGoverflow && iqq!=2) {iqqOld=iqq; iqq=1;}
  else            iqq=iqqOld;

  int cluster=0;
  for (int iMod=0; iMod<4; iMod++) {
    cluster=0;
    for (int ir=0; ir<40; ir++) {
      if (evtEssentials.isCluster0[ir]) { // beginning of cluster or single //
        for (int iCh=0; iCh<4; iCh++) {
          if (HGmax<evtEssentials.miniModule[1+2*iCh+10*iMod+40*cluster]) HGmax=evtEssentials.miniModule[1+2*iCh+10*iMod+40*cluster];
          if (HGmin>evtEssentials.miniModule[1+2*iCh+10*iMod+40*cluster]) HGmin=evtEssentials.miniModule[1+2*iCh+10*iMod+40*cluster];
        }
        cluster++;
      }
    }
  }

  if (MH.isEmpty || iqq==3) continue;

  moduleCanvas->Clear();
  moduleCanvas->Divide(1,4);
  moduleCanvas1->Clear();
  moduleCanvas1->Divide(1,4);
  int iZone=1;
  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, Module[id].ampI[i], Module[id].ampF[i]);
      if (amin >Module[id].ampI[i]) amin =Module[id].ampI[i]; 
      if (amax <Module[id].ampI[i]) amax =Module[id].ampI[i];
      if (aminF>Module[id].ampF[i]) aminF=Module[id].ampF[i]; 
      if (amaxF<Module[id].ampF[i]) amaxF=Module[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);
    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], Module[id].ampI[i-1], time[i], Module[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->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++) {
      myLine->DrawLine(time[is0+i-1], Module[id].ampF[i-1], time[is0+i], Module[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], Module[id2].ampF[i-1], time[is0+i], Module[id2].ampF[i] ); 
      }
    }

    cluster=0;
    for (int ir=0; ir<40; ir++) {
      if (evtEssentials.isCluster0[ir]) { // beginning of cluster or single //
        printf("%10s over=%1d Cluster#%d ->", module_name[iMod], evtEssentials.HGoverflow[iMod], cluster);
        for (int iCh=0; iCh<5; iCh++) {
          printf(" Ch#%d=(%5.1f, %5.1f)", iCh,
                   gain*evtEssentials.miniModule[1+2*iCh+10*iMod+40*cluster],  
                   gain*evtEssentials.miniModule[0+2*iCh+10*iMod+40*cluster]);
        }
        printf("\n");
        cluster++;
      }
    }

    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->DrawCopy();
  for (int i=0; i<40 && 165.*(i+1)<6400.; i++) {
    myLine->DrawLine(165.*i, evtEssentials.GADCevent[i], 165.*(i+1), evtEssentials.GADCevent[i]);
    if (i>0) {
      myLine->DrawLine(165.*i, evtEssentials.GADCevent[i-1], 165.*i, evtEssentials.GADCevent[i]);
    }
  }

// 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., mipCountMAX+0.5);
  moduleCanvas1 ->cd(iZone);
  dummy->SetTitle("S123 trac counting");
  dummy->DrawCopy();
  myLine->SetLineColor(kBlack);
  for (int i=0; i<40 && 165.*(i+1)<6400.; i++) {
    myLine->DrawLine(165.*i, evtEssentials.mipCount[i], 165.*(i+1), evtEssentials.mipCount[i]);
    if (i>0) {
      myLine->DrawLine(165.*i, evtEssentials.mipCount[i-1], 165.*i, evtEssentials.mipCount[i]);
    }
  }
  myLine->SetLineColor(kRed);
  for (int i=0; i<40 && 165.*(i+1)<6400.; i++) {
    myLine->DrawLine(MH.TimeGlobal+165.*i-10.,     evtEssentials.mipCount2[i], 
                     MH.TimeGlobal+165.*(i+1)-10., evtEssentials.mipCount2[i]);
    if (i>0) {
      myLine->DrawLine(MH.TimeGlobal+165.*i-10., evtEssentials.mipCount2[i-1], 
                       MH.TimeGlobal+165.*i-10., evtEssentials.mipCount2[i]);
    }
  }
  myLine->SetLineColor(kBlue);
  for (int i=0; i<40 && 165.*(i+1)<6400.; i++) {
    myLine->DrawLine(MH.TimeGlobal+165.*i+10.,     evtEssentials.mipCount3[i], 
                     MH.TimeGlobal+165.*(i+1)+10., evtEssentials.mipCount3[i]);
    if (i>0) {
      myLine->DrawLine(MH.TimeGlobal+165.*i+10., evtEssentials.mipCount3[i-1], 
                       MH.TimeGlobal+165.*i+10., evtEssentials.mipCount3[i]);
    }
  }
  



  myLine->SetLineColor(kBlack);
  
  moduleCanvas ->Update();
  moduleCanvas1->Update();
  if (iqq < 2) {
    printf("push->");
    cin >> iqq;
    if (iqq==0) break;
  }
}
f->Close();

printf("HGmin=%.1f HGmax=%.1f \n", HGmin, HGmax);
}
}


/*


        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.miniModule[0],   "miniModule[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 {
              }
            }
          }
        }
*/