X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TRD%2FAliTRDCalibraFit.cxx;h=f423656ae8f81d79473ba0218fa6a2417337e62a;hb=811ffdcce3ce8110476abd638c12862aee70fbab;hp=e2c48d16b7aa93fcaf55a7d612b0ac217a3783f0;hpb=053767a494091eeffd93c81d4c088dde743b7fdb;p=u%2Fmrichter%2FAliRoot.git diff --git a/TRD/AliTRDCalibraFit.cxx b/TRD/AliTRDCalibraFit.cxx index e2c48d16b7a..f423656ae8f 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,11 @@ #include #include #include -#include #include #include #include -#include +#include +#include #include "AliLog.h" #include "AliMathBase.h" @@ -76,6 +76,7 @@ #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 +105,6 @@ AliTRDCalibraFit *AliTRDCalibraFit::Instance() return fgInstance; } - //______________________________________________________________________________________ void AliTRDCalibraFit::Terminate() { @@ -121,7 +121,6 @@ void AliTRDCalibraFit::Terminate() } } - //______________________________________________________________________________________ AliTRDCalibraFit::AliTRDCalibraFit() :TObject() @@ -154,6 +153,7 @@ AliTRDCalibraFit::AliTRDCalibraFit() ,fEntriesCurrent(0) ,fCountDet(0) ,fCount(0) + ,fNbDet(0) ,fCalDet(0x0) ,fCalROC(0x0) ,fCalDet2(0x0) @@ -212,6 +212,7 @@ AliTRDCalibraFit::AliTRDCalibraFit(const AliTRDCalibraFit &c) ,fEntriesCurrent(c.fEntriesCurrent) ,fCountDet(c.fCountDet) ,fCount(c.fCount) +,fNbDet(c.fNbDet) ,fCalDet(0x0) ,fCalROC(0x0) ,fCalDet2(0x0) @@ -238,7 +239,7 @@ 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); @@ -296,7 +297,9 @@ AliTRDCalibraFit::~AliTRDCalibraFit() if ( fCalDet ) delete fCalDet; if ( fCalDet2 ) delete fCalDet2; if ( fCalROC ) delete fCalROC; - if ( fCalROC2 ) delete fCalROC2; + if ( fCalROC2 ) delete fCalROC2; + if( fCurrentCoefDetector ) delete [] fCurrentCoefDetector; + if( fCurrentCoefDetector2 ) delete [] fCurrentCoefDetector2; fVectorFit.Delete(); fVectorFit2.Delete(); if (fGeo) { @@ -316,9 +319,20 @@ void AliTRDCalibraFit::Destroy() fgInstance = 0x0; } +} +//_____________________________________________________________________________ +void AliTRDCalibraFit::DestroyDebugStreamer() +{ + // + // Delete DebugStreamer + // + + if ( fDebugStreamer ) delete fDebugStreamer; + fDebugStreamer = 0x0; + } //__________________________________________________________________________________ -void AliTRDCalibraFit::RangeChargeIntegration(Float_t vdrift, Float_t t0, Int_t &begin, Int_t &peak, Int_t &end) +void AliTRDCalibraFit::RangeChargeIntegration(Float_t vdrift, Float_t t0, Int_t &begin, Int_t &peak, Int_t &end) const { // // From the drift velocity and t0 @@ -340,7 +354,7 @@ void AliTRDCalibraFit::RangeChargeIntegration(Float_t vdrift, Float_t t0, Int_t } //____________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 +362,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 @@ -446,8 +461,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 +486,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 @@ -514,10 +527,6 @@ Bool_t AliTRDCalibraFit::AnalyseCH(AliTRDCalibraVector *calvect) } // Fill Infos Fit FillInfosFitCH(idect); - // Memory!!! - if (fDebugLevel != 1) { - delete projch; - } } // Boucle object // Normierungcharge if (fDebugLevel != 1) { @@ -536,7 +545,7 @@ Bool_t AliTRDCalibraFit::AnalyseCH(AliTRDCalibraVector *calvect) return kTRUE; } //________________functions fit Online PH2d____________________________________ -Bool_t AliTRDCalibraFit::AnalysePH(TProfile2D *ph) +Bool_t AliTRDCalibraFit::AnalysePH(const TProfile2D *ph) { // // Take the 1D profiles (average pulse height), projections of the 2D PH @@ -546,26 +555,38 @@ Bool_t AliTRDCalibraFit::AnalysePH(TProfile2D *ph) // // Set the calibration mode - const char *name = ph->GetTitle(); - 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) 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++) { + //printf("idect = %d\n",idect); // Determination of fNnZ, fNnRphi, fNfragZ and fNfragRphi....... UpdatefCountDetAndfCount(idect,1); ReconstructFitRowMinRowMax(idect,1); @@ -585,7 +606,7 @@ Bool_t AliTRDCalibraFit::AnalysePH(TProfile2D *ph) // This detector has not enough statistics or was off if (nentries <= fMinEntries) { //printf("Not enough statistic!\n"); - NotEnoughStatisticPH(idect); + NotEnoughStatisticPH(idect,nentries); if (fDebugLevel != 1) { delete projph; } @@ -598,6 +619,7 @@ Bool_t AliTRDCalibraFit::AnalysePH(TProfile2D *ph) CalculVdriftCoefMean(); CalculT0CoefMean(); //Method choosen + //printf("Method\n"); switch(fMethod) { case 0: FitLagrangePoly((TH1 *) projph); break; @@ -606,7 +628,7 @@ 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); + FillInfosFitPH(idect,nentries); // Memory!!! if (fDebugLevel != 1) { delete projph; @@ -635,8 +657,9 @@ Bool_t AliTRDCalibraFit::AnalysePH(AliTRDCalibraVector *calvect) // // Set the calibration mode - const char *name = calvect->GetNamePH(); - SetModeCalibration(name,1); + //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)) { @@ -657,24 +680,21 @@ Bool_t AliTRDCalibraFit::AnalysePH(AliTRDCalibraVector *calvect) 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); + if(!calvect->GetPHEntries(fCountDet)) { + NotEnoughStatisticPH(idect,fEntriesCurrent); + continue; } + 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 (fEntriesCurrent <= fMinEntries) { //printf("Not enough stat!\n"); - NotEnoughStatisticPH(idect); - if (fDebugLevel != 1) { - if(projph) delete projph; - } + NotEnoughStatisticPH(idect,fEntriesCurrent); continue; } // Statistic of the histos fitted @@ -692,13 +712,9 @@ 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); - // 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)); @@ -712,7 +728,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 +737,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 @@ -799,7 +816,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 +825,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(); @@ -903,8 +921,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 +947,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,10 +976,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) { @@ -986,11 +998,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 +1036,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,10 +1079,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) { @@ -1145,7 +1152,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]; } @@ -1168,7 +1175,7 @@ Bool_t AliTRDCalibraFit::AnalyseLinearFitters(AliTRDCalibraVdriftLinearFit *cali } //____________Functions for seeing if the pad is really okey___________________ //_____________________________________________________________________________ -Int_t AliTRDCalibraFit::GetNumberOfGroupsPRF(const char* nametitle) +Int_t AliTRDCalibraFit::GetNumberOfGroupsPRF(TString nametitle) { // // Get numberofgroupsprf @@ -1184,25 +1191,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; @@ -1210,7 +1217,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() @@ -1224,7 +1231,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 +1247,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 +1342,522 @@ 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))) 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 > -1.5) && (countAll > 15)) coef[(Int_t)(col*rowMax+row)] = meanAll+100.0; + if(ofwhat == 1){ + if((meanDetector[detector] > -1.5) && (countDetector[detector] > 20)) coef[(Int_t)(col*rowMax+row)] = meanDetector[detector]+100.0; + else if((meanSupermodule[sector] > -1.5) && (countSupermodule[sector] > 15)) coef[(Int_t)(col*rowMax+row)] = meanSupermodule[sector]+100.0; + else if((meanAll > -1.5) && (countAll > 15)) coef[(Int_t)(col*rowMax+row)] = meanAll+100.0; + } + } + // 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 +1899,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 @@ -1377,7 +1921,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,7 +1934,10 @@ 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 @@ -1400,7 +1950,7 @@ AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectGain(TObjArray *vectorFit, Double 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 +1970,11 @@ 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]; + // 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 +1982,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 +1996,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 @@ -1479,7 +2036,7 @@ AliTRDCalDet *AliTRDCalibraFit::CreateDetObjectLorentzAngle(TObjArray *vectorFit } //_____________________________________________________________________________ -TObject *AliTRDCalibraFit::CreatePadObjectGain(TObjArray *vectorFit, Double_t scaleFitFactor, AliTRDCalDet *detobject) +TObject *AliTRDCalibraFit::CreatePadObjectGain(const TObjArray *vectorFit, Double_t scaleFitFactor, const AliTRDCalDet *detobject) { // // It Creates the AliTRDCalPad object from AliTRDFitInfo @@ -1511,7 +2068,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 +2084,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 +2131,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 +2168,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 +2180,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 +2212,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 +2250,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 +2495,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; } @@ -2106,13 +2666,13 @@ Bool_t AliTRDCalibraFit::InitFitLinearFitter() 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)); + fCalDet2->SetValue(k,AliTRDCommonParam::Instance()->GetOmegaTau(fCalDet->GetValue(k))); } //printf("test3\n"); } else{ Float_t devalue = 1.5; - Float_t devalue2 = cal->GetOmegaTau(1.5,-fMagneticField); + Float_t devalue2 = AliTRDCommonParam::Instance()->GetOmegaTau(1.5); if(fCalDet) delete fCalDet; if(fCalDet2) delete fCalDet2; //printf("test1\n"); @@ -2152,10 +2712,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 +2739,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 +2822,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 +2842,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,7 +2862,89 @@ 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" @@ -2314,7 +2977,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,6 +2987,105 @@ 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="<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,6 +3135,92 @@ 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="<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); @@ -2448,26 +3296,110 @@ 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 +3407,124 @@ 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 +3537,107 @@ 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; } @@ -2560,7 +3685,7 @@ void AliTRDCalibraFit::FillFillCH(Int_t idect) if ( !fDebugStreamer ) { //debug stream TDirectory *backup = gDirectory; - fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root"); + fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root"); if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer } @@ -2589,7 +3714,7 @@ 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 @@ -2611,7 +3736,7 @@ void AliTRDCalibraFit::FillFillPH(Int_t idect) if ( !fDebugStreamer ) { //debug stream TDirectory *backup = gDirectory; - fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root"); + fDebugStreamer = new TTreeSRedirector("TRDDebugFitPH.root"); if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer } @@ -2633,6 +3758,7 @@ void AliTRDCalibraFit::FillFillPH(Int_t idect) (* fDebugStreamer) << "FillFillPH"<< "detector="<cd(); //we don't want to be cd'd to the debug streamer } @@ -2724,7 +3850,7 @@ void AliTRDCalibraFit::FillFillLinearFitter() if ( !fDebugStreamer ) { //debug stream TDirectory *backup = gDirectory; - fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root"); + fDebugStreamer = new TTreeSRedirector("TRDDebugFitLinearFitter.root"); if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer } @@ -2773,26 +3899,32 @@ Bool_t AliTRDCalibraFit::CalculT0CoefMean() fCurrentCoef2[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++) { 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 +3941,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 +3988,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 @@ -2930,18 +4065,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 +4174,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 +4205,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 +4239,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; @@ -3122,17 +4266,17 @@ void AliTRDCalibraFit::FitPente(TH1* projPH) if (fPhdt0 >= 0.0) { fCurrentCoef2[0] = (fPhdt0 - t0Shift) / widbins; if (fCurrentCoef2[0] < -1.0) { - fCurrentCoef2[0] = fCurrentCoef2[1]; + 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 +4309,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 +4340,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(); @@ -3234,7 +4384,7 @@ 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); + CalculPolynomeLagrange2(x,y,c0,c1,c2,c3,c4); AliInfo("At the limit for beginning!"); break; case 2: @@ -3248,7 +4398,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 +4414,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 +4427,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 +4442,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 +4451,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; @@ -3357,7 +4507,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 +4521,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 +4538,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 +4551,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 +4566,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 +4574,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; @@ -3457,8 +4607,20 @@ void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH) } //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!"); + put = kFALSE; + } + if((projPH->GetBinContent(binmin)+projPH->GetBinError(binmin)) > (projPH->GetBinContent(binmin-1))) { + AliInfo("Too many fluctuations at the end!"); + put = kFALSE; + } + if(TMath::Abs(pente->GetBinContent(binmin+1)) <= 0.0000000000001){ + AliInfo("No entries for the next bin!"); + pente->SetBinContent(binmin,0); + if(pente->GetEntries() > 0) binmin = (Int_t) pente->GetMinimumBin(); + } + x[0] = 0.0; x[1] = 0.0; @@ -3493,18 +4655,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"); + put = 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"); + put = 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 +4697,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 +4719,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 +4737,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"); + put = kFALSE; + } } //pol2 case 2 if(((binmin == (nbins-2)) && ((binmin-1) == TMath::Min(binmax+4, projPH->GetNbinsX()))) || @@ -3571,7 +4756,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,12 +4772,15 @@ 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 "); + put = kFALSE; + } } if((binmin == (nbins-1)) && ((binmin-2) < TMath::Min(binmax+4, projPH->GetNbinsX()))) { put = kFALSE; @@ -3605,14 +4793,14 @@ void AliTRDCalibraFit::FitLagrangePoly(TH1* projPH) AliInfo("For the drift...problem!"); } //pass but should not happen - if((binmin <= (nbins-3)) && (binmin < TMath::Min(binmax+4, projPH->GetNbinsX()))){ + if((binmin <= (nbins-3)) && (binmin < TMath::Min(binmax+6, projPH->GetNbinsX()))){ put = kFALSE; 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; @@ -3628,6 +4816,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; @@ -3645,21 +4836,38 @@ 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) { - fCurrentCoef2[0] = fCurrentCoef2[1]; + 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] < -1.0) { + fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0; + } + } + else { + fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0; + } + } + else{ + fCurrentCoef2[0] = fCurrentCoef2[1] + 100.0; + //printf("Fit failed!\n"); + } } if (fDebugLevel == 1) { @@ -3684,10 +4892,18 @@ 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; projPH->SetDirectory(0); @@ -3764,7 +4980,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 +4988,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 +5010,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 +5073,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++; @@ -3868,7 +5084,7 @@ Double_t AliTRDCalibraFit::FitGausMI(Double_t *arraye, Double_t *arraym, Double_ if ( !fDebugStreamer ) { //debug stream TDirectory *backup = gDirectory; - fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root"); + fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root"); if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer } @@ -3904,7 +5120,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); @@ -3921,7 +5137,7 @@ Double_t AliTRDCalibraFit::FitGausMI(Double_t *arraye, Double_t *arraym, Double_ if ( !fDebugStreamer ) { //debug stream TDirectory *backup = gDirectory; - fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root"); + fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root"); if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer } @@ -4050,7 +5266,7 @@ void AliTRDCalibraFit::FitTnpRange(Double_t *arraye, Double_t *arraym, Double_t if ( !fDebugStreamer ) { //debug stream TDirectory *backup = gDirectory; - fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root"); + fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root"); if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer } @@ -4104,7 +5320,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; @@ -4131,7 +5347,7 @@ void AliTRDCalibraFit::FitTnpRange(Double_t *arraye, Double_t *arraym, Double_t if ( !fDebugStreamer ) { //debug stream TDirectory *backup = gDirectory; - fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root"); + fDebugStreamer = new TTreeSRedirector("TRDDebugFitPRF.root"); if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer } @@ -4422,68 +5638,59 @@ void AliTRDCalibraFit::FitBisCH(TH1* projch, Double_t mean) } } //_____________________________________________________________________________ -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 +5698,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() @@ -4562,7 +5766,7 @@ void AliTRDCalibraFit::NormierungCharge() if ( !fDebugStreamer ) { //debug stream TDirectory *backup = gDirectory; - fDebugStreamer = new TTreeSRedirector("TRDDebugFit.root"); + fDebugStreamer = new TTreeSRedirector("TRDDebugFitCH.root"); if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer } (* fDebugStreamer) << "NormierungCharge"<< @@ -4571,7 +5775,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 +5804,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 +5830,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 +5891,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 +5946,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 +5972,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 +5982,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 +5998,7 @@ 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 @@ -4902,8 +6059,8 @@ Double_t AliTRDCalibraFit::LanGauFun(Double_t *x, Double_t *par) } //_____________________________________________________________________________ -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 { @@ -4941,7 +6098,7 @@ TF1 *AliTRDCalibraFit::LanGauFit(TH1 *his, Double_t *fitrange, Double_t *startva } //_____________________________________________________________________________ -Int_t AliTRDCalibraFit::LanGauPro(Double_t *params, Double_t &maxx, Double_t &fwhm) +Int_t AliTRDCalibraFit::LanGauPro(const Double_t *params, Double_t &maxx, Double_t &fwhm) { // // Function for the fit @@ -5037,7 +6194,7 @@ Int_t AliTRDCalibraFit::LanGauPro(Double_t *params, Double_t &maxx, Double_t &fw return (0); } //_____________________________________________________________________________ -Double_t AliTRDCalibraFit::GausConstant(Double_t *x, Double_t *par) +Double_t AliTRDCalibraFit::GausConstant(const Double_t *x, const Double_t *par) { // // Gaus with identical mean @@ -5048,3 +6205,6 @@ Double_t AliTRDCalibraFit::GausConstant(Double_t *x, Double_t *par) return gauss; } + + +