HarpoClusteringBlocs.cxx

Go to the documentation of this file.

// 
// File HarpoClusteringBlocs.cxx
//
#include "HarpoClusteringBlocs.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 "TMath.h"
#include "TGLabel.h"
#include "TApplication.h"
#include "TSystem.h"

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

ClassImp(HarpoClusteringBlocs);

void   HarpoClusteringBlocs::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   HarpoClusteringBlocs::process()
{
  //  Bool_t drawEvent = kFALSE;
  nEvents++;

  fEvt->GetRecoEvent()->ResetClustersArray();
  TClonesArray* clArray = fEvt->GetRecoEvent()->GetClustersArray();
  fNcl = 0;
  Int_t nClBloc[2][100];

  for(Int_t ndet = 0; ndet<2; ndet++){
    //    Info("process","---------------- PLANE %d -------------",ndet);
    for(Int_t i = 0; i<100; i++)
      nClBloc[ndet][i] = 0;
    for(Int_t i = 0; i<fSizeX; i++){
      for(Int_t j = 0; j<fSizeZ; j++){
        fMapReduced[i][j] = -1;
        fMapUsed[i][j] = -1;
      }
    }
    
    Int_t nClTot = MakeClusters(ndet);
    Int_t nBloc = 0, nClBlocs = 0;
    for(Int_t j=0;j<fSizeZ;j++){ // Channels
      for(Int_t i=0;i<fSizeX;i++){ // Time bins
        if(nBloc>100) break;
        Int_t nCl = GroupClusters(nBloc,i,j, clArray);
        nClBlocs += nCl;
        nClBloc[ndet][nBloc] = nCl;       
        if(nCl>0){
          if(gHarpoDebug>2)
            Info("process","nClBloc[%d][%d] = %d",ndet,nBloc,nClBloc[ndet][nBloc]);
          nBloc++;
        }
        if(nClBlocs >= nClTot) break;
      }
      if(nClBlocs >= nClTot) break;
    }
    if(fNminBlocs>0){
      Int_t iCl = 0;
      for(Int_t j=0;j<fSizeZ;j++){ // Channels
        for(Int_t i=0;i<fSizeX;i++){ // Time bins
          CollateClusters(iCl,i,j, clArray);
        }
      }
    }
    clArray->Compress();
  }

  Int_t nPoints = clArray->GetEntries();
  for(Int_t i = 0; i<nPoints; i++){
    HarpoRecoClusters* cluster = (HarpoRecoClusters*)clArray->At(i);
    if(cluster->GetType() != BLOCCLUSTER) continue;
    Int_t ndet = cluster->GetPlane();
    Int_t nBloc = cluster->GetIdClusterTrack();
    Double_t xSigma = cluster->GetSig();
    Double_t zSigma = cluster->GetMean();
    Double_t x = cluster->GetX();
    Double_t z = cluster->GetZ();
    Double_t q = cluster->GetQ();
    cluster->SetSig(xSigma/q - x*x);
    cluster->SetMean(zSigma/q - z*z);
    cluster->RemoveAllClusterTrack();
    if(nClBloc[ndet][nBloc]<fNminBlocs){
      clArray->RemoveAt(i);
    }
  }
  clArray->Compress();
}

Int_t HarpoClusteringBlocs::GroupClusters(Int_t nBloc, Int_t i, Int_t j, TClonesArray* clArray)
{

  if(i<0) return 0;
  if(j<0) return 0;
  if(i>=fSizeX) return 0;
  if(j>=fSizeZ) return 0;
  if(gHarpoDebug>1)
    Info("GroupClusters","%d %d   %d   %d",i,j,fMapUsed[i][j],fMapReduced[i][j]);
  if(fMapUsed[i][j] == 1) return 0;
  fMapUsed[i][j] = 1;
  Int_t nCl = fMapReduced[i][j];
  if(nCl < 0) return 0;
  if(gHarpoDebug>1)
    Info("GroupClusters","%d   %d %d   %d",nBloc,i,j,nCl);

  HarpoRecoClusters* cluster = (HarpoRecoClusters*)clArray->At(nCl);
  cluster->SetIdClusterTrack(nBloc);
  //  fMapReduced[i][j] = -nCl;
  Int_t n = 1;
  n += GroupClusters(nBloc, i+1, j, clArray);
  n += GroupClusters(nBloc, i, j+1, clArray);
  n += GroupClusters(nBloc, i-1, j, clArray);
  n += GroupClusters(nBloc, i, j-1, clArray);
 
  return n;
}


HarpoRecoClusters* HarpoClusteringBlocs::CollateClusters(Int_t iCl, Int_t i, Int_t j, TClonesArray* clArray)
{

  if(i<0) return 0;
  if(j<0) return 0;
  if(i>=fSizeX) return 0;
  if(j>=fSizeZ) return 0;
  Int_t nCl = fMapReduced[i][j];
  if(nCl < 0) return 0;

  HarpoRecoClusters* cluster = (HarpoRecoClusters*)clArray->At(nCl);
  if(gHarpoDebug>1)
    Info("CollateClusters","%d  %d %d   %d",cluster->GetPlane(),i,j,nCl);
  Double_t x = cluster->GetX();
  Double_t z = cluster->GetZ();
  Double_t qTot = cluster->GetQ();
  Double_t xMean = x*qTot;
  Double_t zMean = z*qTot;
  Double_t qMax = cluster->GetQmax();
  Double_t xSigma = cluster->GetSig();
  Double_t zSigma = cluster->GetMean();
  Int_t ifirst = cluster->GetIfirst();
  Int_t imin = ifirst%NADC;
  Int_t jmin = ifirst/NADC;
  Int_t width = cluster->GetWidth();
  Int_t imax = width%NADC;
  Int_t jmax = width/NADC;
  fMapReduced[i][j] = -1;
  //  Int_t n = 1;
  for(Int_t i1 = -1; i1<2; i1++){
    for(Int_t j1 = -1; j1<2; j1++){
      if(i1 == j1 || i1 == -j1) continue;
      //Info("CollateClusters","%d %d",i1,j1);
      HarpoRecoClusters* cluster1 = CollateClusters(iCl,i+i1, j+j1, clArray);
      if(!cluster1) continue;
      Double_t x1 = cluster1->GetX();
      Double_t z1 = cluster1->GetZ();
      if(TMath::Abs(x1-x)>fRebinX) continue;
      if(TMath::Abs(z1-z)>fRebinZ) continue;
      //if(TMath::Abs(z1-z)/fRebinZ+TMath::Abs(x1-x)/fRebinX>1) continue;
      Double_t q1 = cluster1->GetQ();
      Int_t ifirst1 = cluster1->GetIfirst();
      Int_t imin1 = ifirst1%NADC;
      Int_t jmin1 = ifirst1/NADC;
      Int_t width1 = cluster1->GetWidth();
      Int_t imax1 = width1%NADC;
      Int_t jmax1 = width1/NADC;
      if(imin1<imin) imin = imin1;
      if(jmin1<jmin) jmin = jmin1;
      if(imax1>imax) imax = imax1;
      if(jmax1>jmax) jmax = jmax1;
      qTot += q1;
      xMean += q1*x1;
      zMean += q1*z1;
      if(q1>qMax) qMax = q1;
      xSigma += q1*cluster1->GetSig();
      zSigma += q1*cluster1->GetMean();
      clArray->Remove(cluster1);
    }
  }

  cluster->SetX(xMean/qTot);
  cluster->SetZ(zMean/qTot);
  cluster->SetQ(qTot);
  cluster->SetQmax(qMax);
  cluster->SetSig(xSigma);
  cluster->SetMean(zSigma);
  cluster->SetIfirst(imin + NADC*jmin);
  cluster->SetWidth(imax + NADC*jmax);
  
  return cluster;
}


Int_t HarpoClusteringBlocs::MakeClusters(Int_t ndet)
{

  HarpoDccMap* m = fEvt->GetDccMap(ndet);
  if ( m == NULL ) return 0;

  if(gHarpoDebug>1)
    Info("MakeClusters","%d",ndet);

  Int_t nCl = 0;

  //  Double_t qTotTot = 0;
  //  TGraph* g = new TGraph();
  for(Int_t j=0;j<NCHAN;j+=fRebinX){ // Channels
    for(Int_t i=0;i<NADC;i+=fRebinZ){ // Time bins
      Double_t xMean = 0, zMean = 0;
      Double_t xSigma = 0, zSigma = 0;
      Double_t qTot = 0, qMax = 0;
      for(Int_t j2 = 0; j2<fRebinX; j2++){
        if(j + j2>=NCHAN) break;
        for(Int_t i2 = 0; i2<fRebinZ; i2++){
          if(i+i2>=NADC) break;
          Double_t q = m->GetData(j+j2,i+i2);
          if(q <= -1000) continue;
          if(q>qMax) qMax = q;
          zMean += q*(i+i2);
          xMean += q*(j+j2);
          zSigma += q*(i+i2)*(i+i2);
          xSigma += q*(j+j2)*(j+j2);
          qTot += q;
        }
      }
      if(qTot<=500) continue;
      xMean /= qTot;
      zMean /= qTot;
      xSigma /= qTot;
      zSigma /= qTot;
      xSigma -= xMean*xMean;
      zSigma -= zMean*zMean;
      fMapReduced[j/fRebinX][i/fRebinZ] = fEvt->GetRecoEvent()->GetNcl();
      //      HarpoRecoClusters *cl = new HarpoRecoClusters(BLOCCLUSTER,qTot,qMax,xMean,i,TMath::Sqrt(xSigma),fRebinZ,fRebinX,-1,ndet,0);
      HarpoRecoClusters *cl = new HarpoRecoClusters(BLOCCLUSTER,qTot,qMax,zSigma,i,xSigma,fRebinZ,fRebinX,-1,ndet,0);
      if(gHarpoDebug>2)
        Info("MakeClusters","%d %g %g",ndet,xMean,zMean);
      cl->SetX(xMean);
      cl->SetZ(zMean);
      cl->SetPlane(ndet);
      cl->SetIfirst(i + NADC*j);
      Int_t imax = i + fRebinZ;
      if(imax>=NADC) imax = NADC-1;
      Int_t jmax = j + fRebinX;
      if(jmax>=NADC) jmax = NCHAN-1;
      cl->SetWidth(imax + NADC*jmax);
      fEvt->GetRecoEvent()->AddClusters(cl);
      cl->Delete();
      nCl++;
    }
  }
  
  return nCl;
}






void   HarpoClusteringBlocs::Init()
{
  
  fRebinX = 6;//5;
  fRebinZ = 12;//10;  
  fNcl = 0;
  fNminBlocs = 4;
  
  Long64_t rebinX = 10;
  if( ! gHConfig->Lookup("ClusteringBlocs.rebinX",rebinX ))
    Info("Constructor","Use default RebinX %d",fRebinX);
  else
    fRebinX = rebinX;

  Long64_t rebinZ = 10;
  if( ! gHConfig->Lookup("ClusteringBlocs.rebinZ",rebinZ ))
    Info("Constructor","Use default RebinZ %d",fRebinZ);
  else
    fRebinZ = rebinZ;

  Long64_t nMinBloc = 4;
  if( ! gHConfig->Lookup("ClusteringBlocs.nMinBloc",nMinBloc ))
    Info("Constructor","Use default RebinZ %d",fNminBlocs);
  else
    fRebinZ = nMinBloc;
  
  fSizeX = Int_t(NCHAN/fRebinX)+1;
  fSizeZ = Int_t(NADC/fRebinZ)+1;
  // Int_t fSizeX1 = Int_t(NCHAN/fRebinX/2)+1;
  // Int_t fSizeZ1 = Int_t(NADC/fRebinZ/2)+1;
  // Int_t fSizeX2 = Int_t(NCHAN/fRebinX/4)+1;
  // Int_t fSizeZ2 = Int_t(NADC/fRebinZ/4)+1;
  
  cout << "fSizeX = " << fSizeX << endl;
  cout << "fSizeZ = " << fSizeZ << endl;

  //  return;

  fMapReduced = new Int_t*[fSizeX];
  for(Int_t i = 0; i<fSizeX; i++)
    fMapReduced[i] = new Int_t[fSizeZ];
  fMapUsed = new Int_t*[fSizeX];
  for(Int_t i = 0; i<fSizeX; i++)
    fMapUsed[i] = new Int_t[fSizeZ];
  // fMapReduced1 = new Int_t*[fSizeX1];
  // for(Int_t i = 0; i<fSizeX1; i++)
  //   fMapReduced1[i] = new Int_t[fSizeZ1];
  // fMapReduced2 = new Int_t*[fSizeX2];
  // for(Int_t i = 0; i<fSizeX2; i++)
  //   fMapReduced2[i] = new Int_t[fSizeZ2];
  
 }

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

   OpenHistFile("recoClusteringBlocs");

   // Save histograms here

 }





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

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

  //  std::cout << "HarpoClusteringBlocs::DisplayAnalysis" << std::endl;
  for(Int_t plane = 0; plane<2;plane++){
    TVirtualPad* cP = c->GetPad(plane+1);
    cP->Divide(1,2);
    TVirtualPad* cMaps = cP->GetPad(1);
    //    TVirtualPad* cProj = cP->GetPad(2);
    
    HarpoDccMap* m = fEvt->GetDccMap(plane);
    TH2F* h = new TH2F("h","",512,0,512,288,0,288);
    for(Int_t i=1;i<NADC;i++){ //Time bins
      for(Int_t j=0;j<NCHAN;j++){ // Channels
        Double_t q = m->GetData(j,i);
        h->SetBinContent(i+1,j+1,q);
      }
    } 
    h->SetMinimum(0);
    
    
    MakeNiceHisto(h,cMaps,"colz");
    h->Delete();



    TGraph* gAll = new TGraph();
    TGraph* g[10];
    for(Int_t i = 0; i<10; i++) g[i] = new TGraph();
    TClonesArray* clArray = fEvt->GetRecoEvent()->GetClustersArray();
    Int_t nPoints = clArray->GetEntries();
    // Info("DisplayAnalysis","nPoints = %d",nPoints);
    for(Int_t i = 0; i<nPoints; i++){
      HarpoRecoClusters* cluster = (HarpoRecoClusters*)clArray->At(i);
      if(cluster->GetPlane()!=plane) continue;
      gAll->SetPoint(i,cluster->GetZ(),cluster->GetX());
      if(gHarpoDebug>0)
        Info("DisplayAnalysis","%d %g %g",plane,cluster->GetZ(),cluster->GetX());
      Int_t iTr = cluster->GetIdClusterTrack();
      if(iTr<0 || iTr>=10) continue;
      g[iTr]->SetPoint(g[iTr]->GetN(),cluster->GetZ(),cluster->GetX());
    }
    gAll->SetMarkerStyle(7);
    gAll->Draw("Psame");
    for(Int_t i = 0; i<10; i++){
      g[i]->SetMarkerStyle(6);
      g[i]->SetMarkerColor(i+2);
      g[i]->Draw("Psame");
    }
  }

  c->Update();
}



void HarpoClusteringBlocs::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, "Clustering 2D blocs");

  TGHorizontalFrame* rebinXFrame = new TGHorizontalFrame(f213,xsize,ysize2,kHorizontalFrame);
  TGLabel* rebinXLabel = new TGLabel(rebinXFrame, "X size");
  fChooseRebinX    = new TGNumberEntry(rebinXFrame);
  fChooseRebinX->SetNumber(fRebinX);

  TGHorizontalFrame* rebinZFrame = new TGHorizontalFrame(f213,xsize,ysize2,kHorizontalFrame);
  TGLabel* rebinZLabel = new TGLabel(rebinZFrame, "Z size");
  fChooseRebinZ    = new TGNumberEntry(rebinZFrame);
  fChooseRebinZ->SetNumber(fRebinZ);

  TGHorizontalFrame* nMinBlocsFrame = new TGHorizontalFrame(f213,xsize,ysize2,kHorizontalFrame);
  TGLabel* nMinBlocsLabel = new TGLabel(nMinBlocsFrame, "WmaxCl");
  fChooseNminBlocs    = new TGNumberEntry(nMinBlocsFrame);
  fChooseNminBlocs->SetNumber(fNminBlocs);




  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(rebinXFrame,fLayout3);
  rebinXFrame->AddFrame(rebinXLabel,fLayout1);
  rebinXFrame->AddFrame(fChooseRebinX,fLayout2);

  f213->AddFrame(rebinZFrame,fLayout3);
  rebinZFrame->AddFrame(rebinZLabel,fLayout1);
  rebinZFrame->AddFrame(fChooseRebinZ,fLayout2);

  f213->AddFrame(nMinBlocsFrame,fLayout3);
  nMinBlocsFrame->AddFrame(nMinBlocsLabel,fLayout1);
  nMinBlocsFrame->AddFrame(fChooseNminBlocs,fLayout2);

  f213->AddFrame(setConf,fLayout3);


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


void   HarpoClusteringBlocs::SetConfig()
{
  
  Info("SetConfig","Renewing configuration");
  
  // fDisplay = 1;
  //  Info("SetConfig","###################### (%d)",fDisplay);

  if(!fChooseRebinX) return;

  fRebinX = fChooseRebinX->GetNumber();
  fRebinZ = fChooseRebinZ->GetNumber();
  fNminBlocs = fChooseNminBlocs->GetNumber();


  for(Int_t i = 0; i<fSizeX; i++)
    delete fMapReduced[i];
  delete fMapReduced;
 
  fSizeX = Int_t(NCHAN/fRebinX)+1;
  fSizeZ = Int_t(NADC/fRebinZ)+1;
 
  fMapReduced = new Int_t*[fSizeX];
  for(Int_t i = 0; i<fSizeX; i++)
    fMapReduced[i] = new Int_t[fSizeZ];
  fMapUsed = new Int_t*[fSizeX];
  for(Int_t i = 0; i<fSizeX; i++)
    fMapUsed[i] = new Int_t[fSizeZ];
  
 


  //  gApplication->Terminate(0);
}