Adding more bins in QA (Alis)
[u/mrichter/AliRoot.git] / HMPID / AliHMPIDRecon.cxx
index 44e125f..792da47 100644 (file)
 #include <TH1D.h>            //HoughResponse()
 #include <TClonesArray.h>    //CkovAngle()
 #include <AliESDtrack.h>     //CkovAngle()
+#include <AliESDfriendTrack.h>     //CkovAngle()
 
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-AliHMPIDRecon::AliHMPIDRecon():TTask("RichRec","RichPat")
+AliHMPIDRecon::AliHMPIDRecon():
+  TNamed("RichRec","RichPat"),
+  fPhotCnt(-1),
+  fPhotFlag(0x0),
+  fPhotClusIndex(0x0),    
+  fPhotCkov(0x0),
+  fPhotPhi(0x0),
+  fPhotWei(0x0),
+  fCkovSigma2(0),
+  fIsWEIGHT(kFALSE),
+  fDTheta(0.001),
+  fWindowWidth(0.045),
+  fRingArea(0),
+  fRingAcc(0),
+  fTrkDir(0,0,1),  // Just for test
+  fTrkPos(30,40),  // Just for test
+  fMipPos(0,0),
+  fPc(0,0),
+  fParam(AliHMPIDParam::Instance())
 {
 //..
 //init of data members
 //..
   
-  fPhotCnt  = -1;
-  fPhotFlag = 0x0;
-  fPhotCkov = 0x0;
-  fPhotPhi  = 0x0;
-  fPhotWei  = 0x0;
-  fCkovSigma2 = 0;
-  fIsWEIGHT = kFALSE;
-  fDTheta   = 0.001;
-  fWindowWidth = 0.045;
-  fTrkDir = TVector3(0,0,1); // init just for test
-  fTrkPos = TVector2(30,40); // init just for test
-  
-  AliHMPIDParam *pParam=AliHMPIDParam::Instance();
-  fParam = pParam;
-  
   fParam->SetRefIdx(fParam->MeanIdxRad()); // initialization of ref index to a default one
 }
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@@ -58,78 +62,104 @@ void AliHMPIDRecon::InitVars(Int_t n)
 //..
 //Init some variables
 //..
-  if(n<0) return;
+  if(n<=0) return;
   fPhotFlag = new Int_t[n];
+  fPhotClusIndex  = new Int_t[n];
   fPhotCkov = new Double_t[n];
   fPhotPhi  = new Double_t[n];
   fPhotWei  = new Double_t[n];
 //
 }
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-void AliHMPIDRecon::DeleteVars()
+void AliHMPIDRecon::DeleteVars()const
 {
 //..
 //Delete variables
 //..
-  delete fPhotFlag;
-  delete fPhotCkov;
-  delete fPhotPhi;
-  delete fPhotWei;
+  delete [] fPhotFlag;
+  delete [] fPhotClusIndex;
+  delete [] fPhotCkov;
+  delete [] fPhotPhi;
+  delete [] fPhotWei;
 }
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-void AliHMPIDRecon::CkovAngle(AliESDtrack *pTrk,TClonesArray *pCluLst,Double_t nmean,Double_t qthre)
+void AliHMPIDRecon::CkovAngle(AliESDtrack *pTrk,TClonesArray *pCluLst,Int_t index,Double_t nmean,Float_t xRa,Float_t yRa)
 {
 // Pattern recognition method based on Hough transform
 // Arguments:   pTrk     - track for which Ckov angle is to be found
 //              pCluLst  - list of clusters for this chamber   
 //   Returns:            - track ckov angle, [rad], 
-    
+       
+  const Int_t nMinPhotAcc = 3;                      // Minimum number of photons required to perform the pattern recognition
+  
   Int_t nClusTot = pCluLst->GetEntries();
-  if(nClusTot>fParam->MultCut()) fIsWEIGHT = kTRUE; // offset to take into account bkg in reconstruction
-  else                           fIsWEIGHT = kFALSE;
 
   InitVars(nClusTot);
   
-  Float_t xRa,yRa,th,ph;
-  pTrk->GetHMPIDtrk(xRa,yRa,th,ph);        //initialize this track: th and ph angles at middle of RAD 
+  Float_t xPc,yPc,th,ph;
+  pTrk->GetHMPIDtrk(xPc,yPc,th,ph);        //initialize this track: th and ph angles at middle of RAD 
   SetTrack(xRa,yRa,th,ph);
 
   fParam->SetRefIdx(nmean);
 
-  Float_t dMin=999,mipX=-1,mipY=-1;Int_t chId=-1,mipId=-1,mipQ=-1;                                                                           
-  fPhotCnt=0;                                                      
+  Float_t mipX=-1,mipY=-1;
+  Int_t chId=-1,mipQ=-1,sizeClu = -1;
+  
+  fPhotCnt=0;
+  
+  Int_t nPads = 0;
+  
   for (Int_t iClu=0; iClu<pCluLst->GetEntriesFast();iClu++){//clusters loop
-    AliHMPIDCluster *pClu=(AliHMPIDCluster*)pCluLst->UncheckedAt(iClu);                       //get pointer to current cluster    
+    AliHMPIDCluster *pClu=(AliHMPIDCluster*)pCluLst->UncheckedAt(iClu);                       //get pointer to current cluster
+    nPads+=pClu->Size();    
+    if(iClu == index) {                                                                       // this is the MIP! not a photon candidate: just store mip info
+      mipX = pClu->X();
+      mipY = pClu->Y();
+      mipQ=(Int_t)pClu->Q();
+      sizeClu=pClu->Size();
+      continue;                                                             
+    }
     chId=pClu->Ch();
-    if(pClu->Q()>qthre){                                                                      //charge compartible with MIP clusters      
-      Float_t dX=fPc.X()-pClu->X(),dY=fPc.Y()-pClu->Y(),d =TMath::Sqrt(dX*dX+dY*dY);          //distance between current cluster and intersection point
-      if( d < dMin) {mipId=iClu; dMin=d;mipX=pClu->X();mipY=pClu->Y();mipQ=(Int_t)pClu->Q();} //current cluster is closer, overwrite data for min cluster
-    }else{                                                                                    //charge compatible with photon cluster
-      Double_t thetaCer,phiCer;
-      if(FindPhotCkov(pClu->X(),pClu->Y(),thetaCer,phiCer)){                                  //find ckov angle for this  photon candidate
-        fPhotCkov[fPhotCnt]=thetaCer;                                                         //actual theta Cerenkov (in TRS)
-        fPhotPhi [fPhotCnt]=phiCer;                                                           //actual phi   Cerenkov (in TRS): -pi to come back to "unusual" ref system (X,Y,-Z)
-       //PH        Printf("photon n. %i reconstructed theta = %f",fPhotCnt,fPhotCkov[fPhotCnt]);
-        fPhotCnt++;                                                                           //increment counter of photon candidates
-      }
+    if(pClu->Q()>2*fParam->QCut()) continue;
+    Double_t thetaCer,phiCer;
+    if(FindPhotCkov(pClu->X(),pClu->Y(),thetaCer,phiCer)){                                    //find ckov angle for this  photon candidate
+      fPhotCkov[fPhotCnt]=thetaCer;                                                           //actual theta Cerenkov (in TRS)
+      fPhotPhi [fPhotCnt]=phiCer;
+      fPhotClusIndex[fPhotCnt]=iClu;                                                             //actual phi   Cerenkov (in TRS): -pi to come back to "unusual" ref system (X,Y,-Z)
+      fPhotCnt++;                                                                             //increment counter of photon candidates
     }
   }//clusters loop
-  fMipPos.Set(mipX,mipY);
-  if(fPhotCnt<=3) pTrk->SetHMPIDsignal(kNoPhotAccept);                                        //no reconstruction with <=3 photon candidates
-  Int_t iNacc=FlagPhot(HoughResponse());                                                      //flag photons according to individual theta ckov with respect to most probable
-  pTrk->SetHMPIDmip(mipX,mipY,mipQ,iNacc);                                                    //store mip info 
 
-  if(mipId==-1)              {pTrk->SetHMPIDsignal(kMipQdcCut);  return;}                     //no clusters with QDC more the threshold at all
-  if(dMin>fParam->DistCut()) {pTrk->SetHMPIDsignal(kMipDistCut); return;}                     //closest cluster with enough charge is still too far from intersection
-  pTrk->SetHMPIDcluIdx(chId,mipId);                                                           //set index of cluster
-  if(iNacc<1){
-    pTrk->SetHMPIDsignal(kNoPhotAccept);                                                      //no photon candidates is accepted
-  }
-  else {
-    pTrk->SetHMPIDsignal(FindRingCkov(pCluLst->GetEntries()));                                //find best Theta ckov for ring i.e. track
-    pTrk->SetHMPIDchi2(fCkovSigma2);                                                          //errors squared
+  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
+    pTrk->SetHMPIDsignal(kNoPhotAccept);                                                      //set the appropriate flag
+    return;
   }
+    
+  fMipPos.Set(mipX,mipY);
+    
+//PATTERN RECOGNITION STARTED: 
+  if(fPhotCnt>fParam->MultCut()) fIsWEIGHT = kTRUE; // offset to take into account bkg in reconstruction
+  else                           fIsWEIGHT = kFALSE;
+    
+  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 
 
+  if(iNrec<nMinPhotAcc){
+    pTrk->SetHMPIDsignal(kNoPhotAccept);                                                    //no photon candidates are accepted
+    return;
+  }
+  
+  Int_t occupancy = (Int_t)(1000*(nPads/(6.*80.*48.))); 
+  
+  Double_t thetaC = FindRingCkov(pCluLst->GetEntries());                                    //find the best reconstructed theta Cherenkov
+//    FindRingGeom(thetaC,2);
+  pTrk->SetHMPIDsignal(thetaC+occupancy);                                                   //store theta Cherenkov and chmaber occupancy
+  pTrk->SetHMPIDchi2(fCkovSigma2);                                                          //store experimental ring angular resolution squared
+  
   DeleteVars();
 }//CkovAngle()
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@@ -164,7 +194,7 @@ Bool_t AliHMPIDRecon::FindPhotCkov(Double_t cluX,Double_t cluY,Double_t &thetaCe
     else if(dist<-kTol) ckov2=ckov;                           //cluster @ smaller ckov
     else{                                                     //precision achived: ckov in DRS found
       dirCkov.SetMagThetaPhi(1,ckov,phi);                     //
-      RecPhot(dirCkov,thetaCer,phiCer);                       //find ckov (in TRS:the effective Cherenkov angle!)
+      Lors2Trs(dirCkov,thetaCer,phiCer);                       //find ckov (in TRS:the effective Cherenkov angle!)
       return kTRUE;
     }
   }
@@ -190,16 +220,17 @@ TVector2 AliHMPIDRecon::TraceForward(TVector3 dirCkov)const
   return pos;
 }//TraceForward()
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-void AliHMPIDRecon::RecPhot(TVector3 dirCkov,Double_t &thetaCer,Double_t &phiCer)
+void AliHMPIDRecon::Lors2Trs(TVector3 dirCkov,Double_t &thetaCer,Double_t &phiCer)const
 {
   //Theta Cerenkov reconstruction 
-  // Arguments: (x,y) of initial point in LORS, dirCkov photon vector in LORS
-  //   Returns: thetaCer theta cerenkov reconstructed
+  // Arguments: dirCkov photon vector in LORS
+  //   Returns: thetaCer of photon in TRS
+  //              phiCer of photon in TRS
 //  TVector3 dirTrk;
 //  dirTrk.SetMagThetaPhi(1,fTrkDir.Theta(),fTrkDir.Phi());
 //  Double_t thetaCer = TMath::ACos(dirCkov*dirTrk);
-  TRotation mtheta;   mtheta.RotateY(- fTrkDir.Theta());
-  TRotation mphi;       mphi.RotateZ(- fTrkDir.Phi());
+  TRotation mtheta;   mtheta.RotateY(-fTrkDir.Theta());
+  TRotation mphi;       mphi.RotateZ(-fTrkDir.Phi());
   TRotation mrot=mtheta*mphi;
   TVector3 dirCkovTRS;
   dirCkovTRS=mrot*dirCkov;
@@ -207,21 +238,44 @@ void AliHMPIDRecon::RecPhot(TVector3 dirCkov,Double_t &thetaCer,Double_t &phiCer
   thetaCer= dirCkovTRS.Theta();                                        //actual value of thetaCerenkov of the photon
 }
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-Double_t AliHMPIDRecon::FindRingArea(Double_t ckovAng)const
+void AliHMPIDRecon::Trs2Lors(TVector3 dirCkov,Double_t &thetaCer,Double_t &phiCer)const
+{
+  //Theta Cerenkov reconstruction 
+  // Arguments: dirCkov photon vector in TRS
+  //   Returns: thetaCer of photon in LORS
+  //              phiCer of photon in LORS
+  TRotation mtheta;   mtheta.RotateY(fTrkDir.Theta());
+  TRotation mphi;       mphi.RotateZ(fTrkDir.Phi());
+  TRotation mrot=mphi*mtheta;
+  TVector3 dirCkovLORS;
+  dirCkovLORS=mrot*dirCkov;
+  phiCer  = dirCkovLORS.Phi();                                          //actual value of the phi of the photon
+  thetaCer= dirCkovLORS.Theta();                                        //actual value of thetaCerenkov of the photon
+}
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+void AliHMPIDRecon::FindRingGeom(Double_t ckovAng,Int_t level)
 {
 // Find area covered in the PC acceptance
-// Arguments: ckovAng - cerenkov angle     
+// Arguments: ckovAng - cerenkov angle
+//            level   - precision in finding area and portion of ring accepted (multiple of 50)
 //   Returns: area of the ring in cm^2 for given theta ckov
    
-  const Int_t kN=50;
-  TVector2 pos1;
+  Int_t kN=50*level;
+  Int_t nPoints = 0;
   Double_t area=0;
+  
   Bool_t first=kFALSE;
+  TVector2 pos1;
+  
   for(Int_t i=0;i<kN;i++){
    if(!first) {
-     pos1=TracePhot(ckovAng,Double_t(TMath::TwoPi()*(i+1)/kN));                                     //find a good trace for the first photon
+      pos1=TracePhot(ckovAng,Double_t(TMath::TwoPi()*(i+1)/kN));                                    //find a good trace for the first photon
      if(pos1.X()==-999) continue;                                                                   //no area: open ring                 
-     if(!fParam->IsInside(pos1.X(),pos1.Y(),0)) pos1 = IntWithEdge(fMipPos,pos1);                   // ffind the very first intersection...
+     if(!fParam->IsInside(pos1.X(),pos1.Y(),0)) {
+       pos1 = IntWithEdge(fMipPos,pos1);                                                            // find the very first intersection...
+     } else {
+       if(!AliHMPIDParam::IsInDead(pos1.X(),pos1.Y())) nPoints++;                                   //photon is accepted if not in dead zone
+     }
      first=kTRUE;
      continue;
    }
@@ -229,14 +283,17 @@ Double_t AliHMPIDRecon::FindRingArea(Double_t ckovAng)const
    if(pos2.X()==-999) continue;                                                                     //no area: open ring            
    if(!fParam->IsInside(pos2.X(),pos2.Y(),0)) {
      pos2 = IntWithEdge(fMipPos,pos2);
+   } else {
+     if(!AliHMPIDParam::IsInDead(pos2.X(),pos2.Y())) nPoints++;                                     //photon is accepted if not in dead zone
    }
    area+=TMath::Abs((pos1-fMipPos).X()*(pos2-fMipPos).Y()-(pos1-fMipPos).Y()*(pos2-fMipPos).X());   //add area of the triangle...           
    pos1 = pos2;
   }
-//---  find points from ring
+//---  find area and length of the ring;
+  fRingAcc = (Double_t)nPoints/(Double_t)kN;
   area*=0.5;
-  return area;
-}//FindRingArea()
+  fRingArea = area;
+}//FindRingGeom()
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 TVector2 AliHMPIDRecon::IntWithEdge(TVector2 p1,TVector2 p2)const
 {
@@ -295,7 +352,7 @@ Double_t AliHMPIDRecon::FindRingCkov(Int_t)
       weightThetaCerenkov += fPhotCkov[i]*fPhotWei[i];
       wei += fPhotWei[i];                                                    //collect weight as sum of all candidate weghts   
       
-      sigma2 += 1./Sigma2(fPhotCkov[i],fPhotPhi[i]);
+      sigma2 += 1./fParam->Sigma2(fTrkDir.Theta(),fTrkDir.Phi(),fPhotCkov[i],fPhotPhi[i]);
     }
   }//candidates loop
   
@@ -306,7 +363,7 @@ Double_t AliHMPIDRecon::FindRingCkov(Int_t)
   return weightThetaCerenkov;
 }//FindCkovRing()
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-Int_t AliHMPIDRecon::FlagPhot(Double_t ckov)
+Int_t AliHMPIDRecon::FlagPhot(Double_t ckov,TClonesArray *pCluLst, AliESDtrack *pTrk)
 {
 // Flag photon candidates if their individual ckov angle is inside the window around ckov angle returned by  HoughResponse()
 // Arguments: ckov- value of most probable ckov angle for track as returned by HoughResponse()
@@ -327,27 +384,39 @@ Int_t AliHMPIDRecon::FlagPhot(Double_t ckov)
   Int_t iInsideCnt = 0; //count photons which Theta ckov inside the window
   for(Int_t i=0;i<fPhotCnt;i++){//photon candidates loop
     fPhotFlag[i] = 0;
-    if(fPhotCkov[i] >= tmin && fPhotCkov[i] <= tmax)   { 
-      fPhotFlag[i]=2;    
+    if(fPhotCkov[i] >= tmin && fPhotCkov[i] <= tmax) { 
+      fPhotFlag[i]=2;
+      AddObjectToFriends(pCluLst,i,pTrk);
       iInsideCnt++;
     }
   }
+      
   return iInsideCnt;
+  
 }//FlagPhot()
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+void  AliHMPIDRecon::AddObjectToFriends(TClonesArray *pCluLst, Int_t photonIndex, AliESDtrack *pTrk)
+{
+// Add AliHMPIDcluster object to ESD friends
+    
+  AliHMPIDCluster *pClu=(AliHMPIDCluster*)pCluLst->UncheckedAt(fPhotClusIndex[photonIndex]);     
+  AliHMPIDCluster *pClus = new AliHMPIDCluster(*pClu);
+  pClus->SetChi2(fPhotCkov[photonIndex]);  
+  pTrk->AddCalibObject(pClus);   
+}    
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 TVector2 AliHMPIDRecon::TracePhot(Double_t ckovThe,Double_t ckovPhi)const
 {
 // Trace a single Ckov photon from emission point somewhere in radiator up to photocathode taking into account ref indexes of materials it travereses
 // Arguments: ckovThe,ckovPhi- photon ckov angles in TRS, [rad]
 //   Returns: distance between photon point on PC and track projection  
-  TRotation mtheta;   mtheta.RotateY(fTrkDir.Theta());
-  TRotation mphi;       mphi.RotateZ(fTrkDir.Phi());  
-  TRotation mrot=mphi*mtheta;
-  TVector3  dirCkov,dirCkovTors;   
-
-  dirCkovTors.SetMagThetaPhi(1,ckovThe,ckovPhi);                    //initially photon is directed according to requested ckov angle
-  dirCkov=mrot*dirCkovTors;                                         //now we know photon direction in LORS
-  return TraceForward(dirCkov);
+  
+  Double_t theta,phi;
+  TVector3  dirTRS,dirLORS;   
+  dirTRS.SetMagThetaPhi(1,ckovThe,ckovPhi);                     //photon in TRS
+  Trs2Lors(dirTRS,theta,phi);
+  dirLORS.SetMagThetaPhi(1,theta,phi);                          //photon in LORS
+  return TraceForward(dirLORS);                                 //now foward tracing
 }//TracePhot()
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 void AliHMPIDRecon::Propagate(const TVector3 dir,TVector3 &pos,Double_t z)const
@@ -399,7 +468,11 @@ Double_t AliHMPIDRecon::HoughResponse()
     Double_t weight=1.;
     if(fIsWEIGHT){
       Double_t lowerlimit = ((Double_t)bin)*fDTheta - 0.5*fDTheta;  Double_t upperlimit = ((Double_t)bin)*fDTheta + 0.5*fDTheta;
-      Double_t diffArea = FindRingArea(upperlimit)-FindRingArea(lowerlimit);
+      FindRingGeom(lowerlimit);
+      Double_t areaLow  = GetRingArea();
+      FindRingGeom(upperlimit);
+      Double_t areaHigh = GetRingArea();
+      Double_t diffArea = areaHigh - areaLow;
       if(diffArea>0) weight = 1./diffArea;
     }
     photsw->Fill(angle,weight);
@@ -414,6 +487,7 @@ Double_t AliHMPIDRecon::HoughResponse()
     Double_t sumPhots=phots->Integral(bin1,bin2);
     if(sumPhots<3) continue;                            // if less then 3 photons don't trust to this ring
     Double_t sumPhotsw=photsw->Integral(bin1,bin2);
+    if((Double_t)((i+0.5)*fDTheta)>0.7) continue;
     resultw->Fill((Double_t)((i+0.5)*fDTheta),sumPhotsw);
   } 
 // evaluate the "BEST" theta ckov as the maximum value of histogramm
@@ -424,118 +498,32 @@ Double_t AliHMPIDRecon::HoughResponse()
   return (Double_t)(locMax*fDTheta+0.5*fDTheta); //final most probable track theta ckov   
 }//HoughResponse()
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-Double_t AliHMPIDRecon::Sigma2(Double_t ckovTh, Double_t ckovPh)const
-{
-// Analithical calculation of total error (as a sum of localization, geometrical and chromatic errors) on Cerenkov angle for a given Cerenkov photon 
-// created by a given MIP. Fromulae according to CERN-EP-2000-058 
-// Arguments: Cerenkov and azimuthal angles for Cerenkov photon, [radians]
-//            dip and azimuthal angles for MIP taken at the entrance to radiator, [radians]        
-//            MIP beta
-//   Returns: absolute error on Cerenkov angle, [radians]    
-  
-  TVector3 v(-999,-999,-999);
-  Double_t trkBeta = 1./(TMath::Cos(ckovTh)*fParam->GetRefIdx());
-  
-  if(trkBeta > 1) trkBeta = 1;                 //protection against bad measured thetaCer  
-  if(trkBeta < 0) trkBeta = 0.0001;            //
-
-  v.SetX(SigLoc (ckovTh,ckovPh,trkBeta));
-  v.SetY(SigGeom(ckovTh,ckovPh,trkBeta));
-  v.SetZ(SigCrom(ckovTh,ckovPh,trkBeta));
-
-  return v.Mag2();
-}
-//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-Double_t AliHMPIDRecon::SigLoc(Double_t thetaC, Double_t phiC,Double_t betaM)const
+  Double_t AliHMPIDRecon::FindRingExt(Double_t ckov,Int_t ch,Double_t xPc,Double_t yPc,Double_t thRa,Double_t phRa)
 {
-// Analithical calculation of localization error (due to finite segmentation of PC) on Cerenkov angle for a given Cerenkov photon 
-// created by a given MIP. Fromulae according to CERN-EP-2000-058 
-// Arguments: Cerenkov and azimuthal angles for Cerenkov photon, [radians]
-//            dip and azimuthal angles for MIP taken at the entrance to radiator, [radians]        
-//            MIP beta
-//   Returns: absolute error on Cerenkov angle, [radians]    
-  
-  Double_t phiDelta = phiC - fTrkDir.Phi();
-
-  Double_t sint     = TMath::Sin(fTrkDir.Theta());
-  Double_t cost     = TMath::Cos(fTrkDir.Theta());
-  Double_t sinf     = TMath::Sin(fTrkDir.Phi());
-  Double_t cosf     = TMath::Cos(fTrkDir.Phi());
-  Double_t sinfd    = TMath::Sin(phiDelta);
-  Double_t cosfd    = TMath::Cos(phiDelta);
-  Double_t tantheta = TMath::Tan(thetaC);
-  
-  Double_t alpha =cost-tantheta*cosfd*sint;                                                 // formula (11)
-  Double_t k = 1.-fParam->GetRefIdx()*fParam->GetRefIdx()+alpha*alpha/(betaM*betaM);        // formula (after 8 in the text)
-  if (k<0) return 1e10;
-  Double_t mu =sint*sinf+tantheta*(cost*cosfd*sinf+sinfd*cosf);                             // formula (10)
-  Double_t e  =sint*cosf+tantheta*(cost*cosfd*cosf-sinfd*sinf);                             // formula (9)
-
-  Double_t kk = betaM*TMath::Sqrt(k)/(fParam->GapThick()*alpha);                            // formula (6) and (7)
-  Double_t dtdxc = kk*(k*(cosfd*cosf-cost*sinfd*sinf)-(alpha*mu/(betaM*betaM))*sint*sinfd); // formula (6)           
-  Double_t dtdyc = kk*(k*(cosfd*sinf+cost*sinfd*cosf)+(alpha* e/(betaM*betaM))*sint*sinfd); // formula (7)            pag.4
-
-  Double_t errX = 0.2,errY=0.25;                                                            //end of page 7
-  return  TMath::Sqrt(errX*errX*dtdxc*dtdxc + errY*errY*dtdyc*dtdyc);
-}
-//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-Double_t AliHMPIDRecon::SigCrom(Double_t thetaC, Double_t phiC,Double_t betaM)const
-{
-// Analithical calculation of chromatic error (due to lack of knowledge of Cerenkov photon energy) on Cerenkov angle for a given Cerenkov photon 
-// created by a given MIP. Fromulae according to CERN-EP-2000-058 
-// Arguments: Cerenkov and azimuthal angles for Cerenkov photon, [radians]
-//            dip and azimuthal angles for MIP taken at the entrance to radiator, [radians]        
-//            MIP beta
-//   Returns: absolute error on Cerenkov angle, [radians]    
-  
-  Double_t phiDelta = phiC - fTrkDir.Phi();
-
-  Double_t sint     = TMath::Sin(fTrkDir.Theta());
-  Double_t cost     = TMath::Cos(fTrkDir.Theta());
-  Double_t cosfd    = TMath::Cos(phiDelta);
-  Double_t tantheta = TMath::Tan(thetaC);
-  
-  Double_t alpha =cost-tantheta*cosfd*sint;                                                 // formula (11)
-  Double_t dtdn = cost*fParam->GetRefIdx()*betaM*betaM/(alpha*tantheta);                    // formula (12)
-            
-//  Double_t f = 0.00928*(7.75-5.635)/TMath::Sqrt(12.);
-  Double_t f = 0.0172*(7.75-5.635)/TMath::Sqrt(24.);
-
-  return f*dtdn;
-}//SigCrom()
-//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-Double_t AliHMPIDRecon::SigGeom(Double_t thetaC, Double_t phiC,Double_t betaM)const
-{
-// Analithical calculation of geometric error (due to lack of knowledge of creation point in radiator) on Cerenkov angle for a given Cerenkov photon 
-// created by a given MIP. Formulae according to CERN-EP-2000-058 
-// Arguments: Cerenkov and azimuthal angles for Cerenkov photon, [radians]
-//            dip and azimuthal angles for MIP taken at the entrance to radiator, [radians]        
-//            MIP beta
-//   Returns: absolute error on Cerenkov angle, [radians]    
-
-  Double_t phiDelta = phiC - fTrkDir.Phi();
-
-  Double_t sint     = TMath::Sin(fTrkDir.Theta());
-  Double_t cost     = TMath::Cos(fTrkDir.Theta());
-  Double_t sinf     = TMath::Sin(fTrkDir.Phi());
-  Double_t cosfd    = TMath::Cos(phiDelta);
-  Double_t costheta = TMath::Cos(thetaC);
-  Double_t tantheta = TMath::Tan(thetaC);
-  
-  Double_t alpha =cost-tantheta*cosfd*sint;                                                  // formula (11)
+// To find the acceptance of the ring even from external inputs. 
+//    
+//       
+  Double_t xRa = xPc - (fParam->RadThick()+fParam->WinThick()+fParam->GapThick())*TMath::Cos(phRa)*TMath::Tan(thRa); //just linear extrapolation back to RAD
+  Double_t yRa = yPc - (fParam->RadThick()+fParam->WinThick()+fParam->GapThick())*TMath::Sin(phRa)*TMath::Tan(thRa);
   
-  Double_t k = 1.-fParam->GetRefIdx()*fParam->GetRefIdx()+alpha*alpha/(betaM*betaM);         // formula (after 8 in the text)
-  if (k<0) return 1e10;
-
-  Double_t eTr = 0.5*fParam->RadThick()*betaM*TMath::Sqrt(k)/(fParam->GapThick()*alpha);     // formula (14)
-  Double_t lambda = 1.-sint*sint*sinf*sinf;                                                  // formula (15)
+  Int_t nStep = 500;
+  Int_t nPhi = 0;  
 
-  Double_t c1 = 1./(1.+ eTr*k/(alpha*alpha*costheta*costheta));                              // formula (13.a)
-  Double_t c2 = betaM*TMath::Power(k,1.5)*tantheta*lambda/(fParam->GapThick()*alpha*alpha);  // formula (13.b)
-  Double_t c3 = (1.+eTr*k*betaM*betaM)/((1+eTr)*alpha*alpha);                                // formula (13.c)
-  Double_t c4 = TMath::Sqrt(k)*tantheta*(1-lambda)/(fParam->GapThick()*betaM);               // formula (13.d)
-  Double_t dtdT = c1 * (c2+c3*c4);
-  Double_t trErr = fParam->RadThick()/(TMath::Sqrt(12.)*cost);
-
-  return trErr*dtdT;
-}//SigGeom()
+  Int_t ipc,ipadx,ipady;
+    
+  if(ckov>0){
+    SetTrack(xRa,yRa,thRa,phRa);
+    for(Int_t j=0;j<nStep;j++){
+      TVector2 pos; pos=TracePhot(ckov,j*TMath::TwoPi()/(Double_t)(nStep-1));
+      if(fParam->IsInDead(pos.X(),pos.Y())) continue;
+      fParam->Lors2Pad(pos.X(),pos.Y(),ipc,ipadx,ipady);
+      ipadx+=(ipc%2)*fParam->kPadPcX;
+      ipady+=(ipc/2)*fParam->kPadPcY;
+      if(ipadx<0 || ipady>160 || ipady<0 || ipady>144 || ch<0 || ch>6) continue;
+      if(fParam->IsDeadPad(ipadx,ipady,ch)) continue;
+      nPhi++;
+    }//point loop
+  return ((Double_t)nPhi/(Double_t)nStep); 
+  }//if
+  return -1;
+}