]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG2/RESONANCES/AliRsnFunction.cxx
fixed sig.segv
[u/mrichter/AliRoot.git] / PWG2 / RESONANCES / AliRsnFunction.cxx
index ac10aeddfc5158e637f5e7b46843669a151cca1c..bcbe014c660113c7ae16c73f2cc518327589bcfe 100644 (file)
 // author: A. Pulvirenti             (email: alberto.pulvirenti@ct.infn.it)
 //
 
-#include <Riostream.h>
-
-#include <TH1.h>
-#include <TList.h>
 #include <TString.h>
+#include <TAxis.h>
 
 #include "AliLog.h"
 
-#include "AliRsnMCInfo.h"
 #include "AliRsnDaughter.h"
 #include "AliRsnEvent.h"
 #include "AliRsnPairDef.h"
-#include "AliRsnCut.h"
+#include "AliRsnMother.h"
+#include "AliRsnValue.h"
 
 #include "AliRsnFunction.h"
 
 ClassImp(AliRsnFunction)
 
 //________________________________________________________________________________________
-AliRsnFunction::AliRsnFunction() :
-    fFcnType(kFcnTypes),
-    fRotAngle(0.0),
-    fUseBins(kFALSE),
-    fSkipOutsideInterval(kFALSE),
-    fBins(0),
-    fBinningCut(),
-    fBinningCutType(AliRsnCut::kLastCutType),
-    fHistoDef(0x0)
+AliRsnFunction::AliRsnFunction(Bool_t useTH1) :
+    TObject(),
+    fPairDef(0x0),
+    fAxisList("AliRsnValue", 0),
+    fPair(0x0),
+    fEvent(0x0),
+    fUseTH1(useTH1),
+    fSize(0),
+    fH1(0x0),
+    fHSparse(0x0)
 {
-  //
-  // Constructor.
-  // The histogram data member cannot be passed externally,
-  // its initialization MUST be defined inside the Init() method,
-  // which must be overridden in any derivate implementation.
-  //
-
-  Int_t i;
-  for (i = 0; i < 100; i++)
-  {
-    fHisto[i] = 0x0;
-  }
+//
+// Constructor.
+//
 }
 
 //________________________________________________________________________________________
-AliRsnFunction::AliRsnFunction
-(EFcnType type, AliRsnHistoDef *hd, Bool_t skipOut) :
-    fFcnType(type),
-    fRotAngle(0.0),
-    fUseBins(kFALSE),
-    fSkipOutsideInterval(skipOut),
-    fBins(0),
-    fBinningCut(),
-    fBinningCutType(AliRsnCut::kLastCutType),
-    fHistoDef(hd)
+AliRsnFunction::AliRsnFunction(const AliRsnFunction &copy) :
+    TObject(copy),
+    fPairDef(copy.fPairDef),
+    fAxisList(copy.fAxisList),
+    fPair(copy.fPair),
+    fEvent(copy.fEvent),
+    fUseTH1(copy.fUseTH1),
+    fSize(copy.fSize),
+    fH1(0x0),
+    fHSparse(0x0)
 {
-  //
-  // Constructor.
-  // The histogram data member cannot be passed externally,
-  // its initialization MUST be defined inside the Init() method,
-  // which must be overridden in any derivate implementation.
-  //
-
-  Int_t i;
-  for (i = 0; i < 100; i++)
-  {
-    fHisto[i] = 0x0;
-  }
+//
+// Copy constructor.
+//
 }
 
 //________________________________________________________________________________________
-AliRsnFunction::AliRsnFunction(const AliRsnFunction &copy) :
-    TObject(copy),
-    fFcnType(copy.fFcnType),
-    fRotAngle(copy.fRotAngle),
-    fUseBins(copy.fUseBins),
-    fSkipOutsideInterval(copy.fSkipOutsideInterval),
-    fBins(0),
-    fBinningCut(),
-    fBinningCutType(AliRsnCut::kLastCutType),
-    fHistoDef(copy.fHistoDef)
+const AliRsnFunction& AliRsnFunction::operator=(const AliRsnFunction& copy)
 {
-  //
-  // Copy constructor.
-  // Calls the function to define binning.
-  //
+//
+// Assignment operator.
+//
 
-  Int_t i, n = 100;
-  for (i = 0; i < n; i++)
-  {
-    fHisto[i] = 0x0;
-  }
+  //SetName(copy.GetName());
+  //SetTitle(copy.GetTitle());
 
-  if (fUseBins)
-  {
-    n = copy.fBins.GetSize();
-    Double_t *array = new Double_t[n];
-    for (i = 0; i < n; i++) array[i] = copy.fBins[i];
-    SetBinningCut(copy.fBinningCutType, copy.fBins.GetSize(), array);
-    delete [] array;
-  }
-}
-//________________________________________________________________________________________
-const AliRsnFunction& AliRsnFunction::operator=(const AliRsnFunction& /*copy*/)
-{
-  //
-  // Assignment operator.
-  // Behaves like copy constructor.
-  // Also in this case, the histogram is not copied, and,
-  // if it was present, it is destroyed and will need to be recreated.
-  //
+  fPairDef = copy.fPairDef;
+  fPair = copy.fPair;
+  fEvent = copy.fEvent;
+  fUseTH1 = copy.fUseTH1;
+  fSize = copy.fSize;
+
+  if (fH1) delete fH1;
+  fH1 = 0x0;
+  
+  if (fHSparse) delete fHSparse;
+  fHSparse = 0x0;
 
   return (*this);
 }
+
 //________________________________________________________________________________________
-void AliRsnFunction::Clear(Option_t* /*option*/)
+const char* AliRsnFunction::GetName() const
 {
-  //
-  // Clear arrays and histogram.
-  // For the sake of security, all pointers are also set explicitly to NULL.
-  //
+//
+// Defines the name of this object according to
+// the function type and binning
+//
 
-  Int_t i;
-  for (i = 0; i < 100; i++)
-  {
-    delete fHisto[i];
-    fHisto[i] = 0x0;
-  }
-}
+  TString name("");
 
-//________________________________________________________________________________________
-TList* AliRsnFunction::Init(const char *histoName, const char *histoTitle)
-{
-  //
-  // Initialization function.
-  // By default, it initializes the owned histogram using the method
-  // from AliRsnHistoDef class, giving the same name and title of this.
-  // A user can override this behaviour, if necessary.
-  // Before creating, the HistoDef is checked for proper initialization.
-  //
-
-  Clear();
-
-  Int_t i, ibin, nbins = fHistoDef->GetNBins();
-  Double_t min = fHistoDef->GetMin(), max = fHistoDef->GetMax();
-
-  // list is created and named after the general
-  // settings used for the contained histograms
-  TList *histos = new TList;
-  histos->SetName(Form("%s", GetFcnName().Data()));
-
-  // a general histogram is always added,
-  // which overrides the binning and collects everything
-  fHisto[0] = new TH1D(histoName, histoTitle, nbins, min, max);
-  histos->AddLast(fHisto[0]);
-
-  // if requested a binning w.r. to some cut variable, histograms are added
-  // for that in this part of the method (one per each bin)
-  Char_t hName[255];
-  Char_t hTitle[255];
-  if (fUseBins)
-  {
-    for (ibin = 0, i = 1; ibin < fBins.GetSize() - 1; ibin++, i++)
-    {
-      sprintf(hName, "%s[%.2f-%.2f]", histoName, fBins[ibin], fBins[ibin+1]);
-      sprintf(hTitle, "%s [%.2f-%.2f]", histoTitle, fBins[ibin], fBins[ibin+1]);
-      fHisto[i] = new TH1D(hName, hTitle, nbins, min, max);
-      histos->AddLast(fHisto[i]);
-    }
+  TObjArrayIter next(&fAxisList);
+  AliRsnValue *axis = 0;
+
+  while ((axis = (AliRsnValue*)next())) {
+    if (name.Length() > 1) name += '_';
+    name += axis->GetName();
   }
 
-  // returns the full list at the end
-  return histos;
+  return name.Data();
 }
 
 //________________________________________________________________________________________
-void AliRsnFunction::Init(const char *histoName, const char *histoTitle, TList *histos)
+void AliRsnFunction::AddAxis(AliRsnValue *const axis)
 {
-  //
-  // Initialization function.
-  // By default, it initializes the owned histogram using the method
-  // from AliRsnHistoDef class, giving the same name and title of this.
-  // A user can override this behaviour, if necessary.
-  // Before creating, the HistoDef is checked for proper initialization.
-  //
-
-  Clear();
-
-  Int_t i, ibin, nbins = fHistoDef->GetNBins();
-  Double_t min = fHistoDef->GetMin(), max = fHistoDef->GetMax();
-
-  // list is created and named after the general
-  // settings used for the contained histograms
-  if (!histos)
+  AliDebug(AliLog::kDebug+2,"<-");
+  Int_t size = fAxisList.GetEntries();
+  new(fAxisList[size]) AliRsnValue(*axis);
+  AliDebug(AliLog::kDebug+2,"->");
+  
+  if (fAxisList.GetEntries() > 3)
   {
-    AliError("NULL target list!");
-    return;
-  }
-
-  // a general histogram is always added,
-  // which overrides the binning and collects everything
-  fHisto[0] = new TH1D(histoName, histoTitle, nbins, min, max);
-  histos->AddLast(fHisto[0]);
-
-  // if requested a binning w.r. to some cut variable, histograms are added
-  // for that in this part of the method (one per each bin)
-  Char_t hName[255];
-  Char_t hTitle[255];
-  if (fUseBins)
-  {
-    for (ibin = 0, i = 1; ibin < fBins.GetSize() - 1; ibin++, i++)
-    {
-      sprintf(hName, "%s[%.2f-%.2f]", histoName, fBins[ibin], fBins[ibin+1]);
-      sprintf(hTitle, "%s [%.2f-%.2f]", histoTitle, fBins[ibin], fBins[ibin+1]);
-      fHisto[i] = new TH1D(hName, hTitle, nbins, min, max);
-      histos->AddLast(fHisto[i]);
-    }
+    AliWarning("A TH1-type output cannot add more than 3 axes: switching to THnSparse -- THIS COULD CAUSE VERY LARGE FILES!!!");
+    fUseTH1 = kFALSE;
   }
 }
 
 //________________________________________________________________________________________
-void AliRsnFunction::SetBinningCut
-(AliRsnCut::EType type, Double_t min, Double_t max, Double_t step)
+TH1* AliRsnFunction::CreateHistogram(const char *histoName, const char *histoTitle)
 {
-  //
-  // Set fixed bins
-  //
-
-  fUseBins = kTRUE;
+//
+// 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.
+//
+// This version produces a THnSparseF.
+//
 
-  Int_t i, nBins = (Int_t)((max - min) / step) + 1;
-  fBinningCutType = type;
-  fBins.Set(nBins);
-  for (i = 0; i < nBins; i++)
+  fSize = fAxisList.GetEntries();
+  if (!fSize) {
+    AliError("No axes defined");
+    return 0x0;
+  }
+  else if (fSize < 1 || fSize > 3)
   {
-    fBins[i] = min + (Double_t)i * step;
+    AliError("Too few or too many axes defined");
+    return 0x0;
   }
-}
-
-//________________________________________________________________________________________
-void AliRsnFunction::SetBinningCut
-(AliRsnCut::EType type, Int_t nbins, Double_t *bins)
-{
-  //
-  // Set variable bins
-  //
 
-  fUseBins = kTRUE;
-
-  Int_t i;
-  fBinningCutType = type;
-  fBins.Set(nbins);
-  for (i = 0; i < nbins; i++)
+  // retrieve binnings for main and secondary axes
+  AliRsnValue *fcnAxis;
+  TArrayD      array[3];
+  for (Int_t i = 0; i < fSize; i++) 
   {
-    fBins[i] = bins[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();
+    }
   }
-}
 
-//________________________________________________________________________________________
-TString AliRsnFunction::GetFcnName()
-{
-  //
-  // Return a string which names the function type
-  //
-
-  TString text("Undef");
-
-  switch (fFcnType)
+  // create histogram depending on the number of axes
+  switch (fSize)
   {
-    case kInvMass:
-      text = "IM";
+    case 1:
+      fH1 = new TH1F(histoName, histoTitle, array[0].GetSize() - 1, array[0].GetArray());
       break;
-    case kInvMassMC:
-      text = "IM_MC";
+    case 2:
+      fH1 = new TH2F(histoName, histoTitle, array[0].GetSize() - 1, array[0].GetArray(), array[1].GetSize() - 1, array[1].GetArray());
       break;
-    case kInvMassRotated:
-      text = Form("IMR%.2f", fRotAngle);
+    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;
-    case kResolution:
-      text = "RES";
-      break;
-    case kPtSpectrum:
-      text = "PT";
-    default:
-      AliError("Type not defined");
   }
+  fH1->Sumw2();
 
-  return text;
+  return fH1;
 }
 
 //________________________________________________________________________________________
-TString AliRsnFunction::GetFcnTitle()
+THnSparseF* AliRsnFunction::CreateHistogramSparse(const char *histoName, const char *histoTitle)
 {
-  //
-  // Return a string which names the function type
-  //
+//
+// 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.
+//
+// This version produces a THnSparseF.
+//
 
-  TString text("Undef");
+  fSize = fAxisList.GetEntries();
+  if (!fSize) {
+    AliError("No axes defined");
+    return 0x0;
+  }
+  
+  // 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
+  Double_t     dummyD;
+  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;
+  }
 
-  switch (fFcnType)
+  // create histogram
+  fHSparse = new THnSparseF(histoName, histoTitle, fSize, nbins, &dummyD, &dummyD);
+  fHSparse->Sumw2();
+  
+  // update the various axes using the definitions given in the array of axes here
+  for (Int_t i = 0; i < fSize; i++) 
   {
-    case kInvMass:
-      text = "Invariant mass";
-      break;
-    case kInvMassMC:
-      text = "Invariant mass (MC)";
-      break;
-    case kResolution:
-      text = "Resolution";
-      break;
-    case kPtSpectrum:
-      text = "p_{#perp} distribution";
-    default:
-      AliError("Type not defined");
+    fcnAxis = (AliRsnValue*)fAxisList.At(i);
+    if (!fcnAxis) {
+      AliError("Empty axis: doing unique bin betweeen -100000 and 100000");
+      continue;
+    }
+    fHSparse->SetBinEdges(i, fcnAxis->GetArray().GetArray());
   }
+  
+  delete [] nbins;
 
-  return text;
+  return fHSparse;
 }
 
+
 //________________________________________________________________________________________
-Bool_t AliRsnFunction::Fill(AliRsnPairParticle *pair, AliRsnPairDef *ref, Double_t weight)
+Bool_t AliRsnFunction::Fill()
 {
-  //
-  // Fillse the histogram with data contained in a defined pair.
-  // This method must be overidden by an appropriate definition in each inheriting class.
-  //
+//
+// Fill function histogram with values computed from given input object.
+//
 
-  Double_t value = FcnValue(pair, ref);
-  if (fSkipOutsideInterval)
-  {
-    if (value < fHistoDef->GetMin()) return kFALSE;
-    if (value > fHistoDef->GetMax()) return kFALSE;
-  }
+  AliDebug(AliLog::kDebug +2,"->");
 
-  // fill global histogram
-  if (weight == 0.0) fHisto[0]->Fill(value);
-  else fHisto[0]->Fill(value, weight);
+  Int_t  i;
+  Double_t *values = new Double_t[fSize];
 
-  // if bins are allocated, find right one and fill it
-  if (fUseBins)
+  AliRsnValue *fcnAxis = 0;
+  for (i = 0; i < fSize; i++) {
+    fcnAxis = (AliRsnValue*)fAxisList.At(i);
+    if (!fcnAxis) 
+    {
+      values[i] = 0.0;
+      continue;
+    }
+    if (fcnAxis->Eval(fPair, fPairDef, fEvent)) values[i] = (Double_t)((Float_t)fcnAxis->GetValue());
+  }
+  
+  // fill histogram
+  if (fUseTH1)
   {
-    Int_t i, ibin;
-    for (ibin = 0, i = 1; ibin < fBins.GetSize() - 1; ibin++, i++)
+    // check presence of output histogram
+    if (!fH1) {
+      AliError("Required a TH1 which is not initialized");
+      delete [] values;
+      return kFALSE;
+    }
+    
+    // fill according to dimensions
+    switch (fSize)
     {
-      if (!fHisto[i]) continue;
-      fBinningCut.SetCutValues(fBinningCutType, fBins[ibin], fBins[ibin+1]);
-      if (fBinningCut.IsSelected(AliRsnCut::kPair, pair))
-      {
-        if (weight == 0.0) fHisto[i]->Fill(value);
-        else fHisto[i]->Fill(value, weight);
+      case 1:
+        {
+          TH1F *h1 = (TH1F*)fH1;
+          h1->Fill(values[0]);
+        }
+        break;
+      case 2:
+        {
+          TH2F *h2 = (TH2F*)fH1;
+          h2->Fill(values[0], values[1]);
+        }
+        break;
+      case 3:
+        {
+          TH3F *h3 = (TH3F*)fH1;
+          h3->Fill(values[0], values[1], values[2]);
+        }
         break;
-      }
+      default:
+        AliError(Form("Wrong size : %d", fSize));
+        delete [] values;
+        return kFALSE;
     }
   }
-
-  return kTRUE;
-}
-
-//________________________________________________________________________________________
-Double_t AliRsnFunction::FcnValue(AliRsnPairParticle *pair, AliRsnPairDef *ref)
-{
-  //
-  // This method must be overridden in all inheritin functions.
-  // It computes the value which must be used to fill the histogram.
-  //
-
-  switch (fFcnType)
+  else
   {
-    case kInvMass:
-      return pair->GetInvMass(ref->GetMass(0), ref->GetMass(1));
-    case kInvMassMC:
-      return pair->GetInvMassMC(ref->GetMass(0), ref->GetMass(1));
-    case kInvMassRotated:
-      //AliInfo(Form("*** ROTATION ANGLE = %f ***", fRotAngle));
-      //AliInfo(Form("UNROTATED INV MASS = %f", pair->GetInvMass(ref->GetMass(0), ref->GetMass(1))));
-      //pair->GetDaughter(1)->Print("P");
-      pair->GetDaughter(1)->RotateP(fRotAngle * TMath::DegToRad());
-      pair->ResetPair();
-      //AliInfo(Form("  ROTATED INV MASS = %f", pair->GetInvMass(ref->GetMass(0), ref->GetMass(1))));
-      //pair->GetDaughter(1)->Print("P");
-      return pair->GetInvMass(ref->GetMass(0), ref->GetMass(1));
-    case kResolution:
-      return FcnResolution(pair, ref);
-    case kPtSpectrum:
-      return pair->GetPt();
-    default:
-      AliError("Type not defined");
+    // check presence of output histogram
+    if (!fHSparse) {
+      AliError("Required a THnSparseF which is not initialized");
+      delete [] values;
+      return kFALSE;
+    }
+    
+    fHSparse->Fill(values);
   }
+  
+  delete [] values;
 
-  return 0.0;
-}
-
-//________________________________________________________________________________________
-inline Double_t AliRsnFunction::FcnResolution(AliRsnPairParticle *pair, AliRsnPairDef *ref)
-{
-  //
-  // Invariant mass resolution (compared between reconstructed and montecarlo)
-  //
-
-  Double_t recInvMass = pair->GetInvMass(ref->GetMass(0), ref->GetMass(1));
-  Double_t simInvMass = pair->GetInvMassMC(ref->GetMass(0), ref->GetMass(1));
-
-  return (simInvMass - recInvMass) / simInvMass;
+  AliDebug(AliLog::kDebug +2,"->");
+  return kTRUE;
 }