PWG2/RESONANCES/AliRsnFunction.cxx

   1 //
   2 // Class AliRsnFunction
   3 //
   4 // This class defines a base classe to implement a function
   5 // which uses the internal RSN package event format (AliRsnEvent).
   6 // It contains some default flags which turn out to be useful:
   7 //  - a flag to select only the "true" pairs (tracks from same resonance)
   8 //  - a flag to know if the computation is done over two events (mixing)
   9 //
  10 // Any kind of analysis object should be implemented as inheriting from this
  11 // because the AliRsnAnalyzer which executes the analysis will accept a collection
  12 // of such objects, in order to have a unique format of processing method
  13 //
  14 // The user who implements a kind of computation type should inherit from
  15 // this class and override the virtual functions defined in it, which
  16 // initialize the final output histogram and define how to process data.
  17 //
  18 //
  19 // author: A. Pulvirenti             (email: alberto.pulvirenti@ct.infn.it)
  20 //
  21
  22 #include <TString.h>
  23 #include <TAxis.h>
  24
  25 #include "AliLog.h"
  26
  27 #include "AliRsnDaughter.h"
  28 #include "AliRsnEvent.h"
  29 #include "AliRsnPairDef.h"
  30 #include "AliRsnMother.h"
  31 #include "AliRsnValue.h"
  32
  33 #include "AliRsnFunction.h"
  34
  35 ClassImp(AliRsnFunction)
  36
  37 //________________________________________________________________________________________
  38 AliRsnFunction::AliRsnFunction(Bool_t useTH1) :
  39    TObject(),
  40    fPairDef(0x0),
  41    fAxisList("AliRsnValue", 0),
  42    fPair(0x0),
  43    fEvent(0x0),
  44    fUseTH1(useTH1),
  45    fSize(0),
  46    fH1(0x0),
  47    fHSparse(0x0)
  48 {
  49 //
  50 // Constructor.
  51 //
  52 }
  53
  54 //________________________________________________________________________________________
  55 AliRsnFunction::AliRsnFunction(const AliRsnFunction &copy) :
  56    TObject(copy),
  57    fPairDef(copy.fPairDef),
  58    fAxisList(copy.fAxisList),
  59    fPair(copy.fPair),
  60    fEvent(copy.fEvent),
  61    fUseTH1(copy.fUseTH1),
  62    fSize(copy.fSize),
  63    fH1(0x0),
  64    fHSparse(0x0)
  65 {
  66 //
  67 // Copy constructor.
  68 //
  69 }
  70
  71 //________________________________________________________________________________________
  72 const AliRsnFunction& AliRsnFunction::operator=(const AliRsnFunction& copy)
  73 {
  74 //
  75 // Assignment operator.
  76 //
  77
  78    fPairDef = copy.fPairDef;
  79    fPair = copy.fPair;
  80    fEvent = copy.fEvent;
  81    fUseTH1 = copy.fUseTH1;
  82    fSize = copy.fSize;
  83
  84    if (fH1) delete fH1;
  85    fH1 = 0x0;
  86
  87    if (fHSparse) delete fHSparse;
  88    fHSparse = 0x0;
  89
  90    return (*this);
  91 }
  92
  93 //________________________________________________________________________________________
  94 const char* AliRsnFunction::GetName() const
  95 {
  96 //
  97 // Defines the name of this object according to
  98 // the function type and binning
  99 //
 100
 101    TString name("");
 102
 103    TObjArrayIter next(&fAxisList);
 104    AliRsnValue *axis = 0;
 105
 106    while ((axis = (AliRsnValue*)next())) {
 107       if (name.Length() > 1) name += '_';
 108       name += axis->GetName();
 109    }
 110
 111    return name.Data();
 112 }
 113
 114 //________________________________________________________________________________________
 115 Bool_t AliRsnFunction::AddAxis(AliRsnValue *const axis)
 116 {
 117 //
 118 // Try to add a new axis to this function.
 119 // The operation succeeds only if the related value object
 120 // has a target amon those allowed here (AliRsnMother, AliRsnEvent),
 121 // otherwise the axis is not added.
 122 //
 123 // If more than 3 axes are added, switch to THnSparseF output.
 124 // NOTE: this can cause large files and high memory occupancy.
 125 //
 126
 127    RSNTARGET target = axis->GetTargetType();
 128    if (target != AliRsnTarget::kMother && target != AliRsnTarget::kEvent) {
 129       AliError(Form("Allowed targets are mothers and events; cannot use axis '%s' which has target '%s'", axis->GetName(), axis->GetTargetTypeName()));
 130       return kFALSE;
 131    }
 132
 133    Int_t size = fAxisList.GetEntries();
 134    new(fAxisList[size]) AliRsnValue(*axis);
 135
 136    if (fAxisList.GetEntries() > 3) {
 137       AliWarning("Adding more than 3 axes will produce a THnSparseD output.");
 138       fUseTH1 = kFALSE;
 139    }
 140
 141    return kTRUE;
 142 }
 143
 144 //________________________________________________________________________________________
 145 TH1* AliRsnFunction::CreateHistogram(const char *histoName, const char *histoTitle)
 146 {
 147 //
 148 // Creates and returns the histogram defined using
 149 // arguments fo name and title, and the first histoDef for binning.
 150 // Variable-sized histogram binning is always called, due to use of histoDef,
 151 // even if the bins are equal, since they are defined in this class.
 152 // Eventually present histoDef's in other slots of array (1, 2) are ignored.
 153 //
 154 // This version produces a THnSparseF.
 155 //
 156
 157    fSize = fAxisList.GetEntries();
 158    if (!fSize) {
 159       AliError("No axes defined");
 160       return 0x0;
 161    } else if (fSize > 3) {
 162       AliError("Too many axes defined for a TH1 object");
 163       return 0x0;
 164    }
 165
 166    // retrieve binnings for main and secondary axes
 167    AliRsnValue *fcnAxis;
 168    TArrayD      array[3];
 169    for (Int_t i = 0; i < fSize; i++) {
 170       fcnAxis = (AliRsnValue*)fAxisList.At(i);
 171       if (!fcnAxis) {
 172          AliError("Empty axis");
 173          array[i].Set(2);
 174          array[i][0] = -1E5;
 175          array[i][1] = -1E5;
 176          continue;
 177       } else {
 178          array[i] = fcnAxis->GetArray();
 179       }
 180    }
 181
 182    // create histogram depending on the number of axes
 183    switch (fSize) {
 184       case 1:
 185          fH1 = new TH1F(histoName, histoTitle, array[0].GetSize() - 1, array[0].GetArray());
 186          break;
 187       case 2:
 188          fH1 = new TH2F(histoName, histoTitle, array[0].GetSize() - 1, array[0].GetArray(), array[1].GetSize() - 1, array[1].GetArray());
 189          break;
 190       case 3:
 191          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());
 192          break;
 193    }
 194    fH1->Sumw2();
 195
 196    return fH1;
 197 }
 198
 199 //________________________________________________________________________________________
 200 THnSparseF* AliRsnFunction::CreateHistogramSparse(const char *histoName, const char *histoTitle)
 201 {
 202 //
 203 // Creates and returns the histogram defined using
 204 // arguments fo name and title, and the first histoDef for binning.
 205 // Variable-sized histogram binning is always called, due to use of histoDef,
 206 // even if the bins are equal, since they are defined in this class.
 207 // Eventually present histoDef's in other slots of array (1, 2) are ignored.
 208 //
 209 // This version produces a THnSparseF.
 210 //
 211
 212    fSize = fAxisList.GetEntries();
 213    if (!fSize) {
 214       AliError("No axes defined");
 215       return 0x0;
 216    }
 217
 218    // initialize the array of number of bins for each axis
 219    // taking it from the stored values, while for the bins
 220    // they are set as summied and defined later
 221    Double_t     dummyD;
 222    Int_t       *nbins   = new Int_t[fSize];
 223    AliRsnValue *fcnAxis = 0;
 224    for (Int_t i = 0; i < fSize; i++) {
 225       fcnAxis = (AliRsnValue*)fAxisList.At(i);
 226       if (!fcnAxis) {
 227          nbins[i] = 1;
 228          AliError("Empty axis");
 229          continue;
 230       }
 231       nbins[i] = fcnAxis->GetArray().GetSize() - 1;
 232    }
 233
 234    // create histogram
 235    fHSparse = new THnSparseF(histoName, histoTitle, fSize, nbins, &dummyD, &dummyD);
 236    fHSparse->Sumw2();
 237
 238    // update the various axes using the definitions given in the array of axes here
 239    for (Int_t i = 0; i < fSize; i++) {
 240       fcnAxis = (AliRsnValue*)fAxisList.At(i);
 241       if (!fcnAxis) {
 242          AliError("Empty axis: doing unique bin betweeen -100000 and 100000");
 243          continue;
 244       }
 245       fHSparse->SetBinEdges(i, fcnAxis->GetArray().GetArray());
 246    }
 247
 248    delete [] nbins;
 249
 250    return fHSparse;
 251 }
 252
 253
 254 //________________________________________________________________________________________
 255 Bool_t AliRsnFunction::Fill()
 256 {
 257 //
 258 // Fill function histogram with values computed from given input object.
 259 //
 260
 261    AliDebug(AliLog::kDebug + 2, "->");
 262
 263    Int_t  i;
 264    Double_t *values = new Double_t[fSize];
 265    Bool_t    computeOK = kFALSE;
 266
 267    AliRsnValue *fcnAxis = 0;
 268    for (i = 0; i < fSize; i++) {
 269       values[i] = 0.0;
 270       fcnAxis = (AliRsnValue*)fAxisList.At(i);
 271       if (!fcnAxis) continue;
 272       switch (fcnAxis->GetTargetType()) {
 273          case AliRsnTarget::kMother:
 274             fcnAxis->SetSupportObject(fPairDef);
 275             computeOK = fcnAxis->Eval(fPair);
 276             break;
 277          case AliRsnTarget::kEvent:
 278             computeOK = fcnAxis->Eval(fEvent);
 279             break;
 280          default:
 281             AliError(Form("Allowed targets are mothers and events; cannot use axis '%s' which has target '%s'", fcnAxis->GetName(), fcnAxis->GetTargetTypeName()));
 282             computeOK = kFALSE;
 283       }
 284       if (computeOK) values[i] = ((Float_t)fcnAxis->GetComputedValue());
 285    }
 286
 287    // fill histogram
 288    if (fUseTH1) {
 289       // check presence of output histogram
 290       if (!fH1) {
 291          AliError("Required a TH1 which is not initialized");
 292          delete [] values;
 293          return kFALSE;
 294       }
 295
 296       // fill according to dimensions
 297       switch (fSize) {
 298          case 1: {
 299             TH1F *h1 = (TH1F*)fH1;
 300             h1->Fill(values[0]);
 301          }
 302          break;
 303          case 2: {
 304             TH2F *h2 = (TH2F*)fH1;
 305             h2->Fill(values[0], values[1]);
 306          }
 307          break;
 308          case 3: {
 309             TH3F *h3 = (TH3F*)fH1;
 310             h3->Fill(values[0], values[1], values[2]);
 311          }
 312          break;
 313          default:
 314             AliError(Form("Wrong size : %d", fSize));
 315             delete [] values;
 316             return kFALSE;
 317       }
 318    } else {
 319       // check presence of output histogram
 320       if (!fHSparse) {
 321          AliError("Required a THnSparseF which is not initialized");
 322          delete [] values;
 323          return kFALSE;
 324       }
 325
 326       fHSparse->Fill(values);
 327    }
 328
 329    delete [] values;
 330
 331    AliDebug(AliLog::kDebug + 2, "->");
 332    return kTRUE;
 333 }