1 /**************************************************************************
2 * Copyright(c) 2001-2002, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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 ////////////////////////////////////////////////////////////////////////
27 #include <Riostream.h>
32 #include "AliGausCorr.h"
36 //_______________________________________________________________________
37 AliGausCorr::AliGausCorr():
42 // Default constructor
46 //_______________________________________________________________________
47 AliGausCorr::AliGausCorr(const TMatrixD & vec, Int_t size):
49 fCv(new TMatrixD(fSize,fSize))
52 // Standard constructor
54 for(Int_t j=0;j<fSize;j++){
56 for(Int_t k=0;k<j;k++){
57 accum += (*fCv)(j,k)* (*fCv)(j,k);
59 (*fCv)(j,j)=TMath::Sqrt(TMath::Abs(vec(j,j)-accum));
60 for(Int_t i=j+1;i<fSize;i++){
62 for(Int_t k=0;k<j;k++){
63 accum+=(*fCv)(i,k)* (*fCv)(j,k);
65 (*fCv)(i,j) = (vec(i,j)-accum) / (*fCv)(j,j);
70 //_______________________________________________________________________
71 AliGausCorr::AliGausCorr(const AliGausCorr & tgcorr):
74 fCv(new TMatrixD(fSize,fSize))
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);
84 //_______________________________________________________________________
85 AliGausCorr::~AliGausCorr()
91 //_______________________________________________________________________
92 void AliGausCorr::GetGaussN(TArrayD &vec) const
94 // return fSize correlated gaussian numbers
98 for(Int_t i=0; i<fSize; i++){
105 tmpv[i] = TMath::Sin(x)*TMath::Sqrt(-2*TMath::Log(y));
108 for(Int_t i=0;i<fSize;i++){
110 for(Int_t j=0;j<=i;j++)vec[i] += (*fCv)(i,j)* tmpv[j];
115 //_______________________________________________________________________
116 void AliGausCorr::PrintCv() const
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++){
123 printf("| %12.4f ",(*fCv)(i,j));
126 printf(" %12.4f ",(*fCv)(i,j));
129 printf(" %12.4f | \n",(*fCv)(i,fSize-1));
134 //_______________________________________________________________________
135 AliGausCorr & AliGausCorr::operator=(const AliGausCorr & tgcorr)
137 // Assignment operator
138 if(&tgcorr != this && tgcorr.fSize!=fSize){
140 fSize = tgcorr.fSize;
141 fCv = new TMatrixD(fSize,fSize);
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);