Demo entry 6838797

ana

   

Submitted by sof on Jun 11, 2019 at 20:42
Language: C++. Code size: 5.4 kB.

#include <Riostream.h>
#include <stdlib.h>
#include <TROOT.h>
#include <TSystem.h>
#include "TNtuple.h"
#include "TFile.h"
#include "TTree.h"
#include "TCanvas.h"
#include "TGraph.h"
#include "TGraphErrors.h"
#include "TF1.h"

#ifndef LABO
#define LABO
struct acqParams_t {
	UInt_t		acqTimestamp;
	UShort_t	dppVersion;
	UShort_t	nsPerSample;
	UShort_t	nsPerTimetag;
	UShort_t	numChannels;
	UShort_t	numSamples[32];
};

struct psdParams_t {
	UInt_t		channel;
	UInt_t		threshold;
	UInt_t		pretrigger;
	UInt_t		pregate;
	UInt_t		shortgate;
	UInt_t		longgate;
	UInt_t		numsamples;
};

struct acqEventPSD_t {
    ULong64_t	timetag;
	UInt_t		baseline;
	UShort_t	qshort;
	UShort_t	qlong;
	UShort_t	pur;
	UShort_t	samples[4096];
};
#endif

int AnaBraggDgtz(const char *filename, int ch,  int nsig=0) {

    if (ch<0 || ch>7) return -1;
  
    acqParams_t params;
    acqEventPSD_t inc_data1;

    TFile *infile = new TFile(filename);
    TTree *tree   = (TTree*)infile->Get("acq_tree_0");

    TBranch *br_inc_setup = tree->GetBranch("acq_params");
    br_inc_setup->SetAddress(&params.acqTimestamp);
    br_inc_setup->GetEntry(0);
  
    printf("acqTimestamp = %d \n",params.acqTimestamp);
    printf("nsPerSample = %d \n",params.nsPerSample);
    printf("nsPerTimetag = %d \n",params.nsPerTimetag);
    printf("numChannels = %d \n",params.numChannels);
    printf("numSamples[%d] = %d \n",ch,params.numSamples[ch]);
  
    char st[100];
    sprintf(st,"acq_ch%d",ch);
    TBranch *branch1 = tree->GetBranch(st);
    if (!branch1) {printf("no branch!\n");return -1;}

    branch1->SetAddress(&inc_data1.timetag);
    printf("Entries in data branch: %lld\n", branch1->GetEntries());

    int nev = branch1->GetEntries();

    int blfrom = 0, blto = 600; // regione per il calcolo della baseline
    float vmax; // massimo relativo alla bl
    float bl; //baseline
    float deltaT; //larghezza temporale
    float sum, sum1, sum2, sum3, sum4;
    float somma, somma1, somma2, somma3, somma4;
    float a, b, c, d;
    float delta, DELTA;
    float x_min, x_max;

    char outfilename[200];
    strcpy(outfilename,"anabragg_");
    const char *cc=strrchr(filename,'/');
    if (cc) {cc++; strcat(outfilename,cc);}
    else strcat(outfilename,filename);

    TFile *fout=new TFile(outfilename,"RECREATE")

    TNtuple *nt=new TNtuple("nt","","ev:vmax:baseline:deltaT");

    int maxev=nev;
    if (nsig && nsig<nev) maxev=nsig;

    TGraph * h = new TGraph(0);
    TGraph * f = new TGraph(0);
    TCanvas *cps= NULL;
    TGraph *g = NULL;
    char tit[50];
    sprintf(tit,"canvas_ch%d",ch);
    if (gROOT->FindObject(tit)) cps=(TCanvas*)gROOT->FindObject(tit);
    cps = new TCanvas(tit,tit,0,0,640,480);


// LOOP SUGLI EVENTI
for (int i=0; i<maxev; i++) {

    branch1->GetEntry(i);
    bl=0;
    vmax=0;

    //filtraggio dei segnali
    for(int j=0; j<3999; j++) {
        if(i==0) h->SetPoint(h->GetN(),j,inc_data1.samples[j+1]);
        inc_data1.samples[j+1]=(inc_data1.samples[j]+inc_data1.samples[j+2])/2;
        if(i==0) f->SetPoint(f->GetN(),j,inc_data1.samples[j+1]);
        j+=2;
    }

    // calcolo baseline
    for (int j=blfrom; j<blto; j++) {
        bl += inc_data1.samples[j]; bl /= (blto-blfrom);
    }

    // calcolo vmax
    for (int j=600; j<2500; j++) {
        if ( (inc_data1.samples[j]-bl)>vmax ) vmax = (inc_data1.samples[j]-bl);
    }

    // calcolo larghezza temporale
    for(int j=1300; j<1399; j++) {
        sum1 += pow(j,2);
        sum += j;
        sum3 += inc_data1.samples[j];
        sum4 += (j*inc_data1.samples[j]);
    }
    sum2 = pow(sum,2);
    delta = 100*sum1-sum2;

    a = (1/delta)*(sum1*sum3-sum*sum4);
    b = (1/delta)*(100*sum4-sum*sum3);

    x_min = (bl-a)/b;

    for(int j=1800; j<1899; j++) {
        somma1 += pow(j,2);
        somma += j;
        somma3 += inc_data1.samples[j];
        somma4 += (j*inc_data1.samples[j]);
    }
    somma2 = pow(somma,2);
    DELTA = 100*somma1-somma2;

    c = (1/DELTA)*(somma1*somma3-somma*somma4);
    d = (1/DELTA)*(100*somma4-somma*somma3);

    x_max = (vmax+bl-c)/d;

    deltaT= (x_max-x_min)/250;

    //controlli qualita'
    float dv[500]={0};
    int k=0;
    if (bl<8600 && deltaT>2.8 && deltaT<3.0) {
        float s, s1, s2, s3, s4, m, de;
        // calcolo della derivata temporale
        int u=0, h=125;

        while (k<32) {
            s=0, s1=0, s2=0, s3=0, s4=0, m=0, de=0;
            for (int j=u; j<h; j++) {
                s1 += pow(j,2);
                s += j;
                s3 += inc_data1.samples[j];
                s4 += (j*inc_data1.samples[j]);
            }

        s2 = pow(s,2);
        de = 125*s1-s2;
        m = (1/de)*(125*s4-s*s3);
        dv[k]=m;

        u += 125;
        h += 125;
        ++k;
        }

        g = new TGraph(0);

        if(i<250){
            for(int j=0; j<32; j++) {
                g->SetPoint(g->GetN(), j, dv[j]);
            }

            g->SetLineColor(38);
            g->Draw();
            cps->Update();
            gSystem->ProcessEvents();
        }

        nt->Fill(i,vmax,bl,deltaT);
    }

}

    h->Write("prima");
    f->Write("smooth");
    fout->Write();
    fout->Close();
    infile->Close();
    new TFile(outfilename);

    return 0;
}

This snippet took 0.01 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).