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