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