Interval of intersection doesn't work

David_Fernandes · May 12, 2020, 3:25pm

Please read tips for efficient and successful posting and posting code

_ROOT Version: 6.18.00
_Platform: UBUNTU 18.04

Good Afternoon!

I expected that the code give me only the data in the material G4_WATER, which extended from 0.4 to 22.0, but gives me all the data in all volumes even when I restrict the position as show in the following code, what am I doing wrong?

//the goal is to do a histogram posxEn Th2f from tree

#include "TROOT.h"
#include "TMath.h"
#include "TColor.h"
#include "TCanvas.h"
#include "TFile.h"
#include "TTree.h"
#include "TGraph.h"
#include "TBrowser.h"
#include "TH1.h"
#include "TF1.h"
#include "TH2.h"
#include "TRandom.h"
#include "TStopwatch.h"
#include <iostream>
#include <fstream>
#include <iomanip>
#include <string>
#include <vector>
#include <array>

void posXY(){

//declaration of some counters

Int_t mm=0;
Int_t mm1=0;
Int_t mm2=0;
Int_t mm3=0;
Int_t mm4=0;
Int_t mm5=0;
Int_t mm6=0;
Int_t mm7=0;
Int_t mm8=0;

Double_t r=0;
Double_t dx=0;
Double_t dy=0;
Double_t dz=0;

Double_t max=0;
Double_t volP=0;

Double_t edep=0;
Double_t sumDep=0;
Double_t somaDep=0;

Double_t elost=0;
Double_t sumLos=0;
Double_t somaLos=0;

Double_t ekin=0;
Double_t sumKin=0;
Double_t SumEkin=0;

// declaration of TChain

      TChain chain("testeSA");
      
      chain.Add("testeSA_tree_1000.root");  [testeS_tree_1000.root|attachment](upload://87dqjXu41et80oiWvqxaMOEnvUZ.root) (2.7 MB) 
      
//////// *************** //////////////////////

//Declaration of leaves types
   Int_t           Event_EventID=0, nentries=0;
   Double_t        Event_AccumulatedEnergyDeposited=0,Event_InitialKineticEnergy=0, Event_AccumulatedEnergyLost=0;

   std::vector<double>  *Step_AccumulatedEnergyDeposited=0;
   std::vector<double>  *Step_AccumulatedEnergyLost=0;

   std::vector<double>  *Step_FinalKineticEnergy=0;
   std::vector<double>  *Step_FinalPosX=0;
   std::vector<double>  *Step_FinalPosY=0;
   std::vector<double>  *Step_FinalPosZ=0;
   
   std::vector<double>  *Step_InitialPosX=0;
   std::vector<double>  *Step_InitialPosY=0;
   std::vector<double>  *Step_InitialPosZ=0;
   std::vector<int>     *Step_StepNumber=0;
   std::vector<int>     *Step_TrackID=0;
   
   std::vector<int>     *Track_EventID=0;
   std::vector<int>     *Track_TrackID=0;
   std::vector<double>  *Track_AccumulatedDose=0;
   std::vector<double>  *Track_AccumulatedEnergyDeposited=0;
   std::vector<double>  *Track_AccumulatedEnergyLost=0;
   std::vector<double>  *Track_FinalKineticEnergy=0;

   std::vector<string>  *Step_CreatorProcess=0;
   std::vector<string>  *Step_FinalLogicalVolume=0;
   std::vector<string>  *Step_FinalProcess=0;
   std::vector<string>  *Step_InitialLogicalVolume=0;
   
   std::vector<string>  *Step_InitialMaterial=0;
   std::vector<string>  *Step_FinalMaterial=0;
   std::vector<string>  *Step_InitialProcess=0;
   std::vector<string>  *Step_Particle=0;
   
   std::vector<string>  *Track_FinalLogicalVolume=0;
   std::vector<string>  *Track_InitialMaterial=0;
   std::vector<string>  *Track_FinalMaterial=0;
   std::vector<string>  *Track_Particle=0;
   
// Set branches address 

   chain.SetBranchAddress("Event_EventID",&Event_EventID);
   chain.SetBranchAddress("Event_AccumulatedEnergyDeposited",&Event_AccumulatedEnergyDeposited);
   chain.SetBranchAddress("Event_InitialKineticEnergy",&Event_InitialKineticEnergy);
   chain.SetBranchAddress("Event_AccumulatedEnergyLost",&Event_AccumulatedEnergyLost);
   
   chain.SetBranchAddress("Step_AccumulatedEnergyDeposited",&Step_AccumulatedEnergyDeposited);
   chain.SetBranchAddress("Step_AccumulatedEnergyLost",&Step_AccumulatedEnergyLost);
   chain.SetBranchAddress("Step_FinalKineticEnergy",&Step_FinalKineticEnergy);
   
   chain.SetBranchAddress("Step_FinalPosX",&Step_FinalPosX);
   chain.SetBranchAddress("Step_FinalPosY",&Step_FinalPosY);
   chain.SetBranchAddress("Step_FinalPosZ",&Step_FinalPosZ);
   
   chain.SetBranchAddress("Step_InitialPosX",&Step_InitialPosX);
   chain.SetBranchAddress("Step_InitialPosY",&Step_InitialPosY);
   chain.SetBranchAddress("Step_InitialPosZ",&Step_InitialPosZ);
   chain.SetBranchAddress("Step_CreatorProcess",&Step_CreatorProcess);
   
   chain.SetBranchAddress("Step_FinalLogicalVolume",&Step_FinalLogicalVolume);
   chain.SetBranchAddress("Step_FinalProcess",&Step_FinalProcess);
   chain.SetBranchAddress("Step_InitialLogicalVolume",&Step_InitialLogicalVolume);
   chain.SetBranchAddress("Step_InitialMaterial",&Step_InitialMaterial);
   
   chain.SetBranchAddress("Step_FinalMaterial",&Step_FinalMaterial); 
   chain.SetBranchAddress("Step_InitialProcess",&Step_InitialProcess);
   chain.SetBranchAddress("Step_Particle",&Step_Particle);
   chain.SetBranchAddress("Track_EventID",&Track_EventID);
   
   chain.SetBranchAddress("Step_TrackID",&Step_TrackID);
   chain.SetBranchAddress("Track_TrackID",&Track_TrackID);
   chain.SetBranchAddress("Track_AccumulatedDose",&Track_AccumulatedDose);
   chain.SetBranchAddress("Track_AccumulatedEnergyDeposited",&Track_AccumulatedEnergyDeposited);   
   chain.SetBranchAddress("Track_AccumulatedEnergyLost",&Track_AccumulatedEnergyLost);
   
   chain.SetBranchAddress("Track_FinalKineticEnergy",&Track_FinalKineticEnergy);   
   chain.SetBranchAddress("Track_FinalLogicalVolume",&Track_FinalLogicalVolume);
   chain.SetBranchAddress("Track_InitialMaterial",&Track_InitialMaterial);
   chain.SetBranchAddress("Track_FinalMaterial",&Track_InitialMaterial);
   chain.SetBranchAddress("Track_Particle",&Track_Particle);   
   
// here begins the analysis

//** histogram details
   int nbins = 100;
   gStyle->SetOptStat(111111);
   gStyle->SetPalette(107);
   TCanvas * c1 = new TCanvas("c1", "c1", 800, 600);

   TH2F *pos1 = new TH2F("pos1", "pos1", nbins, -0.60, 0.60, nbins, -0.60, 0.60);

   //TH3F *pos1 = new TH3F("pos1", "pos1", nbins, 0., 0.0455, nbins, -0.30, 0.30, nbins, -0.30, 0.30);
   c1->Divide(1,1);     

//***run for fill the histogram   dx = Step_FinalPosX->at(count)-Step_InitialPosX->at(count); dy = Step_FinalPosY->at(count)-Step_InitialPosY->at(count); dz = Step_FinalPosY->at(count)-Step_InitialPosY->at(count);  cout<<"dx = "<<Step_FinalPosX->at(count)<<" dy = "<<Step_FinalPosY->at(count)<<" dz = "<<Step_FinalPosZ->at(count)<<endl;

   c1->cd(1);
   Long64_t nb = 0, nbytes = 0;
   nentries=chain.GetEntries();
   for(Int_t j=0; j<nentries;j++){ 
		
    	nb = chain.GetEntry(j); nbytes += nb; //run all the trees, branches and leaves  Step_FinalMaterial->at(count)!="control"&&

	cout<<"entries vector = "<<nentries<<endl;
		
	cout<<"Entry = "<<j<<endl;//*/
	
	if(Step_FinalMaterial->size()>0){
	Int_t count = 0;
	while(count<Step_FinalMaterial->size()-1){count++;} //eu preciso do contador para garantir que não pegará aqueles eventos que size=0. && &&Step_FinalMaterial->at(count)!="G4_WATER"	      dx = Step_FinalPosX->at(count)-Step_InitialPosX->at(count); dy = Step_FinalPosY->at(count)-Step_InitialPosY->at(count); dz = Step_FinalPosY->at(count)-Step_InitialPosY->at(count);

	
	     if((0.40<=Step_FinalPosX->at(count)<22.0)&&(0.40<=Step_FinalPosY->at(count)<22.0)&&(0.40<=Step_FinalPosZ->at(count)<22.0)){mm1++;  
	      if(Step_FinalMaterial->at(count)=="G4_WATER"){mm++;
	          cout<<"posX = "<<Step_FinalPosX->at(count)<<endl;
	          
	            for(Int_t l=0;l<Step_FinalMaterial->size();l++){
	                  //dx = dx +  Step_FinalPosX->at(l); 
	               	cout<<"posX = "<<Step_FinalPosX->at(l)<<endl;
	                  edep = edep + Step_AccumulatedEnergyDeposited->at(l);
	               }      
	                  sumDep=edep;
	                  edep=0;
	            }
	         
	                 pos1->Fill(Step_FinalPosX->at(count), Step_FinalPosY->at(count)); 
	           
	      }
	            //pos1->Fill(Step_FinalPosX->at(count), Step_FinalPosY->at(count));
	}
	
    }
cout<<endl;
cout<<"n of Ti = " << mm << endl;

// histogram styles PosX [mm];
pos1 -> SetNameTitle("energy", "PosY versus PosX;  PosX [mm];  PosY [mm]");
pos1 -> SetFillColor(0); 
pos1 -> Draw("COLZ");

// out file
TFile *p = new TFile("posXTi.root", "RECREATE"); // Sam 0.05mm

pos1->Write();

}

testeS_tree_1000.root (2.7 MB)

Wile_E_Coyote · May 12, 2020, 4:32pm

This has nothing to do with ROOT. Google -> C++ logical operators

{
  double x = 0.; // try e.g. 0. and then 10.
  if (0.40 <= x < 22.0) std::cout << "Surpriced?" << std::endl;
  if ((0.40 <= x) && (x < 22.0)) std::cout << "So what!" << std::endl;
}

David_Fernandes · May 12, 2020, 5:45pm

Thanks @Wile_E_Coyote!

Sorry for the doubt!

I just ask because only the second choice that you show works:

  if ((0.40 <= x) && (x < 22.0))

The first doesn’t and I don’t unsderstand why…

Thanks again…

Wile_E_Coyote · May 12, 2020, 6:15pm

Google -> learn C++

eguiraud · May 14, 2020, 11:50am

Hi David,
note that due to operator precedence 0.40 <= x < 22.0 is equivalent to (0.40 <= x) < 22.0, i.e. true/false < 22.0, i.e. always true (true and false correspond to 1.0 and 0.0 when cast to double).

Cheers,
Enrico

David_Fernandes · May 14, 2020, 3:36pm

Thanks @eguiraud this help very much!

God bless you!

system · May 28, 2020, 3:44pm

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.