| 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 |
|
| 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(gRandom->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 |
} |
Parent Directory
| 
Revision Log