[u/mrichter/AliRoot.git] / PWGHF / hfe / AliHFEpid.cxx

/**************************************************************************
* Copyright(c) 1998-1999, 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.                  *
**************************************************************************/
//
// PID Steering Class 
// Interface to the user task
// Several strategies for Electron identification implemented.
// In addition users can combine different detectors or use
// single detector PID
// 
// Authors:
//   Markus Fasel <M.Fasel@gsi.de>
//
#include <TAxis.h>
#include <TClass.h>
#include <TF1.h>
#include <TIterator.h>
#include <TList.h>
#include <TMath.h>
#include <TObjArray.h>
#include <TObjString.h>
#include <TString.h>

#include "AliLog.h"
#include "AliPID.h"
#include "AliVParticle.h"

#include "AliHFEcontainer.h"
#include "AliHFEpid.h"
#include "AliHFEpidQAmanager.h"
#include "AliHFEpidITS.h"
#include "AliHFEpidTPC.h"
#include "AliHFEpidTRD.h"
#include "AliHFEpidTOF.h"
#include "AliHFEpidEMCAL.h"
#include "AliHFEpidMC.h"
#include "AliHFEpidBayes.h"
#include "AliHFEvarManager.h"

ClassImp(AliHFEpid)

const Char_t* AliHFEpid::fgkDetectorName[AliHFEpid::kNdetectorPID + 1] = {
  "MCPID",
  "BAYESPID",
  "ITSPID",
  "TPCPID",
  "TRDPID",
  "TOFPID",
  "EMCALPID",
  "UndefinedPID"
};

//____________________________________________________________
AliHFEpid::AliHFEpid():
  TNamed(),
  fEnabledDetectors(0),
  fNPIDdetectors(0),
  fVarManager(NULL),
  fCommonObjects(NULL)
{
  //
  // Default constructor
  //
  memset(fDetectorPID, 0, sizeof(AliHFEpidBase *) * kNdetectorPID);
  memset(fDetectorOrder, kUndefined, sizeof(UInt_t) * kNdetectorPID);
  memset(fSortedOrder, 0, sizeof(UInt_t) * kNdetectorPID);
}

//____________________________________________________________
AliHFEpid::AliHFEpid(const Char_t *name):
  TNamed(name, ""),
  fEnabledDetectors(0),
  fNPIDdetectors(0),
  fVarManager(NULL),
  fCommonObjects(NULL)
{
  //
  // Default constructor
  // Create PID objects for all detectors
  //
  memset(fDetectorPID, 0, sizeof(AliHFEpidBase *) * kNdetectorPID);
  memset(fDetectorOrder, kUndefined, sizeof(UInt_t) * kNdetectorPID);
  memset(fSortedOrder, 0, sizeof(UInt_t) * kNdetectorPID);

  fDetectorPID[kMCpid] = new AliHFEpidMC("MCPID");
  fDetectorPID[kBAYESpid] = new AliHFEpidBayes("BAYESPID");
  fDetectorPID[kTPCpid] = new AliHFEpidTPC("TPCPID");
  fDetectorPID[kTRDpid] = new AliHFEpidTRD("TRDPID");
  fDetectorPID[kTOFpid] = new AliHFEpidTOF("TOFPID");
  fDetectorPID[kEMCALpid] = new AliHFEpidEMCAL("EMCALPID");

}

//____________________________________________________________
AliHFEpid::AliHFEpid(const AliHFEpid &c):
  TNamed(c),
  fEnabledDetectors(c.fEnabledDetectors),
  fNPIDdetectors(c.fNPIDdetectors),
  fVarManager(c.fVarManager),
  fCommonObjects(NULL)
{
  //
  // Copy Constructor
  //
  memset(fDetectorPID, 0, sizeof(AliHFEpidBase *) * kNdetectorPID);
  c.Copy(*this);
}

//____________________________________________________________
AliHFEpid& AliHFEpid::operator=(const AliHFEpid &c){
  //
  // Assignment operator
  //
  if(&c != this) c.Copy(*this);
  return *this;
}

//____________________________________________________________
AliHFEpid::~AliHFEpid(){
  //
  // Destructor
  //
  for(Int_t idet = 0; idet < kNdetectorPID; idet++)
    if(fDetectorPID[idet]) delete fDetectorPID[idet];
  ClearCommonObjects();
}

//____________________________________________________________
void AliHFEpid::Copy(TObject &o) const{
  //
  // Make copy
  //
  
  TNamed::Copy(o);
  AliHFEpid &target = dynamic_cast<AliHFEpid &>(o);
  target.ClearCommonObjects();

  target.fEnabledDetectors = fEnabledDetectors;
  target.fNPIDdetectors = fNPIDdetectors;
  target.fVarManager = fVarManager;
 
  // Copy detector PIDs
  for(Int_t idet = 0; idet < kNdetectorPID; idet++){
    //Cleanup pointers in case of assignment
    if(target.fDetectorPID[idet])  
      delete target.fDetectorPID[idet];     
    if(fDetectorPID[idet]) 
      target.fDetectorPID[idet] = dynamic_cast<AliHFEpidBase *>(fDetectorPID[idet]->Clone());
  }
  memcpy(target.fDetectorOrder, fDetectorOrder, sizeof(UInt_t) * kNdetectorPID);
  memcpy(target.fSortedOrder, fSortedOrder, sizeof(UInt_t) * kNdetectorPID);
}

//____________________________________________________________
void AliHFEpid::AddCommonObject(TObject * const o){
  //
  // Add common object to the garbage collection
  //
  if(!fCommonObjects) fCommonObjects = new TObjArray;
  fCommonObjects->Add(o);
}

//____________________________________________________________
void AliHFEpid::ClearCommonObjects(){
  //
  // empty garbage collection
  //
  if(fCommonObjects){
    fCommonObjects->Delete();
    delete fCommonObjects;
    fCommonObjects = NULL;
  }
}

//____________________________________________________________
void AliHFEpid::AddDetector(TString detector, UInt_t position){
  //
  // Add Detector in position 
  //
  UInt_t detectorID = kUndefined;
  detector.ToUpper();
  if(!detector.CompareTo("MC")) detectorID = kMCpid;
  else if(!detector.CompareTo("BAYES")) detectorID = kBAYESpid;
  else if(!detector.CompareTo("TPC")) detectorID = kTPCpid;
  else if(!detector.CompareTo("TRD")) detectorID = kTRDpid;
  else if(!detector.CompareTo("TOF")) detectorID = kTOFpid;
  else if(!detector.CompareTo("EMCAL")) detectorID = kEMCALpid;
  else AliError("Detector not available");

  if(detectorID == kUndefined) return;
  if(IsDetectorOn(detectorID)) return;
  SwitchOnDetector(detectorID);
  fDetectorOrder[detectorID] = position;
  fNPIDdetectors++;
}

//____________________________________________________________
Bool_t AliHFEpid::InitializePID(Int_t run){
  //
  // Initializes the PID object
  //

  if(!TestBit(kDetectorsSorted)) SortDetectors();
  Bool_t status = kTRUE;
  for(Int_t idet = 0; idet < kNdetectorPID; idet++){
    if(!IsDetectorOn(idet)) continue;
    if(fDetectorPID[idet]){ 
      status &= fDetectorPID[idet]->InitializePID(run);
      if(HasMCData() && status) fDetectorPID[idet]->SetHasMCData();
    }
  }
  SetBit(kIsInit);
  PrintStatus();
  return status;
}

//____________________________________________________________
Bool_t AliHFEpid::IsSelected(const AliHFEpidObject * const track, AliHFEcontainer *cont, const Char_t *contname, AliHFEpidQAmanager *pidqa){
  //
  // Select Tracks
  //
  Bool_t isSelected = kTRUE;
  AliDebug(1, Form("Particle used for PID, QA available: %s", pidqa ? "Yes" : "No"));
  for(UInt_t idet = 0; idet < fNPIDdetectors; idet++){
    AliDebug(2, Form("Using Detector %s\n", SortedDetectorName(idet)));
    if(TMath::Abs(fDetectorPID[fSortedOrder[idet]]->IsSelected(track, pidqa)) != 11){
      isSelected = kFALSE;
      break;
    }
    AliDebug(2, "Particlae selected by detector");
    if(fVarManager && cont){
      TString reccontname = contname; reccontname += "Reco";
      AliDebug(2, Form("Filling container %s", reccontname.Data()));
      if(fVarManager->IsSignalTrack())
        fVarManager->FillContainerStepname(cont, reccontname.Data(), SortedDetectorName(idet));
      if(HasMCData()){
        TString mccontname = contname; mccontname += "MC";
        AliDebug(2, Form("MC Information available, Filling container %s", mccontname.Data()));
        if(fVarManager->IsSignalTrack()) {
          fVarManager->FillContainerStepname(cont, mccontname.Data(), SortedDetectorName(idet), kTRUE);
	        if(cont->GetCorrelationMatrix("correlationstepafterTOF")){
	          TString tstept("TOFPID"); 
	          if(!tstept.CompareTo(SortedDetectorName(idet))) {
	            fVarManager->FillCorrelationMatrix(cont->GetCorrelationMatrix("correlationstepafterTOF"));
	            //printf("Step %s\n",(const char*) SortedDetectorName(idet));
	          }
	        }
	      }
      }
      // The PID will NOT fill the double counting information
    }
  }
  return isSelected;
}

//____________________________________________________________
void AliHFEpid::SortDetectors(){
  //
  // Make sorted list of detectors
  //
  if(TestBit(kDetectorsSorted)) return; // Don't sort detectors when they are already sorted
  TMath::Sort(static_cast<UInt_t>(kNdetectorPID), fDetectorOrder, fSortedOrder, kFALSE);
  SetBit(kDetectorsSorted);
}

//____________________________________________________________
void AliHFEpid::SetPIDResponse(const AliPIDResponse * const pid){
  //
  // Set ESD PID to the Detector PID objects
  //
  for(Int_t idet = 0; idet < kNdetectorPID; idet++){
    if(fDetectorPID[idet]) fDetectorPID[idet]->SetPIDResponse(pid);
  }
}

//____________________________________________________________
const AliPIDResponse *AliHFEpid::GetPIDResponse() const {
  //
  // Return PID response function
  //
  const AliPIDResponse *response = NULL;
  for(Int_t idet = 0; idet < kNdetectorPID; idet++){
    if(fDetectorPID[idet]){
      response = fDetectorPID[idet]->GetPIDResponse();
      break;
    }
  } 
  return response;
}

//____________________________________________________________
void AliHFEpid::ConfigureTPCasymmetric(Double_t pmin, Double_t pmax, Double_t sigmamin, Double_t sigmamax){
  //
  // TPC alone, symmetric 3 sigma cut and asymmetric sigma cut in the momentum region between 2GeV/c and 10 GeV/c and sigma between -1 and 100
  //
  AliHFEpidTPC *pid = dynamic_cast<AliHFEpidTPC *>(fDetectorPID[kTPCpid]);
  if(pid){
    pid->SetTPCnSigma(3);
    pid->SetAsymmetricTPCsigmaCut(pmin, pmax, sigmamin, sigmamax);
  }
}

//____________________________________________________________
void AliHFEpid::ConfigureTPCrejectionSimple(){
  //
  // TPC alone, symmetric 3 sigma cut and 2 - -100 sigma pion rejection
  //   
  AliHFEpidTPC *pid = dynamic_cast<AliHFEpidTPC *>(fDetectorPID[kTPCpid]);
  if(pid){
    pid->SetTPCnSigma(3);
    pid->SetRejectParticle(AliPID::kPion, 0., -100., 10., 1.);
  }
}

//____________________________________________________________
void AliHFEpid::ConfigureTOF(Float_t TOFCut){
  //
  // Set Number of sigmas for TOF PID
  //
  AliHFEpidTOF *tofpid = dynamic_cast<AliHFEpidTOF *>(fDetectorPID[kTOFpid]);
  if(tofpid) tofpid->SetTOFnSigma(TOFCut);
}

//____________________________________________________________
void AliHFEpid::ConfigureTPCcentralityCut(Int_t centralityBin, const char *lowerCutParam, const Double_t * const params, Float_t upperTPCCut){
  //
  // Cofigure centrality dependent cut function for TPC PID 
  //
  ConfigureTPCcut(centralityBin, lowerCutParam, params, upperTPCCut);
}

//____________________________________________________________
void AliHFEpid::ConfigureTPCdefaultCut(const char *lowerCutParam, const Double_t * const params, Float_t upperTPCCut){
  //
  // Cofigure default cut function for TPC PID 
  //
  ConfigureTPCcut(-1, lowerCutParam, params, upperTPCCut);
}

//____________________________________________________________
void AliHFEpid::ConfigureTPCcut(Int_t centralityBin, const char *lowerCutParam, const Double_t * const params, Float_t upperTPCCut){
  //
  // Cofigure cut function for TPC PID 
  // if no function parameterizaion is given, then the default one (exponential) is chosen
  //
  
  if(HasMCData()) AliInfo("Configuring TPC for MC\n");
  AliHFEpidTPC *tpcpid = dynamic_cast<AliHFEpidTPC *>(fDetectorPID[kTPCpid]);
  //TF1 *upperCut = new TF1("upperCut", "[0] * TMath::Exp([1]*x)", 0, 20);
  TF1 *upperCut = new TF1(Form("upperCut%s", centralityBin <  0 ? "Default" : Form("Bin%d", centralityBin)), "[0]", 0, 20); // Use constant upper cut
  TF1 *lowerCut = new TF1(Form("lowerCut%s", centralityBin <  0 ? "Default" : Form("Bin%d", centralityBin)), lowerCutParam == NULL ? "[0] * TMath::Exp([1]*x) + [2]": lowerCutParam, 0, 20);

  upperCut->SetParameter(0, upperTPCCut); // pp

  if(params){
      for(Int_t ipar = 0; ipar < lowerCut->GetNpar(); ipar++)
      {
	  lowerCut->SetParameter(ipar, params[ipar]);
        //  printf("printout %i %s %f \n", centralityBin, lowerCutParam, params[ipar]);
      }
  } else {
    // Set default parameterization
    if(HasMCData()) lowerCut->SetParameter(0, -2.5);
    else lowerCut->SetParameter(0, -4.03);  //pp
    lowerCut->SetParameter(1, -0.22); // pp
    
    if(HasMCData()) lowerCut->SetParameter(2, -2.2);
    else lowerCut->SetParameter(2, 0.92); //pp
  }


  if(tpcpid){
    tpcpid->SetTPCnSigma(2);
    if(centralityBin < 0){
      tpcpid->SetUpperSigmaCutDefault(upperCut);
      tpcpid->SetLowerSigmaCutDefault(lowerCut);
    } else {
      tpcpid->SetUpperSigmaCutCentrality(upperCut, centralityBin);
      tpcpid->SetLowerSigmaCutCentrality(lowerCut, centralityBin);
    }
  }
  AddCommonObject(upperCut);
  AddCommonObject(lowerCut);
}

//____________________________________________________________
void AliHFEpid::ConfigureBayesDetectorMask(Int_t detmask){
  //
  // Configure detector mask for Bayes PID
  // if no detector mask is set the default mask is chosen
  //
  
  if(HasMCData()) AliInfo("Configuring Bayes for MC\n");
  AliHFEpidBayes *bayespid = dynamic_cast<AliHFEpidBayes *>(fDetectorPID[kBAYESpid]);

  if(bayespid)
  {
      bayespid->SetBayesDetectorMask(detmask);
      printf("detector mask in pid class %i \n",detmask);
  }

}

//____________________________________________________________
void AliHFEpid::ConfigureBayesPIDThreshold(Float_t pidthres){
  //
  // Configure pid threshold for Bayes PID
  // if no threshold is set the default threshold is chosen
  //
  
  if(HasMCData()) AliInfo("Configuring Bayes for MC\n");
  AliHFEpidBayes *bayespid = dynamic_cast<AliHFEpidBayes *>(fDetectorPID[kBAYESpid]);

  if(bayespid)
  {
      bayespid->SetBayesPIDThreshold(pidthres);
      printf("combined pidthreshold %f \n",pidthres);
  }

}

//____________________________________________________________
void AliHFEpid::PrintStatus() const {
  //
  // Print the PID configuration
  //
  printf("\n%s: Printing configuration\n", GetName());
  printf("===============================================\n");
  printf("PID Detectors: \n");
  Int_t npid = 0;
  TString detectors[kNdetectorPID] = {"MC", "BAYES", "ITS", "TPC", "TRD", "TOF", "EMCAL"};
  for(Int_t idet = 0; idet < kNdetectorPID; idet++){
    if(IsDetectorOn(idet)){
      printf("\t%s\n", detectors[idet].Data());
      npid++;
    }
  }
  if(!npid) printf("\tNone\n");
  printf("\n");
}
Commit	Line	Data
809a4336	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, 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	**************************************************************************/
50685501	15	//
	16	// PID Steering Class
	17	// Interface to the user task
9bcfd1ab	18	// Several strategies for Electron identification implemented.
	19	// In addition users can combine different detectors or use
	20	// single detector PID
50685501	21	//
	22	// Authors:
	23	// Markus Fasel <M.Fasel@gsi.de>
	24	//
	25	#include <TAxis.h>
809a4336	26	#include <TClass.h>
faee3b18	27	#include <TF1.h>
809a4336	28	#include <TIterator.h>
809a4336	29	#include <TList.h>
8c1c76e9	30	#include <TMath.h>
809a4336	31	#include <TObjArray.h>
	32	#include <TObjString.h>
	33	#include <TString.h>
	34
75d81601	35	#include "AliLog.h"
75d81601	36	#include "AliPID.h"
3a72645a	37	#include "AliVParticle.h"
809a4336	38
3a72645a	39	#include "AliHFEcontainer.h"
809a4336	40	#include "AliHFEpid.h"
3a72645a	41	#include "AliHFEpidQAmanager.h"
75d81601	42	#include "AliHFEpidITS.h"
809a4336	43	#include "AliHFEpidTPC.h"
	44	#include "AliHFEpidTRD.h"
	45	#include "AliHFEpidTOF.h"
c2690925	46	#include "AliHFEpidEMCAL.h"
809a4336	47	#include "AliHFEpidMC.h"
8c1c76e9	48	#include "AliHFEpidBayes.h"
3a72645a	49	#include "AliHFEvarManager.h"
809a4336	50
	51	ClassImp(AliHFEpid)
	52
3a72645a	53	const Char_t* AliHFEpid::fgkDetectorName[AliHFEpid::kNdetectorPID + 1] = {
3a72645a	54	"MCPID",
8c1c76e9	55	"BAYESPID",
3a72645a	56	"ITSPID",
	57	"TPCPID",
	58	"TRDPID",
	59	"TOFPID",
e3ae862b	60	"EMCALPID",
3a72645a	61	"UndefinedPID"
	62	};
	63
809a4336	64	//____________________________________________________________
809a4336	65	AliHFEpid::AliHFEpid():
e3fc062d	66	TNamed(),
809a4336	67	fEnabledDetectors(0),
3a72645a	68	fNPIDdetectors(0),
3a72645a	69	fVarManager(NULL),
faee3b18	70	fCommonObjects(NULL)
809a4336	71	{
	72	//
	73	// Default constructor
	74	//
	75	memset(fDetectorPID, 0, sizeof(AliHFEpidBase ) kNdetectorPID);
3a72645a	76	memset(fDetectorOrder, kUndefined, sizeof(UInt_t) * kNdetectorPID);
8c1c76e9	77	memset(fSortedOrder, 0, sizeof(UInt_t) * kNdetectorPID);
809a4336	78	}
809a4336	79
e3fc062d	80	//____________________________________________________________
	81	AliHFEpid::AliHFEpid(const Char_t *name):
	82	TNamed(name, ""),
	83	fEnabledDetectors(0),
3a72645a	84	fNPIDdetectors(0),
3a72645a	85	fVarManager(NULL),
e3fc062d	86	fCommonObjects(NULL)
	87	{
	88	//
	89	// Default constructor
	90	// Create PID objects for all detectors
	91	//
	92	memset(fDetectorPID, 0, sizeof(AliHFEpidBase ) kNdetectorPID);
3a72645a	93	memset(fDetectorOrder, kUndefined, sizeof(UInt_t) * kNdetectorPID);
	94	memset(fSortedOrder, 0, sizeof(UInt_t) * kNdetectorPID);
	95
e3fc062d	96	fDetectorPID[kMCpid] = new AliHFEpidMC("MCPID");
8c1c76e9	97	fDetectorPID[kBAYESpid] = new AliHFEpidBayes("BAYESPID");
3a72645a	98	fDetectorPID[kTPCpid] = new AliHFEpidTPC("TPCPID");
e3fc062d	99	fDetectorPID[kTRDpid] = new AliHFEpidTRD("TRDPID");
e3fc062d	100	fDetectorPID[kTOFpid] = new AliHFEpidTOF("TOFPID");
c2690925	101	fDetectorPID[kEMCALpid] = new AliHFEpidEMCAL("EMCALPID");
e3fc062d	102
	103	}
	104
809a4336	105	//____________________________________________________________
809a4336	106	AliHFEpid::AliHFEpid(const AliHFEpid &c):
e3fc062d	107	TNamed(c),
809a4336	108	fEnabledDetectors(c.fEnabledDetectors),
3a72645a	109	fNPIDdetectors(c.fNPIDdetectors),
3a72645a	110	fVarManager(c.fVarManager),
faee3b18	111	fCommonObjects(NULL)
809a4336	112	{
	113	//
	114	// Copy Constructor
	115	//
	116	memset(fDetectorPID, 0, sizeof(AliHFEpidBase ) kNdetectorPID);
e3fc062d	117	c.Copy(*this);
809a4336	118	}
	119
	120	//____________________________________________________________
	121	AliHFEpid& AliHFEpid::operator=(const AliHFEpid &c){
	122	//
	123	// Assignment operator
	124	//
e3fc062d	125	if(&c != this) c.Copy(*this);
e3fc062d	126	return *this;
809a4336	127	}
	128
	129	//____________________________________________________________
	130	AliHFEpid::~AliHFEpid(){
	131	//
	132	// Destructor
	133	//
e3fc062d	134	for(Int_t idet = 0; idet < kNdetectorPID; idet++)
e3fc062d	135	if(fDetectorPID[idet]) delete fDetectorPID[idet];
faee3b18	136	ClearCommonObjects();
	137	}
	138
e3fc062d	139	//____________________________________________________________
	140	void AliHFEpid::Copy(TObject &o) const{
	141	//
	142	// Make copy
	143	//
	144
	145	TNamed::Copy(o);
	146	AliHFEpid &target = dynamic_cast<AliHFEpid &>(o);
	147	target.ClearCommonObjects();
	148
	149	target.fEnabledDetectors = fEnabledDetectors;
3a72645a	150	target.fNPIDdetectors = fNPIDdetectors;
3a72645a	151	target.fVarManager = fVarManager;
e3fc062d	152
	153	// Copy detector PIDs
	154	for(Int_t idet = 0; idet < kNdetectorPID; idet++){
	155	//Cleanup pointers in case of assignment
	156	if(target.fDetectorPID[idet])
	157	delete target.fDetectorPID[idet];
	158	if(fDetectorPID[idet])
	159	target.fDetectorPID[idet] = dynamic_cast<AliHFEpidBase *>(fDetectorPID[idet]->Clone());
	160	}
3a72645a	161	memcpy(target.fDetectorOrder, fDetectorOrder, sizeof(UInt_t) * kNdetectorPID);
3a72645a	162	memcpy(target.fSortedOrder, fSortedOrder, sizeof(UInt_t) * kNdetectorPID);
e3fc062d	163	}
e3fc062d	164
faee3b18	165	//____________________________________________________________
	166	void AliHFEpid::AddCommonObject(TObject * const o){
	167	//
	168	// Add common object to the garbage collection
	169	//
	170	if(!fCommonObjects) fCommonObjects = new TObjArray;
	171	fCommonObjects->Add(o);
	172	}
	173
	174	//____________________________________________________________
	175	void AliHFEpid::ClearCommonObjects(){
	176	//
	177	// empty garbage collection
	178	//
	179	if(fCommonObjects){
	180	fCommonObjects->Delete();
	181	delete fCommonObjects;
	182	fCommonObjects = NULL;
	183	}
809a4336	184	}
	185
	186	//____________________________________________________________
3a72645a	187	void AliHFEpid::AddDetector(TString detector, UInt_t position){
809a4336	188	//
3a72645a	189	// Add Detector in position
809a4336	190	//
3a72645a	191	UInt_t detectorID = kUndefined;
	192	detector.ToUpper();
	193	if(!detector.CompareTo("MC")) detectorID = kMCpid;
8c1c76e9	194	else if(!detector.CompareTo("BAYES")) detectorID = kBAYESpid;
3a72645a	195	else if(!detector.CompareTo("TPC")) detectorID = kTPCpid;
	196	else if(!detector.CompareTo("TRD")) detectorID = kTRDpid;
	197	else if(!detector.CompareTo("TOF")) detectorID = kTOFpid;
e3ae862b	198	else if(!detector.CompareTo("EMCAL")) detectorID = kEMCALpid;
3a72645a	199	else AliError("Detector not available");
e3fc062d	200
3a72645a	201	if(detectorID == kUndefined) return;
	202	if(IsDetectorOn(detectorID)) return;
	203	SwitchOnDetector(detectorID);
	204	fDetectorOrder[detectorID] = position;
	205	fNPIDdetectors++;
	206	}
	207
	208	//____________________________________________________________
8c1c76e9	209	Bool_t AliHFEpid::InitializePID(Int_t run){
3a72645a	210	//
	211	// Initializes the PID object
	212	//
	213
8c1c76e9	214	if(!TestBit(kDetectorsSorted)) SortDetectors();
809a4336	215	Bool_t status = kTRUE;
809a4336	216	for(Int_t idet = 0; idet < kNdetectorPID; idet++){
e3fc062d	217	if(!IsDetectorOn(idet)) continue;
809a4336	218	if(fDetectorPID[idet]){
8c1c76e9	219	status &= fDetectorPID[idet]->InitializePID(run);
809a4336	220	if(HasMCData() && status) fDetectorPID[idet]->SetHasMCData();
	221	}
	222	}
8c1c76e9	223	SetBit(kIsInit);
e3fc062d	224	PrintStatus();
809a4336	225	return status;
	226	}
	227
	228	//____________________________________________________________
e156c3bb	229	Bool_t AliHFEpid::IsSelected(const AliHFEpidObject * const track, AliHFEcontainer cont, const Char_t contname, AliHFEpidQAmanager *pidqa){
809a4336	230	//
3a72645a	231	// Select Tracks
809a4336	232	//
3a72645a	233	Bool_t isSelected = kTRUE;
	234	AliDebug(1, Form("Particle used for PID, QA available: %s", pidqa ? "Yes" : "No"));
	235	for(UInt_t idet = 0; idet < fNPIDdetectors; idet++){
	236	AliDebug(2, Form("Using Detector %s\n", SortedDetectorName(idet)));
	237	if(TMath::Abs(fDetectorPID[fSortedOrder[idet]]->IsSelected(track, pidqa)) != 11){
	238	isSelected = kFALSE;
	239	break;
0792aa82	240	}
3a72645a	241	AliDebug(2, "Particlae selected by detector");
3a72645a	242	if(fVarManager && cont){
bf892a6a	243	TString reccontname = contname; reccontname += "Reco";
bf892a6a	244	AliDebug(2, Form("Filling container %s", reccontname.Data()));
3a72645a	245	if(fVarManager->IsSignalTrack())
bf892a6a	246	fVarManager->FillContainerStepname(cont, reccontname.Data(), SortedDetectorName(idet));
3a72645a	247	if(HasMCData()){
bf892a6a	248	TString mccontname = contname; mccontname += "MC";
	249	AliDebug(2, Form("MC Information available, Filling container %s", mccontname.Data()));
	250	if(fVarManager->IsSignalTrack()) {
	251	fVarManager->FillContainerStepname(cont, mccontname.Data(), SortedDetectorName(idet), kTRUE);
e3ae862b	252	if(cont->GetCorrelationMatrix("correlationstepafterTOF")){
	253	TString tstept("TOFPID");
	254	if(!tstept.CompareTo(SortedDetectorName(idet))) {
	255	fVarManager->FillCorrelationMatrix(cont->GetCorrelationMatrix("correlationstepafterTOF"));
	256	//printf("Step %s\n",(const char*) SortedDetectorName(idet));
	257	}
	258	}
	259	}
3a72645a	260	}
3a72645a	261	// The PID will NOT fill the double counting information
809a4336	262	}
75d81601	263	}
3a72645a	264	return isSelected;
75d81601	265	}
75d81601	266
faee3b18	267	//____________________________________________________________
8c1c76e9	268	void AliHFEpid::SortDetectors(){
	269	//
	270	// Make sorted list of detectors
	271	//
	272	if(TestBit(kDetectorsSorted)) return; // Don't sort detectors when they are already sorted
	273	TMath::Sort(static_cast<UInt_t>(kNdetectorPID), fDetectorOrder, fSortedOrder, kFALSE);
	274	SetBit(kDetectorsSorted);
	275	}
	276
	277	//____________________________________________________________
	278	void AliHFEpid::SetPIDResponse(const AliPIDResponse * const pid){
faee3b18	279	//
	280	// Set ESD PID to the Detector PID objects
	281	//
	282	for(Int_t idet = 0; idet < kNdetectorPID; idet++){
8c1c76e9	283	if(fDetectorPID[idet]) fDetectorPID[idet]->SetPIDResponse(pid);
faee3b18	284	}
	285	}
	286
809a4336	287	//____________________________________________________________
8c1c76e9	288	const AliPIDResponse *AliHFEpid::GetPIDResponse() const {
0792aa82	289	//
8c1c76e9	290	// Return PID response function
0792aa82	291	//
8c1c76e9	292	const AliPIDResponse *response = NULL;
3a72645a	293	for(Int_t idet = 0; idet < kNdetectorPID; idet++){
8c1c76e9	294	if(fDetectorPID[idet]){
	295	response = fDetectorPID[idet]->GetPIDResponse();
	296	break;
	297	}
	298	}
	299	return response;
0792aa82	300	}
	301
	302	//____________________________________________________________
3a72645a	303	void AliHFEpid::ConfigureTPCasymmetric(Double_t pmin, Double_t pmax, Double_t sigmamin, Double_t sigmamax){
0792aa82	304	//
	305	// TPC alone, symmetric 3 sigma cut and asymmetric sigma cut in the momentum region between 2GeV/c and 10 GeV/c and sigma between -1 and 100
	306	//
e3fc062d	307	AliHFEpidTPC pid = dynamic_cast<AliHFEpidTPC >(fDetectorPID[kTPCpid]);
bf892a6a	308	if(pid){
	309	pid->SetTPCnSigma(3);
	310	pid->SetAsymmetricTPCsigmaCut(pmin, pmax, sigmamin, sigmamax);
	311	}
0792aa82	312	}
	313
	314	//____________________________________________________________
3a72645a	315	void AliHFEpid::ConfigureTPCrejectionSimple(){
0792aa82	316	//
	317	// TPC alone, symmetric 3 sigma cut and 2 - -100 sigma pion rejection
	318	//
e3fc062d	319	AliHFEpidTPC pid = dynamic_cast<AliHFEpidTPC >(fDetectorPID[kTPCpid]);
bf892a6a	320	if(pid){
	321	pid->SetTPCnSigma(3);
	322	pid->SetRejectParticle(AliPID::kPion, 0., -100., 10., 1.);
	323	}
0792aa82	324	}
	325
	326	//____________________________________________________________
e156c3bb	327	void AliHFEpid::ConfigureTOF(Float_t TOFCut){
9bcfd1ab	328	//
e156c3bb	329	// Set Number of sigmas for TOF PID
9bcfd1ab	330	//
e3fc062d	331	AliHFEpidTOF tofpid = dynamic_cast<AliHFEpidTOF >(fDetectorPID[kTOFpid]);
c2690925	332	if(tofpid) tofpid->SetTOFnSigma(TOFCut);
e156c3bb	333	}
	334
	335	//____________________________________________________________
	336	void AliHFEpid::ConfigureTPCcentralityCut(Int_t centralityBin, const char lowerCutParam, const Double_t const params, Float_t upperTPCCut){
	337	//
	338	// Cofigure centrality dependent cut function for TPC PID
	339	//
	340	ConfigureTPCcut(centralityBin, lowerCutParam, params, upperTPCCut);
	341	}
3a72645a	342
e156c3bb	343	//____________________________________________________________
	344	void AliHFEpid::ConfigureTPCdefaultCut(const char lowerCutParam, const Double_t const params, Float_t upperTPCCut){
	345	//
	346	// Cofigure default cut function for TPC PID
	347	//
	348	ConfigureTPCcut(-1, lowerCutParam, params, upperTPCCut);
	349	}
	350
	351	//____________________________________________________________
	352	void AliHFEpid::ConfigureTPCcut(Int_t centralityBin, const char lowerCutParam, const Double_t const params, Float_t upperTPCCut){
	353	//
	354	// Cofigure cut function for TPC PID
	355	// if no function parameterizaion is given, then the default one (exponential) is chosen
	356	//
	357
	358	if(HasMCData()) AliInfo("Configuring TPC for MC\n");
	359	AliHFEpidTPC tpcpid = dynamic_cast<AliHFEpidTPC >(fDetectorPID[kTPCpid]);
faee3b18	360	//TF1 upperCut = new TF1("upperCut", "[0] TMath::Exp([1]*x)", 0, 20);
e156c3bb	361	TF1 *upperCut = new TF1(Form("upperCut%s", centralityBin < 0 ? "Default" : Form("Bin%d", centralityBin)), "[0]", 0, 20); // Use constant upper cut
e156c3bb	362	TF1 lowerCut = new TF1(Form("lowerCut%s", centralityBin < 0 ? "Default" : Form("Bin%d", centralityBin)), lowerCutParam == NULL ? "[0] TMath::Exp([1]*x) + [2]": lowerCutParam, 0, 20);
c2690925	363
	364	upperCut->SetParameter(0, upperTPCCut); // pp
	365
e3ae862b	366	if(params){
e156c3bb	367	for(Int_t ipar = 0; ipar < lowerCut->GetNpar(); ipar++)
	368	{
	369	lowerCut->SetParameter(ipar, params[ipar]);
	370	// printf("printout %i %s %f \n", centralityBin, lowerCutParam, params[ipar]);
	371	}
e3ae862b	372	} else {
e3ae862b	373	// Set default parameterization
e156c3bb	374	if(HasMCData()) lowerCut->SetParameter(0, -2.5);
e156c3bb	375	else lowerCut->SetParameter(0, -4.03); //pp
c2690925	376	lowerCut->SetParameter(1, -0.22); // pp
c2690925	377
e3ae862b	378	if(HasMCData()) lowerCut->SetParameter(2, -2.2);
c2690925	379	else lowerCut->SetParameter(2, 0.92); //pp
e3ae862b	380	}
c2690925	381
c2690925	382
bf892a6a	383	if(tpcpid){
bf892a6a	384	tpcpid->SetTPCnSigma(2);
e156c3bb	385	if(centralityBin < 0){
	386	tpcpid->SetUpperSigmaCutDefault(upperCut);
	387	tpcpid->SetLowerSigmaCutDefault(lowerCut);
	388	} else {
	389	tpcpid->SetUpperSigmaCutCentrality(upperCut, centralityBin);
	390	tpcpid->SetLowerSigmaCutCentrality(lowerCut, centralityBin);
	391	}
bf892a6a	392	}
faee3b18	393	AddCommonObject(upperCut);
faee3b18	394	AddCommonObject(lowerCut);
9bcfd1ab	395	}
9bcfd1ab	396
8c1c76e9	397	//____________________________________________________________
	398	void AliHFEpid::ConfigureBayesDetectorMask(Int_t detmask){
	399	//
	400	// Configure detector mask for Bayes PID
	401	// if no detector mask is set the default mask is chosen
	402	//
	403
	404	if(HasMCData()) AliInfo("Configuring Bayes for MC\n");
	405	AliHFEpidBayes bayespid = dynamic_cast<AliHFEpidBayes >(fDetectorPID[kBAYESpid]);
	406
	407	if(bayespid)
	408	{
	409	bayespid->SetBayesDetectorMask(detmask);
	410	printf("detector mask in pid class %i \n",detmask);
	411	}
	412
	413	}
	414
	415	//____________________________________________________________
	416	void AliHFEpid::ConfigureBayesPIDThreshold(Float_t pidthres){
	417	//
	418	// Configure pid threshold for Bayes PID
	419	// if no threshold is set the default threshold is chosen
	420	//
	421
	422	if(HasMCData()) AliInfo("Configuring Bayes for MC\n");
	423	AliHFEpidBayes bayespid = dynamic_cast<AliHFEpidBayes >(fDetectorPID[kBAYESpid]);
	424
	425	if(bayespid)
	426	{
	427	bayespid->SetBayesPIDThreshold(pidthres);
	428	printf("combined pidthreshold %f \n",pidthres);
	429	}
	430
	431	}
	432
e3fc062d	433	//____________________________________________________________
	434	void AliHFEpid::PrintStatus() const {
	435	//
	436	// Print the PID configuration
	437	//
	438	printf("\n%s: Printing configuration\n", GetName());
	439	printf("===============================================\n");
e3fc062d	440	printf("PID Detectors: \n");
e3fc062d	441	Int_t npid = 0;
8c1c76e9	442	TString detectors[kNdetectorPID] = {"MC", "BAYES", "ITS", "TPC", "TRD", "TOF", "EMCAL"};
e3fc062d	443	for(Int_t idet = 0; idet < kNdetectorPID; idet++){
	444	if(IsDetectorOn(idet)){
	445	printf("\t%s\n", detectors[idet].Data());
	446	npid++;
	447	}
	448	}
	449	if(!npid) printf("\tNone\n");
	450	printf("\n");
	451	}