index 5740d4e..aa558c0 100644 (file)
@@ -23,7 +23,9 @@

#include <AliESDtrack.h>       //FindPid()
+#include <TRandom.h>           //Resolution()

//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
//..
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-void AliHMPIDPid::FindPid(AliESDtrack *pTrk,Double_t *prob)
+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)

-  AliPID ppp; //needed
-  Double_t h[AliPID::kSPECIES];
+  AliPID *pPid = new AliPID();
+  Double_t thetaCerExp = pTrk->GetHMPIDsignal();                                                                           //  measured thetaCherenkov

-  if(pTrk->GetHMPIDsignal()<=0){//HMPID does not find anything reasonable for this track, assign 0.2 for all species
-    for(Int_t iPart=0;iPart<AliPID::kSPECIES;iPart++) prob[iPart]=1.0/AliPID::kSPECIES;
-    return;
+  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;
+    delete pPid ; pPid=0x0; return;
}
+
+  Double_t p[3] = {0}, pmod = 0;
+  if(pTrk->GetOuterHmpPxPyPz(p))  pmod = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);  // Momentum of the charged particle
+
+  else {
+    for(Int_t iPart=0;iPart<nsp;iPart++) prob[iPart]=1.0/(Float_t)nsp;
+    delete pPid ; pPid=0x0; return;
+  }
+
+  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

-  Double_t pmod = pTrk->GetP();
-  Double_t hTot=0;
-
-  for(Int_t iPart=0;iPart<AliPID::kSPECIES;iPart++){
-    Double_t mass = AliPID::ParticleMass(iPart);
-    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;
+      delete pPid ; pPid=0x0; delete [] h; 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<AliPID::kSPECIES;iPart++) {//species loop to assign probabilities
-    if(hTot>hMin) prob[iPart]=h[iPart]/hTot;
-    else prob[iPart]=1.0/AliPID::kSPECIES;            //all theoretical values are far away from experemental one
+  for(Int_t iPart=0;iPart<nsp;iPart++) {//species loop to assign probabilities
+
+    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);
+  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;