KVDataQualityAuditReporting_INDRA.cpp

#include "KVDataQualityAuditReporting_INDRA.h"
#include "KVDataSetManager.h"
#include "KVINDRA.h"
#include <KVGeoDNTrajectory.h>
#include <TColor.h>
#include <TLatex.h>
#include <TLegend.h>
#include <TMarker.h>
#include "TLegend.h"
 
ClassImp(KVDataQualityAuditReporting_INDRA)
 
 
 
void KVDataQualityAuditReporting_INDRA::make_canvas(canvas_t style)
{
   // Make an A4-size canvas
   Double_t w = 297 * 6;
   Double_t h = 210 * 6;
   if (style == canvas_t::kPortrait) std::swap(w, h);
   myCanvas = ::new TCanvas("c", "c", w, h);
   myCanvas->SetWindowSize(w + (w - myCanvas->GetWw()), h + (h - myCanvas->GetWh()));
}
 
 
 
 
void KVDataQualityAuditReporting_INDRA::do_report()
{
   if (!gMultiDetArray) {
      auto dsm = new KVDataSetManager;
      dsm->Init();
      KVMultiDetArray::MakeMultiDetector(dataset_name);
   }
 
   make_canvas();
 
   TString pdf_file = Form("%s_DataQualityAudit_%s.pdf", dataset_name.Data(), fAudit->GetName());
   auto first_page = pdf_file + "(";
   auto last_page = pdf_file + ")";
   current_page = first_page;
   // INDRA Z identification bilan
   int nx(2), ny(3);
   myCanvas->Divide(nx, ny);
   int pad = 1;
   std::vector<TString> indra_id_types = {"PHOS_R_L", "CI_SI", "SI_CSI", "CI_CSI", "CI_SI75", "SI75_SILI", indra_csi_idtype};
   for (auto id : indra_id_types) {
      for (int ring = 1; ring <= 17; ++ring) {
         if (id != "PHOS_R_L" && ring == 1) continue;
         if (id == "PHOS_R_L" && ring > 1) break;
         if ((id == "CI_SI" || id == "SI_CSI") && ring > 9) break;
         if ((id == "CI_CSI" || id.Contains("SI75")) && ring < 10) continue;
 
         myCanvas->cd(pad)->SetFillStyle(4000);//transparent pad
         INDRA_ring_reporting_Z(ring, id);
         ++pad;
         if (pad > nx * ny) {
            draw_sidebar_legend();
            myCanvas->Print(current_page, Form("Title:INDRA Z %s ID Quality by Ring", id.Data()));
            if (current_page == first_page) current_page = pdf_file;
            myCanvas->Clear();
            myCanvas->Divide(nx, ny);
            pad = 1;
         }
      }
      if (pad != 1) {
         draw_sidebar_legend();
         myCanvas->Print(current_page, Form("Title:INDRA Z %s ID Quality by Ring", id.Data()));
         if (current_page == first_page) current_page = pdf_file;
         myCanvas->Clear();
         myCanvas->Divide(nx, ny);
         pad = 1;
      }
      nx = 4;
      ny = 3;
      myCanvas->Clear();
      myCanvas->Divide(nx, ny);
      pad = 1;
      for (int ring = 1; ring <= 17; ++ring) {
         INDRA_ring_Z_threshold_vs_Z(ring, id, pad, nx, ny);
      }
      if (pad != 1) {
         myCanvas->Print(current_page, Form("Title:INDRA Z thresh. vs Z %s by Ring", id.Data()));
         myCanvas->Clear();
         myCanvas->Divide(nx, ny);
         pad = 1;
      }
      if (id == indra_csi_idtype) {
         for (int ring = 2; ring <= 17; ++ring) {
            INDRA_ring_mean_A_vs_Z(ring, id, pad, nx, ny);
         }
      }
      else if (id == "SI75_SILI") {
         for (int ring = 10; ring <= 17; ++ring) {
            INDRA_ring_mean_A_vs_Z(ring, id, pad, nx, ny);
         }
      }
      if (pad != 1) {
         myCanvas->Print(current_page, Form("Title:INDRA <A> vs Z %s by Ring", id.Data()));
         myCanvas->Clear();
         myCanvas->Divide(nx, ny);
         pad = 1;
      }
   }
 
   myCanvas->Clear();
   myCanvas->Print(last_page, "Title:Last page");
}
 
 
 
 
void KVDataQualityAuditReporting_INDRA::INDRA_ring_reporting_Z(int ring, const TString& idtype)
{
   if (!gIndra) {
      Error("INDRA_ring_reporting", "You need to build INDRA first...");
      return;
   }
   TList tels;
   fill_list_of_telescopes(idtype, tels, ring);
 
   auto get_module_number = [](const KVIDTelescope * idt) {
      return (int)dynamic_cast<KVINDRADetector*>(idt->GetDetectors()->Last())->GetModuleNumber();
   };
   auto zmean = fReport.get_mean_Z_for_telescopes(&tels, get_module_number);
   auto zmax = fReport.get_max_Z_for_telescopes(&tels, get_module_number, 24);
   auto zmax_iso = fReport.get_max_Z_with_isotopes_for_telescopes(&tels, get_module_number, 30);
   zmax_iso->SetMarkerColor(kOrange - 2);
   auto zmin = fReport.get_min_Z_for_telescopes(&tels, get_module_number, 25);
   TMultiGraph* mg = ::new TMultiGraph;
   mg->SetTitle(Form("Ring %d %s Min/Mean/Max Z vs. Module", ring, idtype.Data()));
   mg->Add(zmin);
   mg->Add(zmean);
   mg->Add(zmax);
   mg->Add(zmax_iso);
   mg->Draw("ap");
}
 
 
 
 
void KVDataQualityAuditReporting_INDRA::fill_list_of_telescopes(const TString& idtype, TList& tels, int ring)
{
   for (int mod = 1; mod <= 24; ++mod) {
      TString name;
      if (ring > 1)
         name = Form("%s_%02d%02d", idtype.Data(), ring, mod);
      else
         name = Form("%s_%02d", idtype.Data(), mod);
      auto tel = gIndra->GetIDTelescope(name);
      if (tel) tels.Add(tel);
   }
}
 
 
 
 
void KVDataQualityAuditReporting_INDRA::INDRA_ring_mean_A_vs_Z(int ring, const TString& idtype, int& pad, int nx, int ny)
{
   TList tels;
   fill_list_of_telescopes(idtype, tels, ring);
   int nmods = tels.GetEntries();
   auto mod_set = get_mod_set(nmods);
   int color_step = TColor::GetPalette().GetSize() / (mod_set + 1);
   KVIDTelescope* idt;
   nmods = 0;
   while (tels.GetEntries()) {
      TMultiGraph* mg = ::new TMultiGraph;
      int i = 1;
      while (i <= mod_set) {
         idt = (KVIDTelescope*)tels.Remove(tels.First()); // "pop" the first one in the list
         if (fAudit->HasTelescope(idt->GetName())) {
            auto gr = fReport[idt->GetName()].get_mean_isotopic_mass_by_Z();
            gr->SetMarkerColor(TColor::GetPalette()[color_step * i]);
            gr->SetMarkerStyle(markers[i - 1]);
            gr->SetLineWidth(0);
            mg->Add(gr);
         }
         ++i;
         ++nmods;
      }
      myCanvas->cd(pad);
      ++pad;
      // change title of graph => title of pad
      mg->SetTitle(Form("INDRA <A> vs. Z %s Ring %d [%d-%d]", idtype.Data(), ring, nmods - mod_set + 1, nmods));
      mg->Draw("ap");
      TLegend* leg;
      if (i > 12) {
         leg = gPad->BuildLegend(.11, .89, .61, .69);
         leg->SetNColumns(3);
      }
      else if (i > 6) {
         leg = gPad->BuildLegend(.11, .89, .61, .69);
         leg->SetNColumns(2);
      }
      else {
         leg = gPad->BuildLegend(.11, .89, .61, .69);
      }
      leg->SetBorderSize(0);
      leg->SetFillColorAlpha(kWhite, 1.00);
      if (pad > nx * ny) {
         myCanvas->Print(current_page, Form("Title:INDRA <A> vs Z %s by Ring", idtype.Data()));
         myCanvas->Clear();
         myCanvas->Divide(nx, ny);
         pad = 1;
      }
   }
}
 
 
 
 
int KVDataQualityAuditReporting_INDRA::get_mod_set(int nmods)
{
   // 24 => 4 sets of 6, 16 => 2 sets of 8, 12 => 4 sets of 3, 8 => 1 set of 8
   int mod_set;
   switch (nmods) {
      case 24:
         mod_set = 6;
         break;
      case 16:
         mod_set = 8;
         break;
      case 12:
         mod_set = 3;
         break;
      case 8:
         mod_set = 8;
         break;
      default:
         mod_set = 1;
   }
 
   return mod_set;
}
 
 
 
 
void KVDataQualityAuditReporting_INDRA::INDRA_ring_Z_threshold_vs_Z(int ring, const TString& idtype, int& pad, int nx, int ny)
{
   TList tels;
   fill_list_of_telescopes(idtype, tels, ring);
   int nmods = tels.GetEntries();
   auto mod_set = get_mod_set(nmods);
 
   int color_step = TColor::GetPalette().GetSize() / (mod_set + 1);
   KVIDTelescope* idt;
   nmods = 0;
   while (tels.GetEntries()) {
      TMultiGraph* mg = ::new TMultiGraph;
      int i = 1;
      while (i <= mod_set) {
         idt = (KVIDTelescope*)tels.Remove(tels.First()); // "pop" the first one in the list
         if (fAudit->HasTelescope(idt->GetName())) {
            auto gr = fReport[idt->GetName()].get_element_thresholds_by_Z_mev_per_nuc(TColor::GetPalette()[color_step * i]);
            gr->SetMarkerStyle(markers[i - 1]);
            gr->SetLineWidth(0);
            mg->Add(gr);
            if (gr->GetMean(2) < 0) {
               // threshold = -1 => uncalibrated particles
               std::cout << idt->GetName() << " : particles are UNCALIBRATED";
               // check detector calibrations
               KVGeoDNTrajectory* traj;
               if (idt->GetSize() > 1) traj = (KVGeoDNTrajectory*)idt->GetDetector(2)->GetNode()->GetForwardTrajectories()->First();
               else traj = (KVGeoDNTrajectory*)idt->GetDetector(1)->GetNode()->GetForwardTrajectories()->First();
               int nuncal = 0;
               traj->IterateFrom();
               KVGeoDetectorNode* dn;
               while ((dn = traj->GetNextNode())) nuncal += (!dn->GetDetector()->IsCalibrated());
               if (nuncal == traj->GetN()) std::cout << " ... just like ALL DETECTORS";
               else {
                  std::cout << " ... just like ";
                  traj->IterateFrom();
                  KVGeoDetectorNode* dn;
                  while ((dn = traj->GetNextNode())) {
                     if (!dn->GetDetector()->IsCalibrated()) {
                        std::cout << dn->GetName() << " ";
                     }
                  }
               }
               std::cout << std::endl;
            }
         }
         ++i;
         ++nmods;
      }
      myCanvas->cd(pad);
      ++pad;
      // change title of graph => title of pad
      mg->SetTitle(Form("INDRA Z thresh. [MeV/u] vs. Z %s Ring %d [%d-%d]", idtype.Data(), ring, nmods - mod_set + 1, nmods));
      mg->Draw("ap");
      double x1, x2;
//        if(idtype == "CSI")
//        {
//            return;// no legend - doesn't fit
//            //x1 = .39; x2 = .89;
//        }
//        else
//        {
      x1 = .11;
      x2 = .61;
//        }
      TLegend* leg;
      if (mod_set > 12) {
         leg = gPad->BuildLegend(x1, .89, x2, .69);
         leg->SetNColumns(3);
      }
      else if (mod_set > 6) {
         leg = gPad->BuildLegend(x1, .89, x2, .69);
         leg->SetNColumns(2);
      }
      else {
         leg = gPad->BuildLegend(x1, .89, x2, .69);
      }
      leg->SetBorderSize(0);
      leg->SetFillColorAlpha(kWhite, 1.00);
      if (pad > nx * ny) {
         myCanvas->Print(current_page, Form("Title:INDRA Z thresh. vs Z %s by Ring", idtype.Data()));
         myCanvas->Clear();
         myCanvas->Divide(nx, ny);
         pad = 1;
      }
   }
}
 
 
 
 
void KVDataQualityAuditReporting_INDRA::draw_sidebar_legend()
{
//=========Macro generated from canvas: c1/c1
//=========  (Tue Jul 27 11:21:28 2021) by ROOT version 6.24/02
 
   myCanvas->cd();
   // make sure all pads are transparent
   TIter it(myCanvas->GetListOfPrimitives());
   TObject* obj;
   while ((obj = it())) {
      if (obj->InheritsFrom("TPad")) {
         dynamic_cast<TPad*>(obj)->SetFillStyle(4000);
      }
   }
   TLatex*    tex = ::new TLatex(0.02336825, 0.86604, "Max. Z");
   tex->SetTextSize(0.028);
   tex->SetTextAngle(90);
   tex->SetLineWidth(2);
   tex->Draw();
   TMarker* marker = ::new TMarker(0.01611604, 0.8519389, 4);
 
   Int_t ci;      // for color index setting
   TColor* color; // for color definition with alpha
   ci = TColor::GetColor("#333399");
   marker->SetMarkerColor(ci);
   marker->SetMarkerStyle(4);
   marker->SetMarkerSize(1.7);
   marker->Draw();
   marker = ::new TMarker(0.01772764, 0.4453584, 30);
 
   ci = TColor::GetColor("#ffcc33");
   marker->SetMarkerColor(ci);
   marker->SetMarkerStyle(30);
   marker->SetMarkerSize(1.7);
   marker->Draw();
   tex = ::new TLatex(0.02497985, 0.4629847, "Max. Z with identified A");
   tex->SetTextSize(0.028);
   tex->SetTextAngle(90);
   tex->SetLineWidth(2);
   tex->Draw();
   tex = ::new TLatex(0.02659146, 0.2632197, "Mean Z");
   tex->SetTextSize(0.028);
   tex->SetTextAngle(90);
   tex->SetLineWidth(2);
   tex->Draw();
   tex = ::new TLatex(0.02820306, 0.1045828, "Min. Z");
   tex->SetTextSize(0.028);
   tex->SetTextAngle(90);
   tex->SetLineWidth(2);
   tex->Draw();
   marker = ::new TMarker(0.01853344, 0.2479436, 20);
 
   ci = TColor::GetColor("#333399");
   marker->SetMarkerColor(ci);
   marker->SetMarkerStyle(20);
   marker->SetMarkerSize(1.7);
   marker->Draw();
   marker = ::new TMarker(0.02095085, 0.08578143, 25);
 
   ci = TColor::GetColor("#333399");
   marker->SetMarkerColor(ci);
   marker->SetMarkerStyle(25);
   marker->SetMarkerSize(1.7);
   marker->Draw();
}