Generation of N gaussian numbers through a known covariance matrix. It probably will...
[u/mrichter/AliRoot.git] / 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