/[PAMELA software]/PamVMC/s4/src/PamVMCS4Dig.cxx
ViewVC logotype

Contents of /PamVMC/s4/src/PamVMCS4Dig.cxx

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (show annotations) (download)
Thu Feb 19 16:00:16 2009 UTC (15 years, 9 months ago) by nikolas
Branch: MAIN
Cleaning files before release

1 #include "PamVMCS4Dig.h"
2 #include <TRandom.h>
3
4 ClassImp(PamVMCS4Dig)
5
6 void PamVMCS4Dig::Digitize(){
7
8 cout<<"Starting S4 Digitization...."<<endl;
9
10 //This is simple approximation, reflection + atennuation..
11 //I hope in future we will use for generation and propagation of
12 //optical photons abilitiees of montecarlo...
13
14
15
16
17 Int_t DEBUG=1;
18 // creato: S. Borisov, INFN Roma2 e MEPHI, Sett 2007
19 Int_t i,j,t,NdF,pmt,S4;
20 Float_t E0,X,Y,Z,x,y,z,V[3],q,w,l;
21
22 const Float_t Xs[2] = {-24.1, 24.1};
23 const Float_t Ys[2] = {-24.1, 24.1};
24 const Float_t Zs[2] = {-0.5, 0.5};
25 const Float_t Yp[6] = {-20., -17., -1., 2., 17., 20.}; //PMT positions
26
27 const Float_t E1 = 1e-6; //Energy which we spend to produce 1 photon
28 const Float_t p = 0.1;
29 const Int_t l0 = 500;
30
31 const Int_t S4p = 32; //pedestal
32
33 UShort_t S4v = 0; // ADC count
34
35 Int_t NdFT = 0;
36
37 Float_t Ert = 0.;
38
39 TClonesArray * hc = (TClonesArray*)fhitscolmap.GetValue("S4");
40 PamVMCDetectorHit * hit = 0;
41
42 if(hc){
43 for(i=0;i<hc->GetEntriesFast();i++){
44
45 hit = (PamVMCDetectorHit*)hc->At(i);
46
47 Ert += hit->GetEREL();
48 NdF=Int_t(hit->GetEREL()/E1);
49 NdFT=0;
50 X=hit->GetXIN();
51 Y=hit->GetYIN();
52 Z=(Float_t)(random())/(Float_t)(0x7fffffff)-0.5;
53 for(j=0;j<NdF;j++){
54 q=(Float_t)random()/(Float_t)0x7fffffff;
55 w=(Float_t)random()/(Float_t)0x7fffffff;
56 V[0]=p*cos(6.28318*q);
57 V[1]=p*sin(6.28318*q);
58 V[2]=p*(2.*w-1.);
59 pmt=0;
60 x=X;
61 y=Y;
62 z=Z;
63 while(pmt==0 && (x>Xs[0] && x<Xs[1])&&(y>Ys[0] && y<Ys[1])&&(z>Zs[0] && z<Zs[1])){
64 l=0;
65 while(pmt==0 && (x>Xs[0] && x<Xs[1])&&(y>Ys[0] && y<Ys[1])&&(z>Zs[0] && z<Zs[1])){
66 x+=V[0];
67 y+=V[1];
68 z+=V[2];
69 l+=p;
70 }
71 if((x<Xs[0]+p || x>Xs[1]-p)&&(y>Ys[0]+p && y<Ys[1]-p)&&(z>Zs[0]+p && z<Zs[1]-p)){
72 for(t=0;t<3;t++){
73 if(y>=Yp[2*t] && y<Yp[2*t+1]){
74 if(pmt==0)NdFT++;
75 pmt=1;
76 break;
77 }
78 }
79 if(pmt==1)break;
80 V[0]=-V[0];
81 }
82 q=(Float_t)random()/(Float_t)0x7fffffff;
83 w=1-exp(-l/l0);
84 if(q<w)break;
85 q=(Float_t)random()/(Float_t)0x7fffffff;
86 w=0.5;
87 if(q<w)break;
88 if((x>Xs[0]+p && x<Xs[1]-p)&&(y<Ys[0]+p || y>Ys[1]-p)&&(z>Zs[0]+p && z<Zs[1]-p))V[1]=-V[1];
89 if((x>Xs[0]+p && x<Xs[1]-p)&&(y>Ys[0]+p && y<Ys[1]-p)&&(z<Zs[0]+p || z>Zs[1]-p))V[2]=-V[2];
90 x+=V[0];
91 y+=V[1];
92 z+=V[2];
93 l=0;
94 }
95 }
96 }
97 Ert=Ert/0.002;
98 q=(Float_t)(random())/(Float_t)0x7fffffff;
99 w=0.7;
100 E0=4064./7.;
101 if(Ert<1) S4=0;
102 else S4=(Int_t)(4064.*(1.-exp(-(Ert-1.)/E0)));
103 i=S4/4;
104 if(S4%4==0)
105 S4v=S4+S4p;
106 else if(S4%4==1){
107 if(q<w) S4v=S4-1+S4p;
108 else S4v=S4+1+S4p;
109 } else if(S4%4==2) S4v=S4+S4p;
110 else if(S4%4==3){
111 if(q<w) S4v=S4+1+S4p;
112 else S4v=S4-1+S4p;
113 }
114 if (DEBUG)
115 cout<<"Ert_S4 = " << Ert << " --- S4v = " << S4v << endl;
116
117 fData.clear();
118
119 fData.push_back(S4v);
120 fData.push_back(0xd800);
121 fData.push_back(0x0300);
122 }
123 }

  ViewVC Help
Powered by ViewVC 1.1.23