/*
 * methods for SiPM CAEN readout output
 *
 * First Version Jul/2023
 *
 */

//#include <unistd.h>
#include <fstream>

#include <Riostream.h>
#include <TSystem.h>
#include <TCanvas.h>
#include <TVector3.h>
#include <TStyle.h>
#include <TGraph.h>
#include <TTree.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TF1.h>
#include <TStyle.h>
#include <TMath.h>
#include <TProfile.h>
#include <TString.h>
#include <TText.h>
#include <TLegend.h>
#include <TLegendEntry.h>
#include <TCut.h>
#include <TLine.h>
#include <TObjArray.h>
#include <TObjString.h>
#include <TFile.h>
#include <TChain.h>
#include <TParameter.h>
#include <TEntryList.h>

#include "scintifers.h"

/*
 *
 */

void setLocalStyle() {
  gStyle->SetTitleW(0.4);
  gStyle->SetTitleH(0.07);
  gStyle->SetTitleSize(0.06,"xyzt");
  gStyle->SetTitleXOffset(.8);
  gStyle->SetTitleYOffset(.8);
  gStyle->SetLabelSize(0.04,"XYZ");
  gStyle->SetPadTopMargin(.1);
  gStyle->SetPadLeftMargin(.1);
  gStyle->SetPadRightMargin(.05);
  gStyle->SetPadBottomMargin(.1);
  gStyle->SetOptStat(0);
  gStyle->SetOptFit(0);
}

/*
 *
 */

int singleBarSpe(TChain *ctl,
		 int selch,  // one of the two channels at the scintillator bar ends
		 int thr_gain=8192
		 )
{

  int ch[2];
  ch[0]=selch;
  // get paired channels
  ch[1]=pairedChannel(ch[0]);
  
  TCanvas *c1 = new TCanvas("c1",Form("Chs: %d %d Side: %d, Pos: %d",ch[0],ch[1],c2side(ch[0]), c2idx(ch[0])));
  c1->Divide(3,2);
  c1->Update();

  TCut cpair=Form("((chan==%d)||(chan==%d))",ch[0],ch[1]);
  
  // channels distribution
  TCut chgain = Form("(hgain<%d)",thr_gain);
  TCut clgain = Form("(lgain<%d)",thr_gain);

  ctl->SetMarkerStyle(1);
  ctl->SetLineWidth(1);
  ctl->SetLineColor(18);
  c1->cd(1);
  ctl->Draw("hgain",chgain&&cpair);

  c1->cd(2);
  ctl->Draw("lgain",clgain&&cpair);

  ctl->SetLineWidth(2);
  for (int i=0;i<2;i++) {
    c1->cd(1)->SetLogy();
    ctl->SetMarkerStyle(20+i);
    ctl->SetMarkerColor(2+i);
    ctl->SetLineColor(2+i);
    ctl->Draw("hgain",Form("(hgain<%d)&&(chan==%d)",thr_gain,ch[i]),"same");
    c1->cd(2)->SetLogy();
    ctl->Draw("lgain",Form("(lgain<%d)&&(chan==%d)",thr_gain,ch[i]),"same");
  }
  c1->Update();

  // combine the paired channels

  int novt;
  int chan[64];  // 64*nboards
  int lgain[64];
  int hgain[64];

  ctl->SetBranchAddress("novt",&novt);
  ctl->SetBranchAddress("chan",chan);
  ctl->SetBranchAddress("hgain",hgain);
  ctl->SetBranchAddress("lgain",lgain);

  int nn = ctl->GetEntries();

  ctl->Draw(">>elist",cpair&&(chgain&&clgain),"entrylist");
  TEntryList *elist = (TEntryList*) gDirectory->Get("elist");
  
  int nen = elist->GetN();

  TGraph *g0cor = new TGraph(nen);
  g0cor->SetNameTitle("ghcor",Form("High-Low Gain Correlation %d",ch[0]));
  g0cor->SetMarkerStyle(22);
  g0cor->SetMarkerColor(6);
  TGraph *g1cor = new TGraph(nen);
  g1cor->SetNameTitle("glcor",Form("High-Low Gain Correlation %d",ch[1]));
  g1cor->SetMarkerStyle(23);
  g1cor->SetMarkerColor(9);
  TGraph *gtcor = new TGraph(nen);
  gtcor->SetNameTitle("gtcor","High vs Low Total Gain Correlation");
  gtcor->SetMarkerStyle(34);
  TGraph *gdcor = new TGraph(nen);
  gdcor->SetNameTitle("gdcor","High vs Low Delta Gain Correlation");
  gdcor->SetMarkerStyle(33);
  int nact=0;

  TH1F *hhg = new TH1F("hhg","High Gain Centroid",100,-1.1,1.1);
  hhg->SetLineColor(2);
  TH1F *hlg = new TH1F("hlg","Low Gain Centroid",100,-1.1,1.1);
  hlg->SetLineColor(4);
    
  int hcharge[2]={0,0};
  int lcharge[2]={0,0};

  printf(" Total Entries %d, after cuts: %d\n",nn,nen);

  //  ctl->Reset();
  //  ctl->SetEntryList(elist);
  
  for (int i=0;i<nn;i++) {
    //    int ei = elist->GetEntry(i);    
    ctl->GetEntry(i);
    int nf=0;
    for (int h=0;h<novt;h++) {
      for (int k=0;k<2;k++) {
	if (chan[h]==ch[k]) {
	  if ((hgain[h]>=thr_gain)||(lgain[h]>=thr_gain)) { break; }
	  hcharge[k] = hgain[h];
	  lcharge[k] = lgain[h];
	  nf+=1;
	}
      }
      if (nf==2) {
	g0cor->SetPoint(nact,lcharge[0],hcharge[0]);
	g1cor->SetPoint(nact,lcharge[1],hcharge[1]);
	float htot=((float) (hcharge[0]+hcharge[1]));
	float ltot=((float) (lcharge[0]+lcharge[1]));
	float hdelta = ((float) (hcharge[0]-hcharge[1]));
	float ldelta = ((float) (lcharge[0]-lcharge[1]));
	gtcor->SetPoint(nact, ltot, htot);
	gdcor->SetPoint(nact, ldelta, hdelta);
	if (htot>0) {
	  hhg->Fill(hdelta/htot);
	}
	if (ltot>0) {
	  hlg->Fill(ldelta/ltot);
	}
	nact++;
	break;
      }
    } // end loop on channels in single entry
  } // end loop on selected entries

  for (int i=nact;i<nen;i++) { //remove unfilled points in TGraph
    g0cor->RemovePoint(i);
    g1cor->RemovePoint(i);
    gtcor->RemovePoint(i);
  }
  
  c1->cd(3);
  g0cor->Draw("PAW");
  g1cor->Draw("P");
  c1->cd(4);
  gtcor->Draw("PAW");
  c1->cd(5);
  gdcor->Draw("PAW");

  c1->cd(6);
  hhg->Draw();
  hlg->Draw("same");
  
  c1->Update();
  
  return 0;
  
  // channels correlation

  // weighted centroids

  //
};


/*
 *
 *
 */

TGraph * plotActiveSChannels(TChain *ctl, TCanvas *cc, TString svar="hgain", TCut gcut="hgain<30000") {

  static int count=0;
  
  float scale=0.1;

  TGraph *gmearms = new TGraph(64); // mean and rms of each distribution
    
  int pflag[16]; // 1 if pad has been already used
  for (int i=0;i<16;i++) { pflag[i]=0; }

  TText *text = new TText(0,0,"mean / rms");
  text->SetTextColor(1);
  text->SetTextSize(0.15);

  cc->Divide(2,8);
  cc->Update();

  cc->cd(1)->SetLogy();
  
  ctl->SetMarkerStyle(1);
  ctl->Draw(Form("%s>>h%ddum",svar.Data(),count),gcut);

  TH1F *hdum = (TH1F*) gDirectory->Get(Form("h%ddum",count));
  hdum->SetTitle("High Gain Charge Distribution");
  hdum->SetXTitle("ADC ch");
  hdum->SetMarkerStyle(1);
  hdum->Scale(scale);
  
  for (int i=0;i<64;i++) {
    int is = c2side(i);
    int ip = c2idx(i);

    gmearms->SetPoint(i,-1,-1); // default value
    
    if (is<0) { continue; }
    
    int icd = ip*2 + (is % 2);
    cc->cd(icd+1)->SetLogy();
    if (pflag[icd]==0) {
      hdum->Draw("p");
      pflag[icd]=1;
    }
    
    ctl->SetLineColor(is+1);
    ctl->SetMarkerColor(is+1);
    ctl->SetMarkerStyle(20+is);
    ctl->Draw(Form("%s>>h%ddus%d",svar.Data(),count,i),Form("chan==%d",i),"p,same");
    if (ctl->GetSelectedRows()<=0) { continue; }
    cc->Update();

    TH1F *hf = (TH1F*) gDirectory->Get(Form("h%ddus%d",count,i));
    float mean = hf->GetMean();
    float rms = hf->GetRMS();
    gmearms->SetPoint(i,mean,rms);

    printf(" %d     %f %f\n",i,mean,rms);

    text->SetTextColor(is+1);
    float yn = 0.4 + 0.2*((float) (is/2));
    text->DrawTextNDC(0.6,yn,Form("%.1f +/- %.1f",mean,rms));
    
    count += 1;
    
  }

  return gmearms;
  
}

/* - - - - - - - - - - - - - - - - - -
 * Timing Mode Process methods
 *
 */

/*
 * get pedestal of time over threshold
 * can be used to plot "time over threshold" distributions of each channels;
 * useful tu evaluate equalization of channels responses
 *
 */

vector<float> *getToTPedestal(TTree *tped, // tree with pedestal data (in NULL) use ped=0 and sigma=def_rms
			      vector<int> vich,  // vector of active channels
			      float def_rms // assumed rms where there is no channel data
			      ) {


  vector<float> *vped = new vector<float>[3]; // chananel (float for simplicity), pedestal and rms
  int nchs = (int) vich.size();

  if (tped==NULL) {
    for (int i=0;i<nchs;i++) {
      vped[0].push_back((float) vich[i]);
      vped[1].push_back(0);
      vped[2].push_back(def_rms);
    }
    return vped;
  }
      
  int mincount=500; // at least this number of counts in the channel for pedestal estimation, otherwise use default
  float fscale=3;

  TCanvas *c1 = new TCanvas("c1","Single Channel ToT - Pedestal");
  c1->Divide(2,8,0.0001,0.002); // [x,y][pos]
  c1->Update();

  short pflag[16]; // 1 if first plot drawn (used to overlap plot of channels from both ends of bar)
  for (int i=0;i<16;i++) { pflag[i]=0; }

  c1->cd(1);
  tped->Draw("tot","");
  TH1F *htot = (TH1F*) ((TH1F *) gPad->GetPrimitive("htemp"))->Clone("htot"); // normalize the plots

  htot->Scale(fscale/((float) nchs));
  htot->SetLineColor(10);

  TText *text = new TText(0,0,"mean / rms");
  text->SetTextColor(1);
  text->SetTextSize(0.15);
  
  tped->SetLineWidth(2);
  for (int i=0;i<nchs;i++) {
    int j = (int) vich[i];
    int iside = c2side(j);
    int idx = c2idx(j);
    int cid = (iside%2)+idx*2;
    //    printf("%d %d %d %d\n",j,iside,idx,cid);
    c1->cd(cid+1)->SetLogy();
    if (pflag[cid] == 0) {
      htot->SetTitle(Form("ToT Distribution, axis %d pos %d\n",c2axis(j),idx));
      htot->DrawCopy();
      pflag[cid] += 1;
    }
    tped->SetLineColor(iside/2+1);
    tped->Draw("tot",Form("chan==%d",j),"same");
    int n = tped->GetSelectedRows();
    float mean=0;
    float rms=def_rms;
    if (n>mincount) {
      mean = TMath::Mean(n, tped->GetV1());
      rms = TMath::RMS(n, tped->GetV1());
    }
    
    text->SetTextColor(iside/2+1);
    float yn = 0.5 + 0.2*((float) (iside/2));
    text->DrawTextNDC(0.7,yn,Form("%.1f +/- %.1f (%d)",mean,rms,n));
    // printf(" Ch %d cnt %d : %f %f\n",j, n, mean,rms);
    vped[0].push_back((float) j);
    vped[1].push_back(mean);
    vped[2].push_back(rms);
    
    c1->Update();
    
  };
  c1->Update();

  return vped;
  
};

/*
 *
 */

summaryList* channelTimeProcess(TTree * tl,
				 float thr_toa=100.,   // [.5 ns]  time gate width for coincidence
				 float nsigma=3.,  // [.5 ns] number of sigma for time of threshold width
				 vector<float> *vped=NULL)
{

  //  float thr_tot = 25.; // shaping time is 25 ns
  //  float thr_toa = 100.;
  
  vector<int> vich; // indices of channels acquired

  summaryList *rpl = new summaryList("channelTimeProcess");
  
  tl->Draw("timeus/1e6","","goff"); // s
  double acqtime = TMath::MaxElement(tl->GetSelectedRows(),tl->GetV1()) - TMath::MinElement(tl->GetSelectedRows(),tl->GetV1());

  printf(" Total Acquisition time [s] %f\n",acqtime);
  
  // find unique indeces of acquired channels

  vich = getActiveChannels(tl);

  int nchs = (int) vich.size();
  printf(" Number of active SiPM channels from data: %d\n",nchs);

  //  getToTPedestal(tl, vich, thr_tot);

  int ntref = tl->GetEntries();
  printf(" Total Reference Times (triggers) : %d\n",ntref);

  // find time coincidences
  TCanvas *cv = new TCanvas("cv");
  cv->Divide(2,2);
  cv->Update();

  TH2F *hxy = new TH2F("hxy","x/y coinc Hit map",20,-60,60, 20,-60,60);
  hxy->SetXTitle("x centroid [mm]");
  hxy->SetYTitle("y centroid [mm]");
  
  TH1F *hx = new TH1F("hx","x Hit Map",20,-60,60); //,50,-260,260);
  hx->SetXTitle("x centroid [mm]");
  hx->SetLineColor(1);
  TH1F *hix = new TH1F("hix","Hit Map",8,-0.5,7.5); 
  hix->SetXTitle("axis position index");
  hix->SetLineColor(1);
  
  TH1F *hy = new TH1F("hy","y Hit Map",20,-60,60); //50,-260,260);
  hy->SetXTitle("y centroid [mm]");
  hy->SetLineColor(2);
  TH1F *hiy = new TH1F("hiy","Hit map",8,-0.5,7.5); 
  hiy->SetXTitle("axis position index");
  hiy->SetLineColor(2);
  
  TH1F *hmulti = new TH1F("hmulti","Signal Multiplicity", 20,0,20);
  //  TProfile *hpm = new TProfile("hpm","Signal Multiplicity vs Channel",64,-0.5,63.5,0,50);

  TString scut="";
  TString sand="";
  
  for (uint i=0;i<vped->size();i++) {  // can be improved including the "and" between opposite scinti-bar ends ToT channels
    int ch = (int) vped[0][i];
    float thr = vped[1][i] + nsigma * vped[2][i];
    scut += Form("%s(tot[%d]>%.2f)",sand.Data(),ch,thr);
    sand="||";
  }
  TCut pedcut= scut.Data();
  pedcut.Print();
  
  TCut cbepu="bPulse<6e12"; // low beam pulse cut
  
  for (int itr=0;itr<ntref;itr++) { // loop on reference times (triggers)

    //    tl->Draw("toa:tot:chan",Form("(Entry$==%d)&&(tot>%f)",itr,thr_tot),"goff");
    TCut ecut = Form("(Entry$==%d)",itr);
    
    tl->Draw("toa:tot:chan",ecut&&pedcut,"goff");
    
    int nsig = tl->GetSelectedRows();
    if (nsig<=0) { continue; }
    
    float tamin = TMath::MinElement(nsig, tl->GetV1());
    float tamax = TMath::MaxElement(nsig, tl->GetV1());

    // parameters of the coincidence list for each channel 
    vector<float> vta[2]; // average time of arrival
    vector<float> vxc[2]; // x centroid
    vector<float> vyc[2]; // y centroid
    vector<float> vtq[2]; // total time of overthreshold (sort of total charge)
    vector<int> vns[2]; // number of signals forming the coincidence
    vector<uint32_t> vcb[2]; // each bit corresponds to a firing channel
    
    for (int i=0;i<nsig;i++) { // loop on signals over threshold

      int ch = tl->GetV3()[i]; // channel
      
      //      float thr = vped[0][ch]+vped[1][ch]; // one sigma
      
      float ta = tl->GetV1()[i]; // average time of arrival
      float to = tl->GetV2()[i]; // time over threshold

      float x = c2x(ch);
      float y = c2y(ch);
      int axis = c2axis(ch);
      
      // printf("sig %d %f %f %d %f %f\n",i,ta,to,ch,x,y);

      // hpm->Fill((float) ch,1.);

      int k = axis;
      int idx=-1;
      for (uint j=0;j<vta[k].size();j++) { // loop on stored coincidences of the given axis
	if (TMath::Abs(ta-vta[k][j])<thr_toa) { // 
	  idx = j;
	  break;
	}
      }
      
      if (idx<0) { // new coincidence
	vta[k].push_back(ta);
	vxc[k].push_back(x);
	vyc[k].push_back(y);
	vns[k].push_back(1);
	vtq[k].push_back(to);
	vcb[k].push_back(((uint32_t) 1) << (ch-32*k));
      } else {
	vta[k][idx] = (vta[k][idx]*vns[k][idx]+ta)/(vns[k][idx]+1);
	vns[k][idx] = vns[k][idx]+1;
	vxc[k][idx] = (vxc[k][idx]*vtq[k][idx]+x*to)/(vtq[k][idx]+to);
	vyc[k][idx] = (vyc[k][idx]*vtq[k][idx]+y*to)/(vtq[k][idx]+to);
	vtq[k][idx] = vtq[k][idx]+to;
	vcb[k][idx] = vcb[k][idx] | (((uint32_t) 1) << (ch-32*k));
      }
    } // loop on signals in single event

    // now search coincidences on both ends of scinti-bar
    
    vector<int> vcoi[2]; // x, y coincidence indices
    
    for (int k=0;k<2;k++) { // loop on axes
      
      for (uint j=0;j<vta[k].size();j++) { // loop on coincidences
      
	if (vns[k][j]<2) { // at least 2 SiPM channels firing
	  continue;
	}
	
	uint32_t vb = evalCoincSingle(vcb[k][j]); // coinc on both ends of scinti-bar

	if (vb > 0) {
	  switch (k) {
	  case 0: // "x" axis
	    hx->Fill(vxc[k][j]);
	    for (int h=0;h<8;h++) {
	      if (vb & (1<<h)) {
		hix->Fill(h);
	      }
	    }
	    break;
	  case 1: // "y" axis
	    hy->Fill(vyc[k][j]);
	    for (int h=0;h<8;h++) {
	      if (vb & (1<<h)) {
		hiy->Fill(h);
	      }
	    }
	    break;
	  default:
	    printf("ERROR: trigger %d - axis %d - coinc 0x%x - something wrong !! extra axis ??? \n",itr,k,vcb[k][j]);
	    break;
	  }
	  vcoi[k].push_back(j);
	} // end vb
	
      }
    } // end loop on axes
    
    // x/y coincidences as minimum time difference within thr_toa window

    vector<float> vdtm;
    vector<int> vidx;
    vector<int> vidy;
    for (uint kx=0; kx<vcoi[0].size(); kx++) {
      float dtmin=10*thr_toa;
      int imin=-1;
      for (uint ky=0; ky<vcoi[1].size(); ky++) {
	float dt = TMath::Abs(vta[0][kx]-vta[1][ky]);
	if ((dt<thr_toa) && (dt<dtmin)) {
	  dtmin = dt;
	  imin = ky;
	}
      }
      if (imin>=0) { // time coinc
	int idd=0;
	for (uint h=0;h<vdtm.size(); h++) {
	  if (dtmin<vdtm[h]) {
	    idd = h;
	    break;
	  }
	}
	vdtm.insert(vdtm.begin() + idd, dtmin);
	vidx.insert(vidx.begin() + idd, kx);
	vidy.insert(vidy.begin() + idd, imin);
      }
    } // loop on kx

    for (uint j=0; j<vdtm.size(); j++) { // loop on x/y time coincidence
      int kx = vidx[j];
      int ky = vidy[j];
      if ((vcb[0][kx]>0)&&(vcb[1][ky]>0)) { // avoid multiple counting
	hxy->Fill(vxc[0][kx],vyc[1][ky]);
	vcb[0][kx]=0; // set to zero to avoid multiple counting
	vcb[1][ky]=0;
      }
    }
    
    if ((itr % 100)==0) {
      printf(" progress trigger: %d (%.1f fraction of total)\n",itr,((float) itr)/((float) ntref));
      cv->cd(1);
      hxy->Draw("colz");
      cv->Update();
      cv->cd(2);
      hy->Draw();
      cv->Update();
      cv->cd(3);
      hx->Draw();
      cv->cd(4);
      hiy->Draw();
      hix->Draw("same");
      //      hmulti->Draw();
      //      cv->cd(3);
      //      hpm->Draw();
      cv->Update();
    }

  } // loop on triggers

  // estimate rates
  float scale = 1./((float) ntref);
    
  float totxy = hxy->Integral();
  float totx = hx->Integral();
  float toty = hy->Integral();
  printf("Statistics:     xy      x      y\n");
  printf("Sum of signals: %.4f %.4f %.4f\n",totxy,totx,toty);
  printf("Signal/Trigger: %.4f %.4f %.4f\n",totxy*scale, totx*scale, toty*scale);
  printf("Signal/Time[s]: %.4f %.4f %.4f\n",totxy/acqtime, totx/acqtime, toty/acqtime);

  vector<float> rval;
  rval.push_back(acqtime);
  rval.push_back((float) ntref);
  rpl->addFloat(acqtime,"AcqTime_s");
  rpl->addInt(ntref,"NEntries");
  
  for (int i=0;i<8;i++) {
    float count = hix->GetBinContent(i+1);
    rval.push_back(count);
    rpl->addFloat(count,Form("xCoinc[%d]",i));
  }
  for (int i=0;i<8;i++) {
    float count = hiy->GetBinContent(i+1);
    rval.push_back(count);
    rpl->addFloat(count,Form("yCoinc[%d]",i));
  }

  hxy->Scale(scale);
  hy->Scale(scale);
  hy->SetYTitle("Signals/Triggers");
  hx->Scale(scale);
  hx->SetYTitle("Signals/Triggers");
  cv->cd(1);
  hxy->Draw("colz");
  cv->Update();
  cv->cd(2);
  hy->Draw();
  cv->Update();
  cv->cd(3);
  hx->Draw();
  cv->cd(4);
  hiy->Draw();
  hix->Draw("same");
  cv->Update();
  
  return rpl;
  
};


/*
 * read and plot staircase threshold scan (ScanThr.txt file)
 */

int readThrScan(TString infile) {

  printf(" read input file %s into Ttree\n",infile.Data());
  TTree *tsci = new TTree("tscan","SIPM Threshold Scans");

  tsci->ReadFile(infile,"thr/I:counts[64]/F:tor/F:qor/F");
  
  tsci->Print();

  TCanvas * cc = new TCanvas("cc","Threshold Scan Out");
  cc->Divide(4,4);
  cc->Update();

  TString sss;
  TH1F *hdummy;

  TGraph *gmean = new TGraph(64);
  gmean->SetTitle("ThrScan mean");
  gmean->SetMarkerStyle(20);
  
  TGraph *grms = new TGraph(64);
  grms->SetTitle("ThrScan RMS");
  grms->SetMarkerStyle(20);
  
  for (int i=0;i<64;i++) {
    cc->cd(i / 4 + 1); //->SetLogy(); // for channels / plot
    if ((i%4)==0) { // prepare same plot windows for all channels
      tsci->SetMarkerStyle(1);
      tsci->SetMarkerColor(0);
      tsci->Draw("counts:thr", "counts>0");
    }
    
    int j = i%4;
    tsci->SetMarkerStyle(20+j);
    tsci->SetMarkerColor(1+j);
    tsci->SetLineColor(1+j);
    tsci->Draw(Form("counts[%d]:thr",i),"","lp,same");
    float mean = TMath::Mean(tsci->GetSelectedRows(), tsci->GetV2(), tsci->GetV1());
    float rms = TMath::RMS(tsci->GetSelectedRows(), tsci->GetV2(), tsci->GetV1());

    gmean->SetPoint(i,(float) i, mean);
    grms->SetPoint(i,(float) i, rms);
    
    if ((i%4) == 3) {
      hdummy = (TH1F*) gPad->GetPrimitive("htemp");
      hdummy->SetTitle(Form("channels: %d to %d",i-3,i));
    }
    cc->Update();

  }

  TCanvas *c2 = new TCanvas("c2",Form("ThrScan mean and RMS (file %s)",infile.Data()));
  c2->Divide(2,2);
  c2->cd(1);
  gmean->Draw("PAW");
  c2->cd(3);
  grms->Draw("PAW");
  c2->Update();
  
  return 0;
};

/*
 * Spectroscopy Pedestals
 *
 */
TGraph * getSpePedestal(TChain *ctlist) { //from root file
  
  TCanvas *cc1 = new TCanvas("c1ped","Pedestal Charge distributions");

  TGraph *gpedestal =  plotActiveSChannels(ctlist, cc1, "hgain");
  
  gpedestal->SetNameTitle("Pedestals Mean vs RMS; Mean; RMS");
  gpedestal->SetMarkerStyle(20);

  return gpedestal;
}

/*
 * Read and plot PHA HG or LG files
 * and estimate pedestals
 * return the graph of pedestal-mean and -rms
 */

TGraph* readPHArun(int run, TString inpath, int board=0, TString sgain="HG") {

  TTree *tdummy;

  TGraph *gr[64];
  float xmin,xmax,ymin,ymax;
  
  xmin=ymin = 1e99;
  xmax=ymax=-1;
  int nbin=0;

  TGraph *gmean = new TGraph(64);
  gmean->SetTitle("PHA mean");
  gmean->SetMarkerStyle(20);
  
  TGraph *grms = new TGraph(64);
  grms->SetTitle("PHA RMS");
  grms->SetMarkerStyle(20);

  TGraph *gpedestal = new TGraph(64);
  gpedestal->SetNameTitle("Pedestals Mean vs RMS; Mean; RMS");
  gpedestal->SetMarkerStyle(20);
  
  for (int i=0;i<64;i++) {
    
    tdummy = new TTree("tdummy","dummy container");

    TString infile = Form("%s/Run%d_PHA_%s_%d_%d.txt",inpath.Data(),run,sgain.Data(),board,i);

    if (tdummy->ReadFile(infile,"count/I") == 0) { continue; } // nothing to read
    
    tdummy->Draw("count:Entry$","count>0","goff");
    
    float mean=0;
    float rms=1e99;

    int n = tdummy->GetSelectedRows();

    if (n>0) {
      mean = TMath::Mean(n, tdummy->GetV2(), tdummy->GetV1());
      rms = TMath::RMS(n, tdummy->GetV2(), tdummy->GetV1());
      gr[i] = new TGraph(n, tdummy->GetV2(), tdummy->GetV1());

      float fdu = TMath::MinElement(n,tdummy->GetV2());
      xmin = (xmin>fdu) ? fdu : xmin;

      fdu = TMath::MinElement(n,tdummy->GetV1());
      ymin = (ymin>fdu) ? fdu : ymin;

      fdu = TMath::MaxElement(n,tdummy->GetV2());
      xmax = (xmax<fdu) ? fdu : xmax;

      fdu = TMath::MaxElement(n,tdummy->GetV1());
      ymax = (ymax<fdu) ? fdu : ymax;
    }
    
    gmean->SetPoint(i,(float) i, mean);
    grms->SetPoint(i,(float) i, rms);

    gpedestal->SetPoint(i, mean, rms);
    
    nbin = n;

  }

  if (nbin==0) { return NULL; }
  
  printf("bins %d Min/Max : x %f %f  y %f %f\n",nbin,xmin,xmax, ymin,ymax);

  TH2F *hdummy = new TH2F("hd","PHA distribution",nbin,xmin,xmax, 10, ymin,ymax);

  TCanvas * cc = new TCanvas("cc",Form("PHA distribution (run %d)",run));
  cc->Divide(4,4);
  cc->Update();

  for (int i=0;i<64;i++) {
    cc->cd(i / 4 + 1); //->SetLogy(); // for channels / plot
    if ((i%4)==0) { // prepare same plot windows for all channels
      hdummy->SetTitle(Form("PHA channels: %d to %d",i,i+3));
      hdummy->DrawCopy();
    }

    int j = i%4;
    gr[i]->SetMarkerStyle(20+j);
    gr[i]->SetMarkerColor(1+j);
    gr[i]->SetLineColor(1+j);
    gr[i]->Draw("lp"); 
    
    cc->Update();

  }

  TCanvas *c2 = new TCanvas("c2",Form("PHA mean ad RMS (run %d)",run));
  c2->Divide(2,2);
  c2->cd(1);
  gmean->Draw("PAW");
  c2->cd(4);
  grms->Draw("PAW");
  c2->cd(3);
  gpedestal->Draw("PAW");
  c2->Update();
  
  return gpedestal;
};


/*
 * Read and plot Staircase txt file
 */

int readStaircaseRun(int run, TString inpath, int board=0) {

  TTree *tdummy;

  TGraph *gr[64];

  float xmin,xmax,ymin,ymax;

  xmin=ymin = 1e99;
  xmax=ymax=-1;
  int nbin=0;

  for (int i=0;i<64;i++) {
    
    tdummy = new TTree("tdummy","dummy container");

    TString infile = Form("%s/Run%d_Staircase_%d_%d.txt",inpath.Data(),run,board,i);

    tdummy->ReadFile(infile,"thr/I:count/F");
  
    tdummy->Draw("count:thr","","goff");

    int n = tdummy->GetSelectedRows();

    gr[i] = new TGraph(n, tdummy->GetV2(), tdummy->GetV1());

    float fdu = TMath::MinElement(n,tdummy->GetV2());
    xmin = (xmin>fdu) ? fdu : xmin;

    fdu = TMath::MinElement(n,tdummy->GetV1());
    ymin = (ymin>fdu) ? fdu : ymin;

    fdu = TMath::MaxElement(n,tdummy->GetV2());
    xmax = (xmax<fdu) ? fdu : xmax;

    fdu = TMath::MaxElement(n,tdummy->GetV1());
    ymax = (ymax<fdu) ? fdu : ymax;

    nbin = n;

  }

  printf("bins %d Min/Max : x %f %f  y %f %f\n",nbin,xmin,xmax, ymin,ymax);

  TH2F *hdummy = new TH2F("hd","Thr Staircase",nbin,xmin,xmax, 10, ymin,ymax);

  TCanvas * cc = new TCanvas("cc",Form("Thr Staircase (run %d)",run));
  cc->Divide(4,4);
  cc->Update();


  for (int i=0;i<64;i++) {
    cc->cd(i / 4 + 1); //->SetLogy(); // for channels / plot
    if ((i%4)==0) { // prepare same plot windows for all channels
      hdummy->SetTitle(Form("THRscan channels: %d to %d",i,i+3));
      hdummy->DrawCopy();
    }

    int j = i%4;
    gr[i]->SetMarkerStyle(20+j);
    gr[i]->SetMarkerColor(1+j);
    gr[i]->SetLineColor(1+j);
    gr[i]->Draw("lp"); 
    
    cc->Update();

  }
  
  return 0;
};

/*
 */

int testIt(TString sin="adogr 123    rtdsw     aesed  addfd") {

  TObjArray *a = sin.Tokenize(" ");

  for (int i=0; i<a->GetEntries(); i++) {
    TObjString *obs = (TObjString *) a->At(i);
    TString ss = obs->GetString();
    //ss.ReplaceAll("\n","").ReplaceAll("\t","").ReplaceAll(" ","");
    printf("%s\n",ss.Data());  
  }
  
  return 0;
}

/*
 * Read Spectroscopy Event List file in txt (ascii) format
 * pedestal runs can be: PHA data (single file for the moment) or list file(s)
 */

int readSEventList(TString basepath, // path of the run files
		   vector<int> srun, // signal (beam) run numbers
		   vector<int> prun={-1}, // pedestal run number, -1 for no pedestal, apply default threshold
		   float nsigma=3, // number of sigma for the pedestal RMS (or its default value)
		   int thr_gain=4096,  // max gain value (13 bits adc values, should not be requires) - fix potential bug in ASCII readout Spectroscopy mode
		   int channel=-1, // if >=0 process channel and related paired on the opposite end of the bar
		   int nboard=0) { // number of electronic board

  setLocalStyle();

  //TString inpath = basepath+"/fers/";
  TString inpath = basepath; //tolto fers perchè sembra cerca fers fers
  
  TString parTitle="[Runs:";
  for (uint i=0;i<srun.size();i++) {
    parTitle += Form(" %d",srun[i]);
  }
  parTitle+=" (";
  for (uint i=0;i<prun.size();i++) {
    parTitle += Form(" %d",prun[i]);
  }
  parTitle+=Form(" ) nsig %.1f thr_g %d]",nsigma,thr_gain);
  
  printf("##### Parse text data if needed\n");
  srun.insert(srun.end(), prun.begin(), prun.end()); // for simplicity parse all runs in the same way
  
  checkAndConvertText2Root(inpath, srun);
  
  srun.erase(srun.end()-prun.size(),srun.end()); // remove pedestal runs from list of signal runs
  
  int nchannel=64;

  int pedrun=prun[0]; // TBD
  // read pedestal, if any
  TChain *ctlist;
  TGraph *gped = readPHArun(pedrun, inpath, nboard); // try as PHA ...
  if (gped == NULL) {
    if (prun[0]>=0) { // try to load pedestals from tree(s)
      ctlist = new TChain("tlist");
      for (uint i=0;i<prun.size();i++) {
	TString ifile = inpath + Form("/Run%d.root",prun[i]);
	ctlist->Add(ifile);
      }
      gped = getSpePedestal(ctlist); // try root file
      ctlist->Delete();
    }
    if (gped == NULL) { // no file provided or bad files
      gped = new TGraph(64);
      for (int i=0;i<64;i++) {
	gped->SetPoint(i,0,0); // no pedestal
      }
    }
  }
  
  ctlist = new TChain("tlist");
  for (uint i=0;i<srun.size();i++) {
    TString ifile = inpath + Form("/Run%d.root",srun[i]);
    ctlist->Add(ifile);
  }


  if ((channel>=0)&&(c2side(channel)>=0)) { // process single bar
    singleBarSpe(ctlist, channel,4096);
    return 0;
  }
  
  TCanvas *css = new TCanvas("css",Form("HGain Charge Distr. %s",parTitle.Data()));
  plotActiveSChannels(ctlist, css, "hgain:lgain", Form("(hgain<%d)&&(lgain<%d)",thr_gain,thr_gain));
    
  int counttrg = ctlist->GetEntries();

  TGraph *ghmean = new TGraph(64);
  ghmean->SetTitle("Gain Mean - Pedestal Subtracted %s");
  ghmean->SetMarkerStyle(20);
  
  TGraph *ghrms = new TGraph(64);
  ghrms->SetTitle("Gain RMS - Pedestal Subtracted");
  ghrms->SetMarkerStyle(20);

  TGraph *glmean = new TGraph(64);
  glmean->SetTitle("Low Gain Mean");
  glmean->SetMarkerStyle(22);
  glmean->SetMarkerColor(2);
    
  TGraph *glrms = new TGraph(64);
  glrms->SetTitle("Low Gain RMS");
  glrms->SetMarkerStyle(22);
  glrms->SetMarkerColor(2);
  
  TCanvas *c1 = new TCanvas("c1",Form("Charge Stats %s",parTitle.Data()));
  c1->Divide(1,2);
  c1->Update();

  TCanvas *c1a = new TCanvas("c1a",Form("Gain Stats %s",parTitle.Data()));
  c1a->Divide(2,2);
  c1a->Update();

  TH2F *h2hit = new TH2F("h2hit","Hit map",64,-0.5,63.5,64,-0.5,63.5);

  TCanvas *cnv0 = new TCanvas("cnv0",Form("Centroids %s",parTitle.Data()));
  cnv0->Divide(2,2);
  cnv0->Update();

  cnv0->cd(1);
  ctlist->Draw("xc[0]:xc[2]","qt[0]+qt[2]","colz");
  cnv0->cd(2);
  ctlist->Draw("yc[1]:yc[3]","qt[1]+qt[3]","colz");
  cnv0->cd(3);
  ctlist->Draw("qt[1]+qt[3]:qt[0]+qt[2]","","colz");
  cnv0->cd(4);
  ctlist->Draw("(yc[1]*qt[1]+yc[3]*qt[3])/(qt[1]+qt[3]):(xc[0]*qt[0]+xc[2]*qt[2])/(qt[0]+qt[2])",
	      "((qt[1]+qt[3])>0)&&((qt[0]+qt[2])>0)","colz");
  cnv0->Update();
  
  // process-present data
  int cntu = 0; // count unmasked channels
  for (int i=0;i<64;i++) { // loop on the 64 channels

    if (c2side(i)==-1) { continue; } // channel masked

    double pmean, prms;
    gped->GetPoint(i, pmean, prms);
    
    if (isnan(prms)) {
      printf(" Warning ch %d has rms = NaN, set to 1e99\n",i);
      prms =1e99;
    }
    double sthr = pmean + nsigma*prms;
    
    //    printf(" channel: %d mean: %f +/- %f ->thr: %f\n",i,pmean,prms,sthr);

    c1->cd(1);
    ctlist->Draw("hgain",Form("(hgain>%f)&&(chan==%d)",sthr,i));
    if (ctlist->GetSelectedRows()<=0) { continue; }
    
    TH1F *hf = (TH1F *) gPad->GetPrimitive("htemp");
    c1->Update();
    
    float mean = hf->GetMean();
    float rms = hf->GetRMS();

    //    printf("  high gain mean / rms : %f %f\n",mean,rms);
    
    ghmean->SetPoint(cntu, (float) i, mean);
    ghrms->SetPoint(cntu, (float) i, rms);
    cntu+=1;

    /* not used for the moment
       tlist->Draw(Form("lgain[0][%d]",i),Form("lgain[0][%d]>=0",i),"");
       hf = (TH1F *) gPad->GetPrimitive("htemp");

       if (hf != NULL) { // to be understood

       mean = hf->GetMean();
       rms = hf->GetRMS();
       
       } else {
       mean = 0;
       rms = 0;
       }
       glmean->SetPoint(i, (float) i, mean);
       glrms->SetPoint(i, (float) i, rms);
       printf(" low gain mean / rms : %f %f\n",mean,rms);
    */
  }

  ghmean->Set(cntu);
  ghrms->Set(cntu);

  //  TLegend *leg = new TLegend(0.7,0.65,0.9,0.9);
  //  leg->AddEntry(ghmean,"High Gain");
  // leg->AddEntry(glmean,"Low Gain");
  
  c1->cd(1);
  ghmean->Draw("paw");
  //  glmean->Draw("p");
  // leg->Draw();
  c1->cd(2);
  ghrms->Draw("paw");
  //  glrms->Draw("p");
  
  c1->Update();

  ctlist->SetMarkerStyle(7);
  c1a->cd(1);
  ctlist->Draw("hgain:lgain");
  c1a->cd(2);
  ctlist->Draw("hgain:lgain",Form("(hgain>=%d)||(lgain>=%d)",thr_gain,thr_gain));
  c1a->cd(3);
  ctlist->Draw("hgain:lgain",Form("(hgain<%d)&&(lgain<%d)",thr_gain,thr_gain));
  c1a->cd(4);
  ctlist->Draw("chan",Form("(hgain>=%d)||(lgain>=%d)",thr_gain,thr_gain));
  c1a->Update();

  
  TCanvas *c2 = new TCanvas("c2","versus time");
  c2->Update();
  c2->cd(1);
  ctlist->SetMarkerStyle(47);
  // tlist->SetEstimate(-1);
  ctlist->Draw(Form("Min$(hgain):Max$(hgain)"),Form("hgain>=0"),"goff");

  float ymin = TMath::MinElement(counttrg, ctlist->GetV1()); 
  float ymax = TMath::MaxElement(counttrg, ctlist->GetV2()); 
  ctlist->Draw(Form("timeus/1000."),"","goff");
  float xmax = TMath::MaxElement(ctlist->GetSelectedRows(), ctlist->GetV1()); 
  printf(" max time: %f, gain min max: %f %f\n",xmax, ymin, ymax);

  /*
    TProfile *hf = (TProfile *) gPad->GetPrimitive("htemp");
    hf->SetXTitle("Trigger Time [ms]");
    hf->SetYTitle("High Gain ADC");
    printf(" max/min : %f %f %f\n",hf->GetYmax(), hf->GetYmin(), ymin);
  */

  ctlist->Draw("hgain:timeus/1000.","hgain>=0","prof");

  int ntrig = ctlist->GetSelectedRows();
  float tmin = TMath::MinElement(ntrig, ctlist->GetV2());
  float tmax = TMath::MaxElement(ntrig, ctlist->GetV2());

  float trate = ((float) ntrig)/(tmax-tmin)*1000.; //trigger/sec
  printf(" Time interval: %.2f %.2f ms triggers: %d rate: %f\n",tmin,tmax,ntrig,trate);
						     
  int ach[8]={0,1,4,5,35,36,41,42}; // to be improved
  for (int i=0;i<8;i++) {
    ctlist->SetMarkerStyle(20+i);
    ctlist->SetMarkerColor(2+i);
    ctlist->Draw("hgain:timeus/1000.",Form("(hgain>=0)&&(chan==%d)",ach[i]),"prof,same");
  }

  TH2F *hf2 = (TH2F *) gPad->GetPrimitive("htemp");
  hf2->SetXTitle("Trigger Time [ms]");
  hf2->SetYTitle("High Gain ADC");
  hf2->SetTitle("ADC values vs Time");

  hf2->SetMinimum(ymin);
  hf2->SetMaximum(ymax);
  
  c2->Update();
  
  return 0;
}

/*
 * Plot relevant (hopefully) data and produce x/y hit map
 */

summaryList* plotTimeData(TChain *chainT) {

  // process and display data stored in TTree

  TCut cbepu="bPulse<6e12"; // low beam pulse cut

  // list of "summary" parameters

  summaryList *spl = new summaryList("plotTimeData");

  //  TParameter<float> *pl[20];
  // for (int i=0;i<20;i++) { pl[i] = NULL; }
  
  chainT->LoadTree(0);
  chainT->GetTree()->GetUserInfo()->Print();
  //  int run = daqp->runNumber();
  TList *userlist = (TList *) chainT->GetTree()->GetUserInfo();
  TParameter<int> *par = (TParameter<int>*) userlist->At(0); // first par is run number
  int run=par->GetVal();
  printf("Plot data from run %d\n",run);
  float timebw = ((TParameter<float>*) userlist->FindObject("TimeBinWidth"))->GetVal();
  std::time_t start_time = ((TParameter<std::time_t>*) userlist->FindObject("StartTime"))->GetVal();

  std::cout << "START TIME : " << start_time << std::endl;
  
  TCanvas *cc1 = new TCanvas("cc1",Form("Timing mode correlation (run %d)",run));
  cc1->Divide(1,1);
  cc1->Update();
  cc1->cd(1);
  chainT->Draw(Form("tot*%f:toa*%f",timebw,timebw),"","colz");
  TH2F* hdummy2 = (TH2F*) gPad->GetPrimitive("htemp");
  hdummy2->SetTitle(Form("Timing / Amplitude correlation (run %d)",run));
  hdummy2->SetXTitle("Time of Arrival [ns]");
  hdummy2->SetYTitle("Time over Threshold [ns]"); 
  cc1->Update();

  TCanvas *cc2 = new TCanvas("cc2","Ref Trigger Time Distribution");
  //  cc2->Divide(2,2);
  //  cc2->Update();
  TPad *pad1 = new TPad("pad1", "Pad1", 0.05,0.51,0.98,0.97);
  TPad *pad2 = new TPad("pad2", "Pad2", 0.05,0.02,0.98,0.49);
  pad1->Draw();
  pad2->Draw();
  pad2->cd();
  TPad *pad21 = new TPad("pad21", "Pad21", 0.02,0.05,0.48,0.98);
  TPad *pad22 = new TPad("pad22", "Pad22", 0.52,0.05,0.98,0.98);
  pad21->Draw();
  pad22->Draw();
  
  // cc2->cd(1);
  pad1->cd();
  chainT->Draw("novt:timeus/1000.","","prof");
  int nn = chainT->GetSelectedRows();

  float smean=-1;
  float srms=-1;
  if (nn>0) {
    //    pl[0] = new TParameter<float>("Signal_mean",TMath::Mean(nn,chainT->GetV1())); // signal over thr time mean ~ charge mean
    //    pl[1] = new TParameter<float>("Signal_rms", TMath::StdDev(nn,chainT->GetV1())); // signal over thr. time rms ~ charge rms
    smean= TMath::Mean(nn,chainT->GetV1());
    srms = TMath::StdDev(nn,chainT->GetV1());
  }
  spl->addFloat(smean, "Signal_mean");
  spl->addFloat(srms, "Signal_rms");
  
  float tstart = TMath::MinElement(nn, chainT->GetV2());
  float tstop = TMath::MaxElement(nn, chainT->GetV2());
  float acqtime = (tstop-tstart)/1000.; // [s]
  float scale = 1./acqtime;

  printf(" number of loaded events: %d  acqtime %f [ms] -> event rate %f [Hz]\n",nn,acqtime, ((float) nn)/acqtime);

  hdummy2 = (TH2F*) gPad->GetPrimitive("htemp");
  hdummy2->SetTitle(Form("OverThreshold Signals vs Trigger Time since start of run %d",run));
  hdummy2->SetXTitle("Trigger Time [ms]");
  hdummy2->SetYTitle("Number OverThr"); 
  
  float min=1e99;
  float max=-1e99;
  
  for (int i=1;i<nn;i++) {
    float delta = chainT->GetV2()[i]-chainT->GetV2()[i-1];
    min = (delta < min) ? delta : min;
    max = (delta > max) ? delta : max;
  }
  float dmm = (max-min)/199.;

  printf(" delta time min/max %f %f [ms]\n",min,max);
  
  TH1F *hdt = new TH1F("hdt","Time Distance between triggers",200,min-dmm,max+dmm);
  hdt->SetXTitle("Delta Time [ms]"); 

  for (int i=1;i<nn;i++) {
    float delta = chainT->GetV2()[i]-chainT->GetV2()[i-1];
    hdt->Fill(delta);
  }
  //  cc2->cd(2)->SetLogy();
  pad21->cd()->SetLogy();
  hdt->Draw();

  pad22->cd();
  chainT->Draw("novt"); // number of overtreshold channels / event
  TH1F *hdummy = (TH1F*) gPad->GetPrimitive("htemp");
  hdummy->SetTitle(Form("Acquired number of channels/event (run %d)",run));
  hdummy->SetXTitle("NumChannels/Event"); 
  TF1 *fpe = new TF1("fpe","pol2(0)+expo(3)",2,62); // fit to find minimum
  fpe->SetRange(2,62);
  hdummy->Fit(fpe,"R");
  float ovt_thr = fpe->GetMinimumX(2,60);
  float ovt_low = hdummy->Integral(0,ovt_thr);
  spl->addFloat(ovt_thr, "Count_thr");

  //    pl[2] = new TParameter<float>("Count_Thr",ovt_thr); // per event

  //  pl[3] = new TParameter<float>("Ovt_count_Low", ovt_low);
  //  pl[4] = new TParameter<float>("Ovt_count_High", hdummy->Integral(0,64)-ovt_low);
  cc2->Update();

  TCanvas *cc3 = new TCanvas("cc3","Compare to beam pulses");
  TPad *cc3p1 = new TPad("cc3p1", "Pad31", 0.05,0.51,0.98,0.97);
  TPad *cc3p2 = new TPad("cc3p2", "Pad32", 0.05,0.02,0.98,0.49);
  cc3p1->Draw();
  cc3p2->Draw();
  cc3p1->cd();
  TPad *cc3p11 = new TPad("cc3p11", "Pad11", 0.02,0.05,0.48,0.98);
  TPad *cc3p12 = new TPad("cc3p12", "Pad312", 0.52,0.05,0.98,0.98);
  cc3p11->Draw();
  cc3p12->Draw();

  cc3p11->cd();
  chainT->SetMarkerStyle(1);
  chainT->Draw("novt:bPulse");
  chainT->SetMarkerStyle(7);
  chainT->SetMarkerColor(2);
  chainT->Draw("novt:bPulse",cbepu, "same");
  float mean_novt0 = TMath::Mean(chainT->GetSelectedRows(), chainT->GetV1());
  chainT->SetMarkerColor(1);
  chainT->Draw("novt:bPulse",!cbepu, "same");
  float mean_novt1 = TMath::Mean(chainT->GetSelectedRows(), chainT->GetV1());

  spl->addFloat(mean_novt0, "Mean_novt_low_beam");
  spl->addFloat(mean_novt1, "Mean_novt_high_beam");

  cc3p12->cd();
  chainT->Draw("bDTime");
  //chainT->Draw(Form("tot*%f:bPulse",timebw),Form("%f",scale),"colz");
  chainT->Draw(Form("toa*%f:bPulse",timebw),Form("%f",scale),"colz");
  
  cc3p2->cd();
  chainT->Draw("novt:timeus/1000000");
  chainT->SetMarkerColor(2);
  chainT->Draw("novt:timeus/1000000","bPulse<6e12","same");
  chainT->SetMarkerColor(1);
  cc3->Update();
  
  chainT->Draw("novt",Form("novt<=%f",ovt_thr),"goff");
  int ndu = chainT->GetSelectedRows();
  //  pl[3] = new TParameter<float>("NSignal_Low", ((float) ndu)); // total signals below threshold
  spl->addInt(ndu, "NSignal_Low");

  smean=srms=-1;
  if (ndu>0) {
    //    pl[4] = new TParameter<float>("Signal_mean_Low", TMath::Mean(ndu,chainT->GetV1())); // signal/event
    //    pl[5] = new TParameter<float>("Signal_rms_Low", TMath::StdDev(ndu,chainT->GetV1()));
    smean = TMath::Mean(ndu,chainT->GetV1());
    srms = TMath::StdDev(ndu,chainT->GetV1());
  }
  spl->addFloat(smean, "Signal_mean_Low");
  spl->addFloat(srms, "Signal_rms_Low");

  chainT->Draw("novt",Form("novt>%f",ovt_thr),"goff");
  ndu = chainT->GetSelectedRows();
  //pl[6] = new TParameter<float>("NSignal_high", ((float) ndu));
  spl->addInt(ndu, "NSignal_high");

  smean=srms=-1;
  if (ndu>0) {
    //    pl[7] = new TParameter<float>("Signal_mean_High", TMath::Mean(ndu,chainT->GetV1()));
    //    pl[8] = new TParameter<float>("Signal_rms_High", TMath::StdDev(ndu,chainT->GetV1()));
    smean = TMath::Mean(ndu,chainT->GetV1());
    srms = TMath::StdDev(ndu,chainT->GetV1());
  }
  spl->addFloat(smean, "Signal_mean_High");
  spl->addFloat(srms, "Signal_rms_High");
  
  printf(" scale %f\n",scale);
  TCanvas *cc = new TCanvas("cc",Form("Timing mode (run %d)",run));
  cc->Divide(1,2);
  cc->Update();
  
  cc->cd(1);
  chainT->Draw(Form("tot*%f",timebw),Form("%f",scale));

  hdummy = (TH1F*) gPad->GetPrimitive("htemp");
  hdummy->SetTitle(Form("Signal Time Lenght - Overthreshold (run %d)",run));
  hdummy->SetXTitle("Time Over Threshold [ns]");

  //  pl[9] = new TParameter<float>("ToT_mean",hdummy->GetMean());
  //  pl[10] = new TParameter<float>("ToT_rms",hdummy->GetRMS());
  spl->addFloat(hdummy->GetMean(), "ToT_mean");
  spl->addFloat(hdummy->GetRMS(), "ToT_rms");
    
  cc->cd(2)->SetLogy();
  chainT->Draw(Form("toa*%f",timebw),Form("%f",scale));
  float reftwin = TMath::MaxElement(chainT->GetSelectedRows(), chainT->GetV1());
  //  float reft0 = TMath::MinElement(chainT->GetSelectedRows(), chainT->GetV1());
  printf(" Ref Time Window : %f\n",reftwin);
  hdummy = (TH1F*) gPad->GetPrimitive("htemp");

  // pl[11] = new TParameter<float>("ToA_mean",hdummy->GetMean());
  // pl[12] = new TParameter<float>("ToA_rms",hdummy->GetRMS());
  spl->addFloat(hdummy->GetMean(), "ToA_mean");
  spl->addFloat(hdummy->GetRMS(), "ToA_rms");
  
  hdummy->SetTitle(Form("Signal Arrival Time run %d",run));
  hdummy->SetXTitle("Time of Arrival [ns]"); 
  cc->Update();

  //  pl[13] = new TParameter<float>("AcqTime_s",acqtime);
  //  pl[14] = new TParameter<float>("Events",((float) chainT->GetEntries()));
  spl->addFloat(acqtime, "pAcqTime_s");
  spl->addInt(chainT->GetEntries(), "pEvents");
  /*  
      TList *rpl = new TList();
      for (int i=0;i<20;i++) { if (pl[i]) { rpl->Add(pl[i]); } }
      rpl->Print();
  */
  
  return spl;
  
};

/* = = = = = = = = = = = = =
 *
 */

summaryList* processTEvents(TString basepath="data",
			    vector<int> srun={15}, // signal runs, require "{}" brackets, if <0 look at pedestal only
			    vector<int> prun={-1}, // pedestal runs, can be cumulated (if negative not considered)
			    float totsig=6.6, // default Time Over Threshold sigma  
			    float nsigma=3.,
			    float timegate=100.,
			    TString ntoFile="cube/run_pkup_sall.root")
{

  TString inpath = basepath + "/fers/";

  setLocalStyle();

  summaryList *srList=NULL; // list with run summary parameters
  
  printf("##### Check and parse text data if needed\n");
  srun.insert(srun.end(), prun.begin(), prun.end()); // for simplicity parse all runs (signal/pedestal) in the same way

  checkAndConvertText2Root(basepath, srun, "/fers/Run", ntoFile);
  
  printf("##### Plot signal data and extract active channels\n");
  srun.erase(srun.end()-prun.size(),srun.end()); // remove pedestal runs from list of signal runs

  for(uint i=0;i<srun.size();i++) { printf(" %d\n",srun[i]); }

  TChain *ctlist;
  
  if (srun[0]!=-1) { 
    ctlist = new TChain("tlist");

    for (uint i=0;i<srun.size();i++) {
      TString ifile = inpath + Form("/Run%d.root",srun[i]);
      ctlist->Add(ifile);
    }
    
    srList = plotTimeData(ctlist);
    
  } else { ctlist = NULL; }  // look at the pedestal only
  
  vector<int> ach = getActiveChannels(ctlist); // this requires the signal ttree and not the pedestal ttree

  ctlist->Delete();

  // get pedestals
  
  printf("##### Get pedestal\n");
  vector<float> *ped;
  
  if (prun[0]>=0) { // load pedestals
    ctlist = new TChain("tlist");
    for (uint i=0;i<prun.size();i++) {
      TString ifile = inpath+"/"+ Form("/Run%d.root",prun[i]);
      ctlist->Add(ifile);
    }
      // TFile *fped = new TFile(ifile,"read");
      // TTree *tped = (TTree *) fped->Get("tlist");
    ped = getToTPedestal(ctlist, ach, totsig);
    ctlist->Delete();
  } else {
    ped = getToTPedestal(NULL, ach, totsig);
  }

  // now reopen the signal chain
  printf("##### Now process signal data\n");

  vector<float> agd;
  summaryList *spList=NULL;  
  if (srun[0]!=-1) { 
    ctlist = new TChain("tlist");
    for (uint i=0;i<srun.size();i++) {
      TString ifile = inpath+"/"+ Form("/Run%d.root",srun[i]);
      printf(" %s\n",ifile.Data());
      ctlist->Add(ifile);
    }  

    if (ctlist != NULL) {
      spList = channelTimeProcess(ctlist, timegate, nsigma, ped);
      printf("\n%d ",srun[0]);
      for (uint i=0;i<agd.size();i++) {
	printf("%f ",agd[i]);
      }
      printf("\n");
    }
  }

  // append the two list objects
  if (srList != NULL) {
    spList->merge(srList);
    delete srList;
    /*
      for (int i=0;i<srList->GetSize();i++) {
      TParameter<float> *pp;
      pp = (TParameter<float> *) srList->At(i);
      agd.push_back(pp->GetVal());
      }
    */
  }
  return spList;
  
};

/*
 * analyze multiple runs (adapted to 2310 CERN test)
 */

void runTimingMode(TString ipath="data") {
  
  vector<vector<int>> run = {{15,31},{16},{17,27}, // at least one for each delay and other run variables
			     {18},{19},{20},
			     {21,26},{28},{29},{30},
			     {44},{45},{46},
			     {48},{49,53},{55},{54},
			     {56},{57},{58}};
  
  vector<int> delay={0,500,1000,1500,2000,2500,3000,3500,4000,10000,  // [ns] delay from cube trigger signal
		     500, 2000, 10000,
		     10000,2000,500,0,
		     2000,2000,2000};

  vector<float> tdgap={20,20,20,20,20,20,20,20,20,20,
		       10,10,10,
		       30,30,30,30,
		       20,20,20}; // [cm] target detector distance

  vector<int> target={0,0,0,0,0,0,0,0,0,0, // target type 0: light 0.3mm, 1: light-vertical, 2:medium-vertical 0.5mm, 3:heavy-vertical 0.5mm+reinforced dome
		      0,0,0,
		      0,0,0,0,
		      1,2,3};
  
  vector<int> prun={13,22,33}; // pedestal runs

  vector<float> ctime={50,100,200,350,500}; // coincidence gate window time

  //  vector<vector<float>> vdb;
  vector<string> vdb;
  TString stitle="";
  
  uint nruns = run.size();

  summaryList *sld;

  for (uint i=0;i<nruns; i++) {
    printf("- - - - first signal run %d\n",run[i][0]);
    for (uint j=0;j<ctime.size(); j++) {
      sld = processTEvents(ipath,run[i],prun,6.6,3.,ctime[j]);
      sld->addFloat(ctime[j],"CoincGate_par");
      sld->addFloat(tdgap[i],"TargDetGap_par");
      sld->addInt(delay[i],"gammaDelay_par");
      sld->addInt(target[i],"Target_par");
      sld->addInt(run[i][0],"Run_par");
      vdb.push_back((sld->getValues()).data());
      if (stitle=="") {
	stitle = sld->getNames(1);
      }
      delete sld;
      /*
	v.insert(v.begin(),ctime[j]);
	v.insert(v.begin(),tdgap[i]);
	v.insert(v.begin(),(float) delay[i]);
	vdb.push_back(v);
      */
    }
  }

  printf("----- Summary\n%s\n",stitle.Data());

  for (uint i=0;i<vdb.size();i++) {
    printf("%s\n",vdb[i].data());
  }
  
  /*
  for (uint i=0;i<vdb.size();i++) {
    for (uint j=0;j<vdb[i].size(); j++) {
      printf("%.2f ",vdb[i][j]);
    }
    printf("\n");
  }
  */

};

/*
 * analyse multiple runs in spectroscopy mode
 *
 * NOTE: before running this macro, use parseListFile(...) to convert the Run*_list.txt file into root file
 *
 */

void runSpectMode(int ori=0, // bar orientation/axis: 0:vertical/x or 1:horizontal/y
		  int pos=5, // position along axis (0 to 7)
		  int thr_gain=4096, // threshold on max ADC-charge (to avoid saturated/overflow data) 
		  float nsigma=3., // number of sigma that define the threshold of the signal, relative to pedestal 
		  TString ipath="/home/cisbani/proj/ntof/x17-gustavino/test/2310_cern/data/fers") {
  
  gStyle->SetTitleW(0.8);
  gStyle->SetTitleH(0.1);
  gStyle->SetTitleSize(0.7,"XYZT");
  gStyle->SetTitleXOffset(.7);
  gStyle->SetTitleYOffset(.7);
  gStyle->SetLabelSize(0.05,"XYZ");
  gStyle->SetPadTopMargin(.1);
  gStyle->SetPadLeftMargin(.06);
  gStyle->SetPadRightMargin(.0);
  gStyle->SetPadBottomMargin(.07);
  gStyle->SetOptStat(0);
  gStyle->SetOptFit(0);
  
  TString sori[2]={"V","H"};
  
  int srun[4]={38,40,41,42}; // spect mode cern oct 2023
  int gdelay[4]={0,500,2000,10000};
  int prun=37; // pedestal run

  // get fers/SiPM channels corresponding to orientation and position
  vector<int> chs;
  TString scut=""; // remove outlier
  int cnt=0;
  for (int i=0;i<64;i++) {
    if ((c2side(i)>=0)&&(c2axis(i) == ori)&&(c2idx(i) == pos)) {
      chs.push_back(i);
      printf(" channels %d selected\n",i);
      if (cnt==0) {
	scut = Form("(chan==%d)",i);
      } else {
	scut = scut + Form("||(chan==%d)",i);
      }
      cnt += 1;
    }
  }
  
  // pedestal
  TChain *ctlist = new TChain("tlist");
  TString ifile = ipath + Form("/Run%d.root",prun);
  ctlist->Add(ifile);
  TGraph *gped = getSpePedestal(ctlist); // try root file
  ctlist->Delete();

  TCanvas *cc = new TCanvas("cc","Spectroscopy channels...");
  cc->Divide(4,2);
  cc->Update();

  // signal runs
  for (int i=0;i<4;i++) { // loop on runs 
    int run = srun[i];
    TChain *ctlist = new TChain("tlist");
    TString ifile = ipath+"/"+ Form("/Run%d.root",run);
    ctlist->Add(ifile);

    // get acqtime in s
    ctlist->Draw(Form("timeus/1000000."),"","goff");
    float tacq = TMath::MaxElement(ctlist->GetSelectedRows(), ctlist->GetV1()); 
    printf(" AcqTime %f s\n",tacq);

    float weight = 1./tacq;
    
    TCut ccut = Form("(%f)*((hgain<%d)&&(%s))",weight,thr_gain,scut.Data());
    ccut.Print();
    
    TString sbeam[2]={"(bPulse>=6e12)","(bPulse<6e12)"};
    
    for (uint j=0;j<chs.size();j++) { // loop on bar side channels

      ctlist->SetMarkerStyle(1);
      ctlist->SetMarkerColor(1);
      ctlist->SetLineColor(1);

      cc->cd(i+1+j*4)->SetLogy();
      ctlist->Draw("hgain",ccut,"p");
      TH1F *hdummy = (TH1F*) gPad->GetPrimitive("htemp");
      hdummy->SetTitle(Form("%s-%d #gamma delay %d ns",sori[ori].Data(),pos,gdelay[i]));
    
      for (int k=0;k<2;k++) { // loop on main / parassitic pulses
      
	double pmean,prms;
	gped->GetPoint(chs[j], pmean, prms);
      
	double sthr = pmean + nsigma*prms;

	ctlist->SetMarkerStyle(22+k);
	ctlist->SetMarkerColor(2+k);
	ctlist->SetLineColor(2+k);
	TCut call  = Form("(%f)*((hgain>%f)&&(chan==%d)&&(%s))",weight,sthr,chs[j],sbeam[k].Data());
	call.Print();

	ctlist->Draw("hgain",call,"pl,same");
	
      }
    }
    cc->Update();
  }
  
}