EXP2003 » EXP2003: Express analysis

#include "processData.h"

void processData() {

        Double_t amplification = 0.5;

        TChain *chCal = new TChain("tree");

        chCal->Add("/mnt/data/exp2001/data/analysed/reco/crun09_0001.lmd.reco.root");

        TChain *ch = new TChain("tree");

        ch->Add("/mnt/data/exp2001/data/analysed/reco/run*.root");

        TCanvas *c_raw = new TCanvas("c_raw","raw Data",1800,1000);

        c_raw->Divide(2,3);

        c_raw->cd(1);

        // chCal->Draw("tF5-tF3+68.475 >> h1(71,160,170)","","");

        chCal->Draw("tF5-tF3+68.475 >> h1(500,160,170)","","");

        c_raw->Update();

        hCal = (TH1F*)gPad->GetPrimitive("h1");

        c_raw->cd(2);

        // ch->Draw("tF5-tF3+68.475 >> h2(71,160,170)","","");

        ch->Draw("tF5-tF3+68.475 >> h2(200,160,170)","","");

        // gPad->SetLogy(1);

        c_raw->Update();

        hExp = (TH1F*)gPad->GetPrimitive("h2");        

        c_raw->cd(3);

        gCal = SetGraphs(hCal);

        gCal->Draw("AL*");

        c_raw->Update();

        c_raw->cd(4);

        gExp = SetGraphs(hExp);

        gExp->Draw("AL*");

        c_raw->Update();

        c_raw->cd(5);

        chCal->Draw("(he8.E()-he8.Mag())*1000","","");

        c_raw->Update();

        c_raw->cd(6);

        ch->Draw("(he8.E()-he8.Mag())*1000","","");

        c_raw->Update();

        // c_raw->cd(6);

        // gCFD_Exp = SetCFDGraph(gExp, amplification, 1.5,&tCFD_exp);

        // gCFD_Exp->Draw("AL*");

        // c_raw->Update();

        // cout << tCFD_exp << " " << gCFD_Exp->GetN() << endl;

        // c_raw->cd(5);

        // gCFD_Cal = SetCFDGraph(gCal, amplification, 1.5,&tCFD_cal);

        // gCFD_Cal->Draw("AL*");

        // cout << tCFD_cal << " "  << gCFD_Cal->GetN() << endl;

        // c_raw->Update();

        // c_raw->cd(6);

        // gCFD_Exp = SetCFDGraph(gExp, amplification, 1.5,&tCFD_exp);

        // gCFD_Exp->Draw("AL*");

        // c_raw->Update();

        // cout << tCFD_exp << " " << gCFD_Exp->GetN() << endl;

        // smooth

        TCanvas *c_smooth = new TCanvas("c_smooth","smooth Data",1800,1000);

        c_smooth->Divide(3,4);

        c_smooth->cd(1);

        hCal_smooth = (TH1F*)hCal->Clone();

        hCal_smooth->Smooth(2);

        hCal_smooth->Draw();

        c_smooth->Update();

        c_smooth->cd(2);

        hExp_smooth = (TH1F*)hExp->Clone();

        hExp_smooth->Smooth(2);

        hExp_smooth->Draw();

        c_smooth->Update();

        c_smooth->cd(3);

        hExp_smooth->Draw();

        TH1F *hCal_smooth_shift = shiftHisto(hCal_smooth,hExp_smooth, 0);

        hCal_smooth_shift->SetLineColor(kRed);

        hCal_smooth_shift->Draw("same");

        c_smooth->Update();

        // graphs

        c_smooth->cd(4);

        gCal_smooth = SetGraphs(hCal_smooth);

        gCal_smooth->Draw("AL");

        c_smooth->Update();

        c_smooth->cd(5);

        gExp_smooth = SetGraphs(hExp_smooth);

        gExp_smooth->Draw("AL");

        c_smooth->Update();

        // CFD

        c_smooth->cd(7);

        Double_t timeCal;

        TGraph *gCal_smooth_CFD = SetCFDGraph(gCal_smooth, amplification, 1.5 ,&timeCal);

        cout << "CFDtime on lvl " << amplification << " for Calibration data: " << timeCal << endl;

        gCal_smooth_CFD->Draw("AL*");

        TMarker *mCal_cfd = new TMarker(timeCal, gCal_smooth_CFD->Eval(timeCal), 3);

        mCal_cfd->SetMarkerColor(kRed);

        mCal_cfd->Draw("same");

        TLine *l_zero = new TLine(*(gCal_smooth_CFD->GetX()), 0, *(gCal_smooth_CFD->GetX()+gCal_smooth_CFD->GetN()-1),0);

        l_zero->SetLineColor(kRed);

        l_zero->Draw("same");

        c_smooth->Update();

        c_smooth->cd(8);

        Double_t timeExp;

        TGraph *gExp_smooth_CFD = SetCFDGraph(gExp_smooth, amplification, 1.5 ,&timeExp);

        cout << "CFDtime on lvl " << amplification << " for Experimental data: " << timeExp << endl;

        gExp_smooth_CFD->Draw("AL*");

        TMarker *mExp_cfd = new TMarker(timeExp, gExp_smooth_CFD->Eval(timeExp), 3);

        mExp_cfd->SetMarkerColor(kRed);

        mExp_cfd->Draw("same");

        l_zero->Draw("same");

        c_smooth->Update();

        // CFD

        c_smooth->cd(10);

        Double_t timeCal_low;

        TGraph *gCal_back_CFD = SetCFDGraph_back(gCal_smooth, amplification, 0.5 ,&timeCal_low);

        cout << "Low Energy CFDtime on lvl " << amplification << " for Calibration data: " << timeCal_low << endl;

        gCal_back_CFD->Draw("AL*");

        TMarker *mCal_cfd_back = new TMarker(timeCal_low, gCal_back_CFD->Eval(timeCal_low), 3);

        mCal_cfd_back->SetMarkerColor(kRed);

        mCal_cfd_back->Draw("same");

        l_zero->Draw("same");

        c_smooth->Update();

        c_smooth->cd(11);

        Double_t timeExp_low;

        TGraph *gExp_back_CFD = SetCFDGraph_back(gExp_smooth, amplification, 0.5 ,&timeExp_low);

        cout << "Low Energy CFDtime on lvl " << amplification << " for Experimental data: " << timeExp_low << endl;

        gExp_back_CFD->Draw("AL*");

        TMarker *mExp_cfd_back = new TMarker(timeExp_low, gExp_back_CFD->Eval(timeExp_low), 3);

        mExp_cfd_back->SetMarkerColor(kRed);

        mExp_cfd_back->Draw("same");

        l_zero->Draw("same");

        c_smooth->Update();

        Double_t yDev_cal,xDev_cal;

        getPeregib(gCal_smooth,timeCal_low ,0.1,&yDev_cal,&xDev_cal);

        Double_t yDev_exp,xDev_exp;

        getPeregib(gExp_smooth,timeExp_low ,0.1,&yDev_exp,&xDev_exp);

        // draw TimeStamps

        c_smooth->cd(4);

        TMarker *mCal = new TMarker(timeCal, gCal_smooth->Eval(timeCal), 3);

        mCal->SetMarkerColor(kRed);

        mCal->Draw("same");

        TMarker *mCal_back = new TMarker(timeCal_low, gCal_smooth->Eval(timeCal_low), 3);

        mCal_back->SetMarkerColor(kGreen);

        mCal_back->Draw("same");

        TMarker *mCal_dev = new TMarker(xDev_cal, gCal_smooth->Eval(xDev_cal), 3);

        mCal_dev->SetMarkerColor(kMagenta);

        mCal_dev->Draw("same");

        TLine *lCal = new TLine(*(gCal_smooth->GetX()), hCal_smooth->GetMaximum()/2, *(gCal_smooth->GetX()+gCal_smooth->GetN()-1), hCal_smooth->GetMaximum()/2);

        lCal->SetLineColor(kRed);

        lCal->Draw("same");

        c_smooth->Update();

        c_smooth->cd(5);

        TMarker *mExp = new TMarker(timeExp, gExp_smooth->Eval(timeExp), 3);

        mExp->SetMarkerColor(kRed);

        mExp->Draw("same");

        TMarker *mExp_back = new TMarker(timeExp_low, gExp_smooth->Eval(timeExp_low), 3);

        mExp_back->SetMarkerColor(kGreen);

        mExp_back->Draw("same");

        TMarker *mExp_dev = new TMarker(xDev_exp, gExp_smooth->Eval(xDev_exp), 3);

        mExp_dev->SetMarkerColor(kMagenta);

        mExp_dev->Draw("same");

        TLine *lExp = new TLine(*(gExp_smooth->GetX()), hExp_smooth->GetMaximum()/2, *(gExp_smooth->GetX()+gExp_smooth->GetN()-1), hExp_smooth->GetMaximum()/2);

        lExp->SetLineColor(kRed);

        lExp->Draw("same");

        c_smooth->Update();

        cout << "Peregib points: " << xDev_exp << " " << xDev_cal << endl;

        cout << " Bin width: " << hExp_smooth->GetBinWidth(1) << endl;

        // fit func

        TGraph *g_temp_exp = (TGraph*)gExp_smooth->Clone();

        TGraph *g_temp_exp1 = (TGraph*)gExp_smooth->Clone();

        TF1 *fitf = new TF1("fitf","[0]*x + [1]");

        fitf->SetLineColor(kRed);

        TCanvas *c_temp = new TCanvas("c_temp","title",1800,1000);

        // c_temp->Divide(1,3);

        // c_temp->cd(1);

        // gCal_smooth->Draw();

        // fitf->SetRange(167.89,168.17);

        // gCal_smooth->Fit("fitf","R");

        // c_temp->Update();

        // c_temp->cd(2);

        // gExp_smooth->Draw();

        // fitf->SetRange(167.35,167.7);

        // gExp_smooth->Fit("fitf","R");

        // c_temp->Update();

        c_temp->cd();

        // g_temp_exp->Draw();

        Double_t lin_c, sh_c;

        lin_c = 0.93;

        sh_c = 0.54;

        for (Int_t i=0;i<g_temp_exp1->GetN();i++) {

                g_temp_exp1->SetPoint(i,*(g_temp_exp1->GetX()+i)+sh_c, *(g_temp_exp1->GetY()+i)*lin_c );

        }

        g_temp_exp1->Draw();

        // TH1F *htemp = shiftHisto(hCal_smooth,hExp_smooth, -0.47);

        TGraph *g_temp = shiftGraph(gCal_smooth,gExp_smooth, 0.);

        g_temp->SetLineColor(kRed);

        g_temp->Draw("same");

        c_temp->Update();

return;

        cout << "Calibration data. Bin width: " << hCal->GetBinWidth(1) << endl;

        cout << getLed_low(hCal,10) << " ";

        cout << getLed_low(hCal,100) << " ";

        cout << getLed_low(hCal,300) << " ";

        cout << getLed_low(hCal,500) << endl << endl;

        // cout << getLed_up(hCal,100) << endl;

        cout << "EXP data. Bin width: " << hExp->GetBinWidth(1) << endl;

        cout << getLed_low(hExp,10) << " ";

        cout << getLed_low(hExp,100) << " ";

        cout << getLed_low(hExp,300) << " ";

        cout << getLed_low(hExp,500) << endl;

}

Double_t getLed_low(TH1F *hist, Double_t thresh) {

        for (Int_t i=1;i<hist->GetNbinsX()+1;i++){

                if (hist->GetBinContent(i)>thresh) {

                        return hist->GetBinCenter(i);

                }

        }

        return -1;

}

Double_t getLed_up(TH1F *hist, Double_t thresh) {

        for (Int_t i=hist->GetNbinsX();i>0;i--){

                if (hist->GetBinContent(i)>thresh) {

                        return hist->GetBinCenter(i);

                }

        }

        return -1;

}

TGraph* SetGraphs(TH1F* hist) {

        TGraph *g = new TGraph();

        for (Int_t i=0; i<hist->GetNbinsX(); i++) {

                g->SetPoint(i, hist->GetBinCenter(i+1), hist->GetBinContent(i+1));

        }

        return g;

}

TGraph* SetCFDGraph_static(TGraph *g, Double_t amp, Int_t delay) {

        Double_t stepX = *(g->GetX()+1) - *(g->GetX());

        TGraph *gCFD = new TGraph();

        gCFD->Set(g->GetN() + delay);

        for (Int_t i=0;i<gCFD->GetN();i++) {

                static Double_t yVal;

                static Double_t xVal;

                if (i<delay) {

                        yVal = *(g->GetY()+i)*amp*(-1);

                        xVal = *(g->GetX()+i)-delay*stepX;

                }

                if (i>=delay && i<g->GetN()) {

                        yVal = *(g->GetY()+i)*amp*(-1) + *(g->GetY()+i-delay);

                        xVal = *(g->GetX()+i-delay);

                }

                if (i>=g->GetN()) {

                        yVal = *(g->GetY()+i-delay-1);

                        xVal = *(g->GetX()+i-delay-1);

                }

                gCFD->SetPoint(i, xVal, yVal);

        }

        return gCFD;

}

TGraph* shiftGraph(TGraph *g1, TGraph *g2, Double_t shift) {

        Double_t max1,max2;

        max1 = *(g1->GetY());

        for (Int_t i=0;i<g1->GetN();i++) {

                if (max1<*(g1->GetY()+i)) max1 = *(g1->GetY()+i);

        }

        max2 = *(g2->GetY());

        for (Int_t i=0;i<g2->GetN();i++) {

                if (max2<*(g2->GetY()+i)) max2 = *(g2->GetY()+i);

        }

        Double_t amp = max2/max1;

        TGraph *gTemp = new TGraph();

        gTemp->Set(g1->GetN());

        for (Int_t i=0;i<gTemp->GetN();i++) {

                gTemp->SetPoint(i,*(g1->GetX()+i)+shift, *(g1->GetY()+i)*amp );

        }

        return gTemp;

}

TH1F* shiftHisto(TH1F *h1, TH1F *h2, Double_t shift) {

        Int_t binShift = shift/h1->GetBinWidth(1);

        Double_t amp = h2->GetMaximum()/h1->GetMaximum();

        TH1F *htemp = (TH1F*)h1->Clone();

        for(Int_t i=1;i<htemp->GetNbinsX()+1;i++) {

                if(i+binShift<=htemp->GetNbinsX()) htemp->SetBinContent(i,h1->GetBinContent(i-binShift)*amp);

                else {

                        htemp->SetBinContent(i,0);

                }

        }

        return htemp;

}

TGraph* SetCFDGraph_back(TGraph *g, Double_t amp, Double_t delay,Double_t* CFDtime) {

        const Double_t timeStep = 0.05;

        Double_t mytime = *(g->GetX());

        Double_t timeMin  = mytime;

        Double_t timeMax  = mytime;

        Double_t valMin  = *(g->GetY());

        Double_t valMax  = *(g->GetY());

        Int_t nBins = ( *(g->GetX()+g->GetN()-1) - *(g->GetX()) + delay)/timeStep;

        TGraph *gCFD = new TGraph();

        gCFD->Set(nBins);

        for (Int_t i=0;i<gCFD->GetN();i++) {

                static Double_t yVal;

                if ( mytime<delay+*(g->GetX()) ) {

                        yVal = g->Eval(mytime);

                }

                // if ( mytime>=*(g->GetX()) && mytime< *(g->GetX()+g->GetN()-1)+delay ) {

                if ( mytime>=delay+*(g->GetX()) && mytime< *(g->GetX()+g->GetN()-1) ) {

                        yVal = g->Eval(mytime-delay)*amp*(-1) + g->Eval(mytime);

                }

                if ( mytime>= *(g->GetX()+g->GetN()-1) ) {

                        yVal = g->Eval(mytime-delay)*amp*(-1);

                }

                if (yVal>valMax) {

                        valMax = yVal;

                        timeMax = mytime;

                }

                if (yVal<valMin) {

                        valMin = yVal;

                        timeMin = mytime;

                }

                gCFD->SetPoint(i, mytime, yVal);

                mytime += timeStep;

        }

        mytime = timeMin;

        while( (gCFD->Eval(mytime) <= 0) && (mytime >= timeMax) /*&& (time >= TminCFD)*/ ) {

                *CFDtime = mytime;

                mytime -= timeStep;

        }

        return gCFD;

}

TGraph* SetCFDGraph(TGraph *g, Double_t amp, Double_t delay,Double_t* CFDtime) {

        const Double_t timeStep = 0.05;

        Double_t mytime = *(g->GetX()) - delay;

        Double_t timeMin  = mytime;

        Double_t timeMax  = mytime;

        Double_t valMin  = *(g->GetY());

        Double_t valMax  = *(g->GetY());

        Int_t nBins = ( *(g->GetX()+g->GetN()-1) - *(g->GetX()) + delay)/timeStep;

        TGraph *gCFD = new TGraph();

        gCFD->Set(nBins);

        for (Int_t i=0;i<gCFD->GetN();i++) {

                static Double_t yVal;

                if ( mytime<*(g->GetX()) ) {

                        yVal = g->Eval(mytime+delay)*amp*(-1);

                }

                if ( mytime>=*(g->GetX()) && mytime< *(g->GetX()+g->GetN()-1)-delay ) {

                        yVal = g->Eval(mytime+delay)*amp*(-1) + g->Eval(mytime);

                }

                if ( mytime>= *(g->GetX()+g->GetN()-1)-delay ) {

                        yVal = g->Eval(mytime);

                }

                if (yVal>valMax) {

                        valMax = yVal;

                        timeMax = mytime;

                }

                if (yVal<valMin) {

                        valMin = yVal;

                        timeMin = mytime;

                }

                gCFD->SetPoint(i, mytime, yVal);

                mytime += timeStep;

        }

        mytime = timeMin;

        while( (gCFD->Eval(mytime) <= 0) && (mytime <= timeMax) /*&& (time >= TminCFD)*/ ) {

                *CFDtime = mytime;

                mytime += timeStep;

        }

        return gCFD;

}

void  getPeregib(TGraph *g, Double_t x, Double_t range,Double_t *maxDev,Double_t *XmaxDev) {

        const Double_t timeStep = 0.05;

        Double_t mytime = x - range;

        Double_t preValue = g->Eval(mytime-timeStep);

        *maxDev = fabs( (g->Eval(mytime)-g->Eval(mytime-timeStep))/timeStep );

        *XmaxDev = mytime;

        while (mytime<x+range) {

                static Double_t value;

                value = g->Eval(mytime);

                static Double_t dev; 

                dev = (value - preValue)/timeStep;

                if (dev>*maxDev) {

                        *maxDev = dev;

                        *XmaxDev = mytime;

                } 

                mytime += timeStep;

                preValue = value;

        }

        return 0;

}