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