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