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