HarpoAnalyseVertex.cxx

Go to the documentation of this file.

// 
// File HarpoAnalyseVertex.cxx
//
#include "HarpoAnalyseVertex.h"
#include "HarpoConfig.h"
#include "HarpoDetSet.h"
#include "HarpoDebug.h"
#include "HarpoDccEvent.h"
#include "Pmm2Event.h"
#include "HarpoEvent.h"
#include "HarpoRecoEvent.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(HarpoAnalyseVertex);

void   HarpoAnalyseVertex::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   HarpoAnalyseVertex::process()
{
  //  Bool_t drawEvent = kFALSE;
  nEvents++;
  if (!gHDetSet->isExist(XDCC)) return;
  if (!gHDetSet->isExist(YDCC)) return;
  // Process RAW data
  HarpoDccMap* mX = fEvt->GetDccMap(XDCC);
  HarpoDccMap* mY = fEvt->GetDccMap(YDCC);
  if ( mX == NULL ) return;
  if ( mY == NULL ) return;

  Int_t i2 = 0;
  Double_t qTotX3 = 0, qTotY3 = 0;
  Double_t meanX3 = 0, meanY3 = 0;
  Double_t meanXevt = 0, meanYevt = 0;

  Int_t triggerType = -1;
  if(gHDetSet->isExist(PMM2)) {
    Pmm2Event *anaEvt = static_cast<Pmm2Event *>(fEvt->GetDetEvent(PMM2)); 
    triggerType = anaEvt->GetTriggerType();
  }
  Float_t data[7];
  for(Int_t i = 0; i<7; i++) data[i] = -1000;
  data[0] = triggerType;

  for(Int_t i=0;i<NADC;i++){ //Time bins
    Double_t qTotX = 0, qTotY = 0;
    Double_t meanX = 0, meanY = 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 > 0){
        qTotX += qX;
        meanX += j*qX;
        qTotX3 += qX;
        meanX3 += j*qX;
      }
      if(qY > 0){
        qTotY += qY;
        meanY += j*qY;
        qTotY3 += qY;
        meanY3 += j*qY;
      }
    }
    if(qTotX*qTotY == 0) continue;
    meanX /= qTotX;
    meanY /= qTotY;
    if(i2 == 2){
      meanXevt = meanX3/qTotX3;
      meanYevt = meanY3/qTotY3;
      //gVertex->SetPoint(gVertex->GetN(),i-1,meanXevt,meanYevt);
      data[1] = i-1;
      data[2] = meanXevt;
      data[3] = meanYevt;
      if(i>=150) hVertex2->Fill(meanXevt,meanYevt);
      if(TMath::Abs(meanXevt-144)>10) break;
      if(TMath::Abs(meanYevt-144)>10) break;
    }
    i2++;
    if(i<150 && i2>2) break;
    //    if(i2>kNhist) break;
    if(i2 == 3){
      //      gVertex2->SetPoint(gVertex2->GetN(),i-1,meanXevt,meanYevt);
      data[4] = i-1;
      data[5] = meanXevt;
      data[6] = meanYevt;

      Double_t xMin = 400, yMin = 400, xMax = -1, yMax = -1;
      for(Int_t J=0;J<NCHAN;J++){ // Channels
        Double_t qXch = 0, qYch = 0;
        for(Int_t I=i-1;I<i+100 && I<NADC;I++){ //Time bins
          Double_t qX = mX->GetData(J,I);
          Double_t qY = mY->GetData(J,I);
          if(qX > 0)
            qXch += qX;
          if(qY > 0)
            qYch += qY;
        }
        if(qXch>1000 && J<xMin) xMin = J;
        if(qYch>1000 && J<yMin) yMin = J;
        if(qXch>1000 && J>xMax) xMax = J;
        if(qYch>1000 && J>yMax) yMax = J;
      }
      //std::cout << xMin << " " << yMin << " " << xMax << " " <<yMax << std::endl;
      Double_t widthX = xMax-xMin+1;
      Double_t widthY = yMax-yMin+1;
      hWidthXvsY->Fill(widthX,widthY);
      if(widthX > 1 && widthY > 1){
        hWidthRatio->Fill(widthX/widthY);
        hWidthRatioLog->Fill(TMath::Log(widthX/widthY));
      }
    }
    // for(Int_t jX=0;jX<NCHAN;jX++){ // Channels
    //   Double_t qX = mX->GetData(jX,i);
    //   if(qX <= 0) continue;
    //   qX /= qTotX;
    //   for(Int_t jY=0;jY<NCHAN;jY++){ // Channels
    //  Double_t qY = mY->GetData(jY,i);
    //  if(qY <= 0) continue;
    //  qY /= qTotY;
    //  if(i2-1<kNhist)
    //    hVertex[i2-1]->Fill(jX-meanX,jY-meanY,qX*qY);
    //   }
    // }


  }


  fNtuple->Fill(data);


}

 void   HarpoAnalyseVertex::Init()
 {
   
   fNtuple = new TNtuple("fNtuple","vertex data","trigger:T:X:Y:T2:X2:Y2");
   fNtuple->SetDirectory(0);

   hVertex2 = new TH2F("hVertex",";X;Y",576,0,288,576,0,288);

   hWidthRatio = new TH1F("hWidthRatio","",1000,0,10);
   hWidthRatioLog = new TH1F("hWidthRatioLog","",1000,-2,2);

   hWidthXvsY = new TH2F("hWidthXvsY",";w_{X};w_{Y};",288,0,288,288,0,288);
 }

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

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

   // Save histograms here
   hVertex2->Write();


   fNtuple->Write("fNtuple");
   hWidthRatio->Write();
   hWidthRatioLog->Write();
   hWidthXvsY->Write();

 
 }