Updating angular resolution calculation and correction factor as a function of chambe...
authorgvolpe <gvolpe@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 23 Nov 2012 16:21:14 +0000 (16:21 +0000)
committergvolpe <gvolpe@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 23 Nov 2012 16:21:14 +0000 (16:21 +0000)
HMPID/AliHMPIDParam.cxx
HMPID/AliHMPIDParam.h
HMPID/AliHMPIDPid.cxx
HMPID/AliHMPIDPid.h
HMPID/AliHMPIDRecon.cxx

index 31f55fa..38ed89d 100644 (file)
@@ -367,3 +367,18 @@ Double_t AliHMPIDParam::SigGeom(Double_t trkTheta,Double_t trkPhi,Double_t theta
   return trErr*dtdT;
 }//SigGeom()
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+Double_t AliHMPIDParam::SigmaCorrFact  (Int_t iPart, Double_t occupancy) 
+{
+  Double_t corr = 1.0;
+                                                                                                             
+  switch(iPart) {
+    case 0: corr = 0.115*occupancy + 1.166; break; 
+    case 1: corr = 0.115*occupancy + 1.166; break;
+    case 2: corr = 0.115*occupancy + 1.166; break;
+    case 3: corr = 0.065*occupancy + 1.137; break;
+    case 4: corr = 0.048*occupancy + 1.202; break;
+  }
+                                                                                                                           
+ return corr; 
+}
+
index 6aa31a7..69756b6 100644 (file)
@@ -143,7 +143,8 @@ public:
   Double_t SigGeom        (Double_t trkTheta,Double_t trkPhi,Double_t ckovTh,Double_t ckovPh,Double_t beta);//error due to unknown photon origin
   Double_t SigCrom        (Double_t trkTheta,Double_t trkPhi,Double_t ckovTh,Double_t ckovPh,Double_t beta);//error due to unknonw photon energy
   Double_t Sigma2         (Double_t trkTheta,Double_t trkPhi,Double_t ckovTh,Double_t ckovPh              );//photon candidate sigma^2
-  static Double_t SigmaCorrFact  (Double_t occupancy                                                             ) {return 0.109*occupancy + 1.15;}//correction facotor for theoretical resolution
+  
+  static Double_t SigmaCorrFact(Int_t iPart, Double_t occupancy                                         );//correction factor for theoretical resolution
 
   //Mathieson Getters
   
index 4a66750..c653c7c 100644 (file)
@@ -72,7 +72,7 @@ void AliHMPIDPid::FindPid(AliESDtrack *pTrk,Double_t nmean,Int_t nsp,Double_t *p
     Double_t cosThetaTh = TMath::Sqrt(mass*mass+pmod*pmod)/(nmean*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);
+    Double_t sigmaRing = Resolution(iPart,thetaCerTh,pTrk);
     
     if(sigmaRing==0) {
       for(Int_t jPart=0;jPart<nsp;jPart++) prob[jPart]=1.0/(Float_t)nsp;
@@ -96,13 +96,16 @@ void AliHMPIDPid::FindPid(AliESDtrack *pTrk,Double_t nmean,Int_t nsp,Double_t *p
   delete pPid ; pPid=0x0;
 }
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-Double_t AliHMPIDPid::Resolution(Double_t thetaCerTh, AliESDtrack *pTrk)
+Double_t AliHMPIDPid::Resolution(Int_t iPart, Double_t thetaCerTh, AliESDtrack *pTrk)
 {
   AliHMPIDParam *pParam = AliHMPIDParam::Instance();
       
   AliHMPIDRecon rec;
   Float_t xPc,yPc,thRa,phRa;
+  Float_t x, y;
+  Int_t q, nph;
   pTrk->GetHMPIDtrk(xPc,yPc,thRa,phRa);
+  pTrk->GetHMPIDmip(x,y,q,nph);
   
   Double_t xRa = xPc - (pParam->RadThick()+pParam->WinThick()+pParam->GapThick())*TMath::Cos(phRa)*TMath::Tan(thRa); //just linear extrapolation back to RAD
   Double_t yRa = yPc - (pParam->RadThick()+pParam->WinThick()+pParam->GapThick())*TMath::Sin(phRa)*TMath::Tan(thRa); //just linear extrapolation back to RAD
@@ -110,18 +113,24 @@ Double_t AliHMPIDPid::Resolution(Double_t thetaCerTh, AliESDtrack *pTrk)
   
   Double_t occupancy = pTrk->GetHMPIDoccupancy();
   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);
+  Int_t    nPhotsTh = (Int_t)(12.*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;
+    
+    Int_t nPhots = 0; 
+    if(nph<nPhotsTh+TMath::Sqrt(nPhotsTh) && nph>nPhotsTh-TMath::Sqrt(nPhotsTh)) nPhots = nph;
+    else nPhots = gRandom->Poisson(nPhotsTh);
+  
     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;
+      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;
@@ -130,5 +139,6 @@ Double_t AliHMPIDPid::Resolution(Double_t thetaCerTh, AliESDtrack *pTrk)
     }      
     if(invSigma!=0) sigmatot += 1./TMath::Sqrt(invSigma);  
   }
-  return (sigmatot/nTrks)*pParam->SigmaCorrFact(occupancy);
+  
+  return (sigmatot/nTrks)*pParam->SigmaCorrFact(iPart,occupancy);
 }
index b8a4bb9..ef070bd 100644 (file)
@@ -23,7 +23,7 @@ public :
     virtual ~AliHMPIDPid() {;} //dtor
     
     void     FindPid(AliESDtrack *pESD,Double_t nmean,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
+    Double_t Resolution(Int_t iPart, Double_t thetaCerTh, AliESDtrack *pTrk);   //Find the sigma for a given ThetaCerTh
 
 //
 protected:
index 6f05d3f..792da47 100644 (file)
@@ -133,7 +133,7 @@ void AliHMPIDRecon::CkovAngle(AliESDtrack *pTrk,TClonesArray *pCluLst,Int_t inde
   pTrk->SetHMPIDmip(mipX,mipY,mipQ,fPhotCnt);                                                 //store mip info in any case 
   pTrk->SetHMPIDcluIdx(chId,index+1000*sizeClu);                                              //set index of cluster
   
-  if(fPhotCnt<=nMinPhotAcc) {                                                                 //no reconstruction with <=3 photon candidates
+  if(fPhotCnt<nMinPhotAcc) {                                                                 //no reconstruction with <=3 photon candidates
     pTrk->SetHMPIDsignal(kNoPhotAccept);                                                      //set the appropriate flag
     return;
   }
@@ -146,10 +146,10 @@ void AliHMPIDRecon::CkovAngle(AliESDtrack *pTrk,TClonesArray *pCluLst,Int_t inde
     
   Int_t iNrec=FlagPhot(HoughResponse(),pCluLst,pTrk);                                                      //flag photons according to individual theta ckov with respect to most probable
   
-  pTrk->SetHMPIDmip(mipX,mipY,mipQ,iNrec);                                                    //store mip info 
+  pTrk->SetHMPIDmip(mipX,mipY,mipQ,iNrec);                                                  //store mip info 
 
-  if(iNrec<1){
-    pTrk->SetHMPIDsignal(kNoPhotAccept);                                                      //no photon candidates are accepted
+  if(iNrec<nMinPhotAcc){
+    pTrk->SetHMPIDsignal(kNoPhotAccept);                                                    //no photon candidates are accepted
     return;
   }
   
@@ -157,8 +157,8 @@ void AliHMPIDRecon::CkovAngle(AliESDtrack *pTrk,TClonesArray *pCluLst,Int_t inde
   
   Double_t thetaC = FindRingCkov(pCluLst->GetEntries());                                    //find the best reconstructed theta Cherenkov
 //    FindRingGeom(thetaC,2);
-  pTrk->SetHMPIDsignal(thetaC+occupancy);                                                             //store theta Cherenkov
-  pTrk->SetHMPIDchi2(fCkovSigma2);                                                              //store errors squared
+  pTrk->SetHMPIDsignal(thetaC+occupancy);                                                   //store theta Cherenkov and chmaber occupancy
+  pTrk->SetHMPIDchi2(fCkovSigma2);                                                          //store experimental ring angular resolution squared
   
   DeleteVars();
 }//CkovAngle()