{
// 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] ) ;
+
}
class TClonesArray;
#include "TParticle.h"
+#include "AliESDtrack.h"
// --- Standard library ---
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
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]) ;
}
//________________________________________________________________________
// 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) ;
+ }
}
}
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
};
}
//____________________________________________________________________________
-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
// ----------------------
// 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 ;
}
// --- AliRoot header files ---
#include "AliPHOSFastRecParticle.h"
-#include "AliESDtrack.h"
class TParticle ;
#include "TVector3.h"
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; }
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
};
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff