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