#include <TVectorD.h>
#include <TROOT.h>
#include <TString.h>
+#include <TLine.h>
#include "AliLog.h"
#include "AliMathBase.h"
#include "AliTRDCalibraMode.h"
#include "AliTRDCalibraVector.h"
#include "AliTRDCalibraVdriftLinearFit.h"
+#include "AliTRDCalibraExbAltFit.h"
#include "AliTRDcalibDB.h"
#include "AliTRDgeometry.h"
#include "AliTRDpadPlane.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
{
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
{
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;
fDebugStreamer = 0x0;
return kTRUE;
+}
+//______________________________________________________________________________________
+Bool_t AliTRDCalibraFit::AnalyseExbAltFit(AliTRDCalibraExbAltFit *calivdli)
+{
+ //
+ // The linear method
+ //
+
+ fStatisticMean = 0.0;
+ fNumberFit = 0;
+ fNumberFitSuccess = 0;
+ fNumberEnt = 0;
+ if(!InitFitExbAlt()) return kFALSE;
+
+
+ for(Int_t idet = 0; idet < 540; idet++){
+
+
+ //printf("detector number %d\n",idet);
+
+ // Take the result
+ TVectorD param(3);
+ TVectorD error(3);
+ fEntriesCurrent = 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];
+ fNumberEnt++;
+ }
+ //printf("Number of entries %d\n",fEntriesCurrent);
+ // Nothing found or not enough statistic
+ if((!heree) || (!here) || (fEntriesCurrent <= fMinEntries)) {
+ NotEnoughStatisticExbAlt();
+ continue;
+ }
+ //param.Print();
+ //error.Print();
+ //Statistics
+ fNumberFit++;
+ fStatisticMean += fEntriesCurrent;
+
+ // Statistics
+ fNumberFitSuccess ++;
+
+ // Put the fCurrentCoef
+ if(TMath::Abs(param[2])>0.0001){
+ fCurrentCoef2[0] = -param[1]/2/param[2];
+ fCurrentCoefE2 = 0;//error[1];
+ }else{
+ fCurrentCoef2[0] = 100;
+ fCurrentCoefE2 = 0;//error[1];
+ }
+
+ // Fill
+ FillInfosFitExbAlt();
+
+ }
+ // 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));
+ }
+ else {
+ AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt));
+ }
+ delete fDebugStreamer;
+ fDebugStreamer = 0x0;
+ return kTRUE;
+
}
//____________Functions fit Online CH2d________________________________________
void AliTRDCalibraFit::AnalyseLinearFittersAllTogether(AliTRDCalibraVdriftLinearFit *calivdli, Double_t &vdriftoverall, Double_t &exboverall)
return object;
+}
+//_____________________________________________________________________________
+AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectExbAlt(const TObjArray *vectorFit)
+{
+ //
+ // It creates the AliTRDCalDet object from the AliTRDFitInfo2
+ // It takes the min value of the coefficients per detector
+ // This object has to be written in the database
+ //
+
+ // Create the DetObject
+ AliTRDCalDet *object = new AliTRDCalDet("tan(lorentzangle)","tan(lorentzangle) (detector value)");
+
+
+ Int_t loop = (Int_t) vectorFit->GetEntriesFast();
+ if(loop != 540) AliInfo("The Vector Fit is not complete!");
+ Int_t detector = -1;
+ Float_t value = 0.0;
+
+ for (Int_t k = 0; k < loop; k++) {
+ detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
+ /*
+ Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
+ Int_t colMax = fGeo->GetColMax(GetLayer(detector));
+ Float_t min = 100.0;
+ for (Int_t row = 0; row < rowMax; row++) {
+ for (Int_t col = 0; col < colMax; col++) {
+ value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
+ mean += -TMath::Abs(value);
+ count++;
+ } // Col
+ } // Row
+ if(count > 0) mean = mean/count;
+ */
+ value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
+ //if(value > 70.0) value = value-100.0;
+ object->SetValue(detector,value);
+ }
+
+ return object;
+
}
//_____________________________________________________________________________
TObject *AliTRDCalibraFit::CreatePadObjectGain(const TObjArray *vectorFit, Double_t scaleFitFactor, const AliTRDCalDet *detobject)
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++) {
return kTRUE;
}
//____________Functions for initialising the AliTRDCalibraFit in the code_________
+Bool_t AliTRDCalibraFit::InitFitExbAlt()
+{
+ //
+ // Init the fCalDet, fVectorFit fCurrentCoefDetector
+ //
+
+ gDirectory = gROOT;
+
+ if( fCurrentCoefDetector2 ) delete [] fCurrentCoefDetector2;
+ fCurrentCoefDetector2 = new Float_t[2304];
+ for (Int_t k = 0; k < 2304; k++) {
+ fCurrentCoefDetector2[k] = 0.0;
+ }
+
+ return kTRUE;
+}
+//____________Functions for initialising the AliTRDCalibraFit in the code_________
void AliTRDCalibraFit::InitfCountDetAndfCount(Int_t i)
{
//
return kTRUE;
}
+//____________Functions for initialising the AliTRDCalibraFit in the code_________
+Bool_t AliTRDCalibraFit::NotEnoughStatisticExbAlt()
+{
+ //
+ // For the case where there are not enough entries in the histograms
+ // of the calibration group, the value present in the choosen database
+ // will be put. A negativ sign enables to know that a fit was not possible.
+ //
+
+ Int_t factor = 0;
+ if(GetStack(fCountDet) == 2) factor = 1728;
+ else factor = 2304;
+
+
+ // Fill the fCurrentCoefDetector
+ for (Int_t k = 0; k < factor; k++) {
+ fCurrentCoefDetector2[k] = 100.0;
+ }
+
+ fCurrentCoef2[0] = 100.0;
+ fCurrentCoefE2 = 0.0;
+
+ FillFillExbAlt();
+
+ return kTRUE;
+}
+
//____________Functions for initialising the AliTRDCalibraFit in the code_________
Bool_t AliTRDCalibraFit::FillInfosFitCH(Int_t idect)
{
return kTRUE;
+}
+//____________Functions for initialising the AliTRDCalibraFit in the code_________
+Bool_t AliTRDCalibraFit::FillInfosFitExbAlt()
+{
+ //
+ // Fill the coefficients found with the fits or other
+ // methods from the Fit functions
+ //
+
+ Int_t factor = 0;
+ if(GetStack(fCountDet) == 2) factor = 1728;
+ else factor = 2304;
+
+ // Pointer to the branch
+ for (Int_t k = 0; k < factor; k++) {
+ fCurrentCoefDetector2[k] = fCurrentCoef2[0];
+ }
+
+ FillFillExbAlt();
+
+ return kTRUE;
+
}
//________________________________________________________________________________
void AliTRDCalibraFit::FillFillCH(Int_t idect)
//
// 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++) {
"\n";
}
+}
+//________________________________________________________________________________
+void AliTRDCalibraFit::FillFillExbAlt()
+{
+ //
+ // DebugStream and fVectorFit
+ //
+
+ // End of one detector
+ FillVectorFit2();
+
+
+ // Reset
+ for (Int_t k = 0; k < 2304; k++) {
+ fCurrentCoefDetector2[k] = 0.0;
+ }
+
+
+ if(fDebugLevel > 1){
+
+ if ( !fDebugStreamer ) {
+ //debug stream
+ TDirectory *backup = gDirectory;
+ fDebugStreamer = new TTreeSRedirector("TRDDebugFitExbAlt.root");
+ if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
+ }
+
+ //Debug: comparaison of the different methods (okey for first time but not for iterative procedure)
+ AliTRDpadPlane *padplane = fGeo->GetPadPlane(GetLayer(fCountDet),GetStack(fCountDet));
+ Float_t rowmd = (padplane->GetRow0()+padplane->GetRowEnd())/2.;
+ Float_t r = AliTRDgeometry::GetTime0(GetLayer(fCountDet));
+ Float_t tiltangle = padplane->GetTiltingAngle();
+ Int_t detector = fCountDet;
+ Int_t stack = GetStack(fCountDet);
+ Int_t layer = GetLayer(fCountDet);
+ Float_t vf = fCurrentCoef2[0];
+ Float_t vfE = fCurrentCoefE2;
+
+ (* fDebugStreamer) << "FillFillLinearFitter"<<
+ "detector="<<detector<<
+ "stack="<<stack<<
+ "layer="<<layer<<
+ "rowmd="<<rowmd<<
+ "r="<<r<<
+ "tiltangle="<<tiltangle<<
+ "vf="<<vf<<
+ "vfE="<<vfE<<
+ "\n";
+ }
+
}
//
//____________Calcul Coef Mean_________________________________________________
y[1] = pentea->GetBinContent(binmax-1);
y[2] = pentea->GetBinContent(binmax);
CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
- AliInfo("At the limit for beginning!");
+ //AliInfo("At the limit for beginning!");
break;
case 2:
minnn = pentea->GetBinCenter(binmax-2);
if (binmin <= 1) {
binmin = 2;
put = 1;
- AliInfo("Put the binmax from 1 to 2 to enable the fit");
+ //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!");
+ //AliInfo("Too many fluctuations at the end!");
put = kFALSE;
}
if((projPH->GetBinContent(binmin)+projPH->GetBinError(binmin)) > (projPH->GetBinContent(binmin-1))) {
- AliInfo("Too many fluctuations at the end!");
+ //AliInfo("Too many fluctuations at the end!");
put = kFALSE;
}
if(TMath::Abs(pente->GetBinContent(binmin+1)) <= 0.0000000000001){
- AliInfo("No entries for the next bin!");
+ //AliInfo("No entries for the next bin!");
pente->SetBinContent(binmin,0);
if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
}
(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");
+ //AliInfo("polynome 4 false 2");
put = kFALSE;
}
// poly 3
}
if((binmin == (nbins-1)) && ((binmin-2) < TMath::Min(binmax+4, projPH->GetNbinsX()))) {
put = kFALSE;
- AliInfo("At the limit for the drift and not usable!");
+ //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;
- AliInfo("For the drift...problem!");
+ //AliInfo("For the drift...problem!");
}
//pass but should not happen
if((binmin <= (nbins-3)) && (binmin < TMath::Min(binmax+6, projPH->GetNbinsX()))){
put = kFALSE;
- AliInfo("For the drift...problem!");
+ //AliInfo("For the drift...problem!");
}
if(put) {
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;
+ Double_t chisqrl = 0.0;
+
+ projch->Fit("landau","WWQ+",""
+ ,(Double_t) mean/fBeginFitCharge
+ ,lastvalue);
+ 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;
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;
+ his->Fit(funname,"WWQRB0"); // Fit within specified range, use ParLimits, do not plot
-
- // 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;
- }
-
- 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