Added new cuts for (TPC pid & TOF match) and (TPC pid & TOF veto) for K* pA analysis
[u/mrichter/AliRoot.git] / PWGLF / RESONANCES / AliRsnMiniOutput.cxx
CommitLineData
03d23846 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
61f275d1 8//
03d23846 9
10#include "Riostream.h"
11
12#include "TH1.h"
13#include "TH2.h"
14#include "TH3.h"
15#include "TList.h"
45aa62b9 16#include "TMath.h"
03d23846 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
31ClassImp(AliRsnMiniOutput)
32
33//__________________________________________________________________________________________________
34AliRsnMiniOutput::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(),
a2455d2a 45 fList(0x0),
46 fSel1(0),
74d60285 47 fSel2(0),
48 fCheckHistRange(kTRUE)
03d23846 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//__________________________________________________________________________________________________
60AliRsnMiniOutput::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(),
a2455d2a 71 fList(0x0),
72 fSel1(0),
74d60285 73 fSel2(0),
74 fCheckHistRange(kTRUE)
03d23846 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//__________________________________________________________________________________________________
86AliRsnMiniOutput::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(),
a2455d2a 97 fList(0x0),
98 fSel1(0),
74d60285 99 fSel2(0),
100 fCheckHistRange(kTRUE)
03d23846 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:
6aaeb33c 107// -- "HIST" --> common histogram (up to 3 dimensions)
108// -- "SPARSE" --> sparse histogram
61f275d1 109//
110// Computation:
6aaeb33c 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)
03d23846 118//
119
120 TString input;
61f275d1 121
03d23846 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));
61f275d1 131
03d23846 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;
6aaeb33c 141 else if (!input.CompareTo("ROTATE1"))
142 fComputation = kTrackPairRotated1;
143 else if (!input.CompareTo("ROTATE2"))
144 fComputation = kTrackPairRotated2;
03d23846 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));
61f275d1 151
03d23846 152 fCutID[0] = fCutID[1] = -1;
153 fDaughter[0] = fDaughter[1] = AliRsnDaughter::kUnknown;
154 fCharge[0] = fCharge[1] = 0;
155}
156
157//__________________________________________________________________________________________________
158AliRsnMiniOutput::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(),
a2455d2a 169 fList(copy.fList),
170 fSel1(0),
74d60285 171 fSel2(0),
172 fCheckHistRange(copy.fCheckHistRange)
03d23846 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//__________________________________________________________________________________________________
61f275d1 187AliRsnMiniOutput &AliRsnMiniOutput::operator=(const AliRsnMiniOutput &copy)
03d23846 188{
189//
190// Assignment operator
191//
61f275d1 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);
74d60285 213 fCheckHistRange = copy.fCheckHistRange;
61f275d1 214
03d23846 215 return (*this);
216}
217
218
219//__________________________________________________________________________________________________
220void AliRsnMiniOutput::AddAxis(Int_t i, Int_t nbins, Double_t min, Double_t max)
221{
222//
223// Create a new axis reference
224//
225
61f275d1 226 Int_t size = fAxes.GetEntries();
03d23846 227 new (fAxes[size]) AliRsnMiniAxis(i, nbins, min, max);
228}
229
230//__________________________________________________________________________________________________
231void AliRsnMiniOutput::AddAxis(Int_t i, Double_t min, Double_t max, Double_t step)
232{
233//
234// Create a new axis reference
235//
236
61f275d1 237 Int_t size = fAxes.GetEntries();
03d23846 238 new (fAxes[size]) AliRsnMiniAxis(i, min, max, step);
239}
240
241//__________________________________________________________________________________________________
242void AliRsnMiniOutput::AddAxis(Int_t i, Int_t nbins, Double_t *values)
243{
244//
245// Create a new axis reference
246//
247
61f275d1 248 Int_t size = fAxes.GetEntries();
03d23846 249 new (fAxes[size]) AliRsnMiniAxis(i, nbins, values);
250}
251
252//__________________________________________________________________________________________________
253Bool_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 }
61f275d1 263
03d23846 264 fList = list;
45aa62b9 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 }
03d23846 270
271 switch (fOutputType) {
272 case kHistogram:
d573d2fb 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 }
03d23846 280 return kTRUE;
281 case kHistogramSparse:
6aaeb33c 282 CreateHistogramSparse(Form("%s_%s", prefix, GetName()));
03d23846 283 return kTRUE;
284 default:
285 AliError("Wrong output histogram definition");
286 return kFALSE;
287 }
288}
289
290//__________________________________________________________________________________________________
291void 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;
61f275d1 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
03d23846 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 }
61f275d1 319
03d23846 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//________________________________________________________________________________________
329void 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
d573d2fb 336 Int_t size = fAxes.GetEntries();
337 AliInfo(Form("Sparse histogram name = '%s', with %d axes", name, size));
61f275d1 338
03d23846 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
03d23846 342 Int_t i, *nbins = new Int_t[size];
343 for (i = 0; i < size; i++) {
61f275d1 344 AliRsnMiniAxis *axis = (AliRsnMiniAxis *)fAxes[i];
03d23846 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++) {
61f275d1 353 AliRsnMiniAxis *axis = (AliRsnMiniAxis *)fAxes[i];
03d23846 354 h1->GetAxis(i)->Set(nbins[i], axis->BinArray());
355 }
356
357 // clear heap
358 delete [] nbins;
61f275d1 359
03d23846 360 // add to list
361 if (h1 && fList) {
362 h1->Sumw2();
363 fList->Add(h1);
364 fOutputID = fList->IndexOf(h1);
365 }
366}
367
368//________________________________________________________________________________________
45aa62b9 369Bool_t AliRsnMiniOutput::FillEvent(AliRsnMiniEvent *event, TClonesArray *valueList)
03d23846 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
45aa62b9 381 // compute & fill
382 ComputeValues(event, valueList);
383 FillHistogram();
384 return kTRUE;
03d23846 385}
386
387//________________________________________________________________________________________
45aa62b9 388Bool_t AliRsnMiniOutput::FillMother(const AliRsnMiniPair *pair, AliRsnMiniEvent *event, TClonesArray *valueList)
03d23846 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 }
61f275d1 399
03d23846 400 // copy passed pair info
401 fPair = (*pair);
61f275d1 402
03d23846 403 // check pair against cuts
404 if (fPairCuts) if (!fPairCuts->IsSelected(&fPair)) return kFALSE;
405
45aa62b9 406 // compute & fill
407 ComputeValues(event, valueList);
408 FillHistogram();
409 return kTRUE;
03d23846 410}
411
412//________________________________________________________________________________________
d573d2fb 413Int_t AliRsnMiniOutput::FillPair(AliRsnMiniEvent *event1, AliRsnMiniEvent *event2, TClonesArray *valueList, Bool_t refFirst)
03d23846 414{
415//
45aa62b9 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.
d573d2fb 419// Last argument tells if the reference event for event-based values is the first or the second.
03d23846 420//
421
422 // check computation type
6aaeb33c 423 Bool_t okComp = kFALSE;
45aa62b9 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;
6aaeb33c 429 if (!okComp) {
45aa62b9 430 AliError(Form("[%s] This method can be called only for pair-based computations", GetName()));
03d23846 431 return kFALSE;
432 }
61f275d1 433
45aa62b9 434 // loop variables
435 Int_t i1, i2, start, nadded = 0;
45aa62b9 436 AliRsnMiniParticle *p1, *p2;
61f275d1 437
45aa62b9 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)
7196ee4f 440 //Bool_t sameCriteria = ((fCharge[0] == fCharge[1]) && (fCutID[0] == fCutID[1]));
441 Bool_t sameCriteria = ((fCharge[0] == fCharge[1]) && (fDaughter[0] == fDaughter[1]));
45aa62b9 442 Bool_t sameEvent = (event1->ID() == event2->ID());
61f275d1 443
9e3a9020 444 TString selList1 = "";
445 TString selList2 = "";
a2455d2a 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()));
9e3a9020 452 if (!n1 || !n2) {
d573d2fb 453 AliDebugClass(1, "No pairs to mix");
454 return 0;
455 }
61f275d1 456
45aa62b9 457 // external loop
458 for (i1 = 0; i1 < n1; i1++) {
a2455d2a 459 p1 = event1->GetParticle(fSel1[i1]);
9e3a9020 460 //p1 = event1->GetParticle(i1);
461 //if (p1->Charge() != fCharge[0]) continue;
462 //if (!p1->HasCutBit(fCutID[0])) continue;
45aa62b9 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);
7196ee4f 469 AliDebugClass(2, Form("Start point = %d", start));
45aa62b9 470 // internal loop
471 for (i2 = start; i2 < n2; i2++) {
a2455d2a 472 p2 = event2->GetParticle(fSel2[i2]);
9e3a9020 473 //p2 = event2->GetParticle(i2);
474 //if (p2->Charge() != fCharge[1]) continue;
475 //if (!p2->HasCutBit(fCutID[1])) continue;
45aa62b9 476 // avoid to mix a particle with itself
477 if (sameEvent && (p1->Index() == p2->Index())) {
d573d2fb 478 AliDebugClass(2, "Skipping same index");
45aa62b9 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;
55de729a 490 } else if (fPair.MotherPDG() != fMotherPDG) {
45aa62b9 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) {
9e7b94f5 502 if (!fPairCuts->IsSelected(&fPair)) continue;
45aa62b9 503 }
504 // get computed values & fill histogram
9e3a9020 505 nadded++;
61f275d1 506 if (refFirst) ComputeValues(event1, valueList); else ComputeValues(event2, valueList);
45aa62b9 507 FillHistogram();
45aa62b9 508 } // end internal loop
509 } // end external loop
61f275d1 510
d573d2fb 511 AliDebugClass(1, Form("Pairs added in total = %4d", nadded));
45aa62b9 512 return nadded;
03d23846 513}
514
515//________________________________________________________________________________________
45aa62b9 516void AliRsnMiniOutput::ComputeValues(AliRsnMiniEvent *event, TClonesArray *valueList)
03d23846 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();
61f275d1 528
03d23846 529 for (i = 0; i < size; i++) {
530 fComputed[i] = 1E20;
61f275d1 531 AliRsnMiniAxis *axis = (AliRsnMiniAxis *)fAxes[i];
03d23846 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 }
61f275d1 541 AliRsnMiniValue *val = (AliRsnMiniValue *)valueList->At(ival);
03d23846 542 if (!val) {
543 AliError(Form("Value in position #%d is NULL", ival));
544 continue;
61f275d1 545 }
03d23846 546 // if none of the above exit points is taken, compute value
547 fComputed[i] = val->Eval(&fPair, event);
548 }
03d23846 549}
550
551//________________________________________________________________________________________
45aa62b9 552void AliRsnMiniOutput::FillHistogram()
03d23846 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");
45aa62b9 563 return;
03d23846 564 }
565 TObject *obj = fList->At(fOutputID);
566
567 if (obj->InheritsFrom(TH1F::Class())) {
61f275d1 568 ((TH1F *)obj)->Fill(fComputed[0]);
03d23846 569 } else if (obj->InheritsFrom(TH2F::Class())) {
61f275d1 570 ((TH2F *)obj)->Fill(fComputed[0], fComputed[1]);
03d23846 571 } else if (obj->InheritsFrom(TH3F::Class())) {
61f275d1 572 ((TH3F *)obj)->Fill(fComputed[0], fComputed[1], fComputed[2]);
03d23846 573 } else if (obj->InheritsFrom(THnSparseF::Class())) {
74d60285 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());
03d23846 581 } else {
582 AliError("No output initialized");
03d23846 583 }
03d23846 584}