#include <AliESDInputHandler.h>
#include <AliAnalysisManager.h>
#include <AliTrackerBase.h>
-#include <AliTRDPIDReference.h>
+#include <AliTRDPIDResponseObject.h>
#include <AliTRDPIDResponse.h>
#include "AliTRDCalChamberStatus.h"
#include <AliTender.h>
// Get new references
TIter refs(arr);
TObject *o = NULL;
- AliTRDPIDReference *ref = NULL;
+ AliTRDPIDResponseObject *ref = NULL;
while((o = refs())){
- if(!TString(o->IsA()->GetName()).CompareTo("AliTRDPIDReference")){
- ref = dynamic_cast<AliTRDPIDReference *>(o);
+ if(!TString(o->IsA()->GetName()).CompareTo("AliTRDPIDResponseObject")){
+ ref = dynamic_cast<AliTRDPIDResponseObject *>(o);
break;
}
}
if(ref){
- fESDpid->GetTRDResponse().Load(ref);
+ fESDpid->GetTRDResponse().SetPIDResponseObject(ref);
fHasReferences = kTRUE;
AliDebug(1, "Reference distributions loaded into the PID Response");
} else {
STEERBase/AliITSPidParams.cxx
STEERBase/AliTRDPIDReference.cxx
STEERBase/AliTRDPIDParams.cxx
- STEERBase/AliTOFPIDParams.cxx
+ STEERBase/AliTRDPIDResponseObject.cxx
+ STEERBase/AliTRDTKDInterpolator.cxx
STEERBase/AliPIDResponse.cxx
STEERBase/AliITSPIDResponse.cxx
STEERBase/AliTPCPIDResponse.cxx
STEERBase/AliTOFPIDResponse.cxx
+ STEERBase/AliTOFPIDParams.cxx
STEERBase/AliTRDPIDResponse.cxx
STEERBase/AliEMCALPIDResponse.cxx
STEERBase/AliPIDCombined.cxx
#include <TF1.h>
#include <TSpline.h>
#include <TFile.h>
+#include <TArrayF.h>
#include <AliVEvent.h>
#include <AliVTrack.h>
#include <AliLog.h>
#include <AliPID.h>
#include <AliOADBContainer.h>
-#include <AliTRDPIDParams.h>
-#include <AliTRDPIDReference.h>
+#include <AliTRDPIDResponseObject.h>
#include <AliTOFPIDParams.h>
#include "AliPIDResponse.h"
fOldRun(0),
fArrPidResponseMaster(0x0),
fResolutionCorrection(0x0),
-fTRDPIDParams(0x0),
-fTRDPIDReference(0x0),
+fTRDPIDResponseObject(0x0),
fTOFtail(1.1),
fTOFPIDParams(0x0),
fEMCALPIDParams(0x0),
//
// dtor
//
- delete fArrPidResponseMaster;
- delete fTRDPIDParams;
- delete fTRDPIDReference;
+ delete fArrPidResponseMaster;
+ delete fTRDPIDResponseObject;
if (fTOFPIDParams) delete fTOFPIDParams;
}
fOldRun(0),
fArrPidResponseMaster(0x0),
fResolutionCorrection(0x0),
-fTRDPIDParams(0x0),
-fTRDPIDReference(0x0),
+fTRDPIDResponseObject(0x0),
fTOFtail(1.1),
fTOFPIDParams(0x0),
fEMCALPIDParams(0x0),
fOldRun=0;
fArrPidResponseMaster=0x0;
fResolutionCorrection=0x0;
- fTRDPIDParams=0x0;
- fTRDPIDReference=0x0;
+ fTRDPIDResponseObject=0x0;
fEMCALPIDParams=0x0;
memset(fTRDslicesForPID,0,sizeof(UInt_t)*2);
fTOFtail=1.1;
//
// Load the TRD pid params and references from the OADB
//
- if(fTRDPIDParams) return;
+ if(fTRDPIDResponseObject) return;
AliOADBContainer contParams("contParams");
- Int_t statusPars = contParams.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDParams.root", fOADBPath.Data()), "AliTRDPIDParams");
- if(statusPars){
- AliError("Failed initializing PID Params from OADB");
+ Int_t statusResponse = contParams.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDResponse.root", fOADBPath.Data()), "AliTRDPIDResponseObject");
+ if(statusResponse){
+ AliError("Failed initializing PID Response Object from OADB");
} else {
- AliInfo(Form("Loading TRD Params from %s/COMMON/PID/data/TRDPIDParams.root", fOADBPath.Data()));
- fTRDPIDParams = dynamic_cast<AliTRDPIDParams *>(contParams.GetObject(fRun));
- if(!fTRDPIDParams){
- AliError(Form("TRD Params not found in run %d", fRun));
+ AliInfo(Form("Loading TRD Response from %s/COMMON/PID/data/TRDPIDResponse.root", fOADBPath.Data()));
+ fTRDPIDResponseObject = dynamic_cast<AliTRDPIDResponseObject *>(contParams.GetObject(fRun));
+ if(!fTRDPIDResponseObject){
+ AliError(Form("TRD Response not found in run %d", fRun));
}
}
-
+ /*
AliOADBContainer contRefs("contRefs");
Int_t statusRefs = contRefs.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDReferenceLQ1D.root", fOADBPath.Data()), "AliTRDPIDReference");
if(statusRefs){
if(!fTRDPIDReference){
AliError(Form("TRD References not found in OADB Container for run %d", fRun));
}
- }
+ }
+ */
}
//______________________________________________________________________________
//
// Set PID Params and references to the TRD PID response
//
- fTRDResponse.SetPIDParams(fTRDPIDParams);
- fTRDResponse.Load(fTRDPIDReference);
- if(fLHCperiod == "LHC10b" || fLHCperiod == "LHC10c" || fLHCperiod == "LHC10d" || fLHCperiod == "LHC10e"){
- fTRDslicesForPID[0] = 0;
- fTRDslicesForPID[1] = 7;
- }
+ fTRDResponse.SetPIDResponseObject(fTRDPIDResponseObject);
+ if(fLHCperiod == "LHC10b" || fLHCperiod == "LHC10c" || fLHCperiod == "LHC10d" || fLHCperiod == "LHC10e"){
+ fTRDslicesForPID[0] = 0;
+ fTRDslicesForPID[1] = 7;
+ }
}
//______________________________________________________________________________
class AliVEvent;
class TF1;
+class AliTRDPIDResponseObject;
class AliPIDResponse : public TNamed {
public:
TObjArray *fArrPidResponseMaster; //! TPC pid splines
TF1 *fResolutionCorrection; //! TPC resolution correction
- AliTRDPIDParams *fTRDPIDParams; //! TRD PID Params
- AliTRDPIDReference *fTRDPIDReference; //! TRD PID References
+ AliTRDPIDResponseObject *fTRDPIDResponseObject; //! TRD PID Response Object
UInt_t fTRDslicesForPID[2]; //! TRD PID slices
Float_t fTOFtail; //! TOF tail effect used in TOF probability
fEntries = new TSortedList;
}
+//____________________________________________________________
+AliTRDPIDParams::AliTRDPIDParams(const AliTRDPIDParams &ref):
+TNamed(ref),
+fEntries(NULL)
+{
+ //
+ // Copy constructor
+ //
+
+ fEntries=(TSortedList*)ref.fEntries->Clone();
+}
+
+//____________________________________________________________
+AliTRDPIDParams::AliTRDPIDParams(const AliTRDPIDParams &ref):
+TNamed(ref),
+fEntries(NULL)
+{
+ //
+ // Copy constructor
+ //
+
+ fEntries=(TSortedList*)ref.fEntries->Clone();
+}
+
//____________________________________________________________
AliTRDPIDParams::~AliTRDPIDParams(){
//
public:
AliTRDPIDParams();
AliTRDPIDParams(const char *name);
+ AliTRDPIDParams(const AliTRDPIDParams &);
virtual ~AliTRDPIDParams();
virtual void Print(Option_t *) const;
private:
class AliTRDPIDThresholds : public TObject{
- public:
- AliTRDPIDThresholds();
- AliTRDPIDThresholds(Int_t nTracklets, Double_t effMin, Double_t effMax, Double_t *params = NULL);
- AliTRDPIDThresholds(Int_t nTracklets, Double_t eff, Double_t *params = NULL);
- AliTRDPIDThresholds(const AliTRDPIDThresholds &);
- AliTRDPIDThresholds &operator=(const AliTRDPIDThresholds &);
- virtual ~AliTRDPIDThresholds() {}
-
+ public:
+ AliTRDPIDThresholds();
+ AliTRDPIDThresholds(Int_t nTracklets, Double_t effMin, Double_t effMax, Double_t *params = NULL);
+ AliTRDPIDThresholds(Int_t nTracklets, Double_t eff, Double_t *params = NULL);
+ AliTRDPIDThresholds(const AliTRDPIDThresholds &);
+ AliTRDPIDThresholds &operator=(const AliTRDPIDThresholds &);
+ virtual ~AliTRDPIDThresholds() {}
+
Int_t GetNTracklets() const { return fNTracklets; }
- Double_t GetElectronEfficiency(Int_t step = 0) const { if(step == 0) return fEfficiency[0]; else return fEfficiency[1]; }
- const Double_t *GetThresholdParams() const { return fParams; }
+ Double_t GetElectronEfficiency(Int_t step = 0) const { if(step == 0) return fEfficiency[0]; else return fEfficiency[1]; }
+ const Double_t *GetThresholdParams() const { return fParams; }
- virtual Bool_t IsSortable() const { return kTRUE; }
- virtual Int_t Compare(const TObject *ref) const;
+ virtual Bool_t IsSortable() const { return kTRUE; }
+ virtual Int_t Compare(const TObject *ref) const;
virtual Bool_t IsEqual(const TObject *ref) const;
- private:
- Int_t fNTracklets; //
- Double_t fEfficiency[2]; //
- Double_t fParams[4]; //
+ private:
+ Int_t fNTracklets; //
+ Double_t fEfficiency[2]; //
+ Double_t fParams[4]; //
- ClassDef(AliTRDPIDThresholds, 1);
+ ClassDef(AliTRDPIDThresholds, 1);
};
+ AliTRDPIDParams &operator=(const AliTRDPIDParams &);
+
static const Double_t kVerySmall;
- AliTRDPIDParams(const AliTRDPIDParams &);
- AliTRDPIDParams &operator=(const AliTRDPIDParams &);
-
TSortedList *fEntries; //
ClassDef(AliTRDPIDParams, 1);
//____________________________________________________________
AliTRDPIDReference::AliTRDPIDReference(const AliTRDPIDReference &ref):
TNamed(ref),
- fRefContainer(ref.fRefContainer),
+ fRefContainer(NULL),
fMomentumBins(ref.fMomentumBins)
{
- //
- // Copy constructor
- // Only copies poiters, object is not the owner of the references
- //
- SetBit(kIsOwner, kFALSE);
+ //
+ // Copy constructor
+ //
+ fRefContainer = new TObjArray(fMomentumBins.GetSize() * AliPID::kSPECIES);
+ fRefContainer->SetOwner();
+
+ for(Int_t ip = 0; ip < GetNumberOfMomentumBins(); ip++){
+ for(Int_t is = 0; is < 5; is++){
+ Int_t ent=is * fMomentumBins.GetSize() + ip;
+ TObject *obj=ref.fRefContainer->At(ent);
+ if(obj){
+ fRefContainer->AddAt(obj->Clone(),ent);
+ }
+ }
+ }
+
+ SetBit(kIsOwner, kTRUE);
}
//____________________________________________________________
return;
}
AliDebug(1, Form("Adding object with address %p to position %d", ref, spec * fMomentumBins.GetSize() + pbin));
- fRefContainer->AddAt(ref, spec * fMomentumBins.GetSize() + pbin);
+ fRefContainer->AddAt(ref->Clone(), spec * fMomentumBins.GetSize() + pbin);
}
//____________________________________________________________
#include "AliLog.h"
-#include "AliTRDPIDParams.h"
-#include "AliTRDPIDReference.h"
+#include "AliTRDPIDResponseObject.h"
+//#include "AliTRDPIDParams.h"
+//#include "AliTRDPIDReference.h"
#include "AliTRDPIDResponse.h"
+#include "AliTRDTKDInterpolator.h"
ClassImp(AliTRDPIDResponse)
//____________________________________________________________
AliTRDPIDResponse::AliTRDPIDResponse():
TObject()
- ,fkPIDReference(NULL)
- ,fkPIDParams(NULL)
+ ,fkPIDResponseObject(NULL)
,fGainNormalisationFactor(1.)
,fPIDmethod(kLQ1D)
{
//____________________________________________________________
AliTRDPIDResponse::AliTRDPIDResponse(const AliTRDPIDResponse &ref):
TObject(ref)
- ,fkPIDReference(ref.fkPIDReference)
- ,fkPIDParams(ref.fkPIDParams)
+ ,fkPIDResponseObject(NULL)
,fGainNormalisationFactor(ref.fGainNormalisationFactor)
,fPIDmethod(ref.fPIDmethod)
{
// Make copy
TObject::operator=(ref);
fGainNormalisationFactor = ref.fGainNormalisationFactor;
- fkPIDReference = ref.fkPIDReference;
- fkPIDParams = ref.fkPIDParams;
+ fkPIDResponseObject = ref.fkPIDResponseObject;
fPIDmethod = ref.fPIDmethod;
return *this;
//
// Destructor
//
- if(IsOwner()) delete fkPIDReference;
+ if(IsOwner()) delete fkPIDResponseObject;
}
//____________________________________________________________
-Bool_t AliTRDPIDResponse::Load(const Char_t * filename, const Char_t *refName){
+Bool_t AliTRDPIDResponse::Load(const Char_t * filename){
//
// Load References into the toolkit
//
}
gROOT->cd();
- fkPIDReference = dynamic_cast<const AliTRDPIDReference *>(in->Get(refName)->Clone());
+ fkPIDResponseObject = dynamic_cast<const AliTRDPIDResponseObject *>(in->Get("TRDPIDResponse")->Clone());
in->Close(); delete in;
owd->cd();
SetBit(kIsOwner, kTRUE);
- AliDebug(2, Form("Successfully loaded References for %d Momentum bins", fkPIDReference->GetNumberOfMomentumBins()));
+ AliDebug(2, Form("Successfully loaded References for %d Momentum bins", fkPIDResponseObject->GetNumberOfMomentumBins()));
return kTRUE;
}
// true if calculation success
//
- if(!fkPIDReference){
- AliWarning("Missing reference data. PID calculation not possible.");
- return kFALSE;
- }
+ if(!fkPIDResponseObject){
+ AliWarning("Missing reference data. PID calculation not possible.");
+ return kFALSE;
+ }
- for(Int_t is(AliPID::kSPECIES); is--;) prob[is]=.2;
- Double_t prLayer[AliPID::kSPECIES];
- Double_t dE[10], s(0.);
- for(Int_t il(kNlayer); il--;){
- memset(prLayer, 0, AliPID::kSPECIES*sizeof(Double_t));
- if(!CookdEdx(n, &dedx[il*n], &dE[0])) continue;
+ for(Int_t is(AliPID::kSPECIES); is--;) prob[is]=.2;
+ Double_t prLayer[AliPID::kSPECIES];
+ Double_t dE[10], s(0.);
+ for(Int_t il(kNlayer); il--;){
+ memset(prLayer, 0, AliPID::kSPECIES*sizeof(Double_t));
+ if(!CookdEdx(n, &dedx[il*n], &dE[0])) continue;
- s=0.;
- for(Int_t is(AliPID::kSPECIES); is--;){
- if((dE[0] > 0.) && (p[il] > 0.)) prLayer[is] = GetProbabilitySingleLayer(is, p[il], dE[0]);
- AliDebug(3, Form("Probability for Species %d in Layer %d: %f", is, il, prLayer[is]));
- s+=prLayer[is];
+ s=0.;
+ Bool_t filled=kTRUE;
+ for(Int_t is(AliPID::kSPECIES); is--;){
+ if((dE[0] > 0.) && (p[il] > 0.)) prLayer[is] = GetProbabilitySingleLayer(is, p[il], &dE[0]);
+ AliDebug(3, Form("Probability for Species %d in Layer %d: %f", is, il, prLayer[is]));
+ if(prLayer[is]<1.e-30){
+ AliDebug(2, Form("Null for species %d species prob layer[%d].",is,il));
+ filled=kFALSE;
+ break;
+ }
+ s+=prLayer[is];
+ }
+ if(!filled){
+ continue;
+ }
+ if(s<1.e-30){
+ AliDebug(2, Form("Null all species prob layer[%d].", il));
+ continue;
+ }
+ for(Int_t is(AliPID::kSPECIES); is--;){
+ if(kNorm) prLayer[is] /= s;
+ prob[is] *= prLayer[is];
+ }
}
+ if(!kNorm) return kTRUE;
+
+ s=0.;
+ for(Int_t is(AliPID::kSPECIES); is--;) s+=prob[is];
if(s<1.e-30){
- AliDebug(2, Form("Null all species prob layer[%d].", il));
- continue;
+ AliDebug(2, "Null total prob.");
+ return kFALSE;
}
- for(Int_t is(AliPID::kSPECIES); is--;){
- if(kNorm) prLayer[is] /= s;
- prob[is] *= prLayer[is];
- }
- }
- if(!kNorm) return kTRUE;
-
- s=0.;
- for(Int_t is(AliPID::kSPECIES); is--;) s+=prob[is];
- if(s<1.e-30){
- AliDebug(2, "Null total prob.");
- return kFALSE;
- }
- for(Int_t is(AliPID::kSPECIES); is--;) prob[is]/=s;
- return kTRUE;
+ for(Int_t is(AliPID::kSPECIES); is--;) prob[is]/=s;
+ return kTRUE;
}
//____________________________________________________________
-Double_t AliTRDPIDResponse::GetProbabilitySingleLayer(Int_t species, Double_t plocal, Double_t dEdx) const {
+Double_t AliTRDPIDResponse::GetProbabilitySingleLayer(Int_t species, Double_t plocal, Double_t *dEdx) const {
//
// Get the non-normalized probability for a certain particle species as coming
// from the reference histogram
// Interpolation between momentum bins
//
AliDebug(1, Form("Make Probability for Species %d with Momentum %f", species, plocal));
- Float_t pLower, pUpper;
+
Double_t probLayer = 0.;
- TH1 *refUpper = dynamic_cast<TH1 *>(fkPIDReference->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper)),
- *refLower = dynamic_cast<TH1 *>(fkPIDReference->GetLowerReference((AliPID::EParticleType)species, plocal, pLower));
- // Do Interpolation exept for underflow and overflow
- if(refLower && refUpper){
- Double_t probLower = refLower->GetBinContent(refLower->GetXaxis()->FindBin(dEdx));
- Double_t probUpper = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx));
-
- probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower);
- } else if(refLower){
- // underflow
- probLayer = refLower->GetBinContent(refLower->GetXaxis()->FindBin(dEdx));
- } else if(refUpper){
- // overflow
- probLayer = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx));
- } else {
- AliError("No references available");
+ Float_t pLower, pUpper;
+
+ switch(fPIDmethod){
+ case kNN: // NN
+ break;
+ case kLQ2D: // 2D LQ
+ {
+ Double_t error;
+ Double_t point[kNslicesLQ2D];
+ for(Int_t idim=0;idim<kNslicesLQ2D;idim++){point[idim]=dEdx[idim];}
+
+ AliTRDTKDInterpolator *refLower = dynamic_cast<AliTRDTKDInterpolator*>(fkPIDResponseObject->GetLowerReference((AliPID::EParticleType)species, plocal, pLower,kLQ2D));
+
+ if(refLower){
+ refLower->Eval(point,probLayer,error);
+ }
+ else {
+ AliError("No references available");
+ }
+ }
+ break;
+ case kLQ1D: // 1D LQ
+ {
+ TH1 *refUpper = dynamic_cast<TH1 *>(fkPIDResponseObject->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper,kLQ1D)),
+ *refLower = dynamic_cast<TH1 *>(fkPIDResponseObject->GetLowerReference((AliPID::EParticleType)species, plocal, pLower,kLQ1D));
+ // Do Interpolation exept for underflow and overflow
+ if(refLower && refUpper){
+ Double_t probLower = refLower->GetBinContent(refLower->GetXaxis()->FindBin(dEdx[0]));
+ Double_t probUpper = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx[0]));
+
+ probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower);
+ } else if(refLower){
+ // underflow
+ probLayer = refLower->GetBinContent(refLower->GetXaxis()->FindBin(dEdx[0]));
+ } else if(refUpper){
+ // overflow
+ probLayer = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx[0]));
+ } else {
+ AliError("No references available");
+ }
+ AliDebug(1, Form("Probability %f", probLayer));
+ }
+ break;
+ default:
+ break;
}
- AliDebug(1, Form("Probability %f", probLayer));
return probLayer;
}
//
// Make Deep Copy of the Reference Histograms
//
- if(!fkPIDReference || IsOwner()) return;
- const AliTRDPIDReference *tmp = fkPIDReference;
- fkPIDReference = dynamic_cast<const AliTRDPIDReference *>(tmp->Clone());
- SetBit(kIsOwner, kTRUE);
+ if(!fkPIDResponseObject || IsOwner()) return;
+ const AliTRDPIDResponseObject *tmp = fkPIDResponseObject;
+ fkPIDResponseObject = dynamic_cast<const AliTRDPIDResponseObject *>(tmp->Clone());
+ SetBit(kIsOwner, kTRUE);
}
//____________________________________________________________
//
switch(fPIDmethod){
case kNN: // NN
- break;
- case kLQ2D: // 2D LQ
- break;
- case kLQ1D: // 1D LQ
- out[0]= 0.;
- for(Int_t islice = 0; islice < nSlice; islice++)
- if(in[islice] > 0) out[0] += in[islice] * fGainNormalisationFactor; // Protect against negative values for slices having no dE/dx information
- if(out[0] < 1e-6) return kFALSE;
- break;
+ break;
+ case kLQ2D: // 2D LQ
+ out[0]=0;
+ out[1]=0;
+ for(Int_t islice = 0; islice < nSlice; islice++){
+ if(islice<fkPIDResponseObject->GetNSlicesQ0())out[0]+= in[islice];
+ else out[1]+= in[islice];
+ }
+ break;
+ case kLQ1D: // 1D LQ
+ out[0]= 0.;
+ for(Int_t islice = 0; islice < nSlice; islice++)
+ if(in[islice] > 0) out[0] += in[islice] * fGainNormalisationFactor; // Protect against negative values for slices having no dE/dx information
+ if(out[0] < 1e-6) return kFALSE;
+ break;
+
default:
- return kFALSE;
+ return kFALSE;
}
return kTRUE;
}
//____________________________________________________________
Bool_t AliTRDPIDResponse::IdentifiedAsElectron(Int_t nTracklets, const Double_t *like, Double_t p, Double_t level) const {
- //
- // Check whether particle is identified as electron assuming a certain electron efficiency level
- // Only electron and pion hypothesis is taken into account
- //
- // Inputs:
- // Number of tracklets
- // Likelihood values
- // Momentum
- // Electron efficiency level
- //
- // If the function fails when the params are not accessible, the function returns true
- //
- if(!fkPIDParams){
+ //
+ // Check whether particle is identified as electron assuming a certain electron efficiency level
+ // Only electron and pion hypothesis is taken into account
+ //
+ // Inputs:
+ // Number of tracklets
+ // Likelihood values
+ // Momentum
+ // Electron efficiency level
+ //
+ // If the function fails when the params are not accessible, the function returns true
+ //
+ if(!fkPIDResponseObject){
AliError("No PID Param object available");
return kTRUE;
}
Double_t probEle = like[AliPID::kElectron]/(like[AliPID::kElectron] + like[AliPID::kPion]);
Double_t params[4];
- if(!fkPIDParams->GetThresholdParameters(nTracklets, level, params)){
+ if(!fkPIDResponseObject->GetThresholdParameters(nTracklets, level, params,fPIDmethod)){
AliError("No Params found for the given configuration");
return kTRUE;
}
class TObjArray;
class AliVTrack;
-class AliTRDPIDParams;
-class AliTRDPIDReference;
+class AliTRDPIDResponseObject;
+
class AliTRDPIDResponse : public TObject {
public:
enum ETRDPIDResponseStatus {
kIsOwner = BIT(14)
};
enum ETRDPIDResponseDef {
- kNlayer = 6
- ,kNPBins = 6
+ kNlayer = 6
+ ,kNPBins = 6
};
enum ETRDPIDMethod {
kNN = 0,
kLQ2D = 1,
kLQ1D = 2
};
+ enum ETRDPIDNMethod {
+ kNMethod=3
+ };
enum ETRDNslices {
kNslicesLQ1D = 1,
kNslicesLQ2D = 2,
void SetOwner();
void SetPIDmethod(ETRDPIDMethod m) {fPIDmethod=m;}
void SetGainNormalisationFactor(Double_t gainFactor) { fGainNormalisationFactor = gainFactor; }
- void SetPIDParams(const AliTRDPIDParams * params) { fkPIDParams = params; }
- Bool_t Load(const Char_t *filename = NULL, const Char_t *refName = "RefTRDLQ1D");
- Bool_t Load(const AliTRDPIDReference *ref) { fkPIDReference = ref; return kTRUE; }
+ Bool_t SetPIDResponseObject(const AliTRDPIDResponseObject * obj) { fkPIDResponseObject = obj; return kTRUE;}
+
+ Bool_t Load(const Char_t *filename = NULL);
Bool_t IdentifiedAsElectron(Int_t nTracklets, const Double_t *like, Double_t p, Double_t level) const;
private:
Bool_t CookdEdx(Int_t nSlice, const Double_t * const in, Double_t *out) const;
- Double_t GetProbabilitySingleLayer(Int_t species, Double_t plocal, Double_t dEdx) const;
+ Double_t GetProbabilitySingleLayer(Int_t species, Double_t plocal, Double_t *dEdx) const;
- const AliTRDPIDReference *fkPIDReference; // PID References
- const AliTRDPIDParams *fkPIDParams; // PID Params
+ const AliTRDPIDResponseObject *fkPIDResponseObject; // PID References and Params
Double_t fGainNormalisationFactor; // Gain normalisation factor
ETRDPIDMethod fPIDmethod; // PID method selector
return slices;
}
#endif
+
--- /dev/null
+/**************************************************************************
+ * 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. *
+ **************************************************************************/
+//
+// Container class for the reference distributions for TRD PID
+// The class contains the reference distributions and the momentum steps
+// the references are taken at. Mapping is done inside. To derive references,
+// the functions GetUpperReference and GetLowerReference return the next
+// reference distribution object and the momentum step above respectively below
+// the tracklet momentum.
+//
+// Authors:
+// Markus Fasel <M.Fasel@gsi.de>
+// Daniel Lohner <Daniel.Lohner@cern.ch>
+
+#include "AliLog.h"
+
+#include "AliTRDPIDResponseObject.h"
+
+#ifndef AliTRDPIDREFERENCE_H
+#include "AliTRDPIDReference.h"
+#endif
+
+#ifndef AliTRDPIDPARAMS_H
+#include "AliTRDPIDParams.h"
+#endif
+
+
+ClassImp(AliTRDPIDResponseObject)
+
+//____________________________________________________________
+AliTRDPIDResponseObject::AliTRDPIDResponseObject():
+ TNamed(),
+ fNSlicesQ0(4)
+{
+ //
+ // Dummy constructor
+ //
+ SetBit(kIsOwner, kTRUE);
+
+ for(Int_t method=0;method<AliTRDPIDResponse::kNMethod;method++){
+ fkPIDParams[method]=NULL;
+ fkPIDReference[method]=NULL;
+ }
+}
+
+//____________________________________________________________
+AliTRDPIDResponseObject::AliTRDPIDResponseObject(const Char_t *name):
+TNamed(name, "TRD PID Response Object"),
+fNSlicesQ0(4)
+{
+ //
+ // Default constructor
+ //
+ SetBit(kIsOwner, kTRUE);
+
+ for(Int_t method=0;method<AliTRDPIDResponse::kNMethod;method++){
+ fkPIDParams[method]=NULL;
+ fkPIDReference[method]=NULL;
+ }
+}
+
+//____________________________________________________________
+AliTRDPIDResponseObject::AliTRDPIDResponseObject(const AliTRDPIDResponseObject &ref):
+TNamed(ref),
+fNSlicesQ0(ref.fNSlicesQ0)
+{
+ //
+ // Copy constructor
+ // Only copies pointers, object is not the owner of the references
+ //
+ SetBit(kIsOwner, kFALSE);
+
+ for(Int_t method=0;method<AliTRDPIDResponse::kNMethod;method++){
+ fkPIDParams[method]=ref.fkPIDParams[method]; // new Object is not owner, copy only pointer
+ fkPIDReference[method]=ref.fkPIDReference[method]; // new Object is not owner, copy only pointer
+ }
+}
+//____________________________________________________________
+AliTRDPIDResponseObject &AliTRDPIDResponseObject::operator=(const AliTRDPIDResponseObject &ref){
+ //
+ // Assginment operator
+ // Only copies poiters, object is not the owner of the references
+ //
+ if(this != &ref){
+ TNamed::operator=(ref);
+ fNSlicesQ0=ref.fNSlicesQ0;
+ for(Int_t method=0;method<AliTRDPIDResponse::kNMethod;method++){
+ if(TestBit(kIsOwner) && fkPIDParams[method])delete fkPIDParams[method];
+ if(TestBit(kIsOwner) && fkPIDReference[method])delete fkPIDReference[method];
+
+ fkPIDParams[method]=ref.fkPIDParams[method]; // new Object is not owner, copy only pointer
+ fkPIDReference[method]=ref.fkPIDReference[method]; // new Object is not owner, copy only pointer
+ }
+ SetBit(kIsOwner, kFALSE);
+ }
+ return *this;
+}
+
+//____________________________________________________________
+AliTRDPIDResponseObject::~AliTRDPIDResponseObject(){
+ //
+ // Destructor
+ // references are deleted if the object is the owner
+ //
+ if(fkPIDParams && TestBit(kIsOwner)) delete fkPIDParams;
+ if(fkPIDReference && TestBit(kIsOwner)) delete fkPIDReference;
+
+}
+
+//____________________________________________________________
+void AliTRDPIDResponseObject::SetPIDParams(AliTRDPIDParams *params,AliTRDPIDResponse::ETRDPIDMethod method){
+
+ if(Int_t(method)>=Int_t(AliTRDPIDResponse::kNMethod)||Int_t(method)<0){
+ AliError("Method does not exist");
+ return;
+ }
+ if(fkPIDParams[method]){
+ delete fkPIDParams[method];
+ fkPIDParams[method]=NULL;
+ }
+
+ fkPIDParams[method]=new AliTRDPIDParams(*params);
+}
+
+//____________________________________________________________
+void AliTRDPIDResponseObject::SetPIDReference(AliTRDPIDReference *reference,AliTRDPIDResponse::ETRDPIDMethod method){
+
+ if(Int_t(method)>=Int_t(AliTRDPIDResponse::kNMethod)||Int_t(method)<0){
+ AliError("Method does not exist");
+ return;
+ }
+ if(fkPIDReference[method]){
+ delete fkPIDReference[method];
+ fkPIDReference[method]=NULL;
+ }
+ fkPIDReference[method]=new AliTRDPIDReference(*reference);
+}
+
+//____________________________________________________________
+TObject *AliTRDPIDResponseObject::GetUpperReference(AliPID::EParticleType spec, Float_t p, Float_t &pUpper,AliTRDPIDResponse::ETRDPIDMethod method) const{
+
+ if(Int_t(method)>=Int_t(AliTRDPIDResponse::kNMethod)||Int_t(method)<0){
+ AliError("Method does not exist");
+ return NULL;
+ }
+
+ if(fkPIDReference[method]){
+ return fkPIDReference[method]->GetUpperReference(spec,p,pUpper);
+ }
+ return NULL;
+}
+
+
+//____________________________________________________________
+TObject *AliTRDPIDResponseObject::GetLowerReference(AliPID::EParticleType spec, Float_t p, Float_t &pLower,AliTRDPIDResponse::ETRDPIDMethod method) const{
+
+ if(Int_t(method)>=Int_t(AliTRDPIDResponse::kNMethod)||Int_t(method)<0){
+ AliError("Method does not exist");
+ return NULL;
+ }
+
+ if(fkPIDReference[method]){
+ return fkPIDReference[method]->GetLowerReference(spec,p,pLower);
+ }
+ return NULL;
+}
+
+//____________________________________________________________
+Bool_t AliTRDPIDResponseObject::GetThresholdParameters(Int_t ntracklets, Double_t efficiency, Double_t *params,AliTRDPIDResponse::ETRDPIDMethod method) const{
+
+ if(Int_t(method)>=Int_t(AliTRDPIDResponse::kNMethod)||Int_t(method)<0){
+ AliError("Method does not exist");
+ return kFALSE;
+ }
+
+ if(fkPIDParams[method]){
+ return fkPIDParams[method]->GetThresholdParameters(ntracklets,efficiency,params);
+ }
+ return kFALSE;
+}
+
+//____________________________________________________________
+Int_t AliTRDPIDResponseObject::GetNumberOfMomentumBins(AliTRDPIDResponse::ETRDPIDMethod method) const{
+
+ if(Int_t(method)>=Int_t(AliTRDPIDResponse::kNMethod)||Int_t(method)<0){
+ AliError("Method does not exist");
+ return 0;
+ }
+
+ if(fkPIDReference[method]){
+ return fkPIDReference[method]->GetNumberOfMomentumBins();
+ }
+ return 0;
+}
+
+//____________________________________________________________
+void AliTRDPIDResponseObject::Print(const Option_t* opt) const{
+ //
+ // Print content of the PID object
+ //
+ printf("Content of AliTRDPIDResponseObject \n\n");
+
+ for(Int_t method=0;method<AliTRDPIDResponse::kNMethod;method++){
+ if(fkPIDReference[method])fkPIDReference[method]->Print(opt);
+ if(fkPIDParams[method])printf("+ Threshold Parameters \n");
+ }
+}
--- /dev/null
+/**************************************************************************
+* 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. *
+**************************************************************************/
+//
+// Container for TRD PID Response Objects stored in the OADB
+//
+#ifndef ALITRDPIDRESPONSEOBJECT_H
+#define ALITRDPIDRESPONSEOBJECT_H
+
+#ifndef ROOT_TNamed
+#include <TNamed.h>
+#endif
+
+#ifndef AliTRDPIDRESPONSE_H
+#include "AliTRDPIDResponse.h"
+#endif
+
+
+class AliTRDPIDParams;
+class AliTRDPIDReference;
+class AliTRDPIDResponse;
+
+class AliTRDPIDResponseObject : public TNamed{
+public:
+ enum ETRDPIDResponseObjectStatus {
+ kIsOwner = BIT(14)
+ };
+
+ AliTRDPIDResponseObject();
+ AliTRDPIDResponseObject(const char *name);
+ AliTRDPIDResponseObject(const AliTRDPIDResponseObject &ref);
+ AliTRDPIDResponseObject &operator=(const AliTRDPIDResponseObject &ref);
+ virtual ~AliTRDPIDResponseObject();
+
+ virtual void Print(Option_t *opt) const;
+
+ void SetPIDParams(AliTRDPIDParams *params,AliTRDPIDResponse::ETRDPIDMethod method=AliTRDPIDResponse::kLQ1D);
+ void SetPIDReference(AliTRDPIDReference *params,AliTRDPIDResponse::ETRDPIDMethod method=AliTRDPIDResponse::kLQ1D);
+
+ // Derive reference
+ TObject *GetLowerReference(AliPID::EParticleType spec, Float_t p, Float_t &pLower,AliTRDPIDResponse::ETRDPIDMethod method=AliTRDPIDResponse::kLQ1D) const;
+ TObject *GetUpperReference(AliPID::EParticleType spec, Float_t p, Float_t &pUpper,AliTRDPIDResponse::ETRDPIDMethod method=AliTRDPIDResponse::kLQ1D) const;
+
+ Int_t GetNumberOfMomentumBins(AliTRDPIDResponse::ETRDPIDMethod method=AliTRDPIDResponse::kLQ1D) const;
+
+ // Derive threshold params
+ Bool_t GetThresholdParameters(Int_t ntracklets, Double_t efficiency, Double_t *params,AliTRDPIDResponse::ETRDPIDMethod method=AliTRDPIDResponse::kLQ1D) const;
+
+ // Number of SlicesQ0
+ Int_t GetNSlicesQ0() const{return fNSlicesQ0;}
+ void SetNSlicesQ0(Int_t nsl){fNSlicesQ0=nsl;}
+
+private:
+
+ const AliTRDPIDParams *fkPIDParams[AliTRDPIDResponse::kNMethod]; // Contains Thresholds
+ const AliTRDPIDReference *fkPIDReference[AliTRDPIDResponse::kNMethod]; // Contains References
+ Int_t fNSlicesQ0; // Number of Slices for Q0
+
+ ClassDef(AliTRDPIDResponseObject, 1);
+};
+#endif
--- /dev/null
+#include "AliTRDTKDInterpolator.h"
+
+#include "TClonesArray.h"
+#include "TLinearFitter.h"
+#include "TMath.h"
+#include "TRandom.h"
+
+#include "iostream"
+using namespace std;
+
+ClassImp(AliTRDTKDInterpolator)
+
+//_________________________________________________________________
+AliTRDTKDInterpolator::AliTRDTKDInterpolator() :
+TKDTreeIF(),
+fNDataNodes(0),
+fNodes(NULL),
+fLambda(0),
+fNPointsI(0),
+fUseHelperNodes(kFALSE),
+fUseWeights(kFALSE),
+fPDFMode(kInterpolation)
+{
+ // default constructor
+}
+
+//_________________________________________________________________
+AliTRDTKDInterpolator::AliTRDTKDInterpolator(Int_t npoints, Int_t ndim, UInt_t bsize, Float_t **data) :
+TKDTreeIF(npoints, ndim, bsize, data),
+fNDataNodes(0),
+fNodes(NULL),
+fLambda(1 + ndim + (ndim*(ndim+1)>>1)),
+fNPointsI(100),
+fUseHelperNodes(kFALSE),
+fUseWeights(kFALSE),
+fPDFMode(kInterpolation)
+{
+}
+
+//_________________________________________________________________
+AliTRDTKDInterpolator::~AliTRDTKDInterpolator()
+{
+ if(fNodes){
+ fNodes->Delete();
+ delete fNodes;
+ fNodes=NULL;
+ }
+}
+//_________________________________________________________________
+AliTRDTKDInterpolator::AliTRDTKDInterpolator(const AliTRDTKDInterpolator &ref):
+TKDTreeIF(),
+fNDataNodes(ref.fNDataNodes),
+fNodes(ref.fNodes),
+fLambda(ref.fLambda),
+fNPointsI(ref.fNPointsI),
+fUseHelperNodes(ref.fUseHelperNodes),
+fUseWeights(ref.fUseWeights),
+fPDFMode(ref.fPDFMode)
+{
+ // Copy constructor
+}
+
+//____________________________________________________________
+
+AliTRDTKDInterpolator &AliTRDTKDInterpolator::operator=(const AliTRDTKDInterpolator &ref){
+ //
+ // Assignment operator
+ //
+ if(this == &ref) return *this;
+
+ // Make copy
+ TObject::operator=(ref);
+
+ return *this;
+}
+
+//_________________________________________________________________
+Bool_t AliTRDTKDInterpolator::Build()
+{
+ TKDTreeIF::Build();
+ if(!fBoundaries) MakeBoundaries();
+
+ // allocate interpolation nodes
+ fNDataNodes = fNPoints/fBucketSize + ((fNPoints%fBucketSize)?1:0);/*TKDTreeIF::GetNTNodes();*/
+
+ if(fNodes){
+ Warning("AliTRDTKDInterpolator::Build()", "Data already allocated.");
+ fNodes->Delete();
+ } else {
+ fNodes = new TClonesArray("AliTRDTKDInterpolator::AliTRDTKDNodeInfo", fNDataNodes);
+ fNodes->SetOwner();
+ }
+ for(int in=0; in<fNDataNodes; in++) new ((*fNodes)[in]) AliTRDTKDNodeInfo(fNDim);
+
+ // Set Interpolator nodes
+
+ for(int inode=0, tnode = fNNodes; inode<fNDataNodes-1; inode++, tnode++){
+ AliTRDTKDNodeInfo *node =GetNodeInfo(inode);
+ memcpy(node->fBounds,GetBoundary(tnode),2*fNDim*sizeof(Float_t));
+ node->fVal[0] = Float_t(fBucketSize)/fNPoints;
+ for(int idim=0; idim<fNDim; idim++) node->fVal[0] /= (node->fBounds[2*idim+1] - node->fBounds[2*idim]);
+ node->fVal[1] = node->fVal[0]/TMath::Sqrt(float(fBucketSize));
+
+ Int_t *indexPoints = GetPointsIndexes(tnode);
+ // loop points in this terminal node
+ for(int idim=0; idim<fNDim; idim++){
+ node->fData[idim] = 0.;
+ for(int ip = 0; ip<fBucketSize; ip++) node->fData[idim] += fData[idim][indexPoints[ip]];
+ node->fData[idim] /= fBucketSize;
+ }
+ }
+
+ // Analyze last (incomplete) terminal node
+
+ Int_t counts = fNPoints%fBucketSize;
+ counts = counts ? counts : fBucketSize;
+ Int_t inode = fNDataNodes - 1, tnode = inode + fNNodes;
+ AliTRDTKDNodeInfo *ftnode = GetNodeInfo(inode);
+ ftnode->fVal[0] = Float_t(counts)/fNPoints;
+ memcpy(ftnode->fBounds,GetBoundary(tnode),2*fNDim*sizeof(Float_t));
+
+ for(int idim=0; idim<fNDim; idim++){
+ Float_t dx = ftnode->fBounds[2*idim+1]-ftnode->fBounds[2*idim];
+ if(dx < 1.e-30){
+ Warning("AliTRDTKDInterpolator::Build()", "Terminal bucket index[%d] too narrow on the %d dimension.", inode, idim);
+ continue;
+ }
+ ftnode->fVal[0] /= (ftnode->fBounds[2*idim+1] - ftnode->fBounds[2*idim]);
+ }
+ ftnode->fVal[1] = ftnode->fVal[0]/TMath::Sqrt(float(counts));
+
+ // loop points in this terminal node
+ Int_t *indexPoints = GetPointsIndexes(tnode);
+ for(int idim=0; idim<fNDim; idim++){
+ ftnode->fData[idim] = 0.;
+ for(int ip = 0; ip<counts; ip++) ftnode->fData[idim] += fData[idim][indexPoints[ip]];
+ ftnode->fData[idim] /= counts;
+ }
+
+ delete [] fBoundaries;
+ fBoundaries = NULL;
+ // Add Helper Nodes
+ if(fUseHelperNodes){BuildBoundaryNodes();}
+
+ if(fNPointsI>GetNTNodes()){fNPointsI=GetNTNodes();}
+
+ BuildInterpolation();
+
+ return kTRUE;
+}
+
+
+//_________________________________________________________________
+Bool_t AliTRDTKDInterpolator::Eval(const Double_t *point, Double_t &result, Double_t &error)
+{
+ Float_t pointF[fNDim]; // local Float_t conversion for "point"
+ for(int idim=0; idim<fNDim; idim++) pointF[idim] = (Float_t)point[idim];
+ Int_t nodeIndex = GetNodeIndex(pointF);
+ if(nodeIndex<0){
+ Error("AliTRDTKDInterpolator::Eval()", "Can not retrive node for data point.");
+ result = 0.;
+ error = 1.E10;
+ return kFALSE;
+ }
+ AliTRDTKDNodeInfo *node =GetNodeInfo(nodeIndex);
+
+ switch(fPDFMode){
+
+ case kInterpolation:
+ return node->CookPDF(point, result, error);
+ case kMinError:
+ node->CookPDF(point, result, error);
+ if(error<node->fVal[1]){
+ return kTRUE;
+ }
+ error=node->fVal[1];
+ result=node->fVal[0];
+ return kTRUE;
+ case kNodeVal:
+ error=node->fVal[1];
+ result=node->fVal[0];
+ return kTRUE;
+ default:
+ return kFALSE;
+ }
+}
+
+//__________________________________________________________________
+void AliTRDTKDInterpolator::Print(const Option_t */*opt*/) const
+{
+ for(Int_t i=GetNTNodes(); i--;){
+ printf("Node %d of %d: \n",i,GetNTNodes());
+ GetNodeInfo(i)->Print();
+ }
+
+}
+
+//__________________________________________________________________
+Int_t AliTRDTKDInterpolator::GetNodeIndex(const Float_t *p)
+{
+ Int_t inode=FindNode(p)-fNDataNodes+1;
+ if(GetNodeInfo(inode)->Has(p))return inode;
+
+ // Search extra nodes
+
+ for(inode=fNDataNodes;inode<GetNTNodes();inode++){
+ if(GetNodeInfo(inode)->Has(p)){return inode;}
+ }
+
+ // Search for nearest neighbor
+ Float_t dist;
+ Float_t closestdist=10000;
+ inode=-1;
+ for(Int_t ii=0;ii<GetNTNodes();ii++){
+ AliTRDTKDNodeInfo *node=GetNodeInfo(ii);
+ dist=0;
+ for(Int_t idim=0;idim<fNDim;idim++){
+ dist+=TMath::Power((node->fData[idim]-p[idim]),2);
+ }
+ dist=TMath::Sqrt(dist);
+ if(dist<closestdist){closestdist=dist;inode=ii;}
+ }
+ return inode;
+}
+
+//_________________________________________________________________
+AliTRDTKDInterpolator::AliTRDTKDNodeInfo* AliTRDTKDInterpolator::GetNodeInfo(Int_t inode) const
+{
+ if(!fNodes || inode >= GetNTNodes()) return NULL;
+ return (AliTRDTKDNodeInfo*)(*fNodes)[inode];
+}
+
+//_________________________________________________________________
+Int_t AliTRDTKDInterpolator::GetNTNodes() const
+{
+ return fNodes?fNodes->GetEntriesFast():0;
+}
+
+//_________________________________________________________________
+Bool_t AliTRDTKDInterpolator::GetRange(Int_t idim,Float_t range[2]) const
+{
+ if(!fNodes) return kFALSE;
+ if(idim<0 || idim>=fNDim){
+ range[0]=0.; range[1]=0.;
+ return kFALSE;
+ }
+ range[0]=1.e10; range[1]=-1.e10;
+ for(Int_t in=GetNTNodes(); in--; ){
+ AliTRDTKDNodeInfo *node = GetNodeInfo(in);
+
+ if(node->fBounds[2*idim]<range[0]) range[0] = node->fBounds[2*idim];
+ if(node->fBounds[2*idim+1]>range[1]) range[1] = node->fBounds[2*idim+1];
+ }
+
+ return kTRUE;
+}
+
+//_________________________________________________________________
+TH2Poly *AliTRDTKDInterpolator::Projection(Int_t xdim,Int_t ydim)
+{
+ Float_t rangex[2],rangey[2];
+ GetRange(xdim,rangex);
+ GetRange(ydim,rangey);
+
+ TH2Poly* h2 = new TH2Poly("hTKDnodes","hTKDnodes", rangex[0],rangex[1],rangey[0],rangey[1]);
+ h2->GetXaxis()->SetTitle(Form("Q_{%d}", xdim));
+ h2->GetYaxis()->SetTitle(Form("Q_{%d}", ydim));
+
+ for(Int_t inode=0;inode<GetNTNodes();inode++){
+
+ AliTRDTKDNodeInfo* node=GetNodeInfo(inode);
+ h2->AddBin(node->fBounds[2*xdim],node->fBounds[2*ydim],node->fBounds[2*xdim+1],node->fBounds[2*ydim+1]);
+ h2->SetBinContent(inode+1, node->fVal[0]);
+ }
+ return h2;
+}
+
+//_________________________________________________________________
+void AliTRDTKDInterpolator::BuildInterpolation()
+{
+
+ // Calculate Interpolation
+
+ Double_t buffer[fLambda];
+
+ Float_t **refPoints = new Float_t*[fNDim];
+ for(int id=0; id<fNDim; id++){
+ refPoints[id] = new Float_t[GetNTNodes()];
+ for(int in=0; in<GetNTNodes(); in++) refPoints[id][in] = GetNodeInfo(in)->fData[id];
+ }
+ TKDTreeIF *KDhelper = new TKDTreeIF(GetNTNodes(), fNDim, 30, refPoints);
+ KDhelper->Build();
+ KDhelper->MakeBoundariesExact();
+
+ Float_t dist[fNPointsI];
+ Int_t ind[fNPointsI];
+
+ TLinearFitter fitter(fLambda, Form("hyp%d", fLambda-1));
+
+ Int_t nodeIndex(0); Float_t param[6], *pp(NULL);
+ nodeIndex=GetNTNodes(); pp=¶m[0];
+ while(nodeIndex--){
+
+ fitter.ClearPoints();
+
+ AliTRDTKDNodeInfo *node = GetNodeInfo(nodeIndex);
+ // find nearest neighbors
+ KDhelper->FindNearestNeighbors(node->fData,fNPointsI, &ind[0], &dist[0]);
+
+ for(int in=0;in<fNPointsI;in++){
+ AliTRDTKDNodeInfo *nnode = GetNodeInfo(ind[in]);
+
+ Float_t w=1; //weight
+ // calculate tri-cubic weighting function
+ if(fUseWeights){
+ Float_t d = dist[in]/dist[fNPointsI-1];
+ Float_t w0 = (1. - d*d*d);
+ w = w0*w0*w0;
+ if(w<1.e-30) continue;
+ }
+ Int_t ipar=0;
+ for(int idim=0; idim<fNDim; idim++){
+ buffer[ipar++] = nnode->fData[idim];
+ for(int jdim=idim; jdim<fNDim; jdim++) buffer[ipar++] = nnode->fData[idim]*nnode->fData[jdim];
+ }
+ fitter.AddPoint(buffer,nnode->fVal[0], nnode->fVal[1]/w);
+
+ // Ensure Boundary Condition
+ for(Int_t kdim=0;kdim<fNDim;kdim++){
+ if(node->fBounds[2*kdim]==0){
+ Float_t zdata[fNDim];
+ memcpy(&zdata[0],node->fData,fNDim*sizeof(Float_t));
+ zdata[kdim]=0;
+ ipar=0;
+ for(int idim=0; idim<fNDim; idim++){
+ buffer[ipar++] = zdata[idim];
+ for(int jdim=idim; jdim<fNDim; jdim++) buffer[ipar++] = zdata[idim]*zdata[jdim];
+ }
+ fitter.AddPoint(buffer,0,1);
+ }
+ }
+ }
+
+ fitter.Eval();
+
+ // retrive fitter results
+ TMatrixD cov(fLambda, fLambda);
+ TVectorD par(fLambda);
+ fitter.GetCovarianceMatrix(cov);
+ fitter.GetParameters(par);
+
+ // store results
+ node->Store(&par,&cov);
+ }
+
+ delete KDhelper;
+ for(int id=0; id<fNDim; id++){
+ delete refPoints[id];
+ }
+ delete[] refPoints;
+}
+
+//_________________________________________________________________
+void AliTRDTKDInterpolator::BuildBoundaryNodes(){
+
+ Int_t nnew=0;
+
+ Float_t treebounds[2*fNDim];
+ for(Int_t idim=0;idim<fNDim;idim++){
+ GetRange(idim,&treebounds[2*idim]);
+ }
+
+ for(int inode=0; inode<GetNTNodes(); inode++){
+
+ AliTRDTKDNodeInfo *node=GetNodeInfo(inode);
+
+ for(Int_t vdim=0;vdim<fNDim;vdim++){
+
+ // Try expansion to lower and higher values
+ for(Int_t iter=0;iter<2;iter++){
+ if(node->fBounds[2*vdim+iter]==treebounds[2*vdim+iter]){
+
+ // Add new Node
+ new ((*fNodes)[GetNTNodes()]) AliTRDTKDNodeInfo(fNDim);
+
+ AliTRDTKDNodeInfo *newnode = GetNodeInfo(GetNTNodes()-1);
+ if(iter==0)newnode->fBounds[2*vdim+iter]=0;
+ if(iter==1)newnode->fBounds[2*vdim+iter]=2*treebounds[2*vdim+iter];
+ newnode->fBounds[2*vdim+!iter]=node->fBounds[2*vdim+iter];
+ for(Int_t idim=0;idim<fNDim;idim++){
+ if(idim==vdim)continue;
+ newnode->fBounds[2*idim]=node->fBounds[2*idim];
+ newnode->fBounds[2*idim+1]=node->fBounds[2*idim+1];
+ }
+ newnode->fVal[0]=0;
+ newnode->fVal[1]=Float_t(1)/fNPoints;
+ for(int idim=0; idim<fNDim; idim++){
+ newnode->fVal[1] /= (newnode->fBounds[2*idim+1] - newnode->fBounds[2*idim]);
+ newnode->fData[idim]=0.5*(newnode->fBounds[2*idim+1] + newnode->fBounds[2*idim]);
+ }
+ nnew++;
+ }
+ }
+ }
+ }
+ printf("%d Boundary Nodes Added \n",nnew);
+}
+
+//_________________________________________________________________
+AliTRDTKDInterpolator::AliTRDTKDNodeInfo::AliTRDTKDNodeInfo(Int_t ndim):
+ TObject()
+ ,fNDim(ndim)
+ ,fNBounds(2*ndim)
+ ,fNPar(1 + ndim + (ndim*(ndim+1)>>1))
+ ,fNCov(fNPar*fNPar)
+ ,fData(NULL)
+ ,fBounds(NULL)
+ ,fPar(NULL)
+ ,fCov(NULL)
+{
+ // Default constructor
+ fVal[0] = 0.; fVal[1] = 0.;
+ fData=new Float_t[fNDim];
+ fBounds=new Float_t[fNBounds];
+}
+
+//_________________________________________________________________
+AliTRDTKDInterpolator::AliTRDTKDNodeInfo::AliTRDTKDNodeInfo(const AliTRDTKDInterpolator::AliTRDTKDNodeInfo &ref):
+TObject(ref)
+ ,fNDim(ref.fNDim)
+ ,fNBounds(ref.fNBounds)
+ ,fNPar(ref.fNPar)
+ ,fNCov(ref.fNCov)
+ ,fData(NULL)
+ ,fBounds(NULL)
+ ,fPar(NULL)
+ ,fCov(NULL)
+{
+ // Copy constructor
+
+ if(ref.fData){
+ fData = new Float_t[fNDim];
+ memcpy(fData, ref.fData, fNDim*sizeof(Float_t));
+ }
+ if(ref.fBounds){
+ fBounds = new Float_t[2*fNDim];
+ memcpy(fBounds, ref.fBounds, 2*fNDim*sizeof(Float_t));
+ }
+
+ fVal[0] = ref.fVal[0];
+ fVal[1] = ref.fVal[1];
+
+ if(ref.fPar){
+ fPar=new Double_t[fNPar];
+ memcpy(fPar, ref.fPar, fNPar*sizeof(Double_t));
+ }
+ if(ref.fCov){
+ fCov=new Double_t[fNCov];
+ memcpy(fCov, ref.fCov, fNCov*sizeof(Double_t));
+ }
+}
+
+//_________________________________________________________________
+AliTRDTKDInterpolator::AliTRDTKDNodeInfo::~AliTRDTKDNodeInfo()
+{
+ // Destructor
+ if(fData) delete [] fData;
+ if(fBounds) delete [] fBounds;
+ if(fPar) delete [] fPar;
+ if(fCov) delete [] fCov;
+}
+
+//_________________________________________________________________
+AliTRDTKDInterpolator::AliTRDTKDNodeInfo &AliTRDTKDInterpolator::AliTRDTKDNodeInfo::operator=(const AliTRDTKDInterpolator::AliTRDTKDNodeInfo &ref)
+{
+ //
+ // Assignment operator
+ //
+
+ if(&ref==this) return *this;
+ memcpy(fData, ref.fData, fNDim*sizeof(Float_t));
+ memcpy(fBounds, ref.fBounds, 2*fNDim*sizeof(Float_t));
+
+ fVal[0] = ref.fVal[0];
+ fVal[1] = ref.fVal[1];
+
+ if(ref.fPar){
+ fPar=new Double_t[fNPar];
+ memcpy(fPar, ref.fPar, fNPar*sizeof(Double_t));
+ }
+ if(ref.fCov){
+ fCov=new Double_t[fNCov];
+ memcpy(fCov, ref.fCov, fNCov*sizeof(Double_t));
+ }
+ return *this;
+}
+
+
+//_________________________________________________________________
+void AliTRDTKDInterpolator::AliTRDTKDNodeInfo::Print(const Option_t *opt) const
+{
+ // Print the content of the node
+ printf("x [");
+ for(int idim=0; idim<fNDim; idim++) printf("%f ", fData?fData[idim]:0.);
+ printf("] f = [%e +- %e]\n", fVal[0], fVal[1]);
+
+ if(fBounds){
+ printf("range [");
+ for(int idim=0; idim<fNDim; idim++) printf("(%f %f) ", fBounds[2*idim], fBounds[2*idim+1]);
+ printf("]\n");
+ }
+ if(strcmp(opt, "a")!=0) return;
+
+ if(fPar){
+ printf("Fit parameters : \n");
+ for(int ip=0; ip<fNPar; ip++) printf("p%d[%f] ", ip, fPar[ip]);
+ printf("\n");
+ }
+ if(!fCov) return;
+ for(int ip(0), n(0); ip<fNPar; ip++){
+ for(int jp(ip); jp<fNPar; jp++) printf("c(%d %d)[%f] ", ip, jp, fCov[n++]);
+ printf("\n");
+ }
+}
+
+//_________________________________________________________________
+void AliTRDTKDInterpolator::AliTRDTKDNodeInfo::Store(TVectorD const *par, TMatrixD const *cov)
+{
+// Store the parameters and the covariance in the node
+
+ if(!fPar){fPar = new Double_t[fNPar];}
+ for(int ip=0; ip<fNPar; ip++) fPar[ip] = (*par)[ip];
+
+ if(!cov) return;
+ if(!fCov){fCov = new Double_t[fNCov];}
+ for(int ip(0), np(0); ip<fNPar; ip++)
+ for(int jp=ip; jp<fNPar; jp++) fCov[np++] = (*cov)(ip,jp);
+}
+
+//_________________________________________________________________
+Bool_t AliTRDTKDInterpolator::AliTRDTKDNodeInfo::CookPDF(const Double_t *point, Double_t &result, Double_t &error) const
+{
+ // Recalculate the PDF for one node from the results of interpolation (parameters and covariance matrix)
+
+ result =0.; error = 1.;
+ if(!fPar) return kFALSE;
+
+ Double_t fdfdp[fNDim];
+ Int_t ipar = 0;
+ fdfdp[ipar++] = 1.;
+ for(int idim=0; idim<fNDim; idim++){
+ fdfdp[ipar++] = point[idim];
+ for(int jdim=idim; jdim<fNDim; jdim++) fdfdp[ipar++] = point[idim]*point[jdim];
+ }
+
+ // calculate estimation
+ for(int i=0; i<fNPar; i++) result += fdfdp[i]*fPar[i];
+
+ if(!fCov)return kTRUE;
+ // calculate error
+ error=0;
+
+ for(int i(0), n(0); i<fNPar; i++){
+ error += fdfdp[i]*fdfdp[i]*fCov[n++];
+ for(int j(i+1); j<fNPar; j++) error += 2.*fdfdp[i]*fdfdp[j]*fCov[n++];
+ }
+ if(error>0)error = TMath::Sqrt(error);
+ else{error=100;}
+
+ // Boundary condition
+ if(result<0){
+ result=fVal[0];
+ error=fVal[1];
+ }
+
+ return kTRUE;
+}
+
+//_____________________________________________________________________
+Bool_t AliTRDTKDInterpolator::AliTRDTKDNodeInfo::Has(const Float_t *p) const
+{
+ Int_t n(0);
+ for(int id=0; id<fNDim; id++)
+ if(p[id]>=fBounds[2*id] && p[id]<fBounds[2*id+1]) n++;
+ if(n==fNDim) return kTRUE;
+ return kFALSE;
+}
+
--- /dev/null
+#ifndef ROOT_ALITRDTKDINTERPOLATOR_H
+#define ROOT_ALITRDTKDINTERPOLATOR_H
+
+#ifndef ROOT_TKDTree
+#include "TKDTree.h"
+#endif
+
+#ifndef ROOT_TObject
+#include "TObject.h"
+#endif
+
+#include "TVectorD.h"
+#include "TMatrixD.h"
+#include "TH2Poly.h"
+class TClonesArray;
+
+class AliTRDTKDInterpolator : public TKDTreeIF
+{
+
+
+
+
+public:
+ enum TRDTKDMode{
+ kInterpolation=0,
+ kMinError=1,
+ kNodeVal=2
+ };
+
+ // Bucket Object class
+ class AliTRDTKDNodeInfo : public TObject
+ {
+ friend class AliTRDTKDInterpolator;
+ public:
+
+ AliTRDTKDNodeInfo(Int_t ndim = 0);
+ AliTRDTKDNodeInfo(const AliTRDTKDInterpolator::AliTRDTKDNodeInfo &ref);
+ AliTRDTKDInterpolator::AliTRDTKDNodeInfo& operator=(const AliTRDTKDInterpolator::AliTRDTKDNodeInfo &ref);
+ virtual ~AliTRDTKDNodeInfo();
+ Bool_t CookPDF(const Double_t *point, Double_t &result, Double_t &error) const;
+ Bool_t Has(const Float_t *p) const;
+ void Print(const Option_t * = "") const;
+ void Store(TVectorD const *par, TMatrixD const *cov=NULL);
+
+ private:
+ Int_t fNDim; // Dimension of Points
+ Int_t fNBounds; // 2* Dimension of Points
+ Int_t fNPar; // Number of Parameters
+ Int_t fNCov; // Size of Cov Matrix
+ Float_t *fData; //[fNDim] Data Point
+ Float_t *fBounds; //[fNBounds] Boundaries
+ Float_t fVal[2]; //measured value for node
+ Double_t *fPar; //[fNPar] interpolator parameters
+ Double_t *fCov; //[fNCov] interpolator covariance matrix
+
+ ClassDef(AliTRDTKDNodeInfo, 1) // node info for interpolator
+ };
+
+public:
+
+ AliTRDTKDInterpolator();
+ AliTRDTKDInterpolator(Int_t npoints, Int_t ndim, UInt_t bsize, Float_t **data);
+ ~AliTRDTKDInterpolator();
+
+ Bool_t Eval(const Double_t *point, Double_t &result, Double_t &error);
+ void Print(const Option_t *opt="") const;
+
+ TH2Poly * Projection(Int_t xdim,Int_t ydim);
+
+ Int_t GetNDIM() const {return fNDim;}
+ Bool_t GetRange(Int_t idim,Float_t range[2]) const;
+
+ void SetNPointsInterpolation(Int_t np){fNPointsI=np;};
+ Int_t GetNPointsInterpolation(){return fNPointsI;};
+
+ void SetUseWeights(Bool_t k=kTRUE){fUseWeights=k;}
+ void SetPDFMode(TRDTKDMode mod){fPDFMode=mod;}
+
+ Bool_t Build();
+
+ void SetUseHelperNodes(Bool_t k){fUseHelperNodes=k;}
+
+private:
+
+ Int_t GetNodeIndex(const Float_t *p);
+ AliTRDTKDInterpolator::AliTRDTKDNodeInfo* GetNodeInfo(Int_t inode) const;
+ Int_t GetNTNodes() const;
+ void BuildInterpolation();
+ void BuildBoundaryNodes();
+ AliTRDTKDInterpolator(const AliTRDTKDInterpolator &ref);
+ AliTRDTKDInterpolator &operator=(const AliTRDTKDInterpolator &ref);
+
+ Int_t fNDataNodes; // Number of filled nodes (total-zero nodes)
+ TClonesArray *fNodes; //interpolation nodes
+ UChar_t fLambda; // number of parameters in polynom
+ Int_t fNPointsI; // number of points for interpolation
+ Bool_t fUseHelperNodes; // Build Helper nodes to ensure boundary conditions
+ Bool_t fUseWeights; // Use tricubic weights
+ TRDTKDMode fPDFMode; // Mode for PDF calculation
+
+ ClassDef(AliTRDTKDInterpolator, 1) // data interpolator based on KD tree
+};
+
+#endif
#pragma link C++ class AliTRDPIDReference+;
#pragma link C++ class AliTRDPIDParams+;
#pragma link C++ class AliTRDPIDParams::AliTRDPIDThresholds+;
+#pragma link C++ class AliTRDPIDResponseObject+;
+#pragma link C++ class AliTRDTKDInterpolator+;
+#pragma link C++ class AliTRDTKDInterpolator::AliTRDTKDNodeInfo+;
#pragma link C++ class AliITSPidParams+;
#pragma link C++ class AliPIDResponse+;
#pragma link C++ class AliITSPIDResponse+;
Bool_t hasReference = kFALSE;
while ((obj = refs())){
if ((ref = dynamic_cast<AliTRDPIDReference *>(obj))){
- fPIDResponse->Load(ref);
+ //fPIDResponse->Load(ref);
hasReference = kTRUE;
break;
}
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff