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

Diff of /PamelaDigitizer/Digitizer.cxx

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

revision 1.4 by orsi, Wed Oct 31 18:17:58 2007 UTC revision 1.5 by silvio, Wed Nov 28 18:54:31 2007 UTC
# Line 1074  Int_t Digitizer::DigitizeTOF() { Line 1074  Int_t Digitizer::DigitizeTOF() {
1074    Float_t pmGain = 3.5e6;  /* PMT Gain: the same for all PMTs */    Float_t pmGain = 3.5e6;  /* PMT Gain: the same for all PMTs */
1075    Float_t effi=0.21;       /* Efficienza di fotocatodo */    Float_t effi=0.21;       /* Efficienza di fotocatodo */
1076    
1077    Float_t ADC_pC0=-58.1;      //  ADC/pC conversion coefficient 0    //   Float_t ADC_pC0=-58.1;      //  ADC/pC conversion coefficient 0
1078    Float_t ADC_pC1=1.728;      //  ADC/pC conversion coefficient 1    //   Float_t ADC_pC1=1.728;      //  ADC/pC conversion coefficient 1
1079    Float_t ADC_pC2=-4.063e-05; //  ADC/pC conversion coefficient 2    //   Float_t ADC_pC2=-4.063e-05; //  ADC/pC conversion coefficient 2
1080    Float_t ADC_pC3=-5.763e-08; //  ADC/pC conversion coefficient 3    //   Float_t ADC_pC3=-5.763e-08; //  ADC/pC conversion coefficient 3
1081      
1082      // pC < 800
1083      Float_t ADC_pC0A =   -4.437616e+01   ;
1084      Float_t ADC_pC1A =   1.573329e+00    ;
1085      Float_t ADC_pC2A =   2.780518e-04  ;
1086      Float_t ADC_pC3A =  -2.302160e-07 ;
1087      
1088      // pC   > 800:
1089      Float_t ADC_pC0B =    -2.245756e+02  ;
1090      Float_t ADC_pC1B =     2.184156e+00  ;
1091      Float_t ADC_pC2B =     -4.171825e-04  ;
1092      Float_t ADC_pC3B =     3.789715e-08  ;
1093    
1094    Float_t pCthres=40.;     // threshold in charge    Float_t pCthres=40.;     // threshold in charge
1095    Int_t ADClast=4095;      // no signal --> ADC ch=4095    Int_t ADClast=4095;      // no signal --> ADC ch=4095
# Line 1085  Int_t Digitizer::DigitizeTOF() { Line 1097  Int_t Digitizer::DigitizeTOF() {
1097    Int_t ADCtof[48];    Int_t ADCtof[48];
1098    
1099    
1100  // ---- introduce scale factors to tune simul ADC to real data   24-oct DC    // ---- introduce scale factors to tune simul ADC to real data   24-oct DC
1101    Float_t ScaleFact[48]={0.18,0.22,0.35,0.26,0.47,0.35,0.31,0.37,    //   Float_t ScaleFact[48]={0.18,0.22,0.35,0.26,0.47,0.35,0.31,0.37,
1102                           0.44,0.23,0.38,0.60,0.39,0.29,0.40,0.23,    //                     0.44,0.23,0.38,0.60,0.39,0.29,0.40,0.23,
1103                           0.30,0.66,0.22,1.53,0.17,0.55,    //                     0.30,0.66,0.22,1.53,0.17,0.55,
1104                           0.84,0.19,0.21,1.64,0.62,0.13,    //                     0.84,0.19,0.21,1.64,0.62,0.13,
1105                           0.18,0.15,0.10,0.14,0.14,0.14,0.14,0.12,    //                     0.18,0.15,0.10,0.14,0.14,0.14,0.14,0.12,
1106                           0.26,0.18,0.25,0.23,0.20,0.40,    //                     0.26,0.18,0.25,0.23,0.20,0.40,0.19,0.23,0.25,0.23,0.25,0.20};
1107                           0.19,0.23,0.25,0.23,0.25,0.20};    
1108        // new scale factors: WM 30-Oct-07
1109      //   Float_t ScaleFact[48]={0.35,0.41,0.32,0.34,0.58,0.47,0.42,0.44,
1110      //                     0.50,0.34,0.50,0.50,0.51,0.42,0.46,0.25,
1111      //                     0.20,0.38,0.29,0.49,0.24,0.68,
1112      //                     0.30,0.26,0.28,0.79,0.31,0.12,
1113      //                     0.25,0.21,0.14,0.20,
1114      //                     0.16,0.17,0.19,0.18,
1115      //                     0.34,0.27,0.34,0.31,0.25,0.57,
1116      //                     0.24,0.34,0.34,0.32,0.31,0.30};
1117      
1118      Float_t ScaleFact[48]={0.39, 0.49, 0.38, 0.40, 0.65, 0.51, 0.43,
1119                             0.49, 0.58, 0.38, 0.53, 0.57, 0.53, 0.45, 0.49, 0.16,
1120                             0.15, 0.44, 0.28, 0.57, 0.26, 0.72, 0.37, 0.29, 0.30, 0.89,
1121                             0.37, 0.08,  0.27, 0.23, 0.12, 0.22, 0.15, 0.16, 0.21,
1122                             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 };
1123    
1124    for(Int_t i=0; i<48; i++){    for(Int_t i=0; i<48; i++){
1125      QevePmt_pC[i] = 0;      QevePmt_pC[i] = 0;
1126      ADCtof[i]=0;      ADCtof[i]=0;
# Line 1140  Int_t Digitizer::DigitizeTOF() { Line 1167  Int_t Digitizer::DigitizeTOF() {
1167    Float_t thresh=10.; // to be defined better... (Wolfgang)    Float_t thresh=10.; // to be defined better... (Wolfgang)
1168    
1169    // === TDC: simulate timing for each paddle    // === TDC: simulate timing for each paddle
1170    Float_t dt1 = 285.e-12 ;   // single PMT resolution    //Float_t dt1 = 285.e-12 ;   // single PMT resolution
1171    //    Float_t dt1 = 10.e-12 ;   // TEST    Float_t dt1 = 425.e-12 ;   // single PMT resolution (WM, Nov'07)
1172    Float_t tdcres[50],c1_S[50],c2_S[50],c3_S[50];    Float_t tdcres[50],c1_S[50],c2_S[50],c3_S[50];
1173    for(Int_t j=0;j<48;j++)  tdcres[j] = 50.E-12;   // TDC resolution 50 picosec    for(Int_t j=0;j<48;j++)  tdcres[j] = 50.E-12;   // TDC resolution 50 picosec
1174    for(Int_t j=0;j<48;j++)  c1_S[j] = 500.;  // cable length in channels    for(Int_t j=0;j<48;j++)  c1_S[j] = 500.;  // cable length in channels
# Line 1169  Int_t Digitizer::DigitizeTOF() { Line 1196  Int_t Digitizer::DigitizeTOF() {
1196      ipad = Ipaddle[nh]-1;      ipad = Ipaddle[nh]-1;
1197      pmtleft=0;      pmtleft=0;
1198      pmtright=0;      pmtright=0;
1199      
1200      if (ip<6) {      // WM: S12 paddles are "reversed" (Nov'07)
1201        if (ip==2)
1202          if (ipad==0)
1203            ipad=1;
1204          else
1205            ipad=0;
1206    
1207       if (ip<6) {
1208        Paddle2Pmt(ip, ipad, &pmtleft, &pmtright);        Paddle2Pmt(ip, ipad, &pmtleft, &pmtright);
1209    
1210        // DC: evaluates mean position and path inside the paddle        // DC: evaluates mean position and path inside the paddle
# Line 1181  Int_t Digitizer::DigitizeTOF() { Line 1215  Int_t Digitizer::DigitizeTOF() {
1215        if(ip==0 || ip==3 || ip==4)        if(ip==0 || ip==3 || ip==4)
1216          tpos = (Yintof[nh]+Youttof[nh])/2.;          tpos = (Yintof[nh]+Youttof[nh])/2.;
1217        else        else
1218          if(ip==1 || ip==2 || ip==5)   //--- Strip in X per S12,S21,S32          if(ip==1 || ip==2 || ip==5)   //--- Strip in X for S12,S21,S32
1219            tpos = (Xintof[nh]+Xouttof[nh])/2.;            tpos = (Xintof[nh]+Xouttof[nh])/2.;
1220          else //if (ip!=6)          else //if (ip!=6)
1221            printf("*** WARNING TOF: this option should never occur! (ip=%2i, nh=%2i)\n",ip,nh);            printf("*** WARNING TOF: this option should never occur! (ip=%2i, nh=%2i)\n",ip,nh);
# Line 1196  Int_t Digitizer::DigitizeTOF() { Line 1230  Int_t Digitizer::DigitizeTOF() {
1230          cout <<"pmtleft, pmtright "<<pmtleft<<" "<<pmtright<<endl;          cout <<"pmtleft, pmtright "<<pmtleft<<" "<<pmtright<<endl;
1231        }        }
1232                
1233        // constant geometric factor, for the moment        // constant geometric factor, the rest is handled by the scaling factor
1234        FGeo[0] =0.5;        FGeo[0] =0.5;
1235        FGeo[1] =0.5;        FGeo[1] =0.5;
1236        //  FGeo[1] = atan(path[1]/dimes[ip])/6.28318; // fraction of photons toward SX        //  FGeo[1] = atan(path[1]/dimes[ip])/6.28318; // fraction of photons toward left
1237        //  FGeo[2] = atan(path[2]/dimes[ip])/6.28318; // toward DX        //  FGeo[2] = atan(path[2]/dimes[ip])/6.28318; // toward right
1238                
1239                
1240        //  Npho = Poisson(ERELTOF[nh])*Pho_keV*1e6  Poissonian fluctuations to be inserted-DC        //  Npho = Poisson(ERELTOF[nh])*Pho_keV*1e6  Poissonian fluctuations to be inserted-DC
# Line 1209  Int_t Digitizer::DigitizeTOF() { Line 1243  Int_t Digitizer::DigitizeTOF() {
1243        Float_t knorm[2]={0., 0.}; // Donatella        Float_t knorm[2]={0., 0.}; // Donatella
1244        Float_t Atten[2]={0., 0.}; // Donatella        Float_t Atten[2]={0., 0.}; // Donatella
1245        for(Int_t j=0; j<2; j++){        for(Int_t j=0; j<2; j++){
1246          QhitPad_pC[j]= Npho*FGeo[j]*effi*pmGain*echarge*1.E12; // corrected WM          QhitPad_pC[j]= Npho*FGeo[j]*effi*pmGain*echarge*1.E12*ScaleFact[pmtleft+j];
1247          // WM          // WM
1248          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)) +
1249            atte2[pmtleft+j]*exp(lambda2[pmtleft+j]*dimel[ip]/2.*pow(-1,j+1));            atte2[pmtleft+j]*exp(lambda2[pmtleft+j]*dimel[ip]/2.*pow(-1,j+1));
# Line 1233  Int_t Digitizer::DigitizeTOF() { Line 1267  Int_t Digitizer::DigitizeTOF() {
1267        QevePmt_pC[pmtright] += QhitPmt_pC[1];        QevePmt_pC[pmtright] += QhitPmt_pC[1];
1268                
1269        // TDC        // TDC
1270  // WM right and left <->        // WM right and left <->
1271  //      t2 = t2 + fabs(path[0]/veff) + s_l_g[ip]/veff1 ;  // Signal reaches PMT        //      t2 = t2 + fabs(path[0]/veff) + s_l_g[ip]/veff1 ;  // Signal reaches PMT
1272  //      t1 = t1 + fabs(path[1]/veff) + s_l_g[ip]/veff1;        //      t1 = t1 + fabs(path[1]/veff) + s_l_g[ip]/veff1;
1273    
1274        t1 = t1 + fabs(path[0]/veff) + s_l_g[ip]/veff1;        t1 = t1 + fabs(path[0]/veff) + s_l_g[ip]/veff1;
1275        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
1276                
1277        Float_t t1save = t1;        t1 = gRandom->Gaus(t1,dt1); //apply gaussian error  dt
1278        Float_t t2save = t2;        t2 = gRandom->Gaus(t2,dt1); //apply gaussian error  dt
   
 /*  
       TRandom r;  
 // This does not work... WM -  but works in my simulation code ??  
 //      t1 = r.Gaus(t1,dt1);  //apply gaussian error  dt  
 //      t2 = r.Gaus(t2,dt1);  //apply gaussian error  dt  
 */        
         t1 = gRandom->Gaus(t1,dt1); //apply gaussian error  dt  
         t2 = gRandom->Gaus(t2,dt1); //apply gaussian error  dt  
   
 //      cout<<1E12*(t1save-t1)<<" "<<1E12*(t2save-t2)<<endl;  
1279    
1280        t1 = t1 + c1_S[pmtleft] ;  // Signal reaches Discriminator ,TDC starts  to run        t1 = t1 + c1_S[pmtleft] ;  // Signal reaches Discriminator ,TDC starts  to run
1281        t2 = t2 + c1_S[pmtright] ;        t2 = t2 + c1_S[pmtright] ;
# Line 1295  Int_t Digitizer::DigitizeTOF() { Line 1318  Int_t Digitizer::DigitizeTOF() {
1318        
1319    // ======  ADC ======    // ======  ADC ======
1320    
   
1321    for(Int_t i=0; i<48; i++){    for(Int_t i=0; i<48; i++){
1322      if(QevePmt_pC[i] >= pCthres){      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));
1323        ADCtof[i]= (Int_t)(ADC_pC0 + ADC_pC1*QevePmt_pC[i] + ADC_pC2*pow(QevePmt_pC[i],2) + ADC_pC3*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));
1324    } else      if (QevePmt_pC[i] > 2485.) ADCtof[i]= (Int_t)(1758. + 0.54*QevePmt_pC[i]);  //assuming a fictional 0.54 ch/pC above ADCsat
1325        ADCtof[i]= ADClast;        if (ADCtof[i]>ADCsat) ADCtof[i]=ADCsat;
1326  }      if (QevePmt_pC[i] < pCthres)  ADCtof[i]= ADClast;
1327          if (ADCtof[i] < 0) ADCtof[i]=ADClast;    
1328  // ---- introduce scale factors to tune simul ADC to real data   24-oct DC      if (ADCtof[i] > ADClast) ADCtof[i]=ADClast;
   
   for(Int_t i=0; i<48; i++){  
     if(ADCtof[i] != ADClast){  
                              //      printf("%3d, %4d, %4.2f\n",i, ADCtof[i],ScaleFact[i]);  
       ADCtof[i]= Int_t (ADCtof[i]*ScaleFact[i]);  
                              //      printf("%3d, %4d,\n",i, ADCtof[i]);  
     }  
1329    }    }
1330    
1331    for(Int_t i=0; i<48; i++){  //   for(Int_t i=0; i<48; i++){
1332      if(ADCtof[i] != ADClast){  //     if(QevePmt_pC[i] >= pCthres){
1333        if(ADCtof[i]> ADCsat) ADCtof[i]=ADCsat;  //       ADCtof[i]= (Int_t)(ADC_pC0 + ADC_pC1*QevePmt_pC[i] + ADC_pC2*pow(QevePmt_pC[i],2) + ADC_pC3*pow(QevePmt_pC[i],3));
1334        else if(ADCtof[i]< 0) ADCtof[i]=ADClast;      //   } else
1335    }  //       ADCtof[i]= ADClast;  
1336    }      // }
1337    // ======  build  TDC coincidence  ======    
1338    // // ---- introduce scale factors to tune simul ADC to real data   24-oct DC
1339    
1340    //   for(Int_t i=0; i<48; i++){
1341    //     if(ADCtof[i] != ADClast){
1342    //                              //      printf("%3d, %4d, %4.2f\n",i, ADCtof[i],ScaleFact[i]);
1343    //       ADCtof[i]= Int_t (ADCtof[i]*ScaleFact[i]);
1344    //                              //      printf("%3d, %4d,\n",i, ADCtof[i]);
1345    //     }
1346    //   }
1347    
1348    //   for(Int_t i=0; i<48; i++){
1349    //     if(ADCtof[i] != ADClast){
1350    //       if(ADCtof[i]> ADCsat) ADCtof[i]=ADCsat;
1351    //       else if(ADCtof[i]< 0) ADCtof[i]=ADClast;    
1352    //   }
1353    //   }    
1354    
1355    
1356    // ======  build  TDC coincidence  ======
1357    
1358    Float_t t_coinc = 0;    Float_t t_coinc = 0;
1359    Int_t ilast = 100;    Int_t ilast = 100;
# Line 1337  Int_t Digitizer::DigitizeTOF() { Line 1371  Int_t Digitizer::DigitizeTOF() {
1371      if (tdc[ii] != 0) tdc1[ii] = t_coinc - tdcpmt[ii];  // test 2      if (tdc[ii] != 0) tdc1[ii] = t_coinc - tdcpmt[ii];  // test 2
1372      tdc1[ii] = tdc1[ii]/tdcres[ii];                     // divide by TDC resolution      tdc1[ii] = tdc1[ii]/tdcres[ii];                     // divide by TDC resolution
1373      if (tdc[ii] != 0) tdc1[ii] = tdc1[ii] + c3_S[ii];  // add cable length c3      if (tdc[ii] != 0) tdc1[ii] = tdc1[ii] + c3_S[ii];  // add cable length c3
   
1374    } // missing parenthesis inserted! (Silvio)    } // missing parenthesis inserted! (Silvio)
1375    
1376    for(Int_t i=0; i<48; i++){    for(Int_t i=0; i<48; i++){
# Line 1358  Int_t Digitizer::DigitizeTOF() { Line 1391  Int_t Digitizer::DigitizeTOF() {
1391    
1392  //------ use channelmap  18-oct WM  //------ use channelmap  18-oct WM
1393    
1394  Int_t  channelmap[] =  {3,21,11,29,19,45,27,37,36,28,44,20,5,12,13,4,    Int_t  channelmap[] =  {3,21,11,29,19,45,27,37,36,28,44,20,5,12,13,4,
1395                          6,47,14,39,22,31,30,23,38,15,46,7,0,33,16,24,                          6,47,14,39,22,31,30,23,38,15,46,7,0,33,16,24,
1396                          8,41,32,40,25,17,34,9,42,1,2,10,18,26,35,43};                          8,41,32,40,25,17,34,9,42,1,2,10,18,26,35,43};
1397    
# Line 1366  Int_t  channelmap[] =  {3,21,11,29,19,45 Line 1399  Int_t  channelmap[] =  {3,21,11,29,19,45
1399    Int_t TDChelp[48];    Int_t TDChelp[48];
1400    
1401    for(Int_t i=0; i<48; i++){    for(Int_t i=0; i<48; i++){
1402        Int_t ii=channelmap[i];      Int_t ii=channelmap[i];
1403        ADChelp[ii]= ADCtof[i];      ADChelp[ii]= ADCtof[i];
1404        TDChelp[ii]= TDCint[i];      TDChelp[ii]= TDCint[i];
1405                            }    }
   
   for(Int_t i=0; i<48; i++){  
       ADCtof[i]= ADChelp[i];  
       TDCint[i]= TDChelp[i];  
                           }  
   
   
 /*  
 //--- fake data ------------------------  
   for(Int_t i=0; i<48; i++){  
       ADCtof[i]= 100 + 10*i;  
       TDCint[i]= 800 + 10*i;  
 //      cout<<i<<" "<<ADCtof[i]<<" "<<TDCint[i]<<endl;  
                           }  
 */  
1406    
 /*  
1407    for(Int_t i=0; i<48; i++){    for(Int_t i=0; i<48; i++){
1408     if (((ADCtof[i]>0)&&(ADCtof[i]<4095)) || ((TDCint[i]>0)&&(TDCint[i]<4095)))  cout<<i<<" "<<ADCtof[i]<<" "<<TDCint[i]<<endl;      ADCtof[i]= ADChelp[i];
1409                            }      TDCint[i]= TDChelp[i];
1410  */    }
1411    
1412    
1413  // ======  write fDataTof  =======  // ======  write fDataTof  =======
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.5
  ViewVC Help
Powered by ViewVC 1.1.23