PWGHF/hfe/AliHFEInclusiveSpectrum.cxx

   1
   2 /**************************************************************************
   3 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   4 *                                                                        *
   5 * Author: The ALICE Off-line Project.                                    *
   6 * Contributors are mentioned in the code where appropriate.              *
   7 *                                                                        *
   8 * Permission to use, copy, modify and distribute this software and its   *
   9 * documentation strictly for non-commercial purposes is hereby granted   *
  10 * without fee, provided that the above copyright notice appears in all   *
  11 * copies and that both the copyright notice and this permission notice   *
  12 * appear in the supporting documentation. The authors make no claims     *
  13 * about the suitability of this software for any purpose. It is          *
  14 * provided "as is" without express or implied warranty.                  *
  15 **************************************************************************/
  16 //
  17 // Class for spectrum correction
  18 // Subtraction of hadronic background, Unfolding of the data and
  19 // Renormalization done here
  20 // The following containers have to be set:
  21 //  - Correction framework container for real data
  22 //  - Correction framework container for MC (Efficiency Map)
  23 //  - Correction framework container for background coming from data
  24 //  - Correction framework container for background coming from MC
  25 //
  26 //  Author:
  27 //            Raphaelle Bailhache <R.Bailhache@gsi.de>
  28 //            Markus Fasel <M.Fasel@gsi.de>
  29 //
  30
  31 #include <TArrayD.h>
  32 #include <TH1.h>
  33 #include <TList.h>
  34 #include <TObjArray.h>
  35 #include <TROOT.h>
  36 #include <TCanvas.h>
  37 #include <TLegend.h>
  38 #include <TStyle.h>
  39 #include <TMath.h>
  40 #include <TAxis.h>
  41 #include <TGraphErrors.h>
  42 #include <TFile.h>
  43 #include <TPad.h>
  44 #include <TH2D.h>
  45 #include <TF1.h>
  46
  47 #include "AliPID.h"
  48 #include "AliCFContainer.h"
  49 #include "AliCFDataGrid.h"
  50 #include "AliCFEffGrid.h"
  51 #include "AliCFGridSparse.h"
  52 #include "AliCFUnfolding.h"
  53 #include "AliLog.h"
  54
  55 #include "AliHFEInclusiveSpectrum.h"
  56 #include "AliHFEInclusiveSpectrumQA.h"
  57 #include "AliHFEcuts.h"
  58 #include "AliHFEcontainer.h"
  59 #include "AliHFEtools.h"
  60
  61 ClassImp(AliHFEInclusiveSpectrum)
  62
  63 //____________________________________________________________
  64 AliHFEInclusiveSpectrum::AliHFEInclusiveSpectrum(const char *name):
  65   AliHFECorrectSpectrumBase(name),
  66   fQA(NULL)
  67 {
  68   //
  69   // Default constructor
  70   //
  71
  72   fQA = new AliHFEInclusiveSpectrumQA("AliHFEInclusiveSpectrumQA");
  73
  74 }
  75 //____________________________________________________________
  76 AliHFEInclusiveSpectrum::AliHFEInclusiveSpectrum(const AliHFEInclusiveSpectrum &ref):
  77   AliHFECorrectSpectrumBase(ref),
  78   fQA(ref.fQA)
  79 {
  80   //
  81   // Copy constructor
  82   //
  83   ref.Copy(*this);
  84
  85 }
  86 //____________________________________________________________
  87 AliHFEInclusiveSpectrum &AliHFEInclusiveSpectrum::operator=(const AliHFEInclusiveSpectrum &ref){
  88   //
  89   // Assignment operator
  90   //
  91   if(this == &ref)
  92     ref.Copy(*this);
  93   return *this;
  94 }
  95 //____________________________________________________________
  96 void AliHFEInclusiveSpectrum::Copy(TObject &o) const {
  97   //
  98   // Copy into object o
  99   //
 100   AliHFEInclusiveSpectrum &target = dynamic_cast<AliHFEInclusiveSpectrum &>(o);
 101   target.fQA = fQA;
 102
 103
 104 }
 105 //____________________________________________________________
 106 AliHFEInclusiveSpectrum::~AliHFEInclusiveSpectrum(){
 107   //
 108   // Destructor
 109   //
 110   if(fQA) delete fQA;
 111
 112 }
 113 //____________________________________________________________
 114 Bool_t AliHFEInclusiveSpectrum::Init(const AliHFEcontainer *datahfecontainer, const AliHFEcontainer *mchfecontainer, const AliHFEcontainer *v0hfecontainer, const AliHFEcontainer */*bghfecontainer*/){
 115   //
 116   // Init what we need for the correction:
 117   //
 118   // Raw spectrum, hadron contamination
 119   // MC efficiency maps, correlation matrix
 120   // V0 efficiency if wanted
 121   //
 122   // This for a given dimension.
 123   //
 124   //
 125
 126   Int_t kNdim = 3;
 127
 128   Int_t dims[kNdim];
 129   switch(fNbDimensions){
 130   case 1:   dims[0] = 0;
 131     break;
 132   case 2:   for(Int_t i = 0; i < 2; i++) dims[i] = i;
 133     break;
 134   case 3:   for(Int_t i = 0; i < 3; i++) dims[i] = i;
 135     break;
 136   default:
 137     AliError("Container with this number of dimensions not foreseen (yet)");
 138     return kFALSE;
 139   };
 140
 141   //
 142   // Data container: raw spectrum + hadron contamination
 143   //
 144   AliCFContainer *datacontainer = datahfecontainer->GetCFContainer("recTrackContReco");
 145   AliCFContainer *contaminationcontainer = datahfecontainer->GetCFContainer("hadronicBackground");
 146   if((!datacontainer) || (!contaminationcontainer)) return kFALSE;
 147   AliCFContainer *datacontainerD = GetSlicedContainer(datacontainer, fNbDimensions, dims, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
 148   AliCFContainer *contaminationcontainerD = GetSlicedContainer(contaminationcontainer, fNbDimensions, dims, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
 149   if((!datacontainerD) || (!contaminationcontainerD)) return kFALSE;
 150   SetContainer(datacontainerD,AliHFECorrectSpectrumBase::kDataContainer);
 151   SetContainer(contaminationcontainerD,AliHFECorrectSpectrumBase::kBackgroundData);
 152   // QA
 153   Int_t dimqa = datacontainer->GetNVar();
 154   Int_t dimsqa[dimqa];
 155   for(Int_t i = 0; i < dimqa; i++) dimsqa[i] = i;
 156   AliCFContainer *datacontainerDQA = GetSlicedContainer(datacontainer, dimqa, dimsqa, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
 157   fQA->AddResultAt(datacontainerDQA,AliHFEInclusiveSpectrumQA::kDataProjection);
 158
 159   //
 160   // MC container: ESD/MC efficiency maps + MC/MC efficiency maps
 161   //
 162   AliCFContainer *mccontaineresd = 0x0;
 163   AliCFContainer *mccontainermc = 0x0;
 164   mccontaineresd = mchfecontainer->MakeMergedCFContainer("sumesd","sumesd","MCTrackCont:recTrackContReco");
 165   mccontainermc = mchfecontainer->MakeMergedCFContainer("summc","summc","MCTrackCont:recTrackContMC");
 166   if((!mccontaineresd) || (!mccontainermc)) return kFALSE;
 167   Int_t source = -1;
 168   AliCFContainer *mccontaineresdD = GetSlicedContainer(mccontaineresd, fNbDimensions, dims, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
 169   AliCFContainer *mccontainermcD = GetSlicedContainer(mccontainermc, fNbDimensions, dims, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
 170   if((!mccontaineresdD) || (!mccontainermcD)) return kFALSE;
 171   SetContainer(mccontainermcD,AliHFECorrectSpectrumBase::kMCContainerMC);
 172   SetContainer(mccontaineresdD,AliHFECorrectSpectrumBase::kMCContainerESD);
 173
 174   //
 175   // Correlation matrix
 176   //
 177   THnSparseF *mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepbeforePID");
 178   if(!mccorrelation) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID");
 179   if(!mccorrelation) return kFALSE;
 180   THnSparseF *mccorrelationD = GetSlicedCorrelation(mccorrelation, fNbDimensions, dims,fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
 181   if(!mccorrelationD) {
 182     printf("No correlation\n");
 183     return kFALSE;
 184   }
 185   SetCorrelation(mccorrelationD);
 186   // QA
 187   fQA->AddResultAt(mccorrelation,AliHFEInclusiveSpectrumQA::kCMProjection);
 188
 189   //
 190   // V0 container Electron, pt eta phi
 191   //
 192   if(v0hfecontainer) {
 193     AliCFContainer *containerV0 = v0hfecontainer->GetCFContainer("taggedTrackContainerReco");
 194     if(!containerV0) return kFALSE;
 195     AliCFContainer *containerV0Electron = GetSlicedContainer(containerV0, fNbDimensions, dims, AliPID::kElectron,fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
 196     if(!containerV0Electron) return kFALSE;
 197     SetContainer(containerV0Electron,AliHFECorrectSpectrumBase::kDataContainerV0);
 198   }
 199
 200   //
 201   fQA->DrawProjections();
 202
 203
 204   return kTRUE;
 205 }
 206 //____________________________________________________________
 207 Bool_t AliHFEInclusiveSpectrum::Correct(Bool_t subtractcontamination){
 208   //
 209   // Correct the spectrum for efficiency and unfolding
 210   // with both method and compare
 211   //
 212
 213   gStyle->SetPalette(1);
 214   gStyle->SetOptStat(1111);
 215   gStyle->SetPadBorderMode(0);
 216   gStyle->SetCanvasColor(10);
 217   gStyle->SetPadLeftMargin(0.13);
 218   gStyle->SetPadRightMargin(0.13);
 219
 220   printf("Steps are: stepdata %d, stepMC %d, steptrue %d, stepV0after %d, stepV0before %d\n",fStepData,fStepMC,fStepTrue,fStepAfterCutsV0,fStepBeforeCutsV0);
 221
 222   ///////////////////////////
 223   // Check initialization
 224   ///////////////////////////
 225
 226   if((!GetContainer(kDataContainer)) || (!GetContainer(kMCContainerMC)) || (!GetContainer(kMCContainerESD))){
 227     AliInfo("You have to init before");
 228     return kFALSE;
 229   }
 230
 231   if((fStepTrue < 0) && (fStepMC < 0) && (fStepData < 0)) {
 232     AliInfo("You have to set the steps before: SetMCTruthStep, SetMCEffStep, SetStepToCorrect");
 233     return kFALSE;
 234   }
 235
 236   SetStepGuessedUnfolding(AliHFEcuts::kStepRecKineITSTPC + AliHFEcuts::kNcutStepsMCTrack);
 237
 238   AliCFDataGrid *dataGridAfterFirstSteps = 0x0;
 239   //////////////////////////////////
 240   // Subtract hadron background
 241   /////////////////////////////////
 242   AliCFDataGrid *dataspectrumaftersubstraction = 0x0;
 243   if(subtractcontamination) {
 244     dataspectrumaftersubstraction = SubtractBackground();
 245     dataGridAfterFirstSteps = dataspectrumaftersubstraction;
 246   }
 247
 248   ////////////////////////////////////////////////
 249   // Correct for TPC efficiency from V0 if any
 250   ///////////////////////////////////////////////
 251   AliCFDataGrid *dataspectrumafterV0efficiencycorrection = 0x0;
 252   AliCFContainer *dataContainerV0 = GetContainer(kDataContainerV0);
 253   if(dataContainerV0){
 254     dataspectrumafterV0efficiencycorrection = CorrectV0Efficiency(dataspectrumaftersubstraction);
 255     dataGridAfterFirstSteps = dataspectrumafterV0efficiencycorrection;
 256   }
 257
 258   //////////////////////////////////////////////////////////////////////////////
 259   // Correct for efficiency parametrized (if TPC efficiency is parametrized)
 260   /////////////////////////////////////////////////////////////////////////////
 261   AliCFDataGrid *dataspectrumafterefficiencyparametrizedcorrection = 0x0;
 262   if(fEfficiencyFunction){
 263     dataspectrumafterefficiencyparametrizedcorrection = CorrectParametrizedEfficiency(dataGridAfterFirstSteps);
 264     dataGridAfterFirstSteps = dataspectrumafterefficiencyparametrizedcorrection;
 265   }
 266
 267   ///////////////
 268   // Unfold
 269   //////////////
 270   THnSparse *correctedspectrum = Unfold(dataGridAfterFirstSteps);
 271   if(!correctedspectrum){
 272     AliError("No corrected spectrum\n");
 273     return kFALSE;
 274   }
 275
 276   /////////////////////
 277   // Simply correct
 278   ////////////////////
 279   AliCFDataGrid *alltogetherCorrection = CorrectForEfficiency(dataGridAfterFirstSteps);
 280
 281   // QA final results
 282   TGraphErrors* correctedspectrumD = Normalize(correctedspectrum);
 283   TGraphErrors* alltogetherspectrumD = Normalize(alltogetherCorrection);
 284   fQA->AddResultAt(correctedspectrumD,AliHFEInclusiveSpectrumQA::kFinalResultUnfolded);
 285   fQA->AddResultAt(alltogetherspectrumD,AliHFEInclusiveSpectrumQA::kFinalResultDirectEfficiency);
 286   fQA->DrawResult();
 287
 288   return kTRUE;
 289 }
 290
 291 //____________________________________________________________
 292 AliCFDataGrid* AliHFEInclusiveSpectrum::SubtractBackground(){
 293   //
 294   // Apply background subtraction
 295   //
 296
 297   // Raw spectrum
 298   AliCFContainer *dataContainer = GetContainer(kDataContainer);
 299   if(!dataContainer){
 300     AliError("Data Container not available");
 301     return NULL;
 302   }
 303   printf("Step data: %d\n",fStepData);
 304   AliCFDataGrid *spectrumSubtracted = new AliCFDataGrid("spectrumSubtracted", "Data Grid for spectrum after Background subtraction", *dataContainer,fStepData);
 305
 306   AliCFDataGrid *dataspectrumbeforesubstraction = (AliCFDataGrid *) ((AliCFDataGrid *)GetSpectrum(GetContainer(kDataContainer),fStepData))->Clone();
 307   dataspectrumbeforesubstraction->SetName("dataspectrumbeforesubstraction");
 308
 309
 310   // Background Estimate
 311   AliCFContainer *backgroundContainer = GetContainer(kBackgroundData);
 312   if(!backgroundContainer){
 313     AliError("MC background container not found");
 314     return NULL;
 315   }
 316
 317   Int_t stepbackground = 1;
 318   AliCFDataGrid *backgroundGrid = new AliCFDataGrid("ContaminationGrid","ContaminationGrid",*backgroundContainer,stepbackground);
 319
 320   // Subtract
 321   spectrumSubtracted->Add(backgroundGrid,-1.0);
 322
 323   // QA
 324   TH1D *subtractedspectrum = (TH1D *) spectrumSubtracted->Project(0);
 325   CorrectFromTheWidth(subtractedspectrum);
 326   TH1D *rawspectrum = (TH1D *) dataspectrumbeforesubstraction->Project(0);
 327   CorrectFromTheWidth(rawspectrum);
 328   fQA->AddResultAt(subtractedspectrum,AliHFEInclusiveSpectrumQA::kAfterSC);
 329   fQA->AddResultAt(rawspectrum,AliHFEInclusiveSpectrumQA::kBeforeSC);
 330   fQA->DrawSubtractContamination();
 331
 332   return spectrumSubtracted;
 333 }
 334
 335 //____________________________________________________________
 336 AliCFDataGrid *AliHFEInclusiveSpectrum::CorrectParametrizedEfficiency(AliCFDataGrid* const bgsubpectrum){
 337
 338   //
 339   // Apply TPC pid efficiency correction from parametrisation
 340   //
 341
 342   // Data in the right format
 343   AliCFDataGrid *dataGrid = 0x0;
 344   if(bgsubpectrum) {
 345     dataGrid = bgsubpectrum;
 346   }
 347   else {
 348
 349     AliCFContainer *dataContainer = GetContainer(kDataContainer);
 350     if(!dataContainer){
 351       AliError("Data Container not available");
 352       return NULL;
 353     }
 354     dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
 355   }
 356   AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
 357   result->SetName("ParametrizedEfficiencyBefore");
 358   THnSparse *h = result->GetGrid();
 359   Int_t nbdimensions = h->GetNdimensions();
 360   //printf("CorrectParametrizedEfficiency::We have dimensions %d\n",nbdimensions);
 361   AliCFContainer *dataContainer = GetContainer(kDataContainer);
 362   if(!dataContainer){
 363     AliError("Data Container not available");
 364     return NULL;
 365   }
 366   AliCFContainer *dataContainerbis = (AliCFContainer *) dataContainer->Clone();
 367   dataContainerbis->Add(dataContainerbis,-1.0);
 368
 369
 370   Int_t* coord = new Int_t[nbdimensions];
 371   memset(coord, 0, sizeof(Int_t) * nbdimensions);
 372   Double_t* points = new Double_t[nbdimensions];
 373
 374   ULong64_t nEntries = h->GetNbins();
 375   for (ULong64_t i = 0; i < nEntries; ++i) {
 376
 377     Double_t value = h->GetBinContent(i, coord);
 378     //Double_t valuecontainer = dataContainerbis->GetBinContent(coord,fStepData);
 379     //printf("Value %f, and valuecontainer %f\n",value,valuecontainer);
 380
 381     // Get the bin co-ordinates given an coord
 382     for (Int_t j = 0; j < nbdimensions; ++j)
 383       points[j] = h->GetAxis(j)->GetBinCenter(coord[j]);
 384
 385     if (!fEfficiencyFunction->IsInside(points))
 386          continue;
 387     TF1::RejectPoint(kFALSE);
 388
 389     // Evaulate function at points
 390     Double_t valueEfficiency = fEfficiencyFunction->EvalPar(points, NULL);
 391     //printf("Value efficiency is %f\n",valueEfficiency);
 392
 393     if(valueEfficiency > 0.0) {
 394       h->SetBinContent(coord,value/valueEfficiency);
 395       dataContainerbis->SetBinContent(coord,fStepData,value/valueEfficiency);
 396     }
 397     Double_t error = h->GetBinError(i);
 398     h->SetBinError(coord,error/valueEfficiency);
 399     dataContainerbis->SetBinError(coord,fStepData,error/valueEfficiency);
 400
 401
 402   }
 403
 404   delete[] coord;
 405   delete[] points;
 406
 407   AliCFDataGrid *resultt = new AliCFDataGrid("spectrumEfficiencyParametrized", "Data Grid for spectrum after Efficiency parametrized", *dataContainerbis,fStepData);
 408
 409  // QA
 410   TH1D *afterE = (TH1D *) resultt->Project(0);
 411   CorrectFromTheWidth(afterE);
 412   TH1D *beforeE = (TH1D *) dataGrid->Project(0);
 413   CorrectFromTheWidth(beforeE);
 414   fQA->AddResultAt(afterE,AliHFEInclusiveSpectrumQA::kAfterPE);
 415   fQA->AddResultAt(beforeE,AliHFEInclusiveSpectrumQA::kBeforePE);
 416   fQA->AddResultAt(fEfficiencyFunction,AliHFEInclusiveSpectrumQA::kPEfficiency);
 417   fQA->DrawCorrectWithEfficiency(AliHFEInclusiveSpectrumQA::kParametrized);
 418
 419   return resultt;
 420
 421 }
 422 //____________________________________________________________
 423 AliCFDataGrid *AliHFEInclusiveSpectrum::CorrectV0Efficiency(AliCFDataGrid* const bgsubpectrum){
 424
 425   //
 426   // Apply TPC pid efficiency correction from V0
 427   //
 428
 429   AliCFContainer *v0Container = GetContainer(kDataContainerV0);
 430   if(!v0Container){
 431     AliError("V0 Container not available");
 432     return NULL;
 433   }
 434
 435   // Efficiency
 436   AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*v0Container);
 437   efficiencyD->CalculateEfficiency(fStepAfterCutsV0,fStepBeforeCutsV0);
 438
 439   // Data in the right format
 440   AliCFDataGrid *dataGrid = 0x0;
 441   if(bgsubpectrum) {
 442     dataGrid = bgsubpectrum;
 443   }
 444   else {
 445     AliCFContainer *dataContainer = GetContainer(kDataContainer);
 446     if(!dataContainer){
 447       AliError("Data Container not available");
 448       return NULL;
 449     }
 450     dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
 451   }
 452
 453   // Correct
 454   AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
 455   result->ApplyEffCorrection(*efficiencyD);
 456
 457   // QA
 458   TH1D *afterE = (TH1D *) result->Project(0);
 459   CorrectFromTheWidth(afterE);
 460   TH1D *beforeE = (TH1D *) dataGrid->Project(0);
 461   CorrectFromTheWidth(beforeE);
 462   TH1D* efficiencyDproj = (TH1D *) efficiencyD->Project(0);
 463   fQA->AddResultAt(afterE,AliHFEInclusiveSpectrumQA::kAfterV0);
 464   fQA->AddResultAt(beforeE,AliHFEInclusiveSpectrumQA::kBeforeV0);
 465   fQA->AddResultAt(efficiencyDproj,AliHFEInclusiveSpectrumQA::kV0Efficiency);
 466   fQA->DrawCorrectWithEfficiency(AliHFEInclusiveSpectrumQA::kV0);
 467
 468
 469   return result;
 470
 471 }
 472 //____________________________________________________________
 473 THnSparse *AliHFEInclusiveSpectrum::Unfold(AliCFDataGrid* const bgsubpectrum){
 474
 475   //
 476   // Return the unfolded spectrum
 477   //
 478
 479   AliCFContainer *mcContainer = GetContainer(kMCContainerMC);
 480   if(!mcContainer){
 481     AliError("MC Container not available");
 482     return NULL;
 483   }
 484
 485   if(!fCorrelation){
 486     AliError("No Correlation map available");
 487     return NULL;
 488   }
 489
 490   // Data
 491   AliCFDataGrid *dataGrid = 0x0;
 492   if(bgsubpectrum) {
 493     dataGrid = bgsubpectrum;
 494   }
 495   else {
 496
 497     AliCFContainer *dataContainer = GetContainer(kDataContainer);
 498     if(!dataContainer){
 499       AliError("Data Container not available");
 500       return NULL;
 501     }
 502
 503     dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
 504   }
 505
 506   // Guessed
 507   AliCFDataGrid* guessedGrid = new AliCFDataGrid("guessed","",*mcContainer, fStepGuessedUnfolding);
 508   THnSparse* guessedTHnSparse = ((AliCFGridSparse*)guessedGrid->GetData())->GetGrid();
 509
 510   // Efficiency
 511   AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
 512   efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
 513
 514   // Unfold
 515
 516   AliCFUnfolding unfolding("unfolding","",fNbDimensions,fCorrelation,efficiencyD->GetGrid(),dataGrid->GetGrid(),guessedTHnSparse,1.e-06,0,fNumberOfIterations);
 517   if(fSetSmoothing) unfolding.UseSmoothing();
 518   unfolding.Unfold();
 519
 520   // Results
 521   THnSparse* result = unfolding.GetUnfolded();
 522   THnSparse* residual = unfolding.GetEstMeasured();
 523
 524   // QA
 525   TH1D *residualh = (TH1D *) residual->Projection(0);
 526   TH1D *beforeE = (TH1D *) dataGrid->Project(0);
 527   TH1D* efficiencyDproj = (TH1D *) efficiencyD->Project(0);
 528   TH1D *afterE = (TH1D *) result->Projection(0);
 529   CorrectFromTheWidth(residualh);
 530   CorrectFromTheWidth(beforeE);
 531   CorrectFromTheWidth(afterE);
 532   fQA->AddResultAt(residualh,AliHFEInclusiveSpectrumQA::kResidualU);
 533   fQA->AddResultAt(afterE,AliHFEInclusiveSpectrumQA::kAfterU);
 534   fQA->AddResultAt(beforeE,AliHFEInclusiveSpectrumQA::kBeforeU);
 535   fQA->AddResultAt(efficiencyDproj,AliHFEInclusiveSpectrumQA::kUEfficiency);
 536   fQA->DrawUnfolding();
 537
 538   return (THnSparse *) result->Clone();
 539
 540 }
 541 //____________________________________________________________
 542 AliCFDataGrid *AliHFEInclusiveSpectrum::CorrectForEfficiency(AliCFDataGrid* const bgsubpectrum){
 543
 544   //
 545   // Apply unfolding and efficiency correction together to bgsubspectrum
 546   //
 547
 548   AliCFContainer *mcContainer = GetContainer(kMCContainerESD);
 549   if(!mcContainer){
 550     AliError("MC Container not available");
 551     return NULL;
 552   }
 553
 554   // Efficiency
 555   AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
 556   efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
 557
 558   // Data in the right format
 559   AliCFDataGrid *dataGrid = 0x0;
 560   if(bgsubpectrum) {
 561     dataGrid = bgsubpectrum;
 562   }
 563   else {
 564
 565     AliCFContainer *dataContainer = GetContainer(kDataContainer);
 566     if(!dataContainer){
 567       AliError("Data Container not available");
 568       return NULL;
 569     }
 570
 571     dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
 572   }
 573
 574   // Correct
 575   AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
 576   result->ApplyEffCorrection(*efficiencyD);
 577
 578   // QA
 579   TH1D *afterE = (TH1D *) result->Project(0);
 580   CorrectFromTheWidth(afterE);
 581   TH1D *beforeE = (TH1D *) dataGrid->Project(0);
 582   CorrectFromTheWidth(beforeE);
 583   TH1D* efficiencyDproj = (TH1D *) efficiencyD->Project(0);
 584   fQA->AddResultAt(afterE,AliHFEInclusiveSpectrumQA::kAfterMCE);
 585   fQA->AddResultAt(beforeE,AliHFEInclusiveSpectrumQA::kBeforeMCE);
 586   fQA->AddResultAt(efficiencyDproj,AliHFEInclusiveSpectrumQA::kMCEfficiency);
 587   fQA->DrawCorrectWithEfficiency(AliHFEInclusiveSpectrumQA::kMC);
 588
 589
 590   return result;
 591
 592 }
 593 //____________________________________________________________
 594 void AliHFEInclusiveSpectrum::WriteResults(const char *filename)
 595 {
 596   //
 597   // Write results
 598   //
 599
 600   AliCFContainer *dataContainer = GetContainer(kDataContainer);
 601   AliCFContainer *mcContainer = GetContainer(kMCContainerMC);
 602   TObject *unfolded = 0x0;
 603   TObject *correctedspectrum = 0x0;
 604   if(fQA) {
 605     unfolded = fQA->GetResult(AliHFEInclusiveSpectrumQA::kFinalResultUnfolded);
 606     correctedspectrum = fQA->GetResult(AliHFEInclusiveSpectrumQA::kFinalResultDirectEfficiency);
 607   }
 608
 609   TFile *file = TFile::Open(filename,"recreate");
 610   if(dataContainer) dataContainer->Write("data");
 611   if(mcContainer) mcContainer->Write("mcefficiency");
 612   if(fCorrelation) fCorrelation->Write("correlationmatrix");
 613   if(unfolded) unfolded->Write("unfoldedspectrum");
 614   if(correctedspectrum) correctedspectrum->Write("correctedspectrum");
 615   if(fQA) fQA->Write("QAResults");
 616   file->Close();
 617
 618 }
 619