#include <TH1D.h> //HoughResponse()
#include <TClonesArray.h> //CkovAngle()
#include <AliESDtrack.h> //CkovAngle()
+#include <AliESDfriendTrack.h> //CkovAngle()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-AliHMPIDRecon::AliHMPIDRecon():TTask("RichRec","RichPat")
+AliHMPIDRecon::AliHMPIDRecon():
+ TTask("RichRec","RichPat"),
+ fPhotCnt(-1),
+ fPhotFlag(0x0),
+ fPhotClusIndex(0x0),
+ fPhotCkov(0x0),
+ fPhotPhi(0x0),
+ fPhotWei(0x0),
+ fCkovSigma2(0),
+ fIsWEIGHT(kFALSE),
+ fDTheta(0.001),
+ fWindowWidth(0.045),
+ fRingArea(0),
+ fRingAcc(0),
+ fTrkDir(0,0,1), // Just for test
+ fTrkPos(30,40), // Just for test
+ fMipPos(0,0),
+ fPc(0,0),
+ fParam(AliHMPIDParam::Instance())
{
//..
//init of data members
//..
- fPhotCnt = -1;
- fPhotFlag = 0x0;
- fPhotCkov = 0x0;
- fPhotPhi = 0x0;
- fPhotWei = 0x0;
- fCkovSigma2 = 0;
- fIsWEIGHT = kFALSE;
- fDTheta = 0.001;
- fWindowWidth = 0.045;
- fTrkDir = TVector3(0,0,1); // init just for test
- fTrkPos = TVector2(30,40); // init just for test
-
- AliHMPIDParam *pParam=AliHMPIDParam::Instance();
- fParam = pParam;
-
fParam->SetRefIdx(fParam->MeanIdxRad()); // initialization of ref index to a default one
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
//..
//Init some variables
//..
- if(n<0) return;
+ if(n<=0) return;
fPhotFlag = new Int_t[n];
+ fPhotClusIndex = new Int_t[n];
fPhotCkov = new Double_t[n];
fPhotPhi = new Double_t[n];
fPhotWei = new Double_t[n];
//..
//Delete variables
//..
- delete fPhotFlag;
- delete fPhotCkov;
- delete fPhotPhi;
- delete fPhotWei;
+ delete [] fPhotFlag;
+ delete [] fPhotClusIndex;
+ delete [] fPhotCkov;
+ delete [] fPhotPhi;
+ delete [] fPhotWei;
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-void AliHMPIDRecon::CkovAngle(AliESDtrack *pTrk,TClonesArray *pCluLst,Double_t nmean,Double_t qthre)
+void AliHMPIDRecon::CkovAngle(AliESDtrack *pTrk,TClonesArray *pCluLst,Int_t index,Double_t nmean,Float_t xRa,Float_t yRa)
{
// Pattern recognition method based on Hough transform
// Arguments: pTrk - track for which Ckov angle is to be found
// pCluLst - list of clusters for this chamber
// Returns: - track ckov angle, [rad],
+
+ AliESDfriendTrack *pFriendTrk = (AliESDfriendTrack*)pTrk->GetFriendTrack();
+
+ const Int_t nMinPhotAcc = 3; // Minimum number of photons required to perform the pattern recognition
+
Int_t nClusTot = pCluLst->GetEntries();
if(nClusTot>fParam->MultCut()) fIsWEIGHT = kTRUE; // offset to take into account bkg in reconstruction
else fIsWEIGHT = kFALSE;
InitVars(nClusTot);
- Float_t xRa,yRa,th,ph;
- pTrk->GetHMPIDtrk(xRa,yRa,th,ph); //initialize this track: th and ph angles at middle of RAD
+ Float_t xPc,yPc,th,ph;
+ pTrk->GetHMPIDtrk(xPc,yPc,th,ph); //initialize this track: th and ph angles at middle of RAD
SetTrack(xRa,yRa,th,ph);
fParam->SetRefIdx(nmean);
- Float_t dMin=999,mipX=-1,mipY=-1;Int_t chId=-1,mipId=-1,mipQ=-1;
- fPhotCnt=0;
+ Float_t mipX=-1,mipY=-1;
+ Int_t chId=-1,mipQ=-1,sizeClu = -1;
+
+ fPhotCnt=0;
+
for (Int_t iClu=0; iClu<pCluLst->GetEntriesFast();iClu++){//clusters loop
AliHMPIDCluster *pClu=(AliHMPIDCluster*)pCluLst->UncheckedAt(iClu); //get pointer to current cluster
+ 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
- }
+ 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
+
+ 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);
- if(fPhotCnt<=3) pTrk->SetHMPIDsignal(kNoPhotAccept); //no reconstruction with <=3 photon candidates
- Int_t iNrec=FlagPhot(HoughResponse()); //flag photons according to individual theta ckov with respect to most probable
+
+ TClonesArray *pPhotCluLst = pFriendTrk->GetHmpPhotClus();
+
+//PATTERN RECOGNITION STARTED:
+
+ Int_t iNrec=FlagPhot(HoughResponse(),pCluLst,pPhotCluLst); //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(iNrec<1){
pTrk->SetHMPIDsignal(kNoPhotAccept); //no photon candidates are accepted
+ return;
}
- else {
- Double_t thetaC = FindRingCkov(pCluLst->GetEntries()); //find the best reconstructed theta Cherenkov
+
+ 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
- }
-
+ pTrk->SetHMPIDsignal(thetaC); //store theta Cherenkov
+ pTrk->SetHMPIDchi2(fCkovSigma2); //store errors squared
+
DeleteVars();
}//CkovAngle()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
return weightThetaCerenkov;
}//FindCkovRing()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-Int_t AliHMPIDRecon::FlagPhot(Double_t ckov)
+Int_t AliHMPIDRecon::FlagPhot(Double_t ckov,TClonesArray *pCluLst, TClonesArray *pPhotCluLst)
{
// 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()
// Photon Flag: Flag = 0 initial set;
// Flag = 1 good candidate (charge compatible with photon);
// Flag = 2 photon used for the ring;
+ Int_t *PhotIndex = new Int_t[fPhotCnt];
Int_t steps = (Int_t)((ckov )/ fDTheta); //how many times we need to have fDTheta to fill the distance between 0 and thetaCkovHough
fPhotFlag[i] = 0;
if(fPhotCkov[i] >= tmin && fPhotCkov[i] <= tmax) {
fPhotFlag[i]=2;
+ PhotIndex[iInsideCnt]=fPhotClusIndex[i];
iInsideCnt++;
}
}
+
+ Int_t nPhot = 0;
+
+ for (Int_t iClu=0; iClu<pCluLst->GetEntriesFast();iClu++){//clusters loop
+ AliHMPIDCluster *pClu=(AliHMPIDCluster*)pCluLst->UncheckedAt(iClu); //get pointer to current cluster
+ for(Int_t j=0; j<iInsideCnt; j++){
+ if(iClu==PhotIndex[j]) {
+ new ((*pPhotCluLst)[nPhot++]) AliHMPIDCluster(*pClu); //add this raw cluster
+ }
+ }
+ }
+
+ delete [] PhotIndex;
+
return iInsideCnt;
+
}//FlagPhot()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
TVector2 AliHMPIDRecon::TracePhot(Double_t ckovThe,Double_t ckovPhi)const
return (Double_t)(locMax*fDTheta+0.5*fDTheta); //final most probable track theta ckov
}//HoughResponse()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ Double_t AliHMPIDRecon::FindRingExt(Double_t ckov,Int_t ch,Double_t xPc,Double_t yPc,Double_t thRa,Double_t phRa)
+{
+// 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);
+
+ Int_t nStep = 500;
+ Int_t nPhi = 0;
+
+ Int_t ipc,ipadx,ipady;
+
+ if(ckov>0){
+ SetTrack(xRa,yRa,thRa,phRa);
+ for(Int_t j=0;j<nStep;j++){
+ TVector2 pos; pos=TracePhot(ckov,j*TMath::TwoPi()/(Double_t)(nStep-1));
+ if(fParam->IsInDead(pos.X(),pos.Y())) continue;
+ fParam->Lors2Pad(pos.X(),pos.Y(),ipc,ipadx,ipady);
+ ipadx+=(ipc%2)*fParam->kPadPcX;
+ ipady+=(ipc/2)*fParam->kPadPcY;
+ if(fParam->IsDeadPad(ipadx,ipady,ch)) continue;
+ nPhi++;
+ }//point loop
+ return ((Double_t)nPhi/(Double_t)nStep);
+ }//if
+ return -1;
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff