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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

1 nikolas 1.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