//////////////////////////////////// // KRATTA Array Event Analyser // for the Asy-Eos experiment // TKRATAnalyser IMPLEMENTATION // July 2012 // revison 7/2012 // E.d.F ver 1.0 // sebastian.kupny@uj.edu.pl //////////////////////////////////// #include "TKRATAnalyser.h" ClassImp ( TKRATAnalyser ); /// %%%%%%%%%%%%%%%%%%%%%%%%%%% TKRATAnalyser %%%%%%%%%%%%%%%%%%%%%%%%%% //______________________________________________________________________ TKRATAnalyser::TKRATAnalyser( void ) { // } //______________________________________________________________________ TKRATAnalyser::~TKRATAnalyser() { } //______________________________________________________________________ Bool_t TKRATAnalyser::isPileUp ( Int_t * A_array, Int_t A_arraySize ) { Bool_t result = kFALSE; TH1F * rHist = new TH1F("tempHist", "tempHist", A_arraySize, 0, A_arraySize ); for ( int bin = 0; bin < A_arraySize ; bin++) { rHist->SetBinContent( bin, A_array[bin] ); } result = isPileUp( rHist ); if (result ) rHist->SaveAs("Histogram_with_pileup.root"); delete rHist; return result; }; //______________________________________________________________________ Bool_t TKRATAnalyser::isPileUp ( TH1F * A_hist ) { Bool_t result = false; TF1 * f1; TH1F * hpxDistr = new TH1F("hpxDistr", "", A_hist->GetNbinsX (), 1, A_hist->GetNbinsX ()); Double_t aSlope = 0.; Int_t BinsToFit = 30; ///Iteracyjnie fitowane sa proste oraz wyznaczane wspolczynniki kierunkowe for ( int kk = BinsToFit; kk < A_hist->GetNbinsX () - BinsToFit ; kk++){ //< Na ktorym binie ma zakonczyc sie analiza. Bo jest 1024 aSlope = 0; A_hist->Fit("pol1", "Q", "", kk - BinsToFit/2, kk + BinsToFit/2); f1 = A_hist->GetFunction("pol1"); aSlope = f1 -> GetParameter(1); hpxDistr->SetBinContent((Double_t)kk, aSlope ); } //hpxDistr->SaveAs("temp_pilep_test.root"); Int_t regist[3] = {0}; ///Okreslenie poziomu od ktorego wartosc pochodnej sygnalu ///jest akceptowana. Badamy na wszystkich trzech poziomach: Double_t threshold[3] = {0.9, 1., 1.1}; ///Detekcja narastajacego zbocza dla trzech roznych poziomow for(int j = 0; j < 3; j++){ Bool_t beforeAbove = false; Bool_t nowAbove = false; for(int i = 0; i < A_hist->GetNbinsX (); i++){ if (hpxDistr -> GetBinContent(i) >= threshold[j]){ nowAbove = true; }else{ nowAbove = false; } if(nowAbove && !beforeAbove){ //jezeli teraz jest wieksze a wczesnie bylo mniejsze regist[j]++; } beforeAbove = nowAbove; } } if(regist[0] > 1 && regist[1] > 1 && regist[2] > 1){ result = true; } delete hpxDistr; return result; } //______________________________________________________________________ ///This function is an overlay for proper function Double_t TKRATAnalyser::getInterpolatedValue(Double_t x, Int_t x_first, Int_t * TabY, Int_t ArraySize) { if(ArraySize <= 0) { return 0.0; } Int_t TabSize = ArraySize; Int_t * TabX = new Int_t[TabSize]; Double_t * TabYd = new Double_t[TabSize]; for ( Int_t i = 0; i < TabSize; i++) { TabX[i]= x_first + i; TabYd[i] = (Double_t) TabY[i]; } Double_t result = getInterpolatedValue ( x, TabX, TabYd, TabSize); delete TabX; return result; } //______________________________________________________________________ ///This function is an overlay for proper function Double_t TKRATAnalyser::getInterpolatedValue(Double_t x, Int_t x_first, Int_t * TabY, Int_t ArraySize, Double_t &A_p0, Double_t &A_p1, Double_t &A_p2) { if(ArraySize <= 0) { return 0.0; } Int_t TabSize = ArraySize; Int_t * TabX = new Int_t[TabSize]; Double_t * TabYd = new Double_t[TabSize]; for ( Int_t i = 0; i < TabSize; i++) { TabX[i]= x_first + i; TabYd[i] = (Double_t) TabY[i]; } Double_t result = getInterpolatedValue ( x, TabX, TabYd, TabSize, A_p0, A_p1, A_p2); delete TabX; return result; } //______________________________________________________________________ ///This function is an overlay for proper function Double_t TKRATAnalyser::getInterpolatedValue(Double_t x, Int_t x_first, Float_t * TabY, Int_t ArraySize, Double_t &A_p0, Double_t &A_p1, Double_t &A_p2) { if(ArraySize <= 0) { return 0.0; } Int_t TabSize = ArraySize; Int_t * TabX = new Int_t[TabSize]; Double_t * TabYd = new Double_t[TabSize]; for ( Int_t i = 0; i < TabSize; i++) { TabX[i]= x_first + i; TabYd[i] = (Double_t) TabY[i]; } Double_t result = getInterpolatedValue ( x, TabX, TabYd, TabSize, A_p0, A_p1, A_p2); delete TabX; return result; } //______________________________________________________________________ ///This function is an overlay for proper function Double_t TKRATAnalyser::getInterpolatedValue(Double_t x, Int_t x1, Double_t y1, Int_t x2, Double_t y2, Int_t x3, Double_t y3) { Int_t TabSize = 3; Int_t * TabX = new Int_t[TabSize]; Double_t * TabY = new Double_t[TabSize]; TabX[0]= x1; TabY[0]= y1; TabX[1]= x2; TabY[1]= y2; TabX[2]= x3; TabY[2]= y3; Double_t result = getInterpolatedValue ( x, TabX, TabY, TabSize); delete TabX; delete TabY; return result; } //____________________________________________________________________________________ ///This function returns interpolation result for given arrays and polynomial 2-nd order interpolation Double_t TKRATAnalyser::getInterpolatedValue(Double_t x, Int_t * ArrayX, Double_t * ArrayY, Int_t ArraySize) { Double_t result; Double_t p0, p1, p2; result = getInterpolatedValue(x,ArrayX, ArrayY, ArraySize, p0, p1, p2 ); return result; } //____________________________________________________________________________________ ///This function returns interpolation result for given arrays and polynomial 2-nd order interpolation Double_t TKRATAnalyser::getInterpolatedValue(Double_t x, Int_t * ArrayX, Double_t * ArrayY, Int_t ArraySize, Double_t &A_p0, Double_t &A_p1 , Double_t &A_p2 ) { Int_t verboseLevel = 0; ///Allowed values:1,2,3,... ///1. Validate input data: if(ArraySize <= 0) { return 0.0; } ///2. Get the smallest and highest x value Int_t Xmin = ArrayX[0]; Int_t Xmax = ArrayX[0]; for(int i = 0; i < ArraySize; i++) { if(Xmin > ArrayX[i]){Xmin =ArrayX[i]; } if(Xmax < ArrayX[i]){Xmax =ArrayX[i]; } } ///3. Creating histogram to get access to fit procedure TH1F *histo = new TH1F("histo","Fitting", ArraySize +1 , Xmin , Xmax); for(int i = 0; i < ArraySize; i++) { histo->Fill(ArrayX[i], ArrayY[i]); //histo->SetBinContent(ArrayX[i], ArrayY[i]); if( verboseLevel > 0 ) cout << "("<< ArrayX[i] << "," << ArrayY[i] << "), "; } if( verboseLevel > 0 ) cout << endl; ///4. Fit polynomial of order 2. TFitResultPtr fitResults; fitResults = histo->Fit("pol2", "SFNWQ"); ///WNSFM, Q-quiet ,"SFNWQ" ///5. Getting results Double_t p0 = fitResults->Parameter(0); Double_t p1 = fitResults->Parameter(1); Double_t p2 = fitResults->Parameter(2); A_p0 = p0; A_p1 = p1; A_p2 = p2; if( verboseLevel > 0 ) cout << "Parameters: p0="<< p0 << ", p1="<< p1 << ",p2=" << p2 << endl; Double_t y = p2*x*x + p1*x + p0; //fitResults->Print("V"); // print full information of fit including covariance matrix if(verboseLevel >0) cout << "Result: ("<< x << "," << y << ")"<Draw(); ///6. Cleaning delete histo; ///7. Returning result return y; } //____________________________________________________________________________________ vector TKRATAnalyser::getPrecentsInDumbWay ( TH1F * A_hist, Int_t x0, Double_t baseline ) { vector points; Int_t xMax = A_hist->GetMaximumBin(); Int_t yMax = A_hist->GetMaximum(); Double_t level[5] = {10., 30., 50., 70., 90.}; Double_t height = (yMax - baseline); Int_t curX; Double_t curY; for ( int lvl = 0; lvl < 5; lvl++) { for ( int i = 0; i < xMax; i++) { curX = x0 + i; curY = A_hist -> GetBinContent( curX ); cout << " Point " << curX << ", "<< curY << endl; if ( curY >= baseline + height * ( level [lvl] /100.)) { cout << "Found ["<<(level[lvl] ) << "] for y=" << curY << " x=" << curX << endl; points.push_back( curX ); points.push_back( curY ); break; } } } return points; } //____________________________________________________________________________________ vector TKRATAnalyser::getPrecents ( TH1F * A_hist, Int_t x0, Double_t baseline ) { vector points; //points.push_back(1); Int_t xMax = A_hist->GetMaximumBin(); Int_t yMax = A_hist->GetMaximum(); Double_t level[5] = {10., 30., 50., 70., 90.}; Double_t height = (yMax - baseline); Int_t curX; Double_t curY; Double_t searchedX; Double_t searchedY; Int_t tabX[10]; Double_t tabY[10]; for ( int lvl = 0; lvl < 5; lvl++) { for ( int i = 0; i < xMax; i++) { curX = x0 + i; curY = A_hist -> GetBinContent( curX ); ///cout << " Point " << curX << ", "<< curY << endl; searchedY = baseline + height * ( level [lvl] /100.); if ( curY >= searchedY) { cout << "Znaleziono przyblizone dla "<<(level[lvl] ) << " % x=" << curX << ": y=" << curY << endl; for(int ti = 0; ti < 10; ti++) { tabX[ti] = curX - 5 + ti; tabY[ti] = A_hist -> GetBinContent( tabX[ti] ); } Double_t a,b,c; getInterpolatedValue( curX, tabX,tabY, 10, c, b, a ); ///Jezeli wierzcholek jest na dole to prawy punkt Double_t x1; Double_t x2; Double_t delta = b*b - 4* (a * (c - searchedY)); if (a != 0 && delta>= 0){ x1 = (-b - TMath:: Sqrt( delta ) ) /(2 * a); x2 = (-b + TMath:: Sqrt( delta ) ) /(2 * a); } else{ x1 = x2 = curX; } ///Jezeli wierzcholek jest do gory to lewy punkt ///lub wybrac wierzcholek blizej binu cout << " -odleglosci od curX: dla x1=" << x1 <<" : " << TMath::Abs(curX - x1) << " i dla x2=" << x2 << " : " << TMath::Abs(curX - x2) << endl; if ( TMath::Abs(curX - x1) < TMath::Abs(curX - x2)) { searchedX = x1; } else { searchedX = x2; } cout << " -Wybrany: " << searchedX << " dla " << searchedY << endl; points.push_back( searchedX ); points.push_back( searchedY ); break; } } } return points; } //______________________________________________________________________ void TKRATAnalyser::test ( void ) { return; }