HarpoAnalysePedestalShift.cxx

Go to the documentation of this file.

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

#include "TFile.h"
#include "TStyle.h"
#include "TCanvas.h"
#include "TLatex.h"
#include "TGraph.h"
#include "TF1.h"
#include "TMath.h"
#include "TSystem.h"
#include "TRootEmbeddedCanvas.h"
#include "TGListBox.h"
#include "TGLabel.h"

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

ClassImp(HarpoAnalysePedestalShift);

void   HarpoAnalysePedestalShift::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   HarpoAnalysePedestalShift::process()
{
  //  Bool_t drawEvent = kFALSE;
  nEvents++;
  
  Double_t QchMax[2];
  QchMax[0] = 0;
  QchMax[1] = 0;
  Double_t qtot[2];
  qtot[0] = 0;
  qtot[1] = 0;
  // 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;
    hMapRaw[ndet]->Reset();

    for(Int_t j=0;j<NCHAN;j++){ // Channels
      Double_t qch = 0;
      for(Int_t i=0;i<NADC;i++){ // Time bins
        Double_t q = m->GetData(j,i);
        if(qch>4e4 && q<=0)
          m->SetData(j,i,-2000);
        if(q <= 0) continue;
        qch += q;
        hMapRaw[ndet]->Fill(i,j,qch);
        hQraw[ndet]->Fill(q);
      }
      qtot[ndet] += qch;
      hQch[ndet]->Fill(qch);
      if(qch>QchMax[ndet]) QchMax[ndet] = qch;
      if(qch>QchMax[ndet]) QchMax[ndet] = qch;
      Int_t index = 0;
      for(Int_t i=NADC-1;i>=0;i--){
        Double_t q = m->GetRawData(j,i);
        if(q<=0) continue;
        if(qch>5e4) hQsat[ndet]->Fill(index,q);
        else hQnosat[ndet]->Fill(index,q);
        if(index == 0) hQlast[ndet]->Fill(qch,q);
        index++;
      }
      index = 0;
      for(Int_t i=0;i<NADC;i++){
        Double_t q = m->GetRawData(j,i);
        if(q<=0) continue;
        if(qch>5e4) hQsat1[ndet]->Fill(index,q);
        else hQnosat1[ndet]->Fill(index,q);
        if(index == 0) hQfirst[ndet]->Fill(qch,q);
        index++;
      }
    }
    if(QchMax[ndet]>2e4 && gHarpoDebug>1)
      Info("process","Map %d saturated",ndet);
  }


  //std::cout << QchMax[1] << "  " << qtot[1] << "  " <<  QchMax[0] << "  " << qtot[0] << std::endl;
  if(qtot[0]<=0) return;
  if(qtot[1]<=0) return;
  //  if(QchMax[1]<45e3){
    hRatioVsQchX->Fill(QchMax[0],qtot[0]/qtot[1]);
    hRatioVsQchXVsQtotX->Fill(QchMax[0],qtot[0],qtot[0]/qtot[1]);
    //  }
    // if(QchMax[0]<45e3){
    hRatioVsQchY->Fill(QchMax[1],qtot[0]/qtot[1]);
    hRatioVsQchYVsQtotY->Fill(QchMax[1],qtot[1],qtot[0]/qtot[1]);
    //}
  //std::cout << QchMax[1] << "  " << qtot[1] << "  " <<  QchMax[0] << "  " << qtot[0] << "  " << qtot[0]/qtot[1] << std::endl;

  
}

void   HarpoAnalysePedestalShift::Init()
{
  
  // Initialise histograms here
  hMapRaw[0] = new TH2F("hMapRawX","map;Time;Channel",511,0,511,288,0,288);
  hMapRaw[1] = new TH2F("hMapRawY","map",511,0,511,288,0,288);
  
  hPmm2 = new TH2F("hPmm2","map",16,0,16,1024,0,1024);
  hAngle = new TH1F("hAngle","map",200,0,TMath::Pi());
  
  hQraw[0] = new TH1F("hQrawX",";Q_{raw}",4096,0,4096);
  hQraw[1] = new TH1F("hQrawY",";Q_{raw}",4096,0,4096);
  
  
  hQsat[0] = new TH2F("hQsatX",";index;Q;",511,0,511,4096,0,4096);
  hQsat[1] = new TH2F("hQsatY",";index;Q;",511,0,511,4096,0,4096);
  hQnosat[0] = new TH2F("hQnosatX",";index;Q;",511,0,511,4096,0,4096);
  hQnosat[1] = new TH2F("hQnosatY",";index;Q;",511,0,511,4096,0,4096);

  hQsat1[0] = new TH2F("hQsat1X",";index;Q;",511,0,511,4096,0,4096);
  hQsat1[1] = new TH2F("hQsat1Y",";index;Q;",511,0,511,4096,0,4096);
  hQnosat1[0] = new TH2F("hQnosat1X",";index;Q;",511,0,511,4096,0,4096);
  hQnosat1[1] = new TH2F("hQnosat1Y",";index;Q;",511,0,511,4096,0,4096);

  


   Int_t fNbins = 2000;
   Double_t fQmin = 100, fQmax = 1e7;
   
   Long64_t Nbins;
   if( ! gHConfig->Lookup("pedshift.nbins",Nbins ))
     Info("Constructor","Use default fnbinr %d",fNbins);
   else
     fNbins = Nbins;

   Double_t qmin;
   if( ! gHConfig->Lookup("pedshift.qmin",qmin ))
     Info("Constructor","Use default Qmin %.3g",fQmin);
   else
     fQmin = qmin;

   Double_t qmax;
   if( ! gHConfig->Lookup("pedshift.qmax",qmax ))
     Info("Constructor","Use default Qmax %.3g",fQmax);
   else
     fQmax = qmax;



  const Int_t nbinsQcl = fNbins;
  Double_t xminQcl = fQmin;
  Double_t xmaxQcl = fQmax;
  Double_t logxminQcl = TMath::Log(xminQcl);
  Double_t logxmaxQcl = TMath::Log(xmaxQcl);
  Double_t binwidthQcl = (logxmaxQcl-logxminQcl)/nbinsQcl;
  Double_t xbinsQcl[nbinsQcl+1];
  xbinsQcl[0] = xminQcl;
  for (Int_t i=1;i<=nbinsQcl;i++)
    xbinsQcl[i] = TMath::Exp(logxminQcl+i*binwidthQcl);

  const Int_t nbinsQcl2 = 200;
  Double_t xminQcl2 = 1e2;
  Double_t xmaxQcl2 = 1e7;
  Double_t logxminQcl2 = TMath::Log(xminQcl2);
  Double_t logxmaxQcl2 = TMath::Log(xmaxQcl2);
  Double_t binwidthQcl2 = (logxmaxQcl2-logxminQcl2)/nbinsQcl2;
  Double_t xbinsQcl2[nbinsQcl2+1];
  xbinsQcl2[0] = xminQcl2;
  for (Int_t i=1;i<=nbinsQcl2;i++)
    xbinsQcl2[i] = TMath::Exp(logxminQcl2+i*binwidthQcl2);

  hQcl[0] = new TH1F("hQclX",";Q_{cl}",nbinsQcl,xbinsQcl);
  hQcl[1] = new TH1F("hQclY",";Q_{cl}",nbinsQcl,xbinsQcl);


  hQch[0] = new TH1F("hQchX",";Q_{ch}",nbinsQcl,xbinsQcl);
  hQch[1] = new TH1F("hQchY",";Q_{ch}",nbinsQcl,xbinsQcl);

  hQfirst[0] = new TH2F("hQfirstX",";Q_{first}",nbinsQcl,xbinsQcl,4096,0,4096);
  hQfirst[1] = new TH2F("hQfirstY",";Q_{first}",nbinsQcl,xbinsQcl,4096,0,4096);
  hQlast[0] = new TH2F("hQlastX",";Q_{last}",nbinsQcl,xbinsQcl,4096,0,4096);
  hQlast[1] = new TH2F("hQlastY",";Q_{last}",nbinsQcl,xbinsQcl,4096,0,4096);

  hRatioVsQchX = new TH2F("hRatioVsQchX",";Q_{channelX}^{max};Q_{x}/Q_{y}",nbinsQcl,xbinsQcl,500,0,5);
  hRatioVsQchY = new TH2F("hRatioVsQchY",";Q_{channelY}^{max};Q_{x}/Q_{y}",nbinsQcl,xbinsQcl,500,0,5);


  const Int_t nbinsRatio = 150;
  Double_t xminRatio = 0;
  Double_t xmaxRatio = 3;
  Double_t xbinsRatio[nbinsRatio+1];
  for (Int_t i=1;i<=nbinsRatio;i++)
    xbinsRatio[i] = xminRatio + i*1./nbinsRatio*(xmaxRatio - xminRatio);


  hRatioVsQchXVsQtotX = new TH3F("hRatioVsQchXVsQtotX",";Q_{channelX}^{max};Q_{X}^{tot};Q_{x}/Q_{y}",nbinsQcl2,xbinsQcl2,nbinsQcl2,xbinsQcl2,nbinsRatio,xbinsRatio);
  hRatioVsQchYVsQtotY = new TH3F("hRatioVsQchYVsQtotY",";Q_{channelY}^{max};Q_{Y}^{tot};Q_{x}/Q_{y}",nbinsQcl2,xbinsQcl2,nbinsQcl2,xbinsQcl2,nbinsRatio,xbinsRatio);
 }

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

   /* TFile *hf = */ OpenHistFile("PedestalShift");

   // Save histograms here
   hMapRaw[0]->Write();
   hMapRaw[1]->Write();
   hQraw[0]->Write();
   hQraw[1]->Write();
   hQsat[0]->Write();
   hQsat[1]->Write();
   hQnosat[0]->Write();
   hQnosat[1]->Write();
   hQsat1[0]->Write();
   hQsat1[1]->Write();
   hQnosat1[0]->Write();
   hQnosat1[1]->Write();
   hQcl[0]->Write();
   hQcl[1]->Write();
   hQch[0]->Write();
   hQch[1]->Write();
   hQfirst[0]->Write();
   hQfirst[1]->Write();
   hQlast[0]->Write();
   hQlast[1]->Write();
   hPmm2->Write();
   hAngle->Write();
   hRatioVsQchX->Write();
   hRatioVsQchY->Write();
   hRatioVsQchXVsQtotX->Write();
   hRatioVsQchYVsQtotY->Write();
 
 }


void HarpoAnalysePedestalShift::DisplayAnalysis(TRootEmbeddedCanvas* ecTab, TGListBox* infobox)
{
  // in hrecomonitorgui: fills the canvas on the "Analysis" tab
  TCanvas* c = ecTab->GetCanvas();
  c->GetPad(1)->Delete();
  c->GetPad(2)->Delete();
  c->Divide(2);

  for(Int_t plane = 0; plane<2;plane++){
    TVirtualPad* cMap = c->GetPad(plane+1);
    cMap->Divide(1,1);
    hQch[plane]->SetMinimum(0.5);
    // MakeNiceHisto(hQch[plane],cMap->GetPad(2),"");
    // cMap->GetPad(2)->SetLogx();
    // cMap->GetPad(2)->SetLogy();
    hMapRaw[plane]->SetMaximum(100000);
    MakeNiceHisto(hMapRaw[plane],cMap->GetPad(1),"colz");
    Double_t qmax = 0;
    infobox->NewEntry(Form("Qchmax = %g",qmax));
  }

  c->Update();
}


void HarpoAnalysePedestalShift::ConfigFrame(TGMainFrame* /* fMain */, Int_t /* id */)
{
  // // in hrecomonitorgui: Creates a popup window for analysis configuration
  // // You must define SetConfig() properly

  // UInt_t xsize = 200;
  // UInt_t ysize2 = 20;
  // UInt_t ysize = 10*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* frame = new TGVerticalFrame(main,xsize,ysize,kVerticalFrame);

  // // Object layout options
  // 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);

  // //________ DO NOT MODIFY ABOVE _____________________________




  // // // Title of the analysis
  // // TGLabel* fAnalysisLabel = new TGLabel(frame, "Template analysis");

  // // // Create a line for choosing value of parameter
  // // TGHorizontalFrame* nBinsFrame = new TGHorizontalFrame(frame,xsize,ysize2,kHorizontalFrame);
  // // TGLabel* nBinsLabel = new TGLabel(nBinsFrame, "fNbins");
  // // fChooseNbins    = new TGNumberEntry(nBinsFrame);
  // // fChooseNbins->SetNumber(fNbins);

  // // main->AddFrame(frame,fLayout4);
  // // frame->AddFrame(fAnalysisLabel,fLayout3);
  // // frame->AddFrame(nBinsFrame,fLayout3);
  // // nBinsFrame->AddFrame(nBinsLabel,fLayout1);
  // // nBinsFrame->AddFrame(fChooseNbins,fLayout2);



  // //________ DO NOT MODIFY BELOW _____________________________
  // // Button to validate configuration
  // TGTextButton* setConf = new TGTextButton(frame,"Save Config",id);
  // setConf->Associate(fMain);

  // frame->AddFrame(setConf,fLayout3);

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



void   HarpoAnalysePedestalShift::SetConfig()
{
  // Update the configuration according to the values in the popup window

  //  if(!fChooseNbins) return;

  //  fNbins = fChooseNbins->GetNumber();
}