Coverity fix
[u/mrichter/AliRoot.git] / PWG2 / 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
8//
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(),
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//__________________________________________________________________________________________________
57AliRsnMiniOutput::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//__________________________________________________________________________________________________
80AliRsnMiniOutput::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:
6aaeb33c 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)
03d23846 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;
6aaeb33c 132 else if (!input.CompareTo("ROTATE1"))
133 fComputation = kTrackPairRotated1;
134 else if (!input.CompareTo("ROTATE2"))
135 fComputation = kTrackPairRotated2;
03d23846 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//__________________________________________________________________________________________________
149AliRsnMiniOutput::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//__________________________________________________________________________________________________
175AliRsnMiniOutput& 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//__________________________________________________________________________________________________
203void 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//__________________________________________________________________________________________________
214void 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//__________________________________________________________________________________________________
225void 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//__________________________________________________________________________________________________
236Bool_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;
45aa62b9 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 }
03d23846 253
254 switch (fOutputType) {
255 case kHistogram:
45aa62b9 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 }
03d23846 263 return kTRUE;
264 case kHistogramSparse:
6aaeb33c 265 CreateHistogramSparse(Form("%s_%s", prefix, GetName()));
03d23846 266 return kTRUE;
267 default:
268 AliError("Wrong output histogram definition");
269 return kFALSE;
270 }
271}
272
273//__________________________________________________________________________________________________
274void 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//________________________________________________________________________________________
312void 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 // retrieve binnings and sizes of all axes
320 // since the check for null values is done in Init(),
321 // we assume that here they must all be well defined
322 Int_t size = fAxes.GetEntries();
323 Int_t i, *nbins = new Int_t[size];
324 for (i = 0; i < size; i++) {
325 AliRsnMiniAxis *axis = (AliRsnMiniAxis*)fAxes[i];
326 nbins[i] = axis->NBins();
327 }
328
329 // create fHSparseogram
330 THnSparseF *h1 = new THnSparseF(name, "", size, nbins);
331
332 // update the various axes using the definitions given in the array of axes here
333 for (i = 0; i < size; i++) {
334 AliRsnMiniAxis *axis = (AliRsnMiniAxis*)fAxes[i];
335 h1->GetAxis(i)->Set(nbins[i], axis->BinArray());
336 }
337
338 // clear heap
339 delete [] nbins;
340
341 // add to list
342 if (h1 && fList) {
343 h1->Sumw2();
344 fList->Add(h1);
345 fOutputID = fList->IndexOf(h1);
346 }
347}
348
349//________________________________________________________________________________________
45aa62b9 350Bool_t AliRsnMiniOutput::FillEvent(AliRsnMiniEvent *event, TClonesArray *valueList)
03d23846 351{
352//
353// Compute values for event-based computations (does not use the pair)
354//
355
356 // check computation type
357 if (fComputation != kEventOnly) {
358 AliError("This method can be called only for event-based computations");
359 return kFALSE;
360 }
361
45aa62b9 362 // compute & fill
363 ComputeValues(event, valueList);
364 FillHistogram();
365 return kTRUE;
03d23846 366}
367
368//________________________________________________________________________________________
45aa62b9 369Bool_t AliRsnMiniOutput::FillMother(const AliRsnMiniPair *pair, AliRsnMiniEvent *event, TClonesArray *valueList)
03d23846 370{
371//
372// Compute values for mother-based computations
373//
374
375 // check computation type
376 if (fComputation != kMother) {
377 AliError("This method can be called only for mother-based computations");
378 return kFALSE;
379 }
380
381 // copy passed pair info
382 fPair = (*pair);
383
384 // check pair against cuts
385 if (fPairCuts) if (!fPairCuts->IsSelected(&fPair)) return kFALSE;
386
45aa62b9 387 // compute & fill
388 ComputeValues(event, valueList);
389 FillHistogram();
390 return kTRUE;
03d23846 391}
392
393//________________________________________________________________________________________
45aa62b9 394Int_t AliRsnMiniOutput::FillPair(AliRsnMiniEvent *event1, AliRsnMiniEvent *event2, TClonesArray *valueList)
03d23846 395{
396//
45aa62b9 397// Loops on the passed mini-event, and for each pair of particles
398// which satisfy the charge and cut requirements defined here, add an entry.
399// Returns the number of successful fillings.
03d23846 400//
401
402 // check computation type
6aaeb33c 403 Bool_t okComp = kFALSE;
45aa62b9 404 if (fComputation == kTrackPair) okComp = kTRUE;
405 if (fComputation == kTrackPairMix) okComp = kTRUE;
406 if (fComputation == kTrackPairRotated1) okComp = kTRUE;
407 if (fComputation == kTrackPairRotated2) okComp = kTRUE;
408 if (fComputation == kTruePair) okComp = kTRUE;
6aaeb33c 409 if (!okComp) {
45aa62b9 410 AliError(Form("[%s] This method can be called only for pair-based computations", GetName()));
03d23846 411 return kFALSE;
412 }
413
45aa62b9 414 // loop variables
415 Int_t i1, i2, start, nadded = 0;
416 Int_t n1 = event1->Particles().GetEntriesFast();
417 Int_t n2 = event2->Particles().GetEntriesFast();
418 AliRsnMiniParticle *p1, *p2;
03d23846 419
45aa62b9 420 // it is necessary to know if criteria for the two daughters are the same
421 // and if the two events are the same or not (mixing)
422 Bool_t sameCriteria = ((fCharge[0] == fCharge[1]) && (fCutID[0] == fCutID[1]));
423 Bool_t sameEvent = (event1->ID() == event2->ID());
6aaeb33c 424
45aa62b9 425 //Int_t nsel1 = event1->CountParticles(fCharge[0], fCutID[0]);
426 //Int_t nsel2 = event2->CountParticles(fCharge[1], fCutID[1]);
427 //cout << "Charge #1 = " << fCharge[0] << " cut bit #1 = " << fCutID[0] << ", tracks = " << nsel1 << endl;
428 //cout << "Charge #2 = " << fCharge[1] << " cut bit #2 = " << fCutID[1] << ", tracks = " << nsel2 << endl;
429 //if (!nsel1 || !nsel2) {
430 // cout << "Nothing to mix" << endl;
431 // return 0;
432 //}
03d23846 433
45aa62b9 434 // external loop
435 for (i1 = 0; i1 < n1; i1++) {
436 p1 = event1->GetParticle(i1);
437 if (p1->Charge() != fCharge[0]) continue;
438 if (!p1->HasCutBit(fCutID[0])) continue;
439 // define starting point for inner loop
440 // if daughter selection criteria (charge, cuts) are the same
441 // and the two events coincide, internal loop must start from
442 // the first track *after* current one;
443 // otherwise it starts from the beginning
444 start = ((sameEvent && sameCriteria) ? i1 + 1 : 0);
445 // internal loop
446 for (i2 = start; i2 < n2; i2++) {
447 p2 = event2->GetParticle(i2);
448 if (p2->Charge() != fCharge[1]) continue;
449 if (!p2->HasCutBit(fCutID[1])) continue;
450 //cout << "Matching: " << i1 << ' ' << i2 << endl;
451 // avoid to mix a particle with itself
452 if (sameEvent && (p1->Index() == p2->Index())) {
453 //cout << "-- skipping same index" << endl;
454 continue;
455 }
456 // sum momenta
457 fPair.Fill(p1, p2, GetMass(0), GetMass(1), fMotherMass);
458 // do rotation if needed
459 if (fComputation == kTrackPairRotated1) fPair.InvertP(kTRUE);
460 if (fComputation == kTrackPairRotated2) fPair.InvertP(kFALSE);
461 // if required, check that this is a true pair
462 if (fComputation == kTruePair) {
463 if (fPair.Mother() < 0) {
464 continue;
465 } else if (TMath::Abs(fPair.MotherPDG()) != fMotherPDG) {
466 continue;
467 }
468 Bool_t decayMatch = kFALSE;
469 if (p1->PDGAbs() == AliRsnDaughter::SpeciesPDG(fDaughter[0]) && p2->PDGAbs() == AliRsnDaughter::SpeciesPDG(fDaughter[1]))
470 decayMatch = kTRUE;
471 if (p2->PDGAbs() == AliRsnDaughter::SpeciesPDG(fDaughter[0]) && p1->PDGAbs() == AliRsnDaughter::SpeciesPDG(fDaughter[1]))
472 decayMatch = kTRUE;
473 if (!decayMatch) continue;
474 }
475 // check pair against cuts
476 if (fPairCuts) {
477 if (!fPairCuts->IsSelected(&fPair)) {
478 //cout << "-- cuts not passed: Y = " << TMath::Abs(fPair.Y(0)) << endl;
479 continue;
480 }
481 }
482 // get computed values & fill histogram
483 ComputeValues(event1, valueList);
484 FillHistogram();
485 nadded++;
486 } // end internal loop
487 } // end external loop
03d23846 488
45aa62b9 489 return nadded;
03d23846 490}
491
492//________________________________________________________________________________________
45aa62b9 493void AliRsnMiniOutput::ComputeValues(AliRsnMiniEvent *event, TClonesArray *valueList)
03d23846 494{
495//
496// Using the arguments and the internal 'fPair' data member,
497// compute all values to be stored in the histogram
498//
499
500 // check size of computed array
501 Int_t size = fAxes.GetEntries();
502 if (fComputed.GetSize() != size) fComputed.Set(size);
503
504 Int_t i, ival, nval = valueList->GetEntries();
505
506 for (i = 0; i < size; i++) {
507 fComputed[i] = 1E20;
508 AliRsnMiniAxis *axis = (AliRsnMiniAxis*)fAxes[i];
509 if (!axis) {
510 AliError("Null axis");
511 continue;
512 }
513 ival = axis->GetValueID();
514 if (ival < 0 || ival >= nval) {
515 AliError(Form("Required value #%d, while maximum is %d", ival, nval));
516 continue;
517 }
518 AliRsnMiniValue *val = (AliRsnMiniValue*)valueList->At(ival);
519 if (!val) {
520 AliError(Form("Value in position #%d is NULL", ival));
521 continue;
522 }
523 // if none of the above exit points is taken, compute value
524 fComputed[i] = val->Eval(&fPair, event);
525 }
03d23846 526}
527
528//________________________________________________________________________________________
45aa62b9 529void AliRsnMiniOutput::FillHistogram()
03d23846 530{
531//
532// Fills the internal histogram using the current values stored in the
533// 'fComputed' array, in the order as they are stored, up to the max
534// dimension of the initialized histogram itself.
535//
536
537 // retrieve object from list
538 if (!fList) {
539 AliError("List pointer is NULL");
45aa62b9 540 return;
03d23846 541 }
542 TObject *obj = fList->At(fOutputID);
543
544 if (obj->InheritsFrom(TH1F::Class())) {
545 ((TH1F*)obj)->Fill(fComputed[0]);
546 } else if (obj->InheritsFrom(TH2F::Class())) {
547 ((TH2F*)obj)->Fill(fComputed[0], fComputed[1]);
548 } else if (obj->InheritsFrom(TH3F::Class())) {
549 ((TH3F*)obj)->Fill(fComputed[0], fComputed[1], fComputed[2]);
550 } else if (obj->InheritsFrom(THnSparseF::Class())) {
551 ((THnSparseF*)obj)->Fill(fComputed.GetArray());
552 } else {
553 AliError("No output initialized");
03d23846 554 }
03d23846 555}