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