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 ClassImp(AliAnaChargedJetResponseMaker)
  16
  17 AliAnaChargedJetResponseMaker::AliAnaChargedJetResponseMaker():
  18   fDebug(kFALSE),
  19   fResolutionType(kParam),
  20   fDeltaPt(0x0),
  21   fhDeltaPt(0x0),
  22   fDimensions(1),
  23   fDimRec(0),
  24   fDimGen(1),
  25   fPtMin(-999),
  26   fPtMax(-999),
  27   fNbins(0),
  28   fBinArrayPtRec(0x0),
  29   fPtMeasured(0x0),
  30   fEffFlat(0),
  31   fEfficiency(0x0),
  32   fEfficiencyFine(0x0),
  33   fResponseMatrix(0x0),
  34   fResponseMatrixFine(0x0),
  35   fPtMinUnfolded(0.),
  36   fPtMaxUnfolded(0.),
  37   fPtMaxUnfoldedHigh(-1.),
  38   fBinWidthFactorUnfolded(2),
  39   fSkipBinsUnfolded(0),
  40   fExtraBinsUnfolded(5),
  41   fbVariableBinning(kFALSE),
  42   fPtMaxUnfVarBinning(0),
  43   f1MergeFunction(0x0),
  44   fFineFrac(10),
  45   fbCalcErrors(kFALSE)
  46 {;}
  47
  48 AliAnaChargedJetResponseMaker::~AliAnaChargedJetResponseMaker() {
  49   //destructor
  50
  51 }
  52
  53 //--------------------------------------------------------------------------------------------------------------------------------------------------
  54 void AliAnaChargedJetResponseMaker::SetMeasuredSpectrum(TH1D *hPtMeasured)
  55 {
  56   //
  57   // Set measured spectrum in THnSparse format
  58   // This defines the minimum and maximum pT on the reconstructed axis of the response matrix
  59   //
  60   if(fDebug) printf(">>AliAnaChargedJetResponseMaker::SetMeasuredSpectrum \n");
  61
  62   fNbins = hPtMeasured->GetXaxis()->GetNbins();
  63   fPtMin = hPtMeasured->GetXaxis()->GetXmin();
  64   fPtMax = hPtMeasured->GetXaxis()->GetXmax();
  65
  66   if(fDebug) printf("fNbins: %d  fPtMin: %f  fPtMax: %f \n",fNbins,fPtMin,fPtMax);
  67
  68   if(fBinArrayPtRec) delete fBinArrayPtRec;
  69   fBinArrayPtRec = new Double_t[fNbins+1];
  70   for(int j = 0; j<fNbins; j++) {
  71     fBinArrayPtRec[j] = hPtMeasured->GetXaxis()->GetBinLowEdge(j+1);
  72   }
  73   fBinArrayPtRec[fNbins] = hPtMeasured->GetXaxis()->GetBinUpEdge(fNbins);
  74
  75
  76   Int_t nbins[fDimensions];
  77   Double_t xmin[fDimensions];
  78   Double_t xmax[fDimensions];
  79   for(int dim = 0; dim<fDimensions; dim++) {
  80     nbins[dim] = fNbins;
  81     xmin[dim]  = fPtMin;
  82     xmax[dim]  = fPtMax;
  83   }
  84
  85   if(fPtMeasured) delete fPtMeasured;
  86   fPtMeasured = new THnSparseD("fPtMeasured","Measured pT spectrum; p_{T}^{rec}",fDimensions,nbins,xmin,xmax);
  87   fPtMeasured->Sumw2();
  88   fPtMeasured->GetAxis(0)->SetTitle("p_{T}^{rec}");
  89   fPtMeasured->SetBinEdges(0,fBinArrayPtRec);
  90
  91   //Fill
  92   if(fDebug) printf("fill measured THnSparse\n");
  93   if(fNbins!=hPtMeasured->GetNbinsX())
  94     printf("WARNING: nbins not correct \t %d vs %d !!!\n",fNbins,hPtMeasured->GetNbinsX());
  95
  96   int bin[1] = {0};
  97   bin[0] = 0;
  98   for(int i = hPtMeasured->FindBin(fPtMin); i<hPtMeasured->FindBin(fPtMax); i++) {
  99     double pT[1];
 100     pT[0]= hPtMeasured->GetBinCenter(i);
 101     fPtMeasured->SetBinContent(bin,(Double_t)hPtMeasured->GetBinContent(i));
 102     fPtMeasured->SetBinError(bin,(Double_t)hPtMeasured->GetBinError(i));
 103     bin[0]++;
 104   }
 105
 106   if(fDebug) printf("fPtMeasured->GetNbins(): %d \n",(int)fPtMeasured->GetNbins());
 107
 108 }
 109
 110 //--------------------------------------------------------------------------------------------------------------------------------------------------
 111 void AliAnaChargedJetResponseMaker::SetFlatEfficiency(Double_t eff) {
 112
 113   fEffFlat = eff;
 114
 115   Int_t nbins[fDimensions];
 116   Double_t xmin[fDimensions];
 117   Double_t xmax[fDimensions];
 118   for(int dim = 0; dim<fDimensions; dim++) {
 119     nbins[dim] = fNbins;
 120     xmin[dim] = fPtMin;
 121     xmax[dim] = fPtMax;
 122   }
 123
 124   if(fEfficiency) delete fEfficiency;
 125   fEfficiency = new THnSparseD("fEfficiency","Efficiency - no resolution effects; p_{T}^{gen}",fDimensions,nbins,xmin,xmax);
 126   fEfficiency->Sumw2();
 127   fEfficiency->GetAxis(0)->SetTitle("p_{T}^{gen}");
 128   fEfficiency->SetBinEdges(0,fBinArrayPtRec);
 129
 130   for(int i=0; i<fNbins; i++) {
 131     int bin[1];
 132     bin[0] = i;
 133     fEfficiency->SetBinContent(bin,fEffFlat);
 134     fEfficiency->SetBinError(bin,0);
 135   }
 136
 137 }
 138
 139 //--------------------------------------------------------------------------------------------------------------------------------------------------
 140 void AliAnaChargedJetResponseMaker::SetEfficiency(TGraphErrors *grEff)
 141 {
 142   Int_t nbins[fDimensions];
 143   Double_t xmin[fDimensions];
 144   Double_t xmax[fDimensions];
 145   for(int dim = 0; dim<fDimensions; dim++) {
 146     nbins[dim] = fNbins;
 147     xmin[dim] = fPtMin;
 148     xmax[dim] = fPtMax;
 149   }
 150
 151   if(fEfficiency) delete fEfficiency;
 152   fEfficiency = new THnSparseD("fEfficiency","Efficiency - no resolution effects; p_{T}^{gen}",fDimensions,nbins,xmin,xmax);
 153   fEfficiency->Sumw2();
 154   fEfficiency->GetAxis(0)->SetTitle("p_{T}^{gen}");
 155   fEfficiency->SetBinEdges(0,fBinArrayPtRec);
 156
 157   double pT[1];
 158   double yield = 0.;
 159   double error = 0.;
 160   double dummy = 0.;
 161   int nbinsgr = grEff->GetN();
 162
 163   for(int i=0; i<nbinsgr; i++) {
 164     grEff->GetPoint(i,dummy,yield);
 165     pT[0] = dummy;
 166     error = grEff->GetErrorY(i);
 167
 168     fEfficiency->Fill(pT,yield);
 169     fEfficiency->SetBinError(i,error);
 170
 171   }
 172
 173 }
 174
 175 //--------------------------------------------------------------------------------------------------------------------------------------------------
 176 void AliAnaChargedJetResponseMaker::MakeResponseMatrixJetsFineMerged(Int_t skipBins, Int_t binWidthFactor, Int_t extraBins, Bool_t bVariableBinning, Double_t ptmax)
 177 {
 178   //
 179   // Make jet response matrix
 180   //
 181
 182   if(!fPtMeasured) return;
 183   if(fResponseMatrix)     { delete fResponseMatrix; }
 184   if(fResponseMatrixFine) { delete fResponseMatrixFine; }
 185
 186   SetSkipBinsUnfolded(skipBins);
 187   SetBinWidthFactorUnfolded(binWidthFactor);
 188   SetExtraBinsUnfolded(extraBins);
 189   SetVariableBinning(bVariableBinning,ptmax);
 190
 191   InitializeResponseMatrix();
 192   InitializeEfficiency();
 193
 194   InitializeResponseMatrixFine();
 195   InitializeEfficiencyFine();
 196
 197   FillResponseMatrixFineAndMerge();
 198
 199 }
 200
 201 //--------------------------------------------------------------------------------------------------------------------------------------------------
 202 void AliAnaChargedJetResponseMaker::InitializeResponseMatrix() {
 203   //
 204   //Define bin width of RM to be used for unfolding
 205   //
 206
 207   Int_t nbins[fDimensions*2];
 208   nbins[fDimRec] = fNbins;
 209
 210   double binWidthMeas = (double)((fPtMax-fPtMin)/fNbins);
 211   double binWidthUnf  = (double)fBinWidthFactorUnfolded*binWidthMeas;
 212   double binWidthUnfLowPt = -1.;
 213   if(fbVariableBinning)
 214     binWidthUnfLowPt = binWidthUnf*0.5;
 215
 216   if(fExtraBinsUnfolded>0) {
 217     fPtMaxUnfolded = fPtMax+(double)(fExtraBinsUnfolded)*binWidthUnf;
 218     nbins[fDimGen]+=fExtraBinsUnfolded;
 219   }
 220
 221   printf("fPtMinMeas: %f  fPtMaxMeas: %f\n",fPtMin,fPtMax);
 222   printf("binWidthMeas: %f  binWidthUnf: %f   fBinWidthFactorUnfolded: %d\n",binWidthMeas,binWidthUnf,fBinWidthFactorUnfolded);
 223   printf("binWidthUnfLowPt: %f\n",binWidthUnfLowPt);
 224
 225   int tmp = round((fPtMaxUnfolded/binWidthUnf)); //nbins which fit between 0 and fPtMaxUnfolded
 226   tmp = tmp - fSkipBinsUnfolded;
 227   fPtMinUnfolded = fPtMaxUnfolded-(double)(tmp)*binWidthUnf; //set ptmin unfolded
 228   fPtMaxUnfolded = fPtMinUnfolded+(double)(tmp)*binWidthUnf; //set ptmax unfolded
 229   nbins[fDimGen] = (int)((fPtMaxUnfolded-fPtMinUnfolded)/binWidthUnf); //adjust nbins to bin width
 230
 231   if(fbVariableBinning) {
 232     tmp = (int)(fPtMaxUnfVarBinning/binWidthUnfLowPt);
 233     tmp = tmp - fSkipBinsUnfolded;
 234     fPtMinUnfolded = fPtMaxUnfVarBinning-(double)(tmp)*binWidthUnfLowPt;
 235     //Redefine also number of bins on generated axis in case of variable binning
 236     nbins[fDimGen] = (int)((fPtMaxUnfVarBinning-fPtMinUnfolded)/binWidthUnfLowPt)+(int)((fPtMaxUnfolded-fPtMaxUnfVarBinning)/binWidthUnf);
 237   }
 238
 239   Double_t binWidth[2];
 240   binWidth[fDimRec] = (double)((fPtMax-fPtMin)/nbins[fDimRec]);
 241   binWidth[fDimGen] = (double)((fPtMaxUnfolded-fPtMinUnfolded)/nbins[fDimGen]);
 242
 243   Double_t xmin[fDimensions*2];
 244   Double_t xmax[fDimensions*2];
 245   xmin[fDimRec] = fPtMin;
 246   xmax[fDimRec] = fPtMax;
 247   xmin[fDimGen] = fPtMinUnfolded;
 248   xmax[fDimGen] = fPtMaxUnfolded;
 249
 250   printf("xmin[fDimRec]: %f  xmin[fDimGen]: %f \n",xmin[fDimRec],xmin[fDimGen]);
 251   printf("xmax[fDimRec]: %f  xmax[fDimGen]: %f \n",xmax[fDimRec],xmax[fDimGen]);
 252   printf("nbins[fDimRec]: %d  nbins[fDimGen]: %d \n",nbins[fDimRec],nbins[fDimGen]);
 253   if(!fbVariableBinning) printf("binWidth[fDimRec]: %f  binWidth[fDimGen]: %f \n",binWidth[fDimRec],binWidth[fDimGen]);
 254
 255   Double_t binArrayPt0[nbins[fDimRec]+1];
 256   Double_t binArrayPt1[nbins[fDimGen]+1];
 257   Double_t xmaxGen = TMath::Max(xmax[fDimGen],fPtMaxUnfoldedHigh);
 258
 259   //Define bin limits reconstructed/measured axis
 260   for(int i=0; i<nbins[fDimRec]; i++) {
 261     binArrayPt0[i] = xmin[fDimRec]+(double)i*binWidth[fDimRec];
 262   }
 263   binArrayPt0[nbins[fDimRec]]= xmax[fDimRec];
 264
 265   //Define bin limits generated/unfolded axis
 266   for(int i=0; i<nbins[fDimGen]; i++) {
 267     if(!fbVariableBinning) binArrayPt1[i] = xmin[fDimGen]+(double)i*binWidth[fDimGen];
 268     if(fbVariableBinning) {
 269       double test = xmin[fDimGen]+(double)i*binWidthUnfLowPt;
 270       if(test<=fPtMaxUnfVarBinning) binArrayPt1[i] = test;
 271       else binArrayPt1[i] = binArrayPt1[i-1]+binWidthUnf;
 272     }
 273     //printf("RM. i = %d \t binArrayPt[i] = %f \n",i,binArrayPt1[i]);
 274   }
 275   binArrayPt1[nbins[fDimGen]]= xmaxGen;
 276
 277
 278   // Response matrix : dimensions must be 2N = 2 x (number of variables)
 279   // dimensions 0 ->  N-1 must be filled with reconstructed values
 280   // dimensions N -> 2N-1 must be filled with generated values
 281   fResponseMatrix = new THnSparseD("fResponseMatrix","Response Matrix pTMC vs pTrec",fDimensions*2,nbins,xmin,xmax);
 282   fResponseMatrix->Sumw2();
 283   fResponseMatrix->GetAxis(fDimRec)->SetTitle("p_{T}^{rec}");
 284   fResponseMatrix->GetAxis(fDimGen)->SetTitle("p_{T}^{gen}");
 285
 286   fResponseMatrix->SetBinEdges(fDimRec,binArrayPt0);
 287   fResponseMatrix->SetBinEdges(fDimGen,binArrayPt1);
 288
 289
 290 }
 291
 292 //--------------------------------------------------------------------------------------------------------------------------------------------------
 293 void AliAnaChargedJetResponseMaker::InitializeEfficiency() {
 294   //
 295   // Define dimensions of efficiency THnSparse
 296   //
 297
 298   if(!fResponseMatrix) {
 299     printf("AliAnaChargedJetResponseMaker::InitializeEfficiency()\n");
 300     printf("Can not define dimensions efficiency without fResponseMatrix dimensions. Aborting \n");
 301     return;
 302   }
 303
 304   TAxis *genAxis = fResponseMatrix->GetAxis(fDimGen);
 305
 306   Int_t nbinsEff[fDimensions];
 307   Double_t xminEff[fDimensions];
 308   Double_t xmaxEff[fDimensions];
 309
 310   for(int dim = 0; dim<fDimensions; dim++) {
 311     nbinsEff[dim] = genAxis->GetNbins();
 312     xminEff[dim]  = genAxis->GetXmin();
 313     xmaxEff[dim]  = genAxis->GetXmax();
 314   }
 315
 316   if(fEfficiency) delete fEfficiency;
 317   fEfficiency = new THnSparseD("fEfficiency","Efficiency - no resolution effects; p_{T}^{gen}",fDimensions,nbinsEff,xminEff,xmaxEff);
 318   fEfficiency->Sumw2();
 319   fEfficiency->GetAxis(0)->SetTitle("p_{T}^{gen}");
 320
 321   const Double_t *binArrayPt = genAxis->GetXbins()->GetArray();
 322   fEfficiency->SetBinEdges(0,binArrayPt);
 323
 324 }
 325
 326 //--------------------------------------------------------------------------------------------------------------------------------------------------
 327 void AliAnaChargedJetResponseMaker::InitializeResponseMatrixFine() {
 328   //
 329   //Initialize fine response matrix
 330   //
 331
 332   Int_t nbinsFine[fDimensions*2];
 333   Double_t xminFine[fDimensions*2];
 334   Double_t xmaxFine[fDimensions*2];
 335   Double_t pTarrayFine[fDimensions*2];
 336
 337   for(int i=0; i<fDimensions*2; i++) {
 338
 339     if(i==0) nbinsFine[i] = fResponseMatrix->GetAxis(fDimRec)->GetNbins()*fFineFrac; // dimension 0 is reconstructed value
 340     else     nbinsFine[i] = fResponseMatrix->GetAxis(fDimRec)->GetNbins()*fFineFrac; // dimension 1 is generated value
 341
 342     xminFine[i] = TMath::Min(fPtMin,0.);
 343     xmaxFine[i] = fResponseMatrix->GetAxis(fDimGen)->GetXmax()+40.;
 344
 345     pTarrayFine[i] = 0.;
 346   }
 347
 348   Double_t binWidth[2];
 349   binWidth[fDimRec] = fResponseMatrix->GetAxis(fDimRec)->GetBinWidth(1);
 350
 351   Double_t binWidthFine[2];
 352   binWidthFine[fDimRec] = binWidth[fDimRec]/((double)fFineFrac);
 353   binWidthFine[fDimGen] = binWidthFine[fDimRec]*(double)fBinWidthFactorUnfolded;
 354
 355   nbinsFine[fDimRec] = (int)((xmaxFine[fDimRec]-xminFine[fDimRec])/binWidthFine[fDimRec]); //adjust nbins to bin width
 356   nbinsFine[fDimGen] = (int)((xmaxFine[fDimGen]-xminFine[fDimGen])/binWidthFine[fDimGen]); //adjust nbins to bin width
 357
 358   printf("xminFine[fDimRec]: %f  xminFine[fDimGen]: %f \n",xminFine[fDimRec],xminFine[fDimGen]);
 359   printf("xmaxFine[fDimRec]: %f  xmaxFine[fDimGen]: %f \n",xmaxFine[fDimRec],xmaxFine[fDimGen]);
 360   printf("nbinsFine[fDimRec]: %d  nbinsFine[fDimGen]: %d \n",nbinsFine[fDimRec],nbinsFine[fDimGen]);
 361   printf("binWidthFine[fDimRec]: %f  binWidthFine[fDimGen]: %f \n",binWidthFine[fDimRec],binWidthFine[fDimGen]);
 362
 363
 364   Double_t binArrayPt0Fine[nbinsFine[fDimRec]+1];
 365   Double_t binArrayPt1Fine[nbinsFine[fDimGen]+1];
 366
 367   for(int i=0; i<nbinsFine[fDimRec]; i++) {
 368     binArrayPt0Fine[i] = xminFine[fDimRec]+(double)i*binWidthFine[fDimRec];
 369     //    printf("RM. i = %d \t binArrayPtFine[i] = %f \n",i,binArrayPt0Fine[i]);
 370   }
 371   binArrayPt0Fine[nbinsFine[fDimRec]]= xmaxFine[fDimRec];
 372
 373   for(int i=0; i<nbinsFine[fDimGen]; i++) {
 374     binArrayPt1Fine[i] = xminFine[fDimGen]+(double)i*binWidthFine[fDimGen];
 375     //    printf("RM. i = %d \t binArrayPtFine[i] = %f \n",i,binArrayPt1Fine[i]);
 376   }
 377   binArrayPt1Fine[nbinsFine[fDimGen]]= xmaxFine[fDimGen];
 378
 379   // Response matrix : dimensions must be 2N = 2 x (number of variables)
 380   // dimensions 0 ->  N-1 must be filled with reconstructed values
 381   // dimensions N -> 2N-1 must be filled with generated values
 382   fResponseMatrixFine = new THnSparseD("fResponseMatrixFine","Response Matrix pTMC vs pTrec",fDimensions*2,nbinsFine,xminFine,xmaxFine);
 383   fResponseMatrixFine->Sumw2();
 384   fResponseMatrixFine->GetAxis(fDimRec)->SetTitle("p_{T}^{rec}");
 385   fResponseMatrixFine->GetAxis(fDimGen)->SetTitle("p_{T}^{gen}");
 386
 387   fResponseMatrixFine->SetBinEdges(fDimRec,binArrayPt0Fine);
 388   fResponseMatrixFine->SetBinEdges(fDimGen,binArrayPt1Fine);
 389
 390 }
 391
 392 //--------------------------------------------------------------------------------------------------------------------------------------------------
 393 void AliAnaChargedJetResponseMaker::InitializeEfficiencyFine() {
 394   //
 395   // Define dimensions of efficiency THnSparse
 396   //
 397
 398   if(!fResponseMatrixFine) {
 399     printf("AliAnaChargedJetResponseMaker::InitializeEfficiencyFine()\n");
 400     printf("Can not define dimensions efficiency without fResponseMatrixFine dimensions. Aborting \n");
 401     return;
 402   }
 403
 404   TAxis *genAxis = fResponseMatrixFine->GetAxis(fDimGen);
 405
 406   Int_t nbinsEff[fDimensions];
 407   Double_t xminEff[fDimensions];
 408   Double_t xmaxEff[fDimensions];
 409
 410   for(int dim = 0; dim<fDimensions; dim++) {
 411     nbinsEff[dim] = genAxis->GetNbins();
 412     xminEff[dim]  = genAxis->GetXmin();
 413     xmaxEff[dim]  = genAxis->GetXmax();
 414   }
 415
 416   if(fEfficiencyFine) delete fEfficiencyFine;
 417   fEfficiencyFine = new THnSparseD("fEfficiencyFine","EfficiencyFine - no resolution effects; p_{T}^{gen}",fDimensions,nbinsEff,xminEff,xmaxEff);
 418   fEfficiencyFine->Sumw2();
 419   fEfficiencyFine->GetAxis(0)->SetTitle("p_{T}^{gen}");
 420
 421   const Double_t *binArrayPt = genAxis->GetXbins()->GetArray();
 422   fEfficiencyFine->SetBinEdges(0,binArrayPt);
 423
 424 }
 425
 426 //--------------------------------------------------------------------------------------------------------------------------------------------------
 427 void AliAnaChargedJetResponseMaker::FillResponseMatrixFineAndMerge() {
 428   //
 429   // Fill fine response matrix
 430   //
 431
 432   TAxis *genAxis = fResponseMatrixFine->GetAxis(fDimGen);
 433   TAxis *recAxis = fResponseMatrixFine->GetAxis(fDimRec);
 434
 435   Int_t nbinsFine[fDimensions*2];
 436   Double_t xminFine[fDimensions*2];
 437   Double_t xmaxFine[fDimensions*2];
 438   Double_t pTarrayFine[fDimensions*2];
 439
 440   nbinsFine[fDimGen] = genAxis->GetNbins();
 441   nbinsFine[fDimRec] = recAxis->GetNbins();
 442   xminFine[fDimGen]  = genAxis->GetXmin();
 443   xminFine[fDimRec]  = recAxis->GetXmin();
 444   xmaxFine[fDimGen]  = genAxis->GetXmax();
 445   xmaxFine[fDimRec]  = recAxis->GetXmax();
 446
 447   double sumyield = 0.;
 448   double sumyielderror2 = 0.;
 449
 450   int ipt[2]    = {0.,0.};
 451   int iptMerged[2]    = {0.,0.};
 452   int ierror[2] = {0.,0.};
 453
 454   Int_t check = 0;
 455   Double_t pTgen, pTrec;
 456   Double_t yield = 0.;
 457   Double_t error = 0.;
 458
 459   const int nng = fResponseMatrix->GetAxis(fDimGen)->GetNbins();
 460   const int nnr = fResponseMatrix->GetAxis(fDimRec)->GetNbins();
 461   Double_t errorArray[nng][nnr];
 462   for(int iig =0; iig<nng; iig++) {
 463     for(int iir =0; iir<nnr; iir++) {
 464       errorArray[iig][iir] = 0.;
 465     }
 466   }
 467
 468   for(int iy=1; iy<=nbinsFine[fDimGen]; iy++) { //gen
 469     pTgen = fResponseMatrixFine->GetAxis(fDimGen)->GetBinCenter(iy);
 470     pTarrayFine[fDimGen] = pTgen;
 471     ierror[fDimGen]=iy;
 472     sumyield = 0.;
 473     check = 0;
 474
 475     for(int ix=1; ix<=nbinsFine[fDimRec]; ix++) { //rec
 476       pTrec = fResponseMatrixFine->GetAxis(fDimRec)->GetBinCenter(ix);
 477       Double_t width = fResponseMatrixFine->GetAxis(fDimRec)->GetBinWidth(ix);
 478       if(fResolutionType==kParam) {
 479         yield = fDeltaPt->Eval(pTrec-pTgen);
 480         error = 0.;
 481       }
 482       else if(fResolutionType==kResiduals) {
 483         yield = fhDeltaPt->Interpolate(pTrec-pTgen);
 484         error = 0.;
 485       }
 486       else if(fResolutionType==kResidualsErr) {
 487         Double_t deltaPt = pTrec-pTgen;
 488         int bin = fhDeltaPt->FindBin(deltaPt);
 489         yield = fhDeltaPt->GetBinContent(bin);
 490         error = fhDeltaPt->GetBinError(bin);
 491       }
 492       yield=yield*width;
 493       error=error*width;
 494       //avoid trouble with empty bins in the high pT tail of deltaPt distribution
 495       if(check==0 && yield>0. && pTrec>pTgen) check=1;
 496       if(check==1 && yield==0.) ix=nbinsFine[fDimRec];
 497
 498       sumyield+=yield;
 499       sumyielderror2 += error*error;
 500
 501       pTarrayFine[fDimRec] = pTrec;
 502       ierror[fDimRec]  = ix;
 503       fResponseMatrixFine->Fill(pTarrayFine,yield);
 504       if(fbCalcErrors) fResponseMatrixFine->SetBinError(ierror,error);
 505
 506     }// ix (dimRec) loop
 507
 508     //Normalize to total number of counts =1
 509     if(sumyield>1) {
 510       ipt[fDimGen]=iy;
 511       for(int ix=1; ix<=nbinsFine[fDimRec]; ix++) {
 512         ipt[fDimRec]=ix;
 513         fResponseMatrixFine->SetBinContent(ipt,fResponseMatrixFine->GetBinContent(ipt)/sumyield);
 514         if(fResolutionType==kResidualsErr && fbCalcErrors) {
 515           Double_t A = 1./sumyield*fResponseMatrix->GetBinError(ipt);
 516           Double_t B = -1.*fResponseMatrix->GetBinContent(ipt)/sumyield/sumyield*TMath::Sqrt(sumyielderror2);
 517           Double_t tmp2 = A*A + B*B;
 518           fResponseMatrix->SetBinError(ipt,TMath::Sqrt(tmp2));
 519         }
 520
 521       }
 522     }
 523
 524     int bin[1];
 525     bin[0] = iy;
 526     fEfficiencyFine->SetBinContent(bin,sumyield);
 527     if(fbCalcErrors) fEfficiencyFine->SetBinError(bin,TMath::Sqrt(sumyielderror2));
 528
 529     //fill merged response matrix
 530
 531     //find bin in fine RM correspoinding to mimimum/maximum bin of merged RM on rec axis
 532     int ixMin = fResponseMatrixFine->GetAxis(fDimRec)->FindBin(fResponseMatrix->GetAxis(fDimRec)->GetXmin());
 533     int ixMax = fResponseMatrixFine->GetAxis(fDimRec)->FindBin(fResponseMatrix->GetAxis(fDimRec)->GetXmax());
 534
 535     if(fResponseMatrixFine->GetAxis(fDimGen)->GetBinLowEdge(iy) >= fResponseMatrix->GetAxis(fDimGen)->GetXmin()) {
 536       ipt[fDimGen]=iy;
 537       iptMerged[fDimGen]=fResponseMatrix->GetAxis(fDimGen)->FindBin(pTgen);
 538
 539       Double_t weight = 1.;
 540       if(f1MergeFunction) {
 541         Double_t loEdge = fResponseMatrix->GetAxis(fDimGen)->GetBinLowEdge(iptMerged[fDimGen]);
 542         Double_t upEdge = fResponseMatrix->GetAxis(fDimGen)->GetBinUpEdge(iptMerged[fDimGen]);
 543         Float_t powInteg = f1MergeFunction->Integral(loEdge,upEdge);
 544         //printf("loEdge = %f  upEdge = %f  powInteg = %f\n",loEdge,upEdge,powInteg);
 545         if(powInteg>0.)
 546           weight = f1MergeFunction->Integral(fResponseMatrixFine->GetAxis(fDimGen)->GetBinLowEdge(iy),fResponseMatrixFine->GetAxis(fDimGen)->GetBinUpEdge(iy))/powInteg;
 547         //      printf("weight: %f \n", weight );
 548       } else {
 549         weight = 1./((double)fFineFrac);
 550       }
 551
 552       for(int ix=ixMin; ix<=ixMax; ix++) {                    //loop reconstructed axis
 553         pTrec = fResponseMatrixFine->GetAxis(fDimRec)->GetBinCenter(ix);
 554         ipt[fDimRec]=ix;
 555         iptMerged[fDimRec]=fResponseMatrix->GetAxis(fDimRec)->FindBin(pTrec);
 556
 557         fResponseMatrix->AddBinContent(iptMerged,fResponseMatrixFine->GetBinContent(ipt)*weight);
 558         if(fbCalcErrors) errorArray[iptMerged[fDimGen]-1][iptMerged[fDimRec]-1] += fResponseMatrixFine->GetBinError(ipt)*fResponseMatrixFine->GetBinError(ipt)*weight*weight;
 559       }
 560
 561    }//loop gen axis fine response matrix
 562
 563   } // iy (dimGen) loop
 564
 565   //Fill Efficiency corresponding to merged response matrix
 566   // FillEfficiency(errorArray,fFineFrac);
 567
 568   for(int iy=1; iy<=fResponseMatrix->GetAxis(fDimGen)->GetNbins(); iy++) { //gen
 569     sumyield = 0.;
 570     ipt[fDimGen] = iy;
 571
 572     for(int ix=1; ix<=fResponseMatrix->GetAxis(fDimRec)->GetNbins(); ix++) { //rec
 573       ipt[fDimRec] = ix;
 574       sumyield += fResponseMatrix->GetBinContent(ipt);
 575
 576       if(fbCalcErrors) fResponseMatrix->SetBinError(ipt,TMath::Sqrt(errorArray[iy][ix]));
 577     }
 578     printf("RM: pTgen: %f \t sumyield: %f \n",fResponseMatrix->GetAxis(fDimGen)->GetBinCenter(iy),sumyield);
 579     int bin[1];
 580     bin[0] = iy;
 581     fEfficiency->SetBinContent(bin,sumyield);
 582     if(fbCalcErrors) fEfficiency->SetBinError(bin,0);
 583   }
 584
 585   if(fDebug) printf("fResponseMatrixFine->GetNbins(): %d \n",(int)fResponseMatrixFine->GetNbins());
 586   if(fDebug) printf("fResponseMatrix->GetNbins(): %d \n",(int)fResponseMatrix->GetNbins());
 587
 588   if(fDebug) printf("Done constructing response matrix\n");
 589
 590 }
 591
 592
 593 //--------------------------------------------------------------------------------------------------------------------------------------------------
 594 TH1D* AliAnaChargedJetResponseMaker::MultiplyResponseGenerated(TH1 *hGen, TH2 *hResponse,TH1 *hEfficiency,Bool_t bDrawSlices) {
 595
 596   //
 597   // Multiply hGen with response matrix to obtain refolded spectrum
 598   //
 599
 600   //For response
 601   //x-axis: generated
 602   //y-axis: reconstructed
 603   if(fDebug>0) cout << "nbins hGen: " << hGen->GetNbinsX() << "\t nbins hResponseGen: " << hResponse->GetXaxis()->GetNbins() << "\t nbins hResponseRec: " << hResponse->GetYaxis()->GetNbins()  << endl;
 604
 605   TH1D *hRec = hResponse->ProjectionY("hRec");
 606   hRec->Sumw2();
 607   hRec->Reset();
 608   hRec->SetTitle("hRec");
 609   hRec->SetName("hRec");
 610
 611   for(int irec=1; irec<=hRec->GetNbinsX(); irec++)
 612     hRec->SetBinContent(irec,0);
 613
 614   TH1D *hRecGenBin = 0x0;
 615   TCanvas *cSlices = 0x0;
 616   if(bDrawSlices) {
 617     cSlices = new TCanvas("cSlices","cSlices:Slices",400,400);
 618     gPad->SetLogy();
 619   }
 620
 621   Double_t yieldMC = 0.;
 622   Double_t yieldMCerror = 0.;
 623   Double_t sumYield = 0.;
 624   const Int_t nbinsRec = hRec->GetNbinsX();
 625   Double_t sumError2[nbinsRec+1];
 626   for(int i=0; i<=nbinsRec; i++) sumError2[i] = 0.;
 627   Double_t eff = 0.;
 628
 629   for(int igen=1; igen<=hGen->GetNbinsX(); igen++) {
 630     //get pTMC
 631     sumYield = 0.;
 632     if(fEffFlat>0.)
 633       eff = fEffFlat;
 634     else
 635       eff = hEfficiency->GetBinContent(igen);
 636     yieldMC = hGen->GetBinContent(igen)*eff;
 637     //    yieldMCerror = hGen->GetBinError(igen);
 638
 639     if(bDrawSlices) {
 640       hRecGenBin = hResponse->ProjectionY(Form("hRecGenBin%d",igen));
 641       hRecGenBin->Sumw2();
 642       hRecGenBin->Reset();
 643       hRecGenBin->SetTitle(Form("hRecGenBin%d",igen));
 644       hRecGenBin->SetName(Form("hRecGenBin%d",igen));
 645     }
 646
 647     for(int irec=1; irec<=hRec->GetNbinsX(); irec++) {
 648       hRec->AddBinContent(irec,yieldMC*hResponse->GetBinContent(igen,irec));
 649       sumYield+=hResponse->GetBinContent(igen,irec);
 650       Double_t A = yieldMC*hResponse->GetBinError(igen,irec);
 651       Double_t B = hResponse->GetBinContent(igen,irec)*yieldMCerror;
 652       Double_t tmp2 = A*A + B*B;
 653       sumError2[irec-1] += tmp2 ;
 654
 655       if(bDrawSlices)
 656         hRecGenBin->SetBinContent(irec,yieldMC*hResponse->GetBinContent(igen,irec));
 657
 658     }
 659     if(bDrawSlices) {
 660       cSlices->cd();
 661       hRecGenBin->SetLineColor(igen);
 662       if(igen==1) hRecGenBin->DrawCopy();
 663       else hRecGenBin->DrawCopy("same");
 664     }
 665
 666     if(hRecGenBin) delete hRecGenBin;
 667
 668     cout << "igen: " << igen << "\tpTMC: " << hGen->GetXaxis()->GetBinCenter(igen) << "\teff:" << eff << "\tsumYield: " << sumYield << endl;
 669   }
 670
 671   for(int i=0; i<=nbinsRec; i++) hRec->SetBinError(i+1,TMath::Sqrt(sumError2[i]));
 672
 673
 674   return hRec;
 675 }
 676
 677 //--------------------------------------------------------------------------------------------------------------------------------------------------
 678 TH1D* AliAnaChargedJetResponseMaker::MultiplyResponseGenerated(TF1 *fGen, TH2 *hResponse,TH1 *hEfficiency) {
 679
 680   //For response
 681   //x-axis: generated
 682   //y-axis: reconstructed
 683
 684   if(fDebug>0) printf(">>AliAnaChargedJetResponseMaker::MultiplyResponseGenerated(TF1 *fGen, TH2 *hResponse,TH1 *hEfficiency)");
 685
 686   TH1D *hRec = hResponse->ProjectionY("hRec");
 687   hRec->Sumw2();
 688   hRec->Reset();
 689   hRec->SetTitle("hRec");
 690   hRec->SetName("hRec");
 691
 692   //  TH1D *hRec = new TH1D("hRec","hRec",hResponse->GetNbinsY(),hResponse->GetYaxis()->GetXmin(),hResponse->GetYaxis()->GetXmax());
 693
 694   for(int irec=1; irec<=hRec->GetNbinsX(); irec++)
 695     hRec->SetBinContent(irec,0);
 696
 697   Double_t yieldMC = 0.;
 698   Double_t sumYield = 0.;
 699   Double_t eff = 0.;
 700   for(int igen=1; igen<=hResponse->GetNbinsX(); igen++) {
 701     //get pTMC
 702     sumYield = 0.;
 703     double pTMC = hResponse->GetXaxis()->GetBinCenter(igen);
 704     int binEff = hEfficiency->FindBin(pTMC);
 705     if(fEffFlat>0.)
 706       eff = fEffFlat;
 707     else
 708       eff = hEfficiency->GetBinContent(binEff);
 709     yieldMC = fGen->Eval(pTMC)*eff;
 710     for(int irec=1; irec<=hResponse->GetNbinsY(); irec++) {
 711       hRec->AddBinContent(irec,yieldMC*hResponse->GetBinContent(igen,irec));
 712       sumYield+=hResponse->GetBinContent(igen,irec);
 713     }
 714     cout << "igen: " << igen << "\tpTMC: " << pTMC << "\tsumYield: " << sumYield << endl;
 715   }
 716
 717   return hRec;
 718 }
 719
 720 //--------------------------------------------------------------------------------------------------------------------------------------------------
 721 Double_t AliAnaChargedJetResponseMaker::InterpolateFast(TGraph *gr, Double_t x) {
 722
 723   Double_t ipmax = gr->GetN()-1.;
 724   Double_t x2,y2,xold,yold;
 725
 726   Double_t xmin,ymin,xmax, ymax;
 727   gr->GetPoint(0,xmin,ymin);
 728   gr->GetPoint(gr->GetN()-1,xmax,ymax);
 729   if(x<xmin || x>xmax) return 0;
 730
 731   Double_t ip = ipmax/2.;
 732   Double_t ipold = 0.;
 733   gr->GetPoint((int)(ip),x2,y2);
 734
 735   Int_t i = 0;
 736
 737   if(x2>x) {
 738     while(x2>x) {
 739       if(i==0) ipold=0.;
 740       ip -= (ip)/2.;
 741       gr->GetPoint((int)(ip),x2,y2);
 742       if(x2>x){}
 743       else ipold = ip;
 744       i++;
 745       //      cout << "ipold: " << ipold << "\tip: " << ip << "\tx2: " << x2 << "\ty2: " << y2 << endl;
 746     }
 747   }
 748   else if(x2<x) {
 749     while(x2<x) {
 750       if(i==0) ipold=ipmax;
 751       ip += (ipold-ip)/2.;
 752       gr->GetPoint((int)(ip),x2,y2);
 753       if(x2>x) ipold = ip;
 754       else {}
 755       i++;
 756       //      cout << "ipold: " << ipold << "\tip: " << ip << "\tx2: " << x2 << "\ty2: " << y2 << endl;
 757     }
 758   }
 759
 760   int ip2 = 0;
 761   if(x2>x) {
 762     ip2 = (int)(ip)-1;
 763     gr->GetPoint(ip2,x2,y2);
 764     while(x2>x) {
 765       ip2--;
 766       gr->GetPoint(ip2,x2,y2);
 767     }
 768     gr->GetPoint(ip2+1,xold,yold);
 769   }
 770   else {
 771     ip2 = (int)(ip)+1;
 772     gr->GetPoint(ip2,x2,y2);
 773     while(x2<x) {
 774       ip2++;
 775       gr->GetPoint(ip2,x2,y2);
 776     }
 777     gr->GetPoint(ip2-1,xold,yold);
 778
 779   }
 780   //  cout << "For x=" << x << " interpolate between: " << xold << " and " << x2 << endl;
 781   return ((x-xold)*y2 + (x2-x)*yold) / (x2-xold);
 782
 783 }
 784
 785 //--------------------------------------------------------------------------------------------------------------------------------------------------
 786 Double_t AliAnaChargedJetResponseMaker::InterpolateFast(TH1 *h, Double_t x) {
 787
 788   // Double_t ipmax = gr->GetN()-1.;
 789   Double_t ipmax = (double)h->GetNbinsX();
 790   Double_t x2,y2,xold,yold;
 791
 792   Double_t xmin = h->GetXaxis()->GetXmin();
 793   Double_t xmax = h->GetXaxis()->GetXmax();
 794   if(x<xmin || x>xmax) return 0;
 795
 796   Double_t ip = ipmax/2.;
 797   Double_t ipold = 0.;
 798
 799   x2 = h->GetBinCenter((int)ip);
 800   y2 = h->GetBinContent((int)ip);
 801
 802   Int_t i = 0;
 803
 804   if(x2>x) {
 805     while(x2>x) {
 806       if(i==0) ipold=0.;
 807       ip -= (ip)/2.;
 808       x2 = h->GetBinCenter((int)ip);
 809       y2 = h->GetBinContent((int)ip);
 810       if(x2>x) {}
 811       else ipold = ip;
 812       i++;
 813       //      cout << "ipold: " << ipold << "\tip: " << ip << "\tx2: " << x2 << "\ty2: " << y2 << endl;
 814     }
 815   }
 816   else if(x2<x) {
 817     while(x2<x) {
 818       if(i==0) ipold=ipmax;
 819       ip += (ipold-ip)/2.;
 820       x2 = h->GetBinCenter((int)ip);
 821       y2 = h->GetBinContent((int)ip);
 822       if(x2>x) ipold = ip;
 823       else {}
 824       i++;
 825       //      cout << "ipold: " << ipold << "\tip: " << ip << "\tx2: " << x2 << "\ty2: " << y2 << endl;
 826     }
 827   }
 828
 829   int ip2 = 0;
 830   if(x2>x) {
 831     ip2 = (int)(ip)-1;
 832     x2 = h->GetBinCenter(ip2);
 833     y2 = h->GetBinContent(ip2);
 834     while(x2>x) {
 835       ip2--;
 836       x2 = h->GetBinCenter(ip2);
 837       y2 = h->GetBinContent(ip2);
 838     }
 839     xold = h->GetBinCenter(ip2+1);
 840     yold = h->GetBinContent(ip2+1);
 841   }
 842   else {
 843     ip2 = (int)(ip)+1;
 844     x2 = h->GetBinCenter(ip2);
 845     y2 = h->GetBinContent(ip2);
 846     while(x2<x) {
 847       ip2++;
 848       x2 = h->GetBinCenter(ip2);
 849       y2 = h->GetBinContent(ip2);
 850     }
 851     xold = h->GetBinCenter(ip2-1);
 852     yold = h->GetBinContent(ip2-1);
 853
 854   }
 855   //  cout << "For x=" << x << " interpolate between: " << xold << " and " << x2 << endl;
 856   return ((x-xold)*y2 + (x2-x)*yold) / (x2-xold);
 857
 858 }
 859
 860