1 |
#include "PamVMCNDDig.h" |
2 |
#include <TRandom.h> |
3 |
#include <TMath.h> |
4 |
|
5 |
|
6 |
#define DEBUG 0 |
7 |
|
8 |
using TMath::Log; |
9 |
using TMath::Exp; |
10 |
using TMath::Abs; |
11 |
|
12 |
ClassImp(PamVMCNDDig) |
13 |
|
14 |
void PamVMCNDDig::Digitize(){ |
15 |
|
16 |
if(DEBUG) cout<<"Starting ND Digitization...."<<endl; |
17 |
|
18 |
Int_t NdN = 0; |
19 |
|
20 |
const Double_t Showercut = 10.e-6; // 10 us |
21 |
const Double_t Window = 200.e-6; //200 us |
22 |
const Double_t Tau = 17.e-6; // 17 us. Time when amplitude of signal reduce |
23 |
// in e times |
24 |
const Double_t SigmaTau = 4.e-6; //Sigma of fall time (obtained from |
25 |
//Lebedev's oscillograms) |
26 |
const Double_t Tresh = 130; //According Lebedev th = 130 mV |
27 |
const Double_t keVtomV = Tresh/200.; //We suppose that 130 mV = 200 keV... everything linear |
28 |
|
29 |
|
30 |
Double_t deadtime[36]; |
31 |
Double_t Time,U, Uback; |
32 |
Int_t pos; |
33 |
|
34 |
for (Int_t i = 0; i<36; i++) deadtime[i] = Showercut; |
35 |
|
36 |
TClonesArray * hc = (TClonesArray*)fhitscolmap.GetValue("NDTI"); |
37 |
|
38 |
|
39 |
PamVMCDetectorHit * hit = 0; |
40 |
PamVMCDetectorHit * hit1 = 0; |
41 |
PamVMCDetectorHit *hit0 =0; |
42 |
|
43 |
|
44 |
if(hc){ |
45 |
|
46 |
Int_t HitsNum = hc->GetEntriesFast(); |
47 |
|
48 |
//sorting by time and put results into vector fpHits |
49 |
if (HitsNum > 1){ |
50 |
for(Int_t i=0; i<hc->GetEntriesFast(); i++){ |
51 |
Time = 0.; |
52 |
for(Int_t j=0; j<hc->GetEntriesFast(); j++){ |
53 |
hit = (PamVMCDetectorHit*)hc->At(j); |
54 |
if((hit->GetTOF()>Time) && (!IsInsideVector(hit))){ |
55 |
hit1 = hit; |
56 |
Time = hit->GetTOF(); |
57 |
} |
58 |
} |
59 |
|
60 |
fpHits.push_back(hit1); |
61 |
} |
62 |
} |
63 |
|
64 |
//now we separate hits by tubes |
65 |
|
66 |
for(Int_t i = fpHits.size()-1; i>=0; i--){ |
67 |
hit = (PamVMCDetectorHit*)fpHits.at(i); |
68 |
pos = hit->GetPOS()-1; |
69 |
tubeHits[pos].push_back(hit); |
70 |
} |
71 |
fpHits.clear();// we don't need it |
72 |
|
73 |
//From this moment all hits sorted by time and present in tubes |
74 |
|
75 |
for(Int_t i = 0; i<NTUBES; i++){ |
76 |
if( tubeHits[i].size() > 0 ){ |
77 |
for(UInt_t j = 0; j<tubeHits[i].size(); j++){ |
78 |
Uback = 0.; |
79 |
hit =(PamVMCDetectorHit*)tubeHits[i].at(j); |
80 |
U = keVtomV*(hit->GetEREL())*1.e6; //initial amplitude of this signal; |
81 |
Time = hit->GetTOF(); |
82 |
for(UInt_t k = 1; k<=j; k++){ |
83 |
hit0 = (PamVMCDetectorHit*)tubeHits[i].at(k-1); |
84 |
//hit1 = (PamVMCDetectorHit*)tubeHits[i].at(j); |
85 |
if (tauHits[i].at(k-1)> 0) |
86 |
Uback += keVtomV*(hit0->GetEREL())*1.e6*Exp((hit0->GetTOF()-/*hit1->GetTOF()*/hit->GetTOF())/Abs(tauHits[i].at(k-1))); |
87 |
} |
88 |
|
89 |
tauHits[i].push_back(frandom->Gaus(Tau,SigmaTau)); |
90 |
if (DEBUG){ |
91 |
cout<<"++++++HIT+++++++"<<endl; |
92 |
cout<<"TubeNo="<<hit->GetPOS()<<endl |
93 |
<<"Particle="<<((TParticlePDG*)(TDatabasePDG::Instance()->GetParticle(hit->GetPDG())))->GetName()<<endl |
94 |
<<" Erel="<<hit->GetEREL()*1.0e6<<" (keV)"<<endl |
95 |
<<" Time="<<hit->GetTOF()*1.e6<<" (us)"<<endl |
96 |
<<" Patn="<<hit->GetPATH()*10.<<" (mm)"<<endl |
97 |
<<" Tau="<<tauHits[i].at(j)*1.e6<<" (us)"<<endl |
98 |
<<" U="<<U<<" (mV)"<<endl |
99 |
<<" Uback="<<Uback<<" (mV)"<<endl; |
100 |
if ((Uback < Tresh) && (U+Uback > Tresh) && (Time<Showercut+Window) && (Time>Showercut)) |
101 |
cout<<" Particle detected as NEUTRON"<<endl; |
102 |
} |
103 |
if ((Uback < Tresh) && (U+Uback > Tresh) && (Time<Showercut+Window) && (Time>Showercut)) NdN++; |
104 |
|
105 |
} |
106 |
} |
107 |
} |
108 |
|
109 |
if(DEBUG) |
110 |
cout<<"==== Total nuber of detected neutrons for for event ==="<<NdN<<" ==="<<endl; |
111 |
|
112 |
} |
113 |
|
114 |
|
115 |
|
116 |
//Cleaning of everything... |
117 |
for (Int_t i = 0; i<NTUBES; i++){ |
118 |
tauHits[i].clear(); |
119 |
tubeHits[i].clear(); |
120 |
} |
121 |
|
122 |
fData.clear(); |
123 |
|
124 |
if (DEBUG) cout<<"REGISTERED: "<< NdN <<" neutrons"<<endl; |
125 |
for(Int_t i=0;i<3;i++){ |
126 |
fData.push_back(0x0000); |
127 |
fData.push_back(0x010F); |
128 |
} |
129 |
fData.at(0) = ((0xFF00 & (256*NdN))); |
130 |
} |