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

Diff of /PamelaDigitizer/DigitizeAC.cxx

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

revision 1.1 by pamelats, Wed May 21 09:50:38 2008 UTC revision 1.5 by mocchiut, Fri Nov 20 10:20:20 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::DigitizeAC() {  void Digitizer::DigitizeAC() {
4    // created:  J. Conrad, KTH    // created:  J. Conrad, KTH
5    // modified: S. Orsi, INFN Roma2    // modified: S. Orsi, INFN Roma2
6    // fDataAC[0-63]:   main AC board    // fDataAC[0-63]:   main AC board
7    // fDataAC[64-127]: extra AC board (identical to main board, for now)    // fDataAC[64-127]: extra AC board (identical to main board, for now)
   
8    // We activate all branches. Once the digitization algorithm is determined    // We activate all branches. Once the digitization algorithm is determined
9    // only the branches that involve needed information will be activated    // only the branches that involve needed information will be activated
     
 //   fhBookTree->SetBranchStatus("Ievnt",&Ievnt);//modified by E.Vannuccini 03/08  
 //   fhBookTree->SetBranchStatus("Nthcat",1);  
 //   fhBookTree->SetBranchStatus("Iparcat",1);  
 //   fhBookTree->SetBranchStatus("Icat",1);  
 //   fhBookTree->SetBranchStatus("Xincat",1);  
 //   fhBookTree->SetBranchStatus("Yincat",1);  
 //   fhBookTree->SetBranchStatus("Zincat",1);  
 //   fhBookTree->SetBranchStatus("Xoutcat",1);  
 //   fhBookTree->SetBranchStatus("Youtcat",1);  
 //   fhBookTree->SetBranchStatus("Zoutcat",1);  
 //   fhBookTree->SetBranchStatus("Erelcat",1);  
 //   fhBookTree->SetBranchStatus("Timecat",1);  
 //   fhBookTree->SetBranchStatus("Pathcat",1);  
 //   fhBookTree->SetBranchStatus("P0cat",1);  
 //   fhBookTree->SetBranchStatus("Nthcas",1);  
 //   fhBookTree->SetBranchStatus("Iparcas",1);  
 //   fhBookTree->SetBranchStatus("Icas",1);  
 //   fhBookTree->SetBranchStatus("Xincas",1);  
 //   fhBookTree->SetBranchStatus("Yincas",1);  
 //   fhBookTree->SetBranchStatus("Zincas",1);  
 //   fhBookTree->SetBranchStatus("Xoutcas",1);  
 //   fhBookTree->SetBranchStatus("Youtcas",1);  
 //   fhBookTree->SetBranchStatus("Zoutcas",1);  
 //   fhBookTree->SetBranchStatus("Erelcas",1);  
 //   fhBookTree->SetBranchStatus("Timecas",1);  
 //   fhBookTree->SetBranchStatus("Pathcas",1);  
 //   fhBookTree->SetBranchStatus("P0cas",1);  
 //   fhBookTree->SetBranchStatus("Nthcard",1);  
 //   fhBookTree->SetBranchStatus("Iparcard",1);  
 //   fhBookTree->SetBranchStatus("Icard",1);  
 //   fhBookTree->SetBranchStatus("Xincard",1);  
 //   fhBookTree->SetBranchStatus("Yincard",1);  
 //   fhBookTree->SetBranchStatus("Zincard",1);  
 //   fhBookTree->SetBranchStatus("Xoutcard",1);  
 //   fhBookTree->SetBranchStatus("Youtcard",1);  
 //   fhBookTree->SetBranchStatus("Zoutcard",1);  
 //   fhBookTree->SetBranchStatus("Erelcard",1);  
 //   fhBookTree->SetBranchStatus("Timecard",1);  
 //   fhBookTree->SetBranchStatus("Pathcard",1);  
 //   fhBookTree->SetBranchStatus("P0card",1);  
10    
11      // Threshold: thr=0.8 MeV.
12      Float_t thr = 8e-4;
13      Int_t nReg=6,i;
14    fDataAC[0] = 0xACAC;    fDataAC[0] = 0xACAC;
15    fDataAC[64]= 0xACAC;    fDataAC[64]= 0xACAC;
16    fDataAC[1] = 0xAC11;      fDataAC[1] = 0xAC11;  
17    fDataAC[65] = 0xAC22;      fDataAC[65] = 0xAC22;  
18      
19    // the third word is a status word (dummy: "no errors are present in the AC boards")    // the third word is a status word (dummy: "no errors are present in the AC boards")
20    fDataAC[2] = 0xFFFF; //FFEF?    fDataAC[2] = 0xFFFF; //FFEF?
21    fDataAC[66] = 0xFFFF;    fDataAC[66] = 0xFFFF;
   
   const UInt_t nReg = 6;  
   
22    // FPGA Registers (dummy)    // FPGA Registers (dummy)
23    for (UInt_t i=0; i<=nReg; i++){    for (i=0; i<=nReg; i++){
24      fDataAC[i+4] = 0xFFFF;      fDataAC[i+4] = 0xFFFF;
25      fDataAC[i+68] = 0xFFFF;      fDataAC[i+68] = 0xFFFF;
26    }    }
27    
   // the last word is a CRC  
   // Dummy for the time being, but it might need to be calculated in the end  
   fDataAC[63] = 0xABCD;  
   fDataAC[127] = 0xABCD;  
   
28    // shift registers (moved to the end of the routine)    // shift registers (moved to the end of the routine)
29    
   //Int_t evntLSB=Ievnt%65536;  
   //Int_t evntMSB=(Int_t)(Ievnt/65536);  
30    Int_t evntLSB=(UShort_t)Ievnt;    Int_t evntLSB=(UShort_t)Ievnt;
31    Int_t evntMSB=Ievnt >> 16;    Int_t evntMSB=Ievnt >> 16;
32    
33    // singles counters are dummy    // singles counters are dummy
34    for (UInt_t i=0; i<=15; i++){  //SO Oct '07: // for (UInt_t i=0; i<=16; i++){    for(i=0; i<=15; i++){  //SO Oct '07:
     //     fDataAC[i+26] = 0x0000;    
     //     fDataAC[i+90] = 0x0000;  
35      fDataAC[i+26] = evntLSB;      fDataAC[i+26] = evntLSB;
36      fDataAC[i+90] = evntLSB;      fDataAC[i+90] = evntLSB;
37    };    }
38        for(i=0; i<=7; i++){
   // coincidences are dummy (increment by 1 at each event)  
   // for (UInt_t i=0; i<=7; i++){  
   //    fDataAC[i+42] = 0x0000;  
   //    fDataAC[i+106] = 0x0000;  
   //   }  
   for (UInt_t i=0; i<=7; i++){  
39      fDataAC[i+42] = evntLSB;      fDataAC[i+42] = evntLSB;
40      fDataAC[i+106] = evntLSB;      fDataAC[i+106] = evntLSB;
41    };    }
   
42    // increments for every trigger might be needed at some point.    // increments for every trigger might be needed at some point.
43    // dummy for now    // dummy for now
44    fDataAC[50]  = 0x0000;    fDataAC[50]  = 0x0000;  
45    fDataAC[114] = 0x0000;    fDataAC[114] = 0x0000;
46    
47    // dummy FPGA clock (increment by 1 at each event)    // dummy FPGA clock (increment by 1 at each event)
48    /*        if(Ievnt<=0xFFFF){
     fDataAC[51] = 0x006C;  
     fDataAC[52] = 0x6C6C;  
     fDataAC[115] = 0x006C;  
     fDataAC[116] = 0x6C6C;  
   */  
   if (Ievnt<=0xFFFF) {  
49      fDataAC[51] = 0x0000;      fDataAC[51] = 0x0000;
50      fDataAC[52] = Ievnt;      fDataAC[52] = Ievnt;
51      fDataAC[115] = 0x0000;      fDataAC[115] = 0x0000;
52      fDataAC[116] = Ievnt;      fDataAC[116] = Ievnt;
53    } else {    }
54      else{
55      fDataAC[51] = evntMSB;      fDataAC[51] = evntMSB;
56      fDataAC[52] = evntLSB;      fDataAC[52] = evntLSB;
57      fDataAC[115] = fDataAC[51];      fDataAC[115] = fDataAC[51];
58      fDataAC[116] = fDataAC[52];      fDataAC[116] = fDataAC[52];
59    }    }
   
60    // dummy temperatures    // dummy temperatures
61    fDataAC[53] = 0x0000;    fDataAC[53] = 0x0000;
62    fDataAC[54] = 0x0000;    fDataAC[54] = 0x0000;
63    fDataAC[117] = 0x0000;    fDataAC[117] = 0x0000;
64    fDataAC[118] = 0x0000;    fDataAC[118] = 0x0000;
   
   
65    // dummy DAC thresholds    // dummy DAC thresholds
66    for (UInt_t i=0; i<=7; i++){    for(i=0; i<=7; i++){
67      fDataAC[i+55] = 0x1A13;        fDataAC[i+55] = 0x1A13;  
68      fDataAC[i+119] = 0x1A13;      fDataAC[i+119] = 0x1A13;
69    }    }
70        // In this simpliefied approach we will assume that once a particle releases > 0.5 mip in one of the 12 AC detectors it will fire. We will furthermore assume that both cards read out identical data.
   // In this simpliefied approach we will assume that once  
   // a particle releases > 0.5 mip in one of the 12 AC detectors it  
   // will fire. We will furthermore assume that both cards read out  
   // identical data.  
   
   // If you develop your digitization algorithm, you should start by  
   // identifying the information present in level2 (post-darth-vader)  
   // data.  
   
   Float_t SumEcat[5];  
   Float_t SumEcas[5];  
   Float_t SumEcard[5];  
   for (Int_t k= 0;k<5;k++){  
     SumEcat[k]=0.;  
     SumEcas[k]=0.;  
     SumEcard[k]=0.;  
   };  
71    
72    if (Nthcat>50 || Nthcas>50 || Nthcard>50)    // If you develop your digitization algorithm, you should start by identifying the information present in level2 (post-darth-vader) data.
     printf("*** ERROR AC! NthAC out of range!\n\n");  
73    
74      Float_t SumEcat[4],SumEcas[4],SumEcard[4];
75      for(i=0;i<4;i++){
76        SumEcat[i]=0.;
77        SumEcas[i]=0.;
78        SumEcard[i]=0.;
79      }
80    // energy dependence on position (see file AcFitOutputDistancePmt.C by S.Orsi)    // energy dependence on position (see file AcFitOutputDistancePmt.C by S.Orsi)
81    // based on J.Lundquist's calculations (PhD thesis, page 94)    // based on J.Lundquist's calculations (PhD thesis, page 94)
82    // function: [0]+[1]*atan([2]/(x+1)), where the 3 parameters are:    // function: [0]+[1]*atan([2]/(x+1)), where the 3 parameters are:
# Line 192  void Digitizer::DigitizeAC() { Line 85  void Digitizer::DigitizeAC() {
85    //   9.81177e+00   +- 1.21284e+00    //   9.81177e+00   +- 1.21284e+00
86    // hp: 1 minimum ionising particle at 35cm from the PMT releases 1mip    // hp: 1 minimum ionising particle at 35cm from the PMT releases 1mip
87    
   TF1 *attenAC = new TF1("fAttAC",".825+.64*atan(9.8/x)",0.,45.);  
   
88    // PMT positions: x,y,z: (average position of the 2 PMTs)    // PMT positions: x,y,z: (average position of the 2 PMTs)
89    Float_t posCasPmt[4][3]={{28.308, -17.168, 63.644},    // 1 - CAS CPU: x,y,z    Float_t posCasPmt[4][3]={{28.308, -17.168, 63.644},    // 1 - CAS CPU: x,y,z
90                             {18.893, 24.913, 63.644},     // 2 - CAS DCDC                             {18.893, 24.913, 63.644},     // 2 - CAS DCDC
91                             {-24.307, 17.162, 63.644},    // 3 - CAS VME                               {-24.307, 17.162, 63.644},    // 3 - CAS VME  
92                             {-17.765, -28.300, 63.644}};  // 4 - CAS IPM                               {-17.765, -28.300, 63.644}};  // 4 - CAS IPM  
     
93    Float_t dAC=0.; // distance from PMT    Float_t dAC=0.; // distance from PMT
     
   // look in CAT  
   //  for (UInt_t k= 0;k<50;k++){  
   for (Int_t k= 0;k<Nthcat;k++){  
     if (Erelcat[k] > 0)  
       SumEcat[Icat[k]] += Erelcat[k];  
   };  
   
   // look in CAS  
   for (Int_t k= 0;k<Nthcas;k++){  
     if (Erelcas[k] >0) {  
       dAC=sqrt(pow((Xincas[k]+Xoutcas[k])/2 - posCasPmt[Icas[k]-1][0],2) + pow((Yincas[k]+Youtcas[k])/2 - posCasPmt[Icas[k]-1][1],2) + pow((Zincas[k]+Zoutcas[k])/2 - posCasPmt[Icas[k]-1][2],2));  
       SumEcas[Icas[k]] += Erelcas[k]*attenAC->Eval(dAC);  
     }  
   };  
94    
95    // look in CARD    if(Nthcat>*ncat){
96    for (Int_t k= 0;k<Nthcard;k++){      cout<<"*** ERROR AC! Nthcat=  "<<Nthcat<<" out of range! "<<endl;
97      if (Erelcard[k] >0)      for(i=0;i<4;i++)SumEcat[i]=2*thr;
98        SumEcard[Icard[k]] += Erelcard[k];    }
99    };    else{
100        for(i=0;i<Nthcat;i++){
101          if(Icat[i]>0 && Icat[i]<5)SumEcat[Icat[i]-1]+=Erelcat[i];
102        }
103      }
104      if(Nthcas>*ncas){
105        cout<<"*** ERROR AC! Nthcas=  "<<Nthcas<<" out of range!"<<endl;
106        for(i=0;i<4;i++)SumEcas[i]=2*thr;
107      }
108      else{
109        for (i=0;i<Nthcas;i++){
110          if(Icas[i]>0 && Icas[i]<5){
111            dAC=sqrt(pow((Xincas[i]+Xoutcas[i])/2 - posCasPmt[Icas[i]-1][0],2) + pow((Yincas[i]+Youtcas[i])/2 - posCasPmt[Icas[i]-1][1],2) + pow((Zincas[i]+Zoutcas[i])/2 - posCasPmt[Icas[i]-1][2],2));
112            SumEcas[Icas[i]-1] += Erelcas[i]*attenAC->Eval(dAC);
113          }
114        }
115      }
116      if(Nthcard>*ncar){
117        cout<<"*** ERROR AC! Nthcard= "<<Nthcard<<" out of range!"<<endl;
118        for(i=0;i<4;i++)SumEcard[i]=2*thr;
119      }
120      else{
121        for(Int_t k= 0;k<Nthcard;k++){
122          if(Icard[i]>0 && Icard[i]<5)SumEcard[Icard[k]-1] += Erelcard[k];
123        }
124      }
125    
126    // channel mapping              Hit Map    // channel mapping              Hit Map
127    // 1 CARD4                      0          LSB    // 1 CARD4                      0          LSB
# Line 243  void Digitizer::DigitizeAC() { Line 143  void Digitizer::DigitizeAC() {
143    
144    // In the first version only the hit-map is filled, not the SR.    // In the first version only the hit-map is filled, not the SR.
145    
   // Threshold: 0.8 MeV.  
   
   Float_t thr = 8e-4;  
   
146    fDataAC[3] = 0x0000;    fDataAC[3] = 0x0000;
147    
148    if (SumEcas[0] > thr)  fDataAC[3]  = 0x0004;    if (SumEcas[0] > thr)  fDataAC[3]  = 0x0004;
# Line 279  void Digitizer::DigitizeAC() { Line 175  void Digitizer::DigitizeAC() {
175    //    printf("%0x  ",fDataAC[i]);      //    printf("%0x  ",fDataAC[i]);  
176    //    if ((i+1)%8 ==0) cout << endl;    //    if ((i+1)%8 ==0) cout << endl;
177    //   }    //   }
178    
179      // the last word is a CRC
180      fDataAC[63] = EvaluateCrcAC (fDataAC, 0);
181      fDataAC[127]= EvaluateCrcAC (fDataAC, 1);
182    };
183    
184    
185    
186    
187    
188    UShort_t Digitizer::EvaluateCrcAC( UShort_t *pAC , Bool_t flag )
189    {
190      // Flag=0for AC1 ; flag=1 for AC2
191      UShort_t check=0;
192      
193      UInt_t l=0;
194      for (UInt_t i=0; i<63; i++)
195        {
196          l = 64*flag + i ;
197          //printf("ACAC  i=%d   l=%d   chech=%d \n",i,l,check);
198          check = crcAC(check, pAC[l]);
199        }
200      return check;
201      
202    };
203    
204    
205    UShort_t Digitizer::crcAC( UShort_t old , UShort_t neww )
206    {
207      union crc_data {
208        short word;
209        struct bit_field {
210          unsigned b0:1;
211          unsigned b1:1;
212          unsigned b2:1;
213          unsigned b3:1;
214          unsigned b4:1;
215          unsigned b5:1;
216          unsigned b6:1;
217          unsigned b7:1;
218          unsigned b8:1;
219          unsigned b9:1;
220          unsigned b10:1;
221          unsigned b11:1;
222          unsigned b12:1;
223          unsigned b13:1;
224          unsigned b14:1;
225          unsigned b15:1;
226        } bit;
227      } c,d,r;
228    
229      c.word = old;
230      d.word = neww;
231      r.word = 0;
232    
233      r.bit.b0 = c.bit.b0 ^ c.bit.b4 ^ c.bit.b8 ^ c.bit.b11 ^ c.bit.b12 ^
234                 d.bit.b0 ^ d.bit.b4 ^ d.bit.b8 ^ d.bit.b11 ^ d.bit.b12;
235    
236      r.bit.b1 = c.bit.b1 ^ c.bit.b5 ^ c.bit.b9 ^ c.bit.b12 ^ c.bit.b13 ^
237                 d.bit.b1 ^ d.bit.b5 ^ d.bit.b9 ^ d.bit.b12 ^ d.bit.b13;
238    
239      r.bit.b2 = c.bit.b2 ^ c.bit.b6 ^ c.bit.b10 ^ c.bit.b13 ^ c.bit.b14 ^
240                 d.bit.b2 ^ d.bit.b6 ^ d.bit.b10 ^ d.bit.b13 ^ d.bit.b14;
241    
242      r.bit.b3 = c.bit.b3 ^ c.bit.b7 ^ c.bit.b11 ^ c.bit.b14 ^ c.bit.b15 ^
243                 d.bit.b3 ^ d.bit.b7 ^ d.bit.b11 ^ d.bit.b14 ^ d.bit.b15;
244    
245      r.bit.b4 = c.bit.b4 ^ c.bit.b8 ^ c.bit.b12 ^ c.bit.b15 ^
246                 d.bit.b4 ^ d.bit.b8 ^ d.bit.b12 ^ d.bit.b15;
247    
248      r.bit.b5 = c.bit.b0 ^ c.bit.b4 ^ c.bit.b5 ^ c.bit.b8 ^ c.bit.b9 ^
249                 c.bit.b11 ^ c.bit.b12 ^ c.bit.b13 ^
250                 d.bit.b0 ^ d.bit.b4 ^ d.bit.b5 ^ d.bit.b8 ^ d.bit.b9 ^
251                 d.bit.b11 ^ d.bit.b12 ^ d.bit.b13;
252    
253      r.bit.b6 = c.bit.b1 ^ c.bit.b5 ^ c.bit.b6 ^ c.bit.b9 ^ c.bit.b10 ^
254                 c.bit.b12 ^ c.bit.b13 ^ c.bit.b14 ^
255                 d.bit.b1 ^ d.bit.b5 ^ d.bit.b6 ^ d.bit.b9 ^ d.bit.b10 ^
256                 d.bit.b12 ^ d.bit.b13 ^ d.bit.b14;
257    
258      r.bit.b7 = c.bit.b2 ^ c.bit.b6 ^ c.bit.b7 ^ c.bit.b10 ^ c.bit.b11 ^
259                 c.bit.b13 ^ c.bit.b14 ^ c.bit.b15 ^
260                 d.bit.b2 ^ d.bit.b6 ^ d.bit.b7 ^ d.bit.b10 ^ d.bit.b11 ^
261                 d.bit.b13 ^ d.bit.b14 ^ d.bit.b15;
262    
263      r.bit.b8 = c.bit.b3 ^ c.bit.b7 ^ c.bit.b8 ^ c.bit.b11 ^ c.bit.b12 ^
264                 c.bit.b14 ^ c.bit.b15 ^
265                 d.bit.b3 ^ d.bit.b7 ^ d.bit.b8 ^ d.bit.b11 ^ d.bit.b12 ^
266                 d.bit.b14 ^ d.bit.b15;
267    
268      r.bit.b9 = c.bit.b4 ^ c.bit.b8 ^ c.bit.b9 ^ c.bit.b12 ^ c.bit.b13 ^
269                 c.bit.b15 ^
270                 d.bit.b4 ^ d.bit.b8 ^ d.bit.b9 ^ d.bit.b12 ^ d.bit.b13 ^
271                 d.bit.b15;
272    
273      r.bit.b10 = c.bit.b5 ^ c.bit.b9 ^ c.bit.b10 ^ c.bit.b13 ^ c.bit.b14 ^
274                  d.bit.b5 ^ d.bit.b9 ^ d.bit.b10 ^ d.bit.b13 ^ d.bit.b14;
275    
276      r.bit.b11 = c.bit.b6 ^ c.bit.b10 ^ c.bit.b11 ^ c.bit.b14 ^ c.bit.b15 ^
277                  d.bit.b6 ^ d.bit.b10 ^ d.bit.b11 ^ d.bit.b14 ^ d.bit.b15;
278    
279      r.bit.b12 = c.bit.b0 ^ c.bit.b4 ^ c.bit.b7 ^ c.bit.b8 ^ c.bit.b15 ^
280                  d.bit.b0 ^ d.bit.b4 ^ d.bit.b7 ^ d.bit.b8 ^ d.bit.b15;
281    
282      r.bit.b13 = c.bit.b1 ^ c.bit.b5 ^ c.bit.b8 ^ c.bit.b9 ^
283                  d.bit.b1 ^ d.bit.b5 ^ d.bit.b8 ^ d.bit.b9;
284    
285      r.bit.b14 = c.bit.b2 ^ c.bit.b6 ^ c.bit.b9 ^ c.bit.b10 ^
286                  d.bit.b2 ^ d.bit.b6 ^ d.bit.b9 ^ d.bit.b10;
287    
288      r.bit.b15 = c.bit.b3 ^ c.bit.b7 ^ c.bit.b10 ^ c.bit.b11 ^
289                  d.bit.b3 ^ d.bit.b7 ^ d.bit.b10 ^ d.bit.b11;
290    
291      return r.word;
292    
293  };  };
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.5
  ViewVC Help
Powered by ViewVC 1.1.23