HarpoAnalyseBaselineFluct.cxx

Go to the documentation of this file.

// 
// File HarpoAnalyseBaselineFluct.cxx
//
#include "HarpoAnalyseBaselineFluct.h"
#include "HarpoConfig.h"
#include "HarpoDetSet.h"
#include "HarpoDebug.h"
#include "HarpoDccEvent.h"
#include "HarpoEvent.h"
#include "HarpoRecoEvent.h"

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

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

ClassImp(HarpoAnalyseBaselineFluct);

void   HarpoAnalyseBaselineFluct::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   HarpoAnalyseBaselineFluct::process()
{
  //   Bool_t drawEvent = kFALSE;
  nEvents++;
  // Process Reconstructed data (clusters, tracks, matched tracks)
  HarpoRecoEvent* reco = fEvt->GetRecoEvent();
  
  //std::cout << "PROCESSING RAW DATA" << std::endl;
  Double_t xStart[2] = {0, 0};
  Double_t tStart[2] = {0, 0};

  Int_t nNoiseChan[2];
  // Process RAW data
  //  Double_t sigQ[NADC][2];
  for(UInt_t ndet=0;ndet<2;ndet++) {
    if (!gHDetSet->isExist(ndet)) continue;
    HarpoDccMap* m = fEvt->GetDccMap(ndet);
    if ( m == NULL ) continue;
    //    Double_t limNch = 30, limSig = 0.06, limDelta = 5;
    // Double_t limNch = 30, limSig = 20, limDelta = 5;
    // Double_t limTruncL = 0.1, limTruncH = 0.7;
    Double_t meanNoise = 0, qNoise = 0; 
    Int_t nChan = 0;
    TH2F* hNoise = new TH2F("hNoise","",511,0,511,4096,-300,3796);
    TH2F* hNoiseTot = new TH2F("hNoiseTot","",511,0,511,1000,-3000,1e5);
    TH2F* hSignal = new TH2F("hSignal","",511,0,511,4096,-300,3796);
    Double_t xNoiseMin = 511, xNoiseMax = 1, qNoiseMax = 1;
    nNoiseChan[ndet] = 0;
    Double_t qStart = 0;
    Int_t nbinsStart = 0;
    for(Int_t i=0;i<NADC;i++){ //Time bins

      Double_t qTot = 0;
      Double_t sigmaQ = 0;
      TArrayD* arr = new TArrayD(NCHAN);
      for(Int_t j=0;j<NCHAN;j++){ // Channels
        Double_t q = m->GetRawData(j,i);
        //if(q<400) m->SetRawData(j,i,-1); // FIXME
        //continue; // FIXME
        if(q < 0) continue; // NO DATA
        arr->AddAt(q,j);
        qTot += q;
        sigmaQ += q*q;
        nChan++;

        if(nbinsStart<4){
          xStart[ndet] += j*q;
          tStart[ndet] += i*q;
          qStart += q;
        }
      }
      if(qTot>0) nbinsStart++;

      // Calculate truncated mean and sigma for that time bin
      Double_t truncSig = 0, truncMean = 0;
      TruncSigma(arr,truncSig,truncMean,fLimTruncL,fLimTruncH);
      delete arr;

      if(truncMean){
        hSigma2[ndet]->Fill(truncSig/truncMean);
      }

      Double_t qNoiseT = 0;
      if(truncSig){
        //      Info("process","truncSig/truncMean = %f/%f = %f", truncSig, truncMean, truncSig/truncMean);
        if(nChan>fLimNch && truncSig<fLimSig) nNoiseChan[ndet]++;
        hSigma[ndet]->Fill(truncSig);
        for(Int_t j=0;j<NCHAN;j++){ // Channels
          Double_t q = m->GetRawData(j,i);
          Double_t q2 = m->GetData(j,i);
          if( q < 0 ) continue;  // NO DATA
          hSigmaVsMean[ndet]->Fill(truncSig,TMath::Abs(q-truncMean)/truncSig);
          //      if(nChan>fLimNch && truncSig/truncMean<fLimSig && TMath::Abs(q-truncMean)<fLimDelta*truncSig){
          if(nChan>fLimNch && truncSig<fLimSig && TMath::Abs(q-truncMean)<fLimDelta*truncSig){
            hSigmaVsMeanCut[ndet]->Fill(truncSig,TMath::Abs(q-truncMean)/truncSig);
            meanNoise += i*q2;
            qNoise += q2;
            qNoiseT += q2;
            //Info("process","%g",q);
            hNoise->Fill(i,q2);
            if(xNoiseMin>i) xNoiseMin = i;
            if(xNoiseMax<i) xNoiseMax = i;
            if(qNoiseMax<q2) qNoiseMax = q2;
            q = 0;
            //      Info("process","Remove pixel %d %d",j,i);
            m->SetRawData(j,i,-1);
          }else{
            hSignal->Fill(i,q2);
          }
        }
        hNoiseTot->Fill(i,qNoiseT);
      }


    }

    if(qStart>0){
      xStart[ndet] /= qStart;
      tStart[ndet] /= qStart;
    }

    //    if(nNoiseChan[ndet] > 4){
    if(qNoise > 250){
      TProfile* pNoise = hNoise->ProfileX(Form("pNoise%d",ndet));
      
      Double_t mean = 0, qTot = 0;
      Int_t binmax = pNoise->GetMaximumBin();
      for(Int_t bin = binmax-2; bin<binmax+3; bin++){
        Double_t q = pNoise->GetBinContent(bin);
        mean += bin*q;
        qTot += q;
      }
      pNoise->Delete();
   
      Double_t qSignal = 0;
      if(qTot>0){
        for(Int_t i = 0; i<511;i++){
          Double_t qSigTot = 0, qSigMax = 0;
          for(Int_t j = 0; j<4096; j++){
            hNoiseShape[ndet]->Fill(i-mean/qTot,hNoise->GetYaxis()->GetBinCenter(j+1),hNoise->GetBinContent(i+1,j+1));
            hNoiseShapeTot[ndet]->Fill(i-mean/qTot,hNoiseTot->GetYaxis()->GetBinCenter(j+1),hNoiseTot->GetBinContent(i+1,j+1));
            hNoiseShapeNorm[ndet]->Fill(i-mean/qTot,hNoise->GetYaxis()->GetBinCenter(j+1)/TMath::Abs(qTot),hNoise->GetBinContent(i+1,j+1));
            //      hSignalShape[ndet]->Fill(i-mean/qTot,hNoise->GetYaxis()->GetBinCenter(j+1),hSignal->GetBinContent(i+1,j+1));
            Double_t qSig = hNoise->GetYaxis()->GetBinCenter(j+1);
            qSigTot += qSig*hSignal->GetBinContent(i+1,j+1);
            if(qSig>qSigMax && hSignal->GetBinContent(i+1,j+1)>0) qSigMax = qSig;
          }
          hSignalShape[ndet]->Fill(i-mean/qTot,qSigTot*0.01);
          hSignalShape2[ndet]->Fill(i-mean/qTot,qSigMax*0.5);
          if(TMath::Abs(i-mean/qTot+13)<10) qSignal += qSigTot;
        }
      }
      hSignalVsNoise[ndet]->Fill(qNoise,qSignal);
    }

    hNoise->Delete();
    hNoiseTot->Delete();
    hSignal->Delete();

   
    for(Int_t j=0;j<NCHAN;j++){ // Channels
      Int_t nPix = 0;
      for(Int_t i=0;i<NADC;i++){ //Time bins
        Double_t q = m->GetRawData(j,i);
        if(q < 0 && nPix == 0) continue;
        if(q < 0){
          if(nPix<10){
            if(gHarpoDebug>1)
              Info("process","Remove small cluster on channel %d, t = %d",j,i);
            for(Int_t iTmp = i - nPix; iTmp<i; iTmp++)
              m->SetRawData(j,iTmp,-1);
          }
          nPix = 0;
        }else nPix++;
      }
    }
  }
  
  reco->SetXstart(xStart[0]);
  reco->SetYstart(xStart[1]);
  reco->SetTstart(0.5*(tStart[0]+tStart[1]));

  if(nEvents%100 == 0)
    std::cout << "Event " << nEvents << " done" << std::endl;
}

 void   HarpoAnalyseBaselineFluct::Init()
 {

   fLimNch = 30;
   fLimSig = 20;
   fLimDelta = 5;
   fLimTruncL = 0.1;
   fLimTruncH = 0.7;

   const Int_t nbinsQ = 100;
   Double_t xminQ = 250;
   Double_t xmaxQ = 1e7;
   Double_t logxminQ = TMath::Log(xminQ);
   Double_t logxmaxQ = TMath::Log(xmaxQ);
   Double_t binwidthQ = (logxmaxQ-logxminQ)/nbinsQ;
   Double_t xbinsQ[nbinsQ+1];
   xbinsQ[0] = xminQ;
   for (Int_t i=1;i<=nbinsQ;i++)
     xbinsQ[i] = TMath::Exp(logxminQ+i*binwidthQ);
   hSignalVsNoise[0] = new TH2F("hSignalVsNoiseX",";Q_{noise} [ADC];Q_{signal} [ADC]",nbinsQ,xbinsQ,nbinsQ,xbinsQ);
   hSignalVsNoise[1] = new TH2F("hSignalVsNoiseY",";Q_{noise} [ADC];Q_{signal} [ADC]",nbinsQ,xbinsQ,nbinsQ,xbinsQ);
  
   hNoiseShape[0] = new TH2F("hNoiseShapeX",";t-t_{noise} [30ns];q [ADCch]",1000,-250,250,4096,-300,3796);
   hNoiseShape[1] = new TH2F("hNoiseShapeY",";t-t_{noise} [30ns];q [ADCch]",1000,-250,250,4096,-300,3796);
   hNoiseShapeTot[0] = new TH2F("hNoiseShapeTotX",";t-t_{noise} [30ns];q [ADCch]",1000,-250,250,1000,-3000,1e5);
   hNoiseShapeTot[1] = new TH2F("hNoiseShapeTotY",";t-t_{noise} [30ns];q [ADCch]",1000,-250,250,1000,-3000,1e5);
   hNoiseShapeNorm[0] = new TH2F("hNoiseShapeNormX",";t-t_{noise} [30ns];q_{norm}",1000,-250,250,500,-1,1);
   hNoiseShapeNorm[1] = new TH2F("hNoiseShapeNormY",";t-t_{noise} [30ns];q_{norm}",1000,-250,250,500,-1,1);
   hSignalShape[0] = new TH2F("hSignalShapeX",";t-t_{noise} [30ns];q_{tot} [ADCch]",1000,-250,250,4096,-300,3796);
   hSignalShape[1] = new TH2F("hSignalShapeY",";t-t_{noise} [30ns];q_{tot} [ADCch]",1000,-250,250,4096,-300,3796);
   hSignalShape2[0] = new TH2F("hSignalShape2X",";t-t_{noise} [30ns];q_{max} [ADCch]",1000,-250,250,4096,-300,3796);
   hSignalShape2[1] = new TH2F("hSignalShape2Y",";t-t_{noise} [30ns];q_{max} [ADCch]",1000,-250,250,4096,-300,3796);

   hSigma[0] = new TH1F("hSigmaX",";#sigma_{trunc}",500,0,200);
   hSigma[1] = new TH1F("hSigmaY",";#sigma_{trunc}",500,0,200);
   hSigmaVsMean[0] = new TH2F("hSigmaVsMeanX",";#sigma_{trunc};(q-#mu_{trunc})/#sigma_{trunc}",500,0,200,500,0,50);
   hSigmaVsMean[1] = new TH2F("hSigmaVsMeanY",";#sigma_{trunc};(q-#mu_{trunc})/#sigma_{trunc}",500,0,200,500,0,50);
   hSigmaVsMeanCut[0] = new TH2F("hSigmaVsMeanCutX",";#sigma_{trunc};(q-#mu_{trunc})/#sigma_{trunc}",500,0,200,500,0,50);
   hSigmaVsMeanCut[1] = new TH2F("hSigmaVsMeanCutY",";#sigma_{trunc};(q-#mu_{trunc})/#sigma_{trunc}",500,0,200,500,0,50);
   hSigmaVsTruncH[0] = new TH2F("hSigmaVsTruncHX",";#sigma_{trunc};truncH",500,0,200,20,0,100);
   hSigmaVsTruncH[1] = new TH2F("hSigmaVsTruncHY",";#sigma_{trunc};truncH",500,0,200,20,0,100);
   hSigma2[0] = new TH1F("hSigma2X",";#sigma_{trunc}/#mu_{trunc}",500,0,5);
   hSigma2[1] = new TH1F("hSigma2Y",";#sigma_{trunc}/#mu_{trunc}",500,0,5);

 }

void   HarpoAnalyseBaselineFluct::Save(char * /* mode */)
{
 
   OpenHistFile("baseline");

  
   // Save Run header (contains all DB info about the run)
   // if(fRunHeader != NULL)
   //   fRunHeader->Write("fRunHeader");
  hSigma[0]->Write();
  hSigma[1]->Write();
  hSigmaVsMean[0]->Write();
  hSigmaVsMean[1]->Write();
  hSigmaVsMeanCut[0]->Write();
  hSigmaVsMeanCut[1]->Write();
  hSigma2[0]->Write();
  hSigma2[1]->Write();
  
  hNoiseShape[0]->Write();
  hNoiseShape[1]->Write();
  hNoiseShapeTot[0]->Write();
  hNoiseShapeTot[1]->Write();
  hNoiseShapeNorm[0]->Write();
  hNoiseShapeNorm[1]->Write();
  hSignalShape[0]->Write();
  hSignalShape[1]->Write();
  hSignalShape2[0]->Write();
  hSignalShape2[1]->Write();

  hSigmaVsTruncH[0]->Write();
  hSigmaVsTruncH[1]->Write();

  hSignalVsNoise[0]->Write(); 
  hSignalVsNoise[1]->Write(); 
  
}



Double_t HarpoAnalyseBaselineFluct::TruncSigma(TArrayD* vect, Double_t &truncS, Double_t &truncM, 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();
  //  Info("TruncSigma","size = %d",size);
  Double_t truncMean = 0;
  Double_t truncSigma = 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]);
    truncSigma += vect->At(index[tind])*vect->At(index[tind]);
  }

  delete[] index;

  //  if(tHigh-tLow == 0) return 0;
  if(tHigh-tLow < 15) return 0;

  truncMean /= (tHigh-tLow);
  truncSigma /= (tHigh-tLow);
  truncSigma -= truncMean*truncMean;

  truncM = truncMean;
  truncS = TMath::Sqrt(truncSigma);

  return truncSigma;
}