HarpoSelectorVertex.cxx

Go to the documentation of this file.

// 
// File HarpoSelectorVertex.cxx
//
#include "HarpoSelectorVertex.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 "HarpoRecoMonitorGui.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 "TVector3.h"

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

using namespace std;


ClassImp(HarpoSelectorVertex);

void   HarpoSelectorVertex::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   HarpoSelectorVertex::process()
{
  //  Bool_t drawEvent = kFALSE;
  nEvents++;
  Int_t test = eventselection();

  fEvt->GetHeader()->SetEvAnaStatus(test);
  Info("process","********** EvAnaStatus = %d *********",test);

  return;
}

Int_t   HarpoSelectorVertex::eventselection()
{
  if (!gHDetSet->isExist(XDCC)) return 0;
  if (!gHDetSet->isExist(YDCC)) return 0;
  if(!fEvt->GetRecoEvent()) return 0;


  TVector3 pgamma = TVector3(fDx,fDy,1).Unit();
  TVector3 ppola = (TVector3(0,1,0).Cross(pgamma)).Unit();

  HarpoRecoEvent* reco = fEvt->GetRecoEvent();

  //Double_t xStart[2] = {-1,-1};
  Double_t zStart[2] = {1000,1000};
  TClonesArray* clArray = reco->GetClustersArray();
  Int_t nCl = clArray->GetEntries();
  for(Int_t iCl = 0; iCl<nCl; iCl++){
    HarpoRecoClusters* cl = (HarpoRecoClusters*)clArray->At(iCl);
    Int_t plane = cl->GetPlane();
    if(cl->GetZ()<zStart[plane]){
      //  xStart[plane] = cl->GetX();
      zStart[plane] = cl->GetZ();
    }
  }

  // if(gHDetSet->isExist(PMM2)){
  //   Pmm2Event *anaEvt = static_cast<Pmm2Event *>(fEvt->GetDetEvent(PMM2));
  // }

  Double_t fTstart = 90;
  Double_t fTend = 430;
  Double_t fThr = 1e4;
  for(Int_t plane = 0; plane<2; plane++){
    if (!gHDetSet->isExist(plane)) continue;
    // Process RAW data
    HarpoDccMap* m = fEvt->GetDccMap(plane);
    if ( m == NULL ) continue;
    //Int_t i2 = 0;
    Double_t qTotIn = 0;
    Double_t qTotBefore = 0;
    Double_t qTotAfter = 0;
    for(Int_t j=0;j<NCHAN;j++){ // Channels
      for(Int_t i=0;i<NADC;i++){ //Time bins
        Double_t q = m->GetData(j,i);
        if(q<0) continue;
        if(i< fTstart            ) qTotBefore += q;
        if(              i >fTend) qTotAfter += q;
        if(i>=fTstart && i<=fTend) qTotIn += q;
      }
    }
    if(qTotBefore>fThr || qTotAfter>fThr)
      return -2;
  }


  TClonesArray* vArray = reco->GetVertex3DArray();
  Int_t nV = vArray->GetEntries();
  for(Int_t iV = 0; iV<nV; iV++){
    HarpoRecoVertex3D* vertex = (HarpoRecoVertex3D*)vArray->At(iV);

    TVector3 p1 = vertex->GetP1Vertex3D();
    TVector3 p2 = vertex->GetP2Vertex3D();
    Double_t x0 = vertex->GetXvertex3D();
    Double_t y0 = vertex->GetYvertex3D();
    Double_t z0 = vertex->GetZvertex3D();

    Double_t u1u2 = (p1.Unit()).Dot(p2.Unit());
    Double_t thetapm=TMath::ACos(u1u2);

    TVector3 ppola2 = ppola.Cross(pgamma);
    TVector3 pPlus = TVector3(p1.Dot(pgamma.Unit()),p1.Dot(ppola.Unit()),p1.Dot(ppola2.Unit()));
    TVector3 pMinus = TVector3(p2.Dot(pgamma.Unit()),p2.Dot(ppola.Unit()),p2.Dot(ppola2.Unit()));
    Double_t omegapm = atan2((pPlus.Px()*pMinus.Py() - pMinus.Px()*pPlus.Py()),
                                 (pPlus.Px()*pMinus.Pz() - pMinus.Px()*pPlus.Pz()));
    
    TVector3 psum = p1 + p2;
    if(!(z0>120 && z0<350 && !(z0>255 && z0<270) && x0>5 && x0<283 && y0>5 && y0<283))
      continue;
    if(TMath::Abs(x0-142-fDx*z0)>20) continue;
    if(TMath::Abs(y0-150-fDy*z0)>20) continue;
    Info("eventselection","X: %g, Y: %g, Z: %g",x0,y0,z0);
    if(thetapm == 0)
      Info("eventselection","omega (thetapm == 0): %g %g %g",x0,y0,z0);

    if(thetapm>=fThetaMax) continue;
    if(thetapm<=fThetaMin) continue;

    Info("eventselection","omega: %g, theta: %g (%g, %g, %g)",omegapm, thetapm,p1.Mag(),p2.Mag(),u1u2);

    return 0;
  }
 

  return -2;


}

void   HarpoSelectorVertex::Init()
{
  
  // Initialise histograms here

  fThetaMin = 0;
  fThetaMax = 0.2;
  fDx = 0;
  fDy = 0;
  fChooseDx = 0;

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

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

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




 }

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

   //   OpenHistFile("eventselector");

   // Save histograms here

 }


void HarpoSelectorVertex::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* DxFrame = new TGHorizontalFrame(frame,xsize,ysize2,kHorizontalFrame);
  TGLabel* DxLabel = new TGLabel(DxFrame, "fDx");
  fChooseDx    = new TGNumberEntry(DxFrame);
  fChooseDx->SetNumber(fDx);
  TGHorizontalFrame* DyFrame = new TGHorizontalFrame(frame,xsize,ysize2,kHorizontalFrame);
  TGLabel* DyLabel = new TGLabel(DyFrame, "fDy");
  fChooseDy    = new TGNumberEntry(DyFrame);
  fChooseDy->SetNumber(fDy);
  TGHorizontalFrame* ThetaMinFrame = new TGHorizontalFrame(frame,xsize,ysize2,kHorizontalFrame);
  TGLabel* ThetaMinLabel = new TGLabel(ThetaMinFrame, "fThetaMin");
  fChooseThetaMin    = new TGNumberEntry(ThetaMinFrame);
  fChooseThetaMin->SetNumber(fThetaMin);
  TGHorizontalFrame* ThetaMaxFrame = new TGHorizontalFrame(frame,xsize,ysize2,kHorizontalFrame);
  TGLabel* ThetaMaxLabel = new TGLabel(ThetaMaxFrame, "fThetaMax");
  fChooseThetaMax    = new TGNumberEntry(ThetaMaxFrame);
  fChooseThetaMax->SetNumber(fThetaMax);

  main->AddFrame(frame,fLayout4);
  frame->AddFrame(fAnalysisLabel,fLayout3);
  frame->AddFrame(DxFrame,fLayout3);
  DxFrame->AddFrame(DxLabel,fLayout1);
  DxFrame->AddFrame(fChooseDx,fLayout2);
  frame->AddFrame(DyFrame,fLayout3);
  DyFrame->AddFrame(DyLabel,fLayout1);
  DyFrame->AddFrame(fChooseDy,fLayout2);
  frame->AddFrame(ThetaMinFrame,fLayout3);
  ThetaMinFrame->AddFrame(ThetaMinLabel,fLayout1);
  ThetaMinFrame->AddFrame(fChooseThetaMin,fLayout2);
  frame->AddFrame(ThetaMaxFrame,fLayout3);
  ThetaMaxFrame->AddFrame(ThetaMaxLabel,fLayout1);
  ThetaMaxFrame->AddFrame(fChooseThetaMax,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   HarpoSelectorVertex::SetConfig()
{
  // Update the configuration according to the values in the popup window

  if(!fChooseDx) return;

  fDx = fChooseDx->GetNumber();
  fDy = fChooseDy->GetNumber();
  fThetaMin = fChooseThetaMin->GetNumber();
  fThetaMax = fChooseThetaMax->GetNumber();
}