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