Subversion Repositories f9daq

Rev

Rev 218 | Blame | Compare with Previous | Last modification | View Log | RSS feed

  1. #include <stdlib.h>
  2. #include <stdio.h>
  3. #include <TFile.h>
  4. #include <TH1F.h>
  5. #include <TH1D.h>
  6. #include <TH2F.h>
  7. #include <TStyle.h>
  8. #include <TCanvas.h>
  9. #include <TLegend.h>
  10. // Program za analizo podatkov zajetih z AITSiPMDAQ
  11.  
  12. int aitana( const char *fname="a1.dat", int range=5, float cut=300, int *offset = NULL)
  13. {
  14.         int debug=0;
  15.        //int off[4]={507,521,497,459};
  16.                 //int off[4]={198,203,190,150}; /// run 1
  17.                 int off[4]={120,120,120,120}; /// run 6
  18.                 if (offset) for (int i=0;i<4;i++){off[i]=offset[i];};
  19.                 gStyle->SetOptStat(0);
  20.  
  21.         char rootname[0xFF];
  22.                 sprintf(rootname, "%s.root", fname);
  23.         TFile *froot = new  TFile(rootname,"RECREATE");
  24.                 char histname[128];
  25.             TLegend * leg = new TLegend(0.8, 0.8, 1, 1);
  26.                 TH1F *hadc[4];
  27.                 TH1F *hadcsum;
  28.                 TH1F *hadcsumx;
  29.                 TH1F *hadcsumy;
  30.                 TH2F *hxy,*hxy1,*hxy2;
  31.                 TH2F *h2d=NULL;
  32.                 TH2F *h2dx=NULL;
  33.                 TH2F *h2dy=NULL;
  34.                 TH2F *hxcorr=NULL;
  35.                 TH2F *hycorr=NULL;
  36.                
  37.                
  38.             for (int i=0; i<4; i++) {
  39.          sprintf(histname,"adc%d",i);
  40.          hadc[i]=new TH1F(histname,histname,200,0.5,200*range+0.5);
  41.                 }
  42.                 hadcsum=new TH1F("adcsum","adcsum",200,0.5,8*200*range+0.5);
  43.                 hadcsumx=new TH1F("adcsumx","adcsumx",200,0.5,200*range+0.5);
  44.                 hadcsumy=new TH1F("adcsumy","adcsumy",200,0.5,200*range+0.5);
  45.         hxy=new TH2F("hxy","Center of gravity;ycog(a.u.);xcog(a.u.)",200,-0.1,1.1,200,-0.1,1.1);
  46.         hxy1=new TH2F("hxy1","Center of gravity for position 1;ycog(a.u.);xcog(a.u.)",200,-0.1,1.1,200,-0.1,1.1);
  47.         hxy2=new TH2F("hxy2","Center of gravity for position 2;ycog(a.u.);xcog(a.u.)",200,-0.1,1.1,200,-0.1,1.1);
  48.                 FILE *fin=fopen(fname,"rb");
  49.         if (fin==NULL) {
  50.                 printf("Error opening file %s\n",fname);
  51.                 return 0;
  52.         }
  53.        
  54.         unsigned short hdr[2];
  55.                 int pos[7]={0,0,0,0,0,0,0};
  56.                 int poshdr[9]={0,0,0,0,0,0,0,0,0};
  57.         float b[4], sum=0;
  58.                 unsigned short a[4];
  59.                 unsigned short buf[10000];
  60.                 float posx=0, posy=0;
  61.         while(!feof(fin)) {
  62.                
  63.                 int stat=fread(hdr,1,4,fin);
  64.                
  65.                 //if (hdr[0]!= 1 ) printf("#%d %d\n",hdr[0], hdr[1] );
  66.                
  67.                                 stat=fread(buf,1,hdr[1],fin);
  68.                                 if (stat!=hdr[1]){
  69.                                        
  70.                                         printf("READ ERROR!\n");
  71.                                         break;
  72.                                        
  73.                                 }
  74.                 switch (hdr[0]) {
  75.                                        
  76.                                         case 0x1:
  77.                                                    
  78.                                                    for (int i=0;i<hdr[1]/4/sizeof(unsigned short);i++){
  79.                                                            unsigned short *adc=&buf[i*4];
  80.                                                            sum=0;                                                          
  81.                                                            for (int j=0;j<4;j++){
  82.                                                                  a[j]=adc[j]&0xFFF;
  83.                                                                  //sum+=a[j];
  84.                                                                  hadc[j]->Fill(a[j]);  
  85.                                  b[j]=a[j]-off[j];
  86.                                                                  if (b[j]< 0) b[j]=0;
  87.                                  sum+=b[j];                                                              
  88.                                                            }
  89.                                                            hadcsum->Fill(sum);
  90.                                                            hadcsumx->Fill(a[0]+a[1]);
  91.                                                            hadcsumy->Fill(a[2]+a[3]);
  92.                                float y = (b[0]+b[1])? b[0]/(b[0]+b[1]) : 0;                                                    
  93.                                float x = (b[2]+b[3])? b[2]/(b[2]+b[3]) : 0;
  94.                                                            if (h2d) h2d->Fill(posx,posy,sum);
  95.                                                            if (h2dx) h2dx->Fill(posx,posy,b[0]+b[1]);
  96.                                                            if (h2dy) h2dy->Fill(posx,posy,b[2]+b[3]);
  97.                                                            
  98.                                                            if (sum>cut) {
  99.                                                                    hxy->Fill(y,x);
  100.                                    if (pos[3]==poshdr[6]/4 && pos[4]==poshdr[7]/4) hxy1->Fill(y,x);
  101.                                    if (pos[3]==poshdr[6]/4*3 && pos[4]==poshdr[7]/4*3) hxy2->Fill(y,x);
  102.                                                                    if (hxcorr) hxcorr->Fill(pos[0],y);
  103.                                                                    if (hycorr) hycorr->Fill(pos[1],x);
  104.                                                            }       
  105.                                                            
  106.                                                    }
  107.                                                    break;
  108.                                                 case 0x2:{
  109.                            int *p= (int *) &buf[0];    
  110.                                                    for (int k=0;k<7;k++) pos[k]=p[k];
  111.                            printf("Position: %d %d\n",pos[3], pos[4] );
  112.                                                    posx = poshdr[0]+pos[3]* poshdr[3];
  113.                                                    posy = poshdr[1]+pos[4]* poshdr[4];
  114.                         }                                                  
  115.                                                    break;
  116.                                                 case 0x3:{
  117.                                                 int *p= (int *) &buf[0];
  118.                         for (int k=0;k<9;k++) poshdr[k]=p[k];                                          
  119.                                                      
  120.                                 h2d=new TH2F("h2d","Position dependency of sum;x(step);y(step)",
  121.                                                 poshdr[6],poshdr[0]-poshdr[3]*0.5 ,poshdr[0]+poshdr[6]*(poshdr[3]-0.5),
  122.                                                 poshdr[7],poshdr[1]-poshdr[4]*0.5 ,poshdr[1]+poshdr[7]*(poshdr[4]-0.5)   );
  123.                
  124.                                                 h2dx=new TH2F("h2dx","Position dependency of xsum;x(step);y(step)",
  125.                                                 poshdr[6],poshdr[0]-poshdr[3]*0.5 ,poshdr[0]+poshdr[6]*(poshdr[3]-0.5),
  126.                                                 poshdr[7],poshdr[1]-poshdr[4]*0.5 ,poshdr[1]+poshdr[7]*(poshdr[4]-0.5)   );
  127.                                                
  128.                                                 h2dy=new TH2F("h2dy","Position dependency of ysum;x(step);y(step)",
  129.                                                 poshdr[6],poshdr[0]-poshdr[3]*0.5 ,poshdr[0]+poshdr[6]*(poshdr[3]-0.5),
  130.                                                 poshdr[7],poshdr[1]-poshdr[4]*0.5 ,poshdr[1]+poshdr[7]*(poshdr[4]-0.5)   );
  131.                                                
  132.                                    
  133.                                 hxcorr=new TH2F("hxcorr","Correlation stage position vs cog;x(step);xcog(a.u.)",
  134.                                             poshdr[6],poshdr[0]-poshdr[3]*0.5 ,poshdr[0]+poshdr[6]*(poshdr[3]-0.5), 100 , 0, 1);
  135.                  
  136.                                 hycorr=new TH2F("hycorr","Correlation stage position vs cog;y(step);ycog(a.u.)",
  137.                                                 poshdr[7],poshdr[1]-poshdr[4]*0.5 ,poshdr[1]+poshdr[7]*(poshdr[4]-0.5), 100 , 0, 1);
  138.                                                 }
  139.                                                 break;
  140.                                                 default: break;
  141.                                 }
  142.         }                              
  143.  
  144.         TCanvas* c= new TCanvas("c","ADC",750,50,700,700);
  145.                 sprintf(rootname, "%s.pdf", fname);
  146.                 c->Print(TString(rootname) + "[","pdf");
  147.                
  148.         c->Divide(2,3);
  149.         for (int i=0; i<4; i++) {
  150.                   c->cd(i+1)->SetLogy();
  151.                   hadc[i]->DrawCopy();
  152.             }
  153.                 c->cd(5)->SetLogy();hadcsumx->DrawCopy();
  154.                 c->cd(6)->SetLogy();hadcsumy->DrawCopy();
  155.                
  156.                 /*
  157.                 c->cd(7);hadcsumx->DrawCopy();
  158.                 hadc[0]->DrawCopy("same");
  159.                 hadc[1]->DrawCopy("same");
  160.                 c->cd(8);hadcsumy->DrawCopy();
  161.                 hadc[2]->DrawCopy("same");
  162.                 hadc[3]->DrawCopy("same");
  163.                 */
  164.                 c->Print(rootname,"pdf");
  165.  
  166.                
  167.                 c= new TCanvas("c1","ait reconstruction",750,50,700,700);
  168.                 c->Divide(2,4);
  169.                 if (!h2d){
  170.               c->Clear();
  171.                   c->Divide(2,3);
  172.                 }      
  173.        
  174.                
  175.                 hxy->SetMinimum(-1);
  176.                 hxy->SetMaximum(1000);
  177.                 c->cd(1)->SetLogz(); hxy->DrawCopy("colz");
  178.         c->cd(3)->SetLogz(); hxy1->DrawCopy("colz");
  179.         c->cd(4)->SetLogz(); hxy2->DrawCopy("colz");
  180.                 c->cd(2)->SetLogy(); hadcsum->DrawCopy();
  181.                
  182.                 if (h2d){
  183.                   hxcorr->SetMinimum(-1);
  184.                   hycorr->SetMinimum(-1);
  185.                   c->cd(7); hxcorr->DrawCopy("colz");
  186.                   c->cd(8); hycorr->DrawCopy("colz");
  187.                 }  
  188.                 TH1D *px = hxy->ProjectionX("_px",0,-1);
  189.                 c->cd(5); px->DrawCopy();
  190.          
  191.             TH1D *py = hxy->ProjectionY("_py",0,-1);
  192.                 c->cd(6); py->DrawCopy();
  193.                
  194.                 c->Modified();
  195.                 c->Update();
  196.                 c->Print(rootname,"pdf");
  197.                
  198.                
  199.                 if (h2d){
  200.                   c= new TCanvas("c2","position dependence sum, xsum, ysum",750,50,700,700);
  201.                   c->Divide(1,3);
  202.                   c->cd(1); h2d->DrawCopy("colz");
  203.                   c->cd(2); h2dx->DrawCopy("colz");
  204.                   c->cd(3); h2dy->DrawCopy("colz");
  205.                   c->Modified();
  206.                   c->Update();
  207.                   c->Print(rootname,"pdf");
  208.                 }                
  209.                
  210.                 c->Print(TString(rootname) + "]","pdf");
  211.                 froot->Write();
  212.                 froot->Close();
  213.                 if (fin) fclose(fin);
  214.        
  215.         return 0;
  216. }
  217.