nd/src/PamVMCNDDig.cxx

#include "PamVMCNDDig.h"
#include <TRandom.h>
#include <TMath.h>


#define DEBUG 0

using TMath::Log;
using TMath::Exp;
using TMath::Abs;

ClassImp(PamVMCNDDig)

void PamVMCNDDig::Digitize(){

  cout<<"Starting ND Digitization...."<<endl;

  Int_t NdN = 0;

  const Double_t Showercut = 10.e-6; // 10 us
  const Double_t Window = 200.e-6; //200 us
  const Double_t Tau = 17.e-6; // 17 us. Time when amplitude of signal reduce
                               // in e times
  const Double_t SigmaTau = 4.e-6; //Sigma of fall time (obtained from 
                                   //Lebedev's oscillograms)
  const Double_t Tresh = 130; //According Lebedev th = 130 mV
  const Double_t keVtomV = Tresh/200.; //We suppose that 130 mV = 200 keV... everything linear
  

  Double_t deadtime[36];
  Double_t Time,U, Uback;
  Int_t pos;

  for (Int_t i = 0; i<36; i++) deadtime[i] = Showercut;

  TClonesArray * hc = (TClonesArray*)fhitscolmap.GetValue("NDTI");

  
  PamVMCDetectorHit * hit = 0;
  PamVMCDetectorHit * hit1 = 0;
  PamVMCDetectorHit *hit0 =0;
   

  if(hc){
    
    Int_t HitsNum = hc->GetEntriesFast();

    //sorting by time and put results into vector fpHits
    if (HitsNum > 1){
      for(Int_t i=0; i<hc->GetEntriesFast(); i++){
        Time = 0.;
        for(Int_t j=0; j<hc->GetEntriesFast(); j++){
          hit = (PamVMCDetectorHit*)hc->At(j);
          if((hit->GetTOF()>Time) && (!IsInsideVector(hit))){
            hit1 = hit;
            Time = hit->GetTOF();
          }
        }
        
        fpHits.push_back(hit1);
      }
    }
    
    //now we separate hits by tubes

    for(Int_t i = fpHits.size()-1; i>=0; i--){
      hit = (PamVMCDetectorHit*)fpHits.at(i);
      pos = hit->GetPOS()-1;
      tubeHits[pos].push_back(hit);
    }
    fpHits.clear();// we don't need it

    //From this moment all hits sorted by time and present in tubes

    for(Int_t i = 0; i<NTUBES; i++){
      if( tubeHits[i].size() > 0 ){
        for(UInt_t j = 0; j<tubeHits[i].size(); j++){
            Uback = 0.;
            hit =(PamVMCDetectorHit*)tubeHits[i].at(j);
            U = keVtomV*(hit->GetEREL())*1.e6; //initial amplitude of this signal;
            Time = hit->GetTOF();
            for(UInt_t k = 1; k<=j; k++){
              hit0 = (PamVMCDetectorHit*)tubeHits[i].at(k-1);
              //hit1 = (PamVMCDetectorHit*)tubeHits[i].at(j);
              if (tauHits[i].at(k-1)> 0)
                Uback += keVtomV*(hit0->GetEREL())*1.e6*Exp((hit0->GetTOF()-/*hit1->GetTOF()*/hit->GetTOF())/Abs(tauHits[i].at(k-1)));    
            }
              
            tauHits[i].push_back(frandom->Gaus(Tau,SigmaTau));
            if (DEBUG){
              cout<<"++++++HIT+++++++"<<endl;
              cout<<"TubeNo="<<hit->GetPOS()<<endl
                  <<"Particle="<<((TParticlePDG*)(TDatabasePDG::Instance()->GetParticle(hit->GetPDG())))->GetName()<<endl
                  <<"    Erel="<<hit->GetEREL()*1.0e6<<" (keV)"<<endl
                  <<"    Time="<<hit->GetTOF()*1.e6<<" (us)"<<endl
                  <<"    Patn="<<hit->GetPATH()*10.<<" (mm)"<<endl
                  <<"     Tau="<<tauHits[i].at(j)*1.e6<<" (us)"<<endl
                  <<"       U="<<U<<" (mV)"<<endl
                  <<"    Uback="<<Uback<<" (mV)"<<endl;
              if ((Uback < Tresh) && (U+Uback > Tresh) && (Time<Showercut+Window) && (Time>Showercut))
                cout<<" Particle detected as NEUTRON"<<endl;
            }
            if ((Uback < Tresh) && (U+Uback > Tresh) && (Time<Showercut+Window) && (Time>Showercut)) NdN++;

          } 
      }
    }

   if(DEBUG)
     cout<<"==== Total nuber of detected neutrons for for event ==="<<NdN<<" ==="<<endl;

      }  
    
          
  //Cleaning of everything...
  for (Int_t i = 0; i<NTUBES; i++){
    tauHits[i].clear();
    tubeHits[i].clear();
  }
  
  fData.clear();
  
  if (DEBUG) cout<<"REGISTERED: "<< NdN <<" neutrons"<<endl;
  for(Int_t i=0;i<3;i++){
    fData.push_back(0x0000);
    fData.push_back(0x010F);
  }
  fData.at(0) = ((0xFF00 & (256*NdN)));
}
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.e6Exp((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	}