]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG2/RESONANCES/AliRsnFunction.cxx
Add new version of macros for RSN analysis
[u/mrichter/AliRoot.git] / PWG2 / RESONANCES / AliRsnFunction.cxx
index d417ea27c350244152a358dc0745a0afdda19537..99064783c6abb0089b33d3f0c79806a0bcb68516 100644 (file)
 // author: A. Pulvirenti             (email: alberto.pulvirenti@ct.infn.it)
 //
 
-#include <Riostream.h>
 #include <TString.h>
+#include <TAxis.h>
 
 #include "AliLog.h"
 
-#include "AliRsnDaughter.h"
-#include "AliRsnEvent.h"
-#include "AliRsnPairDef.h"
-#include "AliRsnCut.h"
+#include "AliRsnValue.h"
 
 #include "AliRsnFunction.h"
 
 ClassImp(AliRsnFunction)
 
 //________________________________________________________________________________________
-AliRsnFunction::AliRsnFunction() :
-  TNamed(),
-  fFcnType(kFcnTypes),
-  fPairDef(0x0),
-  fTrack(0x0),
-  fPair(0x0),
-  fEvent(0x0),
-  fHistogram(0x0)
+AliRsnFunction::AliRsnFunction(Bool_t useTH1) :
+   TObject(),
+   fAxisList("AliRsnValue", 0),
+   fUseTH1(useTH1),
+   fSize(0),
+   fH1(0x0),
+   fHSparse(0x0),
+   fValues(0)
 {
 //
-// Constructor for 1D functions.
-// Requires only the binning of the output function,
-// which is stored as 'main' histoDef in fHistoDef[0]
+// Constructor.
 //
-
-  fBinType[0] = kNoBins;
-  fBinType[1] = kNoBins;
-
-  fHistoDef[0] = 0x0;
-  fHistoDef[1] = 0x0;
-  fHistoDef[2] = 0x0;
-}
-
-//________________________________________________________________________________________
-AliRsnFunction::AliRsnFunction
-(EFcnType fcnType, AliRsnHistoDef *hd) :
-  TNamed(),
-  fFcnType(fcnType),
-  fPairDef(0x0),
-  fTrack(0x0),
-  fPair(0x0),
-  fEvent(0x0),
-  fHistogram(0x0)
-{
-//
-// Constructor for 1D functions.
-// Requires only the binning of the output function,
-// which is stored as 'main' histoDef in fHistoDef[0]
-//
-
-  fBinType[0] = kNoBins;
-  fBinType[1] = kNoBins;
-
-  fHistoDef[0] = hd;
-  fHistoDef[1] = 0x0;
-  fHistoDef[2] = 0x0;
-
-  DefineName();
 }
 
-//________________________________________________________________________________________
-AliRsnFunction::AliRsnFunction
-(EFcnType fcnType, EBinType binType, AliRsnHistoDef *hdMain, AliRsnHistoDef *hdBin) :
-  TNamed(),
-  fFcnType(fcnType),
-  fPairDef(0x0),
-  fTrack(0x0),
-  fPair(0x0),
-  fEvent(0x0),
-  fHistogram(0x0)
-{
-//
-// Constructor for 2D functions.
-// Requires the binning of the output function,
-// which is stored as 'main' histoDef in fHistoDef[0],
-// and a definition for a secondary binning, stored in fHistoDef[1]
-//
-
-  fBinType[0] = binType;
-  fBinType[1] = kNoBins;
-
-  fHistoDef[0] = hdMain;
-  fHistoDef[1] = hdBin;
-  fHistoDef[2] = 0x0;
-
-  DefineName();
-}
-
-//________________________________________________________________________________________
-AliRsnFunction::AliRsnFunction
-(EFcnType fcnType, EBinType binType1, EBinType binType2,
- AliRsnHistoDef *hdMain, AliRsnHistoDef *hdBin1, AliRsnHistoDef *hdBin2) :
-  fFcnType(fcnType),
-  fPairDef(0x0),
-  fTrack(0x0),
-  fPair(0x0),
-  fEvent(0x0),
-  fHistogram(0x0)
-{
-//
-// Constructor for 3D functions.
-// Requires the binning of the output function,
-// which is stored as 'main' histoDef in fHistoDef[0],
-// and a definition for two secondary binnings, stored in fHistoDef[1,2]
-//
-
-  fBinType[0] = binType1;
-  fBinType[1] = binType2;
-
-  fHistoDef[0] = hdMain;
-  fHistoDef[1] = hdBin1;
-  fHistoDef[2] = hdBin2;
-
-  DefineName();
-}    
-
 //________________________________________________________________________________________
 AliRsnFunction::AliRsnFunction(const AliRsnFunction &copy) :
-  TNamed(copy),
-  fFcnType(copy.fFcnType),
-  fPairDef(copy.fPairDef),
-  fTrack(copy.fTrack),
-  fPair(copy.fPair),
-  fEvent(copy.fEvent),
-  fHistogram(0x0)
+   TObject(copy),
+   fAxisList(copy.fAxisList),
+   fUseTH1(copy.fUseTH1),
+   fSize(copy.fSize),
+   fH1(0x0),
+   fHSparse(0x0),
+   fValues(copy.fValues)
 {
 //
 // Copy constructor.
 //
-
-  Int_t i;
-  for (i = 0; i < 3; i++) {
-    fHistoDef[i] = copy.fHistoDef[i];
-    if (i < 2) fBinType[i] = copy.fBinType[i];
-  }
-
-  DefineName();
 }
 
 //________________________________________________________________________________________
@@ -170,165 +67,62 @@ const AliRsnFunction& AliRsnFunction::operator=(const AliRsnFunction& copy)
 // Assignment operator.
 //
 
-  SetName(copy.GetName());
-  SetTitle(copy.GetTitle());
-
-  fFcnType = copy.fFcnType;
-
-  fPairDef = copy.fPairDef;
-
-  Int_t i;
-  for (i = 0; i < 3; i++) {
-    fHistoDef[i] = copy.fHistoDef[i];
-    if (i < 2) fBinType[i] = copy.fBinType[i];
-  }
-
-  fTrack = copy.fTrack;
-  fPair = copy.fPair;
-  fEvent = copy.fEvent;
-
-  if (fHistogram) delete fHistogram;
-  fHistogram = 0x0;
+   fAxisList = copy.fAxisList;
+   fUseTH1 = copy.fUseTH1;
+   fSize = copy.fSize;
+   fValues = copy.fValues;
 
-  DefineName();
+   if (fH1) delete fH1;
+   if (fHSparse) delete fHSparse;
 
-  return (*this);
+   return (*this);
 }
 
 //________________________________________________________________________________________
-const char* AliRsnFunction::FcnName()
+const char* AliRsnFunction::GetName() const
 {
 //
 // Defines the name of this object according to
 // the function type and binning
 //
 
-  switch (fFcnType)
-  {
-    case kTrackPt:
-      return  "TRKPT";
-      break;
-    case kTrackEta:
-      return  "TRKETA";
-      break;
-    case kInvMass:
-      return  "IM";
-      break;
-    case kInvMassMC:
-      return  "IMMC";
-      break;
-    case kResolution:
-      return  "RES";
-      break;
-    case kPairPt:
-      return  "PT";
-      break;
-    case kPairEta:
-      return  "ETA";
-      break;
-    case kEventMult:
-      return  "MULT";
-      break;
-    default:
-      return  "UNDEF";
-  }
-}
+   TString name("");
 
-//________________________________________________________________________________________
-void AliRsnFunction::DefineName()
-{
-//
-// Defines the name of this object according to
-// the function type and binning
-//
+   TObjArrayIter next(&fAxisList);
+   AliRsnValue *axis = 0;
 
-  Int_t  dim = CheckDim();
-
-  switch (dim)
-  {
-    case 1:
-      SetName(FcnName());
-      break;
-    case 2:
-      SetName(Form("%s_%s", FcnName(), BinName(fBinType[0])));
-      break;
-    case 3:
-      SetName(Form("%s_%s_%s", FcnName(), BinName(fBinType[0]), BinName(fBinType[1])));
-      break;
-    default:
-      SetName("UNDEF");
-  }
+   while ((axis = (AliRsnValue*)next())) {
+      if (name.Length() > 1) name += '_';
+      name += axis->GetName();
+   }
+
+   return name.Data();
 }
 
 //________________________________________________________________________________________
-Double_t AliRsnFunction::Eval()
+Bool_t AliRsnFunction::AddAxis(AliRsnValue *const axis)
 {
 //
-// Compute value for functions with 'event' argument type
+// Try to add a new axis to this function.
+// The operation succeeds only if the related value object
+// has a target compatible with the function type:
+// -- 'single'     functions, for tracks/events: AliRsnDaughter or AliRsnEvent
+// -- 'not single' functions, for pairs/events : AliRsnMother or AliRsnEvent
+// otherwise the axis is not added.
+//
+// If more than 3 axes are added, switch to THnSparseF output.
+// NOTE: this can cause large files and high memory occupancy.
 //
 
-  Double_t value;
-
-  switch (fFcnType)
-  {
-    case kTrackPt:
-      return fTrack->Pt();
-    case kTrackEta:
-      return fTrack->Eta();
-    case kInvMass:
-      return fPair->GetInvMass(fPairDef->GetMass(0), fPairDef->GetMass(1));
-    case kInvMassMC:
-      return fPair->GetInvMassMC(fPairDef->GetMass(0), fPairDef->GetMass(1));
-    case kResolution:
-      value  = fPair->GetInvMass(fPairDef->GetMass(0), fPairDef->GetMass(1));
-      value -= fPair->GetInvMassMC(fPairDef->GetMass(0), fPairDef->GetMass(1));
-      value /= fPair->GetInvMassMC(fPairDef->GetMass(0), fPairDef->GetMass(1));
-      return value;
-    case kPairPt:
-      return fPair->GetPt();
-    case kPairEta:
-      return fPair->GetEta();
-    case kEventMult:
-      return fEvent->GetMultiplicity();
-    default:
-      AliWarning("Function type not supported");
-      return -999.0;
-  }
-}
+   Int_t size = fAxisList.GetEntries();
+   new (fAxisList[size]) AliRsnValue(*axis);
 
-//_________________________________________________________________________________________________
-Bool_t AliRsnFunction::CheckInput(Option_t *option)
-{
-//
-// Checks if the argument type is coherent with
-// the function type required
-//
+   if (fAxisList.GetEntries() > 3) {
+      AliWarning("Adding more than 3 axes will produce a THnSparseD output.");
+      fUseTH1 = kFALSE;
+   }
 
-  TString opt(option);
-  opt.ToUpper();
-
-  if (opt.Contains("TRACK")) {
-    if (!fTrack) {
-      AliError("Input track object is NULL");
-      return kFALSE;
-    }
-  }
-
-  if (opt.Contains("PAIR")) {
-    if (!fPair) {
-      AliError("Input pair object is NULL");
-      return kFALSE;
-    }
-  }
-
-  if (opt.Contains("EVENT")) {
-    if (!fEvent) {
-      AliError("Input event object is NULL");
-      return kFALSE;
-    }
-  }
-
-  return kTRUE;
+   return kTRUE;
 }
 
 //________________________________________________________________________________________
@@ -341,174 +135,182 @@ TH1* AliRsnFunction::CreateHistogram(const char *histoName, const char *histoTit
 // even if the bins are equal, since they are defined in this class.
 // Eventually present histoDef's in other slots of array (1, 2) are ignored.
 //
-
-  // first binning is required
-  if (!fHistoDef[0]) return 0;
-
-  // retrieve binnings for main and secondary axes
-  Int_t    i, nbins[3] = {0, 0, 0};
-  Double_t min[3] = {0., 0., 0.}, max[3] = {0., 0., 0.};
-  for (i = 0; i < 3; i++)
-  {
-    if (fHistoDef[i])
-    {
-      nbins[i] = fHistoDef[i]->GetNBins();
-      min[i] = fHistoDef[i]->GetMin();
-      max[i] = fHistoDef[i]->GetMax();
-    }
-  }
-    
-  // define the kind of output according to the number of histoDefs
-  if (fHistogram) delete fHistogram;
-  if (!nbins[1] && !nbins[2]) {
-    fHistogram = new TH1D(histoName, histoTitle, nbins[0], min[0], max[0]);
-    fHistogram->SetXTitle(FcnName());
-  }
-  else if (nbins[1] > 0 && !nbins[2]) {
-    fHistogram = new TH2D(histoName, histoTitle, nbins[0], min[0], max[0], nbins[1], min[1], max[1]);
-    fHistogram->SetXTitle(FcnName());
-    fHistogram->SetYTitle(BinName(fBinType[0]));
-  }
-  else {
-    fHistogram = new TH3D(histoName, histoTitle, nbins[0], min[0], max[0], nbins[1], min[1], max[1], nbins[2], min[2], max[2]);
-    fHistogram->SetXTitle(FcnName());
-    fHistogram->SetYTitle(BinName(fBinType[0]));
-    fHistogram->SetZTitle(BinName(fBinType[1]));
-  }
-  
-  fHistogram->Sumw2();
-
-  return fHistogram;
+// This version produces a THnSparseF.
+//
+
+   fSize = fAxisList.GetEntries();
+   if (!fSize) {
+      AliError("No axes defined");
+      return 0x0;
+   } else if (fSize > 3) {
+      AliError("Too many axes defined for a TH1 object");
+      return 0x0;
+   }
+   
+   // initialize the size of values container
+   fValues.Set(fSize);
+
+   // retrieve binnings for main and secondary axes
+   AliRsnValue *fcnAxis;
+   TArrayD      array[3];
+   for (Int_t i = 0; i < fSize; i++) {
+      fcnAxis = (AliRsnValue*)fAxisList.At(i);
+      if (!fcnAxis) {
+         AliError("Empty axis");
+         array[i].Set(2);
+         array[i][0] = -1E5;
+         array[i][1] = -1E5;
+         continue;
+      } else {
+         array[i] = fcnAxis->GetArray();
+      }
+   }
+
+   // create histogram depending on the number of axes
+   switch (fSize) {
+      case 1:
+         fH1 = new TH1F(histoName, histoTitle, array[0].GetSize() - 1, array[0].GetArray());
+         break;
+      case 2:
+         fH1 = new TH2F(histoName, histoTitle, array[0].GetSize() - 1, array[0].GetArray(), array[1].GetSize() - 1, array[1].GetArray());
+         break;
+      case 3:
+         fH1 = new TH3F(histoName, histoTitle, array[0].GetSize() - 1, array[0].GetArray(), array[1].GetSize() - 1, array[1].GetArray(), array[2].GetSize() - 1, array[2].GetArray());
+         break;
+   }
+   fH1->Sumw2();
+
+   return fH1;
 }
 
 //________________________________________________________________________________________
-Bool_t AliRsnFunction::Fill()
+THnSparseF* AliRsnFunction::CreateHistogramSparse(const char *histoName, const char *histoTitle)
 {
 //
-// Fill function histogram with values computed from given input object.
+// Creates and returns the histogram defined using
+// arguments fo name and title, and the first histoDef for binning.
+// Variable-sized histogram binning is always called, due to use of histoDef,
+// even if the bins are equal, since they are defined in this class.
+// Eventually present histoDef's in other slots of array (1, 2) are ignored.
 //
-  AliDebug(AliLog::kDebug +2,"->");
-
-  // checks coherence between fcn type and passed argument
-  switch (fFcnType) {
-    case kTrackPt:
-    case kTrackEta:
-      if (!CheckInput("TRACK")) return kFALSE;
-      break;
-    case kInvMass:
-    case kInvMassMC:
-    case kResolution:
-    case kPairPt:
-    case kPairEta:
-      if (!CheckInput("PAIR")) return kFALSE;
-      break;
-    case kEventMult:
-      if (!CheckInput("EVENT")) return kFALSE;
-      break;
-    default:
-      AliError(Form("Input type %d not defined", (Int_t)fFcnType));
-      return kFALSE;
-  }
-
-  // check presence of output histogram
-  if (!fHistogram) {
-    AliError("Histogram is not yet initialized");
-    return kFALSE;
-  }
-
-  // compute value and stores into histogram
-  Int_t    dim = CheckDim();
-  Double_t mainValue, binValue[2];
-
-  TH1D *h1 = dynamic_cast<TH1D*>(fHistogram);
-  TH2D *h2 = dynamic_cast<TH2D*>(fHistogram);
-  TH3D *h3 = dynamic_cast<TH3D*>(fHistogram);
-  
-  mainValue = Eval();
-
-  switch (dim)
-  {
-    case 1:
-      if (h1) h1->Fill(mainValue);
-      break;
-    case 2:
-      binValue[0] = BinValue(fBinType[0]);
-      if (h2) h2->Fill(mainValue, binValue[0]);
-      break;
-    case 3:
-      binValue[0] = BinValue(fBinType[0]);
-      binValue[1] = BinValue(fBinType[1]);
-      if (h3) h3->Fill(mainValue, binValue[0], binValue[1]);
-      break;
-    default:
-      AliError("Wrong number of dimensions in the histogram. Check HD initialization");
-      return kFALSE;
-  }
-
-  AliDebug(AliLog::kDebug +2,"->");
-  return kTRUE;
+// This version produces a THnSparseF.
+//
+
+   fSize = fAxisList.GetEntries();
+   if (!fSize) {
+      AliError("No axes defined");
+      return 0x0;
+   }
+   
+   // initialize the size of values container
+   fValues.Set(fSize);
+
+   // initialize the array of number of bins for each axis
+   // taking it from the stored values, while for the bins
+   // they are set as summied and defined later
+   Int_t       *nbins   = new Int_t[fSize];
+   AliRsnValue *fcnAxis = 0;
+   for (Int_t i = 0; i < fSize; i++) {
+      fcnAxis = (AliRsnValue*)fAxisList.At(i);
+      if (!fcnAxis) {
+         nbins[i] = 1;
+         AliError("Empty axis");
+         continue;
+      }
+      nbins[i] = fcnAxis->GetArray().GetSize() - 1;
+   }
+
+   // create histogram
+   fHSparse = new THnSparseF(histoName, histoTitle, fSize, nbins);
+   fHSparse->Sumw2();
+
+   // update the various axes using the definitions given in the array of axes here
+   for (Int_t i = 0; i < fSize; i++) {
+      fcnAxis = (AliRsnValue*)fAxisList.At(i);
+      if (!fcnAxis) {
+         AliError("Empty axis: doing unique bin betweeen -100000 and 100000");
+         continue;
+      }
+      TAxis* axis = fHSparse->GetAxis(i);
+      axis->Set(nbins[i], fcnAxis->GetArray().GetArray());
+   }
+
+   delete [] nbins;
+
+   return fHSparse;
 }
 
-//________________________________________________________________________________________
-Double_t AliRsnFunction::BinValue(EBinType binType)
-{
-//
-// Computes the value for binning from the argument.
-// For each kind of binning type, the object is expected
-// to be of a given type, otherwise an error is raised.
-//
-
-  // checks coherence between bin type and passed argument
-  switch (binType) {
-    case kBinPairPt:
-      if (!CheckInput("PAIR")) return 0.0;
-      return fPair->GetPt();
-    case kBinPairEta:
-      if (!CheckInput("PAIR")) return 0.0;
-      return fPair->GetEta();
-    case kBinEventMult:
-      if (!CheckInput("EVENT")) return 0.0;
-      return fEvent->GetMultiplicity();
-    default:
-      AliError(Form("%s: Binning type not defined", GetName()));
-      return 0.0;
-  }
-}
 
 //________________________________________________________________________________________
-Int_t AliRsnFunction::CheckDim()
+Bool_t AliRsnFunction::Fill(TObject *object)
 {
 //
-// Checks number of dimensions.
-// Makes sure that eventual binnings are coherent and well defined
-//
-
-  if (!fHistoDef[0]) return 0;
-  if (fHistoDef[0] && !fHistoDef[1] && !fHistoDef[2]) return 1;
-  if (fHistoDef[0] && fHistoDef[1] && !fHistoDef[2]) return 2;
-
-  return 3;
+// Fill function histogram using the passed object
+//
+
+   AliDebug(AliLog::kDebug + 2, "->");
+
+   // loop on axes and try to compute values
+   // using this object or, as an alternative
+   // its reference event
+   Int_t  i;
+   Bool_t globalOK = kTRUE, computeOK;
+   AliRsnValue *fcnAxis = 0;
+   for (i = 0; i < fSize; i++) {
+      fValues[i] = 0.0;
+      computeOK = kFALSE;
+      fcnAxis = (AliRsnValue*)fAxisList.At(i);
+      if (fcnAxis) {
+         computeOK = fcnAxis->Eval(object);
+         if (computeOK) fValues[i] = ((Float_t)fcnAxis->GetComputedValue());
+      }
+      if (!computeOK) globalOK = kFALSE;
+   }
+   
+   // if even one of the computations has failes, the histograms are not filled
+   if (!globalOK) return kFALSE;
+
+   // fill histogram
+   if (fUseTH1) {
+      
+      // check presence of output histogram
+      if (!fH1) {
+         AliError("Required a TH1 which is not initialized");
+         return kFALSE;
+      }
+
+      // fill according to dimensions
+      switch (fSize) {
+         case 1: {
+            TH1F *h1 = (TH1F*)fH1;
+            h1->Fill(fValues[0]);
+         }
+         break;
+         case 2: {
+            TH2F *h2 = (TH2F*)fH1;
+            h2->Fill(fValues[0], fValues[1]);
+         }
+         break;
+         case 3: {
+            TH3F *h3 = (TH3F*)fH1;
+            h3->Fill(fValues[0], fValues[1], fValues[2]);
+         }
+         break;
+         default:
+            AliError(Form("Wrong size : %d", fSize));
+            return kFALSE;
+      }
+   } else {
+      
+      // check presence of output histogram
+      if (!fHSparse) {
+         AliError("Required a THnSparseF which is not initialized");
+         return kFALSE;
+      }
+
+      fHSparse->Fill(fValues.GetArray());
+   }
+
+   AliDebug(AliLog::kDebug + 2, "->");
+   return kTRUE;
 }
-
-//________________________________________________________________________________________
-const char* AliRsnFunction::BinName(EBinType binType)
-{
-//
-// Defines the name of binning
-//
-
-  switch (binType)
-  {
-    case kBinPairPt:
-      return "PT";
-      break;
-    case kBinPairEta:
-      return "ETA";
-      break;
-    case kBinEventMult:
-      return "MULT";
-      break;
-    default:
-      return "UNDEF";
-  }
-}
\ No newline at end of file