| 1 | #ifndef MARS_MUnfold
|
|---|
| 2 | #define MARS_MUnfold
|
|---|
| 3 |
|
|---|
| 4 | #ifndef ROOT_TH1
|
|---|
| 5 | #include <TH1.h>
|
|---|
| 6 | #endif
|
|---|
| 7 |
|
|---|
| 8 | #ifndef ROOT_TMatrixD
|
|---|
| 9 | #include <TMatrixD.h>
|
|---|
| 10 | #endif
|
|---|
| 11 |
|
|---|
| 12 | #ifndef MARS_MTask
|
|---|
| 13 | #include "MTask.h"
|
|---|
| 14 | #endif
|
|---|
| 15 |
|
|---|
| 16 | class TH2;
|
|---|
| 17 | class TH3;
|
|---|
| 18 | class TString;
|
|---|
| 19 | class TVectorD;
|
|---|
| 20 | class TH2D;
|
|---|
| 21 | class TH3D;
|
|---|
| 22 |
|
|---|
| 23 | class MUnfold : public MTask
|
|---|
| 24 | {
|
|---|
| 25 |
|
|---|
| 26 | private:
|
|---|
| 27 |
|
|---|
| 28 | TH1 *DrawMatrixClone(const TMatrixD &m, Option_t *opt);
|
|---|
| 29 | TH1 *DrawMatrixColClone(const TMatrixD &m, Option_t *opt, Int_t col);
|
|---|
| 30 | void PrintTH3Content(const TH3 &hist);
|
|---|
| 31 | void PrintTH3Error(const TH3 &hist);
|
|---|
| 32 | void PrintTH2Content(const TH2 &hist);
|
|---|
| 33 | void PrintTH2Error(const TH2 &hist);
|
|---|
| 34 | void PrintTH1Content(const TH1 &hist);
|
|---|
| 35 | void PrintTH1Error(const TH1 &hist);
|
|---|
| 36 | void CopyCol(TMatrixD &m, const TH1 &h, Int_t col);
|
|---|
| 37 | void CopyCol(TH1 &h, const TMatrixD &m, Int_t col);
|
|---|
| 38 | void CopyH2M(TMatrixD &m, const TH2 &h);
|
|---|
| 39 | void CopySqr(TMatrixD &m, const TH1 &h);
|
|---|
| 40 | Double_t GetMatrixSumRow(const TMatrixD &m, Int_t row);
|
|---|
| 41 | Double_t GetMatrixSumDiag(const TMatrixD &m);
|
|---|
| 42 | Double_t GetMatrixSumCol(const TMatrixD &m, Int_t col);
|
|---|
| 43 | Double_t GetMatrixSum(const TMatrixD &m);
|
|---|
| 44 |
|
|---|
| 45 | public:
|
|---|
| 46 | TString bintitle;
|
|---|
| 47 |
|
|---|
| 48 | UInt_t fNa; // Number of bins in the distribution to be unfolded
|
|---|
| 49 | UInt_t fNb; // Number of bins in the unfolded distribution
|
|---|
| 50 |
|
|---|
| 51 | TMatrixD fMigrat; // migration matrix (fNa, fNb)
|
|---|
| 52 | TMatrixD fMigraterr2;// error**2 of migration matrix (fNa, fNb)
|
|---|
| 53 |
|
|---|
| 54 | TMatrixD fMigOrig; // original migration matrix (fNa, fNb)
|
|---|
| 55 | TMatrixD fMigOrigerr2;// error**2 oforiginal migr. matrix (fNa, fNb)
|
|---|
| 56 |
|
|---|
| 57 | TMatrixD fMigSmoo; // smoothed migration matrix M (fNa, fNb)
|
|---|
| 58 | TMatrixD fMigSmooerr2;// error**2 of smoothed migr. matrix (fNa, fNb)
|
|---|
| 59 | TMatrixD fMigChi2; // chi2 contributions for smoothing (fNa, fNb)
|
|---|
| 60 |
|
|---|
| 61 | TMatrixD fVa; // distribution to be unfolded (fNa)
|
|---|
| 62 | TMatrixD fVacov; // error matrix of fVa (fNa, fNa)
|
|---|
| 63 | TMatrixD fVacovInv; // inverse of fVacov (fNa, fNa)
|
|---|
| 64 | Double_t fSpurVacov; // Spur of fVacov
|
|---|
| 65 |
|
|---|
| 66 | // UInt_t fVaevents; // total number of events
|
|---|
| 67 | UInt_t fVapoints; // number of significant measurements
|
|---|
| 68 |
|
|---|
| 69 | TMatrixD fVb; // unfolded distribution (fNb)
|
|---|
| 70 | TMatrixD fVbcov; // error matrix of fVb (fNb, fNb)
|
|---|
| 71 |
|
|---|
| 72 | TMatrixD fVEps; // prior distribution (fNb)
|
|---|
| 73 | TMatrixDColumn fVEps0;
|
|---|
| 74 |
|
|---|
| 75 | Double_t fW; // weight
|
|---|
| 76 | Double_t fWbest; // best weight
|
|---|
| 77 | Int_t ixbest;
|
|---|
| 78 |
|
|---|
| 79 | TMatrixD fResult; // unfolded distribution and errors (fNb, 5)
|
|---|
| 80 | TMatrixD fChi2; // chisquared contribution (fNa, 1)
|
|---|
| 81 |
|
|---|
| 82 | Double_t fChisq; // total chisquared
|
|---|
| 83 | Double_t fNdf; // number of degrees of freedom
|
|---|
| 84 | Double_t fProb; // chisquared probability
|
|---|
| 85 |
|
|---|
| 86 | TMatrixD G; // G = M * M(transposed) (fNa, fNa)
|
|---|
| 87 | TVectorD EigenValue; // vector of eigenvalues lambda of G (fNa)
|
|---|
| 88 | TMatrixD Eigen; // matrix of eigen vectors of G (fNa, fNa)
|
|---|
| 89 | Double_t RankG; // rank of G
|
|---|
| 90 | Double_t tau; // 1 / lambda_max
|
|---|
| 91 | Double_t EpsLambda;
|
|---|
| 92 |
|
|---|
| 93 | // quantities stored for each weight :
|
|---|
| 94 | TVectorD SpSig; // Spur of covariance matrix of fVbcov
|
|---|
| 95 | TVectorD SpAR; // effective rank of G^tilde
|
|---|
| 96 | TVectorD chisq; // chi squared (measures agreement between
|
|---|
| 97 | // fVa and the folded fVb)
|
|---|
| 98 | TVectorD SecDer; // regularization term = sum of (2nd der.)**2
|
|---|
| 99 | TVectorD ZerDer; // regularization term = sum of (fVb)**2
|
|---|
| 100 | TVectorD Entrop; // regularization term = reduced cross-entropy
|
|---|
| 101 | TVectorD DAR2; //
|
|---|
| 102 | TVectorD Dsqbar; //
|
|---|
| 103 |
|
|---|
| 104 | Double_t SpurAR;
|
|---|
| 105 | Double_t SpurSigma;
|
|---|
| 106 | Double_t SecDeriv;
|
|---|
| 107 | Double_t ZerDeriv;
|
|---|
| 108 | Double_t Entropy;
|
|---|
| 109 | Double_t DiffAR2;
|
|---|
| 110 | Double_t Chisq;
|
|---|
| 111 | Double_t D2bar;
|
|---|
| 112 |
|
|---|
| 113 | TMatrixD Chi2;
|
|---|
| 114 |
|
|---|
| 115 | //
|
|---|
| 116 |
|
|---|
| 117 | // plots versus weight
|
|---|
| 118 | Int_t Nix;
|
|---|
| 119 | Double_t xmin;
|
|---|
| 120 | Double_t xmax;
|
|---|
| 121 | Double_t dlogx;
|
|---|
| 122 |
|
|---|
| 123 | TH1D *hBchisq;
|
|---|
| 124 | TH1D *hBSpAR;
|
|---|
| 125 | TH1D *hBDSpAR;
|
|---|
| 126 | TH1D *hBSpSig;
|
|---|
| 127 | TH1D *hBDSpSig;
|
|---|
| 128 | TH1D *hBSecDeriv;
|
|---|
| 129 | TH1D *hBDSecDeriv;
|
|---|
| 130 | TH1D *hBZerDeriv;
|
|---|
| 131 | TH1D *hBDZerDeriv;
|
|---|
| 132 | TH1D *hBEntropy;
|
|---|
| 133 | TH1D *hBDEntropy;
|
|---|
| 134 | TH1D *hBDAR2;
|
|---|
| 135 | TH1D *hBD2bar;
|
|---|
| 136 |
|
|---|
| 137 | //
|
|---|
| 138 | TH1D *hEigen;
|
|---|
| 139 |
|
|---|
| 140 | // plots for the best solution
|
|---|
| 141 | TH2D *fhmig;
|
|---|
| 142 | TH2D *shmig;
|
|---|
| 143 | TH2D *shmigChi2;
|
|---|
| 144 |
|
|---|
| 145 | TH1D *fhb0;
|
|---|
| 146 |
|
|---|
| 147 | TH1D *fha;
|
|---|
| 148 |
|
|---|
| 149 | TH1D *hprior;
|
|---|
| 150 |
|
|---|
| 151 | TH1D *hb;
|
|---|
| 152 |
|
|---|
| 153 | Double_t CalcSpurSigma(TMatrixD &T, Double_t norm=1);
|
|---|
| 154 |
|
|---|
| 155 | // -----------------------------------------------------------------------
|
|---|
| 156 | //
|
|---|
| 157 | // Constructor
|
|---|
| 158 | // copy histograms into matrices
|
|---|
| 159 | //
|
|---|
| 160 | MUnfold(TH1D &ha, TH2D &hacov, TH2D &hmig);
|
|---|
| 161 | ~MUnfold();
|
|---|
| 162 |
|
|---|
| 163 | void SetPriorConstant();
|
|---|
| 164 | Bool_t SetPriorRebin(TH1D &ha);
|
|---|
| 165 | Bool_t SetPriorInput(TH1D &hpr);
|
|---|
| 166 | Bool_t SetPriorPower(Double_t gamma);
|
|---|
| 167 | Bool_t SetInitialWeight(Double_t &weight);
|
|---|
| 168 | void PrintResults();
|
|---|
| 169 | Bool_t Schmelling(TH1D &hb0);
|
|---|
| 170 | void SchmellCore(Int_t full, Double_t &xiter, TMatrixD &gamma,
|
|---|
| 171 | TMatrixD &dgamma, Double_t &dga2);
|
|---|
| 172 | Bool_t SmoothMigrationMatrix(TH2D &hmigorig);
|
|---|
| 173 | Bool_t Tikhonov2(TH1D &hb0);
|
|---|
| 174 | Bool_t Bertero(TH1D &hb0);
|
|---|
| 175 | Bool_t BertCore(Double_t &xiter);
|
|---|
| 176 | Bool_t CalculateG();
|
|---|
| 177 | Bool_t SelectBestWeight();
|
|---|
| 178 | Bool_t DrawPlots();
|
|---|
| 179 | Bool_t CallMinuit(
|
|---|
| 180 | void (*fcnx)(Int_t &, Double_t *, Double_t &, Double_t *, Int_t),
|
|---|
| 181 | UInt_t npar, char name[20][100],
|
|---|
| 182 | Double_t vinit[20], Double_t step[20],
|
|---|
| 183 | Double_t limlo[20], Double_t limup[20], Int_t fix[20]);
|
|---|
| 184 | TMatrixD &GetVb() { return fVb; }
|
|---|
| 185 | TMatrixD &GetVbcov() { return fVbcov; }
|
|---|
| 186 | TMatrixD &GetResult() { return fResult; }
|
|---|
| 187 | TMatrixD &GetChi2() { return fChi2; }
|
|---|
| 188 | Double_t &GetChisq() { return fChisq; }
|
|---|
| 189 | Double_t &GetNdf() { return fNdf; }
|
|---|
| 190 | Double_t &GetProb() { return fProb; }
|
|---|
| 191 | TMatrixD &GetMigSmoo() { return fMigSmoo; }
|
|---|
| 192 | TMatrixD &GetMigSmooerr2() { return fMigSmooerr2; }
|
|---|
| 193 | TMatrixD &GetMigChi2() { return fMigChi2; }
|
|---|
| 194 |
|
|---|
| 195 | ClassDef(MUnfold, 1)
|
|---|
| 196 |
|
|---|
| 197 | };
|
|---|
| 198 |
|
|---|
| 199 | #endif
|
|---|
