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