Using TTreeReader vs RDataFrame

I eventually went back to using TTreeReader because of a problem with using RDataFrame. In the following script, I comment out the Foreach function call, whicih is the part that seems to be causing the problem.

RDataFrame Foreach

double count_passed_events_v2(TString file_prefix, int num_jobs){

  // initialize TChain and loop through all subdirectories, adding root files to chain.
  TChain chain("passTree");

  for(int i=1; i<=num_jobs; i++){

    TString file_name ("");
    file_name.Form(file_prefix + "/run%d/IceFinal_%d_passTree1.root", i ,i);

    chain.Add(file_name); 
  }

  RDataFrame df(chain);
  
  double weighted_passed_count=0.;
  // lambda expression that computes the weighted count 
  auto count = [&weighted_passed_count] (const double &w1, const double &w2, const double &w3){
    weighted_passed_count += w1/(w2*w3);
  };
  
  // df.Foreach(count, {"event.neutrino.path.weight",
  //                     "event.loop.positionWeight",
  //                     "event.loop.directionWeight"});
  return weighted_passed_count;

}

As soon as this is uncommented, the script can only be run once. That is, it will work fine the first time when I type root effar.C where effar.C calls count_passed_events_v2, but in the ROOT interactive mode, if I try to execute the script again using .x effar.C, this script leads to a segmentation fault. It seems like the seg-fault has to do with the lambda expression, but I am not sure.

TODO

Ask about this on ROOT Forum.

Below is the full script in which, instead of RDataFrame, I used TTreeReader, which does not seg-fault.

TTreeReader reader.Next()

double count_passed_events(TString file_prefix, int num_jobs){

  // initialize TChain and loop through all subdirectories, adding root files to chain.
  TChain *my_chain = new TChain("passTree");

  for(int i=1; i<=num_jobs; i++){

    TString file_name ("");
    file_name.Form(file_prefix + "/run%d/IceFinal_%d_passTree1.root", i ,i);

    my_chain->Add(file_name); 
  }
  
  TTreeReader reader(my_chain);

  TTreeReaderValue<double_t> w1(reader, "event.neutrino.path.weight");
  TTreeReaderValue<double_t> w2(reader, "event.loop.positionWeight");
  TTreeReaderValue<double_t> w3(reader, "event.loop.directionWeight");

  double weighted_passed_event_count;
  while (reader.Next()) {
    weighted_passed_event_count += *w1/(*w2 * *w3); 
  }

  delete my_chain;
  return weighted_passed_event_count;
}

void effar(){
  
  const double ice_volume = 26859227.0607;
  const double prefactor = ice_volume * 4 * M_PI;
  const int e21_interaction_length = 141324./1e3;
  const int e20_interaction_length = 267669./1e3;
  const int e19_interaction_length = 528717./1e3;
  const int e18_interaction_length = 1.10053e+06/1e3;

  int num_jobs;
  int nnt;
  int energy;
  TString file_path("");

  nnt = 100;
  energy = 21;
  file_path.Form("../2023-12-21_nnt_%d_energy_%d",nnt, energy);
  num_jobs = count_files(file_path);
  double vol21 = prefactor * count_passed_events(file_path, num_jobs) / (nnt * num_jobs);
  
  nnt = 100;
  energy = 20;
  file_path.Form("../2023-12-21_nnt_%d_energy_%d",nnt, energy);
  num_jobs = count_files(file_path);
  double vol20 = prefactor * count_passed_events(file_path, num_jobs) / (nnt * num_jobs);

  nnt = 400;
  energy = 19;
  file_path.Form("../2023-12-21_nnt_%d_energy_%d",nnt, energy);
  num_jobs = count_files(file_path);
  double vol19 = prefactor * count_passed_events(file_path, num_jobs) / (nnt * num_jobs);

  nnt = 1000;
  energy = 18;
  file_path.Form("../2023-12-21_nnt_%d_energy_%d",nnt, energy);
  num_jobs = count_files(file_path);
  double vol18 = prefactor * count_passed_events(file_path, num_jobs) / (nnt * num_jobs);


  // --------------   Data  --------------  //
  double effar[4] = {vol21/(e21_interaction_length), vol20/(e20_interaction_length),
                     vol19/(e19_interaction_length), vol18/(e18_interaction_length) };
  double e[4] = {1e21, 1e20, 1e19, 1e18};
  
  double effar_old[8] = {0.0030019, 0.2665, 4.45, 18.89, 41.3415, 77.749,122.65, 169.992};
  double e_old[8] = {3.17e+17, 1.01e+18, 3.21e+18, 1.009e+19, 3.205e+19, 1.017e+20, 1.19e+20, 1.013e+21};

  // --------------  Drawing --------------  //
  
  TCanvas * c = new TCanvas("","",1400,900);
  c->SetLogy();
  c->SetLogx();

  auto mg = new TMultiGraph();
  mg->SetTitle("Effective Area Comparison");

  TGraph *gr = new TGraph(8, e_old, effar_old);
  gr->SetMarkerStyle(kFullCircle);
  gr->SetMarkerColorAlpha(kBlue,0.35);
  mg->Add(gr);

  TGraph *gr_2 = new TGraph(4, e, effar);
  gr_2->SetMarkerStyle(kFullCross);
  gr_2->SetMarkerSize(1.5);
  gr_2->SetMarkerColorAlpha(46,0.35);
  mg->Add(gr_2);
  

  mg->GetXaxis()->SetTitle("energy [eV]");
  mg->GetXaxis()->SetTitleOffset(1.1);
  mg->GetYaxis()->SetTitle("effective area [km^{2} sr]");
  mg->Draw("ALP");

  TLegend * leg = new TLegend(0.2,0.75,0.4,0.9);
  leg->AddEntry(gr, "previous data");
  leg->AddEntry(gr_2, "data from my test run");
  leg->Draw();

  TString of_name = "../img/effective_area_plot.png";
  c->SaveAs(of_name);
}

It might be just in my head, but I think using TTreeReader is slightly faster.

TODO

Time this.