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