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

Annotation of /trieste/pamVMC/nd/src/PamVMCNDDig.cxx

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.2 - (hide annotations) (download)
Tue Mar 24 14:04:03 2009 UTC (15 years, 11 months ago) by pizzolot
Branch: MAIN
CVS Tags: HEAD
Changes since 1.1: +3 -2 lines
setting of random seed; new distribution for primary generation

1 pamelats 1.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     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 pizzolot 1.2
75 pamelats 1.1 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 pizzolot 1.2
90     tauHits[i].push_back(rnd->Gaus(Tau,SigmaTau));
91 pamelats 1.1 if (DEBUG){
92     cout<<"++++++HIT+++++++"<<endl;
93     cout<<"TubeNo="<<hit->GetPOS()<<endl
94     <<"Particle="<<((TParticlePDG*)(TDatabasePDG::Instance()->GetParticle(hit->GetPDG())))->GetName()<<endl
95     <<" Erel="<<hit->GetEREL()*1.0e6<<" (keV)"<<endl
96     <<" Time="<<hit->GetTOF()*1.e6<<" (us)"<<endl
97     <<" Patn="<<hit->GetPATH()*10.<<" (mm)"<<endl
98     <<" Tau="<<tauHits[i].at(j)*1.e6<<" (us)"<<endl
99     <<" U="<<U<<" (mV)"<<endl
100     <<" Uback="<<Uback<<" (mV)"<<endl;
101     if ((Uback < Tresh) && (U+Uback > Tresh) && (Time<Showercut+Window) && (Time>Showercut))
102     cout<<" Particle detected as NEUTRON"<<endl;
103     }
104     if ((Uback < Tresh) && (U+Uback > Tresh) && (Time<Showercut+Window) && (Time>Showercut)) NdN++;
105    
106     }
107     }
108     }
109    
110     if(DEBUG)
111     cout<<"==== Total nuber of detected neutrons for for event ==="<<NdN<<" ==="<<endl;
112    
113     }
114    
115    
116    
117     //Cleaning of everything...
118     for (Int_t i = 0; i<NTUBES; i++){
119     tauHits[i].clear();
120     tubeHits[i].clear();
121     }
122    
123     fData.clear();
124    
125     if (DEBUG) cout<<"REGISTERED: "<< NdN <<" neutrons"<<endl;
126     for(Int_t i=0;i<3;i++){
127     fData.push_back(0x0000);
128     fData.push_back(0x010F);
129     }
130     fData.at(0) = ((0xFF00 & (256*NdN)));
131     }

  ViewVC Help
Powered by ViewVC 1.1.23