/[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.2 by orsi, Fri Sep 28 10:46:23 2007 UTC revision 1.3 by orsi, Thu Oct 11 11:29:21 2007 UTC
# Line 260  void Digitizer::Loop() { Line 260  void Digitizer::Loop() {
260      DigitizeAC();      DigitizeAC();
261      DigitizeCALO();      DigitizeCALO();
262      DigitizeTrack();      DigitizeTrack();
263      //DigitizeS4();      DigitizeS4();
264      DigitizeND();      DigitizeND();
265        //        //
266      // Add padding to 64 bits      // Add padding to 64 bits
# Line 269  void Digitizer::Loop() { Line 269  void Digitizer::Loop() {
269  //  //
270      // Create CPU header, we need packet type (0x10 = physics data) and packet length.      // 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;      UInt_t length=2*(fCALOlength+fACbuffer+fTracklength+fNDbuffer+fS4buffer)+fPadding+fTOFbuffer+fTRIGGERbuffer;
273      UInt_t length=2*(fCALOlength+fACbuffer+fTracklength+fNDbuffer)+fPadding+fTOFbuffer+fTRIGGERbuffer;      //UInt_t length=2*(fCALOlength+fACbuffer+fTracklength+fNDbuffer)+fPadding+fTOFbuffer+fTRIGGERbuffer;
274      DigitizePSCU(length,0x10);      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();
# Line 1037  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 1096  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 1126  Int_t Digitizer::DigitizeTOF() { Line 1128  Int_t Digitizer::DigitizeTOF() {
1128    Float_t  tdc[48],tdc1[48],tdcpmt[48];    Float_t  tdc[48],tdc1[48],tdcpmt[48];
1129    for(Int_t i=0; i<48; i++)    for(Int_t i=0; i<48; i++)
1130      tdcpmt[i] = 1000.;      tdcpmt[i] = 1000.;
1131    Float_t thresh=1.; // to be defined better... (Wolfgang)    Float_t thresh=10.; // to be defined better... (Wolfgang)
1132    
1133      // === TDC: simulate timing for each paddle      // === TDC: simulate timing for each paddle
1134      Float_t dt1 = 285.e-12 ;   // single PMT resolution      Float_t dt1 = 285.e-12 ;   // single PMT resolution
# Line 1139  Int_t Digitizer::DigitizeTOF() { Line 1141  Int_t Digitizer::DigitizeTOF() {
1141      for(Int_t j=0;j<48;j++)  c2_S[j] = c2_S[j]*tdcres[j];      for(Int_t j=0;j<48;j++)  c2_S[j] = c2_S[j]*tdcres[j];
1142       // ih = 0 + i1;  // not used?? (Silvio)       // ih = 0 + i1;  // not used?? (Silvio)
1143    
1144    /* **********************************  inizio loop sugli hit */    /* **********************************  start loop over hits */
1145        
1146    for(Int_t nh=0; nh<Nthtof; nh++){    for(Int_t nh=0; nh<Nthtof; nh++){
1147            
# Line 1158  Int_t Digitizer::DigitizeTOF() { Line 1160  Int_t Digitizer::DigitizeTOF() {
1160      t1 = Timetof[nh] ;  // Start      t1 = Timetof[nh] ;  // Start
1161      t2 = Timetof[nh] ;      t2 = Timetof[nh] ;
1162    
1163      // Donatella          // Donatella: redefinition plane and pad for vectors in C
     // ridefiniz. piano e pad per i vettori 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      // TDC        t2 = r.Gaus(t2,dt1);  //apply gaussian error  dt
1245      t2 = t2 + fabs(path[0]/veff) + s_l_g[ip]/veff1 ;  // Signal reaches PMT        
1246      t1 = t1 + fabs(path[1]/veff) + s_l_g[ip]/veff1;        t1 = t1 + c1_S[pmtleft] ;  // Signal reaches Discriminator ,TDC starts  to run
1247              t2 = t2 + c1_S[pmtright] ;
1248      TRandom r;        
1249      t1 = r.Gaus(t1,dt1);  //apply gaussian error  dt        // check if signal is above threshold
1250      t2 = r.Gaus(t2,dt1);  //apply gaussian error  dt        // then check if tdcpmt is already filled by another hit...
1251              // only re-fill if time is smaller
1252      t1 = t1 + c1_S[pmtleft] ;  // Signal reaches Discriminator ,TDC starts  to run        
1253      t2 = t2 + c1_S[pmtright] ;        if (QhitPmt_pC[0] > thresh) {
1254                if (tdcpmt[pmtleft] == 1000.) {  // fill for the first time
     // check if signal is above threshold  
     // then check if tdcpmt is already filled by another hit...  
     // only re-fill if time is smaller  
       
     if (QhitPmt_pC[0] > thresh)  
       if (tdcpmt[pmtleft] < 1000.) // is already filled!  
         if (t1 <  tdcpmt[pmtleft]) {  
1255            tdcpmt[pmtleft] = t1;            tdcpmt[pmtleft] = t1;
1256            t1 = t1 + c2_S[pmtleft] ;  // Signal reaches Coincidence            tdc[pmtleft] = t1 + c2_S[pmtleft] ;  // Signal reaches Coincidence
1257            tdc[pmtleft] = t1;          }
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          }          }
       
     if (QhitPmt_pC[1] > thresh)  
1270          if (tdcpmt[pmtright] < 1000.)  // is already filled!          if (tdcpmt[pmtright] < 1000.)  // is already filled!
1271            if (t2 <  tdcpmt[pmtright]) {            if (t2 <  tdcpmt[pmtright]) {
1272              tdcpmt[pmtright] = t2;              tdcpmt[pmtright] = t2;
1273              t2 = t2 + c2_S[pmtright] ;              t2 = t2 + c2_S[pmtright] ;
1274              tdc[pmtright] = t2;              tdc[pmtright] = t2;
1275            }                  }      
1276              }
   } // ****************************************       end loop over hits  
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 ======    // ======  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.){
# Line 1257  Int_t Digitizer::DigitizeTOF() { Line 1289  Int_t Digitizer::DigitizeTOF() {
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  ======    // ======  build  TDC coincidence  ======
# Line 1284  Int_t Digitizer::DigitizeTOF() { Line 1316  Int_t Digitizer::DigitizeTOF() {
1316    for(Int_t i=0; i<48; i++){    for(Int_t i=0; i<48; i++){
1317      if(tdc1[i] != 0.){      if(tdc1[i] != 0.){
1318        TDCint[i]=(Int_t)tdc1[i];        TDCint[i]=(Int_t)tdc1[i];
1319          if (DEBUG)
1320            cout<<i<<" "<<TDCint[i]<<endl;
1321        //ADC[i]= ADC_pC * QevePmt_pC[i] + ADCoffset;        //ADC[i]= ADC_pC * QevePmt_pC[i] + ADCoffset;
1322        //if(ADC[i]> ADClast) ADC[i]=ADClast;        //if(ADC[i]> ADClast) ADC[i]=ADClast;
1323      } else      } else
1324        TDCint[i]= TDClast;        TDCint[i]= TDClast;
1325    }    }
1326    
1327      if (DEBUG)
1328        cout<<"-----------"<<endl;
1329    
1330  // ======  write fDataTof  =======  // ======  write fDataTof  =======
1331    UChar_t tofBin;    UChar_t tofBin;
1332    for (Int_t j=0; j < 12; j++){    for (Int_t j=0; j < 12; j++){
# Line 1373  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 1402  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 1452  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 1501  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 1515  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 1579  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 1586  void Digitizer::DigitizeAC() { Line 1655  void Digitizer::DigitizeAC() {
1655  };  };
1656    
1657    
   
1658  void Digitizer::DigitizeS4(){  void Digitizer::DigitizeS4(){
1659      Int_t DEBUG=0;
1660    // creato:  S. Borisov, INFN Roma2 e MEPHI, Sett 2007    // creato:  S. Borisov, INFN Roma2 e MEPHI, Sett 2007
1661    TString ciao,modo="ns";    TString ciao,modo="ns";
1662    Int_t i,j,t,NdF,pmt,NdFT,S4,S4v=0,S4p=32;    Int_t i,j,t,NdF,pmt,NdFT,S4,S4v=0,S4p=32;
1663    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.;    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    Xs[0]=-24.1;    Xs[0]=-24.1;
1665    Xs[1]=24.1;    Xs[1]=24.1;
1666    Ys[0]=-24.1;    Ys[0]=-24.1;
# Line 1621  void Digitizer::DigitizeS4(){ Line 1690  void Digitizer::DigitizeS4(){
1690    for(i=0;i<Nthtof;i++){    for(i=0;i<Nthtof;i++){
1691      if(Ipltof[i]!=6) continue;      if(Ipltof[i]!=6) continue;
1692      Ert+=Ereltof[i];      Ert+=Ereltof[i];
1693                    
       
1694      if(modo=="ns") continue;      if(modo=="ns") continue;
1695      NdF=Int_t(Ereltof[i]/E1);      NdF=Int_t(Ereltof[i]/E1);
1696      NdFT=0;      NdFT=0;
1697      X=Xintof[i];      X=Xintof[i];
1698      Y=Yintof[i];      Y=Yintof[i];
1699      Z=((Float_t)random()/(Float_t)0x7fffffff)-0.5;      Z=(Float_t)(random())/(Float_t)(0x7fffffff)-0.5;
1700      //cout<<"XYZ "<<X<<"  "<<Y<<"  "<<Z<<endl;      //cout<<"XYZ "<<X<<"  "<<Y<<"  "<<Z<<endl;
1701      for(j=0;j<NdF;j++){      for(j=0;j<NdF;j++){
1702        q=(Float_t)random()/(Float_t)0x7fffffff;        q=(Float_t)random()/(Float_t)0x7fffffff;
# Line 1681  void Digitizer::DigitizeS4(){ Line 1749  void Digitizer::DigitizeS4(){
1749      }      }
1750    }    }
1751    Ert=Ert/0.002;    Ert=Ert/0.002;
1752    q=(Float_t)(random())/(Float_t)(0x7fffffff);    q=(Float_t)(random())/(Float_t)0x7fffffff;
1753    w=0.7;    w=0.7;
1754    //E0=Float_t(4064)/7;    //E0=(Float_t)(4064./7.);
1755    E0=4064./7.;    E0=4064./7.;
1756    S4=(Int_t)(4064.*(1.-exp(-int(Ert)/E0)));    if(Ert<1) S4=0;
1757    //S4=Ert*7;    else S4=(Int_t)(4064.*(1.-exp(-(Ert-1.)/E0)));
1758    i=S4/4;    i=S4/4;
1759    if(S4%4==0)    if(S4%4==0)
1760      S4v=S4+S4p;      S4v=S4+S4p;
1761    else if(S4%4==1) {    else if(S4%4==1){
1762      if(q<w) S4v=S4-1+S4p;      if(q<w) S4v=S4-1+S4p;
1763      else S4v=S4+1+S4p;      else S4v=S4+1+S4p;
1764    } else if(S4%4==2)    } else if(S4%4==2) S4v=S4+S4p;
     S4v=S4+S4p;  
1765    else if(S4%4==3){    else if(S4%4==3){
1766      if(q<w) S4v=S4+1+S4p;      if(q<w) S4v=S4+1+S4p;
1767      else S4v=S4-1+S4p;      else S4v=S4-1+S4p;
1768    }    }
1769      if (DEBUG)
1770    cout << "Ert= " <<Ert<<"; S4v= "<<S4v<<"; S4= "<<S4<<endl;      cout<<"Ert_S4 = " << Ert << " --- S4v = " << S4v << endl;
1771    fDataS4[0]=S4v;//0xf028;    fDataS4[0]=S4v;//0xf028;
1772    fDataS4[1]=0xd800;    fDataS4[1]=0xd800;
1773    fDataS4[2]=0x0300;    fDataS4[2]=0x0300;
1774    //  cout<<"  PMT  "<<NdFT<<"  "<<NdF<<endl;    //cout<<"  PMT  "<<NdFT<<"  "<<NdF<<endl;
1775    //cin>>ciao;    //cin>>ciao;
1776  }  }
1777    
# Line 1733  void Digitizer::DigitizeND(){ Line 1800  void Digitizer::DigitizeND(){
1800        NdN++;        NdN++;
1801      }      }
1802    }    }
1803    NdN=100;    //NdN=100; //only for debug
1804    
1805    for(i=0;i<3;i++){    for(i=0;i<3;i++){
1806      fDataND[2*i]=0x0000;      fDataND[2*i]=0x0000;
1807      fDataND[2*i+1]=0x010F;      fDataND[2*i+1]=0x010F;
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.3
  ViewVC Help
Powered by ViewVC 1.1.23