Better method to find area and length of a ring (needed to normalize n. of photons)
authordibari <dibari@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 14 Nov 2007 08:28:51 +0000 (08:28 +0000)
committerdibari <dibari@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 14 Nov 2007 08:28:51 +0000 (08:28 +0000)
HMPID/AliHMPIDRecon.cxx
HMPID/AliHMPIDRecon.h

index 44e125f..cf23978 100644 (file)
@@ -116,18 +116,20 @@ void AliHMPIDRecon::CkovAngle(AliESDtrack *pTrk,TClonesArray *pCluLst,Double_t n
   }//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 
+  Int_t iNrec=FlagPhot(HoughResponse());                                                      //flag photons according to individual theta ckov with respect to most probable
+  pTrk->SetHMPIDmip(mipX,mipY,mipQ,iNrec);                                                    //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
+  if(iNrec<1){
+    pTrk->SetHMPIDsignal(kNoPhotAccept);                                                      //no photon candidates are accepted
   }
   else {
-    pTrk->SetHMPIDsignal(FindRingCkov(pCluLst->GetEntries()));                                //find best Theta ckov for ring i.e. track
-    pTrk->SetHMPIDchi2(fCkovSigma2);                                                          //errors squared
+    Double_t thetaC = FindRingCkov(pCluLst->GetEntries());                                    //find the best reconstructed theta Cherenkov
+//    FindRingGeom(thetaC,2);
+    pTrk->SetHMPIDsignal(thetaC);                                                             //store theta Cherenkov
+    pTrk->SetHMPIDchi2(fCkovSigma2);                                                          //store errors squared
   }
 
   DeleteVars();
@@ -164,7 +166,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 +192,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 +210,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 +255,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
 {
@@ -340,14 +369,13 @@ 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 +427,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);
index e7d5c5b..ba0a12c 100644 (file)
@@ -31,7 +31,7 @@ public :
   void     CkovAngle    (AliESDtrack *pTrk,TClonesArray *pCluLst,Double_t nmean, Double_t qthre);  //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 ckovAng                                                   )const;//estimated area of delta ring in cm^2 to weight Hough Transform
+  void     FindRingGeom (Double_t ckovAng,Int_t level=1                                     );     //estimated area of ring in cm^2 and portion accepted by geometry
   TVector2 IntWithEdge  (TVector2 p1,TVector2 p2                                            )const;//find intercection between plane and lines of 2 thetaC
   Int_t    FlagPhot     (Double_t ckov                                                      );     //is photon ckov near most probable track ckov
   Double_t HoughResponse(                                                                   );     //most probable track ckov angle
@@ -39,9 +39,14 @@ public :
   void     Refract      (      TVector3 &dir,                    Double_t n1,    Double_t n2)const;//refract photon on the boundary
   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     Lors2Trs     (TVector3 dirCkov,Double_t &thetaCer,Double_t &phiCer               )const;//LORS to TRS 
+  void     Trs2Lors     (TVector3 dirCkov,Double_t &thetaCer,Double_t &phiCer               )const;//TRS to LORS
   TVector2 GetMip       (                                                                   ) 
                         {return fMipPos;}                                                          //mip coordinates
+  Double_t GetRingArea  (                                                                   )
+                        {return fRingArea;}                                                        //area of the current ring in cm^2 
+  Double_t GetRingAcc   (                                                                   )
+                        {return fRingAcc;}                                                         //portion of the ring ([0,1]) accepted by geometry.To scale n. of photons 
   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                                          )
@@ -66,18 +71,20 @@ protected:
   Float_t   fDTheta;                            // Step for sliding window
   Float_t   fWindowWidth;                       // Hough width of sliding window
   
-  TVector3  fTrkDir;                            //track direction in LORS at RAD
-  TVector2  fTrkPos;                            //track positon in LORS at RAD
-  TVector2  fMipPos;                            //mip positon for a given track
-  TVector2  fPc;                                //track position at PC
+  Double_t  fRingArea;                          // area of a given ring
+  Double_t  fRingAcc;                           // fraction of the ring accepted by geometry
+  TVector3  fTrkDir;                            // track direction in LORS at RAD
+  TVector2  fTrkPos;                            // track positon in LORS at RAD
+  TVector2  fMipPos;                            // mip positon for a given track
+  TVector2  fPc;                                // track position at PC
   
-  AliHMPIDParam *fParam;                        //Pointer to AliHMPIDParam
+  AliHMPIDParam *fParam;                        // Pointer to AliHMPIDParam
   
 private:
   AliHMPIDRecon(const AliHMPIDRecon& r);              //dummy copy constructor
   AliHMPIDRecon &operator=(const AliHMPIDRecon& r);   //dummy assignment operator
 //
-  ClassDef(AliHMPIDRecon,0)
+  ClassDef(AliHMPIDRecon,1)
 };
 
 #endif // #ifdef AliHMPIDRecon_cxx