# Emittance calculation

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

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.

1 Like

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

Hi @faca87,

Cheers,
J.

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;
float xp;
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-xm)*(xp-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-xm)*(xp-xpm);
}
xxpm=xxpmN/nentries;
``````

I want to calculate <xx’>, but I don’t know how to say to root that x must be the entries of `PositionDirection.x` and xp 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

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

1 Like

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 and xp.
Can you help me to fix these lines please?

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;
tree.GetEntry(i); // this loads the contents of the branch into the associated variable `x`
``````

Cheers,
Enrico

Thank you @eguiraud
So, I wrote:

``````float x;
float xp;
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++) {
ts.GetEntry(jentry);
xxpmN+= (x-xm)*(xp-xpm);
}
xxpm=xxpmN/nentries;
``````

but I get errors due to these lines

Here the updated macro

scorelemma.cpp (4.3 KB)

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

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);
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)

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);
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;
``````

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

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

Alright, all good then!

1 Like

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