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