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