read_root.C

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();
}