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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1.1.1 - (hide 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 formato 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     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