HarpoAnalyseTemplate.cxx

Go to the documentation of this file.

// 
// File HarpoAnalyseTemplate.cxx
//
#include "HarpoAnalyseTemplate.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 "TpcSimMCEvent.h"
#include "MakeNiceHisto.h"
#include "HarpoTools.h"

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

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

ClassImp(HarpoAnalyseTemplate)

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

  // Process RAW data
  for(Int_t ndet=0;ndet<2;ndet++) { // Detectors 0 and 1 = X and Y planes
    if (!gHDetSet->isExist(ndet)) continue;
    HarpoDccMap* m = fEvt->GetDccMap(ndet);
    if ( m == NULL ) continue;
    for(Int_t j=0;j<NCHAN;j++){ // Channels
      //Double_t qch = 0;
      for(Int_t i=0;i<NADC;i++){ // Time bins
        Double_t q = m->GetData(j,i);
        if(q <= -1000) continue;
        hQraw[ndet]->Fill(q);
      }
    }
  }
  
  if(gHDetSet->isExist(PMM2)) { // PMm2 (trigger) data
    Pmm2Event *anaEvt = static_cast<Pmm2Event *>(fEvt->GetDetEvent(PMM2));
    //    Int_t triggerType = anaEvt->GetTriggerType(); // Trigger line (with mesh, without, traversing tracks, ...)
    
    int imes;
    int esize = anaEvt->GetHeader()->eventSize;   // number of PMs hit
    Pmm2MesVect * pm = anaEvt->GetMesurements();  // data for each included PM
    for(imes=0;imes<esize;imes++) {
      int ch    = pm->at(imes).getChNum(); // channel
      int q     = pm->at(imes).getCharge(); // charge
      hPmm2->Fill(ch,q);
    } 
  }

  Double_t data[5] = {-1,-1,-1,-1,-1};


  //???????????????????????????????
  if(fEvt->GetRecoEvent()){  // Reconstructed data (clusters, tracks, matched tracks)
    HarpoRecoEvent* reco = fEvt->GetRecoEvent();
    // Clusters in the event (both X and Y directions)
    TClonesArray* clArray = reco->GetClustersArray();
    Int_t nCl = clArray->GetEntries();
    Int_t nClX = 0, nClY = 0;
    for(Int_t icl = 0; icl<nCl; icl++){
      HarpoRecoClusters* cluster = (HarpoRecoClusters*)clArray->At(icl);
      Double_t q = cluster->GetQ();
      Int_t plane = cluster->GetPlane();
      if(plane==XPLANE) nClX++;
      if(plane==YPLANE) nClY++;
      hQ[plane]->Fill(q);
    }
    // Tracks in the event (both X and Y directions)
    TClonesArray* trArray = reco->GetTracksArray();
    Int_t nTr = trArray->GetEntries();
    Int_t nTrX = 0, nTrY = 0;
    for(Int_t itr = 0; itr<nTr; itr++){
      HarpoRecoTracks* track = (HarpoRecoTracks*)trArray->At(itr);
      //      Double_t id = track->GetNtrack();
      Double_t angle = track->GetAngleTrack();
      hAngle->Fill(angle);
      Int_t plane = track->GetPlane();
      if(plane==XPLANE) nTrX++;
      if(plane==YPLANE) nTrY++;
    }
    // 2D Vertexes in the event (both X and Y directions)
    TClonesArray* vArray = reco->GetVertexArray();
    Int_t nV = vArray->GetEntries();
    Int_t nVx = 0, nVy = 0;
    for(Int_t iV = 0; iV<nV; iV++){
      HarpoRecoVertex* vertex = (HarpoRecoVertex*)vArray->At(iV);
      Int_t plane = vertex->GetPlane();
      if(plane==XPLANE) nVx++;
      if(plane==YPLANE) nVy++;
    }
    data[0] = nTrX;
    data[1] = nTrY;
    data[2] = nClX;
    data[3] = nClY;
  }




   if (gHDetSet->isExist(SIMDET)){
     HarpoSimEvent *simEvt = (HarpoSimEvent *)(fEvt->GetDetEvent(SIMDET));
     Int_t nIonTr = simEvt->GetNtracks();
     for(Int_t i = 0; i<nIonTr; i++){
       TpcSimIonisationTrack* tr = simEvt->GetIonisationTrack(i);
       Int_t nPoints = tr->GetNpoints();
       for(Int_t j = 0; j<nPoints; j++){
         TpcSimIonisationPoint* point = tr->GetPoint(j);
         hNeSim->Fill(point->GetNe());
       }
     }
     TpcSimMCEvent* mcEvent = simEvt->GetMCEvent();
     Int_t nMCtr = mcEvent->GetNtracks();
     data[4] = nMCtr;
     for(Int_t i = 0; i<nMCtr; i++){
       TpcSimMCTrack* track = mcEvent->GetTrack(i);
       hPsim->Fill(track->GetPx(),track->GetPy(),track->GetPz());
     }
   }


   Double_t energy = getDoubleVal(run,"Ephoton",dburl);
   cout << "Energy = " << energy << endl;


   fNtuple->Fill(data);

}

void   HarpoAnalyseTemplate::Init()
{
  
  fChooseNbins = 0;

  // Initialise histograms here
   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 qmin;
   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;


   hQ[XPLANE] = HistLog("hQx","",fNbins,fQmin,fQmax);
   hQ[YPLANE] = HistLog("hQy","",fNbins,fQmin,fQmax);
   hNeSim = new TH1F("hNeSim","",500,0,500);
   hAngle = new TH1F("hAngle","",500,-10,10);
   hQraw[XPLANE] = new TH1F("hQrawX","",500,-500,5500);
   hQraw[YPLANE] = new TH1F("hQrawY","",500,-500,5500);
   hPmm2 = new TH2F("hPmm2","",16,0,16,128,0,4096);
   hPsim = new TH3F("hPsim","",20,-1,1,20,-1,1,20,-1,1);


   fNtuple = new TNtupleD("fNtuple","","nTrX:nTrY:nClX:nClY:nTrsim");

 }

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

   OpenHistFile("template");

   // Save histograms here
   hQ[XPLANE]->Write();
   hQ[YPLANE]->Write();
   hQraw[XPLANE]->Write();
   hQraw[YPLANE]->Write();
   hNeSim->Write();
   hAngle->Write();
   hPmm2->Write();
   hPsim->Write();

 }



void HarpoAnalyseTemplate::DisplayAnalysis(TRootEmbeddedCanvas* ecTab, TGListBox* /* infobox */ )
{
  // in hrecomonitorgui: fills the canvas on the "Analysis" tab
  TCanvas* c = ecTab->GetCanvas();
  c->GetPad(1)->Delete();
  c->GetPad(2)->Delete();
  c->Divide(2);

  for(Int_t plane = 0; plane<2;plane++){
    TVirtualPad* cMap = c->GetPad(plane+1);
    HarpoDccMap* m = fEvt->GetDccMap(plane);
    TH2F* h = new TH2F("h","",512,0,512,288,0,288);
    for(Int_t i=1;i<NADC;i++){ //Time bins
      for(Int_t j=0;j<NCHAN;j++){ // Channels
        Double_t q = m->GetData(j,i);
        h->SetBinContent(i+1,j+1,q);
      }
    } 
    h->SetMinimum(0);
    MakeNiceHisto(h,cMap,"colz");
    h->Delete();
  }

  c->Update();
}


void HarpoAnalyseTemplate::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* nBinsFrame = new TGHorizontalFrame(frame,xsize,ysize2,kHorizontalFrame);
  TGLabel* nBinsLabel = new TGLabel(nBinsFrame, "fNbins");
  fChooseNbins    = new TGNumberEntry(nBinsFrame);
  fChooseNbins->SetNumber(fNbins);

  main->AddFrame(frame,fLayout4);
  frame->AddFrame(fAnalysisLabel,fLayout3);
  frame->AddFrame(nBinsFrame,fLayout3);
  nBinsFrame->AddFrame(nBinsLabel,fLayout1);
  nBinsFrame->AddFrame(fChooseNbins,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   HarpoAnalyseTemplate::SetConfig()
{
  // Update the configuration according to the values in the popup window

  if(!fChooseNbins) return;

  fNbins = fChooseNbins->GetNumber();
}