ExtractRF.C

// Macro for Rejection Factor (by EMcal cluster numbers)

// Input : AnalysisResults.root generated by C. Jahnke's task
// Structure
// * PWGDQ_dielectron_EMCal           - TDirectoryFile*
// - cjahnke_QA_[trigger_index]_cent0 - TList
// -- [Cut_Definition]                - THashList (histogram classed)
// --- Track_ev1+/-                   - THashList (histograms)
// ---- EMCal_E                       - TH1 (for rejection factor)
// --- Pair_ev1+/-_ev1+/-
// --- Track_Legs_ev1+/-_ev1+/-
// --- Event - ONLY in RAW
// ---- VtxZ                          - TH1 (for event numbers)

#include "TFile.h"
#include "TString.h"

// Trigger_index : MB=100, EG1=3, EG2=4
// Cut definition: TPC, EMCal, EMCal_loose, EMCal_strict, RAW
Int_t GetEventNumber(TFile* resultFile, Int_t trigger_index){
  TString rootDir = "PWGDQ_dielectron_EMCal";
  TDirectoryFile* rootDirObj = (TDirectoryFile*)(resultFile->Get(rootDir));
  TString triggerDir = Form("cjahnke_QA_%d_cent0", trigger_index);
  TList* triggerDirObj = (TList*)(rootDirObj->FindKey(triggerDir)->ReadObj());
  TH1* vtxz = (TH1*)(triggerDirObj
    ->FindObject("RAW")
    ->FindObject("Event")
    ->FindObject("VtxZ"));
  if(!vtxz){
    cout << "[X] ERROR - Object not found : " << rootDir << "/" << triggerDir << "/RAW/Event/VtxZ" << endl;
    return 0;
  }
  else
    return vtxz->GetEntries();
  return 0;
}
TH1* GetEMCalE(TFile* resultFile, Int_t trigger_index, TString cutDef){

  TString rootDir = "PWGDQ_dielectron_EMCal";
  TDirectoryFile* rootDirObj = (TDirectoryFile*)(resultFile->Get(rootDir));
  TString triggerDir = Form("cjahnke_QA_%d_cent0", trigger_index);
  TList* triggerDirObj = (TList*)(rootDirObj->FindKey(triggerDir)->ReadObj());

  TH1* hEMCalE = (TH1*)(triggerDirObj
    ->FindObject(cutDef)
    ->FindObject("Track_ev1+")->FindObject("EMCal_E"));

  TH1* hEMCalE_minus = (TH1*)(triggerDirObj
    ->FindObject(cutDef)
    ->FindObject("Track_ev1-")->FindObject("EMCal_E"));
  
	if(!hEMCalE && !hEMCalE_minus){
    cout << "[X] ERROR - Object not found : " << rootDir << "/" << triggerDir << "/" + cutDef << "/Track_ev1+/EMCal_E" << endl;
    return NULL;
  }
  else{
		hEMCalE->Add(hEMCalE_minus);
    return hEMCalE;
	}

  return NULL;
}
int ExtractRF(TFile* resultFile, TString cutDef = "RAW"){
  
  if(!resultFile){
    cout << "[X] ERROR - Result file is NULL" << endl;
    return -1;
  }
 
  // N_event for normalization
  Int_t nMB = GetEventNumber(resultFile, 100);
  Int_t nEG1 = GetEventNumber(resultFile, 3);
  Int_t nEG2 = GetEventNumber(resultFile, 4);

  // EMCal_E histogram
  TH1* hMB = GetEMCalE(resultFile, 100, cutDef);
  TH1* hEG1 = GetEMCalE(resultFile, 3, cutDef);
  TH1* hEG2 = GetEMCalE(resultFile, 4, cutDef);

  if(!nMB || !nEG1 || !nEG2 || !hMB || !hEG1 || !hEG2){
    return -1;
  }
  // Global style
  gStyle->SetOptFit(0000);
  gStyle->SetOptStat(0);

  // Function for event normalization
  TF1* fcn = new TF1("fNev","[0]",0, 50);

  // Normalized EMCal cluster energy
  fcn->SetParameter(0, nMB * hMB->GetBinWidth(1));
  hMB->SetName("hMB");
  hMB->SetLineColor(kRed);
  hMB->Divide(fcn);
  
  fcn->SetParameter(0, nEG1 * hEG1->GetBinWidth(1));
  hEG1->SetName("hEG1");
  hEG1->SetLineColor(kBlue);
  hEG1->Divide(fcn);

  fcn->SetParameter(0, nEG2 * hEG2->GetBinWidth(1));
  hEG2->SetName("hEG2");
  hEG2->SetLineColor(kGreen);
  hEG2->Divide(fcn);

  TCanvas* cEMCalE = new  TCanvas("c1", "EMCal cluster energy", 800, 600);
  cEMCalE->SetLogy(1);
  hMB->Draw();
  hMB->GetYaxis()->SetRangeUser(1e-8, 1);
  hMB->GetYaxis()->SetTitle("Number of EMCal clusters / N_{event}");
  hMB->GetYaxis()->SetTimeOffset(1.20);
  hEG1->Draw("same");
  hEG2->Draw("same");

  auto yLgd= new TLegend(0.13,0.60,0.49, 0.88, "", "brNDC");
  yLgd->SetName("eclsLgd");
  yLgd->SetBorderSize(0);
  yLgd->SetTextAlign(12);
  yLgd->SetTextFont(42);
  yLgd->SetTextSize(0.04);
  yLgd->AddEntry(hMB,  Form("MB   (kINT7)  - %.1f M", nMB/1e6));
  yLgd->AddEntry(hEG1, Form("EG1 (7 GeV) - %.1f M", nEG1/1e6));
  yLgd->AddEntry(hEG2, Form("EG2 (5 GeV) - %.1f M", nEG2/1e6));
  yLgd->Draw("same");

  // Calculate Rejection Factor
  TH1* rf1 = (TH1*)(hEG1->Clone());
  rf1->SetName("hRF1");
  rf1->SetMarkerStyle(2);
  rf1->SetMarkerColor(kBlue);
  rf1->Divide(hMB);
  
  TH1* rf2 = (TH1*)(hEG2->Clone());
  rf2->SetName("hRF2");
  rf2->SetMarkerStyle(4);
  rf2->SetMarkerColor(kGreen);
  rf2->Divide(hMB);

  cEMCalE = new  TCanvas("c2", "Rejection Factor", 800, 600);
  cEMCalE->SetLogy(1);
  cEMCalE->SetGridx(1);
  cEMCalE->SetGridy(1);

  rf1->Draw("EP");
  rf1->GetYaxis()->SetMoreLogLabels(kTRUE);
  rf1->GetYaxis()->SetTitle("Rejection Factor");
  rf1->GetYaxis()->SetTimeOffset(1.50);
  rf2->Draw("same EP");

  yLgd= new TLegend(0.13,0.60,0.49, 0.88, "", "brNDC");
  yLgd->SetName("rfLgd");
  yLgd->SetBorderSize(0);
  yLgd->SetTextAlign(12);
  yLgd->SetTextFont(42);
  yLgd->SetTextSize(0.04);
  yLgd->AddEntry(rf1, "EG1 / MB");
  yLgd->AddEntry(rf2, "EG2 / MB");
  yLgd->Draw("same");

  return 0;
}

int ExtractRF(TString resultPath  = "AnalysisResults.root", TString cutDef = "RAW"){

  TFile* resultFile = new TFile(resultPath);

  if(resultFile->IsOpen())
    return ExtractRF(resultFile, cutDef);
  else{
    cout << "[X] ERROR - File not found : " << resultPath << endl;
    return -1;
  }

  return 0;
}