#include <TList.h>
#include <TObjArray.h>
#include <TVectorD.h>
#include <TString.h>
#include <TObjString.h>
#include <AliCFContainer.h>
#include <AliAnalysisFilter.h>
#include <AliAnalysisCuts.h>
#include <AliVParticle.h>
#include <AliLog.h>
#include "AliDielectronCF.h"
#include "AliDielectronMC.h"
#include "AliDielectronPair.h"
#include "AliDielectronSignalMC.h"
ClassImp(AliDielectronCF)
AliDielectronCF::AliDielectronCF() :
TNamed("DielectronCF","DielectronCF"),
fUsedVars(new TBits(AliDielectronVarManager::kNMaxValues)),
fNSteps(0),
fNVars(0),
fVarBinLimits(0x0),
fNVarsLeg(0),
fVarBinLimitsLeg(0x0),
fNCuts(0),
fValues(0x0),
fIsMCTruth(0x0),
fStepForMCtruth(kFALSE),
fStepForNoCutsMCmotherPid(kFALSE),
fStepForAfterAllCuts(kTRUE),
fStepForPreFilter(kFALSE),
fStepsForEachCut(kFALSE),
fStepsForCutsIncreasing(kFALSE),
fStepsForSignal(kTRUE),
fStepsForBackground(kFALSE),
fStepsForMCtruthOnly(kFALSE),
fNStepMasks(0),
fPdgMother(-1),
fSignalsMC(0x0),
fCfContainer(0x0),
fHasMC(kFALSE),
fNAddSteps(0)
{
for (Int_t i=0; i<AliDielectronVarManager::kNMaxValues; ++i){
fVariables[i]=0;
fVariablesLeg[i]=0;
}
for (Int_t i=0; i<kNmaxAddSteps; ++i){
fStepMasks[i]=0xFFFFFF;
}
}
AliDielectronCF::AliDielectronCF(const char* name, const char* title) :
TNamed(name, title),
fUsedVars(new TBits(AliDielectronVarManager::kNMaxValues)),
fNSteps(0),
fNVars(0),
fVarBinLimits(0x0),
fNVarsLeg(0),
fVarBinLimitsLeg(0x0),
fNCuts(0),
fValues(0x0),
fIsMCTruth(0x0),
fStepForMCtruth(kFALSE),
fStepForNoCutsMCmotherPid(kFALSE),
fStepForAfterAllCuts(kTRUE),
fStepForPreFilter(kFALSE),
fStepsForEachCut(kFALSE),
fStepsForCutsIncreasing(kFALSE),
fStepsForSignal(kTRUE),
fStepsForBackground(kFALSE),
fStepsForMCtruthOnly(kFALSE),
fNStepMasks(0),
fPdgMother(-1),
fSignalsMC(0x0),
fCfContainer(0x0),
fHasMC(kFALSE),
fNAddSteps(0)
{
for (Int_t i=0; i<AliDielectronVarManager::kNMaxValues; ++i){
fVariables[i]=0;
fVariablesLeg[i]=0;
}
for (Int_t i=0; i<kNmaxAddSteps; ++i){
fStepMasks[i]=0xFFFFFF;
}
}
AliDielectronCF::~AliDielectronCF()
{
if (fUsedVars) delete fUsedVars;
if (fValues) delete [] fValues;
if (fIsMCTruth) delete [] fIsMCTruth;
if (fVarBinLimits) delete fVarBinLimits;
if (fVarBinLimitsLeg) delete fVarBinLimitsLeg;
}
void AliDielectronCF::AddVariable(AliDielectronVarManager::ValueTypes type, Int_t nbins,
Double_t min, Double_t max, Bool_t leg, Bool_t log)
{
TVectorD *binLimits=0x0;
if (!log) binLimits=MakeLinBinning(nbins,min,max);
else binLimits=MakeLogBinning(nbins,min,max);
AddVariable(type,binLimits,leg);
}
void AliDielectronCF::AddVariable(AliDielectronVarManager::ValueTypes type, const char* binLimitStr, Bool_t leg)
{
TString limits(binLimitStr);
if (limits.IsNull()){
AliError(Form("Bin Limit string is empty, cannot add the variable '%s'",AliDielectronVarManager::GetValueName(type)));
return;
}
TObjArray *arr=limits.Tokenize(",");
Int_t nLimits=arr->GetEntries();
if (nLimits<2){
AliError(Form("Need at leas 2 bin limits, cannot add the variable '%s'",AliDielectronVarManager::GetValueName(type)));
delete arr;
return;
}
TVectorD *binLimits=new TVectorD(nLimits);
for (Int_t iLim=0; iLim<nLimits; ++iLim){
(*binLimits)[iLim]=(static_cast<TObjString*>(arr->At(iLim)))->GetString().Atof();
}
delete arr;
AddVariable(type,binLimits,leg);
}
void AliDielectronCF::AddVariable(AliDielectronVarManager::ValueTypes type, TVectorD *binLimits, Bool_t leg)
{
if (!leg){
if (!fVarBinLimits){
fVarBinLimits=new TObjArray;
fVarBinLimits->SetOwner();
}
fVarBinLimits->Add(binLimits);
fVariables[fNVars] = (UInt_t)type;
++fNVars;
} else {
if (!fVarBinLimitsLeg){
fVarBinLimitsLeg=new TObjArray;
fVarBinLimitsLeg->SetOwner();
}
fVarBinLimitsLeg->Add(binLimits);
fVariablesLeg[fNVarsLeg] = (UInt_t)type;
fUsedVars->SetBitNumber(type,kTRUE);
++fNVarsLeg;
}
}
void AliDielectronCF::InitialiseContainer(const AliAnalysisFilter& filter)
{
fNCuts=filter.GetCuts()->GetEntries();
fHasMC=AliDielectronMC::Instance()->HasMC();
fNAddSteps=1;
if (fHasMC){
if (fStepsForSignal && fSignalsMC) fNAddSteps+=fSignalsMC->GetEntries();
if (fStepsForBackground) ++fNAddSteps;
if (fStepsForMCtruthOnly) --fNAddSteps;
} else {
fStepForMCtruth=kFALSE;
fStepForNoCutsMCmotherPid=kFALSE;
fStepsForSignal=kFALSE;
fStepsForBackground=kFALSE;
}
if (fStepsForCutsIncreasing) fStepForAfterAllCuts=kFALSE;
if (fStepsForEachCut&&fNCuts==1) fStepForAfterAllCuts=kFALSE;
fNSteps=0;
if (fStepForMCtruth && fSignalsMC) fNSteps+=fSignalsMC->GetEntries();
if (fStepForNoCutsMCmotherPid && fSignalsMC) fNSteps+=fSignalsMC->GetEntries();
if (fStepForAfterAllCuts) fNSteps+=fNAddSteps;
if (fStepsForEachCut) fNSteps+=(fNAddSteps*fNCuts);
if (fStepsForCutsIncreasing) fNSteps+=(fNAddSteps*fNCuts);
fNSteps+=(fNAddSteps*fNStepMasks);
if (fStepForPreFilter) fNSteps+=fNAddSteps;
Int_t *nbins=new Int_t[fNVars+2*fNVarsLeg];
for (Int_t i=0;i<fNVars;++i) {
Int_t nBins=(static_cast<TVectorD*>(fVarBinLimits->At(i)))->GetNrows()-1;
nbins[i]=nBins;
}
for (Int_t i=0;i<fNVarsLeg;++i){
Int_t nBins=(static_cast<TVectorD*>(fVarBinLimitsLeg->At(i)))->GetNrows()-1;
nbins[i+fNVars]=nBins;
nbins[i+fNVars+fNVarsLeg]=nBins;
}
fCfContainer = new AliCFContainer(GetName(), GetTitle(), fNSteps, fNVars+2*fNVarsLeg, nbins);
delete [] nbins;
for (Int_t iVar=0; iVar<fNVars; iVar++) {
UInt_t type=fVariables[iVar];
Double_t *binLim = (static_cast<TVectorD*>(fVarBinLimits->At(iVar)))->GetMatrixArray();
fCfContainer->SetBinLimits(iVar, binLim);
fCfContainer->SetVarTitle(iVar, AliDielectronVarManager::GetValueName(type));
}
for (Int_t iVar=0; iVar<fNVarsLeg; iVar++) {
UInt_t type=fVariablesLeg[iVar];
Double_t *binLim=(static_cast<TVectorD*>(fVarBinLimitsLeg->At(iVar)))->GetMatrixArray();
fCfContainer->SetBinLimits(iVar+fNVars, binLim);
fCfContainer->SetVarTitle(iVar+fNVars, Form("Leg1_%s",AliDielectronVarManager::GetValueName(type)));
fCfContainer->SetBinLimits(iVar+fNVars+fNVarsLeg, binLim);
fCfContainer->SetVarTitle(iVar+fNVars+fNVarsLeg, Form("Leg2_%s",AliDielectronVarManager::GetValueName(type)));
}
fValues = new Double_t[fNVars+2*fNVarsLeg];
if (fHasMC && fSignalsMC && fSignalsMC->GetEntries()>0) fIsMCTruth=new Bool_t[fSignalsMC->GetEntries()];
Int_t step=0;
if(fStepForMCtruth && fSignalsMC) {
for(Int_t i=0; i<fSignalsMC->GetEntries(); i++)
fCfContainer->SetStepTitle(step++, Form("MC truth (Signal: %s)", fSignalsMC->At(i)->GetTitle()));
}
if (fStepForNoCutsMCmotherPid && fSignalsMC){
for(Int_t i=0; i<fSignalsMC->GetEntries(); i++)
fCfContainer->SetStepTitle(step++,Form("No cuts (Signal: %s)",fSignalsMC->At(i)->GetTitle()));
}
TString cutName;
if (fStepsForEachCut){
for (Int_t iCut=0; iCut<fNCuts;++iCut) {
cutName=filter.GetCuts()->At(iCut)->GetName();
if (!fStepsForMCtruthOnly) {
fCfContainer->SetStepTitle(step++, cutName.Data());
}
if (fHasMC){
if (fStepsForSignal && fSignalsMC) {
for(Int_t i=0; i<fSignalsMC->GetEntries(); i++) {
fCfContainer->SetStepTitle(step++, Form("%s (Signal: %s)", cutName.Data(), fSignalsMC->At(i)->GetTitle()));
}
}
if (fStepsForBackground)
fCfContainer->SetStepTitle(step++, (cutName+" (Background)").Data());
}
}
}
if (fStepsForCutsIncreasing){
cutName="";
for (Int_t iCut=0; iCut<fNCuts;++iCut) {
if (!cutName.IsNull()) cutName+="&";
cutName+=filter.GetCuts()->At(iCut)->GetName();
if (!fStepsForMCtruthOnly) {
fCfContainer->SetStepTitle(step++, cutName.Data());
}
if (fHasMC){
if (fStepsForSignal && fSignalsMC)
for(Int_t i=0; i<fSignalsMC->GetEntries(); i++) {
fCfContainer->SetStepTitle(step++, Form("%s (Signal: %s)", cutName.Data(), fSignalsMC->At(i)->GetTitle()));
}
if (fStepsForBackground)
fCfContainer->SetStepTitle(step++, (cutName+" (Background)").Data());
}
}
}
for (UInt_t iComb=0; iComb<fNStepMasks; ++iComb){
cutName="";
UInt_t mask=fStepMasks[iComb];
for (Int_t iCut=0; iCut<fNCuts;++iCut) {
if (mask&(1<<iCut)){
if (cutName.IsNull()){
cutName=filter.GetCuts()->At(iCut)->GetName();
}else{
cutName+="&";
cutName+=filter.GetCuts()->At(iCut)->GetName();
}
}
}
fCfContainer->SetStepTitle(step++, cutName.Data());
if (fHasMC){
if (fStepsForSignal && fSignalsMC)
for(Int_t i=0; i<fSignalsMC->GetEntries(); i++) {
fCfContainer->SetStepTitle(step++, Form("%s (Signal: %s)", cutName.Data(), fSignalsMC->At(i)->GetTitle()));
}
if (fStepsForBackground)
fCfContainer->SetStepTitle(step++, (cutName+" (Background)").Data());
}
}
if (fStepForAfterAllCuts){
cutName="No pair cuts";
if (filter.GetCuts()->At(0)){
cutName=filter.GetCuts()->At(0)->GetName();
for (Int_t iCut=1; iCut<fNCuts;++iCut) {
cutName+="&";
cutName+=filter.GetCuts()->At(iCut)->GetName();
}
}
if (!fStepsForMCtruthOnly) {
fCfContainer->SetStepTitle(step++, cutName.Data());
}
if (fHasMC){
if (fStepsForSignal && fSignalsMC)
for(Int_t i=0; i<fSignalsMC->GetEntries(); i++) {
fCfContainer->SetStepTitle(step++, Form("%s (Signal: %s)", cutName.Data(), fSignalsMC->At(i)->GetTitle()));
}
if (fStepsForBackground)
fCfContainer->SetStepTitle(step++, (cutName+" (Background)").Data());
}
}
if (fStepForPreFilter){
cutName="PreFilter";
fCfContainer->SetStepTitle(step++, cutName.Data());
if (fHasMC){
if (fStepsForSignal && fSignalsMC)
for(Int_t i=0; i<fSignalsMC->GetEntries(); i++) {
fCfContainer->SetStepTitle(step++, Form("%s (Signal %s)", cutName.Data(), fSignalsMC->At(i)->GetTitle()));
}
if (fStepsForBackground)
fCfContainer->SetStepTitle(step++, (cutName+" (Background)").Data());
}
}
if (step!=fNSteps) {
AliError(Form("Something went wrong in the naming of the steps!!! (%d != %d)",step,fNSteps));
}
}
void AliDielectronCF::Fill(UInt_t mask, const AliDielectronPair *particle)
{
if(fIsMCTruth) {
for(Int_t i=0; i<fSignalsMC->GetEntries(); i++) fIsMCTruth[i]=kFALSE;
}
Bool_t isMixedPair=(particle->GetType()>2&&particle->GetType()<10);
Bool_t isBackground = kFALSE;
if(fIsMCTruth && !isMixedPair) {
for(Int_t i=0; i<fSignalsMC->GetEntries(); ++i) {
fIsMCTruth[i] = AliDielectronMC::Instance()->IsMCTruth(particle, (AliDielectronSignalMC*)fSignalsMC->At(i));
isBackground = (isBackground || fIsMCTruth[i]);
}
isBackground = !isBackground;
}
Double_t valuesPair[AliDielectronVarManager::kNMaxValues];
AliDielectronVarManager::SetFillMap(fUsedVars);
AliDielectronVarManager::Fill(particle,valuesPair);
for (Int_t iVar=0; iVar<fNVars; ++iVar){
Int_t var=fVariables[iVar];
fValues[iVar]=valuesPair[var];
}
if (fNVarsLeg>0){
Double_t valuesLeg1[AliDielectronVarManager::kNMaxValues]={0};
AliDielectronVarManager::Fill(particle->GetFirstDaughterP(),valuesLeg1);
Double_t valuesLeg2[AliDielectronVarManager::kNMaxValues]={0};
AliDielectronVarManager::Fill(particle->GetSecondDaughterP(),valuesLeg2);
for (Int_t iVar=0; iVar<fNVarsLeg; ++iVar){
Int_t var=fVariablesLeg[iVar];
fValues[iVar+fNVars]=valuesLeg1[var];
fValues[iVar+fNVars+fNVarsLeg]=valuesLeg2[var];
}
}
UInt_t selectedMask=(1<<fNCuts)-1;
Int_t step=0;
if (fStepForMCtruth && fIsMCTruth) step+=fSignalsMC->GetEntries();
if (fStepForNoCutsMCmotherPid && fIsMCTruth){
for(Int_t i=0; i<fSignalsMC->GetEntries(); ++i) {
if(fIsMCTruth[i]) {
fCfContainer->Fill(fValues,step);
}
++step;
}
}
if (fStepsForEachCut){
for (Int_t iCut=0; iCut<fNCuts;++iCut) {
UInt_t cutMask=1<<iCut;
if ((mask&cutMask)==cutMask) {
if(!fStepsForMCtruthOnly) {
fCfContainer->Fill(fValues,step);
++step;
}
if (fHasMC){
if ( fStepsForSignal && fIsMCTruth){
for(Int_t i=0; i<fSignalsMC->GetEntries(); ++i) {
if(fIsMCTruth[i]) {
fCfContainer->Fill(fValues,step);
}
++step;
}
}
if ( fStepsForBackground ){
if (isBackground) fCfContainer->Fill(fValues,step);
++step;
}
}
} else {
step+=fNAddSteps;
}
}
}
if (fStepsForCutsIncreasing&&fNCuts>2){
for (Int_t iCut=0; iCut<fNCuts;++iCut) {
UInt_t cutMask=(1<<(iCut+1))-1;
if ((mask&cutMask)==cutMask) {
if(!fStepsForMCtruthOnly) {
fCfContainer->Fill(fValues,step);
++step;
}
if (fHasMC){
if ( fStepsForSignal && fIsMCTruth){
for(Int_t i=0; i<fSignalsMC->GetEntries(); ++i) {
if(fIsMCTruth[i]) {
fCfContainer->Fill(fValues,step);
}
++step;
}
}
if ( fStepsForBackground ){
if (isBackground) fCfContainer->Fill(fValues,step);
++step;
}
}
} else {
step+=fNAddSteps;
}
}
}
for (UInt_t iComb=0; iComb<fNStepMasks; ++iComb){
UInt_t userMask=fStepMasks[iComb];
if ((mask&userMask)==userMask) {
if(!fStepsForMCtruthOnly) {
fCfContainer->Fill(fValues,step);
++step;
}
if (fHasMC){
if ( fStepsForSignal && fIsMCTruth){
for(Int_t i=0; i<fSignalsMC->GetEntries(); ++i) {
if(fIsMCTruth[i]) {
fCfContainer->Fill(fValues,step);
}
++step;
}
}
if ( fStepsForBackground ){
if (isBackground) fCfContainer->Fill(fValues,step);
++step;
}
}
} else {
step+=fNAddSteps;
}
}
if (fStepForAfterAllCuts){
if (mask == selectedMask){
if(!fStepsForMCtruthOnly) {
fCfContainer->Fill(fValues,step);
++step;
}
if (fHasMC){
if ( fStepsForSignal && fIsMCTruth){
for(Int_t i=0; i<fSignalsMC->GetEntries(); ++i) {
if(fIsMCTruth[i]) {
fCfContainer->Fill(fValues,step);
}
++step;
}
}
if ( fStepsForBackground ){
if (isBackground) fCfContainer->Fill(fValues,step);
++step;
}
}
} else {
step+=fNAddSteps;
}
}
if (fStepForPreFilter) {
if (mask&(1<<fNCuts)) {
if(!fStepsForMCtruthOnly) {
fCfContainer->Fill(fValues,step);
++step;
}
if (fHasMC){
if ( fStepsForSignal && fIsMCTruth){
for(Int_t i=0; i<fSignalsMC->GetEntries(); ++i) {
if(fIsMCTruth[i]) {
fCfContainer->Fill(fValues,step);
}
++step;
}
}
if ( fStepsForBackground ){
if (isBackground) fCfContainer->Fill(fValues,step);
++step;
}
}
}
else {
step+=fNAddSteps;
}
}
if (step!=fNSteps) {
AliError("Something went wrong in the step filling!!!");
}
}
void AliDielectronCF::FillMC(const TObject *particle)
{
if (!fStepForMCtruth) return;
Double_t valuesPair[AliDielectronVarManager::kNMaxValues];
AliDielectronVarManager::SetFillMap(fUsedVars);
AliDielectronVarManager::Fill(particle,valuesPair);
AliVParticle *d1=0x0;
AliVParticle *d2=0x0;
AliDielectronMC::Instance()->GetDaughters(particle,d1,d2);
valuesPair[AliDielectronVarManager::kPairType]=1;
for (Int_t iVar=0; iVar<fNVars; ++iVar){
Int_t var=fVariables[iVar];
fValues[iVar]=valuesPair[var];
}
if (fNVarsLeg>0){
Double_t valuesLeg1[AliDielectronVarManager::kNMaxValues];
Double_t valuesLeg2[AliDielectronVarManager::kNMaxValues];
if (d1->Pt()>d2->Pt()){
AliDielectronVarManager::Fill(d1,valuesLeg1);
AliDielectronVarManager::Fill(d2,valuesLeg2);
} else {
AliDielectronVarManager::Fill(d2,valuesLeg1);
AliDielectronVarManager::Fill(d1,valuesLeg2);
}
for (Int_t iVar=0; iVar<fNVarsLeg; ++iVar){
Int_t var=fVariablesLeg[iVar];
fValues[iVar+fNVars]=valuesLeg1[var];
fValues[iVar+fNVars+fNVarsLeg]=valuesLeg2[var];
}
}
fCfContainer->Fill(fValues,0);
}
void AliDielectronCF::FillMC(Int_t label1, Int_t label2, Int_t nSignal) {
if (!fStepForMCtruth) 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);
if (fNVarsLeg>0){
Double_t valuesLeg1[AliDielectronVarManager::kNMaxValues];
Double_t valuesLeg2[AliDielectronVarManager::kNMaxValues];
if (part1->Pt()>part2->Pt()){
AliDielectronVarManager::Fill(part1,valuesLeg1);
AliDielectronVarManager::Fill(part2,valuesLeg2);
} else {
AliDielectronVarManager::Fill(part2,valuesLeg1);
AliDielectronVarManager::Fill(part1,valuesLeg2);
}
for (Int_t iVar=0; iVar<fNVarsLeg; ++iVar){
Int_t var=fVariablesLeg[iVar];
fValues[iVar+fNVars]=valuesLeg1[var];
fValues[iVar+fNVars+fNVarsLeg]=valuesLeg2[var];
}
}
Double_t valuesPair[AliDielectronVarManager::kNMaxValues];
AliDielectronVarManager::Fill(dieMC->GetMCEvent(), valuesPair);
AliDielectronVarManager::FillVarMCParticle2(part1,part2,valuesPair);
if(part1->Charge()*part2->Charge()<0)
valuesPair[AliDielectronVarManager::kPairType]=1;
else if(part1->Charge()>0)
valuesPair[AliDielectronVarManager::kPairType]=0;
else
valuesPair[AliDielectronVarManager::kPairType]=2;
for(Int_t iVar=0; iVar<fNVars; ++iVar){
Int_t var=fVariables[iVar];
fValues[iVar]=valuesPair[var];
}
fCfContainer->Fill(fValues,nSignal);
}
TVectorD* AliDielectronCF::MakeLogBinning(Int_t nbinsX, Double_t xmin, Double_t xmax) const
{
if (xmin<1e-20 || xmax<1e-20){
AliError("For Log binning xmin and xmax must be > 1e-20. Using linear binning instead!");
return MakeLinBinning(nbinsX, xmin, xmax);
}
if (xmax<xmin){
Double_t tmp=xmin;
xmin=xmax;
xmax=tmp;
}
TVectorD *binLim=new TVectorD(nbinsX+1);
Double_t first=xmin;
Double_t last=xmax;
Double_t expMax=TMath::Log(last/first);
for (Int_t i=0; i<nbinsX+1; ++i){
(*binLim)[i]=first*TMath::Exp(expMax/nbinsX*(Double_t)i);
}
return binLim;
}
TVectorD* AliDielectronCF::MakeLinBinning(Int_t nbinsX, Double_t xmin, Double_t xmax) const
{
if (xmax<xmin){
Double_t tmp=xmin;
xmin=xmax;
xmax=tmp;
}
TVectorD *binLim=new TVectorD(nbinsX+1);
Double_t first=xmin;
Double_t last=xmax;
Double_t binWidth=(last-first)/nbinsX;
for (Int_t i=0; i<nbinsX+1; ++i){
(*binLim)[i]=first+binWidth*(Double_t)i;
}
return binLim;
}