Code clean-up (F.Carminati)
[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
25 #include <Riostream.h>
26 #include <TArrayD.h>
27 #include <TMatrixD.h>
28 #include <TRandom.h>
29 #include "AliGausCorr.h"
30
31 ClassImp(AliGausCorr)
32
33 //_______________________________________________________________________
34 AliGausCorr::AliGausCorr():
35   fSize(0),
36   fCv(0)
37 {
38   //
39   // Default constructor
40   //
41 }
42
43 //_______________________________________________________________________
44 AliGausCorr::AliGausCorr(const TMatrixD & vec, Int_t size):
45   fSize(size),
46   fCv(new TMatrixD(fSize,fSize))
47 {
48   //
49   // Standard constructor
50   //
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
67 //_______________________________________________________________________
68 AliGausCorr::AliGausCorr(const AliGausCorr & tgcorr):
69   TObject(tgcorr),
70   fSize(tgcorr.fSize),
71   fCv(new TMatrixD(fSize,fSize))
72 {
73   //
74   // Copy contructor
75   //
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
81 //_______________________________________________________________________
82 AliGausCorr::~AliGausCorr()
83 {
84   // Destructor
85   delete fCv;
86 }
87
88 //_______________________________________________________________________
89 void 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
112 //_______________________________________________________________________
113 void 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
131 //_______________________________________________________________________
132 AliGausCorr & 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 }