+
+
+
+/////////////////
+// Array tools //
+/////////////////
+
+
+Int_t AliTPCCalibViewer::GetBin(Float_t value, Int_t nbins, Double_t binLow, Double_t binUp){
+ // Returns the 'bin' for 'value'
+ // The interval between 'binLow' and 'binUp' is divided into 'nbins' equidistant bins
+ // avoid index out of bounds error: 'if (bin < binLow) bin = binLow' and vice versa
+ /* Begin_Latex
+ GetBin(value) = #frac{nbins - 1}{binUp - binLow} #upoint (value - binLow) +1
+ End_Latex
+ */
+
+ Int_t bin = TMath::Nint( (Float_t)(value - binLow) / (Float_t)(binUp - binLow) * (nbins-1) ) + 1;
+ // avoid index out of bounds:
+ if (value < binLow) bin = 0;
+ if (value > binUp) bin = nbins + 1;
+ return bin;
+
+}
+
+
+Double_t AliTPCCalibViewer::GetLTM(Int_t n, const Double_t *const array, Double_t *const sigma, Double_t fraction){
+ //
+ // returns the LTM and sigma
+ //
+ Double_t *ddata = new Double_t[n];
+ Double_t mean = 0, lsigma = 0;
+ UInt_t nPoints = 0;
+ for (UInt_t i = 0; i < (UInt_t)n; i++) {
+ ddata[nPoints]= array[nPoints];
+ nPoints++;
+ }
+ Int_t hh = TMath::Min(TMath::Nint(fraction * nPoints), Int_t(n));
+ AliMathBase::EvaluateUni(nPoints, ddata, mean, lsigma, hh);
+ if (sigma) *sigma = lsigma;
+ delete [] ddata;
+ return mean;
+}
+
+
+TH1F* AliTPCCalibViewer::SigmaCut(TH1F *const histogram, Float_t mean, Float_t sigma, Float_t sigmaMax, Float_t sigmaStep, Bool_t pm) {
+ //
+ // Creates a cumulative histogram Begin_Latex S(t, #mu, #sigma) End_Latex, where you can see, how much of the data are inside sigma-intervals around the mean value
+ // The data of the distribution Begin_Latex f(x, #mu, #sigma) End_Latex are given in 'histogram'
+ // 'mean' and 'sigma' are Begin_Latex #mu End_Latex and Begin_Latex #sigma End_Latex of the distribution in 'histogram', to be specified by the user
+ // sigmaMax: up to which sigma around the mean/median/LTM the histogram is generated (in units of sigma, Begin_Latex t #sigma End_Latex)
+ // sigmaStep: the binsize of the generated histogram, -1 means, that the maximal reasonable stepsize is used
+ // pm: Decide weather Begin_Latex t > 0 End_Latex (first case) or Begin_Latex t End_Latex arbitrary (secound case)
+ // The actual work is done on the array.
+ /* Begin_Latex
+ f(x, #mu, #sigma) #Rightarrow S(t, #mu, #sigma) = (#int_{#mu}^{#mu + t #sigma} f(x, #mu, #sigma) dx + #int_{#mu}^{#mu - t #sigma} f(x, #mu, #sigma) dx) / (#int_{-#infty}^{+#infty} f(x, #mu, #sigma) dx), for t > 0
+ or
+ f(x, #mu, #sigma) #Rightarrow S(t, #mu, #sigma) = #int_{#mu}^{#mu + t #sigma} f(x, #mu, #sigma) dx / #int_{-#infty}^{+#infty} f(x, #mu, #sigma) dx
+ End_Latex
+ Begin_Macro(source)
+ {
+ Float_t mean = 0;
+ Float_t sigma = 1.5;
+ Float_t sigmaMax = 4;
+ gROOT->SetStyle("Plain");
+ TH1F *distribution = new TH1F("Distrib1", "Distribution f(x, #mu, #sigma)", 1000,-5,5);
+ TRandom rand(23);
+ for (Int_t i = 0; i <50000;i++) distribution->Fill(rand.Gaus(mean, sigma));
+ Float_t *ar = distribution->GetArray();
+
+ TCanvas* macro_example_canvas = new TCanvas("cAliTPCCalibViewer1", "", 350, 350);
+ macro_example_canvas->Divide(0,3);
+ TVirtualPad *pad1 = macro_example_canvas->cd(1);
+ pad1->SetGridy();
+ pad1->SetGridx();
+ distribution->Draw();
+ TVirtualPad *pad2 = macro_example_canvas->cd(2);
+ pad2->SetGridy();
+ pad2->SetGridx();
+
+ TH1F *shist = AliTPCCalibViewer::SigmaCut(distribution, mean, sigma, sigmaMax);
+ shist->SetNameTitle("Cumulative","Cumulative S(t, #mu, #sigma)");
+ shist->Draw();
+ TVirtualPad *pad3 = macro_example_canvas->cd(3);
+ pad3->SetGridy();
+ pad3->SetGridx();
+ TH1F *shistPM = AliTPCCalibViewer::SigmaCut(distribution, mean, sigma, sigmaMax, -1, kTRUE);
+ shistPM->Draw();
+ return macro_example_canvas;
+ }
+ End_Macro
+ */
+
+ Float_t *array = histogram->GetArray();
+ Int_t nbins = histogram->GetXaxis()->GetNbins();
+ Float_t binLow = histogram->GetXaxis()->GetXmin();
+ Float_t binUp = histogram->GetXaxis()->GetXmax();
+ return AliTPCCalibViewer::SigmaCut(nbins, array, mean, sigma, nbins, binLow, binUp, sigmaMax, sigmaStep, pm);
+}
+
+
+TH1F* AliTPCCalibViewer::SigmaCut(Int_t n, const Float_t *array, Float_t mean, Float_t sigma, Int_t nbins, Float_t binLow, Float_t binUp, Float_t sigmaMax, Float_t sigmaStep, Bool_t pm){
+ //
+ // Creates a histogram Begin_Latex S(t, #mu, #sigma) End_Latex, where you can see, how much of the data are inside sigma-intervals around the mean value
+ // The data of the distribution Begin_Latex f(x, #mu, #sigma) End_Latex are given in 'array', 'n' specifies the length of the array
+ // 'mean' and 'sigma' are Begin_Latex #mu End_Latex and Begin_Latex #sigma End_Latex of the distribution in 'array', to be specified by the user
+ // 'nbins': number of bins, 'binLow': first bin, 'binUp': last bin
+ // sigmaMax: up to which sigma around the mean/median/LTM the histogram is generated (in units of sigma, Begin_Latex t #sigma End_Latex)
+ // sigmaStep: the binsize of the generated histogram
+ // Here the actual work is done.
+
+ if (sigma == 0) return 0;
+ Float_t binWidth = (binUp-binLow)/(nbins - 1);
+ if (sigmaStep <= 0) sigmaStep = binWidth;
+ Int_t kbins = (Int_t)(sigmaMax * sigma / sigmaStep) + 1; // + 1 due to overflow bin in histograms
+ if (pm) kbins = 2 * (Int_t)(sigmaMax * sigma / sigmaStep) + 1;
+ Float_t kbinLow = !pm ? 0 : -sigmaMax;
+ Float_t kbinUp = sigmaMax;
+ TH1F *hist = new TH1F("sigmaCutHisto","Cumulative; Multiples of #sigma; Fraction of included data", kbins, kbinLow, kbinUp);
+ hist->SetDirectory(0);
+ hist->Reset();
+
+ // calculate normalization
+ Double_t normalization = 0;
+ for (Int_t i = 0; i <= n; i++) {
+ normalization += array[i];
+ }
+
+ // given units: units from given histogram
+ // sigma units: in units of sigma
+ // iDelta: integrate in interval (mean +- iDelta), given units
+ // x: ofset from mean for integration, given units
+ // hist: needs
+
+// printf("nbins: %i, binLow: %f, binUp: %f \n", nbins, binLow, binUp);
+ // fill histogram
+ for (Float_t iDelta = 0; iDelta <= sigmaMax * sigma; iDelta += sigmaStep) {
+ // integrate array
+ Double_t valueP = array[GetBin(mean, nbins, binLow, binUp)];
+ Double_t valueM = array[GetBin(mean-binWidth, nbins, binLow, binUp)];
+ // add bin of mean value only once to the histogram
+// printf("++ adding bins: ");
+ for (Float_t x = binWidth; x <= iDelta; x += binWidth) {
+ valueP += (mean + x <= binUp) ? array[GetBin(mean + x, nbins, binLow, binUp)] : 0;
+ valueM += (mean-binWidth - x >= binLow) ? array[GetBin(mean-binWidth - x, nbins, binLow, binUp)] : 0;
+// printf("%i, ", GetBin(mean + x, nbins, binLow, binUp));
+ }
+// printf("\n");
+ if (valueP / normalization > 100) printf("+++ Error, value to big: %f, normalization with %f will fail +++ \n", valueP, normalization);
+ if (valueP / normalization > 100) return hist;
+ if (valueM / normalization > 100) printf("+++ Error, value to big: %f, normalization with %f will fail +++ \n", valueM, normalization);
+ if (valueM / normalization > 100) return hist;
+ valueP = (valueP / normalization);
+ valueM = (valueM / normalization);
+ if (pm) {
+ Int_t bin = GetBin(iDelta/sigma, kbins, kbinLow, kbinUp);
+ hist->SetBinContent(bin, valueP);
+ bin = GetBin(-iDelta/sigma, kbins, kbinLow, kbinUp);
+ hist->SetBinContent(bin, valueM);
+ }
+ else { // if (!pm)
+ Int_t bin = GetBin(iDelta/sigma, kbins, kbinLow, kbinUp);
+ hist->SetBinContent(bin, valueP + valueM);
+// printf(" first integration bin: %i, last integration bin in + direction: %i \n", GetBin(mean+binWidth, nbins, binLow, binUp), GetBin(iDelta, nbins, binLow, binUp));
+// printf(" first integration bin: %i, last integration bin in - direction: %i \n", GetBin(mean+binWidth, nbins, binLow, binUp), GetBin(-iDelta, nbins, binLow, binUp));
+// printf(" value: %f, normalization: %f, iDelta: %f, Bin: %i \n", valueP+valueM, normalization, iDelta, bin);
+ }
+ }
+ //hist->SetMaximum(0.7);
+ if (!pm) hist->SetMaximum(1.2);
+ return hist;
+}
+
+
+TH1F* AliTPCCalibViewer::SigmaCut(Int_t /*n*/, const Double_t */*array*/, Double_t /*mean*/, Double_t /*sigma*/, Int_t /*nbins*/, const Double_t */*xbins*/, Double_t /*sigmaMax*/){
+ //
+ // SigmaCut for variable binsize
+ // NOT YET IMPLEMENTED !!!
+ //
+ printf("SigmaCut with variable binsize, Not yet implemented\n");
+
+ return 0;
+}
+
+
+TH1F* AliTPCCalibViewer::Integrate(TH1F *const histogram, Float_t mean, Float_t sigma, Float_t sigmaMax, Float_t sigmaStep){
+ //
+ // Creates an integrated histogram Begin_Latex S(t, #mu, #sigma) End_Latex, out of the input distribution distribution Begin_Latex f(x, #mu, #sigma) End_Latex, given in "histogram"
+ // "mean" and "sigma" are Begin_Latex #mu End_Latex and Begin_Latex #sigma End_Latex of the distribution in "histogram", to be specified by the user
+ // sigmaMax: up to which sigma around the mean/median/LTM you want to integrate
+ // if "igma == 0" and "sigmaMax == 0" the whole histogram is integrated
+ // "sigmaStep": the binsize of the generated histogram, -1 means, that the maximal reasonable stepsize is used
+ // The actual work is done on the array.
+ /* Begin_Latex
+ f(x, #mu, #sigma) #Rightarrow S(t, #mu, #sigma) = #int_{-#infty}^{#mu + t #sigma} f(x, #mu, #sigma) dx / #int_{-#infty}^{+#infty} f(x, #mu, #sigma) dx
+ End_Latex
+ Begin_Macro(source)
+ {
+ Float_t mean = 0;
+ Float_t sigma = 1.5;
+ Float_t sigmaMax = 4;
+ gROOT->SetStyle("Plain");
+ TH1F *distribution = new TH1F("Distrib2", "Distribution f(x, #mu, #sigma)", 1000,-5,5);
+ TRandom rand(23);
+ for (Int_t i = 0; i <50000;i++) distribution->Fill(rand.Gaus(mean, sigma));
+ Float_t *ar = distribution->GetArray();
+
+ TCanvas* macro_example_canvas = new TCanvas("cAliTPCCalibViewer2", "", 350, 350);
+ macro_example_canvas->Divide(0,2);
+ TVirtualPad *pad1 = macro_example_canvas->cd(1);
+ pad1->SetGridy();
+ pad1->SetGridx();
+ distribution->Draw();
+ TVirtualPad *pad2 = macro_example_canvas->cd(2);
+ pad2->SetGridy();
+ pad2->SetGridx();
+ TH1F *shist = AliTPCCalibViewer::Integrate(distribution, mean, sigma, sigmaMax);
+ shist->SetNameTitle("Cumulative","Cumulative S(t, #mu, #sigma)");
+ shist->Draw();
+
+ }
+ End_Macro
+ */
+
+
+ Float_t *array = histogram->GetArray();
+ Int_t nbins = histogram->GetXaxis()->GetNbins();
+ Float_t binLow = histogram->GetXaxis()->GetXmin();
+ Float_t binUp = histogram->GetXaxis()->GetXmax();
+ return AliTPCCalibViewer::Integrate(nbins, array, nbins, binLow, binUp, mean, sigma, sigmaMax, sigmaStep);
+}
+
+
+TH1F* AliTPCCalibViewer::Integrate(Int_t n, const Float_t *const array, Int_t nbins, Float_t binLow, Float_t binUp, Float_t mean, Float_t sigma, Float_t sigmaMax, Float_t sigmaStep){
+ // Creates an integrated histogram Begin_Latex S(t, #mu, #sigma) End_Latex, out of the input distribution distribution Begin_Latex f(x, #mu, #sigma) End_Latex, given in "histogram"
+ // "mean" and "sigma" are Begin_Latex #mu End_Latex and Begin_Latex #sigma End_Latex of the distribution in "histogram", to be specified by the user
+ // sigmaMax: up to which sigma around the mean/median/LTM you want to integrate
+ // if "igma == 0" and "sigmaMax == 0" the whole histogram is integrated
+ // "sigmaStep": the binsize of the generated histogram, -1 means, that the maximal reasonable stepsize is used
+ // Here the actual work is done.
+
+ Bool_t givenUnits = kTRUE;
+ if (sigma != 0 && sigmaMax != 0) givenUnits = kFALSE;
+ if (givenUnits) {
+ sigma = 1;
+ sigmaMax = (binUp - binLow) / 2.;
+ }
+
+ Float_t binWidth = (binUp-binLow)/(nbins - 1);
+ if (sigmaStep <= 0) sigmaStep = binWidth;
+ Int_t kbins = (Int_t)(sigmaMax * sigma / sigmaStep) + 1; // + 1 due to overflow bin in histograms
+ Float_t kbinLow = givenUnits ? binLow : -sigmaMax;
+ Float_t kbinUp = givenUnits ? binUp : sigmaMax;
+ TH1F *hist = 0;
+ if (givenUnits) hist = new TH1F("integratedHisto","Integrated Histogram; Given x; Fraction of included data", kbins, kbinLow, kbinUp);
+ if (!givenUnits) hist = new TH1F("integratedHisto","Integrated Histogram; Multiples of #sigma; Fraction of included data", kbins, kbinLow, kbinUp);
+ hist->SetDirectory(0);
+ hist->Reset();
+
+ // calculate normalization
+ // printf("calculating normalization, integrating from bin 1 to %i \n", n);
+ Double_t normalization = 0;
+ for (Int_t i = 1; i <= n; i++) {
+ normalization += array[i];
+ }
+ // printf("normalization: %f \n", normalization);
+
+ // given units: units from given histogram
+ // sigma units: in units of sigma
+ // iDelta: integrate in interval (mean +- iDelta), given units
+ // x: ofset from mean for integration, given units
+ // hist: needs
+
+ // fill histogram
+ for (Float_t iDelta = mean - sigmaMax * sigma; iDelta <= mean + sigmaMax * sigma; iDelta += sigmaStep) {
+ // integrate array
+ Double_t value = 0;
+ for (Float_t x = mean - sigmaMax * sigma; x <= iDelta; x += binWidth) {
+ value += (x <= binUp && x >= binLow) ? array[GetBin(x, nbins, binLow, binUp)] : 0;
+ }
+ if (value / normalization > 100) printf("+++ Error, value to big: %f, normalization with %f will fail +++ \n", value, normalization);
+ if (value / normalization > 100) return hist;
+ Int_t bin = GetBin(iDelta/sigma, kbins, kbinLow, kbinUp);
+ // printf("first integration bin: %i, last integration bin: %i \n", GetBin(mean - sigmaMax * sigma, nbins, binLow, binUp), GetBin(iDelta, nbins, binLow, binUp));
+ // printf("value: %f, normalization: %f, normalized value: %f, iDelta: %f, Bin: %i \n", value, normalization, value/normalization, iDelta, bin);
+ value = (value / normalization);
+ hist->SetBinContent(bin, value);
+ }
+ return hist;
+}
+
+
+
+
+
+////////////////////////
+// end of Array tools //
+////////////////////////
+
+
+