Subversion Repositories f9daq

// gSystem->AddLinkedLibs(" -lz");

// .L CAENV729Read.cxx+

#include <stdlib.h>

#include <stdio.h>

#include <zlib.h>

#include <TCanvas.h>

#include <TH2F.h>

#include <TGraph.h>

#include <TH1F.h>

#include <TF1.h>

#include <TText.h>

#include <TLine.h>

#include <math.h>

double sqr(double x){ return x*x; }

Double_t myfunction(Double_t *xx, Double_t *p)

   {

      Float_t x =xx[0];

      double s[5],m[5],v[5];

      double a=p[2];

      double b=p[3];

      const int npar=4;

      for (int i=0;i<npar;i++) {

        m[i]=a+b*i;

        v[i]=p[4+2*i];

        s[i]=p[4+2*i+1];

      }

      Double_t f = p[0]+p[1]*x+v[0]*exp(-0.5*sqr( (x-m[0])/s[0]) )

                              +v[1]*exp(-0.5*sqr( (x-m[1])/s[1]) )

                              +v[2]*exp(-0.5*sqr( (x-m[2])/s[2]) )

                              +v[3]*exp(-0.5*sqr( (x-m[3])/s[3]) );

      return f;

   }

int CAENV729Read(char *fname,int opt=1,int start=310, int end=335, int debug=0){

// open file

  TGraph *gr=NULL;

  TCanvas *c = new TCanvas("c1","c1",1500,800);

  c->Divide(2,2);

  int first=1;

  TH2F *h2d;

  TH1F *h1d;

  int maxx=4096;

  int maxy=4096;

  double x[maxx], y[maxx];

  unsigned short hdr[2];

  unsigned short *data= new unsigned short[100000];

  gzFile fp=gzopen(fname,"r");

  if ( fp==NULL ){

    printf("FILE %s cannot be opened\n",fname);

    return 0;

  }

  int cnt=0;

  while (!gzeof(fp)){

    gzread(fp,hdr,4);

    cnt++;

    if (debug) printf("***** %d %d %d\n", cnt, hdr[0], hdr[1]);

    int nb= hdr[1]-4;

    //continue;

    do {

      gzread(fp,hdr,4);

      int nitems = (hdr[1]-4)/2;

      unsigned int id=hdr[0];

      if (debug) printf("\t->> %d %d nitems=%d\n", hdr[0], hdr[1],nitems);

      if (id>4) break;

      gzread(fp,data,2*nitems);

      nb -= hdr[1];

      if (first) {

        TH2F *frame = new TH2F("frame","",5,0,nitems,5,0,maxy);

        h2d=new TH2F("h2d",fname,nitems,-0.5,nitems-0.5,maxy,-0.5,maxy-0.5);

        h1d=new TH1F("h1d","Charge;signal charge(a.u.);N",500,-maxy,5*maxy);

        for (int i=0;i<4;i++){

          c->cd(i+1);  

          frame->DrawCopy();          

        }

        first=0;

      }

      double baseline=0;

      for (int i=0; i<200;i++) baseline += data[i];

      baseline /=200.;

      double charge=0;

      for (int i=0; i<nitems;i++) {

        //printf("ADC %d %d\n", i, data[i]);

        x[i]=i;

        y[i]=data[i];

        if (id==2) y[i]-=baseline;

        if (i>start && i<end) charge += y[i];

        if (id==2) y[i]+=2000;

        if (id==2) h2d->Fill(x[i],y[i]);

      }

      if (id==2) h1d->Fill(charge);

      if (opt==0){

        c->cd(id);

        gr = new TGraph(nitems, x,y);

        gr->Draw("lsame");

      }

    } while (nb>0);

  }

  printf("***** %d Events read\n", cnt);

  gzclose(fp);

  if (opt){

    c->cd(1);

    h2d->GetXaxis()->SetRangeUser(280.,440.); // Set the range

    double maximum = h2d->GetBinContent(h2d->GetMaximumBin());

    h2d->SetMaximum(maximum/4);

    h2d->Draw("colz");

    c->cd(2)->SetLogy();;

    h1d->Draw();

/*

    const int npar=12;

    TF1 *f1 = new TF1("f1", myfunction, h1d->GetXaxis()->GetXmin(), h1d->GetXaxis()->GetXmax(),npar);

    f1->SetParameter(0,0);

    f1->SetParameter(1,0);

    f1->SetParameter(2,0);    // pozicija prvega vrha

    f1->SetParameter(3,5000); // razdalja med vrhoma

    // visine vrhov

    f1->SetParameter(4,550);

    f1->SetParameter(6,300);

    f1->SetParameter(8,100);

    f1->SetParameter(10,40);

    // sigme gausovk

    f1->SetParameter(5,50);

    f1->SetParameter(7,50);

    f1->SetParameter(9,50);

    f1->SetParameter(11,50);

    h1d->Fit(f1, "R");

*/

    c->cd(3);

    TH1F *hintegral = (TH1F*)h1d->Clone("hintegral");

    double sum=0;

    for (int i=hintegral->GetNbinsX();i>0;i--){

       sum+=hintegral->GetBinContent(i);  

       hintegral->SetBinContent(i,sum);

    }

    hintegral->SetTitle("Integrated charge - threshold scan; a.u.; N");

    hintegral->Draw();

    c->cd(4)->SetLogy();

    TH1D *h=h2d->ProjectionY("_py",323,323,"d");

    h->GetXaxis()->SetRangeUser(1900.,2900.); // Set the range

    h->SetTitle("Pulse height @maximum;Pulse height(a.u.);N");

//    double limits[6]={ -2000, 1000,5200,10000,14500,18000 };

    double limits[6]={ 1912, 2050,2280,2490,2700,2900 };

    char strtxt[255];

     printf("%s\t",fname );

    for (int i=0;i<5;i++) {

      sprintf(strtxt,"%g", h->Integral( h->FindBin(limits[i]), h->FindBin(limits[i+1]) )/cnt);

      printf("%g\t",  h->Integral( h->FindBin(limits[i]), h->FindBin(limits[i+1]) )/cnt);  

      TText *txt = new TText(limits[i]+10,200-i*30, strtxt);

      TLine * line= new TLine(limits[i+1],0,limits[i+1],250);

      line->Draw();

      txt->Draw();

    }

    printf(" %g\n", h->Integral(h->FindBin(limits[0]),h->FindBin(limits[5]))/cnt);

  }

  delete data;

  c->Modified();

  c->Update();

  char pngname[255];

  sprintf(pngname,"%s.png",fname);

  c->SaveAs(pngname);

  return 0;

}