simulatanious fit to histograms (RooFit)

From: Arthur E. Snyder <snyder_at_slac.stanford.edu>
Date: Mon, 14 Jul 2008 10:36:03 -0700


Hi Root users,

I'm trying to do a simulatanious fit to three (stat independent histograms) using RooFit. I can't even figure out how to get past the 1st step of setting up an appropriate data set to fit with three histograms of three categories. RooDataHist does not seem to have a constructor that in an obvious way alows me to work from three histograms.

I have examples of unbinned fits, but nothing for binned case. Does anybody perhaps have a simple example of such a fit? It can't be hard. I could write it myself in plain root, but would rather not as I already have my 1d PDF set up in RooFit.

-Art S.

A.E. Snyder, Group EC                        \!c*p?/
SLAC Mail Stop #95                          ((.   .))
Box 4349                                        |
Stanford, Ca, USA, 94309                      '\|/`
e-mail:snyder_at_slac.stanford.edu                 o
phone:650-926-2701                              _
http://www.slac.stanford.edu/~snyder          BaBar
FAX:650-926-2657                          Collaboration



On Fri, 15 Jun 2007, Rene Brun wrote:

> Hello Art,
>
> I do not understand your request. MakeSelector/MakeClass are already
> generating code
> where all data member dimensions are of the type kMaxnel, kMaxnmu, etc. For
> example
> with the file Event.root generated by $ROOTSYS/test/Event, doing
> T->MakeSelector("Art");
> will generate Art.C and Art.h. Art.h having
> const Int_t kMaxfTracks = 612;
> ..
> Float_t fTracks_fPx[kMaxfTracks];
>
> Rene Brun
>
> Arthur E. Snyder wrote:
>>
>> It would be nice if MakeSelector would generate array dimensions that were
>> not hardwired, so that it can more easily be modified to copy with longer
>> arrays ...
>>
>> Something like this:
>>
>> class Selector : public TSelector {
>> public :
>> TTree *fChain; //!pointer to the analyzed TTree or TChain
>>
>> //array sizes for particles
>> enum { MAX_nK=10,
>> MAX_nel=5,
>> MAX_nmu=5,
>> MAX_nch=50,
>> MAX_nne=50,
>> MAX_nks=10,
>> MAX_nmc=100 };
>>
>> // Declaration of leave types
>> Int_t runnumber;
>> Int_t iammc;
>> UChar_t gotbmatch;
>> Int_t lunB;
>> Float_t deltae;
>> Float_t mes;
>> Float_t pxrecoB;
>> Float_t pyrecoB;
>> Float_t pzrecoB;
>> Float_t erecoB;
>> Float_t intpur;
>> Int_t mode;
>> Float_t R2All;
>> Float_t FoxMom2;
>> Float_t W2;
>> Float_t cosTT;
>> Float_t pxTrueRecoB;
>> Float_t pyTrueRecoB;
>> Float_t pzTrueRecoB;
>> Float_t eTrueRecoB;
>> Int_t lundIdRecoB;
>> Int_t uidRecoB;
>> Float_t pxTrueRecoilB;
>> Float_t pyTrueRecoilB;
>> Float_t pzTrueRecoilB;
>> Float_t eTrueRecoilB;
>> Int_t lundIdRecoilB;
>> Int_t uidRecoilB;
>> Int_t nK;
>> Int_t qK[MAX_nK]; //[nK]
>> Int_t uidK[MAX_nK]; //[nK]
>> Float_t pxrecoK[MAX_nK]; //[nK]
>> Float_t pyrecoK[MAX_nK]; //[nK]
>> Float_t pzrecoK[MAX_nK]; //[nK]
>> Float_t erecoK[MAX_nK]; //[nK]
>> Float_t precoK[MAX_nK]; //[nK]
>> Float_t pxK[MAX_nK]; //[nK]
>> Float_t pyK[MAX_nK]; //[nK]
>> Float_t pzK[MAX_nK]; //[nK]
>> Float_t eK[MAX_nK]; //[nK]
>> Float_t pK[MAX_nK]; //[nK]
>> Float_t xmass[MAX_nK]; //[nK]
>> Float_t xmassfit[MAX_nK]; //[nK]
>> Float_t chisq[MAX_nK]; //[nK]
>> Float_t globchisq[MAX_nK]; //[nK]
>> Int_t ndof[MAX_nK]; //[nK]
>> Int_t nel;
>> Int_t el_q[MAX_nel]; //[nel]
>> Int_t el_uid[MAX_nel]; //[nel]
>> Float_t el_pxLab[MAX_nel]; //[nel]
>> Float_t el_pyLab[MAX_nel]; //[nel]
>> Float_t el_pzLab[MAX_nel]; //[nel]
>> Float_t el_pLab[MAX_nel]; //[nel]
>> Float_t el_pxStar[MAX_nel]; //[nel]
>> Float_t el_pyStar[MAX_nel]; //[nel]
>> Float_t el_pzStar[MAX_nel]; //[nel]
>> Float_t el_eStar[MAX_nel]; //[nel]
>> Float_t el_pStar[MAX_nel]; //[nel]
>> Int_t nmu;
>> Int_t mu_q[MAX_nmu]; //[nmu]
>> Int_t mu_uid[MAX_nmu]; //[nmu]
>> Float_t mu_pxLab[MAX_nmu]; //[nmu]
>> Float_t mu_pyLab[MAX_nmu]; //[nmu]
>> Float_t mu_pzLab[MAX_nmu]; //[nmu]
>> Float_t mu_pLab[MAX_nmu]; //[nmu]
>> Float_t mu_pxStar[MAX_nmu]; //[nmu]
>> Float_t mu_pyStar[MAX_nmu]; //[nmu]
>> Float_t mu_pzStar[MAX_nmu]; //[nmu]
>> Float_t mu_eStar[MAX_nmu]; //[nmu]
>> Float_t mu_pStar[MAX_nmu]; //[nmu]
>> Int_t nch;
>> Int_t ch_q[MAX_nch]; //[nch]
>> Int_t ch_uid[MAX_nch]; //[nch]
>>
>> That was dimension for, e.g., mu momenta could be changed just in the enum.
>>
>> This would be easy to implement in TreePlayer::MakeClass where all the
>> needed info has been accumulated.
>>
>> -Art S.
>>
>> A.E. Snyder, Group EC \!c*p?/
>> SLAC Mail Stop #95 ((. .))
>> Box 4349 |
>> Stanford, Ca, USA, 94309 '\|/`
>> e-mail:snyder_at_slac.stanford.edu o
>> phone:650-926-2701 _
>> http://www.slac.stanford.edu/~snyder BaBar
>> FAX:650-926-2657 Collaboration
>>
>>
>>
>> On Tue, 6 Feb 2007, Rene Brun wrote:
>>
>>> Arthur E. Snyder wrote:
>>> > One thing you want to be careful of is if you wish to keep the
>>> covariance
>>> > matrix. If you've done another fit keeping your TF1 won't keep
>>> > convariance. TF1 caches the errors and values but not the full
>>> convariance
>>> > matrix. This could be a problem if, e.g, you wanted to do a systematic
>>> > study by varying fit parameters according to correlated errors. It would
>>> > not be hard to had covariance cache to TF1 and perhaps it will be done
>>> > someday. In the meanwhile if you want to keep covariance you have to
>>> hang
>>> > on to it for yourself.
>>> > Art, Roger,
>>>
>>> See the small script cov.C below showing how to save the covariance matrix
>>> with the fitted histogram.
>>>
>>> Rene Brun
>>>
>>> //macro cov.C
>>> void cov() {
>>> // shows how to save the Fit covariance matrix as an object
>>> // in the histogram list of functions
>>>
>>> covwrite();
>>> covread();
>>> }
>>>
>>> void covwrite() {
>>> // shows how to save the Fit covariance matrix as an object
>>> // in the histogram list of functions
>>>
>>> TH1F *h = new TH1F("h","sho cov matrix",100,-3,3);
>>> h->FillRandom("gaus",5000);
>>> h->Fit("gaus");
>>> TVirtualFitter *fitter = TVirtualFitter::GetFitter();
>>> TMatrixD *matrix = new TMatrixD(3,3,fitter->GetCovarianceMatrix());
>>> h->GetListOfFunctions()->Add(matrix);
>>> h->GetListOfFunctions()->ls();
>>> TFile f("h.root","recreate");
>>> h->Write();
>>> }
>>>
>>> void covread() {
>>> TFile *f = TFile::Open("h.root");
>>> TH1F *h = (TH1F*)f->Get("h");
>>> h->Draw();
>>> TMatrixD *m =
>>> (TMatrixD*)h->GetListOfFunctions()->FindObject("TMatrixT<double>");
>>> m->Print();
>>> }
>>>
>>> > -Art
>>> >
>>> >
>>> > [a] I'm usually just a consumer of root help, but since this one is
>>> > something I knew about I thought I try to help out. A real root expert
>>> > might have more or better answers.
>>> >
>>> > A.E. Snyder, Group EC \!c*p?/
>>> > SLAC Mail Stop #95 ((. .))
>>> > Box 4349 |
>>> > Stanford, Ca, USA, 94309 '\|/`
>>> > e-mail:snyder_at_slac.stanford.edu o
>>> > phone:650-926-2701 _
>>> > http://www.slac.stanford.edu/~snyder BaBar
>>> > FAX:650-926-2657 Collaboration
>>> >
>>> >
>>> >
>>> > On Tue, 6 Feb 2007, Roger Mason wrote:
>>> >
>>> > >> Hello Art,
>>> >>
>>> >> "Arthur E. Snyder" <snyder_at_slac.stanford.edu> writes:
>>> >>
>>> >> >>> Well in c/c++ an array is just a pointer with a certain amount of
>>> space
>>> >>> reseved beyond it. So GetParErrors give you a pointer (i.e. the array
>>> held
>>> >>> by your TF1 object (fParErrors)).
>>> >>> >> I'm still at the stage with c++ where the pointer syntax, combined
>>> >> with object syntax, eludes me. Such was the case here.
>>> >>
>>> >> >>> If you want a copy that will not go away when your TF1 goes away
>>> (to be
>>> >>> safe) then I think you just have to copy it element by element, or use
>>> it
>>> >>> to construct a TMatrix or TVector or something.
>>> >>> >> Well, as long as I have the TF1 in the root file I am OK because I
>>> can
>>> >> always extract the parameters and errors and dump the into a text
>>> >> file. Why do I need this? Because I always tabulate results in
>>> >> drafts of my publications just so the referees can tell me to toss
>>> >> 'em and keep the pictures.
>>> >>
>>> >> >>> If your a fortran old-foggie like me this can be very confusing!
>>> >>>
>>> >>> >> Fortran was my first programming language too. I'm just glad that
>>> >> root is'nt written in Java.
>>> >>
>>> >> Thanks again,
>>> >> Roger
>>> >>
>>> >>
>>> >>
>>>
>
>
Received on Mon Jul 14 2008 - 19:36:12 CEST

This archive was generated by hypermail 2.2.0 : Tue Jul 15 2008 - 23:50:02 CEST