#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;
}
int bin[1] = {0};
bin[0] = 0;
for(int i = hPtMeasured->FindBin(fPtMin); i<hPtMeasured->FindBin(fPtMax); i++) {
- double pT[1];
- pT[0]= hPtMeasured->GetBinCenter(i);
fPtMeasured->SetBinContent(bin,(Double_t)hPtMeasured->GetBinContent(i));
fPtMeasured->SetBinError(bin,(Double_t)hPtMeasured->GetBinError(i));
bin[0]++;
fEffFlat = eff;
return;
- /*
+ /*
Int_t nbins[fDimensions];
Double_t xmin[fDimensions];
Double_t xmax[fDimensions];
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);
printf("binWidthUnfLowPt: %f\n",binWidthUnfLowPt);
- int tmp = round((fPtMaxUnfolded-fPtMinUnfolded)/binWidthUnf); //nbins which fit between 0 and fPtMaxUnfolded
- tmp = tmp - fSkipBinsUnfolded;
- fPtMinUnfolded = fPtMaxUnfolded-(double)(tmp)*binWidthUnf; //set ptmin unfolded
- fPtMaxUnfolded = fPtMinUnfolded+(double)(tmp)*binWidthUnf; //set ptmax unfolded
- nbins[fDimGen] = (int)((fPtMaxUnfolded-fPtMinUnfolded)/binWidthUnf); //adjust nbins to bin width
+ printf("fPtMinUnfolded: %f fPtMaxUnfolded: %f\n",fPtMinUnfolded,fPtMaxUnfolded);
+
if(fbVariableBinning) {
- tmp = (int)(fPtMaxUnfVarBinning/binWidthUnfLowPt);
+ // cout << "fPtMaxUnfVarBinning: " << fPtMaxUnfVarBinning << " \tfPtMinUnfolded: " << fPtMinUnfolded << " binWidthUnfLowPt: " << binWidthUnfLowPt << endl;
+ Int_t tmp = (int)((fPtMaxUnfVarBinning-fPtMinUnfolded)/binWidthUnfLowPt);
tmp = tmp - fSkipBinsUnfolded;
fPtMinUnfolded = fPtMaxUnfVarBinning-(double)(tmp)*binWidthUnfLowPt;
+ //cout << "tmp = " << tmp << " fSkipBinsUnfolded = " << fSkipBinsUnfolded << " fPtMinUnfolded = " << fPtMinUnfolded << endl;
//Redefine also number of bins on generated axis in case of variable binning
nbins[fDimGen] = (int)((fPtMaxUnfVarBinning-fPtMinUnfolded)/binWidthUnfLowPt)+(int)((fPtMaxUnfolded-fPtMaxUnfVarBinning)/binWidthUnf);
}
+ else {
+ int tmp = round((fPtMaxUnfolded-fPtMinUnfolded)/binWidthUnf); //nbins which fit between 0 and fPtMaxUnfolded
+ tmp = tmp - fSkipBinsUnfolded;
+ fPtMinUnfolded = fPtMaxUnfolded-(double)(tmp)*binWidthUnf; //set ptmin unfolded
+ fPtMaxUnfolded = fPtMinUnfolded+(double)(tmp)*binWidthUnf; //set ptmax unfolded
+ nbins[fDimGen] = (int)((fPtMaxUnfolded-fPtMinUnfolded)/binWidthUnf); //adjust nbins to bin width
+ }
+
+ printf(" nbins[fDimGen] = %d nbins[fDimRec] = %d\n",nbins[fDimGen],nbins[fDimRec]);
Double_t binWidth[2];
binWidth[fDimRec] = (double)((fPtMax-fPtMin)/nbins[fDimRec]);
Int_t nbinsFine[fDimensions*2];
Double_t xminFine[fDimensions*2];
Double_t xmaxFine[fDimensions*2];
- Double_t pTarrayFine[fDimensions*2];
+ // Double_t pTarrayFine[fDimensions*2];
nbinsFine[fDimRec] = fResponseMatrix->GetAxis(fDimRec)->GetNbins()*fFineFrac;
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.;
- pTarrayFine[fDimRec] = 0.;
- pTarrayFine[fDimGen] = 0.;
+ xmaxFine[fDimRec] = fResponseMatrix->GetAxis(fDimGen)->GetXmax();//+40.;
+ xmaxFine[fDimGen] = fResponseMatrix->GetAxis(fDimGen)->GetXmax();//+40.;
+ // pTarrayFine[fDimRec] = 0.;
+ // pTarrayFine[fDimGen] = 0.;
Double_t binWidth[2];
binWidth[fDimRec] = fResponseMatrix->GetAxis(fDimRec)->GetBinWidth(1);
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
TAxis *recAxis = fResponseMatrixFine->GetAxis(fDimRec);
Int_t nbinsFine[fDimensions*2];
- Double_t xminFine[fDimensions*2];
- Double_t xmaxFine[fDimensions*2];
+ //Double_t xminFine[fDimensions*2];
+ //Double_t xmaxFine[fDimensions*2];
Double_t pTarrayFine[fDimensions*2];
nbinsFine[fDimGen] = genAxis->GetNbins();
nbinsFine[fDimRec] = recAxis->GetNbins();
- xminFine[fDimGen] = genAxis->GetXmin();
- xminFine[fDimRec] = recAxis->GetXmin();
- xmaxFine[fDimGen] = genAxis->GetXmax();
- xmaxFine[fDimRec] = recAxis->GetXmax();
+ // xminFine[fDimGen] = genAxis->GetXmin();
+ // xminFine[fDimRec] = recAxis->GetXmin();
+ // xmaxFine[fDimGen] = genAxis->GetXmax();
+ // xmaxFine[fDimRec] = recAxis->GetXmax();
pTarrayFine[fDimGen] = 0.;
pTarrayFine[fDimRec] = 0.;
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());
// printf("weight: %f \n", weight );
} else {
weight = 1./((double)fFineFrac);
+ if(fbVariableBinning && pTgen<fPtMaxUnfVarBinning) weight=1./(0.5*(double)fFineFrac);
}
for(int ix=ixMin; ix<=ixMax; ix++) { //loop reconstructed axis
} // iy (dimGen) loop
//Fill Efficiency corresponding to merged response matrix
- // FillEfficiency(errorArray,fFineFrac);
-
for(int iy=1; iy<=fResponseMatrix->GetAxis(fDimGen)->GetNbins(); iy++) { //gen
sumyield = 0.;
ipt[fDimGen] = iy;
}
//--------------------------------------------------------------------------------------------------------------------------------------------------
-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();
- //Forst normalize columns of input
+ 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(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.;
Int_t ixNew = 0;
Int_t iyNew = 0;
- Double_t xval, yval;
+ Double_t xvalLo, xvalUp, yvalLo, yvalUp;
Double_t xmin = hRM2->GetXaxis()->GetXmin();
Double_t ymin = hRM2->GetYaxis()->GetXmin();
Double_t xmax = hRM2->GetXaxis()->GetXmax();
Double_t ymax = hRM2->GetYaxis()->GetXmax();
for(int ix=1; ix<=hRMFine->GetXaxis()->GetNbins(); ix++) {
- xval = hRMFine->GetXaxis()->GetBinCenter(ix);
- if(xval<xmin || xval>xmax) continue;
- ixNew = hRM2->GetXaxis()->FindBin(xval);
+ xvalLo = hRMFine->GetXaxis()->GetBinLowEdge(ix);
+ xvalUp = hRMFine->GetXaxis()->GetBinUpEdge(ix);
+ if(xvalLo<xmin || xvalUp>xmax) continue;
+ ixNew = hRM2->GetXaxis()->FindBin(hRMFine->GetXaxis()->GetBinCenter(ix));
for(int iy=1; iy<=hRMFine->GetYaxis()->GetNbins(); iy++) {
- yval = hRMFine->GetYaxis()->GetBinCenter(iy);
- if(yval<ymin || yval>ymax) continue;
+ yvalLo = hRMFine->GetYaxis()->GetBinLowEdge(iy);
+ yvalUp = hRMFine->GetYaxis()->GetBinUpEdge(iy);
+ if(yvalLo<ymin || yvalUp>ymax) continue;
content = hRMFine->GetBinContent(ix,iy);
- iyNew = hRM2->GetYaxis()->FindBin(yval);
+ iyNew = hRM2->GetYaxis()->FindBin(hRMFine->GetYaxis()->GetBinCenter(iy));
oldcontent = hRM2->GetBinContent(ixNew,iyNew);
- // if(fDebug) cout << "ixNew: " << ixNew << " " << xval << " iyNew: " << iyNew << " " << yval << " content: " << content << " oldContent: " << oldcontent << " newContent: " << oldcontent+content << endl;
+ //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) {
+
+ //
+ // Make combined response matrix (background flucutuations + detector effects)
+ //
+ // hRMDeltaPt is the response matrix for background fluctuations from \delta(p_t) measurement
+ // hRMDetector is the response matrix for detector effects: needs to be a squared matrix with on each axis the same bins as on the generated axis of the bkg fluctuations response matrix
+ //
+ // 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"); //h2RMDeltaPt is the bkg fluctuations RM which has the dimensions we want for the combined response matrix
+ h2ResponseMatrixCombined->SetTitle("h2ResponseMatrixCombined");
+ h2ResponseMatrixCombined->SetName("h2ResponseMatrixCombined");
+
+ // M = RM_deltaPt * RM_DetEffects * T (M=measured T=truth)
+ // Dimensions:
+ // mx1 = mxd * dxt * tx1
+ // m = measured bins
+ // t = truth bins
+ // d = rec for RM_DetEffects and gen for RM_deltaPt
+ // RM_comb = RM_deltaPt * RM_DetEffects (dimensions mxt)
+ // RM_comb(m,t) = Sum_d RM_deltaPt(m,d)*RM_DetEffects(d,t)
+
+ if(fDebug) {
+ 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
+ // if(t==10) cout << "t,m = " << t << "," << m << "\tvalueSum: " << valueSum << endl;
+ h2ResponseMatrixCombined->SetBinContent(t,m,valueSum);
+ } //m-loop
+ }//t-loop
+
+ return h2ResponseMatrixCombined;
+
+}
+
+//--------------------------------------------------------------------------------------------------------------------------------------------------
+TH2* AliAnaChargedJetResponseMaker::GetTransposeResponsMatrix(TH2 *h2RM) {
+ //
+ // Transpose response matrix
+ //
+
+ Printf("AliAnaChargedJetResponseMaker::GetTransposeResponsMatrix");
+
+ //Initialize transposed response matrix h2RMTranspose
+ // TArrayD *arrayX = (TArrayD*)h2RM->GetXaxis()->GetXbins();
+ // for(int i=0; i<arrayX->GetSize(); i++) Printf("i: %d %f", i,arrayX->At(i));
+ // Double_t *xbinsArray = arrayX->GetArray();
+
+ // TArrayD *arrayY = (TArrayD*)h2RM->GetYaxis()->GetXbins();
+ // for(int i=0; i<arrayY->GetSize(); i++) Printf("i: %d %f", i,arrayY->At(i));
+ // Double_t *ybinsArray = arrayY->GetArray();
+
+
+ Double_t *ybinsArray = new Double_t[h2RM->GetNbinsY()+1];
+ for(int i=1; i<=h2RM->GetNbinsY(); i++) {
+ ybinsArray[i-1] = h2RM->GetYaxis()->GetBinLowEdge(i);
+ Printf("i=%d %f %f",i,ybinsArray[i-1],h2RM->GetYaxis()->GetBinLowEdge(i));
+ }
+ ybinsArray[h2RM->GetNbinsY()] = h2RM->GetYaxis()->GetBinUpEdge(h2RM->GetNbinsY());
+
+ Double_t *xbinsArray = new Double_t[h2RM->GetNbinsX()+1];
+ for(int i=1; i<=h2RM->GetNbinsX(); i++)
+ xbinsArray[i-1] = h2RM->GetXaxis()->GetBinLowEdge(i);
+ xbinsArray[h2RM->GetNbinsX()] = h2RM->GetXaxis()->GetBinUpEdge(h2RM->GetNbinsX());
+
+
+ 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++) {
else
eff = hEfficiency->GetBinContent(igen);
yieldMC = hGen->GetBinContent(igen)*eff;
- // yieldMCerror = hGen->GetBinError(igen);
+ yieldMCerror = hGen->GetBinError(igen)*eff;
if(bDrawSlices) {
hRecGenBin = hResponse->ProjectionY(Form("hRecGenBin%d",igen));
for(int irec=1; irec<=hRec->GetNbinsX(); irec++) {
hRec->AddBinContent(irec,yieldMC*hResponse->GetBinContent(igen,irec));
sumYield+=hResponse->GetBinContent(igen,irec);
- Double_t A = yieldMC*hResponse->GetBinError(igen,irec);
+ // Double_t A = yieldMC*hResponse->GetBinError(igen,irec);
Double_t B = hResponse->GetBinContent(igen,irec)*yieldMCerror;
- Double_t tmp2 = A*A + B*B;
- sumError2[irec-1] += tmp2 ;
+ // Double_t tmp2 = A*A + B*B;
+ //sumError2[irec-1] += tmp2 ;
+ sumError2[irec-1] += B*B;
if(bDrawSlices)
hRecGenBin->SetBinContent(irec,yieldMC*hResponse->GetBinContent(igen,irec));
//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