// //
///////////////////////////////////////////////////////////////////////////////
+#include <iostream>
+#include <fstream>
+#include <string>
+#include <TStyle.h>
+
#include "AliTRDCalROC.h"
+#include "TMath.h"
+#include "AliMathBase.h"
+#include "TLinearFitter.h"
+#include "TArrayI.h"
+#include "TH2F.h"
+#include "TH1F.h"
+#include "TArrayF.h"
+#include "TGraph2D.h"
+#include "TGraphDelaunay.h"
+#include "TList.h"
+
+#include "AliTRDCommonParam.h"
+#include "AliTRDpadPlane.h"
+#include "AliLog.h"
ClassImp(AliTRDCalROC)
Int_t iBin = 0;
- if (((AliTRDCalROC &) c).fData) delete [] ((AliTRDCalROC &) c).fData;
- ((AliTRDCalROC &) c).fData = new UShort_t[fNchannels];
+ fData = new UShort_t[fNchannels];
for (iBin = 0; iBin < fNchannels; iBin++) {
- ((AliTRDCalROC &) c).fData[iBin] = fData[iBin];
+ fData[iBin] = ((AliTRDCalROC &) c).fData[iBin];
}
}
for (iBin = 0; iBin < fNchannels; iBin++) {
((AliTRDCalROC &) c).fData[iBin] = fData[iBin];
}
-
+
TObject::Copy(c);
}
-//_____________________________________________________________________________
-void AliTRDCalROC::Scale(Float_t value)
+//___________________________________________________________________________________
+Double_t AliTRDCalROC::GetMean(AliTRDCalROC* outlierROC)
+{
+ //
+ // Calculate the mean
+ //
+
+ Double_t *ddata = new Double_t[fNchannels];
+ Int_t NPoints = 0;
+ for (Int_t i=0;i<fNchannels;i++) {
+ if ((!outlierROC) || (!(outlierROC->GetValue(i)))) {
+ ddata[NPoints]= (Double_t) fData[NPoints]/10000;
+ NPoints++;
+ }
+ }
+ Double_t mean = TMath::Mean(NPoints,ddata);
+ delete [] ddata;
+ return mean;
+}
+
+//_______________________________________________________________________________________
+Double_t AliTRDCalROC::GetMedian(AliTRDCalROC* outlierROC)
+{
+ //
+ // Calculate the median
+ //
+
+ Double_t *ddata = new Double_t[fNchannels];
+ Int_t NPoints = 0;
+ for (Int_t i=0;i<fNchannels;i++) {
+ if ((!outlierROC) || (!(outlierROC->GetValue(i)))) {
+ ddata[NPoints]= (Double_t) fData[NPoints]/10000;
+ NPoints++;
+ }
+ }
+ Double_t mean = TMath::Median(NPoints,ddata);
+ delete [] ddata;
+ return mean;
+}
+
+//____________________________________________________________________________________________
+Double_t AliTRDCalROC::GetRMS(AliTRDCalROC* outlierROC)
+{
+ //
+ // Calculate the RMS
+ //
+
+ Double_t *ddata = new Double_t[fNchannels];
+ Int_t NPoints = 0;
+ for (Int_t i=0;i<fNchannels;i++) {
+ if ((!outlierROC) || (!(outlierROC->GetValue(i)))) {
+ ddata[NPoints]= (Double_t)fData[NPoints]/10000;
+ NPoints++;
+ }
+ }
+ Double_t mean = TMath::RMS(NPoints,ddata);
+ delete [] ddata;
+ return mean;
+}
+
+//______________________________________________________________________________________________
+Double_t AliTRDCalROC::GetLTM(Double_t *sigma, Double_t fraction, AliTRDCalROC* outlierROC)
{
//
- // Scales all values of this ROC with the provided parameter. Is used if ROC defines
- // local variations of a global (or per detector defined) parameter
+ // Calculate LTM mean and sigma
//
- for (Int_t iBin = 0; iBin < fNchannels; iBin++) {
- fData[iBin] = (UShort_t) (value * fData[iBin]);
+ Double_t *ddata = new Double_t[fNchannels];
+ Double_t mean=0, lsigma=0;
+ UInt_t NPoints = 0;
+ for (Int_t i=0;i<fNchannels;i++) {
+ if (!outlierROC || !(outlierROC->GetValue(i))) {
+ ddata[NPoints]= (Double_t) fData[NPoints]/10000;
+ NPoints++;
+ }
}
+ Int_t hh = TMath::Min(TMath::Nint(fraction *NPoints), Int_t(NPoints));
+ AliMathBase::EvaluateUni(NPoints,ddata, mean, lsigma, hh);
+ if (sigma) *sigma=lsigma;
+ delete [] ddata;
+ return mean;
+}
+//___________________________________________________________________________________
+Bool_t AliTRDCalROC::Add(Float_t c1)
+{
+ //
+ // add constant
+ //
+ Bool_t result = kTRUE;
+ for (Int_t idata = 0; idata< fNchannels; idata++) {
+ if(((GetValue(idata)+c1) <= 6.5535) && ((GetValue(idata)+c1) >= 0.0)) SetValue(idata,GetValue(idata)+c1);
+ else {
+ SetValue(idata,0.0);
+ result = kFALSE;
+ }
+ }
+ return result;
+}
+
+//_______________________________________________________________________________________
+Bool_t AliTRDCalROC::Multiply(Float_t c1)
+{
+ //
+ // multiply constant
+ //
+ Bool_t result = kTRUE;
+ if(c1 < 0) return kFALSE;
+ for (Int_t idata = 0; idata< fNchannels; idata++) {
+ if((GetValue(idata)*c1) <= 6.5535) SetValue(idata,GetValue(idata)*c1);
+ else {
+ SetValue(idata,0.0);
+ result = kFALSE;
+ }
+ }
+ return result;
+}
+
+//____________________________________________________________________________________________
+Bool_t AliTRDCalROC::Add(const AliTRDCalROC * roc, Double_t c1)
+{
+ //
+ // add values
+ //
+ Bool_t result = kTRUE;
+ for (Int_t idata = 0; idata< fNchannels; idata++){
+ if(((GetValue(idata)+roc->GetValue(idata)*c1) <= 6.5535) &&
+ ((GetValue(idata)+roc->GetValue(idata)*c1) >= 0.0))
+ SetValue(idata,GetValue(idata)+roc->GetValue(idata)*c1);
+ else {
+ SetValue(idata,0.0);
+ result = kFALSE;
+ }
+ }
+ return result;
+}
+
+//____________________________________________________________________________________________
+Bool_t AliTRDCalROC::Multiply(const AliTRDCalROC* roc)
+{
+ //
+ // multiply values - per by pad
+ //
+ Bool_t result = kTRUE;
+ for (Int_t idata = 0; idata< fNchannels; idata++){
+ if((GetValue(idata)*roc->GetValue(idata)) <= 6.5535)
+ SetValue(idata,GetValue(idata)*roc->GetValue(idata));
+ else {
+ SetValue(idata,0.0);
+ result = kFALSE;
+ }
+ }
+ return result;
+}
+
+//______________________________________________________________________________________________
+Bool_t AliTRDCalROC::Divide(const AliTRDCalROC* roc)
+{
+ //
+ // divide values
+ //
+ Bool_t result = kTRUE;
+ Float_t kEpsilon=0.00000000000000001;
+ for (Int_t idata = 0; idata< fNchannels; idata++){
+ if (TMath::Abs(roc->GetValue(idata))>kEpsilon){
+ if((GetValue(idata)/roc->GetValue(idata)) <= 6.5535)
+ SetValue(idata,GetValue(idata)/roc->GetValue(idata));
+ else {
+ SetValue(idata,0.0);
+ result = kFALSE;
+ }
+ }
+ else result = kFALSE;
+ }
+ return result;
+}
+
+//__________________________________________________________________________________
+TH2F * AliTRDCalROC::MakeHisto2D(Float_t min, Float_t max,Int_t type)
+{
+ //
+ // make 2D histo
+ // type -1 = user defined range
+ // 0 = nsigma cut nsigma=min
+ // 1 = delta cut around median delta=min
+ Float_t kEpsilonr = 0.005;
+ gStyle->SetPalette(1);
+
+ if (type>=0){
+ if (type==0){
+ // nsigma range
+ Float_t mean = GetMean();
+ Float_t sigma = GetRMS();
+ Float_t nsigma = TMath::Abs(min);
+ sigma *= nsigma;
+ if(sigma < kEpsilonr) sigma = kEpsilonr;
+ min = mean-sigma;
+ max = mean+sigma;
+ }
+ if (type==1){
+ // fixed range
+ Float_t mean = GetMedian();
+ Float_t delta = min;
+ min = mean-delta;
+ max = mean+delta;
+ }
+ if (type==2){
+ Double_t sigma;
+ Float_t mean = GetLTM(&sigma,max);
+ sigma*=min;
+ if(sigma < kEpsilonr) sigma = kEpsilonr;
+ min = mean-sigma;
+ max = mean+sigma;
+
+ }
+ }
+
+ char name[1000];
+ sprintf(name,"%s 2D Plane %d Chamber %d",GetTitle(),fPla, fCha);
+ TH2F * his = new TH2F(name,name,fNrows,0, fNrows, fNcols, 0, fNcols);
+ for (Int_t irow=0; irow<fNrows; irow++){
+ for (Int_t icol=0; icol<fNcols; icol++){
+ his->Fill(irow+0.5,icol+0.5,GetValue(icol,irow));
+ }
+ }
+ his->SetStats(0);
+ his->SetMaximum(max);
+ his->SetMinimum(min);
+ return his;
+}
+
+//_______________________________________________________________________________________
+TH1F * AliTRDCalROC::MakeHisto1D(Float_t min, Float_t max,Int_t type)
+{
+ //
+ // make 1D histo
+ // type -1 = user defined range
+ // 0 = nsigma cut nsigma=min
+ // 1 = delta cut around median delta=min
+ Float_t kEpsilonr = 0.005;
+
+
+ if (type>=0){
+ if (type==0){
+ // nsigma range
+ Float_t mean = GetMean();
+ Float_t sigma = GetRMS();
+ Float_t nsigma = TMath::Abs(min);
+ sigma *= nsigma;
+ if(sigma < kEpsilonr) sigma = kEpsilonr;
+ min = mean-sigma;
+ max = mean+sigma;
+ }
+ if (type==1){
+ // fixed range
+ Float_t mean = GetMedian();
+ Float_t delta = min;
+ min = mean-delta;
+ max = mean+delta;
+ }
+ if (type==2){
+ //
+ // LTM mean +- nsigma
+ //
+ Double_t sigma;
+ Float_t mean = GetLTM(&sigma,max);
+ sigma*=min;
+ if(sigma < kEpsilonr) sigma = kEpsilonr;
+ min = mean-sigma;
+ max = mean+sigma;
+ }
+ }
+ char name[1000];
+ sprintf(name,"%s 1D Plane %d Chamber %d",GetTitle(),fPla, fCha);
+ TH1F * his = new TH1F(name,name,100, min,max);
+ for (Int_t irow=0; irow<fNrows; irow++){
+ for (Int_t icol=0; icol<fNcols; icol++){
+ his->Fill(GetValue(icol,irow));
+ }
+ }
+ return his;
}