HarpoAnalyseNoiseSuppression.cxx

Go to the documentation of this file.

// 
// File HarpoAnalyseNoiseSuppression.cxx
//
#include "HarpoAnalyseNoiseSuppression.h"
#include "HarpoConfig.h"
#include "HarpoDetSet.h"
//#include "HarpoDetSet.h"
#include "HarpoDccHeader.h"
#include "HarpoFeminosHeader.h"
#include "HarpoDccEvent.h"
#include "HarpoEvent.h"
#include "Pmm2Event.h"
#include "HarpoSimEvent.h"

#include "TFile.h"
#include "TStyle.h"
#include "TCanvas.h"
#include "TLatex.h"
#include "TGraph.h"
#include "TF1.h"
#include "TLinearFitter.h"
#include "TMath.h"

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

ClassImp(HarpoAnalyseNoiseSuppression);

const int bufsize = 16000;
const int splitlevel = 2;
const int compress = 6;   //default ???



HarpoAnalyseNoiseSuppression::HarpoAnalyseNoiseSuppression()
{
  

 }

void   HarpoAnalyseNoiseSuppression::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   HarpoAnalyseNoiseSuppression::process()
 {
   nEvents++;
   if(nEvents%1000 == 0)
     std::cout << " Processing Event " << nEvents<< std::endl;
   
   if(isZS>3) return;

   TLinearFitter* fFitCh = new TLinearFitter(1,"pol2","D");

   //   Bool_t noZeroSuppr = kFALSE;
   for(Int_t ndet=0;ndet<2;ndet++) {
     if (!gHDetSet->isExist(ndet)){
       continue;
     }
     if(gHarpoDebug>1) Info("process","Processing detector %d", ndet);
     HarpoDccMap * m = fEvt->GetDccMap(ndet);
     std::cout << "Event " << nEvents << ", detector " << ndet << "  " << m << std::endl;
     if ( m == NULL ){ 
       Warning("process","Event %ld, No Data Map for plane %d", nEvents, ndet);
       continue;
     }
     Int_t Qnoise = 0;
     for(Int_t i=0;i<NADC;i++){ //Time bins
       for(Int_t j=NCHAN;j<NALL;j++){ // Channels
         Double_t q = m->GetRawData(j,i);
         Qnoise += q;
       }
     }
     //  if(!((HarpoDccEvent*)fEvt->GetDetEvent(ndet))->IsZeroSuppressed()){
     if(Qnoise<NADC*50){
       //if(gHarpoDebug>2) 
         Info("process","Data ZERO SUPPRESSED (%d)",Qnoise);
       isZS++;
       continue;
     }
     if(gHarpoDebug>2) Info("process","Data FULL (%d)",Qnoise);
     for(Int_t i=0;i<NADC;i++){ //Time bins
       TArrayD* valuesTB = new TArrayD(NALL);
       for(Int_t j=0;j<NALL;j++){ // Channels
         Double_t q = m->GetData(j,i);
         valuesTB->AddAt(q,j);
       }
       Double_t noiseTB = TruncMean(valuesTB,0.16,0.84);
       TGraph* gMinusCh = new TGraph();
       for(Int_t j=0;j<NALL;j++){ // Channels
         Double_t q = m->GetData(j,i);
         q -= noiseTB;
         m->SetData(j,i,q);
         gMinusCh->SetPoint(gMinusCh->GetN(),j,q);
       }
       fFitCh->AssignData(gMinusCh->GetN(), 1, gMinusCh->GetX(), gMinusCh->GetY());
       fFitCh->Eval();
       delete gMinusCh;
       Double_t mean = 0, rms = 0;
       Int_t nnoise = 0;
       for(Int_t j=0;j<NALL;j++){ // Channels
         Double_t q = m->GetData(j,i);
         q -= fFitCh->GetParameter(0) + j*fFitCh->GetParameter(1) + j*j*fFitCh->GetParameter(2);
         m->SetData(j,i,q);
         if(q<100){
           mean += q;
           rms += q*q;
           nnoise++;
         }
       }
       rms = TMath::Sqrt(rms/nnoise - mean*mean/nnoise/nnoise);
       for(Int_t j=0;j<NALL;j++){ // Channels
         Double_t q = m->GetData(j,i);
         if(gHarpoDebug>4) Info("process","%d %d => %g (%g)",j,i,q,fNoiseCut*rms);
         if(q>fNoiseCut*rms){
           m->SetData(j,i,q);
         }else{
           m->SetRawData(j,i,-1);
         }
       }
       delete valuesTB; 
       //((HarpoDccEvent*)fEvt->GetDetEvent(ndet))->SetZeroSuppressed(false);
       //((HarpoDccEvent*)fEvt->GetDetEvent(ndet))->SetZeroSuppressed(false);
     }
   }
   fFitCh->Delete();
   

}

 void   HarpoAnalyseNoiseSuppression::Init()
 {


  fNoiseCut = 5.;
  Double_t noiseCut;
  if ( ! gHConfig->Lookup("NoiseSuppr.NoiseCut",noiseCut) )
    Info("Constructor","Use default Noise Cut %g",fNoiseCut);
  else{
    fNoiseCut = noiseCut;
    Info("Constructor","Use Noise Cut %g",fNoiseCut);
  }

   isZS = 0;

 }

void   HarpoAnalyseNoiseSuppression::Save(char * /* mode */ )
 {
 
 }

Double_t HarpoAnalyseNoiseSuppression::TruncMean(TArrayD* vect, Double_t tl, Double_t th)
{
  //
  // Calculates the truncated mean mean of the non zero value in vect
  // The truncation is done on the 100*tl% lowest and 100*(1-th) highest value
  //
  //  debug(2,"truncated mean");


  Int_t size = vect->GetSize();
  Double_t truncMean = 0;
  Int_t* index = new Int_t[size];
  TMath::Sort(size,vect->GetArray(),index,kFALSE);
  Int_t t = 0, tLow, tHigh;


  //  while(vect->At(index[t])<10&&t<size-1) t++;
  while(vect->At(index[t])==-1000&&t<size-1) t++;
  tLow  = TMath::FloorNint(t + (size - t)*tl);
  tHigh = TMath::FloorNint(t + (size - t)*th);

  for(Int_t tind = tLow; tind<tHigh; tind++)  truncMean += vect->At(index[tind]);
  delete[] index;
  if(tHigh-tLow) {
    return truncMean/(tHigh-tLow);
  }
  else return 0;
}