Adding more bins in QA (Alis)
[u/mrichter/AliRoot.git] / HMPID / AliHMPIDRecon.cxx
index d9b40bd..792da47 100644 (file)
 //////////////////////////////////////////////////////////////////////////
 
 #include "AliHMPIDRecon.h"   //class header
-#include "AliHMPIDParam.h"   //CkovAngle()
 #include "AliHMPIDCluster.h" //CkovAngle()
-#include <TMinuit.h>         //FitEllipse()
 #include <TRotation.h>       //TracePhot()
 #include <TH1D.h>            //HoughResponse()
 #include <TClonesArray.h>    //CkovAngle()
 #include <AliESDtrack.h>     //CkovAngle()
-
-const Double_t AliHMPIDRecon::fgkRadThick=1.5;
-const Double_t AliHMPIDRecon::fgkWinThick=0.5;
-const Double_t AliHMPIDRecon::fgkGapThick=8.0;
-const Double_t AliHMPIDRecon::fgkWinIdx  =1.5787;
-const Double_t AliHMPIDRecon::fgkGapIdx  =1.0005;
+#include <AliESDfriendTrack.h>     //CkovAngle()
 
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-AliHMPIDRecon::AliHMPIDRecon():TTask("RichRec","RichPat"),
-  fRadNmean(1.292),  
+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),
-  fTrkDir(TVector3(0,0,1)),fTrkPos(TVector2(30,40))  
+  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())
 {
-// main ctor
-  for (Int_t i=0; i<3000; i++) {
-    fPhotFlag[i] =  0;
-    fPhotCkov[i] = -1;
-    fPhotPhi [i] = -1;
-    fPhotWei [i] =  0;
-  }
-//hidden algorithm
-  fMipX=fMipY=fThTrkFit=fPhTrkFit=fCkovFit=fMipQ=fRadX=fRadY=-999;
-  fIdxMip=fNClu=0;
-  fCkovSig2=0;
-  for (Int_t i=0; i<100; i++) {
-    fXClu[i] = fYClu[i] = 0;
-    fClCk[i] = kTRUE;
-  }
+//..
+//init of data members
+//..
+  
+  fParam->SetRefIdx(fParam->MeanIdxRad()); // initialization of ref index to a default one
 }
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-void AliHMPIDRecon::CkovAngle(AliESDtrack *pTrk,TClonesArray *pCluLst,Double_t nmean,Double_t qthre)
+void AliHMPIDRecon::InitVars(Int_t n)
+{
+//..
+//Init some variables
+//..
+  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()const
+{
+//..
+//Delete variables
+//..
+  delete [] fPhotFlag;
+  delete [] fPhotClusIndex;
+  delete [] fPhotCkov;
+  delete [] fPhotPhi;
+  delete [] fPhotWei;
+}
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+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], 
-    
-  AliHMPIDParam *pParam=AliHMPIDParam::Instance();
+       
+  const Int_t nMinPhotAcc = 3;                      // Minimum number of photons required to perform the pattern recognition
   
-  if(pCluLst->GetEntries()>pParam->MultCut()) fIsWEIGHT = kTRUE; // offset to take into account bkg in reconstruction
-  else                                        fIsWEIGHT = kFALSE;
+  Int_t nClusTot = pCluLst->GetEntries();
 
-  Float_t xRa,yRa,th,ph;
-  pTrk->GetHMPIDtrk(xRa,yRa,th,ph);        //initialize this track: th and ph angles at middle of RAD 
+  InitVars(nClusTot);
+  
+  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);
 
-  fRadNmean=nmean;
+  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>pParam->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()
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 Bool_t AliHMPIDRecon::FindPhotCkov(Double_t cluX,Double_t cluY,Double_t &thetaCer,Double_t &phiCer)
@@ -125,9 +171,9 @@ Bool_t AliHMPIDRecon::FindPhotCkov(Double_t cluX,Double_t cluY,Double_t &thetaCe
 
   TVector3 dirCkov;
   
-  Double_t zRad= -0.5*fgkRadThick-0.5*fgkWinThick;                   //z position of middle of RAD
+  Double_t zRad= -0.5*fParam->RadThick()-0.5*fParam->WinThick();     //z position of middle of RAD
   TVector3 rad(fTrkPos.X(),fTrkPos.Y(),zRad);                        //impact point at middle of RAD
-  TVector3  pc(cluX,cluY,0.5*fgkWinThick+fgkGapIdx);                 //mip at PC
+  TVector3  pc(cluX,cluY,0.5*fParam->WinThick()+fParam->GapIdx());   //mip at PC
   Double_t cluR = TMath::Sqrt((cluX-fTrkPos.X())*(cluX-fTrkPos.X())+
                               (cluY-fTrkPos.Y())*(cluY-fTrkPos.Y()));//ref. distance impact RAD-CLUSTER   
   Double_t phi=(pc-rad).Phi();                                       //phi of photon
@@ -148,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;
     }
   }
@@ -159,30 +205,32 @@ TVector2 AliHMPIDRecon::TraceForward(TVector3 dirCkov)const
   //Trace forward a photon from (x,y) up to PC
   // Arguments: dirCkov photon vector in LORS
   //   Returns: pos of traced photon at PC
+  
   TVector2 pos(-999,-999);
   Double_t thetaCer = dirCkov.Theta();
-  if(thetaCer > TMath::ASin(1./fRadNmean))  return pos;         //total refraction on WIN-GAP boundary
-  Double_t zRad= -0.5*fgkRadThick-0.5*fgkWinThick;              //z position of middle of RAD
-  TVector3  posCkov(fTrkPos.X(),fTrkPos.Y(),zRad);              //RAD: photon position is track position @ middle of RAD 
-  Propagate(dirCkov,posCkov,           -0.5*fgkWinThick);       //go to RAD-WIN boundary  
-  Refract  (dirCkov,         fRadNmean,fgkWinIdx);              //RAD-WIN refraction
-  Propagate(dirCkov,posCkov,            0.5*fgkWinThick);       //go to WIN-GAP boundary
-  Refract  (dirCkov,         fgkWinIdx,fgkGapIdx);              //WIN-GAP refraction
-  Propagate(dirCkov,posCkov,0.5*fgkWinThick+fgkGapThick);       //go to PC
+  if(thetaCer > TMath::ASin(1./fParam->GetRefIdx())) return pos;          //total refraction on WIN-GAP boundary
+  Double_t zRad= -0.5*fParam->RadThick()-0.5*fParam->WinThick();          //z position of middle of RAD
+  TVector3  posCkov(fTrkPos.X(),fTrkPos.Y(),zRad);                        //RAD: photon position is track position @ middle of RAD 
+  Propagate(dirCkov,posCkov,           -0.5*fParam->WinThick());          //go to RAD-WIN boundary  
+  Refract  (dirCkov,         fParam->GetRefIdx(),fParam->WinIdx());       //RAD-WIN refraction
+  Propagate(dirCkov,posCkov,            0.5*fParam->WinThick());          //go to WIN-GAP boundary
+  Refract  (dirCkov,         fParam->WinIdx(),fParam->GapIdx());          //WIN-GAP refraction
+  Propagate(dirCkov,posCkov,0.5*fParam->WinThick()+fParam->GapThick());   //go to PC
   pos.Set(posCkov.X(),posCkov.Y());
   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;
@@ -190,36 +238,62 @@ 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=100;
-  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(!AliHMPIDParam::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;
    }
    TVector2 pos2=TracePhot(ckovAng,Double_t(TMath::TwoPi()*(i+1)/kN));                              //trace the next photon
    if(pos2.X()==-999) continue;                                                                     //no area: open ring            
-   if(!AliHMPIDParam::IsInside(pos2.X(),pos2.Y(),0)) {
+   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
 {
@@ -228,8 +302,6 @@ TVector2 AliHMPIDRecon::IntWithEdge(TVector2 p1,TVector2 p2)const
 // Arguments: 2 points obtained tracing the photons
 //   Returns: intersection point with detector (PC) edges
 
-  AliHMPIDParam *pParam = AliHMPIDParam::Instance();
-  
   Double_t xmin = (p1.X()<p2.X())? p1.X():p2.X(); 
   Double_t xmax = (p1.X()<p2.X())? p2.X():p1.X(); 
   Double_t ymin = (p1.Y()<p2.Y())? p1.Y():p2.Y(); 
@@ -239,26 +311,22 @@ TVector2 AliHMPIDRecon::IntWithEdge(TVector2 p1,TVector2 p2)const
   TVector2 pint;
   //intersection with low  X
   pint.Set((Double_t)(p1.X() + (0-p1.Y())/m),0.);
-  pint.Print();
-  if(pint.X()>=0 && pint.X()<=pParam->SizeAllX() &&
+  if(pint.X()>=0 && pint.X()<=fParam->SizeAllX() &&
      pint.X()>=xmin && pint.X()<=xmax            &&
      pint.Y()>=ymin && pint.Y()<=ymax) return pint;
   //intersection with high X  
-  pint.Set((Double_t)(p1.X() + (pParam->SizeAllY()-p1.Y())/m),(Double_t)(pParam->SizeAllY()));
-  pint.Print();
-  if(pint.X()>=0 && pint.X()<=pParam->SizeAllX() &&
+  pint.Set((Double_t)(p1.X() + (fParam->SizeAllY()-p1.Y())/m),(Double_t)(fParam->SizeAllY()));
+  if(pint.X()>=0 && pint.X()<=fParam->SizeAllX() &&
      pint.X()>=xmin && pint.X()<=xmax            &&
      pint.Y()>=ymin && pint.Y()<=ymax) return pint;
   //intersection with left Y  
   pint.Set(0.,(Double_t)(p1.Y() + m*(0-p1.X())));
-  pint.Print();
-  if(pint.Y()>=0 && pint.Y()<=pParam->SizeAllY() &&
+  if(pint.Y()>=0 && pint.Y()<=fParam->SizeAllY() &&
      pint.Y()>=ymin && pint.Y()<=ymax            &&
      pint.X()>=xmin && pint.X()<=xmax) return pint;
   //intersection with righ Y  
-  pint.Set((Double_t)(pParam->SizeAllX()),(Double_t)(p1.Y() + m*(pParam->SizeAllX()-p1.X())));
-  pint.Print();
-  if(pint.Y()>=0 && pint.Y()<=pParam->SizeAllY() &&
+  pint.Set((Double_t)(fParam->SizeAllX()),(Double_t)(p1.Y() + m*(fParam->SizeAllX()-p1.X())));
+  if(pint.Y()>=0 && pint.Y()<=fParam->SizeAllY() &&
      pint.Y()>=ymin && pint.Y()<=ymax            &&
      pint.X()>=xmin && pint.X()<=xmax) return pint;
   return p1;
@@ -284,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
   
@@ -295,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()
@@ -316,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 DRS, [rad]    
+// 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
@@ -388,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);
@@ -403,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
@@ -413,395 +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)*fRadNmean);
-  
-  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
-{
-// 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.-fRadNmean*fRadNmean+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)/(fgkGapThick*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*fRadNmean*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)
-  
-  Double_t k = 1.-fRadNmean*fRadNmean+alpha*alpha/(betaM*betaM);                             // formula (after 8 in the text)
-  if (k<0) return 1e10;
-
-  Double_t eTr = 0.5*fgkRadThick*betaM*TMath::Sqrt(k)/(fgkGapThick*alpha);                   // formula (14)
-  Double_t lambda = 1.-sint*sint*sinf*sinf;                                                  // formula (15)
-
-  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/(fgkGapThick*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)/(fgkGapThick*betaM);                      // formula (13.d)
-  Double_t dtdT = c1 * (c2+c3*c4);
-  Double_t trErr = fgkRadThick/(TMath::Sqrt(12.)*cost);
-
-  return trErr*dtdT;
-}//SigGeom()
-//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-//
-// From here HTA....
-//
-//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-Bool_t AliHMPIDRecon::CkovHiddenTrk(AliESDtrack *pTrk,TClonesArray *pCluLst,Double_t nmean, Double_t qthre)
-{
-// Pattern recognition method without any infos from tracking:HTA (Hidden Track Algorithm)...
-// The method finds in the chmber the cluster with the highest charge
-// compatibile with a MIP, then the strategy is applied
-// Arguments:  pTrk     - pointer to ESD track
-//             pCluLs   - list of clusters for a given chamber 
-//             nmean    - mean freon ref. index
-//   Returns:           - 0=ok,1=not fitted 
-  
-  fRadNmean=nmean;
-
-  if(pCluLst->GetEntriesFast()>100) return kFALSE;                                            //boundary check for CluX,CluY...
-  Float_t mipX=-1,mipY=-1;Int_t mipId=-1,mipQ=-1;                                                                           
-  Double_t qRef = 0;
-  Int_t nCh=0;
-  for (Int_t iClu=0;iClu<pCluLst->GetEntriesFast();iClu++){                                   //clusters loop
-    AliHMPIDCluster *pClu=(AliHMPIDCluster*)pCluLst->UncheckedAt(iClu);                       //get pointer to current cluster    
-    nCh = pClu->Ch();
-    fXClu[iClu] = pClu->X();fYClu[iClu] = pClu->Y();                                          //store x,y for fitting procedure
-    fClCk[iClu] = kTRUE;                                                                      //all cluster are accepted at this stage to be reconstructed
-    if(pClu->Q()>qRef){                                                                       //searching the highest charge to select a MIP      
-      qRef = pClu->Q();
-      mipId=iClu; mipX=pClu->X();mipY=pClu->Y();mipQ=(Int_t)pClu->Q();
-    }                                                                                    
-  }//clusters loop
-
-  fNClu = pCluLst->GetEntriesFast();
-  if(qRef>qthre){                                                                     //charge compartible with MIP clusters
-    fIdxMip = mipId;
-    fClCk[mipId] = kFALSE;
-    fMipX = mipX; fMipY=mipY; fMipQ = qRef;
-    if(!DoRecHiddenTrk(pCluLst)) {
-      pTrk->SetHMPIDsignal(kNoPhotAccept);
-      return kFALSE;
-    }                                                                           //Do track and ring reconstruction,if problems returns 1
-    pTrk->SetHMPIDtrk(fRadX,fRadY,fThTrkFit,fPhTrkFit);                                        //store track intersection info
-    pTrk->SetHMPIDmip(fMipX,fMipY,(Int_t)fMipQ,fNClu);                                         //store mip info 
-    pTrk->SetHMPIDcluIdx(nCh,fIdxMip);                                                         //set cham number and index of cluster
-    pTrk->SetHMPIDsignal(fCkovFit);                                                            //find best Theta ckov for ring i.e. track
-    pTrk->SetHMPIDchi2(fCkovSig2);                                                             //errors squared
-//    Printf(" n clusters tot %i accepted %i",pCluLst->GetEntriesFast(),fNClu);
-    return kTRUE;
-  }
-  
-  return kFALSE;
-}//CkovHiddenTrk()
-//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-Bool_t AliHMPIDRecon::DoRecHiddenTrk(TClonesArray *pCluLst)
-{
-// Pattern recognition method without any infos from tracking...
-// First a preclustering filter to avoid part of the noise
-// Then only ellipsed-rings are fitted (no possibility, 
-// for the moment, to reconstruct very inclined tracks)
-// Finally a fitting with (th,ph) free, starting by very close values
-// previously evaluated.
-// Arguments:   none
-//   Returns:   none
-  Double_t phiRec;
-  if(!CluPreFilter(pCluLst)) {return kFALSE;}
-  if(!FitEllipse(phiRec)) {return kFALSE;}
-  Int_t nClTmp1 = pCluLst->GetEntriesFast()-1;  //minus MIP...
-  Int_t nClTmp2 = 0;
-  while(nClTmp1 != nClTmp2){
-    SetNClu(pCluLst->GetEntriesFast());
-    if(!FitFree(phiRec)) {return kFALSE;}
-    nClTmp2 = NClu();
-    if(nClTmp2!=nClTmp1) {nClTmp1=nClTmp2;nClTmp2=0;}
-  }
-  fNClu = nClTmp2;
-  return kTRUE;
-}//DoRecHiddenTrk()
-//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-Bool_t AliHMPIDRecon::CluPreFilter(TClonesArray *pCluLst)
-{
-// Filter of bkg clusters
-// based on elliptical-shapes...
-//
-  if(pCluLst->GetEntriesFast()>50||pCluLst->GetEntriesFast()<4) return kFALSE; 
-  else return kTRUE;
-}
-//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-Bool_t AliHMPIDRecon::FitEllipse(Double_t &phiRec)
-{
-//Fit a set of clusters with an analitical conical section function:
-  //
-  // Ax^2 + B*y^2 + 2Hxy + 2Gx + 2Fy + 1 = 0   ---> conical section
-  //
-  //  H*H - A*B > 0 hyperbola
-  //            < 0 ellipse
-  //            = 0 parabola
-  //
-  // tan 2alfa = 2H/(A-B)  alfa=angle of rotation
-  //
-  // coordinate of the centre of the conical section:
-  //   x = x' + a
-  //   y = y' + b
-  //
-  //       HF - BG
-  //  a = ---------
-  //       AB - H^2
-  //
-  //       HG - AF
-  //  b = --------
-  //       AB - H^2
-  Double_t cA,cB,cF,cG,cH;
-  Double_t aArg=-1;      Int_t iErrFlg;                                                //tmp vars for TMinuit
-
-  if(!gMinuit) gMinuit = new TMinuit(5);                                               //init MINUIT with this number of parameters (5 params)
-  gMinuit->mncler();                                                                   // reset Minuit list of paramters
-  gMinuit->SetObjectFit((TObject*)this);  gMinuit->SetFCN(AliHMPIDRecon::FunMinEl);    //set fit function
-  gMinuit->mnexcm("SET PRI",&aArg,1,iErrFlg);                                          //suspend all printout from TMinuit 
-  gMinuit->mnexcm("SET NOW",&aArg,0,iErrFlg);                                          //suspend all warning printout from TMinuit
-  
-  Double_t d1,d2,d3;
-  TString sName;
-
-  gMinuit->mnparm(0," A ",1,0.01,0,0,iErrFlg);
-  gMinuit->mnparm(1," B ",1,0.01,0,0,iErrFlg);
-  gMinuit->mnparm(2," H ",1,0.01,0,0,iErrFlg);
-  gMinuit->mnparm(3," G ",1,0.01,0,0,iErrFlg);
-  gMinuit->mnparm(4," F ",1,0.01,0,0,iErrFlg);
-
-  gMinuit->mnexcm("SIMPLEX",&aArg,0,iErrFlg);
-  gMinuit->mnexcm("MIGRAD" ,&aArg,0,iErrFlg);   
-  gMinuit->mnpout(0,sName,cA,d1,d2,d3,iErrFlg);
-  gMinuit->mnpout(1,sName,cB,d1,d2,d3,iErrFlg);
-  gMinuit->mnpout(2,sName,cH,d1,d2,d3,iErrFlg);
-  gMinuit->mnpout(3,sName,cG,d1,d2,d3,iErrFlg);
-  gMinuit->mnpout(4,sName,cF,d1,d2,d3,iErrFlg);
-  delete gMinuit;
-
-  Double_t i2 = cA*cB-cH*cH;                                       //quartic invariant : i2 > 0  ellipse, i2 < 0 hyperbola
-  if(i2<=0) return kFALSE;
-  Double_t aX = (cH*cF-cB*cG)/i2;                                  //x centre of the canonical section 
-  Double_t bY = (cH*cG-cA*cF)/i2;                                  //y centre of the canonical section 
-  Double_t alfa1 = TMath::ATan(2*cH/(cA-cB));                      //alpha = angle of rotation of the conical section
-  if(alfa1<0) alfa1+=TMath::Pi(); 
-  alfa1*=0.5;
-//  Double_t alfa2 = alfa1+TMath::Pi();
-  Double_t phiref = TMath::ATan2(bY-fMipY,aX-fMipX);               //evaluate in a unique way the angle of rotation comparing it
-  if(phiref<0) phiref+=TMath::TwoPi();                             //with the vector that points to the centre from the mip 
-  if(i2<0) phiref+=TMath::Pi();
-  if(phiref>TMath::TwoPi()) phiref-=TMath::TwoPi();
-
-//  Printf(" alfa1 %f",alfa1*TMath::RadToDeg());
-//  Printf(" alfa2 %f",alfa2*TMath::RadToDeg());
-//  Printf(" firef %f",phiref*TMath::RadToDeg());
-//  if(TMath::Abs(alfa1-phiref)<TMath::Abs(alfa2-phiref)) phiRec = alfa1; else phiRec = alfa2;  
-  
-//  Printf("FitEllipse: phi reconstructed %f",phiRec*TMath::RadToDeg());
-  phiRec=phiref;
-  return kTRUE;
-//
-}
-//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-Bool_t AliHMPIDRecon::FitFree(Double_t phiRec)
+  Double_t AliHMPIDRecon::FindRingExt(Double_t ckov,Int_t ch,Double_t xPc,Double_t yPc,Double_t thRa,Double_t phRa)
 {
-// Fit performed by minimizing RMS/sqrt(n) of the
-// photons reconstructed. First phi is fixed and theta
-// is fouond, then (th,ph) of the track
-// as free parameters
-// Arguments:    PhiRec phi of the track
-//   Returns:    none
-  Double_t aArg=-1;  Int_t iErrFlg;                                                    //tmp vars for TMinuit
-  if(!gMinuit) gMinuit = new TMinuit(2);                                               //init MINUIT with this number of parameters (5 params)
-  gMinuit->mncler();                                                                   // reset Minuit list of paramters
-  gMinuit->SetObjectFit((TObject*)this);  gMinuit->SetFCN(AliHMPIDRecon::FunMinPhot);  //set fit function
-  gMinuit->mnexcm("SET PRI",&aArg,1,iErrFlg);                                          //suspend all printout from TMinuit 
-  gMinuit->mnexcm("SET NOW",&aArg,0,iErrFlg);                                          //suspend all warning printout from TMinuit
-  
-  Double_t d1,d2,d3;
-  TString sName;
-  Double_t th,ph;
-  
-  gMinuit->mnparm(0," theta ",  0.01,0.01,0,TMath::PiOver2(),iErrFlg);
-  gMinuit->mnparm(1," phi   ",phiRec,0.01,0,TMath::TwoPi()  ,iErrFlg);
-  
-  gMinuit->FixParameter(1);
-  gMinuit->mnexcm("SIMPLEX" ,&aArg,0,iErrFlg);   
-  gMinuit->mnexcm("MIGRAD"  ,&aArg,0,iErrFlg);
-  gMinuit->Release(1);  
-  gMinuit->mnexcm("MIGRAD"  ,&aArg,0,iErrFlg);
-  
-  gMinuit->mnpout(0,sName,th,d1,d2,d3,iErrFlg);
-  gMinuit->mnpout(1,sName,ph,d1,d2,d3,iErrFlg);   
-
-  Double_t outPar[2] = {th,ph}; Double_t g; Double_t f;Int_t flag = 3;
-  gMinuit->Eval(2, &g, f, outPar,flag);  
-
-  SetTrkFit(th,ph);
+// 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);
   
-  return kTRUE;
-}
-//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-Double_t AliHMPIDRecon::FunConSect(Double_t *c,Double_t x,Double_t y)
-{
-  return c[0]*x*x+c[1]*y*y+2*c[2]*x*y+2*c[3]*x+2*c[4]*y+1;
-}
-//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-void AliHMPIDRecon::FunMinEl(Int_t &/* */,Double_t* /* */,Double_t &f,Double_t *par,Int_t /* */)
-{
-  AliHMPIDRecon *pRec=(AliHMPIDRecon*)gMinuit->GetObjectFit();
-  Double_t minFun = 0;
-  Int_t np = pRec->NClu();
-  for(Int_t i=0;i<np;i++) {
-    if(i==pRec->IdxMip()) continue;
-    Double_t el = pRec->FunConSect(par,pRec->XClu(i),pRec->YClu(i));
-    minFun +=el*el;
-  }
-  f = minFun;
-}
-//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-void AliHMPIDRecon::FunMinPhot(Int_t &/* */,Double_t* /* */,Double_t &f,Double_t *par,Int_t iflag)
-{
-  AliHMPIDRecon *pRec=(AliHMPIDRecon*)gMinuit->GetObjectFit();
-  Double_t sizeCh = 0.5*fgkRadThick+fgkWinThick+fgkGapThick;
-  Double_t thTrk = par[0]; 
-  Double_t phTrk = par[1]; 
-  Double_t xrad = pRec->MipX() - sizeCh*TMath::Tan(thTrk)*TMath::Cos(phTrk);
-  Double_t yrad = pRec->MipY() - sizeCh*TMath::Tan(thTrk)*TMath::Sin(phTrk);
-  pRec->SetRadXY(xrad,yrad);
-  pRec->SetTrack(xrad,yrad,thTrk,phTrk);
+  Int_t nStep = 500;
+  Int_t nPhi = 0;  
 
-  Double_t meanCkov =0;
-  Double_t meanCkov2=0;
-  Double_t thetaCer,phiCer;
-  Int_t nClAcc = 0;
-  Int_t nClTot=pRec->NClu();
+  Int_t ipc,ipadx,ipady;
     
-  for(Int_t i=0;i<nClTot;i++) {
-    if(!(pRec->ClCk(i))) continue;
-    pRec->FindPhotCkov(pRec->XClu(i),pRec->YClu(i),thetaCer,phiCer);  
-    meanCkov  += thetaCer;
-    meanCkov2 += thetaCer*thetaCer;
-    nClAcc++;
-  }
-  if(nClAcc==0) {f=999;return;}
-  meanCkov/=nClAcc;
-  Double_t rms = (meanCkov2 - meanCkov*meanCkov*nClAcc)/nClAcc;
-  if(rms<0) Printf(" rms2 = %f, strange!!!",rms);
-  rms = TMath::Sqrt(rms);
-  f = rms/TMath::Sqrt(nClAcc);
-  
-  
-  if(iflag==3) {
-    Printf("FunMinPhot before: photons candidates %i used %i",nClTot,nClAcc);
-    nClAcc = 0;
-    Double_t meanCkov1=0;
-    Double_t meanCkov2=0;
-    for(Int_t i=0;i<nClTot;i++) {
-      if(!(pRec->ClCk(i))) continue;
-      pRec->FindPhotCkov(pRec->XClu(i),pRec->YClu(i),thetaCer,phiCer);  
-      if(TMath::Abs(thetaCer-meanCkov)<2*rms) {
-        meanCkov1 += thetaCer;
-        meanCkov2 += thetaCer*thetaCer;
-        nClAcc++;
-      } else pRec->SetClCk(i,kFALSE);
-    }
-    meanCkov1/=nClAcc;
-    Double_t rms2 = (meanCkov2 - meanCkov*meanCkov*nClAcc)/nClAcc;
-    Printf("FunMinPhot after: photons candidates %i used %i thetaCer %f",nClTot,nClAcc,meanCkov1);
-    pRec->SetCkovFit(meanCkov1);
-    pRec->SetCkovSig2(rms2);
-    pRec->SetNClu(nClAcc);
-  }
-}//FunMinPhot()
-//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-//
-// ended Hidden track algorithm....
-//
-//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+  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;
+}