Changeset 14476 for fact/tools


Ignore:
Timestamp:
10/15/12 14:27:55 (12 years ago)
Author:
Jens Buss
Message:
calculated model attributes
File:
1 edited

Legend:

Unmodified
Added
Removed

  • fact/tools/rootmacros/PulseTemplates/pulse.C

    r14472 r14476  
    1010Pulse::Pulse()
    1111{
     12    InitMembers();
     13}
     14
     15Pulse::Pulse(TString name)
     16{
     17    InitMembers();
     18    mName           = name;
     19}
     20
     21Pulse::Pulse(TString name, TH1* histo)
     22{
     23    InitMembers();
     24    mName           = name;
     25    mHisto          = histo;
     26    Fit(mName, mOptions, mType);
     27}
     28
     29Pulse::Pulse(TString name, TH1* histo, int type)
     30{
     31    InitMembers();
     32    mName           = name;
     33    mHisto          = histo;
     34    mType           = type;
     35    Fit(mName, mOptions, mType);
     36}
     37
     38Pulse::Pulse(TString name, TH1* histo, TString options)
     39{
     40    InitMembers();
     41    mName           = name;
     42    mHisto          = histo;
     43    mOptions        += options;
     44    Fit(mName, mOptions, mType);
     45}
     46
     47Pulse::Pulse(TString name, TH1* histo, TString options, int type)
     48{
     49    InitMembers();
     50    mName           = name;
     51    mHisto          = histo;
     52    mOptions        += options;
     53    mType           = type;
     54    Fit(mName, mOptions, mType);
     55}
     56
     57
     58Pulse::~Pulse()
     59{
     60    delete mModel;
     61}
     62
     63void
     64Pulse::InitMembers()
     65{
     66    mName           = "pulse";
    1267    mHisto          = NULL;
    13     mOptions        = "";
     68    mModel          = NULL;
     69    mOptions        = "S";
    1470    mBsl            = 0;
    1571    mHeight         = 0;
    16     mStart          = 0;
    17     mRising         = 0;
    18     mTauRising      = 0;
    19     mTauFalling     = 0;
     72    mT0             = 0;
     73    mT1             = 0;
     74    mTau1           = 0;
     75    mTau2           = 0;
    2076    mIntegral       = 0;
    2177    mAmplitude      = 0;
    2278    mPhE            = 0;
    23 }
    24 
    25 Pulse::Pulse(TH1* histo)
    26 {
    27     mHisto          = histo;
    28     mOptions        = "";
    29     mBsl            = 0;
    30     mHeight         = 0;
    31     mStart          = 0;
    32     mRising         = 0;
    33     mTauRising      = 0;
    34     mTauFalling     = 0;
    35     mIntegral       = 0;
    36     mAmplitude      = 0;
    37     mPhE            = 0;
    38 }
    39 
    40 Pulse::Pulse(TH1* histo, TString options)
    41 {
    42     mHisto          = histo;
    43     mOptions        = options;
    44     mBsl            = 0;
    45     mHeight         = 0;
    46     mStart          = 0;
    47     mRising         = 0;
    48     mTauRising      = 0;
    49     mTauFalling     = 0;
    50     mIntegral       = 0;
    51     mAmplitude      = 0;
    52     mPhE            = 0;
    53 }
    54 
    55 Pulse::~Pulse()
    56 {
    57 
     79    mType           = 0;
     80    mFitProb        = 0;
     81    mFitNCalls      = 0;
     82    mFitNdf         = 0;
     83    mChi2           = 0;
    5884}
    5985
    … …  
    6187Pulse::Fit(
    6288        TString fitName,
    63         TString fitOptions
     89        TString fitOptions,
     90        int type
    6491        )
    6592{
    … …  
    7097        fitName,
    7198        fitOptions,
     99        type,
    72100        fitMin,
    73101        fitMax
    74     );
     102        );
    75103}
    76104
    77105void
    78106Pulse::Fit(
     107        TString fitName,
     108        TString fitOptions,
     109        int type,
     110        int fitMin,
     111        int fitMax
     112        )
     113{
     114    fitOptions += "R";
     115    if (type == 0)
     116        FitSectionWise(
     117                    "Model0",
     118                    fitOptions,
     119                    fitMin,
     120                    fitMax
     121                    );
     122    else if (type == 1)
     123        FitContious(
     124                    "Model1",
     125                    fitOptions,
     126                    fitMin,
     127                    fitMax
     128                    );
     129    else{
     130        cout << "not a correct type number --> fitting skipped!" << endl;
     131        return;
     132    }
     133
     134    CalculateParameters();
     135    return;
     136}
     137
     138void
     139Pulse::FitSectionWise(
    79140        TString fitName,
    80141        TString fitOptions,
    … …  
    83144        )
    84145{
    85     TF1 fit(fitName, shapeFunc, fitMin, fitMax, 6 );
     146    mFitMin = fitMin;
     147    mFitMax = fitMax;
     148    mModel  = new TF1(fitName, shapeFunc, fitMin, fitMax, 6 );
    86149
    87150    // ======================================================================
    … …  
    92155    int first_bin   = mHisto->GetXaxis()->GetFirst();
    93156
    94     for (int i = 0; i < 10; i++){
     157    for (int i = 0; i < 30; i++){
    95158        bsl += mHisto->GetBinContent(first_bin+0);
    96159    }
    97     bsl = bsl/10;
     160    bsl = bsl/30;
    98161
    99162    double stop     = mHisto->GetMaximumBin(); //pos of max
    100163    double height   = mHisto->GetBinContent(stop);
    101164    double start    = stop-10; //pos 10 slices before maximum
    102     double rising   = (stop - start)/tau;
    103165    // ======================================================================
    104166
    … …  
    107169
    108170
    109     fit.SetParameters(para);
    110     fit.SetParNames("BSL", "Height", "Start", "Rising", "Tau1", "TauFalling");
    111     fit.SetLineColor(kRed);
    112     mHisto->Fit(&fit, fitOptions);
    113 
    114     mBsl            = fit.GetParameter(0);
    115     mHeight         = fit.GetParameter(1);
    116     mStart          = fit.GetParameter(2);
    117     mRising         = fit.GetParameter(3);
    118     mTauRising      = fit.GetParameter(4);
    119     mTauFalling     = fit.GetParameter(5);
    120 }
    121 
    122 void
    123 Pulse::Fit2(
    124         TString fitName,
    125         TString fitOptions
    126         )
    127 {
    128     int     fitMin = mHisto->GetXaxis()->GetFirst();
    129     int     fitMax = mHisto->GetXaxis()->GetLast();
    130 
    131     Fit2(
    132         fitName,
    133         fitOptions,
    134         fitMin,
    135         fitMax
    136     );
    137 }
    138 
    139 void
    140 Pulse::Fit2(
     171    mModel->SetParameters(para);
     172    mModel->SetParNames("BSL", "A0", "t0", "t1", "Tau1", "Tau2");
     173    mModel->SetLineColor(kRed);
     174
     175    mFitResultPtr = mHisto->Fit(mModel, fitOptions);
     176
     177    mBsl            = mModel->GetParameter(0);
     178    mHeight         = mModel->GetParameter(1);
     179    mT0             = mModel->GetParameter(2);
     180    mT1             = mModel->GetParameter(3);
     181    mTau1           = mModel->GetParameter(4);
     182    mTau2           = mModel->GetParameter(5);
     183}
     184
     185void
     186Pulse::FitContious(
    141187        TString fitName,
    142188        TString fitOptions,
    … …  
    145191        )
    146192{
    147     TF1 fit(fitName, shapeFunc2, fitMin, fitMax, 5 );
     193    mFitMin = fitMin;
     194    mFitMax = fitMax;
     195    mModel  = new TF1(fitName, shapeFunc2, fitMin, fitMax, 5 );
    148196
    149197    // ======================================================================
    … …  
    154202    int first_bin   = mHisto->GetXaxis()->GetFirst();
    155203
    156     for (int i = 0; i < 10; i++){
     204    for (int i = 0; i < 30; i++){
    157205        bsl += mHisto->GetBinContent(first_bin+0);
    158206    }
    159     bsl = bsl/10;
     207    bsl = bsl/30;
    160208
    161209    double stop     = mHisto->GetMaximumBin(); //pos of max
    … …  
    166214
    167215    double para[] = {bsl, amplitude, start, tau, tau};
    168     fit.SetParameters(para);
    169     fit.SetParNames("BSL", "Amplitude", "Start", "Tau1", "TauFalling");
    170     fit.SetLineColor(kBlue);
    171     mHisto->Fit(&fit, fitOptions);
    172 
    173     mBsl            = fit.GetParameter(0);
    174     mAmplitude      = fit.GetParameter(1);
    175     mStart          = fit.GetParameter(2);
    176     mTauRising      = fit.GetParameter(3);
    177     mTauFalling     = fit.GetParameter(4);
     216    mModel->SetParameters(para);
     217    mModel->SetParNames("BSL", "A0", "t0", "Tau1", "Tau2");
     218    mModel->SetLineColor(kBlue);
     219
     220    mFitResultPtr = mHisto->Fit(mModel, fitOptions);
     221
     222    mBsl    = mModel->GetParameter(0);
     223    mHeight = mModel->GetParameter(1);
     224    mT0     = mModel->GetParameter(2);
     225    mTau1   = mModel->GetParameter(3);
     226    mTau2   = mModel->GetParameter(4);
    178227}
    179228
    180229double  Pulse::GetBsl(){ return mBsl;}
    181230double  Pulse::GetHeight(){ return mHeight;}
    182 double  Pulse::GetAvalancheStart(){ return mStart;}
    183 double  Pulse::GetAvalancheEnd(){ return mRising;}
    184 double  Pulse::GetTimeConstRising(){ return mTauRising;}
    185 double  Pulse::GetTimeConstFalling(){ return mTauFalling;}
     231double  Pulse::GetAvalancheStart(){ return mT0;}
     232double  Pulse::GetAvalancheEnd(){ return mT1;}
     233double  Pulse::GetTimeConstRising(){ return mTau1;}
     234double  Pulse::GetTimeConstFalling(){ return mTau2;}
    186235double  Pulse::GetIntegral(){ return mIntegral;}
    187236double  Pulse::GetAmplitude(){ return mAmplitude;}
    188237int     Pulse::GetPE(){ return mPhE;}
    189238
     239void
     240Pulse::CalculateParameters()
     241{
     242    mIntegral   = mModel->Integral(mFitMin, mFitMax);
     243    mAmplitude  = mModel->GetMaximum() - mBsl;
     244    mFitProb    = mFitResultPtr->Prob();
     245    mFitNCalls  = mFitResultPtr->NCalls();
     246    mFitNdf     = mFitResultPtr->Ndf();
     247    mChi2       = mFitResultPtr->Chi2();
     248}
     249
     250
    190251//double
    191252int Heaviside(double val)
    192253{
    193     if( val < 0)    return 0;
    194254    if( val >= 0)   return 1;
     255    else return 0;
    195256}
    196257
    … …  
    214275    // helper variables
    215276//    double max      = start + rising * tau1;
    216     double e1       = height * (1 - Exp(-(t[0]-start)/tau1 ) );
     277    double e1       = (1 - Exp(-(t[0]-start)/tau1 ) );
    217278//    double e2       = (-1) * height * (1 - Exp(-(x[0]-max)/tau2 ) );
    218     double e2       = (-1) * height * (1 - Exp(-(t[0]-stop)/tau2 ) );
     279    double e2       = (-1) * (1 - Exp(-(t[0]-stop)/tau2 ) );
    219280    // calculate return value
     281    returnval += Heaviside(t[0]-start)*e1;
     282    returnval += Heaviside(t[0]-stop)*e2;
     283    returnval *= height;
    220284    returnval += bsl;
    221 
    222     if (t[0] > start)
    223         returnval += e1;
    224     if (t[0] > stop)
    225         returnval += e2;
     285//    if (t[0] > start)
     286//        returnval += e1;
     287//    if (t[0] > stop)
     288//        returnval += e2;
    226289
    227290    return returnval;
    228 }
     291
     292
     293}
     294
     295
    229296
    230297/*
    … …  
    282349    return returnval;
    283350}
    284 
Note: See TracChangeset for help on using the changeset viewer.