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