CheckOF.C

Go to the documentation of this file.

 {
 //   gSystem->Load("/opt/local/libexec/boost/1.76/lib/libboost_thread-mt");
 //   gSystem->Load("/opt/local/lib/libcrypto");
 //   gSystem->Load("/opt/local/lib/libssl");
    gROOT->Reset();
    gSystem->Load("libpq");
        gSystem->Load("libqroot");
    gSystem->Load("libthreadwrapper");
    gSystem->Load("libbase");
    gSystem->Load("libcoretools");
    gSystem->Load("libglobalrw");
    gSystem->Load("libqstyle");
    gSystem->Load("libdbbase");
    gSystem->Load("librundatadb");
    gSystem->Load("libdianadata");
    gSystem->Load("libcalderdata");
    gSystem->Load("libcaldertools");
    gSystem->Load("libqstyleloader");
    gSystem->Load("liboptimumfilter");
    using namespace Diana;
        TF1 a; QVector b; 
 
 
    int N = 2048;
    int Ne = 1e3;
    double t0 = 200;
    double taur = 10;
    double taud = 200;
    double taur_dev = 1.0;
    double taud_dev = 1.0;
    double sigma = 1;
    double amp = 100; // in sigma_of
 
    QVector an(N);
    an.Initialize(sigma*sigma*N);
    QVector ap(N);
    ap.Initialize(0);
    QVector ap2(N);
    ap2.Initialize(0);
 
   for(int i = 0; i < N; i++) {
        double t = i-t0;
        if(t < 0) continue;
        ap[i] += -exp(-t/taur) + exp(-t/taud);
   }
   int apmaxindex = ap.GetMaxIndex();
   for(int i = apmaxindex; i > 0; i--) {
        ap2[i] = ap[i];
        ap2[2*apmaxindex-i] = ap[i];
   }
   ap2 = ap;
 
    double max = ap.GetMax()-ap[0];
 
    QOptimumFilter of(ap,an,-1,false);
    QOptimumFilter of2(ap2,an,-1,false);
    double sigma_of = of.GetFilteredNoiseRMS();
    double sigma_of2 = of2.GetFilteredNoiseRMS();
    amp *= sigma_of;
    std::cout<<"Simga OF: "<<sigma_of<<" amp: "<<amp<<std::endl;
    std::cout<<"Simga OF2: "<<sigma_of<<std::endl;;
 
    TH1D* ahisto = new TH1D("ahisto","ahisto",sqrt(Ne),amp-5*sigma_of,amp+5*sigma_of);
    TH1D* lhisto = new TH1D("lhisto","lhisto",sqrt(Ne),0,10);
    TH1D* rhisto = new TH1D("rhisto","rhisto",sqrt(Ne),0,10);
    TH1D* chisto = new TH1D("chisto","chisto",sqrt(Ne),0,2);
    TH1D* jhisto = new TH1D("jhisto","jhisto",240,-30,30);
 
    TH1D* ahisto2 = new TH1D("ahisto2","ahisto2",sqrt(Ne),amp-5*sigma_of,amp+5*sigma_of);
    TH1D* lhisto2 = new TH1D("lhisto2","lhisto2",sqrt(Ne),0,10);
    TH1D* rhisto2 = new TH1D("rhisto2","rhisto2",sqrt(Ne),0,10);
    TH1D* chisto2 = new TH1D("chisto2","chisto2",sqrt(Ne),0,2);
    TH1D* jhisto2 = new TH1D("jhisto2","jhisto2",50*taur,-5*taur,5*taur);
 
    TRandom rand;
    rand.SetSeed(0);
 
 
    int success = 0, success2 = 0;
    for(int e = 0; e < Ne; e++) {
        if(e%100 == 0) cout<<"Event: "<<e<<std::endl;
 
        double shift=rand.Uniform(-4*taur,4*taur);
        QVector signal(N),signal2(N);
        signal.Initialize(0);
        signal2.Initialize(0);
        for(int i = 0; i < N; i++) {
            double t = double(i-t0)-shift;
            signal[i] += rand.Gaus(0,sigma);
            signal2[i] += rand.Gaus(0,sigma);
            if(t < 0) continue;
            signal[i] += amp/max*(-exp(-t/taur) + exp(-t/taud));
            signal2[i] += amp/max*(-exp(-t/taur/taur_dev) + exp(-t/taud/taud_dev));
        }
 
 
 
        //return;
        QError err;
        err = of.Filter(signal);
        if(err != QERR_SUCCESS) continue;
        err = of.SetJitter(QOptimumFilter::J_LOCMAX);
        if(err != QERR_SUCCESS) continue;
        double jitter,chi2,amplitude,integral,left,right;
        err = of.GetInterpolated(jitter,chi2,amplitude,integral,left,right);
        if(err != QERR_SUCCESS) continue;
        ahisto->Fill(amplitude);
        lhisto->Fill(left);
        rhisto->Fill(right);
        chisto->Fill(chi2);
        jhisto->Fill(jitter-shift);
        success++;
 
        err = of2.Filter(signal2);
        if(err != QERR_SUCCESS) continue;
        err = of2.SetJitter(QOptimumFilter::J_FIXED,jitter);
        if(err != QERR_SUCCESS) continue;
        double jitter2,chi22,amplitude2,integral2,left2,right2;
        err = of2.GetInterpolated(jitter2,chi22,amplitude2,integral2,left2,right2);
        if(err != QERR_SUCCESS) continue;
        ahisto2->Fill(amplitude2);
        lhisto2->Fill(left2);
        rhisto2->Fill(right2);
        chisto2->Fill(chi22);
        jhisto2->Fill(jitter2-shift);
        success2++;
 
 
        if(e == 0) {
            TCanvas * cs = new TCanvas("cs","cs",800,600);
            signal.GetGraph()->Draw("AL");
            TGraph* ofwave = of.GetFilteredShifted().GetGraph();
 
            ofwave->SetLineColor(kRed);
            ofwave->Draw("LSAME");
            TGraph* apgraph = (ap.Shift(shift)*(amp/max)).GetGraph();
            apgraph->SetLineColor(kBlue);
            apgraph->Draw("LSAME");
            cs->SaveAs("checkof_wave.pdf");
        }
 
 
 
    }
    std::cout<<"Successfull filtering  : "<<double(success)/Ne*100<<" %"<<std::endl;
    std::cout<<"Successfull filtering 2: "<<double(success2)/Ne*100<<" %"<<std::endl;
 
    TCanvas* c1 = new TCanvas("checkof","checkof",1200,800);
    c1->Divide(3,2);
    c1->cd(1);
    ahisto2->SetLineColor(kRed);
    ahisto2->Draw();
    ahisto->Draw("SAME");
 
    c1->cd(2);
    jhisto2->SetLineColor(kRed);
    jhisto2->Draw();
    jhisto->Draw("SAME");
    c1->cd(3);
    chisto2->SetLineColor(kRed);
    chisto2->Draw();
    chisto->Draw("SAME");
    c1->cd(4);
    lhisto2->SetLineColor(kRed);
    lhisto2->Draw();
    lhisto->Draw("SAME");
    c1->cd(5);
    rhisto2->SetLineColor(kRed);
    rhisto2->Draw();
    rhisto->Draw("SAME");
    c1->SaveAs("checkof.pdf");
 
 
 }