ROCsoft: spadicplugin/go4/TSpadicProc.cxx (r4864/r3723) Source File

Go to the documentation of this file.
00001 #include "TSpadicProc.h"
00002 
00003 #include "Riostream.h"
00004 
00005 #include "TH1.h"
00006 #include "TH2.h"
00007 #include "TROOT.h"
00008 #include "TSystem.h"
00009 #include "TCutG.h"
00010 #include "TGo4Log.h"
00011 #include "TGo4WinCond.h"
00012 #include "TGo4PolyCond.h"
00013 #include "TGo4CondArray.h"
00014 #include "TGo4Picture.h"
00015 #include "TGo4MbsEvent.h"
00016 #include "TGo4FitDataHistogram.h"
00017 #include "TGo4FitParameter.h"
00018 
00019 #include "TSpadicControl.h"
00020 #include "TSpadicModel.h"
00021 
00022 #include "spadic/Message.h"
00023 
00024 #include "base/commons.h"
00025 
00026 //***********************************************************
00027 TSpadicProc::TSpadicProc() :
00028    TGo4EventProcessor()
00029 {
00030    cout << "**** TSpadicProc: Create instance " << endl;
00031 }
00032 //***********************************************************
00033 TSpadicProc::~TSpadicProc()
00034 {
00035    TGo4Log::Info("TSpadicProc: Delete instance");
00036 }
00037 //***********************************************************
00038 // this one is used in standard factory
00039 TSpadicProc::TSpadicProc(const char* name) :
00040    TGo4EventProcessor(name)
00041 {
00042    TGo4Log::Info("TSpadicProc: Create instance %s", name);
00043 
00045    // Creation of parameters (new Go4 wrapper will check if restored from auto save file):
00046    fControl = (TSpadicControl*) (MakeParameter("SpadicControl", "TSpadicControl", "set_SpadicControl.C"));
00047 
00048    // the following block will always use the parameter settings as defined in the macro
00049    // note that you can always produce a new setter macro with the current parameter values
00050    // by typing ".x saveparam.C("*","set")" in the analysis terminal macro line
00051 
00052    if (fControl)
00053       fControl->PrintParameter();
00054 
00055    TString obname;
00056    TString obtitle;
00057    TString foldername;
00058 
00059    UInt_t sid, ch;
00060    for (sid = 1; sid < MAX_SPADIC; sid++)
00061    {
00062       obname.Form("Spadic_%2d/FitGateLeft_%2d", sid, sid);
00063       obtitle.Form("Lower baseline gate %2d", sid);
00064       fBaselineGate1[sid] = MakeWinCond(obname.Data(), 0, 5);
00065       obname.Form("Spadic_%2d/FitGateRight_%2d", sid, sid);
00066       obtitle.Form("Upper baseline gate %2d ", sid);
00067       fBaselineGate2[sid] = MakeWinCond(obname.Data(), 30, 45);
00068       for (ch = 0; ch < N_CHA; ch++)
00069       {
00070          obname.Form("Spadic_%2d/Channel_%2d/PeakHeightGate_%2d_%2d", sid,
00071                ch, sid, ch);
00072          fPeakHeightGate[sid][ch] = MakeWinCond(obname.Data(), 0, 250);
00073       }
00074    }
00075 
00076    for (sid = 1; sid < MAX_SPADIC; sid++)
00077    {
00078       for (ch = 0; ch < N_CHA; ch++)
00079       {
00080          // trace histogram:
00081          obname.Form("Spadic_%2d/Channel_%2d/TRACE_Spadic:%2d_CHAN:%2d",
00082                sid, ch, sid, ch);
00083          obtitle.Form("Trace %2d  %2d", sid, ch);
00084          h_trace[sid][ch] = MakeTH1('I', obname.Data(), obtitle.Data(),
00085                TRACE_SIZE, 0, TRACE_SIZE, "time (40ns)");
00086 
00087          obname.Form("Spadic_%2d/Channel_%2d/FITTRACE_Spadic:%2d_CHAN:%2d",
00088                sid, ch, sid, ch);
00089          obtitle.Form("FitTrace %2d  %2d", sid, ch);
00090          h_fittrace[sid][ch] = MakeTH1('I', obname.Data(), obtitle.Data(),
00091                TRACE_SIZE, 0, TRACE_SIZE, "time (40ns)");
00092 
00093          // peak heightfrom bin
00094          obname.Form("Spadic_%2d/Channel_%2d/PeakHeight_Spadic:%2d_CHAN%2d",
00095                sid, ch, sid, ch);
00096          obtitle.Form("Peak %2d %2d ", sid, ch);
00097          h_peak[sid][ch] = MakeTH1('I', obname.Data(), obtitle.Data(), 300,
00098                0., 300., "pulseheight");
00099 
00100          // baseline from window
00101          obname.Form("Spadic_%2d/Channel_%2d/Baseline_Spadic:%2d_CHAN:%2d",
00102                sid, ch, sid, ch);
00103          obtitle.Form("Baseline %2d %2d ", sid, ch);
00104          h_baseline[sid][ch] = MakeTH1('I', obname.Data(), obtitle.Data(),
00105                300, 0., 300., "pulseheight");
00106 
00107          // mean position from window
00108          obname.Form("Spadic_%2d/Channel_%2d/PulseMean_Spadic:%2d_CHAN:%2d",
00109                sid, ch, sid, ch);
00110          obtitle.Form("Peak position %2d %2d ", sid, ch);
00111          h_meanpos[sid][ch] = MakeTH1('I', obname.Data(), obtitle.Data(),
00112                TRACE_SIZE, 0., TRACE_SIZE, "time (40ns)");
00113 
00114          obname.Form(
00115                "Spadic_%2d/Channel_%2d/NoiseSigma_Spadic:%2d_CHAN:%2d",
00116                sid, ch, sid, ch);
00117          obtitle.Form("Nose sigma %2d %2d ", sid, ch);
00118          h_noisesigma[sid][ch] = MakeTH1('I', obname.Data(), obtitle.Data(),
00119                100, 0., 50., "ampl");
00120 
00121          // pulse maximum position from fit
00122          obname.Form(
00123                "Spadic_%2d/Channel_%2d/FitPulsePos_Spadic:%2d_CHAN:%2d",
00124                sid, ch, sid, ch);
00125          obtitle.Form("Fitted peak position %2d  %2d", sid, ch);
00126          h_fit_pos[sid][ch] = MakeTH1('I', obname.Data(), obtitle.Data(),
00127                TRACE_SIZE, 0., TRACE_SIZE, "time (40ns)");
00128 
00129          // Overview pictures for each channel:
00130          obname.Form("Overview Spadic%2d CHAN %2d", sid, ch);
00131          p_overview[sid][ch] = GetPicture(obname.Data());
00132          if (p_overview[sid][ch] == 0)
00133          {
00134             obtitle.Form("Overview Spadic%2d CHAN %2d", sid, ch);
00135             foldername.Form("Spadic_%2d/", sid);
00136             p_overview[sid][ch] = new TGo4Picture(obname.Data(),
00137                   obtitle.Data());
00138             p_overview[sid][ch]->SetDivision(3, 2);
00139             p_overview[sid][ch]->Pic(0, 0)->AddObject(h_peak[sid][ch]);
00140             p_overview[sid][ch]->Pic(0, 0)->AddCondition(
00141                   fPeakHeightGate[sid][ch]);
00142             p_overview[sid][ch]->Pic(1, 0)->AddObject(h_baseline[sid][ch]);
00143             p_overview[sid][ch]->Pic(1, 0)->SetFillAtt(5, 3001);
00144             p_overview[sid][ch]->Pic(2, 0)->AddObject(h_noisesigma[sid][ch]);
00145             p_overview[sid][ch]->Pic(0, 1)->AddObject(h_meanpos[sid][ch]);
00146             p_overview[sid][ch]->Pic(0, 1)->SetFillAtt(4, 3001);
00147             p_overview[sid][ch]->Pic(1, 1)->AddObject(h_trace[sid][ch]);
00148             p_overview[sid][ch]->Pic(1, 1)->AddObject(h_fittrace[sid][ch]);
00149             p_overview[sid][ch]->Pic(1, 1)->AddCondition(
00150                   fBaselineGate1[sid]);
00151             p_overview[sid][ch]->Pic(1, 1)->AddCondition(
00152                   fBaselineGate2[sid]);
00153             p_overview[sid][ch]->Pic(2, 1)->AddObject(h_fit_pos[sid][ch]);
00154             AddPicture(p_overview[sid][ch], foldername.Data());
00155          }
00156 
00157       }// for channels
00158 
00159       obname.Form("Traces_Spadic%2d", sid);
00160       p_traces[sid] = GetPicture(obname.Data());
00161       if (p_traces[sid] == 0)
00162       {
00163          obtitle.Form("Pulse Traces of Spadic%2d", sid);
00164          foldername.Form("Spadic_%2d/", sid);
00165          p_traces[sid] = new TGo4Picture(obname.Data(), obtitle.Data());
00166          p_traces[sid]->SetDivision(2, 4);
00167          p_traces[sid]->Pic(0, 0)->AddObject(h_trace[sid][0]);
00168          p_traces[sid]->Pic(0, 0)->AddObject(h_fittrace[sid][0]);
00169          p_traces[sid]->Pic(0, 1)->AddObject(h_trace[sid][1]);
00170          p_traces[sid]->Pic(0, 1)->AddObject(h_fittrace[sid][1]);
00171          p_traces[sid]->Pic(0, 2)->AddObject(h_trace[sid][2]);
00172          p_traces[sid]->Pic(0, 2)->AddObject(h_fittrace[sid][2]);
00173          p_traces[sid]->Pic(0, 3)->AddObject(h_trace[sid][3]);
00174          p_traces[sid]->Pic(0, 3)->AddObject(h_fittrace[sid][3]);
00175          p_traces[sid]->Pic(1, 0)->AddObject(h_trace[sid][4]);
00176          p_traces[sid]->Pic(1, 0)->AddObject(h_fittrace[sid][4]);
00177          p_traces[sid]->Pic(1, 1)->AddObject(h_trace[sid][5]);
00178          p_traces[sid]->Pic(1, 1)->AddObject(h_fittrace[sid][5]);
00179          p_traces[sid]->Pic(1, 2)->AddObject(h_trace[sid][6]);
00180          p_traces[sid]->Pic(1, 2)->AddObject(h_fittrace[sid][6]);
00181          p_traces[sid]->Pic(1, 3)->AddObject(h_trace[sid][7]);
00182          p_traces[sid]->Pic(1, 3)->AddObject(h_fittrace[sid][7]);
00183          AddPicture(p_traces[sid], foldername.Data());
00184       }
00185 
00186       obname.Form("RawSpectra_Spadic%2d", sid);
00187       p_spectra[sid] = GetPicture(obname.Data());
00188       if (p_spectra[sid] == 0)
00189       {
00190          obtitle.Form("Peakheight spectrum of Spadic%2d", sid);
00191          foldername.Form("Spadic_%2d/", sid);
00192          p_spectra[sid] = new TGo4Picture(obname.Data(), obtitle.Data());
00193          p_spectra[sid]->SetDivision(2, 4);
00194          p_spectra[sid]->Pic(0, 0)->AddObject(h_peak[sid][0]);
00195          p_spectra[sid]->Pic(0, 1)->AddObject(h_peak[sid][1]);
00196          p_spectra[sid]->Pic(0, 2)->AddObject(h_peak[sid][2]);
00197          p_spectra[sid]->Pic(0, 3)->AddObject(h_peak[sid][3]);
00198          p_spectra[sid]->Pic(1, 0)->AddObject(h_peak[sid][4]);
00199          p_spectra[sid]->Pic(1, 1)->AddObject(h_peak[sid][5]);
00200          p_spectra[sid]->Pic(1, 2)->AddObject(h_peak[sid][6]);
00201          p_spectra[sid]->Pic(1, 3)->AddObject(h_peak[sid][7]);
00202          AddPicture(p_spectra[sid], foldername.Data());
00203       }
00204 
00205    }
00206 
00207    printf("Histograms created \n");
00208    fflush ( stdout);
00209 
00210 }
00211 //-----------------------------------------------------------
00212 // event function
00213 Bool_t TSpadicProc::BuildEvent(TGo4EventElement* target)
00214 {
00215 
00216    fOutput = (TSpadicEvent*) target;
00217    fOutput->SetValid(kFALSE);
00218    fInput = (TGo4MbsEvent*) GetInputEvent();
00219    if (fInput == 0)
00220    {
00221       cout << "AnlProc: no input event !" << endl;
00222       return kFALSE;
00223    }
00224    if (fInput->GetTrigger() > 11)
00225    {
00226       cout << "**** TSpadicProc: Skip trigger event" << endl;
00227       return kFALSE;
00228    }
00229 
00230    if (fControl->fEnableDebug) // debug message id numbers here:
00231    {
00232       cout << "Mbs Event:" << fInput->GetCount() << endl;
00233    }
00234 
00235    fInput->ResetIterator();
00236    TGo4MbsSubEvent* psubevt;
00237    while ((psubevt = fInput->NextSubEvent()) != 0) // loop over subevents
00238    {
00239       if (psubevt->GetProcid() != roc::proc_TRD_Spadic)
00240          continue; // only evaluate spadic events here
00241 
00242       UInt_t sid = psubevt->GetSubcrate(); // susibo id number
00243       if (sid > MAX_SPADIC)
00244       {
00245          cout
00246                << "**** TSpadicProc: Skipping subevent with invalid Susibo id "
00247                << sid << ", maximum allowed is " << MAX_SPADIC << endl;
00248          continue;
00249       }
00250       spadic::Message mess((uint8_t*) psubevt->GetDataField());
00251       if (!mess.CheckMessage())
00252       {
00253          cout
00254                << "**** TSpadicProc: Skipping subevent with invalid message, sid: "
00255                << sid << endl;
00256          continue;
00257       }
00258       if (fControl->fEnableDebug) // debug message id numbers here:
00259       {
00260          cout << " - Message sid:" << sid << ", status:"
00261                << mess.GetStatusNumber();
00262          cout << ", evid:" << mess.GetEventIDNumber() << ", ts:"
00263                << mess.GetTimeStamp() << endl;
00264       }
00265 
00266       for (UInt_t ch = 0; ch < N_CHA; ch++)
00267       {
00268          // generate sample trace histograms:
00269          TH1* tracehis = h_trace[sid][ch]; // the original sample
00270          tracehis->Reset("");
00271          TH1* fittracehis = h_fittrace[sid][ch]; // the fit shape (optional)
00272          for (UInt_t bin = 0; bin < TRACE_SIZE; bin++)
00273          {
00274             uint8_t val = mess.Sample(ch, bin);
00275 //            if (fControl->fEnableDebug) // debug payload
00276 //                  {
00277 //                     cout << "   - val("<<ch<<","<<bin<<")="<< (int) val<<endl;
00278 //                  }
00279             fOutput->fPulse[sid][ch][bin] = val;
00280             tracehis->SetBinContent(bin+1, (int) val);
00281          }
00282          // find peak and fill histogram (only if enabled in control parameter)
00283          Double_t baseline(0), ampl(0), position(0), sigma(0), posfit(0);
00284 
00285          if (DoModelFit(sid, ch, tracehis, fittracehis, baseline, ampl,
00286                position, sigma, posfit))
00287          {
00288 
00289             fOutput->fShapeHeight[sid][ch] = ampl;
00290             fOutput->fShapeCentroid[sid][ch] = position;
00291             fOutput->fShapeBaseline[sid][ch] = baseline;
00292             fOutput->fFitCentroid[sid][ch] = posfit;
00293 
00294             h_peak[sid][ch]->Fill(ampl);
00295             h_meanpos[sid][ch]->Fill(position);
00296             h_baseline[sid][ch]->Fill(baseline);
00297             h_noisesigma[sid][ch]->Fill(sigma);
00298             h_fit_pos[sid][ch]->Fill(posfit);
00299 
00300             fOutput->SetValid(kTRUE); // only store events that were fitted correctly at least in one channel
00301          }
00302 
00303       }// for ch
00304 
00305 
00306    } //while((psubevt = fInput->NextSubEvent())
00307 
00308 
00309    //fOutput->SetValid(kTRUE); // enable this line to store complete events, not only fitted valid JAM
00310    return kTRUE;
00311 }
00312 
00313 Bool_t TSpadicProc::DoModelFit(UInt_t i, UInt_t k, TH1* hist, TH1* fithist,
00314       Double_t& baseline, Double_t& ampl, Double_t& pos, Double_t& sigma,
00315       Double_t& posfit)
00316 {
00317 //   if (!fControl->fEnableFit)
00318 //      return false;
00319    TGo4Fitter fitter("Fitter", TGo4Fitter::ff_chi_square, kTRUE);
00320 
00321    // add histogram to fitter, which should be fitted
00322    TGo4FitDataHistogram* data = fitter.AddH1("data1", hist, kFALSE,
00323          fBaselineGate1[i]->GetXLow(), fBaselineGate1[i]->GetXUp());
00324    data->SetRange(0, fBaselineGate2[i]->GetXLow(), fBaselineGate2[i]->GetXUp());
00325 
00326    // create polynom of first order
00327    fitter.AddPolynomX("data1", "Pol", 1);
00328 
00329    fitter.SetMemoryUsage(2);
00330 
00331    // execute all actions
00332    //fitter.DoActions();
00333    fitter.EstimateAmplitudes(1);
00334 
00335    // draw data, full model and two gaussains on new canvas
00336    //   fitter.Draw("#data1");
00337 
00338    double k0 = fitter.GetParValue("Pol_0.Ampl");
00339    double k1 = fitter.GetParValue("Pol_1.Ampl");
00340 
00341    double maxdiff(0);
00342    int bin = int(fBaselineGate1[i]->GetXUp());
00343    double maxpos = bin;
00344 
00345    baseline = k0 + bin * k1;
00346 
00347    while (bin <= fBaselineGate2[i]->GetXLow())
00348    {
00349       double diff = TMath::Abs(k0 + bin * k1 - hist->GetBinContent(bin + 1));
00350 
00351       if (diff > maxdiff)
00352       {
00353          maxdiff = diff;
00354          maxpos = bin;
00355       }
00356       bin++;
00357       //cout << "bin = " << bin << "  value = " << func << endl;
00358    }
00359 
00360    double sum0(0), sum1(0), sum2(0);
00361    while (bin <= fBaselineGate2[i]->GetXUp())
00362    {
00363       double diff = k0 + bin * k1 - hist->GetBinContent(bin + 1);
00364       sum0 += 1.;
00365       sum1 += diff;
00366       sum2 += diff * diff;
00367 
00368       bin++;
00369    }
00370 
00371    sigma = TMath::Sqrt(sum2 / sum0 - TMath::Power(sum1 / sum0, 2));
00372 
00373    pos = maxpos;
00374    ampl = maxdiff;
00375 
00376    if (!fControl->fEnableFit || !fPeakHeightGate[i][k]->Test(maxdiff))
00377    {
00378       for (int n = 0; n < TRACE_SIZE; n++)
00379          fithist->SetBinContent(n + 1, k0 + n * k1);
00380       return kFALSE;
00381    }
00382 
00383    TGo4Fitter fitter2("Fitter2", TGo4Fitter::ff_chi_square, kTRUE);
00384 
00385    fitter2.AddH1("data1", hist, kFALSE, fBaselineGate1[i]->GetXLow(),
00386          fBaselineGate2[i]->GetXUp());
00387 
00388    fitter2.AddPolynomX("data1", "Pol", 1);
00389 
00390    TSpadicModel* model = new TSpadicModel("Peak");
00391    fitter2.AddModel("data1", model);
00392 
00393    model->SetParsValues(maxdiff * 25, 3., 12., maxpos - 15);
00394 
00395    model->FindPar("N")->SetFixed(kTRUE);
00396    fitter2.SetParValue("Pol_0.Ampl", k0);
00397    fitter2.SetParValue("Pol_1.Ampl", k1);
00398 
00399    fitter2.SetMemoryUsage(2);
00400 
00401    //   fitter2.Print("Pars");
00402 
00403    fitter2.DoActions();
00404 
00405    //   fitter2.Print("Pars");
00406 
00407    double xx = fBaselineGate1[i]->GetXUp();
00408    maxdiff = 0.;
00409    maxpos = 0;
00410    while (xx < fBaselineGate2[i]->GetXLow())
00411    {
00412       xx += 0.1;
00413       double val = fabs(model->Evaluate(xx));
00414       if (val > maxdiff)
00415       {
00416          maxdiff = val;
00417          maxpos = xx;
00418       }
00419    }
00420 
00421    for (int n = 0; n < TRACE_SIZE; n++)
00422       fithist->SetBinContent(n + 1, k0 + n * k1 + model->Evaluate(n));
00423 
00424    //cout << "Amplitude = " << maxdiff << "  maxpos = " << maxpos << endl;
00425 
00426 
00427    ampl = maxdiff;
00428    pos = maxpos;
00429 
00430    // evaluate posfit:
00431    Double_t tau = model->GetParValue("Tau");
00432    Double_t shift = model->GetParValue("Shift");
00433    //cout << "tau = " << tau << "  shift = " << shift << endl;
00434    posfit = tau + shift;
00435    return kTRUE;
00436 }