Coverity fix
[u/mrichter/AliRoot.git] / 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    // 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 //________________________________________________________________________________________
350 Bool_t AliRsnMiniOutput::FillEvent(AliRsnMiniEvent *event, TClonesArray *valueList)
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
362    // compute & fill
363    ComputeValues(event, valueList);
364    FillHistogram();
365    return kTRUE;
366 }
367
368 //________________________________________________________________________________________
369 Bool_t AliRsnMiniOutput::FillMother(const AliRsnMiniPair *pair, AliRsnMiniEvent *event, TClonesArray *valueList)
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
387    // compute & fill
388    ComputeValues(event, valueList);
389    FillHistogram();
390    return kTRUE;
391 }
392
393 //________________________________________________________________________________________
394 Int_t AliRsnMiniOutput::FillPair(AliRsnMiniEvent *event1, AliRsnMiniEvent *event2, TClonesArray *valueList)
395 {
396 //
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.
400 //
401
402    // check computation type
403    Bool_t okComp = kFALSE;
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;
409    if (!okComp) {
410       AliError(Form("[%s] This method can be called only for pair-based computations", GetName()));
411       return kFALSE;
412    }
413    
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;
419    
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());
424    
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    //}
433    
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
488    
489    return nadded;
490 }
491
492 //________________________________________________________________________________________
493 void AliRsnMiniOutput::ComputeValues(AliRsnMiniEvent *event, TClonesArray *valueList)
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    }
526 }
527
528 //________________________________________________________________________________________
529 void AliRsnMiniOutput::FillHistogram()
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");
540       return;
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");
554    }
555 }