PWGLF/RESONANCES/AliRsnMiniOutput.cxx

   1 //
   2 // Mini-Output
   3 // All the definitions needed for building a RSN histogram
   4 // including:
   5 // -- properties of resonance (mass, PDG code if needed)
   6 // -- properties of daughters (assigned mass, charges)
   7 // -- definition of output histogram
   8 //
   9
  10 #include "Riostream.h"
  11
  12 #include "TH1.h"
  13 #include "TH2.h"
  14 #include "TH3.h"
  15 #include "TList.h"
  16 #include "TMath.h"
  17 #include "THnSparse.h"
  18 #include "TString.h"
  19 #include "TClonesArray.h"
  20
  21 #include "AliRsnMiniParticle.h"
  22 #include "AliRsnMiniPair.h"
  23 #include "AliRsnMiniEvent.h"
  24 #include "AliAODEvent.h"
  25
  26 #include "AliLog.h"
  27 #include "AliRsnCutSet.h"
  28 #include "AliRsnMiniAxis.h"
  29 #include "AliRsnMiniOutput.h"
  30 #include "AliRsnMiniValue.h"
  31
  32 ClassImp(AliRsnMiniOutput)
  33
  34 //__________________________________________________________________________________________________
  35 AliRsnMiniOutput::AliRsnMiniOutput() :
  36    TNamed(),
  37    fOutputType(kTypes),
  38    fComputation(kComputations),
  39    fMotherPDG(0),
  40    fMotherMass(0.0),
  41    fPairCuts(0x0),
  42    fOutputID(-1),
  43    fAxes("AliRsnMiniAxis", 0),
  44    fComputed(0),
  45    fPair(),
  46    fList(0x0),
  47    fSel1(0),
  48    fSel2(0),
  49    fMaxNSisters(-1),
  50    fCheckP(kFALSE),
  51    fCheckHistRange(kTRUE)
  52 {
  53 //
  54 // Constructor
  55 //
  56
  57    fCutID[0] = fCutID[1] = -1;
  58    fDaughter[0] = fDaughter[1] = AliRsnDaughter::kUnknown;
  59    fCharge[0] = fCharge[1] = 0;
  60 }
  61
  62 //__________________________________________________________________________________________________
  63 AliRsnMiniOutput::AliRsnMiniOutput(const char *name, EOutputType type, EComputation src) :
  64    TNamed(name, ""),
  65    fOutputType(type),
  66    fComputation(src),
  67    fMotherPDG(0),
  68    fMotherMass(0.0),
  69    fPairCuts(0x0),
  70    fOutputID(-1),
  71    fAxes("AliRsnMiniAxis", 0),
  72    fComputed(0),
  73    fPair(),
  74    fList(0x0),
  75    fSel1(0),
  76    fSel2(0),
  77    fMaxNSisters(-1),
  78    fCheckP(kFALSE),
  79    fCheckHistRange(kTRUE)
  80 {
  81 //
  82 // Constructor
  83 //
  84
  85    fCutID[0] = fCutID[1] = -1;
  86    fDaughter[0] = fDaughter[1] = AliRsnDaughter::kUnknown;
  87    fCharge[0] = fCharge[1] = 0;
  88 }
  89
  90 //__________________________________________________________________________________________________
  91 AliRsnMiniOutput::AliRsnMiniOutput(const char *name, const char *outType, const char *compType) :
  92    TNamed(name, ""),
  93    fOutputType(kTypes),
  94    fComputation(kComputations),
  95    fMotherPDG(0),
  96    fMotherMass(0.0),
  97    fPairCuts(0x0),
  98    fOutputID(-1),
  99    fAxes("AliRsnMiniAxis", 0),
 100    fComputed(0),
 101    fPair(),
 102    fList(0x0),
 103    fSel1(0),
 104    fSel2(0),
 105    fMaxNSisters(-1),
 106    fCheckP(kFALSE),
 107    fCheckHistRange(kTRUE)
 108 {
 109 //
 110 // Constructor, with a more user friendly implementation, where
 111 // the user sets the type of output and computations through conventional strings:
 112 //
 113 // Output:
 114 //    -- "HIST"    --> common histogram (up to 3 dimensions)
 115 //    -- "SPARSE"  --> sparse histogram
 116 //
 117 // Computation:
 118 //    -- "EVENT"   --> event-only computations
 119 //    -- "PAIR"    --> track pair computations (default)
 120 //    -- "MIX"     --> event mixing (like track pair, but different events)
 121 //    -- "ROTATE1" --> rotated background (rotate first track)
 122 //    -- "ROTATE2" --> rotated background (rotate second track)
 123 //    -- "TRUE"    --> true pairs (like track pair, but checking that come from same mother)
 124 //    -- "MOTHER"  --> mother (loop on MC directly for mothers --> denominator of efficiency)
 125 //
 126
 127    TString input;
 128
 129    // understand output type
 130    input = outType;
 131    input.ToUpper();
 132    if (!input.CompareTo("HIST"))
 133       fOutputType = kHistogram;
 134    else if (!input.CompareTo("SPARSE"))
 135       fOutputType = kHistogramSparse;
 136    else
 137       AliWarning(Form("String '%s' does not define a meaningful output type", outType));
 138
 139    // understand computation type
 140    input = compType;
 141    input.ToUpper();
 142    if (!input.CompareTo("EVENT"))
 143       fComputation = kEventOnly;
 144    else if (!input.CompareTo("PAIR"))
 145       fComputation = kTrackPair;
 146    else if (!input.CompareTo("MIX"))
 147       fComputation = kTrackPairMix;
 148    else if (!input.CompareTo("ROTATE1"))
 149       fComputation = kTrackPairRotated1;
 150    else if (!input.CompareTo("ROTATE2"))
 151       fComputation = kTrackPairRotated2;
 152    else if (!input.CompareTo("TRUE"))
 153       fComputation = kTruePair;
 154    else if (!input.CompareTo("MOTHER"))
 155       fComputation = kMother;
 156    else
 157       AliWarning(Form("String '%s' does not define a meaningful computation type", compType));
 158
 159    fCutID[0] = fCutID[1] = -1;
 160    fDaughter[0] = fDaughter[1] = AliRsnDaughter::kUnknown;
 161    fCharge[0] = fCharge[1] = 0;
 162 }
 163
 164 //__________________________________________________________________________________________________
 165 AliRsnMiniOutput::AliRsnMiniOutput(const AliRsnMiniOutput &copy) :
 166    TNamed(copy),
 167    fOutputType(copy.fOutputType),
 168    fComputation(copy.fComputation),
 169    fMotherPDG(copy.fMotherPDG),
 170    fMotherMass(copy.fMotherMass),
 171    fPairCuts(copy.fPairCuts),
 172    fOutputID(copy.fOutputID),
 173    fAxes(copy.fAxes),
 174    fComputed(copy.fComputed),
 175    fPair(),
 176    fList(copy.fList),
 177    fSel1(0),
 178    fSel2(0),
 179    fMaxNSisters(-1),
 180    fCheckP(kFALSE),
 181    fCheckHistRange(copy.fCheckHistRange)
 182 {
 183 //
 184 // Copy constructor
 185 //
 186
 187    Int_t i;
 188    for (i = 0; i < 2; i++) {
 189       fCutID[i] = copy.fCutID[i];
 190       fDaughter[i] = copy.fDaughter[i];
 191       fCharge[i] = copy.fCharge[i];
 192    }
 193 }
 194
 195 //__________________________________________________________________________________________________
 196 AliRsnMiniOutput &AliRsnMiniOutput::operator=(const AliRsnMiniOutput &copy)
 197 {
 198 //
 199 // Assignment operator
 200 //
 201    if (this == &copy)
 202       return *this;
 203    fOutputType = copy.fOutputType;
 204    fComputation = copy.fComputation;
 205    fMotherPDG = copy.fMotherPDG;
 206    fMotherMass = copy.fMotherMass;
 207    fPairCuts = copy.fPairCuts;
 208    fOutputID = copy.fOutputID;
 209    fAxes = copy.fAxes;
 210    fComputed = copy.fComputed;
 211    fList = copy.fList;
 212
 213    Int_t i;
 214    for (i = 0; i < 2; i++) {
 215       fCutID[i] = copy.fCutID[i];
 216       fDaughter[i] = copy.fDaughter[i];
 217       fCharge[i] = copy.fCharge[i];
 218    }
 219
 220    fSel1.Set(0);
 221    fSel2.Set(0);
 222    fMaxNSisters = copy.fMaxNSisters;
 223    fCheckP = copy.fCheckP;
 224    fCheckHistRange = copy.fCheckHistRange;
 225
 226    return (*this);
 227 }
 228
 229
 230 //__________________________________________________________________________________________________
 231 void AliRsnMiniOutput::AddAxis(Int_t i, Int_t nbins, Double_t min, Double_t max)
 232 {
 233 //
 234 // Create a new axis reference
 235 //
 236
 237    Int_t size = fAxes.GetEntries();
 238    new (fAxes[size]) AliRsnMiniAxis(i, nbins, min, max);
 239 }
 240
 241 //__________________________________________________________________________________________________
 242 void AliRsnMiniOutput::AddAxis(Int_t i, Double_t min, Double_t max, Double_t step)
 243 {
 244 //
 245 // Create a new axis reference
 246 //
 247
 248    Int_t size = fAxes.GetEntries();
 249    new (fAxes[size]) AliRsnMiniAxis(i, min, max, step);
 250 }
 251
 252 //__________________________________________________________________________________________________
 253 void AliRsnMiniOutput::AddAxis(Int_t i, Int_t nbins, Double_t *values)
 254 {
 255 //
 256 // Create a new axis reference
 257 //
 258
 259    Int_t size = fAxes.GetEntries();
 260    new (fAxes[size]) AliRsnMiniAxis(i, nbins, values);
 261 }
 262
 263 //__________________________________________________________________________________________________
 264 Bool_t AliRsnMiniOutput::Init(const char *prefix, TList *list)
 265 {
 266 //
 267 // Initialize properly the histogram and add it to the argument list
 268 //
 269
 270    if (!list) {
 271       AliError("Required an output list");
 272       return kFALSE;
 273    }
 274
 275    fList = list;
 276    Int_t size = fAxes.GetEntries();
 277    if (size < 1) {
 278       AliWarning(Form("[%s] Cannot initialize histogram with less than 1 axis", GetName()));
 279       return kFALSE;
 280    }
 281
 282    switch (fOutputType) {
 283       case kHistogram:
 284          if (size <= 3) {
 285             CreateHistogram(Form("%s_%s", prefix, GetName()));
 286          } else {
 287             AliInfo(Form("[%s] Added %d > 3 axes. Creating a sparse histogram", GetName(), size));
 288             fOutputType = kHistogramSparse;
 289             CreateHistogramSparse(Form("%s_%s", prefix, GetName()));
 290          }
 291          return kTRUE;
 292       case kHistogramSparse:
 293          CreateHistogramSparse(Form("%s_%s", prefix, GetName()));
 294          return kTRUE;
 295       default:
 296          AliError("Wrong output histogram definition");
 297          return kFALSE;
 298    }
 299 }
 300
 301 //__________________________________________________________________________________________________
 302 void AliRsnMiniOutput::CreateHistogram(const char *name)
 303 {
 304 //
 305 // Initialize the 'default' TH1 output object.
 306 // In case one of the expected axes is NULL, the initialization fails.
 307 //
 308
 309    Int_t size = fAxes.GetEntries();
 310    AliInfo(Form("Histogram name = '%s', with %d axes", name, size));
 311
 312    // we expect to have maximum 3 axes in this case
 313    AliRsnMiniAxis *xAxis = 0x0, *yAxis = 0x0, *zAxis = 0x0;
 314    if (size >= 1) xAxis = (AliRsnMiniAxis *)fAxes[0];
 315    if (size >= 2) yAxis = (AliRsnMiniAxis *)fAxes[1];
 316    if (size >= 3) zAxis = (AliRsnMiniAxis *)fAxes[2];
 317
 318    // create histogram depending on the number of axes
 319    TH1 *h1 = 0x0;
 320    if (xAxis && yAxis && zAxis) {
 321       h1 = new TH3F(name, "", xAxis->NBins(), xAxis->BinArray(), yAxis->NBins(), yAxis->BinArray(), zAxis->NBins(), zAxis->BinArray());
 322    } else if (xAxis && yAxis) {
 323       h1 = new TH2F(name, "", xAxis->NBins(), xAxis->BinArray(), yAxis->NBins(), yAxis->BinArray());
 324    } else if (xAxis) {
 325       h1 = new TH1F(name, "", xAxis->NBins(), xAxis->BinArray());
 326    } else {
 327       AliError("No axis was initialized");
 328       return;
 329    }
 330
 331    // switch the correct computation of errors
 332    if (h1 && fList) {
 333       h1->Sumw2();
 334       fList->Add(h1);
 335       fOutputID = fList->IndexOf(h1);
 336    }
 337 }
 338
 339 //________________________________________________________________________________________
 340 void AliRsnMiniOutput::CreateHistogramSparse(const char *name)
 341 {
 342 //
 343 // Initialize the THnSparse output object.
 344 // In case one of the expected axes is NULL, the initialization fails.
 345 //
 346
 347    Int_t size = fAxes.GetEntries();
 348    AliInfo(Form("Sparse histogram name = '%s', with %d axes", name, size));
 349
 350    // retrieve binnings and sizes of all axes
 351    // since the check for null values is done in Init(),
 352    // we assume that here they must all be well defined
 353    Int_t i, *nbins = new Int_t[size];
 354    for (i = 0; i < size; i++) {
 355       AliRsnMiniAxis *axis = (AliRsnMiniAxis *)fAxes[i];
 356       nbins[i] = axis->NBins();
 357    }
 358
 359    // create fHSparseogram
 360    THnSparseF *h1 = new THnSparseF(name, "", size, nbins);
 361
 362    // update the various axes using the definitions given in the array of axes here
 363    for (i = 0; i < size; i++) {
 364       AliRsnMiniAxis *axis = (AliRsnMiniAxis *)fAxes[i];
 365       h1->GetAxis(i)->Set(nbins[i], axis->BinArray());
 366    }
 367
 368    // clear heap
 369    delete [] nbins;
 370
 371    // add to list
 372    if (h1 && fList) {
 373       h1->Sumw2();
 374       fList->Add(h1);
 375       fOutputID = fList->IndexOf(h1);
 376    }
 377 }
 378
 379 //________________________________________________________________________________________
 380 Bool_t AliRsnMiniOutput::FillEvent(AliRsnMiniEvent *event, TClonesArray *valueList)
 381 {
 382 //
 383 // Compute values for event-based computations (does not use the pair)
 384 //
 385
 386    // check computation type
 387    if (fComputation != kEventOnly) {
 388       AliError("This method can be called only for event-based computations");
 389       return kFALSE;
 390    }
 391
 392    // compute & fill
 393    ComputeValues(event, valueList);
 394    FillHistogram();
 395    return kTRUE;
 396 }
 397
 398 //________________________________________________________________________________________
 399 Bool_t AliRsnMiniOutput::FillMother(const AliRsnMiniPair *pair, AliRsnMiniEvent *event, TClonesArray *valueList)
 400 {
 401 //
 402 // Compute values for mother-based computations
 403 //
 404
 405    // check computation type
 406    if (fComputation != kMother) {
 407       AliError("This method can be called only for mother-based computations");
 408       return kFALSE;
 409    }
 410
 411    // copy passed pair info
 412    fPair = (*pair);
 413
 414    // check pair against cuts
 415    if (fPairCuts) if (!fPairCuts->IsSelected(&fPair)) return kFALSE;
 416
 417    // compute & fill
 418    ComputeValues(event, valueList);
 419    FillHistogram();
 420    return kTRUE;
 421 }
 422
 423 //________________________________________________________________________________________
 424 Int_t AliRsnMiniOutput::FillPair(AliRsnMiniEvent *event1, AliRsnMiniEvent *event2, TClonesArray *valueList, Bool_t refFirst)
 425 {
 426 //
 427 // Loops on the passed mini-event, and for each pair of particles
 428 // which satisfy the charge and cut requirements defined here, add an entry.
 429 // Returns the number of successful fillings.
 430 // Last argument tells if the reference event for event-based values is the first or the second.
 431 //
 432
 433    // check computation type
 434    Bool_t okComp = kFALSE;
 435    if (fComputation == kTrackPair)         okComp = kTRUE;
 436    if (fComputation == kTrackPairMix)      okComp = kTRUE;
 437    if (fComputation == kTrackPairRotated1) okComp = kTRUE;
 438    if (fComputation == kTrackPairRotated2) okComp = kTRUE;
 439    if (fComputation == kTruePair)          okComp = kTRUE;
 440    if (!okComp) {
 441       AliError(Form("[%s] This method can be called only for pair-based computations", GetName()));
 442       return kFALSE;
 443    }
 444
 445    // loop variables
 446    Int_t i1, i2, start, nadded = 0;
 447    AliRsnMiniParticle *p1, *p2;
 448
 449    // it is necessary to know if criteria for the two daughters are the same
 450    // and if the two events are the same or not (mixing)
 451    //Bool_t sameCriteria = ((fCharge[0] == fCharge[1]) && (fCutID[0] == fCutID[1]));
 452    Bool_t sameCriteria = ((fCharge[0] == fCharge[1]) && (fDaughter[0] == fDaughter[1]));
 453    Bool_t sameEvent = (event1->ID() == event2->ID());
 454
 455    TString selList1  = "";
 456    TString selList2  = "";
 457    Int_t   n1 = event1->CountParticles(fSel1, fCharge[0], fCutID[0]);
 458    Int_t   n2 = event2->CountParticles(fSel2, fCharge[1], fCutID[1]);
 459    for (i1 = 0; i1 < n1; i1++) selList1.Append(Form("%d ", fSel1[i1]));
 460    for (i2 = 0; i2 < n2; i2++) selList2.Append(Form("%d ", fSel2[i2]));
 461    AliDebugClass(1, Form("[%10s] Part #1: [%s] -- evID %6d -- charge = %c -- cut ID = %d --> %4d tracks (%s)", GetName(), (event1 == event2 ? "def" : "mix"), event1->ID(), fCharge[0], fCutID[0], n1, selList1.Data()));
 462    AliDebugClass(1, Form("[%10s] Part #2: [%s] -- evID %6d -- charge = %c -- cut ID = %d --> %4d tracks (%s)", GetName(), (event1 == event2 ? "def" : "mix"), event2->ID(), fCharge[1], fCutID[1], n2, selList2.Data()));
 463    if (!n1 || !n2) {
 464       AliDebugClass(1, "No pairs to mix");
 465       return 0;
 466    }
 467
 468    // external loop
 469    for (i1 = 0; i1 < n1; i1++) {
 470       p1 = event1->GetParticle(fSel1[i1]);
 471       //p1 = event1->GetParticle(i1);
 472       //if (p1->Charge() != fCharge[0]) continue;
 473       //if (!p1->HasCutBit(fCutID[0])) continue;
 474       // define starting point for inner loop
 475       // if daughter selection criteria (charge, cuts) are the same
 476       // and the two events coincide, internal loop must start from
 477       // the first track *after* current one;
 478       // otherwise it starts from the beginning
 479       start = ((sameEvent && sameCriteria) ? i1 + 1 : 0);
 480       AliDebugClass(2, Form("Start point = %d", start));
 481       // internal loop
 482       for (i2 = start; i2 < n2; i2++) {
 483          p2 = event2->GetParticle(fSel2[i2]);
 484          //p2 = event2->GetParticle(i2);
 485          //if (p2->Charge() != fCharge[1]) continue;
 486          //if (!p2->HasCutBit(fCutID[1])) continue;
 487          // avoid to mix a particle with itself
 488          if (sameEvent && (p1->Index() == p2->Index())) {
 489             AliDebugClass(2, "Skipping same index");
 490             continue;
 491          }
 492          // sum momenta
 493          fPair.Fill(p1, p2, GetMass(0), GetMass(1), fMotherMass);
 494          // do rotation if needed
 495          if (fComputation == kTrackPairRotated1) fPair.InvertP(kTRUE);
 496          if (fComputation == kTrackPairRotated2) fPair.InvertP(kFALSE);
 497          // if required, check that this is a true pair
 498          if (fComputation == kTruePair) {
 499             if (fPair.Mother() < 0)  {
 500                continue;
 501             } else if (fPair.MotherPDG() != fMotherPDG) {
 502                continue;
 503             }
 504             Bool_t decayMatch = kFALSE;
 505             if (p1->PDGAbs() == AliRsnDaughter::SpeciesPDG(fDaughter[0]) && p2->PDGAbs() == AliRsnDaughter::SpeciesPDG(fDaughter[1]))
 506                decayMatch = kTRUE;
 507             if (p2->PDGAbs() == AliRsnDaughter::SpeciesPDG(fDaughter[0]) && p1->PDGAbs() == AliRsnDaughter::SpeciesPDG(fDaughter[1]))
 508                decayMatch = kTRUE;
 509             if (!decayMatch) continue;
 510             if ( (fMaxNSisters>0) && (p1->NTotSisters()==p2->NTotSisters()) && (p1->NTotSisters()>fMaxNSisters)) continue;
 511             if ( fCheckP &&(TMath::Abs(fPair.PmotherX()-(p1->Px(1)+p2->Px(1)))/(TMath::Abs(fPair.PmotherX())+1.e-13)) > 0.00001 &&
 512                           (TMath::Abs(fPair.PmotherY()-(p1->Py(1)+p2->Py(1)))/(TMath::Abs(fPair.PmotherY())+1.e-13)) > 0.00001 &&
 513                           (TMath::Abs(fPair.PmotherZ()-(p1->Pz(1)+p2->Pz(1)))/(TMath::Abs(fPair.PmotherZ())+1.e-13)) > 0.00001 ) continue;
 514          }
 515          // check pair against cuts
 516          if (fPairCuts) {
 517             if (!fPairCuts->IsSelected(&fPair)) continue;
 518          }
 519          // get computed values & fill histogram
 520          nadded++;
 521          if (refFirst) ComputeValues(event1, valueList); else ComputeValues(event2, valueList);
 522          FillHistogram();
 523       } // end internal loop
 524    } // end external loop
 525
 526    AliDebugClass(1, Form("Pairs added in total = %4d", nadded));
 527    return nadded;
 528 }
 529
 530 //________________________________________________________________________________________
 531 void AliRsnMiniOutput::ComputeValues(AliRsnMiniEvent *event, TClonesArray *valueList)
 532 {
 533 //
 534 // Using the arguments and the internal 'fPair' data member,
 535 // compute all values to be stored in the histogram
 536 //
 537
 538    // check size of computed array
 539    Int_t size = fAxes.GetEntries();
 540    if (fComputed.GetSize() != size) fComputed.Set(size);
 541
 542    Int_t i, ival, nval = valueList->GetEntries();
 543
 544    for (i = 0; i < size; i++) {
 545       fComputed[i] = 1E20;
 546       AliRsnMiniAxis *axis = (AliRsnMiniAxis *)fAxes[i];
 547       if (!axis) {
 548          AliError("Null axis");
 549          continue;
 550       }
 551       ival = axis->GetValueID();
 552       if (ival < 0 || ival >= nval) {
 553          AliError(Form("Required value #%d, while maximum is %d", ival, nval));
 554          continue;
 555       }
 556       AliRsnMiniValue *val = (AliRsnMiniValue *)valueList->At(ival);
 557       if (!val) {
 558          AliError(Form("Value in position #%d is NULL", ival));
 559          continue;
 560       }
 561       // if none of the above exit points is taken, compute value
 562       fComputed[i] = val->Eval(&fPair, event);
 563    }
 564 }
 565
 566 //________________________________________________________________________________________
 567 void AliRsnMiniOutput::FillHistogram()
 568 {
 569 //
 570 // Fills the internal histogram using the current values stored in the
 571 // 'fComputed' array, in the order as they are stored, up to the max
 572 // dimension of the initialized histogram itself.
 573 //
 574
 575    // retrieve object from list
 576    if (!fList) {
 577       AliError("List pointer is NULL");
 578       return;
 579    }
 580    TObject *obj = fList->At(fOutputID);
 581
 582    if (obj->InheritsFrom(TH1F::Class())) {
 583       ((TH1F *)obj)->Fill(fComputed[0]);
 584    } else if (obj->InheritsFrom(TH2F::Class())) {
 585       ((TH2F *)obj)->Fill(fComputed[0], fComputed[1]);
 586    } else if (obj->InheritsFrom(TH3F::Class())) {
 587       ((TH3F *)obj)->Fill(fComputed[0], fComputed[1], fComputed[2]);
 588    } else if (obj->InheritsFrom(THnSparseF::Class())) {
 589       THnSparseF *h = (THnSparseF *)obj;
 590       if (fCheckHistRange) {
 591          for (Int_t iAxis = 0; iAxis<h->GetNdimensions(); iAxis++) {
 592             if (fComputed.At(iAxis)>h->GetAxis(iAxis)->GetXmax() || fComputed.At(iAxis)<h->GetAxis(iAxis)->GetXmin()) return;
 593          }
 594       }
 595       h->Fill(fComputed.GetArray());
 596    } else {
 597       AliError("No output initialized");
 598    }
 599 }