STEER/AliGausCorr.cxx

   1 /**************************************************************************
   2  * Copyright(c) 2001-2002, 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 ////////////////////////////////////////////////////////////////////////
  19 // Class used to generate correlated gaussian numbers with mean
  20 // zero and known covariance matrix.
  21 // Adapted from the Fortran code in Cernlib V122 (corset, corgen)
  22 // F. James, Monte Carlo theory and practice,
  23 // Rep. Prog. Phys. 43 (1980) 1145-1189.
  24 // M.Masera 15.03.2001 9:30 - modified on 26.02.2002 17:40
  25 ////////////////////////////////////////////////////////////////////////
  26
  27 #include <Riostream.h>
  28 #include <TArrayD.h>
  29 #include <TMath.h>
  30 #include <TRandom.h>
  31
  32 #include "AliGausCorr.h"
  33
  34 ClassImp(AliGausCorr)
  35
  36 //_______________________________________________________________________
  37 AliGausCorr::AliGausCorr():
  38   fSize(0),
  39   fCv(0)
  40 {
  41   //
  42   // Default constructor
  43   //
  44 }
  45
  46 //_______________________________________________________________________
  47 AliGausCorr::AliGausCorr(const TMatrixD & vec, Int_t size):
  48   fSize(size),
  49   fCv(new TMatrixD(fSize,fSize))
  50 {
  51   //
  52   // Standard constructor
  53   //
  54   for(Int_t j=0;j<fSize;j++){
  55     double accum = 0;
  56     for(Int_t k=0;k<j;k++){
  57       accum += (*fCv)(j,k)* (*fCv)(j,k);
  58     }
  59     (*fCv)(j,j)=TMath::Sqrt(TMath::Abs(vec(j,j)-accum));
  60     for(Int_t i=j+1;i<fSize;i++){
  61       accum = 0;
  62       for(Int_t k=0;k<j;k++){
  63         accum+=(*fCv)(i,k)* (*fCv)(j,k);
  64       }
  65       (*fCv)(i,j) = (vec(i,j)-accum) / (*fCv)(j,j);
  66     }
  67   }
  68 }
  69
  70 //_______________________________________________________________________
  71 AliGausCorr::AliGausCorr(const AliGausCorr & tgcorr):
  72   TObject(tgcorr),
  73   fSize(tgcorr.fSize),
  74   fCv(new TMatrixD(fSize,fSize))
  75 {
  76   //
  77   // Copy contructor
  78   //
  79   for(Int_t i=0;i<fSize;i++){
  80     for(Int_t j=0;j<fSize;j++)(*fCv)(i,j)=(*tgcorr.fCv)(i,j);
  81   }
  82 }
  83
  84 //_______________________________________________________________________
  85 AliGausCorr::~AliGausCorr()
  86 {
  87   // Destructor
  88   delete fCv;
  89 }
  90
  91 //_______________________________________________________________________
  92 void AliGausCorr::GetGaussN(TArrayD &vec) const
  93 {
  94   // return fSize correlated gaussian numbers
  95
  96   TArrayD tmpv(fSize);
  97
  98   for(Int_t i=0; i<fSize; i++){
  99     Double_t x, y, z;
 100     do {
 101       y = gRandom->Rndm();
 102     } while (!y);
 103     z = gRandom->Rndm();
 104     x = z * 6.283185;
 105     tmpv[i] = TMath::Sin(x)*TMath::Sqrt(-2*TMath::Log(y));
 106   }
 107
 108   for(Int_t i=0;i<fSize;i++){
 109     vec[i]=0;
 110     for(Int_t j=0;j<=i;j++)vec[i] += (*fCv)(i,j)* tmpv[j];
 111   }
 112
 113 }
 114
 115 //_______________________________________________________________________
 116 void AliGausCorr::PrintCv() const
 117 {
 118   // Printout of the "square root" cov. matrix
 119   printf("\n AliGausCorr - triangular matrix \n");
 120   for(Int_t i=0; i<fSize;i++){
 121     for(Int_t j=0;j<(fSize-1);j++){
 122       if(j==0){
 123         printf("| %12.4f ",(*fCv)(i,j));
 124       }
 125       else {
 126         printf(" %12.4f ",(*fCv)(i,j));
 127       }
 128     }
 129     printf(" %12.4f | \n",(*fCv)(i,fSize-1));
 130   }
 131   printf("\n");
 132 }
 133
 134 //_______________________________________________________________________
 135 AliGausCorr & AliGausCorr::operator=(const AliGausCorr & tgcorr)
 136 {
 137   // Assignment operator
 138   if(&tgcorr != this && tgcorr.fSize!=fSize){
 139     if(fCv)delete fCv;
 140     fSize = tgcorr.fSize;
 141     fCv = new TMatrixD(fSize,fSize);
 142   }
 143   if(&tgcorr != this){
 144     for(Int_t i=0;i<fSize;i++){
 145       for(Int_t j=0;j<fSize;j++)(*fCv)(i,j)=(*tgcorr.fCv)(i,j);
 146     }
 147   }
 148
 149   return *this;
 150 }