Pattern recognition reorganized (for display,monitoring...)
authordibari <dibari@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 27 Mar 2007 12:32:33 +0000 (12:32 +0000)
committerdibari <dibari@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 27 Mar 2007 12:32:33 +0000 (12:32 +0000)
HMPID/AliHMPIDDigit.h
HMPID/AliHMPIDParam.h
HMPID/AliHMPIDRecon.cxx
HMPID/AliHMPIDRecon.h
HMPID/AliHMPIDTracker.cxx
HMPID/AliHMPIDTracker.h
HMPID/Hdisp.C

index b1fadfa..1466025 100644 (file)
@@ -36,7 +36,7 @@ public:
          void    AddTidOffset(Int_t offset                   )     {for (Int_t i=0; i<3; i++) if (fTracks[i]>0) fTracks[i]+=offset;  } //needed for merging
          Int_t   Ch          (                               )const{return A2C(fPad);                                                } //chamber number
   static Bool_t  IsOverTh    (Float_t q                      )     {return q >= fSigmas;                                             } //is digit over threshold????
-  static Bool_t  IsInside    (Float_t x,Float_t y            )     {return x>0&&y>0&&x<SizeAllX()&&y<SizeAllY();                     } //is point inside chamber boundary?
+  static Bool_t  IsInside    (Float_t x,Float_t y,Float_t margin=0){return x>-margin&&y>-margin&&x<SizeAllX()+margin&&y<SizeAllY()+margin;} //is point inside chamber boundary?
          Float_t LorsX       (                               )const{return LorsX(A2P(fPad),A2X(fPad));                               } //center of the pad x, [cm]
   static Float_t LorsX       (Int_t pc,Int_t padx            )     {return (padx    +0.5)*SizePadX()+(pc  %2)*(SizePcX()+SizeDead());} //center of the pad x, [cm]
          Float_t LorsY       (                               )const{return LorsY(A2P(fPad),A2Y(fPad));                               } //center of the pad y, [cm]
index e3c9ba5..3f51cc4 100644 (file)
@@ -22,8 +22,11 @@ public:
          void     Print(Option_t *opt="") const;                                         //print current parametrization
   static inline AliHMPIDParam* Instance();                                //pointer to AliHMPIDParam singleton
   
-                Double_t   MeanIdxRad              () {return 1.29204;}//???????????
-                Double_t   MeanIdxWin              () {return 1.57819;}//???????????
+                Double_t   MeanIdxRad              () {return 1.29204;}   //<--TEMPORAR--> to be removed in future  Mean ref index C6F14
+                Double_t   MeanIdxWin              () {return 1.57819;}   //<--TEMPORAR--> to be removed in future. Mean ref index quartz
+                Float_t    DistCut                 () {return 1.0;}       //<--TEMPORAR--> to be removed in future. Cut for MIP-TRACK residual 
+                Float_t    QCut                    () {return 100;}       //<--TEMPORAR--> to be removed in future. Separation PHOTON-MIP charge 
+                Float_t    MultCut                 () {return 200;}       //<--TEMPORAR--> to be removed in future. Multiplicity cut to activate WEIGHT procedure 
   static        Int_t      Stack(Int_t evt=-1,Int_t tid=-1);              //Print stack info for event and tid
   static        Int_t      StackCount(Int_t pid,Int_t evt);               //Counts stack particles of given sort in given event  
   static        void       IdealPosition(Int_t iCh,TGeoHMatrix *m);       //ideal position of given chamber 
index 7ec9a93..062ade9 100644 (file)
 // for single chamber                                                   //
 //////////////////////////////////////////////////////////////////////////
 
-#include "AliHMPIDRecon.h"  //class header
+#include "AliHMPIDRecon.h"   //class header
+#include "AliHMPIDParam.h"   //CkovAngle()
 #include "AliHMPIDCluster.h" //CkovAngle()
-#include <TRotation.h>     //TracePhoton()
-#include <TH1D.h>          //HoughResponse()
-#include <TClonesArray.h>  //CkovAngle()
-#include <AliESDtrack.h>   //CkovAngle()
+#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;
@@ -34,9 +35,6 @@ const Double_t AliHMPIDRecon::fgkGapThick=8.0;
 const Double_t AliHMPIDRecon::fgkWinIdx  =1.5787;
 const Double_t AliHMPIDRecon::fgkGapIdx  =1.0005;
 
-Double_t xRad;
-Double_t yRad;
-
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 AliHMPIDRecon::AliHMPIDRecon():TTask("RichRec","RichPat"),
   fRadNmean(1.292),  
@@ -56,94 +54,144 @@ AliHMPIDRecon::AliHMPIDRecon():TTask("RichRec","RichPat"),
   }
 }
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-void AliHMPIDRecon::CkovAngle(Double_t xRa,Double_t yRa,AliESDtrack *pTrk,TClonesArray *pCluLst,Double_t nmean)
+void AliHMPIDRecon::CkovAngle(AliESDtrack *pTrk,TClonesArray *pCluLst,Double_t nmean)
 {
 // 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();
   
-  if(pCluLst->GetEntries()>200) fIsWEIGHT = kTRUE; // offset to take into account bkg in reconstruction
-  else                          fIsWEIGHT = kFALSE;
+  if(pCluLst->GetEntries()>pParam->MultCut()) fIsWEIGHT = kTRUE; // offset to take into account bkg in reconstruction
+  else                                        fIsWEIGHT = kFALSE;
 
-  // Photon Flag:  Flag = 0 initial set; Flag = 1 good candidate (charge compatible with photon); Flag = 2 photon used for the ring;
-  Float_t xPc,yPc,th,ph;      pTrk->GetHMPIDtrk(xPc,yPc,th,ph);  SetTrack(xPc,yPc,th,ph); //initialize this track            
+  Float_t xRa,yRa,th,ph;       
+  pTrk->GetHMPIDtrk(xRa,yRa,th,ph);        //initialize this track: th and ph angles at middle of RAD 
+  
+  ph+=TMath::Pi();                         // right XYZ local orientation
+  SetTrack(xRa,yRa,th,ph);
+  
   fRadNmean=nmean;
 
-  xRad=xRa;yRad=yRa;  
-
   Float_t dMin=999,mipX=-1,mipY=-1;Int_t chId=-1,mipId=-1,mipQ=-1;                                                                           
   fPhotCnt=0;                                                      
   for (Int_t iClu=0; iClu<pCluLst->GetEntriesFast();iClu++){//clusters loop
     AliHMPIDCluster *pClu=(AliHMPIDCluster*)pCluLst->UncheckedAt(iClu);                       //get pointer to current cluster    
     chId=pClu->Ch();
-    if(pClu->Q()>100){                                                                        //charge compartible with MIP clusters      
-      Float_t dX=xPc-pClu->X(),dY=yPc-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 compartible with photon cluster
-      fPhotCkov[fPhotCnt]=FindPhotCkov(pClu->X(),pClu->Y());                                  //find ckov angle for this  photon candidate
-      fPhotCnt++;         //increment counter of photon candidates
+    if(pClu->Q()>pParam->QCut()){                                                             //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-TMath::Pi();                                               //actual phi   Cerenkov (in TRS): -pi to come back to "unusual" ref system (X,Y,-Z)
+        fPhotCnt++;                                                                           //increment counter of photon candidates
+      }
     }
   }//clusters loop
-  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>1)    {pTrk->SetHMPIDsignal   (kMipDistCut); return;}                 //closest cluster with enough charge is still too far from intersection
-                 pTrk->SetHMPIDcluIdx(chId,mipId); 
-  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
+  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);                              //error squared 
+  pTrk->SetHMPIDchi2(fCkovSigma2);                                                            //errors squared 
 
 }//ThetaCerenkov()
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-Double_t AliHMPIDRecon::FindPhotCkov(Double_t cluX,Double_t cluY)
+Bool_t AliHMPIDRecon::FindPhotCkov(Double_t cluX,Double_t cluY,Double_t &thetaCer,Double_t &phiCer)
 {
 // Finds Cerenkov angle  for this photon candidate
 // Arguments: cluX,cluY - position of cadidate's cluster  
-//   Returns: Cerenkov angle 
+// Returns: Cerenkov angle 
 
-  TVector2 pos(cluX,cluY); Double_t cluR=(pos-fTrkPos).Mod();  Double_t phi=FindPhotPhi(cluX,cluY);      
-  Double_t ckov1=0,ckov2=0.75;
-  const Double_t kTol=0.05; 
+  TVector3 dirCkov;
+  
+  Double_t zRad=-0.5*AliHMPIDRecon::fgkRadThick
+                    -AliHMPIDRecon::fgkWinThick
+                    -AliHMPIDRecon::fgkGapThick;                    //z position of middle of RAD
+
+  TVector3 rad(fTrkPos.X(),fTrkPos.Y(),zRad);                       //impact point at middle of RAD
+  TVector3  pc(cluX,cluY,0);                                        //mip at PC: z=0 @ 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
+    
+  Double_t ckov1=0,ckov2=0.75+fTrkDir.Theta();                      //start to find theta cerenkov in DRS
+  const Double_t kTol=0.01; 
   Int_t iIterCnt = 0;
   while(1){
-    if(iIterCnt>=50) return -1;
+    if(iIterCnt>=50) return kFALSE;
     Double_t ckov=0.5*(ckov1+ckov2);
-    Double_t dist=cluR-TracePhot(xRad,yRad,ckov,phi,pos); iIterCnt++;   //get distance between trial point and cluster position
+    dirCkov.SetMagThetaPhi(1,ckov,phi);
+    TVector2 posC=TraceForward(dirCkov);                      //trace photon with actual angles
+    Double_t dist=cluR-(posC-fTrkPos).Mod();                  //get distance between trial point and cluster position
+    if(posC.X()==-999) dist = - 999;                          //total reflection problem
+    iIterCnt++;                                               //counter step
     if     (dist> kTol) ckov1=ckov;                           //cluster @ larger ckov 
     else if(dist<-kTol) ckov2=ckov;                           //cluster @ smaller ckov
-    else return ckov;                                         //precision achived
+    else{                                                     //precision achived: ckov in DRS found
+      dirCkov.SetMagThetaPhi(1,ckov,phi);                     //
+      RecPhot(dirCkov,thetaCer,phiCer);                       //find ckov (in TRS:the effective Cherenkov angle!)
+      return kTRUE;
+    }
   }
 }//FindPhotTheta()
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-Double_t AliHMPIDRecon::FindPhotPhi(Double_t cluX,Double_t cluY)
+TVector2 AliHMPIDRecon::TraceForward(TVector3 dirCkov)const
 {
-// Finds phi angle og photon candidate by considering the cluster's position  of this candudate w.r.t track position
-  
-//  Double_t emiss=0; 
-//  return fPhotPhi[fPhotCnt]=TMath::ATan2(cluY-fTrkPos.Y()-emiss*TMath::Tan(fTrkDir.Theta())*TMath::Sin(fTrkDir.Phi()),
-//                                         cluX-fTrkPos.X()-emiss*TMath::Tan(fTrkDir.Theta())*TMath::Cos(fTrkDir.Phi()));
-  return fPhotPhi[fPhotCnt]=TMath::ATan2(cluY-yRad,cluX-xRad)-(TMath::Pi()+fTrkDir.Phi());
+  //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);
+  if(dirCkov.Theta() > TMath::ASin(1./fRadNmean))  return pos;      //total refraction on WIN-GAP boundary
+  Double_t zRad=-0.5*AliHMPIDRecon::fgkRadThick
+                    -AliHMPIDRecon::fgkWinThick
+                    -AliHMPIDRecon::fgkGapThick;                    //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,-fgkWinThick-fgkGapThick);              //go to RAD-WIN boundary  
+  Refract  (dirCkov,         fRadNmean,fgkWinIdx    );              //RAD-WIN refraction
+  Propagate(dirCkov,posCkov,            -fgkGapThick);              //go to WIN-GAP boundary
+  Refract  (dirCkov,         fgkWinIdx,fgkGapIdx    );              //WIN-GAP refraction
+  Propagate(dirCkov,posCkov,                       0);              //go to PC
+  pos.Set(posCkov.X(),posCkov.Y());
+  return pos;
+}//TraceForward()
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+void AliHMPIDRecon::RecPhot(TVector3 dirCkov,Double_t &thetaCer,Double_t &phiCer)
+{
+  //Theta Cerenkov reconstruction 
+  // Arguments: (x,y) of initial point in LORS, dirCkov photon vector in LORS
+  //   Returns: thetaCer theta cerenkov reconstructed
+//  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 mrot=mtheta*mphi;
+  TVector3 dirCkovTRS;
+  dirCkovTRS=mrot*dirCkov;
+  phiCer  = dirCkovTRS.Phi();                                          //actual value of the phi of the photon
+  thetaCer= dirCkovTRS.Theta();                                        //actual value of thetaCerenkov of the photon
 }
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 Double_t AliHMPIDRecon::FindRingArea(Double_t ckovAng)const
 {
 // Find area inside the cerenkov ring which lays inside PCs
-// Arguments: ckovThe - cernkov    
+// Arguments: ckovAng - cerenkov angle    
 //   Returns: area of the ring in cm^2 for given theta ckov
    
-  
-  TVector2 pos1,pos2;
-  
   const Int_t kN=100;
   Double_t area=0;
   for(Int_t i=0;i<kN;i++){
-    TracePhot(xRad,yRad,ckovAng,Double_t(TMath::TwoPi()*i    /kN),pos1);//trace this photon 
-    TracePhot(xRad,yRad,ckovAng,Double_t(TMath::TwoPi()*(i+1)/kN),pos2);//trace this photon 
-    area+=(pos1-fTrkPos)*(pos2-fTrkPos);      
+    TVector2 pos1=TracePhot(ckovAng,Double_t(TMath::TwoPi()*i    /kN));//trace this photon 
+    TVector2 pos2=TracePhot(ckovAng,Double_t(TMath::TwoPi()*(i+1)/kN));//trace the next photon 
+    area+=(pos1-fTrkPos)*(pos2-fTrkPos);                               //add area of the triangle... 
   }
   return area;
 }//FindRingArea()
@@ -163,10 +211,10 @@ Double_t AliHMPIDRecon::FindRingCkov(Int_t)
   
   for(Int_t i=0;i<fPhotCnt;i++){//candidates loop
     if(fPhotFlag[i] == 2){
-      if(fPhotCkov[i]<=0) continue;//?????????????????Flag photos = 2 may imply CkovEta = 0?????????????? 
-      if(fPhotCkov[i]<ckovMin) ckovMin=fPhotCkov[i];  //find max and min Theta ckov from all candidates within probable window
+      if(fPhotCkov[i]<ckovMin) ckovMin=fPhotCkov[i];                         //find max and min Theta ckov from all candidates within probable window
       if(fPhotCkov[i]>ckovMax) ckovMax=fPhotCkov[i]; 
-      weightThetaCerenkov += fPhotCkov[i]*fPhotWei[i];   wei += fPhotWei[i];                 //collect weight as sum of all candidate weghts   
+      weightThetaCerenkov += fPhotCkov[i]*fPhotWei[i];
+      wei += fPhotWei[i];                                                    //collect weight as sum of all candidate weghts   
       
       sigma2 += 1./Sigma2(fPhotCkov[i],fPhotPhi[i]);
     }
@@ -185,6 +233,9 @@ Int_t AliHMPIDRecon::FlagPhot(Double_t ckov)
 // Arguments: ckov- value of most probable ckov angle for track as returned by HoughResponse()
 //   Returns: number of photon candidates happened to be inside the window
 
+// Photon Flag:  Flag = 0 initial set; 
+//               Flag = 1 good candidate (charge compatible with photon); 
+//               Flag = 2 photon used for the ring;
   
   Int_t steps = (Int_t)((ckov )/ fDTheta); //how many times we need to have fDTheta to fill the distance between 0  and thetaCkovHough
 
@@ -204,30 +255,22 @@ Int_t AliHMPIDRecon::FlagPhot(Double_t ckov)
   return iInsideCnt;
 }//FlagPhot()
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-Double_t AliHMPIDRecon::TracePhot(Double_t x,Double_t y,Double_t ckovThe,Double_t ckovPhi,TVector2 &pos)const
+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, [rad]  (warning: not photon theta and phi)     
+// Arguments: ckovThe,ckovPhi- photon ckov angles in DRS, [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  posCkov(fTrkPos.X(),fTrkPos.Y(),-0.5*fgkRadThick-fgkWinThick-fgkGapThick);   //RAD: photon position is track position @ middle of RAD 
-  TVector3  posCkov(x,y,-0.5*fgkRadThick-fgkWinThick-fgkGapThick);                                                  //RAD: photon position is track position @ middle of RAD 
-  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
-//  dirCkov.SetPhi(ckovPhi);  
-  if(dirCkov.Theta() > TMath::ASin(1./fRadNmean)) return -999;//total refraction on WIN-GAP boundary
-  Propagate(dirCkov,posCkov,-fgkWinThick-fgkGapThick);                           //go to RAD-WIN boundary  
-  Refract  (dirCkov,         fRadNmean,fgkWinIdx    );              //RAD-WIN refraction
-  Propagate(dirCkov,posCkov,            -fgkGapThick);              //go to WIN-GAP boundary
-  Refract  (dirCkov,         fgkWinIdx,fgkGapIdx    );              //WIN-GAP refraction
-  Propagate(dirCkov,posCkov,                       0);   //go to PC
-  pos.Set(posCkov.X(),posCkov.Y());
-  return (pos-fTrkPos).Mod();
-}//TracePhoton()
+  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);
+}//TracePhot()
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-void AliHMPIDRecon::Propagate(const TVector3 &dir,TVector3 &pos,Double_t z)const
+void AliHMPIDRecon::Propagate(const TVector3 dir,TVector3 &pos,Double_t z)const
 {
 // Finds an intersection point between a line and XY plane shifted along Z.
 // Arguments:  dir,pos   - vector along the line and any point of the line
index 1f60a45..f2060c7 100644 (file)
@@ -26,17 +26,21 @@ public :
     virtual ~AliHMPIDRecon()                                                          {}
 
   
-  void     CkovAngle    (Double_t xRa,Double_t yRa,AliESDtrack *pTrk,TClonesArray *pCluLst,Double_t nmean); //reconstructed Theta Cerenkov
-  Double_t FindPhotCkov (Double_t cluX,Double_t cluY                                        );     //find ckov angle for single photon candidate
-  Double_t FindPhotPhi  (Double_t cluX,Double_t cluY                                        );     //find phi angle for single photon candidate
+  void     CkovAngle    (AliESDtrack *pTrk,TClonesArray *pCluLst,Double_t nmean             );     //reconstructed Theta Cerenkov
+  Bool_t   FindPhotCkov (Double_t cluX,Double_t cluY,Double_t &thetaCer,Double_t &phiCer    );     //find ckov angle for single photon candidate
   Double_t FindRingCkov (Int_t iNclus                                                       );     //best ckov for ring formed by found photon candidates
   Double_t FindRingArea (Double_t ckov                                                      )const;//estimated area of ring in cm^2
   Int_t    FlagPhot     (Double_t ckov                                                      );     //is photon ckov near most probable track ckov
   Double_t HoughResponse(                                                                   );     //most probable track ckov angle
-  void     Propagate    (const TVector3 &dir,      TVector3 &pos,Double_t z                 )const;//propagate photon alogn the line  
+  void     Propagate    (const TVector3  dir,      TVector3 &pos,Double_t z                 )const;//propagate photon alogn the line  
   void     Refract      (      TVector3 &dir,                    Double_t n1,    Double_t n2)const;//refract photon on the boundary
-  Double_t TracePhot    (Double_t x,Double_t y,Double_t ckovTh,Double_t ckovPh,TVector2 &pos)const;//trace photon created by track to PC 
-  void     SetTrack     (Double_t x,Double_t y,Double_t theta,Double_t phi                  ){fTrkDir.SetMagThetaPhi(1,theta,phi);  fTrkPos.Set(x,y);}//set track
+  TVector2 TracePhot    (Double_t ckovTh,Double_t ckovPh                                    )const;//trace photon created by track to PC 
+  TVector2 TraceForward (TVector3 dirCkov                                                   )const;//tracing forward a photon from (x,y) to PC
+  void     RecPhot      (TVector3 dirCkov,Double_t &thetaCer,Double_t &phiCer               );     //theta,phi cerenkov reconstructed
+  void     SetTrack     (Double_t xRad,Double_t yRad,Double_t theta,Double_t phi            )
+                                {fTrkDir.SetMagThetaPhi(1,theta,phi);  fTrkPos.Set(xRad,yRad);}    //set track parameter at RAD
+  void     SetImpPC     (Double_t xPc,Double_t yPc                                          )
+                                {fPc.Set(xPc,yPc);}                                                //set track impact to PC 
   Double_t SigLoc       (Double_t ckovTh,Double_t ckovPh,Double_t beta                      )const;//error due to cathode segmetation
   Double_t SigGeom      (Double_t ckovTh,Double_t ckovPh,Double_t beta                      )const;//error due to unknown photon origin
   Double_t SigCrom      (Double_t ckovTh,Double_t ckovPh,Double_t beta                      )const;//error due to unknonw photon energy
@@ -61,8 +65,9 @@ protected:
   Float_t fDTheta;                            // Step for sliding window
   Float_t fWindowWidth;                       // Hough width of sliding window
   
-  TVector3 fTrkDir;                           //track direction in LORS
-  TVector2 fTrkPos;                           //track positon in LORS at the middle of radiator
+  TVector3 fTrkDir;                           //track direction in LORS at RAD
+  TVector2 fTrkPos;                           //track positon in LORS at RAD
+  TVector2 fPc;                               //track position at PC
   ClassDef(AliHMPIDRecon,0)
 };
 
index e35f860..f6f7223 100644 (file)
@@ -55,47 +55,47 @@ Int_t AliHMPIDTracker::PropagateBack(AliESD *pEsd)
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 Int_t AliHMPIDTracker::Recon(AliESD *pEsd,TObjArray *pCluAll,TObjArray *pNmean)
 {
-// Interface callback methode invoked by AliRecontruction::RunTracking() during tracking after TOF. It's done just once per event
+// Interface callback method invoked by AliRecontruction::RunTracking() during tracking after TOF. It's done just once per event
 // Arguments: pEsd - pointer to Event Summary Data class instance which contains a list of tracks
 //   Returns: error code, 0 if no errors   
   Int_t iNtracks=pEsd->GetNumberOfTracks();  AliDebugClass(1,Form("Start with %i tracks",iNtracks));
   
-  
   AliHMPIDRecon recon;                                                                       //instance of reconstruction class, nothing important in ctor
-  Float_t xRa,yRa;
+  Double_t xPc,yPc;
   for(Int_t iTrk=0;iTrk<iNtracks;iTrk++){                                                        //ESD tracks loop
     AliESDtrack *pTrk = pEsd->GetTrack(iTrk);                                                    //get next reconstructed track    
-    Int_t cham=IntTrkCha(pTrk,xRa,yRa);                                                          //get chamber intersected by thie track 
+    Int_t cham=IntTrkCha(pTrk,xPc,yPc);                                                          //get chamber intersected by this track 
     if(cham<0) continue;                                                                         //no intersection at all, go after next track
-    Double_t nmean=((TF1*)pNmean->At(3*cham))->Eval(pEsd->GetTimeStamp());                       //C6F14 Nmean for this chamber   
-    recon.CkovAngle((Double_t)xRa,(Double_t)yRa,pTrk,(TClonesArray *)pCluAll->At(cham),nmean);                               //search for Cerenkov angle for this track
+    Double_t nmean=((TF1*)pNmean->At(3*cham))->Eval(pEsd->GetTimeStamp());                       //C6F14 Nmean for this chamber
+    recon.SetImpPC(xPc,yPc);                                                                     //store track impact to PC
+    recon.CkovAngle(pTrk,(TClonesArray *)pCluAll->At(cham),nmean);                               //search for Cerenkov angle of this track
   }                                                                                              //ESD tracks loop
   AliDebugClass(1,"Stop pattern recognition");
   return 0; // error code: 0=no error;
 }//PropagateBack()
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-Int_t AliHMPIDTracker::IntTrkCha(AliESDtrack *pTrk,Float_t &xRa,Float_t &yRa)
+Int_t AliHMPIDTracker::IntTrkCha(AliESDtrack *pTrk,Double_t &xToPc,Double_t &yToPc)
 {
 // Static method to find intersection in between given track and HMPID chambers
 // Arguments: pTrk    - ESD track 
-//            xRa,yRa - track intersection with the middle of radiator, LORS
+//            xPc,yPc - track intersection with PC, LORS
 //   Returns: intersected chamber ID or -1
   
   AliHMPIDParam *pParam=AliHMPIDParam::Instance();
-  Float_t xPc=0,yPc=0,theta=0,phi=0;                                                            //track intersection point and angles, LORS  
-  for(Int_t i=AliHMPIDDigit::kMinCh;i<=AliHMPIDDigit::kMaxCh;i++){                                                                       //chambers loop
+  Float_t xRa=0,yRa=0,xPc=0,yPc=0,theta=0,phi=0;                                                //track intersection point and angles, LORS  
+  for(Int_t i=AliHMPIDDigit::kMinCh;i<=AliHMPIDDigit::kMaxCh;i++){                              //chambers loop
     Double_t p1[3],n1[3]; pParam->Norm(i,n1); pParam->Lors2Mars(i,0,0,p1,AliHMPIDParam::kRad);  //point & norm  for RAD
     Double_t p2[3],n2[3]; pParam->Norm(i,n2); pParam->Lors2Mars(i,0,0,p2,AliHMPIDParam::kPc);   //point & norm  for PC
       
     if(pTrk->Intersect(p1,n1,-GetBz())==kFALSE) continue;                                       //try to intersect track with the middle of radiator
     if(pTrk->Intersect(p2,n2,-GetBz())==kFALSE) continue;                                       //try to intersect track with PC
       
-    pParam->Mars2LorsVec(i,n1,theta,phi);                                                       //track angles
+    pParam->Mars2LorsVec(i,n1,theta,phi);                                                       //track angles at RAD
     pParam->Mars2Lors   (i,p1,xRa,yRa);                                                         //TRKxRAD position
     pParam->Mars2Lors   (i,p2,xPc,yPc);                                                         //TRKxPC position
-      
-    if(AliHMPIDDigit::IsInside(xPc,yPc)==kFALSE) continue;                                      //not in active area  
-    pTrk->SetHMPIDtrk      (xPc,yPc,theta,phi);                                                 //store track intersection info
+    xToPc=(Double_t)xPc;yToPc=(Double_t)yPc;                                                    //conversion float->double only
+    if(AliHMPIDDigit::IsInside(xPc,yPc,pParam->DistCut())==kFALSE) continue;             //not in active area  
+    pTrk->SetHMPIDtrk      (xRa,yRa,theta,phi);                                                 //store track intersection info
     pTrk->SetHMPIDcluIdx   (i,0);
     return i;
   }                                                                                             //chambers loop
index 9016348..4d722ac 100644 (file)
@@ -21,7 +21,7 @@ public:
          Int_t       RefitInward    (AliESD *                   )       {return 0;} //pure virtual from AliTracker 
          void        UnloadClusters (                           )       {         } //pure virtual from AliTracker 
 //private part  
-  static Int_t       IntTrkCha(AliESDtrack *pTrk,Float_t &x,Float_t &y);                    //find track-chamber intersection, retuns chamber ID
+  static Int_t       IntTrkCha(AliESDtrack *pTrk,Double_t &xPc,Double_t &yPc);              //find track-PC intersection, retuns chamber ID
   static Int_t       Recon    (AliESD *pEsd,TObjArray *pCluAll,TObjArray *pNmean=0);        //do actual job, returns status code  
 protected:
   ClassDef(AliHMPIDTracker,0)
index 69e1e5a..86a7f11 100644 (file)
@@ -87,7 +87,7 @@ void SimHits(AliESD *pEsd, TClonesArray *pHits)
   
   AliHMPIDRecon rec;
   
-  Int_t hc=0; TVector2 pos;
+  Int_t hc=0; 
   for(Int_t iTrk=0;iTrk<pEsd->GetNumberOfTracks();iTrk++){//tracks loop
     AliESDtrack *pTrk=pEsd->GetTrack(iTrk);
     Float_t xRa,yRa;
@@ -101,15 +101,10 @@ void SimHits(AliESD *pEsd, TClonesArray *pHits)
     rec.SetTrack(xRa,yRa,theta,phi); 
     
     if(!AliHMPIDDigit::IsInDead(xPc,yPc)) new((*pHits)[hc++]) AliHMPIDHit(ch,200e-9,kProton  ,iTrk,xPc,yPc);                 //mip hit
-    /*
-    for(int i=0;i<4;i++) { 
-      Float_t x=gRandom->Rndm()*130;Float_t y=gRandom->Rndm()*126;
-      if(!AliHMPIDDigit::IsInDead(x,y)) new((*pHits)[hc++]) AliHMPIDHit(ch,7.5e-9,kFeedback,iTrk,x,y); //bkg hits 4 per track
-    }
-    */
     Int_t nPhots = (Int_t)(20.*TMath::Power(TMath::Sin(ckov),2)/TMath::Power(TMath::Sin(TMath::ACos(1./1.292)),2));
     for(int i=0;i<nPhots;i++){
-      rec.TracePhot(ckov,gRandom->Rndm()*TMath::TwoPi(),pos);
+      TVector2 pos;
+      pos=rec.TracePhot(ckov,gRandom->Rndm()*TMath::TwoPi());
       if(!AliHMPIDDigit::IsInDead(pos.X(),pos.Y())) new((*pHits)[hc++]) AliHMPIDHit(ch,7.5e-9,kCerenkov,iTrk,pos.X(),pos.Y());
     }                      //photon hits  
   }//tracks loop    
@@ -162,21 +157,16 @@ void DrawEvt(TClonesArray *pHitLst,TObjArray *pDigLst,TObjArray *pCluLst,AliESD
     Int_t ch=pTrk->GetHMPIDcluIdx();
     if(ch<0) continue; //this track does not hit HMPID
     ch/=1000000; 
-    Float_t th,ph,xPc,yPc; pTrk->GetHMPIDtrk(xPc,yPc,th,ph);  //get info on current track
-    pTxC[ch]->SetNextPoint(xPc,yPc);                          //add new intersection point
+    Float_t th,ph,xRad,yRad; pTrk->GetHMPIDtrk(xRad,yRad,th,ph);//get info on current track
+//    pTxC[ch]->SetNextPoint(xPc,yPc);                          //add new intersection point  TEMPORARLY DISABLED...no more available in ESD!
     
     Float_t ckov=pTrk->GetHMPIDsignal();  Float_t err=TMath::Sqrt(pTrk->GetHMPIDchi2());
     if(ckov>0){
       Printf("theta %f phi %f ckov %f",th*TMath::RadToDeg(),ph*TMath::RadToDeg(),ckov);
-      rec.SetTrack(xPc,yPc,th,ph+TMath::Pi());
-      TVector2 pos;
-       Double_t allGapz=rec.fgkWinThick+0.5*rec.fgkRadThick+rec.fgkGapThick;   //to semplify, the (x,y) are calculted from (xPc,yPc) back to radiator in straight line (Bz=0)
-      TVector3 dir(0,0,-1);TVector3 miprad(xPc,yPc,allGapz);
-      AliHMPIDRecon::Propagate(dir,miprad,0);
-      Double_t xRad=miprad.X();
-      Double_t yRad=miprad.Y();
-      for(int j=0;j<100;j++){ 
-        rec.TracePhot(xRad,yRad,ckov,j*0.0628,pos);
+      rec.SetTrack(xRad,yRad,th,ph+TMath::Pi());
+      for(Int_t j=0;j<100;j++){ 
+        TVector2 pos;
+        pos=rec.TracePhot(ckov,j*0.0628);
        if(!AliHMPIDDigit::IsInDead(pos.X(),pos.Y())) pRin[ch]->SetNextPoint(pos.X(),pos.Y());
       }      
     }