Utils/src/PUReweightingMulti.cc

// $Id: PUReweighting.cc,v 1.1 2011/06/01 12:36:59 mzanetti Exp $


#include "MitPhysics/Utils/interface/PUReweightingMulti.h"
#include <TAxis.h>

ClassImp(mithep::PUReweightingMulti)

using namespace mithep;


PUReweightingMulti::PUReweightingMulti(const THnSparse *sparse, int maxn, double minprob) :
fSparse(sparse),
fNsBins(0),
fNDim(0),
fMaxN(maxn+1),
fMinProb(minprob),
fWeights(0),
fIndices(0),
fCache(0)
{

  Setup();

}

PUReweightingMulti::~PUReweightingMulti() {
  
  if (fWeights) delete [] fWeights;
  if (fIndices) delete [] fIndices;
  if (fCache) delete [] fCache;
  if (fOffsets) delete [] fOffsets;
  if (fIndOffsets) delete [] fIndOffsets;
  if (fFactors) delete [] fFactors;

  
}

void PUReweightingMulti::Setup() {
  
  //cache all values needed to quickly computed probability density as sum over bins of the sparse histogram
  
  fNsBins = fSparse->GetNbins();
  fNDim = fSparse->GetNdimensions();
  
  fOffsets = new Int_t[fNDim];
  fIndOffsets = new Int_t[fNDim];
  fFactors = new Double_t[fNDim];
  
  fWeights = new Double_t[fNsBins];
  fIndices = new UShort_t[fNDim*fNsBins];

  int nbinstotal = 0;
  int *cacheoffset = new int[fNDim];
  for (int idim=0; idim<fNDim; ++idim) {
    cacheoffset[idim] = nbinstotal;
    nbinstotal += fSparse->GetAxis(idim)->GetNbins()+1;
    fIndOffsets[idim] = fNsBins*idim;
  }
  fNCacheBins = nbinstotal;
  fCache = new Double_t[fMaxN*fNCacheBins];
  for (int idim=0; idim<fNDim; ++idim) {
    int nbins = fSparse->GetAxis(idim)->GetNbins()+1;
    for (int ibin=0; ibin<nbins; ++ibin) {
      double lambda = (ibin>0) ? exp(fSparse->GetAxis(idim)->GetBinCenter(ibin)) : 0.0;
      for (int n=0; n<fMaxN; ++n) {
        int globalbin = cacheoffset[idim] + ibin;
        //fCache[fMaxN*globalbin + n] = TMath::Poisson(n,lambda);
        double val = TMath::Poisson(n,lambda);
        fCache[fNCacheBins*n + globalbin] = (val > fMinProb) ? val : 0.0;
      }
    }
  }  
  
  int *coords = new int[fNDim];
  double sumweights = 0;
  for (int i=0; i<fNsBins; ++i) {
    sumweights += fSparse->GetBinContent(i,coords);
  }

  for (int i=0; i<fNsBins; ++i) {
    double weight = fSparse->GetBinContent(i,coords);
    fWeights[i] = weight/sumweights;
    for (int idim=0; idim<fNDim; ++idim) {
      //fIndices[fNDim*i + idim] = cacheoffset[idim] + coords[idim];
      fIndices[fNsBins*idim + i] = cacheoffset[idim] + coords[idim];
    }
  }
  
  delete[] cacheoffset;
  delete[] coords;

}

double PUReweightingMulti::TargetPdf(int *npus, int ndim) const {
    
  assert(ndim<=fNDim);
  
  for (int idim=0; idim<ndim; ++idim) {
    if (npus[idim] >= fMaxN) return 0.0;
    fOffsets[idim] = fNCacheBins*npus[idim];
  }
  
  double totalweight = 0.0;
  for (int i=0; i<fNsBins; ++i) {
    bool nonzero = true;
    for (int idim=0; idim<ndim; ++idim) {
      fFactors[idim] = fCache[fOffsets[idim]+fIndices[fIndOffsets[idim]+i]];
      if ( !(fFactors[idim]>0.0) ) {
        nonzero = false;
        break;
      }
    }
    if (nonzero) {
      double weight = fWeights[i];
      for (int idim=0; idim<ndim; ++idim) {
        weight *= fFactors[idim];
      }
      totalweight += weight;
    }
  }
  
  return totalweight;
}

TH3D *PUReweightingMulti::Target3D() {
 assert(fNDim>=3);
  
 int np[3];
  
 TH3D *target3d = new TH3D("target3d","",fMaxN,-0.5,fMaxN-0.5,fMaxN,-0.5,fMaxN-0.5,fMaxN,-0.5,fMaxN-0.5);
 for (int i=0; i<fMaxN; ++i) {
   np[0] = i;
   for (int j=0; j<fMaxN; ++j) {
     np[1] = j;
     for (int k=0; k<fMaxN; ++k) {
       np[2] = k;
       target3d->Fill(i,j,k,TargetPdf(np,3));
     }
   }
 }
 return target3d;

}

TH3D *PUReweightingMulti::Weights3D(const TH3D *source3d) {
 assert(fNDim>=3);
  
 TH3D *source3dclone = (TH3D*)source3d->Clone();
 source3dclone->Scale(1.0/source3dclone->GetSumOfWeights());
  
 int np[3];
 TH3D *target3d = new TH3D("target3dtemp","",fMaxN,-0.5,fMaxN-0.5,fMaxN,-0.5,fMaxN-0.5,fMaxN,-0.5,fMaxN-0.5);
 for (int i=0; i<fMaxN; ++i) {
   np[0] = i;
   for (int j=0; j<fMaxN; ++j) {
     np[1] = j;
     for (int k=0; k<fMaxN; ++k) {
       np[2] = k;
       if (source3dclone->GetBinContent(source3dclone->FindFixBin(i,j,k))>0.0) {
         target3d->Fill(i,j,k,TargetPdf(np,3));
       }
       else {
         target3d->Fill(i,j,k,0.0);
       }
     }
   }
 }
 target3d->Scale(1.0/target3d->GetSumOfWeights());
 
 TH3D *weights3d = new TH3D(*target3d / *source3dclone);
 delete source3dclone;
 return weights3d;

}
Revision:	1.1
Committed:	Fri Jul 22 18:21:46 2011 UTC (13 years, 9 months ago) by bendavid
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Mit_029c, Mit_029b, Mit_029a, Mit_028a, Mit_028, Mit_027, Mit_027a, Mit_025e, Mit_025d, Mit_025c, Mit_025b, Mit_025a, Mit_025, Mit_025pre2, Mit_024b, Mit_025pre1, Mit_024a, Mit_024, HEAD
Log Message:	add experimental multidimensional pileup reweighting
#	User	Rev	Content
1	bendavid	1.1	// $Id: PUReweighting.cc,v 1.1 2011/06/01 12:36:59 mzanetti Exp $
2
3
4			#include "MitPhysics/Utils/interface/PUReweightingMulti.h"
5			#include <TAxis.h>
6
7			ClassImp(mithep::PUReweightingMulti)
8
9			using namespace mithep;
10
11
12			PUReweightingMulti::PUReweightingMulti(const THnSparse *sparse, int maxn, double minprob) :
13			fSparse(sparse),
14			fNsBins(0),
15			fNDim(0),
16			fMaxN(maxn+1),
17			fMinProb(minprob),
18			fWeights(0),
19			fIndices(0),
20			fCache(0)
21			{
22
23			Setup();
24
25			}
26
27			PUReweightingMulti::~PUReweightingMulti() {
28
29			if (fWeights) delete [] fWeights;
30			if (fIndices) delete [] fIndices;
31			if (fCache) delete [] fCache;
32			if (fOffsets) delete [] fOffsets;
33			if (fIndOffsets) delete [] fIndOffsets;
34			if (fFactors) delete [] fFactors;
35
36
37			}
38
39			void PUReweightingMulti::Setup() {
40
41			//cache all values needed to quickly computed probability density as sum over bins of the sparse histogram
42
43			fNsBins = fSparse->GetNbins();
44			fNDim = fSparse->GetNdimensions();
45
46			fOffsets = new Int_t[fNDim];
47			fIndOffsets = new Int_t[fNDim];
48			fFactors = new Double_t[fNDim];
49
50			fWeights = new Double_t[fNsBins];
51			fIndices = new UShort_t[fNDim*fNsBins];
52
53			int nbinstotal = 0;
54			int *cacheoffset = new int[fNDim];
55			for (int idim=0; idim<fNDim; ++idim) {
56			cacheoffset[idim] = nbinstotal;
57			nbinstotal += fSparse->GetAxis(idim)->GetNbins()+1;
58			fIndOffsets[idim] = fNsBins*idim;
59			}
60			fNCacheBins = nbinstotal;
61			fCache = new Double_t[fMaxN*fNCacheBins];
62			for (int idim=0; idim<fNDim; ++idim) {
63			int nbins = fSparse->GetAxis(idim)->GetNbins()+1;
64			for (int ibin=0; ibin<nbins; ++ibin) {
65			double lambda = (ibin>0) ? exp(fSparse->GetAxis(idim)->GetBinCenter(ibin)) : 0.0;
66			for (int n=0; n<fMaxN; ++n) {
67			int globalbin = cacheoffset[idim] + ibin;
68			//fCache[fMaxN*globalbin + n] = TMath::Poisson(n,lambda);
69			double val = TMath::Poisson(n,lambda);
70			fCache[fNCacheBins*n + globalbin] = (val > fMinProb) ? val : 0.0;
71			}
72			}
73			}
74
75			int *coords = new int[fNDim];
76			double sumweights = 0;
77			for (int i=0; i<fNsBins; ++i) {
78			sumweights += fSparse->GetBinContent(i,coords);
79			}
80
81			for (int i=0; i<fNsBins; ++i) {
82			double weight = fSparse->GetBinContent(i,coords);
83			fWeights[i] = weight/sumweights;
84			for (int idim=0; idim<fNDim; ++idim) {
85			//fIndices[fNDim*i + idim] = cacheoffset[idim] + coords[idim];
86			fIndices[fNsBins*idim + i] = cacheoffset[idim] + coords[idim];
87			}
88			}
89
90			delete[] cacheoffset;
91			delete[] coords;
92
93			}
94
95			double PUReweightingMulti::TargetPdf(int *npus, int ndim) const {
96
97			assert(ndim<=fNDim);
98
99			for (int idim=0; idim<ndim; ++idim) {
100			if (npus[idim] >= fMaxN) return 0.0;
101			fOffsets[idim] = fNCacheBins*npus[idim];
102			}
103
104			double totalweight = 0.0;
105			for (int i=0; i<fNsBins; ++i) {
106			bool nonzero = true;
107			for (int idim=0; idim<ndim; ++idim) {
108			fFactors[idim] = fCache[fOffsets[idim]+fIndices[fIndOffsets[idim]+i]];
109			if ( !(fFactors[idim]>0.0) ) {
110			nonzero = false;
111			break;
112			}
113			}
114			if (nonzero) {
115			double weight = fWeights[i];
116			for (int idim=0; idim<ndim; ++idim) {
117			weight *= fFactors[idim];
118			}
119			totalweight += weight;
120			}
121			}
122
123			return totalweight;
124			}
125
126			TH3D *PUReweightingMulti::Target3D() {
127			assert(fNDim>=3);
128
129			int np[3];
130
131			TH3D *target3d = new TH3D("target3d","",fMaxN,-0.5,fMaxN-0.5,fMaxN,-0.5,fMaxN-0.5,fMaxN,-0.5,fMaxN-0.5);
132			for (int i=0; i<fMaxN; ++i) {
133			np[0] = i;
134			for (int j=0; j<fMaxN; ++j) {
135			np[1] = j;
136			for (int k=0; k<fMaxN; ++k) {
137			np[2] = k;
138			target3d->Fill(i,j,k,TargetPdf(np,3));
139			}
140			}
141			}
142			return target3d;
143
144			}
145
146			TH3D PUReweightingMulti::Weights3D(const TH3D source3d) {
147			assert(fNDim>=3);
148
149			TH3D source3dclone = (TH3D)source3d->Clone();
150			source3dclone->Scale(1.0/source3dclone->GetSumOfWeights());
151
152			int np[3];
153			TH3D *target3d = new TH3D("target3dtemp","",fMaxN,-0.5,fMaxN-0.5,fMaxN,-0.5,fMaxN-0.5,fMaxN,-0.5,fMaxN-0.5);
154			for (int i=0; i<fMaxN; ++i) {
155			np[0] = i;
156			for (int j=0; j<fMaxN; ++j) {
157			np[1] = j;
158			for (int k=0; k<fMaxN; ++k) {
159			np[2] = k;
160			if (source3dclone->GetBinContent(source3dclone->FindFixBin(i,j,k))>0.0) {
161			target3d->Fill(i,j,k,TargetPdf(np,3));
162			}
163			else {
164			target3d->Fill(i,j,k,0.0);
165			}
166			}
167			}
168			}
169			target3d->Scale(1.0/target3d->GetSumOfWeights());
170
171			TH3D weights3d = new TH3D(target3d / *source3dclone);
172			delete source3dclone;
173			return weights3d;
174
175			}