[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliRDHFCuts.cxx

/**************************************************************************
 * Copyright(c) 1998-2010, 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.                  *
 **************************************************************************/

/////////////////////////////////////////////////////////////
//
// Base class for cuts on AOD reconstructed heavy-flavour decay
//
// Author: A.Dainese, andrea.dainese@pd.infn.it
/////////////////////////////////////////////////////////////
#include <Riostream.h>

#include "AliVEvent.h"
#include "AliESDEvent.h"
#include "AliAODEvent.h"
#include "AliVVertex.h"
#include "AliESDVertex.h"
#include "AliAODVertex.h"
#include "AliESDtrack.h"
#include "AliAODTrack.h"
#include "AliESDtrackCuts.h"
#include "AliAODRecoDecayHF.h"
#include "AliRDHFCuts.h"

ClassImp(AliRDHFCuts)

//--------------------------------------------------------------------------
AliRDHFCuts::AliRDHFCuts(const Char_t* name, const Char_t* title) : 
AliAnalysisCuts(name,title),
fMinVtxType(3),
fMinVtxContr(1),
fMaxVtxRedChi2(1e6),
fMinSPDMultiplicity(0),
fTriggerMask(0),
fTrackCuts(0),
fnPtBins(1),
fnPtBinLimits(1),
fPtBinLimits(0),
fnVars(1),
fVarNames(0),
fnVarsForOpt(0),
fVarsForOpt(0),
fGlobalIndex(1),
fCutsRD(0),
fIsUpperCut(0),
fUsePID(kFALSE),
fPidHF(0),
fWhyRejection(0),
fRemoveDaughtersFromPrimary(kFALSE),
fOptPileup(0),
fMinContrPileup(3),
fMinDzPileup(0.6)
{
  //
  // Default Constructor
  //
}
//--------------------------------------------------------------------------
AliRDHFCuts::AliRDHFCuts(const AliRDHFCuts &source) :
  AliAnalysisCuts(source),
  fMinVtxType(source.fMinVtxType),
  fMinVtxContr(source.fMinVtxContr),
  fMaxVtxRedChi2(source.fMaxVtxRedChi2),
  fMinSPDMultiplicity(source.fMinSPDMultiplicity),
  fTriggerMask(source.fTriggerMask),
  fTrackCuts(0),
  fnPtBins(source.fnPtBins),
  fnPtBinLimits(source.fnPtBinLimits),
  fPtBinLimits(0),
  fnVars(source.fnVars),
  fVarNames(0),
  fnVarsForOpt(source.fnVarsForOpt),
  fVarsForOpt(0),
  fGlobalIndex(source.fGlobalIndex),
  fCutsRD(0),
  fIsUpperCut(0),
  fUsePID(source.fUsePID),
  fPidHF(0),
  fWhyRejection(source.fWhyRejection),
  fRemoveDaughtersFromPrimary(source.fRemoveDaughtersFromPrimary),
  fOptPileup(source.fOptPileup),
  fMinContrPileup(source.fMinContrPileup),
  fMinDzPileup(source.fMinDzPileup)  
{
  //
  // Copy constructor
  //
  cout<<"Copy constructor"<<endl;
  if(source.GetTrackCuts()) AddTrackCuts(source.GetTrackCuts());
  if(source.fPtBinLimits) SetPtBins(source.fnPtBinLimits,source.fPtBinLimits);
  if(source.fVarNames) SetVarNames(source.fnVars,source.fVarNames,source.fIsUpperCut);
  if(source.fCutsRD) SetCuts(source.fGlobalIndex,source.fCutsRD);
  if(source.fVarsForOpt) SetVarsForOpt(source.fnVarsForOpt,source.fVarsForOpt);
  if(source.fPidHF) SetPidHF(source.fPidHF);
  PrintAll();

}
//--------------------------------------------------------------------------
AliRDHFCuts &AliRDHFCuts::operator=(const AliRDHFCuts &source)
{
  //
  // assignment operator
  //
  if(&source == this) return *this;

  AliAnalysisCuts::operator=(source);

  fMinVtxType=source.fMinVtxType;
  fMinVtxContr=source.fMinVtxContr;
  fMaxVtxRedChi2=source.fMaxVtxRedChi2;
  fMinSPDMultiplicity=source.fMinSPDMultiplicity;
  fTriggerMask=source.fTriggerMask;
  fnPtBins=source.fnPtBins;
  fnVars=source.fnVars;
  fGlobalIndex=source.fGlobalIndex;
  fnVarsForOpt=source.fnVarsForOpt;
  fUsePID=source.fUsePID;
  SetPidHF(source.GetPidHF());
  fWhyRejection=source.fWhyRejection;
  fRemoveDaughtersFromPrimary=source.fRemoveDaughtersFromPrimary;
  fOptPileup=source.fOptPileup;
  fMinContrPileup=source.fMinContrPileup;
  fMinDzPileup=source.fMinDzPileup;

  if(source.GetTrackCuts()) AddTrackCuts(source.GetTrackCuts());
  if(source.fPtBinLimits) SetPtBins(source.fnPtBinLimits,source.fPtBinLimits);
  if(source.fVarNames) SetVarNames(source.fnVars,source.fVarNames,source.fIsUpperCut);
  if(source.fCutsRD) SetCuts(source.fGlobalIndex,source.fCutsRD);
  if(source.fVarsForOpt) SetVarsForOpt(source.fnVarsForOpt,source.fVarsForOpt);
  PrintAll();

  return *this;
}
//--------------------------------------------------------------------------
AliRDHFCuts::~AliRDHFCuts() {
  //  
  // Default Destructor
  //
  if(fPtBinLimits) {delete [] fPtBinLimits; fPtBinLimits=0;}
  if(fVarNames) {delete [] fVarNames; fVarNames=0;}
  if(fVarsForOpt) {delete [] fVarsForOpt; fVarsForOpt=0;}
  if(fCutsRD) {
    delete [] fCutsRD;
    fCutsRD=0;
  }
  if(fIsUpperCut) {delete [] fIsUpperCut; fIsUpperCut=0;}
  if(fPidHF){ 
    delete fPidHF;
    fPidHF=0;
  }
}
//---------------------------------------------------------------------------
Bool_t AliRDHFCuts::IsEventSelected(AliVEvent *event) {
  //
  // Event selection
  // 
  //if(fTriggerMask && event->GetTriggerMask()!=fTriggerMask) return kFALSE;

  fWhyRejection=0;

  // multiplicity cuts no implemented yet


  const AliVVertex *vertex = event->GetPrimaryVertex();

  if(!vertex) return kFALSE;

  TString title=vertex->GetTitle();
  if(title.Contains("Z") && fMinVtxType>1) return kFALSE; 
  if(title.Contains("3D") && fMinVtxType>2) return kFALSE; 

  if(vertex->GetNContributors()<fMinVtxContr) return kFALSE; 
  if(fOptPileup==kRejectPileupEvent){
    Int_t cutc=(Int_t)fMinContrPileup;
    Double_t cutz=(Double_t)fMinDzPileup;
    if(event->IsPileupFromSPD(cutc,cutz,3.,2.,10.)) {
      fWhyRejection=1;
      return kFALSE;
    }
  }

  return kTRUE;
}
//---------------------------------------------------------------------------
Bool_t AliRDHFCuts::AreDaughtersSelected(AliAODRecoDecayHF *d) const {
  //
  // Daughter tracks selection
  // 
  if(!fTrackCuts) return kTRUE;

  Int_t ndaughters = d->GetNDaughters();
  AliAODVertex *vAOD = d->GetPrimaryVtx();
  Double_t pos[3],cov[6];
  vAOD->GetXYZ(pos);
  vAOD->GetCovarianceMatrix(cov);
  const AliESDVertex vESD(pos,cov,100.,100);

  Bool_t retval=kTRUE;

  for(Int_t idg=0; idg<ndaughters; idg++) {
    AliAODTrack *dgTrack = (AliAODTrack*)d->GetDaughter(idg);
    if(!dgTrack) retval = kFALSE;
    //printf("charge %d\n",dgTrack->Charge());
    if(dgTrack->Charge()==0) continue; // it's not a track, but a V0

    if(!IsDaughterSelected(dgTrack,&vESD,fTrackCuts)) retval = kFALSE;
  }

  return retval;
}
//---------------------------------------------------------------------------
Bool_t AliRDHFCuts::IsDaughterSelected(AliAODTrack *track,const AliESDVertex *primary,AliESDtrackCuts *cuts) const {
  //
  // Convert to ESDtrack, relate to vertex and check cuts
  //
  if(!cuts) return kTRUE;

  Bool_t retval=kTRUE;

  // convert to ESD track here
  AliESDtrack esdTrack(track);
  // needed to calculate the impact parameters
  esdTrack.RelateToVertex(primary,0.,3.); 
  if(!cuts->IsSelected(&esdTrack)) retval = kFALSE;
 
  if(fOptPileup==kRejectTracksFromPileupVertex){
    // to be implemented
    // we need either to have here the AOD Event, 
    // or to have the pileup vertex object
  }
  return retval; 
}
//---------------------------------------------------------------------------
void AliRDHFCuts::SetPtBins(Int_t nPtBinLimits,Float_t *ptBinLimits) {
  // Set the pt bins

  if(fPtBinLimits) {
    delete [] fPtBinLimits;
    fPtBinLimits = NULL;
    printf("Changing the pt bins\n");
  }

  if(nPtBinLimits != fnPtBins+1){
    cout<<"Warning: ptBinLimits dimention "<<nPtBinLimits<<" != nPtBins+1 ("<<fnPtBins+1<<")\nSetting nPtBins to "<<nPtBinLimits-1<<endl;
    SetNPtBins(nPtBinLimits-1);
  }

  fnPtBinLimits = nPtBinLimits;
  SetGlobalIndex();
  //cout<<"Changing also Global Index -> "<<fGlobalIndex<<endl;
  fPtBinLimits = new Float_t[fnPtBinLimits];
  for(Int_t ib=0; ib<nPtBinLimits; ib++) fPtBinLimits[ib]=ptBinLimits[ib];

  return;
}
//---------------------------------------------------------------------------
void AliRDHFCuts::SetVarNames(Int_t nVars,TString *varNames,Bool_t *isUpperCut){
  // Set the variable names

  if(fVarNames) {
    delete [] fVarNames;
    fVarNames = NULL;
    //printf("Changing the variable names\n");
  }
  if(nVars!=fnVars){
    printf("Wrong number of variables: it has to be %d\n",fnVars);
    return;
  }
  //fnVars=nVars;
  fVarNames = new TString[nVars];
  fIsUpperCut = new Bool_t[nVars];
  for(Int_t iv=0; iv<nVars; iv++) {
    fVarNames[iv] = varNames[iv];
    fIsUpperCut[iv] = isUpperCut[iv];
  }

  return;
}
//---------------------------------------------------------------------------
void AliRDHFCuts::SetVarsForOpt(Int_t nVars,Bool_t *forOpt) {
  // Set the variables to be used for cuts optimization

  if(fVarsForOpt) {
    delete [] fVarsForOpt;
    fVarsForOpt = NULL;
    //printf("Changing the variables for cut optimization\n");
  }
  
  if(nVars==0){//!=fnVars) {
    printf("%d not accepted as number of variables: it has to be %d\n",nVars,fnVars);
    return;
  } 
  
  fnVarsForOpt = 0;
  fVarsForOpt = new Bool_t[fnVars];
  for(Int_t iv=0; iv<fnVars; iv++) {
    fVarsForOpt[iv]=forOpt[iv];
    if(fVarsForOpt[iv]) fnVarsForOpt++;
  }

  return;
}
//---------------------------------------------------------------------------
void AliRDHFCuts::SetCuts(Int_t nVars,Int_t nPtBins,Float_t **cutsRD) {
  //
  // store the cuts
  //
  if(nVars!=fnVars) {
    printf("Wrong number of variables: it has to be %d\n",fnVars);
    return;
  } 
  if(nPtBins!=fnPtBins) {
    printf("Wrong number of pt bins: it has to be %d\n",fnPtBins);
    return;
  } 

  if(!fCutsRD)  fCutsRD = new Float_t[fGlobalIndex];
  

  for(Int_t iv=0; iv<fnVars; iv++) {

    for(Int_t ib=0; ib<fnPtBins; ib++) {

      //check
      if(GetGlobalIndex(iv,ib)>=fGlobalIndex) {
	cout<<"Overflow, exit..."<<endl;
	return;
      }

      fCutsRD[GetGlobalIndex(iv,ib)] = cutsRD[iv][ib];

    }
  }
  return;
}
//---------------------------------------------------------------------------
void AliRDHFCuts::SetCuts(Int_t glIndex,Float_t* cutsRDGlob){
  //
  // store the cuts
  //
  if(glIndex != fGlobalIndex){
    cout<<"Wrong array size: it has to be "<<fGlobalIndex<<endl;
    return;
  }
  if(!fCutsRD)  fCutsRD = new Float_t[fGlobalIndex];

  for(Int_t iGl=0;iGl<fGlobalIndex;iGl++){
    fCutsRD[iGl] = cutsRDGlob[iGl];
  }
  return;
}
//---------------------------------------------------------------------------
void AliRDHFCuts::PrintAll() const {
  //
  // print all cuts values
  // 

  printf("Minimum vtx type %d\n",fMinVtxType);
  printf("Minimum vtx contr %d\n",fMinVtxContr);
  printf("Max vtx red chi2 %f\n",fMaxVtxRedChi2);
  printf("Min SPD mult %d\n",fMinSPDMultiplicity);
  printf("Use PID %d\n",(Int_t)fUsePID);
  printf("Remove daughters from vtx %d\n",(Int_t)fRemoveDaughtersFromPrimary);

  if(fVarNames){
    cout<<"Array of variables"<<endl;
    for(Int_t iv=0;iv<fnVars;iv++){
      cout<<fVarNames[iv]<<"\t";
    }
    cout<<endl;
  }
  if(fVarsForOpt){
    cout<<"Array of optimization"<<endl;
    for(Int_t iv=0;iv<fnVars;iv++){
      cout<<fVarsForOpt[iv]<<"\t";
    }
    cout<<endl;
  }
  if(fIsUpperCut){
    cout<<"Array of upper/lower cut"<<endl;
   for(Int_t iv=0;iv<fnVars;iv++){
     cout<<fIsUpperCut[iv]<<"\t";
   }
   cout<<endl;
  }
  if(fPtBinLimits){
    cout<<"Array of ptbin limits"<<endl;
    for(Int_t ib=0;ib<fnPtBinLimits;ib++){
      cout<<fPtBinLimits[ib]<<"\t";
    }
    cout<<endl;
  }
  if(fCutsRD){
    cout<<"Matrix of cuts"<<endl;
   for(Int_t iv=0;iv<fnVars;iv++){
     for(Int_t ib=0;ib<fnPtBins;ib++){
       cout<<"fCutsRD["<<iv<<"]["<<ib<<"] = "<<fCutsRD[GetGlobalIndex(iv,ib)]<<"\t";
     }
     cout<<endl;
   }
   cout<<endl;
  }
  return;
}
//---------------------------------------------------------------------------
void AliRDHFCuts::GetCuts(Float_t**& cutsRD) const{
  //
  // get the cuts
  //

  //cout<<"Give back a "<<fnVars<<"x"<<fnPtBins<<" matrix."<<endl;


  Int_t iv,ib;
  if(!cutsRD) {
    //cout<<"Initialization..."<<endl;
    cutsRD=new Float_t*[fnVars];
    for(iv=0; iv<fnVars; iv++) {
      cutsRD[iv] = new Float_t[fnPtBins];
    }
  }
  
  for(Int_t iGlobal=0; iGlobal<fGlobalIndex; iGlobal++) {
    GetVarPtIndex(iGlobal,iv,ib);
    cutsRD[iv][ib] = fCutsRD[iGlobal];
  }

  return;
}

//---------------------------------------------------------------------------
Int_t AliRDHFCuts::GetGlobalIndex(Int_t iVar,Int_t iPtBin) const{
  //
  // give the global index from variable and pt bin
  //
  return iPtBin*fnVars+iVar;
}

//---------------------------------------------------------------------------
void AliRDHFCuts::GetVarPtIndex(Int_t iGlob, Int_t& iVar, Int_t& iPtBin) const {
  //
  //give the index of the variable and of the pt bin from the global index
  //
  iPtBin=(Int_t)iGlob/fnVars;
  iVar=iGlob%fnVars;

  return;
}

//---------------------------------------------------------------------------
Int_t AliRDHFCuts::PtBin(Double_t pt) const {
  //
  //give the pt bin where the pt lies.
  //
  Int_t ptbin=-1;
  for (Int_t i=0;i<fnPtBins;i++){
    if(pt<fPtBinLimits[i+1]) {
      ptbin=i;
      break;
    }
  }
  return ptbin;
}
//-------------------------------------------------------------------
Float_t AliRDHFCuts::GetCutValue(Int_t iVar,Int_t iPtBin) const {
  // 
  // Give the value of cut set for the variable iVar and the pt bin iPtBin
  //
  if(!fCutsRD){
    cout<<"Cuts not iniziaisez yet"<<endl;
    return 0;
  }
  return fCutsRD[GetGlobalIndex(iVar,iPtBin)];
}
//-------------------------------------------------------------------
Bool_t AliRDHFCuts::CompareCuts(const AliRDHFCuts *obj) const {
  //
  // Compare two cuts objects
  //

  Bool_t areEqual=kTRUE;

  if(fMinVtxType!=obj->fMinVtxType) { printf("Minimum vtx type %d  %d\n",fMinVtxType,obj->fMinVtxType); areEqual=kFALSE;}

  if(fMinVtxContr!=obj->fMinVtxContr) { printf("Minimum vtx contr %d  %d\n",fMinVtxContr,obj->fMinVtxContr); areEqual=kFALSE;}

  if(TMath::Abs(fMaxVtxRedChi2-obj->fMaxVtxRedChi2)>1.e-10) {   printf("Max vtx red chi2 %f  %f\n",fMaxVtxRedChi2,obj->fMaxVtxRedChi2);areEqual=kFALSE;}

  if(fMinSPDMultiplicity!=obj->fMinSPDMultiplicity) {  printf("Min SPD mult %d\n  %d",fMinSPDMultiplicity,obj->fMinSPDMultiplicity);areEqual=kFALSE;}

  if(fUsePID!=obj->fUsePID) { printf("Use PID %d  %d\n",(Int_t)fUsePID,(Int_t)obj->fUsePID); areEqual=kFALSE;}

  if(fRemoveDaughtersFromPrimary!=obj->fRemoveDaughtersFromPrimary) {printf("Remove daughters from vtx %d  %d\n",(Int_t)fRemoveDaughtersFromPrimary,(Int_t)obj->fRemoveDaughtersFromPrimary); areEqual=kFALSE;}
  if(fTrackCuts){
    if(fTrackCuts->GetMinNClusterTPC()!=obj->fTrackCuts->GetMinNClusterTPC()) {printf("MinNClsTPC %d  %d\n",fTrackCuts->GetMinNClusterTPC(),obj->fTrackCuts->GetMinNClusterTPC()); areEqual=kFALSE;}

    if(fTrackCuts->GetMinNClustersITS()!=obj->fTrackCuts->GetMinNClustersITS()) {printf("MinNClsITS %d  %d\n",fTrackCuts->GetMinNClustersITS(),obj->fTrackCuts->GetMinNClustersITS()); areEqual=kFALSE;}

    if(TMath::Abs(fTrackCuts->GetMaxChi2PerClusterTPC()-obj->fTrackCuts->GetMaxChi2PerClusterTPC())>1.e-10) {printf("MaxChi2ClsTPC %f  %f\n",fTrackCuts->GetMaxChi2PerClusterTPC(),obj->fTrackCuts->GetMaxChi2PerClusterTPC()); areEqual=kFALSE;}

    if(fTrackCuts->GetClusterRequirementITS(AliESDtrackCuts::kSPD)!=obj->fTrackCuts->GetClusterRequirementITS(AliESDtrackCuts::kSPD)) {printf("ClusterReq SPD %d  %d\n",fTrackCuts->GetClusterRequirementITS(AliESDtrackCuts::kSPD),obj->fTrackCuts->GetClusterRequirementITS(AliESDtrackCuts::kSPD)); areEqual=kFALSE;}
  }

  if(fCutsRD) {
   for(Int_t iv=0;iv<fnVars;iv++) {
     for(Int_t ib=0;ib<fnPtBins;ib++) {
       if(TMath::Abs(fCutsRD[GetGlobalIndex(iv,ib)]-obj->fCutsRD[GetGlobalIndex(iv,ib)])>1.e-10) {
	 cout<<"fCutsRD["<<iv<<"]["<<ib<<"] = "<<fCutsRD[GetGlobalIndex(iv,ib)]<<"   "<<obj->fCutsRD[GetGlobalIndex(iv,ib)]<<"\n";
	 areEqual=kFALSE;
       }
     }
   }
  }

  return areEqual;
}
Commit	Line	Data
650b3ced	1	/**************************************************************************
	2	* Copyright(c) 1998-2010, 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	//
	18	// Base class for cuts on AOD reconstructed heavy-flavour decay
	19	//
	20	// Author: A.Dainese, andrea.dainese@pd.infn.it
	21	/////////////////////////////////////////////////////////////
	22	#include <Riostream.h>
	23
	24	#include "AliVEvent.h"
	25	#include "AliESDEvent.h"
	26	#include "AliAODEvent.h"
	27	#include "AliVVertex.h"
	28	#include "AliESDVertex.h"
	29	#include "AliAODVertex.h"
3cbc09cd	30	#include "AliESDtrack.h"
3cbc09cd	31	#include "AliAODTrack.h"
650b3ced	32	#include "AliESDtrackCuts.h"
	33	#include "AliAODRecoDecayHF.h"
	34	#include "AliRDHFCuts.h"
	35
	36	ClassImp(AliRDHFCuts)
	37
	38	//--------------------------------------------------------------------------
	39	AliRDHFCuts::AliRDHFCuts(const Char_t* name, const Char_t* title) :
	40	AliAnalysisCuts(name,title),
	41	fMinVtxType(3),
	42	fMinVtxContr(1),
	43	fMaxVtxRedChi2(1e6),
	44	fMinSPDMultiplicity(0),
	45	fTriggerMask(0),
	46	fTrackCuts(0),
	47	fnPtBins(1),
	48	fnPtBinLimits(1),
	49	fPtBinLimits(0),
	50	fnVars(1),
	51	fVarNames(0),
	52	fnVarsForOpt(0),
	53	fVarsForOpt(0),
	54	fGlobalIndex(1),
	55	fCutsRD(0),
c96634a2	56	fIsUpperCut(0),
fa144fd0	57	fUsePID(kFALSE),
fa144fd0	58	fPidHF(0),
918b6831	59	fWhyRejection(0),
e4f44447	60	fRemoveDaughtersFromPrimary(kFALSE),
	61	fOptPileup(0),
	62	fMinContrPileup(3),
	63	fMinDzPileup(0.6)
650b3ced	64	{
	65	//
	66	// Default Constructor
	67	//
	68	}
	69	//--------------------------------------------------------------------------
	70	AliRDHFCuts::AliRDHFCuts(const AliRDHFCuts &source) :
	71	AliAnalysisCuts(source),
	72	fMinVtxType(source.fMinVtxType),
	73	fMinVtxContr(source.fMinVtxContr),
	74	fMaxVtxRedChi2(source.fMaxVtxRedChi2),
	75	fMinSPDMultiplicity(source.fMinSPDMultiplicity),
	76	fTriggerMask(source.fTriggerMask),
	77	fTrackCuts(0),
	78	fnPtBins(source.fnPtBins),
	79	fnPtBinLimits(source.fnPtBinLimits),
	80	fPtBinLimits(0),
	81	fnVars(source.fnVars),
	82	fVarNames(0),
	83	fnVarsForOpt(source.fnVarsForOpt),
	84	fVarsForOpt(0),
	85	fGlobalIndex(source.fGlobalIndex),
	86	fCutsRD(0),
c96634a2	87	fIsUpperCut(0),
fa144fd0	88	fUsePID(source.fUsePID),
fa144fd0	89	fPidHF(0),
918b6831	90	fWhyRejection(source.fWhyRejection),
e4f44447	91	fRemoveDaughtersFromPrimary(source.fRemoveDaughtersFromPrimary),
	92	fOptPileup(source.fOptPileup),
	93	fMinContrPileup(source.fMinContrPileup),
	94	fMinDzPileup(source.fMinDzPileup)
650b3ced	95	{
	96	//
	97	// Copy constructor
	98	//
	99	cout<<"Copy constructor"<<endl;
	100	if(source.GetTrackCuts()) AddTrackCuts(source.GetTrackCuts());
650b3ced	101	if(source.fPtBinLimits) SetPtBins(source.fnPtBinLimits,source.fPtBinLimits);
	102	if(source.fVarNames) SetVarNames(source.fnVars,source.fVarNames,source.fIsUpperCut);
	103	if(source.fCutsRD) SetCuts(source.fGlobalIndex,source.fCutsRD);
	104	if(source.fVarsForOpt) SetVarsForOpt(source.fnVarsForOpt,source.fVarsForOpt);
fa144fd0	105	if(source.fPidHF) SetPidHF(source.fPidHF);
650b3ced	106	PrintAll();
650b3ced	107
	108	}
	109	//--------------------------------------------------------------------------
	110	AliRDHFCuts &AliRDHFCuts::operator=(const AliRDHFCuts &source)
	111	{
	112	//
	113	// assignment operator
	114	//
	115	if(&source == this) return *this;
	116
	117	AliAnalysisCuts::operator=(source);
	118
	119	fMinVtxType=source.fMinVtxType;
	120	fMinVtxContr=source.fMinVtxContr;
	121	fMaxVtxRedChi2=source.fMaxVtxRedChi2;
	122	fMinSPDMultiplicity=source.fMinSPDMultiplicity;
	123	fTriggerMask=source.fTriggerMask;
	124	fnPtBins=source.fnPtBins;
	125	fnVars=source.fnVars;
	126	fGlobalIndex=source.fGlobalIndex;
	127	fnVarsForOpt=source.fnVarsForOpt;
c96634a2	128	fUsePID=source.fUsePID;
fa144fd0	129	SetPidHF(source.GetPidHF());
fa144fd0	130	fWhyRejection=source.fWhyRejection;
918b6831	131	fRemoveDaughtersFromPrimary=source.fRemoveDaughtersFromPrimary;
e4f44447	132	fOptPileup=source.fOptPileup;
	133	fMinContrPileup=source.fMinContrPileup;
	134	fMinDzPileup=source.fMinDzPileup;
650b3ced	135
bbbf8c32	136	if(source.GetTrackCuts()) AddTrackCuts(source.GetTrackCuts());
650b3ced	137	if(source.fPtBinLimits) SetPtBins(source.fnPtBinLimits,source.fPtBinLimits);
	138	if(source.fVarNames) SetVarNames(source.fnVars,source.fVarNames,source.fIsUpperCut);
	139	if(source.fCutsRD) SetCuts(source.fGlobalIndex,source.fCutsRD);
	140	if(source.fVarsForOpt) SetVarsForOpt(source.fnVarsForOpt,source.fVarsForOpt);
	141	PrintAll();
	142
650b3ced	143	return *this;
	144	}
	145	//--------------------------------------------------------------------------
	146	AliRDHFCuts::~AliRDHFCuts() {
	147	//
	148	// Default Destructor
	149	//
	150	if(fPtBinLimits) {delete [] fPtBinLimits; fPtBinLimits=0;}
	151	if(fVarNames) {delete [] fVarNames; fVarNames=0;}
	152	if(fVarsForOpt) {delete [] fVarsForOpt; fVarsForOpt=0;}
	153	if(fCutsRD) {
	154	delete [] fCutsRD;
	155	fCutsRD=0;
	156	}
	157	if(fIsUpperCut) {delete [] fIsUpperCut; fIsUpperCut=0;}
fa144fd0	158	if(fPidHF){
	159	delete fPidHF;
	160	fPidHF=0;
	161	}
650b3ced	162	}
650b3ced	163	//---------------------------------------------------------------------------
fb38b6dc	164	Bool_t AliRDHFCuts::IsEventSelected(AliVEvent *event) {
650b3ced	165	//
	166	// Event selection
	167	//
3cbc09cd	168	//if(fTriggerMask && event->GetTriggerMask()!=fTriggerMask) return kFALSE;
650b3ced	169
fb38b6dc	170	fWhyRejection=0;
fb38b6dc	171
650b3ced	172	// multiplicity cuts no implemented yet
650b3ced	173
fb38b6dc	174
fb38b6dc	175
650b3ced	176	const AliVVertex *vertex = event->GetPrimaryVertex();
	177
	178	if(!vertex) return kFALSE;
	179
	180	TString title=vertex->GetTitle();
	181	if(title.Contains("Z") && fMinVtxType>1) return kFALSE;
	182	if(title.Contains("3D") && fMinVtxType>2) return kFALSE;
	183
	184	if(vertex->GetNContributors()<fMinVtxContr) return kFALSE;
e4f44447	185	if(fOptPileup==kRejectPileupEvent){
	186	Int_t cutc=(Int_t)fMinContrPileup;
	187	Double_t cutz=(Double_t)fMinDzPileup;
fb38b6dc	188	if(event->IsPileupFromSPD(cutc,cutz,3.,2.,10.)) {
	189	fWhyRejection=1;
	190	return kFALSE;
	191	}
e4f44447	192	}
650b3ced	193
	194	return kTRUE;
	195	}
	196	//---------------------------------------------------------------------------
3cbc09cd	197	Bool_t AliRDHFCuts::AreDaughtersSelected(AliAODRecoDecayHF *d) const {
	198	//
	199	// Daughter tracks selection
	200	//
85d335a7	201	if(!fTrackCuts) return kTRUE;
85d335a7	202
3cbc09cd	203	Int_t ndaughters = d->GetNDaughters();
85d335a7	204	AliAODVertex *vAOD = d->GetPrimaryVtx();
	205	Double_t pos[3],cov[6];
	206	vAOD->GetXYZ(pos);
	207	vAOD->GetCovarianceMatrix(cov);
da6fefc3	208	const AliESDVertex vESD(pos,cov,100.,100);
85d335a7	209
85d335a7	210	Bool_t retval=kTRUE;
3cbc09cd	211
	212	for(Int_t idg=0; idg<ndaughters; idg++) {
	213	AliAODTrack dgTrack = (AliAODTrack)d->GetDaughter(idg);
85d335a7	214	if(!dgTrack) retval = kFALSE;
3cbc09cd	215	//printf("charge %d\n",dgTrack->Charge());
3cbc09cd	216	if(dgTrack->Charge()==0) continue; // it's not a track, but a V0
85d335a7	217
da6fefc3	218	if(!IsDaughterSelected(dgTrack,&vESD,fTrackCuts)) retval = kFALSE;
da6fefc3	219	}
85d335a7	220
85d335a7	221	return retval;
3cbc09cd	222	}
3cbc09cd	223	//---------------------------------------------------------------------------
da6fefc3	224	Bool_t AliRDHFCuts::IsDaughterSelected(AliAODTrack track,const AliESDVertex primary,AliESDtrackCuts *cuts) const {
	225	//
	226	// Convert to ESDtrack, relate to vertex and check cuts
	227	//
	228	if(!cuts) return kTRUE;
	229
	230	Bool_t retval=kTRUE;
	231
	232	// convert to ESD track here
	233	AliESDtrack esdTrack(track);
	234	// needed to calculate the impact parameters
	235	esdTrack.RelateToVertex(primary,0.,3.);
	236	if(!cuts->IsSelected(&esdTrack)) retval = kFALSE;
	237
e4f44447	238	if(fOptPileup==kRejectTracksFromPileupVertex){
	239	// to be implemented
	240	// we need either to have here the AOD Event,
	241	// or to have the pileup vertex object
	242	}
da6fefc3	243	return retval;
	244	}
	245	//---------------------------------------------------------------------------
650b3ced	246	void AliRDHFCuts::SetPtBins(Int_t nPtBinLimits,Float_t *ptBinLimits) {
	247	// Set the pt bins
	248
	249	if(fPtBinLimits) {
	250	delete [] fPtBinLimits;
	251	fPtBinLimits = NULL;
	252	printf("Changing the pt bins\n");
	253	}
	254
	255	if(nPtBinLimits != fnPtBins+1){
	256	cout<<"Warning: ptBinLimits dimention "<<nPtBinLimits<<" != nPtBins+1 ("<<fnPtBins+1<<")\nSetting nPtBins to "<<nPtBinLimits-1<<endl;
	257	SetNPtBins(nPtBinLimits-1);
	258	}
	259
	260	fnPtBinLimits = nPtBinLimits;
	261	SetGlobalIndex();
6e0903f0	262	//cout<<"Changing also Global Index -> "<<fGlobalIndex<<endl;
650b3ced	263	fPtBinLimits = new Float_t[fnPtBinLimits];
	264	for(Int_t ib=0; ib<nPtBinLimits; ib++) fPtBinLimits[ib]=ptBinLimits[ib];
	265
	266	return;
	267	}
	268	//---------------------------------------------------------------------------
	269	void AliRDHFCuts::SetVarNames(Int_t nVars,TString varNames,Bool_t isUpperCut){
	270	// Set the variable names
	271
	272	if(fVarNames) {
	273	delete [] fVarNames;
	274	fVarNames = NULL;
6e0903f0	275	//printf("Changing the variable names\n");
650b3ced	276	}
	277	if(nVars!=fnVars){
	278	printf("Wrong number of variables: it has to be %d\n",fnVars);
	279	return;
	280	}
	281	//fnVars=nVars;
	282	fVarNames = new TString[nVars];
	283	fIsUpperCut = new Bool_t[nVars];
	284	for(Int_t iv=0; iv<nVars; iv++) {
	285	fVarNames[iv] = varNames[iv];
	286	fIsUpperCut[iv] = isUpperCut[iv];
	287	}
	288
	289	return;
	290	}
	291	//---------------------------------------------------------------------------
	292	void AliRDHFCuts::SetVarsForOpt(Int_t nVars,Bool_t *forOpt) {
	293	// Set the variables to be used for cuts optimization
	294
	295	if(fVarsForOpt) {
	296	delete [] fVarsForOpt;
	297	fVarsForOpt = NULL;
6e0903f0	298	//printf("Changing the variables for cut optimization\n");
650b3ced	299	}
	300
	301	if(nVars==0){//!=fnVars) {
	302	printf("%d not accepted as number of variables: it has to be %d\n",nVars,fnVars);
	303	return;
	304	}
	305
	306	fnVarsForOpt = 0;
	307	fVarsForOpt = new Bool_t[fnVars];
	308	for(Int_t iv=0; iv<fnVars; iv++) {
	309	fVarsForOpt[iv]=forOpt[iv];
	310	if(fVarsForOpt[iv]) fnVarsForOpt++;
	311	}
	312
	313	return;
	314	}
	315	//---------------------------------------------------------------------------
	316	void AliRDHFCuts::SetCuts(Int_t nVars,Int_t nPtBins,Float_t **cutsRD) {
4edfaee0	317	//
	318	// store the cuts
	319	//
650b3ced	320	if(nVars!=fnVars) {
	321	printf("Wrong number of variables: it has to be %d\n",fnVars);
	322	return;
	323	}
	324	if(nPtBins!=fnPtBins) {
	325	printf("Wrong number of pt bins: it has to be %d\n",fnPtBins);
	326	return;
	327	}
	328
	329	if(!fCutsRD) fCutsRD = new Float_t[fGlobalIndex];
	330
4edfaee0	331
650b3ced	332	for(Int_t iv=0; iv<fnVars; iv++) {
	333
	334	for(Int_t ib=0; ib<fnPtBins; ib++) {
	335
650b3ced	336	//check
4edfaee0	337	if(GetGlobalIndex(iv,ib)>=fGlobalIndex) {
650b3ced	338	cout<<"Overflow, exit..."<<endl;
	339	return;
	340	}
4edfaee0	341
	342	fCutsRD[GetGlobalIndex(iv,ib)] = cutsRD[iv][ib];
	343
650b3ced	344	}
650b3ced	345	}
650b3ced	346	return;
	347	}
	348	//---------------------------------------------------------------------------
	349	void AliRDHFCuts::SetCuts(Int_t glIndex,Float_t* cutsRDGlob){
4edfaee0	350	//
	351	// store the cuts
	352	//
650b3ced	353	if(glIndex != fGlobalIndex){
	354	cout<<"Wrong array size: it has to be "<<fGlobalIndex<<endl;
	355	return;
	356	}
	357	if(!fCutsRD) fCutsRD = new Float_t[fGlobalIndex];
	358
	359	for(Int_t iGl=0;iGl<fGlobalIndex;iGl++){
	360	fCutsRD[iGl] = cutsRDGlob[iGl];
	361	}
650b3ced	362	return;
650b3ced	363	}
650b3ced	364	//---------------------------------------------------------------------------
4edfaee0	365	void AliRDHFCuts::PrintAll() const {
	366	//
	367	// print all cuts values
	368	//
0198f78c	369
	370	printf("Minimum vtx type %d\n",fMinVtxType);
	371	printf("Minimum vtx contr %d\n",fMinVtxContr);
	372	printf("Max vtx red chi2 %f\n",fMaxVtxRedChi2);
	373	printf("Min SPD mult %d\n",fMinSPDMultiplicity);
	374	printf("Use PID %d\n",(Int_t)fUsePID);
	375	printf("Remove daughters from vtx %d\n",(Int_t)fRemoveDaughtersFromPrimary);
	376
650b3ced	377	if(fVarNames){
	378	cout<<"Array of variables"<<endl;
	379	for(Int_t iv=0;iv<fnVars;iv++){
	380	cout<<fVarNames[iv]<<"\t";
	381	}
	382	cout<<endl;
	383	}
	384	if(fVarsForOpt){
	385	cout<<"Array of optimization"<<endl;
	386	for(Int_t iv=0;iv<fnVars;iv++){
	387	cout<<fVarsForOpt[iv]<<"\t";
	388	}
	389	cout<<endl;
	390	}
	391	if(fIsUpperCut){
	392	cout<<"Array of upper/lower cut"<<endl;
	393	for(Int_t iv=0;iv<fnVars;iv++){
	394	cout<<fIsUpperCut[iv]<<"\t";
	395	}
	396	cout<<endl;
	397	}
	398	if(fPtBinLimits){
	399	cout<<"Array of ptbin limits"<<endl;
	400	for(Int_t ib=0;ib<fnPtBinLimits;ib++){
	401	cout<<fPtBinLimits[ib]<<"\t";
	402	}
	403	cout<<endl;
	404	}
	405	if(fCutsRD){
	406	cout<<"Matrix of cuts"<<endl;
	407	for(Int_t iv=0;iv<fnVars;iv++){
	408	for(Int_t ib=0;ib<fnPtBins;ib++){
	409	cout<<"fCutsRD["<<iv<<"]["<<ib<<"] = "<<fCutsRD[GetGlobalIndex(iv,ib)]<<"\t";
	410	}
	411	cout<<endl;
	412	}
	413	cout<<endl;
	414	}
4edfaee0	415	return;
650b3ced	416	}
650b3ced	417	//---------------------------------------------------------------------------
650b3ced	418	void AliRDHFCuts::GetCuts(Float_t**& cutsRD) const{
4edfaee0	419	//
650b3ced	420	// get the cuts
4edfaee0	421	//
650b3ced	422
4edfaee0	423	//cout<<"Give back a "<<fnVars<<"x"<<fnPtBins<<" matrix."<<endl;
650b3ced	424
650b3ced	425
4edfaee0	426	Int_t iv,ib;
650b3ced	427	if(!cutsRD) {
4edfaee0	428	//cout<<"Initialization..."<<endl;
650b3ced	429	cutsRD=new Float_t*[fnVars];
4edfaee0	430	for(iv=0; iv<fnVars; iv++) {
650b3ced	431	cutsRD[iv] = new Float_t[fnPtBins];
	432	}
	433	}
	434
4edfaee0	435	for(Int_t iGlobal=0; iGlobal<fGlobalIndex; iGlobal++) {
	436	GetVarPtIndex(iGlobal,iv,ib);
	437	cutsRD[iv][ib] = fCutsRD[iGlobal];
650b3ced	438	}
650b3ced	439
4edfaee0	440	return;
650b3ced	441	}
	442
	443	//---------------------------------------------------------------------------
	444	Int_t AliRDHFCuts::GetGlobalIndex(Int_t iVar,Int_t iPtBin) const{
4edfaee0	445	//
3cbc09cd	446	// give the global index from variable and pt bin
4edfaee0	447	//
650b3ced	448	return iPtBin*fnVars+iVar;
	449	}
	450
	451	//---------------------------------------------------------------------------
	452	void AliRDHFCuts::GetVarPtIndex(Int_t iGlob, Int_t& iVar, Int_t& iPtBin) const {
4edfaee0	453	//
650b3ced	454	//give the index of the variable and of the pt bin from the global index
4edfaee0	455	//
650b3ced	456	iPtBin=(Int_t)iGlob/fnVars;
	457	iVar=iGlob%fnVars;
	458
4edfaee0	459	return;
650b3ced	460	}
	461
	462	//---------------------------------------------------------------------------
	463	Int_t AliRDHFCuts::PtBin(Double_t pt) const {
4edfaee0	464	//
650b3ced	465	//give the pt bin where the pt lies.
4edfaee0	466	//
eaa73e44	467	Int_t ptbin=-1;
650b3ced	468	for (Int_t i=0;i<fnPtBins;i++){
	469	if(pt<fPtBinLimits[i+1]) {
	470	ptbin=i;
	471	break;
	472	}
	473	}
	474	return ptbin;
	475	}
aaf4e30f	476	//-------------------------------------------------------------------
	477	Float_t AliRDHFCuts::GetCutValue(Int_t iVar,Int_t iPtBin) const {
	478	//
	479	// Give the value of cut set for the variable iVar and the pt bin iPtBin
	480	//
	481	if(!fCutsRD){
	482	cout<<"Cuts not iniziaisez yet"<<endl;
	483	return 0;
	484	}
	485	return fCutsRD[GetGlobalIndex(iVar,iPtBin)];
	486	}
0198f78c	487	//-------------------------------------------------------------------
	488	Bool_t AliRDHFCuts::CompareCuts(const AliRDHFCuts *obj) const {
	489	//
	490	// Compare two cuts objects
	491	//
	492
	493	Bool_t areEqual=kTRUE;
	494
	495	if(fMinVtxType!=obj->fMinVtxType) { printf("Minimum vtx type %d %d\n",fMinVtxType,obj->fMinVtxType); areEqual=kFALSE;}
	496
	497	if(fMinVtxContr!=obj->fMinVtxContr) { printf("Minimum vtx contr %d %d\n",fMinVtxContr,obj->fMinVtxContr); areEqual=kFALSE;}
	498
	499	if(TMath::Abs(fMaxVtxRedChi2-obj->fMaxVtxRedChi2)>1.e-10) { printf("Max vtx red chi2 %f %f\n",fMaxVtxRedChi2,obj->fMaxVtxRedChi2);areEqual=kFALSE;}
	500
	501	if(fMinSPDMultiplicity!=obj->fMinSPDMultiplicity) { printf("Min SPD mult %d\n %d",fMinSPDMultiplicity,obj->fMinSPDMultiplicity);areEqual=kFALSE;}
	502
	503	if(fUsePID!=obj->fUsePID) { printf("Use PID %d %d\n",(Int_t)fUsePID,(Int_t)obj->fUsePID); areEqual=kFALSE;}
aaf4e30f	504
0198f78c	505	if(fRemoveDaughtersFromPrimary!=obj->fRemoveDaughtersFromPrimary) {printf("Remove daughters from vtx %d %d\n",(Int_t)fRemoveDaughtersFromPrimary,(Int_t)obj->fRemoveDaughtersFromPrimary); areEqual=kFALSE;}
e6ee9709	506	if(fTrackCuts){
e6ee9709	507	if(fTrackCuts->GetMinNClusterTPC()!=obj->fTrackCuts->GetMinNClusterTPC()) {printf("MinNClsTPC %d %d\n",fTrackCuts->GetMinNClusterTPC(),obj->fTrackCuts->GetMinNClusterTPC()); areEqual=kFALSE;}
0198f78c	508
e6ee9709	509	if(fTrackCuts->GetMinNClustersITS()!=obj->fTrackCuts->GetMinNClustersITS()) {printf("MinNClsITS %d %d\n",fTrackCuts->GetMinNClustersITS(),obj->fTrackCuts->GetMinNClustersITS()); areEqual=kFALSE;}
0198f78c	510
e6ee9709	511	if(TMath::Abs(fTrackCuts->GetMaxChi2PerClusterTPC()-obj->fTrackCuts->GetMaxChi2PerClusterTPC())>1.e-10) {printf("MaxChi2ClsTPC %f %f\n",fTrackCuts->GetMaxChi2PerClusterTPC(),obj->fTrackCuts->GetMaxChi2PerClusterTPC()); areEqual=kFALSE;}
0198f78c	512
e6ee9709	513	if(fTrackCuts->GetClusterRequirementITS(AliESDtrackCuts::kSPD)!=obj->fTrackCuts->GetClusterRequirementITS(AliESDtrackCuts::kSPD)) {printf("ClusterReq SPD %d %d\n",fTrackCuts->GetClusterRequirementITS(AliESDtrackCuts::kSPD),obj->fTrackCuts->GetClusterRequirementITS(AliESDtrackCuts::kSPD)); areEqual=kFALSE;}
e6ee9709	514	}
0198f78c	515
	516	if(fCutsRD) {
	517	for(Int_t iv=0;iv<fnVars;iv++) {
	518	for(Int_t ib=0;ib<fnPtBins;ib++) {
	519	if(TMath::Abs(fCutsRD[GetGlobalIndex(iv,ib)]-obj->fCutsRD[GetGlobalIndex(iv,ib)])>1.e-10) {
	520	cout<<"fCutsRD["<<iv<<"]["<<ib<<"] = "<<fCutsRD[GetGlobalIndex(iv,ib)]<<" "<<obj->fCutsRD[GetGlobalIndex(iv,ib)]<<"\n";
	521	areEqual=kFALSE;
	522	}
	523	}
	524	}
	525	}
	526
	527	return areEqual;
	528	}