HarpoSelectorStraightTracks.cxx

Go to the documentation of this file.

// 
// File HarpoSelectorStraightTracks.cxx
//
#include "HarpoSelectorStraightTracks.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 <cstdlib>
#include <cstring>
#include <cassert>
#include <fstream>
#include <iostream>

using namespace std;


ClassImp(HarpoSelectorStraightTracks);

void   HarpoSelectorStraightTracks::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   HarpoSelectorStraightTracks::process()
{
  Int_t stat = eventselection();
  // if(stat==0) Info("process","Status = %d",stat);
  fEvt->GetHeader()->SetEvAnaStatus(stat);
  if(stat==0) 
    nEvents++;

  return;
}

Int_t   HarpoSelectorStraightTracks::eventselection()
{

  // if(nEvents%100 == 0)
  //   std::cout << "Event " << nEvents << std::endl;
  
  // Process Reconstructed data (clusters, tracks, matched tracks)
  HarpoRecoEvent* reco = fEvt->GetRecoEvent();
  TClonesArray* trArray = reco->GetHoughTracksArray();
  // TClonesArray* matchArray = reco->GetMatchingArray();
  // Int_t nMatch = matchArray->GetEntries();
  //  if(nMatch<=0) return -2;
  //  if(nMatch!=1) return -2;
  Int_t nTrX = 0, nTrY = 0;
  
  const Int_t nTr = trArray->GetEntries();
  Int_t test = 0;
  for(Int_t itr = 0; itr<nTr; itr++){
    HarpoRecoHoughTracks* tr = (HarpoRecoHoughTracks*)trArray->At(itr);
    Int_t plane = tr->GetPlane();
    if(plane == XPLANE) nTrX++;
    else if(plane == YPLANE) nTrY++;
    else continue;
    Double_t sigRho = tr->GetSigrho();
    hSigmaRho[plane]->Fill(sigRho);
    if(sigRho>fSigmaRhoMax){
      if(gHarpoDebug>1)
        Info("eventselection","Sigrho = %g", sigRho);
      test=-2;
    }else
      hSigmaRhoCut[plane]->Fill(sigRho);
    Double_t sigTheta = tr->GetSigtheta();
    hSigmaTheta[plane]->Fill(sigTheta);
    if(sigTheta>fSigmaThetaMax){
      if(gHarpoDebug>1)
        Info("eventselection","Sigtheta = %g", sigTheta);
      test=-2;
    }else
      hSigmaThetaCut[plane]->Fill(sigTheta);
    Double_t angle = tr->GetAngleTrack();
    hAngle[plane]->Fill(angle);
    if(TMath::Abs(angle)<fThetaMin || TMath::Abs(angle)>fThetaMax){
      if(gHarpoDebug>1)
        Info("eventselection","Angle = %g", angle);
      test=-2;
    }else
      hAngleCut[plane]->Fill(angle);
      
      if(gHarpoDebug>1)
        Info("eventselection","Sigrho = %g, Sigtheta = %g, Angle = %g", sigRho, sigTheta, angle);
  }
  

  //   if(!(nTrX == nMatch && nTrY == nMatch))
  if(!(nTrX == 1 && nTrY == 1)){
    if(gHarpoDebug>1)
    Info("eventselection","nTrX = %d nTrY = %d",nTrX,nTrY);
    test=-2;
  }
   
  return test;
}

void   HarpoSelectorStraightTracks::Init()
{
  
  fSigmaRhoMax = 5.;
  fSigmaThetaMax = 0.4;
  fThetaMin = 0.03;
  fThetaMax = 0.3;

   Double_t sigmaRhoMax;
   if( ! gHConfig->Lookup("SelectorStraightTracks.sigmaRhoMax",sigmaRhoMax))
     Info("Constructor","Use default sigmaRhoMax %.3g",fSigmaRhoMax);
   else
     fSigmaRhoMax = sigmaRhoMax;
   Double_t sigmaThetaMax;
   if( ! gHConfig->Lookup("SelectorStraightTracks.sigmaThetaMax",sigmaThetaMax))
     Info("Constructor","Use default sigmaThetaMax %.3g",fSigmaThetaMax);
   else
     fSigmaThetaMax = sigmaThetaMax;
   Double_t thetaMax;
   if( ! gHConfig->Lookup("SelectorStraightTracks.thetaMax",thetaMax))
     Info("Constructor","Use default thetaMax %.3g",fThetaMax);
   else
     fThetaMax = thetaMax;
   Double_t thetaMin;
   if( ! gHConfig->Lookup("SelectorStraightTracks.thetaMin",thetaMin))
     Info("Constructor","Use default thetaMin %.3g",fThetaMin);
   else
     fThetaMin = thetaMin;

  hAngle[0] = new TH1F("hAngleX",";#theta",200,-TMath::Pi(),TMath::Pi());
  hAngleCut[0] = new TH1F("hAngleCutX",";#theta",200,-TMath::Pi(),TMath::Pi());
  hSigmaRho[0] = new TH1F("hSigmaRhoX",";#sigma_{#rho}",200,0,50);
  hSigmaRhoCut[0] = new TH1F("hSigmaRhoCutX",";#sigma_{#rho}",200,0,50);
  hSigmaTheta[0] = new TH1F("hSigmaThetaX",";#sigma_{#theta}",200,0,TMath::Pi());
  hSigmaThetaCut[0] = new TH1F("hSigmaThetaCutX",";#sigma_{#theta}",200,0,TMath::Pi());
  hAngle[1] = new TH1F("hAngleY",";#theta",200,-TMath::Pi(),TMath::Pi());
  hAngleCut[1] = new TH1F("hAngleCutY",";#theta",200,-TMath::Pi(),TMath::Pi());
  hSigmaRho[1] = new TH1F("hSigmaRhoY",";#sigma_{#rho}",200,0,50);
  hSigmaRhoCut[1] = new TH1F("hSigmaRhoCutY",";#sigma_{#rho}",200,0,50);
  hSigmaTheta[1] = new TH1F("hSigmaThetaY",";#sigma_{#theta}",200,0,TMath::Pi());
  hSigmaThetaCut[1] = new TH1F("hSigmaThetaCutY",";#sigma_{#theta}",200,0,TMath::Pi());

 }

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

   Info("Save","N selected events = %ld",nEvents);

   OpenHistFile("selectorStraightTracks");

   // Save histograms here
   for(Int_t plane = 0; plane<2; plane++){
     hAngle[plane]->Write();
     hAngleCut[plane]->Write();
     hSigmaRho[plane]->Write();
     hSigmaRhoCut[plane]->Write();
     hSigmaTheta[plane]->Write();
     hSigmaThetaCut[plane]->Write();
   }

 
 }