X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TRD%2FAliTRDCalibraFit.cxx;h=a56103a85385ffef8f0ede9661f3c94c9f0dd761;hb=6245a86ff30677709a7e0966963b15d35a07bdd1;hp=eb9c4de9a50f68270230f692aa3cf4bd182f0d17;hpb=ba1aa7a77f63f4ac875b5adb60c2b59950601a11;p=u%2Fmrichter%2FAliRoot.git diff --git a/TRD/AliTRDCalibraFit.cxx b/TRD/AliTRDCalibraFit.cxx index eb9c4de9a50..a56103a8538 100644 --- a/TRD/AliTRDCalibraFit.cxx +++ b/TRD/AliTRDCalibraFit.cxx @@ -64,6 +64,7 @@ #include #include #include +#include #include "AliLog.h" #include "AliMathBase.h" @@ -72,6 +73,7 @@ #include "AliTRDCalibraMode.h" #include "AliTRDCalibraVector.h" #include "AliTRDCalibraVdriftLinearFit.h" +#include "AliTRDCalibraExbAltFit.h" #include "AliTRDcalibDB.h" #include "AliTRDgeometry.h" #include "AliTRDpadPlane.h" @@ -128,6 +130,8 @@ AliTRDCalibraFit::AliTRDCalibraFit() ,fNumberOfBinsExpected(0) ,fMethod(0) ,fBeginFitCharge(3.5) + ,fOutliersFitChargeLow(0.03) + ,fOutliersFitChargeHigh(0.80) ,fFitPHPeriode(1) ,fTakeTheMaxPH(kTRUE) ,fT0Shift0(0.124797) @@ -136,6 +140,7 @@ AliTRDCalibraFit::AliTRDCalibraFit() ,fAccCDB(kFALSE) ,fMinEntries(800) ,fRebin(1) + ,fScaleGain(0.02431) ,fNumberFit(0) ,fNumberFitSuccess(0) ,fNumberEnt(0) @@ -158,6 +163,8 @@ AliTRDCalibraFit::AliTRDCalibraFit() ,fCalROC(0x0) ,fCalDet2(0x0) ,fCalROC2(0x0) + ,fCalDetVdriftUsed(0x0) + ,fCalDetExBUsed(0x0) ,fCurrentCoefDetector(0x0) ,fCurrentCoefDetector2(0x0) ,fVectorFit(0) @@ -187,6 +194,8 @@ AliTRDCalibraFit::AliTRDCalibraFit(const AliTRDCalibraFit &c) ,fNumberOfBinsExpected(c.fNumberOfBinsExpected) ,fMethod(c.fMethod) ,fBeginFitCharge(c.fBeginFitCharge) +,fOutliersFitChargeLow(c.fOutliersFitChargeLow) +,fOutliersFitChargeHigh(c.fOutliersFitChargeHigh) ,fFitPHPeriode(c.fFitPHPeriode) ,fTakeTheMaxPH(c.fTakeTheMaxPH) ,fT0Shift0(c.fT0Shift0) @@ -195,6 +204,7 @@ AliTRDCalibraFit::AliTRDCalibraFit(const AliTRDCalibraFit &c) ,fAccCDB(c.fAccCDB) ,fMinEntries(c.fMinEntries) ,fRebin(c.fRebin) +,fScaleGain(c.fScaleGain) ,fNumberFit(c.fNumberFit) ,fNumberFitSuccess(c.fNumberFitSuccess) ,fNumberEnt(c.fNumberEnt) @@ -217,6 +227,8 @@ AliTRDCalibraFit::AliTRDCalibraFit(const AliTRDCalibraFit &c) ,fCalROC(0x0) ,fCalDet2(0x0) ,fCalROC2(0x0) +,fCalDetVdriftUsed(0x0) +,fCalDetExBUsed(0x0) ,fCurrentCoefDetector(0x0) ,fCurrentCoefDetector2(0x0) ,fVectorFit(0) @@ -243,6 +255,9 @@ AliTRDCalibraFit::AliTRDCalibraFit(const AliTRDCalibraFit &c) 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(); @@ -298,6 +313,8 @@ AliTRDCalibraFit::~AliTRDCalibraFit() 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(); @@ -330,28 +347,6 @@ void AliTRDCalibraFit::DestroyDebugStreamer() 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) @@ -424,8 +419,10 @@ Bool_t AliTRDCalibraFit::AnalyseCH(const TH2I *ch) { case 0: FitMeanW((TH1 *) projch, nentries); break; case 1: FitMean((TH1 *) projch, nentries, mean); break; - case 2: FitCH((TH1 *) projch, mean); break; - case 3: FitBisCH((TH1 *) projch, mean); break; + case 2: FitLandau((TH1 *) projch, mean, nentries); break; + case 3: FitCH((TH1 *) projch, mean, nentries); break; + case 4: FitBisCH((TH1 *) projch, mean, nentries); break; + case 5: FitBisCHEx((TH1 *) projch, mean, nentries); break; default: return kFALSE; } // Fill Infos Fit @@ -441,7 +438,7 @@ Bool_t AliTRDCalibraFit::AnalyseCH(const TH2I *ch) } // 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 { @@ -521,8 +518,10 @@ Bool_t AliTRDCalibraFit::AnalyseCH(AliTRDCalibraVector *calvect) { 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 @@ -534,7 +533,7 @@ Bool_t AliTRDCalibraFit::AnalyseCH(AliTRDCalibraVector *calvect) } // 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 { @@ -544,8 +543,57 @@ Bool_t AliTRDCalibraFit::AnalyseCH(AliTRDCalibraVector *calvect) 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 @@ -554,72 +602,108 @@ Bool_t AliTRDCalibraFit::AnalysePH(const TProfile2D *ph) // A first calibration of T0 is also made using the same method // - // Set the calibration mode - //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) Int_t nbins = ph->GetNbinsX();// time Int_t nybins = ph->GetNbinsY();// calibration group - if (!InitFit(nybins,1)) { - return kFALSE; + + // 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) + // - //printf("Init fit\n"); + // 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()) { 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++) { - //printf("idect = %d\n",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(); //Method choosen - //printf("Method\n"); switch(fMethod) { case 0: FitLagrangePoly((TH1 *) projph); break; @@ -628,15 +712,12 @@ Bool_t AliTRDCalibraFit::AnalysePH(const TProfile2D *ph) 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 { @@ -646,64 +727,82 @@ Bool_t AliTRDCalibraFit::AnalysePH(const TProfile2D *ph) 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(); - TString 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; @@ -712,12 +811,15 @@ Bool_t AliTRDCalibraFit::AnalysePH(AliTRDCalibraVector *calvect) 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 { @@ -804,8 +906,8 @@ Bool_t AliTRDCalibraFit::AnalysePRF(const TProfile2D *prf) 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 { @@ -900,8 +1002,8 @@ Bool_t AliTRDCalibraFit::AnalysePRFMarianFit(const TProfile2D *prf) 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 { @@ -979,7 +1081,7 @@ Bool_t AliTRDCalibraFit::AnalysePRF(AliTRDCalibraVector *calvect) } // 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)); @@ -1082,7 +1184,7 @@ Bool_t AliTRDCalibraFit::AnalysePRFMarianFit(AliTRDCalibraVector *calvect) } // 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)); @@ -1113,18 +1215,18 @@ Bool_t AliTRDCalibraFit::AnalyseLinearFitters(AliTRDCalibraVdriftLinearFit *cali // 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; } @@ -1132,16 +1234,17 @@ Bool_t AliTRDCalibraFit::AnalyseLinearFitters(AliTRDCalibraVdriftLinearFit *cali //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 ++; @@ -1163,7 +1266,7 @@ Bool_t AliTRDCalibraFit::AnalyseLinearFitters(AliTRDCalibraVdriftLinearFit *cali } // 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)); @@ -1172,6 +1275,208 @@ Bool_t AliTRDCalibraFit::AnalyseLinearFitters(AliTRDCalibraVdriftLinearFit *cali 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)); + } + delete fDebugStreamer; + 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___________________ //_____________________________________________________________________________ @@ -1429,9 +1734,9 @@ void AliTRDCalibraFit::RemoveOutliers(Int_t type, Bool_t perdetector){ // Initialisation //////////////////////// Double_t meanAll = 0.0; - Double_t rmsAll = 0.0; - Int_t countAll = 0; - //////////// + Double_t rmsAll = 0.0; + Int_t countAll = 0; + //////////// // compute //////////// for (Int_t k = 0; k < loop; k++) { @@ -1561,7 +1866,10 @@ void AliTRDCalibraFit::RemoveOutliers2(Bool_t perdetector){ 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))) coef[(Int_t)(col*rowMax+row)] = 100.0; + 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 } @@ -1823,11 +2131,11 @@ void AliTRDCalibraFit::PutMeanValueOtherVectorFit2(Int_t ofwhat, Bool_t perdetec for (Int_t col = 0; col < colMax; col++) { value = coef[(Int_t)(col*rowMax+row)]; if(value > 70.0) { - if((ofwhat == 0) && (meanAll > -1.5) && (countAll > 15)) 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] > -1.5) && (countDetector[detector] > 20)) coef[(Int_t)(col*rowMax+row)] = meanDetector[detector]+100.0; - else if((meanSupermodule[sector] > -1.5) && (countSupermodule[sector] > 15)) coef[(Int_t)(col*rowMax+row)] = meanSupermodule[sector]+100.0; - else if((meanAll > -1.5) && (countAll > 15)) 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 @@ -1910,6 +2218,7 @@ AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectGain(const TObjArray *vectorFit, // 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!"); @@ -1972,6 +2281,7 @@ AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectT0(const TObjArray *vectorFit, Bo 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; // @@ -2030,7 +2340,49 @@ AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectLorentzAngle(const TObjArray *vec 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; @@ -2526,6 +2878,7 @@ Bool_t AliTRDCalibraFit::InitFitCH() 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; @@ -2570,11 +2923,12 @@ Bool_t AliTRDCalibraFit::InitFitPH() 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; @@ -2622,6 +2976,7 @@ Bool_t AliTRDCalibraFit::InitFitPRF() 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; @@ -2643,6 +2998,8 @@ Bool_t AliTRDCalibraFit::InitFitLinearFitter() 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++) { @@ -2650,45 +3007,27 @@ Bool_t AliTRDCalibraFit::InitFitLinearFitter() 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) { @@ -2866,8 +3205,7 @@ Bool_t AliTRDCalibraFit::NotEnoughStatisticCH(Int_t idect) 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++){ @@ -2948,8 +3286,7 @@ Bool_t AliTRDCalibraFit::NotEnoughStatisticCH(Int_t idect) } 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); @@ -2992,8 +3329,7 @@ Bool_t AliTRDCalibraFit::NotEnoughStatisticPH(Int_t idect,Double_t nentries) 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++){ @@ -3089,8 +3425,7 @@ Bool_t AliTRDCalibraFit::NotEnoughStatisticPH(Int_t idect,Double_t nentries) } 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(); @@ -3140,8 +3475,7 @@ Bool_t AliTRDCalibraFit::NotEnoughStatisticPRF(Int_t idect) 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++){ @@ -3224,8 +3558,7 @@ Bool_t AliTRDCalibraFit::NotEnoughStatisticPRF(Int_t idect) } 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(); @@ -3273,14 +3606,14 @@ Bool_t AliTRDCalibraFit::NotEnoughStatisticLinearFitter() 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(); @@ -3288,6 +3621,33 @@ Bool_t AliTRDCalibraFit::NotEnoughStatisticLinearFitter() 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) { @@ -3300,8 +3660,7 @@ 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++){ @@ -3412,8 +3771,7 @@ Bool_t AliTRDCalibraFit::FillInfosFitPH(Int_t idect,Double_t nentries) 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++){ @@ -3542,8 +3900,7 @@ Bool_t AliTRDCalibraFit::FillInfosFitPRF(Int_t idect) 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++){ @@ -3664,6 +4021,28 @@ Bool_t AliTRDCalibraFit::FillInfosFitLinearFitter() 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) @@ -3673,7 +4052,7 @@ 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++) { @@ -3722,7 +4101,7 @@ void AliTRDCalibraFit::FillFillPH(Int_t idect,Double_t nentries) // // End of one detector - if ((idect == (fCount-1))) { + if (idect == (fCount-1)) { FillVectorFit(); FillVectorFit2(); // Reset @@ -3783,7 +4162,7 @@ void AliTRDCalibraFit::FillFillPRF(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++) { @@ -3828,21 +4207,79 @@ void AliTRDCalibraFit::FillFillPRF(Int_t idect) } //________________________________________________________________________________ -void AliTRDCalibraFit::FillFillLinearFitter() +void AliTRDCalibraFit::FillFillLinearFitter() +{ + // + // DebugStream and fVectorFit + // + + // End of one detector + FillVectorFit(); + FillVectorFit2(); + + + // Reset + for (Int_t k = 0; k < 2304; k++) { + fCurrentCoefDetector[k] = 0.0; + fCurrentCoefDetector2[k] = 0.0; + } + + + if(fDebugLevel > 1){ + + if ( !fDebugStreamer ) { + //debug stream + TDirectory *backup = gDirectory; + fDebugStreamer = new TTreeSRedirector("TRDDebugFitLinearFitter.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 = fCurrentCoef[0]; + Float_t vs = fCurrentCoef[1]; + Float_t vfE = fCurrentCoefE; + Float_t lorentzangler = fCurrentCoef2[0]; + Float_t elorentzangler = fCurrentCoefE2; + Float_t lorentzangles = fCurrentCoef2[1]; + + (* fDebugStreamer) << "FillFillLinearFitter"<< + "detector="<cd(); //we don't want to be cd'd to the debug streamer } @@ -3863,12 +4300,8 @@ void AliTRDCalibraFit::FillFillLinearFitter() Int_t detector = fCountDet; Int_t stack = GetStack(fCountDet); Int_t layer = GetLayer(fCountDet); - Float_t vf = fCurrentCoef[0]; - Float_t vs = fCurrentCoef[1]; - Float_t vfE = fCurrentCoefE; - Float_t lorentzangler = fCurrentCoef2[0]; - Float_t elorentzangler = fCurrentCoefE2; - Float_t lorentzangles = fCurrentCoef2[1]; + Float_t vf = fCurrentCoef2[0]; + Float_t vfE = fCurrentCoefE2; (* fDebugStreamer) << "FillFillLinearFitter"<< "detector="<GetValue(fCountDet); - fCurrentCoef2[1] = fCalDet2->GetValue(fCountDet); + fCurrentCoef[1] = fCalDetVdriftUsed->GetValue(fCountDet); + fCurrentCoef2[1] = fCalDetExBUsed->GetValue(fCountDet); return kTRUE; } @@ -4266,7 +4695,7 @@ void AliTRDCalibraFit::FitPente(TH1* projPH) 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; } } @@ -4358,7 +4787,9 @@ void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH) 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))); @@ -4386,7 +4817,7 @@ void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH) y[1] = pentea->GetBinContent(binmax-1); y[2] = pentea->GetBinContent(binmax); CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4); - AliInfo("At the limit for beginning!"); + //AliInfo("At the limit for beginning!"); break; case 2: minnn = pentea->GetBinCenter(binmax-2); @@ -4469,7 +4900,9 @@ void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH) fPhd[0] = placemaximum; } + ///////////////////////////// // Amplification region + ///////////////////////////// binmax = 0; ju = 0; for (Int_t kbin = 1; kbin < projPH->GetNbinsX(); kbin ++) { @@ -4591,8 +5024,11 @@ void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH) } 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))); @@ -4603,21 +5039,21 @@ void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH) //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(TMath::Abs(pente->GetBinContent(binmin+1)) <= 0.0000000000001){ - AliInfo("No entries for the next bin!"); + //AliInfo("No entries for the next bin!"); pente->SetBinContent(binmin,0); if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin(); } @@ -4661,14 +5097,14 @@ void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH) 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)) && @@ -4742,7 +5178,7 @@ void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH) //richtung + if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1))) { //AliInfo("polynome 2+ false"); - put = kFALSE; + putd = kFALSE; } } //pol2 case 2 @@ -4780,26 +5216,26 @@ void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH) //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(c0,c1,c2,c3,c4); //AliInfo(Form("GetMinimum of the function %f",polynome->GetMinimumX())); @@ -4835,13 +5271,15 @@ void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH) if ((fPhd[2] > fPhd[0]) && (fPhd[2] > fPhd[1]) && (fPhd[1] > fPhd[0]) && - (put)) { + (put) && + (putd)) { fCurrentCoef[0] = (kDrWidth) / (fPhd[2]-fPhd[1]); if(fCurrentCoef[0] > 2.5) fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]); else fNumberFitSuccess++; if (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; } } @@ -4850,16 +5288,17 @@ void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH) } } else { - //printf("Put default %f\n",-TMath::Abs(fCurrentCoef[1])); + ////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; @@ -4875,6 +5314,7 @@ void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH) 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()); @@ -4882,7 +5322,8 @@ void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH) 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(); @@ -4905,7 +5346,7 @@ void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH) //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); } @@ -5430,7 +5871,15 @@ void AliTRDCalibraFit::FitMeanW(TH1 *projch, Double_t nentries) Int_t sumCurrent = 0; for(Int_t k = 0; k 0.95) break; + if (fractionGetBinContent(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); @@ -5444,6 +5893,8 @@ void AliTRDCalibraFit::FitMeanW(TH1 *projch, Double_t nentries) 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); @@ -5497,58 +5948,132 @@ void AliTRDCalibraFit::FitMeanWSm(TH1 *projch, Float_t sumAll) 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 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 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); @@ -5565,25 +6090,40 @@ void AliTRDCalibraFit::FitCH(TH1 *projch, Double_t mean) } //_____________________________________________________________________________ -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 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(); @@ -5592,15 +6132,13 @@ void AliTRDCalibraFit::FitBisCH(TH1* projch, Double_t mean) // 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; @@ -5611,9 +6149,9 @@ void AliTRDCalibraFit::FitBisCH(TH1* projch, Double_t mean) ,&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))) { @@ -5633,6 +6171,110 @@ void AliTRDCalibraFit::FitBisCH(TH1* projch, Double_t mean) 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 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; @@ -6025,14 +6667,77 @@ Double_t AliTRDCalibraFit::LanGauFun(const Double_t *x, const Double_t *par) 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]; @@ -6047,16 +6752,17 @@ Double_t AliTRDCalibraFit::LanGauFun(const Double_t *x, const Double_t *par) 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; } //_____________________________________________________________________________ @@ -6085,7 +6791,7 @@ TF1 *AliTRDCalibraFit::LanGauFit(TH1 *his, const Double_t *fitrange, const Doubl 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++) { @@ -6097,115 +6803,45 @@ TF1 *AliTRDCalibraFit::LanGauFit(TH1 *his, const Double_t *fitrange, const Doubl 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; -}