ac/src/PamVMCAcDig.cxx

#include "PamVMCAcDig.h"
#include <TF1.h>
#include <TMath.h>

using TMath::Power;
using TMath::Sqrt;

ClassImp(PamVMCAcDig)

void PamVMCAcDig::DigitizeAC(Int_t EventNo){

  // cout<<"Digitizing AC... "<<endl;

  fDEBUG = kTRUE;
  
  for (Int_t k= 0;k<5;k++) fSumEcat[k]=fSumEcas[k]=fSumEcard[k]=0.;
  //Zeroing arrays holds energy releases in AC

  UShort_t DataAC[128];
  for(Int_t i=0; i<128; i++) DataAC[i]=0;
    

  // created:  J. Conrad, KTH
  // modified: S. Orsi, INFN Roma2
  // DataAC[0-63]:   main AC board
  // DataAC[64-127]: extra AC board (identical to main board, for now)

  DataAC[0] = 0xACAC;
  DataAC[64]= 0xACAC;
  DataAC[1] = 0xAC11;   
  DataAC[65] = 0xAC22;  

  // the third word is a status word 
  //(dummy: "no errors are present in the AC boards")
  DataAC[2] = 0xFFFF; //FFEF?
  DataAC[66] = 0xFFFF;

  const UInt_t nReg = 6;

  // FPGA Registers (dummy)
  for (UInt_t i=0; i<=nReg; i++){
    DataAC[i+4] = 0xFFFF;
    DataAC[i+68] = 0xFFFF;
  }

  // the last word is a CRC
  // Dummy for the time being, but it might need 
  // to be calculated in the end
  DataAC[63] = 0xABCD;
  DataAC[127] = 0xABCD;

  // shift registers (moved to the end of the routine)

  Int_t evntLSB=(UShort_t)EventNo;
  Int_t evntMSB=EventNo >> 16;

  // singles counters are dummy
  for (UInt_t i=0; i<=15; i++){  
    DataAC[i+26] = evntLSB;
    DataAC[i+90] = evntLSB;
  };
  
  for (UInt_t i=0; i<=7; i++){
    DataAC[i+42] = evntLSB;
    DataAC[i+106] = evntLSB;
  };


  // increments for every trigger might be needed at some point.
  // dummy for now
  DataAC[50]  = 0x0000;
  DataAC[114] = 0x0000;

  // dummy FPGA clock (increment by 1 at each event)
  if (EventNo<=0xFFFF) {
    DataAC[51] = 0x0000;
    DataAC[52] = EventNo;
    DataAC[115] = 0x0000;
    DataAC[116] = EventNo;
  } else {
    DataAC[51] = evntMSB;
    DataAC[52] = evntLSB;
    DataAC[115] = DataAC[51];
    DataAC[116] = DataAC[52];
  }

  // dummy temperatures
  DataAC[53] = 0x0000;
  DataAC[54] = 0x0000;
  DataAC[117] = 0x0000;
  DataAC[118] = 0x0000;


  // dummy DAC thresholds
  for (UInt_t i=0; i<=7; i++){
    DataAC[i+55] = 0x1A13;  
    DataAC[i+119] = 0x1A13;
  }

  FillErel();
  
  // 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. 

  // channel mapping              Hit Map
  // 1 CARD4                      0          LSB
  // 2 CAT2                       0
  // 3 CAS1                       0
  // 4 NC                         0
  // 5 CARD2                      0
  // 6 CAT4                       1 
  // 7 CAS4                       0  
  // 8 NC                         0
  // 9 CARD3                      0
  // 10 CAT3                      0
  // 11 CAS3                      0 
  // 12 NC                        0
  // 13 CARD1                     0
  // 14 CAT1                      0
  // 15 CAS2                      0
  // 16 NC                        0          MSB

  // In the first version only the hit-map is filled, not the SR.

  // Threshold: 0.8 MeV.

  const Float_t thr = 8e-4; //(GeV)

  //Important Note: In simulation model we have
  //CAT as solid plasic plate. So, I set values:

  fSumEcat[1]=fSumEcat[2]=fSumEcat[3]=fSumEcat[0];
  
  DataAC[3] = 0x0000;

  if (fSumEcas[0] > thr)  DataAC[3]  = 0x0004;
  if (fSumEcas[1] > thr)  DataAC[3] += 0x4000;
  if (fSumEcas[2] > thr)  DataAC[3] += 0x0400;
  if (fSumEcas[3] > thr)  DataAC[3] += 0x0040;   

  if (fSumEcat[0] > thr)  DataAC[3] += 0x2000; //We have solid CAT scintillator
  if (fSumEcat[1] > thr)  DataAC[3] += 0x0002; //See note above
  if (fSumEcat[2] > thr)  DataAC[3] += 0x0200; //
  if (fSumEcat[3] > thr)  DataAC[3] += 0x0020; //
  
  if (fSumEcard[0] > thr) DataAC[3] += 0x1000; //C1D1(1)
  if (fSumEcard[1] > thr) DataAC[3] += 0x0010; //C1D1(2)
  if (fSumEcard[2] > thr) DataAC[3] += 0x0100; //C2D1(1)
  if (fSumEcard[3] > thr) DataAC[3] += 0x0001; //C2D1(2)

  DataAC[67] = DataAC[3];

  // shift registers
  // the central bin is equal to the hitmap, all 
  // other bins in the shift register are 0
  for (UInt_t i=0; i<=15; i++){
    DataAC[i+11] = 0x0000;   
    DataAC[i+75] = 0x0000;
  }
  DataAC[18] = DataAC[3];
  DataAC[82] = DataAC[3];

  //++++++ WRITE EVERYTHING TO DIGITIZER'S BUFFER +++//
  fData.clear(); //clearing data buffer

  for(Int_t i= 0; i<128; i++) fData.push_back(DataAC[i]);


  if (fDEBUG){
    cout<<"====== AC hit info ======"<<endl;
    cout<<"CAT:   Erel="<<fSumEcat[0]*1.e6<<" (keV)"<<endl
        <<"CARD1: Erel="<<fSumEcard[0]*1.e6<<" (keV)"<<endl
        <<"CARD2: Erel="<<fSumEcard[1]*1.e6<<" (keV)"<<endl
        <<"CARD3: Erel="<<fSumEcard[2]*1.e6<<" (keV)"<<endl
        <<"CARD4: Erel="<<fSumEcard[3]*1.e6<<" (keV)"<<endl
        <<"CAS1:  Erel="<<fSumEcas[0]*1.e6<<" (keV)"<<endl
        <<"CAS2:  Erel="<<fSumEcas[1]*1.e6<<" (keV)"<<endl
        <<"CAS3:  Erel="<<fSumEcas[2]*1.e6<<" (keV)"<<endl
        <<"CAS4:  Erel="<<fSumEcas[3]*1.e6<<" (keV)"<<endl;
    cout<<"========================="<<endl;

  }
}

void PamVMCAcDig::FillErel(){
  /*
    Detector     KEY          Number of paddles

    CAT          TOP1                  1 (solid volume with id=1)
    CAS          SID1                  4 (1-4)
    CARD         C1D1 & C2D1           2 & 2 (1-2)
  */
  TClonesArray* hc=0;
  const char* keyac[4] = {"TOP1","SID1","C1D1","C2D1"};
  for(Int_t i=0; i<4; i++){
    hc = (TClonesArray*)fhitscolmap.GetValue(keyac[i]);
    if(hc) AddErel(hc, i);
    hc = 0;
  }
}

void PamVMCAcDig::AddErel(TClonesArray* HitColl, Int_t detcode){


  // energy dependence on position (see file AcFitOutputDistancePmt.C by S.Orsi)
  // based on J.Lundquist's calculations (PhD thesis, page 94)
  // function: [0]+[1]*atan([2]/(x+1)), where the 3 parameters are:
  //   8.25470e-01   +- 1.79489e-02
  //   6.41609e-01   +- 2.65846e-02
  //   9.81177e+00   +- 1.21284e+00
  // hp: 1 minimum ionising particle at 35cm from the PMT releases 1mip

  // TF1 *attenAC = new TF1("fAttAC",".825+.64*atan(9.8/x)",0.,45.);

  const Float_t posCasPmt[4][3]={{28.308, -17.168, 63.644},// 1 - CAS CPU: x,y,z
                                {18.893, 24.913, 63.644},// 2 - CAS DCDC 
                                {-24.307, 17.162, 63.644},// 3 - CAS VME  
                                {-17.765, -28.300, 63.644}};// 4 - CAS IPM

  Float_t dx,dy,dz,dAC; // distance from PMT
  Int_t hitsNo=HitColl->GetEntriesFast();
  PamVMCDetectorHit * hit = 0;

  if (hitsNo>50) cout<<"WARNING: PamVMCAcDig::AddErel... "
                     <<"Number of hits in AC is out of range!!!"<<endl;

  switch(detcode){
  case 0: //CAT
    for(Int_t k=0; k<hitsNo; k++){
      hit = (PamVMCDetectorHit*)HitColl->At(k);
      fSumEcat[(hit->GetPOS())-1] += hit->GetEREL();
    }
    break;

  case 1: //CAS
    for(Int_t k=0; k<hitsNo; k++){
      hit = (PamVMCDetectorHit*)HitColl->At(k);
      dx = (hit->GetXIN()+hit->GetXOUT())/2.-posCasPmt[hit->GetPOS()-1][0];
      dy = (hit->GetYIN()+hit->GetYOUT())/2.-posCasPmt[hit->GetPOS()-1][1];
      dz = (hit->GetZIN()+hit->GetZOUT())/2.-posCasPmt[hit->GetPOS()-1][2];
      dAC = Sqrt(dx*dx + dy*dy + dz*dz);
      //fSumEcas[(hit->GetPOS())-1] += (hit->GetEREL())*attenAC->Eval(dAC);
      if ((hit->GetPOS()==3)) fSumEcas[3] += (hit->GetEREL())*fattenAC->Eval(dAC);
      if ((hit->GetPOS()==4)) fSumEcas[1] += (hit->GetEREL())*fattenAC->Eval(dAC);
      if ((hit->GetPOS()==2)) fSumEcas[0] += (hit->GetEREL())*fattenAC->Eval(dAC);
      if ((hit->GetPOS()==1)) fSumEcas[2] += (hit->GetEREL())*fattenAC->Eval(dAC);
    }
    break;

  case 2: //C1D1
    for(Int_t k=0; k<hitsNo; k++){
      hit = (PamVMCDetectorHit*)HitColl->At(k);
      if ((hit->GetPOS()==1)) fSumEcard[3] += hit->GetEREL();
      if ((hit->GetPOS()==2)) fSumEcard[0] += hit->GetEREL(); 
      //fSumEcard[(hit->GetPOS())-1] += hit->GetEREL();
    }
    break;
  case 3: //C2D1
    for(Int_t k=0; k<hitsNo; k++){
      hit = (PamVMCDetectorHit*)HitColl->At(k);
      if ((hit->GetPOS()==1)) fSumEcard[2] += hit->GetEREL();
      if ((hit->GetPOS()==2)) fSumEcard[1] += hit->GetEREL(); 
      //fSumEcard[(hit->GetPOS())+1] += hit->GetEREL();
    }
    break;

  default:
    break;
  }


}
1	pamelats	1.1	#include "PamVMCAcDig.h"
2			#include <TF1.h>
3			#include <TMath.h>
4
5			using TMath::Power;
6			using TMath::Sqrt;
7
8			ClassImp(PamVMCAcDig)
9
10			void PamVMCAcDig::DigitizeAC(Int_t EventNo){
11
12			// cout<<"Digitizing AC... "<<endl;
13
14			fDEBUG = kTRUE;
15
16			for (Int_t k= 0;k<5;k++) fSumEcat[k]=fSumEcas[k]=fSumEcard[k]=0.;
17			//Zeroing arrays holds energy releases in AC
18
19			UShort_t DataAC[128];
20			for(Int_t i=0; i<128; i++) DataAC[i]=0;
21
22
23			// created: J. Conrad, KTH
24			// modified: S. Orsi, INFN Roma2
25			// DataAC[0-63]: main AC board
26			// DataAC[64-127]: extra AC board (identical to main board, for now)
27
28			DataAC[0] = 0xACAC;
29			DataAC[64]= 0xACAC;
30			DataAC[1] = 0xAC11;
31			DataAC[65] = 0xAC22;
32
33			// the third word is a status word
34			//(dummy: "no errors are present in the AC boards")
35			DataAC[2] = 0xFFFF; //FFEF?
36			DataAC[66] = 0xFFFF;
37
38			const UInt_t nReg = 6;
39
40			// FPGA Registers (dummy)
41			for (UInt_t i=0; i<=nReg; i++){
42			DataAC[i+4] = 0xFFFF;
43			DataAC[i+68] = 0xFFFF;
44			}
45
46			// the last word is a CRC
47			// Dummy for the time being, but it might need
48			// to be calculated in the end
49			DataAC[63] = 0xABCD;
50			DataAC[127] = 0xABCD;
51
52			// shift registers (moved to the end of the routine)
53
54			Int_t evntLSB=(UShort_t)EventNo;
55			Int_t evntMSB=EventNo >> 16;
56
57			// singles counters are dummy
58			for (UInt_t i=0; i<=15; i++){
59			DataAC[i+26] = evntLSB;
60			DataAC[i+90] = evntLSB;
61			};
62
63			for (UInt_t i=0; i<=7; i++){
64			DataAC[i+42] = evntLSB;
65			DataAC[i+106] = evntLSB;
66			};
67
68
69			// increments for every trigger might be needed at some point.
70			// dummy for now
71			DataAC[50] = 0x0000;
72			DataAC[114] = 0x0000;
73
74			// dummy FPGA clock (increment by 1 at each event)
75			if (EventNo<=0xFFFF) {
76			DataAC[51] = 0x0000;
77			DataAC[52] = EventNo;
78			DataAC[115] = 0x0000;
79			DataAC[116] = EventNo;
80			} else {
81			DataAC[51] = evntMSB;
82			DataAC[52] = evntLSB;
83			DataAC[115] = DataAC[51];
84			DataAC[116] = DataAC[52];
85			}
86
87			// dummy temperatures
88			DataAC[53] = 0x0000;
89			DataAC[54] = 0x0000;
90			DataAC[117] = 0x0000;
91			DataAC[118] = 0x0000;
92
93
94			// dummy DAC thresholds
95			for (UInt_t i=0; i<=7; i++){
96			DataAC[i+55] = 0x1A13;
97			DataAC[i+119] = 0x1A13;
98			}
99
100			FillErel();
101
102			// In this simpliefied approach we will assume that once
103			// a particle releases > 0.5 mip in one of the 12 AC detectors it
104			// will fire. We will furthermore assume that both cards read out
105			// identical data.
106
107			// If you develop your digitization algorithm, you should start by
108			// identifying the information present in level2 (post-darth-vader)
109			// data.
110
111			// channel mapping Hit Map
112			// 1 CARD4 0 LSB
113			// 2 CAT2 0
114			// 3 CAS1 0
115			// 4 NC 0
116			// 5 CARD2 0
117			// 6 CAT4 1
118			// 7 CAS4 0
119			// 8 NC 0
120			// 9 CARD3 0
121			// 10 CAT3 0
122			// 11 CAS3 0
123			// 12 NC 0
124			// 13 CARD1 0
125			// 14 CAT1 0
126			// 15 CAS2 0
127			// 16 NC 0 MSB
128
129			// In the first version only the hit-map is filled, not the SR.
130
131			// Threshold: 0.8 MeV.
132
133			const Float_t thr = 8e-4; //(GeV)
134
135			//Important Note: In simulation model we have
136			//CAT as solid plasic plate. So, I set values:
137
138			fSumEcat[1]=fSumEcat[2]=fSumEcat[3]=fSumEcat[0];
139
140			DataAC[3] = 0x0000;
141
142			if (fSumEcas[0] > thr) DataAC[3] = 0x0004;
143			if (fSumEcas[1] > thr) DataAC[3] += 0x4000;
144			if (fSumEcas[2] > thr) DataAC[3] += 0x0400;
145			if (fSumEcas[3] > thr) DataAC[3] += 0x0040;
146
147			if (fSumEcat[0] > thr) DataAC[3] += 0x2000; //We have solid CAT scintillator
148			if (fSumEcat[1] > thr) DataAC[3] += 0x0002; //See note above
149			if (fSumEcat[2] > thr) DataAC[3] += 0x0200; //
150			if (fSumEcat[3] > thr) DataAC[3] += 0x0020; //
151
152			if (fSumEcard[0] > thr) DataAC[3] += 0x1000; //C1D1(1)
153			if (fSumEcard[1] > thr) DataAC[3] += 0x0010; //C1D1(2)
154			if (fSumEcard[2] > thr) DataAC[3] += 0x0100; //C2D1(1)
155			if (fSumEcard[3] > thr) DataAC[3] += 0x0001; //C2D1(2)
156
157			DataAC[67] = DataAC[3];
158
159			// shift registers
160			// the central bin is equal to the hitmap, all
161			// other bins in the shift register are 0
162			for (UInt_t i=0; i<=15; i++){
163			DataAC[i+11] = 0x0000;
164			DataAC[i+75] = 0x0000;
165			}
166			DataAC[18] = DataAC[3];
167			DataAC[82] = DataAC[3];
168
169			//++++++ WRITE EVERYTHING TO DIGITIZER'S BUFFER +++//
170			fData.clear(); //clearing data buffer
171
172			for(Int_t i= 0; i<128; i++) fData.push_back(DataAC[i]);
173
174
175			if (fDEBUG){
176			cout<<"====== AC hit info ======"<<endl;
177			cout<<"CAT: Erel="<<fSumEcat[0]*1.e6<<" (keV)"<<endl
178			<<"CARD1: Erel="<<fSumEcard[0]*1.e6<<" (keV)"<<endl
179			<<"CARD2: Erel="<<fSumEcard[1]*1.e6<<" (keV)"<<endl
180			<<"CARD3: Erel="<<fSumEcard[2]*1.e6<<" (keV)"<<endl
181			<<"CARD4: Erel="<<fSumEcard[3]*1.e6<<" (keV)"<<endl
182			<<"CAS1: Erel="<<fSumEcas[0]*1.e6<<" (keV)"<<endl
183			<<"CAS2: Erel="<<fSumEcas[1]*1.e6<<" (keV)"<<endl
184			<<"CAS3: Erel="<<fSumEcas[2]*1.e6<<" (keV)"<<endl
185			<<"CAS4: Erel="<<fSumEcas[3]*1.e6<<" (keV)"<<endl;
186			cout<<"========================="<<endl;
187
188			}
189			}
190
191			void PamVMCAcDig::FillErel(){
192			/*
193			Detector KEY Number of paddles
194
195			CAT TOP1 1 (solid volume with id=1)
196			CAS SID1 4 (1-4)
197			CARD C1D1 & C2D1 2 & 2 (1-2)
198			*/
199			TClonesArray* hc=0;
200			const char* keyac[4] = {"TOP1","SID1","C1D1","C2D1"};
201			for(Int_t i=0; i<4; i++){
202			hc = (TClonesArray*)fhitscolmap.GetValue(keyac[i]);
203			if(hc) AddErel(hc, i);
204			hc = 0;
205			}
206			}
207
208			void PamVMCAcDig::AddErel(TClonesArray* HitColl, Int_t detcode){
209
210
211			// energy dependence on position (see file AcFitOutputDistancePmt.C by S.Orsi)
212			// based on J.Lundquist's calculations (PhD thesis, page 94)
213			// function: [0]+[1]*atan([2]/(x+1)), where the 3 parameters are:
214			// 8.25470e-01 +- 1.79489e-02
215			// 6.41609e-01 +- 2.65846e-02
216			// 9.81177e+00 +- 1.21284e+00
217			// hp: 1 minimum ionising particle at 35cm from the PMT releases 1mip
218
219			// TF1 attenAC = new TF1("fAttAC",".825+.64atan(9.8/x)",0.,45.);
220
221			const Float_t posCasPmt[4][3]={{28.308, -17.168, 63.644},// 1 - CAS CPU: x,y,z
222			{18.893, 24.913, 63.644},// 2 - CAS DCDC
223			{-24.307, 17.162, 63.644},// 3 - CAS VME
224			{-17.765, -28.300, 63.644}};// 4 - CAS IPM
225
226			Float_t dx,dy,dz,dAC; // distance from PMT
227			Int_t hitsNo=HitColl->GetEntriesFast();
228			PamVMCDetectorHit * hit = 0;
229
230			if (hitsNo>50) cout<<"WARNING: PamVMCAcDig::AddErel... "
231			<<"Number of hits in AC is out of range!!!"<<endl;
232
233			switch(detcode){
234			case 0: //CAT
235			for(Int_t k=0; k<hitsNo; k++){
236			hit = (PamVMCDetectorHit*)HitColl->At(k);
237			fSumEcat[(hit->GetPOS())-1] += hit->GetEREL();
238			}
239			break;
240
241			case 1: //CAS
242			for(Int_t k=0; k<hitsNo; k++){
243			hit = (PamVMCDetectorHit*)HitColl->At(k);
244			dx = (hit->GetXIN()+hit->GetXOUT())/2.-posCasPmt[hit->GetPOS()-1][0];
245			dy = (hit->GetYIN()+hit->GetYOUT())/2.-posCasPmt[hit->GetPOS()-1][1];
246			dz = (hit->GetZIN()+hit->GetZOUT())/2.-posCasPmt[hit->GetPOS()-1][2];
247			dAC = Sqrt(dxdx + dydy + dz*dz);
248			//fSumEcas[(hit->GetPOS())-1] += (hit->GetEREL())*attenAC->Eval(dAC);
249			if ((hit->GetPOS()==3)) fSumEcas[3] += (hit->GetEREL())*fattenAC->Eval(dAC);
250			if ((hit->GetPOS()==4)) fSumEcas[1] += (hit->GetEREL())*fattenAC->Eval(dAC);
251			if ((hit->GetPOS()==2)) fSumEcas[0] += (hit->GetEREL())*fattenAC->Eval(dAC);
252			if ((hit->GetPOS()==1)) fSumEcas[2] += (hit->GetEREL())*fattenAC->Eval(dAC);
253			}
254			break;
255
256			case 2: //C1D1
257			for(Int_t k=0; k<hitsNo; k++){
258			hit = (PamVMCDetectorHit*)HitColl->At(k);
259			if ((hit->GetPOS()==1)) fSumEcard[3] += hit->GetEREL();
260			if ((hit->GetPOS()==2)) fSumEcard[0] += hit->GetEREL();
261			//fSumEcard[(hit->GetPOS())-1] += hit->GetEREL();
262			}
263			break;
264			case 3: //C2D1
265			for(Int_t k=0; k<hitsNo; k++){
266			hit = (PamVMCDetectorHit*)HitColl->At(k);
267			if ((hit->GetPOS()==1)) fSumEcard[2] += hit->GetEREL();
268			if ((hit->GetPOS()==2)) fSumEcard[1] += hit->GetEREL();
269			//fSumEcard[(hit->GetPOS())+1] += hit->GetEREL();
270			}
271			break;
272
273			default:
274			break;
275			}
276
277
278
279
280			}