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