#include <TVectorD.h>
#include <TH1.h>
#include <TH1F.h>
#include <TH2.h>
#include <TH3.h>
#include <TProfile.h>
#include <TProfile2D.h>
#include <TProfile3D.h>
#include <THnSparse.h>
#include <TAxis.h>
#include <TString.h>
#include <TObjString.h>
#include <TObjArray.h>
#include <AliVParticle.h>
#include <AliLog.h>
#include "AliDielectron.h"
#include "AliDielectronHelper.h"
#include "AliDielectronMC.h"
#include "AliDielectronPair.h"
#include "AliDielectronSignalMC.h"
#include "AliDielectronHistos.h"
#include "AliDielectronHF.h"
ClassImp(AliDielectronHF)
AliDielectronHF::AliDielectronHF() :
TNamed(),
fUsedVars(new TBits(AliDielectronVarManager::kNMaxValues)),
fArrPairType(),
fPairType(kSeOnlyOS),
fSignalsMC(0x0),
fVarCutType(new TBits(kMaxCuts)),
fAxes(kMaxCuts),
fHasMC(kFALSE),
fStepGenerated(kFALSE),
fEventArray(kFALSE),
fRefObj(1)
{
for (Int_t i=0; i<kMaxCuts; ++i){
fVarCuts[i]=0;
fBinType[i]=kStdBin;
}
fAxes.SetOwner(kTRUE);
fRefObj.SetOwner(kTRUE);
fArrPairType.SetOwner(kTRUE);
}
AliDielectronHF::AliDielectronHF(const char* name, const char* title) :
TNamed(name, title),
fUsedVars(new TBits(AliDielectronVarManager::kNMaxValues)),
fArrPairType(),
fPairType(kSeOnlyOS),
fSignalsMC(0x0),
fVarCutType(new TBits(kMaxCuts)),
fAxes(kMaxCuts),
fHasMC(kFALSE),
fStepGenerated(kFALSE),
fEventArray(kFALSE),
fRefObj(1)
{
for (Int_t i=0; i<kMaxCuts; ++i){
fVarCuts[i]=0;
fBinType[i]=kStdBin;
}
fAxes.SetOwner(kTRUE);
fRefObj.SetOwner(kTRUE);
fArrPairType.SetOwner(kTRUE);
}
AliDielectronHF::~AliDielectronHF()
{
if(fUsedVars) delete fUsedVars;
if(fVarCutType) delete fVarCutType;
fAxes.Delete();
fRefObj.Delete();
fArrPairType.Delete();
}
void AliDielectronHF::UserProfile(const char* histClass, UInt_t valTypeP,
const TVectorD * const binsX,
UInt_t valTypeX, TString option, UInt_t valTypeW)
{
TH1 *hist=0x0;
if(valTypeP==AliDielectronHistos::kNoProfile)
hist=new TH1F("","",binsX->GetNrows()-1,binsX->GetMatrixArray());
else {
TString opt=""; Double_t pmin=0., pmax=0.;
if(!option.IsNull()) {
TObjArray *arr=option.Tokenize(";");
arr->SetOwner();
opt=((TObjString*)arr->At(0))->GetString();
if(arr->GetEntriesFast()>1) pmin=(((TObjString*)arr->At(1))->GetString()).Atof();
if(arr->GetEntriesFast()>2) pmax=(((TObjString*)arr->At(2))->GetString()).Atof();
delete arr;
}
hist=new TProfile("","",binsX->GetNrows()-1,binsX->GetMatrixArray(),pmin,pmax,opt.Data());
}
UInt_t valType[4] = {0};
valType[0]=valTypeX; valType[1]=valTypeP;
AliDielectronHistos::StoreVariables(hist, valType);
hist->SetUniqueID(valTypeW);
for(Int_t i=0; i<4; i++) fUsedVars->SetBitNumber(valType[i],kTRUE);
fUsedVars->SetBitNumber(valTypeW,kTRUE);
AliDielectronHistos::AdaptNameTitle(hist, histClass);
hist->SetName(Form("HF_%s",hist->GetName()));
fRefObj.AddLast(hist);
delete binsX;
}
void AliDielectronHF::UserProfile(const char* histClass, UInt_t valTypeP,
const TVectorD * const binsX, const TVectorD * const binsY,
UInt_t valTypeX, UInt_t valTypeY, TString option, UInt_t valTypeW)
{
TH1 *hist=0x0;
if(valTypeP==AliDielectronHistos::kNoProfile) {
hist=new TH2F("","",
binsX->GetNrows()-1,binsX->GetMatrixArray(),
binsY->GetNrows()-1,binsY->GetMatrixArray());
}
else {
TString opt=""; Double_t pmin=0., pmax=0.;
if(!option.IsNull()) {
TObjArray *arr=option.Tokenize(";");
arr->SetOwner();
opt=((TObjString*)arr->At(0))->GetString();
if(arr->GetEntriesFast()>1) pmin=(((TObjString*)arr->At(1))->GetString()).Atof();
if(arr->GetEntriesFast()>2) pmax=(((TObjString*)arr->At(2))->GetString()).Atof();
delete arr;
}
hist=new TProfile2D("","",
binsX->GetNrows()-1,binsX->GetMatrixArray(),
binsY->GetNrows()-1,binsY->GetMatrixArray());
((TProfile2D*)hist)->BuildOptions(pmin,pmax,opt.Data());
}
UInt_t valType[4] = {0};
valType[0]=valTypeX; valType[1]=valTypeY; valType[2]=valTypeP;
AliDielectronHistos::StoreVariables(hist, valType);
hist->SetUniqueID(valTypeW);
for(Int_t i=0; i<4; i++) fUsedVars->SetBitNumber(valType[i],kTRUE);
fUsedVars->SetBitNumber(valTypeW,kTRUE);
AliDielectronHistos::AdaptNameTitle(hist, histClass);
hist->SetName(Form("HF_%s",hist->GetName()));
fRefObj.AddLast(hist);
delete binsX;
delete binsY;
}
void AliDielectronHF::UserProfile(const char* histClass, UInt_t valTypeP,
const TVectorD * const binsX, const TVectorD * const binsY, const TVectorD * const binsZ,
UInt_t valTypeX, UInt_t valTypeY, UInt_t valTypeZ, TString option, UInt_t valTypeW)
{
TH1 *hist=0x0;
if(valTypeP==AliDielectronHistos::kNoProfile) {
hist=new TH3F("","",
binsX->GetNrows()-1,binsX->GetMatrixArray(),
binsY->GetNrows()-1,binsY->GetMatrixArray(),
binsZ->GetNrows()-1,binsZ->GetMatrixArray());
}
else {
TString opt=""; Double_t pmin=0., pmax=0.;
if(!option.IsNull()) {
TObjArray *arr=option.Tokenize(";");
arr->SetOwner();
opt=((TObjString*)arr->At(0))->GetString();
if(arr->GetEntriesFast()>1) pmin=(((TObjString*)arr->At(1))->GetString()).Atof();
if(arr->GetEntriesFast()>2) pmax=(((TObjString*)arr->At(2))->GetString()).Atof();
delete arr;
}
hist=new TProfile3D("","",
binsX->GetNrows()-1,binsX->GetMatrixArray(),
binsY->GetNrows()-1,binsY->GetMatrixArray(),
binsZ->GetNrows()-1,binsZ->GetMatrixArray());
((TProfile3D*)hist)->BuildOptions(pmin,pmax,opt.Data());
}
UInt_t valType[4] = {0};
valType[0]=valTypeX; valType[1]=valTypeY; valType[2]=valTypeZ; valType[3]=valTypeP;
AliDielectronHistos::StoreVariables(hist, valType);
hist->SetUniqueID(valTypeW);
for(Int_t i=0; i<4; i++) fUsedVars->SetBitNumber(valType[i],kTRUE);
fUsedVars->SetBitNumber(valTypeW,kTRUE);
AliDielectronHistos::AdaptNameTitle(hist, histClass);
hist->SetName(Form("HF_%s",hist->GetName()));
fRefObj.AddLast(hist);
delete binsX;
delete binsY;
delete binsZ;
}
void AliDielectronHF::UserSparse(const char* histClass, Int_t ndim, TObjArray *limits, UInt_t *vars, UInt_t valTypeW)
{
THnSparseF *hist=0;
Int_t bins[ndim];
for(Int_t idim=0 ;idim<ndim; idim++) {
TVectorD *vec = (TVectorD*) limits->At(idim);
bins[idim]=vec->GetNrows()-1;
}
hist=new THnSparseF("",histClass, ndim, bins, 0x0, 0x0);
for(Int_t idim=0 ;idim<ndim; idim++) {
TVectorD *vec = (TVectorD*) limits->At(idim);
hist->SetBinEdges(idim,vec->GetMatrixArray());
}
AliDielectronHistos::StoreVariables(hist, vars);
hist->SetUniqueID(valTypeW);
for(Int_t i=0; i<20; i++) fUsedVars->SetBitNumber(vars[i],kTRUE);
fUsedVars->SetBitNumber(valTypeW,kTRUE);
TString name;
for(Int_t iv=0; iv < ndim; iv++) name+=Form("%s_",AliDielectronVarManager::GetValueName(vars[iv]));
name.Resize(name.Length()-1);
hist->SetName(Form("HF_%s",name.Data()));
fRefObj.AddLast(hist);
delete limits;
}
void AliDielectronHF::AddCutVariable(AliDielectronVarManager::ValueTypes type,
Int_t nbins, Double_t min, Double_t max, Bool_t log, Bool_t leg, EBinType btype)
{
if (fAxes.GetEntriesFast()>=kMaxCuts) return;
TVectorD *binLimits=0x0;
if (!log) binLimits=AliDielectronHelper::MakeLinBinning(nbins,min,max);
else binLimits=AliDielectronHelper::MakeLogBinning(nbins,min,max);
if (!binLimits) return;
Int_t size=fAxes.GetEntriesFast();
fVarCuts[size]=(UShort_t)type;
fVarCutType->SetBitNumber(size,leg);
fAxes.Add(binLimits->Clone());
fBinType[size]=btype;
fUsedVars->SetBitNumber(type,kTRUE);
}
void AliDielectronHF::AddCutVariable(AliDielectronVarManager::ValueTypes type,
const char* binLimitStr, Bool_t leg, EBinType btype)
{
if (fAxes.GetEntriesFast()>=kMaxCuts) return;
TVectorD *binLimits=AliDielectronHelper::MakeArbitraryBinning(binLimitStr);
if (!binLimits) return;
Int_t size=fAxes.GetEntriesFast();
fVarCuts[size]=(UShort_t)type;
fVarCutType->SetBitNumber(size,leg);
fAxes.Add(binLimits);
fBinType[size]=btype;
fUsedVars->SetBitNumber(type,kTRUE);
}
void AliDielectronHF::AddCutVariable(AliDielectronVarManager::ValueTypes type,
TVectorD * binLimits, Bool_t leg, EBinType btype)
{
if (fAxes.GetEntriesFast()>=kMaxCuts) return;
if (!binLimits) return;
Int_t size=fAxes.GetEntriesFast();
fVarCuts[size]=(UShort_t)type;
fVarCutType->SetBitNumber(size,leg);
fAxes.Add(binLimits);
fBinType[size]=btype;
fUsedVars->SetBitNumber(type,kTRUE);
}
void AliDielectronHF::Fill(Int_t label1, Int_t label2, Int_t nSignal)
{
if(!fStepGenerated || fEventArray) return;
AliVParticle* part1 = AliDielectronMC::Instance()->GetMCTrackFromMCEvent(label1);
AliVParticle* part2 = AliDielectronMC::Instance()->GetMCTrackFromMCEvent(label2);
if(!part1 || !part2) return;
AliDielectronMC* dieMC = AliDielectronMC::Instance();
Int_t mLabel1 = dieMC->GetMothersLabel(label1);
Int_t mLabel2 = dieMC->GetMothersLabel(label2);
AliDielectronSignalMC* sigMC = (AliDielectronSignalMC*)fSignalsMC->At(nSignal);
if(sigMC->GetMothersRelation()==AliDielectronSignalMC::kSame && mLabel1!=mLabel2) return;
if(sigMC->GetMothersRelation()==AliDielectronSignalMC::kDifferent && mLabel1==mLabel2) return;
AliDielectronVarManager::SetFillMap(fUsedVars);
Double_t valuesLeg1[AliDielectronVarManager::kNMaxValues];
Double_t valuesLeg2[AliDielectronVarManager::kNMaxValues];
AliDielectronVarManager::Fill(part1,valuesLeg1);
AliDielectronVarManager::Fill(part2,valuesLeg2);
Double_t valuesPair[AliDielectronVarManager::kNMaxValues];
AliDielectronVarManager::Fill(dieMC->GetMCEvent(), valuesPair);
AliDielectronVarManager::FillVarMCParticle2(part1,part2,valuesPair);
Int_t istep=0;
if(fPairType!=kMConly) istep=AliDielectron::kEv1PMRot+1;
if(part1->Charge()*part2->Charge()<0) {
Fill(istep+nSignal+fSignalsMC->GetEntries(), valuesPair, valuesLeg1, valuesLeg2);
}
return;
}
void AliDielectronHF::Fill(Int_t pairIndex, const AliDielectronPair *particle)
{
if(!IsPairTypeSelected(pairIndex) || fEventArray) return;
Double_t valuesPair[AliDielectronVarManager::kNMaxValues];
AliDielectronVarManager::SetFillMap(fUsedVars);
AliDielectronVarManager::Fill(particle,valuesPair);
Double_t valuesLeg1[AliDielectronVarManager::kNMaxValues]={0};
if(fVarCutType->CountBits()) AliDielectronVarManager::Fill(particle->GetFirstDaughterP(),valuesLeg1);
Double_t valuesLeg2[AliDielectronVarManager::kNMaxValues]={0};
if(fVarCutType->CountBits()) AliDielectronVarManager::Fill(particle->GetSecondDaughterP(),valuesLeg2);
Int_t istep = 0;
if(fPairType!=kMConly) istep=AliDielectron::kEv1PMRot+1;
if(fHasMC && fSignalsMC && pairIndex==AliDielectron::kEv1PM) {
for(Int_t i=0; i<fSignalsMC->GetEntries(); i++) {
if(AliDielectronMC::Instance()->IsMCTruth(particle, (AliDielectronSignalMC*)fSignalsMC->At(i)))
Fill(istep+i, valuesPair, valuesLeg1, valuesLeg2);
}
}
if(fPairType==kMConly) return;
Fill(pairIndex, valuesPair, valuesLeg1, valuesLeg2);
return;
}
void AliDielectronHF::Fill(Int_t index, Double_t * const valuesPair, Double_t * const valuesLeg1, Double_t * const valuesLeg2)
{
TObjArray *histArr = static_cast<TObjArray*>(fArrPairType.At(index));
if(!histArr) return;
Int_t size = GetNumberOfBins();
for(Int_t ihist=0; ihist<size; ihist++) {
Int_t sizeAdd = 1;
Bool_t selected = kTRUE;
Int_t nvars = fAxes.GetEntriesFast();
for(Int_t ivar=0; ivar<nvars; ivar++) {
TVectorD *bins = static_cast<TVectorD*>(fAxes.At(ivar));
Int_t nbins = bins->GetNrows()-1;
Int_t ibin = (ihist/sizeAdd)%nbins;
Double_t lowEdge = (*bins)[ibin];
Double_t upEdge = (*bins)[ibin+1];
switch(fBinType[ivar]) {
case kStdBin: upEdge=(*bins)[ibin+1]; break;
case kBinToMax: upEdge=(*bins)[nbins]; break;
case kBinFromMin: lowEdge=(*bins)[0]; break;
case kSymBin: upEdge=(*bins)[nbins-ibin];
if(ibin>=((Double_t)(nbins+1))/2) upEdge=(*bins)[nbins];
break;
}
if(fVarCutType->TestBitNumber(ivar)) {
if( (valuesLeg1[fVarCuts[ivar]] < lowEdge || valuesLeg1[fVarCuts[ivar]] >= upEdge) ||
(valuesLeg2[fVarCuts[ivar]] < lowEdge || valuesLeg2[fVarCuts[ivar]] >= upEdge) ) {
selected=kFALSE;
break;
}
}
else {
if( (valuesPair[fVarCuts[ivar]] < lowEdge || valuesPair[fVarCuts[ivar]] >= upEdge) ) {
selected=kFALSE;
break;
}
}
sizeAdd*=nbins;
}
if(!selected) continue;
TObjArray *tmp = (TObjArray*) histArr->At(ihist);
TString title = tmp->GetName();
AliDebug(10,title.Data());
for(Int_t i=0; i<tmp->GetEntriesFast(); i++) {
AliDielectronHistos::FillValues(tmp->At(i), valuesPair);
}
}
}
void AliDielectronHF::Init()
{
fHasMC=AliDielectronMC::Instance()->HasMC();
Int_t steps = 0;
if(fHasMC) steps=fSignalsMC->GetEntries();
if(fStepGenerated) steps*=2;
if(fEventArray) steps=1;
fArrPairType.SetName(Form("%s_HF",GetName()));
if( (fHasMC && fPairType==kMConly) || fEventArray) fArrPairType.Expand(steps);
else fArrPairType.Expand(AliDielectron::kEv1PMRot+1+steps);
Int_t size = GetNumberOfBins();
AliDebug(10,Form("Creating a histo array with size %d \n",size));
Int_t sizeAdd = 1;
TObjArray *histArr = new TObjArray(0);
if(!histArr) return;
histArr->SetOwner(kTRUE);
histArr->Expand(size);
for(Int_t ihist=0; ihist<size; ihist++) {
histArr->AddAt(fRefObj.Clone(""), ihist);
}
Int_t nvars = fAxes.GetEntriesFast();
for(Int_t ivar=0; ivar<nvars; ivar++) {
TVectorD *bins = static_cast<TVectorD*>(fAxes.At(ivar));
Int_t nbins = bins->GetNrows()-1;
for(Int_t ihist=0; ihist<size; ihist++) {
Int_t ibin = (ihist/sizeAdd)%nbins;
Double_t lowEdge = (*bins)[ibin];
Double_t upEdge = (*bins)[ibin+1];
switch(fBinType[ivar]) {
case kStdBin: upEdge=(*bins)[ibin+1]; break;
case kBinToMax: upEdge=(*bins)[nbins]; break;
case kBinFromMin: lowEdge=(*bins)[0]; break;
case kSymBin: upEdge=(*bins)[nbins-ibin];
if(ibin>=((Double_t)(nbins+1))/2) upEdge=(*bins)[nbins];
break;
}
TObjArray *tmp= (TObjArray*) histArr->At(ihist);
TString title = tmp->GetName();
if(!ivar) title ="";
if( ivar) title+=":";
if(fVarCutType->TestBitNumber(ivar)) title+="Leg";
title+=AliDielectronVarManager::GetValueName(fVarCuts[ivar]);
title+=Form("#%.2f#%.2f",lowEdge,upEdge);
tmp->SetName(title.Data());
AliDebug(10,title.Data());
}
sizeAdd*=nbins;
}
Int_t istep=0;
if(fEventArray) {
fArrPairType[istep]=(TObjArray*)histArr->Clone("Event");
((TObjArray*)fArrPairType[istep])->SetOwner();
}
else {
if(fPairType != kMConly) {
for(istep=0; istep<AliDielectron::kEv1PMRot+1; istep++) {
if(IsPairTypeSelected(istep)) {
fArrPairType[istep]=(TObjArray*)histArr->Clone(AliDielectron::PairClassName(istep));
((TObjArray*)fArrPairType[istep])->SetOwner();
}
else {
fArrPairType[istep]=new TObjArray(0);
((TObjArray*)fArrPairType[istep])->SetOwner();
((TObjArray*)fArrPairType[istep])->SetName(AliDielectron::PairClassName(istep));
}
}
}
if(fHasMC) {
for(Int_t i=0; i<fSignalsMC->GetEntries(); i++) {
TString title = Form("(Signal: %s)",fSignalsMC->At(i)->GetTitle());
fArrPairType[istep+i]=(TObjArray*)histArr->Clone(title.Data());
if(fStepGenerated) {
title+=" MC truth";
fArrPairType[istep+i+fSignalsMC->GetEntries()]=(TObjArray*)histArr->Clone(title.Data());
}
}
}
}
if(histArr) {
delete histArr;
histArr=0;
}
}
Int_t AliDielectronHF::GetNumberOfBins() const
{
Int_t size=1;
for (Int_t i=0; i<fAxes.GetEntriesFast(); ++i)
size*=((static_cast<TVectorD*>(fAxes.At(i)))->GetNrows()-1);
return size;
}
Bool_t AliDielectronHF::IsPairTypeSelected(Int_t itype)
{
Bool_t selected = kFALSE;
if(fPairType==kAll) return kTRUE;
switch(itype) {
case AliDielectron::kEv1PP:
case AliDielectron::kEv1MM:
if(fPairType==kSeAll || fPairType==kSeMeAll || fPairType==kSeReAll ) selected = kTRUE;
break;
case AliDielectron::kEv1PM:
if(fPairType!=kMeOnlyOS) selected = kTRUE;
break;
case AliDielectron::kEv1PEv2P:
case AliDielectron::kEv1MEv2M:
if(fPairType==kMeAll || fPairType==kSeMeAll) selected = kTRUE;
break;
case AliDielectron::kEv1PEv2M:
if(fPairType==kMeAll || fPairType==kSeMeAll) selected = kTRUE;
break;
case AliDielectron::kEv1MEv2P:
if(fPairType==kMeAll || fPairType==kSeMeAll || fPairType==kMeOnlyOS || fPairType==kSeMeOnlyOS) selected = kTRUE;
break;
case AliDielectron::kEv2PP:
case AliDielectron::kEv2PM:
case AliDielectron::kEv2MM:
selected = kFALSE;
break;
case AliDielectron::kEv1PMRot:
if(fPairType==kSeReAll || fPairType==kSeReOnlyOS) selected = kTRUE;
break;
}
return selected;
}