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