/[PAMELA software]/PamelaDigitizer/DigitizeTOF.cxx
ViewVC logotype

Diff of /PamelaDigitizer/DigitizeTOF.cxx

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.2 by pamelats, Fri Jun 6 15:16:01 2008 UTC revision 1.7 by pizzolot, Fri Nov 13 09:08:53 2009 UTC
# Line 1  Line 1 
 #include <sstream>  
 #include <fstream>  
 #include <stdlib.h>  
 #include <stdio.h>  
 #include <string.h>  
 #include <ctype.h>  
 #include <time.h>  
 #include "Riostream.h"  
 #include "TFile.h"  
 #include "TDirectory.h"  
 #include "TTree.h"  
 #include "TLeafI.h"  
 #include "TH1.h"  
 #include "TH2.h"  
 #include "TF1.h"  
 #include "TMath.h"  
 #include "TRandom.h"  
 #include "TSQLServer.h"  
 #include "TSystem.h"  
 #include "CalibTrk1Event.h"  
 #include "CalibTrk2Event.h"  
 //  
1  #include "Digitizer.h"  #include "Digitizer.h"
 #include "CRC.h"  
 //  
 #include <PamelaRun.h>  
 #include <physics/calorimeter/CalorimeterEvent.h>  
 #include <CalibCalPedEvent.h>  
 #include "GLTables.h"  
2    
3  void Digitizer::DigitizeTOF(int np,float *atte1,float *atte2,float *lambda1,float *lambda2){  void Digitizer::DigitizeTOF(int np,float *atte1,float *atte2,float *lambda1,float *lambda2){
4  //fDataTof: 12 x 23 bytes (=276 bytes)  //fDataTof: 12 x 23 bytes (=276 bytes)
# Line 38  void Digitizer::DigitizeTOF(int np,float Line 10  void Digitizer::DigitizeTOF(int np,float
10                    1,1,1,1,2,2,2,3,3,3,3,4,4,4,1,  //30-44                    1,1,1,1,2,2,2,3,3,3,3,4,4,4,1,  //30-44
11                    1,2,0,2,0,0,5,5,5,5,6,6,6,6,7,  //45-59                    1,2,0,2,0,0,5,5,5,5,6,6,6,6,7,  //45-59
12                    3,3,4,4,5,5,6,7,8,9,10,11,12,13,14 };  //60-74                    3,3,4,4,5,5,6,7,8,9,10,11,12,13,14 };  //60-74
13      
14    int Z = cdp[Ipa-1];          int Z = cdp[Ipa-1];      
15    
16    float time_res[8] = {425.,210.,170.,130.,120.,120.,120.,120.};    float time_res[8] = {425.,210.,170.,130.,120.,120.,120.,120.};
17      for(Int_t i=0;i<8;i++)time_res[i]/=1.4;//1.17;1.5;1.3*/
18    Float_t dt1 = 1.e-12*time_res[0];  // single PMT resolution for Z=1  (WM, Nov'07)    Float_t dt1 = 0.;// = 1.e-12*time_res[0];  // single PMT resolution for Z=1  (WM, Nov'07)
19      
20    if ((Z > 1) && (Z < 9)) dt1=1.e-12*time_res[(Z-1)];    if ((Z > 1) && (Z < 9)) dt1=1.e-12*time_res[(Z-1)];
21    if  (Z > 8) dt1=120.e-12;    if  (Z > 8) dt1=120.e-12;
22      
23    
24  // ------ evaluate energy in each pmt: ------  // ------ evaluate energy in each pmt: ------
25  // strip geometry (lenght/width)  // strip geometry (lenght/width)
# Line 63  void Digitizer::DigitizeTOF(int np,float Line 35  void Digitizer::DigitizeTOF(int np,float
35    const Float_t echarge = 1.6e-19; // electron charge    const Float_t echarge = 1.6e-19; // electron charge
36    Float_t Npho=0.;    Float_t Npho=0.;
37    Float_t QevePmt_pC[48];    Float_t QevePmt_pC[48];
38    Float_t QhitPad_pC[2]={0., 0.};    Float_t QhitPad_pC[2]={0.,0.};
39    Float_t QhitPmt_pC[2]={0., 0.};    Float_t QhitPmt_pC[2]={0.,0.};
40    Float_t pmGain = 3.5e6;  /* PMT Gain: the same for all PMTs */    Float_t pmGain = 3.5e6;  /* PMT Gain: the same for all PMTs */
41    Float_t effi=0.21;       /* Efficienza di fotocatodo */    Float_t effi=0.21;       /* Efficienza di fotocatodo */
42  // pC < 800  // pC < 800
43    Float_t ADC_pC0A =   -4.437616e+01   ;    Float_t ADC_pC0A =-4.437616e+01;
44    Float_t ADC_pC1A =   1.573329e+00    ;    Float_t ADC_pC1A = 1.573329e+00;
45    Float_t ADC_pC2A =   2.780518e-04  ;    Float_t ADC_pC2A = 2.780518e-04;
46    Float_t ADC_pC3A =  -2.302160e-07 ;    Float_t ADC_pC3A =-2.302160e-07;
47  // pC   > 800:  // pC   > 800:
48    Float_t ADC_pC0B =    -2.245756e+02  ;    Float_t ADC_pC0B =-2.245756e+02;
49    Float_t ADC_pC1B =     2.184156e+00  ;    Float_t ADC_pC1B = 2.184156e+00;
50    Float_t ADC_pC2B =     -4.171825e-04  ;    Float_t ADC_pC2B =-4.171825e-04;
51    Float_t ADC_pC3B =     3.789715e-08  ;    Float_t ADC_pC3B = 3.789715e-08;
52    
53    Float_t pCthres=40.;     // threshold in charge    Float_t pCthres=40.;     // threshold in charge
54    Int_t ADClast=4095;      // no signal --> ADC ch=4095    Int_t ADClast=4095;      // no signal --> ADC ch=4095
55    Int_t ADCsat=3100;       // saturation value for the ADCs    Int_t ADCsat=3100;       // saturation value for the ADCs
56    Int_t ADCtof[48];    Int_t ADCtof[48];
57    Float_t ScaleFact[48]={0.39, 0.49, 0.38, 0.40, 0.65, 0.51, 0.43,    Float_t ScaleFact[48]={0.39, 0.49, 0.38, 0.40, 0.65, 0.51, 0.43, 0.49,
58                           0.49, 0.58, 0.38, 0.53, 0.57, 0.53, 0.45, 0.49, 0.16,                           0.58, 0.38, 0.53, 0.57, 0.53, 0.45, 0.49, 0.22,
59                           0.15, 0.44, 0.28, 0.57, 0.26, 0.72, 0.37, 0.29, 0.30, 0.89,                           0.21, 0.44, 0.28, 0.57, 0.26, 0.72, 0.37, 0.29,
60                           0.37, 0.08,  0.27, 0.23, 0.12, 0.22, 0.15, 0.16, 0.21,                           0.30, 0.89, 0.37, 0.12, 0.27, 0.23, 0.15, 0.22,
61                           0.19, 0.41, 0.32, 0.39, 0.38, 0.28, 0.66, 0.28, 0.40, 0.39, 0.40, 0.37, 0.35 };                           0.19, 0.20, 0.21, 0.19, 0.41, 0.32, 0.39, 0.38,
62                             0.28, 0.66, 0.28, 0.40, 0.39, 0.40, 0.37, 0.35};//15:0.7--0.95, 16:0.9--1.25, 27:0.9--1.3, 30:0.9--1.15, 32:0.85--1.05, 33:0.85--1.05
63    for(Int_t i=0; i<48; i++){    for(Int_t i=0; i<48; i++){
64      QevePmt_pC[i] = 0;      QevePmt_pC[i] = 0;
65      ADCtof[i]=0;      ADCtof[i]=0;
66    }    }
67    Int_t ip,ipad,pmtleft=0,pmtright=0,*pl,*pr;    Int_t ip,ipad,pmtleft=0,pmtright=0;
   pl = &pmtleft;  
   pr = &pmtright;  
68  // TDC variables:  // TDC variables:
69    Int_t TDClast=4095,TDCint[48];    Int_t TDClast=4095,TDCint[48];
70    Float_t  tdc[48],tdc1[48],tdcpmt[48];    Float_t  tdc[48],tdc1[48],tdcpmt[48];
# Line 114  void Digitizer::DigitizeTOF(int np,float Line 85  void Digitizer::DigitizeTOF(int np,float
85      c2_S[j] = c2_S[j]*tdcres[j];      c2_S[j] = c2_S[j]*tdcres[j];
86    }    }
87    /* **********************************  start loop over hits */    /* **********************************  start loop over hits */
88    if(Nthtof>*ntf)cout<<"NTHTOF > "<<*ntf<<" , event rejected ! "<<Nthtof<<endl;    if(Nthtof>*ntof)cout<<"NTHTOF > "<<*ntof<<" , event rejected ! "<<Nthtof<<endl;
89    else{    else{
90      for(Int_t nh=0; nh<Nthtof; nh++){      for(Int_t nh=0; nh<Nthtof; nh++){
 ////      if(Ipartof[nh]!=Ipa)continue;  
91        Float_t s_l_g[6] = {8.0, 8.0, 20.9, 22.0, 9.8, 8.3 };  // length of the lightguide        Float_t s_l_g[6] = {8.0, 8.0, 20.9, 22.0, 9.8, 8.3 };  // length of the lightguide
92        Float_t t1,t2,veff,veff1,veff0 ;        Float_t t1,t2,veff,veff1,veff0 ;
93        veff0 = 100.*1.0e8 ; // light velocity in the scintillator in m/sec        veff0 = 100.*1.0e8 ; // light velocity in the scintillator in m/sec
# Line 131  void Digitizer::DigitizeTOF(int np,float Line 101  void Digitizer::DigitizeTOF(int np,float
101        pmtleft=0;        pmtleft=0;
102        pmtright=0;        pmtright=0;
103  // WM: S12 paddles are "reversed" (Nov'07)  // WM: S12 paddles are "reversed" (Nov'07)
104        if (ip==2)        if (ip==2){
105          if (ipad==0)          if (ipad==0)
106            ipad=1;            ipad=1;
107          else          else
108            ipad=0;            ipad=0;
109        //    if (ip<6) {        }
110        if ((ip>-1)&&(ip<6)) {  //ToF paddles only, not S4        if ((ip>-1)&&(ip<6)) {  //ToF paddles only, not S4
111          Paddle2Pmt(ip, ipad, &pmtleft, &pmtright);          Paddle2Pmt(ip, ipad, &pmtleft, &pmtright);
112  // DC: evaluates mean position and path inside the paddle  // DC: evaluates mean position and path inside the paddle
# Line 166  void Digitizer::DigitizeTOF(int np,float Line 136  void Digitizer::DigitizeTOF(int np,float
136          for(Int_t j=0; j<2; j++){          for(Int_t j=0; j<2; j++){
137            QhitPad_pC[j]= Npho*FGeo[j]*effi*pmGain*echarge*1.E12*ScaleFact[pmtleft+j];            QhitPad_pC[j]= Npho*FGeo[j]*effi*pmGain*echarge*1.E12*ScaleFact[pmtleft+j];
138            // WM            // WM
139            knorm[j]=atte1[pmtleft+j]*exp(lambda1[pmtleft+j]*dimel[ip]/2.*pow(-1,j+1)) +            knorm[j]=atte1[pmtleft+j]*exp(lambda1[pmtleft+j]*dimel[ip]/2.*pow(-1,j+1)) + atte2[pmtleft+j]*exp(lambda2[pmtleft+j]*dimel[ip]/2.*pow(-1,j+1));
140              atte2[pmtleft+j]*exp(lambda2[pmtleft+j]*dimel[ip]/2.*pow(-1,j+1));            Atten[j]=atte1[pmtleft+j]*exp(tpos*lambda1[pmtleft+j]) + atte2[pmtleft+j]*exp(tpos*lambda2[pmtleft+j]) ;
           Atten[j]=atte1[pmtleft+j]*exp(tpos*lambda1[pmtleft+j]) +  
             atte2[pmtleft+j]*exp(tpos*lambda2[pmtleft+j]) ;  
141            QhitPmt_pC[j]= QhitPad_pC[j]*Atten[j]/knorm[j];            QhitPmt_pC[j]= QhitPad_pC[j]*Atten[j]/knorm[j];
142            if (DEBUG) {            if (DEBUG) {
143              cout<<"pmtleft "<<pmtleft<<" j "<<j<<endl;              cout<<"pmtleft "<<pmtleft<<" j "<<j<<endl;
# Line 180  void Digitizer::DigitizeTOF(int np,float Line 148  void Digitizer::DigitizeTOF(int np,float
148          if(DEBUG)cout<<"Npho "<<Npho<<" QhitPmt_pC "<<QhitPmt_pC[0]<<" "<<QhitPmt_pC[1]<<endl;            if(DEBUG)cout<<"Npho "<<Npho<<" QhitPmt_pC "<<QhitPmt_pC[0]<<" "<<QhitPmt_pC[1]<<endl;  
149          QevePmt_pC[pmtleft]  += QhitPmt_pC[0];          QevePmt_pC[pmtleft]  += QhitPmt_pC[0];
150          QevePmt_pC[pmtright] += QhitPmt_pC[1];          QevePmt_pC[pmtright] += QhitPmt_pC[1];
151  // TDC  //TDC
152  // WM right and left <->          // WM right and left <->
153          t1 = t1 + fabs(path[0]/veff) + s_l_g[ip]/veff1;          t1 = t1 + fabs(path[0]/veff) + s_l_g[ip]/veff1;
154          t2 = t2 + fabs(path[1]/veff) + s_l_g[ip]/veff1 ;  // Signal reaches PMT          t2 = t2 + fabs(path[1]/veff) + s_l_g[ip]/veff1 ;  // Signal reaches PMT
155          t1 = gRandom->Gaus(t1,dt1); //apply gaussian error  dt          t1 = gRandom->Gaus(t1,dt1); //apply gaussian error  dt
# Line 221  void Digitizer::DigitizeTOF(int np,float Line 189  void Digitizer::DigitizeTOF(int np,float
189    } // NTHTOF < 200    } // NTHTOF < 200
190  // ======  ADC ======  // ======  ADC ======
191    for(Int_t i=0; i<48; i++){    for(Int_t i=0; i<48; i++){
192      if (QevePmt_pC[i] < 800.)  ADCtof[i]= (Int_t)(ADC_pC0A + ADC_pC1A*QevePmt_pC[i] + ADC_pC2A*pow(QevePmt_pC[i],2) + ADC_pC3A*pow(QevePmt_pC[i],3));      if (QevePmt_pC[i] <= 800.)  ADCtof[i]= (Int_t)(ADC_pC0A + ADC_pC1A*QevePmt_pC[i] + ADC_pC2A*pow(QevePmt_pC[i],2) + ADC_pC3A*pow(QevePmt_pC[i],3));
193      if (QevePmt_pC[i] > 800.)  ADCtof[i]= (Int_t)(ADC_pC0B + ADC_pC1B*QevePmt_pC[i] + ADC_pC2B*pow(QevePmt_pC[i],2) + ADC_pC3B*pow(QevePmt_pC[i],3));      if (QevePmt_pC[i] > 800.)  ADCtof[i]= (Int_t)(ADC_pC0B + ADC_pC1B*QevePmt_pC[i] + ADC_pC2B*pow(QevePmt_pC[i],2) + ADC_pC3B*pow(QevePmt_pC[i],3));
194      if (QevePmt_pC[i] > 2485.) ADCtof[i]= (Int_t)(1758. + 0.54*QevePmt_pC[i]);  //assuming a fictional 0.54 ch/pC above ADCsat      if (QevePmt_pC[i] > 2485.) ADCtof[i]= (Int_t)(1758. + 0.54*QevePmt_pC[i]);  //assuming a fictional 0.54 ch/pC above ADCsat
195      if (ADCtof[i]>ADCsat) ADCtof[i]=ADCsat;      if (ADCtof[i]>ADCsat) ADCtof[i]=ADCsat;
# Line 235  void Digitizer::DigitizeTOF(int np,float Line 203  void Digitizer::DigitizeTOF(int np,float
203    
204  // ======  build  TDC coincidence  ======  // ======  build  TDC coincidence  ======
205    
206      //
207      for(Int_t i=0; i<48; i++) {
208        if((tdcpmt[i] - c1_S[i]) > 1e-7) {
209          tdcpmt[i] = 0.;
210          tdc[i] = 0.;
211        }
212      }// cycle to introduce a window for tdc
213    
214    
215    Float_t t_coinc = 0;    Float_t t_coinc = 0;
216    Int_t ilast = 100;    Int_t ilast = 100;
217    for (Int_t ii=0; ii<48;ii++)    for (Int_t ii=0; ii<48;ii++)
# Line 405  UChar_t Digitizer::EvaluateCrcTof(UChar_ Line 382  UChar_t Digitizer::EvaluateCrcTof(UChar_
382    UChar_t crcTof=0x00;    UChar_t crcTof=0x00;
383    UChar_t *pc=&crcTof, *pc2;    UChar_t *pc=&crcTof, *pc2;
384    pc2=pTof;    pc2=pTof;
385    for (Int_t jp=0; jp < 23; jp++){    for (Int_t jp=0; jp < 22; jp++){ // cecilia: fixed 23->22
386      //crcTof = crc8(...)      //crcTof = crc8(...)
387      Crc8Tof(pc2++,pc);      Crc8Tof(pc2++,pc);
388      //    printf("%2i --- %x\n",jp,crcTof);      //    printf("%2i --- %x\n",jp,crcTof);
# Line 474  void Digitizer::Crc8Tof(UChar_t *oldCRC, Line 451  void Digitizer::Crc8Tof(UChar_t *oldCRC,
451    //return r.word;    //return r.word;
452  };  };
453    
 //void Digitizer::Paddle2Pmt(Int_t plane, Int_t paddle, Int_t* &pmtleft, Int_t* &pmtright){  
454  void Digitizer::Paddle2Pmt(Int_t plane, Int_t paddle, Int_t *pl, Int_t *pr){  void Digitizer::Paddle2Pmt(Int_t plane, Int_t paddle, Int_t *pl, Int_t *pr){
455    //* @param plane    (0 - 5)    //* @param plane    (0 - 5)
456    //* @param paddle   (plane=0, paddle = 0,...5)    //* @param paddle   (plane=0, paddle = 0,...5)
# Line 485  void Digitizer::Paddle2Pmt(Int_t plane, Line 461  void Digitizer::Paddle2Pmt(Int_t plane,
461    //    //
462    Int_t somma=0;    Int_t somma=0;
463    Int_t np=plane;    Int_t np=plane;
464    for(Int_t j=0; j<np; j++)    for(Int_t j=0; j<np; j++)somma+=pads[j];
     somma+=pads[j];  
465    padid=paddle+somma;    padid=paddle+somma;
466    *pl = padid*2;    *pl = padid*2;
467    //  *pr = *pr + 1;    //  *pr = *pr + 1;
# Line 498  void Digitizer::LoadTOFCalib(int np,floa Line 473  void Digitizer::LoadTOFCalib(int np,floa
473    Int_t error = 0,temp=0;    Int_t error = 0,temp=0;
474    GL_PARAM *glparam = new GL_PARAM();    GL_PARAM *glparam = new GL_PARAM();
475    error = glparam->Query_GL_PARAM(3,202,fDbc);    error = glparam->Query_GL_PARAM(3,202,fDbc);
476    calfile.str("");    if(!error){
477    calfile << glparam->PATH.Data() << "/";      calfile.str("");
478    calfile << glparam->NAME.Data();      calfile << glparam->PATH.Data() << "/";
479    printf("\n Using TOF calibration file: \n %s\n",calfile.str().c_str());      calfile << glparam->NAME.Data();
480    ifstream fileTriggerCalib;      printf("\n Using TOF calibration file: \n %s\n",calfile.str().c_str());
481    fileTriggerCalib.open(calfile.str().c_str());      ifstream fileTriggerCalib;
482    if(!fileTriggerCalib)printf("debug: no trigger calib file!\n");      fileTriggerCalib.open(calfile.str().c_str());
483    // correct readout WM Oct '07      for(Int_t i=0; i<np; i++){
484    for(Int_t i=0; i<np; i++){        fileTriggerCalib >> temp;
485      fileTriggerCalib >> temp;        fileTriggerCalib >> atte1[i];
486      fileTriggerCalib >> atte1[i];        fileTriggerCalib >> lambda1[i];
487      fileTriggerCalib >> lambda1[i];        fileTriggerCalib >> atte2[i];
488      fileTriggerCalib >> atte2[i];        fileTriggerCalib >> lambda2[i];
489      fileTriggerCalib >> lambda2[i];        fileTriggerCalib >> temp;
490      fileTriggerCalib >> temp;      }
491        fileTriggerCalib.close();
492      }
493      else{
494        cout<<endl<<"                        *********** ATTENTION ***********"<<endl;
495        cout<<endl<<"                           TOF: NO trigger calib file!"<<endl<<endl;
496        cout<<endl<<"                        TOF digitized data will be wrong!"<<endl<<endl;
497        for(Int_t i=0; i<np; i++){
498          atte1[i]=0.;
499          lambda1[i]=0.;
500          atte2[i]=0.;
501          lambda2[i]=0.;
502        }
503    }    }
   fileTriggerCalib.close();  
504    //end tof calib    //end tof calib
505  }  }
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.7
  ViewVC Help
Powered by ViewVC 1.1.23