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     correctedspectrumD = Normalize(dataGridAfterFirstSteps);
 301     alltogetherspectrumD = Normalize(dataGridAfterFirstSteps);
 302   }
 303
 304   // QA final results
 305
 306   fQA->AddResultAt(correctedspectrumD,AliHFEInclusiveSpectrumQA::kFinalResultUnfolded);
 307   fQA->AddResultAt(alltogetherspectrumD,AliHFEInclusiveSpectrumQA::kFinalResultDirectEfficiency);
 308   fQA->DrawResult();
 309
 310   return kTRUE;
 311 }
 312
 313 //____________________________________________________________
 314 AliCFDataGrid* AliHFEInclusiveSpectrum::SubtractBackground(){
 315   //
 316   // Apply background subtraction
 317   //
 318
 319   // Raw spectrum
 320   AliCFContainer *dataContainer = GetContainer(kDataContainer);
 321   if(!dataContainer){
 322     AliError("Data Container not available");
 323     return NULL;
 324   }
 325   printf("Step data: %d\n",fStepData);
 326   AliCFDataGrid *spectrumSubtracted = new AliCFDataGrid("spectrumSubtracted", "Data Grid for spectrum after Background subtraction", *dataContainer,fStepData);
 327
 328   AliCFDataGrid *dataspectrumbeforesubstraction = (AliCFDataGrid *) ((AliCFDataGrid *)GetSpectrum(GetContainer(kDataContainer),fStepData))->Clone();
 329   dataspectrumbeforesubstraction->SetName("dataspectrumbeforesubstraction");
 330
 331
 332   // Background Estimate
 333   AliCFContainer *backgroundContainer = GetContainer(kBackgroundData);
 334   if(!backgroundContainer){
 335     AliError("MC background container not found");
 336     return NULL;
 337   }
 338
 339   Int_t stepbackground = 1;
 340   AliCFDataGrid *backgroundGrid = new AliCFDataGrid("ContaminationGrid","ContaminationGrid",*backgroundContainer,stepbackground);
 341
 342   // Subtract
 343   spectrumSubtracted->Add(backgroundGrid,-1.0);
 344
 345   // QA
 346   TH1D *subtractedspectrum = (TH1D *) spectrumSubtracted->Project(0);
 347   AliHFEtools::NormaliseBinWidth(subtractedspectrum);
 348   TH1D *rawspectrum = (TH1D *) dataspectrumbeforesubstraction->Project(0);
 349   AliHFEtools::NormaliseBinWidth(rawspectrum);
 350   fQA->AddResultAt(subtractedspectrum,AliHFEInclusiveSpectrumQA::kAfterSC);
 351   fQA->AddResultAt(rawspectrum,AliHFEInclusiveSpectrumQA::kBeforeSC);
 352   fQA->DrawSubtractContamination();
 353
 354   if(fNbDimensions>=2) {
 355     fQA->AddResultAt((TObject *) spectrumSubtracted,AliHFEInclusiveSpectrumQA::kAfterSCND);
 356     fQA->AddResultAt((TObject *) dataspectrumbeforesubstraction,AliHFEInclusiveSpectrumQA::kBeforeSCND);
 357     fQA->AddResultAt((TObject *) backgroundGrid,AliHFEInclusiveSpectrumQA::kHadronContaminationND);
 358     fQA->DrawSubtractContaminationND();
 359   }
 360
 361
 362   return spectrumSubtracted;
 363 }
 364
 365 //____________________________________________________________
 366 AliCFDataGrid* AliHFEInclusiveSpectrum::SubtractPhotonicBackground(){
 367   //
 368   // Apply Photonic background subtraction
 369   //
 370
 371   printf("Photonic Background Subtraction \n");
 372
 373   // Raw spectrum
 374   AliCFContainer *dataContainer = GetContainer(kDataContainer);
 375   if(!dataContainer){
 376     AliError("Data Container not available");
 377     return NULL;
 378   }
 379   printf("Step data: %d\n",fStepData);
 380   AliCFDataGrid *spectrumPhotonicSubtracted = new AliCFDataGrid("spectrumPhotonicSubtracted", "Data Grid for spectrum after Photonic Background subtraction", *dataContainer,fStepData);
 381
 382   AliCFDataGrid *dataSpectrumBeforePhotonicSubstraction = (AliCFDataGrid *) ((AliCFDataGrid *)GetSpectrum(GetContainer(kDataContainer),fStepData))->Clone();
 383   dataSpectrumBeforePhotonicSubstraction->SetName("dataSpectrumBeforePhotonicSubstraction");
 384
 385
 386   // Background Estimate
 387   AliCFContainer *photonicContainer = GetContainer(kPhotonicBackground);
 388   if(!photonicContainer){
 389     AliError("Photonic background container not found");
 390     return NULL;
 391   }
 392
 393   Int_t stepbackground = 0;
 394   AliCFDataGrid *photonicGrid = new AliCFDataGrid("ContaminationGrid","ContaminationGrid",*photonicContainer,stepbackground);
 395
 396   // Subtract
 397   spectrumPhotonicSubtracted->Add(photonicGrid,-1.0);
 398
 399   // QA
 400   TH1D *photonicsubtractedspectrum = (TH1D *) spectrumPhotonicSubtracted->Project(0);
 401   AliHFEtools::NormaliseBinWidth(photonicsubtractedspectrum);
 402   TH1D *newrawspectrum = (TH1D *) dataSpectrumBeforePhotonicSubstraction->Project(0);
 403   AliHFEtools::NormaliseBinWidth(newrawspectrum);
 404   fQA->AddResultAt(photonicsubtractedspectrum,AliHFEInclusiveSpectrumQA::kAfterSPB);
 405   fQA->AddResultAt(newrawspectrum,AliHFEInclusiveSpectrumQA::kBeforeSPB);
 406   fQA->DrawSubtractPhotonicBackground();
 407
 408   return spectrumPhotonicSubtracted;
 409 }
 410
 411
 412 //____________________________________________________________
 413 AliCFDataGrid *AliHFEInclusiveSpectrum::CorrectParametrizedEfficiency(AliCFDataGrid* const bgsubpectrum){
 414
 415   //
 416   // Apply TPC pid efficiency correction from parametrisation
 417   //
 418
 419   // Data in the right format
 420   AliCFDataGrid *dataGrid = 0x0;
 421   if(bgsubpectrum) {
 422     dataGrid = bgsubpectrum;
 423   }
 424   else {
 425
 426     AliCFContainer *dataContainer = GetContainer(kDataContainer);
 427     if(!dataContainer){
 428       AliError("Data Container not available");
 429       return NULL;
 430     }
 431     dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
 432   }
 433   AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
 434   result->SetName("ParametrizedEfficiencyBefore");
 435   THnSparse *h = result->GetGrid();
 436   Int_t nbdimensions = h->GetNdimensions();
 437   //printf("CorrectParametrizedEfficiency::We have dimensions %d\n",nbdimensions);
 438   AliCFContainer *dataContainer = GetContainer(kDataContainer);
 439   if(!dataContainer){
 440     AliError("Data Container not available");
 441     return NULL;
 442   }
 443   AliCFContainer *dataContainerbis = (AliCFContainer *) dataContainer->Clone();
 444   dataContainerbis->Add(dataContainerbis,-1.0);
 445
 446
 447   Int_t* coord = new Int_t[nbdimensions];
 448   memset(coord, 0, sizeof(Int_t) * nbdimensions);
 449   Double_t* points = new Double_t[nbdimensions];
 450
 451   ULong64_t nEntries = h->GetNbins();
 452   for (ULong64_t i = 0; i < nEntries; ++i) {
 453
 454     Double_t value = h->GetBinContent(i, coord);
 455     //Double_t valuecontainer = dataContainerbis->GetBinContent(coord,fStepData);
 456     //printf("Value %f, and valuecontainer %f\n",value,valuecontainer);
 457
 458     // Get the bin co-ordinates given an coord
 459     for (Int_t j = 0; j < nbdimensions; ++j)
 460       points[j] = h->GetAxis(j)->GetBinCenter(coord[j]);
 461
 462     if (!fEfficiencyFunction->IsInside(points))
 463          continue;
 464     TF1::RejectPoint(kFALSE);
 465
 466     // Evaulate function at points
 467     Double_t valueEfficiency = fEfficiencyFunction->EvalPar(points, NULL);
 468     //printf("Value efficiency is %f\n",valueEfficiency);
 469
 470     if(valueEfficiency > 0.0) {
 471       h->SetBinContent(coord,value/valueEfficiency);
 472       dataContainerbis->SetBinContent(coord,fStepData,value/valueEfficiency);
 473       Double_t error = h->GetBinError(i);
 474       h->SetBinError(coord,error/valueEfficiency);
 475       dataContainerbis->SetBinError(coord,fStepData,error/valueEfficiency);
 476     }
 477
 478   }
 479
 480   delete[] coord;
 481   delete[] points;
 482
 483   AliCFDataGrid *resultt = new AliCFDataGrid("spectrumEfficiencyParametrized", "Data Grid for spectrum after Efficiency parametrized", *dataContainerbis,fStepData);
 484
 485  // QA
 486   TH1D *afterE = (TH1D *) resultt->Project(0);
 487   AliHFEtools::NormaliseBinWidth(afterE);
 488   TH1D *beforeE = (TH1D *) dataGrid->Project(0);
 489   AliHFEtools::NormaliseBinWidth(beforeE);
 490   fQA->AddResultAt(afterE,AliHFEInclusiveSpectrumQA::kAfterPE);
 491   fQA->AddResultAt(beforeE,AliHFEInclusiveSpectrumQA::kBeforePE);
 492   fQA->AddResultAt(fEfficiencyFunction,AliHFEInclusiveSpectrumQA::kPEfficiency);
 493   fQA->DrawCorrectWithEfficiency(AliHFEInclusiveSpectrumQA::kParametrized);
 494
 495   if(fNbDimensions>=2) {
 496     fQA->AddResultAt((TObject *) resultt,AliHFEInclusiveSpectrumQA::kAfterPEND);
 497     fQA->AddResultAt((TObject *) dataGrid,AliHFEInclusiveSpectrumQA::kBeforePEND);
 498     fQA->AddResultAt((TObject *) fEfficiencyFunction,AliHFEInclusiveSpectrumQA::kPEfficiencyND);
 499     fQA->DrawCorrectWithEfficiencyND(AliHFEInclusiveSpectrumQA::kParametrized);
 500   }
 501
 502   return resultt;
 503
 504 }
 505 //____________________________________________________________
 506 AliCFDataGrid *AliHFEInclusiveSpectrum::CorrectV0Efficiency(AliCFDataGrid* const bgsubpectrum){
 507
 508   //
 509   // Apply TPC pid efficiency correction from V0
 510   //
 511
 512   AliCFContainer *v0Container = GetContainer(kDataContainerV0);
 513   if(!v0Container){
 514     AliError("V0 Container not available");
 515     return NULL;
 516   }
 517
 518   // Efficiency
 519   AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*v0Container);
 520   efficiencyD->CalculateEfficiency(fStepAfterCutsV0,fStepBeforeCutsV0);
 521
 522   // Data in the right format
 523   AliCFDataGrid *dataGrid = 0x0;
 524   if(bgsubpectrum) {
 525     dataGrid = bgsubpectrum;
 526   }
 527   else {
 528     AliCFContainer *dataContainer = GetContainer(kDataContainer);
 529     if(!dataContainer){
 530       AliError("Data Container not available");
 531       return NULL;
 532     }
 533     dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
 534   }
 535
 536   // Correct
 537   AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
 538   result->ApplyEffCorrection(*efficiencyD);
 539
 540   // QA
 541   TH1D *afterE = (TH1D *) result->Project(0);
 542   AliHFEtools::NormaliseBinWidth(afterE);
 543   TH1D *beforeE = (TH1D *) dataGrid->Project(0);
 544   AliHFEtools::NormaliseBinWidth(beforeE);
 545   TH1D* efficiencyDproj = (TH1D *) efficiencyD->Project(0);
 546   fQA->AddResultAt(afterE,AliHFEInclusiveSpectrumQA::kAfterV0);
 547   fQA->AddResultAt(beforeE,AliHFEInclusiveSpectrumQA::kBeforeV0);
 548   fQA->AddResultAt(efficiencyDproj,AliHFEInclusiveSpectrumQA::kV0Efficiency);
 549   fQA->DrawCorrectWithEfficiency(AliHFEInclusiveSpectrumQA::kV0);
 550
 551   if(fNbDimensions>=2) {
 552     fQA->AddResultAt((TObject *) result,AliHFEInclusiveSpectrumQA::kAfterV0ND);
 553     fQA->AddResultAt((TObject *) dataGrid,AliHFEInclusiveSpectrumQA::kBeforeV0ND);
 554     fQA->AddResultAt((TObject *) efficiencyD,AliHFEInclusiveSpectrumQA::kV0EfficiencyND);
 555     fQA->DrawCorrectWithEfficiencyND(AliHFEInclusiveSpectrumQA::kV0);
 556   }
 557
 558
 559   return result;
 560
 561 }
 562 //____________________________________________________________
 563 THnSparse *AliHFEInclusiveSpectrum::Unfold(AliCFDataGrid* const bgsubpectrum){
 564
 565   //
 566   // Return the unfolded spectrum
 567   //
 568
 569   AliCFContainer *mcContainer = GetContainer(kMCContainerMC);
 570   if(!mcContainer){
 571     AliError("MC Container not available");
 572     return NULL;
 573   }
 574
 575   if(!fCorrelation){
 576     AliError("No Correlation map available");
 577     return NULL;
 578   }
 579
 580   // Data
 581   AliCFDataGrid *dataGrid = 0x0;
 582   if(bgsubpectrum) {
 583     dataGrid = bgsubpectrum;
 584   }
 585   else {
 586
 587     AliCFContainer *dataContainer = GetContainer(kDataContainer);
 588     if(!dataContainer){
 589       AliError("Data Container not available");
 590       return NULL;
 591     }
 592
 593     dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
 594   }
 595
 596   // Guessed
 597   AliCFDataGrid* guessedGrid = new AliCFDataGrid("guessed","",*mcContainer, fStepGuessedUnfolding);
 598   THnSparse* guessedTHnSparse = ((AliCFGridSparse*)guessedGrid->GetData())->GetGrid();
 599
 600   // Efficiency
 601   AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
 602   efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
 603
 604   // Unfold
 605
 606   AliCFUnfolding unfolding("unfolding","",fNbDimensions,fCorrelation,efficiencyD->GetGrid(),dataGrid->GetGrid(),guessedTHnSparse,1.e-06,0,fNumberOfIterations);
 607   if(fSetSmoothing) unfolding.UseSmoothing();
 608   unfolding.Unfold();
 609
 610   // Results
 611   THnSparse* result = unfolding.GetUnfolded();
 612   THnSparse* residual = unfolding.GetEstMeasured();
 613
 614   // QA
 615   TH1D *residualh = (TH1D *) residual->Projection(0);
 616   TH1D *beforeE = (TH1D *) dataGrid->Project(0);
 617   TH1D* efficiencyDproj = (TH1D *) efficiencyD->Project(0);
 618   TH1D *afterE = (TH1D *) result->Projection(0);
 619   AliHFEtools::NormaliseBinWidth(residualh);
 620   AliHFEtools::NormaliseBinWidth(beforeE);
 621   AliHFEtools::NormaliseBinWidth(afterE);
 622   fQA->AddResultAt(residualh,AliHFEInclusiveSpectrumQA::kResidualU);
 623   fQA->AddResultAt(afterE,AliHFEInclusiveSpectrumQA::kAfterU);
 624   fQA->AddResultAt(beforeE,AliHFEInclusiveSpectrumQA::kBeforeU);
 625   fQA->AddResultAt(efficiencyDproj,AliHFEInclusiveSpectrumQA::kUEfficiency);
 626   fQA->DrawUnfolding();
 627
 628   return (THnSparse *) result->Clone();
 629
 630 }
 631 //____________________________________________________________
 632 AliCFDataGrid *AliHFEInclusiveSpectrum::CorrectForEfficiency(AliCFDataGrid* const bgsubpectrum){
 633
 634   //
 635   // Apply unfolding and efficiency correction together to bgsubspectrum
 636   //
 637
 638   AliCFContainer *mcContainer = GetContainer(kMCContainerESD);
 639   if(!mcContainer){
 640     AliError("MC Container not available");
 641     return NULL;
 642   }
 643
 644   // Efficiency
 645   AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
 646   efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
 647
 648   // Data in the right format
 649   AliCFDataGrid *dataGrid = 0x0;
 650   if(bgsubpectrum) {
 651     dataGrid = bgsubpectrum;
 652   }
 653   else {
 654
 655     AliCFContainer *dataContainer = GetContainer(kDataContainer);
 656     if(!dataContainer){
 657       AliError("Data Container not available");
 658       return NULL;
 659     }
 660
 661     dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
 662   }
 663
 664   // Correct
 665   AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
 666   result->ApplyEffCorrection(*efficiencyD);
 667
 668   // QA
 669   TH1D *afterE = (TH1D *) result->Project(0);
 670   AliHFEtools::NormaliseBinWidth(afterE);
 671   TH1D *beforeE = (TH1D *) dataGrid->Project(0);
 672   AliHFEtools::NormaliseBinWidth(beforeE);
 673   TH1D* efficiencyDproj = (TH1D *) efficiencyD->Project(0);
 674   fQA->AddResultAt(afterE,AliHFEInclusiveSpectrumQA::kAfterMCE);
 675   fQA->AddResultAt(beforeE,AliHFEInclusiveSpectrumQA::kBeforeMCE);
 676   fQA->AddResultAt(efficiencyDproj,AliHFEInclusiveSpectrumQA::kMCEfficiency);
 677   fQA->DrawCorrectWithEfficiency(AliHFEInclusiveSpectrumQA::kMC);
 678
 679   if(fNbDimensions>=2) {
 680     fQA->AddResultAt((TObject *) result,AliHFEInclusiveSpectrumQA::kAfterMCEND);
 681     fQA->AddResultAt((TObject *) dataGrid,AliHFEInclusiveSpectrumQA::kBeforeMCEND);
 682     fQA->AddResultAt((TObject *) efficiencyD,AliHFEInclusiveSpectrumQA::kMCEfficiencyND);
 683     fQA->DrawCorrectWithEfficiencyND(AliHFEInclusiveSpectrumQA::kMC);
 684   }
 685
 686   return result;
 687
 688 }
 689 //____________________________________________________________
 690 void AliHFEInclusiveSpectrum::WriteResults(const char *filename)
 691 {
 692   //
 693   // Write results
 694   //
 695
 696   AliCFContainer *dataContainer = GetContainer(kDataContainer);
 697   AliCFContainer *mcContainer = GetContainer(kMCContainerMC);
 698   TObject *unfolded = 0x0;
 699   TObject *correctedspectrum = 0x0;
 700   if(fQA) {
 701     unfolded = fQA->GetResult(AliHFEInclusiveSpectrumQA::kFinalResultUnfolded);
 702     correctedspectrum = fQA->GetResult(AliHFEInclusiveSpectrumQA::kFinalResultDirectEfficiency);
 703   }
 704
 705   TFile *file = TFile::Open(filename,"recreate");
 706   if(dataContainer) dataContainer->Write("data");
 707   if(mcContainer) mcContainer->Write("mcefficiency");
 708   if(fCorrelation) fCorrelation->Write("correlationmatrix");
 709   if(unfolded) unfolded->Write("unfoldedspectrum");
 710   if(correctedspectrum) correctedspectrum->Write("correctedspectrum");
 711   if(fQA) fQA->Write("QAResults");
 712   file->Close();
 713
 714 }
 715