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