X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;ds=sidebyside;f=TRD%2FAliTRDCalibraFit.cxx;h=a56103a85385ffef8f0ede9661f3c94c9f0dd761;hb=86a261ed28339b6b10815ba383313b0ec6bc300c;hp=9773ffb98741d87ffcdd996a091aae7a10f535a1;hpb=4aad967c887383f2cea5fa479136f4ca8df9e016;p=u%2Fmrichter%2FAliRoot.git diff --git a/TRD/AliTRDCalibraFit.cxx b/TRD/AliTRDCalibraFit.cxx index 9773ffb9874..a56103a8538 100644 --- a/TRD/AliTRDCalibraFit.cxx +++ b/TRD/AliTRDCalibraFit.cxx @@ -28,13 +28,13 @@ // previous calibration procedure. // The function SetDebug enables the user to see: // _fDebug = 0: nothing, only the values are written in the tree if wanted -// _fDebug = 1: a comparaison of the coefficients found and the default values +// _fDebug = 1: only the fit of the choosen calibration group fFitVoir (SetFitVoir) +// _fDebug = 2: a comparaison of the coefficients found and the default values // in the choosen database. // fCoef , histogram of the coefs as function of the calibration group number // fDelta , histogram of the relative difference of the coef with the default // value in the database as function of the calibration group number // fError , dirstribution of this relative difference -// _fDebug = 2: only the fit of the choosen calibration group fFitVoir (SetFitVoir) // _fDebug = 3: The coefficients in the choosen detector fDet (SetDet) as function of the // pad row and col number // _fDebug = 4; The coeffcicients in the choosen detector fDet (SetDet) like in the 3 but with @@ -59,11 +59,12 @@ #include #include #include -#include #include #include #include -#include +#include +#include +#include #include "AliLog.h" #include "AliMathBase.h" @@ -72,10 +73,12 @@ #include "AliTRDCalibraMode.h" #include "AliTRDCalibraVector.h" #include "AliTRDCalibraVdriftLinearFit.h" +#include "AliTRDCalibraExbAltFit.h" #include "AliTRDcalibDB.h" #include "AliTRDgeometry.h" #include "AliTRDpadPlane.h" #include "AliTRDgeometry.h" +#include "AliTRDCommonParam.h" #include "./Cal/AliTRDCalROC.h" #include "./Cal/AliTRDCalPad.h" #include "./Cal/AliTRDCalDet.h" @@ -104,7 +107,6 @@ AliTRDCalibraFit *AliTRDCalibraFit::Instance() return fgInstance; } - //______________________________________________________________________________________ void AliTRDCalibraFit::Terminate() { @@ -121,7 +123,6 @@ void AliTRDCalibraFit::Terminate() } } - //______________________________________________________________________________________ AliTRDCalibraFit::AliTRDCalibraFit() :TObject() @@ -129,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) @@ -137,6 +140,7 @@ AliTRDCalibraFit::AliTRDCalibraFit() ,fAccCDB(kFALSE) ,fMinEntries(800) ,fRebin(1) + ,fScaleGain(0.02431) ,fNumberFit(0) ,fNumberFitSuccess(0) ,fNumberEnt(0) @@ -154,10 +158,13 @@ AliTRDCalibraFit::AliTRDCalibraFit() ,fEntriesCurrent(0) ,fCountDet(0) ,fCount(0) + ,fNbDet(0) ,fCalDet(0x0) ,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) @@ -212,10 +222,13 @@ AliTRDCalibraFit::AliTRDCalibraFit(const AliTRDCalibraFit &c) ,fEntriesCurrent(c.fEntriesCurrent) ,fCountDet(c.fCountDet) ,fCount(c.fCount) +,fNbDet(c.fNbDet) ,fCalDet(0x0) ,fCalROC(0x0) ,fCalDet2(0x0) ,fCalROC2(0x0) +,fCalDetVdriftUsed(0x0) +,fCalDetExBUsed(0x0) ,fCurrentCoefDetector(0x0) ,fCurrentCoefDetector2(0x0) ,fVectorFit(0) @@ -238,10 +251,13 @@ AliTRDCalibraFit::AliTRDCalibraFit(const AliTRDCalibraFit &c) } if(c.fCalDet) fCalDet = new AliTRDCalDet(*c.fCalDet); if(c.fCalDet2) fCalDet2 = new AliTRDCalDet(*c.fCalDet2); - + 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(); @@ -296,7 +312,11 @@ AliTRDCalibraFit::~AliTRDCalibraFit() if ( fCalDet ) delete fCalDet; if ( fCalDet2 ) delete fCalDet2; if ( fCalROC ) delete fCalROC; - if ( fCalROC2 ) delete fCalROC2; + if ( fCalROC2 ) delete fCalROC2; + if ( fCalDetVdriftUsed) delete fCalDetVdriftUsed; + if ( fCalDetExBUsed) delete fCalDetExBUsed; + if( fCurrentCoefDetector ) delete [] fCurrentCoefDetector; + if( fCurrentCoefDetector2 ) delete [] fCurrentCoefDetector2; fVectorFit.Delete(); fVectorFit2.Delete(); if (fGeo) { @@ -317,30 +337,19 @@ void AliTRDCalibraFit::Destroy() } } -//__________________________________________________________________________________ -void AliTRDCalibraFit::RangeChargeIntegration(Float_t vdrift, Float_t t0, Int_t &begin, Int_t &peak, Int_t &end) +//_____________________________________________________________________________ +void AliTRDCalibraFit::DestroyDebugStreamer() { // - // From the drift velocity and t0 - // return the position of the peak and maximum negative slope + // Delete DebugStreamer // - - 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); + if ( fDebugStreamer ) delete fDebugStreamer; + fDebugStreamer = 0x0; + } //____________Functions fit Online CH2d________________________________________ -Bool_t AliTRDCalibraFit::AnalyseCH(TH2I *ch) +Bool_t AliTRDCalibraFit::AnalyseCH(const TH2I *ch) { // // Fit the 1D histos, projections of the 2D ch on the Xaxis, for each @@ -348,8 +357,9 @@ Bool_t AliTRDCalibraFit::AnalyseCH(TH2I *ch) // // Set the calibration mode - const char *name = ch->GetTitle(); - SetModeCalibration(name,0); + //const char *name = ch->GetTitle(); + TString name = ch->GetTitle(); + if(!SetModeCalibration(name,0)) return kFALSE; // Number of Ybins (detectors or groups of pads) Int_t nbins = ch->GetNbinsX();// charge @@ -409,8 +419,10 @@ Bool_t AliTRDCalibraFit::AnalyseCH(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 @@ -426,7 +438,7 @@ Bool_t AliTRDCalibraFit::AnalyseCH(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 { @@ -446,8 +458,9 @@ Bool_t AliTRDCalibraFit::AnalyseCH(AliTRDCalibraVector *calvect) // // Set the calibraMode - const char *name = calvect->GetNameCH(); - SetModeCalibration(name,0); + //const char *name = calvect->GetNameCH(); + TString name = calvect->GetNameCH(); + if(!SetModeCalibration(name,0)) return kFALSE; // Number of Xbins (detectors or groups of pads) if (!InitFit((432*calvect->GetDetCha0(0)+108*calvect->GetDetCha2(0)),0)) { @@ -470,34 +483,31 @@ Bool_t AliTRDCalibraFit::AnalyseCH(AliTRDCalibraVector *calvect) // Take the histo Double_t nentries = 0.0; Double_t mean = 0.0; - TH1F *projch = 0x0; - Bool_t something = kTRUE; - if(!calvect->GetCHEntries(fCountDet)) something = kFALSE; - if(something){ - TString tname("CH"); - tname += idect; - projch = calvect->ConvertVectorCHHisto(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0)))),(const char *) tname); - projch->SetDirectory(0); - for (Int_t k = 0; k < calvect->GetNumberBinCharge(); k++) { - nentries += projch->GetBinContent(k+1); - mean += projch->GetBinCenter(k+1)*projch->GetBinContent(k+1); - } - if (nentries > 0) { - fNumberEnt++; - mean /= nentries; - } - //printf("The number of entries for the group %d is %f\n",idect,nentries); - // Rebin - if (fRebin > 1) { - projch = ReBin((TH1F *) projch); - } + if(!calvect->GetCHEntries(fCountDet)) { + NotEnoughStatisticCH(idect); + continue; + } + + TString tname("CH"); + tname += idect; + TH1F *projch = calvect->ConvertVectorCHHisto(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0)))),(const char *) tname); + projch->SetDirectory(0); + for (Int_t k = 0; k < calvect->GetNumberBinCharge(); k++) { + nentries += projch->GetBinContent(k+1); + mean += projch->GetBinCenter(k+1)*projch->GetBinContent(k+1); + } + if (nentries > 0) { + fNumberEnt++; + mean /= nentries; + } + //printf("The number of entries for the group %d is %f\n",idect,nentries); + // Rebin + if (fRebin > 1) { + projch = ReBin((TH1F *) projch); } // This detector has not enough statistics or was not found in VectorCH if (nentries <= fMinEntries) { NotEnoughStatisticCH(idect); - if (fDebugLevel != 1) { - if(projch) delete projch; - } continue; } // Statistic of the histos fitted @@ -508,16 +518,14 @@ 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 FillInfosFitCH(idect); - // Memory!!! - if (fDebugLevel != 1) { - delete projch; - } } // Boucle object // Normierungcharge if (fDebugLevel != 1) { @@ -525,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 { @@ -535,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(TProfile2D *ph) +Double_t AliTRDCalibraFit::AnalysePHAllTogether(const TProfile2D *ph) { // // Take the 1D profiles (average pulse height), projections of the 2D PH @@ -545,14 +602,67 @@ Bool_t AliTRDCalibraFit::AnalysePH(TProfile2D *ph) // A first calibration of T0 is also made using the same method // - // Set the calibration mode - const char *name = ph->GetTitle(); - SetModeCalibration(name,1); - // 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()) { @@ -566,34 +676,30 @@ Bool_t AliTRDCalibraFit::AnalysePH(TProfile2D *ph) 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); - 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(); @@ -606,15 +712,12 @@ Bool_t AliTRDCalibraFit::AnalysePH(TProfile2D *ph) default: return kFALSE; } // Fill the tree if end of a detector or only the pointer to the branch!!! - FillInfosFitPH(idect); - // 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 { @@ -624,66 +727,82 @@ Bool_t AliTRDCalibraFit::AnalysePH(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(); - SetModeCalibration(name,1); + //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 - TH1F *projph = 0x0; - fEntriesCurrent = 0; - Bool_t something = kTRUE; - if(!calvect->GetPHEntries(fCountDet)) something = kFALSE; - if(something){ - TString tname("PH"); - tname += idect; - projph = CorrectTheError((TGraphErrors *) (calvect->ConvertVectorPHTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) tname))); - projph->SetDirectory(0); + 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); } - //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); + if (nentries <= fMinEntries) { + //printf("Not enough statistic!\n"); + NotEnoughStatisticPH(idect,nentries); if (fDebugLevel != 1) { - if(projph) delete projph; + 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; @@ -692,16 +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); + 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 { @@ -712,7 +830,7 @@ Bool_t AliTRDCalibraFit::AnalysePH(AliTRDCalibraVector *calvect) return kTRUE; } //____________Functions fit Online PRF2d_______________________________________ -Bool_t AliTRDCalibraFit::AnalysePRF(TProfile2D *prf) +Bool_t AliTRDCalibraFit::AnalysePRF(const TProfile2D *prf) { // // Take the 1D profiles (pad response function), projections of the 2D PRF @@ -721,8 +839,9 @@ Bool_t AliTRDCalibraFit::AnalysePRF(TProfile2D *prf) // // Set the calibration mode - const char *name = prf->GetTitle(); - SetModeCalibration(name,2); + //const char *name = prf->GetTitle(); + TString name = prf->GetTitle(); + if(!SetModeCalibration(name,2)) return kFALSE; // Number of Ybins (detectors or groups of pads) Int_t nybins = prf->GetNbinsY();// calibration groups @@ -787,8 +906,8 @@ Bool_t AliTRDCalibraFit::AnalysePRF(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 { @@ -799,7 +918,7 @@ Bool_t AliTRDCalibraFit::AnalysePRF(TProfile2D *prf) return kTRUE; } //____________Functions fit Online PRF2d_______________________________________ -Bool_t AliTRDCalibraFit::AnalysePRFMarianFit(TProfile2D *prf) +Bool_t AliTRDCalibraFit::AnalysePRFMarianFit(const TProfile2D *prf) { // // Take the 1D profiles (pad response function), projections of the 2D PRF @@ -808,8 +927,9 @@ Bool_t AliTRDCalibraFit::AnalysePRFMarianFit(TProfile2D *prf) // // Set the calibration mode - const char *name = prf->GetTitle(); - SetModeCalibration(name,2); + //const char *name = prf->GetTitle(); + TString name = prf->GetTitle(); + if(!SetModeCalibration(name,2)) return kFALSE; // Number of Ybins (detectors or groups of pads) TAxis *xprf = prf->GetXaxis(); @@ -882,8 +1002,8 @@ Bool_t AliTRDCalibraFit::AnalysePRFMarianFit(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 { @@ -903,8 +1023,9 @@ Bool_t AliTRDCalibraFit::AnalysePRF(AliTRDCalibraVector *calvect) // // Set the calibra mode - const char *name = calvect->GetNamePRF(); - SetModeCalibration(name,2); + //const char *name = calvect->GetNamePRF(); + TString name = calvect->GetNamePRF(); + if(!SetModeCalibration(name,2)) return kFALSE; //printf("test0 %s\n",name); // Number of Xbins (detectors or groups of pads) @@ -928,22 +1049,19 @@ Bool_t AliTRDCalibraFit::AnalysePRF(AliTRDCalibraVector *calvect) UpdatefCountDetAndfCount(idect,2); ReconstructFitRowMinRowMax(idect,2); // Take the histo - TH1F *projprf = 0x0; fEntriesCurrent = 0; - Bool_t something = kTRUE; - if(!calvect->GetPRFEntries(fCountDet)) something = kFALSE; - if(something){ - TString tname("PRF"); - tname += idect; - projprf = CorrectTheError((TGraphErrors *) (calvect->ConvertVectorPRFTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) tname))); - projprf->SetDirectory(0); + if(!calvect->GetPRFEntries(fCountDet)) { + NotEnoughStatisticPRF(idect); + continue; } + TString tname("PRF"); + tname += idect; + TH1F *projprf = calvect->CorrectTheError((TGraphErrors *) (calvect->ConvertVectorPRFTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) tname)),fEntriesCurrent); + projprf->SetDirectory(0); + if(fEntriesCurrent > 0) fNumberEnt++; // This detector has not enough statistics or was off if (fEntriesCurrent <= fMinEntries) { NotEnoughStatisticPRF(idect); - if (fDebugLevel != 1) { - if(projprf) delete projprf; - } continue; } // Statistic of the histos fitted @@ -960,14 +1078,10 @@ Bool_t AliTRDCalibraFit::AnalysePRF(AliTRDCalibraVector *calvect) } // Fill the tree if end of a detector or only the pointer to the branch!!! FillInfosFitPRF(idect); - // Memory!!! - if (fDebugLevel != 1) { - delete projprf; - } } // 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)); @@ -986,11 +1100,12 @@ Bool_t AliTRDCalibraFit::AnalysePRFMarianFit(AliTRDCalibraVector *calvect) // // Set the calibra mode - const char *name = calvect->GetNamePRF(); - SetModeCalibration(name,2); + //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); - printf("test1 %d\n",nbg); + //printf("test1 %d\n",nbg); if(nbg == -1) return kFALSE; if(nbg > 0) fMethod = 1; else fMethod = 0; @@ -1023,35 +1138,33 @@ Bool_t AliTRDCalibraFit::AnalysePRFMarianFit(AliTRDCalibraVector *calvect) UpdatefCountDetAndfCount(idect,2); ReconstructFitRowMinRowMax(idect,2); // Take the histo - TGraphErrors *projprftree = 0x0; fEntriesCurrent = 0; - Bool_t something = kTRUE; - if(!calvect->GetPRFEntries(fCountDet)) something = kFALSE; - if(something){ - TString tname("PRF"); - tname += idect; - projprftree = calvect->ConvertVectorPRFTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) tname); - nbins = projprftree->GetN(); - arrayx = (Double_t *)projprftree->GetX(); - arraye = (Double_t *)projprftree->GetEX(); - arraym = (Double_t *)projprftree->GetY(); - arrayme = (Double_t *)projprftree->GetEY(); - Float_t step = arrayx[1]-arrayx[0]; - lowedge = arrayx[0] - step/2.0; - upedge = arrayx[(nbins-1)] + step/2.0; - //printf("nbins est %d\n",nbins); - for(Int_t k = 0; k < nbins; k++){ - fEntriesCurrent += (Int_t)arraye[k]; - //printf("for %d we have %f, %f\n",k,arraye[k],((projprftree->GetEX())[k])); - if(arraye[k]>0.0) arrayme[k] = TMath::Sqrt(TMath::Abs(arrayme[k]-arraym[k]*arraym[k])/arraye[k]); - } - if(fEntriesCurrent > 0) fNumberEnt++; + if(!calvect->GetPRFEntries(fCountDet)) { + NotEnoughStatisticPRF(idect); + continue; + } + TString tname("PRF"); + tname += idect; + TGraphErrors *projprftree = calvect->ConvertVectorPRFTGraphErrors(fCountDet,(idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1)))),(const char *) tname); + nbins = projprftree->GetN(); + arrayx = (Double_t *)projprftree->GetX(); + arraye = (Double_t *)projprftree->GetEX(); + arraym = (Double_t *)projprftree->GetY(); + arrayme = (Double_t *)projprftree->GetEY(); + Float_t step = arrayx[1]-arrayx[0]; + lowedge = arrayx[0] - step/2.0; + upedge = arrayx[(nbins-1)] + step/2.0; + //printf("nbins est %d\n",nbins); + for(Int_t k = 0; k < nbins; k++){ + fEntriesCurrent += (Int_t)arraye[k]; + //printf("for %d we have %f, %f\n",k,arraye[k],((projprftree->GetEX())[k])); + if(arraye[k]>0.0) arrayme[k] = TMath::Sqrt(TMath::Abs(arrayme[k]-arraym[k]*arraym[k])/arraye[k]); } + 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 (fEntriesCurrent <= fMinEntries) { NotEnoughStatisticPRF(idect); - if(projprftree) delete projprftree; continue; } // Statistic of the histos fitted @@ -1068,14 +1181,10 @@ Bool_t AliTRDCalibraFit::AnalysePRFMarianFit(AliTRDCalibraVector *calvect) } // Fill the tree if end of a detector or only the pointer to the branch!!! FillInfosFitPRF(idect); - // Memory!!! - if (fDebugLevel != 1) { - delete projprftree; - } } // 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)); @@ -1106,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; } @@ -1125,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 ++; @@ -1145,7 +1255,7 @@ Bool_t AliTRDCalibraFit::AnalyseLinearFitters(AliTRDCalibraVdriftLinearFit *cali 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]; } @@ -1156,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)); @@ -1166,55 +1276,257 @@ Bool_t AliTRDCalibraFit::AnalyseLinearFitters(AliTRDCalibraVdriftLinearFit *cali return kTRUE; } -//____________Functions for seeing if the pad is really okey___________________ -//_____________________________________________________________________________ -Int_t AliTRDCalibraFit::GetNumberOfGroupsPRF(const char* nametitle) +//______________________________________________________________________________________ +Bool_t AliTRDCalibraFit::AnalyseExbAltFit(AliTRDCalibraExbAltFit *calivdli) { // - // Get numberofgroupsprf + // The linear method // + + fStatisticMean = 0.0; + fNumberFit = 0; + fNumberFitSuccess = 0; + fNumberEnt = 0; + if(!InitFitExbAlt()) return kFALSE; + - // Some patterns - const Char_t *pattern0 = "Ngp0"; - const Char_t *pattern1 = "Ngp1"; - const Char_t *pattern2 = "Ngp2"; - const Char_t *pattern3 = "Ngp3"; - const Char_t *pattern4 = "Ngp4"; - const Char_t *pattern5 = "Ngp5"; - const Char_t *pattern6 = "Ngp6"; + for(Int_t idet = 0; idet < 540; idet++){ - // Nrphi mode - if (strstr(nametitle,pattern0)) { - return 0; - } - if (strstr(nametitle,pattern1)) { - return 1; - } - if (strstr(nametitle,pattern2)) { - return 2; - } - if (strstr(nametitle,pattern3)) { - return 3; - } - if (strstr(nametitle,pattern4)) { - return 4; + + //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(); + } - if (strstr(nametitle,pattern5)) { - return 5; + // 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)); } - if (strstr(nametitle,pattern6)){ - return 6; + else { + AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt)); } - else return -1; - - + delete fDebugStreamer; + fDebugStreamer = 0x0; + return kTRUE; + } -//_____________________________________________________________________________ -Bool_t AliTRDCalibraFit::SetModeCalibration(const char *name, Int_t i) +//____________Functions fit Online CH2d________________________________________ +void AliTRDCalibraFit::AnalyseLinearFittersAllTogether(AliTRDCalibraVdriftLinearFit *calivdli, Double_t &vdriftoverall, Double_t &exboverall) { // - // Set fNz[i] and fNrphi[i] of the AliTRDCalibraFit::Instance() - // corresponding to the given name + // 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(TString nametitle) +{ + // + // Get numberofgroupsprf + // + + // Some patterns + const Char_t *pattern0 = "Ngp0"; + const Char_t *pattern1 = "Ngp1"; + const Char_t *pattern2 = "Ngp2"; + const Char_t *pattern3 = "Ngp3"; + const Char_t *pattern4 = "Ngp4"; + const Char_t *pattern5 = "Ngp5"; + const Char_t *pattern6 = "Ngp6"; + + // Nrphi mode + if (strstr(nametitle.Data(),pattern0)) { + return 0; + } + if (strstr(nametitle.Data(),pattern1)) { + return 1; + } + if (strstr(nametitle.Data(),pattern2)) { + return 2; + } + if (strstr(nametitle.Data(),pattern3)) { + return 3; + } + if (strstr(nametitle.Data(),pattern4)) { + return 4; + } + if (strstr(nametitle.Data(),pattern5)) { + return 5; + } + if (strstr(nametitle.Data(),pattern6)){ + return 6; + } + else return -1; + + +} +//_____________________________________________________________________________ +Bool_t AliTRDCalibraFit::SetModeCalibration(TString name, Int_t i) +{ + // + // Set fNz[i] and fNrphi[i] of the AliTRDCalibraFit::Instance() + // corresponding to the given name // if(!SetNzFromTObject(name,i)) return kFALSE; @@ -1224,7 +1536,7 @@ Bool_t AliTRDCalibraFit::SetModeCalibration(const char *name, Int_t i) } //_____________________________________________________________________________ -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() @@ -1240,42 +1552,89 @@ Bool_t AliTRDCalibraFit::SetNrphiFromTObject(const char *name, Int_t i) const Char_t *patternrphi5 = "Nrphi5"; const Char_t *patternrphi6 = "Nrphi6"; + + const Char_t *patternrphi10 = "Nrphi10"; + const Char_t *patternrphi100 = "Nrphi100"; + const Char_t *patternz10 = "Nz10"; + const Char_t *patternz100 = "Nz100"; + // Nrphi mode - if (strstr(name,patternrphi0)) { + if ((strstr(name.Data(),patternrphi100)) && (strstr(name.Data(),patternz100))) { + fCalibraMode->SetAllTogether(i); + fNbDet = 540; + if (fDebugLevel > 1) { + AliInfo(Form("fNbDet %d and 100",fNbDet)); + } + return kTRUE; + } + if ((strstr(name.Data(),patternrphi10)) && (strstr(name.Data(),patternz10))) { + fCalibraMode->SetPerSuperModule(i); + fNbDet = 30; + if (fDebugLevel > 1) { + AliInfo(Form("fNDet %d and 100",fNbDet)); + } + return kTRUE; + } + + 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)); + } return kTRUE; } + if (fDebugLevel > 1) { + AliInfo(Form("fNbDet %d and rest",fNbDet)); + } fCalibraMode->SetNrphi(i ,0); return kFALSE; - + } //_____________________________________________________________________________ -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() @@ -1288,33 +1647,525 @@ Bool_t AliTRDCalibraFit::SetNzFromTObject(const char *name, Int_t i) const Char_t *patternz2 = "Nz2"; const Char_t *patternz3 = "Nz3"; const Char_t *patternz4 = "Nz4"; - - if (strstr(name,patternz0)) { + + const Char_t *patternrphi10 = "Nrphi10"; + const Char_t *patternrphi100 = "Nrphi100"; + const Char_t *patternz10 = "Nz10"; + const Char_t *patternz100 = "Nz100"; + + if ((strstr(name.Data(),patternrphi100)) && (strstr(name.Data(),patternz100))) { + fCalibraMode->SetAllTogether(i); + fNbDet = 540; + if (fDebugLevel > 1) { + AliInfo(Form("fNbDet %d and 100",fNbDet)); + } + return kTRUE; + } + if ((strstr(name.Data(),patternrphi10)) && (strstr(name.Data(),patternz10))) { + fCalibraMode->SetPerSuperModule(i); + fNbDet = 30; + if (fDebugLevel > 1) { + AliInfo(Form("fNbDet %d and 10",fNbDet)); + } + return kTRUE; + } + 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)) { + if (strstr(name.Data(),patternz3)) { fCalibraMode->SetNz(i ,3); + if (fDebugLevel > 1) { + AliInfo(Form("fNbDet %d and 3",fNbDet)); + } return kTRUE; } - if (strstr(name,patternz4)) { + 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++; + } + } + 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 + } + } + 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 + } +} +//______________________________________________________________________ +void AliTRDCalibraFit::PutMeanValueOtherVectorFit(Int_t ofwhat, Bool_t perdetector){ + // + // ofwhat is equaled to 0: mean value of all passing detectors + // ofwhat is equaled to 1: mean value of the detector, otherwise supermodule, otherwise all + // + + 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 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++) { + detector = ((AliTRDFitInfo *) fVectorFit.At(k))->GetDetector(); + sector = GetSector(detector); + 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++; + } + } + 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) { + rmsDetector[detector] += value*value; + meanDetector[detector] += value; + countDetector[detector]++; + rmsSupermodule[sector] += value*value; + meanSupermodule[sector] += value; + countSupermodule[sector]++; + rmsAll += value*value; + meanAll += value; + countAll++; + } + + } // Col + } // Row + } + } + 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]; + 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]; + 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++) { + 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(); + + 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) { + 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) && (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 + if(fDebugLevel > 1){ + + if ( !fDebugStreamer ) { + //debug stream + TDirectory *backup = gDirectory; + fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root"); + if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer + } + + Float_t coefnow = coef[(Int_t)(col*rowMax+row)]; + + (* fDebugStreamer) << "PutMeanValueOtherVectorFit"<< + "detector="<GetDetector(); + sector = GetSector(detector); + 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++; + } + } + 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) { + rmsDetector[detector] += value*value; + meanDetector[detector] += value; + countDetector[detector]++; + rmsSupermodule[sector] += value*value; + meanSupermodule[sector] += value; + countSupermodule[sector]++; + rmsAll += value*value; + meanAll += value; + countAll++; + } + + } // Col + } // Row + } + } + 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]; + 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]; + 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++) { + 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(); + + 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) { + if((ofwhat == 0) && (meanAll > -3.0) && (countAll > 15)) coef[(Int_t)(col*rowMax+row)] = meanAll+100.0; + if(ofwhat == 1){ + 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 + if(fDebugLevel > 1){ + + if ( !fDebugStreamer ) { + //debug stream + TDirectory *backup = gDirectory; + fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root"); + if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer + } + + Float_t coefnow = coef[(Int_t)(col*rowMax+row)]; + + (* fDebugStreamer) << "PutMeanValueOtherVectorFit2"<< + "detector="<At(k))->GetDetector(); Float_t mean = 0.0; @@ -1355,7 +2207,7 @@ AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectVdrift(TObjArray *vectorFit, Bool return object; } //_____________________________________________________________________________ -AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectGain(TObjArray *vectorFit, Double_t scaleFitFactor, Bool_t perdetector) +AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectGain(const TObjArray *vectorFit, Bool_t meanOtherBefore, Double_t scaleFitFactor, Bool_t perdetector) { // // It creates the AliTRDCalDet object from the AliTRDFitInfo @@ -1366,6 +2218,7 @@ AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectGain(TObjArray *vectorFit, Double // 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!"); @@ -1377,7 +2230,10 @@ AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectGain(TObjArray *vectorFit, Double Float_t mean = 0.0; if(perdetector){ value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0]; - if(value > 0) value = value*scaleFitFactor; + if(!meanOtherBefore){ + if(value > 0) value = value*scaleFitFactor; + } + else value = value*scaleFitFactor; mean = TMath::Abs(value); } else{ @@ -1387,20 +2243,24 @@ AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectGain(TObjArray *vectorFit, Double for (Int_t row = 0; row < rowMax; row++) { for (Int_t col = 0; col < colMax; col++) { value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)]; - if(value > 0) value = value*scaleFitFactor; + if(!meanOtherBefore) { + if(value > 0) value = value*scaleFitFactor; + } + else value = value*scaleFitFactor; mean += TMath::Abs(value); count++; } // Col } // Row if(count > 0) mean = mean/count; } + if(mean < 0.1) mean = 0.1; object->SetValue(detector,mean); } return object; } //_____________________________________________________________________________ -AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectT0(TObjArray *vectorFit, Bool_t perdetector) +AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectT0(const TObjArray *vectorFit, Bool_t perdetector) { // // It creates the AliTRDCalDet object from the AliTRDFitInfo2 @@ -1420,7 +2280,12 @@ AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectT0(TObjArray *vectorFit, Bool_t p detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector(); Float_t min = 100.0; if(perdetector){ - min = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0]; + 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; + // + min = value; } else{ Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector)); @@ -1428,6 +2293,9 @@ AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectT0(TObjArray *vectorFit, Bool_t p for (Int_t row = 0; row < rowMax; row++) { for (Int_t col = 0; col < colMax; col++) { value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)]; + // check successful + if(value > 70.0) value = value-100.0; + // if(min > value) min = value; } // Col } // Row @@ -1439,7 +2307,7 @@ AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectT0(TObjArray *vectorFit, Bool_t p } //_____________________________________________________________________________ -AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectLorentzAngle(TObjArray *vectorFit) +AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectLorentzAngle(const TObjArray *vectorFit) { // // It creates the AliTRDCalDet object from the AliTRDFitInfo2 @@ -1472,14 +2340,56 @@ AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectLorentzAngle(TObjArray *vectorFit 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; } //_____________________________________________________________________________ -TObject *AliTRDCalibraFit::CreatePadObjectGain(TObjArray *vectorFit, Double_t scaleFitFactor, AliTRDCalDet *detobject) +AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectExbAlt(const TObjArray *vectorFit) +{ + // + // It creates the AliTRDCalDet object from the AliTRDFitInfo2 + // It takes the min value of the coefficients per detector + // This object has to be written in the database + // + + // Create the DetObject + AliTRDCalDet *object = new AliTRDCalDet("tan(lorentzangle)","tan(lorentzangle) (detector value)"); + + + Int_t loop = (Int_t) vectorFit->GetEntriesFast(); + if(loop != 540) AliInfo("The Vector Fit is not complete!"); + Int_t detector = -1; + Float_t value = 0.0; + + for (Int_t k = 0; k < loop; k++) { + detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector(); + /* + Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector)); + Int_t colMax = fGeo->GetColMax(GetLayer(detector)); + Float_t min = 100.0; + for (Int_t row = 0; row < rowMax; row++) { + for (Int_t col = 0; col < colMax; col++) { + value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[(Int_t)(col*rowMax+row)]; + mean += -TMath::Abs(value); + count++; + } // Col + } // Row + if(count > 0) mean = mean/count; + */ + value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0]; + //if(value > 70.0) value = value-100.0; + object->SetValue(detector,value); + } + + return object; + +} +//_____________________________________________________________________________ +TObject *AliTRDCalibraFit::CreatePadObjectGain(const TObjArray *vectorFit, Double_t scaleFitFactor, const AliTRDCalDet *detobject) { // // It Creates the AliTRDCalPad object from AliTRDFitInfo @@ -1511,7 +2421,7 @@ TObject *AliTRDCalibraFit::CreatePadObjectGain(TObjArray *vectorFit, Double_t sc 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++) { @@ -1527,7 +2437,7 @@ TObject *AliTRDCalibraFit::CreatePadObjectGain(TObjArray *vectorFit, Double_t sc return object; } //_____________________________________________________________________________ -TObject *AliTRDCalibraFit::CreatePadObjectVdrift(TObjArray *vectorFit, AliTRDCalDet *detobject) +TObject *AliTRDCalibraFit::CreatePadObjectVdrift(const TObjArray *vectorFit, const AliTRDCalDet *detobject) { // // It Creates the AliTRDCalPad object from AliTRDFitInfo @@ -1574,7 +2484,7 @@ TObject *AliTRDCalibraFit::CreatePadObjectVdrift(TObjArray *vectorFit, AliTRDCal } //_____________________________________________________________________________ -TObject *AliTRDCalibraFit::CreatePadObjectT0(TObjArray *vectorFit, AliTRDCalDet *detobject) +TObject *AliTRDCalibraFit::CreatePadObjectT0(const TObjArray *vectorFit, const AliTRDCalDet *detobject) { // // It Creates the AliTRDCalPad object from AliTRDFitInfo2 @@ -1611,6 +2521,9 @@ TObject *AliTRDCalibraFit::CreatePadObjectT0(TObjArray *vectorFit, AliTRDCalDet for (Int_t row = 0; row < rowMax; row++) { for (Int_t col = 0; col < colMax; col++) { value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[(Int_t)(col*rowMax+row)]; + // check successful + if(value > 70.0) value = value - 100.0; + // calROC->SetValue(col,row,value-min); } // Col } // Row @@ -1620,7 +2533,7 @@ TObject *AliTRDCalibraFit::CreatePadObjectT0(TObjArray *vectorFit, AliTRDCalDet } //_____________________________________________________________________________ -TObject *AliTRDCalibraFit::CreatePadObjectPRF(TObjArray *vectorFit) +TObject *AliTRDCalibraFit::CreatePadObjectPRF(const TObjArray *vectorFit) { // // It Creates the AliTRDCalPad object from AliTRDFitInfo @@ -1652,7 +2565,7 @@ TObject *AliTRDCalibraFit::CreatePadObjectPRF(TObjArray *vectorFit) } //_____________________________________________________________________________ -AliTRDCalDet *AliTRDCalibraFit::MakeOutliersStatDet(TObjArray *vectorFit, const char *name, Double_t &mean) +AliTRDCalDet *AliTRDCalibraFit::MakeOutliersStatDet(const TObjArray *vectorFit, const char *name, Double_t &mean) { // // It Creates the AliTRDCalDet object from AliTRDFitInfo @@ -1690,7 +2603,7 @@ AliTRDCalDet *AliTRDCalibraFit::MakeOutliersStatDet(TObjArray *vectorFit, const return object; } //_____________________________________________________________________________ -TObject *AliTRDCalibraFit::MakeOutliersStatPad(TObjArray *vectorFit, const char *name, Double_t &mean) +TObject *AliTRDCalibraFit::MakeOutliersStatPad(const TObjArray *vectorFit, const char *name, Double_t &mean) { // // It Creates the AliTRDCalPad object from AliTRDFitInfo @@ -1935,7 +2848,7 @@ Bool_t AliTRDCalibraFit::InitFit(Int_t nbins, Int_t i) // Quick verification that we have the good pad calibration mode! if (fNumberOfBinsExpected != nbins) { - AliInfo("It doesn't correspond to the mode of pad group calibration!"); + AliInfo(Form("It doesn't correspond to the mode of pad group calibration: expected %d and seen %d!",fNumberOfBinsExpected,nbins)); return kFALSE; } @@ -1965,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; @@ -2009,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; @@ -2061,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; @@ -2082,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++) { @@ -2089,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,cal->GetOmegaTau(fCalDet->GetValue(k),-fMagneticField)); - } - //printf("test3\n"); - } - else{ - Float_t devalue = 1.5; - Float_t devalue2 = cal->GetOmegaTau(1.5,-fMagneticField); - 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) { @@ -2152,10 +3052,16 @@ 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 @@ -2173,6 +3079,17 @@ void AliTRDCalibraFit::CalculNumberOfBinsExpected(Int_t i) // fNumberOfBinsExpected = 0; + // All + if((fCalibraMode->GetNz(i) == 100) && (fCalibraMode->GetNrphi(i) == 100)){ + fNumberOfBinsExpected = 1; + return; + } + // Per supermodule + if((fCalibraMode->GetNz(i) == 10) && (fCalibraMode->GetNrphi(i) == 10)){ + fNumberOfBinsExpected = 18; + return; + } + // More fCalibraMode->ModePadCalibration(2,i); fCalibraMode->ModePadFragmentation(0,2,0,i); fCalibraMode->SetDetChamb2(i); @@ -2245,18 +3162,19 @@ void AliTRDCalibraFit::UpdatefCountDetAndfCount(Int_t idect, Int_t i) // Doesn't matter for 2 // if (fCount == idect) { - // On en est au detector - fCountDet += 1; - // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi - fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),i); - fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet) + // On en est au detector (or first detector in the group) + fCountDet += 1; + AliDebug(2,Form("We are at the detector %d\n",fCountDet)); + // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi + fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),i); + fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet) ,(Int_t) GetStack(fCountDet) ,(Int_t) GetSector(fCountDet),i); - // Set for the next detector - fCount += fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i); - // calib objects - SetCalROC(i); - } + // Set for the next detector + fCount += fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i); + // calib objects + SetCalROC(i); + } } //____________Functions for initialising the AliTRDCalibraFit in the code_________ void AliTRDCalibraFit::ReconstructFitRowMinRowMax(Int_t idect, Int_t i) @@ -2264,10 +3182,13 @@ void AliTRDCalibraFit::ReconstructFitRowMinRowMax(Int_t idect, Int_t i) // // Reconstruct the min pad row, max pad row, min pad col and // max pad col of the calibration group for the Fit functions + // idect is the calibration group inside the detector // if (fDebugLevel != 1) { fCalibraMode->ReconstructionRowPadGroup((Int_t) (idect-(fCount-(fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i)))),i); } + AliDebug(2,Form("AliTRDCalibraFit::ReconstructFitRowMinRowMax: the local calibration group is %d",idect-(fCount-(fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i))))); + AliDebug(2,Form("AliTRDCalibraFit::ReconstructFitRowMinRowMax: the number of group per detector is %d",fCalibraMode->GetNfragZ(i)*fCalibraMode->GetNfragRphi(i))); } //____________Functions for initialising the AliTRDCalibraFit in the code_________ Bool_t AliTRDCalibraFit::NotEnoughStatisticCH(Int_t idect) @@ -2281,11 +3202,91 @@ Bool_t AliTRDCalibraFit::NotEnoughStatisticCH(Int_t idect) if (fDebugLevel == 1) { AliInfo("The element has not enough statistic to be fitted"); } - + 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)); + // loop over detectors + for(Int_t det = firstdetector; det < lastdetector; det++){ + + //Set the calibration object again + fCountDet = det; + SetCalROC(0); + + // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi + // Put them at 1 + fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),0); + fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet) + ,(Int_t) GetStack(fCountDet) + ,(Int_t) GetSector(fCountDet),0); + // Reconstruct row min row max + ReconstructFitRowMinRowMax(idect,0); + + // Calcul the coef from the database choosen for the detector + CalculChargeCoefMean(kFALSE); + + //stack 2, not stack 2 + Int_t factor = 0; + if(GetStack(fCountDet) == 2) factor = 12; + else factor = 16; + + // Fill the fCurrentCoefDetector with negative value to say: not fitted + for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) { + for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) { + fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]); + } + } + + //Put default value negative + fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]); + fCurrentCoefE = 0.0; + + // Fill the stuff + FillVectorFit(); + // Debug + if(fDebugLevel > 1){ + + if ( !fDebugStreamer ) { + //debug stream + TDirectory *backup = gDirectory; + fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root"); + if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer + } + + Int_t detector = fCountDet; + Int_t caligroup = idect; + Short_t rowmin = fCalibraMode->GetRowMin(0); + Short_t rowmax = fCalibraMode->GetRowMax(0); + Short_t colmin = fCalibraMode->GetColMin(0); + Short_t colmax = fCalibraMode->GetColMax(0); + Float_t gf = fCurrentCoef[0]; + Float_t gfs = fCurrentCoef[1]; + Float_t gfE = fCurrentCoefE; + + (*fDebugStreamer) << "FillFillCH" << + "detector=" << detector << + "caligroup=" << caligroup << + "rowmin=" << rowmin << + "rowmax=" << rowmax << + "colmin=" << colmin << + "colmax=" << colmax << + "gf=" << gf << + "gfs=" << gfs << + "gfE=" << gfE << + "\n"; + + } + // Reset + for (Int_t k = 0; k < 2304; k++) { + fCurrentCoefDetector[k] = 0.0; + } + + }// loop detector + AliDebug(2,Form("Check the count now: fCountDet %d",fCountDet)); + } 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); @@ -2314,7 +3315,7 @@ Bool_t AliTRDCalibraFit::NotEnoughStatisticCH(Int_t idect) //____________Functions for initialising the AliTRDCalibraFit in the code_________ -Bool_t AliTRDCalibraFit::NotEnoughStatisticPH(Int_t idect) +Bool_t AliTRDCalibraFit::NotEnoughStatisticPH(Int_t idect,Double_t nentries) { // // For the case where there are not enough entries in the histograms @@ -2324,10 +3325,107 @@ Bool_t AliTRDCalibraFit::NotEnoughStatisticPH(Int_t idect) if (fDebugLevel == 1) { AliInfo("The element has not enough statistic to be fitted"); } + 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)); + // loop over detectors + for(Int_t det = firstdetector; det < lastdetector; det++){ + + //Set the calibration object again + fCountDet = det; + SetCalROC(1); + + // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi + // Put them at 1 + fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),1); + fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet) + ,(Int_t) GetStack(fCountDet) + ,(Int_t) GetSector(fCountDet),1); + // Reconstruct row min row max + ReconstructFitRowMinRowMax(idect,1); + + // Calcul the coef from the database choosen for the detector + CalculVdriftCoefMean(); + CalculT0CoefMean(); + + //stack 2, not stack 2 + Int_t factor = 0; + if(GetStack(fCountDet) == 2) factor = 12; + else factor = 16; + + // Fill the fCurrentCoefDetector with negative value to say: not fitted + for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) { + for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) { + fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]); + fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[1] + 100.0; + } + } + + //Put default value negative + fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]); + fCurrentCoefE = 0.0; + fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0; + fCurrentCoefE2 = 0.0; + + // Fill the stuff + FillVectorFit(); + FillVectorFit2(); + // Debug + if(fDebugLevel > 1){ + + if ( !fDebugStreamer ) { + //debug stream + TDirectory *backup = gDirectory; + fDebugStreamer = new TTreeSRedirector("TRDDebugFitPH.root"); + if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer + } + + + Int_t detector = fCountDet; + Int_t caligroup = idect; + Short_t rowmin = fCalibraMode->GetRowMin(1); + Short_t rowmax = fCalibraMode->GetRowMax(1); + Short_t colmin = fCalibraMode->GetColMin(1); + Short_t colmax = fCalibraMode->GetColMax(1); + Float_t vf = fCurrentCoef[0]; + Float_t vs = fCurrentCoef[1]; + Float_t vfE = fCurrentCoefE; + Float_t t0f = fCurrentCoef2[0]; + Float_t t0s = fCurrentCoef2[1]; + Float_t t0E = fCurrentCoefE2; + + + + (* fDebugStreamer) << "FillFillPH"<< + "detector="<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(); @@ -2342,22 +3440,22 @@ Bool_t AliTRDCalibraFit::NotEnoughStatisticPH(Int_t idect) for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) { for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) { fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]); - fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[1]; + fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[1] + 100.0; } } // Put the default value fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]); fCurrentCoefE = 0.0; - fCurrentCoef2[0] = fCurrentCoef2[1]; + fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0; fCurrentCoefE2 = 0.0; - FillFillPH(idect); + FillFillPH(idect,nentries); } return kTRUE; - + } @@ -2373,10 +3471,94 @@ Bool_t AliTRDCalibraFit::NotEnoughStatisticPRF(Int_t idect) if (fDebugLevel == 1) { AliInfo("The element has not enough statistic to be fitted"); } + 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)); + + // loop over detectors + for(Int_t det = firstdetector; det < lastdetector; det++){ + + //Set the calibration object again + fCountDet = det; + SetCalROC(2); + + // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi + // Put them at 1 + fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),2); + fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet) + ,(Int_t) GetStack(fCountDet) + ,(Int_t) GetSector(fCountDet),2); + // Reconstruct row min row max + ReconstructFitRowMinRowMax(idect,2); + + // Calcul the coef from the database choosen for the detector + CalculPRFCoefMean(); + + //stack 2, not stack 2 + Int_t factor = 0; + if(GetStack(fCountDet) == 2) factor = 12; + else factor = 16; + + // Fill the fCurrentCoefDetector with negative value to say: not fitted + for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) { + for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) { + fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]); + } + } + + //Put default value negative + fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]); + fCurrentCoefE = 0.0; + + // Fill the stuff + FillVectorFit(); + // Debug + if(fDebugLevel > 1){ + + if ( !fDebugStreamer ) { + //debug stream + TDirectory *backup = gDirectory; + fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root"); + if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer + } + + Int_t detector = fCountDet; + Int_t layer = GetLayer(fCountDet); + Int_t caligroup = idect; + Short_t rowmin = fCalibraMode->GetRowMin(2); + Short_t rowmax = fCalibraMode->GetRowMax(2); + Short_t colmin = fCalibraMode->GetColMin(2); + Short_t colmax = fCalibraMode->GetColMax(2); + Float_t widf = fCurrentCoef[0]; + Float_t wids = fCurrentCoef[1]; + Float_t widfE = fCurrentCoefE; + + (* fDebugStreamer) << "FillFillPRF"<< + "detector="<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(); @@ -2389,12 +3571,12 @@ Bool_t AliTRDCalibraFit::NotEnoughStatisticPRF(Int_t idect) // Fill the fCurrentCoefDetector for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) { for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) { - fCurrentCoefDetector[(Int_t)(j*factor+k)] = -fCurrentCoef[1]; + fCurrentCoefDetector[(Int_t)(j*factor+k)] = -TMath::Abs(fCurrentCoef[1]); } } // Put the default value - fCurrentCoef[0] = -fCurrentCoef[1]; + fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]); fCurrentCoefE = 0.0; FillFillPRF(idect); @@ -2424,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(); @@ -2439,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) { @@ -2448,26 +3657,109 @@ Bool_t AliTRDCalibraFit::FillInfosFitCH(Int_t idect) // if (fDebugLevel != 1) { - - Int_t factor = 0; - if(GetStack(fCountDet) == 2) factor = 12; - else factor = 16; - - for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) { - for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) { - fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0]; + 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)); + // loop over detectors + for(Int_t det = firstdetector; det < lastdetector; det++){ + + //Set the calibration object again + fCountDet = det; + SetCalROC(0); + + // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi + // Put them at 1 + fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),0); + fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet) + ,(Int_t) GetStack(fCountDet) + ,(Int_t) GetSector(fCountDet),0); + // Reconstruct row min row max + ReconstructFitRowMinRowMax(idect,0); + + // Calcul the coef from the database choosen for the detector + if(fCurrentCoef[0] < 0.0) CalculChargeCoefMean(kFALSE); + else CalculChargeCoefMean(kTRUE); + + //stack 2, not stack 2 + Int_t factor = 0; + if(GetStack(fCountDet) == 2) factor = 12; + else factor = 16; + + // Fill the fCurrentCoefDetector with negative value to say: not fitted + Double_t coeftoput = 1.0; + if(fCurrentCoef[0] < 0.0) coeftoput = - TMath::Abs(fCurrentCoef[1]); + else coeftoput = fCurrentCoef[0]; + for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) { + for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) { + fCurrentCoefDetector[(Int_t)(j*factor+k)] = coeftoput; + } + } + + // Fill the stuff + FillVectorFit(); + // Debug + if(fDebugLevel > 1){ + + if ( !fDebugStreamer ) { + //debug stream + TDirectory *backup = gDirectory; + fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root"); + if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer + } + + Int_t detector = fCountDet; + Int_t caligroup = idect; + Short_t rowmin = fCalibraMode->GetRowMin(0); + Short_t rowmax = fCalibraMode->GetRowMax(0); + Short_t colmin = fCalibraMode->GetColMin(0); + Short_t colmax = fCalibraMode->GetColMax(0); + Float_t gf = fCurrentCoef[0]; + Float_t gfs = fCurrentCoef[1]; + Float_t gfE = fCurrentCoefE; + + (*fDebugStreamer) << "FillFillCH" << + "detector=" << detector << + "caligroup=" << caligroup << + "rowmin=" << rowmin << + "rowmax=" << rowmax << + "colmin=" << colmin << + "colmax=" << colmax << + "gf=" << gf << + "gfs=" << gfs << + "gfE=" << gfE << + "\n"; + + } + // Reset + for (Int_t k = 0; k < 2304; k++) { + fCurrentCoefDetector[k] = 0.0; + } + + }// loop detector + //printf("Check the count now: fCountDet %d\n",fCountDet); + } + else{ + + Int_t factor = 0; + if(GetStack(fCountDet) == 2) factor = 12; + else factor = 16; + + for (Int_t k = fCalibraMode->GetRowMin(0); k < fCalibraMode->GetRowMax(0); k++) { + for (Int_t j = fCalibraMode->GetColMin(0); j < fCalibraMode->GetColMax(0); j++) { + fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0]; + } } + + FillFillCH(idect); } - - FillFillCH(idect); - } return kTRUE; } //____________Functions for initialising the AliTRDCalibraFit in the code_________ -Bool_t AliTRDCalibraFit::FillInfosFitPH(Int_t idect) +Bool_t AliTRDCalibraFit::FillInfosFitPH(Int_t idect,Double_t nentries) { // // Fill the coefficients found with the fits or other @@ -2475,18 +3767,123 @@ Bool_t AliTRDCalibraFit::FillInfosFitPH(Int_t idect) // if (fDebugLevel != 1) { - - Int_t factor = 0; - if(GetStack(fCountDet) == 2) factor = 12; - else factor = 16; - - for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) { - for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) { - fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0]; - fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[0]; - } - } - FillFillPH(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)); + + // loop over detectors + for(Int_t det = firstdetector; det < lastdetector; det++){ + + //Set the calibration object again + fCountDet = det; + SetCalROC(1); + + // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi + // Put them at 1 + fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),1); + fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet) + ,(Int_t) GetStack(fCountDet) + ,(Int_t) GetSector(fCountDet),1); + // Reconstruct row min row max + ReconstructFitRowMinRowMax(idect,1); + + // Calcul the coef from the database choosen for the detector + CalculVdriftCoefMean(); + CalculT0CoefMean(); + + //stack 2, not stack 2 + Int_t factor = 0; + if(GetStack(fCountDet) == 2) factor = 12; + else factor = 16; + + // Fill the fCurrentCoefDetector with negative value to say: not fitted + Double_t coeftoput = 1.5; + Double_t coeftoput2 = 0.0; + + if(fCurrentCoef[0] < 0.0) coeftoput = - TMath::Abs(fCurrentCoef[1]); + else coeftoput = fCurrentCoef[0]; + + if(fCurrentCoef2[0] > 70.0) coeftoput2 = fCurrentCoef2[1] + 100.0; + else coeftoput2 = fCurrentCoef2[0]; + + for (Int_t k = fCalibraMode->GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) { + for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) { + fCurrentCoefDetector[(Int_t)(j*factor+k)] = coeftoput; + fCurrentCoefDetector2[(Int_t)(j*factor+k)] = coeftoput2; + } + } + + // Fill the stuff + FillVectorFit(); + FillVectorFit2(); + // Debug + if(fDebugLevel > 1){ + + if ( !fDebugStreamer ) { + //debug stream + TDirectory *backup = gDirectory; + fDebugStreamer = new TTreeSRedirector("TRDDebugFitPH.root"); + if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer + } + + + Int_t detector = fCountDet; + Int_t caligroup = idect; + Short_t rowmin = fCalibraMode->GetRowMin(1); + Short_t rowmax = fCalibraMode->GetRowMax(1); + Short_t colmin = fCalibraMode->GetColMin(1); + Short_t colmax = fCalibraMode->GetColMax(1); + Float_t vf = fCurrentCoef[0]; + Float_t vs = fCurrentCoef[1]; + Float_t vfE = fCurrentCoefE; + Float_t t0f = fCurrentCoef2[0]; + Float_t t0s = fCurrentCoef2[1]; + Float_t t0E = fCurrentCoefE2; + + + + (* fDebugStreamer) << "FillFillPH"<< + "detector="<GetRowMin(1); k < fCalibraMode->GetRowMax(1); k++) { + for (Int_t j = fCalibraMode->GetColMin(1); j < fCalibraMode->GetColMax(1); j++) { + fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0]; + fCurrentCoefDetector2[(Int_t)(j*factor+k)] = fCurrentCoef2[0]; + } + } + + FillFillPH(idect,nentries); + } } return kTRUE; } @@ -2499,20 +3896,106 @@ Bool_t AliTRDCalibraFit::FillInfosFitPRF(Int_t idect) // if (fDebugLevel != 1) { - - Int_t factor = 0; - if(GetStack(fCountDet) == 2) factor = 12; - else factor = 16; + if (fNbDet > 0){ - // Pointer to the branch - for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) { - for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) { - fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[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)); + + // loop over detectors + for(Int_t det = firstdetector; det < lastdetector; det++){ + + //Set the calibration object again + fCountDet = det; + SetCalROC(2); + + // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi + // Put them at 1 + fCalibraMode->ModePadCalibration((Int_t) GetStack(fCountDet),2); + fCalibraMode->ModePadFragmentation((Int_t) GetLayer(fCountDet) + ,(Int_t) GetStack(fCountDet) + ,(Int_t) GetSector(fCountDet),2); + // Reconstruct row min row max + ReconstructFitRowMinRowMax(idect,2); + + // Calcul the coef from the database choosen for the detector + CalculPRFCoefMean(); + + //stack 2, not stack 2 + Int_t factor = 0; + if(GetStack(fCountDet) == 2) factor = 12; + else factor = 16; + + // Fill the fCurrentCoefDetector with negative value to say: not fitted + Double_t coeftoput = 1.0; + if(fCurrentCoef[0] < 0.0) coeftoput = - TMath::Abs(fCurrentCoef[1]); + else coeftoput = fCurrentCoef[0]; + for (Int_t k = fCalibraMode->GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) { + for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) { + fCurrentCoefDetector[(Int_t)(j*factor+k)] = coeftoput; + } + } + + // Fill the stuff + FillVectorFit(); + // Debug + if(fDebugLevel > 1){ + + if ( !fDebugStreamer ) { + //debug stream + TDirectory *backup = gDirectory; + fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root"); + if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer + } + + Int_t detector = fCountDet; + Int_t layer = GetLayer(fCountDet); + Int_t caligroup = idect; + Short_t rowmin = fCalibraMode->GetRowMin(2); + Short_t rowmax = fCalibraMode->GetRowMax(2); + Short_t colmin = fCalibraMode->GetColMin(2); + Short_t colmax = fCalibraMode->GetColMax(2); + Float_t widf = fCurrentCoef[0]; + Float_t wids = fCurrentCoef[1]; + Float_t widfE = fCurrentCoefE; + + (* fDebugStreamer) << "FillFillPRF"<< + "detector="<GetRowMin(2); k < fCalibraMode->GetRowMax(2); k++) { + for (Int_t j = fCalibraMode->GetColMin(2); j < fCalibraMode->GetColMax(2); j++) { + fCurrentCoefDetector[(Int_t)(j*factor+k)] = fCurrentCoef[0]; + } } + FillFillPRF(idect); } - FillFillPRF(idect); } - + return kTRUE; } @@ -2538,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) @@ -2547,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++) { @@ -2589,14 +4094,14 @@ void AliTRDCalibraFit::FillFillCH(Int_t idect) } } //________________________________________________________________________________ -void AliTRDCalibraFit::FillFillPH(Int_t idect) +void AliTRDCalibraFit::FillFillPH(Int_t idect,Double_t nentries) { // // DebugStream and fVectorFit and fVectorFit2 // // End of one detector - if ((idect == (fCount-1))) { + if (idect == (fCount-1)) { FillVectorFit(); FillVectorFit2(); // Reset @@ -2633,6 +4138,7 @@ void AliTRDCalibraFit::FillFillPH(Int_t idect) (* fDebugStreamer) << "FillFillPH"<< "detector="< 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="<GetNz(1) > 0) || - (fCalibraMode->GetNrphi(1) > 0)) { + if(((fCalibraMode->GetNz(1) > 0) || + (fCalibraMode->GetNrphi(1) > 0)) && ((fCalibraMode->GetNz(1) != 10) && (fCalibraMode->GetNz(1) != 100))) { + for (Int_t row = fCalibraMode->GetRowMin(1); row < fCalibraMode->GetRowMax(1); row++) { for (Int_t col = fCalibraMode->GetColMin(1); col < fCalibraMode->GetColMax(1); col++) { fCurrentCoef2[1] += (Float_t) (fCalROC2->GetValue(col,row)+fCalDet2->GetValue(fCountDet)); } } + fCurrentCoef2[1] = fCurrentCoef2[1] / ((fCalibraMode->GetColMax(1)-fCalibraMode->GetColMin(1))*(fCalibraMode->GetRowMax(1)-fCalibraMode->GetRowMin(1))); + } else { + if(!fAccCDB){ fCurrentCoef2[1] = fCalDet2->GetValue(fCountDet); } else{ + for(Int_t row = 0; row < fGeo->GetRowMax(GetLayer(fCountDet),GetStack(fCountDet),GetSector(fCountDet)); row++){ for(Int_t col = 0; col < fGeo->GetColMax(GetLayer(fCountDet)); col++){ fCurrentCoef2[1] += (Float_t) (fCalROC2->GetValue(col,row)+fCalDet2->GetValue(fCountDet)); } } fCurrentCoef2[1] = fCurrentCoef2[1] / ((fGeo->GetRowMax(GetLayer(fCountDet),GetStack(fCountDet),GetSector(fCountDet)))*(fGeo->GetColMax(GetLayer(fCountDet)))); + } } } @@ -2809,8 +4371,8 @@ Bool_t AliTRDCalibraFit::CalculChargeCoefMean(Bool_t vrai) fCurrentCoef[1] = 0.0; if(fDebugLevel != 1){ - if ((fCalibraMode->GetNz(0) > 0) || - (fCalibraMode->GetNrphi(0) > 0)) { + if (((fCalibraMode->GetNz(0) > 0) || + (fCalibraMode->GetNrphi(0) > 0)) && ((fCalibraMode->GetNz(0) != 10) && (fCalibraMode->GetNz(0) != 100))) { for (Int_t row = fCalibraMode->GetRowMin(0); row < fCalibraMode->GetRowMax(0); row++) { for (Int_t col = fCalibraMode->GetColMin(0); col < fCalibraMode->GetColMax(0); col++) { fCurrentCoef[1] += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet)); @@ -2856,14 +4418,17 @@ Bool_t AliTRDCalibraFit::CalculVdriftCoefMean() fCurrentCoef[1] = 0.0; if(fDebugLevel != 1){ - if ((fCalibraMode->GetNz(1) > 0) || - (fCalibraMode->GetNrphi(1) > 0)) { + if (((fCalibraMode->GetNz(1) > 0) || + (fCalibraMode->GetNrphi(1) > 0)) && ((fCalibraMode->GetNz(1) != 10) && (fCalibraMode->GetNz(1) != 100))) { + for (Int_t row = fCalibraMode->GetRowMin(1); row < fCalibraMode->GetRowMax(1); row++) { for (Int_t col = fCalibraMode->GetColMin(1); col < fCalibraMode->GetColMax(1); col++) { fCurrentCoef[1] += (Float_t) (fCalROC->GetValue(col,row)*fCalDet->GetValue(fCountDet)); } } + fCurrentCoef[1] = fCurrentCoef[1] / ((fCalibraMode->GetColMax(1)-fCalibraMode->GetColMin(1))*(fCalibraMode->GetRowMax(1)-fCalibraMode->GetRowMin(1))); + } else { //per detectors @@ -2879,8 +4444,8 @@ Bool_t AliTRDCalibraFit::CalculVdriftLorentzCoef() // 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; } @@ -2930,18 +4495,27 @@ void AliTRDCalibraFit::SetCalROC(Int_t i) switch (i) { case 0: - if(fCalROC) delete fCalROC; - fCalROC = new AliTRDCalROC(*(cal->GetGainFactorROC(fCountDet))); + if( fCalROC ){ + fCalROC->~AliTRDCalROC(); + new(fCalROC) AliTRDCalROC(*(cal->GetGainFactorROC(fCountDet))); + }else fCalROC = new AliTRDCalROC(*(cal->GetGainFactorROC(fCountDet))); break; case 1: - if(fCalROC) delete fCalROC; - if(fCalROC2) delete fCalROC2; - fCalROC = new AliTRDCalROC(*(cal->GetVdriftROC(fCountDet))); - fCalROC2 = new AliTRDCalROC(*(cal->GetT0ROC(fCountDet))); + if( fCalROC ){ + fCalROC->~AliTRDCalROC(); + new(fCalROC) AliTRDCalROC(*(cal->GetVdriftROC(fCountDet))); + }else fCalROC = new AliTRDCalROC(*(cal->GetVdriftROC(fCountDet))); + if( fCalROC2 ){ + fCalROC2->~AliTRDCalROC(); + new(fCalROC2) AliTRDCalROC(*(cal->GetT0ROC(fCountDet))); + }else fCalROC2 = new AliTRDCalROC(*(cal->GetT0ROC(fCountDet))); break; case 2: - if(fCalROC) delete fCalROC; - fCalROC = new AliTRDCalROC(*(cal->GetPRFROC(fCountDet))); break; + if( fCalROC ){ + fCalROC->~AliTRDCalROC(); + new(fCalROC) AliTRDCalROC(*(cal->GetPRFROC(fCountDet))); + }else fCalROC = new AliTRDCalROC(*(cal->GetPRFROC(fCountDet))); + break; default: return; } } @@ -3030,11 +4604,11 @@ void AliTRDCalibraFit::FitPente(TH1* projPH) 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]; } } @@ -3061,10 +4635,10 @@ void AliTRDCalibraFit::FitPente(TH1* projPH) 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){ @@ -3095,10 +4669,10 @@ void AliTRDCalibraFit::FitPente(TH1* projPH) 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; @@ -3121,18 +4695,18 @@ void AliTRDCalibraFit::FitPente(TH1* projPH) if (fPhdt0 >= 0.0) { fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins; - if (fCurrentCoef2[0] < -1.0) { - fCurrentCoef2[0] = fCurrentCoef2[1]; + if (fCurrentCoef2[0] < -3.0) { + fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0; } } else { - fCurrentCoef2[0] = fCurrentCoef2[1]; + fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0; } } else { fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]); - fCurrentCoef2[0] = fCurrentCoef2[1]; + fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0; } if (fDebugLevel == 1) { @@ -3165,13 +4739,15 @@ void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH) // // 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; @@ -3194,7 +4770,11 @@ void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH) 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(); @@ -3207,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))); @@ -3234,8 +4816,8 @@ void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH) 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); @@ -3248,7 +4830,7 @@ void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH) 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){ @@ -3264,7 +4846,7 @@ void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH) 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); @@ -3277,7 +4859,7 @@ void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH) 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); @@ -3292,7 +4874,7 @@ void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH) 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; @@ -3301,7 +4883,7 @@ void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH) 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; @@ -3318,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 ++) { @@ -3357,7 +4941,7 @@ void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH) 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: @@ -3371,7 +4955,7 @@ void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH) 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) @@ -3388,7 +4972,7 @@ void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH) 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); @@ -3401,7 +4985,7 @@ void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH) 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); @@ -3416,7 +5000,7 @@ void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH) 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; @@ -3424,7 +5008,7 @@ void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH) 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; @@ -3440,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))); @@ -3452,13 +5039,25 @@ 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))) put = kFALSE; - if((projPH->GetBinContent(binmin)+projPH->GetBinError(binmin)) > (projPH->GetBinContent(binmin-1))) put = kFALSE; + if((projPH->GetBinContent(binmin)-projPH->GetBinError(binmin)) < (projPH->GetBinContent(binmin+1))) { + //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!"); + putd = kFALSE; + } + 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(); + } + x[0] = 0.0; x[1] = 0.0; @@ -3493,18 +5092,32 @@ void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH) 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))) put = kFALSE; + (pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) { + //AliInfo("polynome 4 false 1"); + 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))) put = kFALSE; + (pente->GetBinContent(binmin-3) <= pente->GetBinContent(binmin-2))) { + //AliInfo("polynome 4 false 2"); + putd = kFALSE; + } + // poly 3 if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1)) && - (pente->GetBinContent(binmin-2) > pente->GetBinContent(binmin-1))) case1 = kTRUE; + (pente->GetBinContent(binmin-2) > pente->GetBinContent(binmin-1))) { + //AliInfo("polynome 4 case 1"); + case1 = kTRUE; + } if((pente->GetBinContent(binmin+2) > pente->GetBinContent(binmin+1)) && - (pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) case4 = kTRUE; + (pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) { + //AliInfo("polynome 4 case 4"); + case4 = kTRUE; + } + } //case binmin = nbins-2 //pol3 case 1 @@ -3521,10 +5134,13 @@ void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH) 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))) case2 = kTRUE; + if((pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) { + //AliInfo("polynome 3- case 2"); + case2 = kTRUE; + } } //pol3 case 4 if(((binmin <= (nbins-3)) && ((binmin-1) == TMath::Min(binmax+4, projPH->GetNbinsX()))) || @@ -3540,9 +5156,12 @@ void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH) 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))) case2 = kTRUE; + if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1))) { + //AliInfo("polynome 3+ case 2"); + case2 = kTRUE; + } } //pol2 case 5 if((binmin <= (nbins-3)) && (binmin == TMath::Min(binmax+4, projPH->GetNbinsX()))){ @@ -3555,9 +5174,12 @@ void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH) 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))) put = kFALSE; + if((pente->GetBinContent(binmin+2) <= pente->GetBinContent(binmin+1))) { + //AliInfo("polynome 2+ false"); + putd = kFALSE; + } } //pol2 case 2 if(((binmin == (nbins-2)) && ((binmin-1) == TMath::Min(binmax+4, projPH->GetNbinsX()))) || @@ -3571,7 +5193,7 @@ void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH) 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 } @@ -3587,32 +5209,35 @@ void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH) 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))) put = kFALSE; + if((pente->GetBinContent(binmin-2) <= pente->GetBinContent(binmin-1))) { + //AliInfo("polynome 2- false "); + 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+4, projPH->GetNbinsX()))){ - put = kFALSE; - AliInfo("For the drift...problem!"); + if((binmin <= (nbins-3)) && (binmin < TMath::Min(binmax+6, projPH->GetNbinsX()))){ + 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; @@ -3628,6 +5253,9 @@ void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH) } fPhd[2] = placeminimum; } + //printf("La fin %d\n",((Int_t)(fPhd[2]*10.0))+2); + if((((Int_t)(fPhd[2]*10.0))+2) >= projPH->GetNbinsX()) fPhd[2] = 0.0; + if(((((Int_t)(fPhd[2]*10.0))+2) < projPH->GetNbinsX()) && (projPH->GetBinContent(((Int_t)(fPhd[2]*10.0))+2)==0)) fPhd[2] = 0.0; Float_t fPhdt0 = 0.0; Float_t t0Shift = 0.0; @@ -3643,29 +5271,50 @@ 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]); - 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) { - fCurrentCoef2[0] = fCurrentCoef2[1]; + //printf("Value of timeoffset %f\n",fCurrentCoef2[0]); + if (fCurrentCoef2[0] < -3.0) { + fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0; } } else { - fCurrentCoef2[0] = fCurrentCoef2[1]; + fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0; } } else { + ////printf("Put default %f\n",-TMath::Abs(fCurrentCoef[1])); fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]); - fCurrentCoef2[0] = fCurrentCoef2[1]; - //printf("Fit failed!\n"); + + if((fPhd[1] > fPhd[0]) && + (put)) { + if (fPhdt0 >= 0.0) { + fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins; + if (fCurrentCoef2[0] < -3.0) { + fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0; + } + else fNumberFitSuccess++; + } + else { + fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0; + } + } + else{ + fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0; + //printf("Fit failed!\n"); + } } if (fDebugLevel == 1) { 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()); @@ -3673,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(); @@ -3684,11 +5334,19 @@ void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH) else { delete pentea; delete pente; - delete polynome; - delete polynomea; - delete polynomeb; + if(polynome) delete polynome; + if(polynomea) delete polynomea; + if(polynomeb) delete polynomeb; + //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); } @@ -3764,7 +5422,7 @@ void AliTRDCalibraFit::FitPH(TH1* projPH, Int_t idect) } else { fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]); - fCurrentCoef2[0] = fCurrentCoef2[1]; + fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0; } } @@ -3772,7 +5430,7 @@ void AliTRDCalibraFit::FitPH(TH1* projPH, Int_t idect) // Put the default value fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]); - fCurrentCoef2[0] = fCurrentCoef2[1]; + fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0; } if (fDebugLevel != 1) { @@ -3794,7 +5452,7 @@ Bool_t AliTRDCalibraFit::FitPRFGausMI(Double_t *arraye, Double_t *arraym, Double 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; } @@ -3857,7 +5515,7 @@ Double_t AliTRDCalibraFit::FitGausMI(Double_t *arraye, Double_t *arraym, Double_ 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++; @@ -3904,7 +5562,7 @@ Double_t AliTRDCalibraFit::FitGausMI(Double_t *arraye, Double_t *arraym, Double_ 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); @@ -4104,7 +5762,7 @@ void AliTRDCalibraFit::FitTnpRange(Double_t *arraye, Double_t *arraym, Double_t 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; @@ -4213,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); @@ -4227,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); @@ -4280,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); @@ -4348,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(); @@ -4375,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; @@ -4394,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))) { @@ -4416,74 +6171,169 @@ 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; } } //_____________________________________________________________________________ -Double_t *AliTRDCalibraFit::CalculPolynomeLagrange2(Double_t *x, 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(Double_t *x, 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(Double_t *x, 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])); @@ -4491,33 +6341,30 @@ Double_t *AliTRDCalibraFit::CalculPolynomeLagrange4(Double_t *x, Double_t *y) co 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() @@ -4571,7 +6418,7 @@ void AliTRDCalibraFit::NormierungCharge() } } //_____________________________________________________________________________ -TH1I *AliTRDCalibraFit::ReBin(TH1I *hist) const +TH1I *AliTRDCalibraFit::ReBin(const TH1I *hist) const { // // Rebin of the 1D histo for the gain calibration if needed. @@ -4600,7 +6447,7 @@ TH1I *AliTRDCalibraFit::ReBin(TH1I *hist) const } //_____________________________________________________________________________ -TH1F *AliTRDCalibraFit::ReBin(TH1F *hist) const +TH1F *AliTRDCalibraFit::ReBin(const TH1F *hist) const { // // Rebin of the 1D histo for the gain calibration if needed @@ -4626,53 +6473,6 @@ TH1F *AliTRDCalibraFit::ReBin(TH1F *hist) const return rehist; -} - -//_____________________________________________________________________________ -TH1F *AliTRDCalibraFit::CorrectTheError(TGraphErrors *hist) -{ - // - // In the case of the vectors method the trees contains TGraphErrors for PH and PRF - // to be able to add them after - // We convert it to a TH1F to be able to applied the same fit function method - // After having called this function you can not add the statistics anymore - // - - TH1F *rehist = 0x0; - - Int_t nbins = hist->GetN(); - Double_t *x = hist->GetX(); - Double_t *entries = hist->GetEX(); - Double_t *mean = hist->GetY(); - Double_t *square = hist->GetEY(); - fEntriesCurrent = 0; - - if (nbins < 2) { - return rehist; - } - - Double_t step = x[1] - x[0]; - Double_t minvalue = x[0] - step/2; - Double_t maxvalue = x[(nbins-1)] + step/2; - - rehist = new TH1F("projcorrecterror","",nbins,minvalue,maxvalue); - - for (Int_t k = 0; k < nbins; k++) { - rehist->SetBinContent(k+1,mean[k]); - if (entries[k] > 0.0) { - fEntriesCurrent += (Int_t) entries[k]; - Double_t d = TMath::Abs(square[k] - (mean[k]*mean[k])); - rehist->SetBinError(k+1,TMath::Sqrt(d/entries[k])); - } - else { - rehist->SetBinError(k+1,0.0); - } - } - - if(fEntriesCurrent > 0) fNumberEnt++; - - return rehist; - } // //____________Some basic geometry function_____________________________________ @@ -4734,7 +6534,7 @@ void AliTRDCalibraFit::ResetVectorFit() // //_____________________________________________________________________________ -Double_t AliTRDCalibraFit::PH(Double_t *x, Double_t *par) +Double_t AliTRDCalibraFit::PH(const Double_t *x, const Double_t *par) { // // Function for the fit @@ -4789,7 +6589,7 @@ Double_t AliTRDCalibraFit::PH(Double_t *x, Double_t *par) } //_____________________________________________________________________________ -Double_t AliTRDCalibraFit::AsymmGauss(Double_t *x, Double_t *par) +Double_t AliTRDCalibraFit::AsymmGauss(const Double_t *x, const Double_t *par) { // // Function for the fit @@ -4815,9 +6615,9 @@ Double_t AliTRDCalibraFit::AsymmGauss(Double_t *x, Double_t *par) 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); @@ -4825,7 +6625,7 @@ Double_t AliTRDCalibraFit::AsymmGauss(Double_t *x, Double_t *par) } //_____________________________________________________________________________ -Double_t AliTRDCalibraFit::FuncLandauGaus(Double_t *x, Double_t *par) +Double_t AliTRDCalibraFit::FuncLandauGaus(const Double_t *x, const Double_t *par) { // // Sum Landau + Gaus with identical mean @@ -4841,7 +6641,69 @@ Double_t AliTRDCalibraFit::FuncLandauGaus(Double_t *x, Double_t *par) } //_____________________________________________________________________________ -Double_t AliTRDCalibraFit::LanGauFun(Double_t *x, Double_t *par) +Double_t AliTRDCalibraFit::LanGauFun(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 + // + // In the Landau distribution (represented by the CERNLIB approximation), + // the maximum is located at x=-0.22278298 with the location parameter=0. + // This shift is corrected within this function, so that the actual + // maximum is identical to the MP parameter. + // + + // Numeric constants + Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2) + Double_t mpshift = -0.22278298; // Landau maximum location + + // Control constants + Double_t np = 100.0; // Number of convolution steps + Double_t sc = 5.0; // Convolution extends to +-sc Gaussian sigmas + + // Variables + Double_t xx = 0.0; + Double_t mpc = 0.0; + Double_t fland = 0.0; + Double_t sum = 0.0; + Double_t xlow = 0.0; + Double_t xupp = 0.0; + Double_t step = 0.0; + Double_t i = 0.0; + + // MP shift correction + mpc = par[1] - mpshift * par[0]; + + // Range of convolution integral + xlow = x[0] - sc * par[3]; + xupp = x[0] + sc * par[3]; + + step = (xupp - xlow) / np; + + // Convolution integral of Landau and Gaussian by sum + for (i = 1.0; i <= np/2; i++) { + + xx = xlow + (i-.5) * step; + if(par[0] > 0.0) fland = TMath::Landau(xx,mpc,par[0]) / 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 @@ -4851,6 +6713,7 @@ Double_t AliTRDCalibraFit::LanGauFun(Double_t *x, Double_t *par) // 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. @@ -4867,14 +6730,14 @@ Double_t AliTRDCalibraFit::LanGauFun(Double_t *x, 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; - 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]; @@ -4889,21 +6752,22 @@ Double_t AliTRDCalibraFit::LanGauFun(Double_t *x, 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; } //_____________________________________________________________________________ -TF1 *AliTRDCalibraFit::LanGauFit(TH1 *his, Double_t *fitrange, Double_t *startvalues - , Double_t *parlimitslo, Double_t *parlimitshi +TF1 *AliTRDCalibraFit::LanGauFit(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 { @@ -4927,7 +6791,7 @@ TF1 *AliTRDCalibraFit::LanGauFit(TH1 *his, Double_t *fitrange, Double_t *startva 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++) { @@ -4939,112 +6803,45 @@ TF1 *AliTRDCalibraFit::LanGauFit(TH1 *his, Double_t *fitrange, Double_t *startva return (ffit); // Return fit function } - //_____________________________________________________________________________ -Int_t AliTRDCalibraFit::LanGauPro(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(Double_t *x, Double_t *par) -{ - // - // Gaus with identical mean - // - Double_t gauss = par[0] * TMath::Gaus(x[0],0.0,par[1])+par[2]; - - return gauss; -} + +