TPC/AliTPCCalROC.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16
  17 ///////////////////////////////////////////////////////////////////////////////
  18 //                                                                           //
  19 //  Calibration base class for a single ROC                                  //
  20 //  Contains one float value per pad                                         //
  21 //     mapping of the pads taken form AliTPCROC                              //
  22 //                                                                           //
  23 ///////////////////////////////////////////////////////////////////////////////
  24
  25 #include "AliTPCCalROC.h"
  26 #include "TMath.h"
  27 #include "TClass.h"
  28 #include "TFile.h"
  29 #include "TH1F.h"
  30 #include "TH2F.h"
  31 #include "AliMathBase.h"
  32 ClassImp(AliTPCCalROC)
  33
  34
  35 //_____________________________________________________________________________
  36 AliTPCCalROC::AliTPCCalROC()
  37              :TObject(),
  38               fSector(0),
  39               fNChannels(0),
  40               fNRows(0),
  41               fIndexes(0),
  42               fData(0)
  43 {
  44   //
  45   // Default constructor
  46   //
  47
  48 }
  49
  50 //_____________________________________________________________________________
  51 AliTPCCalROC::AliTPCCalROC(UInt_t  sector)
  52              :TObject(),
  53               fSector(0),
  54               fNChannels(0),
  55               fNRows(0),
  56               fIndexes(0),
  57               fData(0)
  58 {
  59   //
  60   // Constructor that initializes a given sector
  61   //
  62   fSector = sector;
  63   fNChannels    =  AliTPCROC::Instance()->GetNChannels(fSector);
  64   fNRows        =  AliTPCROC::Instance()->GetNRows(fSector);
  65   fIndexes      =  AliTPCROC::Instance()->GetRowIndexes(fSector);
  66   fData = new Float_t[fNChannels];
  67   for (UInt_t  idata = 0; idata< fNChannels; idata++) fData[idata] = 0.;
  68 }
  69
  70 //_____________________________________________________________________________
  71 AliTPCCalROC::AliTPCCalROC(const AliTPCCalROC &c)
  72              :TObject(c),
  73               fSector(0),
  74               fNChannels(0),
  75               fNRows(0),
  76               fIndexes(0),
  77               fData(0)
  78 {
  79   //
  80   // AliTPCCalROC copy constructor
  81   //
  82   fSector = c.fSector;
  83   fNChannels    =  AliTPCROC::Instance()->GetNChannels(fSector);
  84   fNRows        =  AliTPCROC::Instance()->GetNRows(fSector);
  85   fIndexes      =  AliTPCROC::Instance()->GetRowIndexes(fSector);
  86   //
  87   fData   = new Float_t[fNChannels];
  88   for (UInt_t  idata = 0; idata< fNChannels; idata++) fData[idata] = c.fData[idata];
  89 }
  90 //____________________________________________________________________________
  91 AliTPCCalROC & AliTPCCalROC::operator =(const AliTPCCalROC & param)
  92 {
  93   //
  94   // assignment operator - dummy
  95   //
  96   fData=param.fData;
  97   return (*this);
  98 }
  99
 100
 101 //_____________________________________________________________________________
 102 AliTPCCalROC::~AliTPCCalROC()
 103 {
 104   //
 105   // AliTPCCalROC destructor
 106   //
 107   if (fData) {
 108     delete [] fData;
 109     fData = 0;
 110   }
 111 }
 112
 113
 114
 115 void AliTPCCalROC::Streamer(TBuffer &R__b)
 116 {
 117    // Stream an object of class AliTPCCalROC.
 118    if (R__b.IsReading()) {
 119       AliTPCCalROC::Class()->ReadBuffer(R__b, this);
 120       fIndexes =  AliTPCROC::Instance()->GetRowIndexes(fSector);
 121    } else {
 122       AliTPCCalROC::Class()->WriteBuffer(R__b,this);
 123    }
 124 }
 125
 126 //  //
 127 //   // algebra
 128 //   void Add(Float_t c1);
 129 //   void Multiply(Float_t c1);
 130 //   void Add(const AliTPCCalROC * roc, Double_t c1 = 1);
 131 //   void Divide(const AliTPCCalROC * roc);
 132
 133 void AliTPCCalROC::Add(Float_t c1){
 134   //
 135   // add constant
 136   //
 137   for (UInt_t  idata = 0; idata< fNChannels; idata++) fData[idata]+=c1;
 138 }
 139 void AliTPCCalROC::Multiply(Float_t c1){
 140   //
 141   // add constant
 142   //
 143   for (UInt_t  idata = 0; idata< fNChannels; idata++) fData[idata]*=c1;
 144 }
 145
 146 void AliTPCCalROC::Add(const AliTPCCalROC * roc, Double_t c1){
 147   //
 148   // add values
 149   //
 150   for (UInt_t  idata = 0; idata< fNChannels; idata++){
 151     fData[idata]+=roc->fData[idata]*c1;
 152   }
 153 }
 154
 155
 156 void AliTPCCalROC::Multiply(const AliTPCCalROC*  roc) {
 157   //
 158   // multiply values - per by pad
 159   //
 160   for (UInt_t  idata = 0; idata< fNChannels; idata++){
 161     fData[idata]*=roc->fData[idata];
 162   }
 163 }
 164
 165
 166 void AliTPCCalROC::Divide(const AliTPCCalROC*  roc) {
 167   //
 168   // divide values
 169   //
 170   Float_t kEpsilon=0.00000000000000001;
 171   for (UInt_t  idata = 0; idata< fNChannels; idata++){
 172     if (TMath::Abs(roc->fData[idata])>kEpsilon)
 173       fData[idata]/=roc->fData[idata];
 174   }
 175 }
 176
 177
 178
 179
 180 Double_t AliTPCCalROC::GetLTM(Double_t *sigma, Double_t fraction){
 181   //
 182   //  Calculate LTM mean and sigma
 183   //
 184   Double_t *ddata = new Double_t[fNChannels];
 185   Double_t mean=0, lsigma=0;
 186   Int_t hh = TMath::Min(TMath::Nint(fraction *fNChannels), Int_t(fNChannels));
 187   for (UInt_t i=0;i<fNChannels;i++) ddata[i]= fData[i];
 188   AliMathBase::EvaluateUni(UInt_t(fNChannels),ddata, mean, lsigma, hh);
 189   if (sigma) *sigma=lsigma;
 190   delete [] ddata;
 191   return mean;
 192 }
 193
 194 TH1F * AliTPCCalROC::MakeHisto1D(Float_t min, Float_t max,Int_t type){
 195   //
 196   // make 1D histo
 197   // type -1 = user defined range
 198   //       0 = nsigma cut nsigma=min
 199   //       1 = delta cut around median delta=min
 200   if (type>=0){
 201     if (type==0){
 202       // nsigma range
 203       Float_t mean  = GetMean();
 204       Float_t sigma = GetRMS();
 205       Float_t nsigma = TMath::Abs(min);
 206       min = mean-nsigma*sigma;
 207       max = mean+nsigma*sigma;
 208     }
 209     if (type==1){
 210       // fixed range
 211       Float_t mean   = GetMedian();
 212       Float_t  delta = min;
 213       min = mean-delta;
 214       max = mean+delta;
 215     }
 216     if (type==2){
 217       //
 218       // LTM mean +- nsigma
 219       //
 220       Double_t sigma;
 221       Float_t mean  = GetLTM(&sigma,max);
 222       sigma*=min;
 223       min = mean-sigma;
 224       max = mean+sigma;
 225     }
 226   }
 227   char  name[1000];
 228   sprintf(name,"%s ROC 1D%d",GetTitle(),fSector);
 229   TH1F * his = new TH1F(name,name,100, min,max);
 230   for (UInt_t irow=0; irow<fNRows; irow++){
 231     UInt_t npads = (Int_t)GetNPads(irow);
 232     for (UInt_t ipad=0; ipad<=npads; ipad++){
 233       his->Fill(GetValue(irow,ipad));
 234     }
 235   }
 236   return his;
 237 }
 238
 239
 240
 241 TH2F * AliTPCCalROC::MakeHisto2D(Float_t min, Float_t max,Int_t type){
 242   //
 243   // make 2D histo
 244   // type -1 = user defined range
 245   //       0 = nsigma cut nsigma=min
 246   //       1 = delta cut around median delta=min
 247   if (type>=0){
 248     if (type==0){
 249       // nsigma range
 250       Float_t mean  = GetMean();
 251       Float_t sigma = GetRMS();
 252       Float_t nsigma = TMath::Abs(min);
 253       min = mean-nsigma*sigma;
 254       max = mean+nsigma*sigma;
 255     }
 256     if (type==1){
 257       // fixed range
 258       Float_t mean   = GetMedian();
 259       Float_t  delta = min;
 260       min = mean-delta;
 261       max = mean+delta;
 262     }
 263     if (type==2){
 264       Double_t sigma;
 265       Float_t mean  = GetLTM(&sigma,max);
 266       sigma*=min;
 267       min = mean-sigma;
 268       max = mean+sigma;
 269
 270     }
 271   }
 272   UInt_t maxPad = 0;
 273   for (UInt_t irow=0; irow<fNRows; irow++){
 274     if (GetNPads(irow)>maxPad) maxPad = GetNPads(irow);
 275   }
 276   char  name[1000];
 277   sprintf(name,"%s ROC%d",GetTitle(),fSector);
 278   TH2F * his = new TH2F(name,name,fNRows+10,-5, fNRows+5, maxPad+10, -(Int_t(maxPad/2))-5, maxPad/2+5);
 279   for (UInt_t irow=0; irow<fNRows; irow++){
 280     UInt_t npads = (Int_t)GetNPads(irow);
 281     for (UInt_t ipad=0; ipad<=npads; ipad++){
 282       his->Fill(irow+0.5,Int_t(ipad)-Int_t(npads/2)+0.5,GetValue(irow,ipad));
 283     }
 284   }
 285   his->SetMaximum(max);
 286   his->SetMinimum(min);
 287   return his;
 288 }
 289
 290 TH2F * AliTPCCalROC::MakeHistoOutliers(Float_t delta, Float_t fraction, Int_t type){
 291   //
 292   // Make Histogram with outliers
 293   // mode = 0 - sigma cut used
 294   // mode = 1 - absolute cut used
 295   // fraction - fraction of values used to define sigma
 296   // delta - in mode 0 - nsigma cut
 297   //            mode 1 - delta cut
 298   Double_t sigma;
 299   Float_t mean  = GetLTM(&sigma,fraction);
 300   if (type==0) delta*=sigma;
 301   UInt_t maxPad = 0;
 302   for (UInt_t irow=0; irow<fNRows; irow++){
 303     if (GetNPads(irow)>maxPad) maxPad = GetNPads(irow);
 304   }
 305
 306   char  name[1000];
 307   sprintf(name,"%s ROC Outliers%d",GetTitle(),fSector);
 308   TH2F * his = new TH2F(name,name,fNRows+10,-5, fNRows+5, maxPad+10, -(Int_t(maxPad/2))-5, maxPad/2+5);
 309   for (UInt_t irow=0; irow<fNRows; irow++){
 310     UInt_t npads = (Int_t)GetNPads(irow);
 311     for (UInt_t ipad=0; ipad<=npads; ipad++){
 312       if (TMath::Abs(GetValue(irow,ipad)-mean)>delta)
 313         his->Fill(irow+0.5,Int_t(ipad)-Int_t(npads/2)+0.5,1);
 314     }
 315   }
 316   return his;
 317 }
 318
 319
 320
 321 void AliTPCCalROC::Draw(Option_t* opt){
 322   //
 323   // create histogram with values and draw it
 324   //
 325   TH1 * his=0;
 326   TString option=opt;
 327   option.ToUpper();
 328   if (option.Contains("1D")){
 329     his = MakeHisto1D();
 330   }
 331   else{
 332     his = MakeHisto2D();
 333   }
 334   his->Draw(option);
 335 }
 336
 337
 338
 339
 340
 341 void AliTPCCalROC::Test(){
 342   //
 343   // example function to show functionality and tes AliTPCCalROC
 344   //
 345   AliTPCCalROC  roc0(0);
 346   for (UInt_t irow = 0; irow <roc0.GetNrows(); irow++){
 347     for (UInt_t ipad = 0; ipad <roc0.GetNPads(irow); ipad++){
 348       Float_t value  = irow+ipad/1000.;
 349       roc0.SetValue(irow,ipad,value);
 350     }
 351   }
 352   //
 353   AliTPCCalROC roc1(roc0);
 354   for (UInt_t irow = 0; irow <roc1.GetNrows(); irow++){
 355     for (UInt_t ipad = 0; ipad <roc1.GetNPads(irow); ipad++){
 356       Float_t value  = irow+ipad/1000.;
 357       if (roc1.GetValue(irow,ipad)!=value){
 358         printf("Read/Write error\trow=%d\tpad=%d\n",irow,ipad);
 359       }
 360     }
 361   }
 362   TFile f("calcTest.root","recreate");
 363   roc0.Write("Roc0");
 364   AliTPCCalROC * roc2 = (AliTPCCalROC*)f.Get("Roc0");
 365   f.Close();
 366   //
 367   for (UInt_t irow = 0; irow <roc0.GetNrows(); irow++){
 368     if (roc0.GetNPads(irow)!=roc2->GetNPads(irow))
 369       printf("NPads - Read/Write error\trow=%d\n",irow);
 370     for (UInt_t ipad = 0; ipad <roc1.GetNPads(irow); ipad++){
 371       Float_t value  = irow+ipad/1000.;
 372       if (roc2->GetValue(irow,ipad)!=value){
 373         printf("Read/Write error\trow=%d\tpad=%d\n",irow,ipad);
 374       }
 375     }
 376   }
 377 }
 378