#include "TCanvas.h"
#include "TF1.h"
#include "THnSparse.h"
+#include "TArrayD.h"
#define round(x) ((x)>=0?(long)((x)+0.5):(long)((x)-0.5))
fResolutionType(kParam),
fDeltaPt(0x0),
fhDeltaPt(0x0),
+ fh1MeasuredTruncated(0x0),
+ fh2DetectorResponse(0x0),
+ fh2DetectorResponseRebin(0x0),
+ fhEfficiencyDet(0x0),
+ fh2ResponseMatrixCombinedFineFull(0x0),
+ fh2ResponseMatrixCombinedFull(0x0),
+ fh2ResponseMatrixCombined(0x0),
+ fhEfficiencyCombined(0x0),
fDimensions(1),
fDimRec(0),
fDimGen(1),
fResolutionType(obj.fResolutionType),
fDeltaPt(obj.fDeltaPt),
fhDeltaPt(obj.fhDeltaPt),
+ fh1MeasuredTruncated(obj.fh1MeasuredTruncated),
+ fh2DetectorResponse(obj.fh2DetectorResponse),
+ fh2DetectorResponseRebin(obj.fh2DetectorResponseRebin),
+ fhEfficiencyDet(obj.fhEfficiencyDet),
+ fh2ResponseMatrixCombinedFineFull(obj.fh2ResponseMatrixCombinedFineFull),
+ fh2ResponseMatrixCombinedFull(obj.fh2ResponseMatrixCombinedFull),
+ fh2ResponseMatrixCombined(obj.fh2ResponseMatrixCombined),
+ fhEfficiencyCombined(obj.fhEfficiencyCombined),
fDimensions(obj.fDimensions),
fDimRec(obj.fDimRec),
fDimGen(obj.fDimGen),
// Assignment
if(&other == this) return *this;
AliAnaChargedJetResponseMaker::operator=(other);
- fDebug = other.fDebug;
- fResolutionType = other.fResolutionType;
- fDeltaPt = other.fDeltaPt;
- fhDeltaPt = other.fhDeltaPt;
- fDimensions = other.fDimensions;
- fDimRec = other.fDimRec;
- fDimGen = other.fDimGen;
- fPtMin = other.fPtMin;
- fPtMax = other.fPtMax;
- fNbins = other.fNbins;
- fBinArrayPtRec = other.fBinArrayPtRec;
- fPtMeasured = other.fPtMeasured;
- fEffFlat = other.fEffFlat;
- fEfficiency = other.fEfficiency;
- fEfficiencyFine = other.fEfficiencyFine;
- fResponseMatrix = other.fResponseMatrix;
- fResponseMatrixFine = other.fResponseMatrixFine;
- fPtMinUnfolded = other.fPtMinUnfolded;
- fPtMaxUnfolded = other.fPtMaxUnfolded;
- fPtMaxUnfoldedHigh = other.fPtMaxUnfoldedHigh;
- fBinWidthFactorUnfolded = other.fBinWidthFactorUnfolded;
- fSkipBinsUnfolded = other.fSkipBinsUnfolded;
- fExtraBinsUnfolded = other.fExtraBinsUnfolded;
- fbVariableBinning = other.fbVariableBinning;
- fPtMaxUnfVarBinning = other.fPtMaxUnfVarBinning;
- f1MergeFunction = other.f1MergeFunction;
- fFineFrac = other.fFineFrac;
- fbCalcErrors = other.fbCalcErrors;
+ fDebug = other.fDebug;
+ fResolutionType = other.fResolutionType;
+ fDeltaPt = other.fDeltaPt;
+ fhDeltaPt = other.fhDeltaPt;
+ fh1MeasuredTruncated = other.fh1MeasuredTruncated;
+ fh2DetectorResponse = other.fh2DetectorResponse;
+ fh2DetectorResponseRebin = other.fh2DetectorResponseRebin;
+ fhEfficiencyDet = other.fhEfficiencyDet;
+ fh2ResponseMatrixCombinedFineFull = other.fh2ResponseMatrixCombinedFineFull;
+ fh2ResponseMatrixCombinedFull = other.fh2ResponseMatrixCombinedFull;
+ fh2ResponseMatrixCombined = other.fh2ResponseMatrixCombined;
+ fhEfficiencyCombined = other.fhEfficiencyCombined;
+ fDimensions = other.fDimensions;
+ fDimRec = other.fDimRec;
+ fDimGen = other.fDimGen;
+ fPtMin = other.fPtMin;
+ fPtMax = other.fPtMax;
+ fNbins = other.fNbins;
+ fBinArrayPtRec = other.fBinArrayPtRec;
+ fPtMeasured = other.fPtMeasured;
+ fEffFlat = other.fEffFlat;
+ fEfficiency = other.fEfficiency;
+ fEfficiencyFine = other.fEfficiencyFine;
+ fResponseMatrix = other.fResponseMatrix;
+ fResponseMatrixFine = other.fResponseMatrixFine;
+ fPtMinUnfolded = other.fPtMinUnfolded;
+ fPtMaxUnfolded = other.fPtMaxUnfolded;
+ fPtMaxUnfoldedHigh = other.fPtMaxUnfoldedHigh;
+ fBinWidthFactorUnfolded = other.fBinWidthFactorUnfolded;
+ fSkipBinsUnfolded = other.fSkipBinsUnfolded;
+ fExtraBinsUnfolded = other.fExtraBinsUnfolded;
+ fbVariableBinning = other.fbVariableBinning;
+ fPtMaxUnfVarBinning = other.fPtMaxUnfVarBinning;
+ f1MergeFunction = other.f1MergeFunction;
+ fFineFrac = other.fFineFrac;
+ fbCalcErrors = other.fbCalcErrors;
return *this;
}
if(fbVariableBinning)
binWidthUnfLowPt = binWidthUnf*0.5;
- if(fExtraBinsUnfolded>0) {
- fPtMaxUnfolded = fPtMax+(double)(fExtraBinsUnfolded)*binWidthUnf;
- nbins[fDimGen]+=fExtraBinsUnfolded;
- }
+ fPtMaxUnfolded = fPtMax+(double)(fExtraBinsUnfolded)*binWidthUnf;
+ nbins[fDimGen]+=fExtraBinsUnfolded;
+
printf("fPtMinMeas: %f fPtMaxMeas: %f\n",fPtMin,fPtMax);
printf("binWidthMeas: %f binWidthUnf: %f fBinWidthFactorUnfolded: %d\n",binWidthMeas,binWidthUnf,fBinWidthFactorUnfolded);
nbinsFine[fDimGen] = fResponseMatrix->GetAxis(fDimRec)->GetNbins()*fFineFrac;
xminFine[fDimRec] = TMath::Min(fPtMin,0.);
xminFine[fDimGen] = TMath::Min(fPtMin,0.);
- xmaxFine[fDimRec] = fResponseMatrix->GetAxis(fDimGen)->GetXmax()+40.;
- xmaxFine[fDimGen] = fResponseMatrix->GetAxis(fDimGen)->GetXmax()+40.;
+ xmaxFine[fDimRec] = fResponseMatrix->GetAxis(fDimGen)->GetXmax();//+40.;
+ xmaxFine[fDimGen] = fResponseMatrix->GetAxis(fDimGen)->GetXmax();//+40.;
pTarrayFine[fDimRec] = 0.;
pTarrayFine[fDimGen] = 0.;
Double_t binWidthFine[2];
binWidthFine[fDimRec] = binWidth[fDimRec]/((double)fFineFrac);
binWidthFine[fDimGen] = binWidthFine[fDimRec]*(double)fBinWidthFactorUnfolded;
+ // binWidthFine[fDimRec] = binWidthFine[fDimGen];
nbinsFine[fDimRec] = (int)((xmaxFine[fDimRec]-xminFine[fDimRec])/binWidthFine[fDimRec]); //adjust nbins to bin width
nbinsFine[fDimGen] = (int)((xmaxFine[fDimGen]-xminFine[fDimGen])/binWidthFine[fDimGen]); //adjust nbins to bin width
double sumyield = 0.;
double sumyielderror2 = 0.;
- int ipt[2] = {0.,0.};
- int iptMerged[2] = {0.,0.};
- int ierror[2] = {0.,0.};
+ int ipt[2] = {0,0};
+ int iptMerged[2] = {0,0};
+ int ierror[2] = {0,0};
Int_t check = 0;
Double_t pTgen, pTrec;
//fill merged response matrix
- //find bin in fine RM correspoinding to mimimum/maximum bin of merged RM on rec axis
+ //find bin in fine RM corresponding to mimimum/maximum bin of merged RM on rec axis
int ixMin = fResponseMatrixFine->GetAxis(fDimRec)->FindBin(fResponseMatrix->GetAxis(fDimRec)->GetXmin());
int ixMax = fResponseMatrixFine->GetAxis(fDimRec)->FindBin(fResponseMatrix->GetAxis(fDimRec)->GetXmax());
}
//--------------------------------------------------------------------------------------------------------------------------------------------------
-TH2* AliAnaChargedJetResponseMaker::MakeResponseMatrixRebin(TH2 *hRMFine, TH2 *hRM) {
+TH2* AliAnaChargedJetResponseMaker::MakeResponseMatrixRebin(TH2 *hRMFine, TH2 *hRM, Bool_t useFunctionWeight) {
//
// Rebin matrix hRMFine to dimensions of hRM
- // function returns matrix in TH2D format with dimensions from hRM
+ // function returns matrix in TH2D format (hRM2) with dimensions from hRM
//
TH2 *hRM2 = (TH2*)hRM->Clone("hRM2");
hRM2->Reset();
+ if(useFunctionWeight) cout << "Use function to do weighting" << endl;
+
//First normalize columns of input
const Int_t nbinsNorm = hRM2->GetNbinsX();
- cout << "nbinsNorm: " << nbinsNorm << endl;
+ if(fDebug) cout << "nbinsNorm: " << nbinsNorm << endl;
TArrayF *normVector = new TArrayF(nbinsNorm);
Int_t binxLowFine = hRMFine->GetXaxis()->FindBin(xLow);
Int_t binxUpFine = hRMFine->GetXaxis()->FindBin(xUp)-1;
//cout << "xLowFine: " << hRMFine->GetXaxis()->GetBinLowEdge(binxLowFine) << "\txUpFine: " << hRMFine->GetXaxis()->GetBinUpEdge(binxUpFine) << endl;
- normVector->SetAt(hRMFine->Integral(binxLowFine,binxUpFine,1,hRMFine->GetYaxis()->GetNbins()),ix-1);
- // if(fDebug) cout << "ix norm: " << normVector->At(ix-1) << endl;
+ if(useFunctionWeight)
+ normVector->SetAt(f1MergeFunction->Integral(xLow,xUp),ix-1);
+ else
+ normVector->SetAt(hRMFine->Integral(binxLowFine,binxUpFine,1,hRMFine->GetYaxis()->GetNbins()),ix-1);
+ //if(fDebug) cout << "ix norm: " << normVector->At(ix-1) << endl;
}
Double_t content, oldcontent = 0.;
//if(fDebug) cout << "ixNew: " << ixNew << " " << xvalLo << " iyNew: " << iyNew << " " << yvalLo << " content: " << content << " oldContent: " << oldcontent << " newContent: " << oldcontent+content << endl;
Double_t weight = 1.;
- if(normVector->At(ixNew-1)>0.) weight = 1./normVector->At(ixNew-1);
+ if(normVector->At(ixNew-1)>0.) {
+ if(useFunctionWeight)
+ weight = f1MergeFunction->Integral(xvalLo,xvalUp)/normVector->At(ixNew-1);
+ else
+ weight = 1./normVector->At(ixNew-1);
+ }
hRM2->SetBinContent(ixNew,iyNew,oldcontent+content*weight);
}
}
}
+//--------------------------------------------------------------------------------------------------------------------------------------------------
+TH2* AliAnaChargedJetResponseMaker::CreateTruncated2DHisto(TH2 *h2, Double_t xmin, Double_t xmax, Double_t ymin, Double_t ymax) {
+
+ //
+ // Limit axis range of 2D histogram
+ //
+
+ Int_t binMinXh2 = h2->GetXaxis()->FindBin(xmin);
+ if(h2->GetXaxis()->GetBinLowEdge(binMinXh2) < xmin ) binMinXh2++;
+ if(h2->GetXaxis()->GetBinLowEdge(binMinXh2) > xmin ) binMinXh2--;
+
+ Int_t binMinYh2 = h2->GetYaxis()->FindBin(ymin);
+ if(h2->GetYaxis()->GetBinLowEdge(binMinYh2) < ymin ) binMinYh2++;
+ if(h2->GetYaxis()->GetBinLowEdge(binMinYh2) > ymin ) binMinYh2--;
+
+ Int_t binMaxXh2 = h2->GetXaxis()->FindBin(xmax);
+ if(h2->GetXaxis()->GetBinUpEdge(binMaxXh2) < xmax ) binMaxXh2++;
+ if(h2->GetXaxis()->GetBinUpEdge(binMaxXh2) > xmax ) binMaxXh2--;
+
+ Int_t binMaxYh2 = h2->GetYaxis()->FindBin(ymax);
+ if(h2->GetYaxis()->GetBinUpEdge(binMaxYh2) < ymax ) binMaxYh2++;
+ if(h2->GetYaxis()->GetBinUpEdge(binMaxYh2) > ymax ) binMaxYh2--;
+
+ Int_t nbinsX = binMaxXh2-binMinXh2+1;
+ Int_t nbinsY = binMaxYh2-binMinYh2+1;
+
+ Double_t *binsX = new Double_t[nbinsX+1];
+ Double_t *binsY = new Double_t[nbinsY+1];
+
+ for(int ix=1; ix<=nbinsX; ix++)
+ binsX[ix-1] = h2->GetXaxis()->GetBinLowEdge(binMinXh2+ix-1);
+ binsX[nbinsX] = h2->GetXaxis()->GetBinUpEdge(binMaxXh2);
+
+ for(int iy=1; iy<=nbinsY; iy++)
+ binsY[iy-1] = h2->GetYaxis()->GetBinLowEdge(binMinYh2+iy-1);
+ binsY[nbinsY] = h2->GetYaxis()->GetBinUpEdge(binMaxYh2);
+
+ TH2 *h2Lim = new TH2D("h2Lim","h2Lim",nbinsX,binsX,nbinsY,binsY);
+
+ for(int ix=1; ix<=nbinsX; ix++) {
+ // cout << "ix: " << ix << " " << binsX[ix] << endl;
+ for(int iy=1; iy<=nbinsY; iy++) {
+ // cout << "ix: " << ix << " " << binsX[ix] << "\tiy: " << iy << " " << binsY[iy] << endl;
+
+ double content = h2->GetBinContent(binMinXh2+ix-1,binMinYh2+iy-1);
+ double error = h2->GetBinContent(binMinXh2+ix-1,binMinYh2+iy-1);
+ h2Lim->SetBinContent(ix,iy,content);
+ if(fbCalcErrors) h2Lim->SetBinError(ix,iy,error);
+
+ }
+ }
+
+
+
+ return h2Lim;
+}
+
+//--------------------------------------------------------------------------------------------------------------------------------------------------
+TH2* AliAnaChargedJetResponseMaker::TruncateAxisRangeResponseMatrix(TH2 *hRMOrig, Double_t xmin, Double_t xmax, Double_t ymin, Double_t ymax) {
+
+ //
+ // Limit axis range of response matrix without changing normalization
+ //
+
+ //TH2 *hRMLimit
+ //TH2 *hRMLimit2 = (TH2*)hRMLimit->Clone("hRMLimit2");
+
+ TH2 *hRMLimit2 = CreateTruncated2DHisto(hRMOrig, xmin, xmax, ymin, ymax);
+ hRMLimit2->Reset();
+
+ double binCent[2] = {0.,0.};
+ double content = 0.;
+ double error = 0.;
+ Int_t binOrig[2] = {0};
+ for(int ix=1; ix<=hRMLimit2->GetXaxis()->GetNbins(); ix++) {
+ binCent[0] = hRMLimit2->GetXaxis()->GetBinCenter(ix);
+ binOrig[0] = hRMOrig->GetXaxis()->FindBin(binCent[0]);
+ for(int iy=1; iy<=hRMLimit2->GetYaxis()->GetNbins(); iy++) {
+ binCent[1] = hRMLimit2->GetYaxis()->GetBinCenter(iy);
+ binOrig[1] = hRMOrig->GetYaxis()->FindBin(binCent[1]);
+
+ content = hRMOrig->GetBinContent(binOrig[0],binOrig[1]);
+ error = hRMOrig->GetBinError(binOrig[0],binOrig[1]);
+
+ hRMLimit2->SetBinContent(ix,iy,content);
+ if(fbCalcErrors) hRMLimit2->SetBinError(ix,iy,error);
+
+ }
+ }
+
+
+ return hRMLimit2;
+
+}
+
+//--------------------------------------------------------------------------------------------------------------------------------------------------
+Bool_t AliAnaChargedJetResponseMaker::CheckInputForCombinedResponse() {
+
+ Bool_t check = kTRUE;
+
+ if(!fhDeltaPt) check = kFALSE;
+ if(!fh2DetectorResponse) check = kFALSE;
+ if(!fhEfficiencyDet) check = kFALSE;
+ if(!fh1MeasuredTruncated) check = kFALSE;
+ if(fPtMin<0. || fPtMax<0.) check = kFALSE;
+
+ return check;
+}
+
+//--------------------------------------------------------------------------------------------------------------------------------------------------
+void AliAnaChargedJetResponseMaker::MakeResponseMatrixCombined(Int_t skipBins, Int_t binWidthFactor, Int_t extraBins, Bool_t bVariableBinning, Double_t ptmin) {
+
+ //Check if all input histos are there otherwise break
+ if(!CheckInputForCombinedResponse()) {
+ printf("Not all input available..... \n");
+ return;
+ }
+
+ // Make response bkg fluctuations
+ MakeResponseMatrixJetsFineMerged(skipBins,binWidthFactor,extraBins,bVariableBinning,ptmin);
+
+ //Get TH2D with dimensions we want final RM
+ TH2D *h2ResponseMatrix = fResponseMatrix->Projection(0,1,"E");
+ h2ResponseMatrix->Reset();
+
+ //Get fine binned response matrix from bkg fluctuations
+ TH2D *h2ResponseMatrixFine = fResponseMatrixFine->Projection(0,1,"E");
+
+ // Rebin response detector effects
+ const TArrayD *arrayX = h2ResponseMatrixFine->GetXaxis()->GetXbins();
+ const Double_t *xbinsArray = arrayX->GetArray();
+
+ TH2D *h2ResponseDetTmp = new TH2D("h2ResponseDetTmp","h2ResponseDetTmp",h2ResponseMatrixFine->GetNbinsX(),xbinsArray,h2ResponseMatrixFine->GetNbinsX(),xbinsArray);
+
+ fh2DetectorResponseRebin = (TH2D*)MakeResponseMatrixRebin(fh2DetectorResponse,h2ResponseDetTmp,kFALSE);
+ fh2DetectorResponseRebin->SetName("fh2DetectorResponseRebin");
+ fh2DetectorResponseRebin->SetTitle("fh2DetectorResponseRebin");
+ fh2DetectorResponseRebin->SetXTitle("p_{T}^{jet} gen");
+ fh2DetectorResponseRebin->SetYTitle("p_{T}^{jet} rec");
+
+ // Make full combined response matrix fine binning (background flucutuations + detector effects)
+ fh2ResponseMatrixCombinedFineFull = (TH2D*)MultiplityResponseMatrices(h2ResponseMatrixFine,fh2DetectorResponseRebin);
+ fh2ResponseMatrixCombinedFineFull->SetName("fh2ResponseMatrixCombinedFineFull");
+ fh2ResponseMatrixCombinedFineFull->SetTitle("fh2ResponseMatrixCombinedFineFull");
+
+ // Rebin full combined response matrix with weighting
+ fh2ResponseMatrixCombinedFull = (TH2D*)MakeResponseMatrixRebin(fh2ResponseMatrixCombinedFineFull,h2ResponseMatrix,kTRUE);
+ fh2ResponseMatrixCombinedFull->SetName("fh2ResponseMatrixCombinedFull");
+ fh2ResponseMatrixCombinedFull->SetTitle("fh2ResponseMatrixCombinedFull");
+ fh2ResponseMatrixCombinedFull->SetXTitle("#it{p}_{T,gen}^{ch jet} (GeV/#it{c})");
+ fh2ResponseMatrixCombinedFull->SetYTitle("#it{p}_{T,rec}^{ch jet} (GeV/#it{c})");
+
+ fh2ResponseMatrixCombined = (TH2D*)TruncateAxisRangeResponseMatrix(fh2ResponseMatrixCombinedFull, h2ResponseMatrix->GetXaxis()->GetXmin(),h2ResponseMatrix->GetXaxis()->GetXmax(),fh1MeasuredTruncated->GetXaxis()->GetXmin(),fh1MeasuredTruncated->GetXaxis()->GetXmax());
+ fh2ResponseMatrixCombined->SetTitle("fh2ResponseMatrixCombined");
+ fh2ResponseMatrixCombined->SetName("fh2ResponseMatrixCombined");
+
+ fhEfficiencyCombined = (TH1D*)fh2ResponseMatrixCombined->ProjectionX("fhEfficiencyCombined");
+ fhEfficiencyCombined->Reset();
+
+ for(int i=1; i<=fh2ResponseMatrixCombined->GetNbinsX(); i++) {
+ Double_t sumyield = 0.;
+ for(int j=1; j<=fh2ResponseMatrixCombined->GetYaxis()->GetNbins(); j++) {
+ sumyield+=fh2ResponseMatrixCombined->GetBinContent(i,j);
+ }
+ Double_t kCounter = 0.;
+ Double_t kPtLowEdge = fh2ResponseMatrixCombined->GetXaxis()->GetBinLowEdge(i);
+ Double_t kPtUpEdge = fh2ResponseMatrixCombined->GetXaxis()->GetBinUpEdge(i);
+ Double_t kJetEffDet = 0.;
+
+ for(int k=1; k<=fhEfficiencyDet->GetNbinsX(); k++) {
+ if(fhEfficiencyDet->GetXaxis()->GetBinLowEdge(k)>=kPtLowEdge && fhEfficiencyDet->GetXaxis()->GetBinUpEdge(k)<=kPtUpEdge) {
+ kJetEffDet+=fhEfficiencyDet->GetBinContent(k);
+ kCounter+=1.;
+ }
+ }
+ kJetEffDet = kJetEffDet/kCounter;
+
+ if(kJetEffDet==0.) kJetEffDet=1.;
+
+ fhEfficiencyCombined->SetBinContent(i,sumyield*kJetEffDet);
+
+ }
+
+}
+
//--------------------------------------------------------------------------------------------------------------------------------------------------
TH2* AliAnaChargedJetResponseMaker::MultiplityResponseMatrices(TH2 *h2RMDeltaPt, TH2 *h2RMDetector) {
// Function assumes that generated/unfolded axis is x-axis and reconstructed is on y-axis on both respone matrices
- TH2D *h2ResponseMatrixCombined = (TH2D*)h2RMDeltaPt->Clone("h2ResponseMatrixCombined"); //h2ResponseMatrix is the bkg fluctuations RM which has the dimensions we want for the combined response matrix
+ TH2D *h2ResponseMatrixCombined = (TH2D*)h2RMDeltaPt->Clone("h2ResponseMatrixCombined"); //h2RMDeltaPt is the bkg fluctuations RM which has the dimensions we want for the combined response matrix
h2ResponseMatrixCombined->SetTitle("h2ResponseMatrixCombined");
h2ResponseMatrixCombined->SetName("h2ResponseMatrixCombined");
// RM_comb(m,t) = Sum_d RM_deltaPt(m,d)*RM_DetEffects(d,t)
if(fDebug) {
- printf("Nt=%d",h2ResponseMatrixCombined->GetNbinsX());
- printf("Nm=%d",h2ResponseMatrixCombined->GetNbinsY());
- printf("Nd=%d",h2RMDetector->GetNbinsX());
+ printf("Nt=%d\n",h2ResponseMatrixCombined->GetNbinsX());
+ printf("Nm=%d\n",h2ResponseMatrixCombined->GetNbinsY());
+ printf("Nd=%d\n",h2RMDeltaPt->GetNbinsX());
}
+
for(Int_t t=1; t<=h2ResponseMatrixCombined->GetNbinsX();t++) {
for(Int_t m=1; m<=h2ResponseMatrixCombined->GetNbinsY();m++) {
Double_t valueSum = 0.;
for(Int_t d=1; d<=h2RMDeltaPt->GetNbinsX();d++) {
valueSum += h2RMDeltaPt->GetBinContent(d,m) * h2RMDetector->GetBinContent(t,d);
+ // if(t==10 && m==10) cout << "sum m,d=" << m << "," << d << endl;
}//d-loop
- // cout << "t,m = " << t << "," << m << endl;
+ // if(t==10) cout << "t,m = " << t << "," << m << "\tvalueSum: " << valueSum << endl;
h2ResponseMatrixCombined->SetBinContent(t,m,valueSum);
} //m-loop
}//t-loop
}
+//--------------------------------------------------------------------------------------------------------------------------------------------------
+TH2* AliAnaChargedJetResponseMaker::GetTransposeResponsMatrix(TH2 *h2RM) {
+ //
+ // Transpose response matrix
+ //
+
+ //Initialize transposed response matrix h2RMTranspose
+ TArrayD *arrayX = (TArrayD*)h2RM->GetXaxis()->GetXbins();
+ for(int i=0; i<arrayX->GetSize(); i++) cout << "i: " << arrayX->At(i) << endl;
+ Double_t *xbinsArray = arrayX->GetArray();
+
+ TArrayD *arrayY = (TArrayD*)h2RM->GetYaxis()->GetXbins();
+ for(int i=0; i<arrayY->GetSize(); i++) cout << "i: " << arrayY->At(i) << endl;
+ Double_t *ybinsArray = arrayY->GetArray();
+
+ TH2D *h2RMTranspose = new TH2D("h2RMTranspose","h2RMTranspose",h2RM->GetNbinsY(),ybinsArray,h2RM->GetNbinsX(),xbinsArray);
+
+ //Fill tranposed response matrix
+ for (Int_t ibin = 1; ibin <= h2RMTranspose->GetNbinsX(); ibin++) {
+ for (Int_t jbin = 1; jbin <= h2RMTranspose->GetNbinsY(); jbin++) {
+ h2RMTranspose->SetBinContent(ibin,jbin, h2RM->GetBinContent(jbin,ibin));
+ }
+ }
+
+
+ return h2RMTranspose;
+
+}
+
+//--------------------------------------------------------------------------------------------------------------------------------------------------
+TH2* AliAnaChargedJetResponseMaker::NormalizeResponsMatrixYaxisWithPrior(TH2 *h2RM, TH1 *hPrior) {
+ //
+ // Normalize such that the Y projection is the prior
+ //
+
+ // TH1D *hProjRespY = (TH1D*)h2RM->ProjectionY("hProjRespY");
+ double intPrior = hPrior->Integral();//"width");
+ for (Int_t jbin = 1; jbin <= h2RM->GetNbinsY(); jbin++) {
+ // double corr = hPrior->GetBinContent(jbin)/hProjRespY->GetBinContent(jbin);
+ for (Int_t ibin = 1; ibin <= h2RM->GetNbinsX(); ibin++) {
+ double content = h2RM->GetBinContent(ibin,jbin);
+ // h2RM->SetBinContent(ibin,jbin,content*corr);
+ h2RM->SetBinContent(ibin,jbin,hPrior->GetBinContent(jbin)/hPrior->GetBinWidth(jbin)/intPrior*content);
+ }
+ }
+
+ return h2RM;
+}
+
//--------------------------------------------------------------------------------------------------------------------------------------------------
TH1D* AliAnaChargedJetResponseMaker::MultiplyResponseGenerated(TH1 *hGen, TH2 *hResponse,TH1 *hEfficiency,Bool_t bDrawSlices) {
//
if(!hEfficiency) {
- printf("Efficiency must be given. In case efficiency is 1 use SetFlatEfficiency(1.) before calling function. Aborting!");
- return 0;
+ // printf("Efficiency must be given. In case efficiency is 1 use SetFlatEfficiency(1.) before calling function. Aborting!");
+ // return 0;
+ printf("Setting efficiency to 1 \n");
+ hEfficiency = (TH1D*)hGen->Clone("hEfficiency");
+ hEfficiency->Reset();
+ for(int i=1; i<=hEfficiency->GetNbinsX(); i++) hEfficiency->SetBinContent(i,1.);
}
//For response
Double_t sumYield = 0.;
const Int_t nbinsRec = hRec->GetNbinsX();
Double_t sumError2[nbinsRec+1];
+ for(int irec=0; irec<=hRec->GetNbinsX(); irec++)
+ sumError2[irec]=0.;
Double_t eff = 0.;
for(int igen=1; igen<=hGen->GetNbinsX(); igen++) {
//x-axis: generated
//y-axis: reconstructed
- if(fDebug>0) printf(">>AliAnaChargedJetResponseMaker::MultiplyResponseGenerated(TF1 *fGen, TH2 *hResponse,TH1 *hEfficiency)");
-
- if(!hEfficiency) {
- printf("Efficiency must be given. In case efficiency is 1 use SetFlatEfficiency(1.) before calling function. Aborting!");
- return 0;
- }
+ if(fDebug>0) printf(">>AliAnaChargedJetResponseMaker::MultiplyResponseGenerated(TF1 *fGen, TH2 *hResponse,TH1 *hEfficiency)\n");
TH1D *hRec = hResponse->ProjectionY("hRec");
hRec->Sumw2();
//get pTMC
sumYield = 0.;
double pTMC = hResponse->GetXaxis()->GetBinCenter(igen);
- int binEff = hEfficiency->FindBin(pTMC);
- if(fEffFlat>0.)
- eff = fEffFlat;
- else
+ if(hEfficiency) {
+ int binEff = hEfficiency->FindBin(pTMC);
eff = hEfficiency->GetBinContent(binEff);
- yieldMC = fGen->Eval(pTMC)*eff;
+ }
+ else eff = 1.;
+ // yieldMC = fGen->Eval(pTMC)*eff;
+ yieldMC = fGen->Integral(hResponse->GetXaxis()->GetBinLowEdge(igen),hResponse->GetXaxis()->GetBinUpEdge(igen))*eff;
for(int irec=1; irec<=hResponse->GetNbinsY(); irec++) {
hRec->AddBinContent(irec,yieldMC*hResponse->GetBinContent(igen,irec));
sumYield+=hResponse->GetBinContent(igen,irec);
}
- cout << "igen: " << igen << "\tpTMC: " << pTMC << "\tsumYield: " << sumYield << endl;
+ // cout << "igen: " << igen << "\tpTMC: " << pTMC << "\tsumYield: " << sumYield << endl;
+ // cout << "yieldMC: " << yieldMC << endl;
+
}
return hRec;
}
+//--------------------------------------------------------------------------------------------------------------------------------------------------
+/* static */ Double_t AliAnaChargedJetResponseMaker::GetBetaPerDOFValue(Int_t betaColl, Int_t betaOpt) {
+
+ Float_t betaPerDOFoptions[6] = {0.};
+ //define beta per degree of freedom (DOF=nbins in measured)
+ if(betaColl==0) {
+ betaPerDOFoptions[0] = 9.1e-9;
+ betaPerDOFoptions[1] = 2.25e-8;
+ betaPerDOFoptions[2] = 4.55e-8;
+ betaPerDOFoptions[3] = 9.1e-8;
+ betaPerDOFoptions[4] = 2.25e-7;
+ betaPerDOFoptions[5] = 4.55e-7;
+ }
+ if(betaColl==1) {
+ betaPerDOFoptions[0] = 9.1e-7;
+ betaPerDOFoptions[1] = 2.25e-6;
+ betaPerDOFoptions[2] = 4.55e-6;
+ betaPerDOFoptions[3] = 9.1e-6;
+ betaPerDOFoptions[4] = 2.25e-5;
+ betaPerDOFoptions[5] = 4.55e-5;
+ }
+ if(betaColl==2) {
+ betaPerDOFoptions[0] = 9.1e-5;
+ betaPerDOFoptions[1] = 2.25e-4;
+ betaPerDOFoptions[2] = 4.55e-4;
+ betaPerDOFoptions[3] = 9.1e-4;
+ betaPerDOFoptions[4] = 2.25e-3;
+ betaPerDOFoptions[5] = 4.55e-3;
+ }
+ if(betaColl==3) {
+ betaPerDOFoptions[0] = 9.1e-3;
+ betaPerDOFoptions[1] = 2.25e-2;
+ betaPerDOFoptions[2] = 4.55e-2;
+ betaPerDOFoptions[3] = 9.1e-2;
+ betaPerDOFoptions[4] = 2.25e-1;
+ betaPerDOFoptions[5] = 4.55e-1;
+ }
+ if(betaColl==4) {
+ betaPerDOFoptions[0] = 9.1e-1;
+ betaPerDOFoptions[1] = 2.25;
+ betaPerDOFoptions[2] = 4.55;
+ betaPerDOFoptions[3] = 9.1;
+ betaPerDOFoptions[4] = 22.5;
+ betaPerDOFoptions[5] = 45.5;
+ }
+
+ return betaPerDOFoptions[betaOpt];
+
+}
+
//--------------------------------------------------------------------------------------------------------------------------------------------------
Double_t AliAnaChargedJetResponseMaker::InterpolateFast(TGraph *gr, Double_t x) {
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff