#*********************************************************************************************** # File name: AnaTEST.py # Creator: T.T. Boehlen # Date: 10.03.2011 # Last Mod.: 10.03.2011, TTB #*********************************************************************************************** ### is meant to test read-out of the root file. # import ####################################################################################### import os,sys import math from Numeric import * from array import array import numpy as nP #import pylab as lP import scipy as sP #import ROOT from ROOT import TCanvas, TFile, TProfile, TH1F, TH2F, TF1,TGraphErrors,TLatex from ROOT import gROOT,gStyle,gSystem, Double,gPad,TGraph,TNtuple from ROOT import TPad, TPaveLabel, TPaveText, TLegend, TStyle,TBox #gBenchmark, gRandom # constants #################################################################################### #mass nucleon=(proton+neutron)/2 times c2 in MeV gMass_n=(938.272+939.566)/2.0 #speed of light m/s gCCC=299792458.0 # global settings ############################################################################## # enable running from another directory than the one where py-file resides workdir = os.path.dirname( sys.argv[0] ) if workdir: os.chdir( workdir ) # running mode ################################################################################# gInpPathInput='./' gInpFileNameEnd="_Out" # gInpRunNb="" gInpFileNameEnd="" gInpFileNameStart="MC_Evt1k_test" # # read data#### ################################################################################ print 'READ ROOT FILE ',gInpPathInput+gInpFileNameStart+gInpRunNb+gInpFileNameEnd+'.root',' ... ', fSimFile = TFile( gInpPathInput+gInpFileNameStart+gInpRunNb+gInpFileNameEnd+'.root' ) print 'DONE' #fSimFile.ls() fSimTree=fSimFile.Get("EventTree") fEventNb = fSimTree.GetEntries() #particle_common tr_n=fSimTree.FindLeaf("trn") tr_id=fSimTree.FindLeaf("trid") tr_paid=fSimTree.FindLeaf("trpaid") tr_type=fSimTree.FindLeaf("trtype") tr_cha=fSimTree.FindLeaf("trcha") tr_bar=fSimTree.FindLeaf("trbar") tr_mass=fSimTree.FindLeaf("trmass") tr_reg=fSimTree.FindLeaf("trreg") tr_ipx=fSimTree.FindLeaf("tripx") tr_ipy=fSimTree.FindLeaf("tripy") tr_ipz=fSimTree.FindLeaf("tripz") tr_ix=fSimTree.FindLeaf("trix") tr_iy=fSimTree.FindLeaf("triy") tr_iz=fSimTree.FindLeaf("triz") tr_fx=fSimTree.FindLeaf("trfx") tr_fy=fSimTree.FindLeaf("trfy") tr_fz=fSimTree.FindLeaf("trfz") tr_dead=fSimTree.FindLeaf("trdead") tr_gen=fSimTree.FindLeaf("trgen") ke_n=fSimTree.FindLeaf("ken") ke_id=fSimTree.FindLeaf("keid") ke_reg=fSimTree.FindLeaf("kereg") ke_inpx=fSimTree.FindLeaf("keinpx") ke_inpy=fSimTree.FindLeaf("keinpy") ke_inpz=fSimTree.FindLeaf("keinpz") # prepare data ################################################################################# mycount=0 #################################################################################################################### #read entries from ROOT tree in lists for iev in range(0,fEventNb): fSimTree.GetEntry(iev) cur_ke_n=ke_n.GetValue()#loop over KENTROS # print '### Event ',iev,' #############' #init for event for cur_i in range(0,int(cur_ke_n)): cur_ke_id=int(ke_id.GetValue(cur_i))-1 cur_ke_reg=ke_reg.GetValue(cur_i) cur_ke_inpx=ke_inpx.GetValue(cur_i) cur_ke_inpy=ke_inpy.GetValue(cur_i) cur_ke_inpz=ke_inpz.GetValue(cur_i) cur_tr_type=tr_type.GetValue(cur_ke_id) cur_tr_mass=tr_mass.GetValue(cur_ke_id) cur_tr_cha=tr_cha.GetValue(cur_ke_id) cur_tr_bar=tr_bar.GetValue(cur_ke_id) cur_tr_ipx=tr_ipx.GetValue(cur_ke_id) cur_tr_ipy=tr_ipy.GetValue(cur_ke_id) cur_tr_ipz=tr_ipz.GetValue(cur_ke_id) #calculate cur_tr_p =sqrt(cur_tr_ipx**2+cur_tr_ipy**2+cur_tr_ipz**2) # if curPART_P==0.0: curPART_P=1.0 cur_tr_T =(sqrt(cur_tr_p**2+cur_tr_mass**2)-cur_tr_mass)*1000.#convert to MeV if cur_tr_bar!=0: cur_tr_Tn =cur_tr_T/cur_tr_bar else: cur_tr_Tn =cur_tr_T cur_ke_p =sqrt(cur_ke_inpx**2+cur_ke_inpy**2+cur_ke_inpz**2) # if curPART_P==0.0: curPART_P=1.0 cur_ke_T =(sqrt(cur_ke_p**2+cur_tr_mass**2)-cur_tr_mass)*1000.#convert to MeV if cur_tr_bar!=0: cur_ke_Tn =cur_ke_T/cur_tr_bar else: cur_ke_Tn =cur_ke_T #look only for ions if (cur_tr_type==1)or(cur_tr_type<-1): # if cur_tr_T ok!) # if cur_tr_Tn>cur_ke_Tn+10:#look for particles which lost more than 10MeV/n if abs(cur_tr_Tn-cur_ke_Tn)<.01:#look for particles which have the same energy # if (iev==970): mycount=mycount+1 print 'Part:',cur_i,' Ty:',cur_tr_type,' Q:',cur_tr_cha,' A:',cur_tr_bar,' keID:',cur_ke_id,' Reg:',cur_ke_reg,' P_tr:',cur_tr_p,' P_ke',cur_ke_p,' Tn_tr:',cur_tr_Tn,' Tn_ke',cur_ke_Tn,' How many time?:',mycount # # curPART_AngX=math.degrees(arccos(curPART_PX/curPART_P)) # curPART_AngY=math.degrees(arccos(curPART_PY/curPART_P)) # curPART_AngZ=math.degrees(arccos(curPART_PZ/curPART_P)) # wait ######################################################################################### ## wait for input to keep the GUI (which lives on a ROOT event dispatcher) alive if __name__ == '__main__': rep = '' while not rep in [ 'q', 'Q' ]: rep = raw_input( 'enter "q" to quit: ' ) if 1 < len(rep): rep = rep[0] ################################################################################################ # EOF ########################################################################################## ################################################################################################