[u/mrichter/AliRoot.git] / PWG3 / dielectron / AliDielectronCF.cxx

/*************************************************************************
* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
*                                                                        *
* Author: The ALICE Off-line Project.                                    *
* Contributors are mentioned in the code where appropriate.              *
*                                                                        *
* Permission to use, copy, modify and distribute this software and its   *
* documentation strictly for non-commercial purposes is hereby granted   *
* without fee, provided that the above copyright notice appears in all   *
* copies and that both the copyright notice and this permission notice   *
* appear in the supporting documentation. The authors make no claims     *
* about the suitability of this software for any purpose. It is          *
* provided "as is" without express or implied warranty.                  *
**************************************************************************/

///////////////////////////////////////////////////////////////////////////
//       Dielectron Correction framework manager                         //
//                                                                       //
/*


*/
//                                                                       //
///////////////////////////////////////////////////////////////////////////

#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"

ClassImp(AliDielectronCF)

AliDielectronCF::AliDielectronCF() :
  TNamed("DielectronCF","DielectronCF"),
  fNSteps(0),
  fNVars(0),
  fVarBinLimits(0x0),
  fNVarsLeg(0),
  fVarBinLimitsLeg(0x0),
  fNCuts(0),
  fValues(0x0),
  fStepForMCtruth(kFALSE),
  fStepForNoCutsMCmotherPid(kFALSE),
  fStepForAfterAllCuts(kTRUE),
  fStepsForEachCut(kFALSE),
  fStepsForCutsIncreasing(kFALSE),
  fStepsForSignal(kTRUE),
  fStepsForBackground(kFALSE),
  fNStepMasks(0),
  fPdgMother(-1),
  fCfContainer(0x0),
  fHasMC(kFALSE),
  fNAddSteps(0)
{
  //
  // Default constructor
  //
  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),
  fNSteps(0),
  fNVars(0),
  fVarBinLimits(0x0),
  fNVarsLeg(0),
  fVarBinLimitsLeg(0x0),
  fNCuts(0),
  fValues(0x0),
  fStepForMCtruth(kFALSE),
  fStepForNoCutsMCmotherPid(kFALSE),
  fStepForAfterAllCuts(kTRUE),
  fStepsForEachCut(kFALSE),
  fStepsForCutsIncreasing(kFALSE),
  fStepsForSignal(kTRUE),
  fStepsForBackground(kFALSE),
  fNStepMasks(0),
  fPdgMother(-1),
  fCfContainer(0x0),
  fHasMC(kFALSE),
  fNAddSteps(0)
{
  //
  // Named constructor
  //
  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()
{
  //
  // Destructor
  //
  if (fValues) delete [] fValues;
  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)
{
  //
  // Add a variable to the CF configuration
  // if leg is true it will add the variables of the leg
  // if log is true log binning will be created
  //

  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/*=kFALSE*/)
{
  //
  // Add a variable to the CF configuration
  // specify arbitrary binning in a string.
  // Bin limits need to be separated by a ","
  //
  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/*=kFALSE*/)
{
  //
  // Add variable with the binning given in the TVectorD
  //
  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;
    ++fNVarsLeg;
  }
}

//________________________________________________________________
void AliDielectronCF::InitialiseContainer(const AliAnalysisFilter& filter)
{
  //
  // Initialise container based on the cuts in the analysis filter
  //

  fNCuts=filter.GetCuts()->GetEntries();

  fHasMC=AliDielectronMC::Instance()->HasMC();
  fNAddSteps=1;
  if (fHasMC){
    if (fStepsForSignal) ++fNAddSteps;
    if (fStepsForBackground) ++fNAddSteps;
  } else {
    //if 
    fStepForMCtruth=kFALSE;
    fStepForNoCutsMCmotherPid=kFALSE;
    fStepsForSignal=kFALSE;
    fStepsForBackground=kFALSE;
  }
    
  fNSteps=0;
  if (fStepForMCtruth) ++fNSteps;
  if (fStepForNoCutsMCmotherPid) ++fNSteps;
  if (fStepForAfterAllCuts) fNSteps+=fNAddSteps;

  if (fStepsForEachCut&&fNCuts>1)        fNSteps+=(fNAddSteps*fNCuts);     //one step for each cut + Signal (MC)
  if (fStepsForCutsIncreasing&&fNCuts>2) fNSteps+=(fNAddSteps*(fNCuts-2)); //one step for the increasing cuts + Signal (MC)
                                                      // e.g. cut2&cut3, cut2&cut3&cut4
  fNSteps+=(fNAddSteps*fNStepMasks);                            // cuts for the additional cut masks
  // create the container
  
  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;
  
  // initialize the variables and their bin limits
  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));
  }
  
  // initialize the variables and their bin limits for the Legs
  for (Int_t iVar=0; iVar<fNVarsLeg; iVar++) {
    UInt_t type=fVariablesLeg[iVar];
    Double_t *binLim=(static_cast<TVectorD*>(fVarBinLimitsLeg->At(iVar)))->GetMatrixArray();

    //Leg1
    fCfContainer->SetBinLimits(iVar+fNVars, binLim);
    fCfContainer->SetVarTitle(iVar+fNVars, Form("Leg1_%s",AliDielectronVarManager::GetValueName(type)));
    
    //Leg2
    fCfContainer->SetBinLimits(iVar+fNVars+fNVarsLeg, binLim);
    fCfContainer->SetVarTitle(iVar+fNVars+fNVarsLeg, Form("Leg2_%s",AliDielectronVarManager::GetValueName(type)));
  }

  // array for storing values
  fValues = new Double_t[fNVars+2*fNVarsLeg];
  
  //=================//
  // Set step titles //
  //=================//
  Int_t step=0;

  //Pure MC truth
  if (fStepForMCtruth){
    fCfContainer->SetStepTitle(step++,"MC truth");
  }

  //before cuts (MC truth)
  if (fStepForNoCutsMCmotherPid){
    fCfContainer->SetStepTitle(step++,"No cuts (Signal)");
  }
  
  TString cutName;
  //Steps for each of the cuts
  if (fStepsForEachCut&&fNCuts>1){
    for (Int_t iCut=0; iCut<fNCuts;++iCut) {
      cutName=filter.GetCuts()->At(iCut)->GetName(); //TODO: User GetTitle???
      
      fCfContainer->SetStepTitle(step++, cutName.Data()); //Step for the cut
      
      if (fHasMC){
        if (fStepsForSignal)
          fCfContainer->SetStepTitle(step++, (cutName+" (Signal)").Data()); //Step for the cut with MC truth
        if (fStepsForBackground)
          fCfContainer->SetStepTitle(step++, (cutName+" (Background)").Data()); //Step for the cut with MC truth
      }
    }
  }

  //Steps for increasing cut match
  if (fStepsForCutsIncreasing&&fNCuts>2){
    cutName=filter.GetCuts()->At(0)->GetName(); //TODO: User GetTitle???
    for (Int_t iCut=1; iCut<fNCuts-1;++iCut) {
      cutName+="&";
      cutName+=filter.GetCuts()->At(iCut)->GetName();
      
      fCfContainer->SetStepTitle(step++, cutName.Data()); //Step for the cut
      
      if (fHasMC){
        if (fStepsForSignal)
          fCfContainer->SetStepTitle(step++, (cutName+" (Signal)").Data()); //Step for the cut with MC truth
        if (fStepsForBackground)
          fCfContainer->SetStepTitle(step++, (cutName+" (Background)").Data()); //Step for the cut with MC truth
      }
    }
  }

  //Steps of user defined cut combinations
  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()); //Step for the cut
    
    if (fHasMC){
      if (fStepsForSignal)
        fCfContainer->SetStepTitle(step++, (cutName+" (Signal)").Data()); //Step for the cut with MC truth
      if (fStepsForBackground)
        fCfContainer->SetStepTitle(step++, (cutName+" (Background)").Data()); //Step for the cut with MC truth
    }
  }

  //After All cuts
  if (fStepForAfterAllCuts){
    cutName="No pair cuts";
    if (filter.GetCuts()->At(0)){
      cutName=filter.GetCuts()->At(0)->GetName(); //TODO: User GetTitle???
      for (Int_t iCut=1; iCut<fNCuts;++iCut) {
        cutName+="&";
        cutName+=filter.GetCuts()->At(iCut)->GetName();
      }
    }
    fCfContainer->SetStepTitle(step++, cutName.Data()); //Step for the cut
    if (fHasMC){
      if (fStepsForSignal)
        fCfContainer->SetStepTitle(step++, (cutName+" (Signal)").Data()); //Step for the cut with MC truth
      if (fStepsForBackground)
        fCfContainer->SetStepTitle(step++, (cutName+" (Background)").Data()); //Step for the cut with MC truth
    }
  }

  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)
{
  //
  // Fill the containers
  //

  Bool_t isMCTruth=kFALSE;
  if (fHasMC) isMCTruth=AliDielectronMC::Instance()->IsMotherPdg(particle,fPdgMother);

  Double_t valuesPair[AliDielectronVarManager::kNMaxValues];
  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];
    AliDielectronVarManager::Fill(particle->GetFirstDaughter(),valuesLeg1);
    Double_t valuesLeg2[AliDielectronVarManager::kNMaxValues];
    AliDielectronVarManager::Fill(particle->GetSecondDaughter(),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;

  //============//
  // Fill steps //
  //============//
  // step 0 would be full MC truth and is handled in FillMC
  Int_t step=0;
  if (fStepForMCtruth) ++step;
  
  //No cuts (MC truth)
  if (fStepForNoCutsMCmotherPid){
    if (isMCTruth) fCfContainer->Fill(fValues,step);
    ++step;
  }
  
  //Steps for each of the cuts
  if (fStepsForEachCut&&fNCuts>1){
    for (Int_t iCut=0; iCut<fNCuts;++iCut) {
      if (mask&(1<<iCut)) {
        fCfContainer->Fill(fValues,step);
        ++step;

        if (fHasMC){
          if ( fStepsForSignal){
            if (isMCTruth) fCfContainer->Fill(fValues,step);
            ++step;
          }
          if ( fStepsForBackground ){
            if (!isMCTruth) fCfContainer->Fill(fValues,step);
            ++step;
          }
        }
      } else {
        step+=fNAddSteps;
      }
    }
  }

  //Steps for increasing cut match
  if (fStepsForCutsIncreasing&&fNCuts>2){
    for (Int_t iCut=1; iCut<fNCuts-1;++iCut) {
      if (mask&(1<<((iCut+1)-1))) {
        fCfContainer->Fill(fValues,step);
        ++step;
        
        if (fHasMC){
          if ( fStepsForSignal){
            if (isMCTruth) fCfContainer->Fill(fValues,step);
            ++step;
          }
          if ( fStepsForBackground ){
            if (!isMCTruth) fCfContainer->Fill(fValues,step);
            ++step;
          }
        }
      } else {
        step+=fNAddSteps;
      }
    }
  }

  //Steps of user defined cut combinations
  for (UInt_t iComb=0; iComb<fNStepMasks; ++iComb){
    UInt_t userMask=fStepMasks[iComb];
    if (mask&userMask) {
      fCfContainer->Fill(fValues,step);
      ++step;
      
      if (fHasMC){
        if ( fStepsForSignal){
          if (isMCTruth) fCfContainer->Fill(fValues,step);
          ++step;
        }
        if ( fStepsForBackground ){
          if (!isMCTruth) fCfContainer->Fill(fValues,step);
          ++step;
        }
      }
    } else {
      step+=fNAddSteps;
    }
  }
  
  //All cuts
  if (fStepForAfterAllCuts){
    if (mask == selectedMask){
      fCfContainer->Fill(fValues,step);
      ++step;
      
      if (fHasMC){
        if ( fStepsForSignal){
          if (isMCTruth) fCfContainer->Fill(fValues,step);
          ++step;
        }
        if ( fStepsForBackground ){
          if (!isMCTruth) 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)
{
  //
  // fill MC part of the Container
  //
  if (!fStepForMCtruth) return;
  
  Double_t valuesPair[AliDielectronVarManager::kNMaxValues];
  AliDielectronVarManager::Fill(particle,valuesPair);

  AliVParticle *d1=0x0;
  AliVParticle *d2=0x0;
  AliDielectronMC::Instance()->GetDaughters(particle,d1,d2);
  
  valuesPair[AliDielectronVarManager::kThetaHE]=AliDielectronPair::ThetaPhiCM(d1,d2,kTRUE,kTRUE);
  valuesPair[AliDielectronVarManager::kPhiHE]=AliDielectronPair::ThetaPhiCM(d1,d2,kTRUE,kFALSE);
  valuesPair[AliDielectronVarManager::kThetaCS]=AliDielectronPair::ThetaPhiCM(d1,d2,kFALSE,kTRUE);
  valuesPair[AliDielectronVarManager::kPhiCS]=AliDielectronPair::ThetaPhiCM(d1,d2,kFALSE,kFALSE);
  
  //TODO: temporary solution, set manually the pair type to 1: unlikesign SE
  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);
}

//_____________________________________________________________________________
TVectorD* AliDielectronCF::MakeLogBinning(Int_t nbinsX, Double_t xmin, Double_t xmax) const
{
  //
  // Make logarithmic binning
  // the user has to delete the array afterwards!!!
  //
  
  //check limits
  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
{
  //
  // Make logarithmic binning
  // the user has to delete the array afterwards!!!
  //
  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;
}
Commit	Line	Data
b2a297fa	1	/*************************************************************************
	2	* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	///////////////////////////////////////////////////////////////////////////
	17	// Dielectron Correction framework manager //
	18	// //
	19	/*
	20
	21
	22
	23
	24
	25
	26
	27
	28
	29	*/
	30	// //
	31	///////////////////////////////////////////////////////////////////////////
	32
	33	#include <TList.h>
61d106d3	34	#include <TObjArray.h>
61d106d3	35	#include <TVectorD.h>
3505bfad	36	#include <TString.h>
3505bfad	37	#include <TObjString.h>
b2a297fa	38
	39	#include <AliCFContainer.h>
	40	#include <AliAnalysisFilter.h>
	41	#include <AliAnalysisCuts.h>
6551594b	42	#include <AliVParticle.h>
b2a297fa	43	#include <AliLog.h>
	44
	45	#include "AliDielectronCF.h"
	46	#include "AliDielectronMC.h"
6551594b	47	#include "AliDielectronPair.h"
b2a297fa	48
	49	ClassImp(AliDielectronCF)
	50
	51	AliDielectronCF::AliDielectronCF() :
	52	TNamed("DielectronCF","DielectronCF"),
	53	fNSteps(0),
	54	fNVars(0),
61d106d3	55	fVarBinLimits(0x0),
6551594b	56	fNVarsLeg(0),
3505bfad	57	fVarBinLimitsLeg(0x0),
b2a297fa	58	fNCuts(0),
6551594b	59	fValues(0x0),
e123f993	60	fStepForMCtruth(kFALSE),
	61	fStepForNoCutsMCmotherPid(kFALSE),
	62	fStepForAfterAllCuts(kTRUE),
	63	fStepsForEachCut(kFALSE),
b2a297fa	64	fStepsForCutsIncreasing(kFALSE),
572b0139	65	fStepsForSignal(kTRUE),
572b0139	66	fStepsForBackground(kFALSE),
b2a297fa	67	fNStepMasks(0),
b2a297fa	68	fPdgMother(-1),
572b0139	69	fCfContainer(0x0),
	70	fHasMC(kFALSE),
	71	fNAddSteps(0)
b2a297fa	72	{
	73	//
	74	// Default constructor
	75	//
	76	for (Int_t i=0; i<AliDielectronVarManager::kNMaxValues; ++i){
	77	fVariables[i]=0;
61d106d3	78	fVariablesLeg[i]=0;
b2a297fa	79	}
	80
	81	for (Int_t i=0; i<kNmaxAddSteps; ++i){
e123f993	82	fStepMasks[i]=0xFFFFFF;
b2a297fa	83	}
	84	}
	85
	86	//________________________________________________________________
	87	AliDielectronCF::AliDielectronCF(const char* name, const char* title) :
	88	TNamed(name, title),
	89	fNSteps(0),
	90	fNVars(0),
61d106d3	91	fVarBinLimits(0x0),
6551594b	92	fNVarsLeg(0),
3505bfad	93	fVarBinLimitsLeg(0x0),
b2a297fa	94	fNCuts(0),
6551594b	95	fValues(0x0),
e123f993	96	fStepForMCtruth(kFALSE),
	97	fStepForNoCutsMCmotherPid(kFALSE),
	98	fStepForAfterAllCuts(kTRUE),
	99	fStepsForEachCut(kFALSE),
b2a297fa	100	fStepsForCutsIncreasing(kFALSE),
572b0139	101	fStepsForSignal(kTRUE),
572b0139	102	fStepsForBackground(kFALSE),
b2a297fa	103	fNStepMasks(0),
b2a297fa	104	fPdgMother(-1),
572b0139	105	fCfContainer(0x0),
	106	fHasMC(kFALSE),
	107	fNAddSteps(0)
b2a297fa	108	{
	109	//
	110	// Named constructor
	111	//
	112	for (Int_t i=0; i<AliDielectronVarManager::kNMaxValues; ++i){
	113	fVariables[i]=0;
61d106d3	114	fVariablesLeg[i]=0;
b2a297fa	115	}
	116
	117	for (Int_t i=0; i<kNmaxAddSteps; ++i){
e123f993	118	fStepMasks[i]=0xFFFFFF;
b2a297fa	119	}
	120	}
	121
	122	//________________________________________________________________
	123	AliDielectronCF::~AliDielectronCF()
	124	{
	125	//
	126	// Destructor
	127	//
6551594b	128	if (fValues) delete [] fValues;
3505bfad	129	if (fVarBinLimits) delete fVarBinLimits;
3505bfad	130	if (fVarBinLimitsLeg) delete fVarBinLimitsLeg;
b2a297fa	131	}
	132
	133	//________________________________________________________________
3505bfad	134	void AliDielectronCF::AddVariable(AliDielectronVarManager::ValueTypes type, Int_t nbins,
3505bfad	135	Double_t min, Double_t max, Bool_t leg, Bool_t log)
b2a297fa	136	{
	137	//
	138	// Add a variable to the CF configuration
3505bfad	139	// if leg is true it will add the variables of the leg
3505bfad	140	// if log is true log binning will be created
b2a297fa	141	//
6551594b	142
3505bfad	143	TVectorD *binLimits=0x0;
	144	if (!log) binLimits=MakeLinBinning(nbins,min,max);
	145	else binLimits=MakeLogBinning(nbins,min,max);
	146	AddVariable(type,binLimits,leg);
	147	}
	148
	149	//________________________________________________________________
	150	void AliDielectronCF::AddVariable(AliDielectronVarManager::ValueTypes type, const char* binLimitStr, Bool_t leg/=kFALSE/)
	151	{
	152	//
	153	// Add a variable to the CF configuration
	154	// specify arbitrary binning in a string.
	155	// Bin limits need to be separated by a ","
	156	//
	157	TString limits(binLimitStr);
	158	if (limits.IsNull()){
	159	AliError(Form("Bin Limit string is empty, cannot add the variable '%s'",AliDielectronVarManager::GetValueName(type)));
	160	return;
	161	}
	162
	163	TObjArray *arr=limits.Tokenize(",");
	164	Int_t nLimits=arr->GetEntries();
	165	if (nLimits<2){
	166	AliError(Form("Need at leas 2 bin limits, cannot add the variable '%s'",AliDielectronVarManager::GetValueName(type)));
	167	delete arr;
	168	return;
	169	}
	170
	171	TVectorD *binLimits=new TVectorD(nLimits);
	172	for (Int_t iLim=0; iLim<nLimits; ++iLim){
	173	(binLimits)[iLim]=(static_cast<TObjString>(arr->At(iLim)))->GetString().Atof();
	174	}
	175
	176	delete arr;
	177	AddVariable(type,binLimits,leg);
	178	}
	179
	180	//________________________________________________________________
	181	void AliDielectronCF::AddVariable(AliDielectronVarManager::ValueTypes type, TVectorD binLimits, Bool_t leg/=kFALSE*/)
	182	{
	183	//
	184	// Add variable with the binning given in the TVectorD
	185	//
6551594b	186	if (!leg){
3505bfad	187	if (!fVarBinLimits){
	188	fVarBinLimits=new TObjArray;
	189	fVarBinLimits->SetOwner();
61d106d3	190	}
3505bfad	191	fVarBinLimits->Add(binLimits);
6551594b	192	fVariables[fNVars] = (UInt_t)type;
6551594b	193	++fNVars;
6551594b	194	} else {
3505bfad	195	if (!fVarBinLimitsLeg){
	196	fVarBinLimitsLeg=new TObjArray;
	197	fVarBinLimitsLeg->SetOwner();
	198	}
	199	fVarBinLimitsLeg->Add(binLimits);
6551594b	200	fVariablesLeg[fNVarsLeg] = (UInt_t)type;
6551594b	201	++fNVarsLeg;
6551594b	202	}
b2a297fa	203	}
	204
	205	//________________________________________________________________
	206	void AliDielectronCF::InitialiseContainer(const AliAnalysisFilter& filter)
	207	{
	208	//
	209	// Initialise container based on the cuts in the analysis filter
	210	//
	211
	212	fNCuts=filter.GetCuts()->GetEntries();
e123f993	213
572b0139	214	fHasMC=AliDielectronMC::Instance()->HasMC();
	215	fNAddSteps=1;
	216	if (fHasMC){
	217	if (fStepsForSignal) ++fNAddSteps;
	218	if (fStepsForBackground) ++fNAddSteps;
	219	} else {
	220	//if
	221	fStepForMCtruth=kFALSE;
	222	fStepForNoCutsMCmotherPid=kFALSE;
	223	fStepsForSignal=kFALSE;
	224	fStepsForBackground=kFALSE;
	225	}
	226
e123f993	227	fNSteps=0;
	228	if (fStepForMCtruth) ++fNSteps;
	229	if (fStepForNoCutsMCmotherPid) ++fNSteps;
572b0139	230	if (fStepForAfterAllCuts) fNSteps+=fNAddSteps;
	231
	232	if (fStepsForEachCut&&fNCuts>1) fNSteps+=(fNAddSteps*fNCuts); //one step for each cut + Signal (MC)
	233	if (fStepsForCutsIncreasing&&fNCuts>2) fNSteps+=(fNAddSteps*(fNCuts-2)); //one step for the increasing cuts + Signal (MC)
b2a297fa	234	// e.g. cut2&cut3, cut2&cut3&cut4
572b0139	235	fNSteps+=(fNAddSteps*fNStepMasks); // cuts for the additional cut masks
b2a297fa	236	// create the container
3505bfad	237
6551594b	238	Int_t nbins=new Int_t[fNVars+2fNVarsLeg];
3505bfad	239	for (Int_t i=0;i<fNVars;++i) {
	240	Int_t nBins=(static_cast<TVectorD*>(fVarBinLimits->At(i)))->GetNrows()-1;
	241	nbins[i]=nBins;
	242	}
	243	for (Int_t i=0;i<fNVarsLeg;++i){
	244	Int_t nBins=(static_cast<TVectorD*>(fVarBinLimitsLeg->At(i)))->GetNrows()-1;
	245	nbins[i+fNVars]=nBins;
	246	nbins[i+fNVars+fNVarsLeg]=nBins;
	247	}
6551594b	248
	249	fCfContainer = new AliCFContainer(GetName(), GetTitle(), fNSteps, fNVars+2*fNVarsLeg, nbins);
	250	delete [] nbins;
b2a297fa	251
	252	// initialize the variables and their bin limits
	253	for (Int_t iVar=0; iVar<fNVars; iVar++) {
	254	UInt_t type=fVariables[iVar];
3505bfad	255	Double_t binLim = (static_cast<TVectorD>(fVarBinLimits->At(iVar)))->GetMatrixArray();
3505bfad	256
b2a297fa	257	fCfContainer->SetBinLimits(iVar, binLim);
b2a297fa	258	fCfContainer->SetVarTitle(iVar, AliDielectronVarManager::GetValueName(type));
6551594b	259	}
	260
	261	// initialize the variables and their bin limits for the Legs
	262	for (Int_t iVar=0; iVar<fNVarsLeg; iVar++) {
	263	UInt_t type=fVariablesLeg[iVar];
3505bfad	264	Double_t binLim=(static_cast<TVectorD>(fVarBinLimitsLeg->At(iVar)))->GetMatrixArray();
6551594b	265
	266	//Leg1
	267	fCfContainer->SetBinLimits(iVar+fNVars, binLim);
	268	fCfContainer->SetVarTitle(iVar+fNVars, Form("Leg1_%s",AliDielectronVarManager::GetValueName(type)));
	269
	270	//Leg2
	271	fCfContainer->SetBinLimits(iVar+fNVars+fNVarsLeg, binLim);
	272	fCfContainer->SetVarTitle(iVar+fNVars+fNVarsLeg, Form("Leg2_%s",AliDielectronVarManager::GetValueName(type)));
b2a297fa	273	}
6551594b	274
	275	// array for storing values
	276	fValues = new Double_t[fNVars+2*fNVarsLeg];
b2a297fa	277
	278	//=================//
	279	// Set step titles //
	280	//=================//
	281	Int_t step=0;
	282
	283	//Pure MC truth
e123f993	284	if (fStepForMCtruth){
8df8e382	285	fCfContainer->SetStepTitle(step++,"MC truth");
e123f993	286	}
b2a297fa	287
b2a297fa	288	//before cuts (MC truth)
e123f993	289	if (fStepForNoCutsMCmotherPid){
a655b716	290	fCfContainer->SetStepTitle(step++,"No cuts (Signal)");
e123f993	291	}
b2a297fa	292
e123f993	293	TString cutName;
b2a297fa	294	//Steps for each of the cuts
e123f993	295	if (fStepsForEachCut&&fNCuts>1){
b2a297fa	296	for (Int_t iCut=0; iCut<fNCuts;++iCut) {
37e9382d	297	cutName=filter.GetCuts()->At(iCut)->GetName(); //TODO: User GetTitle???
572b0139	298
b2a297fa	299	fCfContainer->SetStepTitle(step++, cutName.Data()); //Step for the cut
572b0139	300
	301	if (fHasMC){
	302	if (fStepsForSignal)
	303	fCfContainer->SetStepTitle(step++, (cutName+" (Signal)").Data()); //Step for the cut with MC truth
	304	if (fStepsForBackground)
	305	fCfContainer->SetStepTitle(step++, (cutName+" (Background)").Data()); //Step for the cut with MC truth
	306	}
b2a297fa	307	}
	308	}
	309
	310	//Steps for increasing cut match
78091935	311	if (fStepsForCutsIncreasing&&fNCuts>2){
37e9382d	312	cutName=filter.GetCuts()->At(0)->GetName(); //TODO: User GetTitle???
b2a297fa	313	for (Int_t iCut=1; iCut<fNCuts-1;++iCut) {
	314	cutName+="&";
	315	cutName+=filter.GetCuts()->At(iCut)->GetName();
572b0139	316
b2a297fa	317	fCfContainer->SetStepTitle(step++, cutName.Data()); //Step for the cut
572b0139	318
	319	if (fHasMC){
	320	if (fStepsForSignal)
	321	fCfContainer->SetStepTitle(step++, (cutName+" (Signal)").Data()); //Step for the cut with MC truth
	322	if (fStepsForBackground)
	323	fCfContainer->SetStepTitle(step++, (cutName+" (Background)").Data()); //Step for the cut with MC truth
	324	}
b2a297fa	325	}
	326	}
	327
	328	//Steps of user defined cut combinations
	329	for (UInt_t iComb=0; iComb<fNStepMasks; ++iComb){
37e9382d	330	cutName="";
b2a297fa	331	UInt_t mask=fStepMasks[iComb];
	332	for (Int_t iCut=0; iCut<fNCuts;++iCut) {
	333	if (mask&(1<<iCut)){
	334	if (cutName.IsNull()){
	335	cutName=filter.GetCuts()->At(iCut)->GetName();
	336	}else{
	337	cutName+="&";
	338	cutName+=filter.GetCuts()->At(iCut)->GetName();
	339	}
	340	}
	341	}
572b0139	342
b2a297fa	343	fCfContainer->SetStepTitle(step++, cutName.Data()); //Step for the cut
572b0139	344
	345	if (fHasMC){
	346	if (fStepsForSignal)
	347	fCfContainer->SetStepTitle(step++, (cutName+" (Signal)").Data()); //Step for the cut with MC truth
	348	if (fStepsForBackground)
	349	fCfContainer->SetStepTitle(step++, (cutName+" (Background)").Data()); //Step for the cut with MC truth
	350	}
	351	}
	352
	353	//After All cuts
	354	if (fStepForAfterAllCuts){
	355	cutName="No pair cuts";
	356	if (filter.GetCuts()->At(0)){
	357	cutName=filter.GetCuts()->At(0)->GetName(); //TODO: User GetTitle???
	358	for (Int_t iCut=1; iCut<fNCuts;++iCut) {
	359	cutName+="&";
	360	cutName+=filter.GetCuts()->At(iCut)->GetName();
	361	}
572b0139	362	}
8df8e382	363	fCfContainer->SetStepTitle(step++, cutName.Data()); //Step for the cut
572b0139	364	if (fHasMC){
	365	if (fStepsForSignal)
	366	fCfContainer->SetStepTitle(step++, (cutName+" (Signal)").Data()); //Step for the cut with MC truth
	367	if (fStepsForBackground)
	368	fCfContainer->SetStepTitle(step++, (cutName+" (Background)").Data()); //Step for the cut with MC truth
	369	}
b2a297fa	370	}
	371
	372	if (step!=fNSteps) {
8df8e382	373	AliError(Form("Something went wrong in the naming of the steps!!! (%d != %d)",step,fNSteps));
b2a297fa	374	}
	375	}
	376
	377	//________________________________________________________________
6551594b	378	void AliDielectronCF::Fill(UInt_t mask, const AliDielectronPair *particle)
b2a297fa	379	{
	380	//
	381	// Fill the containers
	382	//
	383
	384	Bool_t isMCTruth=kFALSE;
572b0139	385	if (fHasMC) isMCTruth=AliDielectronMC::Instance()->IsMotherPdg(particle,fPdgMother);
b2a297fa	386
6551594b	387	Double_t valuesPair[AliDielectronVarManager::kNMaxValues];
6551594b	388	AliDielectronVarManager::Fill(particle,valuesPair);
b2a297fa	389
b2a297fa	390	for (Int_t iVar=0; iVar<fNVars; ++iVar){
b2a297fa	391	Int_t var=fVariables[iVar];
6551594b	392	fValues[iVar]=valuesPair[var];
b2a297fa	393	}
b2a297fa	394
6551594b	395	if (fNVarsLeg>0){
	396	Double_t valuesLeg1[AliDielectronVarManager::kNMaxValues];
	397	AliDielectronVarManager::Fill(particle->GetFirstDaughter(),valuesLeg1);
	398	Double_t valuesLeg2[AliDielectronVarManager::kNMaxValues];
	399	AliDielectronVarManager::Fill(particle->GetSecondDaughter(),valuesLeg2);
	400
	401	for (Int_t iVar=0; iVar<fNVarsLeg; ++iVar){
	402	Int_t var=fVariablesLeg[iVar];
	403	fValues[iVar+fNVars]=valuesLeg1[var];
	404	fValues[iVar+fNVars+fNVarsLeg]=valuesLeg2[var];
	405	}
	406	}
	407
b2a297fa	408	UInt_t selectedMask=(1<<fNCuts)-1;
	409
	410	//============//
	411	// Fill steps //
	412	//============//
	413	// step 0 would be full MC truth and is handled in FillMC
e123f993	414	Int_t step=0;
e123f993	415	if (fStepForMCtruth) ++step;
8df8e382	416
b2a297fa	417	//No cuts (MC truth)
e123f993	418	if (fStepForNoCutsMCmotherPid){
6551594b	419	if (isMCTruth) fCfContainer->Fill(fValues,step);
b2a297fa	420	++step;
b2a297fa	421	}
e123f993	422
b2a297fa	423	//Steps for each of the cuts
e123f993	424	if (fStepsForEachCut&&fNCuts>1){
b2a297fa	425	for (Int_t iCut=0; iCut<fNCuts;++iCut) {
b2a297fa	426	if (mask&(1<<iCut)) {
6551594b	427	fCfContainer->Fill(fValues,step);
b2a297fa	428	++step;
572b0139	429
	430	if (fHasMC){
	431	if ( fStepsForSignal){
	432	if (isMCTruth) fCfContainer->Fill(fValues,step);
	433	++step;
	434	}
	435	if ( fStepsForBackground ){
	436	if (!isMCTruth) fCfContainer->Fill(fValues,step);
	437	++step;
	438	}
	439	}
b2a297fa	440	} else {
572b0139	441	step+=fNAddSteps;
b2a297fa	442	}
	443	}
	444	}
	445
	446	//Steps for increasing cut match
78091935	447	if (fStepsForCutsIncreasing&&fNCuts>2){
b2a297fa	448	for (Int_t iCut=1; iCut<fNCuts-1;++iCut) {
b2a297fa	449	if (mask&(1<<((iCut+1)-1))) {
6551594b	450	fCfContainer->Fill(fValues,step);
b2a297fa	451	++step;
572b0139	452
	453	if (fHasMC){
	454	if ( fStepsForSignal){
	455	if (isMCTruth) fCfContainer->Fill(fValues,step);
	456	++step;
	457	}
	458	if ( fStepsForBackground ){
	459	if (!isMCTruth) fCfContainer->Fill(fValues,step);
	460	++step;
	461	}
	462	}
b2a297fa	463	} else {
572b0139	464	step+=fNAddSteps;
b2a297fa	465	}
	466	}
	467	}
	468
	469	//Steps of user defined cut combinations
	470	for (UInt_t iComb=0; iComb<fNStepMasks; ++iComb){
	471	UInt_t userMask=fStepMasks[iComb];
	472	if (mask&userMask) {
6551594b	473	fCfContainer->Fill(fValues,step);
b2a297fa	474	++step;
572b0139	475
	476	if (fHasMC){
	477	if ( fStepsForSignal){
	478	if (isMCTruth) fCfContainer->Fill(fValues,step);
	479	++step;
	480	}
	481	if ( fStepsForBackground ){
	482	if (!isMCTruth) fCfContainer->Fill(fValues,step);
	483	++step;
	484	}
	485	}
	486	} else {
	487	step+=fNAddSteps;
	488	}
	489	}
	490
	491	//All cuts
	492	if (fStepForAfterAllCuts){
	493	if (mask == selectedMask){
	494	fCfContainer->Fill(fValues,step);
b2a297fa	495	++step;
572b0139	496
	497	if (fHasMC){
	498	if ( fStepsForSignal){
	499	if (isMCTruth) fCfContainer->Fill(fValues,step);
	500	++step;
	501	}
	502	if ( fStepsForBackground ){
	503	if (!isMCTruth) fCfContainer->Fill(fValues,step);
	504	++step;
	505	}
	506	}
b2a297fa	507	} else {
572b0139	508	step+=fNAddSteps;
b2a297fa	509	}
	510	}
	511
572b0139	512	if (step!=fNSteps) {
	513	AliError("Something went wrong in the step filling!!!");
	514	}
	515
b2a297fa	516	}
	517
	518	//________________________________________________________________
	519	void AliDielectronCF::FillMC(const TObject *particle)
	520	{
	521	//
	522	// fill MC part of the Container
	523	//
e123f993	524	if (!fStepForMCtruth) return;
e123f993	525
6551594b	526	Double_t valuesPair[AliDielectronVarManager::kNMaxValues];
6551594b	527	AliDielectronVarManager::Fill(particle,valuesPair);
8df8e382	528
	529	AliVParticle *d1=0x0;
	530	AliVParticle *d2=0x0;
	531	AliDielectronMC::Instance()->GetDaughters(particle,d1,d2);
	532
	533	valuesPair[AliDielectronVarManager::kThetaHE]=AliDielectronPair::ThetaPhiCM(d1,d2,kTRUE,kTRUE);
	534	valuesPair[AliDielectronVarManager::kPhiHE]=AliDielectronPair::ThetaPhiCM(d1,d2,kTRUE,kFALSE);
	535	valuesPair[AliDielectronVarManager::kThetaCS]=AliDielectronPair::ThetaPhiCM(d1,d2,kFALSE,kTRUE);
	536	valuesPair[AliDielectronVarManager::kPhiCS]=AliDielectronPair::ThetaPhiCM(d1,d2,kFALSE,kFALSE);
6551594b	537
a655b716	538	//TODO: temporary solution, set manually the pair type to 1: unlikesign SE
6551594b	539	valuesPair[AliDielectronVarManager::kPairType]=1;
b2a297fa	540
b2a297fa	541	for (Int_t iVar=0; iVar<fNVars; ++iVar){
b2a297fa	542	Int_t var=fVariables[iVar];
6551594b	543	fValues[iVar]=valuesPair[var];
b2a297fa	544	}
6551594b	545
6551594b	546	if (fNVarsLeg>0){
6551594b	547	Double_t valuesLeg1[AliDielectronVarManager::kNMaxValues];
6551594b	548	Double_t valuesLeg2[AliDielectronVarManager::kNMaxValues];
a655b716	549	if (d1->Pt()>d2->Pt()){
	550	AliDielectronVarManager::Fill(d1,valuesLeg1);
	551	AliDielectronVarManager::Fill(d2,valuesLeg2);
	552	} else {
	553	AliDielectronVarManager::Fill(d2,valuesLeg1);
	554	AliDielectronVarManager::Fill(d1,valuesLeg2);
	555	}
6551594b	556
	557	for (Int_t iVar=0; iVar<fNVarsLeg; ++iVar){
	558	Int_t var=fVariablesLeg[iVar];
	559	fValues[iVar+fNVars]=valuesLeg1[var];
	560	fValues[iVar+fNVars+fNVarsLeg]=valuesLeg2[var];
	561	}
	562	}
	563
	564	fCfContainer->Fill(fValues,0);
b2a297fa	565	}
b2a297fa	566
3505bfad	567	//_____________________________________________________________________________
	568	TVectorD* AliDielectronCF::MakeLogBinning(Int_t nbinsX, Double_t xmin, Double_t xmax) const
	569	{
	570	//
	571	// Make logarithmic binning
	572	// the user has to delete the array afterwards!!!
	573	//
	574
	575	//check limits
	576	if (xmin<1e-20 \|\| xmax<1e-20){
	577	AliError("For Log binning xmin and xmax must be > 1e-20. Using linear binning instead!");
	578	return MakeLinBinning(nbinsX, xmin, xmax);
	579	}
	580	if (xmax<xmin){
	581	Double_t tmp=xmin;
	582	xmin=xmax;
	583	xmax=tmp;
	584	}
	585	TVectorD *binLim=new TVectorD(nbinsX+1);
	586	Double_t first=xmin;
	587	Double_t last=xmax;
	588	Double_t expMax=TMath::Log(last/first);
	589	for (Int_t i=0; i<nbinsX+1; ++i){
	590	(binLim)[i]=firstTMath::Exp(expMax/nbinsX*(Double_t)i);
	591	}
	592	return binLim;
	593	}
	594
	595	//_____________________________________________________________________________
	596	TVectorD* AliDielectronCF::MakeLinBinning(Int_t nbinsX, Double_t xmin, Double_t xmax) const
	597	{
	598	//
	599	// Make logarithmic binning
	600	// the user has to delete the array afterwards!!!
	601	//
	602	if (xmax<xmin){
	603	Double_t tmp=xmin;
	604	xmin=xmax;
	605	xmax=tmp;
	606	}
	607	TVectorD *binLim=new TVectorD(nbinsX+1);
	608	Double_t first=xmin;
	609	Double_t last=xmax;
	610	Double_t binWidth=(last-first)/nbinsX;
	611	for (Int_t i=0; i<nbinsX+1; ++i){
	612	(binLim)[i]=first+binWidth(Double_t)i;
	613	}
	614	return binLim;
	615	}
	616