Analyse chain root 6.18.00


Please read tips for efficient and successful posting and posting code

ROOT Version: 6.18.00
Platform: Not Provided
Compiler: UBUNTU 18.04.04


Good morning!

I have a doubt: I’m using the code to analyse the chain of the tree generated in GAMOS6.1.0 but when I use the code the energy lost, that is supposed give equal or less than Event_InitialKineticEnergy, give more than Event_InitialKineticEnergy. I think that maybe is something in my code using wrong. Can you help with this? Bellow I place the code and in attached the root file used.

#include "TROOT.h"
#include "TMath.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>

void EnEmit(){  /// fecha linha 175

// declarations

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;
Int_t mm9=0;

Double_t eLost = 0;
Double_t eLos = 0;
Double_t edep = 0;
Double_t ed = 0;
Double_t ekin = 0;
Double_t ek = 0;

// read tree

        TChain chain("silverSPH");
        
        chain.Add("silverSPH_tree_1000.root");
        
// declarations of tree variables 

        Int_t nent=0, EventID=0;
	Double_t EventAEnerDep=0.0, EventInKinEn=0.0;

	std::vector<int> *StepNumber=0;
	std::vector<double> *StepFinKinEn=0;
	std::vector<double> *StepInKinEn=0;
	std::vector<double>  *StepAcEnDep=0;
   	std::vector<double>  *StepAcEnLost=0;

   	std::vector<double>  *StepFinPosX=0;
   	std::vector<double>  *StepFinPosY=0;
   	std::vector<double>  *StepFinPosZ=0;

	std::vector<string> *StepFinLogVol=0;//=0;
	std::vector<string> *StepParticle=0;
	std::vector<string> *StepInMaterial=0;//=0;
	
// read variables of the tree

//numbers
  chain.SetBranchAddress("Event_EventID",&EventID);
  chain.SetBranchAddress("Event_AccumulatedEnergyDeposited",&EventAEnerDep);
    chain.SetBranchAddress("Event_InitialKineticEnergy",&EventInKinEn);

// vector<int>
    chain.SetBranchAddress("Step_StepNumber",&StepNumber);
    
// vector<double>
    chain.SetBranchAddress("Step_FinalKineticEnergy",&StepFinKinEn);
    chain.SetBranchAddress("Step_InitialKineticEnergy",&StepInKinEn);
    chain.SetBranchAddress("Step_AccumulatedEnergyDeposited",&StepAcEnDep);
    chain.SetBranchAddress("Step_AccumulatedEnergyLost",&StepAcEnLost);
    chain.SetBranchAddress("Step_FinalPosX",&StepFinPosX);
    chain.SetBranchAddress("Step_FinalPosY",&StepFinPosY);
    chain.SetBranchAddress("Step_FinalPosZ",&StepFinPosZ);

// vector<string>
    chain.SetBranchAddress("Step_FinalLogicalVolume",&StepFinLogVol);
    chain.SetBranchAddress("Step_Particle",&StepParticle);
    chain.SetBranchAddress("Step_InitialMaterial",&StepInMaterial);
    
// fill histogram and load tree
  int nbins = 5000;
  gStyle->SetOptStat(1111111);
  TCanvas * c1 = new TCanvas("c1", "c1", 800, 600);
  //TH2F *pos = new TH2F("pos", "pos", nbins, 0., 0.0455, nbins, -40.0, 40.0); 
  TH2F *pos1 = new TH2F("pos", "pos", nbins, 0., 0.0455, nbins, -15.0, 15.0); 
  
  c1->Divide(1,3);

  c1->cd(1);
  
	Long64_t nb = 0, nbytes = 0;
	nent=chain.GetEntries();
	
     for(Int_t i = 0; i < nent; i++){
        
        nb = chain.GetEntry(i); nb += nbytes;
        
        cout<<"size vector = " << nent << endl;
        cout<<"entry = " << i << endl;   
        cout<<endl;
                        
        if(StepFinLogVol->size() != 0){
            Int_t cont = 0, cont1 = 0;
            
            while(cont < StepFinLogVol->size()-1 && StepFinLogVol->at(cont) != "control"){cont++;}
            
                if(StepFinLogVol->at(cont) == "control"){ mm4++;
                    for(Int_t j = 0; j < StepFinLogVol->size(); j++){
                        eLost = eLost + StepAcEnLost->at(j);
                        edep = edep + StepAcEnDep->at(j);
                        ekin = ekin + StepFinKinEn->at(j);
                    } eLos = eLost;
                      eLost = 0;
                      
                      ed = edep;
                      edep = 0;
                
                      ek = ekin;
                      ekin = 0;
                      
                      pos1->Fill(StepFinKinEn->at(cont), StepFinPosX->at(cont1));
                
                }
             	cout<<endl;
             	
             	if(StepFinLogVol->at(cont) == "silver"){mm++;}
             	if(StepFinLogVol->at(cont) == "I-125_SOURCE"){mm1++;}
             	if(StepFinLogVol->at(cont) == "air"){mm2++;}
             	if(StepFinLogVol->at(cont1) == "seed"){mm3++;}
             	if(StepFinLogVol->at(cont1) == "world"){mm5++;}
            	if(StepFinLogVol->at(cont1) == "OutOfWorld"){mm6++;}
             	
        
        }  
     
     }

cout<< "entry in silver =  "<< mm<<endl;
cout<< "entry in iodine =  "<< mm1<<endl;
cout<< "entry in air =  "<< mm2<<endl;
cout<< "entry in seed =  "<< mm3<<endl;
cout<<"###"<<endl;
cout<< "entry in controlo =  "<< mm4<<endl;
cout<< "entry in world =  "<< mm5<<endl; //contador contador
cout<< "entry in OutOfWorld =  "<< mm6<<endl;


cout<< "kin energy in control=  "<< ek<<endl;
cout<< "edep in control =  "<< ed<<endl;
cout<< "Elost in control =  "<< eLos<<endl;

pos1 -> SetNameTitle("energy", "PositionX versus EKin 0.05mm Ti;EKin;posX"); //Kinetic Energy Energy Lost
pos1 -> SetFillColor(0); 
pos1 -> Draw("COLZ");
c1->cd(2);
pos1->ProjectionX("PROJECTION IN X")->Draw();
c1->cd(3);
pos1->ProjectionY("PROJECTION IN Y")->Draw();

TFile *p = new TFile("EnEmitSil1.root", "RECREATE"); // Sam 0.05mm
  
        
}

Thanks

silverSPH_tree_1000.root (1.2 MB)

Running your macro I get this:

entry in silver =  4827
entry in iodine =  0
entry in air =  0
entry in seed =  0
###
entry in controlo =  173
entry in world =  0
entry in OutOfWorld =  0
kin energy in control=  0.119667
edep in control =  0.0122147
Elost in control =  0.0482441

What is wrong ?

This is exactly what I get too. The wrong is that the value of Elost in control will not assumed values above of 0.03549MeV, given that this is the energy that the particle was generated. Actually the value kin energy in control= 0.119667 was wrong too for the same reason.

Thanks

Well, the eLos evaluation is done by this code in your macro:

                if(StepFinLogVol->at(cont) == "control"){ mm4++;
                    for(Int_t j = 0; j < StepFinLogVol->size(); j++){
                        eLost = eLost + StepAcEnLost->at(j);
                        edep = edep + StepAcEnDep->at(j);
                        ekin = ekin + StepFinKinEn->at(j);
                    } eLos = eLost;
                      eLost = 0;

Programatically there is nothing wrong. I guess you should check the logic of that code…

It is also possible that you misinterpret the quantities that your tree provides. For example, simulations may report only that part of particle’s “energy loss” which is not “used” by its explicitly generated secondary particles (and so, in order to get the “total” energy loss, you would need to sum up the mother and its daughters).
I guess you would need to talk to your colleagues.

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