X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TRD%2FAliTRDCalibraFit.cxx;h=7a34ac65d32438f1dc9f8ae5fdb35c54e8df5c05;hb=5812a0645a533d7abb928eaf71eecded4912f21c;hp=47ffc9d5472a30acf2b02d81c4ef3f9bf457b9a4;hpb=979b168ff690d343d13ffa7df2c324e1d29617c1;p=u%2Fmrichter%2FAliRoot.git diff --git a/TRD/AliTRDCalibraFit.cxx b/TRD/AliTRDCalibraFit.cxx index 47ffc9d5472..7a34ac65d32 100644 --- a/TRD/AliTRDCalibraFit.cxx +++ b/TRD/AliTRDCalibraFit.cxx @@ -59,11 +59,11 @@ #include #include #include -#include #include #include #include -#include +#include +#include #include "AliLog.h" #include "AliMathBase.h" @@ -72,10 +72,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" @@ -135,6 +137,7 @@ AliTRDCalibraFit::AliTRDCalibraFit() ,fAccCDB(kFALSE) ,fMinEntries(800) ,fRebin(1) + ,fScaleGain(0.02431) ,fNumberFit(0) ,fNumberFitSuccess(0) ,fNumberEnt(0) @@ -157,6 +160,8 @@ AliTRDCalibraFit::AliTRDCalibraFit() ,fCalROC(0x0) ,fCalDet2(0x0) ,fCalROC2(0x0) + ,fCalDetVdriftUsed(0x0) + ,fCalDetExBUsed(0x0) ,fCurrentCoefDetector(0x0) ,fCurrentCoefDetector2(0x0) ,fVectorFit(0) @@ -194,6 +199,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) @@ -216,6 +222,8 @@ AliTRDCalibraFit::AliTRDCalibraFit(const AliTRDCalibraFit &c) ,fCalROC(0x0) ,fCalDet2(0x0) ,fCalROC2(0x0) +,fCalDetVdriftUsed(0x0) +,fCalDetExBUsed(0x0) ,fCurrentCoefDetector(0x0) ,fCurrentCoefDetector2(0x0) ,fVectorFit(0) @@ -242,6 +250,9 @@ AliTRDCalibraFit::AliTRDCalibraFit(const AliTRDCalibraFit &c) if(c.fCalROC) fCalROC = new AliTRDCalROC(*c.fCalROC); if(c.fCalROC2) fCalROC = new AliTRDCalROC(*c.fCalROC2); + if(c.fCalDetVdriftUsed) fCalDetVdriftUsed = new AliTRDCalDet(*c.fCalDetVdriftUsed); + if(c.fCalDetExBUsed) fCalDetExBUsed = new AliTRDCalDet(*c.fCalDetExBUsed); + fVectorFit.SetName(c.fVectorFit.GetName()); for(Int_t k = 0; k < c.fVectorFit.GetEntriesFast(); k++){ AliTRDFitInfo *fitInfo = new AliTRDFitInfo(); @@ -297,6 +308,8 @@ AliTRDCalibraFit::~AliTRDCalibraFit() if ( fCalDet2 ) delete fCalDet2; if ( fCalROC ) delete fCalROC; if ( fCalROC2 ) delete fCalROC2; + if ( fCalDetVdriftUsed) delete fCalDetVdriftUsed; + if ( fCalDetExBUsed) delete fCalDetExBUsed; if( fCurrentCoefDetector ) delete [] fCurrentCoefDetector; if( fCurrentCoefDetector2 ) delete [] fCurrentCoefDetector2; fVectorFit.Delete(); @@ -361,7 +374,8 @@ Bool_t AliTRDCalibraFit::AnalyseCH(const TH2I *ch) // // Set the calibration mode - const char *name = ch->GetTitle(); + //const char *name = ch->GetTitle(); + TString name = ch->GetTitle(); if(!SetModeCalibration(name,0)) return kFALSE; // Number of Ybins (detectors or groups of pads) @@ -459,7 +473,8 @@ Bool_t AliTRDCalibraFit::AnalyseCH(AliTRDCalibraVector *calvect) // // Set the calibraMode - const char *name = calvect->GetNameCH(); + //const char *name = calvect->GetNameCH(); + TString name = calvect->GetNameCH(); if(!SetModeCalibration(name,0)) return kFALSE; // Number of Xbins (detectors or groups of pads) @@ -483,34 +498,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 @@ -527,10 +539,6 @@ Bool_t AliTRDCalibraFit::AnalyseCH(AliTRDCalibraVector *calvect) } // Fill Infos Fit FillInfosFitCH(idect); - // Memory!!! - if (fDebugLevel != 1) { - delete projch; - } } // Boucle object // Normierungcharge if (fDebugLevel != 1) { @@ -548,8 +556,55 @@ 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: FitCH((TH1 *) projch, mean); break; + case 3: FitBisCH((TH1 *) projch, mean); break; + default: return -100.0; + } + delete fDebugStreamer; + fDebugStreamer = 0x0; + + if(fCurrentCoef[0] > 0.0) return fCurrentCoef[0]; + else return -100.0; + +} //________________functions fit Online PH2d____________________________________ -Bool_t AliTRDCalibraFit::AnalysePH(const TProfile2D *ph) +Double_t AliTRDCalibraFit::AnalysePHAllTogether(const TProfile2D *ph) { // // Take the 1D profiles (average pulse height), projections of the 2D PH @@ -558,14 +613,67 @@ Bool_t AliTRDCalibraFit::AnalysePH(const TProfile2D *ph) // A first calibration of T0 is also made using the same method // - // Set the calibration mode - const char *name = ph->GetTitle(); - if(!SetModeCalibration(name,1)) return kFALSE; - // Number of Xbins (detectors or groups of pads) Int_t nbins = ph->GetNbinsX();// time Int_t nybins = ph->GetNbinsY();// calibration group - if (!InitFit(nybins,1)) { + + // Take the histo + TH1D *projph = (TH1D *) ph->ProjectionX("projph",1,nybins+1,(Option_t *) "e"); + projph->SetDirectory(0); + // Number of entries for this calibration group + Double_t nentries = 0; + for(Int_t idect = 0; idect < nybins; idect++){ + for (Int_t k = 0; k < nbins; k++) { + Int_t binnb = (nbins+2)*(idect+1)+(k+1); + nentries += ph->GetBinEntries(binnb); + } + } + //printf("AnalysePHAllTogether:: the number of entries is %f\n",nentries); + // This detector has not enough statistics or was off + if (nentries <= fMinEntries) { + AliInfo(Form("There are %d entries for all together, it is not enough (%d)",(Int_t) nentries, fMinEntries)); + if (fDebugLevel != 1) { + delete projph; + } + return -100.0; + } + //Method choosen + //printf("Method\n"); + switch(fMethod) + { + case 0: FitLagrangePoly((TH1 *) projph); break; + case 1: FitPente((TH1 *) projph); break; + case 2: FitPH((TH1 *) projph,0); break; + default: return -100.0; + } + // Memory!!! + if (fDebugLevel != 1) { + delete projph; + } + delete fDebugStreamer; + fDebugStreamer = 0x0; + + if(fCurrentCoef[0] > 0.0) return fCurrentCoef[0]; + else return -100.0; + +} +//____________Functions fit Online PH2d________________________________________ +Bool_t AliTRDCalibraFit::AnalysePH(AliTRDCalibraVector *calvect) +{ + // + // Reconstruct the average pulse height from the vectorPH for each + // calibration group + // Reconstruct a drift velocity + // A first calibration of T0 is also made using the same method (slope method) + // + + // Set the calibration mode + //const char *name = calvect->GetNamePH(); + TString name = calvect->GetNamePH(); + if(!SetModeCalibration(name,1)) return kFALSE; + + // Number of Xbins (detectors or groups of pads) + if (!InitFit((432*calvect->GetDetCha0(1)+108*calvect->GetDetCha2(1)),1)) { return kFALSE; } if (!InitFitPH()) { @@ -579,34 +687,30 @@ Bool_t AliTRDCalibraFit::AnalysePH(const 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,nentries); - if (fDebugLevel != 1) { - delete projph; - } + if (fEntriesCurrent <= fMinEntries) { + //printf("Not enough stat!\n"); + NotEnoughStatisticPH(idect,fEntriesCurrent); continue; } - // Statistics of the histos fitted + // Statistic of the histos fitted fNumberFit++; - fStatisticMean += nentries; + fStatisticMean += fEntriesCurrent; // Calcul of "real" coef CalculVdriftCoefMean(); CalculT0CoefMean(); @@ -619,12 +723,9 @@ Bool_t AliTRDCalibraFit::AnalysePH(const TProfile2D *ph) default: return kFALSE; } // Fill the tree if end of a detector or only the pointer to the branch!!! - FillInfosFitPH(idect,nentries); - // Memory!!! - if (fDebugLevel != 1) { - delete projph; - } + FillInfosFitPH(idect,fEntriesCurrent); } // Boucle object + // Mean Statistic if (fNumberFit > 0) { AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess)); @@ -637,66 +738,82 @@ Bool_t AliTRDCalibraFit::AnalysePH(const TProfile2D *ph) fDebugStreamer = 0x0; return kTRUE; } -//____________Functions fit Online PH2d________________________________________ -Bool_t AliTRDCalibraFit::AnalysePH(AliTRDCalibraVector *calvect) +//________________functions fit Online PH2d____________________________________ +Bool_t AliTRDCalibraFit::AnalysePH(const TProfile2D *ph) { // - // Reconstruct the average pulse height from the vectorPH for each - // calibration group + // Take the 1D profiles (average pulse height), projections of the 2D PH + // on the Xaxis, for each calibration group // Reconstruct a drift velocity - // A first calibration of T0 is also made using the same method (slope method) + // A first calibration of T0 is also made using the same method // // Set the calibration mode - const char *name = calvect->GetNamePH(); + //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,fEntriesCurrent); + 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; @@ -705,13 +822,12 @@ Bool_t AliTRDCalibraFit::AnalysePH(AliTRDCalibraVector *calvect) default: return kFALSE; } // Fill the tree if end of a detector or only the pointer to the branch!!! - FillInfosFitPH(idect,fEntriesCurrent); + FillInfosFitPH(idect,nentries); // Memory!!! if (fDebugLevel != 1) { delete projph; } } // Boucle object - // Mean Statistic if (fNumberFit > 0) { AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess)); @@ -734,7 +850,8 @@ Bool_t AliTRDCalibraFit::AnalysePRF(const TProfile2D *prf) // // Set the calibration mode - const char *name = prf->GetTitle(); + //const char *name = prf->GetTitle(); + TString name = prf->GetTitle(); if(!SetModeCalibration(name,2)) return kFALSE; // Number of Ybins (detectors or groups of pads) @@ -821,7 +938,8 @@ Bool_t AliTRDCalibraFit::AnalysePRFMarianFit(const TProfile2D *prf) // // Set the calibration mode - const char *name = prf->GetTitle(); + //const char *name = prf->GetTitle(); + TString name = prf->GetTitle(); if(!SetModeCalibration(name,2)) return kFALSE; // Number of Ybins (detectors or groups of pads) @@ -916,7 +1034,8 @@ Bool_t AliTRDCalibraFit::AnalysePRF(AliTRDCalibraVector *calvect) // // Set the calibra mode - const char *name = calvect->GetNamePRF(); + //const char *name = calvect->GetNamePRF(); + TString name = calvect->GetNamePRF(); if(!SetModeCalibration(name,2)) return kFALSE; //printf("test0 %s\n",name); @@ -941,22 +1060,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 @@ -973,10 +1089,6 @@ 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) { @@ -999,7 +1111,8 @@ Bool_t AliTRDCalibraFit::AnalysePRFMarianFit(AliTRDCalibraVector *calvect) // // Set the calibra mode - const char *name = calvect->GetNamePRF(); + //const char *name = calvect->GetNamePRF(); + TString name = calvect->GetNamePRF(); if(!SetModeCalibration(name,2)) return kFALSE; //printf("test0 %s\n",name); Int_t nbg = GetNumberOfGroupsPRF((const char *)name); @@ -1036,35 +1149,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 @@ -1081,10 +1192,6 @@ 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) { @@ -1148,6 +1255,7 @@ Bool_t AliTRDCalibraFit::AnalyseLinearFitters(AliTRDCalibraVdriftLinearFit *cali // 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 ++; @@ -1158,7 +1266,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]; } @@ -1178,10 +1286,212 @@ Bool_t AliTRDCalibraFit::AnalyseLinearFitters(AliTRDCalibraVdriftLinearFit *cali fDebugStreamer = 0x0; return kTRUE; +} +//______________________________________________________________________________________ +Bool_t AliTRDCalibraFit::AnalyseExbAltFit(AliTRDCalibraExbAltFit *calivdli) +{ + // + // The linear method + // + + fStatisticMean = 0.0; + fNumberFit = 0; + fNumberFitSuccess = 0; + fNumberEnt = 0; + if(!InitFitExbAlt()) return kFALSE; + + + for(Int_t idet = 0; idet < 540; idet++){ + + + //printf("detector number %d\n",idet); + + // Take the result + TVectorD param(3); + TVectorD error(3); + fEntriesCurrent = 0; + fCountDet = idet; + Bool_t here = calivdli->GetParam(idet,¶m); + Bool_t heree = calivdli->GetError(idet,&error); + //printf("here %d and heree %d\n",here, heree); + if(heree) { + fEntriesCurrent = (Int_t) error[2]; + fNumberEnt++; + } + //printf("Number of entries %d\n",fEntriesCurrent); + // Nothing found or not enough statistic + if((!heree) || (!here) || (fEntriesCurrent <= fMinEntries)) { + NotEnoughStatisticExbAlt(); + continue; + } + //param.Print(); + //error.Print(); + //Statistics + fNumberFit++; + fStatisticMean += fEntriesCurrent; + + // Statistics + fNumberFitSuccess ++; + + // Put the fCurrentCoef + if(TMath::Abs(param[2])>0.0001){ + fCurrentCoef2[0] = -param[1]/2/param[2]; + fCurrentCoefE2 = 0;//error[1]; + }else{ + fCurrentCoef2[0] = 100; + fCurrentCoefE2 = 0;//error[1]; + } + + // Fill + FillInfosFitExbAlt(); + + } + // Mean Statistics + if (fNumberFit > 0) { + AliInfo(Form("There are %d with at least one entries. %d fits have been proceeded (sucessfully or not...). There is a mean statistic of: %d over these fitted histograms and %d successfulled fits",fNumberEnt, fNumberFit, (Int_t) fStatisticMean/fNumberFit,fNumberFitSuccess)); + } + else { + AliInfo(Form("There are %d with at least one entries. There is no fit!",fNumberEnt)); + } + delete fDebugStreamer; + fDebugStreamer = 0x0; + return kTRUE; + +} +//____________Functions fit Online CH2d________________________________________ +void AliTRDCalibraFit::AnalyseLinearFittersAllTogether(AliTRDCalibraVdriftLinearFit *calivdli, Double_t &vdriftoverall, Double_t &exboverall) +{ + // + // 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 + + Int_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++; + } + else { + if(entries< 1198){ + linearfitter.AddPoint(&x,y); + entries++; + } + } + } + + } + } + } + + //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.0) fCurrentCoef2[0] = (par[0]+meanvdriftused*meanexbused)/fCurrentCoef[0]; + else fCurrentCoef2[0] = 100.0; + + } + else { + + fCurrentCoef[0] = -100.0; + fCurrentCoef2[0] = 100.0; + + } + + + } + else { + + fCurrentCoef[0] = -100.0; + fCurrentCoef2[0] = 100.0; + + } + + vdriftoverall = fCurrentCoef[0]; + exboverall = fCurrentCoef2[0]; + + + delete linearfitterhisto; + delete fDebugStreamer; + fDebugStreamer = 0x0; + } //____________Functions for seeing if the pad is really okey___________________ //_____________________________________________________________________________ -Int_t AliTRDCalibraFit::GetNumberOfGroupsPRF(const char* nametitle) +Int_t AliTRDCalibraFit::GetNumberOfGroupsPRF(TString nametitle) { // // Get numberofgroupsprf @@ -1197,25 +1507,25 @@ Int_t AliTRDCalibraFit::GetNumberOfGroupsPRF(const char* nametitle) const Char_t *pattern6 = "Ngp6"; // Nrphi mode - if (strstr(nametitle,pattern0)) { + if (strstr(nametitle.Data(),pattern0)) { return 0; } - if (strstr(nametitle,pattern1)) { + if (strstr(nametitle.Data(),pattern1)) { return 1; } - if (strstr(nametitle,pattern2)) { + if (strstr(nametitle.Data(),pattern2)) { return 2; } - if (strstr(nametitle,pattern3)) { + if (strstr(nametitle.Data(),pattern3)) { return 3; } - if (strstr(nametitle,pattern4)) { + if (strstr(nametitle.Data(),pattern4)) { return 4; } - if (strstr(nametitle,pattern5)) { + if (strstr(nametitle.Data(),pattern5)) { return 5; } - if (strstr(nametitle,pattern6)){ + if (strstr(nametitle.Data(),pattern6)){ return 6; } else return -1; @@ -1223,7 +1533,7 @@ Int_t AliTRDCalibraFit::GetNumberOfGroupsPRF(const char* nametitle) } //_____________________________________________________________________________ -Bool_t AliTRDCalibraFit::SetModeCalibration(const char *name, Int_t i) +Bool_t AliTRDCalibraFit::SetModeCalibration(TString name, Int_t i) { // // Set fNz[i] and fNrphi[i] of the AliTRDCalibraFit::Instance() @@ -1237,7 +1547,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() @@ -1260,7 +1570,7 @@ Bool_t AliTRDCalibraFit::SetNrphiFromTObject(const char *name, Int_t i) const Char_t *patternz100 = "Nz100"; // Nrphi mode - if ((strstr(name,patternrphi100)) && (strstr(name,patternz100))) { + if ((strstr(name.Data(),patternrphi100)) && (strstr(name.Data(),patternz100))) { fCalibraMode->SetAllTogether(i); fNbDet = 540; if (fDebugLevel > 1) { @@ -1268,7 +1578,7 @@ Bool_t AliTRDCalibraFit::SetNrphiFromTObject(const char *name, Int_t i) } return kTRUE; } - if ((strstr(name,patternrphi10)) && (strstr(name,patternz10))) { + if ((strstr(name.Data(),patternrphi10)) && (strstr(name.Data(),patternz10))) { fCalibraMode->SetPerSuperModule(i); fNbDet = 30; if (fDebugLevel > 1) { @@ -1277,49 +1587,49 @@ Bool_t AliTRDCalibraFit::SetNrphiFromTObject(const char *name, Int_t i) return kTRUE; } - if (strstr(name,patternrphi0)) { + if (strstr(name.Data(),patternrphi0)) { fCalibraMode->SetNrphi(i ,0); if (fDebugLevel > 1) { AliInfo(Form("fNbDet %d and 0",fNbDet)); } return kTRUE; } - if (strstr(name,patternrphi1)) { + if (strstr(name.Data(),patternrphi1)) { fCalibraMode->SetNrphi(i, 1); if (fDebugLevel > 1) { AliInfo(Form("fNbDet %d and 1",fNbDet)); } return kTRUE; } - if (strstr(name,patternrphi2)) { + if (strstr(name.Data(),patternrphi2)) { fCalibraMode->SetNrphi(i, 2); if (fDebugLevel > 1) { AliInfo(Form("fNbDet %d and 2",fNbDet)); } return kTRUE; } - if (strstr(name,patternrphi3)) { + if (strstr(name.Data(),patternrphi3)) { fCalibraMode->SetNrphi(i, 3); if (fDebugLevel > 1) { AliInfo(Form("fNbDet %d and 3",fNbDet)); } return kTRUE; } - if (strstr(name,patternrphi4)) { + if (strstr(name.Data(),patternrphi4)) { fCalibraMode->SetNrphi(i, 4); if (fDebugLevel > 1) { AliInfo(Form("fNbDet %d and 4",fNbDet)); } return kTRUE; } - if (strstr(name,patternrphi5)) { + if (strstr(name.Data(),patternrphi5)) { fCalibraMode->SetNrphi(i, 5); if (fDebugLevel > 1) { AliInfo(Form("fNbDet %d and 5",fNbDet)); } return kTRUE; } - if (strstr(name,patternrphi6)) { + if (strstr(name.Data(),patternrphi6)) { fCalibraMode->SetNrphi(i, 6); if (fDebugLevel > 1) { AliInfo(Form("fNbDet %d and 6",fNbDet)); @@ -1335,7 +1645,7 @@ Bool_t AliTRDCalibraFit::SetNrphiFromTObject(const char *name, Int_t i) } //_____________________________________________________________________________ -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() @@ -1354,7 +1664,7 @@ Bool_t AliTRDCalibraFit::SetNzFromTObject(const char *name, Int_t i) const Char_t *patternz10 = "Nz10"; const Char_t *patternz100 = "Nz100"; - if ((strstr(name,patternrphi100)) && (strstr(name,patternz100))) { + if ((strstr(name.Data(),patternrphi100)) && (strstr(name.Data(),patternz100))) { fCalibraMode->SetAllTogether(i); fNbDet = 540; if (fDebugLevel > 1) { @@ -1362,7 +1672,7 @@ Bool_t AliTRDCalibraFit::SetNzFromTObject(const char *name, Int_t i) } return kTRUE; } - if ((strstr(name,patternrphi10)) && (strstr(name,patternz10))) { + if ((strstr(name.Data(),patternrphi10)) && (strstr(name.Data(),patternz10))) { fCalibraMode->SetPerSuperModule(i); fNbDet = 30; if (fDebugLevel > 1) { @@ -1370,47 +1680,210 @@ Bool_t AliTRDCalibraFit::SetNzFromTObject(const char *name, Int_t i) } return kTRUE; } - if (strstr(name,patternz0)) { + if (strstr(name.Data(),patternz0)) { fCalibraMode->SetNz(i, 0); if (fDebugLevel > 1) { AliInfo(Form("fNbDet %d and 0",fNbDet)); } return kTRUE; } - if (strstr(name,patternz1)) { + if (strstr(name.Data(),patternz1)) { fCalibraMode->SetNz(i ,1); if (fDebugLevel > 1) { AliInfo(Form("fNbDet %d and 1",fNbDet)); } return kTRUE; } - if (strstr(name,patternz2)) { + if (strstr(name.Data(),patternz2)) { fCalibraMode->SetNz(i ,2); if (fDebugLevel > 1) { AliInfo(Form("fNbDet %d and 2",fNbDet)); } return kTRUE; } - if (strstr(name,patternz3)) { + 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)); + 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 } - fCalibraMode->SetNz(i ,0); - return kFALSE; } //______________________________________________________________________ void AliTRDCalibraFit::PutMeanValueOtherVectorFit(Int_t ofwhat, Bool_t perdetector){ @@ -1436,17 +1909,23 @@ void AliTRDCalibraFit::PutMeanValueOtherVectorFit(Int_t ofwhat, Bool_t perdetect 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++) { @@ -1455,10 +1934,13 @@ void AliTRDCalibraFit::PutMeanValueOtherVectorFit(Int_t ofwhat, Bool_t perdetect 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++; } @@ -1470,10 +1952,13 @@ void AliTRDCalibraFit::PutMeanValueOtherVectorFit(Int_t ofwhat, Bool_t perdetect 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++; } @@ -1482,13 +1967,24 @@ void AliTRDCalibraFit::PutMeanValueOtherVectorFit(Int_t ofwhat, Bool_t perdetect } // Row } } - if(countAll > 0) meanAll = meanAll/countAll; + if(countAll > 0) { + meanAll = meanAll/countAll; + rmsAll = TMath::Abs(rmsAll/countAll - (meanAll*meanAll)); + } for(Int_t sm = 0; sm < 18; sm++){ - if(countSupermodule[sm] > 0) meanSupermodule[sm] = meanSupermodule[sm]/countSupermodule[sm]; + if(countSupermodule[sm] > 0) { + meanSupermodule[sm] = meanSupermodule[sm]/countSupermodule[sm]; + rmsSupermodule[sm] = TMath::Abs(rmsSupermodule[sm]/countSupermodule[sm] - (meanSupermodule[sm]*meanSupermodule[sm])); + } } for(Int_t det = 0; det < 540; det++){ - if(countDetector[det] > 0) meanDetector[det] = meanDetector[det]/countDetector[det]; + if(countDetector[det] > 0) { + meanDetector[det] = meanDetector[det]/countDetector[det]; + rmsDetector[det] = TMath::Abs(rmsDetector[det]/countDetector[det] - (meanDetector[det]*meanDetector[det])); + } } + //printf("Put mean value, meanAll %f, rmsAll %f\n",meanAll,rmsAll); + /////////////////////////////////////////////// // Put the mean value for the no-fitted ///////////////////////////////////////////// for (Int_t k = 0; k < loop; k++) { @@ -1502,11 +1998,11 @@ void AliTRDCalibraFit::PutMeanValueOtherVectorFit(Int_t ofwhat, Bool_t perdetect for (Int_t col = 0; col < colMax; col++) { value = coef[(Int_t)(col*rowMax+row)]; if(value < 0.0) { - if((ofwhat == 0) && (meanAll > 0.0)) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanAll); + if((ofwhat == 0) && (meanAll > 0.0) && (countAll > 15)) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanAll); if(ofwhat == 1){ - if(meanDetector[detector] > 0.0) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanDetector[detector]); - else if(meanSupermodule[sector] > 0.0) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanSupermodule[sector]); - else if(meanAll > 0.0) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanAll); + if((meanDetector[detector] > 0.0) && (countDetector[detector] > 20)) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanDetector[detector]); + else if((meanSupermodule[sector] > 0.0) && (countSupermodule[sector] > 15)) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanSupermodule[sector]); + else if((meanAll > 0.0) && (countAll > 15)) coef[(Int_t)(col*rowMax+row)] = -TMath::Abs(meanAll); } } // Debug @@ -1533,7 +2029,6 @@ void AliTRDCalibraFit::PutMeanValueOtherVectorFit(Int_t ofwhat, Bool_t perdetect } // Col } // Row } - } //______________________________________________________________________ void AliTRDCalibraFit::PutMeanValueOtherVectorFit2(Int_t ofwhat, Bool_t perdetector){ @@ -1557,16 +2052,21 @@ void AliTRDCalibraFit::PutMeanValueOtherVectorFit2(Int_t ofwhat, Bool_t perdetec // Initialisation //////////////////////// Double_t meanAll = 0.0; + Double_t rmsAll = 0.0; Double_t meanSupermodule[18]; + Double_t rmsSupermodule[18]; Double_t meanDetector[540]; + Double_t rmsDetector[540]; Int_t countAll = 0; Int_t countSupermodule[18]; Int_t countDetector[540]; for(Int_t sm = 0; sm < 18; sm++){ + rmsSupermodule[sm] = 0.0; meanSupermodule[sm] = 0.0; countSupermodule[sm] = 0; } for(Int_t det = 0; det < 540; det++){ + rmsDetector[det] = 0.0; meanDetector[det] = 0.0; countDetector[det] = 0; } @@ -1578,11 +2078,14 @@ void AliTRDCalibraFit::PutMeanValueOtherVectorFit2(Int_t ofwhat, Bool_t perdetec 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++; } } @@ -1593,10 +2096,13 @@ void AliTRDCalibraFit::PutMeanValueOtherVectorFit2(Int_t ofwhat, Bool_t perdetec 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++; } @@ -1605,13 +2111,24 @@ void AliTRDCalibraFit::PutMeanValueOtherVectorFit2(Int_t ofwhat, Bool_t perdetec } // Row } } - if(countAll > 0) meanAll = meanAll/countAll; + if(countAll > 0) { + meanAll = meanAll/countAll; + rmsAll = TMath::Abs(rmsAll/countAll - (meanAll*meanAll)); + } for(Int_t sm = 0; sm < 18; sm++){ - if(countSupermodule[sm] > 0) meanSupermodule[sm] = meanSupermodule[sm]/countSupermodule[sm]; + if(countSupermodule[sm] > 0) { + meanSupermodule[sm] = meanSupermodule[sm]/countSupermodule[sm]; + rmsSupermodule[sm] = TMath::Abs(rmsSupermodule[sm]/countSupermodule[sm] - (meanSupermodule[sm]*meanSupermodule[sm])); + } } for(Int_t det = 0; det < 540; det++){ - if(countDetector[det] > 0) meanDetector[det] = meanDetector[det]/countDetector[det]; + if(countDetector[det] > 0) { + meanDetector[det] = meanDetector[det]/countDetector[det]; + rmsDetector[det] = TMath::Abs(rmsDetector[det]/countDetector[det] - (meanDetector[det]*meanDetector[det])); + } } + //printf("Put mean value 2: meanAll %f, rmsAll %f\n",meanAll,rmsAll); + //////////////////////////////////////////// // Put the mean value for the no-fitted ///////////////////////////////////////////// for (Int_t k = 0; k < loop; k++) { @@ -1625,11 +2142,11 @@ void AliTRDCalibraFit::PutMeanValueOtherVectorFit2(Int_t ofwhat, Bool_t perdetec for (Int_t col = 0; col < colMax; col++) { value = coef[(Int_t)(col*rowMax+row)]; if(value > 70.0) { - if((ofwhat == 0) && (meanAll > 0.0)) coef[(Int_t)(col*rowMax+row)] = meanAll+100.0; + if((ofwhat == 0) && (meanAll > -3.0) && (countAll > 15)) coef[(Int_t)(col*rowMax+row)] = meanAll+100.0; if(ofwhat == 1){ - if(meanDetector[detector] > 0.0) coef[(Int_t)(col*rowMax+row)] = meanDetector[detector]+100.0; - else if(meanSupermodule[sector] > 0.0) coef[(Int_t)(col*rowMax+row)] = meanSupermodule[sector]+100.0; - else if(meanAll > 0.0) coef[(Int_t)(col*rowMax+row)] = meanAll+100.0; + if((meanDetector[detector] > -3.0) && (countDetector[detector] > 20)) coef[(Int_t)(col*rowMax+row)] = meanDetector[detector]+100.0; + else if((meanSupermodule[sector] > -3.0) && (countSupermodule[sector] > 15)) coef[(Int_t)(col*rowMax+row)] = meanSupermodule[sector]+100.0; + else if((meanAll > -3.0) && (countAll > 15)) coef[(Int_t)(col*rowMax+row)] = meanAll+100.0; } } // Debug @@ -1712,6 +2229,7 @@ AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectGain(const TObjArray *vectorFit, // Create the DetObject AliTRDCalDet *object = new AliTRDCalDet("ChamberGainFactor","GainFactor (detector value)"); + fScaleGain = scaleFitFactor; Int_t loop = (Int_t) vectorFit->GetEntriesFast(); if(loop != 540) AliInfo("The Vector Fit is not complete!"); @@ -1746,6 +2264,7 @@ AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectGain(const TObjArray *vectorFit, } // Row if(count > 0) mean = mean/count; } + if(mean < 0.1) mean = 0.1; object->SetValue(detector,mean); } @@ -1773,6 +2292,7 @@ AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectT0(const TObjArray *vectorFit, Bo Float_t min = 100.0; if(perdetector){ value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0]; + //printf("Create det object %f for %d\n",value,k); // check successful if(value > 70.0) value = value-100.0; // @@ -1831,7 +2351,49 @@ AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectLorentzAngle(const TObjArray *vec if(count > 0) mean = mean/count; */ value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0]; - object->SetValue(detector,-TMath::Abs(value)); + if(value > 70.0) value = value-100.0; + object->SetValue(detector,value); + } + + return object; + +} +//_____________________________________________________________________________ +AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectExbAlt(const TObjArray *vectorFit) +{ + // + // It creates the AliTRDCalDet object from the AliTRDFitInfo2 + // It takes the min value of the coefficients per detector + // This object has to be written in the database + // + + // Create the DetObject + AliTRDCalDet *object = new AliTRDCalDet("tan(lorentzangle)","tan(lorentzangle) (detector value)"); + + + Int_t loop = (Int_t) vectorFit->GetEntriesFast(); + if(loop != 540) AliInfo("The Vector Fit is not complete!"); + Int_t detector = -1; + Float_t value = 0.0; + + for (Int_t k = 0; k < loop; k++) { + detector = ((AliTRDFitInfo *) vectorFit->At(k))->GetDetector(); + /* + Int_t rowMax = fGeo->GetRowMax(GetLayer(detector),GetStack(detector),GetSector(detector)); + Int_t colMax = fGeo->GetColMax(GetLayer(detector)); + Float_t min = 100.0; + for (Int_t row = 0; row < rowMax; row++) { + for (Int_t col = 0; col < colMax; col++) { + value = ((AliTRDFitInfo *) fVectorFit2.At(k))->GetCoef()[(Int_t)(col*rowMax+row)]; + mean += -TMath::Abs(value); + count++; + } // Col + } // Row + if(count > 0) mean = mean/count; + */ + value = ((AliTRDFitInfo *) vectorFit->At(k))->GetCoef()[0]; + //if(value > 70.0) value = value-100.0; + object->SetValue(detector,value); } return object; @@ -1870,7 +2432,7 @@ TObject *AliTRDCalibraFit::CreatePadObjectGain(const TObjArray *vectorFit, Doubl 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++) { @@ -2327,6 +2889,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; @@ -2371,11 +2934,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; @@ -2423,6 +2987,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; @@ -2444,6 +3009,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++) { @@ -2451,45 +3018,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) { @@ -2514,10 +3063,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 @@ -2661,8 +3216,7 @@ Bool_t AliTRDCalibraFit::NotEnoughStatisticCH(Int_t idect) else if (fNbDet > 0){ Int_t firstdetector = fCountDet; Int_t lastdetector = fCountDet+fNbDet; - AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted" - ,idect,firstdetector,lastdetector)); + //AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted",idect,firstdetector,lastdetector)); // loop over detectors for(Int_t det = firstdetector; det < lastdetector; det++){ @@ -2743,8 +3297,7 @@ Bool_t AliTRDCalibraFit::NotEnoughStatisticCH(Int_t idect) } else { - AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted" - ,idect-(fCount-(fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))),fCountDet)); +//AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted",idect-(fCount-(fCalibraMode->GetNfragZ(0)*fCalibraMode->GetNfragRphi(0))),fCountDet)); // Calcul the coef from the database choosen CalculChargeCoefMean(kFALSE); @@ -2787,8 +3340,7 @@ Bool_t AliTRDCalibraFit::NotEnoughStatisticPH(Int_t idect,Double_t nentries) Int_t firstdetector = fCountDet; Int_t lastdetector = fCountDet+fNbDet; - AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted" - ,idect,firstdetector,lastdetector)); +//AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted",idect,firstdetector,lastdetector)); // loop over detectors for(Int_t det = firstdetector; det < lastdetector; det++){ @@ -2884,8 +3436,7 @@ Bool_t AliTRDCalibraFit::NotEnoughStatisticPH(Int_t idect,Double_t nentries) } else { - AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted" - ,idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))),fCountDet)); +//AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted",idect-(fCount-(fCalibraMode->GetNfragZ(1)*fCalibraMode->GetNfragRphi(1))),fCountDet)); CalculVdriftCoefMean(); CalculT0CoefMean(); @@ -2935,8 +3486,7 @@ Bool_t AliTRDCalibraFit::NotEnoughStatisticPRF(Int_t idect) Int_t firstdetector = fCountDet; Int_t lastdetector = fCountDet+fNbDet; - AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted" - ,idect,firstdetector,lastdetector)); +// AliInfo(Form("The element %d containing the detectors %d to %d has not enough statistic to be fitted",idect,firstdetector,lastdetector)); // loop over detectors for(Int_t det = firstdetector; det < lastdetector; det++){ @@ -3019,8 +3569,7 @@ Bool_t AliTRDCalibraFit::NotEnoughStatisticPRF(Int_t idect) } else { - AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted" - ,idect-(fCount-(fCalibraMode->GetNfragZ(2)*fCalibraMode->GetNfragRphi(2))),fCountDet)); +// AliInfo(Form("The element %d in this detector %d has not enough statistic to be fitted",idect-(fCount-(fCalibraMode->GetNfragZ(2)*fCalibraMode->GetNfragRphi(2))),fCountDet)); CalculPRFCoefMean(); @@ -3068,14 +3617,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(); @@ -3083,6 +3632,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) { @@ -3095,8 +3671,7 @@ Bool_t AliTRDCalibraFit::FillInfosFitCH(Int_t idect) if (fNbDet > 0){ Int_t firstdetector = fCountDet; Int_t lastdetector = fCountDet+fNbDet; - AliInfo(Form("The element %d containing the detectors %d to %d has been fitted" - ,idect,firstdetector,lastdetector)); + // AliInfo(Form("The element %d containing the detectors %d to %d has been fitted",idect,firstdetector,lastdetector)); // loop over detectors for(Int_t det = firstdetector; det < lastdetector; det++){ @@ -3207,8 +3782,7 @@ Bool_t AliTRDCalibraFit::FillInfosFitPH(Int_t idect,Double_t nentries) Int_t firstdetector = fCountDet; Int_t lastdetector = fCountDet+fNbDet; - AliInfo(Form("The element %d containing the detectors %d to %d has been fitted" - ,idect,firstdetector,lastdetector)); +// AliInfo(Form("The element %d containing the detectors %d to %d has been fitted",idect,firstdetector,lastdetector)); // loop over detectors for(Int_t det = firstdetector; det < lastdetector; det++){ @@ -3337,8 +3911,7 @@ Bool_t AliTRDCalibraFit::FillInfosFitPRF(Int_t idect) Int_t firstdetector = fCountDet; Int_t lastdetector = fCountDet+fNbDet; - AliInfo(Form("The element %d containing the detectors %d to %d has been fitted" - ,idect,firstdetector,lastdetector)); +// AliInfo(Form("The element %d containing the detectors %d to %d has been fitted",idect,firstdetector,lastdetector)); // loop over detectors for(Int_t det = firstdetector; det < lastdetector; det++){ @@ -3459,6 +4032,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) @@ -3681,6 +4276,56 @@ void AliTRDCalibraFit::FillFillLinearFitter() "\n"; } +} +//________________________________________________________________________________ +void AliTRDCalibraFit::FillFillExbAlt() +{ + // + // DebugStream and fVectorFit + // + + // End of one detector + FillVectorFit2(); + + + // Reset + for (Int_t k = 0; k < 2304; k++) { + fCurrentCoefDetector2[k] = 0.0; + } + + + if(fDebugLevel > 1){ + + if ( !fDebugStreamer ) { + //debug stream + TDirectory *backup = gDirectory; + fDebugStreamer = new TTreeSRedirector("TRDDebugFitExbAlt.root"); + if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer + } + + //Debug: comparaison of the different methods (okey for first time but not for iterative procedure) + AliTRDpadPlane *padplane = fGeo->GetPadPlane(GetLayer(fCountDet),GetStack(fCountDet)); + Float_t rowmd = (padplane->GetRow0()+padplane->GetRowEnd())/2.; + Float_t r = AliTRDgeometry::GetTime0(GetLayer(fCountDet)); + Float_t tiltangle = padplane->GetTiltingAngle(); + Int_t detector = fCountDet; + Int_t stack = GetStack(fCountDet); + Int_t layer = GetLayer(fCountDet); + Float_t vf = fCurrentCoef2[0]; + Float_t vfE = fCurrentCoefE2; + + (* fDebugStreamer) << "FillFillLinearFitter"<< + "detector="<GetValue(fCountDet); - fCurrentCoef2[1] = fCalDet2->GetValue(fCountDet); + fCurrentCoef[1] = fCalDetVdriftUsed->GetValue(fCountDet); + fCurrentCoef2[1] = fCalDetExBUsed->GetValue(fCountDet); return kTRUE; } @@ -3970,11 +4615,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]; } } @@ -4001,10 +4646,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){ @@ -4035,10 +4680,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; @@ -4061,7 +4706,7 @@ void AliTRDCalibraFit::FitPente(TH1* projPH) if (fPhdt0 >= 0.0) { fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins; - if (fCurrentCoef2[0] < -1.0) { + if (fCurrentCoef2[0] < -3.0) { fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0; } } @@ -4105,13 +4750,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; @@ -4134,7 +4781,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(); @@ -4174,8 +4825,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); @@ -4188,7 +4839,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){ @@ -4204,7 +4855,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); @@ -4217,7 +4868,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); @@ -4232,7 +4883,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; @@ -4241,7 +4892,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; @@ -4297,7 +4948,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: @@ -4311,7 +4962,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) @@ -4328,7 +4979,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); @@ -4341,7 +4992,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); @@ -4356,7 +5007,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; @@ -4364,7 +5015,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; @@ -4393,20 +5044,20 @@ void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH) if (binmin <= 1) { binmin = 2; put = 1; - AliInfo("Put the binmax from 1 to 2 to enable the fit"); + //AliInfo("Put the binmax from 1 to 2 to enable the fit"); } //check if((projPH->GetBinContent(binmin)-projPH->GetBinError(binmin)) < (projPH->GetBinContent(binmin+1))) { - AliInfo("Too many fluctuations at the end!"); + //AliInfo("Too many fluctuations at the end!"); put = kFALSE; } if((projPH->GetBinContent(binmin)+projPH->GetBinError(binmin)) > (projPH->GetBinContent(binmin-1))) { - AliInfo("Too many fluctuations at the end!"); + //AliInfo("Too many fluctuations at the end!"); put = kFALSE; } - if(pente->GetBinContent(binmin+1)==0){ - AliInfo("No entries for the next bin!"); + if(TMath::Abs(pente->GetBinContent(binmin+1)) <= 0.0000000000001){ + //AliInfo("No entries for the next bin!"); pente->SetBinContent(binmin,0); if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin(); } @@ -4445,7 +5096,7 @@ 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))) { @@ -4456,7 +5107,7 @@ void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH) (pente->GetBinContent(binmin+3) <= pente->GetBinContent(binmin+2)) && ((binmin-3) >= TMath::Min(binmax+4, projPH->GetNbinsX())) && (pente->GetBinContent(binmin-3) <= pente->GetBinContent(binmin-2))) { - AliInfo("polynome 4 false 2"); + //AliInfo("polynome 4 false 2"); put = kFALSE; } // poly 3 @@ -4487,7 +5138,7 @@ 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))) { @@ -4509,7 +5160,7 @@ 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))) { //AliInfo("polynome 3+ case 2"); @@ -4527,7 +5178,7 @@ 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))) { //AliInfo("polynome 2+ false"); @@ -4546,7 +5197,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 } @@ -4562,7 +5213,7 @@ 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 @@ -4574,23 +5225,23 @@ void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH) } if((binmin == (nbins-1)) && ((binmin-2) < TMath::Min(binmax+4, projPH->GetNbinsX()))) { put = kFALSE; - AliInfo("At the limit for the drift and not usable!"); + //AliInfo("At the limit for the drift and not usable!"); } //pass if((binmin == (nbins-2)) && ((binmin-1) < TMath::Min(binmax+4, projPH->GetNbinsX()))){ put = kFALSE; - AliInfo("For the drift...problem!"); + //AliInfo("For the drift...problem!"); } //pass but should not happen if((binmin <= (nbins-3)) && (binmin < TMath::Min(binmax+6, projPH->GetNbinsX()))){ put = kFALSE; - AliInfo("For the drift...problem!"); + //AliInfo("For the drift...problem!"); } if(put) { 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; @@ -4626,10 +5277,12 @@ void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH) (fPhd[1] > fPhd[0]) && (put)) { fCurrentCoef[0] = (kDrWidth) / (fPhd[2]-fPhd[1]); - fNumberFitSuccess++; + if(fCurrentCoef[0] > 2.5) fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]); + else fNumberFitSuccess++; if (fPhdt0 >= 0.0) { fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins; - if (fCurrentCoef2[0] < -1.0) { + //printf("Value of timeoffset %f\n",fCurrentCoef2[0]); + if (fCurrentCoef2[0] < -3.0) { fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0; } } @@ -4638,11 +5291,14 @@ void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH) } } else { + ////printf("Put default %f\n",-TMath::Abs(fCurrentCoef[1])); + fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]); + if((fPhd[1] > fPhd[0]) && (put)) { if (fPhdt0 >= 0.0) { fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins; - if (fCurrentCoef2[0] < -1.0) { + if (fCurrentCoef2[0] < -3.0) { fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0; } } @@ -4651,7 +5307,6 @@ void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH) } } else{ - fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]); fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0; //printf("Fit failed!\n"); } @@ -4677,14 +5332,13 @@ void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH) pente->Draw(); } else { - if(pentea) delete pentea; - if(pente) delete pente; + delete pentea; + delete pente; if(polynome) delete polynome; if(polynomea) delete polynomea; if(polynomeb) delete polynomeb; - if(x) delete [] x; - if(y) delete [] y; - if(c) delete [] c; + //if(x) delete [] x; + //if(y) delete [] y; if(line) delete line; } @@ -4692,7 +5346,7 @@ void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH) //Provisoire //if(fCurrentCoef[0] > 1.7) fCurrentCoef[0] = -TMath::Abs(fCurrentCoef[1]); //if((fCurrentCoef2[0] > 2.6) || (fCurrentCoef2[0] < 2.1)) fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0; - + //printf("Value of timeoffset final %f\n",fCurrentCoef2[0]); projPH->SetDirectory(0); } @@ -4798,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; } @@ -4861,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++; @@ -4908,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); @@ -5108,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; @@ -5426,68 +6080,59 @@ void AliTRDCalibraFit::FitBisCH(TH1* projch, Double_t mean) } } //_____________________________________________________________________________ -Double_t *AliTRDCalibraFit::CalculPolynomeLagrange2(const Double_t *x, const Double_t *y) const +void AliTRDCalibraFit::CalculPolynomeLagrange2(const Double_t *x, const Double_t *y, Double_t &c0, Double_t &c1, Double_t &c2, Double_t &c3, Double_t &c4) const { // // Calcul the coefficients of the polynome passant par ces trois points de degre 2 // - Double_t *c = new Double_t[5]; Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2])); Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2])); Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1])); - c[4] = 0.0; - c[3] = 0.0; - c[2] = x0+x1+x2; - c[1] = -(x0*(x[1]+x[2])+x1*(x[0]+x[2])+x2*(x[0]+x[1])); - c[0] = x0*x[1]*x[2]+x1*x[0]*x[2]+x2*x[0]*x[1]; - - return c; - + c4 = 0.0; + c3 = 0.0; + c2 = x0+x1+x2; + c1 = -(x0*(x[1]+x[2])+x1*(x[0]+x[2])+x2*(x[0]+x[1])); + c0 = x0*x[1]*x[2]+x1*x[0]*x[2]+x2*x[0]*x[1]; } //_____________________________________________________________________________ -Double_t *AliTRDCalibraFit::CalculPolynomeLagrange3(const Double_t *x, const Double_t *y) const +void AliTRDCalibraFit::CalculPolynomeLagrange3(const Double_t *x, const Double_t *y, Double_t &c0, Double_t &c1, Double_t &c2, Double_t &c3, Double_t &c4) const { // // Calcul the coefficients of the polynome passant par ces quatre points de degre 3 // - Double_t *c = new Double_t[5]; Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2])*(x[0]-x[3])); Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2])*(x[1]-x[3])); Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1])*(x[2]-x[3])); Double_t x3 = y[3]/((x[3]-x[0])*(x[3]-x[1])*(x[3]-x[2])); - c[4] = 0.0; - c[3] = x0+x1+x2+x3; - c[2] = -(x0*(x[1]+x[2]+x[3]) + c4 = 0.0; + c3 = x0+x1+x2+x3; + c2 = -(x0*(x[1]+x[2]+x[3]) +x1*(x[0]+x[2]+x[3]) +x2*(x[0]+x[1]+x[3]) +x3*(x[0]+x[1]+x[2])); - c[1] = (x0*(x[1]*x[2]+x[1]*x[3]+x[2]*x[3]) + c1 = (x0*(x[1]*x[2]+x[1]*x[3]+x[2]*x[3]) +x1*(x[0]*x[2]+x[0]*x[3]+x[2]*x[3]) +x2*(x[0]*x[1]+x[0]*x[3]+x[1]*x[3]) +x3*(x[0]*x[1]+x[0]*x[2]+x[1]*x[2])); - c[0] = -(x0*x[1]*x[2]*x[3] + c0 = -(x0*x[1]*x[2]*x[3] +x1*x[0]*x[2]*x[3] +x2*x[0]*x[1]*x[3] +x3*x[0]*x[1]*x[2]); - return c; - - } //_____________________________________________________________________________ -Double_t *AliTRDCalibraFit::CalculPolynomeLagrange4(const Double_t *x, const Double_t *y) const +void AliTRDCalibraFit::CalculPolynomeLagrange4(const Double_t *x, const Double_t *y, Double_t &c0, Double_t &c1, Double_t &c2, Double_t &c3, Double_t &c4) const { // // Calcul the coefficients of the polynome passant par ces cinqs points de degre 4 // - Double_t *c = new Double_t[5]; Double_t x0 = y[0]/((x[0]-x[1])*(x[0]-x[2])*(x[0]-x[3])*(x[0]-x[4])); Double_t x1 = y[1]/((x[1]-x[0])*(x[1]-x[2])*(x[1]-x[3])*(x[1]-x[4])); Double_t x2 = y[2]/((x[2]-x[0])*(x[2]-x[1])*(x[2]-x[3])*(x[2]-x[4])); @@ -5495,33 +6140,30 @@ Double_t *AliTRDCalibraFit::CalculPolynomeLagrange4(const Double_t *x, const Dou 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() @@ -5630,53 +6272,6 @@ TH1F *AliTRDCalibraFit::ReBin(const TH1F *hist) const return rehist; -} - -//_____________________________________________________________________________ -TH1F *AliTRDCalibraFit::CorrectTheError(const 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_____________________________________ @@ -5819,9 +6414,9 @@ Double_t AliTRDCalibraFit::AsymmGauss(const Double_t *x, const 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); @@ -6052,3 +6647,6 @@ Double_t AliTRDCalibraFit::GausConstant(const Double_t *x, const Double_t *par) return gauss; } + + +