source: JEM-EUSO/esaf_cc_at_lal/packages/reconstruction/modules/shower/energy/src/EnergyFitUtils.cc @ 114

Last change on this file since 114 was 114, checked in by moretto, 11 years ago

actual version of ESAF at CCin2p3
File size: 5.8 KB
Line 
1// $Id: EnergyFitUtils.cc 2975 2011-08-03 05:32:03Z fenu $
2// Author: Dmitry.Naumov   2005/09/26
3
4/*****************************************************************************
5 * ESAF: Euso Simulation and Analysis Framework                              *
6 *                                                                           *
7 *  Id: EnergyFitUtils                                                       *
8 *  Package: Energy                                                          *
9 *  Coordinator: Dmitry.Naumov                                               *
10 *                                                                           *
11 *****************************************************************************/
12
13//_____________________________________________________________________________
14//
15// EnergyFitUtils
16//
17//
18#include "EnergyFitUtils.hh"
19#include "EnergyModule.hh"
20#include "EOpticsResponse.hh"
21#include "OpticsResponseBuilder.hh"
22#include "RecoShowerTrack.hh"
23#include "RecoEvent.hh"
24#include "RecoPixelData.hh"
25
26#include <TMinuit.h>
27#include <TF1.h>
28#include <TProfile.h>
29using namespace TMath;
30using namespace sou;
31
32//_____________________________________________________________________________
33Double_t Saturation(Int_t pe) {
34    //
35    //
36    //
37    Double_t dead_time = EnergyModule::Get()->GetFeDeadTime();
38    Double_t gtu       = EnergyModule::Get()->GetRecoEvent()->GetHeader().GetDetector()->GetElectronics()->GetGtuLength();
39    return pe*Exp(-pe*dead_time/gtu);
40}
41
42//_____________________________________________________________________________
43Int_t UnfoldSaturation(Int_t pe) {
44    //
45    //
46    //
47    Int_t guess = pe;
48    Double_t chi2min=10000;
49    for (Int_t i=pe; i<=250; i++) {
50           Double_t chi2 = TMath::Abs(Saturation(i)-pe);
51           if (chi2 < chi2min) {
52               chi2min = chi2;
53               guess = i;
54           }
55    }
56    return guess;
57}
58//_____________________________________________________________________________
59Double_t Fit_Energy(Double_t *x, Double_t *par) {
60
61    Float_t time    = x[0];
62    Float_t energy  = par[0];
63    Float_t tmax    = par[1];
64    Float_t sigma   = par[2];
65
66    EnergyModule *fEnergyModule = (EnergyModule*)gMinuit->GetObjectFit();
67    if (!fEnergyModule) {
68      TF1::RejectPoint();
69      return 0;
70    }
71
72    Float_t Phi   = fEnergyModule->PhiFoV(time);
73    Float_t Theta = fEnergyModule->ThetaFoV(time);
74    Float_t Efficacy = fEnergyModule->GetEfficacy(Theta,Phi);
75
76    if (Theta*RadToDeg() > 30.) {
77      TF1::RejectPoint();
78      return 0;
79    }
80
81    if (time >= *(fEnergyModule->GetFCtime()+1) && 
82        time <= *(fEnergyModule->GetFCtime()+3)) {
83      TF1::RejectPoint();
84      return 0;
85    }
86
87
88    GtuPixels_t GtuPixels = fEnergyModule->GetMyGtuPixels();
89    Pixels_t OneGtuPixels = GtuPixels[Int_t(time)];
90    Pixels_t::iterator itPix;
91    Float_t background(0);
92    for (itPix = OneGtuPixels.begin(); itPix != OneGtuPixels.end(); itPix++) {
93         background += fEnergyModule->BgrPerPixelPerGtu((*itPix)->GetPixelId(),time);
94    }
95    Float_t          fQuantumEff = fEnergyModule->GetQuantumEff();
96    RecoShowerTrack *fTrack      = fEnergyModule->GetRecoShowerTrack();
97
98    const RecoShowerStep& s = fTrack->GetStepThetaPhi(Theta,Phi);
99    Float_t fAttenuation  = s.GetAttenuation();
100    Float_t fYield        = s.GetFluoYield();
101    Float_t fShowerLength = s.GetLength();
102    Float_t fRmax         = (s.GetDir()).Mag();
103
104    // Use Gauss
105    Double_t var = Power((time-tmax)/sigma,2);
106    Float_t pe = Exp(-0.5*var)*energy*1.e14*fQuantumEff*Efficacy*fAttenuation*fYield*fShowerLength/(1450.*4.e0*Pi()*fRmax*fRmax) +
107                 background;
108    return pe;
109}
110
111
112//_____________________________________________________________________________
113Double_t TimeDistr(Double_t *x, Double_t *par) {
114
115    Float_t time    = x[0];
116    Float_t energy  = par[0];
117    Float_t tmax    = par[1];
118    Float_t sigma   = par[2];
119
120    EnergyModule *fEnergyModule = (EnergyModule*)gMinuit->GetObjectFit();
121    if (!fEnergyModule) {
122      TF1::RejectPoint();
123      return 0;
124    }
125
126    Float_t Phi   = fEnergyModule->PhiFoV(time);
127    Float_t Theta = fEnergyModule->ThetaFoV(time);
128
129    if (Theta*RadToDeg() > 30.) {
130      TF1::RejectPoint();
131      return 0;
132    }
133
134    Float_t          fQuantumEff = fEnergyModule->GetQuantumEff();
135    RecoShowerTrack *fTrack      = fEnergyModule->GetRecoShowerTrack();
136
137    const RecoShowerStep& s = fTrack->GetStepThetaPhi(Theta,Phi);
138    Float_t fAttenuation    = s.GetAttenuation();
139    Float_t fYield          = s.GetFluoYield();
140    Float_t fShowerLength   = s.GetLength();
141    Float_t fRmax           = (s.GetDir()).Mag();
142
143    // Use Gauss
144    Double_t var = Power((time-tmax)/sigma,2);
145    Float_t pe = Exp(-0.5*var)*energy*1.e14*fQuantumEff*fAttenuation*fYield*fShowerLength/(1450.*4.e0*Pi()*fRmax*fRmax);
146    return pe;
147}
148
149//_____________________________________________________________________________
150Double_t Background(Double_t *x, Double_t *par) {
151    //
152    //
153    //
154
155    Float_t time    = x[0];
156    EnergyModule *fEnergyModule = (EnergyModule*)gMinuit->GetObjectFit();
157    if (!fEnergyModule) {
158      TF1::RejectPoint();
159      return 0;
160    }
161
162    Float_t Theta = fEnergyModule->ThetaFoV(time);
163
164    if (Theta*RadToDeg() > 30.) {
165      TF1::RejectPoint();
166      return 0;
167    }
168
169    GtuPixels_t GtuPixels = fEnergyModule->GetMyGtuPixels();
170    Pixels_t OneGtuPixels = GtuPixels[Int_t(time)];
171    Pixels_t::iterator itPix;
172    Float_t background(0);
173    for (itPix = OneGtuPixels.begin(); itPix != OneGtuPixels.end(); itPix++) {
174         background += fEnergyModule->BgrPerPixelPerGtu((*itPix)->GetPixelId(),time);
175    }
176
177    return background;
178}
179
180
181Double_t FunctionFit(Double_t *v, Double_t *par){
182  Double_t f_test=par[0]/1.45e9*exp((v[0]/37.15)-(par[1]/37.15)-2*(v[0]/37.15)*log(2*(v[0]/37.15)/((v[0]/37.15)+(par[1]/37.15))));
183  return f_test;
184}
185
Note: See TracBrowser for help on using the repository browser.