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
1419     THnSparse* hnResponse = 0;
1420     if(type == kFlagPt) hnResponse = fhnResponsePt[i];
1421     else if(type == kFlagZ)  hnResponse = fhnResponseZ[i];
1422     else if(type == kFlagXi) hnResponse = fhnResponseXi[i];
1423
1424
1425     TH1F* hPrior = 0;
1426     if(type == kFlagPt && fh1FFTrackPtPrior[i]  && ((TString(fh1FFTrackPtPrior[i]->GetName())).Length() > 0) ) hPrior = fh1FFTrackPtPrior[i];
1427     else if(type == kFlagZ  && fh1FFZPrior[i]   && ((TString(fh1FFZPrior[i]->GetName())).Length() > 0)       ) hPrior = fh1FFZPrior[i];
1428     else if(type == kFlagXi && fh1FFXiPrior[i]  && ((TString(fh1FFXiPrior[i]->GetName())).Length() > 0)      ) hPrior = fh1FFXiPrior[i];
1429
1430
1431     TString histNameTHn = hist->GetName();
1432     histNameTHn.ReplaceAll("TH1","THn");
1433
1434     TString priorNameTHn;
1435     if(hPrior){
1436       priorNameTHn = hPrior->GetName();
1437       priorNameTHn.ReplaceAll("TH1","THn");
1438     }
1439
1440     TString histNameBackFolded = hist->GetName();
1441     histNameBackFolded.Append("_backfold");
1442
1443     TString histNameRatioFolded = hist->GetName();
1444     histNameRatioFolded.ReplaceAll("fh1FF","hRatioFF");
1445     histNameRatioFolded.Append("_unfold");
1446
1447     TString histNameRatioBackFolded = hist->GetName();
1448     histNameRatioBackFolded.ReplaceAll("fh1FF","hRatioFF");
1449     histNameRatioBackFolded.Append("_backfold");
1450
1451     THnSparse* hnHist           = TH1toSparse(hist,histNameTHn,hist->GetTitle());
1452     THnSparse* hnFlatEfficiency = TH1toSparse(hist,"fhnEfficiency","eff",kTRUE); // could optionally also use real eff
1453     THnSparse* hnPrior          = 0;
1454     if(hPrior) hnPrior = TH1toSparse(hPrior,priorNameTHn,hPrior->GetTitle());
1455
1456     THnSparse* hnUnfolded
1457       = Unfold(hnHist,hnResponse,hnFlatEfficiency,nIter,useCorrelatedErrors,hnPrior);
1458
1459     TH1F* hUnfolded = (TH1F*) hnUnfolded->Projection(0);
1460     hUnfolded->SetNameTitle(hist->GetName(),hist->GetTitle());
1461
1462     if(type == kFlagPt) fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,hUnfolded,0,0);
1463     if(type == kFlagZ)  fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,0,hUnfolded,0);
1464     if(type == kFlagXi) fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,0,0,hUnfolded);
1465
1466     // backfolding: apply response matrix to unfolded spectrum
1467     TH1F* hBackFolded = ApplyResponse(hUnfolded,hnResponse);
1468     hBackFolded->SetNameTitle(histNameBackFolded,hist->GetTitle());
1469
1470     if(type == kFlagPt) fh1FFTrackPtBackFolded[i] = hBackFolded;
1471     if(type == kFlagZ)  fh1FFZBackFolded[i]       = hBackFolded;
1472     if(type == kFlagXi) fh1FFXiBackFolded[i]      = hBackFolded;
1473
1474     // ratio unfolded to original histo
1475     TH1F* hRatioUnfolded = (TH1F*) hUnfolded->Clone(histNameRatioFolded);
1476     hRatioUnfolded->Reset();
1477     hRatioUnfolded->Divide(hUnfolded,hist,1,1,"B");
1478
1479     if(type == kFlagPt) fh1FFRatioTrackPtFolded[i] = hRatioUnfolded;
1480     if(type == kFlagZ)  fh1FFRatioZFolded[i]       = hRatioUnfolded;
1481     if(type == kFlagXi) fh1FFRatioXiFolded[i]      = hRatioUnfolded;
1482
1483
1484     // ratio backfolded to original histo
1485     TH1F* hRatioBackFolded = (TH1F*) hBackFolded->Clone(histNameRatioBackFolded);
1486     hRatioBackFolded->Reset();
1487     hRatioBackFolded->Divide(hBackFolded,hist,1,1,"B");
1488
1489     if(type == kFlagPt) fh1FFRatioTrackPtBackFolded[i] = hRatioBackFolded;
1490     if(type == kFlagZ)  fh1FFRatioZBackFolded[i]       = hRatioBackFolded;
1491     if(type == kFlagXi) fh1FFRatioXiBackFolded[i]      = hRatioBackFolded;
1492
1493     delete hnHist;
1494     delete hnFlatEfficiency;
1495
1496  }
1497 }
1498
1499 //_____________________________________________________________________________________________________
1500 void AliFragmentationFunctionCorrections::UnfoldPt(const Int_t nIter, const Bool_t useCorrelatedErrors)
1501 {
1502
1503   Int_t type = kFlagPt;
1504   UnfoldHistos(nIter, useCorrelatedErrors, type);
1505 }
1506
1507 //_____________________________________________________________________________________________________
1508 void AliFragmentationFunctionCorrections::UnfoldZ(const Int_t nIter, const Bool_t useCorrelatedErrors)
1509 {
1510
1511   Int_t type = kFlagZ;
1512   UnfoldHistos(nIter, useCorrelatedErrors, type);
1513 }
1514
1515 //_____________________________________________________________________________________________________
1516 void AliFragmentationFunctionCorrections::UnfoldXi(const Int_t nIter, const Bool_t useCorrelatedErrors)
1517 {
1518
1519   Int_t type = kFlagXi;
1520   UnfoldHistos(nIter, useCorrelatedErrors, type);
1521 }
1522
1523 //______________________________________________________________________________________________
1524 TH1F* AliFragmentationFunctionCorrections::ApplyResponse(const TH1F* hist, THnSparse* hnResponse)
1525 {
1526   // apply (multiply) response matrix to hist
1527
1528   TH1F* hOut = new TH1F(*hist);
1529   hOut->Reset();
1530
1531   const Int_t axisM = 0;
1532   const Int_t axisT = 1;
1533
1534   Int_t nBinsM = hnResponse->GetAxis(axisM)->GetNbins();
1535   Int_t nBinsT = hnResponse->GetAxis(axisT)->GetNbins();
1536
1537   // response matrix normalization
1538   // do it in this function and not when reading response, since for 'proper' normalization errors are difficult to assign
1539   // so for unfolding proper we leave it to CORRFW ...
1540
1541   Double_t normFacResponse[nBinsT];
1542
1543   for(Int_t iT=1; iT<=nBinsT; iT++){
1544
1545     Double_t sumResp = 0;
1546
1547     for(Int_t iM=1; iM<=nBinsM; iM++){
1548
1549       Double_t coordM = hnResponse->GetAxis(axisM)->GetBinCenter(iM);
1550       Double_t coordT = hnResponse->GetAxis(axisT)->GetBinCenter(iT);
1551
1552       Double_t binCoord[] = {coordM,coordT};
1553
1554       Long64_t binIndex = hnResponse->GetBin(binCoord);
1555
1556       Double_t resp = hnResponse->GetBinContent(binIndex);
1557
1558       sumResp += resp;
1559     }
1560
1561     normFacResponse[iT] = 0;
1562     if(sumResp) normFacResponse[iT] = 1/sumResp;
1563   }
1564
1565
1566
1567   for(Int_t iM=1; iM<=nBinsM; iM++){
1568
1569     Double_t contM   = 0;
1570
1571     for(Int_t iT=1; iT<=nBinsT; iT++){
1572
1573       Double_t contT = hist->GetBinContent(iT);
1574
1575       Double_t coordM = hnResponse->GetAxis(axisM)->GetBinCenter(iM);
1576       Double_t coordT = hnResponse->GetAxis(axisT)->GetBinCenter(iT);
1577
1578       Double_t binCoord[] = {coordM,coordT};
1579
1580       Long64_t binIndex = hnResponse->GetBin(binCoord);
1581
1582       Double_t resp = hnResponse->GetBinContent(binIndex);
1583
1584       contM   += resp*normFacResponse[iT]*contT;
1585     }
1586
1587     hOut->SetBinContent(iM,contM);
1588   }
1589
1590   return hOut;
1591 }
1592
1593 //_______________________________________________________________________________________________________
1594 void AliFragmentationFunctionCorrections::ReadEfficiency(TString strfile, TString strdir, TString strlist)
1595 {
1596   // read reconstruction efficiency from file
1597   // argument strlist optional - read from directory strdir if not specified
1598
1599   // temporary histos to hold histos from file
1600   TH1F* hEffPt[fNJetPtSlices];
1601   TH1F* hEffZ[fNJetPtSlices];
1602   TH1F* hEffXi[fNJetPtSlices];
1603
1604   for(Int_t i=0; i<fNJetPtSlices; i++) hEffPt[i] = 0;
1605   for(Int_t i=0; i<fNJetPtSlices; i++) hEffZ[i]  = 0;
1606   for(Int_t i=0; i<fNJetPtSlices; i++) hEffXi[i] = 0;
1607
1608   TFile f(strfile,"READ");
1609
1610   if(!f.IsOpen()){
1611     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
1612     return;
1613   }
1614
1615   if(fDebug>0) Printf("%s:%d -- read efficiencies from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
1616
1617   if(strdir && strdir.Length()) gDirectory->cd(strdir);
1618
1619   TList* list = 0;
1620
1621   if(strlist && strlist.Length()){
1622
1623     if(!(list = (TList*) gDirectory->Get(strlist))){
1624       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
1625       return;
1626     }
1627   }
1628
1629   for(Int_t i=0; i<fNJetPtSlices; i++){
1630
1631     Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
1632     Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
1633
1634     TString strNameEffPt(Form("hEffPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
1635     TString strNameEffZ(Form("hEffZ_%02d_%02d",jetPtLoLim,jetPtUpLim));
1636     TString strNameEffXi(Form("hEffXi_%02d_%02d",jetPtLoLim,jetPtUpLim));
1637
1638
1639     if(list){
1640       hEffPt[i] = (TH1F*) list->FindObject(strNameEffPt);
1641       hEffZ[i]  = (TH1F*) list->FindObject(strNameEffZ);
1642       hEffXi[i] = (TH1F*) list->FindObject(strNameEffXi);
1643     }
1644     else{
1645       hEffPt[i] = (TH1F*) gDirectory->Get(strNameEffPt);
1646       hEffZ[i]  = (TH1F*) gDirectory->Get(strNameEffZ);
1647       hEffXi[i] = (TH1F*) gDirectory->Get(strNameEffXi);
1648     }
1649
1650     if(!hEffPt[i]){
1651       Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffPt.Data());
1652     }
1653
1654     if(!hEffZ[i]){
1655       Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffZ.Data());
1656     }
1657
1658     if(!hEffXi[i]){
1659       Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffXi.Data());
1660     }
1661
1662
1663     if(fNHistoBinsPt[i]) hEffPt[i] = (TH1F*) hEffPt[i]->Rebin(fNHistoBinsPt[i],strNameEffPt+"_rebin",fHistoBinsPt[i]->GetArray());
1664     if(fNHistoBinsZ[i])  hEffZ[i]  = (TH1F*) hEffZ[i]->Rebin(fNHistoBinsZ[i],strNameEffZ+"_rebin",fHistoBinsZ[i]->GetArray());
1665     if(fNHistoBinsXi[i]) hEffXi[i] = (TH1F*) hEffXi[i]->Rebin(fNHistoBinsXi[i],strNameEffXi+"_rebin",fHistoBinsXi[i]->GetArray());
1666
1667     if(hEffPt[i]) hEffPt[i]->SetDirectory(0);
1668     if(hEffZ[i])  hEffZ[i]->SetDirectory(0);
1669     if(hEffXi[i]) hEffXi[i]->SetDirectory(0);
1670
1671   } // jet slices loop
1672
1673   f.Close();
1674
1675   for(Int_t i=0; i<fNJetPtSlices; i++){ // 2nd loop: need to close input file before placing histos
1676     if(hEffPt[i]) new(fh1EffPt[i]) TH1F(*hEffPt[i]);
1677     if(hEffZ[i])  new(fh1EffZ[i])  TH1F(*hEffZ[i]);
1678     if(hEffXi[i]) new(fh1EffXi[i]) TH1F(*hEffXi[i]);
1679   }
1680 }
1681
1682 //___________________________________________________________________________________________________________
1683 void AliFragmentationFunctionCorrections::ReadBgrEfficiency(TString strfile, TString strdir, TString strlist)
1684 {
1685   // read bgr eff from file
1686   // argument strlist optional - read from directory strdir if not specified
1687
1688   TH1F* hEffPtBgr[fNJetPtSlices];
1689   TH1F* hEffZBgr [fNJetPtSlices];
1690   TH1F* hEffXiBgr[fNJetPtSlices];
1691
1692   for(Int_t i=0; i<fNJetPtSlices; i++) hEffPtBgr[i] = 0;
1693   for(Int_t i=0; i<fNJetPtSlices; i++) hEffZBgr[i]  = 0;
1694   for(Int_t i=0; i<fNJetPtSlices; i++) hEffXiBgr[i] = 0;
1695
1696
1697   TFile f(strfile,"READ");
1698
1699   if(!f.IsOpen()){
1700     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
1701     return;
1702   }
1703
1704   if(fDebug>0) Printf("%s:%d -- read bgr efficiencies from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
1705
1706   if(strdir && strdir.Length()) gDirectory->cd(strdir);
1707
1708   TList* list = 0;
1709
1710   if(strlist && strlist.Length()){
1711
1712     if(!(list = (TList*) gDirectory->Get(strlist))){
1713       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
1714       return;
1715     }
1716   }
1717
1718   for(Int_t i=0; i<fNJetPtSlices; i++){
1719
1720     Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
1721     Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
1722
1723     TString strNameEffPtBgr(Form("hEffPtBgr%02dto%02d",jetPtLoLim,jetPtUpLim));
1724     TString strNameEffZBgr(Form("hEffZBgr%02dto%02d",jetPtLoLim,jetPtUpLim));
1725     TString strNameEffXiBgr(Form("hEffXiBgr%02dto%02d",jetPtLoLim,jetPtUpLim));
1726
1727
1728     if(list){
1729       hEffPtBgr[i] = (TH1F*) list->FindObject(strNameEffPtBgr);
1730       hEffZBgr[i]  = (TH1F*) list->FindObject(strNameEffZBgr);
1731       hEffXiBgr[i] = (TH1F*) list->FindObject(strNameEffXiBgr);
1732     }
1733     else{
1734       hEffPtBgr[i] = (TH1F*) gDirectory->Get(strNameEffPtBgr);
1735       hEffZBgr[i]  = (TH1F*) gDirectory->Get(strNameEffZBgr);
1736       hEffXiBgr[i] = (TH1F*) gDirectory->Get(strNameEffXiBgr);
1737     }
1738
1739     if(!hEffPtBgr[i]){
1740       Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffPtBgr.Data());
1741     }
1742
1743     if(!hEffZBgr[i]){
1744       Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffZBgr.Data());
1745     }
1746
1747     if(!hEffXiBgr[i]){
1748       Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffXiBgr.Data());
1749     }
1750
1751
1752     if(fNHistoBinsPt[i]) hEffPtBgr[i] = (TH1F*) hEffPtBgr[i]->Rebin(fNHistoBinsPt[i],strNameEffPtBgr+"_rebin",fHistoBinsPt[i]->GetArray());
1753     if(fNHistoBinsZ[i])  hEffZBgr[i]  = (TH1F*) hEffZBgr[i]->Rebin(fNHistoBinsZ[i],strNameEffZBgr+"_rebin",fHistoBinsZ[i]->GetArray());
1754     if(fNHistoBinsXi[i]) hEffXiBgr[i] = (TH1F*) hEffXiBgr[i]->Rebin(fNHistoBinsXi[i],strNameEffXiBgr+"_rebin",fHistoBinsXi[i]->GetArray());
1755
1756     if(hEffPtBgr[i]) hEffPtBgr[i]->SetDirectory(0);
1757     if(hEffZBgr[i])  hEffZBgr[i]->SetDirectory(0);
1758     if(hEffXiBgr[i]) hEffXiBgr[i]->SetDirectory(0);
1759
1760   } // jet slices loop
1761
1762   f.Close();
1763
1764   for(Int_t i=0; i<fNJetPtSlices; i++){ // 2nd loop: need to close input file before placing histos
1765     if(hEffPtBgr[i]) new(fh1EffBgrPt[i]) TH1F(*hEffPtBgr[i]);
1766     if(hEffZBgr[i])  new(fh1EffBgrZ[i])  TH1F(*hEffZBgr[i]);
1767     if(hEffXiBgr[i]) new(fh1EffBgrXi[i]) TH1F(*hEffXiBgr[i]);
1768   }
1769 }
1770
1771 // ________________________________________________
1772 void AliFragmentationFunctionCorrections::EffCorr()
1773 {
1774   // apply efficiency correction
1775
1776   AddCorrectionLevel("eff");
1777
1778   for(Int_t i=0; i<fNJetPtSlices; i++){
1779
1780     TH1F* histPt = fCorrFF[fNCorrectionLevels-2]->GetTrackPt(i);
1781     TH1F* histZ  = fCorrFF[fNCorrectionLevels-2]->GetZ(i);
1782     TH1F* histXi = fCorrFF[fNCorrectionLevels-2]->GetXi(i);
1783
1784     TString histNamePt = histPt->GetName();
1785     TString histNameZ  = histZ->GetName();
1786     TString histNameXi = histXi->GetName();
1787
1788
1789     TH1F* hFFTrackPtEffCorr = (TH1F*) histPt->Clone(histNamePt);
1790     hFFTrackPtEffCorr->Divide(histPt,fh1EffPt[i],1,1,"");
1791
1792     TH1F* hFFZEffCorr = (TH1F*) histZ->Clone(histNameZ);
1793     hFFZEffCorr->Divide(histZ,fh1EffZ[i],1,1,"");
1794
1795     TH1F* hFFXiEffCorr = (TH1F*) histXi->Clone(histNameXi);
1796     hFFXiEffCorr->Divide(histXi,fh1EffXi[i],1,1,"");
1797
1798     fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,hFFTrackPtEffCorr,hFFZEffCorr,hFFXiEffCorr);
1799   }
1800 }
1801
1802 //___________________________________________________
1803 void AliFragmentationFunctionCorrections::EffCorrBgr()
1804 {
1805   // apply efficiency correction to bgr distributions
1806
1807   AddCorrectionLevelBgr("eff");
1808
1809   Printf("%s:%d -- apply efficiency correction, corrLevel %d",(char*)__FILE__,__LINE__,fNCorrectionLevels-1);
1810
1811   for(Int_t i=0; i<fNJetPtSlices; i++){
1812
1813     TH1F* histPt = fCorrBgr[fNCorrectionLevelsBgr-2]->GetTrackPt(i);
1814     TH1F* histZ  = fCorrBgr[fNCorrectionLevelsBgr-2]->GetZ(i);
1815     TH1F* histXi = fCorrBgr[fNCorrectionLevelsBgr-2]->GetXi(i);
1816
1817     TString histNamePt = histPt->GetName();
1818     TString histNameZ  = histZ->GetName();
1819     TString histNameXi = histXi->GetName();
1820
1821
1822     TH1F* hFFTrackPtEffCorr = (TH1F*) histPt->Clone(histNamePt);
1823     hFFTrackPtEffCorr->Divide(histPt,fh1EffPt[i],1,1,"");
1824
1825     TH1F* hFFZEffCorr = (TH1F*) histZ->Clone(histNameZ);
1826     hFFZEffCorr->Divide(histZ,fh1EffZ[i],1,1,"");
1827
1828     TH1F* hFFXiEffCorr = (TH1F*) histXi->Clone(histNameXi);
1829     hFFXiEffCorr->Divide(histXi,fh1EffXi[i],1,1,"");
1830
1831     fCorrBgr[fNCorrectionLevelsBgr-1]->AddCorrHistos(i,hFFTrackPtEffCorr,hFFZEffCorr,hFFXiEffCorr);
1832   }
1833 }
1834
1835 //______________________________________________________________________
1836 void AliFragmentationFunctionCorrections::XiShift(const Int_t corrLevel)
1837 {
1838   // re-evaluate jet energy after FF corrections from dN/dpt distribution
1839   // 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'
1840   // argument corrlevel: which level of correction to be corrected/shifted to
1841
1842   if(corrLevel>=fNCorrectionLevels){
1843     Printf(" calc xi shift: corrLevel exceeded - do nothing");
1844     return;
1845   }
1846
1847   Double_t* jetPtUncorr = new Double_t[fNJetPtSlices];
1848   Double_t* jetPtCorr   = new Double_t[fNJetPtSlices];
1849   Double_t* deltaXi     = new Double_t[fNJetPtSlices];
1850
1851   for(Int_t i=0; i<fNJetPtSlices; i++){
1852
1853     TH1F* histPtRaw = fCorrFF[0]->GetTrackPt(i);
1854     TH1F* histPt    = fCorrFF[corrLevel]->GetTrackPt(i);
1855
1856     Double_t ptUncorr = 0;
1857     Double_t ptCorr   = 0;
1858
1859     for(Int_t bin = 1; bin<=histPtRaw->GetNbinsX(); bin++){
1860
1861       Double_t cont   = histPtRaw->GetBinContent(bin);
1862       Double_t width  = histPtRaw->GetBinWidth(bin);
1863       Double_t meanPt = histPtRaw->GetBinCenter(bin);
1864
1865       ptUncorr += meanPt*cont*width;
1866     }
1867
1868     for(Int_t bin = 1; bin<=histPt->GetNbinsX(); bin++){
1869
1870       Double_t cont   = histPt->GetBinContent(bin);
1871       Double_t width  = histPt->GetBinWidth(bin);
1872       Double_t meanPt = histPt->GetBinCenter(bin);
1873
1874       ptCorr += meanPt*cont*width;
1875     }
1876
1877     jetPtUncorr[i] = ptUncorr;
1878     jetPtCorr[i]   = ptCorr;
1879   }
1880
1881   // calc dXi from dN/dpt distribution :
1882   // sum over track pt for raw and corrected FF is equivalent to raw/corrected jet pt
1883
1884   for(Int_t i=0; i<fNJetPtSlices; i++){
1885
1886     Float_t jetPtLoLim = fJetPtSlices->At(i);
1887     Float_t jetPtUpLim = fJetPtSlices->At(i+1);
1888
1889     Double_t meanJetPt = 0.5*(jetPtUpLim+jetPtLoLim);
1890
1891     Double_t ptUncorr = jetPtUncorr[i];
1892     Double_t ptCorr   = jetPtCorr[i];
1893
1894     Double_t dXi = TMath::Log(ptCorr/ptUncorr);
1895
1896     Printf(" calc xi shift: jet pt slice %d, mean jet pt %f, ptUncorr %f, ptCorr %f, ratio corr/uncorr %f, dXi %f "
1897            ,i,meanJetPt,ptUncorr,ptCorr,ptCorr/ptUncorr,dXi);
1898
1899     deltaXi[i] = dXi;
1900   }
1901
1902   // book & fill new dN/dxi histos
1903
1904   for(Int_t i=0; i<fNJetPtSlices; i++){
1905
1906     TH1F* histXi = fCorrFF[corrLevel]->GetXi(i);
1907
1908     Double_t dXi = deltaXi[i];
1909
1910     Int_t nBins  = histXi->GetNbinsX();
1911     const Double_t* binsVec = histXi->GetXaxis()->GetXbins()->GetArray();
1912     Float_t binsVecNew[nBins+1];
1913
1914     TString strName = histXi->GetName();
1915     strName.Append("_shift");
1916     TString strTit  = histXi->GetTitle();
1917
1918     TH1F* hXiCorr;
1919
1920     // create shifted histo ...
1921
1922     if(binsVec){ // binsVec only neq NULL if histo was rebinned before
1923
1924       for(Int_t bin=0; bin<nBins+1; bin++) binsVecNew[bin] = binsVec[bin] + dXi;
1925       hXiCorr = new TH1F(strName,strTit,nBins,binsVecNew);
1926     }
1927     else{ // uniform bin size
1928
1929       Double_t xMin  = histXi->GetXaxis()->GetXmin();
1930       Double_t xMax  = histXi->GetXaxis()->GetXmax();
1931
1932       xMin += dXi;
1933       xMax += dXi;
1934
1935       hXiCorr = new TH1F(strName,strTit,nBins,xMin,xMax);
1936     }
1937
1938     // ... and fill
1939
1940     for(Int_t bin=1; bin<nBins+1; bin++){
1941       Double_t cont = histXi->GetBinContent(bin);
1942       Double_t err  = histXi->GetBinError(bin);
1943
1944       hXiCorr->SetBinContent(bin,cont);
1945       hXiCorr->SetBinError(bin,err);
1946     }
1947
1948     new(fh1FFXiShift[i]) TH1F(*hXiCorr);
1949     delete hXiCorr;
1950   }
1951
1952   delete[] jetPtUncorr;
1953   delete[] jetPtCorr;
1954   delete[] deltaXi;
1955 }
1956
1957 //_____________________________________________________
1958 void AliFragmentationFunctionCorrections::SubtractBgr()
1959 {
1960   // subtract bgr distribution from FF
1961   // the current corr level is used for both
1962
1963   AddCorrectionLevel("bgrSub");
1964
1965   for(Int_t i=0; i<fNJetPtSlices; i++){
1966
1967     TH1F* histPt = fCorrFF[fNCorrectionLevels-2]->GetTrackPt(i);
1968     TH1F* histZ  = fCorrFF[fNCorrectionLevels-2]->GetZ(i);
1969     TH1F* histXi = fCorrFF[fNCorrectionLevels-2]->GetXi(i);
1970
1971     TH1F* histPtBgr = fCorrBgr[fNCorrectionLevelsBgr-1]->GetTrackPt(i);
1972     TH1F* histZBgr  = fCorrBgr[fNCorrectionLevelsBgr-1]->GetZ(i);
1973     TH1F* histXiBgr = fCorrBgr[fNCorrectionLevelsBgr-1]->GetXi(i);
1974
1975     TString histNamePt = histPt->GetName();
1976     TString histNameZ  = histZ->GetName();
1977     TString histNameXi = histXi->GetName();
1978
1979
1980     TH1F* hFFTrackPtBgrSub = (TH1F*) histPt->Clone(histNamePt);
1981     hFFTrackPtBgrSub->Add(histPtBgr,-1);
1982
1983     TH1F* hFFZBgrSub =  (TH1F*) histZ->Clone(histNameZ);
1984     hFFZBgrSub->Add(histZBgr,-1);
1985
1986     TH1F* hFFXiBgrSub = (TH1F*) histXi->Clone(histNameXi);
1987     hFFXiBgrSub->Add(histXiBgr,-1);
1988
1989     fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,hFFTrackPtBgrSub,hFFZBgrSub,hFFXiBgrSub);
1990   }
1991 }
1992
1993 //________________________________________________________________________________________________________________
1994 void AliFragmentationFunctionCorrections::WriteSingleTrackEff(TString strInfile, TString strID, TString strOutfile,
1995                                                               Bool_t updateOutfile, TString strOutDir,TString strPostfix)
1996 {
1997   // read task ouput from MC and write single track eff - standard dir/list
1998
1999   TString strdir  = "PWG4_FragmentationFunction_" + strID;
2000   TString strlist = "fracfunc_" + strID;
2001
2002   WriteSingleTrackEff(strInfile,strdir,strlist,strOutfile,updateOutfile,strOutDir,strPostfix);
2003 }
2004
2005 //___________________________________________________________________________________________________________________________________
2006 void AliFragmentationFunctionCorrections::WriteSingleTrackEff(TString strInfile, TString strdir, TString strlist,
2007                                                               TString strOutfile, Bool_t updateOutfile, TString strOutDir,TString strPostfix)
2008 {
2009   // read task output from MC and write single track eff as function of pt, eta, phi
2010   // argument strlist optional - read from directory strdir if not specified
2011   // write eff to file stroutfile - by default only 'update' file (don't overwrite)
2012
2013
2014   TH1D* hdNdptTracksMCPrimGen;
2015   TH2D* hdNdetadphiTracksMCPrimGen;
2016
2017   TH1D* hdNdptTracksMCPrimRec;
2018   TH2D* hdNdetadphiTracksMCPrimRec;
2019
2020
2021   TFile f(strInfile,"READ");
2022
2023   if(!f.IsOpen()){
2024     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2025     return;
2026   }
2027
2028   if(fDebug>0) Printf("%s:%d -- writeSingleTrackEff: open task ouput file %s ",(char*)__FILE__,__LINE__,strInfile.Data());
2029
2030   if(strdir && strdir.Length()){
2031     if(fDebug>0) Printf("%s:%d -- writeSingleTrackEff: change dir to %s",(char*)__FILE__,__LINE__,strdir.Data());
2032     gDirectory->cd(strdir);
2033   }
2034
2035   TList* list = 0;
2036
2037   if(strlist && strlist.Length()){
2038
2039     if(!(list = (TList*) gDirectory->Get(strlist))){
2040       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2041       return;
2042     }
2043   }
2044
2045
2046   TString hnamePtRecEffGen = "fh1TrackQAPtRecEffGen";
2047   if(strPostfix.Length()) hnamePtRecEffGen.Form("fh1TrackQAPtRecEffGen_%s",strPostfix.Data());
2048
2049   TString hnamePtRecEffRec = "fh1TrackQAPtRecEffRec";
2050   if(strPostfix.Length()) hnamePtRecEffRec.Form("fh1TrackQAPtRecEffRec_%s",strPostfix.Data());
2051
2052   TString hnameEtaPhiRecEffGen = "fh2TrackQAEtaPhiRecEffGen";
2053   if(strPostfix.Length()) hnameEtaPhiRecEffGen.Form("fh2TrackQAEtaPhiRecEffGen_%s",strPostfix.Data());
2054
2055   TString hnameEtaPhiRecEffRec = "fh2TrackQAEtaPhiRecEffRec";
2056   if(strPostfix.Length()) hnameEtaPhiRecEffRec.Form("fh2TrackQAEtaPhiRecEffRec_%s",strPostfix.Data());
2057
2058
2059   if(list){
2060     hdNdptTracksMCPrimGen       = (TH1D*) list->FindObject(hnamePtRecEffGen);
2061     hdNdetadphiTracksMCPrimGen  = (TH2D*) list->FindObject(hnameEtaPhiRecEffGen);
2062
2063     hdNdptTracksMCPrimRec       = (TH1D*) list->FindObject(hnamePtRecEffRec);
2064     hdNdetadphiTracksMCPrimRec  = (TH2D*) list->FindObject(hnameEtaPhiRecEffRec);
2065   }
2066   else{
2067     hdNdptTracksMCPrimGen       = (TH1D*) gDirectory->Get(hnamePtRecEffGen);
2068     hdNdetadphiTracksMCPrimGen  = (TH2D*) gDirectory->Get(hnameEtaPhiRecEffGen);
2069
2070     hdNdptTracksMCPrimRec       = (TH1D*) gDirectory->Get(hnamePtRecEffRec);
2071     hdNdetadphiTracksMCPrimRec  = (TH2D*) gDirectory->Get(hnameEtaPhiRecEffRec);
2072   }
2073
2074   hdNdptTracksMCPrimGen->SetDirectory(0);
2075   hdNdetadphiTracksMCPrimGen->SetDirectory(0);
2076   hdNdptTracksMCPrimRec->SetDirectory(0);
2077   hdNdetadphiTracksMCPrimRec->SetDirectory(0);
2078
2079   f.Close();
2080
2081   // projections: dN/deta, dN/dphi
2082
2083   TH1D* hdNdetaTracksMCPrimGen = (TH1D*) hdNdetadphiTracksMCPrimGen->ProjectionX("hdNdetaTracksMcPrimGen");
2084   TH1D* hdNdphiTracksMCPrimGen = (TH1D*) hdNdetadphiTracksMCPrimGen->ProjectionY("hdNdphiTracksMcPrimGen");
2085
2086   TH1D* hdNdetaTracksMCPrimRec = (TH1D*) hdNdetadphiTracksMCPrimRec->ProjectionX("hdNdetaTracksMcPrimRec");
2087   TH1D* hdNdphiTracksMCPrimRec = (TH1D*) hdNdetadphiTracksMCPrimRec->ProjectionY("hdNdphiTracksMcPrimRec");
2088
2089   // rebin
2090
2091   TString strNamePtGen = "hTrackPtGenPrim";
2092   TString strNamePtRec = "hTrackPtRecPrim";
2093
2094   if(fNHistoBinsSinglePt) hdNdptTracksMCPrimGen = (TH1D*) hdNdptTracksMCPrimGen->Rebin(fNHistoBinsSinglePt,strNamePtGen,fHistoBinsSinglePt->GetArray());
2095   if(fNHistoBinsSinglePt) hdNdptTracksMCPrimRec = (TH1D*) hdNdptTracksMCPrimRec->Rebin(fNHistoBinsSinglePt,strNamePtRec,fHistoBinsSinglePt->GetArray());
2096
2097     // efficiency: divide
2098
2099   TString hNameTrackEffPt = "hSingleTrackEffPt";
2100   if(strPostfix.Length()) hNameTrackEffPt.Form("hSingleTrackEffPt_%s",strPostfix.Data());
2101
2102   TString hNameTrackEffEta = "hSingleTrackEffEta";
2103   if(strPostfix.Length()) hNameTrackEffEta.Form("hSingleTrackEffEta_%s",strPostfix.Data());
2104
2105   TString hNameTrackEffPhi = "hSingleTrackEffPhi";
2106   if(strPostfix.Length()) hNameTrackEffPhi.Form("hSingleTrackEffPhi_%s",strPostfix.Data());
2107
2108
2109   TH1F* hSingleTrackEffPt = (TH1F*) hdNdptTracksMCPrimRec->Clone(hNameTrackEffPt);
2110   hSingleTrackEffPt->Divide(hdNdptTracksMCPrimRec,hdNdptTracksMCPrimGen,1,1,"B"); // binominal errors
2111
2112   TH1F* hSingleTrackEffEta = (TH1F*) hdNdetaTracksMCPrimRec->Clone(hNameTrackEffEta);
2113   hSingleTrackEffEta->Divide(hdNdetaTracksMCPrimRec,hdNdetaTracksMCPrimGen,1,1,"B"); // binominal errors
2114
2115   TH1F* hSingleTrackEffPhi = (TH1F*) hdNdphiTracksMCPrimRec->Clone(hNameTrackEffPhi);
2116   hSingleTrackEffPhi->Divide(hdNdphiTracksMCPrimRec,hdNdphiTracksMCPrimGen,1,1,"B"); // binominal errors
2117
2118
2119   TString outfileOption = "RECREATE";
2120   if(updateOutfile)  outfileOption = "UPDATE";
2121
2122   TFile out(strOutfile,outfileOption);
2123
2124   if(!out.IsOpen()){
2125     Printf("%s:%d -- error opening efficiency output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
2126     return;
2127   }
2128
2129   if(fDebug>0) Printf("%s:%d -- write efficiency to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
2130
2131   if(strOutDir && strOutDir.Length()){
2132
2133     TDirectory* dir;
2134     if((dir = ((TDirectory*) gDirectory->Get(strOutDir)))) dir->cd();
2135     else{
2136       dir = out.mkdir(strOutDir);
2137       dir->cd();
2138     }
2139   }
2140
2141   hSingleTrackEffPt->Write();
2142   hSingleTrackEffEta->Write();
2143   hSingleTrackEffPhi->Write();
2144
2145   out.Close();
2146 }
2147
2148 //________________________________________________________________________________________________________________
2149 void AliFragmentationFunctionCorrections::WriteSingleTrackSecCorr(TString strInfile, TString strID, TString strOutfile,
2150                                                                   Bool_t updateOutfile, TString strOutDir)
2151 {
2152   // read task ouput from MC and write single track eff - standard dir/list
2153
2154   TString strdir  = "PWG4_FragmentationFunction_" + strID;
2155   TString strlist = "fracfunc_" + strID;
2156
2157   WriteSingleTrackSecCorr(strInfile,strdir,strlist,strOutfile,updateOutfile,strOutDir);
2158 }
2159
2160 //___________________________________________________________________________________________________________________________________
2161 void AliFragmentationFunctionCorrections::WriteSingleTrackSecCorr(TString strInfile, TString strdir, TString strlist,
2162                                                                   TString strOutfile, Bool_t updateOutfile, TString strOutDir)
2163 {
2164   // read task output from MC and write single track secondaries contamination correction as function of pt, eta, phi
2165   // argument strlist optional - read from directory strdir if not specified
2166   // write corr factor to file stroutfile - by default only 'update' file (don't overwrite)
2167
2168   TH1D* hdNdptTracksMCPrimRec;
2169   TH2D* hdNdetadphiTracksMCPrimRec;
2170
2171   TH1D* hdNdptTracksMCSecRec;
2172   TH2D* hdNdetadphiTracksMCSecRec;
2173
2174   TFile f(strInfile,"READ");
2175
2176   if(!f.IsOpen()){
2177     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2178     return;
2179   }
2180
2181   if(fDebug>0) Printf("%s:%d -- writeSingleTrackEff: open task ouput file %s",(char*)__FILE__,__LINE__,strInfile.Data());
2182
2183   if(strdir && strdir.Length()) gDirectory->cd(strdir);
2184
2185   TList* list = 0;
2186
2187   if(strlist && strlist.Length()){
2188
2189     if(!(list = (TList*) gDirectory->Get(strlist))){
2190       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2191       return;
2192     }
2193   }
2194
2195
2196   if(list){
2197     hdNdptTracksMCPrimRec       = (TH1D*) list->FindObject("fh1TrackQAPtRecEffGen");
2198     hdNdetadphiTracksMCPrimRec  = (TH2D*) list->FindObject("fh2TrackQAEtaPhiRecEffGen");
2199
2200     hdNdptTracksMCSecRec       = (TH1D*) list->FindObject("fh1TrackQAPtSecRec");
2201     hdNdetadphiTracksMCSecRec  = (TH2D*) list->FindObject("fh2TrackQAEtaPhiSecRec");
2202   }
2203   else{
2204     hdNdptTracksMCPrimRec       = (TH1D*) gDirectory->Get("fh1TrackQAPtRecEffGen");
2205     hdNdetadphiTracksMCPrimRec  = (TH2D*) gDirectory->Get("fh2TrackQAEtaPhiRecEffGen");
2206
2207     hdNdptTracksMCSecRec       = (TH1D*) gDirectory->Get("fh1TrackQAPtSecRec");
2208     hdNdetadphiTracksMCSecRec  = (TH2D*) gDirectory->Get("fh2TrackQAEtaPhiSecRec");
2209   }
2210
2211   hdNdptTracksMCPrimRec->SetDirectory(0);
2212   hdNdetadphiTracksMCPrimRec->SetDirectory(0);
2213   hdNdptTracksMCSecRec->SetDirectory(0);
2214   hdNdetadphiTracksMCSecRec->SetDirectory(0);
2215
2216   f.Close();
2217
2218   // projections: dN/deta, dN/dphi
2219
2220   TH1D* hdNdetaTracksMCPrimRec = (TH1D*) hdNdetadphiTracksMCPrimRec->ProjectionX("hdNdetaTracksMcPrimRec");
2221   TH1D* hdNdphiTracksMCPrimRec = (TH1D*) hdNdetadphiTracksMCPrimRec->ProjectionY("hdNdphiTracksMcPrimRec");
2222
2223   TH1D* hdNdetaTracksMCSecRec = (TH1D*) hdNdetadphiTracksMCSecRec->ProjectionX("hdNdetaTracksMcSecRec");
2224   TH1D* hdNdphiTracksMCSecRec = (TH1D*) hdNdetadphiTracksMCSecRec->ProjectionY("hdNdphiTracksMcSecRec");
2225
2226
2227   // rebin
2228
2229   TString strNamePtPrim = "hTrackPtPrim";
2230   TString strNamePtSec  = "hTrackPtSec";
2231
2232   if(fNHistoBinsSinglePt) hdNdptTracksMCPrimRec = (TH1D*) hdNdptTracksMCPrimRec->Rebin(fNHistoBinsSinglePt,strNamePtPrim,fHistoBinsSinglePt->GetArray());
2233   if(fNHistoBinsSinglePt) hdNdptTracksMCSecRec  = (TH1D*) hdNdptTracksMCSecRec->Rebin(fNHistoBinsSinglePt,strNamePtSec,fHistoBinsSinglePt->GetArray());
2234
2235
2236   // secondary correction factor: divide prim/(prim+sec)
2237
2238   TH1F* hSingleTrackSecCorrPt = (TH1F*) hdNdptTracksMCSecRec->Clone("hSingleTrackSecCorrPt");
2239   TH1F* hSumPrimSecPt = (TH1F*) hdNdptTracksMCSecRec->Clone("hSumPrimSecPt");
2240   hSumPrimSecPt->Add(hdNdptTracksMCPrimRec);
2241   hSingleTrackSecCorrPt->Divide(hdNdptTracksMCPrimRec,hSumPrimSecPt,1,1,"B"); // binominal errors
2242
2243   TH1F* hSingleTrackSecCorrEta = (TH1F*) hdNdetaTracksMCSecRec->Clone("hSingleTrackSecCorrEta");
2244   TH1F* hSumPrimSecEta = (TH1F*) hdNdetaTracksMCSecRec->Clone("hSumPrimSecEta");
2245   hSumPrimSecEta->Add(hdNdetaTracksMCPrimRec);
2246   hSingleTrackSecCorrEta->Divide(hdNdetaTracksMCPrimRec,hSumPrimSecEta,1,1,"B"); // binominal errors
2247
2248   TH1F* hSingleTrackSecCorrPhi = (TH1F*) hdNdphiTracksMCSecRec->Clone("hSingleTrackSecCorrPhi");
2249   TH1F* hSumPrimSecPhi = (TH1F*) hdNdphiTracksMCSecRec->Clone("hSumPrimSecPhi");
2250   hSumPrimSecPhi->Add(hdNdphiTracksMCPrimRec);
2251   hSingleTrackSecCorrPhi->Divide(hdNdphiTracksMCPrimRec,hSumPrimSecPhi,1,1,"B"); // binominal errors
2252
2253
2254   TString outfileOption = "RECREATE";
2255   if(updateOutfile)  outfileOption = "UPDATE";
2256
2257   TFile out(strOutfile,outfileOption);
2258
2259   if(!out.IsOpen()){
2260     Printf("%s:%d -- error opening secCorr output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
2261     return;
2262   }
2263
2264   if(fDebug>0) Printf("%s:%d -- write secCorr to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
2265
2266   if(strOutDir && strOutDir.Length()){
2267
2268     TDirectory* dir;
2269     if((dir = ((TDirectory*) gDirectory->Get(strOutDir)))) dir->cd();
2270     else{
2271       dir = out.mkdir(strOutDir);
2272       dir->cd();
2273     }
2274   }
2275
2276   hdNdptTracksMCSecRec->Write();
2277   hdNdetaTracksMCSecRec->Write();
2278   hdNdphiTracksMCSecRec->Write();
2279
2280   hSingleTrackSecCorrPt->Write();
2281   hSingleTrackSecCorrEta->Write();
2282   hSingleTrackSecCorrPhi->Write();
2283
2284   out.Close();
2285 }
2286
2287 //________________________________________________________________________________________________________________
2288 void AliFragmentationFunctionCorrections::WriteSingleResponse(TString strInfile, TString strID, TString strOutfile,
2289                                                               Bool_t updateOutfile, TString strOutDir)
2290 {
2291   // read task ouput from MC and write single track eff - standard dir/list
2292
2293   TString strdir  = "PWG4_FragmentationFunction_" + strID;
2294   TString strlist = "fracfunc_" + strID;
2295
2296   WriteSingleResponse(strInfile,strdir,strlist,strOutfile,updateOutfile,strOutDir);
2297 }
2298
2299
2300 //_____________________________________________________________________________________________________________________________________
2301 void AliFragmentationFunctionCorrections::WriteSingleResponse(TString strInfile, TString strdir, TString strlist,
2302                                                               TString strOutfile, Bool_t updateOutfile, TString strOutDir)
2303 {
2304   // read 2d THnSparse response matrices in pt from file
2305   // project TH2
2306   // write to strOutfile
2307
2308   THnSparse* hnResponseSinglePt;
2309   TH2F*      h2ResponseSinglePt;
2310
2311   TFile f(strInfile,"READ");
2312
2313   if(!f.IsOpen()){
2314     Printf("%s:%d -- error opening file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2315     return;
2316   }
2317
2318   if(fDebug>0) Printf("%s:%d -- read response matrices from file %s ",(char*)__FILE__,__LINE__,strInfile.Data());
2319
2320   if(strdir && strdir.Length()) gDirectory->cd(strdir);
2321
2322   TList* list = 0;
2323
2324   if(strlist && strlist.Length()){
2325
2326     if(!(list = (TList*) gDirectory->Get(strlist))){
2327       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2328       return;
2329     }
2330   }
2331
2332   if(list) hnResponseSinglePt = (THnSparse*) list->FindObject("fhnResponseSinglePt");
2333   else     hnResponseSinglePt = (THnSparse*) gDirectory->Get("fhnResponseSinglePt");
2334
2335
2336   if(!hnResponseSinglePt){
2337     Printf("%s:%d -- error retrieving response matrix fhnResponseSinglePt",(char*)__FILE__,__LINE__);
2338     return;
2339   }
2340
2341   f.Close();
2342
2343
2344   // project 2d histo
2345   h2ResponseSinglePt = (TH2F*) hnResponseSinglePt->Projection(1,0);// note convention: yDim,xDim
2346   h2ResponseSinglePt->SetNameTitle("h2ResponseSinglePt","");
2347
2348
2349   // write
2350
2351   TString outfileOption = "RECREATE";
2352   if(updateOutfile)  outfileOption = "UPDATE";
2353
2354   TFile out(strOutfile,outfileOption);
2355
2356   if(!out.IsOpen()){
2357     Printf("%s:%d -- error opening response matrix output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
2358     return;
2359   }
2360
2361   if(fDebug>0) Printf("%s:%d -- write response matrices to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
2362
2363   if(strOutDir && strOutDir.Length()){
2364
2365     TDirectory* dir;
2366     if((dir = ((TDirectory*)  gDirectory->Get(strOutDir)))) dir->cd();
2367     else{
2368       dir = out.mkdir(strOutDir);
2369       dir->cd();
2370     }
2371   }
2372
2373   hnResponseSinglePt->Write();
2374   h2ResponseSinglePt->Write();
2375
2376   out.Close();
2377 }
2378
2379 //________________________________________________________________________________________________________________
2380 void AliFragmentationFunctionCorrections::WriteJetTrackEff(TString strInfile, TString strID, TString strOutfile,
2381                                                            Bool_t updateOutfile)
2382 {
2383   // read task ouput from MC and write single track eff - standard dir/list
2384
2385   TString strdir  = "PWG4_FragmentationFunction_" + strID;
2386   TString strlist = "fracfunc_" + strID;
2387
2388   WriteJetTrackEff(strInfile,strdir,strlist,strOutfile,updateOutfile);
2389 }
2390
2391 //___________________________________________________________________________________________________________________________________
2392 void AliFragmentationFunctionCorrections::WriteJetTrackEff(TString strInfile, TString strdir, TString strlist,
2393                                                            TString strOutfile, Bool_t updateOutfile)
2394 {
2395   // read task output from MC and write track eff in jet pt slices as function of pt, z, xi
2396   // argument strlist optional - read from directory strdir if not specified
2397   // write eff to file strOutfile - by default only 'update' file (don't overwrite)
2398
2399   TH1F* hEffPt[fNJetPtSlices];
2400   TH1F* hEffXi[fNJetPtSlices];
2401   TH1F* hEffZ[fNJetPtSlices];
2402
2403   TH1F* hdNdptTracksMCPrimGen[fNJetPtSlices];
2404   TH1F* hdNdxiMCPrimGen[fNJetPtSlices];
2405   TH1F* hdNdzMCPrimGen[fNJetPtSlices];
2406
2407   TH1F* hdNdptTracksMCPrimRec[fNJetPtSlices];
2408   TH1F* hdNdxiMCPrimRec[fNJetPtSlices];
2409   TH1F* hdNdzMCPrimRec[fNJetPtSlices];
2410
2411
2412   TH1F* fh1FFJetPtRecEffGen;
2413
2414   TH2F* fh2FFTrackPtRecEffGen;
2415   TH2F* fh2FFZRecEffGen;
2416   TH2F* fh2FFXiRecEffGen;
2417
2418   TH2F* fh2FFTrackPtRecEffRec;
2419   TH2F* fh2FFZRecEffRec;
2420   TH2F* fh2FFXiRecEffRec;
2421
2422
2423   TFile f(strInfile,"READ");
2424
2425   if(!f.IsOpen()){
2426     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2427     return;
2428   }
2429
2430   if(fDebug>0) Printf("%s:%d -- writeJetTrackEff: open task ouput file %s",(char*)__FILE__,__LINE__,strInfile.Data());
2431
2432   if(strdir && strdir.Length()) gDirectory->cd(strdir);
2433
2434   TList* list = 0;
2435
2436   if(strlist && strlist.Length()){
2437
2438     if(!(list = (TList*) gDirectory->Get(strlist))){
2439       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2440       return;
2441     }
2442   }
2443
2444   if(list){
2445     fh1FFJetPtRecEffGen    = (TH1F*) list->FindObject("fh1FFJetPtRecEffGen");
2446
2447     fh2FFTrackPtRecEffGen  = (TH2F*) list->FindObject("fh2FFTrackPtRecEffGen");
2448     fh2FFZRecEffGen        = (TH2F*) list->FindObject("fh2FFZRecEffGen");
2449     fh2FFXiRecEffGen       = (TH2F*) list->FindObject("fh2FFXiRecEffGen");
2450
2451     fh2FFTrackPtRecEffRec  = (TH2F*) list->FindObject("fh2FFTrackPtRecEffRec");
2452     fh2FFZRecEffRec        = (TH2F*) list->FindObject("fh2FFZRecEffRec");
2453     fh2FFXiRecEffRec       = (TH2F*) list->FindObject("fh2FFXiRecEffRec");
2454   }
2455   else{
2456     fh1FFJetPtRecEffGen    = (TH1F*) gDirectory->FindObject("fh1FFJetPtRecEffGen");
2457
2458     fh2FFTrackPtRecEffGen  = (TH2F*) gDirectory->FindObject("fh2FFTrackPtRecEffGen");
2459     fh2FFZRecEffGen        = (TH2F*) gDirectory->FindObject("fh2FFZRecEffGen");
2460     fh2FFXiRecEffGen       = (TH2F*) gDirectory->FindObject("fh2FFXiRecEffGen");
2461
2462     fh2FFTrackPtRecEffRec  = (TH2F*) gDirectory->FindObject("fh2FFTrackPtRecEffRec");
2463     fh2FFZRecEffRec        = (TH2F*) gDirectory->FindObject("fh2FFZRecEffRec");
2464     fh2FFXiRecEffRec       = (TH2F*) gDirectory->FindObject("fh2FFXiRecEffRec");
2465   }
2466
2467   fh1FFJetPtRecEffGen->SetDirectory(0);
2468
2469   fh2FFTrackPtRecEffGen->SetDirectory(0);
2470   fh2FFZRecEffGen->SetDirectory(0);
2471   fh2FFXiRecEffGen->SetDirectory(0);
2472
2473   fh2FFTrackPtRecEffRec->SetDirectory(0);
2474   fh2FFZRecEffRec->SetDirectory(0);
2475   fh2FFXiRecEffRec->SetDirectory(0);
2476
2477   f.Close();
2478
2479
2480   // projections: FF for generated and reconstructed primaries
2481
2482   for(Int_t i=0; i<fNJetPtSlices; i++){
2483
2484     Float_t jetPtLoLim = fJetPtSlices->At(i);
2485     Float_t jetPtUpLim = fJetPtSlices->At(i+1);
2486
2487     Int_t binLo = static_cast<Int_t>(fh2FFTrackPtRecEffGen->GetXaxis()->FindBin(jetPtLoLim));
2488     Int_t binUp = static_cast<Int_t>(fh2FFTrackPtRecEffGen->GetXaxis()->FindBin(jetPtUpLim))-1;
2489
2490     if(binUp > fh2FFTrackPtRecEffGen->GetNbinsX()){
2491       Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim);
2492       return;
2493     }
2494
2495     TString strNameFFPtGen(Form("fh1FFTrackPtGenPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2496     TString strNameFFZGen(Form("fh1FFZGenPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2497     TString strNameFFXiGen(Form("fh1FFXiGenPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2498
2499     TString strNameFFPtRec(Form("fh1FFTrackPtRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2500     TString strNameFFZRec(Form("fh1FFZRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2501     TString strNameFFXiRec(Form("fh1FFXiRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2502
2503     // project
2504     // appendix 'unbinned' to avoid histos with same name after rebinning
2505
2506     hdNdptTracksMCPrimGen[i] = (TH1F*) fh2FFTrackPtRecEffGen->ProjectionY(strNameFFPtGen+"_unBinned",binLo,binUp,"o"); // option "o": original axis range
2507     hdNdzMCPrimGen[i]        = (TH1F*) fh2FFZRecEffGen->ProjectionY(strNameFFZGen+"_unBinned",binLo,binUp,"o");
2508     hdNdxiMCPrimGen[i]       = (TH1F*) fh2FFXiRecEffGen->ProjectionY(strNameFFXiGen+"_unBinned",binLo,binUp,"o");
2509
2510     hdNdptTracksMCPrimRec[i] = (TH1F*) fh2FFTrackPtRecEffRec->ProjectionY(strNameFFPtRec+"_unBinned",binLo,binUp,"o"); // option "o": original axis range
2511     hdNdzMCPrimRec[i]        = (TH1F*) fh2FFZRecEffRec->ProjectionY(strNameFFZRec+"_unBinned",binLo,binUp,"o");
2512     hdNdxiMCPrimRec[i]       = (TH1F*) fh2FFXiRecEffRec->ProjectionY(strNameFFXiRec+"_unBinned",binLo,binUp,"o");
2513
2514     // rebin
2515
2516     if(fNHistoBinsPt[i]) hdNdptTracksMCPrimGen[i] = (TH1F*) hdNdptTracksMCPrimGen[i]->Rebin(fNHistoBinsPt[i],strNameFFPtGen,fHistoBinsPt[i]->GetArray());
2517     if(fNHistoBinsZ[i])  hdNdzMCPrimGen[i]  = (TH1F*) hdNdzMCPrimGen[i]->Rebin(fNHistoBinsZ[i],strNameFFZGen,fHistoBinsZ[i]->GetArray());
2518     if(fNHistoBinsXi[i]) hdNdxiMCPrimGen[i] = (TH1F*) hdNdxiMCPrimGen[i]->Rebin(fNHistoBinsXi[i],strNameFFXiGen,fHistoBinsXi[i]->GetArray());
2519
2520     if(fNHistoBinsPt[i]) hdNdptTracksMCPrimRec[i] = (TH1F*) hdNdptTracksMCPrimRec[i]->Rebin(fNHistoBinsPt[i],strNameFFPtRec,fHistoBinsPt[i]->GetArray());
2521     if(fNHistoBinsZ[i])  hdNdzMCPrimRec[i]  = (TH1F*) hdNdzMCPrimRec[i]->Rebin(fNHistoBinsZ[i],strNameFFZRec,fHistoBinsZ[i]->GetArray());
2522     if(fNHistoBinsXi[i]) hdNdxiMCPrimRec[i] = (TH1F*) hdNdxiMCPrimRec[i]->Rebin(fNHistoBinsXi[i],strNameFFXiRec,fHistoBinsXi[i]->GetArray());
2523
2524     hdNdptTracksMCPrimGen[i]->SetNameTitle(strNameFFPtGen,"");
2525     hdNdzMCPrimGen[i]->SetNameTitle(strNameFFZGen,"");
2526     hdNdxiMCPrimGen[i]->SetNameTitle(strNameFFXiGen,"");
2527
2528     hdNdptTracksMCPrimRec[i]->SetNameTitle(strNameFFPtRec,"");
2529     hdNdzMCPrimRec[i]->SetNameTitle(strNameFFZRec,"");
2530     hdNdxiMCPrimRec[i]->SetNameTitle(strNameFFXiRec,"");
2531
2532     // normalize
2533
2534     Double_t nJetsBin = fh1FFJetPtRecEffGen->Integral(binLo,binUp);
2535
2536     NormalizeTH1(hdNdptTracksMCPrimGen[i],nJetsBin);
2537     NormalizeTH1(hdNdzMCPrimGen[i],nJetsBin);
2538     NormalizeTH1(hdNdxiMCPrimGen[i],nJetsBin);
2539
2540     NormalizeTH1(hdNdptTracksMCPrimRec[i],nJetsBin);
2541     NormalizeTH1(hdNdzMCPrimRec[i],nJetsBin);
2542     NormalizeTH1(hdNdxiMCPrimRec[i],nJetsBin);
2543
2544     // divide rec/gen : efficiency
2545
2546     TString strNameEffPt(Form("hEffPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2547     TString strNameEffZ(Form("hEffZ_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2548     TString strNameEffXi(Form("hEffXi_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2549
2550     hEffPt[i] = (TH1F*) hdNdptTracksMCPrimRec[i]->Clone(strNameEffPt);
2551     hEffPt[i]->Divide(hdNdptTracksMCPrimRec[i],hdNdptTracksMCPrimGen[i],1,1,"B"); // binominal errors
2552
2553     hEffXi[i] = (TH1F*) hdNdxiMCPrimRec[i]->Clone(strNameEffXi);
2554     hEffXi[i]->Divide(hdNdxiMCPrimRec[i],hdNdxiMCPrimGen[i],1,1,"B"); // binominal errors
2555
2556     hEffZ[i] = (TH1F*) hdNdzMCPrimRec[i]->Clone(strNameEffZ);
2557     hEffZ[i]->Divide(hdNdzMCPrimRec[i],hdNdzMCPrimGen[i],1,1,"B"); // binominal errors
2558   }
2559
2560   // write
2561
2562   TString outfileOption = "RECREATE";
2563   if(updateOutfile)  outfileOption = "UPDATE";
2564
2565   TFile out(strOutfile,outfileOption);
2566
2567   if(!out.IsOpen()){
2568     Printf("%s:%d -- error opening efficiency output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
2569     return;
2570   }
2571
2572   if(fDebug>0) Printf("%s:%d -- write efficiency to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
2573
2574 //   if(strdir && strdir.Length()){
2575 //     TDirectory* dir = out.mkdir(strdir);
2576 //     dir->cd();
2577 //   }
2578
2579   for(Int_t i=0; i<fNJetPtSlices; i++){
2580
2581     hdNdptTracksMCPrimGen[i]->Write();
2582     hdNdxiMCPrimGen[i]->Write();
2583     hdNdzMCPrimGen[i]->Write();
2584
2585     hdNdptTracksMCPrimRec[i]->Write();
2586     hdNdxiMCPrimRec[i]->Write();
2587     hdNdzMCPrimRec[i]->Write();
2588
2589     hEffPt[i]->Write();
2590     hEffXi[i]->Write();
2591     hEffZ[i]->Write();
2592   }
2593
2594   out.Close();
2595
2596 }
2597
2598 //________________________________________________________________________________________________________________
2599 void AliFragmentationFunctionCorrections::WriteJetSecCorr(TString strInfile, TString strID, TString strOutfile,
2600                                                            Bool_t updateOutfile)
2601 {
2602   // read task ouput from MC and write secondary correction - standard dir/list
2603
2604   TString strdir  = "PWG4_FragmentationFunction_" + strID;
2605   TString strlist = "fracfunc_" + strID;
2606
2607   WriteJetSecCorr(strInfile,strdir,strlist,strOutfile,updateOutfile);
2608 }
2609
2610 //___________________________________________________________________________________________________________________________________
2611 void AliFragmentationFunctionCorrections::WriteJetSecCorr(TString strInfile, TString strdir, TString strlist,
2612                                                            TString strOutfile, Bool_t updateOutfile)
2613 {
2614   // read task output from MC and write secondary correction in jet pt slices as function of pt, z, xi
2615   // argument strlist optional - read from directory strdir if not specified
2616   // write eff to file strOutfile - by default only 'update' file (don't overwrite)
2617
2618   TH1F* hSecCorrPt[fNJetPtSlices];
2619   TH1F* hSecCorrXi[fNJetPtSlices];
2620   TH1F* hSecCorrZ[fNJetPtSlices];
2621
2622   TH1F* hdNdptTracksMCPrimRec[fNJetPtSlices];
2623   TH1F* hdNdxiMCPrimRec[fNJetPtSlices];
2624   TH1F* hdNdzMCPrimRec[fNJetPtSlices];
2625
2626   TH1F* hdNdptTracksMCSecRec[fNJetPtSlices];
2627   TH1F* hdNdxiMCSecRec[fNJetPtSlices];
2628   TH1F* hdNdzMCSecRec[fNJetPtSlices];
2629
2630   TH1F* fh1FFJetPtRecEffGen;
2631
2632   TH2F* fh2FFTrackPtRecEffRec;
2633   TH2F* fh2FFZRecEffRec;
2634   TH2F* fh2FFXiRecEffRec;
2635
2636   TH2F* fh2FFTrackPtSecRec;
2637   TH2F* fh2FFZSecRec;
2638   TH2F* fh2FFXiSecRec;
2639
2640   TFile f(strInfile,"READ");
2641
2642   if(!f.IsOpen()){
2643     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2644     return;
2645   }
2646
2647   if(fDebug>0) Printf("%s:%d -- writeJetTrackEff: open task ouput file %s",(char*)__FILE__,__LINE__,strInfile.Data());
2648
2649   if(strdir && strdir.Length()) gDirectory->cd(strdir);
2650
2651   TList* list = 0;
2652
2653   if(strlist && strlist.Length()){
2654
2655     if(!(list = (TList*) gDirectory->Get(strlist))){
2656       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2657       return;
2658     }
2659   }
2660
2661   if(list){
2662     fh1FFJetPtRecEffGen    = (TH1F*) list->FindObject("fh1FFJetPtRecEffGen");
2663
2664     fh2FFTrackPtRecEffRec  = (TH2F*) list->FindObject("fh2FFTrackPtRecEffRec");
2665     fh2FFZRecEffRec        = (TH2F*) list->FindObject("fh2FFZRecEffRec");
2666     fh2FFXiRecEffRec       = (TH2F*) list->FindObject("fh2FFXiRecEffRec");
2667
2668     fh2FFTrackPtSecRec  = (TH2F*) list->FindObject("fh2FFTrackPtSecRec");
2669     fh2FFZSecRec        = (TH2F*) list->FindObject("fh2FFZSecRec");
2670     fh2FFXiSecRec       = (TH2F*) list->FindObject("fh2FFXiSecRec");
2671   }
2672   else{
2673     fh1FFJetPtRecEffGen    = (TH1F*) gDirectory->FindObject("fh1FFJetPtRecEffGen");
2674
2675     fh2FFTrackPtRecEffRec  = (TH2F*) gDirectory->FindObject("fh2FFTrackPtRecEffRec");
2676     fh2FFZRecEffRec        = (TH2F*) gDirectory->FindObject("fh2FFZRecEffRec");
2677     fh2FFXiRecEffRec       = (TH2F*) gDirectory->FindObject("fh2FFXiRecEffRec");
2678
2679     fh2FFTrackPtSecRec  = (TH2F*) gDirectory->FindObject("fh2FFTrackPtSecRec");
2680     fh2FFZSecRec        = (TH2F*) gDirectory->FindObject("fh2FFZSecRec");
2681     fh2FFXiSecRec       = (TH2F*) gDirectory->FindObject("fh2FFXiSecRec");
2682   }
2683
2684   fh1FFJetPtRecEffGen->SetDirectory(0);
2685
2686   fh2FFTrackPtRecEffRec->SetDirectory(0);
2687   fh2FFZRecEffRec->SetDirectory(0);
2688   fh2FFXiRecEffRec->SetDirectory(0);
2689
2690   fh2FFTrackPtSecRec->SetDirectory(0);
2691   fh2FFZSecRec->SetDirectory(0);
2692   fh2FFXiSecRec->SetDirectory(0);
2693
2694   f.Close();
2695
2696
2697   // projections: FF for generated and reconstructed primaries
2698
2699   for(Int_t i=0; i<fNJetPtSlices; i++){
2700
2701     Float_t jetPtLoLim = fJetPtSlices->At(i);
2702     Float_t jetPtUpLim = fJetPtSlices->At(i+1);
2703
2704     Int_t binLo = static_cast<Int_t>(fh2FFTrackPtRecEffRec->GetXaxis()->FindBin(jetPtLoLim));
2705     Int_t binUp = static_cast<Int_t>(fh2FFTrackPtRecEffRec->GetXaxis()->FindBin(jetPtUpLim))-1;
2706
2707     if(binUp > fh2FFTrackPtRecEffRec->GetNbinsX()){
2708       Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim);
2709       return;
2710     }
2711
2712     TString strNameFFPtPrimRec(Form("fh1FFTrackPtRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2713     TString strNameFFZPrimRec(Form("fh1FFZRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2714     TString strNameFFXiPrimRec(Form("fh1FFXiRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2715
2716     TString strNameFFPtSecRec(Form("fh1FFTrackPtRecSec_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2717     TString strNameFFZSecRec(Form("fh1FFZRecSec_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2718     TString strNameFFXiSecRec(Form("fh1FFXiRecSec_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2719
2720     // project
2721     // appendix 'unbinned' to avoid histos with same name after rebinning
2722
2723     hdNdptTracksMCPrimRec[i] = (TH1F*) fh2FFTrackPtRecEffRec->ProjectionY(strNameFFPtPrimRec+"_unBinned",binLo,binUp,"o"); // option "o": original axis range
2724     hdNdzMCPrimRec[i]        = (TH1F*) fh2FFZRecEffRec->ProjectionY(strNameFFZPrimRec+"_unBinned",binLo,binUp,"o");
2725     hdNdxiMCPrimRec[i]       = (TH1F*) fh2FFXiRecEffRec->ProjectionY(strNameFFXiPrimRec+"_unBinned",binLo,binUp,"o");
2726
2727     hdNdptTracksMCSecRec[i]  = (TH1F*) fh2FFTrackPtSecRec->ProjectionY(strNameFFPtSecRec+"_unBinned",binLo,binUp,"o"); // option "o": original axis range
2728     hdNdzMCSecRec[i]         = (TH1F*) fh2FFZSecRec->ProjectionY(strNameFFZSecRec+"_unBinned",binLo,binUp,"o");
2729     hdNdxiMCSecRec[i]        = (TH1F*) fh2FFXiSecRec->ProjectionY(strNameFFXiSecRec+"_unBinned",binLo,binUp,"o");
2730
2731     // rebin
2732
2733     if(fNHistoBinsPt[i]) hdNdptTracksMCPrimRec[i] = (TH1F*) hdNdptTracksMCPrimRec[i]->Rebin(fNHistoBinsPt[i],strNameFFPtPrimRec,fHistoBinsPt[i]->GetArray());
2734     if(fNHistoBinsZ[i])  hdNdzMCPrimRec[i]        = (TH1F*) hdNdzMCPrimRec[i]->Rebin(fNHistoBinsZ[i],strNameFFZPrimRec,fHistoBinsZ[i]->GetArray());
2735     if(fNHistoBinsXi[i]) hdNdxiMCPrimRec[i]       = (TH1F*) hdNdxiMCPrimRec[i]->Rebin(fNHistoBinsXi[i],strNameFFXiPrimRec,fHistoBinsXi[i]->GetArray());
2736
2737     if(fNHistoBinsPt[i]) hdNdptTracksMCSecRec[i] = (TH1F*) hdNdptTracksMCSecRec[i]->Rebin(fNHistoBinsPt[i],strNameFFPtSecRec,fHistoBinsPt[i]->GetArray());
2738     if(fNHistoBinsZ[i])  hdNdzMCSecRec[i]        = (TH1F*) hdNdzMCSecRec[i]->Rebin(fNHistoBinsZ[i],strNameFFZSecRec,fHistoBinsZ[i]->GetArray());
2739     if(fNHistoBinsXi[i]) hdNdxiMCSecRec[i]       = (TH1F*) hdNdxiMCSecRec[i]->Rebin(fNHistoBinsXi[i],strNameFFXiSecRec,fHistoBinsXi[i]->GetArray());
2740
2741     hdNdptTracksMCPrimRec[i]->SetNameTitle(strNameFFPtPrimRec,"");
2742     hdNdzMCPrimRec[i]->SetNameTitle(strNameFFZPrimRec,"");
2743     hdNdxiMCPrimRec[i]->SetNameTitle(strNameFFXiPrimRec,"");
2744
2745     hdNdptTracksMCSecRec[i]->SetNameTitle(strNameFFPtSecRec,"");
2746     hdNdzMCSecRec[i]->SetNameTitle(strNameFFZSecRec,"");
2747     hdNdxiMCSecRec[i]->SetNameTitle(strNameFFXiSecRec,"");
2748
2749     // normalize
2750
2751     Double_t nJetsBin = fh1FFJetPtRecEffGen->Integral(binLo,binUp);
2752
2753     NormalizeTH1(hdNdptTracksMCPrimRec[i],nJetsBin);
2754     NormalizeTH1(hdNdzMCPrimRec[i],nJetsBin);
2755     NormalizeTH1(hdNdxiMCPrimRec[i],nJetsBin);
2756
2757     NormalizeTH1(hdNdptTracksMCSecRec[i],nJetsBin);
2758     NormalizeTH1(hdNdzMCSecRec[i],nJetsBin);
2759     NormalizeTH1(hdNdxiMCSecRec[i],nJetsBin);
2760
2761     // divide rec/gen : efficiency
2762
2763     TString strNameSecCorrPt(Form("hSecCorrPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2764     TString strNameSecCorrZ(Form("hSecCorrZ_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2765     TString strNameSecCorrXi(Form("hSecCorrXi_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2766
2767     hSecCorrPt[i] = (TH1F*) hdNdptTracksMCSecRec[i]->Clone(strNameSecCorrPt);
2768     TH1F* hSumPrimSecPt = (TH1F*) hdNdptTracksMCSecRec[i]->Clone("hSumPrimSecPt");
2769     hSumPrimSecPt->Add(hdNdptTracksMCPrimRec[i]);
2770     hSecCorrPt[i]->Divide(hdNdptTracksMCPrimRec[i],hSumPrimSecPt,1,1,"B"); // binominal errors
2771
2772     hSecCorrXi[i] = (TH1F*) hdNdxiMCSecRec[i]->Clone(strNameSecCorrXi);
2773     TH1F* hSumPrimSecXi = (TH1F*) hdNdxiMCSecRec[i]->Clone("hSumPrimSecXi");
2774     hSumPrimSecXi->Add(hdNdxiMCPrimRec[i]);
2775     hSecCorrXi[i]->Divide(hdNdxiMCPrimRec[i],hSumPrimSecXi,1,1,"B"); // binominal errors
2776
2777     hSecCorrZ[i] = (TH1F*) hdNdzMCSecRec[i]->Clone(strNameSecCorrZ);
2778     TH1F* hSumPrimSecZ = (TH1F*) hdNdzMCSecRec[i]->Clone("hSumPrimSecZ");
2779     hSumPrimSecZ->Add(hdNdzMCPrimRec[i]);
2780     hSecCorrZ[i]->Divide(hdNdzMCPrimRec[i],hSumPrimSecZ,1,1,"B"); // binominal errors
2781   }
2782
2783   // write
2784
2785   TString outfileOption = "RECREATE";
2786   if(updateOutfile)  outfileOption = "UPDATE";
2787
2788   TFile out(strOutfile,outfileOption);
2789
2790   if(!out.IsOpen()){
2791     Printf("%s:%d -- error opening efficiency output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
2792     return;
2793   }
2794
2795   if(fDebug>0) Printf("%s:%d -- write efficiency to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
2796
2797   for(Int_t i=0; i<fNJetPtSlices; i++){
2798
2799     //  hdNdptTracksMCSecRec[i]->Write();
2800     //  hdNdxiMCSecRec[i]->Write();
2801     //  hdNdzMCSecRec[i]->Write();
2802
2803     hSecCorrPt[i]->Write();
2804     hSecCorrXi[i]->Write();
2805     hSecCorrZ[i]->Write();
2806   }
2807
2808   out.Close();
2809 }
2810
2811 //________________________________________________________________________________________________________________
2812 void AliFragmentationFunctionCorrections::WriteJetResponse(TString strInfile, TString strID, TString strOutfile,
2813                                                            Bool_t updateOutfile)
2814 {
2815   // read task ouput from MC and write single track eff - standard dir/list
2816
2817   TString strdir  = "PWG4_FragmentationFunction_" + strID;
2818   TString strlist = "fracfunc_" + strID;
2819
2820   WriteJetResponse(strInfile,strdir,strlist,strOutfile,updateOutfile);
2821 }
2822
2823 //_____________________________________________________________________________________________________________________________________
2824 void AliFragmentationFunctionCorrections::WriteJetResponse(TString strInfile, TString strdir, TString strlist,
2825                                                            TString strOutfile, Bool_t updateOutfile)
2826 {
2827   // read 3d THnSparse response matrices in pt,z,xi vs jet pt from file
2828   // project THnSparse + TH2 in jet pt slices
2829   // write to strOutfile
2830
2831   THnSparse* hn3ResponseJetPt;
2832   THnSparse* hn3ResponseJetZ;
2833   THnSparse* hn3ResponseJetXi;
2834
2835   // 2D response matrices
2836
2837   THnSparse* hnResponsePt[fNJetPtSlices];
2838   THnSparse* hnResponseZ[fNJetPtSlices];
2839   THnSparse* hnResponseXi[fNJetPtSlices];
2840
2841   TH2F* h2ResponsePt[fNJetPtSlices];
2842   TH2F* h2ResponseZ[fNJetPtSlices];
2843   TH2F* h2ResponseXi[fNJetPtSlices];
2844
2845   // 1D projections on gen pt / rec pt axes
2846
2847   TH1F* h1FFPtRec[fNJetPtSlices];
2848   TH1F* h1FFZRec[fNJetPtSlices];
2849   TH1F* h1FFXiRec[fNJetPtSlices];
2850
2851   TH1F* h1FFPtGen[fNJetPtSlices];
2852   TH1F* h1FFZGen[fNJetPtSlices];
2853   TH1F* h1FFXiGen[fNJetPtSlices];
2854
2855   TH1F* h1RatioPt[fNJetPtSlices];
2856   TH1F* h1RatioZ[fNJetPtSlices];
2857   TH1F* h1RatioXi[fNJetPtSlices];
2858
2859
2860
2861   TFile f(strInfile,"READ");
2862
2863   if(!f.IsOpen()){
2864     Printf("%s:%d -- error opening file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2865     return;
2866   }
2867
2868   if(fDebug>0) Printf("%s:%d -- read response matrices from file %s ",(char*)__FILE__,__LINE__,strInfile.Data());
2869
2870   if(strdir && strdir.Length()) gDirectory->cd(strdir);
2871
2872   TList* list = 0;
2873
2874   if(strlist && strlist.Length()){
2875
2876     if(!(list = (TList*) gDirectory->Get(strlist))){
2877       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2878       return;
2879     }
2880   }
2881
2882   if(list){
2883     hn3ResponseJetPt = (THnSparse*) list->FindObject("fhnResponseJetTrackPt");
2884     hn3ResponseJetZ  = (THnSparse*) list->FindObject("fhnResponseJetZ");
2885     hn3ResponseJetXi = (THnSparse*) list->FindObject("fhnResponseJetXi");
2886   }
2887   else{
2888     hn3ResponseJetPt = (THnSparse*) gDirectory->Get("fhnResponseJetTrackPt");
2889     hn3ResponseJetZ  = (THnSparse*) gDirectory->Get("fhnResponseJetZ");
2890     hn3ResponseJetXi = (THnSparse*) gDirectory->Get("fhnResponseJetXi");
2891   }
2892
2893
2894   if(!hn3ResponseJetPt){
2895     Printf("%s:%d -- error retrieving response matrix fhnResponseJetTrackPt",(char*)__FILE__,__LINE__);
2896     return;
2897   }
2898
2899   if(!hn3ResponseJetZ){
2900     Printf("%s:%d -- error retrieving response matrix fhnResponseJetZ",(char*)__FILE__,__LINE__);
2901     return;
2902   }
2903
2904   if(!hn3ResponseJetXi){
2905     Printf("%s:%d -- error retrieving response matrix fhnResponseJetXi",(char*)__FILE__,__LINE__);
2906     return;
2907   }
2908
2909   f.Close();
2910
2911   // axes
2912
2913   Int_t axisJetPtTHn3 = -1;
2914   Int_t axisGenPtTHn3 = -1;
2915   Int_t axisRecPtTHn3 = -1;
2916
2917   for(Int_t i=0; i<hn3ResponseJetPt->GetNdimensions(); i++){
2918
2919     TString title = hn3ResponseJetPt->GetAxis(i)->GetTitle();
2920
2921     if(title.Contains("jet p_{T}")) axisJetPtTHn3 = i;
2922     if(title.Contains("gen p_{T}")) axisGenPtTHn3 = i;
2923     if(title.Contains("rec p_{T}")) axisRecPtTHn3 = i;
2924   }
2925
2926   if(axisJetPtTHn3 == -1){
2927     Printf("%s:%d -- error axisJetPtTHn3",(char*)__FILE__,__LINE__);
2928     return;
2929   }
2930
2931   if(axisGenPtTHn3 == -1){
2932     Printf("%s:%d -- error axisGenPtTHn3",(char*)__FILE__,__LINE__);
2933     return;
2934   }
2935
2936   if(axisRecPtTHn3 == -1){
2937     Printf("%s:%d -- error axisRecPtTHn3",(char*)__FILE__,__LINE__);
2938     return;
2939   }
2940
2941   for(Int_t i=0; i<fNJetPtSlices; i++){
2942
2943     Float_t jetPtLoLim = fJetPtSlices->At(i);
2944     Float_t jetPtUpLim = fJetPtSlices->At(i+1);
2945
2946     Int_t binLo = static_cast<Int_t>(hn3ResponseJetPt->GetAxis(axisJetPtTHn3)->FindBin(jetPtLoLim));
2947     Int_t binUp = static_cast<Int_t>(hn3ResponseJetPt->GetAxis(axisJetPtTHn3)->FindBin(jetPtUpLim))-1;
2948
2949     if(binUp > hn3ResponseJetPt->GetAxis(axisJetPtTHn3)->GetNbins()){
2950       Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim);
2951       return;
2952     }
2953
2954     TString strNameRespPt(Form("hnResponsePt_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2955     TString strNameRespZ(Form("hnResponseZ_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2956     TString strNameRespXi(Form("hnResponseXi_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2957
2958     TString strNameTH2RespPt(Form("h2ResponsePt_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2959     TString strNameTH2RespZ(Form("h2ResponseZ_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2960     TString strNameTH2RespXi(Form("h2ResponseXi_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2961
2962     TString strNameRecPt(Form("h1FFTrackPtRecPrim_recPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2963     TString strNameRecZ(Form("h1FFZRecPrim_recPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2964     TString strNameRecXi(Form("h1FFXiRecPrim_recPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2965
2966     TString strNameGenPt(Form("h1FFTrackPtRecPrim_genPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2967     TString strNameGenZ(Form("h1FFZRecPrim_genPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2968     TString strNameGenXi(Form("h1FFXiRecPrim_genPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2969
2970     TString strNameRatioPt(Form("h1RatioTrackPtRecPrim_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2971     TString strNameRatioZ(Form("h1RatioZRecPrim_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2972     TString strNameRatioXi(Form("h1RatioXiRecPrim_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2973
2974
2975     // 2D projections in jet pt range
2976
2977     hn3ResponseJetPt->GetAxis(axisJetPtTHn3)->SetRange(binLo,binUp);
2978     hn3ResponseJetZ->GetAxis(axisJetPtTHn3)->SetRange(binLo,binUp);
2979     hn3ResponseJetXi->GetAxis(axisJetPtTHn3)->SetRange(binLo,binUp);
2980
2981     Int_t axesProj[2] = {axisRecPtTHn3,axisGenPtTHn3};
2982
2983     hnResponsePt[i] = hn3ResponseJetPt->Projection(2,axesProj);
2984     hnResponseZ[i]  = hn3ResponseJetZ->Projection(2,axesProj);
2985     hnResponseXi[i] = hn3ResponseJetXi->Projection(2,axesProj);
2986
2987     hnResponsePt[i]->SetNameTitle(strNameRespPt,"");
2988     hnResponseZ[i]->SetNameTitle(strNameRespZ,"");
2989     hnResponseXi[i]->SetNameTitle(strNameRespXi,"");
2990
2991     h2ResponsePt[i] = (TH2F*) hnResponsePt[i]->Projection(1,0);// note convention: yDim,xDim
2992     h2ResponseZ[i]  = (TH2F*) hnResponseZ[i]->Projection(1,0); // note convention: yDim,xDim
2993     h2ResponseXi[i] = (TH2F*) hnResponseXi[i]->Projection(1,0);// note convention: yDim,xDim
2994
2995     h2ResponsePt[i]->SetNameTitle(strNameTH2RespPt,"");
2996     h2ResponseZ[i]->SetNameTitle(strNameTH2RespZ,"");
2997     h2ResponseXi[i]->SetNameTitle(strNameTH2RespXi,"");
2998
2999
3000     // 1D projections
3001
3002     Int_t axisGenPtTHn2 = -1;
3003     Int_t axisRecPtTHn2 = -1;
3004
3005     for(Int_t d=0; d<hnResponsePt[i]->GetNdimensions(); d++){
3006
3007       TString title = hnResponsePt[i]->GetAxis(d)->GetTitle();
3008
3009       if(title.Contains("gen p_{T}")) axisGenPtTHn2 = d;
3010       if(title.Contains("rec p_{T}")) axisRecPtTHn2 = d;
3011     }
3012
3013
3014     if(axisGenPtTHn2 == -1){
3015       Printf("%s:%d -- error axisGenPtTHn2",(char*)__FILE__,__LINE__);
3016       return;
3017     }
3018
3019     if(axisRecPtTHn2 == -1){
3020       Printf("%s:%d -- error axisRecPtTHn2",(char*)__FILE__,__LINE__);
3021       return;
3022     }
3023
3024
3025     h1FFPtRec[i] = (TH1F*) hnResponsePt[i]->Projection(axisRecPtTHn2);// note convention: yDim,xDim
3026     h1FFZRec[i]  = (TH1F*) hnResponseZ[i]->Projection(axisRecPtTHn2);// note convention: yDim,xDim
3027     h1FFXiRec[i] = (TH1F*) hnResponseXi[i]->Projection(axisRecPtTHn2);// note convention: yDim,xDim
3028
3029     h1FFPtRec[i]->SetNameTitle(strNameRecPt,"");
3030     h1FFZRec[i]->SetNameTitle(strNameRecZ,"");
3031     h1FFXiRec[i]->SetNameTitle(strNameRecXi,"");
3032
3033     h1FFPtGen[i] = (TH1F*) hnResponsePt[i]->Projection(axisGenPtTHn2);// note convention: yDim,xDim
3034     h1FFZGen[i]  = (TH1F*) hnResponseZ[i]->Projection(axisGenPtTHn2);// note convention: yDim,xDim
3035     h1FFXiGen[i] = (TH1F*) hnResponseXi[i]->Projection(axisGenPtTHn2);// note convention: yDim,xDim
3036
3037     h1FFPtGen[i]->SetNameTitle(strNameGenPt,"");
3038     h1FFZGen[i]->SetNameTitle(strNameGenZ,"");
3039     h1FFXiGen[i]->SetNameTitle(strNameGenXi,"");
3040
3041     // normalize 1D projections
3042
3043     if(fNHistoBinsPt[i]) h1FFPtRec[i] = (TH1F*) h1FFPtRec[i]->Rebin(fNHistoBinsPt[i],"",fHistoBinsPt[i]->GetArray());
3044     if(fNHistoBinsZ[i])  h1FFZRec[i]  = (TH1F*) h1FFZRec[i]->Rebin(fNHistoBinsZ[i],"",fHistoBinsZ[i]->GetArray());
3045     if(fNHistoBinsXi[i]) h1FFXiRec[i] = (TH1F*) h1FFXiRec[i]->Rebin(fNHistoBinsXi[i],"",fHistoBinsXi[i]->GetArray());
3046
3047     if(fNHistoBinsPt[i]) h1FFPtGen[i] = (TH1F*) h1FFPtGen[i]->Rebin(fNHistoBinsPt[i],"",fHistoBinsPt[i]->GetArray());
3048     if(fNHistoBinsZ[i])  h1FFZGen[i]  = (TH1F*) h1FFZGen[i]->Rebin(fNHistoBinsZ[i],"",fHistoBinsZ[i]->GetArray());
3049     if(fNHistoBinsXi[i]) h1FFXiGen[i] = (TH1F*) h1FFXiGen[i]->Rebin(fNHistoBinsXi[i],"",fHistoBinsXi[i]->GetArray());
3050
3051     NormalizeTH1(h1FFPtRec[i],fNJets->At(i));
3052     NormalizeTH1(h1FFZRec[i],fNJets->At(i));
3053     NormalizeTH1(h1FFXiRec[i],fNJets->At(i));
3054
3055     NormalizeTH1(h1FFPtGen[i],fNJets->At(i));
3056     NormalizeTH1(h1FFZGen[i],fNJets->At(i));
3057     NormalizeTH1(h1FFXiGen[i],fNJets->At(i));
3058
3059     // ratios 1D projections
3060
3061     h1RatioPt[i] = (TH1F*) h1FFPtRec[i]->Clone(strNameRatioPt);
3062     h1RatioPt[i]->Reset();
3063     h1RatioPt[i]->Divide(h1FFPtRec[i],h1FFPtGen[i],1,1,"B");
3064
3065     h1RatioZ[i] = (TH1F*) h1FFZRec[i]->Clone(strNameRatioZ);
3066     h1RatioZ[i]->Reset();
3067     h1RatioZ[i]->Divide(h1FFZRec[i],h1FFZGen[i],1,1,"B");
3068
3069     h1RatioXi[i] = (TH1F*) h1FFXiRec[i]->Clone(strNameRatioXi);
3070     h1RatioXi[i]->Reset();
3071     h1RatioXi[i]->Divide(h1FFXiRec[i],h1FFXiGen[i],1,1,"B");
3072   }
3073
3074
3075   // write
3076
3077   TString outfileOption = "RECREATE";
3078   if(updateOutfile)  outfileOption = "UPDATE";
3079
3080   TFile out(strOutfile,outfileOption);
3081
3082   if(!out.IsOpen()){
3083     Printf("%s:%d -- error opening response matrix output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
3084     return;
3085   }
3086
3087   if(fDebug>0) Printf("%s:%d -- write response matrices to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
3088
3089   //if(strdir && strdir.Length()){
3090   //  TDirectory* dir = out.mkdir(strdir);
3091   //  dir->cd();
3092   //}
3093
3094   for(Int_t i=0; i<fNJetPtSlices; i++){
3095
3096     hnResponsePt[i]->Write();
3097     hnResponseXi[i]->Write();
3098     hnResponseZ[i]->Write();
3099
3100     h2ResponsePt[i]->Write();
3101     h2ResponseXi[i]->Write();
3102     h2ResponseZ[i]->Write();
3103
3104     h1FFPtRec[i]->Write();
3105     h1FFZRec[i]->Write();
3106     h1FFXiRec[i]->Write();
3107
3108     h1FFPtGen[i]->Write();
3109     h1FFZGen[i]->Write();
3110     h1FFXiGen[i]->Write();
3111
3112     h1RatioPt[i]->Write();
3113     h1RatioZ[i]->Write();
3114     h1RatioXi[i]->Write();
3115
3116   }
3117
3118   out.Close();
3119 }
3120
3121 //______________________________________________________________________________________________________
3122 void AliFragmentationFunctionCorrections::ReadResponse(TString strfile, TString strdir, TString strlist)
3123 {
3124   // read response matrices from file
3125   // argument strlist optional - read from directory strdir if not specified
3126   // note: THnSparse are not rebinned
3127
3128   THnSparse* hRespPt[fNJetPtSlices];
3129   THnSparse* hRespZ[fNJetPtSlices];
3130   THnSparse* hRespXi[fNJetPtSlices];
3131
3132   for(Int_t i=0; i<fNJetPtSlices; i++) hRespPt[i] = 0;
3133   for(Int_t i=0; i<fNJetPtSlices; i++) hRespZ[i]  = 0;
3134   for(Int_t i=0; i<fNJetPtSlices; i++) hRespXi[i] = 0;
3135
3136   TFile f(strfile,"READ");
3137
3138   if(!f.IsOpen()){
3139     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
3140     return;
3141   }
3142
3143   if(fDebug>0) Printf("%s:%d -- read response matrices from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
3144
3145   if(strdir && strdir.Length()) gDirectory->cd(strdir);
3146
3147   TList* list = 0;
3148
3149   if(strlist && strlist.Length()){
3150
3151     if(!(list = (TList*) gDirectory->Get(strlist))){
3152       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
3153       return;
3154     }
3155   }
3156
3157   for(Int_t i=0; i<fNJetPtSlices; i++){
3158
3159     Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
3160     Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
3161
3162     TString strNameRespPt(Form("hnResponsePt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3163     TString strNameRespZ(Form("hnResponseZ_%02d_%02d",jetPtLoLim,jetPtUpLim));
3164     TString strNameRespXi(Form("hnResponseXi_%02d_%02d",jetPtLoLim,jetPtUpLim));
3165
3166     if(list){
3167       hRespPt[i] = (THnSparse*) list->FindObject(strNameRespPt);
3168       hRespZ[i]  = (THnSparse*) list->FindObject(strNameRespZ);
3169       hRespXi[i] = (THnSparse*) list->FindObject(strNameRespXi);
3170     }
3171     else{
3172       hRespPt[i] = (THnSparse*) gDirectory->Get(strNameRespPt);
3173       hRespZ[i]  = (THnSparse*) gDirectory->Get(strNameRespZ);
3174       hRespXi[i] = (THnSparse*) gDirectory->Get(strNameRespXi);
3175     }
3176
3177     if(!hRespPt[i]){
3178       Printf("%s:%d -- error retrieving response %s", (char*)__FILE__,__LINE__,strNameRespPt.Data());
3179     }
3180
3181     if(!hRespZ[i]){
3182       Printf("%s:%d -- error retrieving response %s", (char*)__FILE__,__LINE__,strNameRespZ.Data());
3183     }
3184
3185     if(!hRespXi[i]){
3186       Printf("%s:%d -- error retrieving response %s", (char*)__FILE__,__LINE__,strNameRespXi.Data());
3187     }
3188
3189     //    if(0){ // can't rebin THnSparse ...
3190     //       if(fNHistoBinsPt[i]) hRespPt[i]->SetBinEdges(0,fHistoBinsPt[i]->GetArray());
3191     //       if(fNHistoBinsPt[i]) hRespPt[i]->SetBinEdges(1,fHistoBinsPt[i]->GetArray());
3192
3193     //       if(fNHistoBinsZ[i])  hRespZ[i]->SetBinEdges(0,fHistoBinsZ[i]->GetArray());
3194     //       if(fNHistoBinsZ[i])  hRespZ[i]->SetBinEdges(1,fHistoBinsZ[i]->GetArray());
3195
3196     //       if(fNHistoBinsXi[i]) hRespXi[i]->SetBinEdges(0,fHistoBinsXi[i]->GetArray());
3197     //       if(fNHistoBinsXi[i]) hRespXi[i]->SetBinEdges(1,fHistoBinsXi[i]->GetArray());
3198     //     }
3199
3200
3201   } // jet slices loop
3202
3203   f.Close();  // THnSparse pointers still valid even if file closed
3204
3205 //   for(Int_t i=0; i<fNJetPtSlices; i++){  // no copy c'tor ...
3206 //     if(hRespPt[i]) new(fhnResponsePt[i]) THnSparseF(*hRespPt[i]);
3207 //     if(hRespZ[i])  new(fhnResponseZ[i])  THnSparseF(*hRespZ[i]);
3208 //     if(hRespXi[i]) new(fhnResponseXi[i]) THnSparseF(*hRespXi[i]);
3209 //   }
3210
3211   for(Int_t i=0; i<fNJetPtSlices; i++){
3212     fhnResponsePt[i] = hRespPt[i];
3213     fhnResponseZ[i]  = hRespZ[i];
3214     fhnResponseXi[i] = hRespXi[i];
3215   }
3216 }
3217
3218 //______________________________________________________________________________________________________________________
3219 void AliFragmentationFunctionCorrections::ReadPriors(TString strfile,const Int_t type)
3220 {
3221   // read priors from file: rec primaries, gen pt dist
3222
3223   if(fDebug>0) Printf("%s:%d -- read priors from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
3224
3225   // temporary histos to store pointers from file
3226   TH1F* hist[fNJetPtSlices];
3227
3228   for(Int_t i=0; i<fNJetPtSlices; i++) hist[i] = 0;
3229
3230   TFile f(strfile,"READ");
3231
3232   if(!f.IsOpen()){
3233     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
3234     return;
3235   }
3236
3237   for(Int_t i=0; i<fNJetPtSlices; i++){
3238
3239     Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
3240     Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
3241
3242     TString strName;
3243
3244     if(type == kFlagPt) strName.Form("h1FFTrackPtRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim);
3245     if(type == kFlagZ)  strName.Form("h1FFZRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim);
3246     if(type == kFlagXi) strName.Form("h1FFXiRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim);
3247
3248     hist[i] = (TH1F*) gDirectory->Get(strName);
3249
3250     if(!hist[i]){
3251       Printf("%s:%d -- error retrieving prior %s", (char*)__FILE__,__LINE__,strName.Data());
3252     }
3253
3254
3255     //if(fNHistoBinsPt[i]) hist[i] = (TH1F*) hist[i]->Rebin(fNHistoBinsPt[i],hist[i]->GetName()+"_rebin",fHistoBinsPt[i]->GetArray());
3256
3257     if(hist[i]) hist[i]->SetDirectory(0);
3258
3259   } // jet slices loop
3260
3261   f.Close();
3262
3263
3264   for(Int_t i=0; i<fNJetPtSlices; i++){ // 2nd loop: need to close input file before placing histos
3265     if(hist[i] && type == kFlagPt) new(fh1FFTrackPtPrior[i]) TH1F(*hist[i]);
3266     if(hist[i] && type == kFlagZ)  new(fh1FFZPrior[i])       TH1F(*hist[i]);
3267     if(hist[i] && type == kFlagXi) new(fh1FFXiPrior[i])      TH1F(*hist[i]);
3268   }
3269 }
3270
3271 //_____________________________________________________
3272 // void AliFragmentationFunctionCorrections::RatioRecGen()
3273 // {
3274 //   // create ratio reconstructed over generated FF
3275 //   // use current highest corrLevel
3276
3277 //   Printf("%s:%d -- build ratio rec.gen, corrLevel %d",(char*)__FILE__,__LINE__,fNCorrectionLevels-1);
3278
3279 //   for(Int_t i=0; i<fNJetPtSlices; i++){
3280
3281 //     TH1F* histPtRec =  fCorrFF[fNCorrectionLevels-1]->GetTrackPt(i); // levels -1: latest corr level
3282 //     TH1F* histZRec  =  fCorrFF[fNCorrectionLevels-1]->GetZ(i);       // levels -1: latest corr level
3283 //     TH1F* histXiRec =  fCorrFF[fNCorrectionLevels-1]->GetXi(i);      // levels -1: latest corr level
3284
3285 //     TH1F* histPtGen = fh1FFTrackPtGenPrim[i];
3286 //     TH1F* histZGen  = fh1FFZGenPrim[i];
3287 //     TH1F* histXiGen = fh1FFXiGenPrim[i];
3288
3289 //     TString histNamePt = histPtRec->GetName();
3290 //     TString histNameZ  = histZRec->GetName();
3291 //     TString histNameXi = histXiRec->GetName();
3292
3293 //     histNamePt.ReplaceAll("fh1FF","fh1FFRatioRecGen");
3294 //     histNameZ.ReplaceAll("fh1FF","fh1FFRatioRecGen");
3295 //     histNameXi.ReplaceAll("fh1FF","fh1FFRatioRecGen");
3296
3297 //     // ratio
3298 //     TH1F* hRatioRecGenPt = (TH1F*) histPtRec->Clone(histNamePt);
3299 //     hRatioRecGenPt->Reset();
3300 //     hRatioRecGenPt->Divide(histPtRec,histPtGen,1,1,"B");
3301
3302 //     TH1F* hRatioRecGenZ = (TH1F*) histZRec->Clone(histNameZ);
3303 //     hRatioRecGenZ->Reset();
3304 //     hRatioRecGenZ->Divide(histZRec,histZGen,1,1,"B");
3305
3306 //     TH1F* hRatioRecGenXi = (TH1F*) histXiRec->Clone(histNameXi);
3307 //     hRatioRecGenXi->Reset();
3308 //     hRatioRecGenXi->Divide(histXiRec,histXiGen,1,1,"B");
3309
3310 //     new(fh1FFRatioRecGenPt[i]) TH1F(*hRatioRecGenPt);
3311 //     new(fh1FFRatioRecGenZ[i])  TH1F(*hRatioRecGenZ);
3312 //     new(fh1FFRatioRecGenXi[i]) TH1F(*hRatioRecGenXi);
3313 //   }
3314 // }
3315
3316 // //___________________________________________________________
3317 // void AliFragmentationFunctionCorrections::RatioRecPrimaries()
3318 // {
3319 //   // create ratio reconstructed tracks over reconstructed primaries
3320 //   // use raw FF (corrLevel 0)
3321
3322 //   Printf("%s:%d -- build ratio rec tracks /rec primaries",(char*)__FILE__,__LINE__);
3323
3324 //   for(Int_t i=0; i<fNJetPtSlices; i++){
3325
3326 //     const Int_t corrLevel = 0;
3327
3328 //     TH1F* histPtRec =  fCorrFF[corrLevel]->GetTrackPt(i); // levels -1: latest corr level
3329 //     TH1F* histZRec  =  fCorrFF[corrLevel]->GetZ(i);       // levels -1: latest corr level
3330 //     TH1F* histXiRec =  fCorrFF[corrLevel]->GetXi(i);      // levels -1: latest corr level
3331
3332 //     TH1F* histPtRecPrim = fh1FFTrackPtRecPrim[i];
3333 //     TH1F* histZRecPrim  = fh1FFZRecPrim[i];
3334 //     TH1F* histXiRecPrim = fh1FFXiRecPrim[i];
3335
3336 //     TString histNamePt = histPtRec->GetName();
3337 //     TString histNameZ  = histZRec->GetName();
3338 //     TString histNameXi = histXiRec->GetName();
3339
3340 //     histNamePt.ReplaceAll("fh1FF","fh1FFRatioRecPrim");
3341 //     histNameZ.ReplaceAll("fh1FF","fh1FFRatioRecPrim");
3342 //     histNameXi.ReplaceAll("fh1FF","fh1FFRatioRecPrim");
3343
3344 //     // ratio
3345 //     TH1F* hRatioRecPrimPt = (TH1F*) histPtRec->Clone(histNamePt);
3346 //     hRatioRecPrimPt->Reset();
3347 //     hRatioRecPrimPt->Divide(histPtRec,histPtRecPrim,1,1,"B");
3348
3349 //     TH1F* hRatioRecPrimZ = (TH1F*) histZRec->Clone(histNameZ);
3350 //     hRatioRecPrimZ->Reset();
3351 //     hRatioRecPrimZ->Divide(histZRec,histZRecPrim,1,1,"B");
3352
3353 //     TH1F* hRatioRecPrimXi = (TH1F*) histXiRec->Clone(histNameXi);
3354 //     hRatioRecPrimXi->Reset();
3355 //     hRatioRecPrimXi->Divide(histXiRec,histXiRecPrim,1,1,"B");
3356
3357
3358 //     new(fh1FFRatioRecPrimPt[i]) TH1F(*hRatioRecPrimPt);
3359 //     new(fh1FFRatioRecPrimZ[i])  TH1F(*hRatioRecPrimZ);
3360 //     new(fh1FFRatioRecPrimXi[i]) TH1F(*hRatioRecPrimXi);
3361 //   }
3362 // }
3363
3364 // __________________________________________________________________________________
3365 void AliFragmentationFunctionCorrections::ProjectJetResponseMatrices(TString strOutfile)
3366 {
3367
3368   // project response matrices on both axes:
3369   // FF for rec primaries, in terms of generated and reconstructed momentum
3370   // write FF and ratios to outFile
3371
3372   Printf("%s:%d -- project response matrices, write to %s",(char*)__FILE__,__LINE__,strOutfile.Data());
3373
3374   TH1F* hFFPtRec[fNJetPtSlices];
3375   TH1F* hFFZRec[fNJetPtSlices];
3376   TH1F* hFFXiRec[fNJetPtSlices];
3377
3378   TH1F* hFFPtGen[fNJetPtSlices];
3379   TH1F* hFFZGen[fNJetPtSlices];
3380   TH1F* hFFXiGen[fNJetPtSlices];
3381
3382   TH1F* hRatioPt[fNJetPtSlices];
3383   TH1F* hRatioZ[fNJetPtSlices];
3384   TH1F* hRatioXi[fNJetPtSlices];
3385
3386
3387   Int_t axisGenPt = 1;
3388   Int_t axisRecPt = 0;
3389
3390   for(Int_t i=0; i<fNJetPtSlices; i++){
3391
3392     Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
3393     Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
3394
3395     TString strNameRecPt(Form("h1FFTrackPtRecPrim_recPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3396     TString strNameRecZ(Form("h1FFZRecPrim_recPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3397     TString strNameRecXi(Form("h1FFXiRecPrim_recPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3398
3399     TString strNameGenPt(Form("h1FFTrackPtRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3400     TString strNameGenZ(Form("h1FFZRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3401     TString strNameGenXi(Form("h1FFXiRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3402
3403     TString strNameRatioPt(Form("h1RatioTrackPtRecPrim_%02d_%02d",jetPtLoLim,jetPtUpLim));
3404     TString strNameRatioZ(Form("h1RatioZRecPrim_%02d_%02d",jetPtLoLim,jetPtUpLim));
3405     TString strNameRatioXi(Form("h1RatioXiRecPrim_%02d_%02d",jetPtLoLim,jetPtUpLim));
3406
3407
3408     hFFPtRec[i] = (TH1F*) fhnResponsePt[i]->Projection(axisRecPt);// note convention: yDim,xDim
3409     hFFZRec[i]  = (TH1F*) fhnResponseZ[i]->Projection(axisRecPt);// note convention: yDim,xDim
3410     hFFXiRec[i] = (TH1F*) fhnResponseXi[i]->Projection(axisRecPt);// note convention: yDim,xDim
3411
3412     hFFPtRec[i]->SetNameTitle(strNameRecPt,"");
3413     hFFZRec[i]->SetNameTitle(strNameRecZ,"");
3414     hFFXiRec[i]->SetNameTitle(strNameRecXi,"");
3415
3416
3417     hFFPtGen[i] = (TH1F*) fhnResponsePt[i]->Projection(axisGenPt);// note convention: yDim,xDim
3418     hFFZGen[i]  = (TH1F*) fhnResponseZ[i]->Projection(axisGenPt);// note convention: yDim,xDim
3419     hFFXiGen[i] = (TH1F*) fhnResponseXi[i]->Projection(axisGenPt);// note convention: yDim,xDim
3420
3421     hFFPtGen[i]->SetNameTitle(strNameGenPt,"");
3422     hFFZGen[i]->SetNameTitle(strNameGenZ,"");
3423     hFFXiGen[i]->SetNameTitle(strNameGenXi,"");
3424
3425
3426     if(fNHistoBinsPt[i]) hFFPtRec[i] = (TH1F*) hFFPtRec[i]->Rebin(fNHistoBinsPt[i],"",fHistoBinsPt[i]->GetArray());
3427     if(fNHistoBinsZ[i])  hFFZRec[i]   = (TH1F*) hFFZRec[i]->Rebin(fNHistoBinsZ[i],"",fHistoBinsZ[i]->GetArray());
3428     if(fNHistoBinsXi[i]) hFFXiRec[i]  = (TH1F*) hFFXiRec[i]->Rebin(fNHistoBinsXi[i],"",fHistoBinsXi[i]->GetArray());
3429
3430     if(fNHistoBinsPt[i]) hFFPtGen[i] = (TH1F*) hFFPtGen[i]->Rebin(fNHistoBinsPt[i],"",fHistoBinsPt[i]->GetArray());
3431     if(fNHistoBinsZ[i])  hFFZGen[i]  = (TH1F*) hFFZGen[i]->Rebin(fNHistoBinsZ[i],"",fHistoBinsZ[i]->GetArray());
3432     if(fNHistoBinsXi[i]) hFFXiGen[i] = (TH1F*) hFFXiGen[i]->Rebin(fNHistoBinsXi[i],"",fHistoBinsXi[i]->GetArray());
3433
3434     NormalizeTH1(hFFPtGen[i],fNJets->At(i));
3435     NormalizeTH1(hFFZGen[i],fNJets->At(i));
3436     NormalizeTH1(hFFXiGen[i],fNJets->At(i));
3437
3438     NormalizeTH1(hFFPtRec[i],fNJets->At(i));
3439     NormalizeTH1(hFFZRec[i],fNJets->At(i));
3440     NormalizeTH1(hFFXiRec[i],fNJets->At(i));
3441
3442
3443     hRatioPt[i] = (TH1F*) hFFPtRec[i]->Clone(strNameRatioPt);
3444     hRatioPt[i]->Reset();
3445     hRatioPt[i]->Divide(hFFPtRec[i],hFFPtGen[i],1,1,"B");
3446
3447     hRatioZ[i] = (TH1F*) hFFZRec[i]->Clone(strNameRatioZ);
3448     hRatioZ[i]->Reset();
3449     hRatioZ[i]->Divide(hFFZRec[i],hFFZGen[i],1,1,"B");
3450
3451     hRatioXi[i] = (TH1F*) hFFXiRec[i]->Clone(strNameRatioXi);
3452     hRatioXi[i]->Reset();
3453     hRatioXi[i]->Divide(hFFXiRec[i],hFFXiGen[i],1,1,"B");
3454   }
3455
3456
3457
3458   // write
3459
3460   TFile out(strOutfile,"RECREATE");
3461
3462   if(!out.IsOpen()){
3463     Printf("%s:%d -- error opening response matrix output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
3464     return;
3465   }
3466
3467   for(Int_t i=0; i<fNJetPtSlices; i++){
3468
3469     hFFPtRec[i]->Write();
3470     hFFZRec[i]->Write();
3471     hFFXiRec[i]->Write();
3472
3473     hFFPtGen[i]->Write();
3474     hFFZGen[i]->Write();
3475     hFFXiGen[i]->Write();
3476
3477     hRatioPt[i]->Write();
3478     hRatioZ[i]->Write();
3479     hRatioXi[i]->Write();
3480   }
3481
3482   out.Close();
3483 }
3484
3485 // ____________________________________________________________________________________________________________________________
3486 void AliFragmentationFunctionCorrections::ProjectSingleResponseMatrix(TString strOutfile, Bool_t updateOutfile, TString strOutDir )
3487 {
3488   // project response matrix on both axes:
3489   // pt spec for rec primaries, in terms of generated and reconstructed momentum
3490   // write spec and ratios to outFile
3491
3492   Printf("%s:%d -- project single pt response matrix, write to %s",(char*)__FILE__,__LINE__,strOutfile.Data());
3493
3494   TH1F* hSpecPtRec;
3495   TH1F* hSpecPtGen;
3496   TH1F* hRatioPt;
3497
3498   Int_t axisGenPt = 1;
3499   Int_t axisRecPt = 0;
3500
3501   TString strNameRecPt   = "h1SpecTrackPtRecPrim_recPt";
3502   TString strNameGenPt   = "h1SpecTrackPtRecPrim_genPt";
3503   TString strNameRatioPt = "h1RatioTrackPtRecPrim";
3504
3505   hSpecPtRec = (TH1F*) fhnResponseSinglePt->Projection(axisRecPt);// note convention: yDim,xDim
3506   hSpecPtRec->SetNameTitle(strNameRecPt,"");
3507
3508   hSpecPtGen = (TH1F*) fhnResponseSinglePt->Projection(axisGenPt);// note convention: yDim,xDim
3509   hSpecPtGen->SetNameTitle(strNameGenPt,"");
3510
3511   hRatioPt = (TH1F*) hSpecPtRec->Clone(strNameRatioPt);
3512   hRatioPt->Reset();
3513   hRatioPt->Divide(hSpecPtRec,hSpecPtGen,1,1,"B");
3514
3515   TString outfileOption = "RECREATE";
3516   if(updateOutfile) outfileOption = "UPDATE";
3517
3518   TFile out(strOutfile,outfileOption);
3519
3520   if(!out.IsOpen()){
3521     Printf("%s:%d -- error opening reponse matrix projections output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
3522     return;
3523   }
3524
3525
3526   if(strOutDir && strOutDir.Length()){
3527
3528     TDirectory* dir;
3529     if((dir = ((TDirectory*) gDirectory->Get(strOutDir)))) dir->cd();
3530     else{
3531       dir = out.mkdir(strOutDir);
3532       dir->cd();
3533     }
3534   }
3535
3536   hSpecPtRec->Write();
3537   hSpecPtGen->Write();
3538   hRatioPt->Write();
3539
3540   out.Close();
3541 }
3542
3543
3544 //__________________________________________________________________________________________________________________________________________________________________
3545 void AliFragmentationFunctionCorrections::RebinHisto(const Int_t jetPtSlice, const Int_t nBinsLimits, Double_t* binsLimits, Double_t* binsWidth, const Int_t type)
3546 {
3547   // rebin histo, rescale bins according to new width
3548   // only correct for input histos with equal bin size
3549
3550   // args: jetPtSlice, type, use current corr level
3551
3552   // function takes array of limits and widths (e.g. 1 GeV bins up to 10 GeV, 2 GeV width up to 20 GeV, ...)
3553   // array size of binsLimits: nBinsLimits
3554   // array size of binsWidth: nBinsLimits-1
3555   // binsLimits have to be in increasing order
3556   // if binning undefined for any slice, original binning will be kept
3557
3558   if(!fNJetPtSlices){
3559     Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
3560     return;
3561   }
3562
3563   if(jetPtSlice>=fNJetPtSlices){
3564     Printf("%s:%d -- jetPtSlice %d exceeds max",(char*)__FILE__,__LINE__,jetPtSlice);
3565     return;
3566   }
3567
3568
3569   Double_t binLimitMin = binsLimits[0];
3570   Double_t binLimitMax = binsLimits[nBinsLimits-1];
3571
3572   Double_t binLimit = binLimitMin; // start value
3573
3574   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 ...
3575   TArrayD binsArray(sizeUpperLim);
3576   Int_t nBins = 0;
3577   binsArray.SetAt(binLimitMin,nBins++);
3578
3579   while(binLimit<binLimitMax && nBins<sizeUpperLim){
3580
3581     Int_t currentSlice = -1;
3582     for(Int_t i=0; i<nBinsLimits; i++){
3583       if(binLimit >= 0.999*binsLimits[i]) currentSlice = i; // 0.999 numerical saftey factor
3584     }
3585
3586     Double_t currentBinWidth = binsWidth[currentSlice];
3587     binLimit += currentBinWidth;
3588
3589     binsArray.SetAt(binLimit,nBins++);
3590   }
3591
3592
3593   TH1F* hist = 0;
3594   if(type == kFlagPt)            hist = fCorrFF[fNCorrectionLevels-1]->GetTrackPt(jetPtSlice);
3595   else if(type == kFlagZ)        hist = fCorrFF[fNCorrectionLevels-1]->GetZ(jetPtSlice);
3596   else if(type == kFlagXi)       hist = fCorrFF[fNCorrectionLevels-1]->GetXi(jetPtSlice);
3597   else if(type == kFlagSinglePt) hist = fCorrSinglePt[fNCorrectionLevelsSinglePt-1]->GetTrackPt(0);
3598   else{
3599     Printf("%s%d unknown type",(char*)__FILE__,__LINE__);
3600     return;
3601   }
3602
3603   Double_t binWidthNoRebin = hist->GetBinWidth(1);
3604
3605   Double_t* bins = binsArray.GetArray();
3606
3607   hist = (TH1F*) hist->Rebin(nBins-1,"",bins);
3608
3609   for(Int_t bin=0; bin <= hist->GetNbinsX(); bin++){
3610
3611     Double_t binWidthRebin = hist->GetBinWidth(bin);
3612     Double_t scaleF = binWidthNoRebin / binWidthRebin;
3613
3614     Double_t binCont  = hist->GetBinContent(bin);
3615     Double_t binErr   = hist->GetBinError(bin);
3616
3617     binCont *= scaleF;
3618     binErr  *= scaleF;
3619
3620     hist->SetBinContent(bin,binCont);
3621     hist->SetBinError(bin,binErr);
3622   }
3623
3624
3625
3626   TH1F* temp = new TH1F(*hist);
3627
3628   if(type == kFlagPt)       fCorrFF[fNCorrectionLevels-1]->ReplaceCorrHistos(jetPtSlice,temp,0,0);
3629   if(type == kFlagZ)        fCorrFF[fNCorrectionLevels-1]->ReplaceCorrHistos(jetPtSlice,0,temp,0);
3630   if(type == kFlagXi)       fCorrFF[fNCorrectionLevels-1]->ReplaceCorrHistos(jetPtSlice,0,0,temp);
3631   if(type == kFlagSinglePt) fCorrSinglePt[fNCorrectionLevelsSinglePt-1]->ReplaceCorrHistos(0,temp,0,0);
3632
3633
3634   delete temp;
3635 }
3636 //__________________________________________________________________________________________________________________________________________________________________
3637 void AliFragmentationFunctionCorrections::WriteJetSpecResponse(TString strInfile, TString strdir, TString strlist/*, TString strOutfile*/)
3638 {
3639
3640  if(fDebug>0) Printf("%s:%d -- read jet spectrum response matrix from file %s ",(char*)__FILE__,__LINE__,strInfile.Data());
3641
3642   if(strdir && strdir.Length()) gDirectory->cd(strdir);
3643
3644   TList* list = 0;
3645
3646   if(strlist && strlist.Length()){
3647
3648     if(!(list = (TList*) gDirectory->Get(strlist))){
3649       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
3650       return;
3651     }
3652   }
3653
3654   THnSparse* hn6ResponseJetPt  = (THnSparse*) list->FindObject("fhnCorrelation");
3655
3656   Int_t axis6RecJetPt = 0;
3657   Int_t axis6GenJetPt = 3;
3658
3659   hn6ResponseJetPt->GetAxis(axis6RecJetPt)->SetTitle("rec jet p_{T} (GeV/c)");
3660   hn6ResponseJetPt->GetAxis(axis6GenJetPt)->SetTitle("gen jet p_{T} (GeV/c)");
3661
3662   Int_t nBinsRecPt    = hn6ResponseJetPt->GetAxis(axis6RecJetPt)->GetNbins();
3663   Double_t loLimRecPt = hn6ResponseJetPt->GetAxis(axis6RecJetPt)->GetBinLowEdge(1);
3664   Double_t upLimRecPt = hn6ResponseJetPt->GetAxis(axis6RecJetPt)->GetBinUpEdge(nBinsRecPt);
3665
3666   Int_t nBinsGenPt    = hn6ResponseJetPt->GetAxis(axis6GenJetPt)->GetNbins();
3667   Double_t loLimGenPt = hn6ResponseJetPt->GetAxis(axis6GenJetPt)->GetBinLowEdge(1);
3668   Double_t upLimGenPt = hn6ResponseJetPt->GetAxis(axis6GenJetPt)->GetBinUpEdge(nBinsGenPt);
3669
3670   Int_t nBinsTrackPt = 200;
3671   Int_t loLimTrackPt = 0;
3672   Int_t upLimTrackPt = 200;
3673
3674
3675   Int_t    nBinsResponse[4]  = {nBinsRecPt,nBinsTrackPt,nBinsGenPt,nBinsTrackPt};
3676   Double_t binMinResponse[4] = {loLimRecPt,loLimTrackPt,loLimGenPt,loLimTrackPt};
3677   Double_t binMaxResponse[4] = {upLimRecPt,upLimTrackPt,upLimGenPt,upLimTrackPt};
3678
3679   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)"};
3680
3681   THnSparseD* hn4ResponseTrackPtJetPt =  new THnSparseD("hn4ResponseTrackPtJetPt","",4,nBinsResponse,binMinResponse,binMaxResponse);
3682
3683   for(Int_t i=0; i<4; i++){
3684     hn4ResponseTrackPtJetPt->GetAxis(i)->SetTitle(labelsResponseSinglePt[i]);
3685   }
3686
3687
3688   // fill
3689
3690
3691 }
3692
3693 //_____________________________________________________________________________________________________________________________________
3694 void AliFragmentationFunctionCorrections::ReadSingleTrackEfficiency(TString strfile, TString strdir, TString strlist, TString strname)
3695 {
3696
3697   ReadSingleTrackCorrection(strfile,strdir,strlist,strname,kFlagEfficiency);
3698
3699 }
3700
3701 //_____________________________________________________________________________________________________________________________________
3702 void AliFragmentationFunctionCorrections::ReadSingleTrackResponse(TString strfile, TString strdir, TString strlist, TString strname)
3703 {
3704
3705   ReadSingleTrackCorrection(strfile,strdir,strlist,strname,kFlagResponse);
3706
3707 }
3708
3709 //_____________________________________________________________________________________________________________________________________
3710 void AliFragmentationFunctionCorrections::ReadSingleTrackSecCorr(TString strfile, TString strdir, TString strlist, TString strname)
3711 {
3712
3713   ReadSingleTrackCorrection(strfile,strdir,strlist,strname,kFlagSecondaries);
3714
3715 }
3716
3717 //______________________________________________________________________________________________________________________________________________________
3718 void AliFragmentationFunctionCorrections::ReadSingleTrackCorrection(TString strfile, TString strdir, TString strlist, TString strname, const Int_t type)
3719 {
3720   // read single track correction (pt) from file
3721   // type: efficiency / response / secondaries correction
3722
3723   if(!((type == kFlagEfficiency)  || (type == kFlagResponse) || (type == kFlagSecondaries))){
3724     Printf("%s:%d -- no such correction ",(char*)__FILE__,__LINE__);
3725     return;
3726   }
3727
3728   TFile f(strfile,"READ");
3729
3730   if(!f.IsOpen()){
3731     Printf("%s:%d -- error opening file %s", (char*)__FILE__,__LINE__,strfile.Data());
3732     return;
3733   }
3734
3735   if(fDebug>0 && type==kFlagEfficiency)  Printf("%s:%d -- read single track corr from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
3736   if(fDebug>0 && type==kFlagResponse)    Printf("%s:%d -- read single track response from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
3737   if(fDebug>0 && type==kFlagSecondaries) Printf("%s:%d -- read single track secondaries corr from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
3738
3739   if(strdir && strdir.Length()) gDirectory->cd(strdir);
3740
3741   TList* list = 0;
3742
3743   if(strlist && strlist.Length()){
3744
3745     if(!(list = (TList*) gDirectory->Get(strlist))){
3746       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
3747       return;
3748     }
3749   }
3750
3751   TH1F*      h1CorrHist = 0; // common TObject pointer not possible, need SetDirectory() later
3752   THnSparse* hnCorrHist = 0;
3753
3754   if(type == kFlagEfficiency || type == kFlagSecondaries){
3755
3756     if(list) h1CorrHist = (TH1F*) list->FindObject(strname);
3757     else     h1CorrHist = (TH1F*) gDirectory->Get(strname);
3758
3759     if(!h1CorrHist){
3760       Printf("%s:%d -- error retrieving histo %s", (char*)__FILE__,__LINE__,strname.Data());
3761       return;
3762     }
3763
3764   }
3765   else if(type == kFlagResponse){
3766
3767     if(list) hnCorrHist = (THnSparse*) list->FindObject(strname);
3768     else     hnCorrHist = (THnSparse*) gDirectory->Get(strname);
3769
3770     if(!hnCorrHist){
3771       Printf("%s:%d -- error retrieving histo %s", (char*)__FILE__,__LINE__,strname.Data());
3772       return;
3773     }
3774
3775   }
3776
3777   if(h1CorrHist) h1CorrHist->SetDirectory(0);
3778   //if(hnCorrHist) hnCorrHist->SetDirectory(0);
3779
3780   f.Close();
3781
3782   if(type == kFlagEfficiency)       fh1EffSinglePt      = h1CorrHist;
3783   else if(type == kFlagResponse)    fhnResponseSinglePt = hnCorrHist;
3784   else if(type == kFlagSecondaries) fh1SecCorrSinglePt  = h1CorrHist;
3785
3786 }
3787
3788 //________________________________________________________________________________________________________________
3789 void AliFragmentationFunctionCorrections::ReadRawPtSpec(TString strInfile, TString strID)
3790 {
3791   // read track pt spec from task ouput - standard dir/list
3792
3793   TString strdir  = "PWG4_FragmentationFunction_" + strID;
3794   TString strlist = "fracfunc_" + strID;
3795
3796   ReadRawPtSpec(strInfile,strdir,strlist);
3797 }
3798
3799 //_______________________________________________________________________________________________________
3800 void AliFragmentationFunctionCorrections::ReadRawPtSpec(TString strfile, TString strdir, TString strlist)
3801 {
3802   // get raw pt spectra from file
3803
3804
3805   // book histos
3806   fNCorrectionLevelsSinglePt = 0;
3807   fCorrSinglePt = new AliFragFuncCorrHistos*[fgMaxNCorrectionLevels];
3808   AddCorrectionLevelSinglePt(); // first 'correction' level = raw spectrum
3809
3810   // get raw pt spec from input file, normalize
3811
3812   TFile f(strfile,"READ");
3813
3814   if(!f.IsOpen()){
3815     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
3816     return;
3817   }
3818
3819   if(fDebug>0) Printf("%s:%d -- read raw spectra from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
3820
3821   gDirectory->cd(strdir);
3822
3823   TList* list = 0;
3824
3825   if(!(list = (TList*) gDirectory->Get(strlist))){
3826     Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
3827     return;
3828   }
3829
3830   TString hnameTrackPt("fh1TrackQAPtRecCuts");
3831   TString hnameEvtSel("fh1EvtSelection");
3832
3833   TH1F* fh1TrackPt = (TH1F*) list->FindObject(hnameTrackPt);
3834   TH1F* fh1EvtSel  = (TH1F*) list->FindObject(hnameEvtSel);
3835
3836   if(!fh1TrackPt){ Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameTrackPt.Data()); return; }
3837   if(!fh1EvtSel) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameEvtSel.Data()); return; }
3838
3839
3840   //Float_t nEvents = fh1EvtSel->GetBinContent(fh1EvtSel->FindBin(0));
3841
3842   // evts after physics selection
3843   Float_t nEvents = fh1EvtSel->GetEntries()
3844     - fh1EvtSel->GetBinContent(fh1EvtSel->FindBin(1))   // evts rejected by trigger selection ( = PhysicsSelection
3845     - fh1EvtSel->GetBinContent(fh1EvtSel->FindBin(2));   // evts with wrong centrality class
3846
3847
3848   fh1TrackPt->SetDirectory(0);
3849
3850   f.Close();
3851
3852
3853   NormalizeTH1(fh1TrackPt,nEvents);
3854
3855   // raw FF = corr level 0
3856   fCorrSinglePt[0]->AddCorrHistos(0,fh1TrackPt);
3857 }
3858
3859
3860 //_______________________________________________________________________________________________________
3861 void AliFragmentationFunctionCorrections::ReadRawPtSpecQATask(TString strfile, TString strdir, TString strlist)
3862 {
3863   // get raw pt spectra from file
3864   // for output from Martas QA task
3865
3866
3867   // book histos
3868   fNCorrectionLevelsSinglePt = 0;
3869   fCorrSinglePt = new AliFragFuncCorrHistos*[fgMaxNCorrectionLevels];
3870   AddCorrectionLevelSinglePt(); // first 'correction' level = raw spectrum
3871
3872   // get raw pt spec from input file, normalize
3873
3874   TFile f(strfile,"READ");
3875
3876   if(!f.IsOpen()){
3877     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
3878     return;
3879   }
3880
3881   if(fDebug>0) Printf("%s:%d -- read raw pt spec from QA task output file %s ",(char*)__FILE__,__LINE__,strfile.Data());
3882
3883   gDirectory->cd(strdir);
3884
3885   TList* list = 0;
3886
3887   if(!(list = (TList*) gDirectory->Get(strlist))){
3888     Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
3889     return;
3890   }
3891
3892   TString hnameTrackPt("fPtSel");
3893   TString hnameEvtSel("fNEventAll");
3894
3895   TH1F* fh1TrackPt = (TH1F*) list->FindObject(hnameTrackPt);
3896   TH1F* fh1EvtSel  = (TH1F*) list->FindObject(hnameEvtSel);
3897
3898   if(!fh1TrackPt){ Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameTrackPt.Data()); return; }
3899   if(!fh1EvtSel) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameEvtSel.Data()); return; }
3900
3901
3902   // evts after physics selection
3903   Float_t nEvents = fh1EvtSel->GetEntries();
3904
3905   fh1TrackPt->SetDirectory(0);
3906
3907   f.Close();
3908
3909
3910   NormalizeTH1(fh1TrackPt,nEvents);
3911
3912   // raw FF = corr level 0
3913   fCorrSinglePt[0]->AddCorrHistos(0,fh1TrackPt);
3914 }
3915
3916 // ________________________________________________________
3917 void AliFragmentationFunctionCorrections::EffCorrSinglePt()
3918 {
3919   // apply efficiency correction to inclusive track pt spec
3920
3921   AddCorrectionLevelSinglePt("eff");
3922
3923   TH1F* histPt = fCorrSinglePt[fNCorrectionLevelsSinglePt-2]->GetTrackPt(0);
3924
3925   if(histPt->GetNbinsX() != fh1EffSinglePt->GetNbinsX()) Printf("%s:%d: inconsistency pt spec and eff corr bins - rebin effCorr ...", (char*)__FILE__,__LINE__);
3926
3927   TString histNamePt = histPt->GetName();
3928   TH1F* hTrackPtEffCorr = (TH1F*) histPt->Clone(histNamePt);
3929   hTrackPtEffCorr->Divide(histPt,fh1EffSinglePt,1,1,"");
3930
3931   fCorrSinglePt[fNCorrectionLevelsSinglePt-1]->AddCorrHistos(0,hTrackPtEffCorr);
3932 }
3933
3934 //___________________________________________________________________________________________________________________________
3935 void AliFragmentationFunctionCorrections::UnfoldSinglePt(const Int_t nIter, const Bool_t useCorrelatedErrors)
3936 {
3937   // unfolde inclusive dN/dpt spectra
3938
3939   AddCorrectionLevelSinglePt("unfold");
3940
3941   TH1F* hist = fCorrSinglePt[fNCorrectionLevelsSinglePt-2]->GetTrackPt(0); // level -2: before unfolding, level -1: unfolded
3942   THnSparse* hnResponse = fhnResponseSinglePt;
3943
3944   TString histNameTHn = hist->GetName();
3945   if(histNameTHn.Contains("TH1")) histNameTHn.ReplaceAll("TH1","THn");
3946   if(histNameTHn.Contains("fPt")) histNameTHn.ReplaceAll("fPt","fhnPt");
3947
3948
3949   TString histNameBackFolded = hist->GetName();
3950   histNameBackFolded.Append("_backfold");
3951
3952   TString histNameRatioFolded = hist->GetName();
3953   if(histNameRatioFolded.Contains("fh1")) histNameRatioFolded.ReplaceAll("fh1","hRatio");
3954   if(histNameRatioFolded.Contains("fPt")) histNameRatioFolded.ReplaceAll("fPt","hRatioPt");
3955   histNameRatioFolded.Append("_unfold");
3956
3957   TString histNameRatioBackFolded = hist->GetName();
3958   if(histNameRatioBackFolded.Contains("fh1")) histNameRatioBackFolded.ReplaceAll("fh1","hRatio");
3959   if(histNameRatioBackFolded.Contains("fPt")) histNameRatioBackFolded.ReplaceAll("fPt","hRatioPt");
3960   histNameRatioBackFolded.Append("_backfold");
3961
3962   THnSparse* hnHist           = TH1toSparse(hist,histNameTHn,hist->GetTitle());
3963   THnSparse* hnFlatEfficiency = TH1toSparse(hist,"fhnEfficiency","eff",kTRUE); // could optionally also use real eff
3964
3965   THnSparse* hnUnfolded
3966     = Unfold(hnHist,hnResponse,hnFlatEfficiency,nIter,useCorrelatedErrors);
3967
3968   TH1F* hUnfolded = (TH1F*) hnUnfolded->Projection(0);
3969   hUnfolded->SetNameTitle(hist->GetName(),hist->GetTitle());
3970
3971   fCorrSinglePt[fNCorrectionLevelsSinglePt-1]->AddCorrHistos(0,hUnfolded);
3972
3973   // backfolding: apply response matrix to unfolded spectrum
3974   TH1F* hBackFolded = ApplyResponse(hUnfolded,hnResponse);
3975   hBackFolded->SetNameTitle(histNameBackFolded,hist->GetTitle());
3976
3977   fh1SingleTrackPtBackFolded = hBackFolded;
3978
3979
3980   // ratio unfolded to original histo
3981   TH1F* hRatioUnfolded = (TH1F*) hUnfolded->Clone(histNameRatioFolded);
3982   hRatioUnfolded->Reset();
3983   hRatioUnfolded->Divide(hUnfolded,hist,1,1,"B");
3984
3985   fh1RatioSingleTrackPtFolded = hRatioUnfolded;
3986
3987
3988   // ratio backfolded to original histo
3989   TH1F* hRatioBackFolded = (TH1F*) hBackFolded->Clone(histNameRatioBackFolded);
3990   hRatioBackFolded->Reset();
3991   hRatioBackFolded->Divide(hBackFolded,hist,1,1,"B");
3992
3993   fh1RatioSingleTrackPtBackFolded = hRatioBackFolded;
3994
3995   delete hnHist;
3996   delete hnFlatEfficiency;
3997
3998 }
3999
4000 // ________________________________________________________
4001 void AliFragmentationFunctionCorrections::SecCorrSinglePt()
4002 {
4003   // apply efficiency correction to inclusive track pt spec
4004
4005   AddCorrectionLevelSinglePt("secCorr");
4006
4007   TH1F* histPt = fCorrSinglePt[fNCorrectionLevelsSinglePt-2]->GetTrackPt(0);
4008
4009   if(histPt->GetNbinsX() != fh1SecCorrSinglePt->GetNbinsX())
4010     Printf("%s:%d: inconsistency pt spec and secondaries corr bins - rebin effCorr ...", (char*)__FILE__,__LINE__);
4011
4012   TString histNamePt = histPt->GetName();
4013   TH1F* hTrackPtSecCorr = (TH1F*) histPt->Clone(histNamePt);
4014
4015   hTrackPtSecCorr->Multiply(histPt,fh1SecCorrSinglePt,1,1,"");
4016
4017   fCorrSinglePt[fNCorrectionLevelsSinglePt-1]->AddCorrHistos(0,hTrackPtSecCorr);
4018 }