]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
Add the calibration object algebra by Raphaelle
authorcblume <cblume@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 11 May 2007 11:11:17 +0000 (11:11 +0000)
committercblume <cblume@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 11 May 2007 11:11:17 +0000 (11:11 +0000)
TRD/Cal/AliTRDCalDet.cxx
TRD/Cal/AliTRDCalDet.h
TRD/Cal/AliTRDCalPad.cxx
TRD/Cal/AliTRDCalPad.h
TRD/Cal/AliTRDCalROC.cxx
TRD/Cal/AliTRDCalROC.h

index c78b56b1a595f9cc65fd1dc9d66be325dea29a46..47d8de8ee1797a829bda4a29703d22828e3840fe 100644 (file)
 //                                                                           //
 ///////////////////////////////////////////////////////////////////////////////
 
+#include <TMath.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TStyle.h>
+
 #include "AliTRDCalDet.h"
+#include "AliTRDgeometry.h"
+#include "AliMathBase.h"
+#include "AliTRDpadPlane.h"
 
 ClassImp(AliTRDCalDet)
 
@@ -99,3 +107,385 @@ void AliTRDCalDet::Copy(TObject &c) const
 
 }
 
+//___________________________________________________________________________________
+Double_t AliTRDCalDet::GetMean(AliTRDCalDet* outlierDet) 
+{
+  //
+  // Calculate the mean
+  //
+
+  if (!outlierDet) return TMath::Mean(kNdet,fData);
+  Double_t *ddata = new Double_t[kNdet];
+  Int_t NPoints = 0;
+  for (Int_t i=0;i<kNdet;i++) {
+    if (!(outlierDet->GetValue(i))) {
+      ddata[NPoints]= fData[NPoints];
+      NPoints++;
+    }
+  }
+  Double_t mean = TMath::Mean(NPoints,ddata);
+  delete [] ddata;
+  return mean;
+}
+
+//_______________________________________________________________________________________
+Double_t AliTRDCalDet::GetMedian(AliTRDCalDet* outlierDet) 
+{
+  //
+  // Calculate the median
+  //
+
+  if (!outlierDet) return (Double_t) TMath::Median(kNdet,fData);
+  Double_t *ddata = new Double_t[kNdet];
+  Int_t NPoints = 0;
+  for (Int_t i=0;i<kNdet;i++) {
+    if (!(outlierDet->GetValue(i))) {
+      ddata[NPoints]= fData[NPoints];
+      NPoints++;
+    }
+  }
+  Double_t mean = TMath::Median(NPoints,ddata);
+  delete [] ddata;
+  return mean;
+
+}
+
+//____________________________________________________________________________________________
+Double_t AliTRDCalDet::GetRMS(AliTRDCalDet* outlierDet) 
+{
+  //
+  // Calculate the RMS
+  //
+
+  if (!outlierDet) return TMath::RMS(kNdet,fData);
+  Double_t *ddata = new Double_t[kNdet];
+  Int_t NPoints = 0;
+  for (Int_t i=0;i<kNdet;i++) {
+    if (!(outlierDet->GetValue(i))) {
+      ddata[NPoints]= fData[NPoints];
+      NPoints++;
+    }
+  }
+  Double_t mean = TMath::RMS(NPoints,ddata);
+  delete [] ddata;
+  return mean;
+}
+
+//______________________________________________________________________________________________
+Double_t AliTRDCalDet::GetLTM(Double_t *sigma, Double_t fraction, AliTRDCalDet* outlierDet)
+{
+  //
+  //  Calculate LTM mean and sigma
+  //
+
+  Double_t *ddata = new Double_t[kNdet];
+  Double_t mean=0, lsigma=0;
+  UInt_t NPoints = 0;
+  for (Int_t i=0;i<kNdet;i++) {
+     if (!outlierDet || !(outlierDet->GetValue(i))) {
+        ddata[NPoints]= fData[NPoints];
+        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;
+}
+
+//_________________________________________________________________________
+TH1F * AliTRDCalDet::MakeHisto1Distribution(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
+  //
+
+  if (type>=0){
+    if (type==0){
+      // nsigma range
+      Float_t mean  = GetMean();
+      Float_t sigma = GetRMS();
+      Float_t nsigma = TMath::Abs(min);
+      min = mean-nsigma*sigma;
+      max = mean+nsigma*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;
+      min = mean-sigma;
+      max = mean+sigma;
+    }
+  }
+  char  name[1000];
+  sprintf(name,"%s CalDet 1Distribution",GetTitle());
+  TH1F * his = new TH1F(name,name,100, min,max);
+  for (Int_t idet=0; idet<kNdet; idet++){
+    his->Fill(GetValue(idet));
+  }
+  return his;
+}
+
+//________________________________________________________________________________
+TH1F * AliTRDCalDet::MakeHisto1DAsFunctionOfDet(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
+  //
+
+  if (type>=0){
+    if (type==0){
+      // nsigma range
+      Float_t mean  = GetMean();
+      Float_t sigma = GetRMS();
+      Float_t nsigma = TMath::Abs(min);
+      min = mean-nsigma*sigma;
+      max = mean+nsigma*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;
+      min = mean-sigma;
+      max = mean+sigma;
+
+    }
+  }
+  char  name[1000];
+  sprintf(name,"%s CalDet as function of det",GetTitle());
+  TH1F * his = new TH1F(name,name,kNdet, 0, kNdet);
+  for(Int_t det = 0; det< kNdet; det++){
+    his->Fill(det+0.5,GetValue(det));
+  }
+
+  his->SetMaximum(max);
+  his->SetMinimum(min);
+  return his;
+}
+
+//_____________________________________________________________________________
+TH2F *AliTRDCalDet::MakeHisto2DCh(Int_t ch, Float_t min, Float_t max, Int_t type)
+{
+  //
+  // Make 2D graph
+  // ch    - chamber
+  // type -1 = user defined range
+  //       0 = nsigma cut nsigma=min
+  //       1 = delta cut around median delta=min
+  //
+
+  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);
+      min = mean-nsigma*sigma;
+      max = mean+nsigma*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;
+      min = mean-sigma;
+      max = mean+sigma;
+
+    }
+  }
+    
+  AliTRDgeometry *TRDgeo = new AliTRDgeometry();
+  TRDgeo->Init();
+      
+
+  Double_t poslocal[3]  = {0.0,0.0,0.0};
+  Double_t posglobal[3] = {0.0,0.0,0.0};
+  
+  char  name[1000];
+  sprintf(name,"%s CalDet 2D ch %d",GetTitle(),ch);
+  TH2F * his = new TH2F(name, name, 400,-400.0,400.0,400,-400.0,400.0);
+
+  // Where we begin
+  Int_t offsetch = 6*ch;
+  
+
+  for (Int_t isec = 0; isec < kNsect; isec++){
+    for(Int_t ipl = 0; ipl < kNplan; ipl++){
+      Int_t det   = offsetch+isec*30+ipl;
+      AliTRDpadPlane *padPlane = new AliTRDpadPlane(ipl,ch);
+      for (Int_t icol=0; icol<padPlane->GetNcols(); icol++){
+       poslocal[0] = TRDgeo->GetTime0(ipl);
+       poslocal[2] = padPlane->GetRowPos(0);
+       poslocal[1] = padPlane->GetColPos(icol);
+       TRDgeo->RotateBack(det,poslocal,posglobal);
+       Int_t binx = 1+TMath::Nint((posglobal[0]+400.0)*0.5);
+       Int_t biny = 1+TMath::Nint((posglobal[1]+400.0)*0.5);
+       his->SetBinContent(binx,biny,fData[det]);
+      }
+    }    
+  }
+  his->SetXTitle("x (cm)");
+  his->SetYTitle("y (cm)");
+  his->SetStats(0);
+  his->SetMaximum(max);
+  his->SetMinimum(min);
+  return his;
+}
+
+//_____________________________________________________________________________
+TH2F *AliTRDCalDet::MakeHisto2DSmPl(Int_t sm, Int_t pl, Float_t min, Float_t max, Int_t type)
+{
+  //
+  // Make 2D graph
+  // sm    - supermodule number
+  // pl    - plane number
+  // type -1 = user defined range
+  //       0 = nsigma cut nsigma=min
+  //       1 = delta cut around median delta=min
+  //
+
+  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);
+      min = mean-nsigma*sigma;
+      max = mean+nsigma*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;
+      min = mean-sigma;
+      max = mean+sigma;
+
+    }
+  }
+     
+  AliTRDpadPlane *padPlane0 = new AliTRDpadPlane(pl,0);
+  Double_t row0    = padPlane0->GetRow0();
+  Double_t col0    = padPlane0->GetCol0();
+
+
+  char  name[1000];
+  sprintf(name,"%s CalDet 2D sm %d and pl %d",GetTitle(),sm,pl);
+  TH2F * his = new TH2F( name, name, 5,  -TMath::Abs(row0),  TMath::Abs(row0)
+                                   , 4,-2*TMath::Abs(col0),2*TMath::Abs(col0));
+
+  // Where we begin
+  Int_t offsetsmpl = 30*sm+pl;
+  
+
+  for (Int_t k = 0; k < kNcham; k++){
+    Int_t det = offsetsmpl+k*6;
+    Int_t kb  = kNcham-1-k;
+    his->SetBinContent(kb+1,2,fData[det]);
+    his->SetBinContent(kb+1,3,fData[det]);
+  }
+  his->SetXTitle("z (cm)");
+  his->SetYTitle("y (cm)");
+  his->SetStats(0);
+  his->SetMaximum(max);
+  his->SetMinimum(min);
+  return his;
+}
+
+//_____________________________________________________________________________
+void AliTRDCalDet::Add(Float_t c1)
+{
+  //
+  // Add constant for all detectors
+  //
+
+  for (Int_t idet = 0; idet < kNdet; idet++) {
+     fData[idet] += c1;
+  }
+}
+
+//_____________________________________________________________________________
+void AliTRDCalDet::Multiply(Float_t c1)
+{
+    //
+    // multiply constant for all detectors
+    //
+    for (Int_t idet = 0; idet < kNdet; idet++) {
+      fData[idet] *= c1;
+    }
+}
+
+//_____________________________________________________________________________
+void AliTRDCalDet::Add(const AliTRDCalDet * calDet, Double_t c1)
+{
+    //
+    // add caldet channel by channel multiplied by c1
+    //
+    for (Int_t idet = 0; idet < kNdet; idet++) {
+      fData[idet] += calDet->GetValue(idet)*c1;
+    }
+}
+
+//_____________________________________________________________________________
+void AliTRDCalDet::Multiply(const AliTRDCalDet * calDet)
+{
+    //
+    // multiply caldet channel by channel 
+    //
+    for (Int_t idet = 0; idet < kNdet; idet++) {
+      fData[idet] *= calDet->GetValue(idet);
+    }
+}
+
+//_____________________________________________________________________________
+void AliTRDCalDet::Divide(const AliTRDCalDet * calDet)
+{
+    //
+    // divide caldet channel by channel 
+    //
+ Float_t kEpsilon=0.00000000000000001;
+    for (Int_t idet = 0; idet < kNdet; idet++) {
+      if (TMath::Abs(calDet->GetValue(idet))>kEpsilon){
+       fData[idet] /= calDet->GetValue(idet);
+       }
+    }
+}
index 6164bcb77b5f4c5b18c375524c77990546923b5e..a1542d0b65d14b7b3526cfba1f4f2cb13a435a38 100644 (file)
 ///////////////////////////////////////////////////////////////////////////////
 
 #include "TNamed.h"
-#include "AliTRDgeometry.h"
+#include "TMath.h"
+#include "../AliTRDgeometry.h"
+
+class TH1F;
+class TH2F;
 
 class AliTRDCalDet : public TNamed {
 
@@ -35,7 +39,26 @@ class AliTRDCalDet : public TNamed {
   void          SetValue(Int_t d, Float_t value) { fData[d] = value; };
   void          SetValue(Int_t p, Int_t c, Int_t s, Float_t value) 
                                                  { fData[AliTRDgeometry::GetDetector(p,c,s)] = value; };
-  
+
+  // statistic
+  Double_t GetMean(AliTRDCalDet *outlierDet=0);
+  Double_t GetRMS(AliTRDCalDet *outlierDet=0);
+  Double_t GetMedian(AliTRDCalDet *outlierDet=0);
+  Double_t GetLTM(Double_t *sigma=0, Double_t fraction=0.9, AliTRDCalDet *outlierDet=0);     
+
+  // Plot functions
+  TH1F * MakeHisto1Distribution(Float_t min=4, Float_t max=-4, Int_t type=0);     
+  TH1F * MakeHisto1DAsFunctionOfDet(Float_t min=4, Float_t max=-4, Int_t type=0);
+  TH2F * MakeHisto2DCh(Int_t ch, Float_t min=4, Float_t max=-4, Int_t type=0);  
+  TH2F * MakeHisto2DSmPl(Int_t sm, Int_t pl, Float_t min=4, Float_t max=-4, Int_t type=0); 
+
+  // algebra functions
+  void Add(Float_t c1);
+  void Multiply(Float_t c1);
+  void Add(const AliTRDCalDet * calDet, Double_t c1 = 1);
+  void Multiply(const AliTRDCalDet * calDet);
+  void Divide(const AliTRDCalDet * calDet);
+    
  protected:
 
   Float_t  fData[kNdet];       //[kNdet] Data
index e1d2ff5e6bab3054056082a57c2af6d113abbcec..a9c9b820c0f0888c487590f8bed9026fcaace70f 100644 (file)
 //                                                                           //
 ///////////////////////////////////////////////////////////////////////////////
 
+#include <TMath.h>
+#include <TH2F.h>
+#include <TH1F.h>
+#include <TStyle.h>
+
 #include "AliTRDCalPad.h"
 #include "AliTRDCalROC.h"
 #include "AliTRDCalDet.h"
+#include "AliTRDpadPlane.h"
+#include "AliMathBase.h"
+#include "AliTRDgeometry.h"
 
 ClassImp(AliTRDCalPad)
 
@@ -60,17 +68,20 @@ AliTRDCalPad::AliTRDCalPad(const Text_t *name, const Text_t *title)
 }
 
 //_____________________________________________________________________________
-AliTRDCalPad::AliTRDCalPad(const AliTRDCalPad &c):TNamed(c)
+AliTRDCalPad::AliTRDCalPad(const AliTRDCalPad &c)
+  :TNamed(c)
 {
   //
   // AliTRDCalPad copy constructor
   //
 
-  ((AliTRDCalPad &) c).Copy(*this);
+  for (Int_t idet = 0; idet < kNdet; idet++) {
+    fROC[idet] = new AliTRDCalROC(*((AliTRDCalPad &) c).fROC[idet]);
+  }
 
 }
 
-///_____________________________________________________________________________
+//_____________________________________________________________________________
 AliTRDCalPad::~AliTRDCalPad()
 {
   //
@@ -116,7 +127,7 @@ void AliTRDCalPad::Copy(TObject &c) const
 }
 
 //_____________________________________________________________________________
-void AliTRDCalPad::ScaleROCs(AliTRDCalDet* values)
+Bool_t AliTRDCalPad::ScaleROCs(const AliTRDCalDet* values)
 {
   // 
   // Scales ROCs of this class with the values from the class <values>
@@ -125,13 +136,593 @@ void AliTRDCalPad::ScaleROCs(AliTRDCalDet* values)
   //
   
   if (!values)
-    return;
-  
+    return kFALSE;
+  Bool_t result = kTRUE;
   for (Int_t idet = 0; idet < kNdet; idet++) {
     if (fROC[idet]) { 
-      fROC[idet]->Scale(values->GetValue(idet));
+      if(!fROC[idet]->Multiply(values->GetValue(idet))) result = kFALSE;
+    }
+  }
+  return result;
+}
+
+//_____________________________________________________________________________
+Double_t AliTRDCalPad::GetMeanRMS(Double_t &rms, const AliTRDCalDet *calDet, Int_t type)
+{
+    //
+    // Calculate mean an RMS of all rocs
+    // If calDet correct the CalROC from the detector coefficient
+    // type == 0 for gain and vdrift
+    // type == 1 for t0 
+    //
+  Double_t factor = 0.0;
+  if(type == 0) factor = 1.0;
+  Double_t sum = 0, sum2 = 0, n=0, val=0;
+  for (Int_t idet = 0; idet < kNdet; idet++) {
+    if(calDet) factor = calDet->GetValue(idet);
+    AliTRDCalROC *calRoc = fROC[idet];
+    if ( calRoc ){
+      for (Int_t irow=0; irow<calRoc->GetNrows(); irow++){
+       for (Int_t icol=0; icol<calRoc->GetNcols(); icol++){
+         if(type == 0) val = calRoc->GetValue(icol,irow)*factor;
+         else  val = calRoc->GetValue(icol,irow)+factor;
+         sum+=val;
+         sum2+=val*val;
+         n++;
+       }
+      }
+    }
+  }
+  Double_t n1 = 1./n;
+  Double_t mean = sum*n1;
+  rms  = TMath::Sqrt(TMath::Abs(sum2*n1-mean*mean));
+  return mean;
+}
+
+//_____________________________________________________________________________
+Double_t AliTRDCalPad::GetMean(const AliTRDCalDet *calDet, Int_t type, AliTRDCalPad* outlierPad)
+{
+    //
+    // return mean of the mean of all ROCs
+    // If calDet correct the CalROC from the detector coefficient
+    // type == 0 for gain and vdrift
+    // type == 1 for t0 
+    //
+  Double_t factor = 0.0;
+  if(type == 0) factor = 1.0;
+  Double_t arr[kNdet];
+  Int_t n=0;
+  for (Int_t idet = 0; idet < kNdet; idet++) {
+    if(calDet) factor = calDet->GetValue(idet);
+    AliTRDCalROC *calRoc = fROC[idet];
+    if ( calRoc ){
+      AliTRDCalROC* outlierROC = 0;
+      if (outlierPad) outlierROC = outlierPad->GetCalROC(idet);
+      if(type == 0)  arr[n] = calRoc->GetMean(outlierROC)*factor;
+      else arr[n] = calRoc->GetMean(outlierROC)+factor;
+      n++;
+    }
+  }
+  return TMath::Mean(n,arr);
+}
+
+//_____________________________________________________________________________
+Double_t AliTRDCalPad::GetRMS(const AliTRDCalDet *calDet, Int_t type, AliTRDCalPad* outlierPad)
+{
+    //
+    // return mean of the RMS of all ROCs
+    // If calDet correct the CalROC from the detector coefficient
+    // type == 0 for gain and vdrift
+    // type == 1 for t0 
+    //
+  Double_t factor = 0.0;
+  if(type == 0) factor = 1.0;
+  Double_t arr[kNdet];
+  Int_t n=0;
+  for (Int_t idet = 0; idet < kNdet; idet++) {
+    if(calDet) factor = calDet->GetValue(idet);
+    AliTRDCalROC *calRoc = fROC[idet];
+    if ( calRoc ){
+      AliTRDCalROC* outlierROC = 0;
+      if (outlierPad) outlierROC = outlierPad->GetCalROC(idet);
+      if(type == 0)  arr[n] = calRoc->GetRMS(outlierROC)*factor;
+      else  arr[n] = calRoc->GetRMS(outlierROC);
+      n++;
+    }
+  }
+    return TMath::Mean(n,arr);
+}
+
+//_____________________________________________________________________________
+Double_t AliTRDCalPad::GetMedian(const AliTRDCalDet *calDet, Int_t type
+                               , AliTRDCalPad* outlierPad)
+{
+    //
+    // return mean of the median of all ROCs
+    // If AliTRDCalDet, the correct the CalROC from the detector coefficient
+    // type == 0 for gain and vdrift
+    // type == 1 for t0 
+    //
+  Double_t factor = 0.0;
+  if(type == 0) factor = 1.0;
+  Double_t arr[kNdet];
+  Int_t n=0;
+  for (Int_t idet = 0; idet < kNdet; idet++) {
+    if(calDet) factor = calDet->GetValue(idet);
+    AliTRDCalROC *calRoc = fROC[idet];
+    if ( calRoc ){
+      AliTRDCalROC* outlierROC = 0;
+      if (outlierPad) outlierROC = outlierPad->GetCalROC(idet);
+      if(type == 0)  arr[n] = calRoc->GetMedian(outlierROC)*factor;
+      else  arr[n] = calRoc->GetMedian(outlierROC)+factor;
+      n++;
+    }
+  }
+  return TMath::Mean(n,arr);
+}
+
+//_____________________________________________________________________________
+Double_t AliTRDCalPad::GetLTM(Double_t *sigma, Double_t fraction
+                            , const AliTRDCalDet *calDet, Int_t type
+                            , AliTRDCalPad* outlierPad)
+{
+    //
+    // return mean of the LTM and sigma of all ROCs
+    // If calDet correct the CalROC from the detector coefficient
+    // type == 0 for gain and vdrift
+    // type == 1 for t0 
+    //
+  Double_t factor = 0.0;
+  if(type == 0) factor = 1.0;
+  Double_t arrm[kNdet];
+  Double_t arrs[kNdet];
+  Double_t *sTemp=0x0;
+  Int_t n=0;
+  
+  for (Int_t idet = 0; idet < kNdet; idet++) {
+    if(calDet) factor = calDet->GetValue(idet);
+    AliTRDCalROC *calRoc = fROC[idet];
+    if ( calRoc ){
+      if ( sigma ) sTemp=arrs+n;
+      AliTRDCalROC* outlierROC = 0;
+      if (outlierPad) outlierROC = outlierPad->GetCalROC(idet);
+      if(type == 0)  arrm[n] = calRoc->GetLTM(sTemp,fraction, outlierROC)*factor;
+      else arrm[n] = calRoc->GetLTM(sTemp,fraction, outlierROC)+factor;
+      n++;
+    }
+  }
+  if ( sigma ) *sigma = TMath::Mean(n,arrs);
+  return TMath::Mean(n,arrm);
+}
+
+//_____________________________________________________________________________
+TH1F * AliTRDCalPad::MakeHisto1D(const AliTRDCalDet *calDet, Int_t typedet
+                               , Float_t min, Float_t max,Int_t type)
+{
+  //
+  // make 1D histo
+  // type -1 = user defined range
+  //       0 = nsigma cut nsigma=min
+  // If calDet correct the CalROC from the detector coefficient
+  // typedet == 0 for gain and vdrift
+  // typedet == 1 for t0
+  //
+  Float_t kEpsilonr = 0.005;
+
+  Double_t factor = 0.0;
+  if(typedet == 0) factor = 1.0;
+  if (type>=0){
+    if (type==0){
+      // nsigma range 
+      Float_t mean  = GetMean(calDet,typedet);
+      Float_t sigma = 0.0;
+      if(GetRMS(calDet,typedet) > kEpsilonr) sigma = GetRMS(calDet,typedet);
+      else {
+       Double_t rms = 0.0;
+       sigma = GetMeanRMS(rms,calDet,typedet);
+      }
+      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(calDet,typedet);
+      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,calDet,typedet);
+      sigma*=min;
+      if(sigma < kEpsilonr) sigma = kEpsilonr;
+      min = mean-sigma;
+      max = mean+sigma;
+    }
+  }
+  char  name[1000];
+  sprintf(name,"%s Pad 1D",GetTitle());
+  TH1F * his = new TH1F(name,name,100, min,max);
+    for (Int_t idet = 0; idet < kNdet; idet++) {
+      if(calDet) factor = calDet->GetValue(idet);
+      if (fROC[idet]){
+       for (Int_t irow=0; irow<fROC[idet]->GetNrows(); irow++){
+         for (Int_t icol=0; icol<fROC[idet]->GetNcols(); icol++){
+           if(typedet == 0) his->Fill(fROC[idet]->GetValue(irow,icol)*factor);
+           else his->Fill(fROC[idet]->GetValue(irow,icol)+factor);
+         }
+       }
+      }
+    }
+    return his;
+}
+
+//_____________________________________________________________________________
+TH2F *AliTRDCalPad::MakeHisto2DSmPl(Int_t sm, Int_t pl, const AliTRDCalDet *calDet
+                                  , Int_t typedet, Float_t min, Float_t max,Int_t type)
+{
+  //
+  // Make 2D graph
+  // sm    - supermodule number
+  // pl    - plane number
+  // If calDet correct the CalROC from the detector coefficient
+  // typedet == 0 for gain and vdrift
+  // typedet == 1 for t0
+  //
+  gStyle->SetPalette(1);
+  Double_t factor = 0.0;
+  if(typedet == 0) factor = 1.0;
+
+  Float_t kEpsilon = 0.000000000001;
+  Float_t kEpsilonr = 0.005;
+
+  if (type>=0){
+    if (type==0){
+      // nsigma range
+      Float_t mean  = GetMean(calDet,typedet);
+      Float_t sigma = 0.0;
+      if(GetRMS(calDet,typedet) > kEpsilonr) sigma = GetRMS(calDet,typedet);
+      else {
+       Double_t rms = 0.0;
+       sigma = GetMeanRMS(rms,calDet,typedet);
+      }
+      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(calDet,typedet);
+      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,calDet,typedet);
+      sigma*=min;
+      if(sigma < kEpsilonr) sigma = kEpsilonr;
+      min = mean-sigma;
+      max = mean+sigma;
+    }
+  }
+  
+  AliTRDpadPlane *padPlane0 = new AliTRDpadPlane(pl,0);
+  Double_t row0    = padPlane0->GetRow0();
+  Double_t col0    = padPlane0->GetCol0();
+
+  char  name[1000];
+  sprintf(name,"%s Pad 2D sm %d pl %d",GetTitle(),sm,pl);
+  TH2F * his = new TH2F( name, name, 76,-TMath::Abs(row0),TMath::Abs(row0),144,-TMath::Abs(col0),TMath::Abs(col0));
+
+  // Where we begin
+  Int_t offsetsmpl = 30*sm+pl;
+  
+
+  for (Int_t k = 0; k < kNcham; k++){
+    Int_t det = offsetsmpl+k*6;
+    if(calDet) factor = calDet->GetValue(det);
+    if (fROC[det]){
+      AliTRDCalROC * calRoc = fROC[det];
+      for (Int_t irow=0; irow<calRoc->GetNrows(); irow++){
+       for (Int_t icol=0; icol<calRoc->GetNcols(); icol++){
+         if (TMath::Abs(calRoc->GetValue(icol,irow))>kEpsilon){
+           Int_t binz     = 0;
+           Int_t kb       = kNcham-1-k;
+           if(kb > 2) binz = 16*(kb-1)+12+irow+1;
+           else binz = 16*kb+irow+1; 
+           Int_t biny = icol+1;
+           Float_t value = calRoc->GetValue(icol,irow);
+           if(typedet == 0) his->SetBinContent(binz,biny,value*factor);
+           else his->SetBinContent(binz,biny,value+factor);
+         }
+       }
+      }
+    }
+  }
+  his->SetXTitle("z (cm)");
+  his->SetYTitle("y (cm)");
+  his->SetStats(0);
+  his->SetMaximum(max);
+  his->SetMinimum(min);
+  return his;
+}
+
+//_____________________________________________________________________________
+TH2F *AliTRDCalPad::MakeHisto2DCh(Int_t ch, const AliTRDCalDet *calDet, Int_t typedet
+                                , Float_t min, Float_t max,Int_t type)
+{
+  //
+  // Make 2D graph mean value in z direction
+  // ch    - chamber
+  // If calDet correct the CalROC from the detector coefficient
+  // typedet == 0 for gain and vdrift
+  // typedet == 1 for t0
+  //
+  gStyle->SetPalette(1);
+  Double_t factor = 0.0;
+  if(typedet == 0) factor = 1.0;
+
+  Float_t kEpsilonr = 0.005;
+  
+  if (type>=0){
+    if (type==0){
+      // nsigma range
+      Float_t mean  = GetMean(calDet,typedet);
+      Float_t sigma = 0.0;
+      if(GetRMS(calDet,typedet) > kEpsilonr) sigma = GetRMS(calDet,typedet);
+      else {
+       Double_t rms = 0.0;
+       sigma = GetMeanRMS(rms,calDet,typedet);
+      }
+      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(calDet,typedet);
+      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,calDet,typedet);
+      sigma*=min;
+      if(sigma < kEpsilonr) sigma = kEpsilonr;
+      min = mean-sigma;
+      max = mean+sigma;
+    }
+  }
+
+  AliTRDgeometry *TRDgeo = new AliTRDgeometry();
+  TRDgeo->Init();
+      
+  Float_t kEpsilon = 0.000000000001;
+
+  Double_t poslocal[3]  = {0.0,0.0,0.0};
+  Double_t posglobal[3] = {0.0,0.0,0.0};
+  
+  char  name[1000];
+  sprintf(name,"%s Pad 2D ch %d",GetTitle(),ch);
+  TH2F * his = new TH2F( name, name, 400,-400.0,400.0,400,-400.0,400.0);
+
+  // Where we begin
+  Int_t offsetch = 6*ch;
+  
+
+  for (Int_t isec = 0; isec < kNsect; isec++){
+    for(Int_t ipl = 0; ipl < kNplan; ipl++){
+      Int_t det = offsetch+isec*30+ipl;
+      if(calDet) factor = calDet->GetValue(det);
+      if (fROC[det]){
+       AliTRDCalROC * calRoc = fROC[det];
+       AliTRDpadPlane *padPlane = new AliTRDpadPlane(ipl,ch);
+       for (Int_t icol=0; icol<calRoc->GetNcols(); icol++){
+         poslocal[0] = TRDgeo->GetTime0(ipl);
+         poslocal[2] = padPlane->GetRowPos(0);
+         poslocal[1] = padPlane->GetColPos(icol);
+         TRDgeo->RotateBack(det,poslocal,posglobal);
+         Int_t binx = 1+TMath::Nint((posglobal[0]+400.0)*0.5);
+         Int_t biny = 1+TMath::Nint((posglobal[1]+400.0)*0.5);
+         Float_t value = 0.0;
+         Int_t nb = 0;
+         for (Int_t irow=0; irow<calRoc->GetNrows(); irow++){
+           if (TMath::Abs(calRoc->GetValue(icol,irow))>kEpsilon){
+             value += calRoc->GetValue(icol,irow);
+             nb++;         
+           }
+         }
+         value = value/nb;
+         if(typedet == 0) his->SetBinContent(binx,biny,value*factor);
+         else his->SetBinContent(binx,biny,value+factor);
+       }
+      }
+    }    
+  }
+  his->SetXTitle("x (cm)");
+  his->SetYTitle("y (cm)");
+  his->SetStats(0);
+  his->SetMaximum(max);
+  his->SetMinimum(min);
+  return his;
+}
+
+//_____________________________________________________________________________
+Bool_t AliTRDCalPad::Add(Float_t c1)
+{
+    //
+    // add constant for all channels of all ROCs
+    //
+
+  Bool_t result = kTRUE;
+  for (Int_t idet = 0; idet < kNdet; idet++) {
+    if (fROC[idet]){
+      if(!fROC[idet]->Add(c1)) result = kFALSE;
     }
   }
+  return result;
+}
+
+//_____________________________________________________________________________
+Bool_t AliTRDCalPad::Multiply(Float_t c1)
+{
+    //
+    // multiply constant for all channels of all ROCs
+    //
+  Bool_t result = kTRUE;
+  for (Int_t idet = 0; idet < kNdet; idet++) {
+    if (fROC[idet]){
+      if(!fROC[idet]->Multiply(c1)) result = kFALSE;
+    }
+  }
+  return result;
+}
 
+//_____________________________________________________________________________
+Bool_t AliTRDCalPad::Add(const AliTRDCalPad * pad, Double_t c1
+                       , const AliTRDCalDet * calDet1, const AliTRDCalDet *calDet2
+                       , Int_t type)
+{
+    //
+    // add calpad channel by channel multiplied by c1 - all ROCs
+    // If calDet1 and calDet2, the correct the CalROC from the detector coefficient
+    // then you have calDet1 and the calPad together
+    // calDet2 and pad together
+    // typedet == 0 for gain and vdrift
+    // typedet == 1 for t0
+    //
+  Float_t kEpsilon = 0.000000000001;
+  Double_t factor1 = 0.0;
+  Double_t factor2 = 0.0;
+  if(type == 0) {
+    factor1 = 1.0;
+    factor2 = 1.0;
+  }
+  Bool_t result = kTRUE;
+  for (Int_t idet = 0; idet < kNdet; idet++) {
+    if(calDet1) factor1 = calDet1->GetValue(idet);
+    if(calDet2) factor2 = calDet2->GetValue(idet);
+    if (fROC[idet]){
+      if(type == 0){
+       if(TMath::Abs(factor1) > kEpsilon){
+         if(!fROC[idet]->Add(pad->GetCalROC(idet),c1*factor2/factor1)) result = kFALSE;
+       }
+       else result = kFALSE;
+      }
+      else{
+       AliTRDCalROC *croc = new AliTRDCalROC((const AliTRDCalROC) *pad->GetCalROC(idet));
+       if(!croc->Add(factor2)) result = kFALSE;
+       if(!fROC[idet]->Add(croc,c1)) result = kFALSE;
+      }
+    }   
+  }
+  return result;
 }
 
+//_____________________________________________________________________________
+Bool_t AliTRDCalPad::Multiply(const AliTRDCalPad * pad, const AliTRDCalDet * calDet1
+                            , const AliTRDCalDet *calDet2, Int_t type)
+{
+    //
+    // multiply calpad channel by channel - all ROCs
+    // If calDet1 and calDet2, the correct the CalROC from the detector coefficient
+    // then you have calDet1 and the calPad together
+    // calDet2 and pad together
+    // typedet == 0 for gain and vdrift
+    // typedet == 1 for t0
+    //
+  Float_t kEpsilon = 0.000000000001;
+  Bool_t result = kTRUE;
+  Double_t factor1 = 0.0;
+  Double_t factor2 = 0.0;
+  if(type == 0){
+    factor1 = 1.0;
+    factor2 = 1.0;
+  }
+  for (Int_t idet = 0; idet < kNdet; idet++) {
+    if(calDet1) factor1 = calDet1->GetValue(idet);
+    if(calDet2) factor2 = calDet2->GetValue(idet);
+    if (fROC[idet]){
+      if(type == 0){
+       if(TMath::Abs(factor1) > kEpsilon){
+         AliTRDCalROC *croc = new AliTRDCalROC((const AliTRDCalROC) *pad->GetCalROC(idet));
+         if(!croc->Multiply(factor2)) result = kFALSE;
+         if(!fROC[idet]->Multiply(croc)) result = kFALSE;
+       }
+       else result = kFALSE;
+      }
+      else{
+       AliTRDCalROC *croc2 = new AliTRDCalROC((const AliTRDCalROC) *pad->GetCalROC(idet));
+       if(!croc2->Add(factor2)) result = kFALSE;
+       if(!fROC[idet]->Add(factor1)) result = kFALSE;
+       if(!fROC[idet]->Multiply(croc2)) result = kFALSE;
+       if(!fROC[idet]->Add(-factor1)) result = kFALSE;
+      }
+    }
+  }
+  return result;
+}
+
+//_____________________________________________________________________________
+Bool_t AliTRDCalPad::Divide(const AliTRDCalPad * pad, const AliTRDCalDet * calDet1
+                          , const AliTRDCalDet *calDet2, Int_t type)
+{
+    //
+    // divide calpad channel by channel - all ROCs
+    // If calDet1 and calDet2, the correct the CalROC from the detector coefficient
+    // then you have calDet1 and the calPad together
+    // calDet2 and pad together
+    // typedet == 0 for gain and vdrift
+    // typedet == 1 for t0
+    //
+  Float_t kEpsilon = 0.000000000001;
+  Bool_t result = kTRUE;
+  Double_t factor1 = 0.0;
+  Double_t factor2 = 0.0;
+  if(type == 0){
+    factor1 = 1.0;
+    factor2 = 1.0;
+  }
+  for (Int_t idet = 0; idet < kNdet; idet++) {
+    if(calDet1) factor1 = calDet1->GetValue(idet);
+    if(calDet2) factor2 = calDet2->GetValue(idet);
+    if (fROC[idet]){
+      if(type == 0){
+       if(TMath::Abs(factor1) > kEpsilon){
+         AliTRDCalROC *croc = new AliTRDCalROC((const AliTRDCalROC) *pad->GetCalROC(idet));
+         if(!croc->Multiply(factor2)) result = kFALSE;
+         if(!fROC[idet]->Divide(croc)) result = kFALSE;
+       }
+       else result = kFALSE;
+      }
+      else{
+       AliTRDCalROC *croc2 = new AliTRDCalROC((const AliTRDCalROC) *pad->GetCalROC(idet));
+       if(!croc2->Add(factor2)) result = kFALSE;
+       if(!fROC[idet]->Add(factor1)) result = kFALSE;
+       if(!fROC[idet]->Divide(croc2)) result = kFALSE;
+       if(!fROC[idet]->Add(-factor1)) result = kFALSE;
+      }
+    }
+  }
+  return result;
+}
index 8c1cc3b5a2305595019d001c9eb6dfa6764ab6e2..c365bb42bf8a781a82bca6fc36767bd473078109 100644 (file)
@@ -7,7 +7,7 @@
 
 ///////////////////////////////////////////////////////////////////////////////
 //                                                                           //
-//  TRD calibration class for parameters which are saved per pad                 //
+//  TRD calibration class for parameters which are saved per pad             //
 //                                                                           //
 ///////////////////////////////////////////////////////////////////////////////
 
 
 class AliTRDCalROC;
 class AliTRDCalDet;
+class TH2F;
+class TH1F;
 
-class AliTRDCalPad : public TNamed {
+class AliTRDCalPad : public TNamed 
+{
 
  public:
  
@@ -36,7 +39,26 @@ class AliTRDCalPad : public TNamed {
   AliTRDCalROC       *GetCalROC(Int_t p, Int_t c, Int_t s) const
                                                         { return fROC[GetDet(p,c,s)]; };
   
-  void                ScaleROCs(AliTRDCalDet* values);
+  Bool_t              ScaleROCs(const AliTRDCalDet* values);
+
+  // Statistic
+  Double_t GetMeanRMS(Double_t &rms, const AliTRDCalDet *calDet = 0, Int_t type = 0);
+  Double_t GetMean(const AliTRDCalDet *calDet = 0, Int_t type = 0, AliTRDCalPad* outlierPad = 0);
+  Double_t GetRMS(const AliTRDCalDet *calDet = 0, Int_t type = 0, AliTRDCalPad* outlierPad = 0) ;
+  Double_t GetMedian(const AliTRDCalDet *calDet = 0, Int_t type = 0, AliTRDCalPad* outlierPad = 0) ;
+  Double_t GetLTM(Double_t *sigma=0, Double_t fraction=0.9, const AliTRDCalDet *calDet = 0, Int_t type = 0, AliTRDCalPad* outlierPad = 0);
+
+  // Plot functions
+  TH1F    *MakeHisto1D(const AliTRDCalDet *calDet = 0, Int_t typedet=0, Float_t min=4, Float_t max=-4,Int_t type=0);
+  TH2F    *MakeHisto2DSmPl(Int_t sm, Int_t pl, const AliTRDCalDet *calDet = 0, Int_t typedet=0, Float_t min=4, Float_t max=-4,Int_t type=0);
+  TH2F    *MakeHisto2DCh(Int_t ch, const AliTRDCalDet *calDet = 0, Int_t typedet=0, Float_t min=4, Float_t max=-4,Int_t type=0);
+
+  // Algebra functions
+  Bool_t Add(Float_t c1);
+  Bool_t Multiply(Float_t c1);
+  Bool_t Add(const AliTRDCalPad * pad, Double_t c1 = 1, const AliTRDCalDet * calDet1 = 0, const AliTRDCalDet *calDet2 = 0, Int_t type = 0);
+  Bool_t Multiply(const AliTRDCalPad * pad, const AliTRDCalDet * calDet1 = 0, const AliTRDCalDet *calDet2 = 0, Int_t type = 0);
+  Bool_t Divide(const AliTRDCalPad * pad, const AliTRDCalDet * calDet1 = 0, const AliTRDCalDet *calDet2 = 0, Int_t type = 0);
 
  protected:
 
index 6530c59d3291a32ccd26bc364dc9a89c9eff78c2..9e9545fdda3d9386abde4c83121661501765e1da 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,290 @@ 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)))) {
+         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;
 }
index 6e8f76532b1116c055ea19756cd336ba3d98d0b5..c6b4aa6535a96af41cd3f3f5f43ae786aa30d20f 100644 (file)
@@ -3,7 +3,7 @@
 /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
  * See cxx source for full Copyright notice                               */
 
-/* $Id: AliTRDCalROC.h,v */
+/* $Id$ */
 
 //////////////////////////////////////////////////
 //                                              //
 //////////////////////////////////////////////////
 
 #include <TObject.h>
+#include <TMath.h>
+#include <TLinearFitter.h>
+
+class TArrayI;
+class TArrayF;
+class TH2F;
+class TH1F;
+class TGraph2D;
+class TH2D;
 
 //_____________________________________________________________________________
-class AliTRDCalROC : public TObject {
+class AliTRDCalROC : public TObject 
+{
 
  public:
 
@@ -39,8 +49,24 @@ class AliTRDCalROC : public TObject {
   void          SetValue(Int_t col, Int_t row, Float_t value) 
                                                 { SetValue(GetChannel(col,row), value);    };
 
-  void          Scale(Float_t value);
-
+  // statistic
+  //
+  Double_t GetMean(AliTRDCalROC *outlierROC=0); 
+  Double_t GetRMS(AliTRDCalROC *outlierROC=0);
+  Double_t GetMedian(AliTRDCalROC *outlierROC=0);
+  Double_t GetLTM(Double_t *sigma=0, Double_t fraction=0.9, AliTRDCalROC *outlierROC=0);
+
+  // algebra
+  Bool_t Add(Float_t c1);
+  Bool_t Multiply(Float_t c1);
+  Bool_t Add(const AliTRDCalROC * roc, Double_t c1 = 1);
+  Bool_t Multiply(const AliTRDCalROC * roc);   
+  Bool_t Divide(const AliTRDCalROC * roc);   
+  
+  //Plots
+  TH2F *   MakeHisto2D(Float_t min, Float_t max,Int_t type);
+  TH1F *   MakeHisto1D(Float_t min, Float_t max,Int_t type);
+  
  protected:
 
   Int_t     fPla;              //  Plane number