integralFormSignal.cpp

{
	gSystem->Load("../libData.so");

	TFile *f = new TFile("../data/dataDSR4/Neurad_081216_2NEW.root");
	TTree *tr = (TTree*)f->Get("atree");
	
	const Int_t nP = 1024;

	const Int_t noBranches = 4;
	Double_t ZeroTime[4],mh[4],RisTime[4],decayT[4],tmax[4],timestep,time;
	Int_t iZero[4],PosZero[4],deltaT[4],Eh[4],Bh[4]; 
	TString bName;
	AEvent *aevent[noBranches];	// pointer to the array (of RawEvent class) in which raw data for each channel will be put
	for (Int_t j = 0; j<noBranches; j++) {
		aevent[j] = new AEvent();	//each raw event element is of class RawEvent()
		bName.Form("Ach%d.", j);
		tr->SetBranchAddress(bName.Data(), &aevent[j]);	//read the tree tr with raw data and fill array revent with raw data
	}

	TString hname,Chname;
	TH1F *hist[4];
	TGraph *gr[4];
	for (Int_t i = 0; i < 4; i++) {
		hname.Form("hist%d",i);
		Chname.Form("CHANNEL%d",i+1);
		hist[i] = new TH1F(hname.Data(), Chname.Data(), 2000, -100, 100);
	}

	TH1F *hist01 = new TH1F("hist11", "h1 raw", 2000, -100, 100);

	//TH1(const char* name, const char* title, Int_t nbinsx, Double_t xlow, Double_t xup)

	Long64_t nEntries = tr->GetEntries();	
	TH1F *hSumm = new TH1F("hSumm", "integral signal", nP-1, 0, 199);
	for(Int_t j=0;j<(nP-1);j++) {hSumm->SetBitContent(j,0);}
	timestep = 199/1022; 	//!!!!!!!!!!! не 1023 а 1022. промежутков меньше на 1 чем точек
	//loop over events
	for (Int_t j = 0; j < nEntries; j++) {
		tr->GetEntry(j);
		for(Int_t k=0;k<4;k++) { // loop for channels
			ZeroTime[k] = aevent[k]->GetfCFD(); // get fCFD for channel №k
			iZero[k] = ZeroTime[k]*10; 
			//if(j==0) { PosZero[k] = iZero[k]; }
			//deltaT[k] = iZero[k] - PosZero[k]; // количество точек !!! 
			if(j==0) { PosZero[k] = ZeroTime[k]; }
			deltaT[k] = ZeroTime[k] - PosZero[k]; // время в нс
		}

		gr[0] = new TGraph(*aevent[0]->GetGraphSignal());
		for(Int_t i = 0; i<(nP-1); i++){	
			//if( ((i+deltaT[0])>-1) && ((i+deltaT[0])<nP))	hist[0]->AddBinContent(i+500,aevent[0]->GetOnefAmpPos(i+deltaT[0]));	
			//hist01->AddBinContent(i+500,aevent[0]->GetOnefAmpPos(i));
			if( ((timestep*i + deltaT[0])>=0) && ((timestep*i + deltaT[0])<tmax[0])) { 
				hSumm->AddBinContent( i, fGraphCFD->Eval(timestep*i + deltaT[0]) ); 
			}
		}
	}
		
	//t->Draw("Ach2.fAmpMax: ( Ach0.fTimeCFD - Ach2.fTimeCFD ) >> h4","","col");
	/*for(j=0;j<4;j++) {

		mh[j] = hist[j]->GetBinContent(hist[j]->GetMaximumBin());

		for(i = hist[j].GetMaximumBin(); i>0;i--) {	// finding the 
			if( hist[j]->GetBinContent(i)<0.9*mh[j] ) {
				Eh[j] = i+1;
				break;
			}
		}

		for(i = hist[j].GetMaximumBin(); i>0;i--) {	// finding the 
			if( hist[j]->GetBinContent(i)<0.1*mh[j] ) {
				Bh[j] = i+1;
				break;
			}
		}
		RisTime[j] = (Eh[j] - Bh[j])*0.1; // время нарастания 
	}
*/
	Double_t rMin[4];
	Double_t rMax[4];
	rMin[0]=22.5; rMax[0] = 27.3;
	rMin[1]=26.6; rMax[1] = 31.7;
	rMin[2]=19.3; rMax[2] = 24.;
	rMin[3]= 26.5; rMax[3] = 29.6;

	TF1 *fit1 = new TF1("fit1","[0]*exp(-x*[1])");
	fit1->SetParName(1,"tD");

	/*TCanvas *c1 = new TCanvas("c1","run 071216_8",10,10,1000,600);
	c1->Divide(2,2);
	for(i=0;i<4;i++) {
		c1->cd(i+1);
		hist[i]->GetXaxis()->SetRangeUser(10, 50);
		hist[i]->Draw();
		fit1->SetParLimits(1,0.1,1);
		fit1->SetRange(rMin[i],rMax[i]);
		//hist[i]->Fit("fit1","R");
		//decayT[i] = 1/(fit1->GetParameter(1));
	}
	
	for(i=0;i<4;i++) {
		cout<<decayT[i]<<" that's decay time  of channel number "<<i<<endl;
		cout<<RisTime[i]<< " that's a RiseTime for this channel "<<i<<endl;
	}*/
///////////////////////

	TCanvas *c2 = new TCanvas("c2","RAW",10,10,1000,600);
	//c2->Divide(2,2);
	//c2->cd(1);

	//hist01->GetXaxis()->SetRangeUser(10, 55);
	//hist01->Draw();
	
	//c2->cd(2);
	hSumm->Draw();
	/*
	c2->cd(2);
	hist11->GetXaxis()->SetRangeUser(10, 55);
	hist11->Draw();

	c2->cd(3);
	hist21->GetXaxis()->SetRangeUser(10, 55);
	hist21->Draw();

	c2->cd(4);
	hist31->GetXaxis()->SetRangeUser(10, 55);
	hist31->Draw();*/
}