From 2d1a9b21e4945b9a31708540aa99d2cedebea64a Mon Sep 17 00:00:00 2001 From: dibari Date: Wed, 14 Nov 2007 08:28:51 +0000 Subject: [PATCH] Better method to find area and length of a ring (needed to normalize n. of photons) --- HMPID/AliHMPIDRecon.cxx | 92 +++++++++++++++++++++++++++-------------- HMPID/AliHMPIDRecon.h | 23 +++++++---- 2 files changed, 77 insertions(+), 38 deletions(-) diff --git a/HMPID/AliHMPIDRecon.cxx b/HMPID/AliHMPIDRecon.cxx index 44e125f8677..cf23978fbd0 100644 --- a/HMPID/AliHMPIDRecon.cxx +++ b/HMPID/AliHMPIDRecon.cxx @@ -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;iIsInside(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); diff --git a/HMPID/AliHMPIDRecon.h b/HMPID/AliHMPIDRecon.h index e7d5c5b7a2d..ba0a12c5f70 100644 --- a/HMPID/AliHMPIDRecon.h +++ b/HMPID/AliHMPIDRecon.h @@ -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 -- 2.43.0