/[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.1 by silvio, Thu Sep 13 11:00:53 2007 UTC revision 1.3 by orsi, Thu Oct 11 11:29:21 2007 UTC
# Line 1  Line 1 
1    //               ------ PAMELA Digitizer ------
2    //
3    // Date, release and how-to: see file Pamelagp2Digits.cxx
4    //
5    // NB: Check length physics packet [packet type (0x10 = physics data)]
6  //  //
7  #include <sstream>  #include <sstream>
8  #include <fstream>  #include <fstream>
9  #include <stdlib.h>  #include <stdlib.h>
10    #include <stdio.h>
11  #include <string.h>  #include <string.h>
12  #include <ctype.h>  #include <ctype.h>
13    #include <time.h>
14  #include "Riostream.h"  #include "Riostream.h"
15  #include "TFile.h"  #include "TFile.h"
16  #include "TDirectory.h"  #include "TDirectory.h"
# Line 253  void Digitizer::Loop() { Line 260  void Digitizer::Loop() {
260      DigitizeAC();      DigitizeAC();
261      DigitizeCALO();      DigitizeCALO();
262      DigitizeTrack();      DigitizeTrack();
263      //DigitizeS4();      DigitizeS4();
264      DigitizeND();      DigitizeND();
265      //        //
     // Create CPU header, we need packet type (0x10 = physics data) and packet length.  
     //  
     UInt_t length = (fCALOlength + fACbuffer + fTracklength)*2;  
     DigitizePSCU(length,0x10);  
     //  
266      // Add padding to 64 bits      // Add padding to 64 bits
267      //      //
268      AddPadding();      AddPadding();
269    //
270        // Create CPU header, we need packet type (0x10 = physics data) and packet length.
271      //      //
272        UInt_t length=2*(fCALOlength+fACbuffer+fTracklength+fNDbuffer+fS4buffer)+fPadding+fTOFbuffer+fTRIGGERbuffer;
273        //UInt_t length=2*(fCALOlength+fACbuffer+fTracklength+fNDbuffer)+fPadding+fTOFbuffer+fTRIGGERbuffer;
274        DigitizePSCU(length,0x10);
275      if ( !i%100 ) std::cout << "writing event " << i << endl;      if ( !i%100 ) std::cout << "writing event " << i << endl;
276      WriteData();      WriteData();
277    };    };
# Line 1030  void Digitizer::DigitizeTRIGGER() { Line 1037  void Digitizer::DigitizeTRIGGER() {
1037  Int_t Digitizer::DigitizeTOF() {  Int_t Digitizer::DigitizeTOF() {
1038    //fDataTof: 12 x 23 bytes (=276 bytes)    //fDataTof: 12 x 23 bytes (=276 bytes)
1039    UChar_t *pTof=fDataTof;    UChar_t *pTof=fDataTof;
1040      Bool_t DEBUG=false;
1041    
1042    // --- activate branches:    // --- activate branches:
1043    fhBookTree->SetBranchStatus("Nthtof",1);    fhBookTree->SetBranchStatus("Nthtof",1);
# Line 1046  Int_t Digitizer::DigitizeTOF() { Line 1054  Int_t Digitizer::DigitizeTOF() {
1054    // ------ evaluate energy in each pmt: ------    // ------ evaluate energy in each pmt: ------
1055    // strip geometry (lenght/width)    // strip geometry (lenght/width)
1056    Float_t dimel[6] = {33.0, 40.8 ,18.0, 15.0, 15.0, 18.0};    Float_t dimel[6] = {33.0, 40.8 ,18.0, 15.0, 15.0, 18.0};
1057    Float_t dimes[6] = {5.1, 5.5, 7.5, 9.0, 6.0, 5.0};    //Float_t dimes[6] = {5.1, 5.5, 7.5, 9.0, 6.0, 5.0};
1058        
1059    //  S11 8 paddles  33.0 x 5.1 cm    //  S11 8 paddles  33.0 x 5.1 cm
1060    //  S12 6 paddles  40.8 x 5.5 cm    //  S12 6 paddles  40.8 x 5.5 cm
# Line 1089  Int_t Digitizer::DigitizeTOF() { Line 1097  Int_t Digitizer::DigitizeTOF() {
1097    };    };
1098    Float_t atte1[48],atte2[48],lambda1[48],lambda2[48];    Float_t atte1[48],atte2[48],lambda1[48],lambda2[48];
1099    Int_t temp=0;    Int_t temp=0;
1100      // correct readout WM Oct '07
1101    for(Int_t i=0; i<48; i++){    for(Int_t i=0; i<48; i++){
1102      fileTriggerCalib >> temp;      fileTriggerCalib >> temp;
1103      fileTriggerCalib >> atte1[i];      fileTriggerCalib >> atte1[i];
     fileTriggerCalib >> atte2[i];  
1104      fileTriggerCalib >> lambda1[i];      fileTriggerCalib >> lambda1[i];
1105        fileTriggerCalib >> atte2[i];
1106      fileTriggerCalib >> lambda2[i];      fileTriggerCalib >> lambda2[i];
1107      fileTriggerCalib >> temp;      fileTriggerCalib >> temp;
1108    }    }
# Line 1106  Int_t Digitizer::DigitizeTOF() { Line 1115  Int_t Digitizer::DigitizeTOF() {
1115    //    fine lettura dal file */    //    fine lettura dal file */
1116    
1117    //const Int_t nmax=??; = Nthtof    //const Int_t nmax=??; = Nthtof
1118    Int_t nh, ip, ipad, ipmt;    Int_t ip, ipad;
1119      //Int_t ipmt;
1120    Int_t pmtleft=0, pmtright=0;    Int_t pmtleft=0, pmtright=0;
1121    Int_t *pl, *pr;    Int_t *pl, *pr;
1122    pl = &pmtleft;    pl = &pmtleft;
1123    pr = &pmtright;    pr = &pmtright;
1124    
1125    /* **********************************  inizio loop sugli hit */    // TDC variables:
1126      Int_t TDClast=4095;      // no signal --> ADC ch=4095
1127      Int_t TDCint[48];
1128      Float_t  tdc[48],tdc1[48],tdcpmt[48];
1129      for(Int_t i=0; i<48; i++)
1130        tdcpmt[i] = 1000.;
1131      Float_t thresh=10.; // to be defined better... (Wolfgang)
1132    
1133        // === TDC: simulate timing for each paddle
1134        Float_t dt1 = 285.e-12 ;   // single PMT resolution
1135        Float_t tdcres[50],c1_S[50],c2_S[50],c3_S[50];
1136        for(Int_t j=0;j<48;j++)  tdcres[j] = 50.E-12;   // TDC resolution 50 picosec
1137        for(Int_t j=0;j<48;j++)  c1_S[j] = 500.;  // cable length in channels
1138        for(Int_t j=0;j<48;j++)  c2_S[j] = 0.;
1139        for(Int_t j=0;j<48;j++)  c3_S[j] = 1000.;
1140        for(Int_t j=0;j<48;j++)  c1_S[j] = c1_S[j]*tdcres[j];   // cable length in sec
1141        for(Int_t j=0;j<48;j++)  c2_S[j] = c2_S[j]*tdcres[j];
1142         // ih = 0 + i1;  // not used?? (Silvio)
1143    
1144      /* **********************************  start loop over hits */
1145        
1146    for(Int_t nh=0; nh<Nthtof; nh++){    for(Int_t nh=0; nh<Nthtof; nh++){
1147            
# Line 1122  Int_t Digitizer::DigitizeTOF() { Line 1151  Int_t Digitizer::DigitizeTOF() {
1151        FGeo[j]=0.;        FGeo[j]=0.;
1152      }      }
1153    
1154      // ridefiniz. piano e pad per i vettori in C      Float_t s_l_g[6] = {8.0, 8.0, 20.9, 22.0, 9.8, 8.3 };  // length of the lightguide
1155        Float_t t1,t2,veff,veff1,veff0 ;
1156        veff0 = 100.*1.0e8 ; // light velocity in the scintillator in m/sec
1157        veff1 = 100.*1.5e8; // light velocity in the lightguide in m/sec
1158        veff=veff0;         // signal velocity in the paddle
1159    
1160        t1 = Timetof[nh] ;  // Start
1161        t2 = Timetof[nh] ;
1162    
1163        // Donatella: redefinition plane and pad for vectors in C
1164      ip = Ipltof[nh]-1;      ip = Ipltof[nh]-1;
1165      ipad = Ipaddle[nh]-1;      ipad = Ipaddle[nh]-1;
1166      pmtleft=0;      pmtleft=0;
1167      pmtright=0;      pmtright=0;
1168        
1169      //Paddle2Pmt((Int_t)ip, (Int_t) ipad, (Int_t*) &pmtleft, (Int_t*) &pmtright);      if (ip<6) {
1170      Paddle2Pmt(ip, ipad, &pmtleft, &pmtright);        Paddle2Pmt(ip, ipad, &pmtleft, &pmtright);
     //Paddle2Pmt(ip, ipad, pl, pr);  
1171            
1172      // per avere anche la corrispondenza pmt --> half board e canale        // per avere anche la corrispondenza pmt --> half board e canale
1173      // metodo GetPMTIndex(Int_t ipmt, Int_t &hb, Int_t &ch) // non lo usiamo x ora        // metodo GetPMTIndex(Int_t ipmt, Int_t &hb, Int_t &ch) // non lo usiamo x ora
1174        
1175      /*calcola la pos media e il path all'interno della paddle */        // evaluates mean position and path inside the paddle
1176        
1177      Float_t tpos=0.;        Float_t tpos=0.;
1178      Float_t path[2] = {0., 0.};        Float_t path[2] = {0., 0.};
1179      //--- Strip in Y = S11,S22,S31 ------        //--- Strip in Y = S11,S22,S31 ------
1180      if(ip==0 || ip==3 || ip==4)        if(ip==0 || ip==3 || ip==4)
1181        tpos = (Yintof[nh]+Youttof[nh])/2.;          tpos = (Yintof[nh]+Youttof[nh])/2.;
1182      else        else
1183        if(ip==1 || ip==2 || ip==5)   //--- Strip in X per S12,S21,S32          if(ip==1 || ip==2 || ip==5)   //--- Strip in X per S12,S21,S32
1184          tpos = (Xintof[nh]+Xouttof[nh])/2.;            tpos = (Xintof[nh]+Xouttof[nh])/2.;
1185        else if (ip!=6)          else //if (ip!=6)
1186          printf("*** Warning: 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);
1187      path[0]= tpos + dimel[ip]/2.;        path[0]= tpos + dimel[ip]/2.;
1188      path[1]= dimel[ip]/2.- tpos;        path[1]= dimel[ip]/2.- tpos;
1189    
1190      //  cout <<"Strip N. ="<< ipaddle <<" piano n.= "<< iplane <<" POSIZ = "<< tpos <<"\n";        //  cout <<"Strip N. ="<< ipaddle <<" piano n.= "<< iplane <<" POSIZ = "<< tpos <<"\n";
1191    
1192      /* per il momento metto un fattore geometrico costante*/          if (DEBUG) {
1193      FGeo[0] =0.5;          cout <<" plane "<<ip<<" strip # ="<< ipad <<" tpos  "<< tpos <<"\n";
1194      FGeo[1] =0.5;          cout <<"pmtleft, pmtright "<<pmtleft<<" "<<pmtright<<endl;
1195      //  FGeo[1] = atan(path[1]/dimes[ip])/6.28318; // frazione fotoni verso SX        }
1196      //  FGeo[2] = atan(path[2]/dimes[ip])/6.28318; // e verso DX  
1197          // constant geometric factor, for the moment
1198          FGeo[0] =0.5;
1199          FGeo[1] =0.5;
1200          //  FGeo[1] = atan(path[1]/dimes[ip])/6.28318; // frazione fotoni verso SX
1201          //  FGeo[2] = atan(path[2]/dimes[ip])/6.28318; // e verso DX
1202                
1203      /*  rimando la fluttuazione poissoniana sui fotoni prodotti        /*  rimando la fluttuazione poissoniana sui fotoni prodotti
1204          sto studiando come funziona la funzione:            sto studiando come funziona la funzione:
1205          long int i = sto.Poisson(double x);  */            long int i = sto.Poisson(double x);  */
1206      //  Npho = Poisson(ERELTOF[nh])*Pho_keV*1e6   Eloss in GeV ?        //  Npho = Poisson(ERELTOF[nh])*Pho_keV*1e6   Eloss in GeV ?
1207      Npho = Ereltof[nh]*Pho_keV*10.0e6;  // Eloss in GeV ?        Npho = Ereltof[nh]*Pho_keV*1.0e6;  // Eloss in GeV ?
1208        
1209      Float_t knorm[2]={0., 0.}; // Donatella        Float_t knorm[2]={0., 0.}; // Donatella
1210      Float_t Atten[2]={0., 0.}; // Donatella        Float_t Atten[2]={0., 0.}; // Donatella
1211      for(Int_t j=0; j<2; j++){        for(Int_t j=0; j<2; j++){
1212        QhitPad_pC[j]= Npho*FGeo[j]*effi*pmGain*echarge;          QhitPad_pC[j]= Npho*FGeo[j]*effi*pmGain*echarge*1.E12; // corrected WM
1213        knorm[j]=QhitPad_pC[j]/(atte1[pmtleft+j]*exp((dimel[ip]/2.*pow(-1,j+1))/lambda1[pmtleft+j]) +          /*      knorm[j]=QhitPad_pC[j]/(atte1[pmtleft+j]*exp((dimel[ip]/2.*pow(-1,j+1))/lambda1[pmtleft+j]) +
1214                                atte2[pmtleft+j]*exp((dimel[ip]/2.*pow(-1,j+1))/lambda2[pmtleft+j]));                                  atte2[pmtleft+j]*exp((dimel[ip]/2.*pow(-1,j+1))/lambda2[pmtleft+j]));
1215            Atten[j]=knorm[j]*(atte1[pmtleft+j]*exp(tpos/lambda1[pmtleft+j]) +
1216                               atte2[pmtleft+j]*exp(tpos/lambda2[pmtleft+j]));
1217            QhitPmt_pC[j]= QhitPad_pC[j]*Atten[j];
1218            */
1219            // WM
1220            knorm[j]=atte1[pmtleft+j]*exp(lambda1[pmtleft+j]*dimel[ip]/2.*pow(-1,j+1)) +
1221              atte2[pmtleft+j]*exp(lambda2[pmtleft+j]*dimel[ip]/2.*pow(-1,j+1));
1222            Atten[j]=atte1[pmtleft+j]*exp(tpos*lambda1[pmtleft+j]) +
1223              atte2[pmtleft+j]*exp(tpos*lambda2[pmtleft+j]) ;
1224            QhitPmt_pC[j]= QhitPad_pC[j]*Atten[j]/knorm[j];
1225            if (DEBUG) {
1226              cout<<"pmtleft "<<pmtleft<<" j "<<j<<endl;
1227              cout<<" atte1 "<<atte1[pmtleft+j]<<"lambda1 "<<lambda1[pmtleft+j]<<" atte2 "<<atte2[pmtleft+j]<<"lambda2 "<<lambda2[pmtleft+j] <<endl;    
1228              cout<<j<<" tpos "<<tpos<<" knorm "<<knorm[j]<<" "<<Atten[j]<<" "<<"QhitPmt_pC "<<QhitPmt_pC[j]<<endl;    
1229            }
1230          }
1231          
1232          if (DEBUG)
1233            cout<<"Npho "<<Npho<<" QhitPmt_pC "<<QhitPmt_pC[0]<<" "<<QhitPmt_pC[1]<<endl;  
1234    
1235        Atten[j]=knorm[j]*(atte1[pmtleft+j]*exp(tpos/lambda1[pmtleft+j]) +        QevePmt_pC[pmtleft]  += QhitPmt_pC[0];
1236                           atte2[pmtleft+j]*exp(tpos/lambda2[pmtleft+j]));        QevePmt_pC[pmtright] += QhitPmt_pC[1];
1237                
1238        QhitPmt_pC[j]= QhitPad_pC[j]*Atten[j];        // TDC
1239      }        t2 = t2 + fabs(path[0]/veff) + s_l_g[ip]/veff1 ;  // Signal reaches PMT
1240              t1 = t1 + fabs(path[1]/veff) + s_l_g[ip]/veff1;
1241      QevePmt_pC[pmtleft]  += QhitPmt_pC[0];        
1242      QevePmt_pC[pmtright] += QhitPmt_pC[1];        TRandom r;
1243              t1 = r.Gaus(t1,dt1);  //apply gaussian error  dt
1244    } // ****************************************       fine loop sugli hit        t2 = r.Gaus(t2,dt1);  //apply gaussian error  dt
1245          
1246          t1 = t1 + c1_S[pmtleft] ;  // Signal reaches Discriminator ,TDC starts  to run
1247          t2 = t2 + c1_S[pmtright] ;
1248          
1249          // check if signal is above threshold
1250          // then check if tdcpmt is already filled by another hit...
1251          // only re-fill if time is smaller
1252          
1253          if (QhitPmt_pC[0] > thresh) {
1254            if (tdcpmt[pmtleft] == 1000.) {  // fill for the first time
1255              tdcpmt[pmtleft] = t1;
1256              tdc[pmtleft] = t1 + c2_S[pmtleft] ;  // Signal reaches Coincidence
1257            }
1258            if (tdcpmt[pmtleft] < 1000.) // is already filled!
1259              if (t1 <  tdcpmt[pmtleft]) {
1260                tdcpmt[pmtleft] = t1;
1261                t1 = t1 + c2_S[pmtleft] ;  // Signal reaches Coincidence
1262                tdc[pmtleft] = t1;
1263              }
1264          }      
1265          if (QhitPmt_pC[1] > thresh) {
1266            if (tdcpmt[pmtright] == 1000.) {  // fill for the first time
1267              tdcpmt[pmtright] = t2;
1268              tdc[pmtright] = t2 + c2_S[pmtright] ;  // Signal reaches Coincidence
1269            }
1270            if (tdcpmt[pmtright] < 1000.)  // is already filled!
1271              if (t2 <  tdcpmt[pmtright]) {
1272                tdcpmt[pmtright] = t2;
1273                t2 = t2 + c2_S[pmtright] ;
1274                tdc[pmtright] = t2;
1275              }      
1276          }
1277    
1278          if (DEBUG)
1279            cout<<nh<<" "<<Timetof[nh]<<" "<<t1<<" "<<t2<<endl;
1280    
1281        } // ip < 6
1282    
1283      }; // ****************************************       end loop over hits
1284        
1285      // ======  ADC ======
1286    for(Int_t i=0; i<48; i++){    for(Int_t i=0; i<48; i++){
1287      if(QevePmt_pC[i] != 0.){      if(QevePmt_pC[i] != 0.){
1288        ADCtof[i]= (Int_t)(ADC_pC*QevePmt_pC[i] + ADCoffset);        ADCtof[i]= (Int_t)(ADC_pC*QevePmt_pC[i] + ADCoffset);
1289        if(ADCtof[i]> ADClast) ADCtof[i]=ADClast;        if(ADCtof[i]> ADClast) ADCtof[i]=ADClast;
1290      } else      } else
1291        ADCtof[i]= ADClast;        ADCtof[i]= ADClast;
1292    };    }
1293    
1294        
1295      // ======  build  TDC coincidence  ======
1296    
1297      Float_t t_coinc = 0;
1298      Int_t ilast = 100;
1299      for (Int_t ii=0; ii<48;ii++)
1300        if (tdc[ii] > t_coinc) {
1301          t_coinc = tdc[ii];
1302          ilast = ii;
1303        }
1304      
1305      //     cout<<ilast<<" "<<t_coinc<<endl;
1306      //     At t_coinc  trigger condition is fulfilled
1307      
1308      for (Int_t ii=0; ii<48;ii++){
1309        //      if (tdc[ii] != 0) tdc1[ii] = t_coinc - tdc[ii];   // test 1
1310        if (tdc[ii] != 0) tdc1[ii] = t_coinc - tdcpmt[ii];  // test 2
1311        tdc1[ii] = tdc1[ii]/tdcres[ii];                     // divide by TDC resolution
1312        if (tdc[ii] != 0) tdc1[ii] = tdc1[ii] + c3_S[ii];  // add cable length c3
1313    
1314      } // missing parenthesis inserted! (Silvio)
1315    
1316      for(Int_t i=0; i<48; i++){
1317        if(tdc1[i] != 0.){
1318          TDCint[i]=(Int_t)tdc1[i];
1319          if (DEBUG)
1320            cout<<i<<" "<<TDCint[i]<<endl;
1321          //ADC[i]= ADC_pC * QevePmt_pC[i] + ADCoffset;
1322          //if(ADC[i]> ADClast) ADC[i]=ADClast;
1323        } else
1324          TDCint[i]= TDClast;
1325      }
1326    
1327      if (DEBUG)
1328        cout<<"-----------"<<endl;
1329    
1330    // ======  write fDataTof  =======
1331    UChar_t tofBin;    UChar_t tofBin;
   // --- write fDataTof:  
1332    for (Int_t j=0; j < 12; j++){    for (Int_t j=0; j < 12; j++){
1333      Int_t j12=j*12;      Int_t j12=j*12;
1334      fDataTof[j12+0]=0x00;   // TDC_ID      fDataTof[j12+0]=0x00;   // TDC_ID
# Line 1204  Int_t Digitizer::DigitizeTOF() { Line 1341  Int_t Digitizer::DigitizeTOF() {
1341        fDataTof[jk12+4] = Bin2GrayTof(tofBin,fDataTof[jk12+4]);        fDataTof[jk12+4] = Bin2GrayTof(tofBin,fDataTof[jk12+4]);
1342        tofBin=(UChar_t)(ADCtof[k+4*j]%256);   // ADC# (lsb)        tofBin=(UChar_t)(ADCtof[k+4*j]%256);   // ADC# (lsb)
1343        fDataTof[jk12+5] = Bin2GrayTof(tofBin,fDataTof[jk12+5]);        fDataTof[jk12+5] = Bin2GrayTof(tofBin,fDataTof[jk12+5]);
1344        fDataTof[jk12+6]=0x00;   // TDC# (msb) -- Wolfgang        tofBin=(UChar_t)(TDCint[k+4*j]/256);   // TDC# (msb)
1345        fDataTof[jk12+7]=0x00;   // TDC# (lsb) -- Wolfgang        fDataTof[jk12+6]=Bin2GrayTof(tofBin,fDataTof[jk12+6]);
1346          tofBin=(UChar_t)(TDCint[k+4*j]%256);   // TDC# (lsb)
1347          fDataTof[jk12+7]=Bin2GrayTof(tofBin,fDataTof[jk12+7]);
1348      };      };
1349      fDataTof[j12+20]=0x00;   // TEMP1      fDataTof[j12+20]=0x00;   // TEMP1
1350      fDataTof[j12+21]=0x00;   // TEMP2      fDataTof[j12+21]=0x00;   // TEMP2
# Line 1271  void Digitizer::Paddle2Pmt(Int_t plane, Line 1410  void Digitizer::Paddle2Pmt(Int_t plane,
1410      somma+=pads[j];      somma+=pads[j];
1411    padid=paddle+somma;    padid=paddle+somma;
1412    *pl = padid*2;    *pl = padid*2;
1413    *pr = *pr + 1;    //  *pr = *pr + 1;
1414      *pr = *pl + 1; // WM
1415  };  };
1416    
1417  void Digitizer::DigitizeAC() {  void Digitizer::DigitizeAC() {
1418    // created:  J. Conrad, KTH    // created:  J. Conrad, KTH
1419    // modified: S. Orsi, INFN Roma2    // modified: S. Orsi, INFN Roma2
1420      // fDataAC[0-63]:   main AC board
1421      // fDataAC[64-127]: extra AC board
1422    
1423    fDataAC[0] = 0xACAC;    fDataAC[0] = 0xACAC;
1424    fDataAC[64]= 0xACAC;    fDataAC[64]= 0xACAC;
1425    fDataAC[1] = 0xAC11;   // main card    fDataAC[1] = 0xAC11;  
1426    fDataAC[65] = 0xAC22;   // extra card    fDataAC[65] = 0xAC22;  
1427    
1428    // the third word is a status word (dummy)    // the third word is a status word (dummy: "no errors are present in the AC boards")
1429    fDataAC[2] = 0xFFFF; //FFEF?    fDataAC[2] = 0xFFFF; //FFEF?
1430    fDataAC[66] = 0xFFFF;    fDataAC[66] = 0xFFFF;
1431    
1432    const UInt_t nReg = 6;    const UInt_t nReg = 6;
1433    
1434    // Registers (dummy)    // FPGA Registers (dummy)
1435    for (UInt_t i=0; i<=nReg; i++){    for (UInt_t i=0; i<=nReg; i++){
1436      fDataAC[i+4] = 0xFFFF;      fDataAC[i+4] = 0xFFFF;
1437      fDataAC[i+68] = 0xFFFF;      fDataAC[i+68] = 0xFFFF;
# Line 1300  void Digitizer::DigitizeAC() { Line 1442  void Digitizer::DigitizeAC() {
1442    fDataAC[63] = 0xABCD;    fDataAC[63] = 0xABCD;
1443    fDataAC[127] = 0xABCD;    fDataAC[127] = 0xABCD;
1444    
1445    // shift registers, which one is with respect to PMT, where in    // shift registers (moved to the end of the routine)
   // shift registers is a question of time relative trigger  
   // In level2: hitmap, hitmap-status (synchronised with a trigger),  
   // status  
             
   for (UInt_t i=0; i<=15; i++){  
     fDataAC[i+11] = 0x0000;    
     fDataAC[i+75] = 0x0000;  
   }  
1446    
1447    // singles counters are dummy    Int_t evntLSB=Ievnt%65536;
1448      Int_t evntMSB=(Int_t)(Ievnt/65536);
1449    
1450    for (UInt_t i=0; i<=16; i++){    // singles counters are dummy
1451      fDataAC[i+26] = 0x0000;      for (UInt_t i=0; i<=15; i++){  //SO Oct '07: // for (UInt_t i=0; i<=16; i++){
1452      fDataAC[i+90] = 0x0000;      //     fDataAC[i+26] = 0x0000;  
1453    }      //     fDataAC[i+90] = 0x0000;
1454        fDataAC[i+26] = evntLSB;
1455        fDataAC[i+90] = evntLSB;
1456      };
1457        
1458    // coincidences are dummy    // coincidences are dummy (increment by 1 at each event)
1459      // for (UInt_t i=0; i<=7; i++){
1460      //    fDataAC[i+42] = 0x0000;
1461      //    fDataAC[i+106] = 0x0000;
1462      //   }
1463    for (UInt_t i=0; i<=7; i++){    for (UInt_t i=0; i<=7; i++){
1464      fDataAC[i+42] = 0x0000;      fDataAC[i+42] = evntLSB;
1465      fDataAC[i+106] = 0x0000;      fDataAC[i+106] = evntLSB;
1466    }    };
1467    
1468    // increments for every trigger might be needed at some point.    // increments for every trigger might be needed at some point.
1469    // dummy for now    // dummy for now
1470    fDataAC[50]  = 0x0000;    fDataAC[50]  = 0x0000;
1471    fDataAC[114] = 0x0000;    fDataAC[114] = 0x0000;
1472    
1473    // dummy FPGA clock    // dummy FPGA clock (increment by 1 at each event)
1474      /*    
1475    fDataAC[51] = 0x006C;      fDataAC[51] = 0x006C;
1476    fDataAC[52] = 0x6C6C;      fDataAC[52] = 0x6C6C;
1477    fDataAC[115] = 0x006C;      fDataAC[115] = 0x006C;
1478    fDataAC[116] = 0x6C6C;      fDataAC[116] = 0x6C6C;
1479      */
1480      if (Ievnt<=0xFFFF) {
1481        fDataAC[51] = 0x0000;
1482        fDataAC[52] = Ievnt;
1483        fDataAC[115] = 0x0000;
1484        fDataAC[116] = Ievnt;
1485      } else {
1486        fDataAC[51] = evntMSB;
1487        fDataAC[52] = evntLSB;
1488        fDataAC[115] = fDataAC[51];
1489        fDataAC[116] = fDataAC[52];
1490      }
1491    
1492    // dummy temperatures    // dummy temperatures
1493    fDataAC[53] = 0x0000;    fDataAC[53] = 0x0000;
# Line 1350  void Digitizer::DigitizeAC() { Line 1502  void Digitizer::DigitizeAC() {
1502      fDataAC[i+119] = 0x1A13;      fDataAC[i+119] = 0x1A13;
1503    }    }
1504        
1505    // We activate all branches. Once the digitization algorithm    // We activate all branches. Once the digitization algorithm is determined
1506    // is determined only the branches need to activated which involve needed    // only the branches that involve needed information will be activated
   // information  
1507        
1508      fhBookTree->SetBranchAddress("Ievnt",&Ievnt);
1509    fhBookTree->SetBranchStatus("Nthcat",1);    fhBookTree->SetBranchStatus("Nthcat",1);
1510    fhBookTree->SetBranchStatus("Iparcat",1);    fhBookTree->SetBranchStatus("Iparcat",1);
1511    fhBookTree->SetBranchStatus("Icat",1);    fhBookTree->SetBranchStatus("Icat",1);
# Line 1399  void Digitizer::DigitizeAC() { Line 1551  void Digitizer::DigitizeAC() {
1551    // will fire. We will furthermore assume that both cards read out    // will fire. We will furthermore assume that both cards read out
1552    // identical data.    // identical data.
1553    
1554    // If you develop you digitization algorithm, you should start by    // If you develop your digitization algorithm, you should start by
1555    // identifying the information present in level2 (post-darth-vader)    // identifying the information present in level2 (post-darth-vader)
1556    // data.    // data.
1557    
# Line 1413  void Digitizer::DigitizeAC() { Line 1565  void Digitizer::DigitizeAC() {
1565    };    };
1566    
1567    if (Nthcat>50 || Nthcas>50 || Nthcard>50)    if (Nthcat>50 || Nthcas>50 || Nthcard>50)
1568      printf("Error! NthAC out of range!\n\n");      printf("*** ERROR AC! NthAC out of range!\n\n");
1569    
1570      // energy dependence on position (see file AcFitOutputDistancePmt.C by S.Orsi)
1571      // based on J.Lundquist's calculations (PhD thesis, page 94)
1572      // function: [0]+[1]*atan([2]/(x+1)), where the 3 parameters are:
1573      //   8.25470e-01   +- 1.79489e-02
1574      //   6.41609e-01   +- 2.65846e-02
1575      //   9.81177e+00   +- 1.21284e+00
1576      // hp: 1 minimum ionising particle at 35cm from the PMT releases 1mip
1577      //
1578      // NB: the PMT positions are needed!
1579    
1580    // look in CAT    // look in CAT
1581    //  for (UInt_t k= 0;k<50;k++){    //  for (UInt_t k= 0;k<50;k++){
# Line 1477  void Digitizer::DigitizeAC() { Line 1639  void Digitizer::DigitizeAC() {
1639    
1640    fDataAC[67] = fDataAC[3];    fDataAC[67] = fDataAC[3];
1641    
1642      // shift registers
1643      // the central bin is equal to the hitmap, all other bins in the shift register are 0
1644      for (UInt_t i=0; i<=15; i++){
1645        fDataAC[i+11] = 0x0000;  
1646        fDataAC[i+75] = 0x0000;
1647      }
1648      fDataAC[18] = fDataAC[3];
1649      fDataAC[82] = fDataAC[3];
1650    
1651    //    for (Int_t i=0; i<fACbuffer; i++){    //    for (Int_t i=0; i<fACbuffer; i++){
1652    //    printf("%0x  ",fDataAC[i]);      //    printf("%0x  ",fDataAC[i]);  
1653    //    if ((i+1)%8 ==0) cout << endl;    //    if ((i+1)%8 ==0) cout << endl;
# Line 1484  void Digitizer::DigitizeAC() { Line 1655  void Digitizer::DigitizeAC() {
1655  };  };
1656    
1657    
1658  void Digitizer::DigitizeND(){  void Digitizer::DigitizeS4(){
1659    // creato:  S. Borisov, INFN Roma2 e MEPHI, Sept 2007    Int_t DEBUG=0;
1660    // 4 bytes: 16bit header, 8bit trigPhysics, 16bit up&low background    // creato:  S. Borisov, INFN Roma2 e MEPHI, Sett 2007
1661      TString ciao,modo="ns";
1662    // ND header    Int_t i,j,t,NdF,pmt,NdFT,S4,S4v=0,S4p=32;
1663    fDataND[0] = 0x0000;    Float_t E0,E1=1e-6,Ert,X,Y,Z,x,y,z,V[3],Xs[2],Ys[2],Zs[2],Yp[6],q,w,p=0.1,l,l0=500;
1664    fDataND[1] = 0x000F;    Xs[0]=-24.1;
1665      Xs[1]=24.1;
1666      Ys[0]=-24.1;
1667      Ys[1]=24.1;
1668      Zs[0]=-0.5;
1669      Zs[1]=0.5;
1670      Yp[0]=-20.;
1671      Yp[2]=-1.;
1672      Yp[4]=17.;
1673      for(i=0;i<3;i++)
1674        Yp[2*i+1]=Yp[2*i]+3;
1675      srand(time(NULL));
1676      // --- activate branches:
1677      fhBookTree->SetBranchStatus("Nthtof",1);
1678      fhBookTree->SetBranchStatus("Ipltof",1);
1679      fhBookTree->SetBranchStatus("Ipaddle",1);
1680      
1681      fhBookTree->SetBranchStatus("Xintof",1);
1682      fhBookTree->SetBranchStatus("Yintof",1);
1683      fhBookTree->SetBranchStatus("Xouttof",1);
1684      fhBookTree->SetBranchStatus("Youttof",1);
1685        
1686      fhBookTree->SetBranchStatus("Ereltof",1);
1687      fhBookTree->SetBranchStatus("Timetof",1);
1688      NdFT=0;
1689      Ert=0;
1690      for(i=0;i<Nthtof;i++){
1691        if(Ipltof[i]!=6) continue;
1692        Ert+=Ereltof[i];
1693              
1694        if(modo=="ns") continue;
1695        NdF=Int_t(Ereltof[i]/E1);
1696        NdFT=0;
1697        X=Xintof[i];
1698        Y=Yintof[i];
1699        Z=(Float_t)(random())/(Float_t)(0x7fffffff)-0.5;
1700        //cout<<"XYZ "<<X<<"  "<<Y<<"  "<<Z<<endl;
1701        for(j=0;j<NdF;j++){
1702          q=(Float_t)random()/(Float_t)0x7fffffff;
1703          w=(Float_t)random()/(Float_t)0x7fffffff;
1704          // cout<<"qw "<<q<<" "<<w<<endl;
1705          V[0]=p*cos(6.28318*q);
1706          V[1]=p*sin(6.28318*q);
1707          V[2]=p*(2.*w-1.);
1708          pmt=0;
1709          x=X;
1710          y=Y;
1711          z=Z;
1712          while(pmt==0 && (x>Xs[0] && x<Xs[1])&&(y>Ys[0] && y<Ys[1])&&(z>Zs[0] && z<Zs[1])){
1713            l=0;
1714            while(pmt==0 && (x>Xs[0] && x<Xs[1])&&(y>Ys[0] && y<Ys[1])&&(z>Zs[0] && z<Zs[1])){
1715              x+=V[0];
1716              y+=V[1];
1717              z+=V[2];
1718              l+=p;
1719              //cout<<x<<"  "<<y<<"  "<<z<<"  "<<l<<endl;
1720              //cin>>ciao;
1721            }
1722            if((x<Xs[0]+p || x>Xs[1]-p)&&(y>Ys[0]+p && y<Ys[1]-p)&&(z>Zs[0]+p && z<Zs[1]-p)){
1723              for(t=0;t<3;t++){
1724                if(y>=Yp[2*t] && y<Yp[2*t+1]){
1725                  if(pmt==0)NdFT++;
1726                  pmt=1;
1727                  //cout<<NdFT<<endl;
1728                  break;
1729                }
1730              }
1731              if(pmt==1)break;
1732              V[0]=-V[0];
1733            }
1734            q=(Float_t)random()/(Float_t)0x7fffffff;
1735            w=1-exp(-l/l0);
1736            if(q<w)break;
1737            q=(Float_t)random()/(Float_t)0x7fffffff;
1738            w=0.5;
1739            if(q<w)break;
1740            if((x>Xs[0]+p && x<Xs[1]-p)&&(y<Ys[0]+p || y>Ys[1]-p)&&(z>Zs[0]+p && z<Zs[1]-p))V[1]=-V[1];
1741            if((x>Xs[0]+p && x<Xs[1]-p)&&(y>Ys[0]+p && y<Ys[1]-p)&&(z<Zs[0]+p || z>Zs[1]-p))V[2]=-V[2];
1742            x+=V[0];
1743            y+=V[1];
1744            z+=V[2];
1745            l=0;
1746            //cout<<x<<"  "<<y<<"  "<<z<<"  "<<l<<endl;
1747                    //cin>>ciao;
1748          }
1749        }
1750      }
1751      Ert=Ert/0.002;
1752      q=(Float_t)(random())/(Float_t)0x7fffffff;
1753      w=0.7;
1754      //E0=(Float_t)(4064./7.);
1755      E0=4064./7.;
1756      if(Ert<1) S4=0;
1757      else S4=(Int_t)(4064.*(1.-exp(-(Ert-1.)/E0)));
1758      i=S4/4;
1759      if(S4%4==0)
1760        S4v=S4+S4p;
1761      else if(S4%4==1){
1762        if(q<w) S4v=S4-1+S4p;
1763        else S4v=S4+1+S4p;
1764      } else if(S4%4==2) S4v=S4+S4p;
1765      else if(S4%4==3){
1766        if(q<w) S4v=S4+1+S4p;
1767        else S4v=S4-1+S4p;
1768      }
1769      if (DEBUG)
1770        cout<<"Ert_S4 = " << Ert << " --- S4v = " << S4v << endl;
1771      fDataS4[0]=S4v;//0xf028;
1772      fDataS4[1]=0xd800;
1773      fDataS4[2]=0x0300;
1774      //cout<<"  PMT  "<<NdFT<<"  "<<NdF<<endl;
1775      //cin>>ciao;
1776    }
1777    
1778    
1779    
1780    void Digitizer::DigitizeND(){
1781      // creato:  S. Borisov, INFN Roma2 e MEPHI, Sett 2007
1782      Int_t i=0;
1783      UShort_t NdN=0;
1784    fhBookTree->SetBranchStatus("Nthnd",1);    fhBookTree->SetBranchStatus("Nthnd",1);
1785    fhBookTree->SetBranchStatus("Itubend",1);    fhBookTree->SetBranchStatus("Itubend",1);
1786    fhBookTree->SetBranchStatus("Iparnd",1);      fhBookTree->SetBranchStatus("Iparnd",1);  
# Line 1505  void Digitizer::DigitizeND(){ Line 1794  void Digitizer::DigitizeND(){
1794    fhBookTree->SetBranchStatus("Timend",1);    fhBookTree->SetBranchStatus("Timend",1);
1795    fhBookTree->SetBranchStatus("Pathnd",1);    fhBookTree->SetBranchStatus("Pathnd",1);
1796    fhBookTree->SetBranchStatus("P0nd",1);    fhBookTree->SetBranchStatus("P0nd",1);
1797      //cout<<"n="<<Nthnd<<"  "<<NdN<<"\n";
1798      for(i=0;i<Nthnd;i++){
1799        if(Iparnd[i]==13){
1800          NdN++;
1801        }
1802      }
1803      //NdN=100; //only for debug
1804    
1805    UShort_t NdN=0;    for(i=0;i<3;i++){
1806    for(Int_t i=0;i<Nthnd;i++)      fDataND[2*i]=0x0000;
1807      if(Iparnd[i]==13)      fDataND[2*i+1]=0x010F;
1808        NdN++;          }
1809      fDataND[0]=0xFF00 & (256*NdN);
   NdN=10; // test!  
   fDataND[2]=0x0F00 & (NdN*256);  
   //fDataND[2]=0xFFFF; //test  
   fDataND[2]=0x0000; //background neutrons  
1810  }  }
1811    
1812    
# Line 1530  void Digitizer::DigitizeDummy() { Line 1822  void Digitizer::DigitizeDummy() {
1822      //   printf("%0x  ",fDataDummy[i]);        //   printf("%0x  ",fDataDummy[i]);  
1823      //if ((i+1)%8 ==0) cout << endl;      //if ((i+1)%8 ==0) cout << endl;
1824    }    }
     
   
   
1825  };  };
1826    
1827    
# Line 1559  void Digitizer::WriteData(){ Line 1848  void Digitizer::WriteData(){
1848    swab(fDataTrack,temp,sizeof(UShort_t)*fTracklength);  // WE MUST SWAP THE BYTES!!!    swab(fDataTrack,temp,sizeof(UShort_t)*fTracklength);  // WE MUST SWAP THE BYTES!!!
1849    fOutputfile.write(reinterpret_cast<char*>(temp),sizeof(UShort_t)*fTracklength);    fOutputfile.write(reinterpret_cast<char*>(temp),sizeof(UShort_t)*fTracklength);
1850    fTracklength=0;    fTracklength=0;
1851    // S4     // padding to 64 bytes
   // ...to be done...  
   // ND  
   memset(temp,0,sizeof(UShort_t)*1000000);  
   swab(fDataND,temp,sizeof(UShort_t)*4);  // WE MUST SWAP THE BYTES!!!  
   fOutputfile.write(reinterpret_cast<char*>(temp),sizeof(UShort_t)*4);  
   
   //  
   //  fOutputfile.write(reinterpret_cast<char*>(fDataDummy),sizeof(UShort_t)*fDummybuffer);  
   //  
   // padding to 64 bytes  
1852    //    //
1853    if ( fPadding ){    if ( fPadding ){
1854      fOutputfile.write(reinterpret_cast<char*>(fDataPadding),sizeof(UChar_t)*fPadding);      fOutputfile.write(reinterpret_cast<char*>(fDataPadding),sizeof(UChar_t)*fPadding);
1855    };    };
1856    //    // S4
1857      memset(temp,0,sizeof(UShort_t)*1000000);
1858      swab(fDataS4,temp,sizeof(UShort_t)*fS4buffer);  // WE MUST SWAP THE BYTES!!!
1859      fOutputfile.write(reinterpret_cast<char*>(temp),sizeof(UShort_t)*fS4buffer);
1860      // ND
1861      memset(temp,0,sizeof(UShort_t)*1000000);
1862      swab(fDataND,temp,sizeof(UShort_t)*fNDbuffer);  // WE MUST SWAP THE BYTES!!!
1863      fOutputfile.write(reinterpret_cast<char*>(temp),sizeof(UShort_t)*fNDbuffer);
1864  };  };
1865    
1866    
# Line 1667  Int_t Nstrip; Line 1953  Int_t Nstrip;
1953    
1954    
1955    
     
1956    Float_t ADCfull;    Float_t ADCfull;
1957      Int_t iladd=0;
1958    for (Int_t ix=0; ix<Nstrpx;ix++) {    for (Int_t ix=0; ix<Nstrpx;ix++) {
1959    Iview=Npstripx[ix]*2-1;    Iview=Npstripx[ix]*2-1;
1960    Nstrip=(Int_t)Istripx[ix]-1;    Nstrip=(Int_t)Istripx[ix]-1;
1961    ADCfull=AdcTrack[Iview][Nstrip] += Qstripx[ix]*fMipCor;    if(Nstrip<fNstrips_ladder) iladd=0;
1962      if((Nstrip>=fNstrips_ladder)&&(Nstrip<2*fNstrips_ladder)) iladd=1;
1963      if((Nstrip>=2*fNstrips_ladder)&&(Nstrip<3*fNstrips_ladder)) iladd=2;
1964      ADCfull=AdcTrack[Iview][Nstrip] += Qstripx[ix]*fMipCor[iladd][Iview];
1965    AdcTrack[Iview][Nstrip] *= SaturationTrack(ADCfull);    AdcTrack[Iview][Nstrip] *= SaturationTrack(ADCfull);
1966    
1967    };    };
# Line 1681  Int_t Nstrip; Line 1970  Int_t Nstrip;
1970    for (Int_t iy=0; iy<Nstrpy;iy++) {    for (Int_t iy=0; iy<Nstrpy;iy++) {
1971    Iview=Npstripy[iy]*2-2;    Iview=Npstripy[iy]*2-2;
1972    Nstrip=(Int_t)Istripy[iy]-1;    Nstrip=(Int_t)Istripy[iy]-1;
1973    ADCfull=AdcTrack[Iview][Nstrip] -= Qstripy[iy]*fMipCor;    if(Nstrip<fNstrips_ladder) iladd=0;
1974      if((Nstrip>=fNstrips_ladder)&&(Nstrip<2*fNstrips_ladder)) iladd=1;
1975      if((Nstrip>=2*fNstrips_ladder)&&(Nstrip<3*fNstrips_ladder)) iladd=2;
1976      ADCfull=AdcTrack[Iview][Nstrip] -= Qstripy[iy]*fMipCor[iladd][Iview];
1977    AdcTrack[Iview][Nstrip] *= SaturationTrack(ADCfull);    AdcTrack[Iview][Nstrip] *= SaturationTrack(ADCfull);
1978    
1979    };      };  
# Line 1929  std:: cout << "Entering LoadTrackCalib " Line 2221  std:: cout << "Entering LoadTrackCalib "
2221    
2222  void Digitizer::LoadMipCor() {  void Digitizer::LoadMipCor() {
2223  std:: cout << "Entering LoadMipCor" << endl;  std:: cout << "Entering LoadMipCor" << endl;
2224      Float_t xfactor=1./151.6*1.04;
2225      Float_t yfactor=1./152.1;
2226    
2227      fMipCor[0][0]=140.02*yfactor;
2228      fMipCor[0][1]=140.99*xfactor;
2229      fMipCor[0][2]=134.48*yfactor;
2230      fMipCor[0][3]=144.41*xfactor;
2231      fMipCor[0][4]=140.74*yfactor;
2232      fMipCor[0][5]=142.28*xfactor;
2233      fMipCor[0][6]=134.53*yfactor;
2234      fMipCor[0][7]=140.63*xfactor;
2235      fMipCor[0][8]=135.55*yfactor;
2236      fMipCor[0][9]=138.00*xfactor;
2237      fMipCor[0][10]=154.95*yfactor;
2238      fMipCor[0][11]=158.44*xfactor;
2239      
2240      
2241      fMipCor[1][0]=136.07*yfactor;
2242      fMipCor[1][1]=135.59*xfactor;
2243      fMipCor[1][2]=142.69*yfactor;
2244      fMipCor[1][3]=138.19*xfactor;
2245      fMipCor[1][4]=137.35*yfactor;
2246      fMipCor[1][5]=140.23*xfactor;
2247      fMipCor[1][6]=153.15*yfactor;
2248      fMipCor[1][7]=151.42*xfactor;
2249      fMipCor[1][8]=129.76*yfactor;
2250      fMipCor[1][9]=140.63*xfactor;
2251      fMipCor[1][10]=157.87*yfactor;
2252      fMipCor[1][11]=153.64*xfactor;
2253    
2254      fMipCor[2][0]=134.98*yfactor;
2255      fMipCor[2][1]=143.95*xfactor;
2256      fMipCor[2][2]=140.23*yfactor;
2257      fMipCor[2][3]=138.88*xfactor;
2258      fMipCor[2][4]=137.95*yfactor;
2259      fMipCor[2][5]=134.87*xfactor;
2260      fMipCor[2][6]=157.56*yfactor;
2261      fMipCor[2][7]=157.31*xfactor;
2262      fMipCor[2][8]=141.37*yfactor;
2263      fMipCor[2][9]=143.39*xfactor;
2264      fMipCor[2][10]=156.15*yfactor;
2265      fMipCor[2][11]=158.79*xfactor;
2266    
2267  /*  /*
2268    for (Int_t j=0; j<fNviews;j++) {    for (Int_t j=0; j<fNviews;j++) {
2269      for (Int_t i=0; i<fNstrips_view;i++) {      for (Int_t i=0; i<fNstrips_view;i++) {
# Line 2221  fTracklength=0; Line 2556  fTracklength=0;
2556    
2557  };  };
2558    
 Float_t Digitizer::SaturationTrack(Float_t ADC) {  
2559    
2560  Float_t SatFact=1.;  Float_t Digitizer::SaturationTrack(Float_t ADC) {
2561  return SatFact;    Float_t SatFact=1.;
2562      if(ADC<70.) { SatFact=80./ADC; };
2563      if(ADC>3000.) { SatFact=3000./ADC; };
2564      return SatFact;
2565  };  };
2566    
2567    
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.3
  ViewVC Help
Powered by ViewVC 1.1.23