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 <TRandom.h>
  30
  31 #include "AliGausCorr.h"
  32
  33 ClassImp(AliGausCorr)
  34
  35 //_______________________________________________________________________
  36 AliGausCorr::AliGausCorr():
  37   fSize(0),
  38   fCv(0)
  39 {
  40   //
  41   // Default constructor
  42   //
  43 }
  44
  45 //_______________________________________________________________________
  46 AliGausCorr::AliGausCorr(const TMatrixD & vec, Int_t size):
  47   fSize(size),
  48   fCv(new TMatrixD(fSize,fSize))
  49 {
  50   //
  51   // Standard constructor
  52   //
  53   for(Int_t j=0;j<fSize;j++){
  54     double accum = 0;
  55     for(Int_t k=0;k<j;k++){
  56       accum += (*fCv)(j,k)* (*fCv)(j,k);
  57     }
  58     (*fCv)(j,j)=TMath::Sqrt(TMath::Abs(vec(j,j)-accum));
  59     for(Int_t i=j+1;i<fSize;i++){
  60       accum = 0;
  61       for(Int_t k=0;k<j;k++){
  62         accum+=(*fCv)(i,k)* (*fCv)(j,k);
  63       }
  64       (*fCv)(i,j) = (vec(i,j)-accum) / (*fCv)(j,j);
  65     }
  66   }
  67 }
  68
  69 //_______________________________________________________________________
  70 AliGausCorr::AliGausCorr(const AliGausCorr & tgcorr):
  71   TObject(tgcorr),
  72   fSize(tgcorr.fSize),
  73   fCv(new TMatrixD(fSize,fSize))
  74 {
  75   //
  76   // Copy contructor
  77   //
  78   for(Int_t i=0;i<fSize;i++){
  79     for(Int_t j=0;j<fSize;j++)(*fCv)(i,j)=(*tgcorr.fCv)(i,j);
  80   }
  81 }
  82
  83 //_______________________________________________________________________
  84 AliGausCorr::~AliGausCorr()
  85 {
  86   // Destructor
  87   delete fCv;
  88 }
  89
  90 //_______________________________________________________________________
  91 void AliGausCorr::GetGaussN(TArrayD &vec) const
  92 {
  93   // return fSize correlated gaussian numbers
  94
  95   TArrayD tmpv(fSize);
  96
  97   for(Int_t i=0; i<fSize; i++){
  98     Double_t x, y, z;
  99     do {
 100       y = gRandom->Rndm();
 101     } while (!y);
 102     z = gRandom->Rndm();
 103     x = z * 6.283185;
 104     tmpv[i] = TMath::Sin(x)*TMath::Sqrt(-2*TMath::Log(y));
 105   }
 106
 107   for(Int_t i=0;i<fSize;i++){
 108     vec[i]=0;
 109     for(Int_t j=0;j<=i;j++)vec[i] += (*fCv)(i,j)* tmpv[j];
 110   }
 111
 112 }
 113
 114 //_______________________________________________________________________
 115 void AliGausCorr::PrintCv() const
 116 {
 117   // Printout of the "square root" cov. matrix
 118   printf("\n AliGausCorr - triangular matrix \n");
 119   for(Int_t i=0; i<fSize;i++){
 120     for(Int_t j=0;j<(fSize-1);j++){
 121       if(j==0){
 122         printf("| %12.4f ",(*fCv)(i,j));
 123       }
 124       else {
 125         printf(" %12.4f ",(*fCv)(i,j));
 126       }
 127     }
 128     printf(" %12.4f | \n",(*fCv)(i,fSize-1));
 129   }
 130   printf("\n");
 131 }
 132
 133 //_______________________________________________________________________
 134 AliGausCorr & AliGausCorr::operator=(const AliGausCorr & tgcorr)
 135 {
 136   // Assignment operator
 137   if(&tgcorr != this && tgcorr.fSize!=fSize){
 138     if(fCv)delete fCv;
 139     fSize = tgcorr.fSize;
 140     fCv = new TMatrixD(fSize,fSize);
 141   }
 142   if(&tgcorr != this){
 143     for(Int_t i=0;i<fSize;i++){
 144       for(Int_t j=0;j<fSize;j++)(*fCv)(i,j)=(*tgcorr.fCv)(i,j);
 145     }
 146   }
 147
 148   return *this;
 149 }