Emittance calculation

[faca87]

I simulated by geant muon production in e+e- annihilation e+e- ->mu+ mu-
By the geant simulation, I save in a ROOT file the muon directions x,y,z and I calulcated the angles of produced muons with x and y axes as:

x'=arcos (px/P) and y'=arcos(py/P)

where px and py are the components of the momentum along x and y axes, P is the total momentum.

Then by root I make x-x’ and y-y’ plots, for example, by writing

  TH2F *hemitxxpsp = new TH2F("hemitxxpsp", "", 500, 0., 0.,5000, 0.,0.);
         TString varhemitxxpsp = TString::Format("SecondaryParticleAng.x : PositionDirection.x >> hemitxxpsp");
         ts->Draw(varhemitxxpsp);

I get this one

Now I’ve to calculate the emittance epsilon of produced muons, that is

so I need

1. <x^2>
2. <x’^2>
3. <xx’>

How can I get these values? For example printing them in the stat box?

Thank you

[dastudillo]

In the Stats box, Std Dev x and Std Dev y are the square root of 1 and 2, resp.
You can calculate 3 yourself by looping over the data (tree, I suppose) and using the corresponding formula; and since you are doing that, I suggest you also calculate 1 and 2, to compare against the stats box.

[faca87]

Hi @dastudillo, thank you for your reply.
This is the macro that I wrote to plot data

simlemma.cpp (129.7 KB)

and here the ROOT file https://we.tl/t-ueQZfFmmp7

please, can you help me to modify it to calculate <x^2>, 1. <x’^2> and <xx’> ? Do you have an example of the code?

Thank you

[jalopezg]

Hi @faca87,

I am inviting @couet to this thread, as he may know the answer to your question.

Cheers,
J.

[faca87]

Thank you @jalopezg to ping @couet
Anyway, waiting for @couet or someone else (for example, probably @Wile_E_Coyote @bellenot or @moneta may know the answer), I was trying to start to write the code and this is the updated macro

simlemma.cpp (130.3 KB)

As I wrote in my first message, I make the plots in this way

TH2F *hemitxxpsp = new TH2F("hemitxxpsp", "", 500, 0., 0.,5000, 0.,0.);
         TString varhemitxxpsp = TString::Format("SecondaryParticleAng.x : PositionDirection.x >> hemitxxpsp");
         ts->Draw(varhemitxxpsp);

i.e. my TTree in the Root file is ts and my branches are SecondaryParticleAng.x and PositionDirection.x (for the x-x’ plots).

So, to calculate <xx’> (look the formula in my first message), I started to write this

	TBranch        *b_PositionDirection.x;   //!
		TBranch        *b_SecondaryParticleAng.x;   //!
		float x[3];	
		float xp[3];
		float xm;
		float xpm;
		float xxpmN;
		float xxpm;
		Long64_t nentries = ts->GetEntriesFast();
		xm= hemitxxpsp->GetXaxis()->GetRMS();
		xpm=hemitxxpsp->GetYaxis()->GetRMS();
	
	for(Long64_t jentry=0; jentry<nentries;jentry++) {
			xxpmN+= (x[0]-xm)*(xp[0]-xpm);
		}
		xxpm=xxpmN/nentries;

i.e.

1. Long64_t nentries = ts->GetEntriesFast();
   I want to read how many entries I’ve in ts

		xm= hemitxxpsp->GetXaxis()->GetRMS();
			xpm=hemitxxpsp->GetYaxis()->GetRMS();

I want to calculate the means < x > and <x’>

for(Long64_t jentry=0; jentry<nentries;jentry++) {
				xxpmN+= (x[0]-xm)*(xp[0]-xpm);
			}
			xxpm=xxpmN/nentries;

I want to calculate <xx’>, but I don’t know how to say to root that x[0] must be the entries of PositionDirection.x and xp[0] must be the entries of SecondaryParticleAng.x

Please, do someone know how to finish my code (and eventually fix it if there are some errors)?
Thank you

[eguiraud]

Hi @faca87,
the code you shared is almost 3000 lines, that’s way too much for us to check it and debug it.

Given you have data in a ROOT file, you can read it using TTree, TTreeReader or RDataFrame (see e.g. this post) and perform any calculation you want with the values. The link points e.g. to a TTreeReader tutorial that shows how to read all values in your TTree in a loop.

I hope that helps!
Enrico

[faca87]

Ho @eguiraud yes the macro is very long, but I don’t need your help for all the macro!
I just need your help for the piece of the macro that I wrote in my previous message!
You see they are few linea…in particular I need help to assign the values of SecondaryParticleAng.x and PositionDirection.x (branche of my root file) to the arrays x[0] and xp[0].
Can you help me to fix these lines please?

[eguiraud]

Hi,
yes I got that You need to read the TTree entries into your x and xp variables – that post that I linked above points to various method to do that, and related tutorials.

The short version is:

// extract the TTree from the file 
auto tree = file.Get<TTree>("treename");
// attach the TTree branch to your variable
float x[3];
tree.SetBranchAddress("x", x, "x[3]/I");
// read each entry
tree.GetEntry(i); // this loads the contents of the branch into the associated variable `x`

Cheers,
Enrico

[faca87]

Thank you @eguiraud
So, I wrote:

float x[3];	
			float xp[3];
			float xm;
			float xpm;
			float xxpmN;
			float xxpm;
			Long64_t nentries = ts->GetEntriesFast();
			xm= hemitxxpsp->GetXaxis()->GetRMS();
			xpm=hemitxxpsp->GetYaxis()->GetRMS();
			ts.SetBranchAdress("PositionDirection.x",PositionDirection.x,"x[3]/I");
			ts.SetBranchAdress("SecondaryParticleAng.x",SecondaryParticleAng.x,"xp[3]/I");
		
		for(Long64_t jentry=0; jentry<nentries;jentry++) {
			ts.GetEntry(jentry);
				xxpmN+= (x[3]-xm)*(xp[3]-xpm);
			}
			xxpm=xxpmN/nentries;

but I get errors due to these lines

please, can you check it?

Here the updated macro

scorelemma.cpp (4.3 KB)

[eguiraud]

Hi,
the compiler error says it all really: GetXaxis() and GetYaxis() return a TAxis object and TAxis does not have a GetRMS method. Maybe you meant to call hemitxxpsp->GetRMS(), i.e. TH1::GetRMS?

Similarly for other errors: the compiler is pretty explicit, you need ts-> instead of ts. because ts is a pointer, and PositionDirection.x is not a valid variable, you probably want just x if you want to associate branch "PositionDirection.x" with variable x.

Hope this helps!
Enrico

[faca87]

Thank you @eguiraud

1. hemitxxpsp is a TH2F in which on x axis I plot the position (i.e. the branch PositionDirection.x ) and on y axis I plot the angle (i.e. the branch SecondaryParticleAng.x ). In my formula I need the mean value both of PositionDirection.x and SecondaryParticleAng.x (i.e. < x > and < x’>). For erro I wrote GetRMS() but I need GetMean().
   Given that I need both the means of x and y axis, here Mean from a 2d histogram found that I must use GetMean(1) for the x axis and GetMean(2) for the y one! Did I understand right?

2. I wrote

float x;	
			float xp;
			float xm;
			float xpm;
			float xxpmN;
			float xxpm;
			int jentry;
			int nentries = ts->GetEntries();
			xm= hemitxxpsp->GetMean(1);
			xpm= hemitxxpsp->GetMean(2);
			ts->SetBranchAdress("PositionDirection.x",PositionDirection.x,"x/I");
			ts->SetBranchAdress("SecondaryParticleAng.x",SecondaryParticleAng.x,"xp/I");
		    for(jentry=0; jentry<nentries;jentry++) {
			ts->GetEntry(jentry);
				xxpmN+= (x-xm)*(xp-xpm);
			}
			xxpm=xxpmN/nentries;

but I get these errors

anyway, I wrote a very short macro (just 52 lines), please, can you check it? I think it’s easier for you to help me looking this very shor macro

simlemmared.cpp (1.6 KB)

[faca87]

I’m not sure, but maybe I solved myself in this way

	Double_t x;	
			Double_t xp;
			float xm;
			float xpm;
			float xxpmN;
			float xxpm;
			int jentry;
			int nentries = ts->GetEntries();
			xm= hemitxxpsp->GetMean(1);
			xpm= hemitxxpsp->GetMean(2);
			ts->SetBranchAddress("PositionDirection.x", &x);
			ts->SetBranchAddress("SecondaryParticleAng.x", &xp);
		    for(jentry=0; jentry<nentries;jentry++) {
			ts->GetEntry(jentry);
				xxpmN+= (x-xm)*(xp-xpm);
				std::cout << "entry = " << jentry << "\t x = " << x << " \t xp = " << xp <<  " \t xxpmN = " << xxpmN <<  std::endl;
			}
			xxpm=xxpmN/nentries;

[eguiraud]

That’s ok if PositionDirection.x is a single value, if it’s an array you need x to be an array too.

[faca87]

Hi @eguiraud PositionDirection.x has a value for each entry…i.e. 1 entry= 1 value

[eguiraud]

Alright, all good then!