HarpoVertexing.cxx

Go to the documentation of this file.

// 
// File HarpoVertexing.cxx
//
#include "HarpoVertexing.h"
#include "HarpoConfig.h"
#include "HarpoDetSet.h"
#include "HarpoDebug.h"
#include "HarpoDccEvent.h"
#include "Pmm2Event.h"
#include "HarpoEvent.h"
#include "HarpoRecoEvent.h"
#include "HarpoSimEvent.h"
#include "MakeNiceHisto.h"

#include "TFile.h"
#include "TStyle.h"
#include "TCanvas.h"
#include "TRootEmbeddedCanvas.h"
#include "TLatex.h"
#include "TGraph.h"
#include "TF1.h"
#include "TLine.h"
#include "TMath.h"
#include "TGLabel.h"
#include "TEllipse.h"
#include "TSpectrum.h"
#include "TLinearFitter.h"
#include "TPolyMarker.h"
#include "TSystem.h"

#include "RVersion.h"

#include <cstdlib>
#include <cstring>
#include <cassert>
#include <fstream>
#include <iostream>

ClassImp(HarpoVertexing);



// typedef struct vPoint_t {
//   Double_t x;
//   Double_t y;
//   Double_t w;
// } vPoint_t;






void   HarpoVertexing::print()
 {
   unsigned int nd; // number of detectors
   HarpoEventHeader *hdr = fEvt->GetHeader();

   assert(hdr != NULL);
   hdr->print();
   
   for (nd = 0; nd < gkNDetectors; nd++) {
     //     if (fEvt->isdataExist(nd)) {
       HarpoDccMap *plane = fEvt->GetDccMap(nd);
       if (plane != NULL )plane->print();
   }
 }

void   HarpoVertexing::process()
{
  //  Bool_t drawEvent = kFALSE;
  nEvents++;
  if(!fEvt->GetRecoEvent()) return;


  for(Int_t i = 0; i<20; i++){
    fHistProj2D[i]->Reset();
    fHistProj2Draw[i]->Reset();
    fHistMean[i]->Reset();
    fHistMap[i]->Reset();
  }
  
  // Find vertex candidates
  FindVertexes();
  //RefineVertexPosition();

  // Compute track directions associated to the vertex
  FillVertexAngles();




}


void HarpoVertexing::FillVertexAngles()
{

  // Compute track directions associated to the vertices found by FindVertexes()
  // Completes the HarpoRecoVertex objects in fEvt
  
  for(Int_t i = 0; i<2; i++){
    for(Int_t j = 0; j<10; j++){
      for(Int_t k = 0; k<3; k++){
        fHist[i][j][k] = 0;
        fHist2[i][j][k] = 0;
      }
    }
  }
  
  TClonesArray* vArray = fEvt->GetRecoEvent()->GetVertexArray();
  Int_t nV2 = vArray->GetEntries();
  HarpoDccMap* m[2];
  for(Int_t plane = 0; plane<2;plane++)
    m[plane] = fEvt->GetDccMap(plane);

  // Loop over the vertex candidates found by FindVertexes()
  for(Int_t iV = 0; iV<nV2; iV++){
    HarpoRecoVertex* vertex = (HarpoRecoVertex*)vArray->At(iV);
    Int_t plane = vertex->GetPlane();
    Double_t xV = vertex->GetXvertex();
    Double_t zV = vertex->GetZvertex();
    // Double_t pxV = vertex->GetPxVertex();
    // Double_t pzV = vertex->GetPzVertex();
    // Double_t thV = vertex->GetThetaVertex();

    // Fill a 2D histogram in polar coordinates
    TH2F* hProj2DV_1 = (TH2F*)fHistProj2D[iV]->Clone("hProj2D_1");
    // Use only the data connex with the vertex candidate
    GetConnex(Int_t(zV),Int_t(xV),3,3,m[plane],fHistMap[iV]);
    for(Int_t i=0;i<NADC;i++){ //Time bins
      for(Int_t j=0;j<NCHAN;j++){ // Channels
        Double_t q = fHistMap[iV]->GetBinContent(i+1,j+1);
        if(q<=0) continue;
        Double_t xTmp = j - xV, zTmp = i - zV;
        Double_t rho = TMath::Sqrt(xTmp*xTmp+zTmp*zTmp);
        Double_t theta = TMath::ATan2(xTmp,zTmp);
        hProj2DV_1->Fill(rho,theta,q);
      }
    }
    // Correct the charge to account for the surface change
    for(Int_t bin = 1; bin<=fNbins; bin++){
      Int_t k = Int_t(0.5/TMath::Pi()*fNbinsTheta/hProj2DV_1->GetXaxis()->GetBinCenter(bin));
      for(Int_t i = 0; i<fNbinsTheta; i++){
        Double_t q = hProj2DV_1->GetBinContent(bin,i+1)/(2*k+1.);
        if(q<=0) continue;
        //q *= bin+2;
        for(Int_t j = -k; j<=k; j++){
          Int_t biny = i+1+j;
          // Wrap around the circle
          if(biny<1)           biny =  fNbinsTheta + biny;
          if(biny>fNbinsTheta) biny = -fNbinsTheta + biny;
          fHistProj2Draw[iV]->SetBinContent(bin,biny,fHistProj2Draw[iV]->GetBinContent(bin,biny)+q);
        }
      }
    }
    hProj2DV_1->Delete();
    //    Double_t thetaTmp[fNbinsTheta*fNbins];
    //    Double_t qTmp[fNbinsTheta*fNbins];
    Int_t nTmp = 0;
    // In the polar histogram fHistProj2Draw[iV]
    // Remove the bin which are not connected in a straight line to the vertex
    for(Int_t i = 0; i<fNbinsTheta; i++){ // For each theta bin
      Int_t nZero = 0;
      Double_t qTot = 0;
      for(Int_t bin = 2; bin<=fNbins; bin++){
        Double_t q = fHistProj2Draw[iV]->GetBinContent(bin,i+1);
        if(nZero>1){ // If at least 2 consecutive rho bins are empty
          // Erase the rest of the theta bin
          fHistProj2D[iV]->SetBinContent(bin,i+1,0);
          continue;
        }else{
          fHistProj2D[iV]->SetBinContent(bin,i+1,q);
          qTot+=q;
          //      thetaTmp[nTmp] = fHistProj2D[iV]->GetYaxis()->GetBinCenter(i+1);
          //      qTmp[nTmp] = q;
          nTmp++;
        }
        // Count number of consecutive empty rho bin
        if(q<=10) nZero++;
        else nZero = 0;
      }
      // Fill projection on the theta axis
      fHistMean[iV]->SetBinContent(fNbinsTheta/2+i+1,qTot);
      // Wrap around the circle
      if(i<fNbinsTheta/2)
        fHistMean[iV]->SetBinContent( fNbinsTheta+fNbinsTheta/2+i+1,qTot);
      else
        fHistMean[iV]->SetBinContent(-fNbinsTheta+fNbinsTheta/2+i+1,qTot);

    }
    //    Double_t mu[4000], sigma[4000];
    //    Int_t nPeak = 0;
#if 1
    // Look for peaks in the theta projection fHistMean[iV] (extended to -2pi,2pi)
    TSpectrum *s = new TSpectrum(8);
    Int_t nfound = s->Search(fHistMean[iV],2,"",0.01);
    if(gHarpoDebug>1 || fDisplay)
      Info("process","%d peaks found",nfound);
#if ROOT_VERSION_CODE >= ROOT_VERSION(6,0,0)
    Double_t *xpeaks = s->GetPositionX();
    Double_t *ypeaks = s->GetPositionY();
#else
    Float_t *xpeaks = s->GetPositionX();
    Float_t *ypeaks = s->GetPositionY();
#endif    
    Double_t peaks[nfound];
    Int_t npeaks = 0;
    Int_t* index = new Int_t[nfound];
    // Select the highest peaks
    TMath::Sort(nfound,ypeaks,index,kTRUE);
    for(Int_t i = 0; i<nfound; i++){
      // reject peaks outside (-pi,pi)
      if(xpeaks[index[i]]<-TMath::Pi()) continue;
      if(xpeaks[index[i]]>TMath::Pi()) continue;
      peaks[npeaks] = xpeaks[index[i]];
      if(gHarpoDebug>1)
        Info("process","%d   %g",npeaks,peaks[npeaks]);
      npeaks++;
    }
    delete index;
    s->Delete();
    nfound = npeaks;
    for(Int_t i = 0; i<nfound; i++){
      for(Int_t j = i+1; j<nfound; j++){
        if(TMath::Abs(peaks[i]-peaks[j])>2 &&
           TMath::Abs(peaks[i]-peaks[j]+TMath::Pi()*2)>2 &&
           TMath::Abs(peaks[i]-peaks[j]-TMath::Pi()*2)>2){
          peaks[j] = -1000;
        }
        if(i==0) continue;
        // merge extra peaks if they are close together
        if(TMath::Abs(peaks[i]-peaks[j])<0.05)
          peaks[i] = (peaks[i]+peaks[j])/2.;
        if(TMath::Abs(peaks[i]-peaks[j]+TMath::Pi()*2)<0.05)
          peaks[i] = (peaks[i]+peaks[j])/2.+TMath::Pi();
        if(TMath::Abs(peaks[i]-peaks[j]-TMath::Pi()*2)<0.05)
          peaks[i] = (peaks[i]+peaks[j])/2.+TMath::Pi();
      }
    }
    npeaks = 0;
    for(Int_t i = 0; i<nfound; i++){
      if(peaks[i]<-10) continue;
      peaks[npeaks] = peaks[i];
      npeaks++;
    }
    // Keep only the two highest peaks, inside -pi,pi
    if(npeaks>2) npeaks = 2;
#else
    // Failed attempt to make a better peak search
    //   TArrayD* arr = GetPeaks(nTmp,thetaTmp,qTmp,plane,iV);
    TArrayD* arr = GetPeaks2(fHistMean[iV]);
    Int_t npeaks = 0;
    Double_t peaks[20], ypeaks[20];
    for(Int_t i = 0; i<10; i++){
      if(arr->At(3*i) == -1000) break;
      npeaks++;
      peaks[i] = arr->At(3*i);
      ypeaks[i] = arr->At(3*i+2)/10;
    }
    TPolyMarker * pm =
      (TPolyMarker*)fHistMean[iV]->GetListOfFunctions()->FindObject("TPolyMarker");
    if (pm) {
      fHistMean[iV]->GetListOfFunctions()->Remove(pm);
      delete pm;
    }
    pm = new TPolyMarker(npeaks, peaks, ypeaks);
    fHistMean[iV]->GetListOfFunctions()->Add(pm);
    pm->SetMarkerStyle(23);
    pm->SetMarkerColor(kRed);
    pm->SetMarkerSize(1.3);
    //    if(npeaks>2) npeaks = 2;
#endif
    if(gHarpoDebug>1 || fDisplay)
      Info("process","%d peaks confirmed",npeaks);
    

    vertex->SetChi2Vertex(npeaks);

    if(npeaks==1){ // Single track
      vertex->SetPxVertex(TMath::Sin(peaks[0]));
      vertex->SetPzVertex(TMath::Cos(peaks[0]));
      vertex->SetThetaVertex(0);
    }
    if(npeaks==2){ // Pair vertex
      vertex->SetPxVertex(TMath::Sin(0.5*(peaks[0]+peaks[1])));
      vertex->SetPzVertex(TMath::Cos(0.5*(peaks[0]+peaks[1])));
      Double_t theta = TMath::Abs(peaks[0]-peaks[1]);
      if(theta>2*TMath::Pi()) theta -= 2*TMath::Pi();
      if(theta>TMath::Pi()) theta = 2*TMath::Pi() - theta;
      vertex->SetThetaVertex(theta);
    }
  }
  return;
  
}


void HarpoVertexing::RefineVertexPosition()
{

  TClonesArray* vArray = fEvt->GetRecoEvent()->GetVertexArray();
  Int_t nV2 = vArray->GetEntries();
  HarpoDccMap* m[2];
  for(Int_t plane = 0; plane<2;plane++)
    m[plane] = fEvt->GetDccMap(plane);
  
  for(Int_t iV = 0; iV<nV2; iV++){
    HarpoRecoVertex* vertex = (HarpoRecoVertex*)vArray->At(iV);
    Int_t plane = vertex->GetPlane();
    Double_t xV = vertex->GetXvertex();
    Double_t zV = vertex->GetZvertex();
    // Double_t pxV = vertex->GetPxVertex();
    // Double_t pzV = vertex->GetPzVertex();
    //    Double_t thV = vertex->GetThetaVertex();

    Double_t meanX = 0, meanZ = 0, qTot = 0;
    for(Int_t ch = Int_t(xV - fRadiusX); ch<Int_t(xV + fRadiusX); ch++){
      for(Int_t tb = Int_t(zV - fRadiusZ); tb<Int_t(zV + fRadiusZ); tb++){
        if(ch<0 || ch>=NCHAN || tb<0 || tb>=NADC) continue;
        Double_t q = m[plane]->GetData(ch,tb);
        if(q<=0) continue;
        meanX += ch*q;
        meanZ += tb*q;
        qTot += q;
      }
    }
    if(qTot<=0) continue;

    vertex->SetXvertex(meanX/qTot);
    vertex->SetZvertex(meanZ/qTot);

  }

}





void HarpoVertexing::FindVertexes()
{
  // Find vertex candidates
  // Reconstructed data (clusters, tracks, matched tracks)
  HarpoRecoEvent* reco = fEvt->GetRecoEvent();

#if 1
  // Clusters in the event (both X and Y directions)
  TClonesArray* clArray = reco->GetClustersArray();
  Int_t nCl = clArray->GetEntries();

  Double_t x[2][nCl], z[2][nCl], px[2][nCl], pz[2][nCl], theta[2][nCl];
  Int_t nV[2];
  nV[0] = 0;
  nV[1] = 0;
  for(Int_t i = 0; i<nCl; i++){ // Loop over clusters
    HarpoRecoClusters* cl = (HarpoRecoClusters*)clArray->At(i);
    Double_t q = cl->GetQ();
    if(q<fQmin) continue;
    Int_t plane = cl->GetPlane();
    Double_t xTmp = cl->GetX();
    Double_t zTmp = cl->GetZ();
    Double_t phi0, phi1;
    Double_t frac = GetFraction(plane,xTmp,zTmp,phi0,phi1);
    if(frac>fMinFrac) continue; // reject opening angle too large
    px[plane][nV[plane]] = -TMath::Cos(0.5*(phi0 + phi1));
    pz[plane][nV[plane]] = -TMath::Sin(0.5*(phi0 + phi1));
    theta[plane][nV[plane]] = TMath::Abs(2*TMath::Pi() - phi1 + phi0);
    if(theta[plane][nV[plane]] > +2*TMath::Pi()) 
      theta[plane][nV[plane]] -= 2*TMath::Pi();
    x[plane][nV[plane]] = xTmp;
    z[plane][nV[plane]] = zTmp;
    if(gHarpoDebug>1)
      Info("FindVertex","%d %d   %g %g   %g %g   %g",plane,nV[plane],zTmp, xTmp, pz[plane][nV[plane]], px[plane][nV[plane]], theta[plane][nV[plane]]);
    nV[plane]++;
  }

  for(Int_t plane = 0; plane<2; plane++){ // Remove multiple occurences of same vertex
    Int_t used[nV[plane]];
    for(Int_t i = 0; i<nV[plane]; i++)
      used[i] = 0;
    for(Int_t i = 0; i<nV[plane]; i++){
      if(used[i]) continue;
      used[i] = 1;
      Double_t chi2 = -10;
      Int_t jmin = i;
      //      Double_t thetaMin = 10;
      for(Int_t j = 0; j<nV[plane]; j++){
        if(used[j]) continue;
        if(TMath::Abs(x[plane][i]-x[plane][j])<fRadiusX && TMath::Abs(z[plane][i]-z[plane][j])<fRadiusZ){
          used[j] = 1;
          if(TMath::Abs(theta[plane][j])<TMath::Abs(theta[plane][i])){
            jmin = j;
          }
        }
      }
      HarpoRecoVertex* vertex = new HarpoRecoVertex(plane,z[plane][jmin], x[plane][jmin], pz[plane][jmin], px[plane][jmin], theta[plane][jmin], chi2, -1, -1);
      reco->AddVertex(vertex);
      vertex->Delete();
    }  
  }
#else


  HarpoRecoVertex* vertexX = new HarpoRecoVertex(0,reco->GetTstart()+8, reco->GetXstart(), 1, 0, 0, 0, -1, -1);
  reco->AddVertex(vertexX);
  HarpoRecoVertex* vertexY = new HarpoRecoVertex(1,reco->GetTstart()+8, reco->GetYstart(), 1, 0, 0, 0, -1, -1);
  reco->AddVertex(vertexY);


#endif


}



void HarpoVertexing::FitVertex(TH2F* h, HarpoRecoVertex* vertex)
{
  // Unused

  //  Int_t plane = vertex->GetPlane();
  Double_t xV = vertex->GetXvertex();
  Double_t zV = vertex->GetZvertex();
  Double_t pxV = vertex->GetPxVertex();
  Double_t pzV = vertex->GetPzVertex();
  //  Double_t thV = vertex->GetThetaVertex();
  
  TGraph* g = GetGraph(h,3,xV,zV,pxV,pzV);

  TArrayD* par = FitGraph(g,2);
  delete par;
}


TGraph* HarpoVertexing::GetGraph(TH2F* h, Int_t n, Double_t x0, Double_t z0, Double_t px, Double_t pz)
{
  // Unused

  Int_t nbinsx = h->GetYaxis()->GetNbins();
  Int_t nbinsz = h->GetXaxis()->GetNbins();

  TGraph* g = new TGraph();

  Int_t nP = 0;
  for(Int_t i = 0; i<nbinsz; i+=n){
    for(Int_t j = 0; j<nbinsx; j+=n){
      Double_t x = 0, z = 0, q = 0;
      for(Int_t i2 = 0; i2<n; i2++){
        for(Int_t j2 = 0; j2<n; j2++){
          Double_t val = h->GetBinContent(i+i2+1,j+j2+1);
          z += (i+i2)*val;
          x += (j+j2)*val;
          q += val;
        }
      }
      if(q<=0) continue;
      x /= q;
      z /= q;

      Double_t dx = (x-x0)/fRadiusX;
      Double_t dz = (z-z0)/fRadiusZ;
      if(dx*dx + dz*dz<0.5) continue;

      //      g->SetPoint(nP,(x-x0)*px + (z-z0)*pz,-(z-z0)*px+(x-x0)*pz,q);
      g->SetPoint(nP,(x-x0)*px + (z-z0)*pz,-(z-z0)*px+(x-x0)*pz);
      nP++;
    }
  }

  return g;
}


TArrayD* HarpoVertexing::FitGraph(TGraph* g, const Int_t nFit)
{
  // Unused

  Int_t nP = g->GetN();
  Double_t* x = g->GetX();
  Double_t* y = g->GetY();
  //  Double_t* q = g->GetZ();

  // Double_t w[nFit][nP];
  Double_t p[nFit];
  for(Int_t j = 0; j<nFit; j++)
    p[j] = 1./nFit;
  TArrayD* par = new TArrayD(2*nFit);
  if(nFit == 1){
    par->AddAt(0,0);
    par->AddAt(0,1);
  }else{
    for(Int_t j = 0; j<nFit; j++){
      par->AddAt(0,2*j);
      par->AddAt(-1 + 2*j*1./(nFit-1),2*j+1);
    }
  }

  for(Int_t k = 0; k<5; k++){
    
    Double_t w[nFit][nP];
    for(Int_t i = 0; i<nP; i++){
      Double_t norm = 0;
      for(Int_t j = 0; j<nFit; j++){
        Double_t d = y[i] - x[i]*par->At(2*j+1) - par->At(2*j);
        w[j][i] = p[j]*TMath::Gaus(d,0,1);
        norm += w[j][i];
      }
      for(Int_t j = 0; j<nFit; j++)
        w[j][i] /= norm;
    }
    
    for(Int_t j=0; j<nFit; j++){
      
      Double_t xTmp[nP], yTmp[nP];
      Int_t nTmp = 0;
      for(Int_t i = 0; i<nP; i++){
        Bool_t test = kTRUE;
        for(Int_t J = 0; J<nFit; J++){
          if(w[J][i]>w[j][i]){
            test = kFALSE;
            break;
          }
        }
        if(!test) continue;
        xTmp[nTmp] = x[i];
        yTmp[nTmp] = y[i];
        nTmp++;
      }
      
      if(nTmp<2){
        Info("FitGraph","Too many fits %d",nFit);
        return 0;
      }
      
      p[j] = nTmp*1./nP;
      
      TLinearFitter* fFitCh = new TLinearFitter(1,"pol1","D");
      fFitCh->AssignData(nTmp, 1, xTmp, yTmp);
      fFitCh->Eval();
      
      par->AddAt(fFitCh->GetParameter(0),2*j);
      par->AddAt(fFitCh->GetParameter(1),2*j+1);
      
      Info("FitGraph","Iteration %d, Fit %d, chi2 = %g",k,j,fFitCh->GetChisquare());
      
      fFitCh->Delete();
    }
    
  }
  
  return par;
}


TArrayD* HarpoVertexing::GetPeaks2(TH1F* h)
{
  // Unused

  TArrayD* arr = new TArrayD(30);
  for(Int_t k = 0; k<30; k++)
    arr->AddAt(-1000,k);
  Int_t nbins = h->GetXaxis()->GetNbins();
  Int_t up = 0;
  Double_t noise = 1000;
  Double_t th = 0, q = 0, sig = 0, qold = 0;
  Int_t nCl = 0;
  Int_t istart = 1;
  for(Int_t i = 1; i<=nbins; i++){
    Double_t val = h->GetBinContent(i);
    if(val<noise){
      istart = i;
      break;
    }
  }
  Info("GetPeaks2","istart = %d",istart);

  for(Int_t i = istart; i<=nbins+istart; i++){
    Double_t val, x;
    if(i<=nbins){
      val = h->GetBinContent(i);
      x = h->GetXaxis()->GetBinCenter(i);
    }else{
      val = h->GetBinContent(i-nbins);
      x = h->GetXaxis()->GetBinCenter(i-nbins);
    }
    if(up == 0 && val<noise) continue;
    if((up == 2 && val>qold+noise) || val<100){ // end of cluster
      if(q>0001){
        //      if(gHarpoDebug>1 || fDisplay)
        Info("GetPeaks2","End cluster %g: mu = %g, sig = %g, q = %g, qold = %g",x,th/q,TMath::Sqrt(sig/q-th/q*th/q),q,qold);
        arr->AddAt(th/q,3*nCl);
        arr->AddAt(TMath::Sqrt(sig/q-th/q*th/q),3*nCl+1);
        arr->AddAt(q,3*nCl+2);
        nCl++;
      }
      th = 0; q = 0; sig = 0;
    }
    if(val>qold+noise){
      if(up == 0) Info("GetPeaks2","Start cluster %g (%g %g)",x,val, qold);
      up = 1;
    }
    if(val<qold-noise){ // peak
      if(up == 1) Info("GetPeaks2","Peak %g (%g %g)",x,val, qold);
      up = 2;
    }
    th += x*val;
    sig += x*x*val;
    q += val;
    if(val<noise){
      th = 0; q = 0; sig = 0; up = 0;
    }
    if(TMath::Abs(val-qold)>noise)
      qold = val;
  }

  return arr;
}


TArrayD* HarpoVertexing::GetPeaks(Int_t n, Double_t *x, Double_t *w, Int_t plane, Int_t v)
{
  // Unused
  Int_t nLeftOld = 0;
  TArrayD* arr = new TArrayD(30);
  for(Int_t k = 0; k<30; k++)
    arr->AddAt(-1000,k);
  Int_t i = 0;
  Int_t loop = 0;
  while(1){
    Double_t mu, sig;
    GetPeak(n,x,w,mu,sig,plane,v);
    Int_t nLeft = 0;
    Double_t q = 0;
    Double_t qTot = 0, qTot2 = 0;
    for(Int_t j = 0; j<n; j++){
      qTot += w[j];
      Double_t dist = TMath::Abs(x[j]-mu);
      if(TMath::Abs(x[j]-mu-2*TMath::Pi())<dist) dist = TMath::Abs(x[j]-mu-2*TMath::Pi());
      if(TMath::Abs(x[j]-mu+2*TMath::Pi())<dist) dist = TMath::Abs(x[j]-mu+2*TMath::Pi());
      Double_t frac = dist < 3*sig;
      if(frac<0.01) frac = 0;
      //else Info("","frac = %g",frac);
      q += w[j]*frac;
      w[j] *= 1.-frac;
      qTot2 += w[j];
      if(w[j]<10) nLeft++;
    }
    if(gHarpoDebug>1 || fDisplay){
      Info("GetPeaks","loop %d, qTot = %g, qTot2 = %g, q = %g",loop,qTot,qTot2, q);
      Info("GetPeaks","mu = %g, sig = %g, q = %g, nLeft = %d",mu,sig,q,nLeft);
    }
    if(sig<0.2 && q > 1000){
      arr->AddAt(mu,i);
      i++;
      arr->AddAt(sig,i);
      i++;
      arr->AddAt(q,i);
      i++;
      if(i>=20) break;
    }
    if(nLeft<10) break;
    if(nLeft==nLeftOld) break;
    nLeftOld = nLeft;
    loop++;
  }

  return arr;

}

void HarpoVertexing::GetPeak(Int_t nElem, Double_t *x, Double_t *w, Double_t &muOut, Double_t &sigOut, Int_t plane, Int_t v)
{
  // Unused

  muOut = -100;
  sigOut = -100;

  Double_t xCenter = 0.;
  Double_t xWidth  = 6.5;
  Int_t fNbinr = 100;
  Int_t nloop = 0, maxLoops = 5;
  while(nloop < maxLoops){
    if(gHarpoDebug>1) 
      Info("FindPairNew","loop == %d",nloop);

    TH1F* h = new TH1F("h","",fNbinr,xCenter - xWidth, xCenter + xWidth);
    TH1F* hB = new TH1F("hB","",fNbinr,xCenter - xWidth, xCenter + xWidth);
    TH1F* hW = new TH1F("hW","",fNbinr,xCenter - xWidth, xCenter + xWidth);
    TH1F* h1 = new TH1F("h1","",fNbinr,xCenter - xWidth, xCenter + xWidth);
    TH1F* h2 = new TH1F("h2","",fNbinr,xCenter - xWidth, xCenter + xWidth);
    TH1F* h3 = new TH1F("h3","",fNbinr,xCenter - xWidth, xCenter + xWidth);
    TH1F* h4 = new TH1F("h4","",fNbinr,xCenter - xWidth, xCenter + xWidth);

    for(Int_t i=0;i<nElem;i++){
      if(x[i]<-10) continue;
      if(w[i]<10) continue;
      Double_t theta = x[i];
      Double_t weight= w[i];
      h->Fill(theta);
      hW->Fill(theta, weight);
      h1->Fill(theta, theta*weight);
      h2->Fill(theta, theta*theta*weight);
      h3->Fill(theta, theta*theta*theta*weight);
      h4->Fill(theta, theta*theta*theta*theta*weight);
      
      if(theta>0) theta = theta - 2*TMath::Pi();
      else theta = theta + 2*TMath::Pi();
      h->Fill(theta);
      hW->Fill(theta, weight);
      h1->Fill(theta, theta*weight);
      h2->Fill(theta, theta*theta*weight);
      h3->Fill(theta, theta*theta*theta*weight);
      h4->Fill(theta, theta*theta*theta*theta*weight);
    }
   
    if(h->GetEntries()<3) {
      h->Delete();
      hB->Delete();
      hW->Delete();
      h1->Delete();
      h2->Delete();
      h3->Delete();
      h4->Delete();
      break;
    }
    //    Double_t sRmax = 0, sTmax = 0;

    if(nloop<3){
      //      Info("","%d_%d_%d",plane,v,nloop);
      fHist[plane][v][nloop] = (TH1F*)h->Clone(Form("h%d_%d_%d",plane,v,nloop));
    }
    // Calculate Mean and RMS in histograms
    for (Int_t ir=1;ir<=fNbinr;ir++){
      Double_t n = h->GetBinContent(ir);
      Double_t weight = hW->GetBinContent(ir);
      if(n<=0) continue;
      if(weight<=0) continue;
      Double_t m1Tmp = h1->GetBinContent(ir)/weight;
      Double_t m2Tmp = h2->GetBinContent(ir)/weight;
      Double_t m3Tmp = h3->GetBinContent(ir)/weight;
      Double_t m4Tmp = h4->GetBinContent(ir)/weight;
      m4Tmp = m4Tmp - 4*m1Tmp*m3Tmp + 6*m1Tmp*m1Tmp*m2Tmp - 3*m1Tmp*m1Tmp*m1Tmp*m1Tmp;
      m2Tmp = m2Tmp - m1Tmp*m1Tmp;
      //      Double_t var4 = 1./n*(m4Tmp - (n-3)*m2Tmp*m2Tmp/(n-1));
      //Info("GetPeak","%d  %g  %g %g  (%d)",ir,n,m1Tmp,m2Tmp,nElem);
      // if(m2Tmp>0)
      //        hB->SetBinContent(ir,n*h->GetXaxis()->GetBinWidth(1)/m2Tmp);
      // else
        hB->SetBinContent(ir,n);

    }


    if(nloop<3){
      fHist2[plane][v][nloop] = (TH1F*)h->Clone(Form("hB%d_%d_%d",plane,v,nloop));
    }

    hB->GetXaxis()->SetRangeUser(-TMath::Pi(),TMath::Pi());
    Int_t imax = hB->GetMaximumBin();
    Double_t max = hB->GetXaxis()->GetBinLowEdge(imax);    


    Double_t muPeak = h1->GetBinContent(imax)/hW->GetBinContent(imax);
    Double_t sigmaPeak = TMath::Sqrt(h2->GetBinContent(imax)/hW->GetBinContent(imax) - muPeak*muPeak);
    //Info("GetPeak","%d imax = %d, thmax = %g, max=%g, mu = %g, sigma = %g, w = %g",nloop,imax,max,h->GetBinContent(imax),muPeak,sigmaPeak,h->GetXaxis()->GetBinWidth(1));


    Double_t fSigmaMaxRebin = 0.25, fCutoffHough = 6, kfFacHough = 3;
    Int_t fMinPairs = 3;
    if(sigmaPeak<fSigmaMaxRebin*h->GetXaxis()->GetBinWidth(1) && nloop < maxLoops - 1 && h->GetBinContent(imax)>fMinPairs){
      xCenter = max;
      xWidth /= kfFacHough;
      nloop++;
      h->Delete();
      hB->Delete();
      hW->Delete();
      h1->Delete();
      h2->Delete();
      h3->Delete();
      h4->Delete();
      continue;
    }

    Int_t ioMax = h->GetXaxis()->GetNbins()/2;
    Double_t rapr0 = 0;
    for(Int_t io=1;io<=ioMax;io++){
      Double_t n = hW->Integral(imax-io,imax+io);
      Double_t mu=h1->Integral(imax-io,imax+io)/n;
      Double_t sig=h2->Integral(imax-io,imax+io)/n - mu*mu;
          
      if(sig>0) sig = TMath::Sqrt(sig);
      else sig = 0;
      //Info("GetPeak","n = %g, mu = %g, sigma = %g",n,mu,sig);
    
      Double_t rapr=0;
      if(sig>0) rapr = (io) * h->GetXaxis()->GetBinWidth(1) / sig;
      //Info("GetPeak","%g = %d * %g / %g",rapr,io-1,h->GetXaxis()->GetBinWidth(1),sig);

      if(io==ioMax || (rapr>fCutoffHough && rapr0<fCutoffHough) ){
        muOut = mu;
        if(muOut>TMath::Pi()) muOut -= 2*TMath::Pi();
        if(muOut<-TMath::Pi()) muOut += 2*TMath::Pi();
        sigOut = sig;
        h->Delete();
        hB->Delete();
        hW->Delete();
        h1->Delete();
        h2->Delete();
        h3->Delete();
        h4->Delete();
        return;
      }
      rapr0 = rapr;
    }  
    h->Delete();
    hB->Delete();
    hW->Delete();
    h1->Delete();
    h2->Delete();
    h3->Delete();
    h4->Delete();
  }
}

Double_t HarpoVertexing::GetFraction(Int_t ndet, Double_t x, Double_t z, Double_t &phi0, Double_t &phi1)
{

  // Look for the longest section of an ellipse that does not encounter a region with signal
  // (x,z) is the centre of the ellipse, (fRadiusX,fRadiusZ) their radius
  //
  
  HarpoDccMap* m = fEvt->GetDccMap(ndet);
  Int_t n = 0;
  Int_t l = 0, start = 1000; // start is initialised at 1000 to avoid problems when there is no signal at all along the ellipse
  Bool_t empty = kFALSE;
  Double_t phiStart = 0;
  phi0 = -1; phi1 = -2;
  for(Int_t i=0;i<1000;i++){
    // Loop over points along the ellipse (twice around)
    Double_t phi = i*4*TMath::Pi()/1000;
    Double_t xTmp = x + fRadiusX*TMath::Cos(phi);
    Double_t zTmp = z + fRadiusZ*TMath::Sin(phi);
    // corresponding pixel
    Int_t channel = Int_t(xTmp);
    Int_t timebin = Int_t(zTmp);
    if(channel<0 || channel>=NCHAN) continue;
    if(timebin<0 || timebin>=NADC) continue;
    Int_t q = m->GetData(channel,timebin);
    // Add the neighbouring positive signals
    if(channel>0) q+= TMath::Max(m->GetData(channel-1,timebin),0.);
    if(channel<NCHAN-1) q+= TMath::Max(m->GetData(channel+1,timebin),0.);
    if(timebin>0) q+= TMath::Max(m->GetData(channel,timebin-1),0.);
    if(timebin<NADC-1) q+= TMath::Max(m->GetData(channel,timebin+1),0.);
    if(q>0){
      n++;
      empty = kFALSE;
    }else{
      if(!empty){
        phiStart = phi;
        start = i;
      }
      if(i-start>l){
        l = i-start;
        phi0 = phiStart;
        phi1 = phi;
      }
      empty = kTRUE;
    }
  }

  // The full length of the circle is 1000/2 (we go twice around)
  // l*2./1000 is the longest fraction that does not meet any signal
  // We return the converse, corresponding to an approximate bend of the track in this region
  return 1 - l*2./1000;  
}



void HarpoVertexing::GetConnex(Int_t i, Int_t j, Int_t n, Int_t n0, HarpoDccMap* m, TH2F* h)
{

  // Create a 2D histogram containing the connex region connected to bin i,j
  // Recursive function
  
  if(i<0 || j<0 || i>=NADC || j>=NCHAN) return;
  if(h->GetBinContent(i+1,j+1)>0) return;

  Double_t q = m->GetData(j,i);
  if(q==-2000){
    q = 1;
  }
  if(q<=0){
    //   q=0;
    if(n>0){
      //Info("GetConnex","Gap %d   %d %d   %g",n,i,j,q);
      n--;
      q = 1;
    }
    else return;
  }else n = n0;

  h->SetBinContent(i+1,j+1,q);
  GetConnex(i-1,j  ,n,n0,m,h);
  GetConnex(i+1,j  ,n,n0,m,h);
  GetConnex(i  ,j-1,n,n0,m,h);
  GetConnex(i  ,j+1,n,n0,m,h);

  return;
}


void HarpoVertexing::DisplayAnalysis(TRootEmbeddedCanvas* ecTab, TGListBox* /* infobox */)
{

  TCanvas* c = ecTab->GetCanvas();
  c->GetPad(1)->Delete();
  c->GetPad(2)->Delete();
  c->Divide(2);



  TClonesArray* vArray = fEvt->GetRecoEvent()->GetVertexArray();
  Int_t nV = vArray->GetEntries();
  HarpoDccMap* m[2];
  for(Int_t plane = 0; plane<2;plane++)
    m[plane] = fEvt->GetDccMap(plane);
  
  Int_t iV[2][10];
  Int_t nVp[2];
  nVp[0] = 0;
  nVp[1] = 0;
  //  TGraphErrors* gProjV[2][nV];
  TH1F* hProjV[2][nV];
  TH2F* hProj2DV[2][nV];
  TH2F* hMap2DV[2][nV];
  TH1F* hMeanV[2][nV];
  TH1F* hMean2V[2][nV];
  TH1F* hMean3V[2][nV];
  TH1F* hMean4V[2][nV];
  //  TH1D* hProf[2][nV];
  //  Int_t nbins = 200;//200;
  //  Double_t xMin = -20, xMax = 180;
  Double_t xMin = 0, xMax = 400;
  for(Int_t i = 0; i<nV; i++){
    HarpoRecoVertex* vertex = (HarpoRecoVertex*)vArray->At(i);
    Int_t plane = vertex->GetPlane();
    Double_t xV = vertex->GetXvertex();
    Double_t zV = vertex->GetZvertex();
    Double_t pxV = vertex->GetPxVertex();
    Double_t pzV = vertex->GetPzVertex();
    //    Double_t thV = vertex->GetThetaVertex();

    iV[plane][nVp[plane]] = i;

    hProjV[plane][nVp[plane]] = new TH1F(Form("hProj%d_%d",nVp[plane],plane),"",fNbins,xMin,xMax);
    //    hProj2DV[plane][nVp[plane]] = new TH2F(Form("hProj2D%d_%d",nVp[plane],plane),"",fNbins,xMin,xMax,80,-40,40);
    //    Int_t nbinsTheta = 1000;
    TH2F* hProj2DVtmp = new TH2F("hProj2Dtmp","",fNbins,xMin,xMax,fNbinsTheta,-TMath::Pi(),TMath::Pi());
    hProj2DV[plane][nVp[plane]] = new TH2F(Form("hProj2D%d_%d",nVp[plane],plane),"",fNbins,xMin,xMax,fNbinsTheta,-TMath::Pi(),TMath::Pi());
    hMap2DV[plane][nVp[plane]] = new TH2F(Form("hMap2D%d_%d",nVp[plane],plane),"",512,0,512,288,0,288);
    GetConnex(Int_t(zV),Int_t(xV),3,3,m[plane],hMap2DV[plane][nVp[plane]]);
    //    hMeanV[plane][nVp[plane]] = new TH1F(Form("hMean%d_%d",nVp[plane],plane),"",fNbins,xMin,xMax);
    hMeanV[plane][nVp[plane]] = new TH1F(Form("hMean%d_%d",nVp[plane],plane),"",fNbinsTheta,-TMath::Pi(),TMath::Pi());
    hMean2V[plane][nVp[plane]] = new TH1F(Form("hMean2%d_%d",nVp[plane],plane),"",fNbins,xMin,xMax);
    hMean3V[plane][nVp[plane]] = new TH1F(Form("hMean3%d_%d",nVp[plane],plane),"",fNbins,xMin,xMax);
    hMean4V[plane][nVp[plane]] = new TH1F(Form("hMean4%d_%d",nVp[plane],plane),"",fNbins,xMin,xMax);
    //    hProf[plane][nVp[plane]] = 0;//new TH1F(Form("hProf%d_%d",nVp[plane],plane),"",100,-50,50);
    hProjV[plane][nVp[plane]]->SetLineColor(kRed);
    hMeanV[plane][nVp[plane]]->SetLineColor(kBlack);
    hMean2V[plane][nVp[plane]]->SetLineColor(kBlue);
    hMean4V[plane][nVp[plane]]->SetLineColor(kGreen);
    for(Int_t i=0;i<NADC;i++){ //Time bins
      for(Int_t j=0;j<NCHAN;j++){ // Channels
        //      Double_t q = m[plane]->GetData(j,i);
        Double_t q = hMap2DV[plane][nVp[plane]]->GetBinContent(i+1,j+1);
        if(q<=0) continue;
        Double_t xTmp = j - xV, zTmp = i - zV;
        Double_t a = xTmp*pzV - zTmp*pxV;
        Double_t b = zTmp*pzV+xTmp*pxV;
        Double_t rho = TMath::Sqrt(xTmp*xTmp+zTmp*zTmp);
        Double_t theta = TMath::ATan2(xTmp,zTmp);
        a = theta;
        b = rho;

        hProj2DVtmp->Fill(b,a,q);
        //      hProj2DV[plane][nVp[plane]]->Fill(b,a,q);
        //if(q<=100) continue;
        hProjV[plane][nVp[plane]]->Fill(b,q);
        //hMeanV[plane][nVp[plane]]->Fill(b,q*a);
        hMean2V[plane][nVp[plane]]->Fill(b,q*a*a);
        hMean3V[plane][nVp[plane]]->Fill(b,q*a*a*a);
        hMean4V[plane][nVp[plane]]->Fill(b,q*a*a*a*a);
      }
    } 
    for(Int_t bin = 1; bin<=fNbins; bin++){
      Int_t k = Int_t(0.5/TMath::Pi()*fNbinsTheta/hProj2DV[plane][nVp[plane]]->GetXaxis()->GetBinCenter(bin));
      //if(k<1) break;
      for(Int_t i = 0; i<fNbinsTheta; i++){
        Double_t q = hProj2DVtmp->GetBinContent(bin,i+1);
        if(q<=0) continue;
        for(Int_t j = -k; j<=k; j++)
          hProj2DV[plane][nVp[plane]]->SetBinContent(bin,i+1+j,hProj2DV[plane][nVp[plane]]->GetBinContent(bin,i+1+j)+q);
      }
    }
    hProj2DVtmp->Delete();
    for(Int_t i = 0; i<fNbinsTheta; i++){
      Int_t nZero = 0;
      for(Int_t bin = 1; bin<=fNbins; bin++){
        //if(hProjV[plane][nVp[plane]]->GetXaxis()->GetBinCenter(bin)<10) continue;
        Double_t q = hProj2DV[plane][nVp[plane]]->GetBinContent(bin,i+1);
        if(nZero>1){
          hProj2DV[plane][nVp[plane]]->SetBinContent(bin,i+1,0);
          continue;
        }
        if(q<=0) nZero++;
        else nZero = 0;
      }
    }

    hProjV[plane][nVp[plane]]->SetMinimum(1);
    hProjV[plane][nVp[plane]]->GetXaxis()->SetRangeUser(-10,10);
    Double_t max = hProjV[plane][nVp[plane]]->GetMaximum();
    hProjV[plane][nVp[plane]]->GetXaxis()->SetRangeUser(xMin,xMax);
    Double_t bV;//, sigma2V = 1000;
    Bool_t foundV = kFALSE; //, foundS = kFALSE;
    for(Int_t bin = 1; bin<=fNbins; bin++){
      Double_t q = hProjV[plane][nVp[plane]]->GetBinContent(bin);
      if(q<=0) continue;

      Double_t mean = 0;//hMeanV[plane][nVp[plane]]->GetBinContent(bin)/q;
      Double_t mean2 = hMean2V[plane][nVp[plane]]->GetBinContent(bin)/q;
      Double_t sig2 = mean2-mean*mean;
      Double_t mean3 = hMean3V[plane][nVp[plane]]->GetBinContent(bin)/q;
      Double_t mean4 = hMean4V[plane][nVp[plane]]->GetBinContent(bin)/q;
      Double_t kurt = mean4 - 4*mean*mean3+6*mean*mean*mean2-3*mean*mean*mean*mean;
      kurt /= sig2*sig2;
      kurt -= 3;
      //hMeanV[plane][nVp[plane]]->SetBinContent(bin,mean);
      hMean2V[plane][nVp[plane]]->SetBinContent(bin,sig2);
      hMean4V[plane][nVp[plane]]->SetBinContent(bin,kurt);
      if(q>max*0.5 && !foundV){
        Double_t qTmp = hProjV[plane][nVp[plane]]->GetBinContent(bin-1);
        bV = q*hProjV[plane][nVp[plane]]->GetXaxis()->GetBinCenter(bin)
          + qTmp*hProjV[plane][nVp[plane]]->GetXaxis()->GetBinCenter(bin-1);
        bV /= q+qTmp;
        //      sigma2V = sig2;
        foundV = kTRUE;
      }
    }
    nVp[plane]++;


  }

  //  TH2F* hDummy = new TH2F("hDummy","",100,-100,200,100,-100,100);
  for(Int_t plane = 0; plane<2;plane++){
    TVirtualPad* cP = c->GetPad(plane+1);
    cP->Divide(3,nVp[plane]+1);
    for(Int_t i = 0; i<nVp[plane]; i++){
      // MakeNiceHisto(hMap2DV[plane][i],cP->GetPad(3*i+1),"colz");
      HarpoRecoVertex* vertex = (HarpoRecoVertex*)vArray->At(iV[plane][i]);
      Double_t xV = vertex->GetXvertex();
      Double_t zV = vertex->GetZvertex();
      MakeNiceHisto(fHistProj2Draw[iV[plane][i]],cP->GetPad(3*i+2),"colz");
      MakeNiceHisto(fHistMap[iV[plane][i]],cP->GetPad(3*i+1),"colz");
      // FIXME always 0.0 phi0, phi1
      Double_t phi0 = 0.0, phi1 = 0.0;
      //      Double_t frac = GetFraction(plane,xV,zV,phi0,phi1);
      TEllipse* eStart = new TEllipse(zV,xV,fRadiusZ,fRadiusX,-phi0*180/TMath::Pi()+90,-phi1*180/TMath::Pi()+90);
      eStart->SetLineColor(kGray);
      eStart->SetLineWidth(3);
      eStart->SetFillStyle(0);
      //      eStart->SetFillColor(kGray);
      eStart->Draw();



      Info("","Draw(fHistMean[%d])",iV[plane][i]);
      MakeNiceHisto(fHistMean[iV[plane][i]],cP->GetPad(3*i+3),"colz");
      // cP->GetPad(3*i+3)->SetLogy();
      hMeanV[plane][i]->Delete();

      //     MakeNiceHisto(fHistProj2D[iV[plane][i]],cP->GetPad(3*i+3),"colz");
    }
    if(plane == XPLANE)
      MakeNiceHisto(fHistDiff,cP->GetPad(3*nVp[plane]+1));
  }
  c->Update();

  // c->GetPad(2)->SaveAs("exampleVertex.C");
  // c->GetPad(2)->SaveAs("exampleVertex.png");
  // c->GetPad(2)->SaveAs("exampleVertex.pdf");

  return;





}



void HarpoVertexing::ConfigFrame(TGMainFrame* fMain, Int_t id)
{

  UInt_t xsize = 200;
  UInt_t ysize2 = 20;
  UInt_t ysize = 9*ysize2;
  TGTransientFrame* main = new TGTransientFrame(gClient->GetRoot(), fMain, xsize, ysize);
  main->Connect("CloseWindow()", "HarpoRecoMonitorGui", main, "CloseWindow()");
  main->DontCallClose(); // to avoid double deletions.
  
  // use hierarchical cleaning
  main->SetCleanup(kDeepCleanup);

  TGVerticalFrame* f213 = new TGVerticalFrame(main,xsize,ysize,kVerticalFrame);
  TGLabel* fLabel = new TGLabel(f213, "Vertexing");


  TGHorizontalFrame* RadiusXFrame = new TGHorizontalFrame(f213,xsize,ysize2,kHorizontalFrame);
  TGLabel* RadiusXLabel = new TGLabel(RadiusXFrame, "Radius X");
  fChooseRadiusX    = new TGNumberEntry(RadiusXFrame);
  fChooseRadiusX->SetNumber(fRadiusX);

  TGHorizontalFrame* RadiusZFrame = new TGHorizontalFrame(f213,xsize,ysize2,kHorizontalFrame);
  TGLabel* RadiusZLabel = new TGLabel(RadiusZFrame, "Radius Z");
  fChooseRadiusZ    = new TGNumberEntry(RadiusZFrame);
  fChooseRadiusZ->SetNumber(fRadiusZ);

  TGHorizontalFrame* MinFracFrame = new TGHorizontalFrame(f213,xsize,ysize2,kHorizontalFrame);
  TGLabel* MinFracLabel = new TGLabel(MinFracFrame, "Min Fraction");
  fChooseMinFrac    = new TGNumberEntry(MinFracFrame);
  fChooseMinFrac->SetNumber(fMinFrac);

  TGHorizontalFrame* SeparationFrame = new TGHorizontalFrame(f213,xsize,ysize2,kHorizontalFrame);
  TGLabel* SeparationLabel = new TGLabel(SeparationFrame, "Min Fraction");
  fChooseSeparation    = new TGNumberEntry(SeparationFrame);
  fChooseSeparation->SetNumber(fSeparation);

  TGHorizontalFrame* QminFrame = new TGHorizontalFrame(f213,xsize,ysize2,kHorizontalFrame);
  TGLabel* QminLabel = new TGLabel(QminFrame, "Min Fraction");
  fChooseQmin    = new TGNumberEntry(QminFrame);
  fChooseQmin->SetNumber(fQmin);

  // TGHorizontalFrame* DUMMYFrame = new TGHorizontalFrame(f213,xsize,ysize2,kHorizontalFrame);
  // TGLabel* DUMMYLabel = new TGLabel(DUMMYFrame, "Qx/Qy");
  // fChooseDUMMY    = new TGNumberEntry(DUMMYFrame);
  // fChooseDUMMY->SetNumber(fDUMMY);


  TGTextButton* setConf = new TGTextButton(f213,"Save Config",id);
  setConf->Associate(fMain);

  TGLayoutHints* fLayout1 = new TGLayoutHints(kLHintsLeft ,5,5,5,5);
  TGLayoutHints* fLayout2 = new TGLayoutHints(kLHintsRight ,5,5,5,5);
  TGLayoutHints* fLayout3 = new TGLayoutHints(kLHintsTop | kLHintsExpandX ,5,5,5,5);
  TGLayoutHints* fLayout4 = new TGLayoutHints(kLHintsExpandX | kLHintsExpandY,5,5,5,5);

  main->AddFrame(f213,fLayout4);
  f213->AddFrame(fLabel,fLayout3);

  f213->AddFrame(RadiusXFrame,fLayout3);
  RadiusXFrame->AddFrame(RadiusXLabel,fLayout1);
  RadiusXFrame->AddFrame(fChooseRadiusX,fLayout2);

  f213->AddFrame(RadiusZFrame,fLayout3);
  RadiusZFrame->AddFrame(RadiusZLabel,fLayout1);
  RadiusZFrame->AddFrame(fChooseRadiusZ,fLayout2);

  f213->AddFrame(MinFracFrame,fLayout3);
  MinFracFrame->AddFrame(MinFracLabel,fLayout1);
  MinFracFrame->AddFrame(fChooseMinFrac,fLayout2);

  f213->AddFrame(SeparationFrame,fLayout3);
  SeparationFrame->AddFrame(SeparationLabel,fLayout1);
  SeparationFrame->AddFrame(fChooseSeparation,fLayout2);

  f213->AddFrame(QminFrame,fLayout3);
  QminFrame->AddFrame(QminLabel,fLayout1);
  QminFrame->AddFrame(fChooseQmin,fLayout2);

  // f213->AddFrame(DUMMYFrame,fLayout3);
  // DUMMYFrame->AddFrame(DUMMYLabel,fLayout1);
  // DUMMYFrame->AddFrame(fChooseDUMMY,fLayout2);

  f213->AddFrame(setConf,fLayout3);


  main->MapSubwindows();
  main->MapWindow();
  main->Resize();

}

void HarpoVertexing::SetConfig()
{
  fDisplay = kTRUE;

  if(!fChooseRadiusX) return;

  fQmin =  fChooseQmin->GetNumber();
  fSeparation =  fChooseSeparation->GetNumber();
  fRadiusX =  fChooseRadiusX->GetNumber();
  fRadiusZ =  fChooseRadiusZ->GetNumber();
  fMinFrac =  fChooseMinFrac->GetNumber();

}
void   HarpoVertexing::Init()
{
  fChooseRadiusX = 0;
  fDisplay = kFALSE;

  fQmin = 1000;
  fSeparation = 20;
  // fRadiusX = 15;
  // fRadiusZ = 30;
  fRadiusX = 13;//10;
  fRadiusZ = 25;//20;
  fMinFrac = 0.4;

  fNbins = 200;
  Double_t xMin = 0, xMax = 400;
  fNbinsTheta = 1000;
  for(Int_t i = 0; i<20; i++){
    fHistProj2D[i] = new TH2F(Form("fHistProj2D%d",i),"",fNbins,xMin,xMax,fNbinsTheta,-TMath::Pi(),TMath::Pi());
    fHistProj2Draw[i] = new TH2F(Form("fHistProj2Draw%d",i),"",fNbins,xMin,xMax,fNbinsTheta,-TMath::Pi(),TMath::Pi());
    fHistMean[i]   = new TH1F(Form("fHistMean%d",i),"",2*fNbinsTheta,-2*TMath::Pi(),2*TMath::Pi());
    //fHistMean[i]   = new TH1F(Form("fHistMean%d",i),"",fNbinsTheta,-TMath::Pi(),TMath::Pi());
    fHistMap[i]   = new TH2F(Form("fHistMap%d",i),"",512,0,512,288,0,288);

  }


  fHistDiff = new TH1F("fHistDiff","",fNbinsTheta,-TMath::Pi()/10,TMath::Pi()/10);

 }

void   HarpoVertexing::Save(char * /* mode */)
 {


   OpenHistFile("recovertexing");

   fHistDiff->Write();

   // Save histograms here

 
 }