#include "Riostream.h"
#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include "TList.h"
#include "TMath.h"
#include "THnSparse.h"
#include "TString.h"
#include "TClonesArray.h"
#include "AliRsnMiniParticle.h"
#include "AliRsnMiniPair.h"
#include "AliRsnMiniEvent.h"
#include "AliAODEvent.h"
#include "AliLog.h"
#include "AliRsnCutSet.h"
#include "AliRsnMiniAxis.h"
#include "AliRsnMiniOutput.h"
#include "AliRsnMiniValue.h"
ClassImp(AliRsnMiniOutput)
AliRsnMiniOutput::AliRsnMiniOutput() :
TNamed(),
fOutputType(kTypes),
fComputation(kComputations),
fMotherPDG(0),
fMotherMass(0.0),
fPairCuts(0x0),
fOutputID(-1),
fAxes("AliRsnMiniAxis", 0),
fComputed(0),
fPair(),
fList(0x0),
fSel1(0),
fSel2(0),
fMaxNSisters(-1),
fCheckP(kFALSE),
fCheckFeedDown(kFALSE),
fOriginDselection(kFALSE),
fKeepDfromB(kFALSE),
fKeepDfromBOnly(kFALSE),
fRejectIfNoQuark(kFALSE),
fCheckHistRange(kTRUE)
{
fCutID[0] = fCutID[1] = -1;
fDaughter[0] = fDaughter[1] = AliRsnDaughter::kUnknown;
fCharge[0] = fCharge[1] = 0;
}
AliRsnMiniOutput::AliRsnMiniOutput(const char *name, EOutputType type, EComputation src) :
TNamed(name, ""),
fOutputType(type),
fComputation(src),
fMotherPDG(0),
fMotherMass(0.0),
fPairCuts(0x0),
fOutputID(-1),
fAxes("AliRsnMiniAxis", 0),
fComputed(0),
fPair(),
fList(0x0),
fSel1(0),
fSel2(0),
fMaxNSisters(-1),
fCheckP(kFALSE),
fCheckFeedDown(kFALSE),
fOriginDselection(kFALSE),
fKeepDfromB(kFALSE),
fKeepDfromBOnly(kFALSE),
fRejectIfNoQuark(kFALSE),
fCheckHistRange(kTRUE)
{
fCutID[0] = fCutID[1] = -1;
fDaughter[0] = fDaughter[1] = AliRsnDaughter::kUnknown;
fCharge[0] = fCharge[1] = 0;
}
AliRsnMiniOutput::AliRsnMiniOutput(const char *name, const char *outType, const char *compType) :
TNamed(name, ""),
fOutputType(kTypes),
fComputation(kComputations),
fMotherPDG(0),
fMotherMass(0.0),
fPairCuts(0x0),
fOutputID(-1),
fAxes("AliRsnMiniAxis", 0),
fComputed(0),
fPair(),
fList(0x0),
fSel1(0),
fSel2(0),
fMaxNSisters(-1),
fCheckP(kFALSE),
fCheckFeedDown(kFALSE),
fOriginDselection(kFALSE),
fKeepDfromB(kFALSE),
fKeepDfromBOnly(kFALSE),
fRejectIfNoQuark(kFALSE),
fCheckHistRange(kTRUE)
{
TString input;
input = outType;
input.ToUpper();
if (!input.CompareTo("HIST"))
fOutputType = kHistogram;
else if (!input.CompareTo("SPARSE"))
fOutputType = kHistogramSparse;
else
AliWarning(Form("String '%s' does not define a meaningful output type", outType));
input = compType;
input.ToUpper();
if (!input.CompareTo("EVENT"))
fComputation = kEventOnly;
else if (!input.CompareTo("PAIR"))
fComputation = kTrackPair;
else if (!input.CompareTo("MIX"))
fComputation = kTrackPairMix;
else if (!input.CompareTo("ROTATE1"))
fComputation = kTrackPairRotated1;
else if (!input.CompareTo("ROTATE2"))
fComputation = kTrackPairRotated2;
else if (!input.CompareTo("TRUE"))
fComputation = kTruePair;
else if (!input.CompareTo("MOTHER"))
fComputation = kMother;
else if (!input.CompareTo("MOTHER_IN_ACC"))
fComputation = kMotherInAcc;
else
AliWarning(Form("String '%s' does not define a meaningful computation type", compType));
fCutID[0] = fCutID[1] = -1;
fDaughter[0] = fDaughter[1] = AliRsnDaughter::kUnknown;
fCharge[0] = fCharge[1] = 0;
}
AliRsnMiniOutput::AliRsnMiniOutput(const AliRsnMiniOutput ©) :
TNamed(copy),
fOutputType(copy.fOutputType),
fComputation(copy.fComputation),
fMotherPDG(copy.fMotherPDG),
fMotherMass(copy.fMotherMass),
fPairCuts(copy.fPairCuts),
fOutputID(copy.fOutputID),
fAxes(copy.fAxes),
fComputed(copy.fComputed),
fPair(),
fList(copy.fList),
fSel1(0),
fSel2(0),
fMaxNSisters(-1),
fCheckP(kFALSE),
fCheckFeedDown(kFALSE),
fOriginDselection(kFALSE),
fKeepDfromB(kFALSE),
fKeepDfromBOnly(kFALSE),
fRejectIfNoQuark(kFALSE),
fCheckHistRange(copy.fCheckHistRange)
{
Int_t i;
for (i = 0; i < 2; i++) {
fCutID[i] = copy.fCutID[i];
fDaughter[i] = copy.fDaughter[i];
fCharge[i] = copy.fCharge[i];
}
}
AliRsnMiniOutput &AliRsnMiniOutput::operator=(const AliRsnMiniOutput ©)
{
if (this == ©)
return *this;
fOutputType = copy.fOutputType;
fComputation = copy.fComputation;
fMotherPDG = copy.fMotherPDG;
fMotherMass = copy.fMotherMass;
fPairCuts = copy.fPairCuts;
fOutputID = copy.fOutputID;
fAxes = copy.fAxes;
fComputed = copy.fComputed;
fList = copy.fList;
Int_t i;
for (i = 0; i < 2; i++) {
fCutID[i] = copy.fCutID[i];
fDaughter[i] = copy.fDaughter[i];
fCharge[i] = copy.fCharge[i];
}
fSel1.Set(0);
fSel2.Set(0);
fMaxNSisters = copy.fMaxNSisters;
fCheckP = copy.fCheckP;
fCheckFeedDown = copy.fCheckFeedDown;
fOriginDselection = copy.fOriginDselection;
fKeepDfromB = copy.fOriginDselection;
fKeepDfromBOnly = copy.fKeepDfromBOnly;
fRejectIfNoQuark = copy.fRejectIfNoQuark;
fCheckHistRange = copy.fCheckHistRange;
return (*this);
}
void AliRsnMiniOutput::AddAxis(Int_t i, Int_t nbins, Double_t min, Double_t max)
{
Int_t size = fAxes.GetEntries();
new (fAxes[size]) AliRsnMiniAxis(i, nbins, min, max);
}
void AliRsnMiniOutput::AddAxis(Int_t i, Double_t min, Double_t max, Double_t step)
{
Int_t size = fAxes.GetEntries();
new (fAxes[size]) AliRsnMiniAxis(i, min, max, step);
}
void AliRsnMiniOutput::AddAxis(Int_t i, Int_t nbins, Double_t *values)
{
Int_t size = fAxes.GetEntries();
new (fAxes[size]) AliRsnMiniAxis(i, nbins, values);
}
Bool_t AliRsnMiniOutput::Init(const char *prefix, TList *list)
{
if (!list) {
AliError("Required an output list");
return kFALSE;
}
fList = list;
Int_t size = fAxes.GetEntries();
if (size < 1) {
AliWarning(Form("[%s] Cannot initialize histogram with less than 1 axis", GetName()));
return kFALSE;
}
switch (fOutputType) {
case kHistogram:
if (size <= 3) {
CreateHistogram(Form("%s_%s", prefix, GetName()));
} else {
AliInfo(Form("[%s] Added %d > 3 axes. Creating a sparse histogram", GetName(), size));
fOutputType = kHistogramSparse;
CreateHistogramSparse(Form("%s_%s", prefix, GetName()));
}
return kTRUE;
case kHistogramSparse:
CreateHistogramSparse(Form("%s_%s", prefix, GetName()));
return kTRUE;
default:
AliError("Wrong output histogram definition");
return kFALSE;
}
}
void AliRsnMiniOutput::CreateHistogram(const char *name)
{
Int_t size = fAxes.GetEntries();
AliInfo(Form("Histogram name = '%s', with %d axes", name, size));
AliRsnMiniAxis *xAxis = 0x0, *yAxis = 0x0, *zAxis = 0x0;
if (size >= 1) xAxis = (AliRsnMiniAxis *)fAxes[0];
if (size >= 2) yAxis = (AliRsnMiniAxis *)fAxes[1];
if (size >= 3) zAxis = (AliRsnMiniAxis *)fAxes[2];
TH1 *h1 = 0x0;
if (xAxis && yAxis && zAxis) {
h1 = new TH3F(name, "", xAxis->NBins(), xAxis->BinArray(), yAxis->NBins(), yAxis->BinArray(), zAxis->NBins(), zAxis->BinArray());
} else if (xAxis && yAxis) {
h1 = new TH2F(name, "", xAxis->NBins(), xAxis->BinArray(), yAxis->NBins(), yAxis->BinArray());
} else if (xAxis) {
h1 = new TH1F(name, "", xAxis->NBins(), xAxis->BinArray());
} else {
AliError("No axis was initialized");
return;
}
if (h1 && fList) {
h1->Sumw2();
fList->Add(h1);
fOutputID = fList->IndexOf(h1);
}
}
void AliRsnMiniOutput::CreateHistogramSparse(const char *name)
{
Int_t size = fAxes.GetEntries();
AliInfo(Form("Sparse histogram name = '%s', with %d axes", name, size));
Int_t i, *nbins = new Int_t[size];
for (i = 0; i < size; i++) {
AliRsnMiniAxis *axis = (AliRsnMiniAxis *)fAxes[i];
nbins[i] = axis->NBins();
}
THnSparseF *h1 = new THnSparseF(name, "", size, nbins);
for (i = 0; i < size; i++) {
AliRsnMiniAxis *axis = (AliRsnMiniAxis *)fAxes[i];
h1->GetAxis(i)->Set(nbins[i], axis->BinArray());
}
delete [] nbins;
if (h1 && fList) {
h1->Sumw2();
fList->Add(h1);
fOutputID = fList->IndexOf(h1);
}
}
Bool_t AliRsnMiniOutput::FillEvent(AliRsnMiniEvent *event, TClonesArray *valueList)
{
if (fComputation != kEventOnly) {
AliError("This method can be called only for event-based computations");
return kFALSE;
}
ComputeValues(event, valueList);
FillHistogram();
return kTRUE;
}
Bool_t AliRsnMiniOutput::FillMother(const AliRsnMiniPair *pair, AliRsnMiniEvent *event, TClonesArray *valueList)
{
if (fComputation != kMother) {
AliError("This method can be called only for mother-based computations");
return kFALSE;
}
fPair = (*pair);
if (fPairCuts) if (!fPairCuts->IsSelected(&fPair)) return kFALSE;
ComputeValues(event, valueList);
FillHistogram();
return kTRUE;
}
Bool_t AliRsnMiniOutput::FillMotherInAcceptance(const AliRsnMiniPair *pair, AliRsnMiniEvent *event, TClonesArray *valueList)
{
if (fComputation != kMotherInAcc) {
AliError("This method can be called only for mother-based computations");
return kFALSE;
}
fPair = (*pair);
if (fPairCuts) if (!fPairCuts->IsSelected(&fPair)) return kFALSE;
ComputeValues(event, valueList);
FillHistogram();
return kTRUE;
}
Int_t AliRsnMiniOutput::FillPair(AliRsnMiniEvent *event1, AliRsnMiniEvent *event2, TClonesArray *valueList, Bool_t refFirst)
{
Bool_t okComp = kFALSE;
if (fComputation == kTrackPair) okComp = kTRUE;
if (fComputation == kTrackPairMix) okComp = kTRUE;
if (fComputation == kTrackPairRotated1) okComp = kTRUE;
if (fComputation == kTrackPairRotated2) okComp = kTRUE;
if (fComputation == kTruePair) okComp = kTRUE;
if (!okComp) {
AliError(Form("[%s] This method can be called only for pair-based computations", GetName()));
return kFALSE;
}
Int_t i1, i2, start, nadded = 0;
AliRsnMiniParticle *p1, *p2;
Bool_t sameCriteria = ((fCharge[0] == fCharge[1]) && (fDaughter[0] == fDaughter[1]));
Bool_t sameEvent = (event1->ID() == event2->ID());
TString selList1 = "";
TString selList2 = "";
Int_t n1 = event1->CountParticles(fSel1, fCharge[0], fCutID[0]);
Int_t n2 = event2->CountParticles(fSel2, fCharge[1], fCutID[1]);
for (i1 = 0; i1 < n1; i1++) selList1.Append(Form("%d ", fSel1[i1]));
for (i2 = 0; i2 < n2; i2++) selList2.Append(Form("%d ", fSel2[i2]));
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()));
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()));
if (!n1 || !n2) {
AliDebugClass(1, "No pairs to mix");
return 0;
}
for (i1 = 0; i1 < n1; i1++) {
p1 = event1->GetParticle(fSel1[i1]);
start = ((sameEvent && sameCriteria) ? i1 + 1 : 0);
AliDebugClass(2, Form("Start point = %d", start));
for (i2 = start; i2 < n2; i2++) {
p2 = event2->GetParticle(fSel2[i2]);
if (sameEvent && (p1->Index() == p2->Index())) {
AliDebugClass(2, "Skipping same index");
continue;
}
fPair.Fill(p1, p2, GetMass(0), GetMass(1), fMotherMass);
if (fComputation == kTrackPairRotated1) fPair.InvertP(kTRUE);
if (fComputation == kTrackPairRotated2) fPair.InvertP(kFALSE);
if (fComputation == kTruePair) {
if (fPair.Mother() < 0) {
continue;
} else if (fPair.MotherPDG() != fMotherPDG) {
continue;
}
Bool_t decayMatch = kFALSE;
if (p1->PDGAbs() == AliRsnDaughter::SpeciesPDG(fDaughter[0]) && p2->PDGAbs() == AliRsnDaughter::SpeciesPDG(fDaughter[1]))
decayMatch = kTRUE;
if (p2->PDGAbs() == AliRsnDaughter::SpeciesPDG(fDaughter[0]) && p1->PDGAbs() == AliRsnDaughter::SpeciesPDG(fDaughter[1]))
decayMatch = kTRUE;
if (!decayMatch) continue;
if ( (fMaxNSisters>0) && (p1->NTotSisters()==p2->NTotSisters()) && (p1->NTotSisters()>fMaxNSisters)) continue;
if ( fCheckP &&(TMath::Abs(fPair.PmotherX()-(p1->Px(1)+p2->Px(1)))/(TMath::Abs(fPair.PmotherX())+1.e-13)) > 0.00001 &&
(TMath::Abs(fPair.PmotherY()-(p1->Py(1)+p2->Py(1)))/(TMath::Abs(fPair.PmotherY())+1.e-13)) > 0.00001 &&
(TMath::Abs(fPair.PmotherZ()-(p1->Pz(1)+p2->Pz(1)))/(TMath::Abs(fPair.PmotherZ())+1.e-13)) > 0.00001 ) continue;
if ( fCheckFeedDown ){
Int_t pdgGranma = 0;
Bool_t isFromB=kFALSE;
Bool_t isQuarkFound=kFALSE;
if(fPair.IsFromB() == kTRUE) isFromB = kTRUE;
if(fPair.IsQuarkFound() == kTRUE) isQuarkFound = kTRUE;
if(fRejectIfNoQuark && !isQuarkFound) pdgGranma = -99999;
if(isFromB){
if (!fKeepDfromB) pdgGranma = -9999;
}
else{
if (fKeepDfromBOnly) pdgGranma = -999;
}
if (pdgGranma == -99999){
AliDebug(2,"This particle does not have a quark in his genealogy\n");
continue;
}
if (pdgGranma == -9999){
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");
continue;
}
if (pdgGranma == -999){
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");
continue;
}
}
}
if (fPairCuts) {
if (!fPairCuts->IsSelected(&fPair)) continue;
}
nadded++;
if (refFirst) ComputeValues(event1, valueList); else ComputeValues(event2, valueList);
FillHistogram();
}
}
AliDebugClass(1, Form("Pairs added in total = %4d", nadded));
return nadded;
}
void AliRsnMiniOutput::SetDselection(UShort_t originDselection)
{
fOriginDselection = originDselection;
if (fOriginDselection == 0) {
fKeepDfromB = kFALSE;
fKeepDfromBOnly = kFALSE;
}
if (fOriginDselection == 1) {
fKeepDfromB = kTRUE;
fKeepDfromBOnly = kTRUE;
}
if (fOriginDselection == 2) {
fKeepDfromB = kTRUE;
fKeepDfromBOnly = kFALSE;
}
return;
}
void AliRsnMiniOutput::ComputeValues(AliRsnMiniEvent *event, TClonesArray *valueList)
{
Int_t size = fAxes.GetEntries();
if (fComputed.GetSize() != size) fComputed.Set(size);
Int_t i, ival, nval = valueList->GetEntries();
for (i = 0; i < size; i++) {
fComputed[i] = 1E20;
AliRsnMiniAxis *axis = (AliRsnMiniAxis *)fAxes[i];
if (!axis) {
AliError("Null axis");
continue;
}
ival = axis->GetValueID();
if (ival < 0 || ival >= nval) {
AliError(Form("Required value #%d, while maximum is %d", ival, nval));
continue;
}
AliRsnMiniValue *val = (AliRsnMiniValue *)valueList->At(ival);
if (!val) {
AliError(Form("Value in position #%d is NULL", ival));
continue;
}
fComputed[i] = val->Eval(&fPair, event);
}
}
void AliRsnMiniOutput::FillHistogram()
{
if (!fList) {
AliError("List pointer is NULL");
return;
}
TObject *obj = fList->At(fOutputID);
if (obj->InheritsFrom(TH1F::Class())) {
((TH1F *)obj)->Fill(fComputed[0]);
} else if (obj->InheritsFrom(TH2F::Class())) {
((TH2F *)obj)->Fill(fComputed[0], fComputed[1]);
} else if (obj->InheritsFrom(TH3F::Class())) {
((TH3F *)obj)->Fill(fComputed[0], fComputed[1], fComputed[2]);
} else if (obj->InheritsFrom(THnSparseF::Class())) {
THnSparseF *h = (THnSparseF *)obj;
if (fCheckHistRange) {
for (Int_t iAxis = 0; iAxis<h->GetNdimensions(); iAxis++) {
if (fComputed.At(iAxis)>h->GetAxis(iAxis)->GetXmax() || fComputed.At(iAxis)<h->GetAxis(iAxis)->GetXmin()) return;
}
}
h->Fill(fComputed.GetArray());
} else {
AliError("No output initialized");
}
}