From db0e2c5f7ecb63300fd90881a5a3969dfc91d2aa Mon Sep 17 00:00:00 2001 From: morsch Date: Wed, 27 Jun 2012 15:21:03 +0000 Subject: [PATCH] Updates in order to enable the '2D' PID for the TRD developed by Daniel Lohner. --- .../TenderSupplies/AliTRDTenderSupply.cxx | 10 +- STEER/CMakelibSTEERBase.pkg | 4 +- STEER/STEERBase/AliPIDResponse.cxx | 50 +- STEER/STEERBase/AliPIDResponse.h | 4 +- STEER/STEERBase/AliTRDPIDParams.cxx | 24 + STEER/STEERBase/AliTRDPIDParams.h | 40 +- STEER/STEERBase/AliTRDPIDReference.cxx | 26 +- STEER/STEERBase/AliTRDPIDResponse.cxx | 225 ++++--- STEER/STEERBase/AliTRDPIDResponse.h | 23 +- STEER/STEERBase/AliTRDPIDResponseObject.cxx | 219 +++++++ STEER/STEERBase/AliTRDPIDResponseObject.h | 72 +++ STEER/STEERBase/AliTRDTKDInterpolator.cxx | 588 ++++++++++++++++++ STEER/STEERBase/AliTRDTKDInterpolator.h | 104 ++++ STEER/STEERBaseLinkDef.h | 3 + TRD/AliTRDcalibDB.cxx | 2 +- 15 files changed, 1230 insertions(+), 164 deletions(-) create mode 100644 STEER/STEERBase/AliTRDPIDResponseObject.cxx create mode 100644 STEER/STEERBase/AliTRDPIDResponseObject.h create mode 100644 STEER/STEERBase/AliTRDTKDInterpolator.cxx create mode 100644 STEER/STEERBase/AliTRDTKDInterpolator.h diff --git a/ANALYSIS/TenderSupplies/AliTRDTenderSupply.cxx b/ANALYSIS/TenderSupplies/AliTRDTenderSupply.cxx index d4b3763cd8e..0f798f71be0 100644 --- a/ANALYSIS/TenderSupplies/AliTRDTenderSupply.cxx +++ b/ANALYSIS/TenderSupplies/AliTRDTenderSupply.cxx @@ -47,7 +47,7 @@ #include #include #include -#include +#include #include #include "AliTRDCalChamberStatus.h" #include @@ -276,15 +276,15 @@ void AliTRDTenderSupply::LoadReferences(){ // 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(o); + if(!TString(o->IsA()->GetName()).CompareTo("AliTRDPIDResponseObject")){ + ref = dynamic_cast(o); break; } } if(ref){ - fESDpid->GetTRDResponse().Load(ref); + fESDpid->GetTRDResponse().SetPIDResponseObject(ref); fHasReferences = kTRUE; AliDebug(1, "Reference distributions loaded into the PID Response"); } else { diff --git a/STEER/CMakelibSTEERBase.pkg b/STEER/CMakelibSTEERBase.pkg index d7916e01046..ef4d53ded76 100644 --- a/STEER/CMakelibSTEERBase.pkg +++ b/STEER/CMakelibSTEERBase.pkg @@ -62,11 +62,13 @@ set ( SRCS 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 diff --git a/STEER/STEERBase/AliPIDResponse.cxx b/STEER/STEERBase/AliPIDResponse.cxx index 511475dccb9..57ed3ff64de 100644 --- a/STEER/STEERBase/AliPIDResponse.cxx +++ b/STEER/STEERBase/AliPIDResponse.cxx @@ -29,14 +29,14 @@ #include #include #include +#include #include #include #include #include #include -#include -#include +#include #include #include "AliPIDResponse.h" @@ -67,8 +67,7 @@ fRun(0), fOldRun(0), fArrPidResponseMaster(0x0), fResolutionCorrection(0x0), -fTRDPIDParams(0x0), -fTRDPIDReference(0x0), +fTRDPIDResponseObject(0x0), fTOFtail(1.1), fTOFPIDParams(0x0), fEMCALPIDParams(0x0), @@ -91,9 +90,8 @@ AliPIDResponse::~AliPIDResponse() // // dtor // - delete fArrPidResponseMaster; - delete fTRDPIDParams; - delete fTRDPIDReference; + delete fArrPidResponseMaster; + delete fTRDPIDResponseObject; if (fTOFPIDParams) delete fTOFPIDParams; } @@ -120,8 +118,7 @@ fRun(0), fOldRun(0), fArrPidResponseMaster(0x0), fResolutionCorrection(0x0), -fTRDPIDParams(0x0), -fTRDPIDReference(0x0), +fTRDPIDResponseObject(0x0), fTOFtail(1.1), fTOFPIDParams(0x0), fEMCALPIDParams(0x0), @@ -163,8 +160,7 @@ AliPIDResponse& AliPIDResponse::operator=(const AliPIDResponse &other) fOldRun=0; fArrPidResponseMaster=0x0; fResolutionCorrection=0x0; - fTRDPIDParams=0x0; - fTRDPIDReference=0x0; + fTRDPIDResponseObject=0x0; fEMCALPIDParams=0x0; memset(fTRDslicesForPID,0,sizeof(UInt_t)*2); fTOFtail=1.1; @@ -846,20 +842,20 @@ void AliPIDResponse::SetTRDPidResponseMaster() // // 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(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(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){ @@ -870,7 +866,8 @@ void AliPIDResponse::SetTRDPidResponseMaster() if(!fTRDPIDReference){ AliError(Form("TRD References not found in OADB Container for run %d", fRun)); } - } + } + */ } //______________________________________________________________________________ @@ -878,12 +875,11 @@ void AliPIDResponse::InitializeTRDResponse(){ // // 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; + } } //______________________________________________________________________________ diff --git a/STEER/STEERBase/AliPIDResponse.h b/STEER/STEERBase/AliPIDResponse.h index 6f11e9f646a..d7157305643 100644 --- a/STEER/STEERBase/AliPIDResponse.h +++ b/STEER/STEERBase/AliPIDResponse.h @@ -25,6 +25,7 @@ class AliVEvent; class TF1; +class AliTRDPIDResponseObject; class AliPIDResponse : public TNamed { public: @@ -130,8 +131,7 @@ private: 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 diff --git a/STEER/STEERBase/AliTRDPIDParams.cxx b/STEER/STEERBase/AliTRDPIDParams.cxx index 5f8e60c9740..950ac424f53 100644 --- a/STEER/STEERBase/AliTRDPIDParams.cxx +++ b/STEER/STEERBase/AliTRDPIDParams.cxx @@ -49,6 +49,30 @@ AliTRDPIDParams::AliTRDPIDParams(const char *name) : 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(){ // diff --git a/STEER/STEERBase/AliTRDPIDParams.h b/STEER/STEERBase/AliTRDPIDParams.h index bbadcae79bd..0d46a9f3e38 100644 --- a/STEER/STEERBase/AliTRDPIDParams.h +++ b/STEER/STEERBase/AliTRDPIDParams.h @@ -28,6 +28,7 @@ class AliTRDPIDParams : public TNamed{ public: AliTRDPIDParams(); AliTRDPIDParams(const char *name); + AliTRDPIDParams(const AliTRDPIDParams &); virtual ~AliTRDPIDParams(); virtual void Print(Option_t *) const; @@ -36,35 +37,34 @@ class AliTRDPIDParams : public TNamed{ 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); diff --git a/STEER/STEERBase/AliTRDPIDReference.cxx b/STEER/STEERBase/AliTRDPIDReference.cxx index e3e6379f918..c50707d5fa2 100644 --- a/STEER/STEERBase/AliTRDPIDReference.cxx +++ b/STEER/STEERBase/AliTRDPIDReference.cxx @@ -58,14 +58,26 @@ AliTRDPIDReference::AliTRDPIDReference(const Char_t *name): //____________________________________________________________ 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); } //____________________________________________________________ @@ -124,7 +136,7 @@ void AliTRDPIDReference::AddReference(TObject *ref, AliPID::EParticleType spec, 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); } //____________________________________________________________ diff --git a/STEER/STEERBase/AliTRDPIDResponse.cxx b/STEER/STEERBase/AliTRDPIDResponse.cxx index edff1a18725..fe69bca7258 100644 --- a/STEER/STEERBase/AliTRDPIDResponse.cxx +++ b/STEER/STEERBase/AliTRDPIDResponse.cxx @@ -36,17 +36,18 @@ #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) { @@ -58,8 +59,7 @@ AliTRDPIDResponse::AliTRDPIDResponse(): //____________________________________________________________ AliTRDPIDResponse::AliTRDPIDResponse(const AliTRDPIDResponse &ref): TObject(ref) - ,fkPIDReference(ref.fkPIDReference) - ,fkPIDParams(ref.fkPIDParams) + ,fkPIDResponseObject(NULL) ,fGainNormalisationFactor(ref.fGainNormalisationFactor) ,fPIDmethod(ref.fPIDmethod) { @@ -78,8 +78,7 @@ AliTRDPIDResponse &AliTRDPIDResponse::operator=(const AliTRDPIDResponse &ref){ // Make copy TObject::operator=(ref); fGainNormalisationFactor = ref.fGainNormalisationFactor; - fkPIDReference = ref.fkPIDReference; - fkPIDParams = ref.fkPIDParams; + fkPIDResponseObject = ref.fkPIDResponseObject; fPIDmethod = ref.fPIDmethod; return *this; @@ -90,11 +89,11 @@ AliTRDPIDResponse::~AliTRDPIDResponse(){ // // 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 // @@ -110,11 +109,11 @@ Bool_t AliTRDPIDResponse::Load(const Char_t * filename, const Char_t *refName){ } gROOT->cd(); - fkPIDReference = dynamic_cast(in->Get(refName)->Clone()); + fkPIDResponseObject = dynamic_cast(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; } @@ -139,73 +138,110 @@ Bool_t AliTRDPIDResponse::GetResponse(Int_t n, const Double_t * const dedx, cons // 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(fkPIDReference->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper)), - *refLower = dynamic_cast(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(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(fkPIDResponseObject->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper,kLQ1D)), + *refLower = dynamic_cast(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; } @@ -214,10 +250,10 @@ void AliTRDPIDResponse::SetOwner(){ // // Make Deep Copy of the Reference Histograms // - if(!fkPIDReference || IsOwner()) return; - const AliTRDPIDReference *tmp = fkPIDReference; - fkPIDReference = dynamic_cast(tmp->Clone()); - SetBit(kIsOwner, kTRUE); + if(!fkPIDResponseObject || IsOwner()) return; + const AliTRDPIDResponseObject *tmp = fkPIDResponseObject; + fkPIDResponseObject = dynamic_cast(tmp->Clone()); + SetBit(kIsOwner, kTRUE); } //____________________________________________________________ @@ -228,42 +264,49 @@ Bool_t AliTRDPIDResponse::CookdEdx(Int_t nSlice, const Double_t * const in, Doub // 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(isliceGetNSlicesQ0())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; } diff --git a/STEER/STEERBase/AliTRDPIDResponse.h b/STEER/STEERBase/AliTRDPIDResponse.h index 53d90de4a7e..f1c079ba566 100644 --- a/STEER/STEERBase/AliTRDPIDResponse.h +++ b/STEER/STEERBase/AliTRDPIDResponse.h @@ -30,22 +30,25 @@ 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, @@ -65,19 +68,18 @@ class AliTRDPIDResponse : public TObject { 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 @@ -96,3 +98,4 @@ AliTRDPIDResponse::ETRDNslices AliTRDPIDResponse::GetNumberOfSlices() const { return slices; } #endif + diff --git a/STEER/STEERBase/AliTRDPIDResponseObject.cxx b/STEER/STEERBase/AliTRDPIDResponseObject.cxx new file mode 100644 index 00000000000..566fa908892 --- /dev/null +++ b/STEER/STEERBase/AliTRDPIDResponseObject.cxx @@ -0,0 +1,219 @@ +/************************************************************************** + * 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 +// Daniel Lohner + +#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=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;methodPrint(opt); + if(fkPIDParams[method])printf("+ Threshold Parameters \n"); + } +} diff --git a/STEER/STEERBase/AliTRDPIDResponseObject.h b/STEER/STEERBase/AliTRDPIDResponseObject.h new file mode 100644 index 00000000000..701dc5a9984 --- /dev/null +++ b/STEER/STEERBase/AliTRDPIDResponseObject.h @@ -0,0 +1,72 @@ +/************************************************************************** +* 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 +#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 diff --git a/STEER/STEERBase/AliTRDTKDInterpolator.cxx b/STEER/STEERBase/AliTRDTKDInterpolator.cxx new file mode 100644 index 00000000000..d2842ed9cb1 --- /dev/null +++ b/STEER/STEERBase/AliTRDTKDInterpolator.cxx @@ -0,0 +1,588 @@ +#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; infBounds,GetBoundary(tnode),2*fNDim*sizeof(Float_t)); + node->fVal[0] = Float_t(fBucketSize)/fNPoints; + for(int idim=0; idimfVal[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; idimfData[idim] = 0.; + for(int ip = 0; ipfData[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; idimfBounds[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; idimfData[idim] = 0.; + for(int ip = 0; ipfData[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; idimCookPDF(point, result, error); + case kMinError: + node->CookPDF(point, result, error); + if(errorfVal[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;inodeHas(p)){return inode;} + } + + // Search for nearest neighbor + Float_t dist; + Float_t closestdist=10000; + inode=-1; + for(Int_t ii=0;iifData[idim]-p[idim]),2); + } + dist=TMath::Sqrt(dist); + if(dist= 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]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;inodeAddBin(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; idfData[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;infData[idim]; + for(int jdim=idim; jdimfData[idim]*nnode->fData[jdim]; + } + fitter.AddPoint(buffer,nnode->fVal[0], nnode->fVal[1]/w); + + // Ensure Boundary Condition + for(Int_t kdim=0;kdimfBounds[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; idimStore(&par,&cov); + } + + delete KDhelper; + for(int id=0; idfBounds[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;idimfBounds[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; idimfVal[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; idim0)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=fBounds[2*id] && p[id](obj))){ - fPIDResponse->Load(ref); + //fPIDResponse->Load(ref); hasReference = kTRUE; break; } -- 2.18.4