Dear ROOT guys,
in a 2D distribution, I am trying to find the contour level that will cut away, say, 10% of the data. To achieve this, I start from the FirstContour.C tutorial example that shows how to retrieve the first contour of a histogram contour plot. Now I simply want to repeat this procedure by raising the level of the first contour step by step. But as soon as I repeat the procedure to retrieve the contour and use it to generate a cut in a loop, ROOT produces a segmentation violation (somewhere in TSelectorDraw::ClearFormula it seems).
How would I have to adapt my code so that I can safely run the procedure in a loop?
Many thanks in advance!
Here’s my macro, based on $ROOTSYS/tutorials/hist/FirstContour.C
Note: to run it, you should be in the same directory, or have hsimple.C in the …/ directory!
void FirstContour()
{
//this macro generates a color contour plot by selecting entries
//from an ntuple file.
//The TGraph object corresponding to the first contour line is
//accessed and displayed into a separate canvas.
//Author: Rene Brun
TString dir = gSystem->UnixPathName(TCint::GetCurrentMacroName());
dir.ReplaceAll("FirstContour.C","../hsimple.C");
dir.ReplaceAll("/./","/");
if (!gInterpreter->IsLoaded(dir.Data())) gInterpreter->LoadMacro(dir.Data());
TFile *file = (TFile*)gROOT->ProcessLineFast("hsimple(1)");
if (!file) return;
TTree *ntuple = (TTree*)file->Get("ntuple");
TCanvas *c1 = new TCanvas("c1","Contours",10,10,800,600);
gStyle->SetPalette(1);
TH2D* histo = new TH2D("histo","histo", 50, -5, 5, 50, -5, 5);
Double_t n0 = ntuple->GetEntries();
ntuple->Draw("py:px>>histo","px*px+py*py < 20");
Double_t conts[3];
conts[0] = 20;
conts[1] = 60;
conts[2] = 100;
// for (Int_t i=10;i<20;i+=10) { // run it once: no problem
for (Int_t i=10;i<100;i+=10) { // run it several times: crash!
conts[0] = i;
histo->SetContour(3,conts);
histo->Draw("contz,list");
//we must call Update to force the canvas to be painted. When
//painting the contour plot, the list of contours is generated
//and a reference to it added to the Root list of special objects
c1->Update();
TObjArray *contours =
(TObjArray*)gROOT->GetListOfSpecials()->FindObject("contours");
TList *lcontour1 = (TList*)contours->At(0);
TGraph *gc1 = (TGraph*)lcontour1->First();
gc1->SetMarkerStyle(21);
gc1->Draw("alp");
//We make a TCutG object with the array obtained from this graph
TCutG *cutg = new TCutG("cutg",gc1->GetN(),gc1->GetX(),gc1->GetY());
cutg->SetVarX("px");
cutg->SetVarY("py");
TH2D *histo2 = new TH2D("histo2","histo2", 50, -5, 5, 50, -5, 5);
ntuple->Draw("py:px>>histo2",cutg->GetName());
Double_t ncut = histo2->GetEntries();
if (cutg) delete cutg;
if (histo2) delete histo2;
Double_t frac = ncut/n0;
cout << " i = " << i << " n(cut) = " << ncut << " fraction = " << frac << endl;
}
} 
I am working on it.