CheckOFShape.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; 
 
 
    double fs = 1;
    int N = 2048*fs;
    int Ne = 1e3;
    double t0 = 200*fs;
    double taur = 10*fs;
    double taud = 200*fs;
    double taur_dev = 1;
    double taud_dev = 1;
    double sigma = 1*sqrt(fs);
    double amp = 10;
 
    QVector an(N);
    an.Initialize(sigma*sigma*N);
    QVector ap(N);
    ap.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();
   double max = ap.GetMax()-ap[0];
   ap/=max;
   max = 1;
 
   QOptimumFilter of(ap,an,-1,false);
   QOptimumShapeFilter of2(ap,an,-1,false,2);
   QVector ap1filtered2 = of2.GetOriginalAveragePulseFilteredShifted();
   double cbase = ap1filtered2.GetMean(t0*3./4);
   double cgain = ap1filtered2.GetMax()-ap1filtered2.GetMin();
  // std::cout<<"cgain: "<<cgain<<" cbase: "<<cbase<<std::endl;
   double ofbase = ap1filtered2.GetMean(3./4*t0);
   ap1filtered2 -= ofbase;
   ap1filtered2 /= cgain; 
   int maxpos = ap1filtered2.GetMaxIndex();
   int deltat50 = maxpos;
   int deltat80 = maxpos;
   int deltat20 = maxpos;
   for(int i = maxpos; i > 0; i--) {
        if((ap1filtered2.GetMax()-ap1filtered2[i])<0.8) {
            deltat20 = maxpos-i;
        }
        if((ap1filtered2.GetMax()-ap1filtered2[i])<0.5) {
            deltat50 = maxpos-i;
        }
        if((ap1filtered2.GetMax()-ap1filtered2[i])<0.2) {
            deltat80 = maxpos-i;
        }
 
   }
   /*
   deltat80 += maxpos-ap.GetMaxIndex();
   deltat50 += maxpos-ap.GetMaxIndex();
   deltat20 += maxpos-ap.GetMaxIndex();
   */
    std::cout<<deltat20<<" "<<deltat50<<" "<<deltat80<<std::endl;
 
 
   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_of2<<std::endl;;
 
   TH1D* ahisto = new TH1D("ahisto","ahisto",sqrt(Ne),amp*0.2,amp*1.5);
   TH1D* lhisto = new TH1D("lhisto","lhisto",sqrt(Ne),0,50);
   TH1D* rhisto = new TH1D("rhisto","rhisto",sqrt(Ne),0,50);
   TH1D* chisto = new TH1D("chisto","chisto",sqrt(Ne),0,2);
   TH1D* jhisto = new TH1D("jhisto","jhisto",240,-30,30);
   TH1D* shisto = new TH1D("shisto","shisto",sqrt(Ne),0,amp);
 
   TH1D* ahisto2 = new TH1D("ahisto2","ahisto2",sqrt(Ne),amp*0.2,amp*1.5);
   TH1D* lhisto2 = new TH1D("lhisto2","lhisto2",sqrt(Ne),0,50);
   TH1D* rhisto2 = new TH1D("rhisto2","rhisto2",sqrt(Ne),0,50);
   TH1D* chisto2 = new TH1D("chisto2","chisto2",sqrt(Ne),0,2);
   TH1D* jhisto2 = new TH1D("jhisto2","jhisto2",50*taur,-5*taur,5*taur);
   TH1D* shisto2 = new TH1D("shisto2","shisto2",sqrt(Ne),0,amp);
 
   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);
       QVector noise(N),noise2(N);
       signal.Initialize(0);
       signal2.Initialize(0);
       noise.Initialize(0);
       noise2.Initialize(0);
 
       for(int i = 0; i < N; i++) {
           double t = double(i-t0)-shift;
           noise[i] += rand.Gaus(0,sigma);
           noise2[i] += rand.Gaus(0,sigma);
           if(t < 0) continue;
           signal[i] += (-exp(-t/taur) + exp(-t/taud));
           signal2[i] +=(-exp(-t/taur/taur_dev) + exp(-t/taud/taud_dev));
       }
       signal *= amp/signal.GetMax();
       signal2 *= amp/signal2.GetMax()/1.15;
       signal += noise;
       signal2 += noise2;
      
 
       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);
 
       err = of.Filter(signal2);
       if(err != QERR_SUCCESS) continue;
       err = of.SetJitter(QOptimumFilter::J_LOCMAX);
       if(err != QERR_SUCCESS) continue;
       double jitter2,chi22,amplitude2,integral2,left2,right2;
       err = of.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);
 
       success++;
 
 
       double fac = 1.5;
       err = of2.Filter(signal);
       if(err != QERR_SUCCESS) continue;
       double t1=0,t2=0;
       for(int i = maxpos+jitter; i >0; i--) {
           if(amplitude-of2.GetFilteredShifted()[i] < 0.25*amplitude) t2 = i;
           if(amplitude-of2.GetFilteredShifted()[i] < 0.75*amplitude) t1 = i;
       }
       //shisto->Fill(t2-t1);
       double diff = of2.GetFilteredShifted()[maxpos+jitter-deltat80]-of2.GetFilteredShifted()[maxpos+jitter-deltat20];
       shisto->Fill(diff/cgain);
 
       err = of2.Filter(signal2);
       if(err != QERR_SUCCESS) continue;
 
       t1=0,t2=0;
       for(int i = maxpos+jitter2; i >0; i--) {
           if(amplitude2-of2.GetFilteredShifted()[i] < 0.25*amplitude) t2 = i;
           if(amplitude2-of2.GetFilteredShifted()[i] < 0.75*amplitude) t1 = i;
       }
       //shisto2->Fill(t2-t1);
 
       double diff2 = of2.GetFilteredShifted()[maxpos+jitter2-deltat80]-of2.GetFilteredShifted()[maxpos+jitter2-deltat20];
      shisto2->Fill(diff2/cgain);
 
       success2++;
 
       if(success2 == 1) {
           TCanvas * cs = new TCanvas("cs","cs",800,600);
           signal2.GetGraph()->Draw("AL");
           QVector ofwave_man = of2.GetFilteredShifted();
           ofbase = ofwave_man.GetMean(3./4*t0);
           ofwave_man -= ofbase;
           ofwave_man /= cgain; 
 
           TGraph* ofwave = ofwave_man.GetGraph();
           TGraph* offilter = (of2.GetFilterTD()/of2.GetFilterTD().GetMax()*amp).GetGraph();
 
           ofwave->SetLineColor(kRed);
           ofwave->Draw("LSAME");
           TGraph* apgraph = (ap.Shift(shift)*(amp/max)).GetGraph();
           apgraph->SetLineColor(kBlue);
           apgraph->Draw("LSAME");
           offilter->SetLineColor(kGreen);
           offilter->Draw("LSAME");
           cs->SaveAs("checkof_wave.pdf");
           return;
       }
 
 
 
   }
   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->cd(6);
   shisto2->SetLineColor(kRed);
   shisto2->Draw();
   shisto->Draw("SAME");
   TF1* fg = new TF1("fg","gaus");
   shisto->Fit(fg);
   double m1,s1,m2,s2;
   m1 = fg->GetParameter(1);
   s1 = fg->GetParameter(2);
   shisto2->Fit(fg);
   m2 = fg->GetParameter(1);
   s2 = fg->GetParameter(2);
   double dp = (m1-m2)/sqrt(s1*s1+s2*s2);
   std::cout<<"DP: "<<dp<<std::endl;
 
   c1->SaveAs("checkof.pdf");
 
 
 }