Generation of N gaussian numbers through a known covariance matrix. It probably will...
[u/mrichter/AliRoot.git] / STEER / AliGausCorr.cxx
CommitLineData
740ebff3 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
30ClassImp(AliGausCorr)
31
32//________________________________________________________
33AliGausCorr::AliGausCorr()
34{
35 // Default constructor
36 fSize = 0;
37 fCv = 0;
38}
39
40//________________________________________________________
41AliGausCorr::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//________________________________________________________
63AliGausCorr::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//________________________________________________________
75AliGausCorr::~AliGausCorr()
76{
77 // Destructor
78 delete fCv;
79}
80
81//________________________________________________________
82void 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//________________________________________________________
106void 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//________________________________________________________
125AliGausCorr & 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