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