#include <TAxis.h>
#include <TMath.h>
#include <TDirectory.h>
-#include <TROOT.h>
#include <TTreeStream.h>
#include <TLinearFitter.h>
#include <TVectorD.h>
-#include <TArrayF.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)
,fAccCDB(kFALSE)
,fMinEntries(800)
,fRebin(1)
+ ,fScaleGain(0.02431)
,fNumberFit(0)
,fNumberFitSuccess(0)
,fNumberEnt(0)
,fCalROC(0x0)
,fCalDet2(0x0)
,fCalROC2(0x0)
+ ,fCalDetVdriftUsed(0x0)
+ ,fCalDetExBUsed(0x0)
,fCurrentCoefDetector(0x0)
,fCurrentCoefDetector2(0x0)
,fVectorFit(0)
,fNumberOfBinsExpected(c.fNumberOfBinsExpected)
,fMethod(c.fMethod)
,fBeginFitCharge(c.fBeginFitCharge)
+,fOutliersFitChargeLow(c.fOutliersFitChargeLow)
+,fOutliersFitChargeHigh(c.fOutliersFitChargeHigh)
,fFitPHPeriode(c.fFitPHPeriode)
,fTakeTheMaxPH(c.fTakeTheMaxPH)
,fT0Shift0(c.fT0Shift0)
,fAccCDB(c.fAccCDB)
,fMinEntries(c.fMinEntries)
,fRebin(c.fRebin)
+,fScaleGain(c.fScaleGain)
,fNumberFit(c.fNumberFit)
,fNumberFitSuccess(c.fNumberFitSuccess)
,fNumberEnt(c.fNumberEnt)
,fCalROC(0x0)
,fCalDet2(0x0)
,fCalROC2(0x0)
+,fCalDetVdriftUsed(0x0)
+,fCalDetExBUsed(0x0)
,fCurrentCoefDetector(0x0)
,fCurrentCoefDetector2(0x0)
,fVectorFit(0)
if(c.fCalROC) fCalROC = new AliTRDCalROC(*c.fCalROC);
if(c.fCalROC2) fCalROC = new AliTRDCalROC(*c.fCalROC2);
+ if(c.fCalDetVdriftUsed) fCalDetVdriftUsed = new AliTRDCalDet(*c.fCalDetVdriftUsed);
+ if(c.fCalDetExBUsed) fCalDetExBUsed = new AliTRDCalDet(*c.fCalDetExBUsed);
+
fVectorFit.SetName(c.fVectorFit.GetName());
for(Int_t k = 0; k < c.fVectorFit.GetEntriesFast(); k++){
AliTRDFitInfo *fitInfo = new AliTRDFitInfo();
if ( fCalDet2 ) delete fCalDet2;
if ( fCalROC ) delete fCalROC;
if ( fCalROC2 ) delete fCalROC2;
+ if ( fCalDetVdriftUsed) delete fCalDetVdriftUsed;
+ if ( fCalDetExBUsed) delete fCalDetExBUsed;
if( fCurrentCoefDetector ) delete [] fCurrentCoefDetector;
if( fCurrentCoefDetector2 ) delete [] fCurrentCoefDetector2;
fVectorFit.Delete();
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)
//
// Set the calibration mode
- const char *name = ch->GetTitle();
+ //const char *name = ch->GetTitle();
+ TString name = ch->GetTitle();
if(!SetModeCalibration(name,0)) return kFALSE;
// Number of Ybins (detectors or groups of pads)
{
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 {
//
// Set the calibraMode
- const char *name = calvect->GetNameCH();
+ //const char *name = calvect->GetNameCH();
+ TString name = calvect->GetNameCH();
if(!SetModeCalibration(name,0)) return kFALSE;
// Number of Xbins (detectors or groups of pads)
{
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 {
fDebugStreamer = 0x0;
return kTRUE;
}
+//____________Functions fit Online CH2d________________________________________
+Double_t AliTRDCalibraFit::AnalyseCHAllTogether(const TH2I *ch)
+{
+ //
+ // Fit the 1D histos, projections of the 2D ch on the Xaxis, for each
+ // calibration group normalized the resulted coefficients (to 1 normally)
+ //
+ Int_t nbins = ch->GetNbinsX();// charge
+ Int_t nybins = ch->GetNbinsY();// groups number
+ // Take the histo
+ TH1I *projch = (TH1I *) ch->ProjectionX("projch",1,nybins+1,(Option_t *)"e");
+ projch->SetDirectory(0);
+ // Number of entries for this calibration group
+ Double_t nentries = 0.0;
+ Double_t mean = 0.0;
+ for (Int_t k = 0; k < nbins; k++) {
+ nentries += projch->GetBinContent(k+1);
+ mean += projch->GetBinCenter(k+1)*projch->GetBinContent(k+1);
+ projch->SetBinError(k+1,TMath::Sqrt(projch->GetBinContent(k+1)));
+ }
+ projch->SetEntries(nentries);
+ //printf("The number of entries for the group %d is %f\n",idect,nentries);
+ if (nentries > 0) {
+ mean /= nentries;
+ }
+ // This detector has not enough statistics or was off
+ if (nentries <= fMinEntries) {
+ delete projch;
+ AliInfo(Form("There are %d entries for all together, it is not enough (%d)",(Int_t) nentries, fMinEntries));
+ return -100.0;
+ }
+ //Method choosen
+ switch(fMethod)
+ {
+ case 0: FitMeanW((TH1 *) projch, nentries); break;
+ case 1: FitMean((TH1 *) projch, nentries, 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;
+
+ if(fCurrentCoef[0] > 0.0) return fCurrentCoef[0];
+ else return -100.0;
+
+}
//________________functions fit Online PH2d____________________________________
-Bool_t AliTRDCalibraFit::AnalysePH(const TProfile2D *ph)
+Double_t AliTRDCalibraFit::AnalysePHAllTogether(const TProfile2D *ph)
{
//
// Take the 1D profiles (average pulse height), projections of the 2D PH
// A first calibration of T0 is also made using the same method
//
- // Set the calibration mode
- const char *name = ph->GetTitle();
- if(!SetModeCalibration(name,1)) return kFALSE;
-
// Number of Xbins (detectors or groups of pads)
Int_t nbins = ph->GetNbinsX();// time
Int_t nybins = ph->GetNbinsY();// calibration group
- if (!InitFit(nybins,1)) {
+
+ // Take the histo
+ TH1D *projph = (TH1D *) ph->ProjectionX("projph",1,nybins+1,(Option_t *) "e");
+ projph->SetDirectory(0);
+ // Number of entries for this calibration group
+ Double_t nentries = 0;
+ for(Int_t idect = 0; idect < nybins; idect++){
+ for (Int_t k = 0; k < nbins; k++) {
+ Int_t binnb = (nbins+2)*(idect+1)+(k+1);
+ nentries += ph->GetBinEntries(binnb);
+ }
+ }
+ //printf("AnalysePHAllTogether:: the number of entries is %f\n",nentries);
+ // This detector has not enough statistics or was off
+ if (nentries <= fMinEntries) {
+ AliInfo(Form("There are %d entries for all together, it is not enough (%d)",(Int_t) nentries, fMinEntries));
+ if (fDebugLevel != 1) {
+ delete projph;
+ }
+ return -100.0;
+ }
+ //Method choosen
+ //printf("Method\n");
+ switch(fMethod)
+ {
+ case 0: FitLagrangePoly((TH1 *) projph); break;
+ case 1: FitPente((TH1 *) projph); break;
+ case 2: FitPH((TH1 *) projph,0); break;
+ default: return -100.0;
+ }
+ // Memory!!!
+ if (fDebugLevel != 1) {
+ delete projph;
+ }
+ delete fDebugStreamer;
+ fDebugStreamer = 0x0;
+
+ if(fCurrentCoef[0] > 0.0) return fCurrentCoef[0];
+ else return -100.0;
+
+}
+//____________Functions fit Online PH2d________________________________________
+Bool_t AliTRDCalibraFit::AnalysePH(AliTRDCalibraVector *calvect)
+{
+ //
+ // Reconstruct the average pulse height from the vectorPH for each
+ // calibration group
+ // Reconstruct a drift velocity
+ // A first calibration of T0 is also made using the same method (slope method)
+ //
+
+ // Set the calibration mode
+ //const char *name = calvect->GetNamePH();
+ TString name = calvect->GetNamePH();
+ if(!SetModeCalibration(name,1)) return kFALSE;
+
+ // Number of Xbins (detectors or groups of pads)
+ if (!InitFit((432*calvect->GetDetCha0(1)+108*calvect->GetDetCha2(1)),1)) {
return kFALSE;
}
if (!InitFitPH()) {
InitfCountDetAndfCount(1);
// Beginning of the loop
for (Int_t idect = fDect1; idect < fDect2; idect++) {
- // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi.......
+ // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi...........
UpdatefCountDetAndfCount(idect,1);
ReconstructFitRowMinRowMax(idect,1);
// Take the histo
- TH1D *projph = (TH1D *) ph->ProjectionX("projph",idect+1,idect+1,(Option_t *) "e");
- projph->SetDirectory(0);
- // Number of entries for this calibration group
- Double_t nentries = 0;
- for (Int_t k = 0; k < nbins; k++) {
- Int_t binnb = (nbins+2)*(idect+1)+(k+1);
- nentries += ph->GetBinEntries(binnb);
+ fEntriesCurrent = 0;
+ if(!calvect->GetPHEntries(fCountDet)) {
+ NotEnoughStatisticPH(idect,fEntriesCurrent);
+ continue;
}
- if (nentries > 0) {
- fNumberEnt++;
- }
- //printf("The number of entries for the group %d is %f\n",idect,nentries);
+ TString tname("PH");
+ tname += idect;
+ TH1F *projph = calvect->CorrectTheError((TGraphErrors *) (calvect->ConvertVectorPHTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) tname)),fEntriesCurrent);
+ projph->SetDirectory(0);
+ if(fEntriesCurrent > 0) fNumberEnt++;
+ //printf("The number of entries for the group %d is %d\n",idect,fEntriesCurrent);
// This detector has not enough statistics or was off
- if (nentries <= fMinEntries) {
- //printf("Not enough statistic!\n");
- NotEnoughStatisticPH(idect,nentries);
- if (fDebugLevel != 1) {
- delete projph;
- }
+ if (fEntriesCurrent <= fMinEntries) {
+ //printf("Not enough stat!\n");
+ NotEnoughStatisticPH(idect,fEntriesCurrent);
continue;
}
- // Statistics of the histos fitted
+ // Statistic of the histos fitted
fNumberFit++;
- fStatisticMean += nentries;
+ fStatisticMean += fEntriesCurrent;
// Calcul of "real" coef
CalculVdriftCoefMean();
CalculT0CoefMean();
default: return kFALSE;
}
// Fill the tree if end of a detector or only the pointer to the branch!!!
- FillInfosFitPH(idect,nentries);
- // Memory!!!
- if (fDebugLevel != 1) {
- delete projph;
- }
+ FillInfosFitPH(idect,fEntriesCurrent);
} // 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 {
fDebugStreamer = 0x0;
return kTRUE;
}
-//____________Functions fit Online PH2d________________________________________
-Bool_t AliTRDCalibraFit::AnalysePH(AliTRDCalibraVector *calvect)
+//________________functions fit Online PH2d____________________________________
+Bool_t AliTRDCalibraFit::AnalysePH(const TProfile2D *ph)
{
//
- // Reconstruct the average pulse height from the vectorPH for each
- // calibration group
+ // Take the 1D profiles (average pulse height), projections of the 2D PH
+ // on the Xaxis, for each calibration group
// Reconstruct a drift velocity
- // A first calibration of T0 is also made using the same method (slope method)
+ // A first calibration of T0 is also made using the same method
//
// Set the calibration mode
- const char *name = calvect->GetNamePH();
+ //const char *name = ph->GetTitle();
+ TString name = ph->GetTitle();
if(!SetModeCalibration(name,1)) return kFALSE;
+
+ //printf("Mode calibration set\n");
// Number of Xbins (detectors or groups of pads)
- if (!InitFit((432*calvect->GetDetCha0(1)+108*calvect->GetDetCha2(1)),1)) {
+ Int_t nbins = ph->GetNbinsX();// time
+ Int_t nybins = ph->GetNbinsY();// calibration group
+ if (!InitFit(nybins,1)) {
return kFALSE;
}
+
+ //printf("Init fit\n");
+
if (!InitFitPH()) {
return kFALSE;
}
+
+ //printf("Init fit PH\n");
+
fStatisticMean = 0.0;
fNumberFit = 0;
fNumberFitSuccess = 0;
fNumberEnt = 0;
// Init fCountDet and fCount
InitfCountDetAndfCount(1);
+ //printf("Init Count Det and fCount %d, %d\n",fDect1,fDect2);
+
// Beginning of the loop
for (Int_t idect = fDect1; idect < fDect2; idect++) {
- // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi...........
+ //printf("idect = %d\n",idect);
+ // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi.......
UpdatefCountDetAndfCount(idect,1);
ReconstructFitRowMinRowMax(idect,1);
// Take the histo
- fEntriesCurrent = 0;
- if(!calvect->GetPHEntries(fCountDet)) {
- NotEnoughStatisticPH(idect,fEntriesCurrent);
- continue;
+ TH1D *projph = (TH1D *) ph->ProjectionX("projph",idect+1,idect+1,(Option_t *) "e");
+ projph->SetDirectory(0);
+ // Number of entries for this calibration group
+ Double_t nentries = 0;
+ for (Int_t k = 0; k < nbins; k++) {
+ Int_t binnb = (nbins+2)*(idect+1)+(k+1);
+ nentries += ph->GetBinEntries(binnb);
}
- TString tname("PH");
- tname += idect;
- TH1F *projph = calvect->CorrectTheError((TGraphErrors *) (calvect->ConvertVectorPHTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) tname)),fEntriesCurrent);
- projph->SetDirectory(0);
- if(fEntriesCurrent > 0) fNumberEnt++;
- //printf("The number of entries for the group %d is %d\n",idect,fEntriesCurrent);
+ if (nentries > 0) {
+ fNumberEnt++;
+ }
+ //printf("The number of entries for the group %d is %f\n",idect,nentries);
// This detector has not enough statistics or was off
- if (fEntriesCurrent <= fMinEntries) {
- //printf("Not enough stat!\n");
- NotEnoughStatisticPH(idect,fEntriesCurrent);
+ if (nentries <= fMinEntries) {
+ //printf("Not enough statistic!\n");
+ NotEnoughStatisticPH(idect,nentries);
+ if (fDebugLevel != 1) {
+ delete projph;
+ }
continue;
}
- // Statistic of the histos fitted
+ // Statistics of the histos fitted
fNumberFit++;
- fStatisticMean += fEntriesCurrent;
+ fStatisticMean += nentries;
// Calcul of "real" coef
CalculVdriftCoefMean();
CalculT0CoefMean();
//Method choosen
+ //printf("Method\n");
switch(fMethod)
{
case 0: FitLagrangePoly((TH1 *) projph); break;
default: return kFALSE;
}
// Fill the tree if end of a detector or only the pointer to the branch!!!
- FillInfosFitPH(idect,fEntriesCurrent);
+ FillInfosFitPH(idect,nentries);
+ // Memory!!!
+ if (fDebugLevel != 1) {
+ delete projph;
+ }
} // 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 {
//
// Set the calibration mode
- const char *name = prf->GetTitle();
+ //const char *name = prf->GetTitle();
+ TString name = prf->GetTitle();
if(!SetModeCalibration(name,2)) return kFALSE;
// Number of Ybins (detectors or groups of pads)
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 {
//
// Set the calibration mode
- const char *name = prf->GetTitle();
+ //const char *name = prf->GetTitle();
+ TString name = prf->GetTitle();
if(!SetModeCalibration(name,2)) return kFALSE;
// Number of Ybins (detectors or groups of pads)
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 {
//
// Set the calibra mode
- const char *name = calvect->GetNamePRF();
+ //const char *name = calvect->GetNamePRF();
+ TString name = calvect->GetNamePRF();
if(!SetModeCalibration(name,2)) return kFALSE;
//printf("test0 %s\n",name);
} // 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));
//
// Set the calibra mode
- const char *name = calvect->GetNamePRF();
+ //const char *name = calvect->GetNamePRF();
+ TString name = calvect->GetNamePRF();
if(!SetModeCalibration(name,2)) return kFALSE;
//printf("test0 %s\n",name);
Int_t nbg = GetNumberOfGroupsPRF((const char *)name);
} // 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;
}
// CalculDatabaseVdriftandTan
CalculVdriftLorentzCoef();
+ //printf("AliTRDCalibraFit::AnalyzeVdriftLinearFit detector %d, vdrift %f and %f and exB %f and %f\n",idet,fCalDetVdriftUsed->GetValue(idet),fCurrentCoef[1],fCalDetExBUsed->GetValue(idet),fCurrentCoef2[1]);
// Statistics
fNumberFitSuccess ++;
fCurrentCoef2[0] = (param[0]+fCurrentCoef[1]*fCurrentCoef2[1])/fCurrentCoef[0];
fCurrentCoefE = error[1];
fCurrentCoefE2 = error[0];
- if((fCurrentCoef2[0] != 0.0) && (param[0] != 0.0)){
+ if((TMath::Abs(fCurrentCoef2[0]) > 0.0000001) && (TMath::Abs(param[0]) > 0.0000001)){
fCurrentCoefE2 = (fCurrentCoefE2/param[0]+fCurrentCoefE/fCurrentCoef[0])*fCurrentCoef2[0];
}
}
// 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));
+ }
+ 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);
+ 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) {
+ entriesCurrent = error[2];
+ fNumberEnt++;
+ }
+ //printf("Number of entries %d\n",fEntriesCurrent);
+ // Nothing found or not enough statistic
+ if((!heree) || (!here) || (entriesCurrent <= fMinEntries)) {
+ NotEnoughStatisticExbAlt();
+ continue;
+ }
+ //param.Print();
+ //error.Print();
+ //Statistics
+ fNumberFit++;
+ fStatisticMean += entriesCurrent;
+
+ // 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: %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));
fDebugStreamer = 0x0;
return kTRUE;
+}
+//____________Functions fit Online CH2d________________________________________
+void AliTRDCalibraFit::AnalyseLinearFittersAllTogether(AliTRDCalibraVdriftLinearFit *calivdli, Double_t &vdriftoverall, Double_t &exboverall)
+{
+ //
+ // The linear method
+ //
+
+ // Get the mean vdrift and exb used
+ Double_t meanvdriftused = 0.0;
+ Double_t meanexbused = 0.0;
+ Double_t counterdet = 0.0;
+ if((!fCalDetVdriftUsed) || (!fCalDetExBUsed)) {
+ vdriftoverall = -100.0;
+ exboverall = 100.0;
+ return;
+ }
+
+ // Add histos
+
+ TH2S *linearfitterhisto = 0x0;
+
+ for(Int_t idet = 0; idet < 540; idet++){
+
+ TH2S * u = calivdli->GetLinearFitterHistoForce(idet);
+ Double_t detectorentries = u->Integral();
+ meanvdriftused += fCalDetVdriftUsed->GetValue(idet)*detectorentries;
+ meanexbused += fCalDetExBUsed->GetValue(idet)*detectorentries;
+ counterdet += detectorentries;
+
+ //printf("detectorentries %f\n",detectorentries);
+
+ //printf("AliTRDCalibraFit::AnalyzeVdriftLinearFitsAllTogether detector %d, vdrift %f and exB %f\n",idet,fCalDetVdriftUsed->GetValue(idet),fCalDetExBUsed->GetValue(idet));
+
+ if(idet == 0) linearfitterhisto = u;
+ else linearfitterhisto->Add(u);
+
+ }
+ if(counterdet > 0.0){
+ meanvdriftused = meanvdriftused/counterdet;
+ meanexbused = meanexbused/counterdet;
+ }
+ else {
+ vdriftoverall = -100.0;
+ exboverall = 100.0;
+ return;
+ }
+
+
+ //printf("AliTRDCalibraFit::AnalyzeVdriftLinearFitsAllTogether MEAN vdrift %f and exB %f\n",meanvdriftused,meanexbused);
+
+ // Fit
+
+ Double_t entries = 0;
+ TAxis *xaxis = linearfitterhisto->GetXaxis();
+ TAxis *yaxis = linearfitterhisto->GetYaxis();
+ TLinearFitter linearfitter = TLinearFitter(2,"pol1");
+ //printf("test\n");
+ Double_t integral = linearfitterhisto->Integral();
+ //printf("Integral is %f\n",integral);
+ Bool_t securitybreaking = kFALSE;
+ if(TMath::Abs(integral-1199) < 0.00001) securitybreaking = kTRUE;
+ for(Int_t ibinx = 0; ibinx < linearfitterhisto->GetNbinsX(); ibinx++){
+ for(Int_t ibiny = 0; ibiny < linearfitterhisto->GetNbinsY(); ibiny++){
+ if(linearfitterhisto->GetBinContent(ibinx+1,ibiny+1)>0){
+ Double_t x = xaxis->GetBinCenter(ibinx+1);
+ Double_t y = yaxis->GetBinCenter(ibiny+1);
+
+ for(Int_t k = 0; k < (Int_t)linearfitterhisto->GetBinContent(ibinx+1,ibiny+1); k++){
+ if(!securitybreaking){
+ linearfitter.AddPoint(&x,y);
+ entries = entries+1.;
+ }
+ else {
+ if(entries< 1198.0){
+ linearfitter.AddPoint(&x,y);
+ entries = entries + 1.;
+ }
+ }
+ }
+
+ }
+ }
+ }
+
+ //printf("AnalyseLinearFittersAllTogether::Find %d entries\n",entries);
+ //printf("Minstats %d\n",fMinEntries);
+
+
+
+ // Eval the linear fitter
+ if(entries > fMinEntries){
+ TVectorD par = TVectorD(2);
+ //printf("Fit\n");
+ if((linearfitter.EvalRobust(0.8)==0)) {
+ //printf("Take the param\n");
+ linearfitter.GetParameters(par);
+ //printf("Done\n");
+ //par.Print();
+ //printf("Finish\n");
+ // 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.00001) fCurrentCoef2[0] = (par[0]+meanvdriftused*meanexbused)/fCurrentCoef[0];
+ else fCurrentCoef2[0] = 100.0;
+
+ }
+ else {
+
+ fCurrentCoef[0] = -100.0;
+ fCurrentCoef2[0] = 100.0;
+
+ }
+
+
+ }
+ else {
+
+ fCurrentCoef[0] = -100.0;
+ fCurrentCoef2[0] = 100.0;
+
+ }
+
+ vdriftoverall = fCurrentCoef[0];
+ exboverall = fCurrentCoef2[0];
+
+
+ delete linearfitterhisto;
+ delete fDebugStreamer;
+ fDebugStreamer = 0x0;
+
}
//____________Functions for seeing if the pad is really okey___________________
//_____________________________________________________________________________
-Int_t AliTRDCalibraFit::GetNumberOfGroupsPRF(const char* nametitle)
+Int_t AliTRDCalibraFit::GetNumberOfGroupsPRF(TString nametitle)
{
//
// Get numberofgroupsprf
const Char_t *pattern6 = "Ngp6";
// Nrphi mode
- if (strstr(nametitle,pattern0)) {
+ if (strstr(nametitle.Data(),pattern0)) {
return 0;
}
- if (strstr(nametitle,pattern1)) {
+ if (strstr(nametitle.Data(),pattern1)) {
return 1;
}
- if (strstr(nametitle,pattern2)) {
+ if (strstr(nametitle.Data(),pattern2)) {
return 2;
}
- if (strstr(nametitle,pattern3)) {
+ if (strstr(nametitle.Data(),pattern3)) {
return 3;
}
- if (strstr(nametitle,pattern4)) {
+ if (strstr(nametitle.Data(),pattern4)) {
return 4;
}
- if (strstr(nametitle,pattern5)) {
+ if (strstr(nametitle.Data(),pattern5)) {
return 5;
}
- if (strstr(nametitle,pattern6)){
+ if (strstr(nametitle.Data(),pattern6)){
return 6;
}
else return -1;
}
//_____________________________________________________________________________
-Bool_t AliTRDCalibraFit::SetModeCalibration(const char *name, Int_t i)
+Bool_t AliTRDCalibraFit::SetModeCalibration(TString name, Int_t i)
{
//
// Set fNz[i] and fNrphi[i] of the AliTRDCalibraFit::Instance()
}
//_____________________________________________________________________________
-Bool_t AliTRDCalibraFit::SetNrphiFromTObject(const char *name, Int_t i)
+Bool_t AliTRDCalibraFit::SetNrphiFromTObject(TString name, Int_t i)
{
//
// Set fNrphi[i] of the AliTRDCalibraFit::Instance()
const Char_t *patternz100 = "Nz100";
// Nrphi mode
- if ((strstr(name,patternrphi100)) && (strstr(name,patternz100))) {
+ if ((strstr(name.Data(),patternrphi100)) && (strstr(name.Data(),patternz100))) {
fCalibraMode->SetAllTogether(i);
fNbDet = 540;
if (fDebugLevel > 1) {
}
return kTRUE;
}
- if ((strstr(name,patternrphi10)) && (strstr(name,patternz10))) {
+ if ((strstr(name.Data(),patternrphi10)) && (strstr(name.Data(),patternz10))) {
fCalibraMode->SetPerSuperModule(i);
fNbDet = 30;
if (fDebugLevel > 1) {
return kTRUE;
}
- if (strstr(name,patternrphi0)) {
+ if (strstr(name.Data(),patternrphi0)) {
fCalibraMode->SetNrphi(i ,0);
if (fDebugLevel > 1) {
AliInfo(Form("fNbDet %d and 0",fNbDet));
}
return kTRUE;
}
- if (strstr(name,patternrphi1)) {
+ if (strstr(name.Data(),patternrphi1)) {
fCalibraMode->SetNrphi(i, 1);
if (fDebugLevel > 1) {
AliInfo(Form("fNbDet %d and 1",fNbDet));
}
return kTRUE;
}
- if (strstr(name,patternrphi2)) {
+ if (strstr(name.Data(),patternrphi2)) {
fCalibraMode->SetNrphi(i, 2);
if (fDebugLevel > 1) {
AliInfo(Form("fNbDet %d and 2",fNbDet));
}
return kTRUE;
}
- if (strstr(name,patternrphi3)) {
+ if (strstr(name.Data(),patternrphi3)) {
fCalibraMode->SetNrphi(i, 3);
if (fDebugLevel > 1) {
AliInfo(Form("fNbDet %d and 3",fNbDet));
}
return kTRUE;
}
- if (strstr(name,patternrphi4)) {
+ if (strstr(name.Data(),patternrphi4)) {
fCalibraMode->SetNrphi(i, 4);
if (fDebugLevel > 1) {
AliInfo(Form("fNbDet %d and 4",fNbDet));
}
return kTRUE;
}
- if (strstr(name,patternrphi5)) {
+ if (strstr(name.Data(),patternrphi5)) {
fCalibraMode->SetNrphi(i, 5);
if (fDebugLevel > 1) {
AliInfo(Form("fNbDet %d and 5",fNbDet));
}
return kTRUE;
}
- if (strstr(name,patternrphi6)) {
+ if (strstr(name.Data(),patternrphi6)) {
fCalibraMode->SetNrphi(i, 6);
if (fDebugLevel > 1) {
AliInfo(Form("fNbDet %d and 6",fNbDet));
}
//_____________________________________________________________________________
-Bool_t AliTRDCalibraFit::SetNzFromTObject(const char *name, Int_t i)
+Bool_t AliTRDCalibraFit::SetNzFromTObject(TString name, Int_t i)
{
//
// Set fNz[i] of the AliTRDCalibraFit::Instance()
const Char_t *patternz10 = "Nz10";
const Char_t *patternz100 = "Nz100";
- if ((strstr(name,patternrphi100)) && (strstr(name,patternz100))) {
+ if ((strstr(name.Data(),patternrphi100)) && (strstr(name.Data(),patternz100))) {
fCalibraMode->SetAllTogether(i);
fNbDet = 540;
if (fDebugLevel > 1) {
}
return kTRUE;
}
- if ((strstr(name,patternrphi10)) && (strstr(name,patternz10))) {
+ if ((strstr(name.Data(),patternrphi10)) && (strstr(name.Data(),patternz10))) {
fCalibraMode->SetPerSuperModule(i);
fNbDet = 30;
if (fDebugLevel > 1) {
}
return kTRUE;
}
- if (strstr(name,patternz0)) {
+ if (strstr(name.Data(),patternz0)) {
fCalibraMode->SetNz(i, 0);
if (fDebugLevel > 1) {
AliInfo(Form("fNbDet %d and 0",fNbDet));
}
return kTRUE;
}
- if (strstr(name,patternz1)) {
+ if (strstr(name.Data(),patternz1)) {
fCalibraMode->SetNz(i ,1);
if (fDebugLevel > 1) {
AliInfo(Form("fNbDet %d and 1",fNbDet));
}
return kTRUE;
}
- if (strstr(name,patternz2)) {
+ if (strstr(name.Data(),patternz2)) {
fCalibraMode->SetNz(i ,2);
if (fDebugLevel > 1) {
AliInfo(Form("fNbDet %d and 2",fNbDet));
}
return kTRUE;
}
- if (strstr(name,patternz3)) {
- fCalibraMode->SetNz(i ,3);
- if (fDebugLevel > 1) {
- AliInfo(Form("fNbDet %d and 3",fNbDet));
+ if (strstr(name.Data(),patternz3)) {
+ fCalibraMode->SetNz(i ,3);
+ if (fDebugLevel > 1) {
+ AliInfo(Form("fNbDet %d and 3",fNbDet));
+ }
+ return kTRUE;
+ }
+ if (strstr(name.Data(),patternz4)) {
+ fCalibraMode->SetNz(i ,4);
+ if (fDebugLevel > 1) {
+ AliInfo(Form("fNbDet %d and 4",fNbDet));
+ }
+ return kTRUE;
+ }
+
+ if (fDebugLevel > 1) {
+ AliInfo(Form("fNbDet %d and rest",fNbDet));
+ }
+ fCalibraMode->SetNz(i ,0);
+ return kFALSE;
+}
+//______________________________________________________________________
+void AliTRDCalibraFit::RemoveOutliers(Int_t type, Bool_t perdetector){
+ //
+ // Remove the results too far from the mean value and rms
+ // type: 0 gain, 1 vdrift
+ // perdetector
+ //
+
+ Int_t loop = (Int_t) fVectorFit.GetEntriesFast();
+ if(loop != 540) {
+ AliInfo("The Vector Fit is not complete!");
+ return;
+ }
+ Int_t detector = -1;
+ Int_t sector = -1;
+ Float_t value = 0.0;
+
+ /////////////////////////////////
+ // Calculate the mean values
+ ////////////////////////////////
+ // Initialisation
+ ////////////////////////
+ Double_t meanAll = 0.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);
+ if(perdetector){
+ value = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef()[0];
+ if(value > 0.0) {
+ rmsAll += value*value;
+ meanAll += value;
+ countAll++;
+ }
+ }
+ else {
+ Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
+ Int_t colMax = fGeo->GetColMax(GetLayer(detector));
+ for (Int_t row = 0; row < rowMax; row++) {
+ for (Int_t col = 0; col < colMax; col++) {
+ value = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
+ if(value > 0.0) {
+ rmsAll += value*value;
+ meanAll += value;
+ countAll++;
+ }
+
+ } // Col
+ } // Row
+ }
+ }
+ if(countAll > 0) {
+ meanAll = meanAll/countAll;
+ rmsAll = TMath::Sqrt(TMath::Abs(rmsAll/countAll - (meanAll*meanAll)));
+ }
+ //printf("RemoveOutliers: meanAll %f and rmsAll %f\n",meanAll,rmsAll);
+ /////////////////////////////////////////////////
+ // Remove outliers
+ ////////////////////////////////////////////////
+ Double_t defaultvalue = -1.0;
+ if(type==1) defaultvalue = -1.5;
+ for (Int_t k = 0; k < loop; k++) {
+ detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector();
+ 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();
+
+ // remove the results too far away
+ for (Int_t row = 0; row < rowMax; row++) {
+ for (Int_t col = 0; col < colMax; col++) {
+ value = coef[(Int_t)(col*rowMax+row)];
+ if((value > 0.0) && (rmsAll > 0.0) && (TMath::Abs(value-meanAll) > (2*rmsAll))) {
+ coef[(Int_t)(col*rowMax+row)] = defaultvalue;
+ }
+ } // Col
+ } // Row
+ }
+}
+//______________________________________________________________________
+void AliTRDCalibraFit::RemoveOutliers2(Bool_t perdetector){
+ //
+ // Remove the results too far from the mean and rms
+ // perdetector
+ //
+
+ Int_t loop = (Int_t) fVectorFit2.GetEntriesFast();
+ if(loop != 540) {
+ AliInfo("The Vector Fit is not complete!");
+ return;
+ }
+ Int_t detector = -1;
+ Int_t sector = -1;
+ Float_t value = 0.0;
+
+ /////////////////////////////////
+ // Calculate the mean values
+ ////////////////////////////////
+ // Initialisation
+ ////////////////////////
+ Double_t meanAll = 0.0;
+ Double_t rmsAll = 0.0;
+ Int_t countAll = 0;
+ /////////////
+ // compute
+ ////////////
+ for (Int_t k = 0; k < loop; k++) {
+ detector = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetDetector();
+ sector = GetSector(detector);
+ if(perdetector){
+ value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[0];
+ if(value < 70.0) {
+ meanAll += value;
+ rmsAll += value*value;
+ countAll++;
+ }
}
- return kTRUE;
- }
- if (strstr(name,patternz4)) {
- fCalibraMode->SetNz(i ,4);
- if (fDebugLevel > 1) {
- AliInfo(Form("fNbDet %d and 4",fNbDet));
+ else {
+ Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector));
+ Int_t colMax = fGeo->GetColMax(GetLayer(detector));
+ 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)];
+ if(value < 70.0) {
+ rmsAll += value*value;
+ meanAll += value;
+ countAll++;
+ }
+ } // Col
+ } // Row
}
- return kTRUE;
- }
-
- if (fDebugLevel > 1) {
- AliInfo(Form("fNbDet %d and rest",fNbDet));
+ }
+ if(countAll > 0) {
+ meanAll = meanAll/countAll;
+ rmsAll = TMath::Sqrt(TMath::Abs(rmsAll/countAll - (meanAll*meanAll)));
+ }
+ //printf("Remove outliers 2: meanAll %f, rmsAll %f\n",meanAll,rmsAll);
+ /////////////////////////////////////////////////
+ // Remove outliers
+ ////////////////////////////////////////////////
+ for (Int_t k = 0; k < loop; k++) {
+ detector = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetDetector();
+ 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();
+
+ // remove the results too far away
+ for (Int_t row = 0; row < rowMax; row++) {
+ for (Int_t col = 0; col < colMax; col++) {
+ value = coef[(Int_t)(col*rowMax+row)];
+ if((value < 70.0) && (rmsAll > 0.0) && (TMath::Abs(value-meanAll) > (2.5*rmsAll))) {
+ //printf("value outlier %f\n",value);
+ coef[(Int_t)(col*rowMax+row)] = 100.0;
+ }
+ } // Col
+ } // Row
}
- fCalibraMode->SetNz(i ,0);
- return kFALSE;
}
//______________________________________________________________________
void AliTRDCalibraFit::PutMeanValueOtherVectorFit(Int_t ofwhat, Bool_t perdetector){
Double_t meanAll = 0.0;
Double_t meanSupermodule[18];
Double_t meanDetector[540];
+ Double_t rmsAll = 0.0;
+ Double_t rmsSupermodule[18];
+ Double_t rmsDetector[540];
Int_t countAll = 0;
Int_t countSupermodule[18];
Int_t countDetector[540];
for(Int_t sm = 0; sm < 18; sm++){
+ rmsSupermodule[sm] = 0.0;
meanSupermodule[sm] = 0.0;
countSupermodule[sm] = 0;
}
for(Int_t det = 0; det < 540; det++){
+ rmsDetector[det] = 0.0;
meanDetector[det] = 0.0;
countDetector[det] = 0;
}
+ ////////////
// compute
////////////
for (Int_t k = 0; k < loop; k++) {
if(perdetector){
value = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef()[0];
if(value > 0.0) {
+ rmsDetector[detector] += value*value;
meanDetector[detector] += value;
countDetector[detector]++;
+ rmsSupermodule[sector] += value*value;
meanSupermodule[sector] += value;
countSupermodule[sector]++;
+ rmsAll += value*value;
meanAll += value;
countAll++;
}
for (Int_t col = 0; col < colMax; col++) {
value = ((AliTRDFitInfo *) fVectorFit.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
if(value > 0.0) {
+ rmsDetector[detector] += value*value;
meanDetector[detector] += value;
countDetector[detector]++;
+ rmsSupermodule[sector] += value*value;
meanSupermodule[sector] += value;
countSupermodule[sector]++;
+ rmsAll += value*value;
meanAll += value;
countAll++;
}
} // Row
}
}
- if(countAll > 0) meanAll = meanAll/countAll;
+ if(countAll > 0) {
+ meanAll = meanAll/countAll;
+ rmsAll = TMath::Abs(rmsAll/countAll - (meanAll*meanAll));
+ }
for(Int_t sm = 0; sm < 18; sm++){
- if(countSupermodule[sm] > 0) meanSupermodule[sm] = meanSupermodule[sm]/countSupermodule[sm];
+ if(countSupermodule[sm] > 0) {
+ meanSupermodule[sm] = meanSupermodule[sm]/countSupermodule[sm];
+ rmsSupermodule[sm] = TMath::Abs(rmsSupermodule[sm]/countSupermodule[sm] - (meanSupermodule[sm]*meanSupermodule[sm]));
+ }
}
for(Int_t det = 0; det < 540; det++){
- if(countDetector[det] > 0) meanDetector[det] = meanDetector[det]/countDetector[det];
+ if(countDetector[det] > 0) {
+ meanDetector[det] = meanDetector[det]/countDetector[det];
+ rmsDetector[det] = TMath::Abs(rmsDetector[det]/countDetector[det] - (meanDetector[det]*meanDetector[det]));
+ }
}
+ //printf("Put mean value, meanAll %f, rmsAll %f\n",meanAll,rmsAll);
+ ///////////////////////////////////////////////
// Put the mean value for the no-fitted
/////////////////////////////////////////////
for (Int_t k = 0; k < loop; k++) {
for (Int_t col = 0; col < colMax; col++) {
value = coef[(Int_t)(col*rowMax+row)];
if(value < 0.0) {
- if((ofwhat == 0) && (meanAll > 0.0)) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanAll);
+ if((ofwhat == 0) && (meanAll > 0.0) && (countAll > 15)) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanAll);
if(ofwhat == 1){
- if(meanDetector[detector] > 0.0) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanDetector[detector]);
- else if(meanSupermodule[sector] > 0.0) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanSupermodule[sector]);
- else if(meanAll > 0.0) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanAll);
+ if((meanDetector[detector] > 0.0) && (countDetector[detector] > 20)) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanDetector[detector]);
+ else if((meanSupermodule[sector] > 0.0) && (countSupermodule[sector] > 15)) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanSupermodule[sector]);
+ else if((meanAll > 0.0) && (countAll > 15)) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanAll);
}
}
// Debug
} // Col
} // Row
}
-
}
//______________________________________________________________________
void AliTRDCalibraFit::PutMeanValueOtherVectorFit2(Int_t ofwhat, Bool_t perdetector){
// Initialisation
////////////////////////
Double_t meanAll = 0.0;
+ Double_t rmsAll = 0.0;
Double_t meanSupermodule[18];
+ Double_t rmsSupermodule[18];
Double_t meanDetector[540];
+ Double_t rmsDetector[540];
Int_t countAll = 0;
Int_t countSupermodule[18];
Int_t countDetector[540];
for(Int_t sm = 0; sm < 18; sm++){
+ rmsSupermodule[sm] = 0.0;
meanSupermodule[sm] = 0.0;
countSupermodule[sm] = 0;
}
for(Int_t det = 0; det < 540; det++){
+ rmsDetector[det] = 0.0;
meanDetector[det] = 0.0;
countDetector[det] = 0;
}
if(perdetector){
value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[0];
if(value < 70.0) {
+ rmsDetector[detector] += value*value;
meanDetector[detector] += value;
countDetector[detector]++;
+ rmsSupermodule[sector] += value*value;
meanSupermodule[sector] += value;
countSupermodule[sector]++;
meanAll += value;
+ rmsAll += value*value;
countAll++;
}
}
for (Int_t col = 0; col < colMax; col++) {
value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[(Int_t)(col*rowMax+row)];
if(value < 70.0) {
+ rmsDetector[detector] += value*value;
meanDetector[detector] += value;
countDetector[detector]++;
+ rmsSupermodule[sector] += value*value;
meanSupermodule[sector] += value;
countSupermodule[sector]++;
+ rmsAll += value*value;
meanAll += value;
countAll++;
}
} // Row
}
}
- if(countAll > 0) meanAll = meanAll/countAll;
+ if(countAll > 0) {
+ meanAll = meanAll/countAll;
+ rmsAll = TMath::Abs(rmsAll/countAll - (meanAll*meanAll));
+ }
for(Int_t sm = 0; sm < 18; sm++){
- if(countSupermodule[sm] > 0) meanSupermodule[sm] = meanSupermodule[sm]/countSupermodule[sm];
+ if(countSupermodule[sm] > 0) {
+ meanSupermodule[sm] = meanSupermodule[sm]/countSupermodule[sm];
+ rmsSupermodule[sm] = TMath::Abs(rmsSupermodule[sm]/countSupermodule[sm] - (meanSupermodule[sm]*meanSupermodule[sm]));
+ }
}
for(Int_t det = 0; det < 540; det++){
- if(countDetector[det] > 0) meanDetector[det] = meanDetector[det]/countDetector[det];
+ if(countDetector[det] > 0) {
+ meanDetector[det] = meanDetector[det]/countDetector[det];
+ rmsDetector[det] = TMath::Abs(rmsDetector[det]/countDetector[det] - (meanDetector[det]*meanDetector[det]));
+ }
}
+ //printf("Put mean value 2: meanAll %f, rmsAll %f\n",meanAll,rmsAll);
+ ////////////////////////////////////////////
// Put the mean value for the no-fitted
/////////////////////////////////////////////
for (Int_t k = 0; k < loop; k++) {
for (Int_t col = 0; col < colMax; col++) {
value = coef[(Int_t)(col*rowMax+row)];
if(value > 70.0) {
- if((ofwhat == 0) && (meanAll > 0.0)) coef[(Int_t)(col*rowMax+row)] = meanAll+100.0;
+ if((ofwhat == 0) && (meanAll > -3.0) && (countAll > 15)) coef[(Int_t)(col*rowMax+row)] = meanAll+100.0;
if(ofwhat == 1){
- if(meanDetector[detector] > 0.0) coef[(Int_t)(col*rowMax+row)] = meanDetector[detector]+100.0;
- else if(meanSupermodule[sector] > 0.0) coef[(Int_t)(col*rowMax+row)] = meanSupermodule[sector]+100.0;
- else if(meanAll > 0.0) coef[(Int_t)(col*rowMax+row)] = meanAll+100.0;
+ if((meanDetector[detector] > -3.0) && (countDetector[detector] > 20)) coef[(Int_t)(col*rowMax+row)] = meanDetector[detector]+100.0;
+ else if((meanSupermodule[sector] > -3.0) && (countSupermodule[sector] > 15)) coef[(Int_t)(col*rowMax+row)] = meanSupermodule[sector]+100.0;
+ else if((meanAll > -3.0) && (countAll > 15)) coef[(Int_t)(col*rowMax+row)] = meanAll+100.0;
}
}
// Debug
// Create the DetObject
AliTRDCalDet *object = new AliTRDCalDet("ChamberGainFactor","GainFactor (detector value)");
+ fScaleGain = scaleFitFactor;
Int_t loop = (Int_t) vectorFit->GetEntriesFast();
if(loop != 540) AliInfo("The Vector Fit is not complete!");
} // Row
if(count > 0) mean = mean/count;
}
+ if(mean < 0.1) mean = 0.1;
object->SetValue(detector,mean);
}
Float_t min = 100.0;
if(perdetector){
value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
+ //printf("Create det object %f for %d\n",value,k);
// check successful
if(value > 70.0) value = value-100.0;
//
if(count > 0) mean = mean/count;
*/
value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0];
- object->SetValue(detector,-TMath::Abs(value));
+ if(value > 70.0) value = value-100.0;
+ object->SetValue(detector,value);
+ }
+
+ 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;
detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector();
AliTRDCalROC *calROC = object->GetCalROC(detector);
Float_t mean = detobject->GetValue(detector);
- if(mean == 0) continue;
+ if(TMath::Abs(mean) <= 0.0000000001) continue;
Int_t rowMax = calROC->GetNrows();
Int_t colMax = calROC->GetNcols();
for (Int_t row = 0; row < rowMax; row++) {
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++) {
fCurrentCoefDetector2[k] = 0.0;
}
- //printf("test0\n");
-
- AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
- if (!cal) {
- AliInfo("Could not get calibDB");
- return kFALSE;
- }
+ if((!fCalDetVdriftUsed) || (!fCalDetExBUsed)) return kFALSE;
+
+ return kTRUE;
+}
+//____________Functions for initialising the AliTRDCalibraFit in the code_________
+Bool_t AliTRDCalibraFit::InitFitExbAlt()
+{
+ //
+ // Init the fCalDet, fVectorFit fCurrentCoefDetector
+ //
- //Get the CalDet object
- if(fAccCDB){
- if(fCalDet) delete fCalDet;
- if(fCalDet2) delete fCalDet2;
- fCalDet = new AliTRDCalDet(*(cal->GetVdriftDet()));
- //printf("test1\n");
- fCalDet2 = new AliTRDCalDet("lorentz angle tan","lorentz angle tan (detector value)");
- //printf("test2\n");
- for(Int_t k = 0; k < 540; k++){
- fCalDet2->SetValue(k,AliTRDCommonParam::Instance()->GetOmegaTau(fCalDet->GetValue(k)));
- }
- //printf("test3\n");
- }
- else{
- Float_t devalue = 1.5;
- Float_t devalue2 = AliTRDCommonParam::Instance()->GetOmegaTau(1.5);
- if(fCalDet) delete fCalDet;
- if(fCalDet2) delete fCalDet2;
- //printf("test1\n");
- fCalDet = new AliTRDCalDet("ChamberVdrift","TRD drift velocities (detector value)");
- fCalDet2 = new AliTRDCalDet("lorentz angle tan","lorentz angle tan (detector value)");
- //printf("test2\n");
- for(Int_t k = 0; k < 540; k++){
- fCalDet->SetValue(k,devalue);
- fCalDet2->SetValue(k,devalue2);
- }
- //printf("test3\n");
+ 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)
{
if(fDebugLevel == 1) {
fCountDet = 0;
fCalibraMode->CalculXBins(fCountDet,i);
- while(fCalibraMode->GetXbins(i) <=fFitVoir){
+ if((fCalibraMode->GetNz(i)!=100) && (fCalibraMode->GetNrphi(i)!=100)){
+ while(fCalibraMode->GetXbins(i) <=fFitVoir){
+ fCountDet++;
+ fCalibraMode->CalculXBins(fCountDet,i);
+ //printf("GetXBins %d\n",fCalibraMode->GetXbins(i));
+ }
+ }
+ else {
fCountDet++;
- fCalibraMode->CalculXBins(fCountDet,i);
- }
+ }
fCount = fCalibraMode->GetXbins(i);
fCountDet--;
// Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi
else if (fNbDet > 0){
Int_t firstdetector = fCountDet;
Int_t lastdetector = fCountDet+fNbDet;
- AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted"
- ,idect,firstdetector,lastdetector));
+ //AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted",idect,firstdetector,lastdetector));
// loop over detectors
for(Int_t det = firstdetector; det < lastdetector; det++){
}
else {
- AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted"
- ,idect-(fCount-(fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))),fCountDet));
+//AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted",idect-(fCount-(fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))),fCountDet));
// Calcul the coef from the database choosen
CalculChargeCoefMean(kFALSE);
Int_t firstdetector = fCountDet;
Int_t lastdetector = fCountDet+fNbDet;
- AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted"
- ,idect,firstdetector,lastdetector));
+//AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted",idect,firstdetector,lastdetector));
// loop over detectors
for(Int_t det = firstdetector; det < lastdetector; det++){
}
else {
- AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted"
- ,idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))),fCountDet));
+//AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted",idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))),fCountDet));
CalculVdriftCoefMean();
CalculT0CoefMean();
Int_t firstdetector = fCountDet;
Int_t lastdetector = fCountDet+fNbDet;
- AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted"
- ,idect,firstdetector,lastdetector));
+// AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted",idect,firstdetector,lastdetector));
// loop over detectors
for(Int_t det = firstdetector; det < lastdetector; det++){
}
else {
- AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted"
- ,idect-(fCount-(fCalibraMode->GetNfragZ(2)*fCalibraMode->GetNfragRphi(2))),fCountDet));
+// AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted",idect-(fCount-(fCalibraMode->GetNfragZ(2)*fCalibraMode->GetNfragRphi(2))),fCountDet));
CalculPRFCoefMean();
for (Int_t k = 0; k < factor; k++) {
fCurrentCoefDetector[k] = -TMath::Abs(fCurrentCoef[1]);
// should be negative
- fCurrentCoefDetector2[k] = +TMath::Abs(fCurrentCoef2[1]);
+ fCurrentCoefDetector2[k] = fCurrentCoef2[1]+100.0;
}
- //Put default opposite sign
+ //Put default opposite sign only for vdrift
fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
fCurrentCoefE = 0.0;
- fCurrentCoef2[0] = +TMath::Abs(fCurrentCoef2[1]);
+ fCurrentCoef2[0] = fCurrentCoef2[1]+100.0;
fCurrentCoefE2 = 0.0;
FillFillLinearFitter();
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)
{
if (fNbDet > 0){
Int_t firstdetector = fCountDet;
Int_t lastdetector = fCountDet+fNbDet;
- AliInfo(Form("The element %d containing the detectors %d to %d has been fitted"
- ,idect,firstdetector,lastdetector));
+ // AliInfo(Form("The element %d containing the detectors %d to %d has been fitted",idect,firstdetector,lastdetector));
// loop over detectors
for(Int_t det = firstdetector; det < lastdetector; det++){
Int_t firstdetector = fCountDet;
Int_t lastdetector = fCountDet+fNbDet;
- AliInfo(Form("The element %d containing the detectors %d to %d has been fitted"
- ,idect,firstdetector,lastdetector));
+// AliInfo(Form("The element %d containing the detectors %d to %d has been fitted",idect,firstdetector,lastdetector));
// loop over detectors
for(Int_t det = firstdetector; det < lastdetector; det++){
Int_t firstdetector = fCountDet;
Int_t lastdetector = fCountDet+fNbDet;
- AliInfo(Form("The element %d containing the detectors %d to %d has been fitted"
- ,idect,firstdetector,lastdetector));
+// AliInfo(Form("The element %d containing the detectors %d to %d has been fitted",idect,firstdetector,lastdetector));
// loop over detectors
for(Int_t det = firstdetector; det < lastdetector; det++){
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++) {
"r="<<r<<
"tiltangle="<<tiltangle<<
"vf="<<vf<<
- "vs="<<vs<<
+ "vs="<<vs<<
+ "vfE="<<vfE<<
+ "lorentzangler="<<lorentzangler<<
+ "Elorentzangler="<<elorentzangler<<
+ "lorentzangles="<<lorentzangles<<
+ "\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<<
- "lorentzangler="<<lorentzangler<<
- "Elorentzangler="<<elorentzangler<<
- "lorentzangles="<<lorentzangles<<
"\n";
}
// For the detector fCountDet, mean drift velocity and tan lorentzangle
//
- fCurrentCoef[1] = fCalDet->GetValue(fCountDet);
- fCurrentCoef2[1] = fCalDet2->GetValue(fCountDet);
+ fCurrentCoef[1] = fCalDetVdriftUsed->GetValue(fCountDet);
+ fCurrentCoef2[1] = fCalDetExBUsed->GetValue(fCountDet);
return kTRUE;
}
Float_t l3P2am = pentea->GetFunction("pol2")->GetParameter(2);
Float_t l3P1amE = pentea->GetFunction("pol2")->GetParError(1);
Float_t l3P2amE = pentea->GetFunction("pol2")->GetParError(2);
- if (l3P2am != 0) {
+ if (TMath::Abs(l3P2am) > 0.00000001) {
fPhd[0] = -(l3P1am / (2 * l3P2am));
}
if(!fTakeTheMaxPH){
- if((l3P1am != 0.0) && (l3P2am != 0.0)){
+ if((TMath::Abs(l3P1am) > 0.0000000001) && (TMath::Abs(l3P2am) > 0.00000000001)){
fCurrentCoefE2 = (l3P1amE/l3P1am + l3P2amE/l3P2am)*fPhd[0];
}
}
Float_t l3P2amf = projPH->GetFunction("pol2")->GetParameter(2);
Float_t l3P1amfE = projPH->GetFunction("pol2")->GetParError(1);
Float_t l3P2amfE = projPH->GetFunction("pol2")->GetParError(2);
- if (l3P2amf != 0) {
+ if (TMath::Abs(l3P2amf) > 0.00000000001) {
fPhd[1] = -(l3P1amf / (2 * l3P2amf));
}
- if((l3P1amf != 0.0) && (l3P2amf != 0.0)){
+ if((TMath::Abs(l3P1amf) > 0.0000000001) && (TMath::Abs(l3P2amf) > 0.000000001)){
fCurrentCoefE = (l3P1amfE/l3P1amf + l3P2amfE/l3P2amf)*fPhd[1];
}
if(fTakeTheMaxPH){
Float_t l3P2dr = pente->GetFunction("pol2")->GetParameter(2);
Float_t l3P1drE = pente->GetFunction("pol2")->GetParError(1);
Float_t l3P2drE = pente->GetFunction("pol2")->GetParError(2);
- if (l3P2dr != 0) {
+ if (TMath::Abs(l3P2dr) > 0.00000001) {
fPhd[2] = -(l3P1dr / (2 * l3P2dr));
}
- if((l3P1dr != 0.0) && (l3P2dr != 0.0)){
+ if((TMath::Abs(l3P1dr) > 0.0000000001) && (TMath::Abs(l3P2dr) > 0.00000000001)){
fCurrentCoefE += (l3P1drE/l3P1dr + l3P2drE/l3P2dr)*fPhd[2];
}
Float_t fPhdt0 = 0.0;
if (fPhdt0 >= 0.0) {
fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins;
- if (fCurrentCoef2[0] < -1.0) {
+ if (fCurrentCoef2[0] < -3.0) {
fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
}
}
//
// Slope methode but with polynomes de Lagrange
//
-
+
// Constants
const Float_t kDrWidth = AliTRDgeometry::DrThick();
Int_t binmax = 0;
Int_t binmin = 0;
- Double_t *x = new Double_t[5];
- Double_t *y = new Double_t[5];
+ //Double_t *x = new Double_t[5];
+ //Double_t *y = new Double_t[5];
+ Double_t x[5];
+ Double_t y[5];
x[0] = 0.0;
x[1] = 0.0;
x[2] = 0.0;
TF1 * polynome = 0x0;
TF1 * polynomea = 0x0;
TF1 * polynomeb = 0x0;
- Double_t *c = 0x0;
+ Double_t c0 = 0.0;
+ Double_t c1 = 0.0;
+ Double_t c2 = 0.0;
+ Double_t c3 = 0.0;
+ Double_t c4 = 0.0;
// Some variables
TAxis *xpph = projPH->GetXaxis();
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)));
y[0] = pentea->GetBinContent(binmax-2);
y[1] = pentea->GetBinContent(binmax-1);
y[2] = pentea->GetBinContent(binmax);
- c = CalculPolynomeLagrange2(x,y);
- AliInfo("At the limit for beginning!");
+ CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
+ //AliInfo("At the limit for beginning!");
break;
case 2:
minnn = pentea->GetBinCenter(binmax-2);
y[1] = pentea->GetBinContent(binmax-1);
y[2] = pentea->GetBinContent(binmax);
y[3] = pentea->GetBinContent(binmax+1);
- c = CalculPolynomeLagrange3(x,y);
+ CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
break;
default:
switch(binmax){
y[0] = pentea->GetBinContent(binmax);
y[1] = pentea->GetBinContent(binmax+1);
y[2] = pentea->GetBinContent(binmax+2);
- c = CalculPolynomeLagrange2(x,y);
+ CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
break;
case 2:
minnn = pentea->GetBinCenter(binmax-1);
y[1] = pentea->GetBinContent(binmax);
y[2] = pentea->GetBinContent(binmax+1);
y[3] = pentea->GetBinContent(binmax+2);
- c = CalculPolynomeLagrange3(x,y);
+ CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
break;
default:
minnn = pentea->GetBinCenter(binmax-2);
y[2] = pentea->GetBinContent(binmax);
y[3] = pentea->GetBinContent(binmax+1);
y[4] = pentea->GetBinContent(binmax+2);
- c = CalculPolynomeLagrange4(x,y);
+ CalculPolynomeLagrange4(x,y,c0,c1,c2,c3,c4);
break;
}
break;
if(put) {
polynomeb = new TF1("polb","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",minnn,maxxx);
- polynomeb->SetParameters(c[0],c[1],c[2],c[3],c[4]);
+ polynomeb->SetParameters(c0,c1,c2,c3,c4);
Double_t step = (maxxx-minnn)/10000;
Double_t l = minnn;
fPhd[0] = placemaximum;
}
+ /////////////////////////////
// Amplification region
+ /////////////////////////////
binmax = 0;
ju = 0;
for (Int_t kbin = 1; kbin < projPH->GetNbinsX(); kbin ++) {
y[0] = projPH->GetBinContent(binmax-2);
y[1] = projPH->GetBinContent(binmax-1);
y[2] = projPH->GetBinContent(binmax);
- c = CalculPolynomeLagrange2(x,y);
+ CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
//AliInfo("At the limit for the drift!");
break;
case 1:
y[1] = projPH->GetBinContent(binmax-1);
y[2] = projPH->GetBinContent(binmax);
y[3] = projPH->GetBinContent(binmax+1);
- c = CalculPolynomeLagrange3(x,y);
+ CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
break;
default:
switch(binmax)
y[0] = projPH->GetBinContent(binmax);
y[1] = projPH->GetBinContent(binmax+1);
y[2] = projPH->GetBinContent(binmax+2);
- c = CalculPolynomeLagrange2(x,y);
+ CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
break;
case 2:
minn = projPH->GetBinCenter(binmax-1);
y[1] = projPH->GetBinContent(binmax);
y[2] = projPH->GetBinContent(binmax+1);
y[3] = projPH->GetBinContent(binmax+2);
- c = CalculPolynomeLagrange3(x,y);
+ CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
break;
default:
minn = projPH->GetBinCenter(binmax-2);
y[2] = projPH->GetBinContent(binmax);
y[3] = projPH->GetBinContent(binmax+1);
y[4] = projPH->GetBinContent(binmax+2);
- c = CalculPolynomeLagrange4(x,y);
+ CalculPolynomeLagrange4(x,y,c0,c1,c2,c3,c4);
break;
}
break;
if(put) {
polynomea = new TF1("pola","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",minn,maxx);
- polynomea->SetParameters(c[0],c[1],c[2],c[3],c[4]);
+ polynomea->SetParameters(c0,c1,c2,c3,c4);
Double_t step = (maxx-minn)/1000;
Double_t l = minn;
}
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;
- AliInfo("Put the binmax from 1 to 2 to enable the fit");
+ 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;
+ //AliInfo("Too many fluctuations at the end!");
+ putd = kFALSE;
}
if((projPH->GetBinContent(binmin)+projPH->GetBinError(binmin)) > (projPH->GetBinContent(binmin-1))) {
- AliInfo("Too many fluctuations at the end!");
- put = kFALSE;
+ //AliInfo("Too many fluctuations at the end!");
+ putd = kFALSE;
}
- if(pente->GetBinContent(binmin+1)==0){
- AliInfo("No entries for the next bin!");
+ if(TMath::Abs(pente->GetBinContent(binmin+1)) <= 0.0000000000001){
+ //AliInfo("No entries for the next bin!");
pente->SetBinContent(binmin,0);
if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin();
}
y[3] = pente->GetBinContent(binmin+1);
y[4] = pente->GetBinContent(binmin+2);
//Calcul the polynome de Lagrange
- c = CalculPolynomeLagrange4(x,y);
+ CalculPolynomeLagrange4(x,y,c0,c1,c2,c3,c4);
//richtung +/-
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;
+ //AliInfo("polynome 4 false 2");
+ putd = kFALSE;
}
// poly 3
if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1)) &&
y[2] = pente->GetBinContent(binmin);
y[3] = pente->GetBinContent(binmin+1);
//Calcul the polynome de Lagrange
- c = CalculPolynomeLagrange3(x,y);
+ CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
//richtung +: nothing
//richtung -
if((pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) {
y[2] = pente->GetBinContent(binmin+1);
y[3] = pente->GetBinContent(binmin+2);
//Calcul the polynome de Lagrange
- c = CalculPolynomeLagrange3(x,y);
+ CalculPolynomeLagrange3(x,y,c0,c1,c2,c3,c4);
//richtung +
if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1))) {
//AliInfo("polynome 3+ case 2");
y[1] = pente->GetBinContent(binmin+1);
y[2] = pente->GetBinContent(binmin+2);
//Calcul the polynome de Lagrange
- c = CalculPolynomeLagrange2(x,y);
+ CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
//richtung +
if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1))) {
//AliInfo("polynome 2+ false");
- put = kFALSE;
+ putd = kFALSE;
}
}
//pol2 case 2
y[1] = pente->GetBinContent(binmin);
y[2] = pente->GetBinContent(binmin+1);
//Calcul the polynome de Lagrange
- c = CalculPolynomeLagrange2(x,y);
+ CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
//richtung +: nothing
//richtung -: nothing
}
y[1] = pente->GetBinContent(binmin-1);
y[2] = pente->GetBinContent(binmin);
//Calcul the polynome de Lagrange
- c = CalculPolynomeLagrange2(x,y);
+ CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4);
//AliInfo("At the limit for the drift!");
//fluctuation too big!
//richtung +: nothing
//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;
- AliInfo("At the limit for the drift and not usable!");
+ 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;
- AliInfo("For the drift...problem!");
+ 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;
- AliInfo("For the drift...problem!");
+ 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(c[0],c[1],c[2],c[3],c[4]);
+ polynome->SetParameters(c0,c1,c2,c3,c4);
//AliInfo(Form("GetMinimum of the function %f",polynome->GetMinimumX()));
Double_t step = (max-min)/1000;
Double_t l = min;
if ((fPhd[2] > fPhd[0]) &&
(fPhd[2] > fPhd[1]) &&
(fPhd[1] > fPhd[0]) &&
- (put)) {
+ (put) &&
+ (putd)) {
fCurrentCoef[0] = (kDrWidth) / (fPhd[2]-fPhd[1]);
- fNumberFitSuccess++;
+ if(fCurrentCoef[0] > 2.5) fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
+ else fNumberFitSuccess++;
if (fPhdt0 >= 0.0) {
fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins;
- if (fCurrentCoef2[0] < -1.0) {
+ //printf("Value of timeoffset %f\n",fCurrentCoef2[0]);
+ if (fCurrentCoef2[0] < -3.0) {
fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
}
}
}
}
else {
+ ////printf("Put default %f\n",-TMath::Abs(fCurrentCoef[1]));
+ fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
+
if((fPhd[1] > fPhd[0]) &&
(put)) {
if (fPhdt0 >= 0.0) {
fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins;
- if (fCurrentCoef2[0] < -1.0) {
+ if (fCurrentCoef2[0] < -3.0) {
fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
}
+ else fNumberFitSuccess++;
}
else {
fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
}
}
else{
- fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
//printf("Fit failed!\n");
}
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();
pente->Draw();
}
else {
- if(pentea) delete pentea;
- if(pente) delete pente;
+ delete pentea;
+ delete pente;
if(polynome) delete polynome;
if(polynomea) delete polynomea;
if(polynomeb) delete polynomeb;
- if(x) delete [] x;
- if(y) delete [] y;
- if(c) delete [] c;
+ //if(x) delete [] x;
+ //if(y) delete [] y;
if(line) delete line;
}
//Provisoire
//if(fCurrentCoef[0] > 1.7) fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]);
//if((fCurrentCoef2[0] > 2.6) || (fCurrentCoef2[0] < 2.1)) fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0;
-
+ //printf("Value of timeoffset final %f\n",fCurrentCoef2[0]);
projPH->SetDirectory(0);
}
Double_t ret = FitGausMI(arraye, arraym, arrayme, nBins, xMin, xMax,¶m);
- if(ret == -4){
+ if(TMath::Abs(ret+4) <= 0.000000001){
fCurrentCoef[0] = -fCurrentCoef[1];
return kFALSE;
}
if((valueI - arrayme[ibin]) > 0.0) errorn = TMath::Log((valueI - arrayme[ibin])/valueI);
error = TMath::Max(TMath::Abs(errorm),TMath::Abs(errorn));
}
- if(error == 0.0) continue;
+ if(TMath::Abs(error) < 0.000000001) continue;
val = TMath::Log(Float_t(valueI));
fitter.AddPoint(&xcenter,val,error);
npoints++;
if (!param) param = new TVectorD(3);
- if(par[2] == 0.0) return -4.0;
+ if(TMath::Abs(par[2]) <= 0.000000001) return -4.0;
Double_t x = TMath::Sqrt(TMath::Abs(-2*par[2]));
Double_t deltax = (fitter.GetParError(2))/x;
Double_t errorparam2 = TMath::Abs(deltax)/(x*x);
Double_t ermin0 = 0.0;
//Double_t prfe0 = 0.0;
Double_t prf0 = 0.0;
- if((pars0[2] > 0.0) && (pars0[1] != 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]));
prf0 = pars0[0]+pars0[1]*min0+pars0[2]*min0*min0;
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 *AliTRDCalibraFit::CalculPolynomeLagrange2(const Double_t *x, const Double_t *y) const
+void AliTRDCalibraFit::CalculPolynomeLagrange2(const Double_t *x, const Double_t *y, Double_t &c0, Double_t &c1, Double_t &c2, Double_t &c3, Double_t &c4) const
{
//
// Calcul the coefficients of the polynome passant par ces trois points de degre 2
//
- Double_t *c = new Double_t[5];
Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2]));
Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2]));
Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1]));
- c[4] = 0.0;
- c[3] = 0.0;
- c[2] = x0+x1+x2;
- c[1] = -(x0*(x[1]+x[2])+x1*(x[0]+x[2])+x2*(x[0]+x[1]));
- c[0] = x0*x[1]*x[2]+x1*x[0]*x[2]+x2*x[0]*x[1];
-
- return c;
-
+ c4 = 0.0;
+ c3 = 0.0;
+ c2 = x0+x1+x2;
+ c1 = -(x0*(x[1]+x[2])+x1*(x[0]+x[2])+x2*(x[0]+x[1]));
+ c0 = x0*x[1]*x[2]+x1*x[0]*x[2]+x2*x[0]*x[1];
}
//_____________________________________________________________________________
-Double_t *AliTRDCalibraFit::CalculPolynomeLagrange3(const Double_t *x, const Double_t *y) const
+void AliTRDCalibraFit::CalculPolynomeLagrange3(const Double_t *x, const Double_t *y, Double_t &c0, Double_t &c1, Double_t &c2, Double_t &c3, Double_t &c4) const
{
//
// Calcul the coefficients of the polynome passant par ces quatre points de degre 3
//
- Double_t *c = new Double_t[5];
Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2])*(x[0]-x[3]));
Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2])*(x[1]-x[3]));
Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1])*(x[2]-x[3]));
Double_t x3 = y[3]/((x[3]-x[0])*(x[3]-x[1])*(x[3]-x[2]));
- c[4] = 0.0;
- c[3] = x0+x1+x2+x3;
- c[2] = -(x0*(x[1]+x[2]+x[3])
+ c4 = 0.0;
+ c3 = x0+x1+x2+x3;
+ c2 = -(x0*(x[1]+x[2]+x[3])
+x1*(x[0]+x[2]+x[3])
+x2*(x[0]+x[1]+x[3])
+x3*(x[0]+x[1]+x[2]));
- c[1] = (x0*(x[1]*x[2]+x[1]*x[3]+x[2]*x[3])
+ c1 = (x0*(x[1]*x[2]+x[1]*x[3]+x[2]*x[3])
+x1*(x[0]*x[2]+x[0]*x[3]+x[2]*x[3])
+x2*(x[0]*x[1]+x[0]*x[3]+x[1]*x[3])
+x3*(x[0]*x[1]+x[0]*x[2]+x[1]*x[2]));
- c[0] = -(x0*x[1]*x[2]*x[3]
+ c0 = -(x0*x[1]*x[2]*x[3]
+x1*x[0]*x[2]*x[3]
+x2*x[0]*x[1]*x[3]
+x3*x[0]*x[1]*x[2]);
- return c;
-
-
}
//_____________________________________________________________________________
-Double_t *AliTRDCalibraFit::CalculPolynomeLagrange4(const Double_t *x, const Double_t *y) const
+void AliTRDCalibraFit::CalculPolynomeLagrange4(const Double_t *x, const Double_t *y, Double_t &c0, Double_t &c1, Double_t &c2, Double_t &c3, Double_t &c4) const
{
//
// Calcul the coefficients of the polynome passant par ces cinqs points de degre 4
//
- Double_t *c = new Double_t[5];
Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2])*(x[0]-x[3])*(x[0]-x[4]));
Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2])*(x[1]-x[3])*(x[1]-x[4]));
Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1])*(x[2]-x[3])*(x[2]-x[4]));
Double_t x4 = y[4]/((x[4]-x[0])*(x[4]-x[1])*(x[4]-x[2])*(x[4]-x[3]));
- c[4] = x0+x1+x2+x3+x4;
- c[3] = -(x0*(x[1]+x[2]+x[3]+x[4])
+ c4 = x0+x1+x2+x3+x4;
+ c3 = -(x0*(x[1]+x[2]+x[3]+x[4])
+x1*(x[0]+x[2]+x[3]+x[4])
+x2*(x[0]+x[1]+x[3]+x[4])
+x3*(x[0]+x[1]+x[2]+x[4])
+x4*(x[0]+x[1]+x[2]+x[3]));
- c[2] = (x0*(x[1]*x[2]+x[1]*x[3]+x[1]*x[4]+x[2]*x[3]+x[2]*x[4]+x[3]*x[4])
+ c2 = (x0*(x[1]*x[2]+x[1]*x[3]+x[1]*x[4]+x[2]*x[3]+x[2]*x[4]+x[3]*x[4])
+x1*(x[0]*x[2]+x[0]*x[3]+x[0]*x[4]+x[2]*x[3]+x[2]*x[4]+x[3]*x[4])
+x2*(x[0]*x[1]+x[0]*x[3]+x[0]*x[4]+x[1]*x[3]+x[1]*x[4]+x[3]*x[4])
+x3*(x[0]*x[1]+x[0]*x[2]+x[0]*x[4]+x[1]*x[2]+x[1]*x[4]+x[2]*x[4])
+x4*(x[0]*x[1]+x[0]*x[2]+x[0]*x[3]+x[1]*x[2]+x[1]*x[3]+x[2]*x[3]));
- c[1] = -(x0*(x[1]*x[2]*x[3]+x[1]*x[2]*x[4]+x[1]*x[3]*x[4]+x[2]*x[3]*x[4])
+ c1 = -(x0*(x[1]*x[2]*x[3]+x[1]*x[2]*x[4]+x[1]*x[3]*x[4]+x[2]*x[3]*x[4])
+x1*(x[0]*x[2]*x[3]+x[0]*x[2]*x[4]+x[0]*x[3]*x[4]+x[2]*x[3]*x[4])
+x2*(x[0]*x[1]*x[3]+x[0]*x[1]*x[4]+x[0]*x[3]*x[4]+x[1]*x[3]*x[4])
+x3*(x[0]*x[1]*x[2]+x[0]*x[1]*x[4]+x[0]*x[2]*x[4]+x[1]*x[2]*x[4])
+x4*(x[0]*x[1]*x[2]+x[0]*x[1]*x[3]+x[0]*x[2]*x[3]+x[1]*x[2]*x[3]));
- c[0] = (x0*x[1]*x[2]*x[3]*x[4]
+ c0 = (x0*x[1]*x[2]*x[3]*x[4]
+x1*x[0]*x[2]*x[3]*x[4]
+x2*x[0]*x[1]*x[3]*x[4]
+x3*x[0]*x[1]*x[2]*x[4]
+x4*x[0]*x[1]*x[2]*x[3]);
- return c;
-
-
}
//_____________________________________________________________________________
void AliTRDCalibraFit::NormierungCharge()
Double_t sigma2 = par2save*par2save;
Double_t sqrt2 = TMath::Sqrt(2.0);
Double_t exp1 = par3save * TMath::Exp(-par3save * (dx - 0.5 * par3save * sigma2))
- * (1.0 - TMath::Erf((par3save * sigma2 - dx) / (sqrt2 * par2save)));
+ * (1.0 - AliMathBase::ErfFast((par3save * sigma2 - dx) / (sqrt2 * par2save)));
Double_t exp2 = par5save * TMath::Exp(-par5save * (dx - 0.5 * par5save * sigma2))
- * (1.0 - TMath::Erf((par5save * sigma2 - dx) / (sqrt2 * par2save)));
+ * (1.0 - AliMathBase::ErfFast((par5save * sigma2 - dx) / (sqrt2 * par2save)));
//return par[0]*(exp1+par[4]*exp2);
return par[0] * (exp1 + 1.00124 * exp2);
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 = 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]) / 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]) / 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;
+
+}
+//_____________________________________________________________________________
+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;
- 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])*TMath::Exp(-par[4]*xx) / 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])*TMath::Exp(-par[4]*xx) / 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;
}
//_____________________________________________________________________________
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