Corrected protection.
[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
acd84897 16/* $Id$ */
fb17acd4 17
740ebff3 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////////////////////////////////////////////////////////////////////////
e2afb3b6 26
b16a1b1e 27#include <Riostream.h>
740ebff3 28#include <TArrayD.h>
3010c308 29#include <TMath.h>
740ebff3 30#include <TRandom.h>
116cbefd 31
740ebff3 32#include "AliGausCorr.h"
33
34ClassImp(AliGausCorr)
35
e2afb3b6 36//_______________________________________________________________________
37AliGausCorr::AliGausCorr():
38 fSize(0),
39 fCv(0)
740ebff3 40{
e2afb3b6 41 //
740ebff3 42 // Default constructor
e2afb3b6 43 //
740ebff3 44}
45
e2afb3b6 46//_______________________________________________________________________
47AliGausCorr::AliGausCorr(const TMatrixD & vec, Int_t size):
48 fSize(size),
49 fCv(new TMatrixD(fSize,fSize))
740ebff3 50{
e2afb3b6 51 //
740ebff3 52 // Standard constructor
e2afb3b6 53 //
740ebff3 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
e2afb3b6 70//_______________________________________________________________________
71AliGausCorr::AliGausCorr(const AliGausCorr & tgcorr):
72 TObject(tgcorr),
73 fSize(tgcorr.fSize),
74 fCv(new TMatrixD(fSize,fSize))
740ebff3 75{
e2afb3b6 76 //
740ebff3 77 // Copy contructor
e2afb3b6 78 //
740ebff3 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
e2afb3b6 84//_______________________________________________________________________
740ebff3 85AliGausCorr::~AliGausCorr()
86{
87 // Destructor
88 delete fCv;
89}
90
e2afb3b6 91//_______________________________________________________________________
740ebff3 92void 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
e2afb3b6 115//_______________________________________________________________________
740ebff3 116void 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
e2afb3b6 134//_______________________________________________________________________
740ebff3 135AliGausCorr & AliGausCorr::operator=(const AliGausCorr & tgcorr)
136{
5d8718b8 137 // Assignment operator
740ebff3 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}