From cc1fe362b25017db65870c9f4a5dd2ea20dc9c9c Mon Sep 17 00:00:00 2001 From: schutz Date: Wed, 12 May 2004 15:02:54 +0000 Subject: [PATCH] Bayesian weight for PID calculated for time of flight info --- PHOS/AliPHOSFastRecParticle.cxx | 43 +++++++++------ PHOS/AliPHOSFastRecParticle.h | 3 ++ PHOS/AliPHOSPIDv1.cxx | 96 ++++++++++++++++++++++++++++++--- PHOS/AliPHOSPIDv1.h | 16 ++++-- PHOS/AliPHOSRecParticle.cxx | 36 ++----------- PHOS/AliPHOSRecParticle.h | 5 +- 6 files changed, 138 insertions(+), 61 deletions(-) diff --git a/PHOS/AliPHOSFastRecParticle.cxx b/PHOS/AliPHOSFastRecParticle.cxx index 367fc68c018..f179d7fff42 100644 --- a/PHOS/AliPHOSFastRecParticle.cxx +++ b/PHOS/AliPHOSFastRecParticle.cxx @@ -490,20 +490,31 @@ void AliPHOSFastRecParticle::Print()const { // Print the type, energy and momentum of the reconstructed particle - TString message ; - message = "\n PID bits are %d%d%d %d%d%d %d%d%d %d%d%d" ; - message += ", type is \"%s\"\n" ; - message += " (E,Px,Py,Pz) = (% .3e, % .3e, % .3e, % .3e) GeV\n" ; - Info("Print", message.Data(), - TestPIDBit(0),TestPIDBit(1), - TestPIDBit(2),TestPIDBit(3), - TestPIDBit(4),TestPIDBit(5), - TestPIDBit(6),TestPIDBit(7), - TestPIDBit(8),TestPIDBit(9), - TestPIDBit(10),TestPIDBit(11), - Name().Data(), - Energy(), - Px(), - Py(), - Pz() ); + Info("Print", "-----------------------------") ; + printf("PID bits are %d%d%d %d%d%d %d%d%d %d%d%d", + TestPIDBit(0),TestPIDBit(1), + TestPIDBit(2),TestPIDBit(3), + TestPIDBit(4),TestPIDBit(5), + TestPIDBit(6),TestPIDBit(7), + TestPIDBit(8),TestPIDBit(9), + TestPIDBit(10),TestPIDBit(11)) ; + printf(", type is \"%s\"\n", Name().Data()) ; + printf(" (E,Px,Py,Pz) = (% .3e, % .3e, % .3e, % .3e) GeV\n", + Energy(), + Px(), + Py(), + Pz() ) ; + printf(" TOF = %.3e ns\n", ToF() ) ; + printf(" PID weight: \n" ) ; + printf(" photon -> %f\n", fPID[AliESDtrack::kPhoton] ) ; + printf(" electron -> %f\n", fPID[AliESDtrack::kElectron] ) ; + printf(" Conversion electron -> %f\n", fPID[AliESDtrack::kEleCon] ) ; + printf(" muon -> %f\n", fPID[AliESDtrack::kMuon] ) ; + printf(" neutral pion -> %f\n", fPID[AliESDtrack::kPi0] ) ; + printf(" charged pion -> %f\n", fPID[AliESDtrack::kPion] ) ; + printf(" charged kaon -> %f\n", fPID[AliESDtrack::kKaon] ) ; + printf(" neutral kaon -> %f\n", fPID[AliESDtrack::kKaon0] ) ; + printf(" proton -> %f\n", fPID[AliESDtrack::kProton] ) ; + printf(" neutron -> %f\n", fPID[AliESDtrack::kNeutron] ) ; + } diff --git a/PHOS/AliPHOSFastRecParticle.h b/PHOS/AliPHOSFastRecParticle.h index 6bd9a40196f..954a7ce5358 100644 --- a/PHOS/AliPHOSFastRecParticle.h +++ b/PHOS/AliPHOSFastRecParticle.h @@ -17,6 +17,7 @@ class TClonesArray; #include "TParticle.h" +#include "AliESDtrack.h" // --- Standard library --- @@ -96,6 +97,8 @@ class AliPHOSFastRecParticle : public TParticle { Int_t fIndexInList ; // the index of this RecParticle in the list stored in TreeR (to be set by analysis) Float_t fTof ; // time of fliht Int_t fType ; // particle type obtained by "virtual" reconstruction + Double_t fPID[AliESDtrack::kSPECIESN] ; // PID probability densities + private: ClassDef(AliPHOSFastRecParticle,2) // Reconstructed Particle produced by the fast simulation diff --git a/PHOS/AliPHOSPIDv1.cxx b/PHOS/AliPHOSPIDv1.cxx index 4ed5ee2efb8..da960df4141 100644 --- a/PHOS/AliPHOSPIDv1.cxx +++ b/PHOS/AliPHOSPIDv1.cxx @@ -175,6 +175,32 @@ void AliPHOSPIDv1::InitParameters() fRecParticlesInRun = 0 ; SetParameters() ; // fill the parameters matrix from parameters file SetEventRange(0,-1) ; + // initialisation of response function parameters + // Tof + // Photons + fTphoton[0] = 2.97E-1 ; + fTphoton[1] = 1.55E-8 ; + fTphoton[2] = 5.40E-10 ; + fTFphoton = new TFormula("ToF response to photon" , "gaus") ; + fTFphoton->SetParameters( fTphoton[0], fTphoton[1], fTphoton[2]) ; + // Electrons + fTelectron[0] = 2.73E-1 ; + fTelectron[1] = 1.56E-8 ; + fTelectron[2] = 5.70E-10 ; + fTFelectron = new TFormula("ToF response to electron" , "gaus") ; + fTFelectron->SetParameters( fTelectron[0], fTelectron[1], fTelectron[2]) ; + //Charged Hadrons + fTchargedhadron[0] = 1.58E-1 ; + fTchargedhadron[1] = 1.64E-8 ; + fTchargedhadron[2] = 3.59E-10 ; + fTFchargedhadron = new TFormula("ToF response to charged hadron" , "landau") ; + fTFchargedhadron->SetParameters( fTchargedhadron[0], fTchargedhadron[1], fTchargedhadron[2]) ; + //Neutral Hadrons + fTneutralhadron[0] = 9.62E-1 ; + fTneutralhadron[1] = 1.65E-8 ; + fTneutralhadron[2] = 6.46E-10 ; + fTFneutralhadron = new TFormula("ToF response to neutral hadron" , "landau") ; + fTFneutralhadron->SetParameters( fTneutralhadron[0], fTneutralhadron[1], fTneutralhadron[2]) ; } //________________________________________________________________________ @@ -535,17 +561,71 @@ void AliPHOSPIDv1::MakePID() // construct the PID weight from a Bayesian Method Int_t index ; - Float_t pid1, pid2, pid3, pid4, pid5, pid6 ; - pid1 = pid2 = pid3 = pid4 = pid5 = pid6 = 0 ; + const Int_t kSPECIES = AliESDtrack::kSPECIESN ; + Double_t pid[kSPECIES] = {0., 0., 0., 0., 0., 0.} ; Int_t nparticles = AliPHOSGetter::Instance()->RecParticles()->GetEntriesFast() ; + const Int_t kMAXPARTICLES = 200 ; + if (nparticles >= kMAXPARTICLES) + Error("MakePID", "Change size of MAXPARTICLES") ; + Double_t stof[kSPECIES][kMAXPARTICLES] ; + // make the normalized distribution of pid for this event + // w(pid) in the Bayesian formulation for(index = 0 ; index < nparticles ; index ++) { AliPHOSRecParticle * recpar = AliPHOSGetter::Instance()->RecParticle(index) ; - if (recpar->IsPhoton() || recpar->IsHardPhoton()) pid1++ ; - else if (recpar->IsPi0() || recpar->IsHardPi0()) pid2++ ; - else if (recpar->IsElectron()) pid3++ ; - else if (recpar->IsChargedHadron()) pid4++ ; - else if (recpar->IsNeutralHadron()) pid5++ ; - else if (recpar->IsEleCon()) pid6++ ; + + pid[AliESDtrack::kKaon0] = 0 ; + + if (recpar->IsPhoton() || recpar->IsHardPhoton()) + pid[AliESDtrack::kPhoton]++ ; + + else if (recpar->IsPi0() || recpar->IsHardPi0()) + pid[AliESDtrack::kPi0]++ ; + + else if (recpar->IsElectron()) { + pid[AliESDtrack::kElectron]++ ; + pid[AliESDtrack::kMuon]++ ; + } + + else if (recpar->IsChargedHadron()){ + pid[AliESDtrack::kPion]++ ; + pid[AliESDtrack::kKaon]++ ; + pid[AliESDtrack::kProton]++ ; + } + + else if (recpar->IsNeutralHadron()) + pid[AliESDtrack::kNeutron]++ ; + + else if (recpar->IsEleCon()) + pid[AliESDtrack::kEleCon]++ ; + + // now get the signals probability + // s(pid) in the Bayesian formulation + // Tof + stof[AliESDtrack::kPhoton][index] = fTFphoton->Eval(recpar->ToF()) ; + stof[AliESDtrack::kPi0][index] = fTFphoton->Eval(recpar->ToF()) ; // pi0 are detected via decay photon + stof[AliESDtrack::kElectron][index] = fTFelectron->Eval(recpar->ToF()) ; + stof[AliESDtrack::kPion][index] = fTFchargedhadron->Eval(recpar->ToF()) ; + stof[AliESDtrack::kKaon][index] = fTFchargedhadron->Eval(recpar->ToF()) ; + stof[AliESDtrack::kProton][index] = fTFchargedhadron->Eval(recpar->ToF()) ; + stof[AliESDtrack::kNeutron][index] = fTFneutralhadron->Eval(recpar->ToF()) ; + stof[AliESDtrack::kEleCon][index] = fTFphoton->Eval(recpar->ToF()) ; // a conversion electron has the photon ToF + stof[AliESDtrack::kKaon0][index] = 0 ; // do not know yet what to to with K0 + + } + for (index = 0 ; index < kSPECIES ; index++) + pid[index] /= nparticles ; + + + for(index = 0 ; index < nparticles ; index ++) { + // calculates the Bayesian weight + Int_t jndex ; + Double_t wn = 0.0 ; + for (jndex = 0 ; jndex < kSPECIES ; jndex++) + wn += stof[jndex][index] * pid[jndex] ; + AliPHOSRecParticle * recpar = AliPHOSGetter::Instance()->RecParticle(index) ; + for (jndex = 0 ; jndex < kSPECIES ; jndex++) { + recpar->SetPID(jndex, stof[jndex][index] * pid[jndex] / wn) ; + } } } diff --git a/PHOS/AliPHOSPIDv1.h b/PHOS/AliPHOSPIDv1.h index a04a35dea0f..71fade12814 100644 --- a/PHOS/AliPHOSPIDv1.h +++ b/PHOS/AliPHOSPIDv1.h @@ -106,9 +106,19 @@ private: Double_t *fPPi0 ; //! Principal pi0 eigenvalues Int_t fRecParticlesInRun ; //! Total number of recparticles in one run TMatrix *fParameters; //! Matrix of identification Parameters - - - ClassDef( AliPHOSPIDv1,9) // Particle identifier implementation version 1 + // response function parameters + // ToF + Double_t fTphoton[3] ; // gaussian response for photon + TFormula * fTFphoton ; // the formula + Double_t fTelectron[3] ; // gaussian response for electrons + TFormula * fTFelectron ; // the formula + Double_t fTchargedhadron[3] ; // landau response for charged hadrons + TFormula * fTFchargedhadron ; // the formula + Double_t fTneutralhadron[3] ; // landau response for neutral hadrons + TFormula * fTFneutralhadron ; // the formula + + + ClassDef( AliPHOSPIDv1,10) // Particle identifier implementation version 1 }; diff --git a/PHOS/AliPHOSRecParticle.cxx b/PHOS/AliPHOSRecParticle.cxx index 7f13cb64ffe..60dec667dde 100644 --- a/PHOS/AliPHOSRecParticle.cxx +++ b/PHOS/AliPHOSRecParticle.cxx @@ -161,9 +161,9 @@ const TParticle * AliPHOSRecParticle::GetPrimary(Int_t index) const } //____________________________________________________________________________ -const Double_t * AliPHOSRecParticle::GetPID() +void AliPHOSRecParticle::SetPID(Int_t type, Double_t weight) { - // Get the probability densities that this reconstructed particle + // Set the probability densities that this reconstructed particle // has a type of i: // i particle types // ---------------------- @@ -176,33 +176,7 @@ const Double_t * AliPHOSRecParticle::GetPID() // 6 pi0 at high pt // 7 neutron // 8 K0L - - - if (IsElectron() ) fPID[0] = 1.0; - if (IsChargedHadron()) { - fPID[1] = 0.25; - fPID[2] = 0.25; - fPID[3] = 0.25; - fPID[4] = 0.25; - } - if (IsFastChargedHadron()) { - fPID[1] = 0.33; - fPID[2] = 0.33; - fPID[3] = 0.33; - fPID[4] = 0.00; - } - if (IsSlowChargedHadron()) { - fPID[1] = 0.00; - fPID[2] = 0.00; - fPID[3] = 0.00; - fPID[4] = 1.00; - } - - if (IsPhoton() || IsHardPhoton()) fPID[5] = 1.0; - if (IsHardPi0()) fPID[6] = 1.0; - if (IsFastNeutralHadron()) fPID[7] = 1.0; - if (IsSlowNeutralHadron()) fPID[8] = 1.0; - - if (IsEleCon()) fPID[9] = 1.0; - return fPID; + // 9 Conversion electron + + fPID[type] = weight ; } diff --git a/PHOS/AliPHOSRecParticle.h b/PHOS/AliPHOSRecParticle.h index 32195b88d66..1d603ce3aab 100644 --- a/PHOS/AliPHOSRecParticle.h +++ b/PHOS/AliPHOSRecParticle.h @@ -18,7 +18,6 @@ // --- AliRoot header files --- #include "AliPHOSFastRecParticle.h" -#include "AliESDtrack.h" class TParticle ; #include "TVector3.h" @@ -37,8 +36,9 @@ class AliPHOSRecParticle : public AliPHOSFastRecParticle { TVector3 GetPos() const { return fPos ; } virtual const TParticle * GetPrimary(Int_t index) const ; virtual const TParticle * GetPrimary() const ; - const Double_t *GetPID(); + const Double_t *GetPID() { return fPID ; } void SetDebug() { fDebug = kTRUE ; } + void SetPID(Int_t type, Double_t weight) ; void SetPos(TVector3 pos) { fPos.SetXYZ( pos.X(), pos.Y(), pos.Z() ); } void UnsetDebug() { fDebug = kFALSE ; } void SetTrackSegment(Int_t index){fPHOSTrackSegment = index; } @@ -50,7 +50,6 @@ class AliPHOSRecParticle : public AliPHOSFastRecParticle { Int_t fPHOSTrackSegment ; // pointer to the associated track segment in PHOS Bool_t fDebug ; // to steer debug output TVector3 fPos ; // position in the global alice coordinate system - Double_t fPID[AliESDtrack::kSPECIESN] ; // PID probability densities ClassDef(AliPHOSRecParticle,3) // Reconstructed Particle }; -- 2.39.3