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	#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	}