<commit_before>// This tutorial shows how to optimize a lens shape by using the ROOT Minuite2
// optimization engine.

#include "Math/Functor.h"
#include "Minuit2/Minuit2Minimizer.h"
#include "TCanvas.h"
#include "TGeoBBox.h"
#include "TGeoMatrix.h"
#include "TGeoPhysicalNode.h"
#include "TH2D.h"

#include "AFocalSurface.h"
#include "AGeoAsphericDisk.h"
#include "AGlassCatalog.h"
#include "ALens.h"
#include "AMirror.h"
#include "AOpticsManager.h"
#include "ARayArray.h"
#include "ARayShooter.h"

static const Double_t mm = AOpticsManager::mm();
static const Double_t um = AOpticsManager::um();
static const Double_t nm = AOpticsManager::nm();

static AOpticsManager* gOpticsManager = 0;
static const int kNtheta = 3;
static ARayArray* gArray[kNtheta] = {0, 0, 0};
static double kTheta[kNtheta] = {0., 2.5, 5.0};
static const int kNlambda = 3;
static const double kRadius = 5.0 * mm;

double GetSpotSize(ARayArray* array) {
  // Returns standard deviation of a spot
  TObjArray* focused = array->GetFocused();
  TH2D h("", "", 1, 0, 0, 1, 0, 0);

  for (int i = 0; i < focused->GetLast(); i++) {
    ARay* ray = (ARay*)(*focused)[i];
    double p[4];
    ray->GetLastPoint(p);
    h.Fill(p[0], p[1]);
  }  // i

  double sx = h.GetStdDev(1);
  double sy = h.GetStdDev(2);

  return TMath::Sqrt(sx * sx + sy * sy);
}

int Nphotons(ARayArray* array) {
  int total = 0;
  total += array->GetFocused()->GetLast() + 1;
  total += array->GetStopped()->GetLast() + 1;
  total += array->GetSuspended()->GetLast() + 1;
  total += array->GetAbsorbed()->GetLast() + 1;

  return total;
}

double Func(const double* par) {
  static bool initialized = kFALSE;
  static TGeoRotation rot[kNtheta];
  static TGeoTranslation tr("tr", 0, 0, -3 * mm);

  if (!initialized) {
    for (int i = 0; i < kNtheta; i++) {
      rot[i] = TGeoRotation();
      rot[i].SetAngles(0, kTheta[i], 0.);
    }  // i
    initialized = kTRUE;
  }  // if

  if (gOpticsManager) {
    delete gOpticsManager;
  }  // if
  gOpticsManager = new AOpticsManager("manager", "manager");
  gOpticsManager->DisableFresnelReflection(kTRUE);

  // Make the world
  TGeoBBox* box = new TGeoBBox("box", 100 * mm, 100 * mm, 100 * mm);
  AOpticalComponent* top = new AOpticalComponent("top", box);
  gOpticsManager->SetTopVolume(top);

  AGeoAsphericDisk* disk = new AGeoAsphericDisk("disk", 0 * mm, 0., par[0],
                                                par[1], kRadius, 0. * mm);
  disk->SetConicConstants(0, par[2]);
  double f1 = disk->CalcF1(kRadius);
  double f2 = disk->CalcF2(kRadius);
  if (f2 - f1 < 0) {  // negative edge thickness
    return 1e100;
  }  // if

  ALens* lens = new ALens("lens", disk);
  AGlassCatalog schott("../misc/schottzemax-20180601.agf");
  lens->SetRefractiveIndex(schott.GetRefractiveIndex("N-BK7"));
  gOpticsManager->GetTopVolume()->AddNode(lens, 1);

  double origin[3] = {0, 0, 50 * mm + 1 * um};
  TGeoBBox* box2 = new TGeoBBox("box2", 10 * mm, 10 * mm, 1 * um, origin);
  AFocalSurface* screen = new AFocalSurface("screen", box2);
  top->AddNode(screen, 1);

  AGeoAsphericDisk* disk2 = new AGeoAsphericDisk("disk2", 0 * mm, 0., par[0],
                                                 0., kRadius * 1.2, kRadius);
  AObscuration* obs = new AObscuration("obs", disk2);
  gOpticsManager->GetTopVolume()->AddNode(obs, 1);
<commit_msg>Set the extinction coefficients of N-BK7 to zero (NULL) to avoid minimization failure probably due to internal photon absorption.<commit_after>// This tutorial shows how to optimize a lens shape by using the ROOT Minuite2
// optimization engine.

#include "Math/Functor.h"
#include "Minuit2/Minuit2Minimizer.h"
#include "TCanvas.h"
#include "TGeoBBox.h"
#include "TGeoMatrix.h"
#include "TGeoPhysicalNode.h"
#include "TH2D.h"

#include "AFocalSurface.h"
#include "AGeoAsphericDisk.h"
#include "AGlassCatalog.h"
#include "ALens.h"
#include "AMirror.h"
#include "AOpticsManager.h"
#include "ARayArray.h"
#include "ARayShooter.h"

static const Double_t mm = AOpticsManager::mm();
static const Double_t um = AOpticsManager::um();
static const Double_t nm = AOpticsManager::nm();

static AOpticsManager* gOpticsManager = 0;
static const int kNtheta = 3;
static ARayArray* gArray[kNtheta] = {0, 0, 0};
static double kTheta[kNtheta] = {0., 2.5, 5.0};
static const int kNlambda = 3;
static const double kRadius = 5.0 * mm;

double GetSpotSize(ARayArray* array) {
  // Returns standard deviation of a spot
  TObjArray* focused = array->GetFocused();
  TH2D h("", "", 1, 0, 0, 1, 0, 0);

  for (int i = 0; i < focused->GetLast(); i++) {
    ARay* ray = (ARay*)(*focused)[i];
    double p[4];
    ray->GetLastPoint(p);
    h.Fill(p[0], p[1]);
  }  // i

  double sx = h.GetStdDev(1);
  double sy = h.GetStdDev(2);

  return TMath::Sqrt(sx * sx + sy * sy);
}

int Nphotons(ARayArray* array) {
  int total = 0;
  total += array->GetFocused()->GetLast() + 1;
  total += array->GetStopped()->GetLast() + 1;
  total += array->GetSuspended()->GetLast() + 1;
  total += array->GetAbsorbed()->GetLast() + 1;

  return total;
}

double Func(const double* par) {
  static bool initialized = kFALSE;
  static TGeoRotation rot[kNtheta];
  static TGeoTranslation tr("tr", 0, 0, -3 * mm);

  if (!initialized) {
    for (int i = 0; i < kNtheta; i++) {
      rot[i] = TGeoRotation();
      rot[i].SetAngles(0, kTheta[i], 0.);
    }  // i
    initialized = kTRUE;
  }  // if

  if (gOpticsManager) {
    delete gOpticsManager;
  }  // if
  gOpticsManager = new AOpticsManager("manager", "manager");
  gOpticsManager->DisableFresnelReflection(kTRUE);

  // Make the world
  TGeoBBox* box = new TGeoBBox("box", 100 * mm, 100 * mm, 100 * mm);
  AOpticalComponent* top = new AOpticalComponent("top", box);
  gOpticsManager->SetTopVolume(top);

  AGeoAsphericDisk* disk = new AGeoAsphericDisk("disk", 0 * mm, 0., par[0],
                                                par[1], kRadius, 0. * mm);
  disk->SetConicConstants(0, par[2]);
  double f1 = disk->CalcF1(kRadius);
  double f2 = disk->CalcF2(kRadius);
  if (f2 - f1 < 0) {  // negative edge thickness
    return 1e100;
  }  // if

  ALens* lens = new ALens("lens", disk);
  AGlassCatalog schott("../misc/schottzemax-20180601.agf");
  auto bk7 = schott.GetRefractiveIndex("N-BK7");
  // Temporarily remove the coefficients because the optimization does not
  // converge well when they are non-zero, resulting in random internal
  // absorption in the lens. A better workaround is needed.
  bk7->SetExtinctionCoefficient(0);
  lens->SetRefractiveIndex(bk7);
  gOpticsManager->GetTopVolume()->AddNode(lens, 1);

  double origin[3] = {0, 0, 50 * mm + 1 * um};
  TGeoBBox* box2 = new TGeoBBox("box2", 10 * mm, 10 * mm, 1 * um, origin);
  AFocalSurface* screen = new AFocalSurface("screen", box2);
  top->AddNode(screen, 1);

  AGeoAsphericDisk* disk2 = new AGeoAsphericDisk("disk2", 0 * mm, 0., par[0],
                                                 0., kRadius * 1.2, kRadius);
  AObscuration* obs = new AObscuration("obs", disk2);
  gOpticsManager->GetTopVolume()->AddNode(obs, 1);

{'fraction_non_alphanumeric': 0.15650319829424306, 'fraction_numerical': 0.029992892679459843, 'mean_word_length': 3.5305859626529297, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 48, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '313219', 'n_tokens_mistral': 2759, 'n_tokens_neox': 2569, 'n_words': 672}