PWGJE/AliAnaChargedJetResponseMaker.cxx

   1 #include "AliAnaChargedJetResponseMaker.h"
   2 #include "TGraph.h"
   3 #include "TGraphErrors.h"
   4 #include "TMath.h"
   5 #include "Riostream.h"
   6 #include "TH1.h"
   7 #include "TRandom.h"
   8 #include "TFile.h"
   9 #include "TCanvas.h"
  10 #include "TF1.h"
  11 #include "THnSparse.h"
  12 #include "TArrayD.h"
  13
  14 #define round(x) ((x)>=0?(long)((x)+0.5):(long)((x)-0.5))
  15
  16 using std::cout;
  17 using std::endl;
  18
  19 ClassImp(AliAnaChargedJetResponseMaker)
  20
  21 AliAnaChargedJetResponseMaker::AliAnaChargedJetResponseMaker():
  22   fDebug(kFALSE),
  23   fResolutionType(kParam),
  24   fDeltaPt(0x0),
  25   fhDeltaPt(0x0),
  26   fh1MeasuredTruncated(0x0),
  27   fh2DetectorResponse(0x0),
  28   fh2DetectorResponseRebin(0x0),
  29   fhEfficiencyDet(0x0),
  30   fh2ResponseMatrixCombinedFineFull(0x0),
  31   fh2ResponseMatrixCombinedFull(0x0),
  32   fh2ResponseMatrixCombined(0x0),
  33   fhEfficiencyCombined(0x0),
  34   fDimensions(1),
  35   fDimRec(0),
  36   fDimGen(1),
  37   fPtMin(-999),
  38   fPtMax(-999),
  39   fNbins(0),
  40   fBinArrayPtRec(0x0),
  41   fPtMeasured(0x0),
  42   fEffFlat(0),
  43   fEfficiency(0x0),
  44   fEfficiencyFine(0x0),
  45   fResponseMatrix(0x0),
  46   fResponseMatrixFine(0x0),
  47   fPtMinUnfolded(0.),
  48   fPtMaxUnfolded(0.),
  49   fPtMaxUnfoldedHigh(-1.),
  50   fBinWidthFactorUnfolded(2),
  51   fSkipBinsUnfolded(0),
  52   fExtraBinsUnfolded(5),
  53   fbVariableBinning(kFALSE),
  54   fPtMaxUnfVarBinning(0),
  55   f1MergeFunction(0x0),
  56   fFineFrac(10),
  57   fbCalcErrors(kFALSE)
  58 {;}
  59
  60
  61 //--------------------------------------------------------------------------------------------------------------------------------------------------
  62 AliAnaChargedJetResponseMaker::AliAnaChargedJetResponseMaker(const AliAnaChargedJetResponseMaker& obj):
  63   fDebug(obj.fDebug),
  64   fResolutionType(obj.fResolutionType),
  65   fDeltaPt(obj.fDeltaPt),
  66   fhDeltaPt(obj.fhDeltaPt),
  67   fh1MeasuredTruncated(obj.fh1MeasuredTruncated),
  68   fh2DetectorResponse(obj.fh2DetectorResponse),
  69   fh2DetectorResponseRebin(obj.fh2DetectorResponseRebin),
  70   fhEfficiencyDet(obj.fhEfficiencyDet),
  71   fh2ResponseMatrixCombinedFineFull(obj.fh2ResponseMatrixCombinedFineFull),
  72   fh2ResponseMatrixCombinedFull(obj.fh2ResponseMatrixCombinedFull),
  73   fh2ResponseMatrixCombined(obj.fh2ResponseMatrixCombined),
  74   fhEfficiencyCombined(obj.fhEfficiencyCombined),
  75   fDimensions(obj.fDimensions),
  76   fDimRec(obj.fDimRec),
  77   fDimGen(obj.fDimGen),
  78   fPtMin(obj.fPtMin),
  79   fPtMax(obj.fPtMax),
  80   fNbins(obj.fNbins),
  81   fBinArrayPtRec(obj.fBinArrayPtRec),
  82   fPtMeasured(obj.fPtMeasured),
  83   fEffFlat(obj.fEffFlat),
  84   fEfficiency(obj.fEfficiency),
  85   fEfficiencyFine(obj.fEfficiencyFine),
  86   fResponseMatrix(obj.fResponseMatrix),
  87   fResponseMatrixFine(obj.fResponseMatrixFine),
  88   fPtMinUnfolded(obj.fPtMinUnfolded),
  89   fPtMaxUnfolded(obj.fPtMaxUnfolded),
  90   fPtMaxUnfoldedHigh(obj.fPtMaxUnfoldedHigh),
  91   fBinWidthFactorUnfolded(obj.fBinWidthFactorUnfolded),
  92   fSkipBinsUnfolded(obj.fSkipBinsUnfolded),
  93   fExtraBinsUnfolded(obj.fExtraBinsUnfolded),
  94   fbVariableBinning(obj.fbVariableBinning),
  95   fPtMaxUnfVarBinning(obj.fPtMaxUnfVarBinning),
  96   f1MergeFunction(obj.f1MergeFunction),
  97   fFineFrac(obj.fFineFrac),
  98   fbCalcErrors(obj.fbCalcErrors)
  99 {;}
 100
 101 //--------------------------------------------------------------------------------------------------------------------------------------------------
 102 AliAnaChargedJetResponseMaker& AliAnaChargedJetResponseMaker::operator=(const AliAnaChargedJetResponseMaker& other)
 103 {
 104 // Assignment
 105   if(&other == this) return *this;
 106   AliAnaChargedJetResponseMaker::operator=(other);
 107   fDebug                    = other.fDebug;
 108   fResolutionType           = other.fResolutionType;
 109   fDeltaPt                  = other.fDeltaPt;
 110   fhDeltaPt                 = other.fhDeltaPt;
 111   fh1MeasuredTruncated      = other.fh1MeasuredTruncated;
 112   fh2DetectorResponse       = other.fh2DetectorResponse;
 113   fh2DetectorResponseRebin  = other.fh2DetectorResponseRebin;
 114   fhEfficiencyDet           = other.fhEfficiencyDet;
 115   fh2ResponseMatrixCombinedFineFull = other.fh2ResponseMatrixCombinedFineFull;
 116   fh2ResponseMatrixCombinedFull = other.fh2ResponseMatrixCombinedFull;
 117   fh2ResponseMatrixCombined = other.fh2ResponseMatrixCombined;
 118   fhEfficiencyCombined      = other.fhEfficiencyCombined;
 119   fDimensions               = other.fDimensions;
 120   fDimRec                   = other.fDimRec;
 121   fDimGen                   = other.fDimGen;
 122   fPtMin                    = other.fPtMin;
 123   fPtMax                    = other.fPtMax;
 124   fNbins                    = other.fNbins;
 125   fBinArrayPtRec            = other.fBinArrayPtRec;
 126   fPtMeasured               = other.fPtMeasured;
 127   fEffFlat                  = other.fEffFlat;
 128   fEfficiency               = other.fEfficiency;
 129   fEfficiencyFine           = other.fEfficiencyFine;
 130   fResponseMatrix           = other.fResponseMatrix;
 131   fResponseMatrixFine       = other.fResponseMatrixFine;
 132   fPtMinUnfolded            = other.fPtMinUnfolded;
 133   fPtMaxUnfolded            = other.fPtMaxUnfolded;
 134   fPtMaxUnfoldedHigh        = other.fPtMaxUnfoldedHigh;
 135   fBinWidthFactorUnfolded   = other.fBinWidthFactorUnfolded;
 136   fSkipBinsUnfolded         = other.fSkipBinsUnfolded;
 137   fExtraBinsUnfolded        = other.fExtraBinsUnfolded;
 138   fbVariableBinning         = other.fbVariableBinning;
 139   fPtMaxUnfVarBinning       = other.fPtMaxUnfVarBinning;
 140   f1MergeFunction           = other.f1MergeFunction;
 141   fFineFrac                 = other.fFineFrac;
 142   fbCalcErrors              = other.fbCalcErrors;
 143
 144   return *this;
 145 }
 146
 147 //--------------------------------------------------------------------------------------------------------------------------------------------------
 148 void AliAnaChargedJetResponseMaker::SetMeasuredSpectrum(TH1D *hPtMeasured)
 149 {
 150   //
 151   // Set measured spectrum in THnSparse format
 152   // This defines the minimum and maximum pT on the reconstructed axis of the response matrix
 153   //
 154   if(fDebug) printf(">>AliAnaChargedJetResponseMaker::SetMeasuredSpectrum \n");
 155
 156   fNbins = hPtMeasured->GetXaxis()->GetNbins();
 157   fPtMin = hPtMeasured->GetXaxis()->GetXmin();
 158   fPtMax = hPtMeasured->GetXaxis()->GetXmax();
 159
 160   if(fDebug) printf("fNbins: %d  fPtMin: %f  fPtMax: %f \n",fNbins,fPtMin,fPtMax);
 161
 162   if(fBinArrayPtRec) delete fBinArrayPtRec;
 163   fBinArrayPtRec = new Double_t[fNbins+1];
 164   for(int j = 0; j<fNbins; j++) {
 165     fBinArrayPtRec[j] = hPtMeasured->GetXaxis()->GetBinLowEdge(j+1);
 166   }
 167   fBinArrayPtRec[fNbins] = hPtMeasured->GetXaxis()->GetBinUpEdge(fNbins);
 168
 169
 170   Int_t nbins[fDimensions];
 171   Double_t xmin[fDimensions];
 172   Double_t xmax[fDimensions];
 173   for(int dim = 0; dim<fDimensions; dim++) {
 174     nbins[dim] = fNbins;
 175     xmin[dim]  = fPtMin;
 176     xmax[dim]  = fPtMax;
 177   }
 178
 179   if(fPtMeasured) delete fPtMeasured;
 180   fPtMeasured = new THnSparseD("fPtMeasured","Measured pT spectrum; p_{T}^{rec}",fDimensions,nbins,xmin,xmax);
 181   fPtMeasured->Sumw2();
 182   fPtMeasured->GetAxis(0)->SetTitle("p_{T}^{rec}");
 183   fPtMeasured->SetBinEdges(0,fBinArrayPtRec);
 184
 185   //Fill
 186   if(fDebug) printf("fill measured THnSparse\n");
 187   if(fNbins!=hPtMeasured->GetNbinsX())
 188     printf("WARNING: nbins not correct \t %d vs %d !!!\n",fNbins,hPtMeasured->GetNbinsX());
 189
 190   int bin[1] = {0};
 191   bin[0] = 0;
 192   for(int i = hPtMeasured->FindBin(fPtMin); i<hPtMeasured->FindBin(fPtMax); i++) {
 193     double pT[1];
 194     pT[0]= hPtMeasured->GetBinCenter(i);
 195     fPtMeasured->SetBinContent(bin,(Double_t)hPtMeasured->GetBinContent(i));
 196     fPtMeasured->SetBinError(bin,(Double_t)hPtMeasured->GetBinError(i));
 197     bin[0]++;
 198   }
 199
 200   if(fDebug) printf("fPtMeasured->GetNbins(): %d \n",(int)fPtMeasured->GetNbins());
 201
 202 }
 203
 204 //--------------------------------------------------------------------------------------------------------------------------------------------------
 205 void AliAnaChargedJetResponseMaker::SetFlatEfficiency(Double_t eff) {
 206
 207   //
 208   // Set flat efficiency to value of eff
 209   //
 210
 211   fEffFlat = eff;
 212   return;
 213
 214   /*
 215   Int_t nbins[fDimensions];
 216   Double_t xmin[fDimensions];
 217   Double_t xmax[fDimensions];
 218   for(int dim = 0; dim<fDimensions; dim++) {
 219     nbins[dim] = fNbins;
 220     xmin[dim] = fPtMin;
 221     xmax[dim] = fPtMax;
 222   }
 223
 224   if(fEfficiency) delete fEfficiency;
 225   fEfficiency = new THnSparseD("fEfficiency","Efficiency - no resolution effects; p_{T}^{gen}",fDimensions,nbins,xmin,xmax);
 226   fEfficiency->Sumw2();
 227   fEfficiency->GetAxis(0)->SetTitle("p_{T}^{gen}");
 228   fEfficiency->SetBinEdges(0,fBinArrayPtRec);
 229
 230   for(int i=0; i<fNbins; i++) {
 231     int bin[1];
 232     bin[0] = i;
 233     fEfficiency->SetBinContent(bin,fEffFlat);
 234     fEfficiency->SetBinError(bin,0);
 235   }
 236   */
 237 }
 238
 239 //--------------------------------------------------------------------------------------------------------------------------------------------------
 240 void AliAnaChargedJetResponseMaker::SetEfficiency(TGraphErrors *grEff)
 241 {
 242   //
 243   // Fill fEfficiency (THnSparse) with values from grEff
 244   //
 245
 246   Int_t nbins[fDimensions];
 247   Double_t xmin[fDimensions];
 248   Double_t xmax[fDimensions];
 249   for(int dim = 0; dim<fDimensions; dim++) {
 250     nbins[dim] = fNbins;
 251     xmin[dim] = fPtMin;
 252     xmax[dim] = fPtMax;
 253   }
 254
 255   if(fEfficiency) delete fEfficiency;
 256   fEfficiency = new THnSparseD("fEfficiency","Efficiency - no resolution effects; p_{T}^{gen}",fDimensions,nbins,xmin,xmax);
 257   fEfficiency->Sumw2();
 258   fEfficiency->GetAxis(0)->SetTitle("p_{T}^{gen}");
 259   fEfficiency->SetBinEdges(0,fBinArrayPtRec);
 260
 261   double pT[1];
 262   double yield = 0.;
 263   double error = 0.;
 264   double dummy = 0.;
 265   int nbinsgr = grEff->GetN();
 266
 267   for(int i=0; i<nbinsgr; i++) {
 268     grEff->GetPoint(i,dummy,yield);
 269     pT[0] = dummy;
 270     error = grEff->GetErrorY(i);
 271
 272     fEfficiency->Fill(pT,yield);
 273     fEfficiency->SetBinError(i,error);
 274
 275   }
 276
 277 }
 278
 279 //--------------------------------------------------------------------------------------------------------------------------------------------------
 280 void AliAnaChargedJetResponseMaker::MakeResponseMatrixJetsFineMerged(Int_t skipBins, Int_t binWidthFactor, Int_t extraBins, Bool_t bVariableBinning, Double_t ptmax)
 281 {
 282   //
 283   // Make jet response matrix
 284   //
 285
 286   if(fDebug) printf(">>AliAnaChargedJetResponseMaker::MakeResponseMatrixJetsFineMerged\n");
 287
 288   if(!fPtMeasured) {
 289     if(fDebug) printf("fPtMeasured does not exist. Aborting!!\n");
 290     return;
 291   }
 292   if(fResponseMatrix)     { delete fResponseMatrix; }
 293   if(fResponseMatrixFine) { delete fResponseMatrixFine; }
 294
 295   SetSkipBinsUnfolded(skipBins);
 296   SetBinWidthFactorUnfolded(binWidthFactor);
 297   SetExtraBinsUnfolded(extraBins);
 298   SetVariableBinning(bVariableBinning,ptmax);
 299
 300   InitializeResponseMatrix();
 301   InitializeEfficiency();
 302
 303   InitializeResponseMatrixFine();
 304   InitializeEfficiencyFine();
 305
 306   FillResponseMatrixFineAndMerge();
 307
 308 }
 309
 310 //--------------------------------------------------------------------------------------------------------------------------------------------------
 311 void AliAnaChargedJetResponseMaker::InitializeResponseMatrix() {
 312   //
 313   //Define bin edges of RM to be used for unfolding
 314   //
 315
 316   Int_t nbins[fDimensions*2];
 317   nbins[fDimRec] = fNbins;
 318   nbins[fDimGen] = fNbins;
 319
 320   double binWidthMeas = (double)((fPtMax-fPtMin)/fNbins);
 321   double binWidthUnf  = (double)fBinWidthFactorUnfolded*binWidthMeas;
 322   double binWidthUnfLowPt = -1.;
 323   if(fbVariableBinning)
 324     binWidthUnfLowPt = binWidthUnf*0.5;
 325
 326   fPtMaxUnfolded = fPtMax+(double)(fExtraBinsUnfolded)*binWidthUnf;
 327   nbins[fDimGen]+=fExtraBinsUnfolded;
 328
 329
 330   printf("fPtMinMeas: %f  fPtMaxMeas: %f\n",fPtMin,fPtMax);
 331   printf("binWidthMeas: %f  binWidthUnf: %f   fBinWidthFactorUnfolded: %d\n",binWidthMeas,binWidthUnf,fBinWidthFactorUnfolded);
 332   printf("binWidthUnfLowPt: %f\n",binWidthUnfLowPt);
 333
 334   printf("fPtMinUnfolded: %f  fPtMaxUnfolded: %f\n",fPtMinUnfolded,fPtMaxUnfolded);
 335
 336
 337   if(fbVariableBinning) {
 338     // cout << "fPtMaxUnfVarBinning: " << fPtMaxUnfVarBinning << " \tfPtMinUnfolded: " << fPtMinUnfolded << "  binWidthUnfLowPt: " << binWidthUnfLowPt << endl;
 339     Int_t tmp = (int)((fPtMaxUnfVarBinning-fPtMinUnfolded)/binWidthUnfLowPt);
 340     tmp = tmp - fSkipBinsUnfolded;
 341     fPtMinUnfolded = fPtMaxUnfVarBinning-(double)(tmp)*binWidthUnfLowPt;
 342     //cout << "tmp = " << tmp << "  fSkipBinsUnfolded = " << fSkipBinsUnfolded << " fPtMinUnfolded = " << fPtMinUnfolded << endl;
 343     //Redefine also number of bins on generated axis in case of variable binning
 344     nbins[fDimGen] = (int)((fPtMaxUnfVarBinning-fPtMinUnfolded)/binWidthUnfLowPt)+(int)((fPtMaxUnfolded-fPtMaxUnfVarBinning)/binWidthUnf);
 345   }
 346   else {
 347     int tmp = round((fPtMaxUnfolded-fPtMinUnfolded)/binWidthUnf); //nbins which fit between 0 and fPtMaxUnfolded
 348     tmp = tmp - fSkipBinsUnfolded;
 349     fPtMinUnfolded = fPtMaxUnfolded-(double)(tmp)*binWidthUnf; //set ptmin unfolded
 350     fPtMaxUnfolded = fPtMinUnfolded+(double)(tmp)*binWidthUnf; //set ptmax unfolded
 351     nbins[fDimGen] = (int)((fPtMaxUnfolded-fPtMinUnfolded)/binWidthUnf); //adjust nbins to bin width
 352   }
 353
 354   printf("   nbins[fDimGen] = %d   nbins[fDimRec] = %d\n",nbins[fDimGen],nbins[fDimRec]);
 355
 356   Double_t binWidth[2];
 357   binWidth[fDimRec] = (double)((fPtMax-fPtMin)/nbins[fDimRec]);
 358   binWidth[fDimGen] = (double)((fPtMaxUnfolded-fPtMinUnfolded)/nbins[fDimGen]);
 359
 360   Double_t xmin[fDimensions*2];
 361   Double_t xmax[fDimensions*2];
 362   xmin[fDimRec] = fPtMin;
 363   xmax[fDimRec] = fPtMax;
 364   xmin[fDimGen] = fPtMinUnfolded;
 365   xmax[fDimGen] = fPtMaxUnfolded;
 366
 367   printf("xmin[fDimRec]: %f  xmin[fDimGen]: %f \n",xmin[fDimRec],xmin[fDimGen]);
 368   printf("xmax[fDimRec]: %f  xmax[fDimGen]: %f \n",xmax[fDimRec],xmax[fDimGen]);
 369   printf("nbins[fDimRec]: %d  nbins[fDimGen]: %d \n",nbins[fDimRec],nbins[fDimGen]);
 370   if(!fbVariableBinning) printf("binWidth[fDimRec]: %f  binWidth[fDimGen]: %f \n",binWidth[fDimRec],binWidth[fDimGen]);
 371
 372   Double_t binArrayPt0[nbins[fDimRec]+1];
 373   Double_t binArrayPt1[nbins[fDimGen]+1];
 374   Double_t xmaxGen = TMath::Max(xmax[fDimGen],fPtMaxUnfoldedHigh);
 375
 376   //Define bin limits reconstructed/measured axis
 377   for(int i=0; i<nbins[fDimRec]; i++) {
 378     binArrayPt0[i] = xmin[fDimRec]+(double)i*binWidth[fDimRec];
 379   }
 380   binArrayPt0[nbins[fDimRec]]= xmax[fDimRec];
 381
 382   //Define bin limits generated/unfolded axis
 383   binArrayPt1[0] = fPtMinUnfolded;
 384   for(int i=1; i<nbins[fDimGen]; i++) {
 385     if(fbVariableBinning) {
 386       double test = xmin[fDimGen]+(double)i*binWidthUnfLowPt;
 387       if(test<=fPtMaxUnfVarBinning) binArrayPt1[i] = test;
 388       else binArrayPt1[i] = binArrayPt1[i-1]+binWidthUnf;
 389     }
 390     else
 391       binArrayPt1[i] = xmin[fDimGen]+(double)i*binWidth[fDimGen];
 392     //printf("RM. i = %d \t binArrayPt[i] = %f \n",i,binArrayPt1[i]);
 393   }
 394   binArrayPt1[nbins[fDimGen]]= xmaxGen;
 395
 396
 397   // Response matrix : dimensions must be 2N = 2 x (number of variables)
 398   // dimensions 0 ->  N-1 must be filled with reconstructed values
 399   // dimensions N -> 2N-1 must be filled with generated values
 400   if(fResponseMatrix) delete fResponseMatrix;
 401   fResponseMatrix = new THnSparseD("fResponseMatrix","Response Matrix pTMC vs pTrec",fDimensions*2,nbins,xmin,xmax);
 402   fResponseMatrix->Sumw2();
 403   fResponseMatrix->GetAxis(fDimRec)->SetTitle("p_{T}^{rec}");
 404   fResponseMatrix->GetAxis(fDimGen)->SetTitle("p_{T}^{gen}");
 405
 406   fResponseMatrix->SetBinEdges(fDimRec,binArrayPt0);
 407   fResponseMatrix->SetBinEdges(fDimGen,binArrayPt1);
 408
 409   Int_t bin[2] = {1,1};
 410   for(int i=1; i<fResponseMatrix->GetAxis(0)->GetNbins(); i++) {
 411     bin[0]=i;
 412     for(int j=1; j<fResponseMatrix->GetAxis(1)->GetNbins(); j++) {
 413     bin[1]=j;
 414       fResponseMatrix->SetBinContent(bin,0.);
 415     }
 416   }
 417
 418
 419
 420 }
 421
 422 //--------------------------------------------------------------------------------------------------------------------------------------------------
 423 void AliAnaChargedJetResponseMaker::InitializeEfficiency() {
 424   //
 425   // Define dimensions of efficiency THnSparse
 426   //
 427
 428   if(!fResponseMatrix) {
 429     printf("AliAnaChargedJetResponseMaker::InitializeEfficiency()\n");
 430     printf("Can not define dimensions efficiency without fResponseMatrix dimensions. Aborting \n");
 431     return;
 432   }
 433
 434   TAxis *genAxis = fResponseMatrix->GetAxis(fDimGen);
 435
 436   Int_t nbinsEff[fDimensions];
 437   Double_t xminEff[fDimensions];
 438   Double_t xmaxEff[fDimensions];
 439
 440   for(int dim = 0; dim<fDimensions; dim++) {
 441     nbinsEff[dim] = genAxis->GetNbins();
 442     xminEff[dim]  = genAxis->GetXmin();
 443     xmaxEff[dim]  = genAxis->GetXmax();
 444   }
 445
 446   if(fEfficiency) delete fEfficiency;
 447   fEfficiency = new THnSparseD("fEfficiency","Efficiency - no resolution effects; p_{T}^{gen}",fDimensions,nbinsEff,xminEff,xmaxEff);
 448   fEfficiency->Sumw2();
 449   fEfficiency->GetAxis(0)->SetTitle("p_{T}^{gen}");
 450
 451   const Double_t *binArrayPt = genAxis->GetXbins()->GetArray();
 452   fEfficiency->SetBinEdges(0,binArrayPt);
 453
 454 }
 455
 456 //--------------------------------------------------------------------------------------------------------------------------------------------------
 457 void AliAnaChargedJetResponseMaker::InitializeResponseMatrixFine() {
 458   //
 459   //Initialize fine response matrix
 460   //
 461
 462   Int_t nbinsFine[fDimensions*2];
 463   Double_t xminFine[fDimensions*2];
 464   Double_t xmaxFine[fDimensions*2];
 465   Double_t pTarrayFine[fDimensions*2];
 466
 467   nbinsFine[fDimRec] = fResponseMatrix->GetAxis(fDimRec)->GetNbins()*fFineFrac;
 468   nbinsFine[fDimGen] = fResponseMatrix->GetAxis(fDimRec)->GetNbins()*fFineFrac;
 469   xminFine[fDimRec] = TMath::Min(fPtMin,0.);
 470   xminFine[fDimGen] = TMath::Min(fPtMin,0.);
 471   xmaxFine[fDimRec] = fResponseMatrix->GetAxis(fDimGen)->GetXmax();//+40.;
 472   xmaxFine[fDimGen] = fResponseMatrix->GetAxis(fDimGen)->GetXmax();//+40.;
 473   pTarrayFine[fDimRec] = 0.;
 474   pTarrayFine[fDimGen] = 0.;
 475
 476   Double_t binWidth[2];
 477   binWidth[fDimRec] = fResponseMatrix->GetAxis(fDimRec)->GetBinWidth(1);
 478
 479   Double_t binWidthFine[2];
 480   binWidthFine[fDimRec] = binWidth[fDimRec]/((double)fFineFrac);
 481   binWidthFine[fDimGen] = binWidthFine[fDimRec]*(double)fBinWidthFactorUnfolded;
 482   //  binWidthFine[fDimRec] = binWidthFine[fDimGen];
 483
 484   nbinsFine[fDimRec] = (int)((xmaxFine[fDimRec]-xminFine[fDimRec])/binWidthFine[fDimRec]); //adjust nbins to bin width
 485   nbinsFine[fDimGen] = (int)((xmaxFine[fDimGen]-xminFine[fDimGen])/binWidthFine[fDimGen]); //adjust nbins to bin width
 486
 487   printf("xminFine[fDimRec]: %f  xminFine[fDimGen]: %f \n",xminFine[fDimRec],xminFine[fDimGen]);
 488   printf("xmaxFine[fDimRec]: %f  xmaxFine[fDimGen]: %f \n",xmaxFine[fDimRec],xmaxFine[fDimGen]);
 489   printf("nbinsFine[fDimRec]: %d  nbinsFine[fDimGen]: %d \n",nbinsFine[fDimRec],nbinsFine[fDimGen]);
 490   printf("binWidthFine[fDimRec]: %f  binWidthFine[fDimGen]: %f \n",binWidthFine[fDimRec],binWidthFine[fDimGen]);
 491
 492
 493   Double_t binArrayPt0Fine[nbinsFine[fDimRec]+1];
 494   Double_t binArrayPt1Fine[nbinsFine[fDimGen]+1];
 495
 496   for(int i=0; i<nbinsFine[fDimRec]; i++) {
 497     binArrayPt0Fine[i] = xminFine[fDimRec]+(double)i*binWidthFine[fDimRec];
 498     //    printf("RM. i = %d \t binArrayPtFine[i] = %f \n",i,binArrayPt0Fine[i]);
 499   }
 500   binArrayPt0Fine[nbinsFine[fDimRec]]= xmaxFine[fDimRec];
 501
 502   for(int i=0; i<nbinsFine[fDimGen]; i++) {
 503     binArrayPt1Fine[i] = xminFine[fDimGen]+(double)i*binWidthFine[fDimGen];
 504     //    printf("RM. i = %d \t binArrayPtFine[i] = %f \n",i,binArrayPt1Fine[i]);
 505   }
 506   binArrayPt1Fine[nbinsFine[fDimGen]]= xmaxFine[fDimGen];
 507
 508   // Response matrix : dimensions must be 2N = 2 x (number of variables)
 509   // dimensions 0 ->  N-1 must be filled with reconstructed values
 510   // dimensions N -> 2N-1 must be filled with generated values
 511   if(fResponseMatrixFine) delete fResponseMatrixFine;
 512   fResponseMatrixFine = new THnSparseD("fResponseMatrixFine","Response Matrix pTMC vs pTrec",fDimensions*2,nbinsFine,xminFine,xmaxFine);
 513   fResponseMatrixFine->Sumw2();
 514   fResponseMatrixFine->GetAxis(fDimRec)->SetTitle("p_{T}^{rec}");
 515   fResponseMatrixFine->GetAxis(fDimGen)->SetTitle("p_{T}^{gen}");
 516
 517   fResponseMatrixFine->SetBinEdges(fDimRec,binArrayPt0Fine);
 518   fResponseMatrixFine->SetBinEdges(fDimGen,binArrayPt1Fine);
 519
 520   Int_t bin[2] = {1,1};
 521   for(int i=1; i<fResponseMatrixFine->GetAxis(0)->GetNbins(); i++) {
 522     bin[0]=i;
 523     for(int j=1; j<fResponseMatrixFine->GetAxis(1)->GetNbins(); j++) {
 524     bin[1]=j;
 525       fResponseMatrixFine->SetBinContent(bin,0.);
 526     }
 527   }
 528
 529 }
 530
 531 //--------------------------------------------------------------------------------------------------------------------------------------------------
 532 void AliAnaChargedJetResponseMaker::InitializeEfficiencyFine() {
 533   //
 534   // Define dimensions of efficiency THnSparse
 535   //
 536
 537   if(!fResponseMatrixFine) {
 538     printf("AliAnaChargedJetResponseMaker::InitializeEfficiencyFine()\n");
 539     printf("Can not define dimensions efficiency without fResponseMatrixFine dimensions. Aborting \n");
 540     return;
 541   }
 542
 543   TAxis *genAxis = fResponseMatrixFine->GetAxis(fDimGen);
 544
 545   Int_t nbinsEff[fDimensions];
 546   Double_t xminEff[fDimensions];
 547   Double_t xmaxEff[fDimensions];
 548
 549   for(int dim = 0; dim<fDimensions; dim++) {
 550     nbinsEff[dim] = genAxis->GetNbins();
 551     xminEff[dim]  = genAxis->GetXmin();
 552     xmaxEff[dim]  = genAxis->GetXmax();
 553   }
 554
 555   if(fEfficiencyFine) delete fEfficiencyFine;
 556   fEfficiencyFine = new THnSparseD("fEfficiencyFine","EfficiencyFine - no resolution effects; p_{T}^{gen}",fDimensions,nbinsEff,xminEff,xmaxEff);
 557   fEfficiencyFine->Sumw2();
 558   fEfficiencyFine->GetAxis(0)->SetTitle("p_{T}^{gen}");
 559
 560   const Double_t *binArrayPt = genAxis->GetXbins()->GetArray();
 561   fEfficiencyFine->SetBinEdges(0,binArrayPt);
 562
 563 }
 564
 565 //--------------------------------------------------------------------------------------------------------------------------------------------------
 566 void AliAnaChargedJetResponseMaker::FillResponseMatrixFineAndMerge() {
 567   //
 568   // Fill fine response matrix
 569   //
 570
 571   if(!fResponseMatrix) {
 572     printf("Dimensions of fResponseMatrix have to be defined first. Aborting!");
 573     return;
 574   }
 575
 576   if(!fResponseMatrixFine) {
 577     printf("Dimensions of fResponseMatrixFine have to be defined first. Aborting!");
 578     return;
 579   }
 580
 581   TAxis *genAxis = fResponseMatrixFine->GetAxis(fDimGen);
 582   TAxis *recAxis = fResponseMatrixFine->GetAxis(fDimRec);
 583
 584   Int_t nbinsFine[fDimensions*2];
 585   Double_t xminFine[fDimensions*2];
 586   Double_t xmaxFine[fDimensions*2];
 587   Double_t pTarrayFine[fDimensions*2];
 588
 589   nbinsFine[fDimGen] = genAxis->GetNbins();
 590   nbinsFine[fDimRec] = recAxis->GetNbins();
 591   xminFine[fDimGen]  = genAxis->GetXmin();
 592   xminFine[fDimRec]  = recAxis->GetXmin();
 593   xmaxFine[fDimGen]  = genAxis->GetXmax();
 594   xmaxFine[fDimRec]  = recAxis->GetXmax();
 595   pTarrayFine[fDimGen] = 0.;
 596   pTarrayFine[fDimRec] = 0.;
 597
 598   double sumyield = 0.;
 599   double sumyielderror2 = 0.;
 600
 601   int ipt[2]    = {0,0};
 602   int iptMerged[2]    = {0,0};
 603   int ierror[2] = {0,0};
 604
 605   Int_t check = 0;
 606   Double_t pTgen, pTrec;
 607   Double_t yield = 0.;
 608   Double_t error = 0.;
 609
 610   const int nng = fResponseMatrix->GetAxis(fDimGen)->GetNbins();
 611   const int nnr = fResponseMatrix->GetAxis(fDimRec)->GetNbins();
 612   Double_t errorArray[nng][nnr];
 613   for(int iig =0; iig<nng; iig++) {
 614     for(int iir =0; iir<nnr; iir++) {
 615       errorArray[iig][iir] = 0.;
 616     }
 617   }
 618
 619   for(int iy=1; iy<=nbinsFine[fDimGen]; iy++) { //gen
 620     pTgen = fResponseMatrixFine->GetAxis(fDimGen)->GetBinCenter(iy);
 621     pTarrayFine[fDimGen] = pTgen;
 622     ierror[fDimGen]=iy;
 623     sumyield = 0.;
 624     check = 0;
 625
 626     for(int ix=1; ix<=nbinsFine[fDimRec]; ix++) { //rec
 627       pTrec = fResponseMatrixFine->GetAxis(fDimRec)->GetBinCenter(ix);
 628       Double_t width = fResponseMatrixFine->GetAxis(fDimRec)->GetBinWidth(ix);
 629       if(fResolutionType==kParam) {
 630         yield = fDeltaPt->Eval(pTrec-pTgen);
 631         error = 0.;
 632       }
 633       else if(fResolutionType==kResiduals) {
 634         yield = fhDeltaPt->Interpolate(pTrec-pTgen);
 635         error = 0.;
 636       }
 637       else if(fResolutionType==kResidualsErr) {
 638         Double_t deltaPt = pTrec-pTgen;
 639         int bin = fhDeltaPt->FindBin(deltaPt);
 640         yield = fhDeltaPt->GetBinContent(bin);
 641         error = fhDeltaPt->GetBinError(bin);
 642       }
 643       yield=yield*width;
 644       error=error*width;
 645       //avoid trouble with empty bins in the high pT tail of deltaPt distribution
 646       if(check==0 && yield>0. && pTrec>pTgen) check=1;
 647       if(check==1 && yield==0.) ix=nbinsFine[fDimRec];
 648
 649       sumyield+=yield;
 650       sumyielderror2 += error*error;
 651
 652       pTarrayFine[fDimRec] = pTrec;
 653       ierror[fDimRec]  = ix;
 654       fResponseMatrixFine->Fill(pTarrayFine,yield);
 655       if(fbCalcErrors) fResponseMatrixFine->SetBinError(ierror,error);
 656
 657     }// ix (dimRec) loop
 658
 659     //Normalize to total number of counts =1
 660     if(sumyield>1) {
 661       ipt[fDimGen]=iy;
 662       for(int ix=1; ix<=nbinsFine[fDimRec]; ix++) {
 663         ipt[fDimRec]=ix;
 664         fResponseMatrixFine->SetBinContent(ipt,fResponseMatrixFine->GetBinContent(ipt)/sumyield);
 665         if(fResolutionType==kResidualsErr && fbCalcErrors) {
 666           Double_t A = 1./sumyield*fResponseMatrix->GetBinError(ipt);
 667           Double_t B = -1.*fResponseMatrix->GetBinContent(ipt)/sumyield/sumyield*TMath::Sqrt(sumyielderror2);
 668           Double_t tmp2 = A*A + B*B;
 669           fResponseMatrix->SetBinError(ipt,TMath::Sqrt(tmp2));
 670         }
 671
 672       }
 673     }
 674
 675     int bin[1];
 676     bin[0] = iy;
 677     fEfficiencyFine->SetBinContent(bin,sumyield);
 678     if(fbCalcErrors) fEfficiencyFine->SetBinError(bin,TMath::Sqrt(sumyielderror2));
 679
 680     //fill merged response matrix
 681
 682     //find bin in fine RM corresponding to mimimum/maximum bin of merged RM on rec axis
 683     int ixMin = fResponseMatrixFine->GetAxis(fDimRec)->FindBin(fResponseMatrix->GetAxis(fDimRec)->GetXmin());
 684     int ixMax = fResponseMatrixFine->GetAxis(fDimRec)->FindBin(fResponseMatrix->GetAxis(fDimRec)->GetXmax());
 685
 686     if(fResponseMatrixFine->GetAxis(fDimGen)->GetBinLowEdge(iy) >= fResponseMatrix->GetAxis(fDimGen)->GetXmin()) {
 687       ipt[fDimGen]=iy;
 688       iptMerged[fDimGen]=fResponseMatrix->GetAxis(fDimGen)->FindBin(pTgen);
 689
 690       Double_t weight = 1.;
 691       if(f1MergeFunction) {
 692         Double_t loEdge = fResponseMatrix->GetAxis(fDimGen)->GetBinLowEdge(iptMerged[fDimGen]);
 693         Double_t upEdge = fResponseMatrix->GetAxis(fDimGen)->GetBinUpEdge(iptMerged[fDimGen]);
 694         Float_t powInteg = f1MergeFunction->Integral(loEdge,upEdge);
 695         //printf("loEdge = %f  upEdge = %f  powInteg = %f\n",loEdge,upEdge,powInteg);
 696         if(powInteg>0.)
 697           weight = f1MergeFunction->Integral(fResponseMatrixFine->GetAxis(fDimGen)->GetBinLowEdge(iy),fResponseMatrixFine->GetAxis(fDimGen)->GetBinUpEdge(iy))/powInteg;
 698         //      printf("weight: %f \n", weight );
 699       } else {
 700         weight = 1./((double)fFineFrac);
 701         if(fbVariableBinning && pTgen<fPtMaxUnfVarBinning) weight=1./(0.5*(double)fFineFrac);
 702       }
 703
 704       for(int ix=ixMin; ix<=ixMax; ix++) {                    //loop reconstructed axis
 705         pTrec = fResponseMatrixFine->GetAxis(fDimRec)->GetBinCenter(ix);
 706         ipt[fDimRec]=ix;
 707         iptMerged[fDimRec]=fResponseMatrix->GetAxis(fDimRec)->FindBin(pTrec);
 708
 709         fResponseMatrix->AddBinContent(iptMerged,fResponseMatrixFine->GetBinContent(ipt)*weight);
 710         if(fbCalcErrors) errorArray[iptMerged[fDimGen]-1][iptMerged[fDimRec]-1] += fResponseMatrixFine->GetBinError(ipt)*fResponseMatrixFine->GetBinError(ipt)*weight*weight;
 711       }
 712
 713    }//loop gen axis fine response matrix
 714
 715   } // iy (dimGen) loop
 716
 717   //Fill Efficiency corresponding to merged response matrix
 718   for(int iy=1; iy<=fResponseMatrix->GetAxis(fDimGen)->GetNbins(); iy++) { //gen
 719     sumyield = 0.;
 720     ipt[fDimGen] = iy;
 721
 722     for(int ix=1; ix<=fResponseMatrix->GetAxis(fDimRec)->GetNbins(); ix++) { //rec
 723       ipt[fDimRec] = ix;
 724       sumyield += fResponseMatrix->GetBinContent(ipt);
 725
 726       if(fbCalcErrors) fResponseMatrix->SetBinError(ipt,TMath::Sqrt(errorArray[iy][ix]));
 727     }
 728     printf("RM: pTgen: %f \t sumyield: %f \n",fResponseMatrix->GetAxis(fDimGen)->GetBinCenter(iy),sumyield);
 729     int bin[1];
 730     bin[0] = iy;
 731     fEfficiency->SetBinContent(bin,sumyield);
 732     if(fbCalcErrors) fEfficiency->SetBinError(bin,0);
 733   }
 734
 735   if(fDebug) printf("fResponseMatrixFine->GetNbins(): %d \n",(int)fResponseMatrixFine->GetNbins());
 736   if(fDebug) printf("fResponseMatrix->GetNbins(): %d \n",(int)fResponseMatrix->GetNbins());
 737
 738   if(fDebug) printf("Done constructing response matrix\n");
 739
 740 }
 741
 742 //--------------------------------------------------------------------------------------------------------------------------------------------------
 743 TH2* AliAnaChargedJetResponseMaker::MakeResponseMatrixRebin(TH2 *hRMFine, TH2 *hRM, Bool_t useFunctionWeight) {
 744
 745   //
 746   // Rebin matrix hRMFine to dimensions of hRM
 747   // function returns matrix in TH2D format (hRM2) with dimensions from hRM
 748   //
 749
 750   TH2 *hRM2 = (TH2*)hRM->Clone("hRM2");
 751   hRM2->Reset();
 752
 753   if(useFunctionWeight)  cout << "Use function to do weighting" << endl;
 754
 755   //First normalize columns of input
 756   const Int_t nbinsNorm = hRM2->GetNbinsX();
 757   if(fDebug) cout << "nbinsNorm: " << nbinsNorm << endl;
 758
 759   TArrayF *normVector = new TArrayF(nbinsNorm);
 760
 761   for(int ix=1; ix<=hRM2->GetNbinsX(); ix++) {
 762     Double_t xLow = hRM2->GetXaxis()->GetBinLowEdge(ix);
 763     Double_t xUp = hRM2->GetXaxis()->GetBinUpEdge(ix);
 764     //cout << "xLow: " << xLow << " xUp: " << xUp << "\t center: " << hRM2->GetXaxis()->GetBinCenter(ix) << endl;
 765     Int_t binxLowFine = hRMFine->GetXaxis()->FindBin(xLow);
 766     Int_t binxUpFine = hRMFine->GetXaxis()->FindBin(xUp)-1;
 767     //cout << "xLowFine: " << hRMFine->GetXaxis()->GetBinLowEdge(binxLowFine) << "\txUpFine: " << hRMFine->GetXaxis()->GetBinUpEdge(binxUpFine) << endl;
 768     if(useFunctionWeight)
 769       normVector->SetAt(f1MergeFunction->Integral(xLow,xUp),ix-1);
 770     else
 771       normVector->SetAt(hRMFine->Integral(binxLowFine,binxUpFine,1,hRMFine->GetYaxis()->GetNbins()),ix-1);
 772     //if(fDebug) cout << "ix norm: " << normVector->At(ix-1) << endl;
 773   }
 774
 775   Double_t content, oldcontent = 0.;
 776   Int_t ixNew = 0;
 777   Int_t iyNew = 0;
 778   Double_t xvalLo, xvalUp, yvalLo, yvalUp;
 779   Double_t xmin = hRM2->GetXaxis()->GetXmin();
 780   Double_t ymin = hRM2->GetYaxis()->GetXmin();
 781   Double_t xmax = hRM2->GetXaxis()->GetXmax();
 782   Double_t ymax = hRM2->GetYaxis()->GetXmax();
 783   for(int ix=1; ix<=hRMFine->GetXaxis()->GetNbins(); ix++) {
 784     xvalLo = hRMFine->GetXaxis()->GetBinLowEdge(ix);
 785     xvalUp = hRMFine->GetXaxis()->GetBinUpEdge(ix);
 786     if(xvalLo<xmin || xvalUp>xmax) continue;
 787     ixNew = hRM2->GetXaxis()->FindBin(hRMFine->GetXaxis()->GetBinCenter(ix));
 788     for(int iy=1; iy<=hRMFine->GetYaxis()->GetNbins(); iy++) {
 789       yvalLo = hRMFine->GetYaxis()->GetBinLowEdge(iy);
 790       yvalUp = hRMFine->GetYaxis()->GetBinUpEdge(iy);
 791       if(yvalLo<ymin || yvalUp>ymax) continue;
 792       content = hRMFine->GetBinContent(ix,iy);
 793       iyNew = hRM2->GetYaxis()->FindBin(hRMFine->GetYaxis()->GetBinCenter(iy));
 794       oldcontent = hRM2->GetBinContent(ixNew,iyNew);
 795
 796       //if(fDebug) cout << "ixNew: " << ixNew << " " << xvalLo << " iyNew: " << iyNew << " " << yvalLo << " content: " << content << " oldContent: " << oldcontent << " newContent: " << oldcontent+content << endl;
 797       Double_t weight = 1.;
 798       if(normVector->At(ixNew-1)>0.) {
 799         if(useFunctionWeight)
 800           weight = f1MergeFunction->Integral(xvalLo,xvalUp)/normVector->At(ixNew-1);
 801         else
 802           weight = 1./normVector->At(ixNew-1);
 803       }
 804       hRM2->SetBinContent(ixNew,iyNew,oldcontent+content*weight);
 805     }
 806   }
 807
 808   if(normVector) delete normVector;
 809
 810   return hRM2;
 811
 812 }
 813
 814 //--------------------------------------------------------------------------------------------------------------------------------------------------
 815 TH2* AliAnaChargedJetResponseMaker::CreateTruncated2DHisto(TH2 *h2, Double_t xmin, Double_t xmax, Double_t ymin, Double_t ymax) {
 816
 817   //
 818   // Limit axis range of 2D histogram
 819   //
 820
 821   Int_t binMinXh2 = h2->GetXaxis()->FindBin(xmin);
 822   if(h2->GetXaxis()->GetBinLowEdge(binMinXh2) < xmin ) binMinXh2++;
 823   if(h2->GetXaxis()->GetBinLowEdge(binMinXh2) > xmin ) binMinXh2--;
 824
 825   Int_t binMinYh2 = h2->GetYaxis()->FindBin(ymin);
 826   if(h2->GetYaxis()->GetBinLowEdge(binMinYh2) < ymin ) binMinYh2++;
 827   if(h2->GetYaxis()->GetBinLowEdge(binMinYh2) > ymin ) binMinYh2--;
 828
 829   Int_t binMaxXh2 = h2->GetXaxis()->FindBin(xmax);
 830   if(h2->GetXaxis()->GetBinUpEdge(binMaxXh2) < xmax ) binMaxXh2++;
 831   if(h2->GetXaxis()->GetBinUpEdge(binMaxXh2) > xmax ) binMaxXh2--;
 832
 833   Int_t binMaxYh2 = h2->GetYaxis()->FindBin(ymax);
 834   if(h2->GetYaxis()->GetBinUpEdge(binMaxYh2) < ymax ) binMaxYh2++;
 835   if(h2->GetYaxis()->GetBinUpEdge(binMaxYh2) > ymax ) binMaxYh2--;
 836
 837   Int_t nbinsX = binMaxXh2-binMinXh2+1;
 838   Int_t nbinsY = binMaxYh2-binMinYh2+1;
 839
 840   Double_t *binsX = new Double_t[nbinsX+1];
 841   Double_t *binsY = new Double_t[nbinsY+1];
 842
 843   for(int ix=1; ix<=nbinsX; ix++)
 844     binsX[ix-1] = h2->GetXaxis()->GetBinLowEdge(binMinXh2+ix-1);
 845   binsX[nbinsX] = h2->GetXaxis()->GetBinUpEdge(binMaxXh2);
 846
 847   for(int iy=1; iy<=nbinsY; iy++)
 848     binsY[iy-1] = h2->GetYaxis()->GetBinLowEdge(binMinYh2+iy-1);
 849   binsY[nbinsY] = h2->GetYaxis()->GetBinUpEdge(binMaxYh2);
 850
 851   TH2 *h2Lim = new TH2D("h2Lim","h2Lim",nbinsX,binsX,nbinsY,binsY);
 852
 853   for(int ix=1; ix<=nbinsX; ix++) {
 854     //    cout << "ix: " << ix << "  " << binsX[ix] << endl;
 855     for(int iy=1; iy<=nbinsY; iy++) {
 856       // cout << "ix: " << ix << "  " << binsX[ix] << "\tiy: " << iy << "  " << binsY[iy] << endl;
 857
 858       double content = h2->GetBinContent(binMinXh2+ix-1,binMinYh2+iy-1);
 859       double error = h2->GetBinContent(binMinXh2+ix-1,binMinYh2+iy-1);
 860       h2Lim->SetBinContent(ix,iy,content);
 861       if(fbCalcErrors) h2Lim->SetBinError(ix,iy,error);
 862
 863     }
 864   }
 865
 866
 867
 868   return h2Lim;
 869 }
 870
 871 //--------------------------------------------------------------------------------------------------------------------------------------------------
 872 TH2* AliAnaChargedJetResponseMaker::TruncateAxisRangeResponseMatrix(TH2 *hRMOrig, Double_t xmin, Double_t xmax, Double_t ymin, Double_t ymax) {
 873
 874   //
 875   // Limit axis range of response matrix without changing normalization
 876   //
 877
 878   //TH2 *hRMLimit
 879   //TH2 *hRMLimit2 = (TH2*)hRMLimit->Clone("hRMLimit2");
 880
 881   TH2 *hRMLimit2 = CreateTruncated2DHisto(hRMOrig, xmin, xmax, ymin, ymax);
 882   hRMLimit2->Reset();
 883
 884   double binCent[2] = {0.,0.};
 885   double content = 0.;
 886   double error = 0.;
 887   Int_t binOrig[2] = {0};
 888   for(int ix=1; ix<=hRMLimit2->GetXaxis()->GetNbins(); ix++) {
 889     binCent[0] = hRMLimit2->GetXaxis()->GetBinCenter(ix);
 890     binOrig[0] = hRMOrig->GetXaxis()->FindBin(binCent[0]);
 891     for(int iy=1; iy<=hRMLimit2->GetYaxis()->GetNbins(); iy++) {
 892       binCent[1] = hRMLimit2->GetYaxis()->GetBinCenter(iy);
 893       binOrig[1] = hRMOrig->GetYaxis()->FindBin(binCent[1]);
 894
 895       content = hRMOrig->GetBinContent(binOrig[0],binOrig[1]);
 896       error = hRMOrig->GetBinError(binOrig[0],binOrig[1]);
 897
 898       hRMLimit2->SetBinContent(ix,iy,content);
 899       if(fbCalcErrors) hRMLimit2->SetBinError(ix,iy,error);
 900
 901     }
 902   }
 903
 904
 905   return hRMLimit2;
 906
 907 }
 908
 909 //--------------------------------------------------------------------------------------------------------------------------------------------------
 910 Bool_t AliAnaChargedJetResponseMaker::CheckInputForCombinedResponse() {
 911
 912   Bool_t check = kTRUE;
 913
 914   if(!fhDeltaPt)             check = kFALSE;
 915   if(!fh2DetectorResponse)   check = kFALSE;
 916   if(!fhEfficiencyDet)       check = kFALSE;
 917   if(!fh1MeasuredTruncated)  check = kFALSE;
 918   if(fPtMin<0. || fPtMax<0.) check = kFALSE;
 919
 920   return check;
 921 }
 922
 923 //--------------------------------------------------------------------------------------------------------------------------------------------------
 924 void AliAnaChargedJetResponseMaker::MakeResponseMatrixCombined(Int_t skipBins, Int_t binWidthFactor, Int_t extraBins, Bool_t bVariableBinning, Double_t ptmin) {
 925
 926   //Check if all input histos are there otherwise break
 927   if(!CheckInputForCombinedResponse()) {
 928     printf("Not all input available..... \n");
 929     return;
 930   }
 931
 932   // Make response bkg fluctuations
 933   MakeResponseMatrixJetsFineMerged(skipBins,binWidthFactor,extraBins,bVariableBinning,ptmin);
 934
 935   //Get TH2D with dimensions we want final RM
 936   TH2D *h2ResponseMatrix = fResponseMatrix->Projection(0,1,"E");
 937   h2ResponseMatrix->Reset();
 938
 939   //Get fine binned response matrix from bkg fluctuations
 940   TH2D *h2ResponseMatrixFine = fResponseMatrixFine->Projection(0,1,"E");
 941
 942   // Rebin response detector effects
 943   const TArrayD *arrayX = h2ResponseMatrixFine->GetXaxis()->GetXbins();
 944   const Double_t *xbinsArray = arrayX->GetArray();
 945
 946   TH2D *h2ResponseDetTmp = new TH2D("h2ResponseDetTmp","h2ResponseDetTmp",h2ResponseMatrixFine->GetNbinsX(),xbinsArray,h2ResponseMatrixFine->GetNbinsX(),xbinsArray);
 947
 948   fh2DetectorResponseRebin = (TH2D*)MakeResponseMatrixRebin(fh2DetectorResponse,h2ResponseDetTmp,kFALSE);
 949   fh2DetectorResponseRebin->SetName("fh2DetectorResponseRebin");
 950   fh2DetectorResponseRebin->SetTitle("fh2DetectorResponseRebin");
 951   fh2DetectorResponseRebin->SetXTitle("p_{T}^{jet} gen");
 952   fh2DetectorResponseRebin->SetYTitle("p_{T}^{jet} rec");
 953
 954   // Make full combined response matrix fine binning (background flucutuations + detector effects)
 955   fh2ResponseMatrixCombinedFineFull = (TH2D*)MultiplityResponseMatrices(h2ResponseMatrixFine,fh2DetectorResponseRebin);
 956   fh2ResponseMatrixCombinedFineFull->SetName("fh2ResponseMatrixCombinedFineFull");
 957   fh2ResponseMatrixCombinedFineFull->SetTitle("fh2ResponseMatrixCombinedFineFull");
 958
 959   // Rebin full combined response matrix with weighting
 960   fh2ResponseMatrixCombinedFull = (TH2D*)MakeResponseMatrixRebin(fh2ResponseMatrixCombinedFineFull,h2ResponseMatrix,kTRUE);
 961   fh2ResponseMatrixCombinedFull->SetName("fh2ResponseMatrixCombinedFull");
 962   fh2ResponseMatrixCombinedFull->SetTitle("fh2ResponseMatrixCombinedFull");
 963   fh2ResponseMatrixCombinedFull->SetXTitle("#it{p}_{T,gen}^{ch jet} (GeV/#it{c})");
 964   fh2ResponseMatrixCombinedFull->SetYTitle("#it{p}_{T,rec}^{ch jet} (GeV/#it{c})");
 965
 966   fh2ResponseMatrixCombined = (TH2D*)TruncateAxisRangeResponseMatrix(fh2ResponseMatrixCombinedFull, h2ResponseMatrix->GetXaxis()->GetXmin(),h2ResponseMatrix->GetXaxis()->GetXmax(),fh1MeasuredTruncated->GetXaxis()->GetXmin(),fh1MeasuredTruncated->GetXaxis()->GetXmax());
 967   fh2ResponseMatrixCombined->SetTitle("fh2ResponseMatrixCombined");
 968   fh2ResponseMatrixCombined->SetName("fh2ResponseMatrixCombined");
 969
 970   fhEfficiencyCombined = (TH1D*)fh2ResponseMatrixCombined->ProjectionX("fhEfficiencyCombined");
 971   fhEfficiencyCombined->Reset();
 972
 973   for(int i=1; i<=fh2ResponseMatrixCombined->GetNbinsX(); i++) {
 974     Double_t sumyield = 0.;
 975     for(int j=1; j<=fh2ResponseMatrixCombined->GetYaxis()->GetNbins(); j++) {
 976       sumyield+=fh2ResponseMatrixCombined->GetBinContent(i,j);
 977     }
 978     Double_t kCounter = 0.;
 979     Double_t kPtLowEdge = fh2ResponseMatrixCombined->GetXaxis()->GetBinLowEdge(i);
 980     Double_t kPtUpEdge  = fh2ResponseMatrixCombined->GetXaxis()->GetBinUpEdge(i);
 981     Double_t kJetEffDet = 0.;
 982
 983     for(int k=1; k<=fhEfficiencyDet->GetNbinsX(); k++) {
 984       if(fhEfficiencyDet->GetXaxis()->GetBinLowEdge(k)>=kPtLowEdge && fhEfficiencyDet->GetXaxis()->GetBinUpEdge(k)<=kPtUpEdge) {
 985         kJetEffDet+=fhEfficiencyDet->GetBinContent(k);
 986         kCounter+=1.;
 987       }
 988     }
 989     kJetEffDet = kJetEffDet/kCounter;
 990
 991     if(kJetEffDet==0.) kJetEffDet=1.;
 992
 993     fhEfficiencyCombined->SetBinContent(i,sumyield*kJetEffDet);
 994
 995   }
 996
 997 }
 998
 999 //--------------------------------------------------------------------------------------------------------------------------------------------------
1000 TH2* AliAnaChargedJetResponseMaker::MultiplityResponseMatrices(TH2 *h2RMDeltaPt, TH2 *h2RMDetector) {
1001
1002   //
1003   // Make combined response matrix (background flucutuations + detector effects)
1004   //
1005   // hRMDeltaPt is the response matrix for background fluctuations from \delta(p_t) measurement
1006   // 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
1007   //
1008   // Function assumes that generated/unfolded axis is x-axis and reconstructed is on y-axis on both respone matrices
1009
1010
1011   TH2D *h2ResponseMatrixCombined = (TH2D*)h2RMDeltaPt->Clone("h2ResponseMatrixCombined"); //h2RMDeltaPt is the bkg fluctuations RM which has the dimensions we want for the combined response matrix
1012   h2ResponseMatrixCombined->SetTitle("h2ResponseMatrixCombined");
1013   h2ResponseMatrixCombined->SetName("h2ResponseMatrixCombined");
1014
1015   // M = RM_deltaPt * RM_DetEffects * T   (M=measured T=truth)
1016   // Dimensions:
1017   // mx1 = mxd * dxt * tx1
1018   // m = measured bins
1019   // t = truth bins
1020   // d = rec for RM_DetEffects and gen for RM_deltaPt
1021   // RM_comb = RM_deltaPt * RM_DetEffects (dimensions mxt)
1022   // RM_comb(m,t) = Sum_d RM_deltaPt(m,d)*RM_DetEffects(d,t)
1023
1024   if(fDebug) {
1025     printf("Nt=%d\n",h2ResponseMatrixCombined->GetNbinsX());
1026     printf("Nm=%d\n",h2ResponseMatrixCombined->GetNbinsY());
1027     printf("Nd=%d\n",h2RMDeltaPt->GetNbinsX());
1028   }
1029
1030
1031   for(Int_t t=1; t<=h2ResponseMatrixCombined->GetNbinsX();t++) {
1032     for(Int_t m=1; m<=h2ResponseMatrixCombined->GetNbinsY();m++) {
1033       Double_t valueSum = 0.;
1034       for(Int_t d=1; d<=h2RMDeltaPt->GetNbinsX();d++) {
1035         valueSum += h2RMDeltaPt->GetBinContent(d,m) * h2RMDetector->GetBinContent(t,d);
1036         //      if(t==10 && m==10) cout << "sum m,d=" << m << "," << d << endl;
1037       }//d-loop
1038       //  if(t==10) cout << "t,m = " << t << "," << m << "\tvalueSum: " << valueSum << endl;
1039       h2ResponseMatrixCombined->SetBinContent(t,m,valueSum);
1040     } //m-loop
1041   }//t-loop
1042
1043   return h2ResponseMatrixCombined;
1044
1045 }
1046
1047 //--------------------------------------------------------------------------------------------------------------------------------------------------
1048 TH2* AliAnaChargedJetResponseMaker::GetTransposeResponsMatrix(TH2 *h2RM) {
1049   //
1050   // Transpose response matrix
1051   //
1052
1053   //Initialize transposed response matrix h2RMTranspose
1054   TArrayD *arrayX = (TArrayD*)h2RM->GetXaxis()->GetXbins();
1055   for(int i=0; i<arrayX->GetSize(); i++) cout << "i: " << arrayX->At(i) << endl;
1056   Double_t *xbinsArray = arrayX->GetArray();
1057
1058   TArrayD *arrayY = (TArrayD*)h2RM->GetYaxis()->GetXbins();
1059   for(int i=0; i<arrayY->GetSize(); i++) cout << "i: " << arrayY->At(i) << endl;
1060   Double_t *ybinsArray = arrayY->GetArray();
1061
1062   TH2D *h2RMTranspose = new TH2D("h2RMTranspose","h2RMTranspose",h2RM->GetNbinsY(),ybinsArray,h2RM->GetNbinsX(),xbinsArray);
1063
1064   //Fill tranposed response matrix
1065   for (Int_t ibin = 1; ibin <= h2RMTranspose->GetNbinsX(); ibin++) {
1066     for (Int_t jbin = 1; jbin <= h2RMTranspose->GetNbinsY(); jbin++) {
1067       h2RMTranspose->SetBinContent(ibin,jbin, h2RM->GetBinContent(jbin,ibin));
1068     }
1069   }
1070
1071
1072   return h2RMTranspose;
1073
1074 }
1075
1076 //--------------------------------------------------------------------------------------------------------------------------------------------------
1077 TH2* AliAnaChargedJetResponseMaker::NormalizeResponsMatrixYaxisWithPrior(TH2 *h2RM, TH1 *hPrior) {
1078   //
1079   // Normalize such that the Y projection is the prior
1080   //
1081
1082   //  TH1D *hProjRespY = (TH1D*)h2RM->ProjectionY("hProjRespY");
1083   double intPrior = hPrior->Integral();//"width");
1084   for (Int_t jbin = 1; jbin <= h2RM->GetNbinsY(); jbin++) {
1085     //    double corr = hPrior->GetBinContent(jbin)/hProjRespY->GetBinContent(jbin);
1086       for (Int_t ibin = 1; ibin <= h2RM->GetNbinsX(); ibin++) {
1087         double content = h2RM->GetBinContent(ibin,jbin);
1088         //      h2RM->SetBinContent(ibin,jbin,content*corr);
1089         h2RM->SetBinContent(ibin,jbin,hPrior->GetBinContent(jbin)/hPrior->GetBinWidth(jbin)/intPrior*content);
1090     }
1091   }
1092
1093   return h2RM;
1094 }
1095
1096 //--------------------------------------------------------------------------------------------------------------------------------------------------
1097 TH1D* AliAnaChargedJetResponseMaker::MultiplyResponseGenerated(TH1 *hGen, TH2 *hResponse,TH1 *hEfficiency,Bool_t bDrawSlices) {
1098
1099   //
1100   // Multiply hGen with response matrix to obtain refolded spectrum
1101   // Efficiency must be given. In case efficiency is 1 use SetFlatEfficiency(1.) before calling function
1102   //
1103
1104   if(!hEfficiency) {
1105     //    printf("Efficiency must be given. In case efficiency is 1 use SetFlatEfficiency(1.) before calling function. Aborting!");
1106     //    return 0;
1107     printf("Setting efficiency to 1 \n");
1108     hEfficiency = (TH1D*)hGen->Clone("hEfficiency");
1109     hEfficiency->Reset();
1110     for(int i=1; i<=hEfficiency->GetNbinsX(); i++) hEfficiency->SetBinContent(i,1.);
1111   }
1112
1113   //For response
1114   //x-axis: generated
1115   //y-axis: reconstructed
1116   if(fDebug>0) cout << "nbins hGen: " << hGen->GetNbinsX() << "\t nbins hResponseGen: " << hResponse->GetXaxis()->GetNbins() << "\t nbins hResponseRec: " << hResponse->GetYaxis()->GetNbins()  << endl;
1117
1118   TH1D *hRec = hResponse->ProjectionY("hRec");
1119   hRec->Sumw2();
1120   hRec->Reset();
1121   hRec->SetTitle("hRec");
1122   hRec->SetName("hRec");
1123
1124   for(int irec=1; irec<=hRec->GetNbinsX(); irec++)
1125     hRec->SetBinContent(irec,0);
1126
1127   TH1D *hRecGenBin = 0x0;
1128   TCanvas *cSlices = 0x0;
1129   if(bDrawSlices) {
1130     cSlices = new TCanvas("cSlices","cSlices:Slices",400,400);
1131     gPad->SetLogy();
1132   }
1133
1134   Double_t yieldMC = 0.;
1135   Double_t yieldMCerror = 0.;
1136   Double_t sumYield = 0.;
1137   const Int_t nbinsRec = hRec->GetNbinsX();
1138   Double_t sumError2[nbinsRec+1];
1139   for(int irec=0; irec<=hRec->GetNbinsX(); irec++)
1140     sumError2[irec]=0.;
1141   Double_t eff = 0.;
1142
1143   for(int igen=1; igen<=hGen->GetNbinsX(); igen++) {
1144     //get pTMC
1145     sumYield = 0.;
1146     if(fEffFlat>0.)
1147       eff = fEffFlat;
1148     else
1149       eff = hEfficiency->GetBinContent(igen);
1150     yieldMC = hGen->GetBinContent(igen)*eff;
1151     yieldMCerror = hGen->GetBinError(igen)*eff;
1152
1153     if(bDrawSlices) {
1154       hRecGenBin = hResponse->ProjectionY(Form("hRecGenBin%d",igen));
1155       hRecGenBin->Sumw2();
1156       hRecGenBin->Reset();
1157       hRecGenBin->SetTitle(Form("hRecGenBin%d",igen));
1158       hRecGenBin->SetName(Form("hRecGenBin%d",igen));
1159     }
1160
1161     for(int irec=1; irec<=hRec->GetNbinsX(); irec++) {
1162       hRec->AddBinContent(irec,yieldMC*hResponse->GetBinContent(igen,irec));
1163       sumYield+=hResponse->GetBinContent(igen,irec);
1164       //      Double_t A = yieldMC*hResponse->GetBinError(igen,irec);
1165       Double_t B = hResponse->GetBinContent(igen,irec)*yieldMCerror;
1166       //      Double_t tmp2 = A*A + B*B;
1167       //sumError2[irec-1] += tmp2 ;
1168       sumError2[irec-1] += B*B;
1169
1170       if(bDrawSlices)
1171         hRecGenBin->SetBinContent(irec,yieldMC*hResponse->GetBinContent(igen,irec));
1172
1173     }
1174     if(bDrawSlices) {
1175       cSlices->cd();
1176       hRecGenBin->SetLineColor(igen);
1177       if(igen==1) hRecGenBin->DrawCopy();
1178       else hRecGenBin->DrawCopy("same");
1179     }
1180
1181     if(hRecGenBin) delete hRecGenBin;
1182
1183     cout << "igen: " << igen << "\tpTMC: " << hGen->GetXaxis()->GetBinCenter(igen) << "\teff:" << eff << "\tsumYield: " << sumYield << endl;
1184   }
1185
1186   for(int i=0; i<=nbinsRec; i++) {
1187     if(sumError2[i]>0.)
1188       hRec->SetBinError(i+1,TMath::Sqrt(sumError2[i]));
1189   }
1190
1191
1192   return hRec;
1193 }
1194
1195 //--------------------------------------------------------------------------------------------------------------------------------------------------
1196 TH1D* AliAnaChargedJetResponseMaker::MultiplyResponseGenerated(TF1 *fGen, TH2 *hResponse,TH1 *hEfficiency) {
1197
1198   //
1199   // Multiply fGen function with response matrix to obtain (re)folded spectrum
1200   // Efficiency must be given. In case efficiency is 1 use SetFlatEfficiency(1.) before calling function
1201   //
1202
1203   //For response
1204   //x-axis: generated
1205   //y-axis: reconstructed
1206
1207   if(fDebug>0) printf(">>AliAnaChargedJetResponseMaker::MultiplyResponseGenerated(TF1 *fGen, TH2 *hResponse,TH1 *hEfficiency)\n");
1208
1209   TH1D *hRec = hResponse->ProjectionY("hRec");
1210   hRec->Sumw2();
1211   hRec->Reset();
1212   hRec->SetTitle("hRec");
1213   hRec->SetName("hRec");
1214
1215   //  TH1D *hRec = new TH1D("hRec","hRec",hResponse->GetNbinsY(),hResponse->GetYaxis()->GetXmin(),hResponse->GetYaxis()->GetXmax());
1216
1217   for(int irec=1; irec<=hRec->GetNbinsX(); irec++)
1218     hRec->SetBinContent(irec,0);
1219
1220   Double_t yieldMC = 0.;
1221   Double_t sumYield = 0.;
1222   Double_t eff = 0.;
1223   for(int igen=1; igen<=hResponse->GetNbinsX(); igen++) {
1224     //get pTMC
1225     sumYield = 0.;
1226     double pTMC = hResponse->GetXaxis()->GetBinCenter(igen);
1227     if(hEfficiency) {
1228       int binEff = hEfficiency->FindBin(pTMC);
1229       eff = hEfficiency->GetBinContent(binEff);
1230     }
1231     else eff = 1.;
1232     // yieldMC = fGen->Eval(pTMC)*eff;
1233     yieldMC = fGen->Integral(hResponse->GetXaxis()->GetBinLowEdge(igen),hResponse->GetXaxis()->GetBinUpEdge(igen))*eff;
1234     for(int irec=1; irec<=hResponse->GetNbinsY(); irec++) {
1235       hRec->AddBinContent(irec,yieldMC*hResponse->GetBinContent(igen,irec));
1236       sumYield+=hResponse->GetBinContent(igen,irec);
1237     }
1238     //    cout << "igen: " << igen << "\tpTMC: " << pTMC << "\tsumYield: " << sumYield << endl;
1239     //    cout << "yieldMC: " << yieldMC << endl;
1240
1241   }
1242
1243   return hRec;
1244 }
1245
1246 //--------------------------------------------------------------------------------------------------------------------------------------------------
1247 /* static */ Double_t AliAnaChargedJetResponseMaker::GetBetaPerDOFValue(Int_t betaColl, Int_t betaOpt) {
1248
1249   Float_t betaPerDOFoptions[6] = {0.};
1250   //define beta per degree of freedom (DOF=nbins in measured)
1251   if(betaColl==0) {
1252     betaPerDOFoptions[0] = 9.1e-9;
1253     betaPerDOFoptions[1] = 2.25e-8;
1254     betaPerDOFoptions[2] = 4.55e-8;
1255     betaPerDOFoptions[3] = 9.1e-8;
1256     betaPerDOFoptions[4] = 2.25e-7;
1257     betaPerDOFoptions[5] = 4.55e-7;
1258   }
1259   if(betaColl==1) {
1260     betaPerDOFoptions[0] = 9.1e-7;
1261     betaPerDOFoptions[1] = 2.25e-6;
1262     betaPerDOFoptions[2] = 4.55e-6;
1263     betaPerDOFoptions[3] = 9.1e-6;
1264     betaPerDOFoptions[4] = 2.25e-5;
1265     betaPerDOFoptions[5] = 4.55e-5;
1266   }
1267   if(betaColl==2) {
1268     betaPerDOFoptions[0] = 9.1e-5;
1269     betaPerDOFoptions[1] = 2.25e-4;
1270     betaPerDOFoptions[2] = 4.55e-4;
1271     betaPerDOFoptions[3] = 9.1e-4;
1272     betaPerDOFoptions[4] = 2.25e-3;
1273     betaPerDOFoptions[5] = 4.55e-3;
1274   }
1275   if(betaColl==3) {
1276     betaPerDOFoptions[0] = 9.1e-3;
1277     betaPerDOFoptions[1] = 2.25e-2;
1278     betaPerDOFoptions[2] = 4.55e-2;
1279     betaPerDOFoptions[3] = 9.1e-2;
1280     betaPerDOFoptions[4] = 2.25e-1;
1281     betaPerDOFoptions[5] = 4.55e-1;
1282   }
1283   if(betaColl==4) {
1284     betaPerDOFoptions[0] = 9.1e-1;
1285     betaPerDOFoptions[1] = 2.25;
1286     betaPerDOFoptions[2] = 4.55;
1287     betaPerDOFoptions[3] = 9.1;
1288     betaPerDOFoptions[4] = 22.5;
1289     betaPerDOFoptions[5] = 45.5;
1290   }
1291
1292   return betaPerDOFoptions[betaOpt];
1293
1294 }
1295
1296 //--------------------------------------------------------------------------------------------------------------------------------------------------
1297 Double_t AliAnaChargedJetResponseMaker::InterpolateFast(TGraph *gr, Double_t x) {
1298
1299   Double_t ipmax = gr->GetN()-1.;
1300   Double_t x2,y2,xold,yold;
1301
1302   Double_t xmin,ymin,xmax, ymax;
1303   gr->GetPoint(0,xmin,ymin);
1304   gr->GetPoint(gr->GetN()-1,xmax,ymax);
1305   if(x<xmin || x>xmax) return 0;
1306
1307   Double_t ip = ipmax/2.;
1308   Double_t ipold = 0.;
1309   gr->GetPoint((int)(ip),x2,y2);
1310
1311   Int_t i = 0;
1312
1313   if(x2>x) {
1314     while(x2>x) {
1315       if(i==0) ipold=0.;
1316       ip -= (ip)/2.;
1317       gr->GetPoint((int)(ip),x2,y2);
1318       if(x2>x){}
1319       else ipold = ip;
1320       i++;
1321       //      cout << "ipold: " << ipold << "\tip: " << ip << "\tx2: " << x2 << "\ty2: " << y2 << endl;
1322     }
1323   }
1324   else if(x2<x) {
1325     while(x2<x) {
1326       if(i==0) ipold=ipmax;
1327       ip += (ipold-ip)/2.;
1328       gr->GetPoint((int)(ip),x2,y2);
1329       if(x2>x) ipold = ip;
1330       else {}
1331       i++;
1332       //      cout << "ipold: " << ipold << "\tip: " << ip << "\tx2: " << x2 << "\ty2: " << y2 << endl;
1333     }
1334   }
1335
1336   int ip2 = 0;
1337   if(x2>x) {
1338     ip2 = (int)(ip)-1;
1339     gr->GetPoint(ip2,x2,y2);
1340     while(x2>x) {
1341       ip2--;
1342       gr->GetPoint(ip2,x2,y2);
1343     }
1344     gr->GetPoint(ip2+1,xold,yold);
1345   }
1346   else {
1347     ip2 = (int)(ip)+1;
1348     gr->GetPoint(ip2,x2,y2);
1349     while(x2<x) {
1350       ip2++;
1351       gr->GetPoint(ip2,x2,y2);
1352     }
1353     gr->GetPoint(ip2-1,xold,yold);
1354
1355   }
1356   //  cout << "For x=" << x << " interpolate between: " << xold << " and " << x2 << endl;
1357   return ((x-xold)*y2 + (x2-x)*yold) / (x2-xold);
1358
1359 }
1360
1361 //--------------------------------------------------------------------------------------------------------------------------------------------------
1362 Double_t AliAnaChargedJetResponseMaker::InterpolateFast(TH1 *h, Double_t x) {
1363
1364   // Double_t ipmax = gr->GetN()-1.;
1365   Double_t ipmax = (double)h->GetNbinsX();
1366   Double_t x2,y2,xold,yold;
1367
1368   Double_t xmin = h->GetXaxis()->GetXmin();
1369   Double_t xmax = h->GetXaxis()->GetXmax();
1370   if(x<xmin || x>xmax) return 0;
1371
1372   Double_t ip = ipmax/2.;
1373   Double_t ipold = 0.;
1374
1375   x2 = h->GetBinCenter((int)ip);
1376   y2 = h->GetBinContent((int)ip);
1377
1378   Int_t i = 0;
1379
1380   if(x2>x) {
1381     while(x2>x) {
1382       if(i==0) ipold=0.;
1383       ip -= (ip)/2.;
1384       x2 = h->GetBinCenter((int)ip);
1385       y2 = h->GetBinContent((int)ip);
1386       if(x2>x) {}
1387       else ipold = ip;
1388       i++;
1389       //      cout << "ipold: " << ipold << "\tip: " << ip << "\tx2: " << x2 << "\ty2: " << y2 << endl;
1390     }
1391   }
1392   else if(x2<x) {
1393     while(x2<x) {
1394       if(i==0) ipold=ipmax;
1395       ip += (ipold-ip)/2.;
1396       x2 = h->GetBinCenter((int)ip);
1397       y2 = h->GetBinContent((int)ip);
1398       if(x2>x) ipold = ip;
1399       else {}
1400       i++;
1401       //      cout << "ipold: " << ipold << "\tip: " << ip << "\tx2: " << x2 << "\ty2: " << y2 << endl;
1402     }
1403   }
1404
1405   int ip2 = 0;
1406   if(x2>x) {
1407     ip2 = (int)(ip)-1;
1408     x2 = h->GetBinCenter(ip2);
1409     y2 = h->GetBinContent(ip2);
1410     while(x2>x) {
1411       ip2--;
1412       x2 = h->GetBinCenter(ip2);
1413       y2 = h->GetBinContent(ip2);
1414     }
1415     xold = h->GetBinCenter(ip2+1);
1416     yold = h->GetBinContent(ip2+1);
1417   }
1418   else {
1419     ip2 = (int)(ip)+1;
1420     x2 = h->GetBinCenter(ip2);
1421     y2 = h->GetBinContent(ip2);
1422     while(x2<x) {
1423       ip2++;
1424       x2 = h->GetBinCenter(ip2);
1425       y2 = h->GetBinContent(ip2);
1426     }
1427     xold = h->GetBinCenter(ip2-1);
1428     yold = h->GetBinContent(ip2-1);
1429
1430   }
1431   //  cout << "For x=" << x << " interpolate between: " << xold << " and " << x2 << endl;
1432   return ((x-xold)*y2 + (x2-x)*yold) / (x2-xold);
1433
1434 }
1435
1436