HarpoDccAnalyseNoiseMonitor.cxx

Go to the documentation of this file.

// 
// File HarpoDccAnalyseNoiseMonitor.cxx
//
#include "HarpoDccAnalyseNoiseMonitor.h"
#include "HarpoConfig.h"
//#include "HarpoDetSet.h"
#include "HarpoDccEvent.h"

#include "TROOT.h"
#include "TStyle.h"
#include "TApplication.h"
#include "TFile.h"
#include "TCanvas.h"
#include "TLatex.h"
#include "TGraph.h"
#include "TF1.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TLinearFitter.h"
#include "TMath.h"

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

ClassImp(HarpoDccAnalyseNoiseMonitor)

void   HarpoDccAnalyseNoiseMonitor::print()
{
  //  int nd; // number of detectors
  //HarpoEventHeader *hdr = anaEvt->GetHeader();
  
  // assert(hdr != NULL);
  // hdr->print();
  
  // for (nd = 0; nd < gkNDetectors; nd++) {
  //   //     if (anaEvt->isdataExist(nd)) {
  //   HarpoDccMap *plane = anaEvt->GetDccMap(nd);
  //   if (plane != NULL )plane->print();
  // }
}

void HarpoDccAnalyseNoiseMonitor::process()
{
  nEvents++;
  std::cout << " Processing Event " << nEvents<< std::endl;




  
  TCanvas *c = (TCanvas *)gROOT->GetListOfCanvases()->FindObject("cHarpo");
  if(!c){
    c = new TCanvas("cHarpo","DCC: All Events cumulated",0,0,1000,800);
    c->Clear();
    c->Divide(3,3);//,0.001,0.001);
  }
   
  TCanvas *c2 = (TCanvas *)gROOT->GetListOfCanvases()->FindObject("cHarpo2");
  if(!c2){
    c2 = new TCanvas("cHarpo2","DCC: Current Event",0,0,1000,800);
    c2->Clear();
    c2->Divide(3,3);//,0.001,0.001);
  }
   
  gStyle->SetOptStat(0);
 
  TH1F* h1Ev1Ch = new TH1F("h1Ev1Ch",";time [tb]",511,0,511);
  TH1F* h1EvAllCh = new TH1F("h1EvAllCh",";time [tb]",511,0,511);
  TH2F* h1Ev2D = new TH2F("h1Ev2D",";time [tb]; channel",511,0,511,304,0,304);

  TH1F* h1EvFFT1Ch = new TH1F("h1EvFFT1Ch",";freq []",511,0,511);
  TH1F* h1EvFFTAllCh = new TH1F("h1EvFFTAllCh",";freq []",511,0,511);
  TH2F* h1EvFFT2D = new TH2F("h1EvFFT2D",";freq []; channel",511,0,511,304,0,304);

  TH1F* h1EvPed1Ch = new TH1F("h1EvPed1Ch",";pedesteal [ADC]",500,0,500);
  TH1F* h1EvPedAllCh = new TH1F("h1EvPedAllCh",";pedesteal [ADC]",500,0,500);
  TH2F* h1EvPed2D = new TH2F("h1EvPed2D",";pedesteal [ADC]; channel",500,0,500,304,0,304);
 
  // Int_t npix = 0, nch = 0, ntb = 0;
  
  HarpoDccMap * m = anaEvt->GetMap();
  std::cout << "Event " << nEvents << "  " << m << std::endl;
  if ( m != NULL ) {
    //    for(j=0;j<NCHAN;j++) // Channels
    //    Double_t qTot = 0;
    //TH2F* hFull = new TH2F("hFull",";time [25ns];channel [1mm]",511,0,511,304,0,304);
    
    TH1F* hTmp = new TH1F("hTmp","",511,0,511);
    for(Int_t i=0;i<NADC;i++){ //Time bins
      Double_t qTotT = 0;
      for(Int_t j=0;j<NALL;j++){ // Channels
        //      Double_t q = m->map(j,i);
        Double_t q = m->GetData(j,i);
        h2D->SetBinContent(i+1,j+1,h2D->GetBinContent(i+1,j+1)+q);
        h1Ev2D->SetBinContent(i+1,j+1,q);
        if(j == 100){
          h1Ch->SetBinContent(i+1,h1Ch->GetBinContent(i+1)+q);
          h1Ev1Ch->SetBinContent(i+1,q);
          hPed1Ch->Fill(q);
          h1EvPed1Ch->Fill(q);
        }
        hPed2D->Fill(q,j);
        h1EvPed2D->Fill(q,j);
        qTotT += q;
      }
      hAllCh->SetBinContent(i+1,hAllCh->GetBinContent(i+1)+qTotT);
      h1EvAllCh->SetBinContent(i+1,qTotT);
      hPedAllCh->Fill(qTotT/288);
      h1EvPedAllCh->Fill(qTotT/288);
      hTmp->SetBinContent(i+1,qTotT);
    }
    TH1* hm1 = 0;
    hm1 = hTmp->FFT(hm1, "MAG");
    for(Int_t i=2;i<512;i++){ //Time bins
      hFFTAllCh->SetBinContent(i+1,hFFTAllCh->GetBinContent(i+1)+hm1->GetBinContent(i+1));
      h1EvFFTAllCh->SetBinContent(i+1,hm1->GetBinContent(i+1));
    }
    hTmp->Delete();
    hm1->Delete();

    for(Int_t j=0;j<NALL;j++){ // Channels
      TH1F* hTmp = new TH1F("hTmp","",511,0,511);
      for(Int_t i=0;i<NADC;i++){ //Time bins
        //      Double_t q = m->map(j,i);
        Double_t q = m->GetData(j,i);
        hTmp->SetBinContent(i+1,q);
      }

      TH1* hm = 0;
      hm = hTmp->FFT(hm, "MAG");
      for(Int_t i=2;i<512;i++){ //Time bins
        hFFT2D->SetBinContent(i+1,j+1,hFFT2D->GetBinContent(i+1,j+1)+hm->GetBinContent(i+1));
        h1EvFFT2D->SetBinContent(i+1,j+1,hm->GetBinContent(i+1));
        if(j==100){
          hFFT1Ch->SetBinContent(i+1,hFFT1Ch->GetBinContent(i+1,j+1)+hm->GetBinContent(i+1));
          h1EvFFT1Ch->SetBinContent(i+1,hm->GetBinContent(i+1));
        }
      }
      hTmp->Delete();
      hm->Delete();

    }

    hFFT1Ch->GetXaxis()->SetRangeUser(-1,255);
    hFFTAllCh->GetXaxis()->SetRangeUser(-1,255);
    hFFT2D->GetXaxis()->SetRangeUser(-1,255);
    h1EvFFT1Ch->GetXaxis()->SetRangeUser(-1,255);
    h1EvFFTAllCh->GetXaxis()->SetRangeUser(-1,255);
    h1EvFFT2D->GetXaxis()->SetRangeUser(-1,255);


    c->cd(1);
    h1Ch->DrawCopy();
    c->cd(2);
    hAllCh->DrawCopy();
    c->cd(3);
    h2D->DrawCopy("colz");

    c->cd(4);
    hFFT1Ch->DrawCopy();
    c->cd(5);
    hFFTAllCh->DrawCopy();
    c->cd(6);
    hFFT2D->DrawCopy("colz");

    c->cd(7);
    hPed1Ch->DrawCopy();
    c->cd(8);
    hPedAllCh->DrawCopy();
    c->cd(9);
    hPed2D->DrawCopy("colz");

    c2->cd(1);
    h1Ev1Ch->DrawCopy();
    c2->cd(2);
    h1EvAllCh->DrawCopy();
    c2->cd(3);
    h1Ev2D->DrawCopy("colz");

    c2->cd(4);
    h1EvFFT1Ch->DrawCopy();
    c2->cd(5);
    h1EvFFTAllCh->DrawCopy();
    c2->cd(6);
    h1EvFFT2D->DrawCopy("colz");

    c2->cd(7);
    h1EvPed1Ch->DrawCopy();
    c2->cd(8);
    h1EvPedAllCh->DrawCopy();
    c2->cd(9);
    h1EvPed2D->DrawCopy("colz");

  TFile* file = new TFile("noise1Ev.root","recreate");

  h1Ev1Ch->Write();
  h1EvAllCh->Write();
  h1Ev2D->Write();

  h1EvFFT1Ch->Write();
  h1EvFFTAllCh->Write();
  h1EvFFT2D->Write();

  h1EvPed1Ch->Write();
  h1EvPedAllCh->Write();
  h1EvPed2D->Write();

  file->Close();
    
    h1Ev1Ch->Delete();
    h1EvAllCh->Delete();
    h1Ev2D->Delete();
    h1EvFFT1Ch->Delete();
    h1EvFFTAllCh->Delete();
    h1EvFFT2D->Delete();
    h1EvPed1Ch->Delete();
    h1EvPedAllCh->Delete();
    h1EvPed2D->Delete();



  } else {
    std::cout << "No Data Map for plane " << std::endl;
  }
  
  c->Update();
  c2->Update();

  
  int iscan;
  printf("Event %d done. Waiting for input to continue: ",nEvents);
  if(scanf("%d",&iscan) == 1)
    std::cout << "Advance to next event" << std::endl;
  else
    std::cout << "Scan all events" << std::endl;
    
}

void   HarpoDccAnalyseNoiseMonitor::Init()
{
  // m_hmap[0] = new HarpoMap(0);
  // m_hmap[1] = new HarpoMap(1);

  h1Ch = new TH1F("h1Ch",";time [tb]",511,0,511);
  hAllCh = new TH1F("hAllCh",";time [tb]",511,0,511);
  h2D = new TH2F("h2D",";time [tb]; channel",511,0,511,304,0,304);

  hFFT1Ch = new TH1F("hFFT1Ch",";freq []",511,0,511);
  hFFTAllCh = new TH1F("hFFTAllCh",";freq []",511,0,511);
  hFFT2D = new TH2F("hFFT2D",";freq []; channel",511,0,511,304,0,304);

  hPed1Ch = new TH1F("hPed1Ch",";pedesteal [ADC]",500,0,500);
  hPedAllCh = new TH1F("hPedAllCh",";pedesteal [ADC]",500,0,500);
  hPed2D = new TH2F("hPed2D",";pedesteal [ADC]; channel",500,0,500,304,0,304);

}

void   HarpoDccAnalyseNoiseMonitor::Save(char * /* mode */)
{
  TCanvas *c = (TCanvas *)gROOT->GetListOfCanvases()->FindObject("cHarpo");
  if(c){
    c->SaveAs("noise.png");
    c->SaveAs("noise.C");
  }
  // TString * hstFile = gHConfig->GetHistFile();
  // if ( hstFile == NULL ) {
  //   std::cout << gHConfig->GetProgramName() << " " <<
  //     "No Hist File name gived, use default" << std::endl;
  //   hstFile = new TString("harpohst.root");
  // }
  
  TFile* file = new TFile("noise.root","recreate");

  h1Ch->Write();
  hAllCh->Write();
  h2D->Write();

  hFFT1Ch->Write();
  hFFTAllCh->Write();
  hFFT2D->Write();

  hPed1Ch->Write();
  hPedAllCh->Write();
  hPed2D->Write();

  file->Close();
  // TFile *hf = new TFile(hstFile->Data(),"RECREATE");
  // hf->Close();
}







// Double_t HarpoDccAnalyseNoiseMonitor::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++;
//   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]);
//   if(tHigh-tLow) {
//     return truncMean/(tHigh-tLow);
//   }
//   else return 0;
// }