#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];
//
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-void AliHMPIDRecon::DeleteVars()
+void AliHMPIDRecon::DeleteVars()const
{
//..
//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;
+
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
+
+//PATTERN RECOGNITION STARTED:
+
+ 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
+ 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()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
weightThetaCerenkov += fPhotCkov[i]*fPhotWei[i];
wei += fPhotWei[i]; //collect weight as sum of all candidate weghts
- sigma2 += 1./Sigma2(fPhotCkov[i],fPhotPhi[i]);
+ sigma2 += 1./fParam->Sigma2(fTrkDir.Theta(),fTrkDir.Phi(),fPhotCkov[i],fPhotPhi[i]);
}
}//candidates loop
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()
// 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++;
}
}
+
+ 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]) {
+ AliHMPIDCluster *pClus = new AliHMPIDCluster(*pClu);
+ pTrk->AddCalibObject(pClus);
+ }
+ }
+ }
+
+ 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::Sigma2(Double_t ckovTh, Double_t ckovPh)const
-{
-// Analithical calculation of total error (as a sum of localization, geometrical and chromatic errors) on Cerenkov angle for a given Cerenkov photon
-// created by a given MIP. Fromulae according to CERN-EP-2000-058
-// Arguments: Cerenkov and azimuthal angles for Cerenkov photon, [radians]
-// dip and azimuthal angles for MIP taken at the entrance to radiator, [radians]
-// MIP beta
-// Returns: absolute error on Cerenkov angle, [radians]
-
- TVector3 v(-999,-999,-999);
- Double_t trkBeta = 1./(TMath::Cos(ckovTh)*fParam->GetRefIdx());
-
- if(trkBeta > 1) trkBeta = 1; //protection against bad measured thetaCer
- if(trkBeta < 0) trkBeta = 0.0001; //
-
- v.SetX(SigLoc (ckovTh,ckovPh,trkBeta));
- v.SetY(SigGeom(ckovTh,ckovPh,trkBeta));
- v.SetZ(SigCrom(ckovTh,ckovPh,trkBeta));
-
- return v.Mag2();
-}
-//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-Double_t AliHMPIDRecon::SigLoc(Double_t thetaC, Double_t phiC,Double_t betaM)const
-{
-// Analithical calculation of localization error (due to finite segmentation of PC) on Cerenkov angle for a given Cerenkov photon
-// created by a given MIP. Fromulae according to CERN-EP-2000-058
-// Arguments: Cerenkov and azimuthal angles for Cerenkov photon, [radians]
-// dip and azimuthal angles for MIP taken at the entrance to radiator, [radians]
-// MIP beta
-// Returns: absolute error on Cerenkov angle, [radians]
-
- Double_t phiDelta = phiC - fTrkDir.Phi();
-
- Double_t sint = TMath::Sin(fTrkDir.Theta());
- Double_t cost = TMath::Cos(fTrkDir.Theta());
- Double_t sinf = TMath::Sin(fTrkDir.Phi());
- Double_t cosf = TMath::Cos(fTrkDir.Phi());
- Double_t sinfd = TMath::Sin(phiDelta);
- Double_t cosfd = TMath::Cos(phiDelta);
- Double_t tantheta = TMath::Tan(thetaC);
-
- Double_t alpha =cost-tantheta*cosfd*sint; // formula (11)
- Double_t k = 1.-fParam->GetRefIdx()*fParam->GetRefIdx()+alpha*alpha/(betaM*betaM); // formula (after 8 in the text)
- if (k<0) return 1e10;
- Double_t mu =sint*sinf+tantheta*(cost*cosfd*sinf+sinfd*cosf); // formula (10)
- Double_t e =sint*cosf+tantheta*(cost*cosfd*cosf-sinfd*sinf); // formula (9)
-
- Double_t kk = betaM*TMath::Sqrt(k)/(fParam->GapThick()*alpha); // formula (6) and (7)
- Double_t dtdxc = kk*(k*(cosfd*cosf-cost*sinfd*sinf)-(alpha*mu/(betaM*betaM))*sint*sinfd); // formula (6)
- Double_t dtdyc = kk*(k*(cosfd*sinf+cost*sinfd*cosf)+(alpha* e/(betaM*betaM))*sint*sinfd); // formula (7) pag.4
-
- Double_t errX = 0.2,errY=0.25; //end of page 7
- return TMath::Sqrt(errX*errX*dtdxc*dtdxc + errY*errY*dtdyc*dtdyc);
-}
-//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-Double_t AliHMPIDRecon::SigCrom(Double_t thetaC, Double_t phiC,Double_t betaM)const
+ Double_t AliHMPIDRecon::FindRingExt(Double_t ckov,Int_t ch,Double_t xPc,Double_t yPc,Double_t thRa,Double_t phRa)
{
-// Analithical calculation of chromatic error (due to lack of knowledge of Cerenkov photon energy) on Cerenkov angle for a given Cerenkov photon
-// created by a given MIP. Fromulae according to CERN-EP-2000-058
-// Arguments: Cerenkov and azimuthal angles for Cerenkov photon, [radians]
-// dip and azimuthal angles for MIP taken at the entrance to radiator, [radians]
-// MIP beta
-// Returns: absolute error on Cerenkov angle, [radians]
-
- Double_t phiDelta = phiC - fTrkDir.Phi();
-
- Double_t sint = TMath::Sin(fTrkDir.Theta());
- Double_t cost = TMath::Cos(fTrkDir.Theta());
- Double_t cosfd = TMath::Cos(phiDelta);
- Double_t tantheta = TMath::Tan(thetaC);
-
- Double_t alpha =cost-tantheta*cosfd*sint; // formula (11)
- Double_t dtdn = cost*fParam->GetRefIdx()*betaM*betaM/(alpha*tantheta); // formula (12)
-
-// Double_t f = 0.00928*(7.75-5.635)/TMath::Sqrt(12.);
- Double_t f = 0.0172*(7.75-5.635)/TMath::Sqrt(24.);
-
- return f*dtdn;
-}//SigCrom()
-//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-Double_t AliHMPIDRecon::SigGeom(Double_t thetaC, Double_t phiC,Double_t betaM)const
-{
-// Analithical calculation of geometric error (due to lack of knowledge of creation point in radiator) on Cerenkov angle for a given Cerenkov photon
-// created by a given MIP. Formulae according to CERN-EP-2000-058
-// Arguments: Cerenkov and azimuthal angles for Cerenkov photon, [radians]
-// dip and azimuthal angles for MIP taken at the entrance to radiator, [radians]
-// MIP beta
-// Returns: absolute error on Cerenkov angle, [radians]
-
- Double_t phiDelta = phiC - fTrkDir.Phi();
-
- Double_t sint = TMath::Sin(fTrkDir.Theta());
- Double_t cost = TMath::Cos(fTrkDir.Theta());
- Double_t sinf = TMath::Sin(fTrkDir.Phi());
- Double_t cosfd = TMath::Cos(phiDelta);
- Double_t costheta = TMath::Cos(thetaC);
- Double_t tantheta = TMath::Tan(thetaC);
-
- Double_t alpha =cost-tantheta*cosfd*sint; // formula (11)
+// 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);
- Double_t k = 1.-fParam->GetRefIdx()*fParam->GetRefIdx()+alpha*alpha/(betaM*betaM); // formula (after 8 in the text)
- if (k<0) return 1e10;
+ Int_t nStep = 500;
+ Int_t nPhi = 0;
- Double_t eTr = 0.5*fParam->RadThick()*betaM*TMath::Sqrt(k)/(fParam->GapThick()*alpha); // formula (14)
- Double_t lambda = 1.-sint*sint*sinf*sinf; // formula (15)
-
- Double_t c1 = 1./(1.+ eTr*k/(alpha*alpha*costheta*costheta)); // formula (13.a)
- Double_t c2 = betaM*TMath::Power(k,1.5)*tantheta*lambda/(fParam->GapThick()*alpha*alpha); // formula (13.b)
- Double_t c3 = (1.+eTr*k*betaM*betaM)/((1+eTr)*alpha*alpha); // formula (13.c)
- Double_t c4 = TMath::Sqrt(k)*tantheta*(1-lambda)/(fParam->GapThick()*betaM); // formula (13.d)
- Double_t dtdT = c1 * (c2+c3*c4);
- Double_t trErr = fParam->RadThick()/(TMath::Sqrt(12.)*cost);
-
- return trErr*dtdT;
-}//SigGeom()
+ 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