TPC/AliH2F.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 $Log$
  18 Revision 1.2.4.2  2000/06/26 07:39:42  kowal2
  19 Changes to obey the coding rules
  20
  21 Revision 1.2.4.1  2000/06/25 08:38:41  kowal2
  22 Splitted from AliTPCtracking
  23
  24 Revision 1.2  2000/04/17 09:37:33  kowal2
  25 removed obsolete AliTPCDigitsDisplay.C
  26
  27 Revision 1.1.4.2  2000/04/10 11:32:37  kowal2
  28
  29 "ROOT"-based class with some extra functionality
  30
  31 */
  32 //
  33
  34 //----------------------------------------------------------------------------
  35 //  Author:   Marian Ivanov
  36 //
  37 //  Implementation of class AliH2F
  38 //
  39 //-----------------------------------------------------------------------------
  40
  41 #include "AliH2F.h"
  42 #include "TClonesArray.h"
  43 #include "AliTPC.h"
  44 #include "TRandom.h"
  45 #include "AliCluster.h"
  46 #include "AliTPCClusterFinder.h"
  47 //*KEEP,TMath.
  48
  49 // other include files follow here
  50
  51
  52 ClassImp(AliH2F)
  53 //***********************************************************************
  54 //***********************************************************************
  55 //***********************************************************************
  56 //***********************************************************************
  57 AliH2F::AliH2F():TH2F()
  58 {
  59   //
  60
  61 }
  62 AliH2F::AliH2F(const Text_t *name,const Text_t *title,
  63                        Int_t nbinsx,Axis_t xlow,Axis_t xup
  64                        ,Int_t nbinsy,Axis_t ylow,Axis_t yup):
  65   TH2F(name,title,nbinsx,xlow,xup
  66        ,nbinsy,ylow,yup)
  67 {
  68   //
  69
  70 }
  71
  72 AliH2F::~AliH2F()
  73 {
  74   //
  75 }
  76
  77 AliH2F::AliH2F(const AliH2F &his)
  78 {
  79   //
  80
  81 }
  82
  83 AliH2F & AliH2F::operator = (const AliH2F & his)
  84 {
  85   //
  86   return *this;
  87 }
  88
  89 TClonesArray * AliH2F::FindPeaks(Float_t threshold, Float_t noise)
  90 {
  91   //find peaks and write it in form of AliTPCcluster to array
  92
  93   //firstly we need to create object for cluster finding
  94   //and fill it with contents of histogram
  95   AliTPCClusterFinder cfinder;
  96   cfinder.SetThreshold(threshold);
  97   cfinder.SetNoise(noise);
  98   cfinder.GetHisto(this);
  99   return cfinder.FindPeaks3();
 100 }
 101
 102 void AliH2F::ClearSpectrum()
 103 {
 104   //clera histogram
 105   Int_t dimx =  fXaxis.GetNbins();
 106   Int_t dimy =  fYaxis.GetNbins();
 107   for (Int_t i = 0 ;i<dimx;i++)
 108     for (Int_t j = 0 ;j<dimy;j++)
 109       {
 110         SetCellContent(i,j,0);
 111         SetCellError(i,j,0);
 112       }
 113 }
 114
 115
 116 void AliH2F::AddNoise(Float_t sn)
 117 {
 118   // add gauss noise with sigma sn
 119   Int_t dimx =  fXaxis.GetNbins();
 120   Int_t dimy =  fYaxis.GetNbins();
 121   for (Int_t i = 0 ;i<dimx;i++)
 122     for (Int_t j = 0 ;j<dimy;j++)
 123       {
 124         Float_t noise = gRandom->Gaus(0,sn);
 125         Float_t oldv  =GetCellContent(i,j);
 126         Float_t olds  =GetCellError(i,j);
 127         if (noise >0)
 128           {
 129             SetCellContent(i,j,noise+oldv);
 130             SetCellError(i,j,TMath::Sqrt((noise*noise+olds*olds)));
 131           }
 132       }
 133 }
 134 void AliH2F::AddGauss(Float_t x, Float_t y,
 135                           Float_t sx, Float_t sy, Float_t max)
 136 {
 137   //transform to histogram coordinata
 138   Int_t dimx =  fXaxis.GetNbins();
 139   Int_t dimy =  fYaxis.GetNbins();
 140   Float_t dx =(GetXaxis()->GetXmax()-GetXaxis()->GetXmin())/Float_t(dimx);
 141   Float_t dy =(GetYaxis()->GetXmax()-GetYaxis()->GetXmin())/Float_t(dimy);
 142   //  x=(x-GetXaxis()->GetXmin())/dx;
 143   //y=(y-GetYaxis()->GetXmin())/dy;
 144   sx/=dx;
 145   sy/=dy;
 146
 147
 148   for (Int_t i = 0 ;i<dimx;i++)
 149     for (Int_t j = 0 ;j<dimy;j++)
 150       {
 151         Float_t x2 =GetXaxis()->GetBinCenter(i+1);
 152         Float_t y2 =GetYaxis()->GetBinCenter(j+1);
 153         Float_t dx2 = (x2-x)*(x2-x);
 154         Float_t dy2 = (y2-y)*(y2-y);
 155         Float_t amp =max*exp(-(dx2/(2*sx*sx)+dy2/(2*sy*sy)));
 156         //Float_t oldv  =GetCellContent(i+1,j+1);
 157         //      SetCellContent(i+1,j+1,amp+oldv);
 158         Fill(x2,y2,amp);
 159       }
 160 }
 161
 162 void AliH2F::ClearUnderTh(Int_t threshold)
 163 {
 164   //clear histogram for bin under threshold
 165   Int_t dimx =  fXaxis.GetNbins();
 166   Int_t dimy =  fYaxis.GetNbins();
 167   for (Int_t i = 0 ;i<=dimx;i++)
 168     for (Int_t j = 0 ;j<=dimy;j++)
 169       {
 170         Float_t oldv  =GetCellContent(i,j);
 171         if (oldv <threshold)
 172           SetCellContent(i,j,0);
 173       }
 174 }
 175
 176 void AliH2F::Round()
 177 {
 178   //round float to integer
 179   Int_t dimx =  fXaxis.GetNbins();
 180   Int_t dimy =  fYaxis.GetNbins();
 181   for (Int_t i = 0 ;i<=dimx;i++)
 182     for (Int_t j = 0 ;j<=dimy;j++)
 183       {
 184         Float_t oldv  =GetCellContent(i,j);
 185         oldv=(Int_t)oldv;
 186         SetCellContent(i,j,oldv);
 187       }
 188 }
 189
 190
 191
 192 AliH2F *AliH2F::GetSubrange2d(Float_t xmin, Float_t xmax,
 193                                       Float_t ymin, Float_t ymax)
 194 {
 195   //this function return pointer to the new created
 196   //histogram which is subhistogram of the
 197   //calculate number
 198   //subhistogram range must be inside histogram
 199
 200   if (xmax<=xmin) {
 201     xmin=fXaxis.GetXmin();
 202     xmax=fXaxis.GetXmax();
 203   }
 204   if (ymax<=ymin) {
 205      ymin=fYaxis.GetXmin();
 206      ymax=fYaxis.GetXmax();
 207   }
 208
 209   Int_t nx = Int_t((xmax-xmin)/(fXaxis.GetXmax()-fXaxis.GetXmin())  *
 210                    Float_t(fXaxis.GetNbins()));
 211   Int_t ny = Int_t((ymax-ymin)/(fYaxis.GetXmax()-fYaxis.GetXmin())  *
 212                    Float_t(fYaxis.GetNbins()));
 213   TString  t1 = fName ;
 214   TString  t2 = fTitle ;
 215   t1+="_subrange";
 216   t2+="_subrange";
 217   const Text_t *ktt1 = t1;
 218   const Text_t *ktt2 = t2;
 219
 220   AliH2F * sub = new AliH2F(ktt1,ktt2,nx,xmin,xmax,ny,ymin,ymax);
 221
 222   Int_t i1 = Int_t( Float_t(fXaxis.GetNbins())*(xmin-fXaxis.GetXmin())/
 223                     (fXaxis.GetXmax()-fXaxis.GetXmin()) ) ;
 224   Int_t i2 = Int_t( Float_t(fYaxis.GetNbins())*(ymin-fYaxis.GetXmin())/
 225                     (fYaxis.GetXmax()-fYaxis.GetXmin()) ) ;
 226   for (Int_t i=0;i<nx;i++)
 227     for (Int_t j=0;j<ny;j++)
 228       {
 229         Int_t index1 = GetBin(i1+i,i2+j);
 230         //        Int_t index2 = sub->GetBin(i,j);
 231         Float_t val = GetBinContent(index1);
 232         //        sub->SetBinContent(index2,val);
 233         //        Float_t err = GetBinError(index1);
 234         //sub->SetBinError(index2,GetBinError(index1));
 235         sub->SetCellContent(i,j,val);
 236       }
 237    return sub;
 238 }
 239
 240 TH1F *AliH2F::GetAmplitudes(Float_t zmin, Float_t zmax, Float_t th, Float_t xmin, Float_t xmax,
 241                                       Float_t ymin, Float_t ymax)
 242 {
 243   //this function return pointer to the new created
 244   //histogram which is subhistogram of the
 245   //calculate number
 246   //subhistogram range must be inside histogram
 247
 248   if (xmax<=xmin) {
 249     xmin=fXaxis.GetXmin();
 250     xmax=fXaxis.GetXmax();
 251   }
 252   if (ymax<=ymin) {
 253      ymin=fYaxis.GetXmin();
 254      ymax=fYaxis.GetXmax();
 255   }
 256   Int_t nx = Int_t((xmax-xmin)/(fXaxis.GetXmax()-fXaxis.GetXmin())  *
 257                    Float_t(fXaxis.GetNbins()));
 258   Int_t ny = Int_t((ymax-ymin)/(fYaxis.GetXmax()-fYaxis.GetXmin())  *
 259                    Float_t(fYaxis.GetNbins()));
 260   TString  t1 = fName ;
 261   TString  t2 = fTitle ;
 262   t1+="_amplitudes";
 263   t2+="_amplitudes";
 264   const  Text_t *ktt1 = t1;
 265   const Text_t *ktt2 = t2;
 266
 267   TH1F * h = new TH1F(ktt1,ktt2,100,zmin,zmax);
 268
 269   Int_t i1 = Int_t( Float_t(fXaxis.GetNbins())*(xmin-fXaxis.GetXmin())/
 270                     (fXaxis.GetXmax()-fXaxis.GetXmin()) ) ;
 271   Int_t i2 = Int_t( Float_t(fYaxis.GetNbins())*(ymin-fYaxis.GetXmin())/
 272                     (fYaxis.GetXmax()-fYaxis.GetXmin()) ) ;
 273   for (Int_t i=0;i<nx;i++)
 274     for (Int_t j=0;j<ny;j++)
 275       {
 276         Int_t index1 = GetBin(i1+i,i2+j);
 277         Float_t val = GetBinContent(index1);
 278         if (val>th) h->Fill(val);
 279       }
 280    return h;
 281 }
 282
 283 Float_t   AliH2F::GetOccupancy(Float_t th , Float_t xmin, Float_t xmax,
 284                              Float_t ymin, Float_t ymax)
 285 {
 286   //this function return pointer to the new created
 287   //histogram which is subhistogram of the
 288   //calculate number
 289   //subhistogram range must be inside histogram
 290
 291   if (xmax<=xmin) {
 292     xmin=fXaxis.GetXmin();
 293     xmax=fXaxis.GetXmax();
 294   }
 295   if (ymax<=ymin) {
 296      ymin=fYaxis.GetXmin();
 297      ymax=fYaxis.GetXmax();
 298   }
 299   Int_t nx = Int_t((xmax-xmin)/(fXaxis.GetXmax()-fXaxis.GetXmin())  *
 300                    Float_t(fXaxis.GetNbins()));
 301   Int_t ny = Int_t((ymax-ymin)/(fYaxis.GetXmax()-fYaxis.GetXmin())  *
 302                    Float_t(fYaxis.GetNbins()));
 303
 304   Int_t over =0;
 305   Int_t i1 = Int_t( Float_t(fXaxis.GetNbins())*(xmin-fXaxis.GetXmin())/
 306                     (fXaxis.GetXmax()-fXaxis.GetXmin()) ) ;
 307   Int_t i2 = Int_t( Float_t(fYaxis.GetNbins())*(ymin-fYaxis.GetXmin())/
 308                     (fYaxis.GetXmax()-fYaxis.GetXmin()) ) ;
 309   for (Int_t i=0;i<nx;i++)
 310     for (Int_t j=0;j<ny;j++)
 311       {
 312         Int_t index1 = GetBin(i1+i,i2+j);
 313         Float_t val = GetBinContent(index1);
 314         if (val>th) over++;
 315       }
 316   Int_t  all = nx*ny;
 317   if (all>0)  return Float_t(over)/Float_t(all);
 318   else
 319     return 0;
 320 }