HarpoAnalyseMatching.cxx

Go to the documentation of this file.

// 
// File HarpoAnalyseMatching.cxx
//
#include "HarpoAnalyseMatching.h"
#include "HarpoConfig.h"
#include "HarpoDetSet.h"
#include "HarpoDebug.h"
#include "HarpoDccEvent.h"
#include "Pmm2Event.h"
#include "HarpoEvent.h"
#include "HarpoRecoEvent.h"
#include "HarpoSimEvent.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 <cstdlib>
#include <cstring>
#include <cassert>
#include <fstream>
#include <iostream>

ClassImp(HarpoAnalyseMatching);

void   HarpoAnalyseMatching::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   HarpoAnalyseMatching::process()
{
  //  Bool_t drawEvent = kFALSE;
  nEvents++;
  if(nEvents%500==0) 
    Info("process","%ld events",nEvents);  

  // Process RAW data
  if (!gHDetSet->isExist(XDCC)) return;
  if (!gHDetSet->isExist(YDCC)) return;
  HarpoDccMap* mX = fEvt->GetDccMap(XDCC);
  HarpoDccMap* mY = fEvt->GetDccMap(YDCC);
  if ( mX == NULL ) return;
  if ( mY== NULL ) return;



  Double_t QchXmax = 0, QchYmax = 0;
  //  Int_t satX[NCHAN], satY[NCHAN];
  Int_t nSatTX[NADC], nSatTY[NADC];
  for(Int_t i=0;i<NADC;i++){ // Time bins
    nSatTX[i] = 0;
    nSatTY[i] = 0;
  }
  Int_t nSatX = 0, nSatY = 0, nChX = 0, nChY = 0;
  for(Int_t j=0;j<NCHAN;j++){ // Channels
    Double_t qchX = 0, qchY = 0;
    for(Int_t i=0;i<NADC;i++){ // Time bins
      Double_t qX = mX->GetData(j,i);
      if(qX > -1000)
        qchX += qX;
      Double_t qY = mY->GetData(j,i);
      if(qY > -1000)
        qchY += qY;
      if(qchX>fLimSat)
        nSatTX[i]++;
      if(qchY>fLimSat)
        nSatTY[i]++;
    }
    if(qchX>QchXmax) QchXmax = qchX;
    if(qchY>QchYmax) QchYmax = qchY;

    hQch->Fill(qchX,qchY);

    if(qchX>fLimSat){
      //      satX[j] = 1;
      nSatX++;
    }// else satX[j] = 0;
    if(qchY>fLimSat){
      //      satY[j] = 1;
      nSatY++;
    }// else satY[j] = 0;

    if(qchX) nChX++;
    if(qchY) nChY++;


  }

  //  std::cout << nSatX << "  " << nSatY << "  " << hNsat << std::endl;
  hNsat->Fill(nSatX,nSatY);

  for(Int_t i=0;i<NADC;i++){ // Time bins
    if(!(nSatTX[i]==0 && nSatTY[i]==0))
      hNsat2->Fill(nSatTX[i],nSatTY[i]);
  }
  
  Double_t qXtot = 0, qYtot = 0;
  Double_t qxtot[NCHAN];
  Int_t xstart[NCHAN];
  for(Int_t j=0;j<NCHAN;j++){
    qxtot[j] = 0;
    xstart[j] = 0;
  }
  Int_t n = 0;
  Double_t meanRatio = 0, sigmaRatio = 0;
  TArrayD* arrRatio = new TArrayD(NADC);
  for(Int_t i=0;i<NADC;i++){ // Time bins
    arrRatio->AddAt(-1,i);
    Double_t qXch = 0, qYch = 0;
    for(Int_t j=0;j<NCHAN;j++){ // Channels
      Double_t qX = mX->GetData(j,i);
      Double_t qY = mY->GetData(j,i);
      if(qX > -1000){
        qXch += qX;
        qxtot[j] += qX;
        xstart[j] = 1;
      }else{
        if(xstart[j] == 1)
          xstart[j] = 2;
      }
      if(qY > -1000) qYch += qY;
    }
    qXtot += qXch;
    qYtot += qYch;

    if(qYch>0 && nSatTX[i] <=1 && nSatTY[i] == 0){
      for(Int_t j=0;j<NCHAN;j++){
        Double_t qcorr = qYtot + qXtot - qxtot[j];
        // for(Int_t j2=0;j2<NCHAN;j2++){
        //      if(j2 != j) qcorr -= qxtot[j2];
        // }
        //      if(mX->GetData(j,i)>-100 && qYch>0 && nSatX <=1)
        if(xstart[j] != 2){
          hQxVsQyLog->Fill(qcorr,qxtot[j]);
          hQxVsQy->Fill(qcorr,qxtot[j]);
        }
      }
    }

    if(nSatX<5 && nSatY<5 && (qXch>0 || qYch>0)){
      //     std::cout << i << "  " << nSatX << ", " << nSatY << " => " << qXtot << ", " << qYtot << std::endl; 
      hSaturation[nSatX][nSatY]->Fill(qXtot,qYtot);
    }

    if(nSatTX[i]<5 && nSatTY[i]<5 && (qXch>0 || qYch>0)){
      //     std::cout << i << "  " << nSatX << ", " << nSatY << " => " << qXtot << ", " << qYtot << std::endl; 
      hSaturation2[nSatTX[i]][nSatTY[i]]->Fill(qXtot,qYtot);
    }

    if(qXch <= 0) continue;
    Double_t ratio = qYch/qXch/0.92;
    hRatio->Fill(ratio);

    if(ratio>0)
    // hRatioVsZ->Fill(i,TMath::Log(ratio));
    // hRatioVsNchX->Fill(nChX,TMath::Log(ratio));
    // hRatioVsNchY->Fill(nChY,TMath::Log(ratio));

    ratio = TMath::Log(ratio);
    hRatioVsZ->Fill(i,ratio);
    hRatioVsNchX->Fill(nChX,ratio);
    hRatioVsNchY->Fill(nChY,ratio);

    meanRatio += ratio;
    sigmaRatio += ratio*ratio;
    arrRatio->AddAt(ratio,i);
    n++;


  }


  if(n){
    meanRatio /= n;
    sigmaRatio /= n;
    sigmaRatio -= meanRatio*meanRatio;
      hMeanVsSigmaRatio->Fill(TMath::Sqrt(sigmaRatio),meanRatio);
    //if(meanRatio>0)
    //      hMeanVsSigmaRatio->Fill(TMath::Sqrt(sigmaRatio),TMath::Log(meanRatio));
  }
  Double_t truncM, truncS;
  TruncMean(arrRatio, truncS, truncM, 0.1, 0.9);
  delete arrRatio;
  if(truncM>0)
    hTruncMeanVsSigmaRatio->Fill(truncS,truncM);
  //    hTruncMeanVsSigmaRatio->Fill(truncS,TMath::Log(truncM));
  
  if(qXtot){
    hTest->Fill(QchXmax,QchYmax);
    hRatioTot->Fill(n,TMath::Log(qYtot/qXtot));
    if(QchXmax>fLimSat || QchYmax>fLimSat)
      hRatioTotCut1->Fill(n,TMath::Log(qYtot/qXtot));
    else
      hRatioTotCut2->Fill(n,TMath::Log(qYtot/qXtot));
  }
}

void   HarpoAnalyseMatching::Init()
{
  
  fLimSat = 5e4;


  // Initialise histograms here
  hRatioVsZ = new TH2F("hRatioVsZ",";Z",512,0,512,200,-5,5);
  hRatioVsNchX = new TH2F("hRatioVsNchX",";N_{chX}",288,0,288,200,-5,5);
  hRatioVsNchY = new TH2F("hRatioVsNchY",";N_{chY}",288,0,288,200,-5,5);
  hMeanVsSigmaRatio = new TH2F("hMeanVsSigmaRatio","",200,0,5,200,-5,5);
  hTruncMeanVsSigmaRatio = new TH2F("hTruncMeanVsSigmaRatio","",200,0,5,200,-5,5);
  hRatioTot = new TH2F("hRatioTot","",512,0,512,2000,-5.,5.);
  hRatioTotCut1 = new TH2F("hRatioTotCut1","",512,0,512,2000,-5.,5.);
  hRatioTotCut2 = new TH2F("hRatioTotCut2","",512,0,512,2000,-5.,5.);
  hTest = new TH2F("hTest","",2,0,2*fLimSat,2,0,2*fLimSat);


  const Int_t nbinsQcl = 1000;
  Double_t xminQcl = 10;
  Double_t xmaxQcl = 1e6;
  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);
  hQch = new TH2F("hQch",";Q_{ch X};Q_{ch Y}",nbinsQcl,xbinsQcl,nbinsQcl,xbinsQcl);
  for(Int_t i = 0; i<5; i++){
    for(Int_t j = 0; j<5; j++){
      hSaturation[i][j] = new TH2F(Form("hSaturation%d_%d",i,j),";Q_{X};Q_{Y}",nbinsQcl,xbinsQcl,nbinsQcl,xbinsQcl);
      hSaturation2[i][j] = new TH2F(Form("hSaturation2%d_%d",i,j),";Q_{X};Q_{Y}",nbinsQcl,xbinsQcl,nbinsQcl,xbinsQcl);
    }
  }
  hQxVsQyLogLog = new TH2F("hQxVsQyLog",";Q_{X};Q_{Y}",nbinsQcl,xbinsQcl,nbinsQcl,xbinsQcl);
  hQxVsQyLog = new TH2F("hQxVsQyLog",";Q_{X};Q_{Y}",10000,0,1e6,nbinsQcl,xbinsQcl);
  hQxVsQy = new TH2F("hQxVsQy",";Q_{X};Q_{Y}",2000,0,1e6,1000,0,1e5);

  hNsat = new TH2F("hNsat",";N_{satX};N_{satY}",10,0,10,10,0,10);
  hNsat2 = new TH2F("hNsat2",";N_{satX};N_{satY}",10,0,10,10,0,10);



  hRatio = new TH1F("hRatio",";Q_{x}/Q_{y}",1000,0.5,1.5);

}

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

   /* TFile* hf = */ OpenHistFile("matching");
   // Save histograms here
   hRatioVsZ->Write();
   hRatioVsNchX->Write();
   hRatioVsNchY->Write();
   hMeanVsSigmaRatio->Write();
   hTruncMeanVsSigmaRatio->Write();
   hRatioTot->Write();
   hRatioTotCut1->Write();
   hRatioTotCut2->Write();
   hTest->Write();
   hQxVsQyLogLog->Write();
   hQxVsQyLog->Write();
   hQxVsQy->Write();


   hQch->Write();
   for(Int_t i = 0; i<5; i++){
     for(Int_t j = 0; j<5; j++){
       hSaturation[i][j]->Write();
     }
   }
   for(Int_t i = 0; i<5; i++){
     for(Int_t j = 0; j<5; j++){
       hSaturation2[i][j]->Write();
     }
   }
   hNsat->Write();
   hNsat2->Write();

   hRatio->Write();
 
 }




Double_t HarpoAnalyseMatching::TruncMean(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])<0&&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;

  //Info("TruncMean","%d %d (%d) => %g %g",tLow,tHigh,tHigh-tLow,truncMean,truncSigma);
  //  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;
}