#include <TH1D.h> //HoughResponse()
#include <TClonesArray.h> //CkovAngle()
#include <AliESDtrack.h> //CkovAngle()
-
- Int_t fPhotCnt; // counter of photons candidate
- Int_t *fPhotFlag; // flags of photon candidates
- Double_t *fPhotCkov; // Ckov angles of photon candidates, [rad]
- Double_t *fPhotPhi; // phis of photons candidates, [rad]
- Double_t *fPhotWei; // weigths of photon candidates
- Double_t fCkovSigma2; // sigma2 of the reconstructed ring
-
- Bool_t fIsWEIGHT; // flag to consider weight procedure
- Float_t fDTheta; // Step for sliding window
- Float_t fWindowWidth; // Hough width of sliding window
-
- 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
-
+#include <AliESDfriendTrack.h> //CkovAngle()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
AliHMPIDRecon::AliHMPIDRecon():
- TTask("RichRec","RichPat"),
+ TNamed("RichRec","RichPat"),
fPhotCnt(-1),
fPhotFlag(0x0),
+ fPhotClusIndex(0x0),
fPhotCkov(0x0),
fPhotPhi(0x0),
fPhotWei(0x0),
//..
//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],
-
+
+ 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;
+
+ 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
+
+ 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
- pTrk->SetHMPIDmip(mipX,mipY,mipQ,iNrec); //store mip info
+
+//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(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
+ if(iNrec<nMinPhotAcc){
+ pTrk->SetHMPIDsignal(kNoPhotAccept); //no photon candidates are accepted
+ return;
}
- else {
- Double_t thetaC = FindRingCkov(pCluLst->GetEntries()); //find the best reconstructed theta Cherenkov
+
+ 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); //store theta Cherenkov
- pTrk->SetHMPIDchi2(fCkovSigma2); //store errors squared
- }
-
+ pTrk->SetHMPIDsignal(thetaC+occupancy); //store theta Cherenkov and chmaber occupancy
+ pTrk->SetHMPIDchi2(fCkovSigma2); //store experimental ring angular resolution squared
+
DeleteVars();
}//CkovAngle()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
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()
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
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
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(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;
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff