#include <TVectorD.h>
#include <TROOT.h>
#include <TString.h>
+#include <TLine.h>
#include "AliLog.h"
#include "AliMathBase.h"
,fNumberOfBinsExpected(0)
,fMethod(0)
,fBeginFitCharge(3.5)
+ ,fOutliersFitChargeLow(0.03)
+ ,fOutliersFitChargeHigh(0.80)
,fFitPHPeriode(1)
,fTakeTheMaxPH(kTRUE)
,fT0Shift0(0.124797)
,fNumberOfBinsExpected(c.fNumberOfBinsExpected)
,fMethod(c.fMethod)
,fBeginFitCharge(c.fBeginFitCharge)
+,fOutliersFitChargeLow(c.fOutliersFitChargeLow)
+,fOutliersFitChargeHigh(c.fOutliersFitChargeHigh)
,fFitPHPeriode(c.fFitPHPeriode)
,fTakeTheMaxPH(c.fTakeTheMaxPH)
,fT0Shift0(c.fT0Shift0)
if ( fDebugStreamer ) delete fDebugStreamer;
fDebugStreamer = 0x0;
-}
-//__________________________________________________________________________________
-void AliTRDCalibraFit::RangeChargeIntegration(Float_t vdrift, Float_t t0, Int_t &begin, Int_t &peak, Int_t &end) const
-{
- //
- // From the drift velocity and t0
- // return the position of the peak and maximum negative slope
- //
-
- const Float_t kDrWidth = AliTRDgeometry::DrThick(); // drift region
- Double_t widbins = 0.1; // 0.1 mus
-
- //peak and maxnegslope in mus
- Double_t begind = t0*widbins + fT0Shift0;
- Double_t peakd = t0*widbins + fT0Shift1;
- Double_t maxnegslope = (kDrWidth + vdrift*peakd)/vdrift;
-
- // peak and maxnegslope in timebin
- begin = TMath::Nint(begind*widbins);
- peak = TMath::Nint(peakd*widbins);
- end = TMath::Nint(maxnegslope*widbins);
-
}
//____________Functions fit Online CH2d________________________________________
Bool_t AliTRDCalibraFit::AnalyseCH(const TH2I *ch)
{
case 0: FitMeanW((TH1 *) projch, nentries); break;
case 1: FitMean((TH1 *) projch, nentries, mean); break;
- case 2: FitCH((TH1 *) projch, mean); break;
- case 3: FitBisCH((TH1 *) projch, mean); break;
+ case 2: FitLandau((TH1 *) projch, mean, nentries); break;
+ case 3: FitCH((TH1 *) projch, mean, nentries); break;
+ case 4: FitBisCH((TH1 *) projch, mean, nentries); break;
+ case 5: FitBisCHEx((TH1 *) projch, mean, nentries); break;
default: return kFALSE;
}
// Fill Infos Fit
}
// Mean Statistic
if (fNumberFit > 0) {
- AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit, fNumberFitSuccess));
+ AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %f over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Double_t) (fStatisticMean/fNumberFit), fNumberFitSuccess));
fStatisticMean = fStatisticMean / fNumberFit;
}
else {
{
case 0: FitMeanW((TH1 *) projch, nentries); break;
case 1: FitMean((TH1 *) projch, nentries, mean); break;
- case 2: FitCH((TH1 *) projch, mean); break;
- case 3: FitBisCH((TH1 *) projch, mean); break;
+ case 2: FitLandau((TH1 *) projch, mean, nentries); break;
+ case 3: FitCH((TH1 *) projch, mean, nentries); break;
+ case 4: FitBisCH((TH1 *) projch, mean, nentries); break;
+ case 5: FitBisCHEx((TH1 *) projch, mean, nentries); break;
default: return kFALSE;
}
// Fill Infos Fit
}
// Mean Statistics
if (fNumberFit > 0) {
- AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit, fNumberFitSuccess));
+ AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %f over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Double_t) (fStatisticMean/fNumberFit), fNumberFitSuccess));
fStatisticMean = fStatisticMean / fNumberFit;
}
else {
{
case 0: FitMeanW((TH1 *) projch, nentries); break;
case 1: FitMean((TH1 *) projch, nentries, mean); break;
- case 2: FitCH((TH1 *) projch, mean); break;
- case 3: FitBisCH((TH1 *) projch, mean); break;
+ case 2: FitLandau((TH1 *) projch, mean, nentries); break;
+ case 3: FitCH((TH1 *) projch, mean, nentries); break;
+ case 4: FitBisCH((TH1 *) projch, mean, nentries); break;
+ case 5: FitBisCHEx((TH1 *) projch, mean, nentries); break;
default: return -100.0;
}
delete fDebugStreamer;
// Mean Statistic
if (fNumberFit > 0) {
- AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
+ AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %f over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Double_t) (fStatisticMean/fNumberFit),fNumberFitSuccess));
fStatisticMean = fStatisticMean / fNumberFit;
}
else {
} // Boucle object
// Mean Statistic
if (fNumberFit > 0) {
- AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
+ AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %f over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Double_t) (fStatisticMean/fNumberFit),fNumberFitSuccess));
fStatisticMean = fStatisticMean / fNumberFit;
}
else {
if (fNumberFit > 0) {
AliInfo(Form("There are %d with at least one entries.",fNumberEnt));
AliInfo(Form("%d fits have been proceeded (sucessfully or not...).",fNumberFit));
- AliInfo(Form("There is a mean statistic of: %d over these fitted histograms and %d successfulled fits"
- ,(Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
+ AliInfo(Form("There is a mean statistic of: %f over these fitted histograms and %d successfulled fits"
+ ,(Double_t) (fStatisticMean/fNumberFit),fNumberFitSuccess));
fStatisticMean = fStatisticMean / fNumberFit;
}
else {
if(!SetModeCalibration(name,2)) return kFALSE;
// Number of Ybins (detectors or groups of pads)
- TAxis *xprf = prf->GetXaxis();
- TAxis *yprf = prf->GetYaxis();
+ const TAxis *xprf = prf->GetXaxis();
+ const TAxis *yprf = prf->GetYaxis();
Int_t nybins = yprf->GetNbins();// calibration groups
Int_t nbins = xprf->GetNbins();// bins
Float_t lowedge = (Float_t) xprf->GetBinLowEdge(1);//lowedge in bins
if (fNumberFit > 0) {
AliInfo(Form("There are %d with at least one entries.",fNumberEnt));
AliInfo(Form("%d fits have been proceeded (sucessfully or not...).",fNumberFit));
- AliInfo(Form("There is a mean statistic of: %d over these fitted histograms and %d successfulled fits"
- ,(Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
+ AliInfo(Form("There is a mean statistic of: %f over these fitted histograms and %d successfulled fits"
+ ,(Double_t) (fStatisticMean/fNumberFit),fNumberFitSuccess));
fStatisticMean = fStatisticMean / fNumberFit;
}
else {
} // Boucle object
// Mean Statistics
if (fNumberFit > 0) {
- AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
+ AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %f over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Double_t) (fStatisticMean/fNumberFit),fNumberFitSuccess));
}
else {
AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
} // Boucle object
// Mean Statistics
if (fNumberFit > 0) {
- AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
+ AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %f over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Double_t) (fStatisticMean/fNumberFit),fNumberFitSuccess));
}
else {
AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
// Take the result
TVectorD param(2);
TVectorD error(3);
- fEntriesCurrent = 0;
+ Double_t entriesCurrent = 0;
fCountDet = idet;
Bool_t here = calivdli->GetParam(idet,¶m);
Bool_t heree = calivdli->GetError(idet,&error);
//printf("here %d and heree %d\n",here, heree);
if(heree) {
- fEntriesCurrent = (Int_t) error[2];
+ entriesCurrent = error[2];
fNumberEnt++;
}
//printf("Number of entries %d\n",fEntriesCurrent);
// Nothing found or not enough statistic
- if((!heree) || (!here) || (fEntriesCurrent <= fMinEntries)) {
+ if((!heree) || (!here) || (entriesCurrent <= fMinEntries)) {
NotEnoughStatisticLinearFitter();
continue;
}
//error.Print();
//Statistics
fNumberFit++;
- fStatisticMean += fEntriesCurrent;
+ fStatisticMean += entriesCurrent;
// Check the fit
- if((-(param[1])) <= 0.0) {
+ if((-(param[1])) <= 0.000001) {
NotEnoughStatisticLinearFitter();
continue;
}
}
// Mean Statistics
if (fNumberFit > 0) {
- AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
+ AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %f over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Double_t) (fStatisticMean/fNumberFit),fNumberFitSuccess));
}
else {
AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
// Take the result
TVectorD param(3);
TVectorD error(3);
- fEntriesCurrent = 0;
+ Double_t entriesCurrent = 0;
fCountDet = idet;
Bool_t here = calivdli->GetParam(idet,¶m);
Bool_t heree = calivdli->GetError(idet,&error);
//printf("here %d and heree %d\n",here, heree);
if(heree) {
- fEntriesCurrent = (Int_t) error[2];
+ entriesCurrent = error[2];
fNumberEnt++;
}
//printf("Number of entries %d\n",fEntriesCurrent);
// Nothing found or not enough statistic
- if((!heree) || (!here) || (fEntriesCurrent <= fMinEntries)) {
+ if((!heree) || (!here) || (entriesCurrent <= fMinEntries)) {
NotEnoughStatisticExbAlt();
continue;
}
//error.Print();
//Statistics
fNumberFit++;
- fStatisticMean += fEntriesCurrent;
+ fStatisticMean += entriesCurrent;
// Statistics
fNumberFitSuccess ++;
}
// Mean Statistics
if (fNumberFit > 0) {
- AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess));
+ AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %f over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Double_t) (fStatisticMean/fNumberFit),fNumberFitSuccess));
}
else {
AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
// Fit
- Int_t entries = 0;
+ Double_t entries = 0;
TAxis *xaxis = linearfitterhisto->GetXaxis();
TAxis *yaxis = linearfitterhisto->GetYaxis();
TLinearFitter linearfitter = TLinearFitter(2,"pol1");
for(Int_t k = 0; k < (Int_t)linearfitterhisto->GetBinContent(ibinx+1,ibiny+1); k++){
if(!securitybreaking){
linearfitter.AddPoint(&x,y);
- entries++;
+ entries = entries+1.;
}
else {
- if(entries< 1198){
+ if(entries< 1198.0){
linearfitter.AddPoint(&x,y);
- entries++;
+ entries = entries + 1.;
}
}
}
// Put the fCurrentCoef
fCurrentCoef[0] = -par[1];
// here the database must be the one of the reconstruction for the lorentz angle....
- if(fCurrentCoef[0] > 0.0) fCurrentCoef2[0] = (par[0]+meanvdriftused*meanexbused)/fCurrentCoef[0];
+ if(fCurrentCoef[0] > 0.00001) fCurrentCoef2[0] = (par[0]+meanvdriftused*meanexbused)/fCurrentCoef[0];
else fCurrentCoef2[0] = 100.0;
}
return;
}
Int_t detector = -1;
- Int_t sector = -1;
+ // Coverity
+ //Int_t sector = -1;
Float_t value = 0.0;
/////////////////////////////////
// Initialisation
////////////////////////
Double_t meanAll = 0.0;
- Double_t rmsAll = 0.0;
- Int_t countAll = 0;
- ////////////
+ Double_t rmsAll = 0.0;
+ Int_t countAll = 0;
+ ////////////
// compute
////////////
for (Int_t k = 0; k < loop; k++) {
detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
- sector = GetSector(detector);
+ // Coverity
+ //sector = GetSector(detector);
if(perdetector){
value = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef()[0];
if(value > 0.0) {
if(type==1) defaultvalue = -1.5;
for (Int_t k = 0; k < loop; k++) {
detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
- sector = GetSector(detector);
+ // Coverity
+ //sector = GetSector(detector);
Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
Int_t colMax = fGeo->GetColMax(GetLayer(detector));
Float_t *coef = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef();
return;
}
Int_t detector = -1;
- Int_t sector = -1;
+ // Coverity
+ //Int_t sector = -1;
Float_t value = 0.0;
/////////////////////////////////
////////////
for (Int_t k = 0; k < loop; k++) {
detector = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetDetector();
- sector = GetSector(detector);
+ // Coverity
+ //sector = GetSector(detector);
if(perdetector){
value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[0];
if(value < 70.0) {
////////////////////////////////////////////////
for (Int_t k = 0; k < loop; k++) {
detector = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetDetector();
- sector = GetSector(detector);
+ // Coverity
+ //sector = GetSector(detector);
Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
Int_t colMax = fGeo->GetColMax(GetLayer(detector));
Float_t *coef = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef();
gDirectory = gROOT;
fScaleFitFactor = 0.0;
+ if( fCurrentCoefDetector ) delete [] fCurrentCoefDetector;
fCurrentCoefDetector = new Float_t[2304];
for (Int_t k = 0; k < 2304; k++) {
fCurrentCoefDetector[k] = 0.0;
fVectorFit.SetName("driftvelocitycoefficients");
fVectorFit2.SetName("t0coefficients");
+ if( fCurrentCoefDetector ) delete [] fCurrentCoefDetector;
fCurrentCoefDetector = new Float_t[2304];
for (Int_t k = 0; k < 2304; k++) {
fCurrentCoefDetector[k] = 0.0;
}
-
+ if( fCurrentCoefDetector2 ) delete [] fCurrentCoefDetector2;
fCurrentCoefDetector2 = new Float_t[2304];
for (Int_t k = 0; k < 2304; k++) {
fCurrentCoefDetector2[k] = 0.0;
gDirectory = gROOT;
fVectorFit.SetName("prfwidthcoefficients");
+ if( fCurrentCoefDetector ) delete [] fCurrentCoefDetector;
fCurrentCoefDetector = new Float_t[2304];
for (Int_t k = 0; k < 2304; k++) {
fCurrentCoefDetector[k] = 0.0;
gDirectory = gROOT;
+ if( fCurrentCoefDetector ) delete [] fCurrentCoefDetector;
+ if( fCurrentCoefDetector2 ) delete [] fCurrentCoefDetector2;
fCurrentCoefDetector = new Float_t[2304];
fCurrentCoefDetector2 = new Float_t[2304];
for (Int_t k = 0; k < 2304; k++) {
gDirectory = gROOT;
- fCurrentCoefDetector = new Float_t[2304];
+ if( fCurrentCoefDetector2 ) delete [] fCurrentCoefDetector2;
+ fCurrentCoefDetector2 = new Float_t[2304];
for (Int_t k = 0; k < 2304; k++) {
fCurrentCoefDetector2[k] = 0.0;
}
//
// End of one detector
- if ((idect == (fCount-1))) {
+ if (idect == (fCount-1)) {
FillVectorFit();
// Reset
for (Int_t k = 0; k < 2304; k++) {
//
// End of one detector
- if ((idect == (fCount-1))) {
+ if (idect == (fCount-1)) {
FillVectorFit();
FillVectorFit2();
// Reset
//
// End of one detector
- if ((idect == (fCount-1))) {
+ if (idect == (fCount-1)) {
FillVectorFit();
// Reset
for (Int_t k = 0; k < 2304; k++) {
Bool_t put = kTRUE;
+ ///////////////////////////////
// Beginning of the signal
+ //////////////////////////////
TH1D *pentea = new TH1D("pentea","pentea",projPH->GetNbinsX(),0,(Float_t) limit);
for (Int_t k = 1; k < projPH->GetNbinsX(); k++) {
pentea->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
fPhd[0] = placemaximum;
}
+ /////////////////////////////
// Amplification region
+ /////////////////////////////
binmax = 0;
ju = 0;
for (Int_t kbin = 1; kbin < projPH->GetNbinsX(); kbin ++) {
}
fPhd[1] = placemaximum;
}
-
+
+ //////////////////
// Drift region
+ //////////////////
+ Bool_t putd = kTRUE;
TH1D *pente = new TH1D("pente","pente", projPH->GetNbinsX(),0,(Float_t) limit);
for (Int_t k = TMath::Min(binmax+4, projPH->GetNbinsX()); k < projPH->GetNbinsX(); k++) {
pente->SetBinContent(k,(Double_t) (projPH->GetBinContent(k+1) - projPH->GetBinContent(k)));
//should not happen
if (binmin <= 1) {
binmin = 2;
- put = 1;
+ putd = 1;
//AliInfo("Put the binmax from 1 to 2 to enable the fit");
}
//check
if((projPH->GetBinContent(binmin)-projPH->GetBinError(binmin)) < (projPH->GetBinContent(binmin+1))) {
//AliInfo("Too many fluctuations at the end!");
- put = kFALSE;
+ putd = kFALSE;
}
if((projPH->GetBinContent(binmin)+projPH->GetBinError(binmin)) > (projPH->GetBinContent(binmin-1))) {
//AliInfo("Too many fluctuations at the end!");
- put = kFALSE;
+ putd = kFALSE;
}
if(TMath::Abs(pente->GetBinContent(binmin+1)) <= 0.0000000000001){
//AliInfo("No entries for the next bin!");
if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1)) &&
(pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
//AliInfo("polynome 4 false 1");
- put = kFALSE;
+ putd = kFALSE;
}
if(((binmin+3) <= (nbins-1)) &&
(pente->GetBinContent(binmin+3) <= pente->GetBinContent(binmin+2)) &&
((binmin-3) >= TMath::Min(binmax+4, projPH->GetNbinsX())) &&
(pente->GetBinContent(binmin-3) <= pente->GetBinContent(binmin-2))) {
//AliInfo("polynome 4 false 2");
- put = kFALSE;
+ putd = kFALSE;
}
// poly 3
if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1)) &&
//richtung +
if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1))) {
//AliInfo("polynome 2+ false");
- put = kFALSE;
+ putd = kFALSE;
}
}
//pol2 case 2
//richtung -
if((pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
//AliInfo("polynome 2- false ");
- put = kFALSE;
+ putd = kFALSE;
}
}
if((binmin == (nbins-1)) && ((binmin-2) < TMath::Min(binmax+4, projPH->GetNbinsX()))) {
- put = kFALSE;
+ putd = kFALSE;
//AliInfo("At the limit for the drift and not usable!");
}
//pass
if((binmin == (nbins-2)) && ((binmin-1) < TMath::Min(binmax+4, projPH->GetNbinsX()))){
- put = kFALSE;
+ putd = kFALSE;
//AliInfo("For the drift...problem!");
}
//pass but should not happen
if((binmin <= (nbins-3)) && (binmin < TMath::Min(binmax+6, projPH->GetNbinsX()))){
- put = kFALSE;
+ putd = kFALSE;
//AliInfo("For the drift...problem!");
}
- if(put) {
+ if(putd) {
polynome = new TF1("pol","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",min,max);
polynome->SetParameters(c0,c1,c2,c3,c4);
//AliInfo(Form("GetMinimum of the function %f",polynome->GetMinimumX()));
if ((fPhd[2] > fPhd[0]) &&
(fPhd[2] > fPhd[1]) &&
(fPhd[1] > fPhd[0]) &&
- (put)) {
+ (put) &&
+ (putd)) {
fCurrentCoef[0] = (kDrWidth) / (fPhd[2]-fPhd[1]);
if(fCurrentCoef[0] > 2.5) fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
else fNumberFitSuccess++;
if (fCurrentCoef2[0] < -3.0) {
fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
}
+ else fNumberFitSuccess++;
}
else {
fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
TCanvas *cpentei = new TCanvas("cpentei","cpentei",50,50,600,800);
cpentei->cd();
projPH->Draw();
+ if(polynomea) polynomea->Draw("same");
line->SetLineColor(2);
line->DrawLine(fPhd[0],0,fPhd[0],projPH->GetMaximum());
line->DrawLine(fPhd[1],0,fPhd[1],projPH->GetMaximum());
AliInfo(Form("fPhd[0] (beginning of the signal): %f" ,(Float_t) fPhd[0]));
AliInfo(Form("fPhd[1] (end of the amplification region): %f" ,(Float_t) fPhd[1]));
AliInfo(Form("fPhd[2] (end of the drift region): %f" ,(Float_t) fPhd[2]));
- AliInfo(Form("fVriftCoef[3] (with only the drift region(default)): %f",(Float_t) fCurrentCoef[0]));
+ AliInfo(Form("Vdrift (with only the drift region(default)): %f",(Float_t) fCurrentCoef[0]));
+ AliInfo(Form("Timeoffset (with only the drift region(default)): %f",(Float_t) fCurrentCoef2[0]));
TCanvas *cpentei2 = new TCanvas("cpentei2","cpentei2",50,50,600,800);
cpentei2->cd();
pentea->Draw();
TVectorD pars0;
linearfitter.Eval();
linearfitter.GetParameters(pars0);
- Double_t pointError0 = TMath::Sqrt(linearfitter.GetChisquare()/npoints);
- Double_t errorsx0 = linearfitter.GetParError(2)*pointError0;
+ // Coverity
+ //Double_t pointError0 = TMath::Sqrt(linearfitter.GetChisquare()/npoints);
+ //Double_t errorsx0 = linearfitter.GetParError(2)*pointError0;
Double_t min0 = 0.0;
- Double_t ermin0 = 0.0;
+ // Coverity
+ //Double_t ermin0 = 0.0;
//Double_t prfe0 = 0.0;
Double_t prf0 = 0.0;
if((pars0[2] > 0.000000000001) && (TMath::Abs(pars0[1]) >= 0.000000000001)) {
min0 = -pars0[1]/(2*pars0[2]);
- ermin0 = TMath::Abs(min0*(errorsx0/pars0[2]+linearfitter.GetParError(1)*pointError0/pars0[1]));
+ // Coverity
+ //ermin0 = TMath::Abs(min0*(errorsx0/pars0[2]+linearfitter.GetParError(1)*pointError0/pars0[1]));
prf0 = pars0[0]+pars0[1]*min0+pars0[2]*min0*min0;
if(prf0 > 0.0) {
/*
Int_t sumCurrent = 0;
for(Int_t k = 0; k <nybins; k++){
Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
- if (fraction>0.95) break;
+ if (fraction<fOutliersFitChargeLow) {
+ sumCurrent += (Int_t) projch->GetBinContent(k+1);
+ //printf("Take only after bin %d\n",k);
+ continue;
+ }
+ if (fraction>fOutliersFitChargeHigh) {
+ //printf("Break by the bin %d\n",k);
+ break;
+ }
Double_t weight = a + b*fraction + c*fraction*fraction + d *fraction*fraction*fraction+
e*fraction*fraction*fraction*fraction;
sumw += weight*projch->GetBinContent(k+1)*projch->GetBinCenter(k+1);
TCanvas *cpmeanw = new TCanvas("cpmeanw","cpmeanw",50,50,600,800);
cpmeanw->cd();
projch->Draw();
+ TLine *line = new TLine(fCurrentCoef[0],0.0,fCurrentCoef[0],20000.0);
+ line->Draw("same");
}
fNumberFitSuccess++;
CalculChargeCoefMean(kTRUE);
TCanvas *cpmeanw = new TCanvas("cpmeanw","cpmeanw",50,50,600,800);
cpmeanw->cd();
projch->Draw();
+ TLine *line = new TLine(fCurrentCoef[0],0.0,fCurrentCoef[0],20000.0);
+ line->Draw("same");
}
fNumberFitSuccess++;
}
//_____________________________________________________________________________
-void AliTRDCalibraFit::FitCH(TH1 *projch, Double_t mean)
+void AliTRDCalibraFit::FitLandau(TH1 *projch, Double_t mean, Double_t nentries)
+{
+ //
+ // Fit methode for the gain factor
+ //
+
+
+ //Calcul Range of the fit
+ Double_t lastvalue = 0.0;
+ Float_t sumAll = (Float_t) nentries;
+ Int_t sumCurrent = 0;
+ //printf("There are %d bins\n",nybins);
+ for(Int_t k = 0; k <projch->GetNbinsX(); k++){
+ Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
+ if (fraction>fOutliersFitChargeHigh) {
+ lastvalue = projch->GetBinCenter(k+1);
+ //printf("Break by %f\n",lastvalue);
+ break;
+ }
+ sumCurrent += (Int_t) projch->GetBinContent(k+1);
+ }
+ //
+
+ fCurrentCoef[0] = 0.0;
+ fCurrentCoefE = 0.0;
+ // Coverity
+ //Double_t chisqrl = 0.0;
+
+ projch->Fit("landau","WWQ+",""
+ ,(Double_t) mean/fBeginFitCharge
+ ,lastvalue);
+ // Coverity
+ //chisqrl = projch->GetFunction("landau")->GetChisquare();
+
+ if (fDebugLevel == 1) {
+ TCanvas *cp = new TCanvas("cp","cp",50,50,600,800);
+ cp->cd();
+ projch->Draw();
+ TLine *line = new TLine( projch->GetFunction("landau")->GetParameter(1),0.0,projch->GetFunction("landau")->GetParameter(1),20000.0);
+ line->Draw("same");
+ }
+
+ if ((projch->GetFunction("landau")->GetParameter(1) > 0) && (projch->GetFunction("landau")->GetParError(1) < (0.05*projch->GetFunction("landau")->GetParameter(1)))) {
+ fNumberFitSuccess++;
+ CalculChargeCoefMean(kTRUE);
+ fCurrentCoef[0] = projch->GetFunction("landau")->GetParameter(1);
+ fCurrentCoefE = projch->GetFunction("landau")->GetParError(1);
+ }
+ else {
+ CalculChargeCoefMean(kFALSE);
+ fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
+ }
+
+
+
+}
+//_____________________________________________________________________________
+void AliTRDCalibraFit::FitCH(TH1 *projch, Double_t mean, Double_t nentries)
{
//
// Fit methode for the gain factor
//
+
+ //Calcul Range of the fit
+ Double_t lastvalue = 0.0;
+ Float_t sumAll = (Float_t) nentries;
+ Int_t sumCurrent = 0;
+ //printf("There are %d bins\n",nybins);
+ for(Int_t k = 0; k <projch->GetNbinsX(); k++){
+ Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
+ if (fraction>fOutliersFitChargeHigh) {
+ lastvalue = projch->GetBinCenter(k+1);
+ //printf("Break by %f\n",lastvalue);
+ break;
+ }
+ sumCurrent += (Int_t) projch->GetBinContent(k+1);
+ }
+ //
fCurrentCoef[0] = 0.0;
fCurrentCoefE = 0.0;
Double_t chisqrl = 0.0;
Double_t chisqrg = 0.0;
Double_t chisqr = 0.0;
- TF1 *fLandauGaus = new TF1("fLandauGaus",FuncLandauGaus,0,300,5);
+ TF1 *fLandauGaus = new TF1("fLandauGaus",FuncLandauGaus,(Double_t) mean/fBeginFitCharge,lastvalue,5);
- projch->Fit("landau","0",""
+ projch->Fit("landau","WWQ0",""
,(Double_t) mean/fBeginFitCharge
- ,projch->GetBinCenter(projch->GetNbinsX()));
+ ,lastvalue);
Double_t l3P0 = projch->GetFunction("landau")->GetParameter(0);
Double_t l3P1 = projch->GetFunction("landau")->GetParameter(1);
Double_t l3P2 = projch->GetFunction("landau")->GetParameter(2);
chisqrl = projch->GetFunction("landau")->GetChisquare();
- projch->Fit("gaus","0",""
+ projch->Fit("gaus","WWQ0",""
,(Double_t) mean/fBeginFitCharge
- ,projch->GetBinCenter(projch->GetNbinsX()));
+ ,lastvalue);
Double_t g3P0 = projch->GetFunction("gaus")->GetParameter(0);
Double_t g3P2 = projch->GetFunction("gaus")->GetParameter(2);
chisqrg = projch->GetFunction("gaus")->GetChisquare();
fLandauGaus->SetParameters(l3P0,l3P1,l3P2,g3P0,g3P2);
if (fDebugLevel != 1) {
- projch->Fit("fLandauGaus","0",""
+ projch->Fit("fLandauGaus","WWQ0",""
,(Double_t) mean/fBeginFitCharge
- ,projch->GetBinCenter(projch->GetNbinsX()));
+ ,lastvalue);
chisqr = projch->GetFunction("fLandauGaus")->GetChisquare();
}
else {
TCanvas *cp = new TCanvas("cp","cp",50,50,600,800);
cp->cd();
- projch->Fit("fLandauGaus","+",""
+ projch->Fit("fLandauGaus","WWQ+",""
,(Double_t) mean/fBeginFitCharge
- ,projch->GetBinCenter(projch->GetNbinsX()));
+ ,lastvalue);
chisqr = projch->GetFunction("fLandauGaus")->GetChisquare();
projch->Draw();
fLandauGaus->Draw("same");
+ TLine *line = new TLine(projch->GetFunction("fLandauGaus")->GetParameter(1),0.0,projch->GetFunction("fLandauGaus")->GetParameter(1),20000.0);
+ line->Draw("same");
}
if ((projch->GetFunction("fLandauGaus")->GetParameter(1) > 0) && (projch->GetFunction("fLandauGaus")->GetParError(1) < (0.05*projch->GetFunction("fLandauGaus")->GetParameter(1))) && (chisqr < chisqrl) && (chisqr < chisqrg)) {
- //if ((projch->GetFunction("fLandauGaus")->GetParameter(1) > 0) && (chisqr < chisqrl) && (chisqr < chisqrg)) {
fNumberFitSuccess++;
CalculChargeCoefMean(kTRUE);
fCurrentCoef[0] = projch->GetFunction("fLandauGaus")->GetParameter(1);
}
//_____________________________________________________________________________
-void AliTRDCalibraFit::FitBisCH(TH1* projch, Double_t mean)
+void AliTRDCalibraFit::FitBisCH(TH1* projch, Double_t mean, Double_t nentries)
{
//
// Fit methode for the gain factor more time consuming
//
+ //Calcul Range of the fit
+ Double_t lastvalue = 0.0;
+ Float_t sumAll = (Float_t) nentries;
+ Int_t sumCurrent = 0;
+ //printf("There are %d bins\n",nybins);
+ for(Int_t k = 0; k <projch->GetNbinsX(); k++){
+ Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
+ if (fraction>fOutliersFitChargeHigh) {
+ lastvalue = projch->GetBinCenter(k+1);
+ //printf("Break by %f\n",lastvalue);
+ break;
+ }
+ sumCurrent += (Int_t) projch->GetBinContent(k+1);
+ }
+ //
//Some parameters to initialise
Double_t widthLandau, widthGaus, mPV, integral;
Double_t chisquarel = 0.0;
Double_t chisquareg = 0.0;
- projch->Fit("landau","0M+",""
- ,(Double_t) mean/6
- ,projch->GetBinCenter(projch->GetNbinsX()));
+ projch->Fit("landau","WWQ0M+",""
+ ,(Double_t) mean/fBeginFitCharge
+ ,lastvalue);
widthLandau = projch->GetFunction("landau")->GetParameter(2);
chisquarel = projch->GetFunction("landau")->GetChisquare();
- projch->Fit("gaus","0M+",""
- ,(Double_t) mean/6
- ,projch->GetBinCenter(projch->GetNbinsX()));
+ projch->Fit("gaus","WWQ0M+",""
+ ,(Double_t) mean/fBeginFitCharge
+ ,lastvalue);
widthGaus = projch->GetFunction("gaus")->GetParameter(2);
chisquareg = projch->GetFunction("gaus")->GetChisquare();
// Setting fit range and start values
Double_t fr[2];
- //Double_t sv[4] = { l3P2, fChargeCoef[1], projch->Integral("width"), fG3P2 };
- //Double_t sv[4] = { fL3P2, fChargeCoef[1], fL3P0, fG3P2 };
Double_t sv[4] = { widthLandau, mPV, integral, widthGaus};
Double_t pllo[4] = { 0.001, 0.001, projch->Integral()/3, 0.001};
Double_t plhi[4] = { 300.0, 300.0, 30*projch->Integral(), 300.0};
Double_t fp[4] = { 1.0, 1.0, 1.0, 1.0 };
Double_t fpe[4] = { 1.0, 1.0, 1.0, 1.0 };
- fr[0] = 0.3 * mean;
- fr[1] = 3.0 * mean;
+ fr[0] = mean/fBeginFitCharge;
+ fr[1] = lastvalue;
fCurrentCoef[0] = 0.0;
fCurrentCoefE = 0.0;
,&fp[0],&fpe[0]
,&chisqr,&ndf);
- Double_t projchPeak;
- Double_t projchFWHM;
- LanGauPro(fp,projchPeak,projchFWHM);
+ //Double_t projchPeak;
+ //Double_t projchFWHM;
+ //LanGauPro(fp,projchPeak,projchFWHM);
if ((fp[1] > 0) && ((fpe[1] < (0.05*fp[1])) && (chisqr < chisquarel) && (chisqr < chisquareg))) {
//if ((fp[1] > 0) && ((chisqr < chisquarel) && (chisqr < chisquareg))) {
cpy->cd();
projch->Draw();
fitsnr->Draw("same");
+ TLine *line = new TLine(fp[1],0.0,fp[1],20000.0);
+ line->Draw("same");
+ }
+ else {
+ delete fitsnr;
+ }
+}
+//_____________________________________________________________________________
+void AliTRDCalibraFit::FitBisCHEx(TH1* projch, Double_t mean, Double_t nentries)
+{
+ //
+ // Fit methode for the gain factor more time consuming
+ //
+
+ //Calcul Range of the fit
+ Double_t lastvalue = 0.0;
+ Float_t sumAll = (Float_t) nentries;
+ Int_t sumCurrent = 0;
+ //printf("There are %d bins\n",nybins);
+ for(Int_t k = 0; k <projch->GetNbinsX(); k++){
+ Double_t fraction = Float_t(sumCurrent)/Float_t(sumAll);
+ if (fraction>fOutliersFitChargeHigh) {
+ lastvalue = projch->GetBinCenter(k+1);
+ //printf("Break by %f\n",lastvalue);
+ break;
+ }
+ sumCurrent += (Int_t) projch->GetBinContent(k+1);
+ }
+ //
+
+
+ //Some parameters to initialise
+ Double_t widthLandau, widthGaus, mPV, integral;
+ Double_t chisquarel = 0.0;
+ Double_t chisquareg = 0.0;
+ projch->Fit("landau","WWQM+",""
+ ,(Double_t) mean/fBeginFitCharge
+ ,lastvalue);
+ widthLandau = projch->GetFunction("landau")->GetParameter(2);
+ chisquarel = projch->GetFunction("landau")->GetChisquare();
+ projch->Fit("gaus","WWQM+",""
+ ,(Double_t) mean/fBeginFitCharge
+ ,lastvalue);
+ widthGaus = projch->GetFunction("gaus")->GetParameter(2);
+ chisquareg = projch->GetFunction("gaus")->GetChisquare();
+
+ mPV = (projch->GetFunction("landau")->GetParameter(1))/2;
+ integral = (projch->GetFunction("gaus")->Integral(0.3*mean,3*mean)+projch->GetFunction("landau")->Integral(0.3*mean,3*mean))/2;
+
+ // Setting fit range and start values
+ Double_t fr[2];
+ //Double_t sv[4] = { l3P2, fChargeCoef[1], projch->Integral("width"), fG3P2 };
+ //Double_t sv[4] = { fL3P2, fChargeCoef[1], fL3P0, fG3P2 };
+ Double_t sv[5] = { widthLandau, mPV, integral, widthGaus, 0.0};
+ Double_t pllo[5] = { 0.001, 0.001, projch->Integral()/3, 0.001, 0.0};
+ Double_t plhi[5] = { 300.0, 300.0, 30*projch->Integral(), 300.0, 2.0};
+ Double_t fp[5] = { 1.0, 1.0, 1.0, 1.0, 1.0};
+ Double_t fpe[5] = { 1.0, 1.0, 1.0, 1.0, 1.0};
+ //
+ //fr[0] = 0.3 * mean;
+ //fr[1] = 3.0 * mean;
+ //
+ fr[0] = mean/fBeginFitCharge;
+ fr[1] = lastvalue;
+
+ fCurrentCoef[0] = 0.0;
+ fCurrentCoefE = 0.0;
+
+ Double_t chisqr = 100.0;
+ Int_t ndf = 1;
+
+ TF1 *fitsnr = 0x0;
+
+ if((mPV > 0.0) && (projch->GetFunction("gaus")->GetParameter(1) > 0.0)) {
+ fitsnr = LanGauFitEx(projch,&fr[0],&sv[0]
+ ,&pllo[0],&plhi[0]
+ ,&fp[0],&fpe[0]
+ ,&chisqr,&ndf);
+ }
+
+ //Double_t projchPeak;
+ //Double_t projchFWHM;
+ //LanGauProEx(fp,projchPeak,projchFWHM);
+
+ if ((fp[1] > 0) && ((fpe[1] < (0.05*fp[1])) && (chisqr < chisquarel) && (chisqr < chisquareg))) {
+ //if ((fp[1] > 0) && ((chisqr < chisquarel) && (chisqr < chisquareg))) {
+ fNumberFitSuccess++;
+ CalculChargeCoefMean(kTRUE);
+ fCurrentCoef[0] = fp[1];
+ fCurrentCoefE = fpe[1];
+ //chargeCoefE2 = chisqr;
+ }
+ else {
+ CalculChargeCoefMean(kFALSE);
+ fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
+ }
+ if (fDebugLevel == 1) {
+ AliInfo(Form("fChargeCoef[0]: %f",(Float_t) fCurrentCoef[0]));
+ TCanvas *cpy = new TCanvas("cpy","cpy",50,50,600,800);
+ cpy->cd();
+ projch->Draw();
+ if(fitsnr) fitsnr->Draw("same");
+ TLine *line = new TLine(fp[1],0.0,fp[1],20000.0);
+ line->Draw("same");
}
else {
delete fitsnr;
// you have to choose fRebin, divider of fNumberBinCharge
//
- TAxis *xhist = hist->GetXaxis();
+ const TAxis *xhist = hist->GetXaxis();
TH1I *rehist = new TH1I("projrebin","",(Int_t) xhist->GetNbins()/fRebin
,xhist->GetBinLowEdge(1)
,xhist->GetBinUpEdge(xhist->GetNbins()));
// you have to choose fRebin divider of fNumberBinCharge
//
- TAxis *xhist = hist->GetXaxis();
+ const TAxis *xhist = hist->GetXaxis();
TH1F *rehist = new TH1F("projrebin","",(Int_t) xhist->GetNbins()/fRebin
,xhist->GetBinLowEdge(1)
,xhist->GetBinUpEdge(xhist->GetNbins()));
Double_t sc = 5.0; // Convolution extends to +-sc Gaussian sigmas
// Variables
- Double_t xx;
- Double_t mpc;
- Double_t fland;
+ Double_t xx = 0.0;
+ Double_t mpc = 0.0;
+ Double_t fland = 0.0;
Double_t sum = 0.0;
- Double_t xlow;
- Double_t xupp;
- Double_t step;
- Double_t i;
+ Double_t xlow = 0.0;
+ Double_t xupp = 0.0;
+ Double_t step = 0.0;
+ Double_t i = 0.0;
// MP shift correction
mpc = par[1] - mpshift * par[0];
for (i = 1.0; i <= np/2; i++) {
xx = xlow + (i-.5) * step;
- fland = TMath::Landau(xx,mpc,par[0]) / par[0];
+ if(par[0] > 0.0) fland = TMath::Landau(xx,mpc,par[0]) / par[0];
sum += fland * TMath::Gaus(x[0],xx,par[3]);
xx = xupp - (i-.5) * step;
- fland = TMath::Landau(xx,mpc,par[0]) / par[0];
+ if(par[0] > 0.0) fland = TMath::Landau(xx,mpc,par[0]) / par[0];
sum += fland * TMath::Gaus(x[0],xx,par[3]);
}
- return (par[2] * step * sum * invsq2pi / par[3]);
+ if(par[3] > 0.0) return (par[2] * step * sum * invsq2pi / par[3]);
+ else return 0.0;
+
+}
+//_____________________________________________________________________________
+Double_t AliTRDCalibraFit::LanGauFunEx(const Double_t *x, const Double_t *par)
+{
+ //
+ // Function for the fit
+ //
+ // Fit parameters:
+ // par[0]=Width (scale) parameter of Landau density
+ // par[1]=Most Probable (MP, location) parameter of Landau density
+ // par[2]=Total area (integral -inf to inf, normalization constant)
+ // par[3]=Width (sigma) of convoluted Gaussian function
+ // par[4]=Exponential Slope Parameter
+ //
+ // In the Landau distribution (represented by the CERNLIB approximation),
+ // the maximum is located at x=-0.22278298 with the location parameter=0.
+ // This shift is corrected within this function, so that the actual
+ // maximum is identical to the MP parameter.
+ //
+
+ // Numeric constants
+ Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2)
+ Double_t mpshift = -0.22278298; // Landau maximum location
+
+ // Control constants
+ Double_t np = 100.0; // Number of convolution steps
+ Double_t sc = 5.0; // Convolution extends to +-sc Gaussian sigmas
+
+ // Variables
+ Double_t xx= 0.0;
+ Double_t mpc= 0.0;
+ Double_t fland = 0.0;
+ Double_t sum = 0.0;
+ Double_t xlow= 0.0;
+ Double_t xupp= 0.0;
+ Double_t step= 0.0;
+ Double_t i= 0.0;
+
+ // MP shift correction
+ mpc = par[1] - mpshift * par[0];
+
+ // Range of convolution integral
+ xlow = x[0] - sc * par[3];
+ xupp = x[0] + sc * par[3];
+
+ step = (xupp - xlow) / np;
+
+ // Convolution integral of Landau and Gaussian by sum
+ for (i = 1.0; i <= np/2; i++) {
+
+ xx = xlow + (i-.5) * step;
+ if(par[0] > 0.0) fland = TMath::Landau(xx,mpc,par[0])*TMath::Exp(-par[4]*xx) / par[0];
+ sum += fland * TMath::Gaus(x[0],xx,par[3]);
+
+ xx = xupp - (i-.5) * step;
+ if(par[0] > 0.0) fland = TMath::Landau(xx,mpc,par[0])*TMath::Exp(-par[4]*xx) / par[0];
+ sum += fland * TMath::Gaus(x[0],xx,par[3]);
+
+ }
+
+ if(par[3] > 0.0) return (par[2] * step * sum * invsq2pi / par[3]);
+ else return 0.0;
}
//_____________________________________________________________________________
ffit->SetParLimits(i,parlimitslo[i],parlimitshi[i]);
}
- his->Fit(funname,"RB0"); // Fit within specified range, use ParLimits, do not plot
+ his->Fit(funname,"WWQRB0"); // Fit within specified range, use ParLimits, do not plot
ffit->GetParameters(fitparams); // Obtain fit parameters
for (i = 0; i < 4; i++) {
return (ffit); // Return fit function
}
-
//_____________________________________________________________________________
-Int_t AliTRDCalibraFit::LanGauPro(const Double_t *params, Double_t &maxx, Double_t &fwhm)
+TF1 *AliTRDCalibraFit::LanGauFitEx(TH1 *his, const Double_t *fitrange, const Double_t *startvalues
+ , const Double_t *parlimitslo, const Double_t *parlimitshi
+ , Double_t *fitparams, Double_t *fiterrors
+ , Double_t *chiSqr, Int_t *ndf) const
{
//
// Function for the fit
//
-
- Double_t p;
- Double_t x;
- Double_t fy;
- Double_t fxr;
- Double_t fxl;
- Double_t step;
- Double_t l;
- Double_t lold;
-
- Int_t i = 0;
- Int_t maxcalls = 10000;
- // Search for maximum
- p = params[1] - 0.1 * params[0];
- step = 0.05 * params[0];
- lold = -2.0;
- l = -1.0;
-
- while ((l != lold) && (i < maxcalls)) {
- i++;
- lold = l;
- x = p + step;
- l = LanGauFun(&x,params);
- if (l < lold) {
- step = -step / 10.0;
- }
- p += step;
- }
+ Int_t i;
+ Char_t funname[100];
- if (i == maxcalls) {
- return (-1);
- }
- maxx = x;
- fy = l / 2.0;
+ TF1 *ffitold = (TF1 *) gROOT->GetListOfFunctions()->FindObject(funname);
+ if (ffitold) {
+ delete ffitold;
+ }
- // Search for right x location of fy
- p = maxx + params[0];
- step = params[0];
- lold = -2.0;
- l = -1e300;
- i = 0;
+ TF1 *ffit = new TF1(funname,LanGauFunEx,fitrange[0],fitrange[1],5);
+ ffit->SetParameters(startvalues);
+ ffit->SetParNames("Width","MP","Area","GSigma","Ex");
- while ( (l != lold) && (i < maxcalls) ) {
- i++;
-
- lold = l;
- x = p + step;
- l = TMath::Abs(LanGauFun(&x,params) - fy);
-
- if (l > lold)
- step = -step/10;
-
- p += step;
+ for (i = 0; i < 5; i++) {
+ ffit->SetParLimits(i,parlimitslo[i],parlimitshi[i]);
}
- if (i == maxcalls)
- return (-2);
-
- fxr = x;
-
-
- // Search for left x location of fy
-
- p = maxx - 0.5 * params[0];
- step = -params[0];
- lold = -2.0;
- l = -1.0e300;
- i = 0;
-
- while ((l != lold) && (i < maxcalls)) {
- i++;
- lold = l;
- x = p + step;
- l = TMath::Abs(LanGauFun(&x,params) - fy);
- if (l > lold) {
- step = -step / 10.0;
- }
- p += step;
- }
+ his->Fit(funname,"WWQRB0"); // Fit within specified range, use ParLimits, do not plot
- if (i == maxcalls) {
- return (-3);
+ ffit->GetParameters(fitparams); // Obtain fit parameters
+ for (i = 0; i < 5; i++) {
+ fiterrors[i] = ffit->GetParError(i); // Obtain fit parameter errors
}
+ chiSqr[0] = ffit->GetChisquare(); // Obtain chi^2
+ ndf[0] = ffit->GetNDF(); // Obtain ndf
- fxl = x;
- fwhm = fxr - fxl;
-
- return (0);
+ return (ffit); // Return fit function
+
}
-//_____________________________________________________________________________
-Double_t AliTRDCalibraFit::GausConstant(const Double_t *x, const Double_t *par)
-{
- //
- // Gaus with identical mean
- //
-
- Double_t gauss = par[0] * TMath::Gaus(x[0],0.0,par[1])+par[2];
-
- return gauss;
-}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff