Smearing an energy according to a crystal ball function

Newbie
1

Hi there,
I have a TTree containing decay products, with energies given by pure theory (conservation of energy, momentum etc.), in the proportions predicted by pure theory (QFT decay width).

I would like to create distributions of what those particle energies would be reconstructed as after passing through an crystal calorimeter.

I have a crystal ball function:

double crystalball_function(double x, double alpha, double n, double sigma, double mean) {
  // evaluate the crystal ball function
  if (sigma < 0.)     return 0.;
  double z = (x - mean)/sigma; 
  if (alpha < 0) z = -z; 
  double abs_alpha = std::abs(alpha);
  // double C = n/abs_alpha * 1./(n-1.) * std::exp(-alpha*alpha/2.);
  // double D = std::sqrt(M_PI/2.)*(1.+ROOT::Math::erf(abs_alpha/std::sqrt(2.)));
  // double N = 1./(sigma*(C+D));
  if (z  > - abs_alpha)
    return std::exp(- 0.5 * z * z);
  else {
    //double A = std::pow(n/abs_alpha,n) * std::exp(-0.5*abs_alpha*abs_alpha);
    double nDivAlpha = n/abs_alpha;
    double AA =  std::exp(-0.5*abs_alpha*abs_alpha);
    double B = nDivAlpha -abs_alpha;
    double arg = nDivAlpha/(B-z);
    return AA * std::pow(arg,n);
  }
}

And I get entries from sampleTree1 to write to sampleTree2

  TChain *sampleTree1 = new TChain("sample","");
  sampleTree1->Add("filenames");
  Double_t x,y;
  sampleTree1->SetBranchAddress("x",&x);
  sampleTree1->SetBranchAddress("y",&y);

double alpha=1.19, double n=10, double sigma=3,49, double mean=70;
TTree *sampleTree2 = new TTree("sampleTree2","sampleTree2");
  Double_t truex=0,truey=0,genweight=0, reconx=0,recony=0;
  sampleTree2->Branch("truex",&truex);
  sampleTree2->Branch("truey",&truey);
  sampleTree2->Branch("reconx", &reconx);
  sampleTree2->Branch("recony", &recony);
  double smear1=1,smear2=1,smear3=1,smear4=1;
  for (Long64_t jentry=0; jentry<nentries; jentry++) {
    nb = sampleTree1->GetEntry(jentry);   nbytes += nb;
    truex=x;
    truey=y;
    reconx=crystalBall(x,alpha, n, sigma, mean);
    recony=crystalBall(y,alpha, n, sigma, mean);
    sampledTree2->Fill();
  }

I am confused because when I do this I get reconstructed values that are exactly proportional to the true values–just a straight line plot. Shouldn’t the crystal ball function smear the value? Or do I need to insert a smearing manually in some other way (how?).

Many thanks!

EDIT: my basic question is equivalent to: how would I smear a ttree full of mono-energetic particles so that when I plot the resulting smeared energy, it looks like a crystal ball function?

2

Perhaps @moneta or @StephanH can give a hand here?

3

Hi,

For smearing x you need to generate random values from the Cristal ball with mean = xtrue.
You would need to do :

//outside the TTree loop 
TF1 f("f","crystalball",-10*sigma,10*sigma);
// set the c.b.  parameters (mean=0) and constant=1 
f.SetParameters(1,0,sigma,alpha,n);  

// inside the event loop 

       reconx = truex + f.GetRandom();

Lorenzo

4

Many thanks!

I have a more physics related question, if this is an appropriate place. (It is also more subtle question I think.)

I have two types of particles in my detector.

Particle A is mono-energetic, with energy ~139.57018MeV (if such a particle exists… :-> )

Particle B has an energy distribution that depends on the parameters of theory, and I wish to extract those parameters. So I cannot fit the crystal ball function directly to the more interesting Particle B histogram because that observable depends on the parameters I wish to extract.

My desired process is:
–Fit a crystal ball function to the measured energy histogram of the mono-energetic Particle A and get the fit parameters.

–Then use those (adapted?) fit parameters to smear the true energies of Particle B, to create a simulated distribution. Then compare simulated histograms to measured and extract parameters.

In other words, I want to get the detector response from a particle where the energy does not depend on my parameters of interest. Then apply it to my simulated particles, holding that response constant while I extract parameters.

For the record, I know that both Particles A and B will have response functions very close to a crystal ball. And on physical grounds the response should be very similar. I just cannot fit the crystal ball function directly to the distribution of Particle B, because aspects of the fit (like tail fraction) would change with my theory parameters, confounding the extraction of parameters.

Perhaps my hope is impossible. But we must dream.

Thank you so much for any help you can offer!

EDIT: I realize the actual question is not specified. I would like to know if my procedure is possible: to use the response function for a known particle and apply to another similar particle. Possible answers could be:
–No your dream is impossible. The crystal ball fit parameters are too closely related to the specifics of the particle you fit it against.
–Yes your dream can be done, all you have to do is get the fit parameters from particle A, set the crystal ball function to these parameters, change the constant=1 to constant=someValue, and apply the smearing as moneta indicates.
–Yes your dream can be done, all you have to do is get the fit parameters from particle A, and perform a rain dance.

5

It sounds like what you have in mind is “unfolding”. You could search for “unfolding high energy physics”, and check if that’s what you are after.

Unfolding means to take an observed distribution, and reconstruct the “true” distribution (what you do with particle A). Now, you can either apply the inverse transformation to generate the “observable” distribution of B, or you directly compare B unfolded with B truth.

6

Sort of. I am not unfolding with particle A. I am simply finding the fit parameters. It is possible that once I have the fit parameters from particle A, I could unfold the measured distribution of particle B, rather than folding the true values and comparing to measurement.

The question is: whether I can use the parameters from a fit on particle A on the folding (or unfolding) of particle B. And if so, how this works with crystal ball function. I think, e.g. that I keep the mean=0, but e.g. should I change parameters related to tail fraction before applying to particle B?

Many thanks!

EDIT: I think I see what you mean however. I could use the crystal ball function to fold (i.e. convolve) the crystal ball distribution with the true distribution, rather than changing each each event energy in the tree. But I am still a little unclear about whether I can use the crystal ball function from a mono-energetic particle applied to a particle with varying energy…

EDIT 2: I need to do it in the tree. The variable that I eventually plot is not the energy itself, but a transformation of the energy. The best method is still to change the energy in the tree rather than fold or unfold the histograms after filling. My concern is that the crystal ball will subtract energy that is “lost” by the detector, but the proportion of energy lost in this way–related to the “tail fraction”–is dependent on the energy of the incident particle. For a small energy particle defining reconx by

reconx = truex + f.GetRandom();

could wind up being a negative number, because f has parameters from a particle with large energy.

7

I guess what you could do is

  1. Fit the crystal ball parameters with particle A.
  2. Use the parameters as starting values for the “smearing” of B.
  3. Set the energy-dependent parameter (I would say that this is the mean) to whatever particle B (truth) has.
  4. Throw a random number to generate a point on the crystal ball distribution.
  5. Save this as the smeared energy.
  6. Start at step 2. with the next particle B.

You might need some sanitisation of outputs to e.g. cut off negative energies etc.
You might also need to check which parameters actually depend on the energy. Maybe it’s not only the mean but also other parameters. If you can vary the energy of particle A, you can directly measure this by fitting to distributions generated from multiple energies. If multiple parameters change with the energy, you might have to adjust them simultaneously.