// author: A. Pulvirenti (email: alberto.pulvirenti@ct.infn.it)
//
-#include <Riostream.h>
-
-#include <TH1.h>
-#include <TList.h>
#include <TString.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),
- fNumberOfBinTypes(0),
-// 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, j;
- for (j = 0 ; j < kFcnBinTypes; j++)
- for (i = 0; i < 100; i++)
- fHisto[j][i] = 0x0;
-
+//
+// Constructor.
+//
}
//________________________________________________________________________________________
-AliRsnFunction::AliRsnFunction
-(EFcnType type, AliRsnHistoDef *hd, Bool_t skipOut) :
- fFcnType(type),
- fRotAngle(0.0),
- fUseBins(kFALSE),
- fSkipOutsideInterval(skipOut),
- fNumberOfBinTypes(0),
-// fBins(0),
-// fBinningCut(),
-// fBinningCutType(AliRsnCut::kLastCutType),
- fHistoDef(hd)
+AliRsnFunction::AliRsnFunction(const AliRsnFunction ©) :
+ 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, j;
- for (j = 0 ; j < kFcnBinTypes; j++)
- for (i = 0; i < 100; i++)
- fHisto[j][i] = 0x0;
+//
+// Copy constructor.
+//
}
//________________________________________________________________________________________
-AliRsnFunction::AliRsnFunction(const AliRsnFunction ©) :
- TObject(copy),
- fFcnType(copy.fFcnType),
- fRotAngle(copy.fRotAngle),
- fUseBins(copy.fUseBins),
- fSkipOutsideInterval(copy.fSkipOutsideInterval),
- fNumberOfBinTypes(copy.fNumberOfBinTypes),
-// 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, j, n;
- for (j = 0 ; j < kFcnBinTypes; j++)
- for (i = 0; i < 100; i++)
- fHisto[j][i] = 0x0;
+ //SetName(copy.GetName());
+ //SetTitle(copy.GetTitle());
- if (fUseBins)
- {
- for (i = 0 ; i < kFcnBinTypes; i++){
- if (fNumberOfBinTypes<=i) continue;
- n = copy.fBins[i].GetSize();
- Double_t *array = new Double_t[n];
- for (j = 0; j < n; j++) array[j] = copy.fBins[i][j];
- SetBinningCut(copy.fBinningCutType[i], copy.fBins[i].GetSize(), array,i,kTRUE);
- 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*/)
-{
- //
- // Clear arrays and histogram.
- // For the sake of security, all pointers are also set explicitly to NULL.
- //
-
- Int_t i, j;
- for (j = 0 ; j < kFcnBinTypes; j++)
- for (i = 0; i < 100; i++)
- {
- delete fHisto[j][i];
- fHisto[j][i] = 0x0;
- }
-}
//________________________________________________________________________________________
-TList* AliRsnFunction::Init(const char *histoName, const char *histoTitle)
+const char* AliRsnFunction::GetName() const
{
- //
- // 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][0] = new TH1D(histoName, histoTitle, nbins, min, max);
- fHisto[0][0]->Sumw2();
- histos->AddLast(fHisto[0][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];
- Int_t j;
- if (fUseBins)
- {
- for (j = 0 ; j < kFcnBinTypes; j++){
- if (fNumberOfBinTypes<=j) continue;
- for (ibin = 0, i = 1; ibin < fBins[j].GetSize() - 1; ibin++, i++)
- {
- sprintf(hName, "%s_%d%02d_[%.2f-%.2f]", histoName, j,i,fBins[j][ibin], fBins[j][ibin+1]);
- sprintf(hTitle, "%s [%.2f-%.2f]", histoTitle, fBins[j][ibin], fBins[j][ibin+1]);
-// AliInfo(Form("Adding %s",hName));
- fHisto[j][i] = new TH1D(hName, hTitle, nbins, min, max);
- fHisto[j][i]->Sumw2();
- histos->AddLast(fHisto[j][i]);
- }
- }
- }
+//
+// Defines the name of this object according to
+// the function type and binning
+//
- // returns the full list at the end
- return histos;
-}
+ TString name("");
-//________________________________________________________________________________________
-void AliRsnFunction::Init(const char *histoName, const char *histoTitle, TList *histos)
-{
- //
- // 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)
- {
- AliError("NULL target list!");
- return;
- }
+ TObjArrayIter next(&fAxisList);
+ AliRsnValue *axis = 0;
- // a general histogram is always added,
- // which overrides the binning and collects everything
- fHisto[0][0] = new TH1D(histoName, histoTitle, nbins, min, max);
- histos->AddLast(fHisto[0][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];
- Int_t j;
- if (fUseBins)
- {
- for (j = 0 ; j < kFcnBinTypes; j++){
- if (fNumberOfBinTypes<=j) continue;
- for (ibin = 0, i = 1; ibin < fBins[j].GetSize() - 1; ibin++, i++)
- {
-
- sprintf(hName, "%s_%d_%02d[%.2f-%.2f]", histoName,j,i, fBins[j][ibin], fBins[j][ibin+1]);
- sprintf(hTitle, "%s [%.2f-%.2f]", histoTitle, fBins[j][ibin], fBins[j][ibin+1]);
-// AliInfo(Form("Adding %s",hName));
- fHisto[j][i] = new TH1D(hName, hTitle, nbins, min, max);
- histos->AddLast(fHisto[j][i]);
- }
- }
+ while ((axis = (AliRsnValue*)next())) {
+ if (name.Length() > 1) name += '_';
+ name += axis->GetName();
}
+
+ return name.Data();
}
//________________________________________________________________________________________
-void AliRsnFunction::SetBinningCut
-(AliRsnCut::EType type, Double_t min, Double_t max, Double_t step,Int_t index,Bool_t IsCopyConstructor)
+void AliRsnFunction::AddAxis(AliRsnValue *const axis)
{
- //
- // Set fixed bins
- //
-
- if (index >= kFcnBinTypes) {
- AliError(Form("We support only %d Binning cuts(0-%d). Skipping...",kFcnBinTypes,kFcnBinTypes-1));
- return;
- }
+ AliDebug(AliLog::kDebug+2,"<-");
+ Int_t size = fAxisList.GetEntries();
+ new(fAxisList[size]) AliRsnValue(*axis);
+ AliDebug(AliLog::kDebug+2,"->");
- if (!IsCopyConstructor){
- // TODO if some one sets indexes 0,2,3 it is a bug here(i'll solve it)
- if (index == fNumberOfBinTypes)
- fNumberOfBinTypes++;
- else {
- AliError(Form("Wrong index %d. fUseBins is set to kFALSE",index));
-// fUseBins = kFALSE;
- return;
- }
- }
-
- fUseBins = kTRUE;
- Int_t i, nBins = (Int_t)((max - min) / step) + 1;
- fBinningCutType[index] = type;
- fBins[index].Set(nBins);
- for (i = 0; i < nBins; i++)
+ if (fAxisList.GetEntries() > 3)
{
- fBins[index][i] = min + (Double_t)i * step;
+ 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, Int_t nbins, Double_t *bins,Int_t index,Bool_t IsCopyConstructor)
+TH1* AliRsnFunction::CreateHistogram(const char *histoName, const char *histoTitle)
{
- //
- // Set variable bins
- //
+//
+// 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.
+//
- if (index >= kFcnBinTypes) {
- AliError(Form("We support only %d Binning cuts(0-%d). Skipping...",kFcnBinTypes,kFcnBinTypes-1));
- return;
+ fSize = fAxisList.GetEntries();
+ if (!fSize) {
+ AliError("No axes defined");
+ return 0x0;
}
- if (!IsCopyConstructor){
- // TODO if some one sets indexes 0,2,3 it is a bug here(i'll solve it)
- if (index >= fNumberOfBinTypes)
- fNumberOfBinTypes++;
- else {
- AliError(Form("Wrong index %d (%d). fUseBins is set to kFALSE",index,fNumberOfBinTypes));
- // fUseBins = kFALSE;
- return;
- }
- }
-
- fUseBins = kTRUE;
- Int_t i;
- fBinningCutType[index] = type;
- fBins[index].Set(nbins);
- for (i = 0; i < nbins; i++)
+ else if (fSize < 1 || fSize > 3)
{
- fBins[index][i] = bins[i];
+ AliError("Too few or too many axes defined");
+ return 0x0;
}
-}
-
-//________________________________________________________________________________________
-TString AliRsnFunction::GetFcnName()
-{
- //
- // Return a string which names the function type
- //
- TString text("Undef");
+ Int_t *nbins = new Int_t [fSize];
+ Double_t *min = new Double_t[fSize];
+ Double_t *max = new Double_t[fSize];
+
+ // retrieve binnings for main and secondary axes
+ AliRsnValue *fcnAxis = 0;
+ for (Int_t i = 0; i < fSize; i++) {
+ fcnAxis = (AliRsnValue*)fAxisList.At(i);
+ if (!fcnAxis) {
+ nbins[i] = 0;
+ min[i] = 0.0;
+ max[i] = 0.0;
+ AliError("Empty axis");
+ continue;
+ }
+ nbins[i] = fcnAxis->GetNBins();
+ min[i] = fcnAxis->GetMin();
+ max[i] = fcnAxis->GetMax();
+ }
- switch (fFcnType)
+ // create histogram depending on the number of axes
+ switch (fSize)
{
- case kInvMass:
- text = "IM";
- break;
- case kInvMassMC:
- text = "IM_MC";
- break;
- case kInvMassRotated:
- text = Form("IMR%.2f", fRotAngle);
+ case 1:
+ fH1 = new TH1F(histoName, histoTitle, nbins[0], min[0], max[0]);
break;
- case kResolution:
- text = "RES";
+ case 2:
+ fH1 = new TH2F(histoName, histoTitle, nbins[0], min[0], max[0], nbins[1], min[1], max[1]);
break;
- case kPtSpectrum:
- text = "PT";
+ case 3:
+ fH1 = new TH3F(histoName, histoTitle, nbins[0], min[0], max[0], nbins[1], min[1], max[1], nbins[2], min[2], max[2]);
break;
- case kEtaSpectrum:
- text = "ETA";
- break;
- 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;
+ }
- switch (fFcnType)
- {
- case kInvMass:
- text = "Invariant mass";
- break;
- case kInvMassMC:
- text = "Invariant mass (MC)";
- break;
- case kResolution:
- text = "Resolution";
- break;
- case kPtSpectrum:
- text = "p_{#perp} distribution";
- break;
- case kEtaSpectrum:
- text = "#eta distribution";
- break;
- default:
- AliError("Type not defined");
+ Int_t *nbins = new Int_t [fSize];
+ Double_t *min = new Double_t[fSize];
+ Double_t *max = new Double_t[fSize];
+
+ // retrieve binnings for main and secondary axes
+ AliRsnValue *fcnAxis = 0;
+ for (Int_t i = 0; i < fSize; i++) {
+ fcnAxis = (AliRsnValue*)fAxisList.At(i);
+ if (!fcnAxis) {
+ nbins[i] = 0;
+ min[i] = 0.0;
+ max[i] = 0.0;
+ AliError("Empty axis");
+ continue;
+ }
+ nbins[i] = fcnAxis->GetNBins();
+ min[i] = fcnAxis->GetMin();
+ max[i] = fcnAxis->GetMax();
}
- return text;
+ Int_t size = fAxisList.GetEntries();
+ if (!size) {
+ AliError("No axes defined");
+ return 0x0;
+ }
+
+ // create histogram
+ fHSparse = new THnSparseF(histoName, histoTitle, size, nbins, min, max);
+ fHSparse->Sumw2();
+
+ // clean heap
+ delete [] nbins;
+ delete [] min;
+ delete [] max;
+
+ return fHSparse;
}
+
//________________________________________________________________________________________
-Bool_t AliRsnFunction::Fill(AliRsnPairParticle *pair, AliRsnPairDef *ref)
+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
- fHisto[0][0]->Fill(value);
+ 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] = fcnAxis->GetValue();
+ }
+
+ // fill histogram
+ if (fUseTH1)
{
- Int_t i, j, ibin;
- for (j = 0 ; j < kFcnBinTypes; j++){
- if (fNumberOfBinTypes<=j) continue;
- for (ibin = 0, i = 1; ibin < fBins[j].GetSize() - 1; ibin++, i++)
- {
- if (!fHisto[j][i]) continue;
- fBinningCut[j].SetCutValues(fBinningCutType[j], fBins[j][ibin], fBins[j][ibin+1]);
- if (fBinningCut[j].IsSelected(AliRsnCut::kPair, pair))
- {
- fHisto[j][i]->Fill(value);
- break;
- }
- }
+ // check presence of output histogram
+ if (!fH1) {
+ AliError("Required a TH1 whish is not initialized");
+ return kFALSE;
+ }
+
+ // fill according to dimensions
+ switch (fSize)
+ {
+ 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));
+ 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();
- case kEtaSpectrum:
- return pair->GetEta();
- default:
- AliError("Type not defined");
+ // check presence of output histogram
+ if (!fHSparse) {
+ AliError("Required a THnSparseF which is not initialized");
+ 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;
}