From ac66a50d01b66323854c535ca369db14398ea5a8 Mon Sep 17 00:00:00 2001 From: dibari Date: Thu, 10 Apr 2008 14:42:58 +0000 Subject: [PATCH] Pid optimized with new strategy + minors --- HMPID/AliHMPIDPid.cxx | 80 +++++++++++++++++++++++++++++++++++-------- HMPID/AliHMPIDPid.h | 3 +- HMPID/HESDfromKin.C | 15 +++++++- HMPID/Hdisp.C | 2 +- 4 files changed, 82 insertions(+), 18 deletions(-) diff --git a/HMPID/AliHMPIDPid.cxx b/HMPID/AliHMPIDPid.cxx index 2d1edad64a1..25f8af65e72 100644 --- 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 //FindPid() +#include //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;iPartGetP(); - 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;iPartMeanIdxRad()*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;iPartParticleMass(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;jPartGetHMPIDsignal()-4*TMath::Sqrt(pTrk->GetHMPIDchi2()),pTrk->GetHMPIDsignal(),TMath::Sqrt(pTrk->GetHMPIDchi2()),kTRUE);//5 sigma protection - for(Int_t iPart=0;iParthMin) 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;iTrkRndm()*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 7b48b598a3f..1dd2942f8d0 100644 --- 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 ea63651c7ff..ec1377a5419 100644 --- 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;iTrkGetNumberOfTracks();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 3453ea7afd0..9b39b391bdf 100644 --- 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)); } } -- 2.43.0