Add JE validation

codeJE/Compare.C

@@ -0,0 +1,214 @@
+// Comparison of AliPhysics and O2 histograms
+
+#include "utils_plot.h"
+
+// vectors of histogram specifications
+using VecSpecHis = std::vector<std::tuple<TString, TString, TString, int, bool, bool, TString>>;
+
+// Add histogram specification in the vector.
+void AddHistogram(VecSpecHis& vec, TString label, TString nameRun2, TString nameRun3, int rebin, bool logH, bool logR, TString proj = "x")
+{
+  vec.push_back(std::make_tuple(label, nameRun2, nameRun3, rebin, logH, logR, proj));
+}
+
+Int_t Compare(TString filerun3 = "AnalysisResults_O2.root", TString filerun2 = "AnalysisResults_ALI.root", TString options = "jets", bool doRatio = true)
+{
+  gStyle->SetOptStat(0);
+  gStyle->SetPalette(0);
+  gStyle->SetCanvasColor(0);
+  gStyle->SetFrameFillColor(0);
+
+  TFile* fRun3 = new TFile(filerun3.Data());
+  if (fRun3->IsZombie()) {
+    printf("Failed to open file %s\n", filerun3.Data());
+    return 1;
+  }
+  TFile* fRun2 = new TFile(filerun2.Data());
+  if (fRun2->IsZombie()) {
+    printf("Failed to open file %s\n", filerun2.Data());
+    return 1;
+  }
+
+  TString pathListRun2 = Form("ChJetSpectraAliAnalysisTaskEmcalJetValidation/AliAnalysisTaskEmcalJetValidation"); //
+  TList* lRun2 = nullptr;
+  fRun2->GetObject(pathListRun2.Data(), lRun2);
+  if (!lRun2) {
+    printf("Failed to load list %s from %s\n", pathListRun2.Data(), filerun2.Data());
+    return 1;
+  }
+
+  // Histogram specification: axis label, AliPhysics name, O2Physics path/name, rebin, log scale histogram, log scale ratio, projection axis
+
+  VecSpecHis vecHisEvents;
+
+  VecSpecHis vecHisJets;                                                                                                  //X axis ,                  Hist Name Run2,  Hist  Name in Run3, rebin, logScaleH , logScaleR
+  AddHistogram(vecHisJets, "#it{p}_{T,ch jet} (GeV/#it{c}", "jetPt", "jet-validation-track-collision-qa/jetPt", 1, 1, 0); //
+  AddHistogram(vecHisJets, "#varphi_{jet}", "jetPhi", "jet-validation-track-collision-qa/jetPhi", 1, 0, 0);
+  AddHistogram(vecHisJets, "#eta_{jet}", "jetEta", "jet-validation-track-collision-qa/jetEta", 1, 0, 0);
+  AddHistogram(vecHisJets, "#it{p}_{T,track} (GeV/#it{c}", "jetTrackPt", "jet-validation-track-collision-qa/selectedTrackPt", 1, 1, 0); //
+  AddHistogram(vecHisJets, "#varphi_{track}", "jetTrackPhi", "jet-validation-track-collision-qa/selectedTrackPhi", 1, 0, 0);
+  AddHistogram(vecHisJets, "#eta_{track}", "jetTrackEta", "jet-validation-track-collision-qa/selectedTrackEta", 1, 0, 0);
+
+  // vector of specifications of vectors: name, VecSpecHis, pads X, pads Y
+  std::vector<std::tuple<TString, VecSpecHis, int, int>> vecSpecVecSpec;
+
+  // Add vector specifications in the vector.
+  // if (options.Contains(" events "))
+  //  vecSpecVecSpec.push_back(std::make_tuple("events", vecHisEvents, 4, 2));
+  if (options.Contains("jets"))
+    vecSpecVecSpec.push_back(std::make_tuple("jets", vecHisJets, 3, 2));
+  // if (options.Contains(" skim "))
+  //  vecSpecVecSpec.push_back(std::make_tuple("skim", vecHisSkim, 5, 3));
+
+  // Histogram plot vertical margins
+  Float_t marginHigh = 0.05;
+  Float_t marginLow = 0.05;
+  bool logScaleH = false;
+  // Ratio plot vertical margins
+  Float_t marginRHigh = 0.05;
+  Float_t marginRLow = 0.05;
+  bool logScaleR = false;
+  Float_t yMin, yMax;
+  Int_t nRun2, nRun3, rebin;
+
+  TH1F* hRun2 = nullptr;
+  TH1D* hRun3 = nullptr;
+  TH1F* hRatio = nullptr;
+  TString labelAxis = "";
+  TString nameHisRun2 = "";
+  TString nameHisRun3 = "";
+  TString projAx = "";
+  TCanvas* canHis = nullptr;
+  TCanvas* canRat = nullptr;
+
+  // loop over lists
+  for (const auto& specVecSpec : vecSpecVecSpec) {
+    auto nameSpec = std::get<0>(specVecSpec); // list name
+    auto vecSpec = std::get<1>(specVecSpec);  // list of histogram specs.
+    int nPadsX = std::get<2>(specVecSpec);    // number of horizontal pads
+    int nPadsY = std::get<3>(specVecSpec);    // number of vertical pads
+    Printf("\nProcessing histogram list: %s (%d)", nameSpec.Data(), (int)vecSpec.size());
+    if (nPadsX * nPadsY < vecSpec.size()) {
+      Printf("Not enough pads (%d)", nPadsX * nPadsY);
+      return 1;
+    }
+
+    canHis = new TCanvas(Form("canHis_%s", nameSpec.Data()), Form("Histos_%s", nameSpec.Data()), 3000, 1600);
+    SetCanvas(canHis, nPadsX, nPadsY);
+    if (doRatio) {
+      canRat = new TCanvas(Form("canRat_%s", nameSpec.Data()), Form("Ratios_%s", nameSpec.Data()), 3000, 1600);
+      SetCanvas(canRat, nPadsX, nPadsY);
+    }
+
+    // loop over histograms
+    for (int index = 0; index < vecSpec.size(); index++) {
+      auto spec = vecSpec[index];
+      labelAxis = std::get<0>(spec);
+      nameHisRun2 = std::get<1>(spec);
+      nameHisRun3 = std::get<2>(spec);
+      rebin = std::get<3>(spec);
+      logScaleH = std::get<4>(spec);
+      logScaleR = std::get<5>(spec);
+      projAx = std::get<6>(spec);
+
+      // Get AliPhysics histogram.
+      hRun2 = (TH1F*)lRun2->FindObject(nameHisRun2.Data());
+      if (!hRun2) {
+        printf("Failed to load %s from %s\n", nameHisRun2.Data(), filerun2.Data());
+        return 1;
+      }
+
+      // Get O2 histogram.
+      auto oRun3 = fRun3->Get(nameHisRun3.Data());
+      if (!oRun3) {
+        printf("Failed to load %s from %s\n", nameHisRun3.Data(), filerun3.Data());
+        return 1;
+      }
+
+      if (oRun3->InheritsFrom("TH3")) {
+        if (projAx == "x") {
+          hRun3 = ((TH3D*)oRun3)->ProjectionX();
+        } else if (projAx == "y") {
+          hRun3 = ((TH3D*)oRun3)->ProjectionY();
+        }
+      } else if (oRun3->InheritsFrom("TH2")) {
+        if (projAx == "x") {
+          hRun3 = ((TH2D*)oRun3)->ProjectionX();
+        } else if (projAx == "y") {
+          hRun3 = ((TH2D*)oRun3)->ProjectionY();
+        }
+      } else {
+        hRun3 = (TH1D*)oRun3;
+      }
+
+      Printf("%d (%s, %s): bins: %d, %d, ranges: %g-%g, %g-%g",
+             index, nameHisRun2.Data(), nameHisRun3.Data(),
+             hRun2->GetNbinsX(), hRun3->GetNbinsX(),
+             hRun2->GetXaxis()->GetBinLowEdge(1), hRun2->GetXaxis()->GetBinUpEdge(hRun2->GetNbinsX()),
+             hRun3->GetXaxis()->GetBinLowEdge(1), hRun3->GetXaxis()->GetBinUpEdge(hRun3->GetNbinsX()));
+
+      nRun2 = hRun2->GetEntries();
+      nRun3 = hRun3->GetEntries();
+      // Histograms
+      auto padH = canHis->cd(index + 1);
+
+      Printf("histo name: %s", hRun2->GetName());
+      hRun2->Rebin(rebin);
+      hRun3->Rebin(rebin);
+      hRun2->SetLineColor(1);
+      hRun2->SetLineWidth(2);
+      hRun2->SetMarkerStyle(22);
+      hRun2->SetMarkerColor(1);
+      hRun3->SetLineColor(2);
+      hRun3->SetLineWidth(2);
+      hRun3->SetMarkerStyle(22);
+      hRun3->SetMarkerColor(2);
+      hRun2->SetTitle(" ");
+      hRun2->GetYaxis()->SetTitle("number of entries");
+      hRun2->GetYaxis()->SetMaxDigits(3);
+      yMin = TMath::Min(hRun3->GetMinimum(0), hRun2->GetMinimum(0));
+      yMax = TMath::Max(hRun3->GetMaximum(), hRun2->GetMaximum());
+      SetHistogram(hRun2, yMin, yMax, marginLow, marginHigh, logScaleH);
+      SetPad(padH, logScaleH);
+      hRun2->Draw("");
+      hRun3->Draw("Esame");
+
+      if (std::string(hRun2->GetName()) == "jetPt" || std::string(hRun2->GetName()) == "jetTrackPt") {
+        TLegend* legend = new TLegend(0.2, 0.84, 0.82, 0.92);
+        legend->SetNColumns(2);
+        legend->AddEntry(hRun2, Form("AliPhysics: %d", nRun2), "L");
+        legend->AddEntry(hRun3, Form("O2Physics: %d", nRun3), "L");
+        legend->Draw();
+      } else {
+        TLegend* legend = new TLegend(0.2, 0.92, 0.82, 1.0);
+        legend->SetNColumns(2);
+        legend->SetBorderSize(0);
+        legend->AddEntry(hRun2, Form("AliPhysics: %d", nRun2), "L");
+        legend->AddEntry(hRun3, Form("O2Physics: %d", nRun3), "L");
+        legend->Draw();
+      }
+
+      // Ratio
+      if (doRatio) {
+        auto padR = canRat->cd(index + 1);
+        hRatio = (TH1F*)hRun3->Clone(Form("hRatio%d", index));
+        hRatio->Divide(hRun2);
+        hRatio->SetTitle("");
+        hRatio->GetYaxis()->SetTitle("O2Physics/AliPhysics");
+        yMin = hRatio->GetMinimum(0) * 0.3;
+        yMax = hRatio->GetMaximum() * 1.3;
+        SetHistogram(hRatio, yMin, yMax, marginRLow, marginRHigh, logScaleR);
+        SetPad(padR, logScaleR);
+        hRatio->Draw();
+      }
+    }
+    canHis->SaveAs(Form("comparison_histos_%s.png", nameSpec.Data()));
+    if (doRatio) {
+      canRat->SaveAs(Form("comparison_ratios_%s.png", nameSpec.Data()));
+    }
+    //    delete canHis;
+    delete canRat;
+  }
+
+  return 0;
+}

codeJE/RunJetTaskLocal.C

@@ -0,0 +1,101 @@
+#ifdef __CLING__
+// Tell ROOT where to find AliRoot headers
+R__ADD_INCLUDE_PATH($ALICE_ROOT)
+// Tell ROOT where to find AliPhysics headers
+R__ADD_INCLUDE_PATH($ALICE_PHYSICS)
+
+#include "PWGPP/PilotTrain/AddTaskCDBconnect.C"
+
+#endif
+
+TChain* CreateLocalChain(const char* txtfile);
+
+Long64_t RunJetTaskLocal(TString txtfile = "./list_ali.txt",
+                         TString jsonfilename = "dpl-config_std.json",
+                         Bool_t isMC = kFALSE,
+                         Bool_t useO2Vertexer = kFALSE,
+                         Bool_t useAliEventCuts = kFALSE)
+{
+  // Load common libraries
+  gSystem->Load("libCore.so");
+  gSystem->Load("libTree.so");
+  gSystem->Load("libGeom.so");
+  gSystem->Load("libVMC.so");
+  gSystem->Load("libPhysics.so");
+  gSystem->Load("libSTEERBase");
+  gSystem->Load("libESD");
+  gSystem->Load("libAOD");
+  gSystem->Load("libANALYSIS");
+  gSystem->Load("libANALYSISalice");
+  gSystem->SetIncludePath("-I. -I$ROOTSYS/include -I$ALICE_ROOT -I$ALICE_ROOT/include -I$ALICE_ROOT/ITS -I$ALICE_ROOT/TPC -I$ALICE_ROOT/CONTAINERS -I$ALICE_ROOT/STEER/STEER -I$ALICE_ROOT/STEER/STEERBase -I$ALICE_ROOT/STEER/ESD -I$ALICE_ROOT/STEER/CDB -I$ALICE_ROOT/TRD -I$ALICE_ROOT/macros -I$ALICE_ROOT/ANALYSIS -I$ALICE_PHYSICS -I$ALICE_PHYSICS/include -g");
+
+  AliAnalysisManager* mgr = new AliAnalysisManager("testAnalysis");
+
+  TChain* chainESD = CreateLocalChain(txtfile.Data());
+  if (!chainESD) {
+    Error("CreateLocalChain", "Failed to create chain from file %s", txtfile.Data());
+    return -1;
+  }
+
+  // Create and configure the alien handler plugin
+  AliESDInputHandler* esdH = new AliESDInputHandler();
+  //  esdH->SetNeedField(kTRUE);
+  mgr->SetInputEventHandler(esdH);
+
+  AliMCEventHandler* handler = NULL;
+  if (isMC) {
+    handler = new AliMCEventHandler;
+    handler->SetReadTR(kFALSE);
+    mgr->SetMCtruthEventHandler(handler);
+  }
+
+  // CDBconnect task
+  AliTaskCDBconnect* taskCDB = AddTaskCDBconnect();
+  taskCDB->SetFallBackToRaw(kTRUE);
+
+  // Apply the event selection
+  AliPhysicsSelectionTask* physSelTask = reinterpret_cast<AliPhysicsSelectionTask*>(gInterpreter->ProcessLine(Form(".x %s(%d)", gSystem->ExpandPathName("$ALICE_PHYSICS/OADB/macros/AddTaskPhysicsSelection.C"), isMC)));
+
+  AliAnalysisTaskEmcalJetValidation* taskJet = reinterpret_cast<AliAnalysisTaskEmcalJetValidation*>(gInterpreter->ProcessLine(Form(".x %s(\"\",\"%s\",%d)", gSystem->ExpandPathName("$ALICE_PHYSICS/PWGJE/EMCALJetTasks/macros/AddTaskEmcalJetValidation.C"), jsonfilename.Data(), isMC)));
+  if (useAliEventCuts) {
+    taskJet->SetUseAliEventCuts(useAliEventCuts);
+  }
+  // if (useO2Vertexer) {
+  //   taskJet->SetUseO2Vertexer();
+  // }
+
+  mgr->InitAnalysis();
+  mgr->PrintStatus();
+  return mgr->StartAnalysis("local", chainESD);
+};
+
+TChain* CreateLocalChain(const char* txtfile)
+{
+  // Open the file
+  ifstream in;
+  in.open(txtfile);
+  Int_t count = 0;
+  // Read the input list of files and add them to the chain
+  TString line;
+  TChain* chain = new TChain("esdTree");
+  while (in.good()) {
+    in >> line;
+    if (line.IsNull() || line.BeginsWith("#"))
+      continue;
+    TString esdFile(line);
+    TFile* file = TFile::Open(esdFile);
+    if (file && !file->IsZombie()) {
+      chain->Add(esdFile);
+      file->Close();
+    } else {
+      Error("CreateLocalChain", "Skipping un-openable file: %s", esdFile.Data());
+    }
+  }
+  in.close();
+  if (!chain->GetListOfFiles()->GetEntries()) {
+    Error("CreateLocalChain", "No file from %s could be opened", txtfile);
+    delete chain;
+    return nullptr;
+  }
+  return chain;
+}

codeJE/clean.sh

@@ -0,0 +1,12 @@
+#!/bin/bash
+
+# Script to delete created files
+
+rm -rf \
+AnalysisResults_ALI.root AnalysisResults_O2.root \
+comparison_histos_jets.* comparison_ratios_jets.* \
+./*.log \
+output_* \
+|| { echo "Error: Failed to delete files."; exit 1; }
+
+exit 0