Bayesian weight for PID calculated for time of flight info
authorschutz <schutz@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 12 May 2004 15:02:54 +0000 (15:02 +0000)
committerschutz <schutz@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 12 May 2004 15:02:54 +0000 (15:02 +0000)
PHOS/AliPHOSFastRecParticle.cxx
PHOS/AliPHOSFastRecParticle.h
PHOS/AliPHOSPIDv1.cxx
PHOS/AliPHOSPIDv1.h
PHOS/AliPHOSRecParticle.cxx
PHOS/AliPHOSRecParticle.h

index 367fc68..f179d7f 100644 (file)
@@ -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] ) ; 
+
 }
index 6bd9a40..954a7ce 100644 (file)
@@ -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 
index 4ed5ee2..da960df 100644 (file)
@@ -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) ; 
+    }
   }
 }
 
index a04a35d..71fade1 100644 (file)
@@ -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
 
 };
 
index 7f13cb6..60dec66 100644 (file)
@@ -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 ; 
 }
index 32195b8..1d603ce 100644 (file)
@@ -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
 };