]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TRD/Cal/AliTRDCalROC.cxx
Making online tracklets usable in offline reconstruction
[u/mrichter/AliRoot.git] / TRD / Cal / AliTRDCalROC.cxx
index 6530c59d3291a32ccd26bc364dc9a89c9eff78c2..5320731038b60e37bbf6d904286f5d877b5d4ba4 100644 (file)
 //                                                                           //
 ///////////////////////////////////////////////////////////////////////////////
 
+#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)
 
@@ -151,10 +170,9 @@ AliTRDCalROC::AliTRDCalROC(const AliTRDCalROC &c)
 
   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];
   }
 
 }
@@ -207,21 +225,345 @@ void AliTRDCalROC::Copy(TObject &c) const
   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)))) {
+       if(fData[i] > 0.000000000000001){
+         ddata[NPoints]= (Double_t) fData[i]/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)))) {
+        if(fData[i] > 0.000000000000001){         
+          ddata[NPoints]= (Double_t) fData[i]/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)))) {
+       if(fData[i] > 0.000000000000001){
+         ddata[NPoints]= (Double_t)fData[i]/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)
+{
+  //
+  //  Calculate LTM mean and sigma
+  //
+
+  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))) {
+       if(fData[i] > 0.000000000000001){
+        ddata[NPoints]= (Double_t) fData[i]/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)
 {
   //
-  // 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
+  // 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;
+}
 
-  for (Int_t iBin = 0; iBin < fNchannels; iBin++) {
-    fData[iBin] = (UShort_t) (value * fData[iBin]);
+//_______________________________________________________________________________________
+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;
+}
+//______________________________________________________________________________________________
+Bool_t AliTRDCalROC::Unfold() 
+{
+  //
+  // Compute the mean value per pad col
+  // Divide with this value each pad col
+  // This is for the noise study
+  // Return kFALSE if one or more of the pad col was not normalised
+  //
+  Bool_t result = kTRUE;
+  Float_t kEpsilon=0.00000000000000001;
+  
+  // calcul the mean value per col
+  for(Int_t icol = 0; icol < fNcols; icol++){
+
+    Float_t mean = 0.0;
+    Float_t nb   = 0.0;
+
+    for(Int_t irow = 0; irow < fNrows; irow++){
+      if((GetValue(icol,irow) > 0.06) && (GetValue(icol,irow) < 0.15)){
+       mean += GetValue(icol,irow);
+       nb += 1.0;
+      }
+    }
+    
+    if(nb > kEpsilon) {
+      
+      mean = mean/nb;
+
+      if(mean > kEpsilon){
+       for(Int_t irow = 0; irow < fNrows; irow++){
+         Float_t value = GetValue(icol,irow);
+         SetValue(icol,irow,(Float_t)(value/mean));
+       }    
+      }
+      else result = kFALSE;
+      
+    }
+
+    else result = kFALSE;
+
+    
+  }
+
+  return result;
+}
+
+//__________________________________________________________________________________
+TH2F * AliTRDCalROC::MakeHisto2D(Float_t min, Float_t max,Int_t type,  Float_t mu)
+{
+  //
+  // 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)*mu);
+    }
+  }
+  his->SetStats(0);
+  his->SetMaximum(max);
+  his->SetMinimum(min);
+  return his;
+}
+
+//_______________________________________________________________________________________
+TH1F * AliTRDCalROC::MakeHisto1D(Float_t min, Float_t max,Int_t type,  Float_t mu)
+{
+  //
+  // 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)*mu);
+    }
+  }
+  return his;
 }