/[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.3 by pamelats, Wed Oct 15 14:03:15 2008 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    // the last word is a CRC
28    // Dummy for the time being, but it might need to be calculated in the end    // Dummy for the time being, but it might need to be calculated in the end
29    fDataAC[63] = 0xABCD;    fDataAC[63] = 0xABCD;
# Line 103  void Digitizer::DigitizeAC() { Line 31  void Digitizer::DigitizeAC() {
31    
32    // shift registers (moved to the end of the routine)    // shift registers (moved to the end of the routine)
33    
   //Int_t evntLSB=Ievnt%65536;  
   //Int_t evntMSB=(Int_t)(Ievnt/65536);  
34    Int_t evntLSB=(UShort_t)Ievnt;    Int_t evntLSB=(UShort_t)Ievnt;
35    Int_t evntMSB=Ievnt >> 16;    Int_t evntMSB=Ievnt >> 16;
36    
37    // singles counters are dummy    // singles counters are dummy
38    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;  
39      fDataAC[i+26] = evntLSB;      fDataAC[i+26] = evntLSB;
40      fDataAC[i+90] = evntLSB;      fDataAC[i+90] = evntLSB;
41    };    }
42        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++){  
43      fDataAC[i+42] = evntLSB;      fDataAC[i+42] = evntLSB;
44      fDataAC[i+106] = evntLSB;      fDataAC[i+106] = evntLSB;
45    };    }
   
46    // increments for every trigger might be needed at some point.    // increments for every trigger might be needed at some point.
47    // dummy for now    // dummy for now
48    fDataAC[50]  = 0x0000;    fDataAC[50]  = 0x0000;
49    fDataAC[114] = 0x0000;    fDataAC[114] = 0x0000;
50    
51    // dummy FPGA clock (increment by 1 at each event)    // dummy FPGA clock (increment by 1 at each event)
52    /*        if(Ievnt<=0xFFFF){
     fDataAC[51] = 0x006C;  
     fDataAC[52] = 0x6C6C;  
     fDataAC[115] = 0x006C;  
     fDataAC[116] = 0x6C6C;  
   */  
   if (Ievnt<=0xFFFF) {  
53      fDataAC[51] = 0x0000;      fDataAC[51] = 0x0000;
54      fDataAC[52] = Ievnt;      fDataAC[52] = Ievnt;
55      fDataAC[115] = 0x0000;      fDataAC[115] = 0x0000;
56      fDataAC[116] = Ievnt;      fDataAC[116] = Ievnt;
57    } else {    }
58      else{
59      fDataAC[51] = evntMSB;      fDataAC[51] = evntMSB;
60      fDataAC[52] = evntLSB;      fDataAC[52] = evntLSB;
61      fDataAC[115] = fDataAC[51];      fDataAC[115] = fDataAC[51];
62      fDataAC[116] = fDataAC[52];      fDataAC[116] = fDataAC[52];
63    }    }
   
64    // dummy temperatures    // dummy temperatures
65    fDataAC[53] = 0x0000;    fDataAC[53] = 0x0000;
66    fDataAC[54] = 0x0000;    fDataAC[54] = 0x0000;
67    fDataAC[117] = 0x0000;    fDataAC[117] = 0x0000;
68    fDataAC[118] = 0x0000;    fDataAC[118] = 0x0000;
   
   
69    // dummy DAC thresholds    // dummy DAC thresholds
70    for (UInt_t i=0; i<=7; i++){    for(i=0; i<=7; i++){
71      fDataAC[i+55] = 0x1A13;        fDataAC[i+55] = 0x1A13;  
72      fDataAC[i+119] = 0x1A13;      fDataAC[i+119] = 0x1A13;
73    }    }
74        // 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.;  
   };  
75    
76    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");  
77    
78      Float_t SumEcat[4],SumEcas[4],SumEcard[4];
79      for(i=0;i<4;i++){
80        SumEcat[i]=0.;
81        SumEcas[i]=0.;
82        SumEcard[i]=0.;
83      }
84    // energy dependence on position (see file AcFitOutputDistancePmt.C by S.Orsi)    // energy dependence on position (see file AcFitOutputDistancePmt.C by S.Orsi)
85    // based on J.Lundquist's calculations (PhD thesis, page 94)    // based on J.Lundquist's calculations (PhD thesis, page 94)
86    // 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 89  void Digitizer::DigitizeAC() {
89    //   9.81177e+00   +- 1.21284e+00    //   9.81177e+00   +- 1.21284e+00
90    // hp: 1 minimum ionising particle at 35cm from the PMT releases 1mip    // hp: 1 minimum ionising particle at 35cm from the PMT releases 1mip
91    
   TF1 *attenAC = new TF1("fAttAC",".825+.64*atan(9.8/x)",0.,45.);  
   
92    // PMT positions: x,y,z: (average position of the 2 PMTs)    // PMT positions: x,y,z: (average position of the 2 PMTs)
93    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
94                             {18.893, 24.913, 63.644},     // 2 - CAS DCDC                             {18.893, 24.913, 63.644},     // 2 - CAS DCDC
95                             {-24.307, 17.162, 63.644},    // 3 - CAS VME                               {-24.307, 17.162, 63.644},    // 3 - CAS VME  
96                             {-17.765, -28.300, 63.644}};  // 4 - CAS IPM                               {-17.765, -28.300, 63.644}};  // 4 - CAS IPM  
     
97    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);  
     }  
   };  
98    
99    // look in CARD    if(Nthcat>*ncat){
100    for (Int_t k= 0;k<Nthcard;k++){      cout<<"*** ERROR AC! Nthcat=  "<<Nthcat<<" out of range! "<<endl;
101      if (Erelcard[k] >0)      for(i=0;i<4;i++)SumEcat[i]=2*thr;
102        SumEcard[Icard[k]] += Erelcard[k];    }
103    };    else{
104        for(i=0;i<Nthcat;i++){
105          if(Icat[i]>0 && Icat[i]<5)SumEcat[Icat[i]-1]+=Erelcat[i];
106        }
107      }
108      if(Nthcas>*ncas){
109        cout<<"*** ERROR AC! Nthcas=  "<<Nthcas<<" out of range!"<<endl;
110        for(i=0;i<4;i++)SumEcas[i]=2*thr;
111      }
112      else{
113        for (i=0;i<Nthcas;i++){
114          if(Icas[i]>0 && Icas[i]<5){
115            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));
116            SumEcas[Icas[i]-1] += Erelcas[i]*attenAC->Eval(dAC);
117          }
118        }
119      }
120      if(Nthcard>*ncar){
121        cout<<"*** ERROR AC! Nthcard= "<<Nthcard<<" out of range!"<<endl;
122        for(i=0;i<4;i++)SumEcard[i]=2*thr;
123      }
124      else{
125        for(Int_t k= 0;k<Nthcard;k++){
126          if(Icard[i]>0 && Icard[i]<5)SumEcard[Icard[k]-1] += Erelcard[k];
127        }
128      }
129    
130    // channel mapping              Hit Map    // channel mapping              Hit Map
131    // 1 CARD4                      0          LSB    // 1 CARD4                      0          LSB
# Line 243  void Digitizer::DigitizeAC() { Line 147  void Digitizer::DigitizeAC() {
147    
148    // 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.
149    
   // Threshold: 0.8 MeV.  
   
   Float_t thr = 8e-4;  
   
150    fDataAC[3] = 0x0000;    fDataAC[3] = 0x0000;
151    
152    if (SumEcas[0] > thr)  fDataAC[3]  = 0x0004;    if (SumEcas[0] > thr)  fDataAC[3]  = 0x0004;
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.3
  ViewVC Help
Powered by ViewVC 1.1.23