Ignoring if() statement

ROOT
1

Hi everyone,

I have been working on a root macro for quite some time. The macro is designed to get data produced by a Monte Carlo simulation from a root file and plot certain tests to determine the effects of different geometric elements. My macro has been running great until I decided to make it a bit more “user friendly” by providing some options such as omitting certain histos when they aren’t needed. However, for some reason root is ignoring my logical statements, particularly this one:

1if ((S1>0) && (sYes==1)) { 
2	Neff+=1; 
3	for (int z=0; z<nBins; z++) {
4	  S2[z]=S2[z]/S1;
5	  s2s1->Fill((rtop-(rMin+(z+.05)*binSize)), S2[z]); 
6	} 
7}

The details of the elements are as follows:
s2s1 is a TH1D histogram

Neff is an integer that increments to keep track of the number of events that make the cuts for s2s1

S1 and S2 are doubles determined earlier in the program

sYes is a bool that if true will enable all the necessary parts to create s2s1, but if false will ignore them to make the program run faster

The problem arises when sYes is set to 0. This error occurs:

Error: non class,struct,union object $s2s1 used with . or -> folder/file.C:#:
*** Interpreter error recovered ***

According to the logic, it should not even be executing the line where s2s1 is filled. However, when sYes is set to 1 the whole thing runs fine. I placed cout statements all throughout my program and have found that Neff is not incremented on line 2 of the provided code and that S2[z] does not appear to change. Also, when I leave sYes false and comment out line 5 of the given code the program runs as expected. Can anybody tell me why ROOT is ignoring the if() statement and trying to execute line 5 anyways?

2

Hi,

I cannot reproduce that. Can you simplify your script so much that I can run it (ideally without ROOT files etc) and it still shows the issue?

Cheers, Axel.

3

Well that was only a part of it. Similar if statements in my script posed no problem. I guess my question is whether there are any other cases will ignore an if statement simply because filling a histogram takes priority, or some other similar thing. If this is not the case, here is my entire code:

void plots(int Nev=1000, double binSize=0.5, int rMin=0, char *fName="./comparison_tests/test_internal_gammas/test_internal_gammas.root", double ringsize=2, const bool radYes=0, const bool ravYes=0, const bool sYes=1, const bool centYes=1) {

  //load library definition
  gSystem->Load("./tools/ConfigHandler/lib/libEventData.so");
  gSystem->Load("./tools/ConfigHandler/lib/libConfigHandler.so");

  //open file
  TFile* fdata = new TFile(fName);

  //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  // CONFIGURATION
  //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  TTree* config = (TTree*)(fdata->Get("ConfigOptions"));

  //find the vessel size, livetime, mass, etc.
  char acx[30]="", acy[30]="", acz[30]="";
  //char acxx[30]="", acyy[30]="", aczz[30]="";
  char lar_top_d[30]="", lar_bot_d[30]="", lar_h[30]="";
  char ppd[30]="", pgen[30]="";
  config->SetBranchAddress("active_center_z",acz);
  config->SetBranchAddress("active_center_x",acx);
  config->SetBranchAddress("active_center_y",acy);
  config->SetBranchAddress("lar_top_d",lar_top_d);
  config->SetBranchAddress("lar_bot_d",lar_bot_d);
  config->SetBranchAddress("lar_h",lar_h);
  config->GetEntry(0);
  
  const Double_t active_center_z = atof(acz);
  const Double_t active_center_x = atof(acx);
  const Double_t active_center_y = atof(acy);
  const double rtop = atof(lar_top_d)/2.;
  const double rbot = atof(lar_bot_d)/2.;
  const double h = atof(lar_h);

  const int overBin= (int)(1/binSize+0.5); //1 over binSize
  const int nBins=rtop*overBin+1; //number of bins used to split the radius of the LAr chamber
  const double printBins=(rtop-rMin)*overBin; //number of bins printed 
  const int lessRings = nBins-(1+ringsize*overBin); //determines number of bins to cut off to simulate upper ring
  const double maxR=rtop; //outer edge of detector, or outer radius to be studied for average energy deposits in a radius

  //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  // Prepare Histogram
  //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  if (radYes==1) {
    //radgen plots the average energy deposit made by a particle in each bin (dr)
    TH1D *radgen  = new TH1D("radgen", "Average Energy Deposit in dr (for all events)", printBins, rMin, rtop);
    radgen->GetXaxis()->SetNdivisions(rtop/2);
    radgen->GetXaxis()->CenterTitle();
    radgen->GetXaxis()->SetTitle("Radius in cm");
    radgen->GetYaxis()->CenterTitle();
    radgen->GetYaxis()->SetTitle("Average Energy Deposit in keV");
  }

  if (ravYes==1) {
    //ravgen plots the average accumulation of edeps looping through the radius (for all events)
    TH1D *ravgen  = new TH1D("ravgen", "Average Total Energy Deposit at Radius R", printBins,rMin,rtop);
    ravgen->GetXaxis()->SetNdivisions(rtop/2);
    ravgen->GetXaxis()->CenterTitle();
    ravgen->GetXaxis()->SetTitle("Radius in cm");
    ravgen->GetYaxis()->CenterTitle();
    ravgen->GetYaxis()->SetTitle("Accumulation of Edeps in keV");
  }  

  if (sYes==1) {
    //s2s1 plots the average ratio of s2/s1 with increasing ring sizes
    TH1D *s2s1  = new TH1D("s2s1", "S2/S1 vs. Rrid Width", printBins, rMin, rtop);
    s2s1->GetXaxis()->SetNdivisions(rtop/2);
    s2s1->GetXaxis()->CenterTitle();
    s2s1->GetXaxis()->SetTitle("Grid Width in cm");
    s2s1->GetYaxis()->CenterTitle();
    s2s1->GetYaxis()->SetTitle("Average S2/S21");
  }

  if (centYes==1) {
    //centroid plots the s2/s1 ratio as a function of R_centroid on an event by event basis
    TH2D *centroid = new TH2D("centroid", "s2/s1 vs. R_centroid per Event", 100, 0, rtop, 100, 0, 1.1);
    centroid->GetXaxis()->CenterTitle();
    centroid->GetXaxis()->SetTitle("r centroid in cm");
    centroid->GetYaxis()->CenterTitle();
    centroid->GetYaxis()->SetTitle("s1/s2");
  }

  //get data tree from file
  TTree *Events = (TTree*)(fdata->Get("Events"));
  EventData* evt = 0;
  Events->SetBranchAddress("particle_data",&evt);

  //Define 
  double S1=0;
  double S2[nBins];
  int Neff=0; //number of events used to average s2s1
  double s2=0; //these are used for histogram of s1/s2 vs r_centroid 
  double ratio=0; //this is the ratio of s2/s1 used in centroid

  double En[7][nBins]; //stores edep information as well as particle counts in radial bins
  for(int q=0; q<7; q++) //fill En with zeros
    {
      if (q==0) {
	for (int W=0; W<nBins; W++) {
	  En[q][W]=(W+1)*binSize;
	}
      }
      if (q>0){
	for (int w=0; w<nBins; w++) {
	  En[q][w]=0;
	}
      }
    }
  
  
  //Ensure no repeated events by only processing the number of events in the tree
  if (Nev>Events->GetEntries()){
    std::cout<<"Too Many Events Requested"<<endl<<"Changing Number Of Events To "<<Events->GetEntries()<<" "<<endl;
    Nev=Events->GetEntries();
  }
  cout << "This tree contains " << Events->GetEntries() << " events; only print " << Nev << endl;
  
  //loop through NeV events in the data tree
  for(int event = 0; event < Nev; event++)
    {
      //print every once in a while, just so I can see that the code is running
      if((event%5000)==0)cout << "At event " << event << "." << endl;
      
      //get the next event
      Events->GetEntry(event);
      
      //define variables pertaining to this specific event and set them equal to zero
      double r;
      S1=0;
      for (int k=0; k<nBins; k++) S2[k]=0;
      double EdepActive=0, EdepTot=0;
      double xCentroid=0, yCentroid=0, zCentroid=0;
      double rCentroid=0;                                                                           
      ratio=0;
      
      for(int j=0; j<evt->n_particles; j++) //loop through the events and get info from each
 	{
	  r=sqrt((evt->spawn_x.at(j))*(evt->spawn_x.at(j)) + (evt->spawn_y.at(j))*(evt->spawn_y.at(j)));
	  EdepActive+=evt->edep_active.at(j);
	  xCentroid+=evt->spawn_x.at(j)*evt->edep_active.at(j);
	  yCentroid+=evt->spawn_y.at(j)*evt->edep_active.at(j);
	  zCentroid+=evt->spawn_z.at(j)*evt->edep_active.at(j);
	  
	  for (int k=0; k<nBins; k++) //loops through R in binSize incriments and places edep data in En and S2 arrays
	    {
	      if ((r>=rMin) && (r>=rMin+(k*binSize)) && (r<rMin+(k+1)*binSize)) {
		En[1][k]+=evt->edep_active.at(j); //store energy deposited at this radius
		En[2][k]+=1;                //store number of events that have occured in dr thus far 
	      }
	      
	      if (r<rMin+(k+1)*binSize) {
		En[4][k]+=evt->edep_active.at(j); //accumulates edeptotals as dr increases
		En[5][k]+=1;     //keeps track of the number of events that have occured by this radius
		S2[k]+=evt->edep_active.at(j); //records S2 signal if ringsize was to this bin
	      }
	    }
	}
      
      if(EdepActive>0) { //calculate rCentroid
	xCentroid /= EdepActive;
	yCentroid /= EdepActive;
	zCentroid /= EdepActive;
	rCentroid=sqrt(xCentroid*xCentroid+yCentroid*yCentroid);
      }
      
      S1=S2[nBins-1]; //S1 is total edep in LAr chamber
      
      if ((centYes==1) && (S1!=0) && (EdepActive>0)) {
	ratio=S2[lessRings]/S1;
	centroid->Fill(rCentroid, ratio);
      }
      
      if ((S1>0) && (sYes==1)) { 
	Neff+=1; 
	for (int z=0; z<nBins; z++) {
	  S2[z]=S2[z]/S1;
	  s2s1->Fill((rtop-(rMin+(z+.05)*binSize)), S2[z]); 
	} 
      }   
    }    
  
  if (sYes==1) s2s1->Scale(1./Neff);
  
  if ((radYes==1) || (ravYes==1)) {
    for (int J=0; J<=(nBins-1); J++) { //loop over J and produce the En[3][J] which is the average energy
      //deposited in that bin for all events
      if (radYes==1)  En[3][J]=En[1][J]/Nev;
      if (ravYes==1) En[6][J]=En[4][J]/Nev;
    }

    for (int n=0; n<=((printBins+rMin*overBin)-1); n++) {
      if (radYes==1) radgen->Fill(n*binSize, En[3][n]);
      if (ravYes==1) ravgen->Fill(n*binSize, En[6][n]);
    }
  }

  TFile* out = new TFile("splotter.root", "recreate");
  if (radYes==1) radgen->Write();
  if (ravYes==1) ravgen->Write();
  if (centYes==1)centroid->Write();
  if (sYes==1) s2s1->Write();
 
  if (radYes==1) {
    if (gROOT->FindObject("c1")) delete gROOT->FindObject("c1");
    TCanvas* c1=new TCanvas("c1");
    radgen->SetTitle("Radial Energy Distribution");
    radgen->SetLineColor(1);
    radgen->Draw();
  }

  if (ravYes==1) {
    if (gROOT->FindObject("c2")) delete gROOT->FindObject("c2"); 
    TCanvas* c2=new TCanvas("c2");
    ravgen->SetTitle("Radial Energy Distribution");
    ravgen->SetLineColor(1);
    ravgen->Draw();
  }
  
  if (sYes==1) {
    if (gROOT->FindObject("c3")) delete gROOT->FindObject("c3");
    TCanvas* c3=new TCanvas("c3");
    s2s1->SetLineColor(1);
    s2s1->GetXaxis()->SetRangeUser(0,5);
    s2s1->Draw();
  }

  if (centYes==1) {
    if (gROOT->FindObject("c4")) delete gROOT->FindObject("c4");
    TCanvas* c4=new TCanvas("c4");
    centroid->Draw();
  }
}
4

Hi,

[quote]Error: non class,struct,union object $s2s1 used with . or -> folder/file.C:#:
*** Interpreter error recovered ***[/quote]This is likely caused only by the parsing of the clause (even if it will not be executed) and by the fact that s2s1 in indeed NOT declared (as far as C++ is concerned) at this point in the process. What you have is:if (sYes==1) { //s2s1 plots the average ratio of s2/s1 with increasing ring sizes TH1D *s2s1 = new TH1D("s2s1", "S2/S1 vs. Rrid Width", printBins, rMin, rtop); ... } which means that according to C++ scoping rules, the variable s2s1 is only visible within that if statement. When sYes==1, then the histogram is created and the line that provokes the error message does work properly but it relies on an extension provided by ROOT and CINT for interpreted code where you can use the ‘name’ of a known entity if this entity is held in the current ROOT directory (i.e. the line s2s1->Fill((rtop-(rMin+(z+.05)*binSize)), S2[z]); is turned by CINT into: TH1D *s2s1; gDirecotry->GetObject("s2s1",s2s1); s2s1->Fill((rtop-(rMin+(z+.05)*binSize)), S2[z]);

If you attempted to compiled your code, it would fail with a compilation error.

To solve both, simply do:TH1D *s2s1 = 0; if (sYes==1) { //s2s1 plots the average ratio of s2/s1 with increasing ring sizes s2s1 = new TH1D("s2s1", "S2/S1 vs. Rrid Width", printBins, rMin, rtop); ... }

Cheers,
Philippe.