PWGJE/AliFragmentationFunctionCorrections.cxx

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