/[PAMELA software]/PamVMC_update/nd/src/PamVMCNDDig.cxx
ViewVC logotype

Contents of /PamVMC_update/nd/src/PamVMCNDDig.cxx

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1.1.1 - (show annotations) (download) (vendor branch)
Tue Oct 15 15:51:45 2013 UTC (11 years, 1 month ago) by formato
Branch: MAIN, rel
CVS Tags: reltag, HEAD
Changes since 1.1: +0 -0 lines
PamVMC update

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 }

  ViewVC Help
Powered by ViewVC 1.1.23