HarpoClustering.cxx

Go to the documentation of this file.

// 
// File HarpoClustering.cxx
//
#include "HarpoClustering.h"
#include "HarpoConfig.h"
#include "HarpoDetSet.h"
#include "HarpoDebug.h"
#include "HarpoDccEvent.h"
#include "Pmm2Event.h"
#include "HarpoEvent.h"
#include "HarpoRecoEvent.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 "TSystem.h"
#include "TPolyLine.h"

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

ClassImp(HarpoClustering);

void   HarpoClustering::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   HarpoClustering::process()
{
  // Main process, calls FindCluster for each map



  //  Bool_t drawEvent = kFALSE;
  nEvents++;
  
  fNcl = 0;  
  fEvt->GetRecoEvent()->ResetClustersArray();
 
  // Process RAW data
  for(Int_t ndet=0;ndet<2;ndet++) { // Detectors 0 and 1 = X and Y planes
    if (!gHDetSet->isExist(ndet)) continue;
    HarpoDccMap* m = fEvt->GetDccMap(ndet);
    if ( m == NULL ) continue;
    fTmin = 1000;
    fTmax = -1000;
    fXmin = 1000;
    fXmax = -1000;
    fNCcl = 0;
    fNTcl = 0;
    fQtot = 0;
    fXstart = -1000;
    fTstart = -1000;

    for(Int_t i = 0; i<NALL; i++) fBadCh[i] = -1;
    fBadCh[93] = 2;  // GEM separators
    fBadCh[94] = 2;  // GEM separators
    fBadCh[193] = 2;  // GEM separators
    fBadCh[194] = 2;  // GEM separators

    if(ndet == 0){
      fBadCh[28] = 2;  // GEM sector limits
      fBadCh[61] = 2;  // GEM sector limits
      fBadCh[127] = 2;  // GEM sector limits
      fBadCh[160] = 2;  // GEM sector limits
      fBadCh[226] = 2;  // GEM sector limits
      fBadCh[259] = 2;  // GEM sector limits
    }

    if(fHist[ndet]){
      fHist[ndet]->GetXaxis()->SetRangeUser(0,511);
      fHist[ndet]->GetYaxis()->SetRangeUser(0,304);
      fCanvas[ndet]->Update();
    }
    // MAKE CLUSTERS    
    FindCluster(m,CCLUSTER,ndet);
    FindCluster(m,TCLUSTER,ndet);
    if(fHist[ndet]){
      fHist[ndet]->GetXaxis()->SetRangeUser(0,511);
      fHist[ndet]->GetYaxis()->SetRangeUser(0,304);
      fCanvas[ndet]->Update();
    }

    if(ndet==XPLANE){
      fEvt->GetRecoEvent()->SetTmin(fTmin);
      fEvt->GetRecoEvent()->SetTmax(fTmax);
      fEvt->GetRecoEvent()->SetXmin(fXmin);
      fEvt->GetRecoEvent()->SetXmax(fXmax);
      fEvt->GetRecoEvent()->SetNCclXEvt(fNCcl);
      fEvt->GetRecoEvent()->SetNTclXEvt(fNTcl);
      fEvt->GetRecoEvent()->SetNclXEvt(fNCcl + fNTcl);
      fEvt->GetRecoEvent()->SetQtotXEvt(fQtot);

      fEvt->GetRecoEvent()->SetXstart(fXstart);
      fEvt->GetRecoEvent()->SetTstart(fTstart);
    }
    if(ndet==YPLANE){
      fEvt->GetRecoEvent()->SetTmin(fTmin);
      fEvt->GetRecoEvent()->SetTmax(fTmax);
      fEvt->GetRecoEvent()->SetYmin(fXmin);
      fEvt->GetRecoEvent()->SetYmax(fXmax);
      fEvt->GetRecoEvent()->SetNCclYEvt(fNCcl);
      fEvt->GetRecoEvent()->SetNTclYEvt(fNTcl);
      fEvt->GetRecoEvent()->SetNclYEvt(fNCcl + fNTcl);
      fEvt->GetRecoEvent()->SetQtotYEvt(fQtot);

      fEvt->GetRecoEvent()->SetYstart(fXstart);
      fEvt->GetRecoEvent()->SetTstart(0.5*(fTstart + fEvt->GetRecoEvent()->GetTstart()));
    }

  }
}
 


void HarpoClustering::FindCluster(HarpoDccMap* m, Int_t type, Int_t plane)
{



  Int_t iMax = 0, jMax = 0;
  switch(type){
  case TCLUSTER:
    jMax = NCHAN;
    iMax = NADC;
    break;
  case CCLUSTER:
    jMax = NADC;
    iMax = NCHAN;
    break;
  default:
    Warning("FindCluster","Unknown cluster type %d",type);
    return;
  }

  Double_t xStart = 0, tStart = 0, qStart = 0;
  Int_t nbinsStart = 0;

  for(Int_t i=0;i<iMax;i++){
    // Initialisation
    Double_t qTot = 0, qTotCh = 0, qMax = 0;
    Double_t mean = 0, sigma = 0;
    Int_t width = 0;
    Double_t qold = -1000;
    Int_t flag = 0;
    Int_t up = -1;
    Bool_t noisy = kFALSE;
    for(Int_t j=0;j<jMax;j++){
      Double_t q, noise, noiseThr, thr;
      Double_t xmin, xmax, ymin, ymax, xp, yp;
      if(type == TCLUSTER){
        q = m->GetData(j,i);
        if(q>-1000){
          fQtot += q;
          if(nbinsStart<4){
            xStart += j*q;
            tStart += i*q;
            qStart += q;
            nbinsStart++;
          }
        }
        noise = m->GetSigmaNoise(j);
        thr = noise*fThr;
        // if(fBadCh[j]>0){
        //   if(q==-1000) q=0;
        //   noise = 500*(7*fBadCh[j] - 6) ;
        // }
        if(fBadCh[j]>1){
          //noisy = kTRUE;
          noise = 500*(7*fBadCh[j] - 6) ;
          if(gHarpoDebug>2)
            Info("FindCluster","Dead Channel %d (%d, %d)", j, i, up);
          //continue;
        }
        noiseThr = noise*fNoiseThr;
        xmin = i - 10;
        xmax = i + 10;
        ymin = j - 10;
        ymax = j + 10;
        xp = i;
        yp = j;
      }else{
        q = m->GetData(i,j);
        noise = m->GetSigmaNoise(i);
        noiseThr = noise*fNoiseThr;
        thr = noise*fThr;

        xmin = j - 10;
        xmax = j + 10;
        ymin = i - 10;
        ymax = i + 10;
        xp = j;
        yp = i;
      }
      if(q>10){
        if(xp>fTmax) fTmax = xp;
        if(xp<fTmin) fTmin = xp;
        if(yp>fXmax) fXmax = yp;
        if(yp<fXmin) fXmin = yp;
      }
      if(fHist[plane] != 0 && q>-1000){
        //gSystem->Sleep(1);
        if(xmin<0) xmin = 0;
        if(ymin<0) ymin = 0;
        if(xmax>NADC) xmax = NADC;
        if(ymax>NCHAN) ymax = NCHAN;
        fCanvas[plane]->cd();
        Float_t x[4]; // = {xp,xp,xp+1,xp+1};
        x[0] = xp;
        x[1] = xp;
        x[2] = xp+1;
        x[3] = xp+1;
        Float_t y[4]; // = {yp,yp+1,yp+1,yp};
        y[0] = yp;
        y[1] = yp+1;
        y[2] = yp+1;
        y[3] = yp;
        TPolyLine* triangle = new TPolyLine(4,x,y,"FL");
        triangle->SetLineColor(kBlack);
        triangle->SetFillColor(3+up);
        //triangle->SetFillStyle(3002);
        fCanvas[plane]->cd();
        triangle->Draw("FL");
        // if(plane == 1 && yp == 141)
        //   std::cout << noiseThr << "  " << thr << "    " << type << std::endl;

        // fHist[plane]->GetXaxis()->SetRangeUser(xmin,xmax);
        // fHist[plane]->GetYaxis()->SetRangeUser(ymin,ymax);
        //fCanvas[plane]->Update();
      }

      if(q!=-1000.)
        qTotCh += q;
     // if(q != -1000)
     //         std::cout << "[" << q << "," << up << "]" << std::flush;
      //      if(up == 1 && q>thr && q < qold - noiseThr){ // Peak reached
      if(up == 1 && q < qold - noiseThr){ // Peak reached
        up = 2; // DOWN
      }
                        // DOWN slope
      if((up == 2 && q > qold + noiseThr) || (up != -1 && (q < thr || j == jMax-1))){ // End of cluster
        if(gHarpoDebug>2){
          Info("FindCluster","%d %d (%d, %d)", i, j, type, plane);
          Info("FindCluster","q=%f qold=%f noise=%f fThr=%f", q, qold, noise, fThr);
        }
        // std::cout << std::endl;
        // std::cout << "--------- " << qTot << " -------" << std::endl;
        mean /= qTot;
        sigma /= qTot;
        sigma -= mean*mean;
        if(q > qold + noiseThr || j == jMax-1)
          flag = flag | 1; // Cut After;
        if(noisy){
          flag &= 4;
          if(gHarpoDebug>1)
            Info("FindCluster","Noisy Cluster %d %d %d %g",fNcl,flag,type,mean);
        }// else
        if(gHarpoDebug>2)
          Info("FindCluster","%d: iFirst = %d",fNcl,j-width);
        HarpoRecoClusters *cl = new HarpoRecoClusters(type,qTot,qMax,mean,i,TMath::Sqrt(sigma),j-width,width,-1,plane,flag);
        //fEvt->GetRecoEvent()->SetClusters(cl,fNcl);
        if(type == CCLUSTER){
          cl->SetX(i);
          cl->SetZ(mean);
        }
        if(type == TCLUSTER){
          cl->SetZ(i);
          cl->SetX(mean);
        }
        fEvt->GetRecoEvent()->AddClusters(cl);
        fNcl++;
        if(type == CCLUSTER) fNCcl++;
        if(type == TCLUSTER) fNTcl++;
        cl->Delete();
        qTot = 0;
        qMax = 0;
        mean = 0;
        sigma = 0;
        width = 0;
        if(flag&1)
          flag = 2; // Cut before
        else
          flag = 0;
        up = -1; //kTRUE; // OUT OF CLUSTER
        noisy = kFALSE;
      }
      // if(plane == 1 && yp == 141 && type == 0)
      //        std::cout << q << "  " << qold << "    " << up << std::endl;
      if(q<thr){
        qold = q;
        continue;      
      }
      if(q > qold + noiseThr || (qold<thr && q>thr))
        up = 1;



      if(type == TCLUSTER && fBadCh[j]>1)
        noisy = kTRUE;
      if(j-1>0)
        if(type == TCLUSTER && fBadCh[j-1]>1)
          noisy = kTRUE;
      if(j+1<NCHAN)
        if(type == TCLUSTER && fBadCh[j+1]>1)
          noisy = kTRUE;
      

      if(j == 0) flag = 2; // Cut before
      qTot += q;
      mean += q*j;
      sigma += q*j*j;
      width +=1;
      if(q>qMax) qMax = q;
      qold = q;



    }

    //std::cout << qTotCh << std::endl;
    // if(type == 0 && qTotCh>6.2e4){
    //   //Info("FindCluster","Bad channel %d (%g) [%d %d %d]",i,qTotCh,fBadCh[i-1],fBadCh[i],fBadCh[i+1]);
    //   fBadCh[i-1] += 1;
    //   fBadCh[i] += 1;
    //   fBadCh[i+1] += 1;
    //   //Info("FindCluster","**Bad channel %d (%g) [%d %d %d]",i,qTotCh,fBadCh[i-1],fBadCh[i],fBadCh[i+1]);
    // }
  }


  if(type == TCLUSTER){
    if(qStart>0){
      fXstart = xStart/qStart;
      fTstart = tStart/qStart;
    }else{
      fXstart = -1000;
      fTstart = -1000;
    }
  }


  //  std::cout << fNcl << std::endl;

}

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

  TCanvas* cTab = ecTab->GetCanvas();
  HarpoRecoEvent* reco = fEvt->GetRecoEvent();
  TClonesArray* clArray = reco->GetClustersArray();
  Int_t nCl = clArray->GetEntries();

  TGraph* gCX = new TGraph();
  TGraph* gCY = new TGraph();
  TGraph* gTX = new TGraph();
  TGraph* gTY = new TGraph();

  TGraph* gTXnoise = new TGraph();
  TGraph* gTYnoise = new TGraph();
  gTXnoise->SetMarkerStyle(28);
  gTYnoise->SetMarkerStyle(28);
  // TGraph* gCXnotrack = new TGraph();
  // TGraph* gCYnotrack = new TGraph();
  // TGraph* gTXnotrack = new TGraph();
  // TGraph* gTYnotrack = new TGraph();

  for(Int_t icl = 0; icl<nCl; icl++){
    HarpoRecoClusters* cluster = (HarpoRecoClusters*)clArray->At(icl);
    Int_t ind = cluster->GetIndex();
    Double_t mean = cluster->GetMean();
    Int_t plane = cluster->GetPlane();
    Int_t type = cluster->GetType();
    //Int_t width = cluster->GetWidth();
    //Double_t q = cluster->GetQ();
    //Int_t idcluster = cluster->GetIdClusterTrack();
    if(plane!=XPLANE && plane !=YPLANE) continue;
    // if(type == CCLUSTER && width<10) continue;
    // if(type == TCLUSTER && width>10) continue;
    if(plane==XPLANE && type == CCLUSTER) gCX->SetPoint(gCX->GetN(),mean,ind+0.5);
    if(plane==XPLANE && type == TCLUSTER) gTX->SetPoint(gTX->GetN(),ind+0.5,mean);
    if(plane==YPLANE && type == CCLUSTER) gCY->SetPoint(gCY->GetN(),mean,ind+0.5);
    if(plane==YPLANE && type == TCLUSTER) gTY->SetPoint(gTY->GetN(),ind+0.5,mean);



    if(cluster->GetQuality()<0){
      //      std::cout << "Noisy cluster: " << icl << "  " << cluster->GetQuality() << std::endl;
      if(plane==XPLANE && type == TCLUSTER) gTXnoise->SetPoint(gTXnoise->GetN(),ind+0.5,mean);
      if(plane==YPLANE && type == TCLUSTER) gTYnoise->SetPoint(gTYnoise->GetN(),ind+0.5,mean);
    }
    // if(idcluster>=0) continue; 
    // if(plane==0 && type == 0) gCXnotrack->SetPoint(gCXnotrack->GetN(),mean,ind);
    // if(plane==0 && type == 1) gTXnotrack->SetPoint(gTXnotrack->GetN(),ind,mean);
    // if(plane==1 && type == 0) gCYnotrack->SetPoint(gCYnotrack->GetN(),mean,ind);
    // if(plane==1 && type == 1) gTYnotrack->SetPoint(gTYnotrack->GetN(),ind,mean);
  }

  gCX->SetMarkerStyle(7);
  gTX->SetMarkerStyle(7);
  gCY->SetMarkerStyle(7);
  gTY->SetMarkerStyle(7);
  gCX->SetMarkerSize(0.5);
  gTX->SetMarkerSize(0.5);
  gCY->SetMarkerSize(0.5);
  gTY->SetMarkerSize(0.5);
  gCX->SetMarkerColor(kGray);
  gCY->SetMarkerColor(kGray);
  // gCXnotrack->SetMarkerStyle(7);
  // gTXnotrack->SetMarkerStyle(7);
  // gCYnotrack->SetMarkerStyle(7);
  // gTYnotrack->SetMarkerStyle(7);
  
  std::cout << gCX->GetN() << " " << gCY->GetN() << " " << gTX->GetN() << " " << gTY->GetN() << std::endl;


  cTab->cd(1);
  if(gCX->GetN()>0) gCX->Draw("Psame");
  if(gTX->GetN()>0) gTX->Draw("Psame");
  if(gTXnoise->GetN()>0) gTXnoise->Draw("Psame");
  cTab->cd(2);
  if(gCY->GetN()>0) gCY->Draw("Psame");
  if(gTY->GetN()>0) gTY->Draw("Psame");
  if(gTYnoise->GetN()>0) gTYnoise->Draw("Psame");

  cTab->Update();
}







 void   HarpoClustering::Init()
 {

   for(Int_t i = 0; i<NALL; i++) fBadCh[i] = 0;
 
   fNcl = 0;

   fThr = 3;
   Double_t thr;
   if ( !gHConfig->Lookup("Clustering.Thr",thr) )
     Info("Init","Use default Threshold %g",fThr);
   else{
     fThr = thr;
     Info("Init","Use Threshold %g",fThr);
   }

   fNoiseThr = 3;
   Double_t noiseThr;
   if ( !gHConfig->Lookup("Clustering.NoiseThr",noiseThr) )
     Info("Init","Use default Noise Threshold %g",fNoiseThr);
   else{
     fNoiseThr = noiseThr;
     Info("Init","Use Noise Threshold %g",fNoiseThr);
   }


   fHist[0] = 0;
   fHist[1] = 0;
   fCanvas[0] = 0;
   fCanvas[1] = 0;

   fXstart = -1000;
   fTstart = -1000;
 }

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

   // TString * hstFile = gHConfig->GetHistFile();
   // if ( hstFile == NULL ) {
   //   Info("Save","No Hist File name given, use default");
   //   hstFile = new TString(Form("outputs/anaph%lli.root",gHConfig->GetRunNo()));
   //   Info("Save","No Hist File name given, use default %s",hstFile->Data());
   // }

   // TFile *hf = new TFile(hstFile->Data(),"RECREATE");
   // // Save histograms here
   // hf->Close();
   // printf("fNewFile %s closed\n", hstFile->Data() );
 
 }