Commit 6833af68 authored by Vratislav Chudoba's avatar Vratislav Chudoba

Class AEvent repaired.

parent 4a13def9
void read_root()
{
TFile *f = new TFile("/home/dariak/NeuRad_tests/data/rawDataDSR4/NeuRad_test_07_1.root");
TTree *tr = (TTree*)f->Get("rtree");
tr->SetMakeClass(1);
// TArrayD ampl1;
// TArrayD time1;
// ampl1.Set(ncellMax);
// time1.Set(ncellMax);
const int ncellMax = 1030;
// int ncells;
Double_t ampl1[1023];
Double_t time1[1023];
Double_t ampl1_pos[1023]; //branch for positive signals
Double_t maxAmplitude1;
Double_t timemaxAmplitude1;
Double_t time1_pos[1023];
TFile *fw = new TFile("/home/dariak/NeuRad_tests/data/dataDSR4/analysis_07_1.root", "RECREATE");
TTree *tw = new TTree("drs4analysis", "title of drs4 analysis tree");
tr->SetBranchAddress("amp_ch1", ampl1);
tr->SetBranchAddress("time_ch1", time1);
// tr->SetBranchAddress("ncell", ncells);
tw->Branch("maxAmplitude1", &maxAmplitude1, "maxAmplitude1/D");
tw->Branch("timemaxAmplitude1", &timemaxAmplitude1, "timemaxAmplitude1/D");
tw->Branch("ampl1_pos", ampl1_pos, "ampl1_pos[1023]/D"); //branch for positive signals
tw->Branch("time1_pos", time1_pos, "time1_pos[1023]/D");
const Long64_t nEntries = tr->GetEntries();
std::cout <<"Number of entries: "<<nEntries<<std::endl;
//----for pictures
/* TCanvas *c1 = new TCanvas();
c1->Divide(3,4);
for(Int_t i=12; i<24; i++) {
c1->cd(i-11);
tr->GetEntry(i);
TGraph *gr1 = new TGraph(1023, time1, ampl1);
gr1->SetTitle("Signal shape for one event");
gr1->GetXaxis();
gr1->Draw("Al*");
}
tr->GetEntry(499);
TGraph *gr1 = new TGraph(1023, time1, ampl1);
gr1->SetTitle("Signal shape for one event");
gr1->Draw("Al*");*/
//-----
Double_t maxAmpl1 = 0.;
Double_t timemaxAmpl1 = 0.;
for(Int_t i=0; i<tr->GetEntries(); i++) {
tr->GetEntry(i);
//changing polarity of the signal
for(Int_t j=0; j<1023; j++) {
ampl1_pos[j] = ampl1[j]*(-1.);
time1_pos[j] = time1[j];
}
//find maximum by hand
maxAmpl1 = ampl1_pos[0];
for(Int_t j=0; j<1023; j++) {
if(ampl1_pos[j] > maxAmpl1) {
maxAmpl1 = ampl1_pos[j];
timemaxAmpl1 = time1_pos[j];
}
//cout<<"Time "<<i<<" "<<j<<" "<<time1[j]<<endl;
//cout<<"Amplitude "<<i<<" "<<j<<" "<<ampl1[j]<<endl;
//cout<<endl;
}
maxAmplitude1 = maxAmpl1;
timemaxAmplitude1 = timemaxAmpl1;
std::cout<<"Max amplitude "<<maxAmpl1<<" for entry "<<i<< std::endl;
std::cout<<"Time for max amplitude "<<timemaxAmpl1<<" for entry "<<i<< std::endl;
//fitting 90 percent of rising edge
for(Int_t k=0; k<1023; k++) {
if(ampl1[k] > 0.1*maxAmpl1 && ampl1[k] < 0.9*maxAmpl1) { //we have negative signals, amplitude between 10 and 90 proc from maximum
/*cout<<"rjgnfre"<<endl;
TGraph *gr2 = new TGraph(1023, time1, ampl1);
gr2->SetTitle("stupido");
gr2->Fit("pol1");
gr2->Draw("Al*");*/
}
}
//getting integral
TGraph *gr3 = new TGraph(1023, time1, ampl1_pos);
gr3->Integral(120,180);
std::cout << "INtegral "<<gr3->Integral(120,180)<< std::endl;
tw->Fill();
}
tw->GetEntry(498);
TGraph *gr2 = new TGraph(1023, time1_pos, ampl1_pos);
gr2->SetTitle("Signal shape for one gsgsg event");
gr2->Draw("Al*");
fw->cd();
tw->Write();
fw->Close();
}
...@@ -94,7 +94,7 @@ void AEvent::Reset() { ...@@ -94,7 +94,7 @@ void AEvent::Reset() {
fTimeAmpMax = 0.; fTimeAmpMax = 0.;
fTimeCFD = 0.; fTimeCFD = 0.;
fZeroLevel = 0.; fZeroLevel = 0.;
fCharge = 0.; fChargeCFD = 0.;
} }
void AEvent::SetInputEvent(RawEvent** event) { void AEvent::SetInputEvent(RawEvent** event) {
...@@ -195,17 +195,17 @@ Double_t AEvent::FindZeroLevel(Int_t pmin, Int_t pmax) { ...@@ -195,17 +195,17 @@ Double_t AEvent::FindZeroLevel(Int_t pmin, Int_t pmax) {
void AEvent::SetChargeCFD(Int_t tmin, Int_t tmax) { void AEvent::SetChargeCFD(Int_t tmin, Int_t tmax) {
Double_t integral = 0.; //voltage Double_t integral = 0.; //voltage
Double_t time_sig = 0; //approximate signal duration in seconds Double_t time_sig = 0; //approximate signal duration in seconds
Double_t res = 50.; //resistance 50 Om const Double_t res = 50.; //resistance 50 Om
time_sig = (double)(tmin + tmax)*(1e-9); time_sig = (double)(-tmin + tmax)*(1e-9);
for(Int_t i = 0; i < fNPoints; i++) { for(Int_t i = 0; i < fNPoints; i++) {
if( fTime[i] > (fTimeCFD - tmin) && fTime[i] < (fTimeCFD + tmax) ) { if( fTime[i] > (fTimeCFD + tmin) && fTime[i] < (fTimeCFD + tmax) ) {
integral = integral + fAmpPos[i]; integral = integral + fAmpPos[i];
} }
} }
fCharge = integral*time_sig/res; fChargeCFD = integral*time_sig/res;
// cout<<fCharge<<endl; // cout<<fCharge<<endl;
// printf("\nIntegral is %f , charge is %lf \n", integral, fCharge); // printf("\nIntegral is %f , charge is %lf \n", integral, fCharge);
......
...@@ -42,11 +42,10 @@ private: ...@@ -42,11 +42,10 @@ private:
TArrayD fAmpCFD; //array for CFD amplitudes (attenuated, inversed and delayed) TArrayD fAmpCFD; //array for CFD amplitudes (attenuated, inversed and delayed)
Double_t fTimeCFD; //zero-crossing time Double_t fTimeCFD; //zero-crossing time
Double_t fCharge; //charge of the signal in Coulomb Double_t fChargeCFD; //charge of the signal in Coulomb
TGraph *fGraphSignal; TGraph *fGraphSignal;
TGraph *fGraphCFD; TGraph *fGraphCFD;
TGraph *fGraphZero;
RawEvent *fInputEvent; //! RawEvent *fInputEvent; //!
...@@ -90,7 +89,7 @@ public: ...@@ -90,7 +89,7 @@ public:
//for zero level correction. one parameter fit between pmin and pmax //for zero level correction. one parameter fit between pmin and pmax
//returns fit parameter i.e. number on which amplitude should be corrected //returns fit parameter i.e. number on which amplitude should be corrected
void SetChargeCFD(Int_t tmin = 3, Int_t tmax = 17); void SetChargeCFD(Int_t tmin = -3, Int_t tmax = 17);
//calculates charge of the signal (i.e. its integral //calculates charge of the signal (i.e. its integral
//in range of (tmin,tmax) in ns) //in range of (tmin,tmax) in ns)
//CFD time is taken as a start point //CFD time is taken as a start point
......
...@@ -5,7 +5,7 @@ void testShowCFD() ...@@ -5,7 +5,7 @@ void testShowCFD()
const Long64_t kFirstEvent = 128; const Long64_t kFirstEvent = 128;
TFile fr("../data/dataDSR4/analysis_07_1.root"); TFile fr("../data/dataDSR4/analysis_07_8.root");
TTree *tr = (TTree*)fr.Get("atree"); TTree *tr = (TTree*)fr.Get("atree");
AEvent *revent = new AEvent(); AEvent *revent = new AEvent();
......
...@@ -5,7 +5,8 @@ void testShowGraphs() ...@@ -5,7 +5,8 @@ void testShowGraphs()
const Long64_t kFirstEvent = 128; const Long64_t kFirstEvent = 128;
TFile fr("../data/dataDSR4/analysis_07_1.root"); // TFile fr("../data/dataDSR4/analysis_07_1.root");
TFile fr("../data/dataDSR4/analysis_07_8.root");
TTree *tr = (TTree*)fr.Get("atree"); TTree *tr = (TTree*)fr.Get("atree");
AEvent *revent = new AEvent(); AEvent *revent = new AEvent();
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment