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 */*bghfecontainer*/, const AliHFEcontainer *v0hfecontainer,AliCFContainer *photoniccontainerD){
 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
 153   // Photonic Background
 154   SetContainer(photoniccontainerD,AliHFECorrectSpectrumBase::kPhotonicBackground);
 155
 156   // QA
 157   Int_t dimqa = datacontainer->GetNVar();
 158   Int_t dimsqa[dimqa];
 159   for(Int_t i = 0; i < dimqa; i++) dimsqa[i] = i;
 160   AliCFContainer *datacontainerDQA = GetSlicedContainer(datacontainer, dimqa, dimsqa, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
 161   fQA->AddResultAt(datacontainerDQA,AliHFEInclusiveSpectrumQA::kDataProjection);
 162
 163   //
 164   // MC container: ESD/MC efficiency maps + MC/MC efficiency maps
 165   //
 166   AliCFContainer *mccontaineresd = 0x0;
 167   AliCFContainer *mccontainermc = 0x0;
 168   mccontaineresd = mchfecontainer->MakeMergedCFContainer("sumesd","sumesd","MCTrackCont:recTrackContReco");
 169   mccontainermc = mchfecontainer->MakeMergedCFContainer("summc","summc","MCTrackCont:recTrackContMC");
 170   if((!mccontaineresd) || (!mccontainermc)) return kFALSE;
 171   Int_t source = -1;
 172   AliCFContainer *mccontaineresdD = GetSlicedContainer(mccontaineresd, fNbDimensions, dims, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
 173   AliCFContainer *mccontainermcD = GetSlicedContainer(mccontainermc, fNbDimensions, dims, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
 174   if((!mccontaineresdD) || (!mccontainermcD)) return kFALSE;
 175   SetContainer(mccontainermcD,AliHFECorrectSpectrumBase::kMCContainerMC);
 176   SetContainer(mccontaineresdD,AliHFECorrectSpectrumBase::kMCContainerESD);
 177
 178   //
 179   // Correlation matrix
 180   //
 181   THnSparseF *mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepbeforePID");
 182   if(!mccorrelation) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID");
 183   if(!mccorrelation) return kFALSE;
 184   THnSparseF *mccorrelationD = GetSlicedCorrelation(mccorrelation, fNbDimensions, dims,fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
 185   if(!mccorrelationD) {
 186     printf("No correlation\n");
 187     return kFALSE;
 188   }
 189   SetCorrelation(mccorrelationD);
 190   // QA
 191   fQA->AddResultAt(mccorrelation,AliHFEInclusiveSpectrumQA::kCMProjection);
 192
 193   //
 194   // V0 container Electron, pt eta phi
 195   //
 196   if(v0hfecontainer) {
 197     AliCFContainer *containerV0 = v0hfecontainer->GetCFContainer("taggedTrackContainerReco");
 198     if(!containerV0) return kFALSE;
 199     AliCFContainer *containerV0Electron = GetSlicedContainer(containerV0, fNbDimensions, dims, AliPID::kElectron,fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
 200     if(!containerV0Electron) return kFALSE;
 201     SetContainer(containerV0Electron,AliHFECorrectSpectrumBase::kDataContainerV0);
 202   }
 203
 204   //
 205   fQA->DrawProjections();
 206
 207
 208   return kTRUE;
 209 }
 210 //____________________________________________________________
 211 Bool_t AliHFEInclusiveSpectrum::Correct(Bool_t subtractcontamination, Bool_t subtractphotonic){
 212   //
 213   // Correct the spectrum for efficiency and unfolding
 214   // with both method and compare
 215   //
 216
 217   gStyle->SetPalette(1);
 218   gStyle->SetOptStat(1111);
 219   gStyle->SetPadBorderMode(0);
 220   gStyle->SetCanvasColor(10);
 221   gStyle->SetPadLeftMargin(0.13);
 222   gStyle->SetPadRightMargin(0.13);
 223
 224   printf("Steps are: stepdata %d, stepMC %d, steptrue %d, stepV0after %d, stepV0before %d\n",fStepData,fStepMC,fStepTrue,fStepAfterCutsV0,fStepBeforeCutsV0);
 225
 226   ///////////////////////////
 227   // Check initialization
 228   ///////////////////////////
 229
 230   if((!GetContainer(kDataContainer)) || (!GetContainer(kMCContainerMC)) || (!GetContainer(kMCContainerESD))){
 231     AliInfo("You have to init before");
 232     return kFALSE;
 233   }
 234
 235   if((fStepTrue < 0) && (fStepMC < 0) && (fStepData < 0)) {
 236     AliInfo("You have to set the steps before: SetMCTruthStep, SetMCEffStep, SetStepToCorrect");
 237     return kFALSE;
 238   }
 239
 240   SetStepGuessedUnfolding(AliHFEcuts::kStepRecKineITSTPC + AliHFEcuts::kNcutStepsMCTrack);
 241
 242   AliCFDataGrid *dataGridAfterFirstSteps = 0x0;
 243
 244   //////////////////////////////////
 245   // Subtract hadron background
 246   /////////////////////////////////
 247   AliCFDataGrid *dataspectrumaftersubstraction = 0x0;
 248   if(subtractcontamination) {
 249     dataspectrumaftersubstraction = SubtractBackground();
 250     dataGridAfterFirstSteps = dataspectrumaftersubstraction;
 251   }
 252
 253   //////////////////////////////////
 254   // Subtract Photonic background
 255   /////////////////////////////////
 256   AliCFDataGrid *dataspectrumafterphotonicsubstraction = 0x0;
 257   if(subtractphotonic) {
 258     dataspectrumafterphotonicsubstraction = SubtractPhotonicBackground();
 259     dataGridAfterFirstSteps = dataspectrumafterphotonicsubstraction;
 260   }
 261
 262   ////////////////////////////////////////////////
 263   // Correct for TPC efficiency from V0 if any
 264   ///////////////////////////////////////////////
 265   AliCFDataGrid *dataspectrumafterV0efficiencycorrection = 0x0;
 266   AliCFContainer *dataContainerV0 = GetContainer(kDataContainerV0);
 267   if(dataContainerV0){
 268     dataspectrumafterV0efficiencycorrection = CorrectV0Efficiency(dataspectrumaftersubstraction);
 269     dataGridAfterFirstSteps = dataspectrumafterV0efficiencycorrection;
 270   }
 271
 272   //////////////////////////////////////////////////////////////////////////////
 273   // Correct for efficiency parametrized (if TPC efficiency is parametrized)
 274   /////////////////////////////////////////////////////////////////////////////
 275   AliCFDataGrid *dataspectrumafterefficiencyparametrizedcorrection = 0x0;
 276   if(fEfficiencyFunction){
 277     dataspectrumafterefficiencyparametrizedcorrection = CorrectParametrizedEfficiency(dataGridAfterFirstSteps);
 278     dataGridAfterFirstSteps = dataspectrumafterefficiencyparametrizedcorrection;
 279   }
 280
 281   ///////////////
 282   // Unfold
 283   //////////////
 284   THnSparse *correctedspectrum = Unfold(dataGridAfterFirstSteps);
 285   if(!correctedspectrum){
 286     AliError("No corrected spectrum\n");
 287     return kFALSE;
 288   }
 289
 290   /////////////////////
 291   // Simply correct
 292   ////////////////////
 293   AliCFDataGrid *alltogetherCorrection = CorrectForEfficiency(dataGridAfterFirstSteps);
 294
 295   // QA final results
 296   TGraphErrors* correctedspectrumD = Normalize(correctedspectrum);
 297   TGraphErrors* alltogetherspectrumD = Normalize(alltogetherCorrection);
 298   fQA->AddResultAt(correctedspectrumD,AliHFEInclusiveSpectrumQA::kFinalResultUnfolded);
 299   fQA->AddResultAt(alltogetherspectrumD,AliHFEInclusiveSpectrumQA::kFinalResultDirectEfficiency);
 300   fQA->DrawResult();
 301
 302   return kTRUE;
 303 }
 304
 305 //____________________________________________________________
 306 AliCFDataGrid* AliHFEInclusiveSpectrum::SubtractBackground(){
 307   //
 308   // Apply background subtraction
 309   //
 310
 311   // Raw spectrum
 312   AliCFContainer *dataContainer = GetContainer(kDataContainer);
 313   if(!dataContainer){
 314     AliError("Data Container not available");
 315     return NULL;
 316   }
 317   printf("Step data: %d\n",fStepData);
 318   AliCFDataGrid *spectrumSubtracted = new AliCFDataGrid("spectrumSubtracted", "Data Grid for spectrum after Background subtraction", *dataContainer,fStepData);
 319
 320   AliCFDataGrid *dataspectrumbeforesubstraction = (AliCFDataGrid *) ((AliCFDataGrid *)GetSpectrum(GetContainer(kDataContainer),fStepData))->Clone();
 321   dataspectrumbeforesubstraction->SetName("dataspectrumbeforesubstraction");
 322
 323
 324   // Background Estimate
 325   AliCFContainer *backgroundContainer = GetContainer(kBackgroundData);
 326   if(!backgroundContainer){
 327     AliError("MC background container not found");
 328     return NULL;
 329   }
 330
 331   Int_t stepbackground = 1;
 332   AliCFDataGrid *backgroundGrid = new AliCFDataGrid("ContaminationGrid","ContaminationGrid",*backgroundContainer,stepbackground);
 333
 334   // Subtract
 335   spectrumSubtracted->Add(backgroundGrid,-1.0);
 336
 337   // QA
 338   TH1D *subtractedspectrum = (TH1D *) spectrumSubtracted->Project(0);
 339   CorrectFromTheWidth(subtractedspectrum);
 340   TH1D *rawspectrum = (TH1D *) dataspectrumbeforesubstraction->Project(0);
 341   CorrectFromTheWidth(rawspectrum);
 342   fQA->AddResultAt(subtractedspectrum,AliHFEInclusiveSpectrumQA::kAfterSC);
 343   fQA->AddResultAt(rawspectrum,AliHFEInclusiveSpectrumQA::kBeforeSC);
 344   fQA->DrawSubtractContamination();
 345
 346   return spectrumSubtracted;
 347 }
 348
 349 //____________________________________________________________
 350 AliCFDataGrid* AliHFEInclusiveSpectrum::SubtractPhotonicBackground(){
 351   //
 352   // Apply Photonic background subtraction
 353   //
 354
 355   printf("Photonic Background Subtraction \n");
 356
 357   // Raw spectrum
 358   AliCFContainer *dataContainer = GetContainer(kDataContainer);
 359   if(!dataContainer){
 360     AliError("Data Container not available");
 361     return NULL;
 362   }
 363   printf("Step data: %d\n",fStepData);
 364   AliCFDataGrid *spectrumPhotonicSubtracted = new AliCFDataGrid("spectrumPhotonicSubtracted", "Data Grid for spectrum after Photonic Background subtraction", *dataContainer,fStepData);
 365
 366   AliCFDataGrid *dataSpectrumBeforePhotonicSubstraction = (AliCFDataGrid *) ((AliCFDataGrid *)GetSpectrum(GetContainer(kDataContainer),fStepData))->Clone();
 367   dataSpectrumBeforePhotonicSubstraction->SetName("dataSpectrumBeforePhotonicSubstraction");
 368
 369
 370   // Background Estimate
 371   AliCFContainer *photonicContainer = GetContainer(kPhotonicBackground);
 372   if(!photonicContainer){
 373     AliError("Photonic background container not found");
 374     return NULL;
 375   }
 376
 377   Int_t stepbackground = 0;
 378   AliCFDataGrid *photonicGrid = new AliCFDataGrid("ContaminationGrid","ContaminationGrid",*photonicContainer,stepbackground);
 379
 380   // Subtract
 381   spectrumPhotonicSubtracted->Add(photonicGrid,-1.0);
 382
 383   // QA
 384   TH1D *photonicsubtractedspectrum = (TH1D *) spectrumPhotonicSubtracted->Project(0);
 385   CorrectFromTheWidth(photonicsubtractedspectrum);
 386   TH1D *newrawspectrum = (TH1D *) dataSpectrumBeforePhotonicSubstraction->Project(0);
 387   CorrectFromTheWidth(newrawspectrum);
 388   fQA->AddResultAt(photonicsubtractedspectrum,AliHFEInclusiveSpectrumQA::kAfterSPB);
 389   fQA->AddResultAt(newrawspectrum,AliHFEInclusiveSpectrumQA::kBeforeSPB);
 390   fQA->DrawSubtractPhotonicBackground();
 391
 392   return spectrumPhotonicSubtracted;
 393 }
 394
 395
 396 //____________________________________________________________
 397 AliCFDataGrid *AliHFEInclusiveSpectrum::CorrectParametrizedEfficiency(AliCFDataGrid* const bgsubpectrum){
 398
 399   //
 400   // Apply TPC pid efficiency correction from parametrisation
 401   //
 402
 403   // Data in the right format
 404   AliCFDataGrid *dataGrid = 0x0;
 405   if(bgsubpectrum) {
 406     dataGrid = bgsubpectrum;
 407   }
 408   else {
 409
 410     AliCFContainer *dataContainer = GetContainer(kDataContainer);
 411     if(!dataContainer){
 412       AliError("Data Container not available");
 413       return NULL;
 414     }
 415     dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
 416   }
 417   AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
 418   result->SetName("ParametrizedEfficiencyBefore");
 419   THnSparse *h = result->GetGrid();
 420   Int_t nbdimensions = h->GetNdimensions();
 421   //printf("CorrectParametrizedEfficiency::We have dimensions %d\n",nbdimensions);
 422   AliCFContainer *dataContainer = GetContainer(kDataContainer);
 423   if(!dataContainer){
 424     AliError("Data Container not available");
 425     return NULL;
 426   }
 427   AliCFContainer *dataContainerbis = (AliCFContainer *) dataContainer->Clone();
 428   dataContainerbis->Add(dataContainerbis,-1.0);
 429
 430
 431   Int_t* coord = new Int_t[nbdimensions];
 432   memset(coord, 0, sizeof(Int_t) * nbdimensions);
 433   Double_t* points = new Double_t[nbdimensions];
 434
 435   ULong64_t nEntries = h->GetNbins();
 436   for (ULong64_t i = 0; i < nEntries; ++i) {
 437
 438     Double_t value = h->GetBinContent(i, coord);
 439     //Double_t valuecontainer = dataContainerbis->GetBinContent(coord,fStepData);
 440     //printf("Value %f, and valuecontainer %f\n",value,valuecontainer);
 441
 442     // Get the bin co-ordinates given an coord
 443     for (Int_t j = 0; j < nbdimensions; ++j)
 444       points[j] = h->GetAxis(j)->GetBinCenter(coord[j]);
 445
 446     if (!fEfficiencyFunction->IsInside(points))
 447          continue;
 448     TF1::RejectPoint(kFALSE);
 449
 450     // Evaulate function at points
 451     Double_t valueEfficiency = fEfficiencyFunction->EvalPar(points, NULL);
 452     //printf("Value efficiency is %f\n",valueEfficiency);
 453
 454     if(valueEfficiency > 0.0) {
 455       h->SetBinContent(coord,value/valueEfficiency);
 456       dataContainerbis->SetBinContent(coord,fStepData,value/valueEfficiency);
 457     }
 458     Double_t error = h->GetBinError(i);
 459     h->SetBinError(coord,error/valueEfficiency);
 460     dataContainerbis->SetBinError(coord,fStepData,error/valueEfficiency);
 461
 462
 463   }
 464
 465   delete[] coord;
 466   delete[] points;
 467
 468   AliCFDataGrid *resultt = new AliCFDataGrid("spectrumEfficiencyParametrized", "Data Grid for spectrum after Efficiency parametrized", *dataContainerbis,fStepData);
 469
 470  // QA
 471   TH1D *afterE = (TH1D *) resultt->Project(0);
 472   CorrectFromTheWidth(afterE);
 473   TH1D *beforeE = (TH1D *) dataGrid->Project(0);
 474   CorrectFromTheWidth(beforeE);
 475   fQA->AddResultAt(afterE,AliHFEInclusiveSpectrumQA::kAfterPE);
 476   fQA->AddResultAt(beforeE,AliHFEInclusiveSpectrumQA::kBeforePE);
 477   fQA->AddResultAt(fEfficiencyFunction,AliHFEInclusiveSpectrumQA::kPEfficiency);
 478   fQA->DrawCorrectWithEfficiency(AliHFEInclusiveSpectrumQA::kParametrized);
 479
 480   return resultt;
 481
 482 }
 483 //____________________________________________________________
 484 AliCFDataGrid *AliHFEInclusiveSpectrum::CorrectV0Efficiency(AliCFDataGrid* const bgsubpectrum){
 485
 486   //
 487   // Apply TPC pid efficiency correction from V0
 488   //
 489
 490   AliCFContainer *v0Container = GetContainer(kDataContainerV0);
 491   if(!v0Container){
 492     AliError("V0 Container not available");
 493     return NULL;
 494   }
 495
 496   // Efficiency
 497   AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*v0Container);
 498   efficiencyD->CalculateEfficiency(fStepAfterCutsV0,fStepBeforeCutsV0);
 499
 500   // Data in the right format
 501   AliCFDataGrid *dataGrid = 0x0;
 502   if(bgsubpectrum) {
 503     dataGrid = bgsubpectrum;
 504   }
 505   else {
 506     AliCFContainer *dataContainer = GetContainer(kDataContainer);
 507     if(!dataContainer){
 508       AliError("Data Container not available");
 509       return NULL;
 510     }
 511     dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
 512   }
 513
 514   // Correct
 515   AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
 516   result->ApplyEffCorrection(*efficiencyD);
 517
 518   // QA
 519   TH1D *afterE = (TH1D *) result->Project(0);
 520   CorrectFromTheWidth(afterE);
 521   TH1D *beforeE = (TH1D *) dataGrid->Project(0);
 522   CorrectFromTheWidth(beforeE);
 523   TH1D* efficiencyDproj = (TH1D *) efficiencyD->Project(0);
 524   fQA->AddResultAt(afterE,AliHFEInclusiveSpectrumQA::kAfterV0);
 525   fQA->AddResultAt(beforeE,AliHFEInclusiveSpectrumQA::kBeforeV0);
 526   fQA->AddResultAt(efficiencyDproj,AliHFEInclusiveSpectrumQA::kV0Efficiency);
 527   fQA->DrawCorrectWithEfficiency(AliHFEInclusiveSpectrumQA::kV0);
 528
 529
 530   return result;
 531
 532 }
 533 //____________________________________________________________
 534 THnSparse *AliHFEInclusiveSpectrum::Unfold(AliCFDataGrid* const bgsubpectrum){
 535
 536   //
 537   // Return the unfolded spectrum
 538   //
 539
 540   AliCFContainer *mcContainer = GetContainer(kMCContainerMC);
 541   if(!mcContainer){
 542     AliError("MC Container not available");
 543     return NULL;
 544   }
 545
 546   if(!fCorrelation){
 547     AliError("No Correlation map available");
 548     return NULL;
 549   }
 550
 551   // Data
 552   AliCFDataGrid *dataGrid = 0x0;
 553   if(bgsubpectrum) {
 554     dataGrid = bgsubpectrum;
 555   }
 556   else {
 557
 558     AliCFContainer *dataContainer = GetContainer(kDataContainer);
 559     if(!dataContainer){
 560       AliError("Data Container not available");
 561       return NULL;
 562     }
 563
 564     dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
 565   }
 566
 567   // Guessed
 568   AliCFDataGrid* guessedGrid = new AliCFDataGrid("guessed","",*mcContainer, fStepGuessedUnfolding);
 569   THnSparse* guessedTHnSparse = ((AliCFGridSparse*)guessedGrid->GetData())->GetGrid();
 570
 571   // Efficiency
 572   AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
 573   efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
 574
 575   // Unfold
 576
 577   AliCFUnfolding unfolding("unfolding","",fNbDimensions,fCorrelation,efficiencyD->GetGrid(),dataGrid->GetGrid(),guessedTHnSparse,1.e-06,0,fNumberOfIterations);
 578   if(fSetSmoothing) unfolding.UseSmoothing();
 579   unfolding.Unfold();
 580
 581   // Results
 582   THnSparse* result = unfolding.GetUnfolded();
 583   THnSparse* residual = unfolding.GetEstMeasured();
 584
 585   // QA
 586   TH1D *residualh = (TH1D *) residual->Projection(0);
 587   TH1D *beforeE = (TH1D *) dataGrid->Project(0);
 588   TH1D* efficiencyDproj = (TH1D *) efficiencyD->Project(0);
 589   TH1D *afterE = (TH1D *) result->Projection(0);
 590   CorrectFromTheWidth(residualh);
 591   CorrectFromTheWidth(beforeE);
 592   CorrectFromTheWidth(afterE);
 593   fQA->AddResultAt(residualh,AliHFEInclusiveSpectrumQA::kResidualU);
 594   fQA->AddResultAt(afterE,AliHFEInclusiveSpectrumQA::kAfterU);
 595   fQA->AddResultAt(beforeE,AliHFEInclusiveSpectrumQA::kBeforeU);
 596   fQA->AddResultAt(efficiencyDproj,AliHFEInclusiveSpectrumQA::kUEfficiency);
 597   fQA->DrawUnfolding();
 598
 599   return (THnSparse *) result->Clone();
 600
 601 }
 602 //____________________________________________________________
 603 AliCFDataGrid *AliHFEInclusiveSpectrum::CorrectForEfficiency(AliCFDataGrid* const bgsubpectrum){
 604
 605   //
 606   // Apply unfolding and efficiency correction together to bgsubspectrum
 607   //
 608
 609   AliCFContainer *mcContainer = GetContainer(kMCContainerESD);
 610   if(!mcContainer){
 611     AliError("MC Container not available");
 612     return NULL;
 613   }
 614
 615   // Efficiency
 616   AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
 617   efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
 618
 619   // Data in the right format
 620   AliCFDataGrid *dataGrid = 0x0;
 621   if(bgsubpectrum) {
 622     dataGrid = bgsubpectrum;
 623   }
 624   else {
 625
 626     AliCFContainer *dataContainer = GetContainer(kDataContainer);
 627     if(!dataContainer){
 628       AliError("Data Container not available");
 629       return NULL;
 630     }
 631
 632     dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
 633   }
 634
 635   // Correct
 636   AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
 637   result->ApplyEffCorrection(*efficiencyD);
 638
 639   // QA
 640   TH1D *afterE = (TH1D *) result->Project(0);
 641   CorrectFromTheWidth(afterE);
 642   TH1D *beforeE = (TH1D *) dataGrid->Project(0);
 643   CorrectFromTheWidth(beforeE);
 644   TH1D* efficiencyDproj = (TH1D *) efficiencyD->Project(0);
 645   fQA->AddResultAt(afterE,AliHFEInclusiveSpectrumQA::kAfterMCE);
 646   fQA->AddResultAt(beforeE,AliHFEInclusiveSpectrumQA::kBeforeMCE);
 647   fQA->AddResultAt(efficiencyDproj,AliHFEInclusiveSpectrumQA::kMCEfficiency);
 648   fQA->DrawCorrectWithEfficiency(AliHFEInclusiveSpectrumQA::kMC);
 649
 650   return result;
 651
 652 }
 653 //____________________________________________________________
 654 void AliHFEInclusiveSpectrum::WriteResults(const char *filename)
 655 {
 656   //
 657   // Write results
 658   //
 659
 660   AliCFContainer *dataContainer = GetContainer(kDataContainer);
 661   AliCFContainer *mcContainer = GetContainer(kMCContainerMC);
 662   TObject *unfolded = 0x0;
 663   TObject *correctedspectrum = 0x0;
 664   if(fQA) {
 665     unfolded = fQA->GetResult(AliHFEInclusiveSpectrumQA::kFinalResultUnfolded);
 666     correctedspectrum = fQA->GetResult(AliHFEInclusiveSpectrumQA::kFinalResultDirectEfficiency);
 667   }
 668
 669   TFile *file = TFile::Open(filename,"recreate");
 670   if(dataContainer) dataContainer->Write("data");
 671   if(mcContainer) mcContainer->Write("mcefficiency");
 672   if(fCorrelation) fCorrelation->Write("correlationmatrix");
 673   if(unfolded) unfolded->Write("unfoldedspectrum");
 674   if(correctedspectrum) correctedspectrum->Write("correctedspectrum");
 675   if(fQA) fQA->Write("QAResults");
 676   file->Close();
 677
 678 }
 679