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;
}