sp_evaluation.C

float totalChargeWithBin(TH2F* h2, int ix, int ymin=0, int ymax=6000){
  float maxCharge = 0;
  float totCharge = 0;
  if (ymax>h2->GetNbinsY()) ymax = h2->GetNbinsY();
  for(int iy=ymin; iy<=ymax; iy++){
    float content = h2->GetBinContent(ix, iy);
    if (content>0) totCharge += content;
    if (content>maxCharge) maxCharge = content;
  }

  if(maxCharge<1000 or totCharge<5000) return 0;
  else 
    return totCharge;
}

void sp_evaluation(const char* filename="protodunevd-sim-check-deposplat.root", const char* plane="u"){

  auto file = TFile::Open(filename);

  // cru0
  auto htru0 = (TH2F*)file->Get(Form("h%s_ductor0",plane));
  auto hgau0 = (TH2F*)file->Get(Form("h%s_gauss0",plane));
  auto hdnn0 = (TH2F*)file->Get(Form("h%s_dnnsp0",plane));
  // cru1
  auto htru1 = (TH2F*)file->Get(Form("h%s_ductor1",plane));
  auto hgau1 = (TH2F*)file->Get(Form("h%s_gauss1",plane));
  auto hdnn1 = (TH2F*)file->Get(Form("h%s_dnnsp1",plane));

  vector<TH2F*> htrus {htru0, htru1};
  vector<TH2F*> hgaus {hgau0, hgau1};
  vector<TH2F*> hdnns {hdnn0, hdnn1};

  auto hRatioGau = new TH1F("hRatioGau","",100,0.1,1.9);
  auto hRatioDnn = new TH1F("hRatioDnn","",100,0.1,1.9);
  float ntru=0, ngau=0, ndnn=0; // counter for truth labels, bad gauss channels, bad dnn channels


  for(int ih=0; ih<2; ih++){
      int nx = htrus[ih]->GetNbinsX();
      int ny = htrus[ih]->GetNbinsY();
      for(int ix=1; ix<=nx; ix++){
        int chan = htrus[ih]->GetXaxis()->GetBinCenter(ix);
        float truCharge = totalChargeWithBin(htrus[ih], ix);
        if (truCharge<5000) continue;

        ntru ++;

        int binGau = hgaus[ih]->GetXaxis()->FindBin(chan);
        if (binGau>0 and binGau<=hgaus[ih]->GetNbinsX() ){
          float gauCharge = totalChargeWithBin(hgaus[ih], binGau);
          if (gauCharge/truCharge>0.1 and gauCharge/truCharge<1.9)
            hRatioGau->Fill(gauCharge/truCharge);
          else
            ngau ++;
        }
        else{
          ngau ++;
        }

        int binDnn = hdnns[ih]->GetXaxis()->FindBin(chan);
        if (binDnn>0 and binDnn<=hdnns[ih]->GetNbinsX() ){
          float dnnCharge = totalChargeWithBin(hdnns[ih], binDnn);
          if(dnnCharge/truCharge>0.1 and dnnCharge/truCharge<1.9)
            hRatioDnn->Fill(dnnCharge/truCharge);
          else{
            ndnn ++;
            // cout << "channel: " << chan << endl;
          }
        }
        else{
          ndnn ++;
            // cout << "channel @2: " << chan << endl;
        }
        
      }


  } // end of looping cru


  hRatioGau->SetLineWidth(2);
  hRatioGau->SetLineColor(kBlack);
  hRatioDnn->SetLineWidth(2);
  hRatioDnn->SetLineColor(kBlue);

  hRatioDnn->Draw();
  hRatioGau->Draw("same");


  cout << "=== Bias " << endl;
  cout << "Gau bias (%): " <<  100.0*(hRatioGau->GetMean() - 1.0) << endl;
  cout << "Dnn bias (%): " <<  100.0*(hRatioDnn->GetMean() - 1.0) << endl;

  cout << "=== Resolution " << endl;
  cout << "Gau RMS (%): " <<  100.0*hRatioGau->GetRMS()/hRatioGau->GetMean() << endl;
  cout << "Dnn RMS (%): " <<  100.0*hRatioDnn->GetRMS()/hRatioDnn->GetMean() << endl;

  cout << "=== Inefficiency " << endl;
  cout << "ntru: " << ntru << ", bad ndnn: " << ndnn << ", bad ngau: " << ngau << endl;
  float prob = ngau/ntru; 
  cout << "Gau inefficiency (%): " <<  100.0*prob << " +/- " << 100.0* sqrt(prob*(1-prob)/ntru) << endl;
  prob = ndnn/ntru;
  cout << "Dnn inefficiency (%): " <<  100.0*prob << " +/- " << 100.0* sqrt(prob*(1-prob)/ntru) << endl;
}