search

ytatus94 / Higgsino

The NUHM2 Higgsino study
0 stars 0 forks source link

Skype discussion with Giuliano #33

Open ytatus94 opened 7 years ago

ytatus94 commented 7 years ago

[03.08.17, 12:44:57 PM] Yu-Ting: Hi Giuliano, do you use PyROOT?

[03.08.17, 12:45:37 PM] Giuliano Gustavino: Hi! yes sometimes I do

[03.08.17, 12:45:59 PM] Yu-Ting: ok, then I have a question about using PyROOT to do the fitting

[03.08.17, 12:46:30 PM] Giuliano Gustavino: which kind of fit? histgram fitting?

[03.08.17, 12:47:05 PM] Yu-Ting: yes, here is the tutorial part for using C++ https://root.cern.ch/root/html534/guides/users-guide/FittingHistograms.html#fitting-multiple-sub-ranges for these lines of example

// Get the parameters from the fit
   g1->GetParameters(&par[0]);
   g2->GetParameters(&par[3]);
   g3->GetParameters(&par[6]);

How can I do the same thing in PyROOT because I cannot use &par[x] in PyROOT

[03.08.17, 12:49:33 PM] Giuliano Gustavino: https://root-forum.cern.ch/t/how-can-i-do-a-fit-with-a-custom-function-in-pyroot/20361

[03.08.17, 12:50:07 PM] Yu-Ting: In my case, I want to fit the peak using gaus and to fit the tail using expo. So it is similar to this example. I can fit peak and tail separately but I don't know how to get the parameters and put it into the total function to fit

[03.08.17, 12:50:41 PM] Giuliano Gustavino: you use the GetParameters function and then you have a list with all the parameters and you can set/use whatever you want

[03.08.17, 12:51:25 PM] Yu-Ting: I know I have to use GetParameters() but I cannot use "&" to get the memory address in PyROOT That is my question

[03.08.17, 12:52:44 PM] Giuliano Gustavino: why do you need it?

[03.08.17, 12:52:49 PM] Yu-Ting: In the ROOT tutorial example, there is an array par, and it use GetParameters() to set the parameters and assign to par elements Because it uses address to access the par elements in C++, I don't know how shall I do in PyROOT

**[03.08.17, 12:53:49 PM] Giuliano Gustavino: look at the link that I sent you it uses GetParamters() to have a list of all the parameters

parameters = fit.GetParameters()

and then you can use them w/o using memory address**

[03.08.17, 12:55:10 PM] Yu-Ting: Oh~ I didn't notice there is a scroll bar so I didn't notice the part at the bottom. Thank you very much! I will try it

[03.08.17, 12:55:54 PM] Giuliano Gustavino: np ;)

ytatus94 commented 7 years ago

[04.08.17, 4:13:01 PM] Yu-Ting: Hi Giuliano, I have question about using PyROOT to do fitting again

[04.08.17, 4:13:12 PM] Giuliano Gustavino: wait meeting write me but I will reply as soon as I can ok?

[04.08.17, 4:16:44 PM] Yu-Ting: Our convenor asked me to use user defined function to do fit. And I implement it. But when I use this user defined function, the root complain 

Error in <Fit>: function fit_higgsino_core has illegal number of parameters = 0
Error in <Fit>: function fit_higgsino_core has illegal number of parameters = 0
Error in <Fit>: function fit_higgsino_core has illegal number of parameters = 0

But I did use fit_higgsino_core.SetParameters(1., n1, n2) in my script, why root still use 0 as initial values?

[04.08.17, 4:19:16 PM] Giuliano Gustavino: which is the fit_higgsino_core function?

[04.08.17, 4:19:41 PM] Yu-Ting: fit_higgsino_core is the histogram to be fit oh~ I see, I should use function to set parameters not histogram

[04.08.17, 4:20:36 PM] Giuliano Gustavino: yep

[04.08.17, 4:21:36 PM] Yu-Ting: sorry I didn't notice this

[04.08.17, 4:22:38 PM] Giuliano Gustavino: np ;)

[04.08.17, 4:23:27 PM] Yu-Ting: Wait! fit_higgsino_core is TF1 fit_higgsino_core = ROOT.TF1("fit_higgsino_core", funcMllDistr, 0, higgsino_dm)

def funcMllDistr(x, par):
    '''
    taken from https://arxiv.org/abs/0704.2515
    equation 4, page 9 
    NB there is a typo in the formula (just in writing it in the the paper, not in the work)
    the correct formula is in this code
    '''
    m1 = par[1]
    m2 = par[2]
    mu = m2 - m1
    m = x
    M = m1 + m2
    m_Z = 91
    norm = par[0] * m

    radice = math.sqrt( pow(m, 4) - pow(m, 2)  ( pow(mu, 2) + pow(M, 2) ) + pow(mu  M, 2) )
    normalizzazione = pow(  pow(m, 2) - pow(m_Z, 2), 2)
    var = (norm * radice) / normalizzazione
    var = var  (-2  pow(m, 4) + pow(m, 2)  (2  pow(M, 2) - pow(mu, 2)) + pow((mu * M), 2) )

    if x > math.fabs(m2) - math.fabs(m1):
        var = 1

    return var

[04.08.17, 4:26:16 PM] Giuliano Gustavino: (it is clearly written by an italian) don’t you thinkg that the best is starting from a simpler function? like something that is for example already included among the defaults pdfs in root? like pol0,1,….polN look from them if they have a not null list of parameters if not the problem is in the TF1 definition

[04.08.17, 4:28:24 PM] Yu-Ting: You mean use a simpler function to test and if it works then modify it to my function?

[04.08.17, 4:28:31 PM] Giuliano Gustavino: yep

[04.08.17, 4:29:14 PM] Yu-Ting: ok, i will try The problem is not related my code. I think it is PyROOT problem This is my test

[yushen@lxplus053 pythons]$ python
Python 2.7.13 (default, Mar  8 2017, 14:28:24) 
[GCC 4.9.3] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> def pol2(x, par):
...   return par[0] + par[1]  x + par[2]  x * x
... 
>>> import ROOT
>>> fit = ROOT.TF1("fit_1", pol2, 0, 20)
>>> print fit.GetParameter(0)
0.0
>>> print fit.GetParameter(1)
0.0
>>> print fit.GetParameter(2)
0.0
>>> fit.SetParameters(5, 6, 7)
>>> print fit.GetParameter(0)
0.0
>>> print fit.GetParameter(1)
0.0
>>> print fit.GetParameter(2)
0.0
>>>

After SetParameters(), and I try to GetParameter() but they are still 0.0

[04.08.17, 4:57:53 PM] Giuliano Gustavino: I would say to start from this

import ROOT
from ROOT import TF1, TH1F

x = ( 1.913521, 1.953769, 2.347435, 2.883654, 3.493567,
       4.047560, 4.337210, 4.364347, 4.563004, 5.054247,
       5.194183, 5.380521, 5.303213, 5.384578, 5.563983,
       5.728500, 5.685752, 5.080029, 4.251809, 3.372246,
       2.207432, 1.227541, 0.8597788,0.8220503,0.8046592,
       0.7684097,0.7469761,0.8019787,0.8362375,0.8744895,
       0.9143721,0.9462768,0.9285364,0.8954604,0.8410891,
       0.7853871,0.7100883,0.6938808,0.7363682,0.7032954,
       0.6029015,0.5600163,0.7477068,1.188785, 1.938228,
       2.602717, 3.472962, 4.465014, 5.177035 )

np = len(x)
h = TH1F( 'h', 'Example of several fits in subranges', np, 85, 134 )
h.SetMaximum( 7 )

for i in xrange(np):
    h.SetBinContent( i+1, x[i] )

g1 = TF1( 'g1', 'gaus',  85,  95 )

h.Fit( g1, 'R' )

parameters = g1.GetParameters()
print parameters[0], parameters[1],parameters[2]

[04.08.17, 4:57:57 PM] Giuliano Gustavino: it works

[04.08.17, 4:59:04 PM] Yu-Ting: gaus is default function and I know it works. Because I used gaus + expo to fit yesterday

[04.08.17, 4:59:23 PM] Giuliano Gustavino: so it is what i was suggesting start from a default function

[04.08.17, 4:59:35 PM] Yu-Ting: but our convenor asked me to change the fitting function.

[04.08.17, 4:59:48 PM] Giuliano Gustavino: so the problem is that you are not writing well the function

[04.08.17, 4:59:51 PM] Yu-Ting: and specify the function I have to use. No. I think the problem is not the function itself. I use interactive mode in python and define a pol2 and I set parameters for this function then I print the parameters. but all the parameters show 0 even I already set them

[04.08.17, 5:02:00 PM] Giuliano Gustavino: the same example of before but using `g1 = TF1("parabola","[0]+[1]x+[2]x2", 85, 95 )` works as well**

[04.08.17, 5:02:33 PM] Yu-Ting: then it is very strange... :(

ytatus94 commented 7 years ago

[14.08.17, 11:32:19 AM] Yu-Ting: Hi Giuliano, do you know how to scale the TRUTH3 MC sample to the correctly luminosity?

[14.08.17, 11:33:17 AM] Giuliano Gustavino: Hi what do you mean? why should truth3 be different to the reco MC samples?

[14.08.17, 11:34:27 AM] Yu-Ting: Because reco MC has event weight, lepton weight, jet weight, pileup weight information, but I don't know how to get these information in truth3 MC

[14.08.17, 11:34:45 AM] Giuliano Gustavino: ah well because you cannot have them you have to use only the event weight that is proper of the MC generator the other weights are evaluated by studying the object efficiencies so you cannot assign them to the truth event

[14.08.17, 11:36:44 AM] Yu-Ting: so I only can get the event weight from truth3 samples. Does the "weight = Xsec filter efficiency lump / (total event weight)" only in truth3? and where can I find the event weight?

[14.08.17, 11:52:27 AM] Giuliano Gustavino: yes exaclty do you have a framework which is able to run on the DAOD right? so you run it on the truth3 sample and you can get the event weight as usual

[14.08.17, 11:54:14 AM] Yu-Ting: I don't understand it. I have two framework, one for running reco and one for running truth3. I need to have two because the tree structure in reco and in truth levels are different.