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