3 // All the definitions needed for building a RSN histogram
5 // -- properties of resonance (mass, PDG code if needed)
6 // -- properties of daughters (assigned mass, charges)
7 // -- definition of output histogram
10 #include "Riostream.h"
17 #include "THnSparse.h"
19 #include "TClonesArray.h"
21 #include "AliRsnMiniParticle.h"
22 #include "AliRsnMiniPair.h"
23 #include "AliRsnMiniEvent.h"
24 #include "AliAODEvent.h"
27 #include "AliRsnCutSet.h"
28 #include "AliRsnMiniAxis.h"
29 #include "AliRsnMiniOutput.h"
30 #include "AliRsnMiniValue.h"
32 ClassImp(AliRsnMiniOutput)
34 //__________________________________________________________________________________________________
35 AliRsnMiniOutput::AliRsnMiniOutput() :
38 fComputation(kComputations),
43 fAxes("AliRsnMiniAxis", 0),
51 fCheckFeedDown(kFALSE),
52 fOriginDselection(kFALSE),
54 fKeepDfromBOnly(kFALSE),
55 fRejectIfNoQuark(kFALSE),
56 fCheckHistRange(kTRUE)
62 fCutID[0] = fCutID[1] = -1;
63 fDaughter[0] = fDaughter[1] = AliRsnDaughter::kUnknown;
64 fCharge[0] = fCharge[1] = 0;
67 //__________________________________________________________________________________________________
68 AliRsnMiniOutput::AliRsnMiniOutput(const char *name, EOutputType type, EComputation src) :
76 fAxes("AliRsnMiniAxis", 0),
84 fCheckFeedDown(kFALSE),
85 fOriginDselection(kFALSE),
87 fKeepDfromBOnly(kFALSE),
88 fRejectIfNoQuark(kFALSE),
89 fCheckHistRange(kTRUE)
95 fCutID[0] = fCutID[1] = -1;
96 fDaughter[0] = fDaughter[1] = AliRsnDaughter::kUnknown;
97 fCharge[0] = fCharge[1] = 0;
100 //__________________________________________________________________________________________________
101 AliRsnMiniOutput::AliRsnMiniOutput(const char *name, const char *outType, const char *compType) :
104 fComputation(kComputations),
109 fAxes("AliRsnMiniAxis", 0),
117 fCheckFeedDown(kFALSE),
118 fOriginDselection(kFALSE),
120 fKeepDfromBOnly(kFALSE),
121 fRejectIfNoQuark(kFALSE),
122 fCheckHistRange(kTRUE)
125 // Constructor, with a more user friendly implementation, where
126 // the user sets the type of output and computations through conventional strings:
129 // -- "HIST" --> common histogram (up to 3 dimensions)
130 // -- "SPARSE" --> sparse histogram
133 // -- "EVENT" --> event-only computations
134 // -- "PAIR" --> track pair computations (default)
135 // -- "MIX" --> event mixing (like track pair, but different events)
136 // -- "ROTATE1" --> rotated background (rotate first track)
137 // -- "ROTATE2" --> rotated background (rotate second track)
138 // -- "TRUE" --> true pairs (like track pair, but checking that come from same mother)
139 // -- "MOTHER" --> mother (loop on MC directly for mothers --> denominator of efficiency)
140 // -- "MOTHER_IN_ACC" --> mother (loop on MC directly for mothers (in a defined acceptance interval)--> needed for efficiency calcutation using an enriched sample)
145 // understand output type
148 if (!input.CompareTo("HIST"))
149 fOutputType = kHistogram;
150 else if (!input.CompareTo("SPARSE"))
151 fOutputType = kHistogramSparse;
153 AliWarning(Form("String '%s' does not define a meaningful output type", outType));
155 // understand computation type
158 if (!input.CompareTo("EVENT"))
159 fComputation = kEventOnly;
160 else if (!input.CompareTo("PAIR"))
161 fComputation = kTrackPair;
162 else if (!input.CompareTo("MIX"))
163 fComputation = kTrackPairMix;
164 else if (!input.CompareTo("ROTATE1"))
165 fComputation = kTrackPairRotated1;
166 else if (!input.CompareTo("ROTATE2"))
167 fComputation = kTrackPairRotated2;
168 else if (!input.CompareTo("TRUE"))
169 fComputation = kTruePair;
170 else if (!input.CompareTo("MOTHER"))
171 fComputation = kMother;
172 else if (!input.CompareTo("MOTHER_IN_ACC"))
173 fComputation = kMotherInAcc;
175 AliWarning(Form("String '%s' does not define a meaningful computation type", compType));
177 fCutID[0] = fCutID[1] = -1;
178 fDaughter[0] = fDaughter[1] = AliRsnDaughter::kUnknown;
179 fCharge[0] = fCharge[1] = 0;
182 //__________________________________________________________________________________________________
183 AliRsnMiniOutput::AliRsnMiniOutput(const AliRsnMiniOutput ©) :
185 fOutputType(copy.fOutputType),
186 fComputation(copy.fComputation),
187 fMotherPDG(copy.fMotherPDG),
188 fMotherMass(copy.fMotherMass),
189 fPairCuts(copy.fPairCuts),
190 fOutputID(copy.fOutputID),
192 fComputed(copy.fComputed),
199 fCheckFeedDown(kFALSE),
200 fOriginDselection(kFALSE),
202 fKeepDfromBOnly(kFALSE),
203 fRejectIfNoQuark(kFALSE),
204 fCheckHistRange(copy.fCheckHistRange)
211 for (i = 0; i < 2; i++) {
212 fCutID[i] = copy.fCutID[i];
213 fDaughter[i] = copy.fDaughter[i];
214 fCharge[i] = copy.fCharge[i];
218 //__________________________________________________________________________________________________
219 AliRsnMiniOutput &AliRsnMiniOutput::operator=(const AliRsnMiniOutput ©)
222 // Assignment operator
226 fOutputType = copy.fOutputType;
227 fComputation = copy.fComputation;
228 fMotherPDG = copy.fMotherPDG;
229 fMotherMass = copy.fMotherMass;
230 fPairCuts = copy.fPairCuts;
231 fOutputID = copy.fOutputID;
233 fComputed = copy.fComputed;
237 for (i = 0; i < 2; i++) {
238 fCutID[i] = copy.fCutID[i];
239 fDaughter[i] = copy.fDaughter[i];
240 fCharge[i] = copy.fCharge[i];
245 fMaxNSisters = copy.fMaxNSisters;
246 fCheckP = copy.fCheckP;
247 fCheckFeedDown = copy.fCheckFeedDown;
248 fOriginDselection = copy.fOriginDselection;
249 fKeepDfromB = copy.fOriginDselection;
250 fKeepDfromBOnly = copy.fKeepDfromBOnly;
251 fRejectIfNoQuark = copy.fRejectIfNoQuark;
252 fCheckHistRange = copy.fCheckHistRange;
258 //__________________________________________________________________________________________________
259 void AliRsnMiniOutput::AddAxis(Int_t i, Int_t nbins, Double_t min, Double_t max)
262 // Create a new axis reference
265 Int_t size = fAxes.GetEntries();
266 new (fAxes[size]) AliRsnMiniAxis(i, nbins, min, max);
269 //__________________________________________________________________________________________________
270 void AliRsnMiniOutput::AddAxis(Int_t i, Double_t min, Double_t max, Double_t step)
273 // Create a new axis reference
276 Int_t size = fAxes.GetEntries();
277 new (fAxes[size]) AliRsnMiniAxis(i, min, max, step);
280 //__________________________________________________________________________________________________
281 void AliRsnMiniOutput::AddAxis(Int_t i, Int_t nbins, Double_t *values)
284 // Create a new axis reference
287 Int_t size = fAxes.GetEntries();
288 new (fAxes[size]) AliRsnMiniAxis(i, nbins, values);
291 //__________________________________________________________________________________________________
292 Bool_t AliRsnMiniOutput::Init(const char *prefix, TList *list)
295 // Initialize properly the histogram and add it to the argument list
299 AliError("Required an output list");
304 Int_t size = fAxes.GetEntries();
306 AliWarning(Form("[%s] Cannot initialize histogram with less than 1 axis", GetName()));
310 switch (fOutputType) {
313 CreateHistogram(Form("%s_%s", prefix, GetName()));
315 AliInfo(Form("[%s] Added %d > 3 axes. Creating a sparse histogram", GetName(), size));
316 fOutputType = kHistogramSparse;
317 CreateHistogramSparse(Form("%s_%s", prefix, GetName()));
320 case kHistogramSparse:
321 CreateHistogramSparse(Form("%s_%s", prefix, GetName()));
324 AliError("Wrong output histogram definition");
329 //__________________________________________________________________________________________________
330 void AliRsnMiniOutput::CreateHistogram(const char *name)
333 // Initialize the 'default' TH1 output object.
334 // In case one of the expected axes is NULL, the initialization fails.
337 Int_t size = fAxes.GetEntries();
338 AliInfo(Form("Histogram name = '%s', with %d axes", name, size));
340 // we expect to have maximum 3 axes in this case
341 AliRsnMiniAxis *xAxis = 0x0, *yAxis = 0x0, *zAxis = 0x0;
342 if (size >= 1) xAxis = (AliRsnMiniAxis *)fAxes[0];
343 if (size >= 2) yAxis = (AliRsnMiniAxis *)fAxes[1];
344 if (size >= 3) zAxis = (AliRsnMiniAxis *)fAxes[2];
346 // create histogram depending on the number of axes
348 if (xAxis && yAxis && zAxis) {
349 h1 = new TH3F(name, "", xAxis->NBins(), xAxis->BinArray(), yAxis->NBins(), yAxis->BinArray(), zAxis->NBins(), zAxis->BinArray());
350 } else if (xAxis && yAxis) {
351 h1 = new TH2F(name, "", xAxis->NBins(), xAxis->BinArray(), yAxis->NBins(), yAxis->BinArray());
353 h1 = new TH1F(name, "", xAxis->NBins(), xAxis->BinArray());
355 AliError("No axis was initialized");
359 // switch the correct computation of errors
363 fOutputID = fList->IndexOf(h1);
367 //________________________________________________________________________________________
368 void AliRsnMiniOutput::CreateHistogramSparse(const char *name)
371 // Initialize the THnSparse output object.
372 // In case one of the expected axes is NULL, the initialization fails.
375 Int_t size = fAxes.GetEntries();
376 AliInfo(Form("Sparse histogram name = '%s', with %d axes", name, size));
378 // retrieve binnings and sizes of all axes
379 // since the check for null values is done in Init(),
380 // we assume that here they must all be well defined
381 Int_t i, *nbins = new Int_t[size];
382 for (i = 0; i < size; i++) {
383 AliRsnMiniAxis *axis = (AliRsnMiniAxis *)fAxes[i];
384 nbins[i] = axis->NBins();
387 // create fHSparseogram
388 THnSparseF *h1 = new THnSparseF(name, "", size, nbins);
390 // update the various axes using the definitions given in the array of axes here
391 for (i = 0; i < size; i++) {
392 AliRsnMiniAxis *axis = (AliRsnMiniAxis *)fAxes[i];
393 h1->GetAxis(i)->Set(nbins[i], axis->BinArray());
403 fOutputID = fList->IndexOf(h1);
407 //________________________________________________________________________________________
408 Bool_t AliRsnMiniOutput::FillEvent(AliRsnMiniEvent *event, TClonesArray *valueList)
411 // Compute values for event-based computations (does not use the pair)
414 // check computation type
415 if (fComputation != kEventOnly) {
416 AliError("This method can be called only for event-based computations");
421 ComputeValues(event, valueList);
426 //________________________________________________________________________________________
427 Bool_t AliRsnMiniOutput::FillMother(const AliRsnMiniPair *pair, AliRsnMiniEvent *event, TClonesArray *valueList)
430 // Compute values for mother-based computations
433 // check computation type
434 if (fComputation != kMother) {
435 AliError("This method can be called only for mother-based computations");
439 // copy passed pair info
442 // check pair against cuts
443 if (fPairCuts) if (!fPairCuts->IsSelected(&fPair)) return kFALSE;
446 ComputeValues(event, valueList);
451 //________________________________________________________________________________________
452 Bool_t AliRsnMiniOutput::FillMotherInAcceptance(const AliRsnMiniPair *pair, AliRsnMiniEvent *event, TClonesArray *valueList)
455 // Compute values for mother-based computations
458 // check computation type
459 if (fComputation != kMotherInAcc) {
460 AliError("This method can be called only for mother-based computations");
464 // copy passed pair info
467 // check pair against cuts
468 if (fPairCuts) if (!fPairCuts->IsSelected(&fPair)) return kFALSE;
471 ComputeValues(event, valueList);
476 //________________________________________________________________________________________
477 Int_t AliRsnMiniOutput::FillPair(AliRsnMiniEvent *event1, AliRsnMiniEvent *event2, TClonesArray *valueList, Bool_t refFirst)
480 // Loops on the passed mini-event, and for each pair of particles
481 // which satisfy the charge and cut requirements defined here, add an entry.
482 // Returns the number of successful fillings.
483 // Last argument tells if the reference event for event-based values is the first or the second.
486 // check computation type
487 Bool_t okComp = kFALSE;
488 if (fComputation == kTrackPair) okComp = kTRUE;
489 if (fComputation == kTrackPairMix) okComp = kTRUE;
490 if (fComputation == kTrackPairRotated1) okComp = kTRUE;
491 if (fComputation == kTrackPairRotated2) okComp = kTRUE;
492 if (fComputation == kTruePair) okComp = kTRUE;
494 AliError(Form("[%s] This method can be called only for pair-based computations", GetName()));
499 Int_t i1, i2, start, nadded = 0;
500 AliRsnMiniParticle *p1, *p2;
502 // it is necessary to know if criteria for the two daughters are the same
503 // and if the two events are the same or not (mixing)
504 //Bool_t sameCriteria = ((fCharge[0] == fCharge[1]) && (fCutID[0] == fCutID[1]));
505 Bool_t sameCriteria = ((fCharge[0] == fCharge[1]) && (fDaughter[0] == fDaughter[1]));
506 Bool_t sameEvent = (event1->ID() == event2->ID());
508 TString selList1 = "";
509 TString selList2 = "";
510 Int_t n1 = event1->CountParticles(fSel1, fCharge[0], fCutID[0]);
511 Int_t n2 = event2->CountParticles(fSel2, fCharge[1], fCutID[1]);
512 for (i1 = 0; i1 < n1; i1++) selList1.Append(Form("%d ", fSel1[i1]));
513 for (i2 = 0; i2 < n2; i2++) selList2.Append(Form("%d ", fSel2[i2]));
514 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()));
515 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()));
517 AliDebugClass(1, "No pairs to mix");
522 for (i1 = 0; i1 < n1; i1++) {
523 p1 = event1->GetParticle(fSel1[i1]);
524 //p1 = event1->GetParticle(i1);
525 //if (p1->Charge() != fCharge[0]) continue;
526 //if (!p1->HasCutBit(fCutID[0])) continue;
527 // define starting point for inner loop
528 // if daughter selection criteria (charge, cuts) are the same
529 // and the two events coincide, internal loop must start from
530 // the first track *after* current one;
531 // otherwise it starts from the beginning
532 start = ((sameEvent && sameCriteria) ? i1 + 1 : 0);
533 AliDebugClass(2, Form("Start point = %d", start));
535 for (i2 = start; i2 < n2; i2++) {
536 p2 = event2->GetParticle(fSel2[i2]);
537 //p2 = event2->GetParticle(i2);
538 //if (p2->Charge() != fCharge[1]) continue;
539 //if (!p2->HasCutBit(fCutID[1])) continue;
540 // avoid to mix a particle with itself
541 if (sameEvent && (p1->Index() == p2->Index())) {
542 AliDebugClass(2, "Skipping same index");
546 fPair.Fill(p1, p2, GetMass(0), GetMass(1), fMotherMass);
547 // do rotation if needed
548 if (fComputation == kTrackPairRotated1) fPair.InvertP(kTRUE);
549 if (fComputation == kTrackPairRotated2) fPair.InvertP(kFALSE);
550 // if required, check that this is a true pair
551 if (fComputation == kTruePair) {
552 if (fPair.Mother() < 0) {
554 } else if (fPair.MotherPDG() != fMotherPDG) {
557 Bool_t decayMatch = kFALSE;
558 if (p1->PDGAbs() == AliRsnDaughter::SpeciesPDG(fDaughter[0]) && p2->PDGAbs() == AliRsnDaughter::SpeciesPDG(fDaughter[1]))
560 if (p2->PDGAbs() == AliRsnDaughter::SpeciesPDG(fDaughter[0]) && p1->PDGAbs() == AliRsnDaughter::SpeciesPDG(fDaughter[1]))
562 if (!decayMatch) continue;
563 if ( (fMaxNSisters>0) && (p1->NTotSisters()==p2->NTotSisters()) && (p1->NTotSisters()>fMaxNSisters)) continue;
564 if ( fCheckP &&(TMath::Abs(fPair.PmotherX()-(p1->Px(1)+p2->Px(1)))/(TMath::Abs(fPair.PmotherX())+1.e-13)) > 0.00001 &&
565 (TMath::Abs(fPair.PmotherY()-(p1->Py(1)+p2->Py(1)))/(TMath::Abs(fPair.PmotherY())+1.e-13)) > 0.00001 &&
566 (TMath::Abs(fPair.PmotherZ()-(p1->Pz(1)+p2->Pz(1)))/(TMath::Abs(fPair.PmotherZ())+1.e-13)) > 0.00001 ) continue;
567 if ( fCheckFeedDown ){
569 Bool_t isFromB=kFALSE;
570 Bool_t isQuarkFound=kFALSE;
572 if(fPair.IsFromB() == kTRUE) isFromB = kTRUE;
573 if(fPair.IsQuarkFound() == kTRUE) isQuarkFound = kTRUE;
574 if(fRejectIfNoQuark && !isQuarkFound) pdgGranma = -99999;
576 if (!fKeepDfromB) pdgGranma = -9999; //skip particle if come from a B meson.
579 if (fKeepDfromBOnly) pdgGranma = -999;
581 if (pdgGranma == -99999){
582 AliDebug(2,"This particle does not have a quark in his genealogy\n");
585 if (pdgGranma == -9999){
586 AliDebug(2,"This particle come from a B decay channel but according to the settings of the task, we keep only the prompt charm particles\n");
590 if (pdgGranma == -999){
591 AliDebug(2,"This particle come from a prompt charm particles but according to the settings of the task, we want only the ones coming from B\n");
596 // check pair against cuts
598 if (!fPairCuts->IsSelected(&fPair)) continue;
600 // get computed values & fill histogram
602 if (refFirst) ComputeValues(event1, valueList); else ComputeValues(event2, valueList);
604 } // end internal loop
605 } // end external loop
607 AliDebugClass(1, Form("Pairs added in total = %4d", nadded));
610 //___________________________________________________________
611 void AliRsnMiniOutput::SetDselection(UShort_t originDselection)
613 // setting the way the D0 will be selected
614 // 0 --> only from c quarks
615 // 1 --> only from b quarks
616 // 2 --> from both c quarks and b quarks
618 fOriginDselection = originDselection;
620 if (fOriginDselection == 0) {
621 fKeepDfromB = kFALSE;
622 fKeepDfromBOnly = kFALSE;
625 if (fOriginDselection == 1) {
627 fKeepDfromBOnly = kTRUE;
630 if (fOriginDselection == 2) {
632 fKeepDfromBOnly = kFALSE;
637 //________________________________________________________________________________________
638 void AliRsnMiniOutput::ComputeValues(AliRsnMiniEvent *event, TClonesArray *valueList)
641 // Using the arguments and the internal 'fPair' data member,
642 // compute all values to be stored in the histogram
645 // check size of computed array
646 Int_t size = fAxes.GetEntries();
647 if (fComputed.GetSize() != size) fComputed.Set(size);
649 Int_t i, ival, nval = valueList->GetEntries();
651 for (i = 0; i < size; i++) {
653 AliRsnMiniAxis *axis = (AliRsnMiniAxis *)fAxes[i];
655 AliError("Null axis");
658 ival = axis->GetValueID();
659 if (ival < 0 || ival >= nval) {
660 AliError(Form("Required value #%d, while maximum is %d", ival, nval));
663 AliRsnMiniValue *val = (AliRsnMiniValue *)valueList->At(ival);
665 AliError(Form("Value in position #%d is NULL", ival));
668 // if none of the above exit points is taken, compute value
669 fComputed[i] = val->Eval(&fPair, event);
673 //________________________________________________________________________________________
674 void AliRsnMiniOutput::FillHistogram()
677 // Fills the internal histogram using the current values stored in the
678 // 'fComputed' array, in the order as they are stored, up to the max
679 // dimension of the initialized histogram itself.
682 // retrieve object from list
684 AliError("List pointer is NULL");
687 TObject *obj = fList->At(fOutputID);
689 if (obj->InheritsFrom(TH1F::Class())) {
690 ((TH1F *)obj)->Fill(fComputed[0]);
691 } else if (obj->InheritsFrom(TH2F::Class())) {
692 ((TH2F *)obj)->Fill(fComputed[0], fComputed[1]);
693 } else if (obj->InheritsFrom(TH3F::Class())) {
694 ((TH3F *)obj)->Fill(fComputed[0], fComputed[1], fComputed[2]);
695 } else if (obj->InheritsFrom(THnSparseF::Class())) {
696 THnSparseF *h = (THnSparseF *)obj;
697 if (fCheckHistRange) {
698 for (Int_t iAxis = 0; iAxis<h->GetNdimensions(); iAxis++) {
699 if (fComputed.At(iAxis)>h->GetAxis(iAxis)->GetXmax() || fComputed.At(iAxis)<h->GetAxis(iAxis)->GetXmin()) return;
702 h->Fill(fComputed.GetArray());
704 AliError("No output initialized");