PWGJE/AliFragmentationFunctionCorrections.cxx

   1 // *************************************************************************
   2 // *                                                                       *
   3 // *               corrections to Fragmentation Functions                  *
   4 // *                                                                       *
   5 // *************************************************************************
   6
   7
   8 /**************************************************************************
   9  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
  10  *                                                                        *
  11  * Author: The ALICE Off-line Project.                                    *
  12  * Contributors are mentioned in the code where appropriate.              *
  13  *                                                                        *
  14  * Permission to use, copy, modify and distribute this software and its   *
  15  * documentation strictly for non-commercial purposes is hereby granted   *
  16  * without fee, provided that the above copyright notice appears in all   *
  17  * copies and that both the copyright notice and this permission notice   *
  18  * appear in the supporting documentation. The authors make no claims     *
  19  * about the suitability of this software for any purpose. It is          *
  20  * provided "as is" without express or implied warranty.                  *
  21  **************************************************************************/
  22
  23 /* $Id: */
  24
  25 #include "TMath.h"
  26 #include "TH2F.h"
  27 #include "THnSparse.h"
  28 #include "TFile.h"
  29 #include "TDirectory.h"
  30 #include "AliCFUnfolding.h"
  31 #include "AliFragmentationFunctionCorrections.h"
  32
  33 ///////////////////////////////////////////////////////////////////////////////
  34 //                                                                           //
  35 //  correction class for ouput of AliAnalysisTaskFragmentationFunction       //
  36 //                                                                           //
  37 //  corrections for: reconstruction efficiency, momentum resolution,         //
  38 //                   secondaries, UE / HI background, fluctuations           //
  39 //                   back-correction on jet energy on dN/dxi                 //
  40 //                                                                           //
  41 //  read MC ouput and write out efficiency histos, response matrices etc.    //
  42 //  read measured distributions and apply efficiency, response matrices etc. //
  43 //                                                                           //
  44 //  contact: o.busch@gsi.de                                                  //
  45 //                                                                           //
  46 ///////////////////////////////////////////////////////////////////////////////
  47
  48
  49 ClassImp(AliFragmentationFunctionCorrections)
  50
  51 //________________________________________________________________________
  52 AliFragmentationFunctionCorrections::AliFragmentationFunctionCorrections()
  53    : TObject()
  54    ,fDebug(0)
  55    ,fNJetPtSlices(0)
  56    ,fJetPtSlices(0)
  57    ,fNJets(0)
  58    ,fNJetsBgr(0)
  59    ,fNHistoBinsSinglePt(0)
  60    ,fHistoBinsSinglePt(0)
  61    ,fNHistoBinsPt(0)
  62    ,fNHistoBinsZ(0)
  63    ,fNHistoBinsXi(0)
  64    ,fHistoBinsPt(0)
  65    ,fHistoBinsZ(0)
  66    ,fHistoBinsXi(0)
  67    ,fNCorrectionLevels(0)
  68    ,fCorrFF(0)
  69    ,fNCorrectionLevelsBgr(0)
  70    ,fCorrBgr(0)
  71    ,fNCorrectionLevelsSinglePt(0)
  72    ,fCorrSinglePt(0)
  73    ,fh1FFXiShift(0)
  74    ,fh1EffSinglePt(0)
  75    ,fh1EffPt(0)
  76    ,fh1EffZ(0)
  77    ,fh1EffXi(0)
  78    ,fh1EffBgrPt(0)
  79    ,fh1EffBgrZ(0)
  80    ,fh1EffBgrXi(0)
  81    ,fh1FFTrackPtBackFolded(0)
  82    ,fh1FFZBackFolded(0)
  83    ,fh1FFXiBackFolded(0)
  84    ,fh1FFRatioTrackPtFolded(0)
  85    ,fh1FFRatioZFolded(0)
  86    ,fh1FFRatioXiFolded(0)
  87    ,fh1FFRatioTrackPtBackFolded(0)
  88    ,fh1FFRatioZBackFolded(0)
  89    ,fh1FFRatioXiBackFolded(0)
  90    ,fh1SingleTrackPtBackFolded(0)
  91    ,fh1RatioSingleTrackPtFolded(0)
  92    ,fh1RatioSingleTrackPtBackFolded(0)
  93    ,fhnResponseSinglePt(0)
  94    ,fhnResponsePt(0)
  95    ,fhnResponseZ(0)
  96    ,fhnResponseXi(0)
  97    ,fh1FFTrackPtPrior(0)
  98    ,fh1FFZPrior(0)
  99    ,fh1FFXiPrior(0)
 100    ,fh1SecCorrSinglePt(0)
 101 {
 102    // default constructor
 103 }
 104
 105 //________________________________________________________________________________________________________________________
 106 AliFragmentationFunctionCorrections::AliFragmentationFunctionCorrections(const  AliFragmentationFunctionCorrections &copy)
 107   : TObject()
 108   ,fDebug(copy.fDebug)
 109   ,fNJetPtSlices(copy.fNJetPtSlices)
 110   ,fJetPtSlices(copy.fJetPtSlices)
 111   ,fNJets(copy.fNJets)
 112   ,fNJetsBgr(copy.fNJetsBgr)
 113   ,fNHistoBinsSinglePt(copy.fNHistoBinsSinglePt)
 114   ,fHistoBinsSinglePt(copy.fHistoBinsSinglePt)
 115   ,fNHistoBinsPt(copy.fNHistoBinsPt)
 116   ,fNHistoBinsZ(copy.fNHistoBinsZ)
 117   ,fNHistoBinsXi(copy.fNHistoBinsXi)
 118   ,fHistoBinsPt(copy.fHistoBinsPt)
 119   ,fHistoBinsZ(copy.fHistoBinsZ)
 120   ,fHistoBinsXi(copy.fHistoBinsXi)
 121   ,fNCorrectionLevels(copy.fNCorrectionLevels)
 122   ,fCorrFF(copy.fCorrFF)
 123   ,fNCorrectionLevelsBgr(copy.fNCorrectionLevelsBgr)
 124   ,fCorrBgr(copy.fCorrBgr)
 125   ,fNCorrectionLevelsSinglePt(copy.fNCorrectionLevelsSinglePt)
 126   ,fCorrSinglePt(copy.fCorrSinglePt)
 127   ,fh1FFXiShift(copy.fh1FFXiShift)
 128   ,fh1EffSinglePt(copy.fh1EffSinglePt)
 129   ,fh1EffPt(copy.fh1EffPt)
 130   ,fh1EffZ(copy.fh1EffZ)
 131   ,fh1EffXi(copy.fh1EffXi)
 132   ,fh1EffBgrPt(copy.fh1EffBgrPt)
 133   ,fh1EffBgrZ(copy.fh1EffBgrZ)
 134   ,fh1EffBgrXi(copy.fh1EffBgrXi)
 135   ,fh1FFTrackPtBackFolded(copy.fh1FFTrackPtBackFolded)
 136   ,fh1FFZBackFolded(copy.fh1FFZBackFolded)
 137   ,fh1FFXiBackFolded(copy.fh1FFXiBackFolded)
 138   ,fh1FFRatioTrackPtFolded(copy.fh1FFRatioTrackPtFolded)
 139   ,fh1FFRatioZFolded(copy.fh1FFRatioZFolded)
 140   ,fh1FFRatioXiFolded(copy.fh1FFRatioXiFolded)
 141   ,fh1FFRatioTrackPtBackFolded(copy.fh1FFRatioTrackPtBackFolded)
 142   ,fh1FFRatioZBackFolded(copy.fh1FFRatioZBackFolded)
 143   ,fh1FFRatioXiBackFolded(copy.fh1FFRatioXiBackFolded)
 144   ,fh1SingleTrackPtBackFolded(copy.fh1SingleTrackPtBackFolded)
 145   ,fh1RatioSingleTrackPtFolded(copy.fh1RatioSingleTrackPtFolded)
 146   ,fh1RatioSingleTrackPtBackFolded(copy.fh1RatioSingleTrackPtBackFolded)
 147   ,fhnResponseSinglePt(copy.fhnResponseSinglePt)
 148   ,fhnResponsePt(copy.fhnResponsePt)
 149   ,fhnResponseZ(copy.fhnResponseZ)
 150   ,fhnResponseXi(copy.fhnResponseXi)
 151   ,fh1FFTrackPtPrior(copy.fh1FFTrackPtPrior)
 152   ,fh1FFZPrior(copy.fh1FFZPrior)
 153   ,fh1FFXiPrior(copy.fh1FFXiPrior)
 154   ,fh1SecCorrSinglePt(copy.fh1SecCorrSinglePt)
 155 {
 156   // copy constructor
 157
 158 }
 159
 160 // ______________________________________________________________________________________________________________________________
 161 AliFragmentationFunctionCorrections& AliFragmentationFunctionCorrections::operator=(const AliFragmentationFunctionCorrections& o)
 162 {
 163   // assignment
 164
 165   if(this!=&o){
 166     TObject::operator=(o);
 167     fDebug                          = o.fDebug;
 168     fNJetPtSlices                   = o.fNJetPtSlices;
 169     fJetPtSlices                    = o.fJetPtSlices;
 170     fNJets                          = o.fNJets;
 171     fNJetsBgr                       = o.fNJetsBgr;
 172     fNHistoBinsSinglePt             = o.fNHistoBinsSinglePt;
 173     fHistoBinsSinglePt              = o.fHistoBinsSinglePt;
 174     fNHistoBinsPt                   = o.fNHistoBinsPt;
 175     fNHistoBinsZ                    = o.fNHistoBinsZ;
 176     fNHistoBinsXi                   = o.fNHistoBinsXi;
 177     fHistoBinsPt                    = o.fHistoBinsPt;
 178     fHistoBinsZ                     = o.fHistoBinsZ;
 179     fHistoBinsXi                    = o.fHistoBinsXi;
 180     fNCorrectionLevels              = o.fNCorrectionLevels;
 181     fCorrFF                         = o.fCorrFF;
 182     fNCorrectionLevelsBgr           = o.fNCorrectionLevelsBgr;
 183     fCorrBgr                        = o.fCorrBgr;
 184     fNCorrectionLevelsSinglePt      = o.fNCorrectionLevelsSinglePt;
 185     fCorrSinglePt                   = o.fCorrSinglePt;
 186     fh1FFXiShift                    = o.fh1FFXiShift;
 187     fh1EffSinglePt                  = o.fh1EffSinglePt;
 188     fh1EffPt                        = o.fh1EffPt;
 189     fh1EffZ                         = o.fh1EffZ;
 190     fh1EffXi                        = o.fh1EffXi;
 191     fh1EffBgrPt                     = o.fh1EffBgrPt;
 192     fh1EffBgrZ                      = o.fh1EffBgrZ;
 193     fh1EffBgrXi                     = o.fh1EffBgrXi;
 194     fh1FFTrackPtBackFolded          = o.fh1FFTrackPtBackFolded;
 195     fh1FFZBackFolded                = o.fh1FFZBackFolded;
 196     fh1FFXiBackFolded               = o.fh1FFXiBackFolded;
 197     fh1FFRatioTrackPtFolded         = o.fh1FFRatioTrackPtFolded;
 198     fh1FFRatioZFolded               = o.fh1FFRatioZFolded;
 199     fh1FFRatioXiFolded              = o.fh1FFRatioXiFolded;
 200     fh1FFRatioTrackPtBackFolded     = o.fh1FFRatioTrackPtBackFolded;
 201     fh1FFRatioZBackFolded           = o.fh1FFRatioZBackFolded;
 202     fh1FFRatioXiBackFolded          = o.fh1FFRatioXiBackFolded;
 203     fh1SingleTrackPtBackFolded      = o.fh1SingleTrackPtBackFolded;
 204     fh1RatioSingleTrackPtFolded     = o.fh1RatioSingleTrackPtFolded;
 205     fh1RatioSingleTrackPtBackFolded = o.fh1RatioSingleTrackPtBackFolded;
 206     fhnResponseSinglePt             = o.fhnResponseSinglePt;
 207     fhnResponsePt                   = o.fhnResponsePt;
 208     fhnResponseZ                    = o.fhnResponseZ;
 209     fhnResponseXi                   = o.fhnResponseXi;
 210     fh1FFTrackPtPrior               = o.fh1FFTrackPtPrior;
 211     fh1FFZPrior                     = o.fh1FFZPrior;
 212     fh1FFXiPrior                    = o.fh1FFXiPrior;
 213     fh1SecCorrSinglePt              = o.fh1SecCorrSinglePt;
 214   }
 215
 216   return *this;
 217 }
 218
 219 //_________________________________________________________________________
 220 AliFragmentationFunctionCorrections::~AliFragmentationFunctionCorrections()
 221 {
 222   // destructor
 223
 224   if(fJetPtSlices) delete fJetPtSlices;
 225   if(fNJets)       delete fNJets;
 226   if(fNJetsBgr)    delete fNJetsBgr;
 227
 228   DeleteHistoArray(fh1FFXiShift);
 229
 230   DeleteHistoArray(fh1EffPt);
 231   DeleteHistoArray(fh1EffXi);
 232   DeleteHistoArray(fh1EffZ );
 233
 234   DeleteHistoArray(fh1EffBgrPt);
 235   DeleteHistoArray(fh1EffBgrXi);
 236   DeleteHistoArray(fh1EffBgrZ);
 237
 238   // unfolding
 239
 240   DeleteHistoArray(fh1FFTrackPtBackFolded);
 241   DeleteHistoArray(fh1FFZBackFolded);
 242   DeleteHistoArray(fh1FFXiBackFolded);
 243
 244   DeleteHistoArray(fh1FFRatioTrackPtFolded);
 245   DeleteHistoArray(fh1FFRatioZFolded);
 246   DeleteHistoArray(fh1FFRatioXiFolded);
 247
 248   DeleteHistoArray(fh1FFRatioTrackPtBackFolded);
 249   DeleteHistoArray(fh1FFRatioZBackFolded);
 250   DeleteHistoArray(fh1FFRatioXiBackFolded);
 251
 252   DeleteTHnSparseArray(fhnResponsePt);
 253   DeleteTHnSparseArray(fhnResponseZ);
 254   DeleteTHnSparseArray(fhnResponseXi);
 255
 256   DeleteHistoArray(fh1FFTrackPtPrior);
 257   DeleteHistoArray(fh1FFZPrior);
 258   DeleteHistoArray(fh1FFXiPrior);
 259
 260
 261   // clean up corrected FF
 262
 263   for(Int_t c=0; c<fNCorrectionLevels; c++) delete  fCorrFF[c];
 264   delete[] fCorrFF;
 265
 266   // clean up bgr
 267
 268   for(Int_t c=0; c<fNCorrectionLevelsBgr; c++) delete fCorrBgr[c];
 269   delete[] fCorrBgr;
 270
 271   // clean up inclusive pt
 272   for(Int_t c=0; c<fNCorrectionLevelsSinglePt; c++) delete fCorrSinglePt[c];
 273   delete[] fCorrSinglePt;
 274
 275   delete[] fNHistoBinsPt;
 276   delete[] fNHistoBinsZ;
 277   delete[] fNHistoBinsXi;
 278
 279   // clean up histo bins
 280
 281   if(fHistoBinsSinglePt) delete fHistoBinsSinglePt;
 282
 283   for(Int_t i=0; i<fNJetPtSlices; i++) delete fHistoBinsPt[i];
 284   for(Int_t i=0; i<fNJetPtSlices; i++) delete fHistoBinsZ[i];
 285   for(Int_t i=0; i<fNJetPtSlices; i++) delete fHistoBinsXi[i];
 286
 287   delete[] fHistoBinsPt;
 288   delete[] fHistoBinsZ;
 289   delete[] fHistoBinsXi;
 290 }
 291
 292 //_________________________________________________________________________________
 293 AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::AliFragFuncCorrHistos()
 294   : TObject()
 295   ,fArraySize(0)
 296   ,fh1CorrFFTrackPt(0)
 297   ,fh1CorrFFZ(0)
 298   ,fh1CorrFFXi(0)
 299   ,fCorrLabel(0)
 300 {
 301   // default constructor
 302
 303 }
 304
 305 //__________________________________________________________________________________________________________________
 306 AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::AliFragFuncCorrHistos(const AliFragFuncCorrHistos& copy)
 307   : TObject()
 308   ,fArraySize(copy.fArraySize)
 309   ,fh1CorrFFTrackPt(copy.fh1CorrFFTrackPt)
 310   ,fh1CorrFFZ(copy.fh1CorrFFZ)
 311   ,fh1CorrFFXi(copy.fh1CorrFFXi)
 312   ,fCorrLabel(copy.fCorrLabel)
 313 {
 314   // copy constructor
 315   for(Int_t i=0; i<copy.fArraySize; i++){
 316     fh1CorrFFTrackPt[i] = copy.fh1CorrFFTrackPt[i];
 317     fh1CorrFFZ[i]       = copy.fh1CorrFFZ[i];
 318     fh1CorrFFXi[i]      = copy.fh1CorrFFXi[i];
 319   }
 320 }
 321
 322 //_______________________________________________________________________________________________________________________________________________________________
 323 AliFragmentationFunctionCorrections::AliFragFuncCorrHistos& AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::operator=(const AliFragmentationFunctionCorrections::AliFragFuncCorrHistos& o)
 324 {
 325   // assignment
 326
 327   if(this!=&o){
 328     TObject::operator=(o);
 329     fArraySize       = o.fArraySize;
 330     fCorrLabel       = o.fCorrLabel;
 331
 332     for(Int_t i=0; i<o.fArraySize; i++){
 333       fh1CorrFFTrackPt[i] = o.fh1CorrFFTrackPt[i];
 334       fh1CorrFFZ[i]       = o.fh1CorrFFZ[i];
 335       fh1CorrFFXi[i]      = o.fh1CorrFFXi[i];
 336     }
 337   }
 338
 339   return *this;
 340 }
 341
 342 //__________________________________________________________________________________
 343 AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::~AliFragFuncCorrHistos()
 344 {
 345   // destructor
 346
 347   if(fArraySize){
 348     for(int i=0; i<fArraySize; i++) delete fh1CorrFFTrackPt[i];
 349     for(int i=0; i<fArraySize; i++) delete fh1CorrFFZ[i];
 350     for(int i=0; i<fArraySize; i++) delete fh1CorrFFXi[i];
 351   }
 352
 353   if(fh1CorrFFTrackPt) delete[] fh1CorrFFTrackPt;
 354   if(fh1CorrFFZ)       delete[] fh1CorrFFZ;
 355   if(fh1CorrFFXi)      delete[] fh1CorrFFXi;
 356
 357 }
 358
 359 //___________________________________________________________________________________________________________________
 360 AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::AliFragFuncCorrHistos(const char* label, Int_t arraySize)
 361   : TObject()
 362   ,fArraySize(0)
 363   ,fh1CorrFFTrackPt(0)
 364   ,fh1CorrFFZ(0)
 365   ,fh1CorrFFXi(0)
 366   ,fCorrLabel(label)
 367 {
 368   // constructor
 369
 370   fArraySize = arraySize;
 371   fh1CorrFFTrackPt = new TH1F*[arraySize];
 372   fh1CorrFFZ       = new TH1F*[arraySize];
 373   fh1CorrFFXi      = new TH1F*[arraySize];
 374
 375   for(int i=0; i<arraySize; i++) fh1CorrFFTrackPt[i] = new TH1F;
 376   for(int i=0; i<arraySize; i++) fh1CorrFFZ[i]       = new TH1F;
 377   for(int i=0; i<arraySize; i++) fh1CorrFFXi[i]      = new TH1F;
 378 }
 379
 380 //_______________________________________________________________________________________________________________________________
 381 void AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::AddCorrHistos(Int_t slice,TH1F* histPt,TH1F* histZ,TH1F* histXi)
 382 {
 383   // place histo in array, add corrLabel to histo name
 384
 385   if(slice>=fArraySize){
 386     Printf("%s:%d -- slice > array size", (char*)__FILE__,__LINE__);
 387     return;
 388   }
 389
 390   if(histPt){
 391     TString name = histPt->GetName();
 392     if(fCorrLabel.Length()>0) name += "_"+fCorrLabel;
 393     histPt->SetName(name);
 394
 395     if(!(histPt->GetSumw2())) histPt->Sumw2();
 396
 397     new(fh1CorrFFTrackPt[slice]) TH1F(*histPt);
 398   }
 399
 400   if(histZ){
 401     TString name = histZ->GetName();
 402     if(fCorrLabel.Length()>0) name += "_"+fCorrLabel;
 403     histZ->SetName(name);
 404
 405     if(!(histZ->GetSumw2())) histZ->Sumw2();
 406
 407     new(fh1CorrFFZ[slice]) TH1F(*histZ);
 408   }
 409
 410   if(histXi){
 411     TString name = histXi->GetName();
 412     if(fCorrLabel.Length()>0) name += "_"+fCorrLabel;
 413     histXi->SetName(name);
 414
 415     if(!(histXi->GetSumw2())) histXi->Sumw2();
 416
 417     new(fh1CorrFFXi[slice]) TH1F(*histXi);
 418   }
 419 }
 420
 421 //___________________________________________________________________________________________________________________________________
 422 void AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::ReplaceCorrHistos(Int_t slice,TH1F* histPt,TH1F* histZ,TH1F* histXi)
 423 {
 424   // replace histo in array
 425
 426   if(slice>=fArraySize){
 427     Printf("%s:%d -- slice > array size", (char*)__FILE__,__LINE__);
 428     return;
 429   }
 430
 431   if(histPt){
 432     if(!(histPt->GetSumw2())) histPt->Sumw2();
 433
 434     TString name = histPt->GetName();
 435     histPt->SetName(name);
 436
 437     delete fh1CorrFFTrackPt[slice];
 438     fh1CorrFFTrackPt[slice] =  new TH1F;
 439     new(fh1CorrFFTrackPt[slice]) TH1F(*histPt);
 440   }
 441
 442   if(histZ){
 443     if(!(histZ->GetSumw2())) histZ->Sumw2();
 444
 445     TString name = histZ->GetName();
 446     histZ->SetName(name);
 447
 448     delete fh1CorrFFZ[slice];
 449     fh1CorrFFZ[slice] = new TH1F;
 450     new(fh1CorrFFZ[slice]) TH1F(*histZ);
 451   }
 452
 453   if(histXi){
 454     if(!(histXi->GetSumw2())) histXi->Sumw2();
 455
 456     TString name = histXi->GetName();
 457     histXi->SetName(name);
 458
 459     delete fh1CorrFFXi[slice];
 460     fh1CorrFFXi[slice] = new TH1F;
 461     new(fh1CorrFFXi[slice]) TH1F(*histXi);
 462   }
 463 }
 464
 465 // ___________________________________________________________________________________________
 466 TH1F* AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::GetTrackPt(const Int_t slice)
 467 {
 468   // return pt histo
 469
 470   if(slice>=fArraySize){
 471     Printf("%s:%d -- slice > array size", (char*)__FILE__,__LINE__);
 472     return 0;
 473   }
 474
 475   return fh1CorrFFTrackPt[slice];
 476
 477 }
 478
 479 // ______________________________________________________________________________________
 480 TH1F* AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::GetZ(const Int_t slice)
 481 {
 482   // return z histo
 483
 484   if(slice>=fArraySize){
 485     Printf("%s:%d -- slice > array size", (char*)__FILE__,__LINE__);
 486     return 0;
 487   }
 488
 489   return fh1CorrFFZ[slice];
 490 }
 491
 492 // ________________________________________________________________________________________
 493 TH1F* AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::GetXi(const Int_t slice)
 494 {
 495   // return xi histo
 496
 497   if(slice>=fArraySize){
 498     Printf("%s:%d -- slice > array size", (char*)__FILE__,__LINE__);
 499     return 0;
 500   }
 501
 502   return fh1CorrFFXi[slice];
 503 }
 504
 505 // __________________________________________________________________________
 506 void AliFragmentationFunctionCorrections::DeleteHistoArray(TH1F** hist) const
 507 {
 508   // delete array of TH1
 509
 510   if(hist){
 511     for(Int_t i=0; i<fNJetPtSlices; i++) delete hist[i];
 512     delete[] hist;
 513   }
 514 }
 515
 516 // ____________________________________________________________________________________
 517 void AliFragmentationFunctionCorrections::DeleteTHnSparseArray(THnSparse** hist) const
 518 {
 519   // delete array of THnSparse
 520
 521   if(hist){
 522     for(Int_t i=0; i<fNJetPtSlices; i++) delete hist[i];
 523     delete[] hist;
 524   }
 525 }
 526
 527 // _________________________________________________________
 528 TH1F** AliFragmentationFunctionCorrections::BookHistoArray()
 529 {
 530   // book array of TH1
 531
 532   if(!fNJetPtSlices){
 533     Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
 534     return 0;
 535   }
 536
 537   TH1F** hist = new TH1F*[fNJetPtSlices];
 538   for(Int_t i=0; i<fNJetPtSlices; i++) hist[i] = new TH1F();
 539
 540   return hist;
 541 }
 542
 543 //__________________________________________________________________
 544 THnSparse** AliFragmentationFunctionCorrections::BookTHnSparseArray()
 545 {
 546   // book array of THnSparse
 547
 548   if(!fNJetPtSlices){
 549     Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
 550     return 0;
 551   }
 552
 553   THnSparse** hist = new THnSparse*[fNJetPtSlices];
 554   for(Int_t i=0; i<fNJetPtSlices; i++) hist[i] = new THnSparseF();
 555
 556   return hist;
 557 }
 558
 559 //_____________________________________________________________________________
 560 void AliFragmentationFunctionCorrections::AddCorrectionLevel(const char* label)
 561 {
 562   // increase corr level
 563
 564   if(!fNJetPtSlices){
 565     Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
 566     return;
 567   }
 568
 569   if(fNCorrectionLevels >= fgMaxNCorrectionLevels){
 570     Printf("%s:%d -- max correction level exceeded", (char*)__FILE__,__LINE__);
 571     return;
 572   }
 573
 574   fCorrFF[fNCorrectionLevels] = new AliFragFuncCorrHistos(label,fNJetPtSlices);
 575   fNCorrectionLevels++;
 576 }
 577
 578
 579 //________________________________________________________________________________
 580 void AliFragmentationFunctionCorrections::AddCorrectionLevelBgr(const char* label)
 581 {
 582   // increase corr level for bgr FF
 583
 584   if(!fNJetPtSlices){
 585     if(fDebug>0)  Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
 586     return;
 587   }
 588
 589   if(fNCorrectionLevelsBgr >= fgMaxNCorrectionLevels){
 590     Printf("%s:%d -- max correction level exceeded", (char*)__FILE__,__LINE__);
 591     return;
 592   }
 593
 594   fCorrBgr[fNCorrectionLevelsBgr] = new AliFragFuncCorrHistos(label,fNJetPtSlices);
 595   fNCorrectionLevelsBgr++;
 596 }
 597
 598 //_____________________________________________________________________________________
 599 void AliFragmentationFunctionCorrections::AddCorrectionLevelSinglePt(const char* label)
 600 {
 601   // increase corr level single pt spec
 602
 603   Int_t nJetPtSlicesSingle = 1; // pro forma
 604
 605   if(fNCorrectionLevelsSinglePt >= fgMaxNCorrectionLevels){
 606     Printf("%s:%d -- max correction level exceeded", (char*)__FILE__,__LINE__);
 607     return;
 608   }
 609
 610   fCorrSinglePt[fNCorrectionLevelsSinglePt] = new AliFragFuncCorrHistos(label,nJetPtSlicesSingle);
 611   fNCorrectionLevelsSinglePt++;
 612 }
 613
 614 //_____________________________________________________________________________________________
 615 void AliFragmentationFunctionCorrections::SetJetPtSlices(Float_t* bins, const Int_t nJetPtSlices)
 616 {
 617   // set jet pt slices
 618   // once slices are known, can book all other histos
 619
 620   fNJetPtSlices = nJetPtSlices;
 621
 622   const Int_t size = nJetPtSlices+1;
 623   fJetPtSlices  = new TArrayF(size,bins);
 624
 625   // nJets array
 626
 627   fNJets = new TArrayF(size);
 628   fNJets->Reset(0);
 629
 630   fNJetsBgr = new TArrayF(size);
 631   fNJetsBgr->Reset(0);
 632
 633   // histos
 634
 635   fNCorrectionLevels = 0;
 636   fCorrFF = new AliFragFuncCorrHistos*[fgMaxNCorrectionLevels];
 637   AddCorrectionLevel(); // first 'correction' level = raw FF
 638
 639   fNCorrectionLevelsBgr = 0;
 640   fCorrBgr = new AliFragFuncCorrHistos*[fgMaxNCorrectionLevels];
 641   AddCorrectionLevelBgr(); // first 'correction' level = raw bgr dist
 642
 643   fh1FFXiShift = BookHistoArray();
 644
 645   // eff histos
 646
 647   fh1EffPt = BookHistoArray();
 648   fh1EffXi = BookHistoArray();
 649   fh1EffZ  = BookHistoArray();
 650
 651   fh1EffBgrPt = BookHistoArray();
 652   fh1EffBgrXi = BookHistoArray();
 653   fh1EffBgrZ  = BookHistoArray();
 654
 655
 656   // unfolding
 657
 658   fh1FFTrackPtBackFolded = BookHistoArray();
 659   fh1FFXiBackFolded      = BookHistoArray();
 660   fh1FFZBackFolded       = BookHistoArray();
 661
 662   fh1FFRatioTrackPtFolded = BookHistoArray();
 663   fh1FFRatioXiFolded      = BookHistoArray();
 664   fh1FFRatioZFolded       = BookHistoArray();
 665
 666   fh1FFRatioTrackPtBackFolded = BookHistoArray();
 667   fh1FFRatioXiBackFolded      = BookHistoArray();
 668   fh1FFRatioZBackFolded       = BookHistoArray();
 669
 670   //
 671
 672   fhnResponsePt      = BookTHnSparseArray();
 673   fhnResponseZ       = BookTHnSparseArray();
 674   fhnResponseXi      = BookTHnSparseArray();
 675
 676   fh1FFTrackPtPrior =  BookHistoArray();
 677   fh1FFZPrior       =  BookHistoArray();
 678   fh1FFXiPrior      =  BookHistoArray();
 679
 680
 681   // histos bins
 682
 683   fNHistoBinsPt = new Int_t[fNJetPtSlices];
 684   fNHistoBinsXi = new Int_t[fNJetPtSlices];
 685   fNHistoBinsZ  = new Int_t[fNJetPtSlices];
 686
 687   for(Int_t i=0; i<fNJetPtSlices; i++) fNHistoBinsPt[i] = 0;
 688   for(Int_t i=0; i<fNJetPtSlices; i++) fNHistoBinsXi[i] = 0;
 689   for(Int_t i=0; i<fNJetPtSlices; i++) fNHistoBinsZ[i]  = 0;
 690
 691   fHistoBinsPt = new TArrayD*[fNJetPtSlices];
 692   fHistoBinsXi = new TArrayD*[fNJetPtSlices];
 693   fHistoBinsZ  = new TArrayD*[fNJetPtSlices];
 694
 695   for(Int_t i=0; i<fNJetPtSlices; i++) fHistoBinsPt[i] = new TArrayD(0);
 696   for(Int_t i=0; i<fNJetPtSlices; i++) fHistoBinsXi[i] = new TArrayD(0);
 697   for(Int_t i=0; i<fNJetPtSlices; i++) fHistoBinsZ[i]  = new TArrayD(0);
 698 }
 699
 700 //_____________________________________________________________________________________________________________________________________
 701 void AliFragmentationFunctionCorrections::SetHistoBins(const Int_t jetPtSlice, const Int_t sizeBins, Double_t* bins, const Int_t type)
 702 {
 703   // set histo bins for jet pt slice
 704   // if binning undefined for any slice, original binning will be used
 705
 706   if(!fNJetPtSlices){
 707     Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
 708     return;
 709   }
 710
 711   if(jetPtSlice>=fNJetPtSlices){
 712     Printf("%s:%d -- jetPtSlice %d exceeds max",(char*)__FILE__,__LINE__,jetPtSlice);
 713     return;
 714   }
 715
 716   if(type == kFlagPt){
 717     fNHistoBinsPt[jetPtSlice] = sizeBins-1;
 718     fHistoBinsPt[jetPtSlice]->Set(sizeBins,bins);
 719   }
 720   else if(type==kFlagZ){
 721     fNHistoBinsZ[jetPtSlice] = sizeBins-1;
 722     fHistoBinsZ[jetPtSlice]->Set(sizeBins,bins);
 723   }
 724
 725   else if(type==kFlagXi){
 726     fNHistoBinsXi[jetPtSlice] = sizeBins-1;
 727     fHistoBinsXi[jetPtSlice]->Set(sizeBins,bins);
 728   }
 729 }
 730
 731 //__________________________________________________________________________________________________________________________________________________________________
 732 void AliFragmentationFunctionCorrections::SetHistoBins(const Int_t jetPtSlice, const Int_t nBinsLimits, Double_t* binsLimits, Double_t* binsWidth, const Int_t type)
 733 {
 734   // set histo bins for jet pt slice
 735   // function takes array of limits and widths (e.g. 1 GeV bins up to 10 GeV, 2 GeV width up to 20 GeV, ...)
 736   // array size of binsLimits: nBinsLimits
 737   // array size of binsWidth: nBinsLimits-1
 738   // binsLimits have to be in increasing order
 739   // if binning undefined for any slice, original binning will be used
 740
 741   if(!fNJetPtSlices){
 742     Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
 743     return;
 744   }
 745
 746   if(jetPtSlice>=fNJetPtSlices){
 747     Printf("%s:%d -- jetPtSlice %d exceeds max",(char*)__FILE__,__LINE__,jetPtSlice);
 748     return;
 749   }
 750
 751
 752   Double_t binLimitMin = binsLimits[0];
 753   Double_t binLimitMax = binsLimits[nBinsLimits-1];
 754
 755   Double_t binLimit = binLimitMin; // start value
 756
 757   Int_t sizeUpperLim = 10000; //static_cast<Int_t>(binLimitMax/binsWidth[0])+1;
 758   TArrayD binsArray(sizeUpperLim);
 759   Int_t nBins = 0;
 760   binsArray.SetAt(binLimitMin,nBins++);
 761
 762   while(binLimit<binLimitMax && nBins<sizeUpperLim){
 763
 764     Int_t currentSlice = -1;
 765     for(Int_t i=0; i<nBinsLimits; i++){
 766       if(binLimit >= 0.999*binsLimits[i]) currentSlice = i; // 0.999 numerical saftey factor
 767     }
 768
 769     Double_t currentBinWidth = binsWidth[currentSlice];
 770     binLimit += currentBinWidth;
 771
 772     binsArray.SetAt(binLimit,nBins++);
 773   }
 774
 775   Double_t* bins = binsArray.GetArray();
 776
 777   SetHistoBins(jetPtSlice,nBins,bins,type);
 778 }
 779
 780 //__________________________________________________________________________________________________
 781 void AliFragmentationFunctionCorrections::SetHistoBinsSinglePt(const Int_t sizeBins, Double_t* bins)
 782 {
 783   // set histo bins for inclusive pt spectra
 784   // if binning undefined, original binning will be used
 785
 786   fNHistoBinsSinglePt = sizeBins-1;
 787
 788   fHistoBinsSinglePt = new TArrayD(sizeBins);
 789   fHistoBinsSinglePt->Set(sizeBins,bins);
 790 }
 791
 792 //__________________________________________________________________________________________________________________________________________________________________
 793 void AliFragmentationFunctionCorrections::SetHistoBinsSinglePt(const Int_t nBinsLimits, Double_t* binsLimits, Double_t* binsWidth)
 794 {
 795   // set histo bins for inclusive pt spectra
 796   // function takes array of limits and widths (e.g. 1 GeV bins up to 10 GeV, 2 GeV width up to 20 GeV, ...)
 797   // array size of binsLimits: nBinsLimits
 798   // array size of binsWidth: nBinsLimits-1
 799   // binsLimits have to be in increasing order
 800   // if binning undefined for any slice, original binning will be used
 801
 802
 803   Double_t binLimitMin = binsLimits[0];
 804   Double_t binLimitMax = binsLimits[nBinsLimits-1];
 805
 806   Double_t binLimit = binLimitMin; // start value
 807
 808   Int_t sizeUpperLim = 10000; //static_cast<Int_t>(binLimitMax/binsWidth[0])+1;
 809   TArrayD binsArray(sizeUpperLim);
 810   Int_t nBins = 0;
 811   binsArray.SetAt(binLimitMin,nBins++);
 812
 813   while(binLimit<binLimitMax && nBins<sizeUpperLim){
 814
 815     Int_t currentSlice = -1;
 816     for(Int_t i=0; i<nBinsLimits; i++){
 817       if(binLimit >= 0.999*binsLimits[i]) currentSlice = i; // 0.999 numerical saftey factor
 818     }
 819
 820     Double_t currentBinWidth = binsWidth[currentSlice];
 821     binLimit += currentBinWidth;
 822
 823     binsArray.SetAt(binLimit,nBins++);
 824   }
 825
 826   Double_t* bins = binsArray.GetArray();
 827
 828   SetHistoBinsSinglePt(nBins,bins);
 829 }
 830
 831 //____________________________________________________________________________________
 832 void AliFragmentationFunctionCorrections::NormalizeTH1(TH1* hist, const Float_t nJets)
 833 {
 834   // FF normalization: divide by bin width and normalize to entries/jets
 835   // should also work for TH2, but to be tested !!!
 836
 837   if(nJets == 0){ // nothing to do
 838     if(fDebug>0)  Printf("%s:%d -- normalize: nJets = 0, do nothing", (char*)__FILE__,__LINE__);
 839     return;
 840   }
 841
 842   Int_t nBins = hist->GetNbinsX()*hist->GetNbinsY()*hist->GetNbinsZ();
 843
 844   for(Int_t bin=0; bin<=nBins; bin++){ // include under-/overflow (?)
 845
 846     Double_t binWidth = hist->GetBinWidth(bin);
 847     Double_t binCont  = hist->GetBinContent(bin);
 848     Double_t binErr   = hist->GetBinError(bin);
 849
 850     binCont /= binWidth;
 851     binErr  /= binWidth;
 852
 853     hist->SetBinContent(bin,binCont);
 854     hist->SetBinError(bin,binErr);
 855   }
 856
 857   hist->Scale(1/nJets);
 858 }
 859
 860 //_____________________________________________________
 861 void AliFragmentationFunctionCorrections::NormalizeFF()
 862 {
 863   // normalize FF
 864
 865   if(fNCorrectionLevels>1){
 866     Printf("%s:%d -- FF normalization should be done BEFORE any correction !!!", (char*)__FILE__,__LINE__);
 867   }
 868
 869   for(Int_t i=0; i<fNJetPtSlices; i++){
 870
 871     if(fDebug>0) Printf(" normalizeFF: i %d, nJets %f",i,fNJets->At(i));
 872
 873     NormalizeTH1(fCorrFF[0]->GetTrackPt(i),fNJets->At(i)); // always normalize corr level 0 = raw FF
 874     NormalizeTH1(fCorrFF[0]->GetZ(i),fNJets->At(i));
 875     NormalizeTH1(fCorrFF[0]->GetXi(i),fNJets->At(i));
 876   }
 877 }
 878
 879 //______________________________________________________
 880 void AliFragmentationFunctionCorrections::NormalizeBgr()
 881 {
 882   // normalize bgr/UE FF
 883
 884   if(fNCorrectionLevelsBgr>1){
 885     Printf("%s:%d -- FF normalization should be done BEFORE any correction !!!", (char*)__FILE__,__LINE__);
 886   }
 887
 888   for(Int_t i=0; i<fNJetPtSlices; i++){
 889     NormalizeTH1(fCorrBgr[0]->GetTrackPt(i), fNJetsBgr->At(i)); // always normalize corr level 0 = raw FF
 890     NormalizeTH1(fCorrBgr[0]->GetZ(i), fNJetsBgr->At(i));
 891     NormalizeTH1(fCorrBgr[0]->GetXi(i),fNJetsBgr->At(i));
 892   }
 893
 894 }
 895
 896 //__________________________________________________________________________________________________
 897 void AliFragmentationFunctionCorrections::ReadRawFF(TString strfile, TString strID, TString strFFID)
 898 {
 899   // read raw FF - standard dir/list name
 900
 901   TString strdir  = "PWG4_FragmentationFunction_" + strID;
 902   TString strlist = "fracfunc_" + strID;
 903
 904   ReadRawFF(strfile,strdir,strlist,strFFID);
 905 }
 906
 907 //____________________________________________________________________________________________________________________
 908 void AliFragmentationFunctionCorrections::ReadRawFF(TString strfile, TString strdir, TString strlist, TString strFFID)
 909 {
 910   // get raw FF from input file, project in jet pt slice
 911   // normalization done separately
 912
 913   TFile f(strfile,"READ");
 914
 915   if(!f.IsOpen()){
 916     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
 917     return;
 918   }
 919
 920   if(fDebug>0) Printf("%s:%d -- read FF from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
 921
 922   gDirectory->cd(strdir);
 923
 924   TList* list = 0;
 925
 926   if(!(list = (TList*) gDirectory->Get(strlist))){
 927     Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
 928     return;
 929   }
 930
 931   TString hnameJetPt(Form("fh1FFJetPt%s",strFFID.Data()));
 932   TString hnameTrackPt(Form("fh2FFTrackPt%s",strFFID.Data()));
 933   TString hnameZ(Form("fh2FFZ%s",strFFID.Data()));
 934   TString hnameXi(Form("fh2FFXi%s",strFFID.Data()));
 935
 936   TH1F* fh1FFJetPt   = (TH1F*) list->FindObject(hnameJetPt);
 937   TH2F* fh2FFTrackPt = (TH2F*) list->FindObject(hnameTrackPt);
 938   TH2F* fh2FFZ       = (TH2F*) list->FindObject(hnameZ);
 939   TH2F* fh2FFXi      = (TH2F*) list->FindObject(hnameXi);
 940
 941   if(!fh1FFJetPt)  { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameJetPt.Data());   return; }
 942   if(!fh2FFTrackPt){ Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameTrackPt.Data()); return; }
 943   if(!fh2FFZ)      { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameZ.Data());       return; }
 944   if(!fh2FFXi)     { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameXi.Data());      return; }
 945
 946   fh1FFJetPt->SetDirectory(0);
 947   fh2FFTrackPt->SetDirectory(0);
 948   fh2FFZ->SetDirectory(0);
 949   fh2FFXi->SetDirectory(0);
 950
 951   f.Close();
 952
 953
 954   // nJets per bin
 955
 956   for(Int_t i=0; i<fNJetPtSlices; i++){
 957
 958     Float_t jetPtLoLim = fJetPtSlices->At(i);
 959     Float_t jetPtUpLim = fJetPtSlices->At(i+1);
 960
 961     Int_t binLo = static_cast<Int_t>(fh1FFJetPt->FindBin(jetPtLoLim));
 962     Int_t binUp = static_cast<Int_t>(fh1FFJetPt->FindBin(jetPtUpLim)) - 1;
 963
 964     Float_t nJetsBin = fh1FFJetPt->Integral(binLo,binUp);
 965
 966     fNJets->SetAt(nJetsBin,i);
 967
 968     if(fDebug>0) Printf("jet pt %d to %d: nJets %f",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),fNJets->At(i));
 969   }
 970
 971   // projections: FF
 972
 973   for(Int_t i=0; i<fNJetPtSlices; i++){
 974
 975     Float_t jetPtLoLim = fJetPtSlices->At(i);
 976     Float_t jetPtUpLim = fJetPtSlices->At(i+1);
 977
 978     Int_t binLo = static_cast<Int_t>(fh2FFTrackPt->GetXaxis()->FindBin(jetPtLoLim));
 979     Int_t binUp = static_cast<Int_t>(fh2FFTrackPt->GetXaxis()->FindBin(jetPtUpLim))-1;
 980
 981     if(binUp > fh2FFTrackPt->GetNbinsX()){
 982       Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim);
 983       return;
 984     }
 985
 986     TString strNameFFPt(Form("fh1FFTrackPt%s_%02d_%02d",strFFID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
 987     TString strNameFFZ(Form("fh1FFZ%s_%02d_%02d",strFFID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
 988     TString strNameFFXi(Form("fh1FFXi%s_%02d_%02d",strFFID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
 989
 990     // appendix 'unbinned' to avoid histos with same name after rebinning
 991     TH1F* projPt = (TH1F*) fh2FFTrackPt->ProjectionY(strNameFFPt+"_unBinned",binLo,binUp,"o"); // option "o": original axis range
 992     TH1F* projZ  = (TH1F*) fh2FFZ->ProjectionY(strNameFFZ+"_unBinned",binLo,binUp,"o");
 993     TH1F* projXi = (TH1F*) fh2FFXi->ProjectionY(strNameFFXi+"_unBinned",binLo,binUp,"o");
 994
 995     if(fNHistoBinsPt[i]) projPt = (TH1F*) projPt->Rebin(fNHistoBinsPt[i],strNameFFPt,fHistoBinsPt[i]->GetArray());
 996     if(fNHistoBinsZ[i])  projZ  = (TH1F*) projZ->Rebin(fNHistoBinsZ[i],strNameFFZ,fHistoBinsZ[i]->GetArray());
 997     if(fNHistoBinsXi[i]) projXi = (TH1F*) projXi->Rebin(fNHistoBinsXi[i],strNameFFXi,fHistoBinsXi[i]->GetArray());
 998
 999     projPt->SetNameTitle(strNameFFPt,"");
1000     projZ->SetNameTitle(strNameFFZ,"");
1001     projXi->SetNameTitle(strNameFFXi,"");
1002
1003     // raw FF = corr level 0
1004     fCorrFF[0]->AddCorrHistos(i,projPt,projZ,projXi);
1005   }
1006 }
1007
1008 //_____________________________________________________________________________________________________________________
1009 void AliFragmentationFunctionCorrections::ReadRawBgr(TString strfile, TString strID, TString strBgrID, TString strFFID)
1010 {
1011   // read raw FF - standard dir/list name
1012
1013   TString strdir  = "PWG4_FragmentationFunction_" + strID;
1014   TString strlist = "fracfunc_" + strID;
1015
1016   ReadRawBgr(strfile,strdir,strlist,strBgrID,strFFID);
1017 }
1018
1019 //_______________________________________________________________________________________________________________________________________
1020 void AliFragmentationFunctionCorrections::ReadRawBgr(TString strfile, TString strdir, TString strlist, TString strBgrID, TString strFFID)
1021 {
1022   // get raw FF from input file, project in jet pt slice
1023   // use jet dN/dpt corresponding to strFFID, bgr FF to strBgrID+strFFID
1024   // e.g. "fh1FFJetPtRecCuts", "fh2FFXiBgrPerpRecCuts"
1025   // normalization done separately
1026
1027   TString strID = strBgrID + strFFID;
1028
1029   TFile f(strfile,"READ");
1030
1031   if(!f.IsOpen()){
1032     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
1033     return;
1034   }
1035
1036   if(fDebug>0) Printf("%s:%d -- read Bgr %s from file %s ",(char*)__FILE__,__LINE__,strBgrID.Data(),strfile.Data());
1037
1038   gDirectory->cd(strdir);
1039
1040   TList* list = 0;
1041
1042   if(!(list = (TList*) gDirectory->Get(strlist))){
1043     Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
1044     return;
1045   }
1046
1047   TString hnameNJets = "fh1nRecJetsCuts";
1048   TString hnameJetPt(Form("fh1FFJetPt%s",strFFID.Data())); // not: strID.Data() !!! would not be proper normalization
1049   TString hnameBgrTrackPt(Form("fh2FFTrackPt%s",strID.Data()));
1050   TString hnameBgrZ(Form("fh2FFZ%s",strID.Data()));
1051   TString hnameBgrXi(Form("fh2FFXi%s",strID.Data()));
1052
1053   TH1F* fh1NJets        = (TH1F*) list->FindObject(hnameNJets); // needed for normalization of bgr out of 2 jets
1054   TH1F* fh1FFJetPtBgr   = (TH1F*) list->FindObject(hnameJetPt);
1055   TH2F* fh2FFTrackPtBgr = (TH2F*) list->FindObject(hnameBgrTrackPt);
1056   TH2F* fh2FFZBgr       = (TH2F*) list->FindObject(hnameBgrZ);
1057   TH2F* fh2FFXiBgr      = (TH2F*) list->FindObject(hnameBgrXi);
1058
1059   if(!fh1FFJetPtBgr)  { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameJetPt.Data());      return; }
1060   if(!fh1NJets)       { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameNJets.Data());      return; }
1061   if(!fh2FFTrackPtBgr){ Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrTrackPt.Data()); return; }
1062   if(!fh2FFZBgr)      { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrZ.Data());       return; }
1063   if(!fh2FFXiBgr)     { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrXi.Data());      return; }
1064
1065   fh1FFJetPtBgr->SetDirectory(0);
1066   fh1NJets->SetDirectory(0);
1067   fh2FFTrackPtBgr->SetDirectory(0);
1068   fh2FFZBgr->SetDirectory(0);
1069   fh2FFXiBgr->SetDirectory(0);
1070
1071   f.Close();
1072
1073   // nJets per bin
1074
1075   for(Int_t i=0; i<fNJetPtSlices; i++){
1076
1077     Float_t jetPtLoLim = fJetPtSlices->At(i);
1078     Float_t jetPtUpLim = fJetPtSlices->At(i+1);
1079
1080     Int_t binLo = static_cast<Int_t>(fh1FFJetPtBgr->FindBin(jetPtLoLim));
1081     Int_t binUp = static_cast<Int_t>(fh1FFJetPtBgr->FindBin(jetPtUpLim)) - 1;
1082
1083     Float_t nJetsBin = fh1FFJetPtBgr->Integral(binLo,binUp);
1084     Double_t scaleF = 1;
1085
1086     //if(strBgrID.Contains("Out2Jets")){  // scale by ratio 2 jets events / all events
1087     //  scaleF = fh1NJets->Integral(fh1NJets->FindBin(2),fh1NJets->GetNbinsX())
1088     //  / fh1NJets->Integral(fh1NJets->FindBin(1),fh1NJets->GetNbinsX());}
1089
1090
1091     if(strBgrID.Contains("OutAllJets")){  // scale by ratio >3 jets events / all events
1092       scaleF = fh1NJets->Integral(fh1NJets->FindBin(4),fh1NJets->GetNbinsX())
1093         / fh1NJets->Integral(fh1NJets->FindBin(1),fh1NJets->GetNbinsX());
1094     }
1095
1096     fNJetsBgr->SetAt(nJetsBin*scaleF,i);
1097
1098     if(fDebug>0) Printf("bgr jet pt %d to %d: nJets %f, scaleF %.2f",
1099                         static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),nJetsBin,scaleF);
1100
1101   }
1102
1103   // projections: FF
1104
1105   for(Int_t i=0; i<fNJetPtSlices; i++){
1106
1107     Float_t jetPtLoLim = fJetPtSlices->At(i);
1108     Float_t jetPtUpLim = fJetPtSlices->At(i+1);
1109
1110     Int_t binLo = static_cast<Int_t>(fh2FFTrackPtBgr->GetXaxis()->FindBin(jetPtLoLim));
1111     Int_t binUp = static_cast<Int_t>(fh2FFTrackPtBgr->GetXaxis()->FindBin(jetPtUpLim))-1;
1112
1113     if(binUp > fh2FFTrackPtBgr->GetNbinsX()){
1114       Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim);
1115       return;
1116     }
1117
1118     TString strNameBgrPt(Form("fh1BgrTrackPt%s_%02d_%02d",strID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1119     TString strNameBgrZ(Form("fh1BgrZ%s_%02d_%02d",strID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1120     TString strNameBgrXi(Form("fh1BgrXi%s_%02d_%02d",strID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1121
1122     // appendix 'unbinned' to avoid histos with same name after rebinning
1123     TH1F* projPt = (TH1F*) fh2FFTrackPtBgr->ProjectionY(strNameBgrPt+"_unBinned",binLo,binUp,"o"); // option "o": original axis range
1124     TH1F* projZ  = (TH1F*) fh2FFZBgr->ProjectionY(strNameBgrZ+"_unBinned",binLo,binUp,"o");
1125     TH1F* projXi = (TH1F*) fh2FFXiBgr->ProjectionY(strNameBgrXi+"_unBinned",binLo,binUp,"o");
1126
1127     if(fNHistoBinsPt[i]) projPt = (TH1F*) projPt->Rebin(fNHistoBinsPt[i],strNameBgrPt,fHistoBinsPt[i]->GetArray());
1128     if(fNHistoBinsZ[i])  projZ  = (TH1F*) projZ->Rebin(fNHistoBinsZ[i],strNameBgrZ,fHistoBinsZ[i]->GetArray());
1129     if(fNHistoBinsXi[i]) projXi = (TH1F*) projXi->Rebin(fNHistoBinsXi[i],strNameBgrXi,fHistoBinsXi[i]->GetArray());
1130
1131     projPt->SetNameTitle(strNameBgrPt,"");
1132     projZ->SetNameTitle(strNameBgrZ,"");
1133     projXi->SetNameTitle(strNameBgrXi,"");
1134
1135     // raw bgr = corr level 0
1136     fCorrBgr[0]->AddCorrHistos(i,projPt,projZ,projXi);
1137   }
1138 }
1139
1140 //_____________________________________________________________________________________________________________________
1141 void AliFragmentationFunctionCorrections::ReadRawBgrEmbedding(TString strfile, TString strID, TString strFFID)
1142 {
1143   // read raw FF - standard dir/list name
1144
1145   TString strdir  = "PWG4_FragmentationFunction_" + strID;
1146   TString strlist = "fracfunc_" + strID;
1147
1148   ReadRawBgrEmbedding(strfile,strdir,strlist,strFFID);
1149 }
1150
1151 //_______________________________________________________________________________________________________________________________________
1152 void AliFragmentationFunctionCorrections::ReadRawBgrEmbedding(TString strfile, TString strdir, TString strlist, TString strFFID)
1153 {
1154   // get raw FF from input file, project in jet pt slice
1155   // for embedding, the bgr FF are taken from histos "fh1FFJetPtRecCuts", "fh2FFXiRecCuts"
1156   // normalization done separately
1157
1158   TString strBgrID = "BckgEmbed";
1159   TString strID = strBgrID + strFFID;
1160
1161   TFile f(strfile,"READ");
1162
1163   if(!f.IsOpen()){
1164     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
1165     return;
1166   }
1167
1168   if(fDebug>0) Printf("%s:%d -- read Bgr %s from file %s ",(char*)__FILE__,__LINE__,strFFID.Data(),strfile.Data());
1169
1170   gDirectory->cd(strdir);
1171
1172   TList* list = 0;
1173
1174   if(!(list = (TList*) gDirectory->Get(strlist))){
1175     Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
1176     return;
1177   }
1178
1179   TString hnameNJets = "fh1nRecJetsCuts";
1180   TString hnameJetPt(Form("fh1FFJetPt%s",strFFID.Data()));
1181   TString hnameBgrTrackPt(Form("fh2FFTrackPt%s",strFFID.Data()));
1182   TString hnameBgrZ(Form("fh2FFZ%s",strFFID.Data()));
1183   TString hnameBgrXi(Form("fh2FFXi%s",strFFID.Data()));
1184
1185   TH1F* fh1NJets        = (TH1F*) list->FindObject(hnameNJets); // needed for normalization of bgr out of 2 jets
1186   TH1F* fh1FFJetPtBgr   = (TH1F*) list->FindObject(hnameJetPt);
1187   TH2F* fh2FFTrackPtBgr = (TH2F*) list->FindObject(hnameBgrTrackPt);
1188   TH2F* fh2FFZBgr       = (TH2F*) list->FindObject(hnameBgrZ);
1189   TH2F* fh2FFXiBgr      = (TH2F*) list->FindObject(hnameBgrXi);
1190
1191   if(!fh1FFJetPtBgr)  { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameJetPt.Data());      return; }
1192   if(!fh1NJets)       { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameNJets.Data());      return; }
1193   if(!fh2FFTrackPtBgr){ Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrTrackPt.Data()); return; }
1194   if(!fh2FFZBgr)      { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrZ.Data());       return; }
1195   if(!fh2FFXiBgr)     { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrXi.Data());      return; }
1196
1197   fh1FFJetPtBgr->SetDirectory(0);
1198   fh1NJets->SetDirectory(0);
1199   fh2FFTrackPtBgr->SetDirectory(0);
1200   fh2FFZBgr->SetDirectory(0);
1201   fh2FFXiBgr->SetDirectory(0);
1202
1203   f.Close();
1204
1205   // nJets per bin
1206
1207   for(Int_t i=0; i<fNJetPtSlices; i++){
1208
1209     Float_t jetPtLoLim = fJetPtSlices->At(i);
1210     Float_t jetPtUpLim = fJetPtSlices->At(i+1);
1211
1212     Int_t binLo = static_cast<Int_t>(fh1FFJetPtBgr->FindBin(jetPtLoLim));
1213     Int_t binUp = static_cast<Int_t>(fh1FFJetPtBgr->FindBin(jetPtUpLim)) - 1;
1214
1215     Float_t nJetsBin = fh1FFJetPtBgr->Integral(binLo,binUp);
1216     Double_t scaleF = 1;
1217
1218     fNJetsBgr->SetAt(nJetsBin*scaleF,i);
1219
1220     if(fDebug>0) Printf("bgr jet pt %d to %d: nJets %f, scaleF %.2f",
1221                         static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),nJetsBin,scaleF);
1222
1223   }
1224
1225   // projections: FF
1226
1227   for(Int_t i=0; i<fNJetPtSlices; i++){
1228
1229     Float_t jetPtLoLim = fJetPtSlices->At(i);
1230     Float_t jetPtUpLim = fJetPtSlices->At(i+1);
1231
1232     Int_t binLo = static_cast<Int_t>(fh2FFTrackPtBgr->GetXaxis()->FindBin(jetPtLoLim));
1233     Int_t binUp = static_cast<Int_t>(fh2FFTrackPtBgr->GetXaxis()->FindBin(jetPtUpLim))-1;
1234
1235     if(binUp > fh2FFTrackPtBgr->GetNbinsX()){
1236       Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim);
1237       return;
1238     }
1239
1240     TString strNameBgrPt(Form("fh1BgrTrackPt%s_%02d_%02d",strID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1241     TString strNameBgrZ(Form("fh1BgrZ%s_%02d_%02d",strID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1242     TString strNameBgrXi(Form("fh1BgrXi%s_%02d_%02d",strID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1243
1244     // appendix 'unbinned' to avoid histos with same name after rebinning
1245     TH1F* projPt = (TH1F*) fh2FFTrackPtBgr->ProjectionY(strNameBgrPt+"_unBinned",binLo,binUp,"o"); // option "o": original axis range
1246     TH1F* projZ  = (TH1F*) fh2FFZBgr->ProjectionY(strNameBgrZ+"_unBinned",binLo,binUp,"o");
1247     TH1F* projXi = (TH1F*) fh2FFXiBgr->ProjectionY(strNameBgrXi+"_unBinned",binLo,binUp,"o");
1248
1249     if(fNHistoBinsPt[i]) projPt = (TH1F*) projPt->Rebin(fNHistoBinsPt[i],strNameBgrPt,fHistoBinsPt[i]->GetArray());
1250     if(fNHistoBinsZ[i])  projZ  = (TH1F*) projZ->Rebin(fNHistoBinsZ[i],strNameBgrZ,fHistoBinsZ[i]->GetArray());
1251     if(fNHistoBinsXi[i]) projXi = (TH1F*) projXi->Rebin(fNHistoBinsXi[i],strNameBgrXi,fHistoBinsXi[i]->GetArray());
1252
1253     projPt->SetNameTitle(strNameBgrPt,"");
1254     projZ->SetNameTitle(strNameBgrZ,"");
1255     projXi->SetNameTitle(strNameBgrXi,"");
1256
1257     // raw bgr = corr level 0
1258     fCorrBgr[0]->AddCorrHistos(i,projPt,projZ,projXi);
1259   }
1260 }
1261
1262
1263 //__________________________________________________________________________________________________________
1264 void AliFragmentationFunctionCorrections::WriteOutput(TString strfile, TString strdir, Bool_t updateOutfile)
1265 {
1266   // write histos to file
1267   // skip histos with 0 entries
1268
1269   TString outfileOption = "RECREATE";
1270   if(updateOutfile) outfileOption = "UPDATE";
1271
1272   TFile f(strfile,outfileOption);
1273
1274   if(!f.IsOpen()){
1275     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
1276     return;
1277   }
1278
1279   if(fDebug>0) Printf("%s:%d -- write FF to file %s ",(char*)__FILE__,__LINE__,strfile.Data());
1280
1281   if(strdir && strdir.Length()){
1282     TDirectory* dir = f.mkdir(strdir);
1283     dir->cd();
1284   }
1285
1286   for(Int_t i=0; i<fNJetPtSlices; i++){
1287
1288     for(Int_t c=0; c<fNCorrectionLevels; c++) if(fCorrFF[c]->GetTrackPt(i)->GetEntries()) fCorrFF[c]->GetTrackPt(i)->Write();
1289     for(Int_t c=0; c<fNCorrectionLevels; c++) if(fCorrFF[c]->GetZ(i)->GetEntries())       fCorrFF[c]->GetZ(i)->Write();
1290     for(Int_t c=0; c<fNCorrectionLevels; c++) if(fCorrFF[c]->GetXi(i)->GetEntries())      fCorrFF[c]->GetXi(i)->Write();
1291
1292     if(fh1FFXiShift[i]->GetEntries()) fh1FFXiShift[i]->Write();
1293
1294     for(Int_t c=0; c<fNCorrectionLevelsBgr; c++) if(fCorrBgr[c]->GetTrackPt(i)->GetEntries()) fCorrBgr[c]->GetTrackPt(i)->Write();
1295     for(Int_t c=0; c<fNCorrectionLevelsBgr; c++) if(fCorrBgr[c]->GetZ(i)->GetEntries())       fCorrBgr[c]->GetZ(i)->Write();
1296     for(Int_t c=0; c<fNCorrectionLevelsBgr; c++) if(fCorrBgr[c]->GetXi(i)->GetEntries())      fCorrBgr[c]->GetXi(i)->Write();
1297
1298
1299     if(fh1FFTrackPtBackFolded[i] && fh1FFTrackPtBackFolded[i]->GetEntries()) fh1FFTrackPtBackFolded[i]->Write();
1300     if(fh1FFZBackFolded[i]       && fh1FFZBackFolded[i]->GetEntries())       fh1FFZBackFolded[i]->Write();
1301     if(fh1FFXiBackFolded[i]      && fh1FFXiBackFolded[i]->GetEntries())      fh1FFXiBackFolded[i]->Write();
1302
1303
1304     if(fh1FFRatioTrackPtFolded[i] && fh1FFRatioTrackPtFolded[i]->GetEntries()) fh1FFRatioTrackPtFolded[i]->Write();
1305     if(fh1FFRatioZFolded[i]       && fh1FFRatioZFolded[i]->GetEntries())       fh1FFRatioZFolded[i]->Write();
1306     if(fh1FFRatioXiFolded[i]      && fh1FFRatioXiFolded[i]->GetEntries())      fh1FFRatioXiFolded[i]->Write();
1307
1308     if(fh1FFRatioTrackPtBackFolded[i] && fh1FFRatioTrackPtBackFolded[i]->GetEntries()) fh1FFRatioTrackPtBackFolded[i]->Write();
1309     if(fh1FFRatioZBackFolded[i] && fh1FFRatioZBackFolded[i]->GetEntries())             fh1FFRatioZBackFolded[i]->Write();
1310     if(fh1FFRatioXiBackFolded[i] &&fh1FFRatioXiBackFolded[i]->GetEntries())            fh1FFRatioXiBackFolded[i]->Write();
1311
1312   }
1313
1314   // inclusive track pt
1315
1316   for(Int_t c=0; c<fNCorrectionLevelsSinglePt; c++) if(fCorrSinglePt[c]->GetTrackPt(0)->GetEntries()) fCorrSinglePt[c]->GetTrackPt(0)->Write();
1317   if(fh1SingleTrackPtBackFolded)      fh1SingleTrackPtBackFolded->Write();
1318   if(fh1RatioSingleTrackPtFolded)     fh1RatioSingleTrackPtFolded->Write();
1319   if(fh1RatioSingleTrackPtBackFolded) fh1RatioSingleTrackPtBackFolded->Write();
1320
1321   f.Close();
1322 }
1323
1324 //____________________________________________________________________________________________________________________________________
1325 THnSparse* AliFragmentationFunctionCorrections::TH1toSparse(const TH1F* hist, TString strName, TString strTit, const Bool_t fillConst)
1326 {
1327   // copy 1-dimensional histo to THnSparse
1328   // if fillConst TRUE, create THnSparse with same binning as hist but all entries = 1
1329   // histos with variable bin size are supported
1330
1331   // note: function returns pointer to 'new' THnSparse on heap, object needs to be deleted by user
1332
1333   THnSparseF* fhnHist;
1334
1335   Int_t nBins       = hist->GetXaxis()->GetNbins();
1336   Int_t nBinsVec[1] = { nBins };
1337
1338   const Double_t* binsVec = hist->GetXaxis()->GetXbins()->GetArray();
1339
1340   if(binsVec){ // binsVec only neq NULL if histo was rebinned before
1341
1342     fhnHist = new THnSparseF(strName,strTit,1,nBinsVec/*,binMinVec,binMaxVec*/);
1343     fhnHist->SetBinEdges(0,binsVec);
1344   }
1345   else{ // uniform bin size
1346
1347     Double_t xMin  = hist->GetXaxis()->GetXmin();
1348     Double_t xMax  = hist->GetXaxis()->GetXmax();
1349
1350     Double_t binMinVec[1] = { xMin };
1351     Double_t binMaxVec[1] = { xMax };
1352
1353     fhnHist = new THnSparseF(strName,strTit,1,nBinsVec,binMinVec,binMaxVec);
1354   }
1355
1356
1357   for(Int_t bin=1; bin<=nBins; bin++){
1358
1359     Double_t binCenter = fhnHist->GetAxis(0)->GetBinCenter(bin);
1360
1361     Double_t binCoord[] = {binCenter};
1362     fhnHist->Fill(binCoord,1); // initially need to create the bin
1363
1364     Long64_t binIndex = fhnHist->GetBin(binCoord);
1365
1366     Double_t cont =  hist->GetBinContent(bin);
1367     Double_t err  =  hist->GetBinError(bin);
1368
1369     if(fillConst){
1370       cont = 1;
1371       err  = 0;
1372     }
1373
1374     fhnHist->SetBinContent(binIndex,cont);
1375     fhnHist->SetBinError(binIndex,err);
1376  }
1377
1378  return fhnHist;
1379 }
1380
1381 //______________________________________________________________________________________________________________________________________________
1382 THnSparse* AliFragmentationFunctionCorrections::Unfold(THnSparse* hnHist, const THnSparse* hnResponse, const THnSparse* hnEff, const Int_t nIter,
1383                                                        const Bool_t useCorrelatedErrors, const THnSparse* hnPrior)
1384 {
1385   // unfold input histo
1386
1387   AliCFUnfolding unfolding("unfolding","",1,hnResponse,hnEff,hnHist,hnPrior);  // arg3: nVar; hnEff required, hnPrior not (defaults to 0x0)
1388   unfolding.SetMaxNumberOfIterations(nIter);
1389   // unfolding.SetMaxChi2PerDOF(1.e-07); // OBSOLETE !!!
1390   // if(useSmoothing) unfolding.UseSmoothing();
1391
1392   if(useCorrelatedErrors) unfolding.SetUseCorrelatedErrors();
1393   unfolding.Unfold();
1394
1395   THnSparse* unfolded = unfolding.GetUnfolded();
1396
1397   TString hnameUnf = hnHist->GetName();
1398   hnameUnf.Append("_unf");
1399   unfolded->SetNameTitle(hnameUnf,hnHist->GetTitle());
1400
1401   return unfolded;
1402 }
1403
1404 //___________________________________________________________________________________________________________________________
1405 void AliFragmentationFunctionCorrections::UnfoldHistos(const Int_t nIter, const Bool_t useCorrelatedErrors, const Int_t type)
1406 {
1407   // loop over jet pt slices and unfold dN/dpt spectra
1408
1409   TString labelCorr = fCorrFF[fNCorrectionLevels-1]->GetLabel();
1410   if(!labelCorr.Contains("unfold")) AddCorrectionLevel("unfold");
1411
1412   for(Int_t i=0; i<fNJetPtSlices; i++){
1413
1414     TH1F* hist = 0;
1415     if(type == kFlagPt)      hist = fCorrFF[fNCorrectionLevels-2]->GetTrackPt(i); // level -2: before unfolding, level -1: unfolded
1416     else if(type == kFlagZ)  hist = fCorrFF[fNCorrectionLevels-2]->GetZ(i);       // level -2: before unfolding, level -1: unfolded
1417     else if(type == kFlagXi) hist = fCorrFF[fNCorrectionLevels-2]->GetXi(i);      // level -2: before unfolding, level -1: unfolded
1418     else{
1419       Printf("%s%d unknown type",(char*)__FILE__,__LINE__);
1420       return;
1421     }
1422
1423     THnSparse* hnResponse = 0;
1424     if(type == kFlagPt) hnResponse = fhnResponsePt[i];
1425     else if(type == kFlagZ)  hnResponse = fhnResponseZ[i];
1426     else if(type == kFlagXi) hnResponse = fhnResponseXi[i];
1427     else{
1428       Printf("%s%d unknown type",(char*)__FILE__,__LINE__);
1429       return;
1430     }
1431
1432
1433     TH1F* hPrior = 0;
1434     if(type == kFlagPt && fh1FFTrackPtPrior[i]  && ((TString(fh1FFTrackPtPrior[i]->GetName())).Length() > 0) ) hPrior = fh1FFTrackPtPrior[i];
1435     else if(type == kFlagZ  && fh1FFZPrior[i]        && ((TString(fh1FFZPrior[i]->GetName())).Length() > 0)       ) hPrior = fh1FFZPrior[i];
1436     else if(type == kFlagXi && fh1FFXiPrior[i]       && ((TString(fh1FFXiPrior[i]->GetName())).Length() > 0)      ) hPrior = fh1FFXiPrior[i];
1437     else{
1438       Printf("%s%d unknown type",(char*)__FILE__,__LINE__);
1439       return;
1440     }
1441
1442
1443     TString histNameTHn = hist->GetName();
1444     histNameTHn.ReplaceAll("TH1","THn");
1445
1446     TString priorNameTHn;
1447     if(hPrior){
1448       priorNameTHn = hPrior->GetName();
1449       priorNameTHn.ReplaceAll("TH1","THn");
1450     }
1451
1452     TString histNameBackFolded = hist->GetName();
1453     histNameBackFolded.Append("_backfold");
1454
1455     TString histNameRatioFolded = hist->GetName();
1456     histNameRatioFolded.ReplaceAll("fh1FF","hRatioFF");
1457     histNameRatioFolded.Append("_unfold");
1458
1459     TString histNameRatioBackFolded = hist->GetName();
1460     histNameRatioBackFolded.ReplaceAll("fh1FF","hRatioFF");
1461     histNameRatioBackFolded.Append("_backfold");
1462
1463     THnSparse* hnHist           = TH1toSparse(hist,histNameTHn,hist->GetTitle());
1464     THnSparse* hnFlatEfficiency = TH1toSparse(hist,"fhnEfficiency","eff",kTRUE); // could optionally also use real eff
1465     THnSparse* hnPrior          = 0;
1466     if(hPrior) hnPrior = TH1toSparse(hPrior,priorNameTHn,hPrior->GetTitle());
1467
1468     THnSparse* hnUnfolded
1469       = Unfold(hnHist,hnResponse,hnFlatEfficiency,nIter,useCorrelatedErrors,hnPrior);
1470
1471     TH1F* hUnfolded = (TH1F*) hnUnfolded->Projection(0);
1472     hUnfolded->SetNameTitle(hist->GetName(),hist->GetTitle());
1473
1474     if(type == kFlagPt) fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,hUnfolded,0,0);
1475     if(type == kFlagZ)  fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,0,hUnfolded,0);
1476     if(type == kFlagXi) fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,0,0,hUnfolded);
1477
1478     // backfolding: apply response matrix to unfolded spectrum
1479     TH1F* hBackFolded = ApplyResponse(hUnfolded,hnResponse);
1480     hBackFolded->SetNameTitle(histNameBackFolded,hist->GetTitle());
1481
1482     if(type == kFlagPt) fh1FFTrackPtBackFolded[i] = hBackFolded;
1483     if(type == kFlagZ)  fh1FFZBackFolded[i]       = hBackFolded;
1484     if(type == kFlagXi) fh1FFXiBackFolded[i]      = hBackFolded;
1485
1486     // ratio unfolded to original histo
1487     TH1F* hRatioUnfolded = (TH1F*) hUnfolded->Clone(histNameRatioFolded);
1488     hRatioUnfolded->Reset();
1489     hRatioUnfolded->Divide(hUnfolded,hist,1,1,"B");
1490
1491     if(type == kFlagPt) fh1FFRatioTrackPtFolded[i] = hRatioUnfolded;
1492     if(type == kFlagZ)  fh1FFRatioZFolded[i]       = hRatioUnfolded;
1493     if(type == kFlagXi) fh1FFRatioXiFolded[i]      = hRatioUnfolded;
1494
1495
1496     // ratio backfolded to original histo
1497     TH1F* hRatioBackFolded = (TH1F*) hBackFolded->Clone(histNameRatioBackFolded);
1498     hRatioBackFolded->Reset();
1499     hRatioBackFolded->Divide(hBackFolded,hist,1,1,"B");
1500
1501     if(type == kFlagPt) fh1FFRatioTrackPtBackFolded[i] = hRatioBackFolded;
1502     if(type == kFlagZ)  fh1FFRatioZBackFolded[i]       = hRatioBackFolded;
1503     if(type == kFlagXi) fh1FFRatioXiBackFolded[i]      = hRatioBackFolded;
1504
1505     delete hnHist;
1506     delete hnFlatEfficiency;
1507
1508  }
1509 }
1510
1511 //_____________________________________________________________________________________________________
1512 void AliFragmentationFunctionCorrections::UnfoldPt(const Int_t nIter, const Bool_t useCorrelatedErrors)
1513 {
1514
1515   Int_t type = kFlagPt;
1516   UnfoldHistos(nIter, useCorrelatedErrors, type);
1517 }
1518
1519 //_____________________________________________________________________________________________________
1520 void AliFragmentationFunctionCorrections::UnfoldZ(const Int_t nIter, const Bool_t useCorrelatedErrors)
1521 {
1522
1523   Int_t type = kFlagZ;
1524   UnfoldHistos(nIter, useCorrelatedErrors, type);
1525 }
1526
1527 //_____________________________________________________________________________________________________
1528 void AliFragmentationFunctionCorrections::UnfoldXi(const Int_t nIter, const Bool_t useCorrelatedErrors)
1529 {
1530
1531   Int_t type = kFlagXi;
1532   UnfoldHistos(nIter, useCorrelatedErrors, type);
1533 }
1534
1535 //______________________________________________________________________________________________
1536 TH1F* AliFragmentationFunctionCorrections::ApplyResponse(const TH1F* hist, THnSparse* hnResponse)
1537 {
1538   // apply (multiply) response matrix to hist
1539
1540   TH1F* hOut = new TH1F(*hist);
1541   hOut->Reset();
1542
1543   const Int_t axisM = 0;
1544   const Int_t axisT = 1;
1545
1546   Int_t nBinsM = hnResponse->GetAxis(axisM)->GetNbins();
1547   Int_t nBinsT = hnResponse->GetAxis(axisT)->GetNbins();
1548
1549   // response matrix normalization
1550   // do it in this function and not when reading response, since for 'proper' normalization errors are difficult to assign
1551   // so for unfolding proper we leave it to CORRFW ...
1552
1553   Double_t normFacResponse[nBinsT];
1554
1555   for(Int_t iT=1; iT<=nBinsT; iT++){
1556
1557     Double_t sumResp = 0;
1558
1559     for(Int_t iM=1; iM<=nBinsM; iM++){
1560
1561       Double_t coordM = hnResponse->GetAxis(axisM)->GetBinCenter(iM);
1562       Double_t coordT = hnResponse->GetAxis(axisT)->GetBinCenter(iT);
1563
1564       Double_t binCoord[] = {coordM,coordT};
1565
1566       Long64_t binIndex = hnResponse->GetBin(binCoord);
1567
1568       Double_t resp = hnResponse->GetBinContent(binIndex);
1569
1570       sumResp += resp;
1571     }
1572
1573     normFacResponse[iT] = 0;
1574     if(sumResp) normFacResponse[iT] = 1/sumResp;
1575   }
1576
1577
1578
1579   for(Int_t iM=1; iM<=nBinsM; iM++){
1580
1581     Double_t contM   = 0;
1582
1583     for(Int_t iT=1; iT<=nBinsT; iT++){
1584
1585       Double_t contT = hist->GetBinContent(iT);
1586
1587       Double_t coordM = hnResponse->GetAxis(axisM)->GetBinCenter(iM);
1588       Double_t coordT = hnResponse->GetAxis(axisT)->GetBinCenter(iT);
1589
1590       Double_t binCoord[] = {coordM,coordT};
1591
1592       Long64_t binIndex = hnResponse->GetBin(binCoord);
1593
1594       Double_t resp = hnResponse->GetBinContent(binIndex);
1595
1596       contM   += resp*normFacResponse[iT]*contT;
1597     }
1598
1599     hOut->SetBinContent(iM,contM);
1600   }
1601
1602   return hOut;
1603 }
1604
1605 //_______________________________________________________________________________________________________
1606 void AliFragmentationFunctionCorrections::ReadEfficiency(TString strfile, TString strdir, TString strlist)
1607 {
1608   // read reconstruction efficiency from file
1609   // argument strlist optional - read from directory strdir if not specified
1610
1611   // temporary histos to hold histos from file
1612   TH1F* hEffPt[fNJetPtSlices];
1613   TH1F* hEffZ[fNJetPtSlices];
1614   TH1F* hEffXi[fNJetPtSlices];
1615
1616   for(Int_t i=0; i<fNJetPtSlices; i++) hEffPt[i] = 0;
1617   for(Int_t i=0; i<fNJetPtSlices; i++) hEffZ[i]  = 0;
1618   for(Int_t i=0; i<fNJetPtSlices; i++) hEffXi[i] = 0;
1619
1620   TFile f(strfile,"READ");
1621
1622   if(!f.IsOpen()){
1623     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
1624     return;
1625   }
1626
1627   if(fDebug>0) Printf("%s:%d -- read efficiencies from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
1628
1629   if(strdir && strdir.Length()) gDirectory->cd(strdir);
1630
1631   TList* list = 0;
1632
1633   if(strlist && strlist.Length()){
1634
1635     if(!(list = (TList*) gDirectory->Get(strlist))){
1636       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
1637       return;
1638     }
1639   }
1640
1641   for(Int_t i=0; i<fNJetPtSlices; i++){
1642
1643     Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
1644     Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
1645
1646     TString strNameEffPt(Form("hEffPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
1647     TString strNameEffZ(Form("hEffZ_%02d_%02d",jetPtLoLim,jetPtUpLim));
1648     TString strNameEffXi(Form("hEffXi_%02d_%02d",jetPtLoLim,jetPtUpLim));
1649
1650
1651     if(list){
1652       hEffPt[i] = (TH1F*) list->FindObject(strNameEffPt);
1653       hEffZ[i]  = (TH1F*) list->FindObject(strNameEffZ);
1654       hEffXi[i] = (TH1F*) list->FindObject(strNameEffXi);
1655     }
1656     else{
1657       hEffPt[i] = (TH1F*) gDirectory->Get(strNameEffPt);
1658       hEffZ[i]  = (TH1F*) gDirectory->Get(strNameEffZ);
1659       hEffXi[i] = (TH1F*) gDirectory->Get(strNameEffXi);
1660     }
1661
1662     if(!hEffPt[i]){
1663       Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffPt.Data());
1664     }
1665
1666     if(!hEffZ[i]){
1667       Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffZ.Data());
1668     }
1669
1670     if(!hEffXi[i]){
1671       Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffXi.Data());
1672     }
1673
1674
1675     if(fNHistoBinsPt[i]) hEffPt[i] = (TH1F*) hEffPt[i]->Rebin(fNHistoBinsPt[i],strNameEffPt+"_rebin",fHistoBinsPt[i]->GetArray());
1676     if(fNHistoBinsZ[i])  hEffZ[i]  = (TH1F*) hEffZ[i]->Rebin(fNHistoBinsZ[i],strNameEffZ+"_rebin",fHistoBinsZ[i]->GetArray());
1677     if(fNHistoBinsXi[i]) hEffXi[i] = (TH1F*) hEffXi[i]->Rebin(fNHistoBinsXi[i],strNameEffXi+"_rebin",fHistoBinsXi[i]->GetArray());
1678
1679     if(hEffPt[i]) hEffPt[i]->SetDirectory(0);
1680     if(hEffZ[i])  hEffZ[i]->SetDirectory(0);
1681     if(hEffXi[i]) hEffXi[i]->SetDirectory(0);
1682
1683   } // jet slices loop
1684
1685   f.Close();
1686
1687   for(Int_t i=0; i<fNJetPtSlices; i++){ // 2nd loop: need to close input file before placing histos
1688     if(hEffPt[i]) new(fh1EffPt[i]) TH1F(*hEffPt[i]);
1689     if(hEffZ[i])  new(fh1EffZ[i])  TH1F(*hEffZ[i]);
1690     if(hEffXi[i]) new(fh1EffXi[i]) TH1F(*hEffXi[i]);
1691   }
1692 }
1693
1694 //___________________________________________________________________________________________________________
1695 void AliFragmentationFunctionCorrections::ReadBgrEfficiency(TString strfile, TString strdir, TString strlist)
1696 {
1697   // read bgr eff from file
1698   // argument strlist optional - read from directory strdir if not specified
1699
1700   TH1F* hEffPtBgr[fNJetPtSlices];
1701   TH1F* hEffZBgr [fNJetPtSlices];
1702   TH1F* hEffXiBgr[fNJetPtSlices];
1703
1704   for(Int_t i=0; i<fNJetPtSlices; i++) hEffPtBgr[i] = 0;
1705   for(Int_t i=0; i<fNJetPtSlices; i++) hEffZBgr[i]  = 0;
1706   for(Int_t i=0; i<fNJetPtSlices; i++) hEffXiBgr[i] = 0;
1707
1708
1709   TFile f(strfile,"READ");
1710
1711   if(!f.IsOpen()){
1712     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
1713     return;
1714   }
1715
1716   if(fDebug>0) Printf("%s:%d -- read bgr efficiencies from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
1717
1718   if(strdir && strdir.Length()) gDirectory->cd(strdir);
1719
1720   TList* list = 0;
1721
1722   if(strlist && strlist.Length()){
1723
1724     if(!(list = (TList*) gDirectory->Get(strlist))){
1725       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
1726       return;
1727     }
1728   }
1729
1730   for(Int_t i=0; i<fNJetPtSlices; i++){
1731
1732     Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
1733     Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
1734
1735     TString strNameEffPtBgr(Form("hEffPtBgr%02dto%02d",jetPtLoLim,jetPtUpLim));
1736     TString strNameEffZBgr(Form("hEffZBgr%02dto%02d",jetPtLoLim,jetPtUpLim));
1737     TString strNameEffXiBgr(Form("hEffXiBgr%02dto%02d",jetPtLoLim,jetPtUpLim));
1738
1739
1740     if(list){
1741       hEffPtBgr[i] = (TH1F*) list->FindObject(strNameEffPtBgr);
1742       hEffZBgr[i]  = (TH1F*) list->FindObject(strNameEffZBgr);
1743       hEffXiBgr[i] = (TH1F*) list->FindObject(strNameEffXiBgr);
1744     }
1745     else{
1746       hEffPtBgr[i] = (TH1F*) gDirectory->Get(strNameEffPtBgr);
1747       hEffZBgr[i]  = (TH1F*) gDirectory->Get(strNameEffZBgr);
1748       hEffXiBgr[i] = (TH1F*) gDirectory->Get(strNameEffXiBgr);
1749     }
1750
1751     if(!hEffPtBgr[i]){
1752       Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffPtBgr.Data());
1753     }
1754
1755     if(!hEffZBgr[i]){
1756       Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffZBgr.Data());
1757     }
1758
1759     if(!hEffXiBgr[i]){
1760       Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffXiBgr.Data());
1761     }
1762
1763
1764     if(fNHistoBinsPt[i]) hEffPtBgr[i] = (TH1F*) hEffPtBgr[i]->Rebin(fNHistoBinsPt[i],strNameEffPtBgr+"_rebin",fHistoBinsPt[i]->GetArray());
1765     if(fNHistoBinsZ[i])  hEffZBgr[i]  = (TH1F*) hEffZBgr[i]->Rebin(fNHistoBinsZ[i],strNameEffZBgr+"_rebin",fHistoBinsZ[i]->GetArray());
1766     if(fNHistoBinsXi[i]) hEffXiBgr[i] = (TH1F*) hEffXiBgr[i]->Rebin(fNHistoBinsXi[i],strNameEffXiBgr+"_rebin",fHistoBinsXi[i]->GetArray());
1767
1768     if(hEffPtBgr[i]) hEffPtBgr[i]->SetDirectory(0);
1769     if(hEffZBgr[i])  hEffZBgr[i]->SetDirectory(0);
1770     if(hEffXiBgr[i]) hEffXiBgr[i]->SetDirectory(0);
1771
1772   } // jet slices loop
1773
1774   f.Close();
1775
1776   for(Int_t i=0; i<fNJetPtSlices; i++){ // 2nd loop: need to close input file before placing histos
1777     if(hEffPtBgr[i]) new(fh1EffBgrPt[i]) TH1F(*hEffPtBgr[i]);
1778     if(hEffZBgr[i])  new(fh1EffBgrZ[i])  TH1F(*hEffZBgr[i]);
1779     if(hEffXiBgr[i]) new(fh1EffBgrXi[i]) TH1F(*hEffXiBgr[i]);
1780   }
1781 }
1782
1783 // ________________________________________________
1784 void AliFragmentationFunctionCorrections::EffCorr()
1785 {
1786   // apply efficiency correction
1787
1788   AddCorrectionLevel("eff");
1789
1790   for(Int_t i=0; i<fNJetPtSlices; i++){
1791
1792     TH1F* histPt = fCorrFF[fNCorrectionLevels-2]->GetTrackPt(i);
1793     TH1F* histZ  = fCorrFF[fNCorrectionLevels-2]->GetZ(i);
1794     TH1F* histXi = fCorrFF[fNCorrectionLevels-2]->GetXi(i);
1795
1796     TString histNamePt = histPt->GetName();
1797     TString histNameZ  = histZ->GetName();
1798     TString histNameXi = histXi->GetName();
1799
1800
1801     TH1F* hFFTrackPtEffCorr = (TH1F*) histPt->Clone(histNamePt);
1802     hFFTrackPtEffCorr->Divide(histPt,fh1EffPt[i],1,1,"");
1803
1804     TH1F* hFFZEffCorr = (TH1F*) histZ->Clone(histNameZ);
1805     hFFZEffCorr->Divide(histZ,fh1EffZ[i],1,1,"");
1806
1807     TH1F* hFFXiEffCorr = (TH1F*) histXi->Clone(histNameXi);
1808     hFFXiEffCorr->Divide(histXi,fh1EffXi[i],1,1,"");
1809
1810     fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,hFFTrackPtEffCorr,hFFZEffCorr,hFFXiEffCorr);
1811   }
1812 }
1813
1814 //___________________________________________________
1815 void AliFragmentationFunctionCorrections::EffCorrBgr()
1816 {
1817   // apply efficiency correction to bgr distributions
1818
1819   AddCorrectionLevelBgr("eff");
1820
1821   Printf("%s:%d -- apply efficiency correction, corrLevel %d",(char*)__FILE__,__LINE__,fNCorrectionLevels-1);
1822
1823   for(Int_t i=0; i<fNJetPtSlices; i++){
1824
1825     TH1F* histPt = fCorrBgr[fNCorrectionLevelsBgr-2]->GetTrackPt(i);
1826     TH1F* histZ  = fCorrBgr[fNCorrectionLevelsBgr-2]->GetZ(i);
1827     TH1F* histXi = fCorrBgr[fNCorrectionLevelsBgr-2]->GetXi(i);
1828
1829     TString histNamePt = histPt->GetName();
1830     TString histNameZ  = histZ->GetName();
1831     TString histNameXi = histXi->GetName();
1832
1833
1834     TH1F* hFFTrackPtEffCorr = (TH1F*) histPt->Clone(histNamePt);
1835     hFFTrackPtEffCorr->Divide(histPt,fh1EffPt[i],1,1,"");
1836
1837     TH1F* hFFZEffCorr = (TH1F*) histZ->Clone(histNameZ);
1838     hFFZEffCorr->Divide(histZ,fh1EffZ[i],1,1,"");
1839
1840     TH1F* hFFXiEffCorr = (TH1F*) histXi->Clone(histNameXi);
1841     hFFXiEffCorr->Divide(histXi,fh1EffXi[i],1,1,"");
1842
1843     fCorrBgr[fNCorrectionLevelsBgr-1]->AddCorrHistos(i,hFFTrackPtEffCorr,hFFZEffCorr,hFFXiEffCorr);
1844   }
1845 }
1846
1847 //______________________________________________________________________
1848 void AliFragmentationFunctionCorrections::XiShift(const Int_t corrLevel)
1849 {
1850   // re-evaluate jet energy after FF corrections from dN/dpt distribution
1851   // apply correction (shift) to dN/dxi distribution: xi = ln (pt/E) -> xi' = ln (pt/E') = ln (pt/E x E/E') = xi + ln E/E'
1852   // argument corrlevel: which level of correction to be corrected/shifted to
1853
1854   if(corrLevel>=fNCorrectionLevels){
1855     Printf(" calc xi shift: corrLevel exceeded - do nothing");
1856     return;
1857   }
1858
1859   Double_t* jetPtUncorr = new Double_t[fNJetPtSlices];
1860   Double_t* jetPtCorr   = new Double_t[fNJetPtSlices];
1861   Double_t* deltaXi     = new Double_t[fNJetPtSlices];
1862
1863   for(Int_t i=0; i<fNJetPtSlices; i++){
1864
1865     TH1F* histPtRaw = fCorrFF[0]->GetTrackPt(i);
1866     TH1F* histPt    = fCorrFF[corrLevel]->GetTrackPt(i);
1867
1868     Double_t ptUncorr = 0;
1869     Double_t ptCorr   = 0;
1870
1871     for(Int_t bin = 1; bin<=histPtRaw->GetNbinsX(); bin++){
1872
1873       Double_t cont   = histPtRaw->GetBinContent(bin);
1874       Double_t width  = histPtRaw->GetBinWidth(bin);
1875       Double_t meanPt = histPtRaw->GetBinCenter(bin);
1876
1877       ptUncorr += meanPt*cont*width;
1878     }
1879
1880     for(Int_t bin = 1; bin<=histPt->GetNbinsX(); bin++){
1881
1882       Double_t cont   = histPt->GetBinContent(bin);
1883       Double_t width  = histPt->GetBinWidth(bin);
1884       Double_t meanPt = histPt->GetBinCenter(bin);
1885
1886       ptCorr += meanPt*cont*width;
1887     }
1888
1889     jetPtUncorr[i] = ptUncorr;
1890     jetPtCorr[i]   = ptCorr;
1891   }
1892
1893   // calc dXi from dN/dpt distribution :
1894   // sum over track pt for raw and corrected FF is equivalent to raw/corrected jet pt
1895
1896   for(Int_t i=0; i<fNJetPtSlices; i++){
1897
1898     Float_t jetPtLoLim = fJetPtSlices->At(i);
1899     Float_t jetPtUpLim = fJetPtSlices->At(i+1);
1900
1901     Double_t meanJetPt = 0.5*(jetPtUpLim+jetPtLoLim);
1902
1903     Double_t ptUncorr = jetPtUncorr[i];
1904     Double_t ptCorr   = jetPtCorr[i];
1905
1906     Double_t dXi = TMath::Log(ptCorr/ptUncorr);
1907
1908     Printf(" calc xi shift: jet pt slice %d, mean jet pt %f, ptUncorr %f, ptCorr %f, ratio corr/uncorr %f, dXi %f "
1909            ,i,meanJetPt,ptUncorr,ptCorr,ptCorr/ptUncorr,dXi);
1910
1911     deltaXi[i] = dXi;
1912   }
1913
1914   // book & fill new dN/dxi histos
1915
1916   for(Int_t i=0; i<fNJetPtSlices; i++){
1917
1918     TH1F* histXi = fCorrFF[corrLevel]->GetXi(i);
1919
1920     Double_t dXi = deltaXi[i];
1921
1922     Int_t nBins  = histXi->GetNbinsX();
1923     const Double_t* binsVec = histXi->GetXaxis()->GetXbins()->GetArray();
1924     Float_t binsVecNew[nBins+1];
1925
1926     TString strName = histXi->GetName();
1927     strName.Append("_shift");
1928     TString strTit  = histXi->GetTitle();
1929
1930     TH1F* hXiCorr;
1931
1932     // create shifted histo ...
1933
1934     if(binsVec){ // binsVec only neq NULL if histo was rebinned before
1935
1936       for(Int_t bin=0; bin<nBins+1; bin++) binsVecNew[bin] = binsVec[bin] + dXi;
1937       hXiCorr = new TH1F(strName,strTit,nBins,binsVecNew);
1938     }
1939     else{ // uniform bin size
1940
1941       Double_t xMin  = histXi->GetXaxis()->GetXmin();
1942       Double_t xMax  = histXi->GetXaxis()->GetXmax();
1943
1944       xMin += dXi;
1945       xMax += dXi;
1946
1947       hXiCorr = new TH1F(strName,strTit,nBins,xMin,xMax);
1948     }
1949
1950     // ... and fill
1951
1952     for(Int_t bin=1; bin<nBins+1; bin++){
1953       Double_t cont = histXi->GetBinContent(bin);
1954       Double_t err  = histXi->GetBinError(bin);
1955
1956       hXiCorr->SetBinContent(bin,cont);
1957       hXiCorr->SetBinError(bin,err);
1958     }
1959
1960     new(fh1FFXiShift[i]) TH1F(*hXiCorr);
1961     delete hXiCorr;
1962   }
1963
1964   delete[] jetPtUncorr;
1965   delete[] jetPtCorr;
1966   delete[] deltaXi;
1967 }
1968
1969 //_____________________________________________________
1970 void AliFragmentationFunctionCorrections::SubtractBgr()
1971 {
1972   // subtract bgr distribution from FF
1973   // the current corr level is used for both
1974
1975   AddCorrectionLevel("bgrSub");
1976
1977   for(Int_t i=0; i<fNJetPtSlices; i++){
1978
1979     TH1F* histPt = fCorrFF[fNCorrectionLevels-2]->GetTrackPt(i);
1980     TH1F* histZ  = fCorrFF[fNCorrectionLevels-2]->GetZ(i);
1981     TH1F* histXi = fCorrFF[fNCorrectionLevels-2]->GetXi(i);
1982
1983     TH1F* histPtBgr = fCorrBgr[fNCorrectionLevelsBgr-1]->GetTrackPt(i);
1984     TH1F* histZBgr  = fCorrBgr[fNCorrectionLevelsBgr-1]->GetZ(i);
1985     TH1F* histXiBgr = fCorrBgr[fNCorrectionLevelsBgr-1]->GetXi(i);
1986
1987     TString histNamePt = histPt->GetName();
1988     TString histNameZ  = histZ->GetName();
1989     TString histNameXi = histXi->GetName();
1990
1991
1992     TH1F* hFFTrackPtBgrSub = (TH1F*) histPt->Clone(histNamePt);
1993     hFFTrackPtBgrSub->Add(histPtBgr,-1);
1994
1995     TH1F* hFFZBgrSub =  (TH1F*) histZ->Clone(histNameZ);
1996     hFFZBgrSub->Add(histZBgr,-1);
1997
1998     TH1F* hFFXiBgrSub = (TH1F*) histXi->Clone(histNameXi);
1999     hFFXiBgrSub->Add(histXiBgr,-1);
2000
2001     fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,hFFTrackPtBgrSub,hFFZBgrSub,hFFXiBgrSub);
2002   }
2003 }
2004
2005 //________________________________________________________________________________________________________________
2006 void AliFragmentationFunctionCorrections::WriteSingleTrackEff(TString strInfile, TString strID, TString strOutfile,
2007                                                               Bool_t updateOutfile, TString strOutDir,TString strPostfix)
2008 {
2009   // read task ouput from MC and write single track eff - standard dir/list
2010
2011   TString strdir  = "PWG4_FragmentationFunction_" + strID;
2012   TString strlist = "fracfunc_" + strID;
2013
2014   WriteSingleTrackEff(strInfile,strdir,strlist,strOutfile,updateOutfile,strOutDir,strPostfix);
2015 }
2016
2017 //___________________________________________________________________________________________________________________________________
2018 void AliFragmentationFunctionCorrections::WriteSingleTrackEff(TString strInfile, TString strdir, TString strlist,
2019                                                               TString strOutfile, Bool_t updateOutfile, TString strOutDir,TString strPostfix)
2020 {
2021   // read task output from MC and write single track eff as function of pt, eta, phi
2022   // argument strlist optional - read from directory strdir if not specified
2023   // write eff to file stroutfile - by default only 'update' file (don't overwrite)
2024
2025
2026   TH1D* hdNdptTracksMCPrimGen;
2027   TH2D* hdNdetadphiTracksMCPrimGen;
2028
2029   TH1D* hdNdptTracksMCPrimRec;
2030   TH2D* hdNdetadphiTracksMCPrimRec;
2031
2032
2033   TFile f(strInfile,"READ");
2034
2035   if(!f.IsOpen()){
2036     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2037     return;
2038   }
2039
2040   if(fDebug>0) Printf("%s:%d -- writeSingleTrackEff: open task ouput file %s ",(char*)__FILE__,__LINE__,strInfile.Data());
2041
2042   if(strdir && strdir.Length()){
2043     if(fDebug>0) Printf("%s:%d -- writeSingleTrackEff: change dir to %s",(char*)__FILE__,__LINE__,strdir.Data());
2044     gDirectory->cd(strdir);
2045   }
2046
2047   TList* list = 0;
2048
2049   if(strlist && strlist.Length()){
2050
2051     if(!(list = (TList*) gDirectory->Get(strlist))){
2052       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2053       return;
2054     }
2055   }
2056
2057
2058   TString hnamePtRecEffGen = "fh1TrackQAPtRecEffGen";
2059   if(strPostfix.Length()) hnamePtRecEffGen.Form("fh1TrackQAPtRecEffGen_%s",strPostfix.Data());
2060
2061   TString hnamePtRecEffRec = "fh1TrackQAPtRecEffRec";
2062   if(strPostfix.Length()) hnamePtRecEffRec.Form("fh1TrackQAPtRecEffRec_%s",strPostfix.Data());
2063
2064   TString hnameEtaPhiRecEffGen = "fh2TrackQAEtaPhiRecEffGen";
2065   if(strPostfix.Length()) hnameEtaPhiRecEffGen.Form("fh2TrackQAEtaPhiRecEffGen_%s",strPostfix.Data());
2066
2067   TString hnameEtaPhiRecEffRec = "fh2TrackQAEtaPhiRecEffRec";
2068   if(strPostfix.Length()) hnameEtaPhiRecEffRec.Form("fh2TrackQAEtaPhiRecEffRec_%s",strPostfix.Data());
2069
2070
2071   if(list){
2072     hdNdptTracksMCPrimGen       = (TH1D*) list->FindObject(hnamePtRecEffGen);
2073     hdNdetadphiTracksMCPrimGen  = (TH2D*) list->FindObject(hnameEtaPhiRecEffGen);
2074
2075     hdNdptTracksMCPrimRec       = (TH1D*) list->FindObject(hnamePtRecEffRec);
2076     hdNdetadphiTracksMCPrimRec  = (TH2D*) list->FindObject(hnameEtaPhiRecEffRec);
2077   }
2078   else{
2079     hdNdptTracksMCPrimGen       = (TH1D*) gDirectory->Get(hnamePtRecEffGen);
2080     hdNdetadphiTracksMCPrimGen  = (TH2D*) gDirectory->Get(hnameEtaPhiRecEffGen);
2081
2082     hdNdptTracksMCPrimRec       = (TH1D*) gDirectory->Get(hnamePtRecEffRec);
2083     hdNdetadphiTracksMCPrimRec  = (TH2D*) gDirectory->Get(hnameEtaPhiRecEffRec);
2084   }
2085
2086   hdNdptTracksMCPrimGen->SetDirectory(0);
2087   hdNdetadphiTracksMCPrimGen->SetDirectory(0);
2088   hdNdptTracksMCPrimRec->SetDirectory(0);
2089   hdNdetadphiTracksMCPrimRec->SetDirectory(0);
2090
2091   f.Close();
2092
2093   // projections: dN/deta, dN/dphi
2094
2095   TH1D* hdNdetaTracksMCPrimGen = (TH1D*) hdNdetadphiTracksMCPrimGen->ProjectionX("hdNdetaTracksMcPrimGen");
2096   TH1D* hdNdphiTracksMCPrimGen = (TH1D*) hdNdetadphiTracksMCPrimGen->ProjectionY("hdNdphiTracksMcPrimGen");
2097
2098   TH1D* hdNdetaTracksMCPrimRec = (TH1D*) hdNdetadphiTracksMCPrimRec->ProjectionX("hdNdetaTracksMcPrimRec");
2099   TH1D* hdNdphiTracksMCPrimRec = (TH1D*) hdNdetadphiTracksMCPrimRec->ProjectionY("hdNdphiTracksMcPrimRec");
2100
2101   // rebin
2102
2103   TString strNamePtGen = "hTrackPtGenPrim";
2104   TString strNamePtRec = "hTrackPtRecPrim";
2105
2106   if(fNHistoBinsSinglePt) hdNdptTracksMCPrimGen = (TH1D*) hdNdptTracksMCPrimGen->Rebin(fNHistoBinsSinglePt,strNamePtGen,fHistoBinsSinglePt->GetArray());
2107   if(fNHistoBinsSinglePt) hdNdptTracksMCPrimRec = (TH1D*) hdNdptTracksMCPrimRec->Rebin(fNHistoBinsSinglePt,strNamePtRec,fHistoBinsSinglePt->GetArray());
2108
2109     // efficiency: divide
2110
2111   TString hNameTrackEffPt = "hSingleTrackEffPt";
2112   if(strPostfix.Length()) hNameTrackEffPt.Form("hSingleTrackEffPt_%s",strPostfix.Data());
2113
2114   TString hNameTrackEffEta = "hSingleTrackEffEta";
2115   if(strPostfix.Length()) hNameTrackEffEta.Form("hSingleTrackEffEta_%s",strPostfix.Data());
2116
2117   TString hNameTrackEffPhi = "hSingleTrackEffPhi";
2118   if(strPostfix.Length()) hNameTrackEffPhi.Form("hSingleTrackEffPhi_%s",strPostfix.Data());
2119
2120
2121   TH1F* hSingleTrackEffPt = (TH1F*) hdNdptTracksMCPrimRec->Clone(hNameTrackEffPt);
2122   hSingleTrackEffPt->Divide(hdNdptTracksMCPrimRec,hdNdptTracksMCPrimGen,1,1,"B"); // binominal errors
2123
2124   TH1F* hSingleTrackEffEta = (TH1F*) hdNdetaTracksMCPrimRec->Clone(hNameTrackEffEta);
2125   hSingleTrackEffEta->Divide(hdNdetaTracksMCPrimRec,hdNdetaTracksMCPrimGen,1,1,"B"); // binominal errors
2126
2127   TH1F* hSingleTrackEffPhi = (TH1F*) hdNdphiTracksMCPrimRec->Clone(hNameTrackEffPhi);
2128   hSingleTrackEffPhi->Divide(hdNdphiTracksMCPrimRec,hdNdphiTracksMCPrimGen,1,1,"B"); // binominal errors
2129
2130
2131   TString outfileOption = "RECREATE";
2132   if(updateOutfile)  outfileOption = "UPDATE";
2133
2134   TFile out(strOutfile,outfileOption);
2135
2136   if(!out.IsOpen()){
2137     Printf("%s:%d -- error opening efficiency output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
2138     return;
2139   }
2140
2141   if(fDebug>0) Printf("%s:%d -- write efficiency to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
2142
2143   if(strOutDir && strOutDir.Length()){
2144
2145     TDirectory* dir;
2146     if((dir = ((TDirectory*) gDirectory->Get(strOutDir)))) dir->cd();
2147     else{
2148       dir = out.mkdir(strOutDir);
2149       dir->cd();
2150     }
2151   }
2152
2153   hSingleTrackEffPt->Write();
2154   hSingleTrackEffEta->Write();
2155   hSingleTrackEffPhi->Write();
2156
2157   out.Close();
2158 }
2159
2160 //________________________________________________________________________________________________________________
2161 void AliFragmentationFunctionCorrections::WriteSingleTrackSecCorr(TString strInfile, TString strID, TString strOutfile,
2162                                                                   Bool_t updateOutfile, TString strOutDir)
2163 {
2164   // read task ouput from MC and write single track eff - standard dir/list
2165
2166   TString strdir  = "PWG4_FragmentationFunction_" + strID;
2167   TString strlist = "fracfunc_" + strID;
2168
2169   WriteSingleTrackSecCorr(strInfile,strdir,strlist,strOutfile,updateOutfile,strOutDir);
2170 }
2171
2172 //___________________________________________________________________________________________________________________________________
2173 void AliFragmentationFunctionCorrections::WriteSingleTrackSecCorr(TString strInfile, TString strdir, TString strlist,
2174                                                                   TString strOutfile, Bool_t updateOutfile, TString strOutDir)
2175 {
2176   // read task output from MC and write single track secondaries contamination correction as function of pt, eta, phi
2177   // argument strlist optional - read from directory strdir if not specified
2178   // write corr factor to file stroutfile - by default only 'update' file (don't overwrite)
2179
2180   TH1D* hdNdptTracksMCPrimRec;
2181   TH2D* hdNdetadphiTracksMCPrimRec;
2182
2183   TH1D* hdNdptTracksMCSecRec;
2184   TH2D* hdNdetadphiTracksMCSecRec;
2185
2186   TFile f(strInfile,"READ");
2187
2188   if(!f.IsOpen()){
2189     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2190     return;
2191   }
2192
2193   if(fDebug>0) Printf("%s:%d -- writeSingleTrackEff: open task ouput file %s",(char*)__FILE__,__LINE__,strInfile.Data());
2194
2195   if(strdir && strdir.Length()) gDirectory->cd(strdir);
2196
2197   TList* list = 0;
2198
2199   if(strlist && strlist.Length()){
2200
2201     if(!(list = (TList*) gDirectory->Get(strlist))){
2202       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2203       return;
2204     }
2205   }
2206
2207
2208   if(list){
2209     hdNdptTracksMCPrimRec       = (TH1D*) list->FindObject("fh1TrackQAPtRecEffGen");
2210     hdNdetadphiTracksMCPrimRec  = (TH2D*) list->FindObject("fh2TrackQAEtaPhiRecEffGen");
2211
2212     hdNdptTracksMCSecRec       = (TH1D*) list->FindObject("fh1TrackQAPtSecRec");
2213     hdNdetadphiTracksMCSecRec  = (TH2D*) list->FindObject("fh2TrackQAEtaPhiSecRec");
2214   }
2215   else{
2216     hdNdptTracksMCPrimRec       = (TH1D*) gDirectory->Get("fh1TrackQAPtRecEffGen");
2217     hdNdetadphiTracksMCPrimRec  = (TH2D*) gDirectory->Get("fh2TrackQAEtaPhiRecEffGen");
2218
2219     hdNdptTracksMCSecRec       = (TH1D*) gDirectory->Get("fh1TrackQAPtSecRec");
2220     hdNdetadphiTracksMCSecRec  = (TH2D*) gDirectory->Get("fh2TrackQAEtaPhiSecRec");
2221   }
2222
2223   hdNdptTracksMCPrimRec->SetDirectory(0);
2224   hdNdetadphiTracksMCPrimRec->SetDirectory(0);
2225   hdNdptTracksMCSecRec->SetDirectory(0);
2226   hdNdetadphiTracksMCSecRec->SetDirectory(0);
2227
2228   f.Close();
2229
2230   // projections: dN/deta, dN/dphi
2231
2232   TH1D* hdNdetaTracksMCPrimRec = (TH1D*) hdNdetadphiTracksMCPrimRec->ProjectionX("hdNdetaTracksMcPrimRec");
2233   TH1D* hdNdphiTracksMCPrimRec = (TH1D*) hdNdetadphiTracksMCPrimRec->ProjectionY("hdNdphiTracksMcPrimRec");
2234
2235   TH1D* hdNdetaTracksMCSecRec = (TH1D*) hdNdetadphiTracksMCSecRec->ProjectionX("hdNdetaTracksMcSecRec");
2236   TH1D* hdNdphiTracksMCSecRec = (TH1D*) hdNdetadphiTracksMCSecRec->ProjectionY("hdNdphiTracksMcSecRec");
2237
2238
2239   // rebin
2240
2241   TString strNamePtPrim = "hTrackPtPrim";
2242   TString strNamePtSec  = "hTrackPtSec";
2243
2244   if(fNHistoBinsSinglePt) hdNdptTracksMCPrimRec = (TH1D*) hdNdptTracksMCPrimRec->Rebin(fNHistoBinsSinglePt,strNamePtPrim,fHistoBinsSinglePt->GetArray());
2245   if(fNHistoBinsSinglePt) hdNdptTracksMCSecRec  = (TH1D*) hdNdptTracksMCSecRec->Rebin(fNHistoBinsSinglePt,strNamePtSec,fHistoBinsSinglePt->GetArray());
2246
2247
2248   // secondary correction factor: divide prim/(prim+sec)
2249
2250   TH1F* hSingleTrackSecCorrPt = (TH1F*) hdNdptTracksMCSecRec->Clone("hSingleTrackSecCorrPt");
2251   TH1F* hSumPrimSecPt = (TH1F*) hdNdptTracksMCSecRec->Clone("hSumPrimSecPt");
2252   hSumPrimSecPt->Add(hdNdptTracksMCPrimRec);
2253   hSingleTrackSecCorrPt->Divide(hdNdptTracksMCPrimRec,hSumPrimSecPt,1,1,"B"); // binominal errors
2254
2255   TH1F* hSingleTrackSecCorrEta = (TH1F*) hdNdetaTracksMCSecRec->Clone("hSingleTrackSecCorrEta");
2256   TH1F* hSumPrimSecEta = (TH1F*) hdNdetaTracksMCSecRec->Clone("hSumPrimSecEta");
2257   hSumPrimSecEta->Add(hdNdetaTracksMCPrimRec);
2258   hSingleTrackSecCorrEta->Divide(hdNdetaTracksMCPrimRec,hSumPrimSecEta,1,1,"B"); // binominal errors
2259
2260   TH1F* hSingleTrackSecCorrPhi = (TH1F*) hdNdphiTracksMCSecRec->Clone("hSingleTrackSecCorrPhi");
2261   TH1F* hSumPrimSecPhi = (TH1F*) hdNdphiTracksMCSecRec->Clone("hSumPrimSecPhi");
2262   hSumPrimSecPhi->Add(hdNdphiTracksMCPrimRec);
2263   hSingleTrackSecCorrPhi->Divide(hdNdphiTracksMCPrimRec,hSumPrimSecPhi,1,1,"B"); // binominal errors
2264
2265
2266   TString outfileOption = "RECREATE";
2267   if(updateOutfile)  outfileOption = "UPDATE";
2268
2269   TFile out(strOutfile,outfileOption);
2270
2271   if(!out.IsOpen()){
2272     Printf("%s:%d -- error opening secCorr output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
2273     return;
2274   }
2275
2276   if(fDebug>0) Printf("%s:%d -- write secCorr to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
2277
2278   if(strOutDir && strOutDir.Length()){
2279
2280     TDirectory* dir;
2281     if((dir = ((TDirectory*) gDirectory->Get(strOutDir)))) dir->cd();
2282     else{
2283       dir = out.mkdir(strOutDir);
2284       dir->cd();
2285     }
2286   }
2287
2288   hdNdptTracksMCSecRec->Write();
2289   hdNdetaTracksMCSecRec->Write();
2290   hdNdphiTracksMCSecRec->Write();
2291
2292   hSingleTrackSecCorrPt->Write();
2293   hSingleTrackSecCorrEta->Write();
2294   hSingleTrackSecCorrPhi->Write();
2295
2296   out.Close();
2297 }
2298
2299 //________________________________________________________________________________________________________________
2300 void AliFragmentationFunctionCorrections::WriteSingleResponse(TString strInfile, TString strID, TString strOutfile,
2301                                                               Bool_t updateOutfile, TString strOutDir)
2302 {
2303   // read task ouput from MC and write single track eff - standard dir/list
2304
2305   TString strdir  = "PWG4_FragmentationFunction_" + strID;
2306   TString strlist = "fracfunc_" + strID;
2307
2308   WriteSingleResponse(strInfile,strdir,strlist,strOutfile,updateOutfile,strOutDir);
2309 }
2310
2311
2312 //_____________________________________________________________________________________________________________________________________
2313 void AliFragmentationFunctionCorrections::WriteSingleResponse(TString strInfile, TString strdir, TString strlist,
2314                                                               TString strOutfile, Bool_t updateOutfile, TString strOutDir)
2315 {
2316   // read 2d THnSparse response matrices in pt from file
2317   // project TH2
2318   // write to strOutfile
2319
2320   THnSparse* hnResponseSinglePt;
2321   TH2F*      h2ResponseSinglePt;
2322
2323   TFile f(strInfile,"READ");
2324
2325   if(!f.IsOpen()){
2326     Printf("%s:%d -- error opening file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2327     return;
2328   }
2329
2330   if(fDebug>0) Printf("%s:%d -- read response matrices from file %s ",(char*)__FILE__,__LINE__,strInfile.Data());
2331
2332   if(strdir && strdir.Length()) gDirectory->cd(strdir);
2333
2334   TList* list = 0;
2335
2336   if(strlist && strlist.Length()){
2337
2338     if(!(list = (TList*) gDirectory->Get(strlist))){
2339       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2340       return;
2341     }
2342   }
2343
2344   if(list) hnResponseSinglePt = (THnSparse*) list->FindObject("fhnResponseSinglePt");
2345   else     hnResponseSinglePt = (THnSparse*) gDirectory->Get("fhnResponseSinglePt");
2346
2347
2348   if(!hnResponseSinglePt){
2349     Printf("%s:%d -- error retrieving response matrix fhnResponseSinglePt",(char*)__FILE__,__LINE__);
2350     return;
2351   }
2352
2353   f.Close();
2354
2355
2356   // project 2d histo
2357   h2ResponseSinglePt = (TH2F*) hnResponseSinglePt->Projection(1,0);// note convention: yDim,xDim
2358   h2ResponseSinglePt->SetNameTitle("h2ResponseSinglePt","");
2359
2360
2361   // write
2362
2363   TString outfileOption = "RECREATE";
2364   if(updateOutfile)  outfileOption = "UPDATE";
2365
2366   TFile out(strOutfile,outfileOption);
2367
2368   if(!out.IsOpen()){
2369     Printf("%s:%d -- error opening response matrix output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
2370     return;
2371   }
2372
2373   if(fDebug>0) Printf("%s:%d -- write response matrices to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
2374
2375   if(strOutDir && strOutDir.Length()){
2376
2377     TDirectory* dir;
2378     if((dir = ((TDirectory*)  gDirectory->Get(strOutDir)))) dir->cd();
2379     else{
2380       dir = out.mkdir(strOutDir);
2381       dir->cd();
2382     }
2383   }
2384
2385   hnResponseSinglePt->Write();
2386   h2ResponseSinglePt->Write();
2387
2388   out.Close();
2389 }
2390
2391 //________________________________________________________________________________________________________________
2392 void AliFragmentationFunctionCorrections::WriteJetTrackEff(TString strInfile, TString strID, TString strOutfile,
2393                                                            Bool_t updateOutfile)
2394 {
2395   // read task ouput from MC and write single track eff - standard dir/list
2396
2397   TString strdir  = "PWG4_FragmentationFunction_" + strID;
2398   TString strlist = "fracfunc_" + strID;
2399
2400   WriteJetTrackEff(strInfile,strdir,strlist,strOutfile,updateOutfile);
2401 }
2402
2403 //___________________________________________________________________________________________________________________________________
2404 void AliFragmentationFunctionCorrections::WriteJetTrackEff(TString strInfile, TString strdir, TString strlist,
2405                                                            TString strOutfile, Bool_t updateOutfile)
2406 {
2407   // read task output from MC and write track eff in jet pt slices as function of pt, z, xi
2408   // argument strlist optional - read from directory strdir if not specified
2409   // write eff to file strOutfile - by default only 'update' file (don't overwrite)
2410
2411   TH1F* hEffPt[fNJetPtSlices];
2412   TH1F* hEffXi[fNJetPtSlices];
2413   TH1F* hEffZ[fNJetPtSlices];
2414
2415   TH1F* hdNdptTracksMCPrimGen[fNJetPtSlices];
2416   TH1F* hdNdxiMCPrimGen[fNJetPtSlices];
2417   TH1F* hdNdzMCPrimGen[fNJetPtSlices];
2418
2419   TH1F* hdNdptTracksMCPrimRec[fNJetPtSlices];
2420   TH1F* hdNdxiMCPrimRec[fNJetPtSlices];
2421   TH1F* hdNdzMCPrimRec[fNJetPtSlices];
2422
2423
2424   TH1F* fh1FFJetPtRecEffGen;
2425
2426   TH2F* fh2FFTrackPtRecEffGen;
2427   TH2F* fh2FFZRecEffGen;
2428   TH2F* fh2FFXiRecEffGen;
2429
2430   TH2F* fh2FFTrackPtRecEffRec;
2431   TH2F* fh2FFZRecEffRec;
2432   TH2F* fh2FFXiRecEffRec;
2433
2434
2435   TFile f(strInfile,"READ");
2436
2437   if(!f.IsOpen()){
2438     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2439     return;
2440   }
2441
2442   if(fDebug>0) Printf("%s:%d -- writeJetTrackEff: open task ouput file %s",(char*)__FILE__,__LINE__,strInfile.Data());
2443
2444   if(strdir && strdir.Length()) gDirectory->cd(strdir);
2445
2446   TList* list = 0;
2447
2448   if(strlist && strlist.Length()){
2449
2450     if(!(list = (TList*) gDirectory->Get(strlist))){
2451       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2452       return;
2453     }
2454   }
2455
2456   if(list){
2457     fh1FFJetPtRecEffGen    = (TH1F*) list->FindObject("fh1FFJetPtRecEffGen");
2458
2459     fh2FFTrackPtRecEffGen  = (TH2F*) list->FindObject("fh2FFTrackPtRecEffGen");
2460     fh2FFZRecEffGen        = (TH2F*) list->FindObject("fh2FFZRecEffGen");
2461     fh2FFXiRecEffGen       = (TH2F*) list->FindObject("fh2FFXiRecEffGen");
2462
2463     fh2FFTrackPtRecEffRec  = (TH2F*) list->FindObject("fh2FFTrackPtRecEffRec");
2464     fh2FFZRecEffRec        = (TH2F*) list->FindObject("fh2FFZRecEffRec");
2465     fh2FFXiRecEffRec       = (TH2F*) list->FindObject("fh2FFXiRecEffRec");
2466   }
2467   else{
2468     fh1FFJetPtRecEffGen    = (TH1F*) gDirectory->FindObject("fh1FFJetPtRecEffGen");
2469
2470     fh2FFTrackPtRecEffGen  = (TH2F*) gDirectory->FindObject("fh2FFTrackPtRecEffGen");
2471     fh2FFZRecEffGen        = (TH2F*) gDirectory->FindObject("fh2FFZRecEffGen");
2472     fh2FFXiRecEffGen       = (TH2F*) gDirectory->FindObject("fh2FFXiRecEffGen");
2473
2474     fh2FFTrackPtRecEffRec  = (TH2F*) gDirectory->FindObject("fh2FFTrackPtRecEffRec");
2475     fh2FFZRecEffRec        = (TH2F*) gDirectory->FindObject("fh2FFZRecEffRec");
2476     fh2FFXiRecEffRec       = (TH2F*) gDirectory->FindObject("fh2FFXiRecEffRec");
2477   }
2478
2479   fh1FFJetPtRecEffGen->SetDirectory(0);
2480
2481   fh2FFTrackPtRecEffGen->SetDirectory(0);
2482   fh2FFZRecEffGen->SetDirectory(0);
2483   fh2FFXiRecEffGen->SetDirectory(0);
2484
2485   fh2FFTrackPtRecEffRec->SetDirectory(0);
2486   fh2FFZRecEffRec->SetDirectory(0);
2487   fh2FFXiRecEffRec->SetDirectory(0);
2488
2489   f.Close();
2490
2491
2492   // projections: FF for generated and reconstructed primaries
2493
2494   for(Int_t i=0; i<fNJetPtSlices; i++){
2495
2496     Float_t jetPtLoLim = fJetPtSlices->At(i);
2497     Float_t jetPtUpLim = fJetPtSlices->At(i+1);
2498
2499     Int_t binLo = static_cast<Int_t>(fh2FFTrackPtRecEffGen->GetXaxis()->FindBin(jetPtLoLim));
2500     Int_t binUp = static_cast<Int_t>(fh2FFTrackPtRecEffGen->GetXaxis()->FindBin(jetPtUpLim))-1;
2501
2502     if(binUp > fh2FFTrackPtRecEffGen->GetNbinsX()){
2503       Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim);
2504       return;
2505     }
2506
2507     TString strNameFFPtGen(Form("fh1FFTrackPtGenPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2508     TString strNameFFZGen(Form("fh1FFZGenPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2509     TString strNameFFXiGen(Form("fh1FFXiGenPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2510
2511     TString strNameFFPtRec(Form("fh1FFTrackPtRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2512     TString strNameFFZRec(Form("fh1FFZRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2513     TString strNameFFXiRec(Form("fh1FFXiRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2514
2515     // project
2516     // appendix 'unbinned' to avoid histos with same name after rebinning
2517
2518     hdNdptTracksMCPrimGen[i] = (TH1F*) fh2FFTrackPtRecEffGen->ProjectionY(strNameFFPtGen+"_unBinned",binLo,binUp,"o"); // option "o": original axis range
2519     hdNdzMCPrimGen[i]        = (TH1F*) fh2FFZRecEffGen->ProjectionY(strNameFFZGen+"_unBinned",binLo,binUp,"o");
2520     hdNdxiMCPrimGen[i]       = (TH1F*) fh2FFXiRecEffGen->ProjectionY(strNameFFXiGen+"_unBinned",binLo,binUp,"o");
2521
2522     hdNdptTracksMCPrimRec[i] = (TH1F*) fh2FFTrackPtRecEffRec->ProjectionY(strNameFFPtRec+"_unBinned",binLo,binUp,"o"); // option "o": original axis range
2523     hdNdzMCPrimRec[i]        = (TH1F*) fh2FFZRecEffRec->ProjectionY(strNameFFZRec+"_unBinned",binLo,binUp,"o");
2524     hdNdxiMCPrimRec[i]       = (TH1F*) fh2FFXiRecEffRec->ProjectionY(strNameFFXiRec+"_unBinned",binLo,binUp,"o");
2525
2526     // rebin
2527
2528     if(fNHistoBinsPt[i]) hdNdptTracksMCPrimGen[i] = (TH1F*) hdNdptTracksMCPrimGen[i]->Rebin(fNHistoBinsPt[i],strNameFFPtGen,fHistoBinsPt[i]->GetArray());
2529     if(fNHistoBinsZ[i])  hdNdzMCPrimGen[i]  = (TH1F*) hdNdzMCPrimGen[i]->Rebin(fNHistoBinsZ[i],strNameFFZGen,fHistoBinsZ[i]->GetArray());
2530     if(fNHistoBinsXi[i]) hdNdxiMCPrimGen[i] = (TH1F*) hdNdxiMCPrimGen[i]->Rebin(fNHistoBinsXi[i],strNameFFXiGen,fHistoBinsXi[i]->GetArray());
2531
2532     if(fNHistoBinsPt[i]) hdNdptTracksMCPrimRec[i] = (TH1F*) hdNdptTracksMCPrimRec[i]->Rebin(fNHistoBinsPt[i],strNameFFPtRec,fHistoBinsPt[i]->GetArray());
2533     if(fNHistoBinsZ[i])  hdNdzMCPrimRec[i]  = (TH1F*) hdNdzMCPrimRec[i]->Rebin(fNHistoBinsZ[i],strNameFFZRec,fHistoBinsZ[i]->GetArray());
2534     if(fNHistoBinsXi[i]) hdNdxiMCPrimRec[i] = (TH1F*) hdNdxiMCPrimRec[i]->Rebin(fNHistoBinsXi[i],strNameFFXiRec,fHistoBinsXi[i]->GetArray());
2535
2536     hdNdptTracksMCPrimGen[i]->SetNameTitle(strNameFFPtGen,"");
2537     hdNdzMCPrimGen[i]->SetNameTitle(strNameFFZGen,"");
2538     hdNdxiMCPrimGen[i]->SetNameTitle(strNameFFXiGen,"");
2539
2540     hdNdptTracksMCPrimRec[i]->SetNameTitle(strNameFFPtRec,"");
2541     hdNdzMCPrimRec[i]->SetNameTitle(strNameFFZRec,"");
2542     hdNdxiMCPrimRec[i]->SetNameTitle(strNameFFXiRec,"");
2543
2544     // normalize
2545
2546     Double_t nJetsBin = fh1FFJetPtRecEffGen->Integral(binLo,binUp);
2547
2548     NormalizeTH1(hdNdptTracksMCPrimGen[i],nJetsBin);
2549     NormalizeTH1(hdNdzMCPrimGen[i],nJetsBin);
2550     NormalizeTH1(hdNdxiMCPrimGen[i],nJetsBin);
2551
2552     NormalizeTH1(hdNdptTracksMCPrimRec[i],nJetsBin);
2553     NormalizeTH1(hdNdzMCPrimRec[i],nJetsBin);
2554     NormalizeTH1(hdNdxiMCPrimRec[i],nJetsBin);
2555
2556     // divide rec/gen : efficiency
2557
2558     TString strNameEffPt(Form("hEffPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2559     TString strNameEffZ(Form("hEffZ_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2560     TString strNameEffXi(Form("hEffXi_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2561
2562     hEffPt[i] = (TH1F*) hdNdptTracksMCPrimRec[i]->Clone(strNameEffPt);
2563     hEffPt[i]->Divide(hdNdptTracksMCPrimRec[i],hdNdptTracksMCPrimGen[i],1,1,"B"); // binominal errors
2564
2565     hEffXi[i] = (TH1F*) hdNdxiMCPrimRec[i]->Clone(strNameEffXi);
2566     hEffXi[i]->Divide(hdNdxiMCPrimRec[i],hdNdxiMCPrimGen[i],1,1,"B"); // binominal errors
2567
2568     hEffZ[i] = (TH1F*) hdNdzMCPrimRec[i]->Clone(strNameEffZ);
2569     hEffZ[i]->Divide(hdNdzMCPrimRec[i],hdNdzMCPrimGen[i],1,1,"B"); // binominal errors
2570   }
2571
2572   // write
2573
2574   TString outfileOption = "RECREATE";
2575   if(updateOutfile)  outfileOption = "UPDATE";
2576
2577   TFile out(strOutfile,outfileOption);
2578
2579   if(!out.IsOpen()){
2580     Printf("%s:%d -- error opening efficiency output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
2581     return;
2582   }
2583
2584   if(fDebug>0) Printf("%s:%d -- write efficiency to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
2585
2586 //   if(strdir && strdir.Length()){
2587 //     TDirectory* dir = out.mkdir(strdir);
2588 //     dir->cd();
2589 //   }
2590
2591   for(Int_t i=0; i<fNJetPtSlices; i++){
2592
2593     hdNdptTracksMCPrimGen[i]->Write();
2594     hdNdxiMCPrimGen[i]->Write();
2595     hdNdzMCPrimGen[i]->Write();
2596
2597     hdNdptTracksMCPrimRec[i]->Write();
2598     hdNdxiMCPrimRec[i]->Write();
2599     hdNdzMCPrimRec[i]->Write();
2600
2601     hEffPt[i]->Write();
2602     hEffXi[i]->Write();
2603     hEffZ[i]->Write();
2604   }
2605
2606   out.Close();
2607
2608 }
2609
2610 //________________________________________________________________________________________________________________
2611 void AliFragmentationFunctionCorrections::WriteJetSecCorr(TString strInfile, TString strID, TString strOutfile,
2612                                                            Bool_t updateOutfile)
2613 {
2614   // read task ouput from MC and write secondary correction - standard dir/list
2615
2616   TString strdir  = "PWG4_FragmentationFunction_" + strID;
2617   TString strlist = "fracfunc_" + strID;
2618
2619   WriteJetSecCorr(strInfile,strdir,strlist,strOutfile,updateOutfile);
2620 }
2621
2622 //___________________________________________________________________________________________________________________________________
2623 void AliFragmentationFunctionCorrections::WriteJetSecCorr(TString strInfile, TString strdir, TString strlist,
2624                                                            TString strOutfile, Bool_t updateOutfile)
2625 {
2626   // read task output from MC and write secondary correction in jet pt slices as function of pt, z, xi
2627   // argument strlist optional - read from directory strdir if not specified
2628   // write eff to file strOutfile - by default only 'update' file (don't overwrite)
2629
2630   TH1F* hSecCorrPt[fNJetPtSlices];
2631   TH1F* hSecCorrXi[fNJetPtSlices];
2632   TH1F* hSecCorrZ[fNJetPtSlices];
2633
2634   TH1F* hdNdptTracksMCPrimRec[fNJetPtSlices];
2635   TH1F* hdNdxiMCPrimRec[fNJetPtSlices];
2636   TH1F* hdNdzMCPrimRec[fNJetPtSlices];
2637
2638   TH1F* hdNdptTracksMCSecRec[fNJetPtSlices];
2639   TH1F* hdNdxiMCSecRec[fNJetPtSlices];
2640   TH1F* hdNdzMCSecRec[fNJetPtSlices];
2641
2642   TH1F* fh1FFJetPtRecEffGen;
2643
2644   TH2F* fh2FFTrackPtRecEffRec;
2645   TH2F* fh2FFZRecEffRec;
2646   TH2F* fh2FFXiRecEffRec;
2647
2648   TH2F* fh2FFTrackPtSecRec;
2649   TH2F* fh2FFZSecRec;
2650   TH2F* fh2FFXiSecRec;
2651
2652   TFile f(strInfile,"READ");
2653
2654   if(!f.IsOpen()){
2655     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2656     return;
2657   }
2658
2659   if(fDebug>0) Printf("%s:%d -- writeJetTrackEff: open task ouput file %s",(char*)__FILE__,__LINE__,strInfile.Data());
2660
2661   if(strdir && strdir.Length()) gDirectory->cd(strdir);
2662
2663   TList* list = 0;
2664
2665   if(strlist && strlist.Length()){
2666
2667     if(!(list = (TList*) gDirectory->Get(strlist))){
2668       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2669       return;
2670     }
2671   }
2672
2673   if(list){
2674     fh1FFJetPtRecEffGen    = (TH1F*) list->FindObject("fh1FFJetPtRecEffGen");
2675
2676     fh2FFTrackPtRecEffRec  = (TH2F*) list->FindObject("fh2FFTrackPtRecEffRec");
2677     fh2FFZRecEffRec        = (TH2F*) list->FindObject("fh2FFZRecEffRec");
2678     fh2FFXiRecEffRec       = (TH2F*) list->FindObject("fh2FFXiRecEffRec");
2679
2680     fh2FFTrackPtSecRec  = (TH2F*) list->FindObject("fh2FFTrackPtSecRec");
2681     fh2FFZSecRec        = (TH2F*) list->FindObject("fh2FFZSecRec");
2682     fh2FFXiSecRec       = (TH2F*) list->FindObject("fh2FFXiSecRec");
2683   }
2684   else{
2685     fh1FFJetPtRecEffGen    = (TH1F*) gDirectory->FindObject("fh1FFJetPtRecEffGen");
2686
2687     fh2FFTrackPtRecEffRec  = (TH2F*) gDirectory->FindObject("fh2FFTrackPtRecEffRec");
2688     fh2FFZRecEffRec        = (TH2F*) gDirectory->FindObject("fh2FFZRecEffRec");
2689     fh2FFXiRecEffRec       = (TH2F*) gDirectory->FindObject("fh2FFXiRecEffRec");
2690
2691     fh2FFTrackPtSecRec  = (TH2F*) gDirectory->FindObject("fh2FFTrackPtSecRec");
2692     fh2FFZSecRec        = (TH2F*) gDirectory->FindObject("fh2FFZSecRec");
2693     fh2FFXiSecRec       = (TH2F*) gDirectory->FindObject("fh2FFXiSecRec");
2694   }
2695
2696   fh1FFJetPtRecEffGen->SetDirectory(0);
2697
2698   fh2FFTrackPtRecEffRec->SetDirectory(0);
2699   fh2FFZRecEffRec->SetDirectory(0);
2700   fh2FFXiRecEffRec->SetDirectory(0);
2701
2702   fh2FFTrackPtSecRec->SetDirectory(0);
2703   fh2FFZSecRec->SetDirectory(0);
2704   fh2FFXiSecRec->SetDirectory(0);
2705
2706   f.Close();
2707
2708
2709   // projections: FF for generated and reconstructed primaries
2710
2711   for(Int_t i=0; i<fNJetPtSlices; i++){
2712
2713     Float_t jetPtLoLim = fJetPtSlices->At(i);
2714     Float_t jetPtUpLim = fJetPtSlices->At(i+1);
2715
2716     Int_t binLo = static_cast<Int_t>(fh2FFTrackPtRecEffRec->GetXaxis()->FindBin(jetPtLoLim));
2717     Int_t binUp = static_cast<Int_t>(fh2FFTrackPtRecEffRec->GetXaxis()->FindBin(jetPtUpLim))-1;
2718
2719     if(binUp > fh2FFTrackPtRecEffRec->GetNbinsX()){
2720       Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim);
2721       return;
2722     }
2723
2724     TString strNameFFPtPrimRec(Form("fh1FFTrackPtRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2725     TString strNameFFZPrimRec(Form("fh1FFZRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2726     TString strNameFFXiPrimRec(Form("fh1FFXiRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2727
2728     TString strNameFFPtSecRec(Form("fh1FFTrackPtRecSec_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2729     TString strNameFFZSecRec(Form("fh1FFZRecSec_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2730     TString strNameFFXiSecRec(Form("fh1FFXiRecSec_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2731
2732     // project
2733     // appendix 'unbinned' to avoid histos with same name after rebinning
2734
2735     hdNdptTracksMCPrimRec[i] = (TH1F*) fh2FFTrackPtRecEffRec->ProjectionY(strNameFFPtPrimRec+"_unBinned",binLo,binUp,"o"); // option "o": original axis range
2736     hdNdzMCPrimRec[i]        = (TH1F*) fh2FFZRecEffRec->ProjectionY(strNameFFZPrimRec+"_unBinned",binLo,binUp,"o");
2737     hdNdxiMCPrimRec[i]       = (TH1F*) fh2FFXiRecEffRec->ProjectionY(strNameFFXiPrimRec+"_unBinned",binLo,binUp,"o");
2738
2739     hdNdptTracksMCSecRec[i]  = (TH1F*) fh2FFTrackPtSecRec->ProjectionY(strNameFFPtSecRec+"_unBinned",binLo,binUp,"o"); // option "o": original axis range
2740     hdNdzMCSecRec[i]         = (TH1F*) fh2FFZSecRec->ProjectionY(strNameFFZSecRec+"_unBinned",binLo,binUp,"o");
2741     hdNdxiMCSecRec[i]        = (TH1F*) fh2FFXiSecRec->ProjectionY(strNameFFXiSecRec+"_unBinned",binLo,binUp,"o");
2742
2743     // rebin
2744
2745     if(fNHistoBinsPt[i]) hdNdptTracksMCPrimRec[i] = (TH1F*) hdNdptTracksMCPrimRec[i]->Rebin(fNHistoBinsPt[i],strNameFFPtPrimRec,fHistoBinsPt[i]->GetArray());
2746     if(fNHistoBinsZ[i])  hdNdzMCPrimRec[i]        = (TH1F*) hdNdzMCPrimRec[i]->Rebin(fNHistoBinsZ[i],strNameFFZPrimRec,fHistoBinsZ[i]->GetArray());
2747     if(fNHistoBinsXi[i]) hdNdxiMCPrimRec[i]       = (TH1F*) hdNdxiMCPrimRec[i]->Rebin(fNHistoBinsXi[i],strNameFFXiPrimRec,fHistoBinsXi[i]->GetArray());
2748
2749     if(fNHistoBinsPt[i]) hdNdptTracksMCSecRec[i] = (TH1F*) hdNdptTracksMCSecRec[i]->Rebin(fNHistoBinsPt[i],strNameFFPtSecRec,fHistoBinsPt[i]->GetArray());
2750     if(fNHistoBinsZ[i])  hdNdzMCSecRec[i]        = (TH1F*) hdNdzMCSecRec[i]->Rebin(fNHistoBinsZ[i],strNameFFZSecRec,fHistoBinsZ[i]->GetArray());
2751     if(fNHistoBinsXi[i]) hdNdxiMCSecRec[i]       = (TH1F*) hdNdxiMCSecRec[i]->Rebin(fNHistoBinsXi[i],strNameFFXiSecRec,fHistoBinsXi[i]->GetArray());
2752
2753     hdNdptTracksMCPrimRec[i]->SetNameTitle(strNameFFPtPrimRec,"");
2754     hdNdzMCPrimRec[i]->SetNameTitle(strNameFFZPrimRec,"");
2755     hdNdxiMCPrimRec[i]->SetNameTitle(strNameFFXiPrimRec,"");
2756
2757     hdNdptTracksMCSecRec[i]->SetNameTitle(strNameFFPtSecRec,"");
2758     hdNdzMCSecRec[i]->SetNameTitle(strNameFFZSecRec,"");
2759     hdNdxiMCSecRec[i]->SetNameTitle(strNameFFXiSecRec,"");
2760
2761     // normalize
2762
2763     Double_t nJetsBin = fh1FFJetPtRecEffGen->Integral(binLo,binUp);
2764
2765     NormalizeTH1(hdNdptTracksMCPrimRec[i],nJetsBin);
2766     NormalizeTH1(hdNdzMCPrimRec[i],nJetsBin);
2767     NormalizeTH1(hdNdxiMCPrimRec[i],nJetsBin);
2768
2769     NormalizeTH1(hdNdptTracksMCSecRec[i],nJetsBin);
2770     NormalizeTH1(hdNdzMCSecRec[i],nJetsBin);
2771     NormalizeTH1(hdNdxiMCSecRec[i],nJetsBin);
2772
2773     // divide rec/gen : efficiency
2774
2775     TString strNameSecCorrPt(Form("hSecCorrPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2776     TString strNameSecCorrZ(Form("hSecCorrZ_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2777     TString strNameSecCorrXi(Form("hSecCorrXi_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2778
2779     hSecCorrPt[i] = (TH1F*) hdNdptTracksMCSecRec[i]->Clone(strNameSecCorrPt);
2780     TH1F* hSumPrimSecPt = (TH1F*) hdNdptTracksMCSecRec[i]->Clone("hSumPrimSecPt");
2781     hSumPrimSecPt->Add(hdNdptTracksMCPrimRec[i]);
2782     hSecCorrPt[i]->Divide(hdNdptTracksMCPrimRec[i],hSumPrimSecPt,1,1,"B"); // binominal errors
2783
2784     hSecCorrXi[i] = (TH1F*) hdNdxiMCSecRec[i]->Clone(strNameSecCorrXi);
2785     TH1F* hSumPrimSecXi = (TH1F*) hdNdxiMCSecRec[i]->Clone("hSumPrimSecXi");
2786     hSumPrimSecXi->Add(hdNdxiMCPrimRec[i]);
2787     hSecCorrXi[i]->Divide(hdNdxiMCPrimRec[i],hSumPrimSecXi,1,1,"B"); // binominal errors
2788
2789     hSecCorrZ[i] = (TH1F*) hdNdzMCSecRec[i]->Clone(strNameSecCorrZ);
2790     TH1F* hSumPrimSecZ = (TH1F*) hdNdzMCSecRec[i]->Clone("hSumPrimSecZ");
2791     hSumPrimSecZ->Add(hdNdzMCPrimRec[i]);
2792     hSecCorrZ[i]->Divide(hdNdzMCPrimRec[i],hSumPrimSecZ,1,1,"B"); // binominal errors
2793   }
2794
2795   // write
2796
2797   TString outfileOption = "RECREATE";
2798   if(updateOutfile)  outfileOption = "UPDATE";
2799
2800   TFile out(strOutfile,outfileOption);
2801
2802   if(!out.IsOpen()){
2803     Printf("%s:%d -- error opening efficiency output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
2804     return;
2805   }
2806
2807   if(fDebug>0) Printf("%s:%d -- write efficiency to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
2808
2809   for(Int_t i=0; i<fNJetPtSlices; i++){
2810
2811     //  hdNdptTracksMCSecRec[i]->Write();
2812     //  hdNdxiMCSecRec[i]->Write();
2813     //  hdNdzMCSecRec[i]->Write();
2814
2815     hSecCorrPt[i]->Write();
2816     hSecCorrXi[i]->Write();
2817     hSecCorrZ[i]->Write();
2818   }
2819
2820   out.Close();
2821 }
2822
2823 //________________________________________________________________________________________________________________
2824 void AliFragmentationFunctionCorrections::WriteJetResponse(TString strInfile, TString strID, TString strOutfile,
2825                                                            Bool_t updateOutfile)
2826 {
2827   // read task ouput from MC and write single track eff - standard dir/list
2828
2829   TString strdir  = "PWG4_FragmentationFunction_" + strID;
2830   TString strlist = "fracfunc_" + strID;
2831
2832   WriteJetResponse(strInfile,strdir,strlist,strOutfile,updateOutfile);
2833 }
2834
2835 //_____________________________________________________________________________________________________________________________________
2836 void AliFragmentationFunctionCorrections::WriteJetResponse(TString strInfile, TString strdir, TString strlist,
2837                                                            TString strOutfile, Bool_t updateOutfile)
2838 {
2839   // read 3d THnSparse response matrices in pt,z,xi vs jet pt from file
2840   // project THnSparse + TH2 in jet pt slices
2841   // write to strOutfile
2842
2843   THnSparse* hn3ResponseJetPt;
2844   THnSparse* hn3ResponseJetZ;
2845   THnSparse* hn3ResponseJetXi;
2846
2847   // 2D response matrices
2848
2849   THnSparse* hnResponsePt[fNJetPtSlices];
2850   THnSparse* hnResponseZ[fNJetPtSlices];
2851   THnSparse* hnResponseXi[fNJetPtSlices];
2852
2853   TH2F* h2ResponsePt[fNJetPtSlices];
2854   TH2F* h2ResponseZ[fNJetPtSlices];
2855   TH2F* h2ResponseXi[fNJetPtSlices];
2856
2857   // 1D projections on gen pt / rec pt axes
2858
2859   TH1F* h1FFPtRec[fNJetPtSlices];
2860   TH1F* h1FFZRec[fNJetPtSlices];
2861   TH1F* h1FFXiRec[fNJetPtSlices];
2862
2863   TH1F* h1FFPtGen[fNJetPtSlices];
2864   TH1F* h1FFZGen[fNJetPtSlices];
2865   TH1F* h1FFXiGen[fNJetPtSlices];
2866
2867   TH1F* h1RatioPt[fNJetPtSlices];
2868   TH1F* h1RatioZ[fNJetPtSlices];
2869   TH1F* h1RatioXi[fNJetPtSlices];
2870
2871
2872
2873   TFile f(strInfile,"READ");
2874
2875   if(!f.IsOpen()){
2876     Printf("%s:%d -- error opening file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2877     return;
2878   }
2879
2880   if(fDebug>0) Printf("%s:%d -- read response matrices from file %s ",(char*)__FILE__,__LINE__,strInfile.Data());
2881
2882   if(strdir && strdir.Length()) gDirectory->cd(strdir);
2883
2884   TList* list = 0;
2885
2886   if(strlist && strlist.Length()){
2887
2888     if(!(list = (TList*) gDirectory->Get(strlist))){
2889       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2890       return;
2891     }
2892   }
2893
2894   if(list){
2895     hn3ResponseJetPt = (THnSparse*) list->FindObject("fhnResponseJetTrackPt");
2896     hn3ResponseJetZ  = (THnSparse*) list->FindObject("fhnResponseJetZ");
2897     hn3ResponseJetXi = (THnSparse*) list->FindObject("fhnResponseJetXi");
2898   }
2899   else{
2900     hn3ResponseJetPt = (THnSparse*) gDirectory->Get("fhnResponseJetTrackPt");
2901     hn3ResponseJetZ  = (THnSparse*) gDirectory->Get("fhnResponseJetZ");
2902     hn3ResponseJetXi = (THnSparse*) gDirectory->Get("fhnResponseJetXi");
2903   }
2904
2905
2906   if(!hn3ResponseJetPt){
2907     Printf("%s:%d -- error retrieving response matrix fhnResponseJetTrackPt",(char*)__FILE__,__LINE__);
2908     return;
2909   }
2910
2911   if(!hn3ResponseJetZ){
2912     Printf("%s:%d -- error retrieving response matrix fhnResponseJetZ",(char*)__FILE__,__LINE__);
2913     return;
2914   }
2915
2916   if(!hn3ResponseJetXi){
2917     Printf("%s:%d -- error retrieving response matrix fhnResponseJetXi",(char*)__FILE__,__LINE__);
2918     return;
2919   }
2920
2921   f.Close();
2922
2923   // axes
2924
2925   Int_t axisJetPtTHn3 = -1;
2926   Int_t axisGenPtTHn3 = -1;
2927   Int_t axisRecPtTHn3 = -1;
2928
2929   for(Int_t i=0; i<hn3ResponseJetPt->GetNdimensions(); i++){
2930
2931     TString title = hn3ResponseJetPt->GetAxis(i)->GetTitle();
2932
2933     if(title.Contains("jet p_{T}")) axisJetPtTHn3 = i;
2934     if(title.Contains("gen p_{T}")) axisGenPtTHn3 = i;
2935     if(title.Contains("rec p_{T}")) axisRecPtTHn3 = i;
2936   }
2937
2938   if(axisJetPtTHn3 == -1){
2939     Printf("%s:%d -- error axisJetPtTHn3",(char*)__FILE__,__LINE__);
2940     return;
2941   }
2942
2943   if(axisGenPtTHn3 == -1){
2944     Printf("%s:%d -- error axisGenPtTHn3",(char*)__FILE__,__LINE__);
2945     return;
2946   }
2947
2948   if(axisRecPtTHn3 == -1){
2949     Printf("%s:%d -- error axisRecPtTHn3",(char*)__FILE__,__LINE__);
2950     return;
2951   }
2952
2953   for(Int_t i=0; i<fNJetPtSlices; i++){
2954
2955     Float_t jetPtLoLim = fJetPtSlices->At(i);
2956     Float_t jetPtUpLim = fJetPtSlices->At(i+1);
2957
2958     Int_t binLo = static_cast<Int_t>(hn3ResponseJetPt->GetAxis(axisJetPtTHn3)->FindBin(jetPtLoLim));
2959     Int_t binUp = static_cast<Int_t>(hn3ResponseJetPt->GetAxis(axisJetPtTHn3)->FindBin(jetPtUpLim))-1;
2960
2961     if(binUp > hn3ResponseJetPt->GetAxis(axisJetPtTHn3)->GetNbins()){
2962       Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim);
2963       return;
2964     }
2965
2966     TString strNameRespPt(Form("hnResponsePt_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2967     TString strNameRespZ(Form("hnResponseZ_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2968     TString strNameRespXi(Form("hnResponseXi_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2969
2970     TString strNameTH2RespPt(Form("h2ResponsePt_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2971     TString strNameTH2RespZ(Form("h2ResponseZ_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2972     TString strNameTH2RespXi(Form("h2ResponseXi_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2973
2974     TString strNameRecPt(Form("h1FFTrackPtRecPrim_recPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2975     TString strNameRecZ(Form("h1FFZRecPrim_recPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2976     TString strNameRecXi(Form("h1FFXiRecPrim_recPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2977
2978     TString strNameGenPt(Form("h1FFTrackPtRecPrim_genPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2979     TString strNameGenZ(Form("h1FFZRecPrim_genPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2980     TString strNameGenXi(Form("h1FFXiRecPrim_genPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2981
2982     TString strNameRatioPt(Form("h1RatioTrackPtRecPrim_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2983     TString strNameRatioZ(Form("h1RatioZRecPrim_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2984     TString strNameRatioXi(Form("h1RatioXiRecPrim_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2985
2986
2987     // 2D projections in jet pt range
2988
2989     hn3ResponseJetPt->GetAxis(axisJetPtTHn3)->SetRange(binLo,binUp);
2990     hn3ResponseJetZ->GetAxis(axisJetPtTHn3)->SetRange(binLo,binUp);
2991     hn3ResponseJetXi->GetAxis(axisJetPtTHn3)->SetRange(binLo,binUp);
2992
2993     Int_t axesProj[2] = {axisRecPtTHn3,axisGenPtTHn3};
2994
2995     hnResponsePt[i] = hn3ResponseJetPt->Projection(2,axesProj);
2996     hnResponseZ[i]  = hn3ResponseJetZ->Projection(2,axesProj);
2997     hnResponseXi[i] = hn3ResponseJetXi->Projection(2,axesProj);
2998
2999     hnResponsePt[i]->SetNameTitle(strNameRespPt,"");
3000     hnResponseZ[i]->SetNameTitle(strNameRespZ,"");
3001     hnResponseXi[i]->SetNameTitle(strNameRespXi,"");
3002
3003     h2ResponsePt[i] = (TH2F*) hnResponsePt[i]->Projection(1,0);// note convention: yDim,xDim
3004     h2ResponseZ[i]  = (TH2F*) hnResponseZ[i]->Projection(1,0); // note convention: yDim,xDim
3005     h2ResponseXi[i] = (TH2F*) hnResponseXi[i]->Projection(1,0);// note convention: yDim,xDim
3006
3007     h2ResponsePt[i]->SetNameTitle(strNameTH2RespPt,"");
3008     h2ResponseZ[i]->SetNameTitle(strNameTH2RespZ,"");
3009     h2ResponseXi[i]->SetNameTitle(strNameTH2RespXi,"");
3010
3011
3012     // 1D projections
3013
3014     Int_t axisGenPtTHn2 = -1;
3015     Int_t axisRecPtTHn2 = -1;
3016
3017     for(Int_t d=0; d<hnResponsePt[i]->GetNdimensions(); d++){
3018
3019       TString title = hnResponsePt[i]->GetAxis(d)->GetTitle();
3020
3021       if(title.Contains("gen p_{T}")) axisGenPtTHn2 = d;
3022       if(title.Contains("rec p_{T}")) axisRecPtTHn2 = d;
3023     }
3024
3025
3026     if(axisGenPtTHn2 == -1){
3027       Printf("%s:%d -- error axisGenPtTHn2",(char*)__FILE__,__LINE__);
3028       return;
3029     }
3030
3031     if(axisRecPtTHn2 == -1){
3032       Printf("%s:%d -- error axisRecPtTHn2",(char*)__FILE__,__LINE__);
3033       return;
3034     }
3035
3036
3037     h1FFPtRec[i] = (TH1F*) hnResponsePt[i]->Projection(axisRecPtTHn2);// note convention: yDim,xDim
3038     h1FFZRec[i]  = (TH1F*) hnResponseZ[i]->Projection(axisRecPtTHn2);// note convention: yDim,xDim
3039     h1FFXiRec[i] = (TH1F*) hnResponseXi[i]->Projection(axisRecPtTHn2);// note convention: yDim,xDim
3040
3041     h1FFPtRec[i]->SetNameTitle(strNameRecPt,"");
3042     h1FFZRec[i]->SetNameTitle(strNameRecZ,"");
3043     h1FFXiRec[i]->SetNameTitle(strNameRecXi,"");
3044
3045     h1FFPtGen[i] = (TH1F*) hnResponsePt[i]->Projection(axisGenPtTHn2);// note convention: yDim,xDim
3046     h1FFZGen[i]  = (TH1F*) hnResponseZ[i]->Projection(axisGenPtTHn2);// note convention: yDim,xDim
3047     h1FFXiGen[i] = (TH1F*) hnResponseXi[i]->Projection(axisGenPtTHn2);// note convention: yDim,xDim
3048
3049     h1FFPtGen[i]->SetNameTitle(strNameGenPt,"");
3050     h1FFZGen[i]->SetNameTitle(strNameGenZ,"");
3051     h1FFXiGen[i]->SetNameTitle(strNameGenXi,"");
3052
3053     // normalize 1D projections
3054
3055     if(fNHistoBinsPt[i]) h1FFPtRec[i] = (TH1F*) h1FFPtRec[i]->Rebin(fNHistoBinsPt[i],"",fHistoBinsPt[i]->GetArray());
3056     if(fNHistoBinsZ[i])  h1FFZRec[i]  = (TH1F*) h1FFZRec[i]->Rebin(fNHistoBinsZ[i],"",fHistoBinsZ[i]->GetArray());
3057     if(fNHistoBinsXi[i]) h1FFXiRec[i] = (TH1F*) h1FFXiRec[i]->Rebin(fNHistoBinsXi[i],"",fHistoBinsXi[i]->GetArray());
3058
3059     if(fNHistoBinsPt[i]) h1FFPtGen[i] = (TH1F*) h1FFPtGen[i]->Rebin(fNHistoBinsPt[i],"",fHistoBinsPt[i]->GetArray());
3060     if(fNHistoBinsZ[i])  h1FFZGen[i]  = (TH1F*) h1FFZGen[i]->Rebin(fNHistoBinsZ[i],"",fHistoBinsZ[i]->GetArray());
3061     if(fNHistoBinsXi[i]) h1FFXiGen[i] = (TH1F*) h1FFXiGen[i]->Rebin(fNHistoBinsXi[i],"",fHistoBinsXi[i]->GetArray());
3062
3063     NormalizeTH1(h1FFPtRec[i],fNJets->At(i));
3064     NormalizeTH1(h1FFZRec[i],fNJets->At(i));
3065     NormalizeTH1(h1FFXiRec[i],fNJets->At(i));
3066
3067     NormalizeTH1(h1FFPtGen[i],fNJets->At(i));
3068     NormalizeTH1(h1FFZGen[i],fNJets->At(i));
3069     NormalizeTH1(h1FFXiGen[i],fNJets->At(i));
3070
3071     // ratios 1D projections
3072
3073     h1RatioPt[i] = (TH1F*) h1FFPtRec[i]->Clone(strNameRatioPt);
3074     h1RatioPt[i]->Reset();
3075     h1RatioPt[i]->Divide(h1FFPtRec[i],h1FFPtGen[i],1,1,"B");
3076
3077     h1RatioZ[i] = (TH1F*) h1FFZRec[i]->Clone(strNameRatioZ);
3078     h1RatioZ[i]->Reset();
3079     h1RatioZ[i]->Divide(h1FFZRec[i],h1FFZGen[i],1,1,"B");
3080
3081     h1RatioXi[i] = (TH1F*) h1FFXiRec[i]->Clone(strNameRatioXi);
3082     h1RatioXi[i]->Reset();
3083     h1RatioXi[i]->Divide(h1FFXiRec[i],h1FFXiGen[i],1,1,"B");
3084   }
3085
3086
3087   // write
3088
3089   TString outfileOption = "RECREATE";
3090   if(updateOutfile)  outfileOption = "UPDATE";
3091
3092   TFile out(strOutfile,outfileOption);
3093
3094   if(!out.IsOpen()){
3095     Printf("%s:%d -- error opening response matrix output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
3096     return;
3097   }
3098
3099   if(fDebug>0) Printf("%s:%d -- write response matrices to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
3100
3101   //if(strdir && strdir.Length()){
3102   //  TDirectory* dir = out.mkdir(strdir);
3103   //  dir->cd();
3104   //}
3105
3106   for(Int_t i=0; i<fNJetPtSlices; i++){
3107
3108     hnResponsePt[i]->Write();
3109     hnResponseXi[i]->Write();
3110     hnResponseZ[i]->Write();
3111
3112     h2ResponsePt[i]->Write();
3113     h2ResponseXi[i]->Write();
3114     h2ResponseZ[i]->Write();
3115
3116     h1FFPtRec[i]->Write();
3117     h1FFZRec[i]->Write();
3118     h1FFXiRec[i]->Write();
3119
3120     h1FFPtGen[i]->Write();
3121     h1FFZGen[i]->Write();
3122     h1FFXiGen[i]->Write();
3123
3124     h1RatioPt[i]->Write();
3125     h1RatioZ[i]->Write();
3126     h1RatioXi[i]->Write();
3127
3128   }
3129
3130   out.Close();
3131 }
3132
3133 //______________________________________________________________________________________________________
3134 void AliFragmentationFunctionCorrections::ReadResponse(TString strfile, TString strdir, TString strlist)
3135 {
3136   // read response matrices from file
3137   // argument strlist optional - read from directory strdir if not specified
3138   // note: THnSparse are not rebinned
3139
3140   THnSparse* hRespPt[fNJetPtSlices];
3141   THnSparse* hRespZ[fNJetPtSlices];
3142   THnSparse* hRespXi[fNJetPtSlices];
3143
3144   for(Int_t i=0; i<fNJetPtSlices; i++) hRespPt[i] = 0;
3145   for(Int_t i=0; i<fNJetPtSlices; i++) hRespZ[i]  = 0;
3146   for(Int_t i=0; i<fNJetPtSlices; i++) hRespXi[i] = 0;
3147
3148   TFile f(strfile,"READ");
3149
3150   if(!f.IsOpen()){
3151     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
3152     return;
3153   }
3154
3155   if(fDebug>0) Printf("%s:%d -- read response matrices from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
3156
3157   if(strdir && strdir.Length()) gDirectory->cd(strdir);
3158
3159   TList* list = 0;
3160
3161   if(strlist && strlist.Length()){
3162
3163     if(!(list = (TList*) gDirectory->Get(strlist))){
3164       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
3165       return;
3166     }
3167   }
3168
3169   for(Int_t i=0; i<fNJetPtSlices; i++){
3170
3171     Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
3172     Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
3173
3174     TString strNameRespPt(Form("hnResponsePt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3175     TString strNameRespZ(Form("hnResponseZ_%02d_%02d",jetPtLoLim,jetPtUpLim));
3176     TString strNameRespXi(Form("hnResponseXi_%02d_%02d",jetPtLoLim,jetPtUpLim));
3177
3178     if(list){
3179       hRespPt[i] = (THnSparse*) list->FindObject(strNameRespPt);
3180       hRespZ[i]  = (THnSparse*) list->FindObject(strNameRespZ);
3181       hRespXi[i] = (THnSparse*) list->FindObject(strNameRespXi);
3182     }
3183     else{
3184       hRespPt[i] = (THnSparse*) gDirectory->Get(strNameRespPt);
3185       hRespZ[i]  = (THnSparse*) gDirectory->Get(strNameRespZ);
3186       hRespXi[i] = (THnSparse*) gDirectory->Get(strNameRespXi);
3187     }
3188
3189     if(!hRespPt[i]){
3190       Printf("%s:%d -- error retrieving response %s", (char*)__FILE__,__LINE__,strNameRespPt.Data());
3191     }
3192
3193     if(!hRespZ[i]){
3194       Printf("%s:%d -- error retrieving response %s", (char*)__FILE__,__LINE__,strNameRespZ.Data());
3195     }
3196
3197     if(!hRespXi[i]){
3198       Printf("%s:%d -- error retrieving response %s", (char*)__FILE__,__LINE__,strNameRespXi.Data());
3199     }
3200
3201     //    if(0){ // can't rebin THnSparse ...
3202     //       if(fNHistoBinsPt[i]) hRespPt[i]->SetBinEdges(0,fHistoBinsPt[i]->GetArray());
3203     //       if(fNHistoBinsPt[i]) hRespPt[i]->SetBinEdges(1,fHistoBinsPt[i]->GetArray());
3204
3205     //       if(fNHistoBinsZ[i])  hRespZ[i]->SetBinEdges(0,fHistoBinsZ[i]->GetArray());
3206     //       if(fNHistoBinsZ[i])  hRespZ[i]->SetBinEdges(1,fHistoBinsZ[i]->GetArray());
3207
3208     //       if(fNHistoBinsXi[i]) hRespXi[i]->SetBinEdges(0,fHistoBinsXi[i]->GetArray());
3209     //       if(fNHistoBinsXi[i]) hRespXi[i]->SetBinEdges(1,fHistoBinsXi[i]->GetArray());
3210     //     }
3211
3212
3213   } // jet slices loop
3214
3215   f.Close();  // THnSparse pointers still valid even if file closed
3216
3217 //   for(Int_t i=0; i<fNJetPtSlices; i++){  // no copy c'tor ...
3218 //     if(hRespPt[i]) new(fhnResponsePt[i]) THnSparseF(*hRespPt[i]);
3219 //     if(hRespZ[i])  new(fhnResponseZ[i])  THnSparseF(*hRespZ[i]);
3220 //     if(hRespXi[i]) new(fhnResponseXi[i]) THnSparseF(*hRespXi[i]);
3221 //   }
3222
3223   for(Int_t i=0; i<fNJetPtSlices; i++){
3224     fhnResponsePt[i] = hRespPt[i];
3225     fhnResponseZ[i]  = hRespZ[i];
3226     fhnResponseXi[i] = hRespXi[i];
3227   }
3228 }
3229
3230 //______________________________________________________________________________________________________________________
3231 void AliFragmentationFunctionCorrections::ReadPriors(TString strfile,const Int_t type)
3232 {
3233   // read priors from file: rec primaries, gen pt dist
3234
3235   if(fDebug>0) Printf("%s:%d -- read priors from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
3236
3237   // temporary histos to store pointers from file
3238   TH1F* hist[fNJetPtSlices];
3239
3240   for(Int_t i=0; i<fNJetPtSlices; i++) hist[i] = 0;
3241
3242   TFile f(strfile,"READ");
3243
3244   if(!f.IsOpen()){
3245     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
3246     return;
3247   }
3248
3249   for(Int_t i=0; i<fNJetPtSlices; i++){
3250
3251     Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
3252     Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
3253
3254     TString strName;
3255
3256     if(type == kFlagPt) strName.Form("h1FFTrackPtRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim);
3257     if(type == kFlagZ)  strName.Form("h1FFZRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim);
3258     if(type == kFlagXi) strName.Form("h1FFXiRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim);
3259
3260     hist[i] = (TH1F*) gDirectory->Get(strName);
3261
3262     if(!hist[i]){
3263       Printf("%s:%d -- error retrieving prior %s", (char*)__FILE__,__LINE__,strName.Data());
3264     }
3265
3266
3267     //if(fNHistoBinsPt[i]) hist[i] = (TH1F*) hist[i]->Rebin(fNHistoBinsPt[i],hist[i]->GetName()+"_rebin",fHistoBinsPt[i]->GetArray());
3268
3269     if(hist[i]) hist[i]->SetDirectory(0);
3270
3271   } // jet slices loop
3272
3273   f.Close();
3274
3275
3276   for(Int_t i=0; i<fNJetPtSlices; i++){ // 2nd loop: need to close input file before placing histos
3277     if(hist[i] && type == kFlagPt) new(fh1FFTrackPtPrior[i]) TH1F(*hist[i]);
3278     if(hist[i] && type == kFlagZ)  new(fh1FFZPrior[i])       TH1F(*hist[i]);
3279     if(hist[i] && type == kFlagXi) new(fh1FFXiPrior[i])      TH1F(*hist[i]);
3280   }
3281 }
3282
3283 //_____________________________________________________
3284 // void AliFragmentationFunctionCorrections::RatioRecGen()
3285 // {
3286 //   // create ratio reconstructed over generated FF
3287 //   // use current highest corrLevel
3288
3289 //   Printf("%s:%d -- build ratio rec.gen, corrLevel %d",(char*)__FILE__,__LINE__,fNCorrectionLevels-1);
3290
3291 //   for(Int_t i=0; i<fNJetPtSlices; i++){
3292
3293 //     TH1F* histPtRec =  fCorrFF[fNCorrectionLevels-1]->GetTrackPt(i); // levels -1: latest corr level
3294 //     TH1F* histZRec  =  fCorrFF[fNCorrectionLevels-1]->GetZ(i);       // levels -1: latest corr level
3295 //     TH1F* histXiRec =  fCorrFF[fNCorrectionLevels-1]->GetXi(i);      // levels -1: latest corr level
3296
3297 //     TH1F* histPtGen = fh1FFTrackPtGenPrim[i];
3298 //     TH1F* histZGen  = fh1FFZGenPrim[i];
3299 //     TH1F* histXiGen = fh1FFXiGenPrim[i];
3300
3301 //     TString histNamePt = histPtRec->GetName();
3302 //     TString histNameZ  = histZRec->GetName();
3303 //     TString histNameXi = histXiRec->GetName();
3304
3305 //     histNamePt.ReplaceAll("fh1FF","fh1FFRatioRecGen");
3306 //     histNameZ.ReplaceAll("fh1FF","fh1FFRatioRecGen");
3307 //     histNameXi.ReplaceAll("fh1FF","fh1FFRatioRecGen");
3308
3309 //     // ratio
3310 //     TH1F* hRatioRecGenPt = (TH1F*) histPtRec->Clone(histNamePt);
3311 //     hRatioRecGenPt->Reset();
3312 //     hRatioRecGenPt->Divide(histPtRec,histPtGen,1,1,"B");
3313
3314 //     TH1F* hRatioRecGenZ = (TH1F*) histZRec->Clone(histNameZ);
3315 //     hRatioRecGenZ->Reset();
3316 //     hRatioRecGenZ->Divide(histZRec,histZGen,1,1,"B");
3317
3318 //     TH1F* hRatioRecGenXi = (TH1F*) histXiRec->Clone(histNameXi);
3319 //     hRatioRecGenXi->Reset();
3320 //     hRatioRecGenXi->Divide(histXiRec,histXiGen,1,1,"B");
3321
3322 //     new(fh1FFRatioRecGenPt[i]) TH1F(*hRatioRecGenPt);
3323 //     new(fh1FFRatioRecGenZ[i])  TH1F(*hRatioRecGenZ);
3324 //     new(fh1FFRatioRecGenXi[i]) TH1F(*hRatioRecGenXi);
3325 //   }
3326 // }
3327
3328 // //___________________________________________________________
3329 // void AliFragmentationFunctionCorrections::RatioRecPrimaries()
3330 // {
3331 //   // create ratio reconstructed tracks over reconstructed primaries
3332 //   // use raw FF (corrLevel 0)
3333
3334 //   Printf("%s:%d -- build ratio rec tracks /rec primaries",(char*)__FILE__,__LINE__);
3335
3336 //   for(Int_t i=0; i<fNJetPtSlices; i++){
3337
3338 //     const Int_t corrLevel = 0;
3339
3340 //     TH1F* histPtRec =  fCorrFF[corrLevel]->GetTrackPt(i); // levels -1: latest corr level
3341 //     TH1F* histZRec  =  fCorrFF[corrLevel]->GetZ(i);       // levels -1: latest corr level
3342 //     TH1F* histXiRec =  fCorrFF[corrLevel]->GetXi(i);      // levels -1: latest corr level
3343
3344 //     TH1F* histPtRecPrim = fh1FFTrackPtRecPrim[i];
3345 //     TH1F* histZRecPrim  = fh1FFZRecPrim[i];
3346 //     TH1F* histXiRecPrim = fh1FFXiRecPrim[i];
3347
3348 //     TString histNamePt = histPtRec->GetName();
3349 //     TString histNameZ  = histZRec->GetName();
3350 //     TString histNameXi = histXiRec->GetName();
3351
3352 //     histNamePt.ReplaceAll("fh1FF","fh1FFRatioRecPrim");
3353 //     histNameZ.ReplaceAll("fh1FF","fh1FFRatioRecPrim");
3354 //     histNameXi.ReplaceAll("fh1FF","fh1FFRatioRecPrim");
3355
3356 //     // ratio
3357 //     TH1F* hRatioRecPrimPt = (TH1F*) histPtRec->Clone(histNamePt);
3358 //     hRatioRecPrimPt->Reset();
3359 //     hRatioRecPrimPt->Divide(histPtRec,histPtRecPrim,1,1,"B");
3360
3361 //     TH1F* hRatioRecPrimZ = (TH1F*) histZRec->Clone(histNameZ);
3362 //     hRatioRecPrimZ->Reset();
3363 //     hRatioRecPrimZ->Divide(histZRec,histZRecPrim,1,1,"B");
3364
3365 //     TH1F* hRatioRecPrimXi = (TH1F*) histXiRec->Clone(histNameXi);
3366 //     hRatioRecPrimXi->Reset();
3367 //     hRatioRecPrimXi->Divide(histXiRec,histXiRecPrim,1,1,"B");
3368
3369
3370 //     new(fh1FFRatioRecPrimPt[i]) TH1F(*hRatioRecPrimPt);
3371 //     new(fh1FFRatioRecPrimZ[i])  TH1F(*hRatioRecPrimZ);
3372 //     new(fh1FFRatioRecPrimXi[i]) TH1F(*hRatioRecPrimXi);
3373 //   }
3374 // }
3375
3376 // __________________________________________________________________________________
3377 void AliFragmentationFunctionCorrections::ProjectJetResponseMatrices(TString strOutfile)
3378 {
3379
3380   // project response matrices on both axes:
3381   // FF for rec primaries, in terms of generated and reconstructed momentum
3382   // write FF and ratios to outFile
3383
3384   Printf("%s:%d -- project response matrices, write to %s",(char*)__FILE__,__LINE__,strOutfile.Data());
3385
3386   TH1F* hFFPtRec[fNJetPtSlices];
3387   TH1F* hFFZRec[fNJetPtSlices];
3388   TH1F* hFFXiRec[fNJetPtSlices];
3389
3390   TH1F* hFFPtGen[fNJetPtSlices];
3391   TH1F* hFFZGen[fNJetPtSlices];
3392   TH1F* hFFXiGen[fNJetPtSlices];
3393
3394   TH1F* hRatioPt[fNJetPtSlices];
3395   TH1F* hRatioZ[fNJetPtSlices];
3396   TH1F* hRatioXi[fNJetPtSlices];
3397
3398
3399   Int_t axisGenPt = 1;
3400   Int_t axisRecPt = 0;
3401
3402   for(Int_t i=0; i<fNJetPtSlices; i++){
3403
3404     Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
3405     Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
3406
3407     TString strNameRecPt(Form("h1FFTrackPtRecPrim_recPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3408     TString strNameRecZ(Form("h1FFZRecPrim_recPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3409     TString strNameRecXi(Form("h1FFXiRecPrim_recPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3410
3411     TString strNameGenPt(Form("h1FFTrackPtRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3412     TString strNameGenZ(Form("h1FFZRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3413     TString strNameGenXi(Form("h1FFXiRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3414
3415     TString strNameRatioPt(Form("h1RatioTrackPtRecPrim_%02d_%02d",jetPtLoLim,jetPtUpLim));
3416     TString strNameRatioZ(Form("h1RatioZRecPrim_%02d_%02d",jetPtLoLim,jetPtUpLim));
3417     TString strNameRatioXi(Form("h1RatioXiRecPrim_%02d_%02d",jetPtLoLim,jetPtUpLim));
3418
3419
3420     hFFPtRec[i] = (TH1F*) fhnResponsePt[i]->Projection(axisRecPt);// note convention: yDim,xDim
3421     hFFZRec[i]  = (TH1F*) fhnResponseZ[i]->Projection(axisRecPt);// note convention: yDim,xDim
3422     hFFXiRec[i] = (TH1F*) fhnResponseXi[i]->Projection(axisRecPt);// note convention: yDim,xDim
3423
3424     hFFPtRec[i]->SetNameTitle(strNameRecPt,"");
3425     hFFZRec[i]->SetNameTitle(strNameRecZ,"");
3426     hFFXiRec[i]->SetNameTitle(strNameRecXi,"");
3427
3428
3429     hFFPtGen[i] = (TH1F*) fhnResponsePt[i]->Projection(axisGenPt);// note convention: yDim,xDim
3430     hFFZGen[i]  = (TH1F*) fhnResponseZ[i]->Projection(axisGenPt);// note convention: yDim,xDim
3431     hFFXiGen[i] = (TH1F*) fhnResponseXi[i]->Projection(axisGenPt);// note convention: yDim,xDim
3432
3433     hFFPtGen[i]->SetNameTitle(strNameGenPt,"");
3434     hFFZGen[i]->SetNameTitle(strNameGenZ,"");
3435     hFFXiGen[i]->SetNameTitle(strNameGenXi,"");
3436
3437
3438     if(fNHistoBinsPt[i]) hFFPtRec[i] = (TH1F*) hFFPtRec[i]->Rebin(fNHistoBinsPt[i],"",fHistoBinsPt[i]->GetArray());
3439     if(fNHistoBinsZ[i])  hFFZRec[i]   = (TH1F*) hFFZRec[i]->Rebin(fNHistoBinsZ[i],"",fHistoBinsZ[i]->GetArray());
3440     if(fNHistoBinsXi[i]) hFFXiRec[i]  = (TH1F*) hFFXiRec[i]->Rebin(fNHistoBinsXi[i],"",fHistoBinsXi[i]->GetArray());
3441
3442     if(fNHistoBinsPt[i]) hFFPtGen[i] = (TH1F*) hFFPtGen[i]->Rebin(fNHistoBinsPt[i],"",fHistoBinsPt[i]->GetArray());
3443     if(fNHistoBinsZ[i])  hFFZGen[i]  = (TH1F*) hFFZGen[i]->Rebin(fNHistoBinsZ[i],"",fHistoBinsZ[i]->GetArray());
3444     if(fNHistoBinsXi[i]) hFFXiGen[i] = (TH1F*) hFFXiGen[i]->Rebin(fNHistoBinsXi[i],"",fHistoBinsXi[i]->GetArray());
3445
3446     NormalizeTH1(hFFPtGen[i],fNJets->At(i));
3447     NormalizeTH1(hFFZGen[i],fNJets->At(i));
3448     NormalizeTH1(hFFXiGen[i],fNJets->At(i));
3449
3450     NormalizeTH1(hFFPtRec[i],fNJets->At(i));
3451     NormalizeTH1(hFFZRec[i],fNJets->At(i));
3452     NormalizeTH1(hFFXiRec[i],fNJets->At(i));
3453
3454
3455     hRatioPt[i] = (TH1F*) hFFPtRec[i]->Clone(strNameRatioPt);
3456     hRatioPt[i]->Reset();
3457     hRatioPt[i]->Divide(hFFPtRec[i],hFFPtGen[i],1,1,"B");
3458
3459     hRatioZ[i] = (TH1F*) hFFZRec[i]->Clone(strNameRatioZ);
3460     hRatioZ[i]->Reset();
3461     hRatioZ[i]->Divide(hFFZRec[i],hFFZGen[i],1,1,"B");
3462
3463     hRatioXi[i] = (TH1F*) hFFXiRec[i]->Clone(strNameRatioXi);
3464     hRatioXi[i]->Reset();
3465     hRatioXi[i]->Divide(hFFXiRec[i],hFFXiGen[i],1,1,"B");
3466   }
3467
3468
3469
3470   // write
3471
3472   TFile out(strOutfile,"RECREATE");
3473
3474   if(!out.IsOpen()){
3475     Printf("%s:%d -- error opening response matrix output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
3476     return;
3477   }
3478
3479   for(Int_t i=0; i<fNJetPtSlices; i++){
3480
3481     hFFPtRec[i]->Write();
3482     hFFZRec[i]->Write();
3483     hFFXiRec[i]->Write();
3484
3485     hFFPtGen[i]->Write();
3486     hFFZGen[i]->Write();
3487     hFFXiGen[i]->Write();
3488
3489     hRatioPt[i]->Write();
3490     hRatioZ[i]->Write();
3491     hRatioXi[i]->Write();
3492   }
3493
3494   out.Close();
3495 }
3496
3497 // ____________________________________________________________________________________________________________________________
3498 void AliFragmentationFunctionCorrections::ProjectSingleResponseMatrix(TString strOutfile, Bool_t updateOutfile, TString strOutDir )
3499 {
3500   // project response matrix on both axes:
3501   // pt spec for rec primaries, in terms of generated and reconstructed momentum
3502   // write spec and ratios to outFile
3503
3504   Printf("%s:%d -- project single pt response matrix, write to %s",(char*)__FILE__,__LINE__,strOutfile.Data());
3505
3506   TH1F* hSpecPtRec;
3507   TH1F* hSpecPtGen;
3508   TH1F* hRatioPt;
3509
3510   Int_t axisGenPt = 1;
3511   Int_t axisRecPt = 0;
3512
3513   TString strNameRecPt   = "h1SpecTrackPtRecPrim_recPt";
3514   TString strNameGenPt   = "h1SpecTrackPtRecPrim_genPt";
3515   TString strNameRatioPt = "h1RatioTrackPtRecPrim";
3516
3517   hSpecPtRec = (TH1F*) fhnResponseSinglePt->Projection(axisRecPt);// note convention: yDim,xDim
3518   hSpecPtRec->SetNameTitle(strNameRecPt,"");
3519
3520   hSpecPtGen = (TH1F*) fhnResponseSinglePt->Projection(axisGenPt);// note convention: yDim,xDim
3521   hSpecPtGen->SetNameTitle(strNameGenPt,"");
3522
3523   hRatioPt = (TH1F*) hSpecPtRec->Clone(strNameRatioPt);
3524   hRatioPt->Reset();
3525   hRatioPt->Divide(hSpecPtRec,hSpecPtGen,1,1,"B");
3526
3527   TString outfileOption = "RECREATE";
3528   if(updateOutfile) outfileOption = "UPDATE";
3529
3530   TFile out(strOutfile,outfileOption);
3531
3532   if(!out.IsOpen()){
3533     Printf("%s:%d -- error opening reponse matrix projections output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
3534     return;
3535   }
3536
3537
3538   if(strOutDir && strOutDir.Length()){
3539
3540     TDirectory* dir;
3541     if((dir = ((TDirectory*) gDirectory->Get(strOutDir)))) dir->cd();
3542     else{
3543       dir = out.mkdir(strOutDir);
3544       dir->cd();
3545     }
3546   }
3547
3548   hSpecPtRec->Write();
3549   hSpecPtGen->Write();
3550   hRatioPt->Write();
3551
3552   out.Close();
3553 }
3554
3555
3556 //__________________________________________________________________________________________________________________________________________________________________
3557 void AliFragmentationFunctionCorrections::RebinHisto(const Int_t jetPtSlice, const Int_t nBinsLimits, Double_t* binsLimits, Double_t* binsWidth, const Int_t type)
3558 {
3559   // rebin histo, rescale bins according to new width
3560   // only correct for input histos with equal bin size
3561
3562   // args: jetPtSlice, type, use current corr level
3563
3564   // function takes array of limits and widths (e.g. 1 GeV bins up to 10 GeV, 2 GeV width up to 20 GeV, ...)
3565   // array size of binsLimits: nBinsLimits
3566   // array size of binsWidth: nBinsLimits-1
3567   // binsLimits have to be in increasing order
3568   // if binning undefined for any slice, original binning will be kept
3569
3570   if(!fNJetPtSlices){
3571     Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
3572     return;
3573   }
3574
3575   if(jetPtSlice>=fNJetPtSlices){
3576     Printf("%s:%d -- jetPtSlice %d exceeds max",(char*)__FILE__,__LINE__,jetPtSlice);
3577     return;
3578   }
3579
3580
3581   Double_t binLimitMin = binsLimits[0];
3582   Double_t binLimitMax = binsLimits[nBinsLimits-1];
3583
3584   Double_t binLimit = binLimitMin; // start value
3585
3586   Int_t sizeUpperLim = 1000; //static_cast<Int_t>(binLimitMax/binsWidth[0])+1; -  only if first bin has smallest width, but not case for dN/dxi ...
3587   TArrayD binsArray(sizeUpperLim);
3588   Int_t nBins = 0;
3589   binsArray.SetAt(binLimitMin,nBins++);
3590
3591   while(binLimit<binLimitMax && nBins<sizeUpperLim){
3592
3593     Int_t currentSlice = -1;
3594     for(Int_t i=0; i<nBinsLimits; i++){
3595       if(binLimit >= 0.999*binsLimits[i]) currentSlice = i; // 0.999 numerical saftey factor
3596     }
3597
3598     Double_t currentBinWidth = binsWidth[currentSlice];
3599     binLimit += currentBinWidth;
3600
3601     binsArray.SetAt(binLimit,nBins++);
3602   }
3603
3604
3605   TH1F* hist = 0;
3606   if(type == kFlagPt)       hist = fCorrFF[fNCorrectionLevels-1]->GetTrackPt(jetPtSlice);
3607   else if(type == kFlagZ)        hist = fCorrFF[fNCorrectionLevels-1]->GetZ(jetPtSlice);
3608   else if(type == kFlagXi)       hist = fCorrFF[fNCorrectionLevels-1]->GetXi(jetPtSlice);
3609   else if(type == kFlagSinglePt) hist = fCorrSinglePt[fNCorrectionLevelsSinglePt-1]->GetTrackPt(0);
3610   else{
3611     Printf("%s%d unknown type",(char*)__FILE__,__LINE__);
3612     return;
3613   }
3614
3615   Double_t binWidthNoRebin = hist->GetBinWidth(1);
3616
3617   Double_t* bins = binsArray.GetArray();
3618
3619   hist = (TH1F*) hist->Rebin(nBins-1,"",bins);
3620
3621   for(Int_t bin=0; bin <= hist->GetNbinsX(); bin++){
3622
3623     Double_t binWidthRebin = hist->GetBinWidth(bin);
3624     Double_t scaleF = binWidthNoRebin / binWidthRebin;
3625
3626     Double_t binCont  = hist->GetBinContent(bin);
3627     Double_t binErr   = hist->GetBinError(bin);
3628
3629     binCont *= scaleF;
3630     binErr  *= scaleF;
3631
3632     hist->SetBinContent(bin,binCont);
3633     hist->SetBinError(bin,binErr);
3634   }
3635
3636
3637
3638   TH1F* temp = new TH1F(*hist);
3639
3640   if(type == kFlagPt)       fCorrFF[fNCorrectionLevels-1]->ReplaceCorrHistos(jetPtSlice,temp,0,0);
3641   if(type == kFlagZ)        fCorrFF[fNCorrectionLevels-1]->ReplaceCorrHistos(jetPtSlice,0,temp,0);
3642   if(type == kFlagXi)       fCorrFF[fNCorrectionLevels-1]->ReplaceCorrHistos(jetPtSlice,0,0,temp);
3643   if(type == kFlagSinglePt) fCorrSinglePt[fNCorrectionLevelsSinglePt-1]->ReplaceCorrHistos(0,temp,0,0);
3644
3645
3646   delete temp;
3647 }
3648 //__________________________________________________________________________________________________________________________________________________________________
3649 void AliFragmentationFunctionCorrections::WriteJetSpecResponse(TString strInfile, TString strdir, TString strlist/*, TString strOutfile*/)
3650 {
3651
3652  if(fDebug>0) Printf("%s:%d -- read jet spectrum response matrix from file %s ",(char*)__FILE__,__LINE__,strInfile.Data());
3653
3654   if(strdir && strdir.Length()) gDirectory->cd(strdir);
3655
3656   TList* list = 0;
3657
3658   if(strlist && strlist.Length()){
3659
3660     if(!(list = (TList*) gDirectory->Get(strlist))){
3661       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
3662       return;
3663     }
3664   }
3665
3666   THnSparse* hn6ResponseJetPt;
3667
3668   if(list){
3669     hn6ResponseJetPt = (THnSparse*) list->FindObject("fhnCorrelation");
3670   }
3671   else{
3672     hn6ResponseJetPt = (THnSparse*) list->FindObject("fhnCorrelation");
3673   }
3674
3675   Int_t axis6RecJetPt = 0;
3676   Int_t axis6GenJetPt = 3;
3677
3678   hn6ResponseJetPt->GetAxis(axis6RecJetPt)->SetTitle("rec jet p_{T} (GeV/c)");
3679   hn6ResponseJetPt->GetAxis(axis6GenJetPt)->SetTitle("gen jet p_{T} (GeV/c)");
3680
3681   Int_t nBinsRecPt    = hn6ResponseJetPt->GetAxis(axis6RecJetPt)->GetNbins();
3682   Double_t loLimRecPt = hn6ResponseJetPt->GetAxis(axis6RecJetPt)->GetBinLowEdge(1);
3683   Double_t upLimRecPt = hn6ResponseJetPt->GetAxis(axis6RecJetPt)->GetBinUpEdge(nBinsRecPt);
3684
3685   Int_t nBinsGenPt    = hn6ResponseJetPt->GetAxis(axis6GenJetPt)->GetNbins();
3686   Double_t loLimGenPt = hn6ResponseJetPt->GetAxis(axis6GenJetPt)->GetBinLowEdge(1);
3687   Double_t upLimGenPt = hn6ResponseJetPt->GetAxis(axis6GenJetPt)->GetBinUpEdge(nBinsGenPt);
3688
3689   Int_t nBinsTrackPt = 200;
3690   Int_t loLimTrackPt = 0;
3691   Int_t upLimTrackPt = 200;
3692
3693
3694   Int_t    nBinsResponse[4]  = {nBinsRecPt,nBinsTrackPt,nBinsGenPt,nBinsTrackPt};
3695   Double_t binMinResponse[4] = {loLimRecPt,loLimTrackPt,loLimGenPt,loLimTrackPt};
3696   Double_t binMaxResponse[4] = {upLimRecPt,upLimTrackPt,upLimGenPt,upLimTrackPt};
3697
3698   const char* labelsResponseSinglePt[4] = {"rec jet p_{T} (GeV/c)", "rec track p_{T} (GeV/c)", "gen jet p_{T} (GeV/c)", "gen track p_{T} (GeV/c)"};
3699
3700   THnSparseD* hn4ResponseTrackPtJetPt =  new THnSparseD("hn4ResponseTrackPtJetPt","",4,nBinsResponse,binMinResponse,binMaxResponse);
3701
3702   for(Int_t i=0; i<4; i++){
3703     hn4ResponseTrackPtJetPt->GetAxis(i)->SetTitle(labelsResponseSinglePt[i]);
3704   }
3705
3706
3707   // fill
3708
3709
3710 }
3711
3712 //_____________________________________________________________________________________________________________________________________
3713 void AliFragmentationFunctionCorrections::ReadSingleTrackEfficiency(TString strfile, TString strdir, TString strlist, TString strname)
3714 {
3715
3716   ReadSingleTrackCorrection(strfile,strdir,strlist,strname,kFlagEfficiency);
3717
3718 }
3719
3720 //_____________________________________________________________________________________________________________________________________
3721 void AliFragmentationFunctionCorrections::ReadSingleTrackResponse(TString strfile, TString strdir, TString strlist, TString strname)
3722 {
3723
3724   ReadSingleTrackCorrection(strfile,strdir,strlist,strname,kFlagResponse);
3725
3726 }
3727
3728 //_____________________________________________________________________________________________________________________________________
3729 void AliFragmentationFunctionCorrections::ReadSingleTrackSecCorr(TString strfile, TString strdir, TString strlist, TString strname)
3730 {
3731
3732   ReadSingleTrackCorrection(strfile,strdir,strlist,strname,kFlagSecondaries);
3733
3734 }
3735
3736 //______________________________________________________________________________________________________________________________________________________
3737 void AliFragmentationFunctionCorrections::ReadSingleTrackCorrection(TString strfile, TString strdir, TString strlist, TString strname, const Int_t type)
3738 {
3739   // read single track correction (pt) from file
3740   // type: efficiency / response / secondaries correction
3741
3742   if(!((type == kFlagEfficiency)  || (type == kFlagResponse) || (type == kFlagSecondaries))){
3743     Printf("%s:%d -- no such correction ",(char*)__FILE__,__LINE__);
3744     return;
3745   }
3746
3747   TFile f(strfile,"READ");
3748
3749   if(!f.IsOpen()){
3750     Printf("%s:%d -- error opening file %s", (char*)__FILE__,__LINE__,strfile.Data());
3751     return;
3752   }
3753
3754   if(fDebug>0 && type==kFlagEfficiency)  Printf("%s:%d -- read single track corr from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
3755   if(fDebug>0 && type==kFlagResponse)    Printf("%s:%d -- read single track response from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
3756   if(fDebug>0 && type==kFlagSecondaries) Printf("%s:%d -- read single track secondaries corr from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
3757
3758   if(strdir && strdir.Length()) gDirectory->cd(strdir);
3759
3760   TList* list = 0;
3761
3762   if(strlist && strlist.Length()){
3763
3764     if(!(list = (TList*) gDirectory->Get(strlist))){
3765       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
3766       return;
3767     }
3768   }
3769
3770   TH1F*      h1CorrHist = 0; // common TObject pointer not possible, need SetDirectory() later
3771   THnSparse* hnCorrHist = 0;
3772
3773   if(type == kFlagEfficiency || type == kFlagSecondaries){
3774
3775     if(list) h1CorrHist = (TH1F*) list->FindObject(strname);
3776     else     h1CorrHist = (TH1F*) gDirectory->Get(strname);
3777
3778     if(!h1CorrHist){
3779       Printf("%s:%d -- error retrieving histo %s", (char*)__FILE__,__LINE__,strname.Data());
3780       return;
3781     }
3782
3783   }
3784   else if(type == kFlagResponse){
3785
3786     if(list) hnCorrHist = (THnSparse*) list->FindObject(strname);
3787     else     hnCorrHist = (THnSparse*) gDirectory->Get(strname);
3788
3789     if(!hnCorrHist){
3790       Printf("%s:%d -- error retrieving histo %s", (char*)__FILE__,__LINE__,strname.Data());
3791       return;
3792     }
3793
3794   }
3795
3796   if(h1CorrHist) h1CorrHist->SetDirectory(0);
3797   //if(hnCorrHist) hnCorrHist->SetDirectory(0);
3798
3799   f.Close();
3800
3801   if(type == kFlagEfficiency)       fh1EffSinglePt      = h1CorrHist;
3802   else if(type == kFlagResponse)    fhnResponseSinglePt = hnCorrHist;
3803   else if(type == kFlagSecondaries) fh1SecCorrSinglePt  = h1CorrHist;
3804
3805 }
3806
3807 //________________________________________________________________________________________________________________
3808 void AliFragmentationFunctionCorrections::ReadRawPtSpec(TString strInfile, TString strID)
3809 {
3810   // read track pt spec from task ouput - standard dir/list
3811
3812   TString strdir  = "PWG4_FragmentationFunction_" + strID;
3813   TString strlist = "fracfunc_" + strID;
3814
3815   ReadRawPtSpec(strInfile,strdir,strlist);
3816 }
3817
3818 //_______________________________________________________________________________________________________
3819 void AliFragmentationFunctionCorrections::ReadRawPtSpec(TString strfile, TString strdir, TString strlist)
3820 {
3821   // get raw pt spectra from file
3822
3823
3824   // book histos
3825   fNCorrectionLevelsSinglePt = 0;
3826   fCorrSinglePt = new AliFragFuncCorrHistos*[fgMaxNCorrectionLevels];
3827   AddCorrectionLevelSinglePt(); // first 'correction' level = raw spectrum
3828
3829   // get raw pt spec from input file, normalize
3830
3831   TFile f(strfile,"READ");
3832
3833   if(!f.IsOpen()){
3834     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
3835     return;
3836   }
3837
3838   if(fDebug>0) Printf("%s:%d -- read raw spectra from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
3839
3840   gDirectory->cd(strdir);
3841
3842   TList* list = 0;
3843
3844   if(!(list = (TList*) gDirectory->Get(strlist))){
3845     Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
3846     return;
3847   }
3848
3849   TString hnameTrackPt("fh1TrackQAPtRecCuts");
3850   TString hnameEvtSel("fh1EvtSelection");
3851
3852   TH1F* fh1TrackPt = (TH1F*) list->FindObject(hnameTrackPt);
3853   TH1F* fh1EvtSel  = (TH1F*) list->FindObject(hnameEvtSel);
3854
3855   if(!fh1TrackPt){ Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameTrackPt.Data()); return; }
3856   if(!fh1EvtSel) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameEvtSel.Data()); return; }
3857
3858
3859   //Float_t nEvents = fh1EvtSel->GetBinContent(fh1EvtSel->FindBin(0));
3860
3861   // evts after physics selection
3862   Float_t nEvents = fh1EvtSel->GetEntries()
3863     - fh1EvtSel->GetBinContent(fh1EvtSel->FindBin(1))   // evts rejected by trigger selection ( = PhysicsSelection
3864     - fh1EvtSel->GetBinContent(fh1EvtSel->FindBin(2));   // evts with wrong centrality class
3865
3866
3867   fh1TrackPt->SetDirectory(0);
3868
3869   f.Close();
3870
3871
3872   NormalizeTH1(fh1TrackPt,nEvents);
3873
3874   // raw FF = corr level 0
3875   fCorrSinglePt[0]->AddCorrHistos(0,fh1TrackPt);
3876 }
3877
3878
3879 //_______________________________________________________________________________________________________
3880 void AliFragmentationFunctionCorrections::ReadRawPtSpecQATask(TString strfile, TString strdir, TString strlist)
3881 {
3882   // get raw pt spectra from file
3883   // for output from Martas QA task
3884
3885
3886   // book histos
3887   fNCorrectionLevelsSinglePt = 0;
3888   fCorrSinglePt = new AliFragFuncCorrHistos*[fgMaxNCorrectionLevels];
3889   AddCorrectionLevelSinglePt(); // first 'correction' level = raw spectrum
3890
3891   // get raw pt spec from input file, normalize
3892
3893   TFile f(strfile,"READ");
3894
3895   if(!f.IsOpen()){
3896     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
3897     return;
3898   }
3899
3900   if(fDebug>0) Printf("%s:%d -- read raw pt spec from QA task output file %s ",(char*)__FILE__,__LINE__,strfile.Data());
3901
3902   gDirectory->cd(strdir);
3903
3904   TList* list = 0;
3905
3906   if(!(list = (TList*) gDirectory->Get(strlist))){
3907     Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
3908     return;
3909   }
3910
3911   TString hnameTrackPt("fPtSel");
3912   TString hnameEvtSel("fNEventAll");
3913
3914   TH1F* fh1TrackPt = (TH1F*) list->FindObject(hnameTrackPt);
3915   TH1F* fh1EvtSel  = (TH1F*) list->FindObject(hnameEvtSel);
3916
3917   if(!fh1TrackPt){ Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameTrackPt.Data()); return; }
3918   if(!fh1EvtSel) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameEvtSel.Data()); return; }
3919
3920
3921   // evts after physics selection
3922   Float_t nEvents = fh1EvtSel->GetEntries();
3923
3924   fh1TrackPt->SetDirectory(0);
3925
3926   f.Close();
3927
3928
3929   NormalizeTH1(fh1TrackPt,nEvents);
3930
3931   // raw FF = corr level 0
3932   fCorrSinglePt[0]->AddCorrHistos(0,fh1TrackPt);
3933 }
3934
3935 // ________________________________________________________
3936 void AliFragmentationFunctionCorrections::EffCorrSinglePt()
3937 {
3938   // apply efficiency correction to inclusive track pt spec
3939
3940   AddCorrectionLevelSinglePt("eff");
3941
3942   TH1F* histPt = fCorrSinglePt[fNCorrectionLevelsSinglePt-2]->GetTrackPt(0);
3943
3944   if(histPt->GetNbinsX() != fh1EffSinglePt->GetNbinsX()) Printf("%s:%d: inconsistency pt spec and eff corr bins - rebin effCorr ...", (char*)__FILE__,__LINE__);
3945
3946   TString histNamePt = histPt->GetName();
3947   TH1F* hTrackPtEffCorr = (TH1F*) histPt->Clone(histNamePt);
3948   hTrackPtEffCorr->Divide(histPt,fh1EffSinglePt,1,1,"");
3949
3950   fCorrSinglePt[fNCorrectionLevelsSinglePt-1]->AddCorrHistos(0,hTrackPtEffCorr);
3951 }
3952
3953 //___________________________________________________________________________________________________________________________
3954 void AliFragmentationFunctionCorrections::UnfoldSinglePt(const Int_t nIter, const Bool_t useCorrelatedErrors)
3955 {
3956   // unfolde inclusive dN/dpt spectra
3957
3958   AddCorrectionLevelSinglePt("unfold");
3959
3960   TH1F* hist = fCorrSinglePt[fNCorrectionLevelsSinglePt-2]->GetTrackPt(0); // level -2: before unfolding, level -1: unfolded
3961   THnSparse* hnResponse = fhnResponseSinglePt;
3962
3963   TString histNameTHn = hist->GetName();
3964   if(histNameTHn.Contains("TH1")) histNameTHn.ReplaceAll("TH1","THn");
3965   if(histNameTHn.Contains("fPt")) histNameTHn.ReplaceAll("fPt","fhnPt");
3966
3967
3968   TString histNameBackFolded = hist->GetName();
3969   histNameBackFolded.Append("_backfold");
3970
3971   TString histNameRatioFolded = hist->GetName();
3972   if(histNameRatioFolded.Contains("fh1")) histNameRatioFolded.ReplaceAll("fh1","hRatio");
3973   if(histNameRatioFolded.Contains("fPt")) histNameRatioFolded.ReplaceAll("fPt","hRatioPt");
3974   histNameRatioFolded.Append("_unfold");
3975
3976   TString histNameRatioBackFolded = hist->GetName();
3977   if(histNameRatioBackFolded.Contains("fh1")) histNameRatioBackFolded.ReplaceAll("fh1","hRatio");
3978   if(histNameRatioBackFolded.Contains("fPt")) histNameRatioBackFolded.ReplaceAll("fPt","hRatioPt");
3979   histNameRatioBackFolded.Append("_backfold");
3980
3981   THnSparse* hnHist           = TH1toSparse(hist,histNameTHn,hist->GetTitle());
3982   THnSparse* hnFlatEfficiency = TH1toSparse(hist,"fhnEfficiency","eff",kTRUE); // could optionally also use real eff
3983
3984   THnSparse* hnUnfolded
3985     = Unfold(hnHist,hnResponse,hnFlatEfficiency,nIter,useCorrelatedErrors);
3986
3987   TH1F* hUnfolded = (TH1F*) hnUnfolded->Projection(0);
3988   hUnfolded->SetNameTitle(hist->GetName(),hist->GetTitle());
3989
3990   fCorrSinglePt[fNCorrectionLevelsSinglePt-1]->AddCorrHistos(0,hUnfolded);
3991
3992   // backfolding: apply response matrix to unfolded spectrum
3993   TH1F* hBackFolded = ApplyResponse(hUnfolded,hnResponse);
3994   hBackFolded->SetNameTitle(histNameBackFolded,hist->GetTitle());
3995
3996   fh1SingleTrackPtBackFolded = hBackFolded;
3997
3998
3999   // ratio unfolded to original histo
4000   TH1F* hRatioUnfolded = (TH1F*) hUnfolded->Clone(histNameRatioFolded);
4001   hRatioUnfolded->Reset();
4002   hRatioUnfolded->Divide(hUnfolded,hist,1,1,"B");
4003
4004   fh1RatioSingleTrackPtFolded = hRatioUnfolded;
4005
4006
4007   // ratio backfolded to original histo
4008   TH1F* hRatioBackFolded = (TH1F*) hBackFolded->Clone(histNameRatioBackFolded);
4009   hRatioBackFolded->Reset();
4010   hRatioBackFolded->Divide(hBackFolded,hist,1,1,"B");
4011
4012   fh1RatioSingleTrackPtBackFolded = hRatioBackFolded;
4013
4014   delete hnHist;
4015   delete hnFlatEfficiency;
4016
4017 }
4018
4019 // ________________________________________________________
4020 void AliFragmentationFunctionCorrections::SecCorrSinglePt()
4021 {
4022   // apply efficiency correction to inclusive track pt spec
4023
4024   AddCorrectionLevelSinglePt("secCorr");
4025
4026   TH1F* histPt = fCorrSinglePt[fNCorrectionLevelsSinglePt-2]->GetTrackPt(0);
4027
4028   if(histPt->GetNbinsX() != fh1SecCorrSinglePt->GetNbinsX())
4029     Printf("%s:%d: inconsistency pt spec and secondaries corr bins - rebin effCorr ...", (char*)__FILE__,__LINE__);
4030
4031   TString histNamePt = histPt->GetName();
4032   TH1F* hTrackPtSecCorr = (TH1F*) histPt->Clone(histNamePt);
4033
4034   hTrackPtSecCorr->Multiply(histPt,fh1SecCorrSinglePt,1,1,"");
4035
4036   fCorrSinglePt[fNCorrectionLevelsSinglePt-1]->AddCorrHistos(0,hTrackPtSecCorr);
4037 }