Subversion Repositories f9daq

Rev

Blame | Last modification | View Log | RSS feed

  1. // gSystem->AddLinkedLibs(" -lz");
  2. // .L CAENV729Read.cxx+
  3. #include <stdlib.h>
  4. #include <stdio.h>
  5. #include <zlib.h>
  6. #include <TCanvas.h>
  7. #include <TH2F.h>
  8. #include <TGraph.h>
  9. #include <TH1F.h>
  10. #include <TF1.h>
  11. #include <TText.h>
  12. #include <TLine.h>
  13. #include <math.h>
  14.  
  15. double sqr(double x){ return x*x; }
  16.  
  17. Double_t myfunction(Double_t *xx, Double_t *p)
  18.    {
  19.       Float_t x =xx[0];
  20.       double s[5],m[5],v[5];
  21.       double a=p[2];
  22.       double b=p[3];
  23.       const int npar=4;
  24.       for (int i=0;i<npar;i++) {
  25.         m[i]=a+b*i;
  26.         v[i]=p[4+2*i];
  27.         s[i]=p[4+2*i+1];
  28.       }
  29.      
  30.  
  31.       Double_t f = p[0]+p[1]*x+v[0]*exp(-0.5*sqr( (x-m[0])/s[0]) )
  32.                               +v[1]*exp(-0.5*sqr( (x-m[1])/s[1]) )
  33.                               +v[2]*exp(-0.5*sqr( (x-m[2])/s[2]) )
  34.                               +v[3]*exp(-0.5*sqr( (x-m[3])/s[3]) );
  35.       return f;
  36.    }
  37.  
  38. int CAENV729Read(char *fname,int opt=1,int start=310, int end=335, int debug=0){
  39.  
  40. // open file
  41.   TGraph *gr=NULL;
  42.   TCanvas *c = new TCanvas("c1","c1",1500,800);
  43.   c->Divide(2,2);
  44.   int first=1;
  45.  
  46.   TH2F *h2d;
  47.   TH1F *h1d;
  48.   int maxx=4096;
  49.   int maxy=4096;
  50.   double x[maxx], y[maxx];
  51.  
  52.  
  53.   unsigned short hdr[2];
  54.   unsigned short *data= new unsigned short[100000];
  55.   gzFile fp=gzopen(fname,"r");
  56.   if ( fp==NULL ){
  57.     printf("FILE %s cannot be opened\n",fname);
  58.     return 0;
  59.   }
  60.   int cnt=0;
  61.   while (!gzeof(fp)){
  62.     gzread(fp,hdr,4);
  63.     cnt++;
  64.     if (debug) printf("***** %d %d %d\n", cnt, hdr[0], hdr[1]);
  65.     int nb= hdr[1]-4;
  66.     //continue;
  67.     do {
  68.       gzread(fp,hdr,4);
  69.       int nitems = (hdr[1]-4)/2;
  70.       unsigned int id=hdr[0];
  71.       if (debug) printf("\t->> %d %d nitems=%d\n", hdr[0], hdr[1],nitems);
  72.       if (id>4) break;
  73.       gzread(fp,data,2*nitems);
  74.       nb -= hdr[1];
  75.      
  76.       if (first) {
  77.         TH2F *frame = new TH2F("frame","",5,0,nitems,5,0,maxy);
  78.         h2d=new TH2F("h2d",fname,nitems,-0.5,nitems-0.5,maxy,-0.5,maxy-0.5);
  79.         h1d=new TH1F("h1d","Charge;signal charge(a.u.);N",500,-maxy,5*maxy);
  80.         for (int i=0;i<4;i++){
  81.           c->cd(i+1);  
  82.           frame->DrawCopy();          
  83.         }
  84.         first=0;
  85.       }
  86.       double baseline=0;
  87.       for (int i=0; i<200;i++) baseline += data[i];
  88.       baseline /=200.;
  89.       double charge=0;
  90.       for (int i=0; i<nitems;i++) {
  91.         //printf("ADC %d %d\n", i, data[i]);
  92.         x[i]=i;
  93.         y[i]=data[i];
  94.         if (id==2) y[i]-=baseline;
  95.         if (i>start && i<end) charge += y[i];
  96.         if (id==2) y[i]+=2000;
  97.         if (id==2) h2d->Fill(x[i],y[i]);
  98.        
  99.       }
  100.       if (id==2) h1d->Fill(charge);
  101.       if (opt==0){
  102.         c->cd(id);
  103.         gr = new TGraph(nitems, x,y);
  104.         gr->Draw("lsame");
  105.       }
  106.     } while (nb>0);
  107.   }
  108.   printf("***** %d Events read\n", cnt);
  109.   gzclose(fp);
  110.   if (opt){
  111.  
  112.     c->cd(1);
  113.     h2d->GetXaxis()->SetRangeUser(280.,440.); // Set the range
  114.     double maximum = h2d->GetBinContent(h2d->GetMaximumBin());
  115.     h2d->SetMaximum(maximum/4);
  116.     h2d->Draw("colz");
  117.  
  118.     c->cd(2)->SetLogy();;
  119.     h1d->Draw();
  120. /*
  121.     const int npar=12;
  122.     TF1 *f1 = new TF1("f1", myfunction, h1d->GetXaxis()->GetXmin(), h1d->GetXaxis()->GetXmax(),npar);
  123.     f1->SetParameter(0,0);
  124.     f1->SetParameter(1,0);
  125.    
  126.     f1->SetParameter(2,0);    // pozicija prvega vrha
  127.     f1->SetParameter(3,5000); // razdalja med vrhoma
  128.     // visine vrhov
  129.     f1->SetParameter(4,550);
  130.     f1->SetParameter(6,300);
  131.     f1->SetParameter(8,100);
  132.     f1->SetParameter(10,40);
  133.  
  134.     // sigme gausovk
  135.     f1->SetParameter(5,50);
  136.     f1->SetParameter(7,50);
  137.     f1->SetParameter(9,50);
  138.     f1->SetParameter(11,50);
  139.  
  140.     h1d->Fit(f1, "R");
  141. */
  142.  
  143.  
  144.     c->cd(3);
  145.  
  146.     TH1F *hintegral = (TH1F*)h1d->Clone("hintegral");
  147.     double sum=0;
  148.     for (int i=hintegral->GetNbinsX();i>0;i--){
  149.        sum+=hintegral->GetBinContent(i);  
  150.        hintegral->SetBinContent(i,sum);
  151.     }
  152.     hintegral->SetTitle("Integrated charge - threshold scan; a.u.; N");
  153.     hintegral->Draw();
  154.  
  155.     c->cd(4)->SetLogy();
  156.     TH1D *h=h2d->ProjectionY("_py",323,323,"d");
  157.     h->GetXaxis()->SetRangeUser(1900.,2900.); // Set the range
  158.     h->SetTitle("Pulse height @maximum;Pulse height(a.u.);N");
  159.  
  160. //    double limits[6]={ -2000, 1000,5200,10000,14500,18000 };
  161.     double limits[6]={ 1912, 2050,2280,2490,2700,2900 };
  162.     char strtxt[255];
  163.      printf("%s\t",fname );
  164.     for (int i=0;i<5;i++) {
  165.       sprintf(strtxt,"%g", h->Integral( h->FindBin(limits[i]), h->FindBin(limits[i+1]) )/cnt);
  166.       printf("%g\t",  h->Integral( h->FindBin(limits[i]), h->FindBin(limits[i+1]) )/cnt);  
  167.       TText *txt = new TText(limits[i]+10,200-i*30, strtxt);
  168.       TLine * line= new TLine(limits[i+1],0,limits[i+1],250);
  169.       line->Draw();
  170.       txt->Draw();
  171.     }
  172.     printf(" %g\n", h->Integral(h->FindBin(limits[0]),h->FindBin(limits[5]))/cnt);
  173.  
  174.    
  175.   }
  176.  
  177.   delete data;
  178.   c->Modified();
  179.   c->Update();
  180.   char pngname[255];
  181.   sprintf(pngname,"%s.png",fname);
  182.   c->SaveAs(pngname);
  183.   return 0;
  184. }
  185.