Pid optimized with new strategy + minors

diff --git a/HMPID/AliHMPIDPid.cxx b/HMPID/AliHMPIDPid.cxx
index 2d1edad64a1eb8134c3996d345761bfb866ee279..25f8af65e72cebde5e31c1da55031144a0ae95eb 100644 (file)
--- a/HMPID/AliHMPIDPid.cxx
+++ b/HMPID/AliHMPIDPid.cxx
@@ -23,7 +23,9 @@
 
 #include "AliHMPIDPid.h"       //class header
 #include "AliHMPIDParam.h"     //class header
+#include "AliHMPIDRecon.h"     //class header
 #include <AliESDtrack.h>       //FindPid()
+#include <TRandom.h>           //Resolution()
 
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 AliHMPIDPid::AliHMPIDPid():TTask("HMPIDrec","HMPIDPid")
@@ -35,35 +37,83 @@ AliHMPIDPid::AliHMPIDPid():TTask("HMPIDrec","HMPIDPid")
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 void AliHMPIDPid::FindPid(AliESDtrack *pTrk,Int_t nsp,Double_t *prob)
 {
-// Calculates probability to be a electron-muon-pion-kaon-proton
+// Calculates probability to be a electron-muon-pion-kaon-proton with the "amplitude" method
 // from the given Cerenkov angle and momentum assuming no initial particle composition
 // (i.e. apriory probability to be the particle of the given sort is the same for all sorts)
 
-  if(pTrk->GetHMPIDsignal()<=0){//HMPID does not find anything reasonable for this track, assign 0.2 for all species
+  AliPID *pPid = new AliPID();
+  Double_t thetaCerExp = pTrk->GetHMPIDsignal();                                                                           //  measured thetaCherenkov
+  
+  if(thetaCerExp<=0){                                         //HMPID does not find anything reasonable for this track, assign 0.2 for all species
     for(Int_t iPart=0;iPart<nsp;iPart++) prob[iPart]=1.0/(Float_t)nsp;
     return;
   } 
 
-  Double_t pmod = pTrk->GetP();
-  Double_t hTot=0;
-  Double_t *h = new Double_t [nsp];
+  Double_t pmod = pTrk->GetP();                  // Momentum of the charged particle
+  Double_t hTot=0;                               // Initialize the total height of the amplitude method
+  Double_t *h = new Double_t [nsp];              // number of charged particles to be considered
 
-  for(Int_t iPart=0;iPart<nsp;iPart++){
-    Double_t mass = AliPID::ParticleMass(iPart);
-    Double_t cosThetaTh = TMath::Sqrt(mass*mass+pmod*pmod)/(AliHMPIDParam::Instance()->MeanIdxRad()*pmod);
-    if(cosThetaTh<1) //calculate the height of theoretical theta ckov on the gaus of experimental one
-      h[iPart] =TMath::Gaus(TMath::ACos(cosThetaTh),pTrk->GetHMPIDsignal(),TMath::Sqrt(pTrk->GetHMPIDchi2()),kTRUE);      
-    else             //beta < 1/ref. idx. => no light at all  
-      h[iPart] =0 ;       
+  Bool_t desert = kTRUE;                                                                                                     //  Flag to evaluate if ThetaC is far ("desert") from the given Gaussians
+  
+  for(Int_t iPart=0;iPart<nsp;iPart++){                                                                                      //  for each particle
+    
+    h[iPart] = 0;                                                                                                            //  reset the height
+    Double_t mass = pPid->ParticleMass(iPart);                                                                             //  with the given mass
+    Double_t cosThetaTh = TMath::Sqrt(mass*mass+pmod*pmod)/(AliHMPIDParam::Instance()->MeanIdxRad()*pmod);                   //  evaluate the theor. Theta Cherenkov
+    if(cosThetaTh>1) continue;                                                                                               //  no light emitted, zero height
+    Double_t thetaCerTh = TMath::ACos(cosThetaTh);                                                                           //  theoretical Theta Cherenkov
+    Double_t sigmaRing = Resolution(thetaCerTh,pTrk);
+    
+    if(sigmaRing==0) {
+      for(Int_t jPart=0;jPart<nsp;jPart++) prob[jPart]=1.0/(Float_t)nsp;
+      return;
+    } 
+    
+    if(TMath::Abs(thetaCerExp-thetaCerTh)<4*sigmaRing) desert = kFALSE;                                                                //   
+    h[iPart] =TMath::Gaus(thetaCerTh,thetaCerExp,sigmaRing,kTRUE);
     hTot    +=h[iPart]; //total height of all theoretical heights for normalization
+    
   }//species loop
 
-  Double_t hMin=TMath::Gaus(pTrk->GetHMPIDsignal()-4*TMath::Sqrt(pTrk->GetHMPIDchi2()),pTrk->GetHMPIDsignal(),TMath::Sqrt(pTrk->GetHMPIDchi2()),kTRUE);//5 sigma protection
-
   for(Int_t iPart=0;iPart<nsp;iPart++) {//species loop to assign probabilities
-    if(hTot>hMin) prob[iPart]=h[iPart]/hTot;
+    
+    if(!desert) prob[iPart]=h[iPart]/hTot;
     else prob[iPart]=1.0/(Float_t)nsp;            //all theoretical values are far away from experemental one
+    
   }
+  
   delete [] h;
+  delete pPid ; pPid=0x0;
 }
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+Double_t AliHMPIDPid::Resolution(Double_t thetaCerTh, AliESDtrack *pTrk)
+{
+  AliHMPIDParam *pParam = AliHMPIDParam::Instance();
+      
+  AliHMPIDRecon rec;
+  Float_t xRa,yRa,thRa,phRa;
+  pTrk->GetHMPIDtrk(xRa,yRa,thRa,phRa);
+  rec.SetTrack(xRa,yRa,thRa,phRa);
+  Double_t thetaMax = TMath::ACos(1./pParam->MeanIdxRad());
+  Int_t nPhots = (Int_t)(21.*TMath::Sin(thetaCerTh)*TMath::Sin(thetaCerTh)/(TMath::Sin(thetaMax)*TMath::Sin(thetaMax))+0.01);
+
+  Double_t sigmatot = 0;
+  Int_t nTrks = 20;
+  for(Int_t iTrk=0;iTrk<nTrks;iTrk++) {
+    Double_t invSigma = 0;
+    Int_t nPhotsAcc = 0;
+    for(Int_t j=0;j<nPhots;j++){
+      Double_t phi = gRandom->Rndm()*TMath::TwoPi();
+      TVector2 pos; pos=rec.TracePhot(thetaCerTh,phi);
+      if(!pParam->IsInside(pos.X(),pos.Y())) continue;
+      if(pParam->IsInDead(pos.X(),pos.Y())) continue;
+      Double_t sigma2 = pParam->Sigma2(thRa,phRa,thetaCerTh,phi);//photon candidate sigma^2
+      if(sigma2!=0) {
+        invSigma += 1./sigma2;
+        nPhotsAcc++;
+      }
+    }      
+    if(invSigma!=0) sigmatot += 1./TMath::Sqrt(invSigma);  
+  }
+  return sigmatot/nTrks;
+}

diff --git a/HMPID/AliHMPIDPid.h b/HMPID/AliHMPIDPid.h
index 7b48b598a3f04983b4b529c63aaabacc705b6a0e..1dd2942f8d05156c189f41a3a69de92cd6738090 100644 (file)
--- a/HMPID/AliHMPIDPid.h
+++ b/HMPID/AliHMPIDPid.h
@@ -23,7 +23,8 @@ public :
              AliHMPIDPid();    //ctor
     virtual ~AliHMPIDPid() {;} //dtor
     
-    void FindPid(AliESDtrack *pESD,Int_t nsp,Double_t *prob);  //Find PID for tracks
+    void     FindPid(AliESDtrack *pESD,Int_t nsp,Double_t *prob);  //Find PID for tracks
+    Double_t Resolution(Double_t thetaCerTh, AliESDtrack *pTrk);   //Find the sigma for a given ThetaCerTh
 
 //
 protected:

diff --git a/HMPID/HESDfromKin.C b/HMPID/HESDfromKin.C
index ea63651c7ffd715ea67cc4cc5bf8f0f5c5d2490f..ec1377a5419a5824549d47d1cbd260649aaf229a 100644 (file)
--- a/HMPID/HESDfromKin.C
+++ b/HMPID/HESDfromKin.C
@@ -68,10 +68,23 @@ void SimEsd(AliLoader *pHL,AliESDEvent *pEsd)
       trk.SetHMPIDcluIdx   (iCh,99999);                                                          //chamber not found, mip not yet considered
       trk.SetHMPIDtrk(xRa,yRa,thRa,phRa);                                                        //store initial infos
       pEsd->AddTrack(&trk);
-      AliHMPIDTracker::Recon(pEsd,pH->CluLst(),pNmean,pQthre);
     }// track loop
+    
+    AliHMPIDTracker::Recon(pEsd,pH->CluLst(),pNmean,pQthre);
+// Perform PID
+    
+    for(Int_t iTrk=0;iTrk<pEsd->GetNumberOfTracks();iTrk++){                                       //ESD tracks loop
+      AliESDtrack *pTrk = pEsd->GetTrack(iTrk);                                                    //get reconstructed track    
+      AliHMPIDPid pID;
+      Double_t prob[5];
+      pID.FindPid(pTrk,5,prob);
+      trk.SetHMPIDpid(prob);
+//      Printf(" Prob e- %6.2f mu %6.2f pi %6.2f k %6.2f p %6.2f",prob[0]*100,prob[1]*100,prob[2]*100,prob[3]*100,prob[4]*100);
+    }
+    
     gEsdTr->Fill();
     pEsd->Reset();
+    
   }// event loop
   Printf("Events processed %i",iEvt);
   gAL->UnloadHeader();  gAL->UnloadKinematics();

diff --git a/HMPID/Hdisp.C b/HMPID/Hdisp.C
index 3453ea7afd0e283372db5968154b369bd27b3350..9b39b391bdfb43fc4c85369d679366d8f8018d19 100644 (file)
--- a/HMPID/Hdisp.C
+++ b/HMPID/Hdisp.C
@@ -578,7 +578,7 @@ void DisplayInfo(Int_t evt,Int_t px, Int_t py)
     pTrk->GetHMPIDpid(prob);
     if(ckov>0){
       Float_t x,y;Int_t q,nacc;   pTrk->GetHMPIDmip(x,y,q,nacc);
-      text3="";text3.Append(Form("Theta Cherenkov %5.3f with %i photons |  e %5.1f%% | u %5.1f%% | K %5.1f%% | pi %5.1f%% | p %5.1f%%",
+      text3="";text3.Append(Form("Theta Cherenkov %5.3f with %i photons |  e %5.1f%% | u %5.1f%% | pi %5.1f%% | K %5.1f%% | p %5.1f%%",
           pTrk->GetHMPIDsignal(),nacc,prob[0]*100,prob[1]*100,prob[2]*100,prob[3]*100,prob[4]*100));
     }      
   }