X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=HMPID%2FAliHMPIDRecon.cxx;h=b3aa7c550f1ffd6da3a84958950f249733fc3182;hb=30b5dcd4c5d84955fa64f219666e60bfee4a9fbb;hp=820e818dce401f2a9cee40371753cea5cd3ba706;hpb=43d3333b8b92069a40de8f98ca477729166e258f;p=u%2Fmrichter%2FAliRoot.git diff --git a/HMPID/AliHMPIDRecon.cxx b/HMPID/AliHMPIDRecon.cxx index 820e818dce4..b3aa7c550f1 100644 --- a/HMPID/AliHMPIDRecon.cxx +++ b/HMPID/AliHMPIDRecon.cxx @@ -27,12 +27,14 @@ #include //HoughResponse() #include //CkovAngle() #include //CkovAngle() +#include //CkovAngle() //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ AliHMPIDRecon::AliHMPIDRecon(): TTask("RichRec","RichPat"), fPhotCnt(-1), fPhotFlag(0x0), + fPhotClusIndex(0x0), fPhotCkov(0x0), fPhotPhi(0x0), fPhotWei(0x0), @@ -60,8 +62,9 @@ void AliHMPIDRecon::InitVars(Int_t n) //.. //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]; @@ -73,89 +76,84 @@ 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; - Int_t sizeClu = -1; + Float_t mipX=-1,mipY=-1; + Int_t chId=-1,mipQ=-1,sizeClu = -1; + fPhotCnt=0; + for (Int_t iClu=0; iCluGetEntriesFast();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();sizeClu=pClu->Size();} //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->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 - pTrk->SetHMPIDcluIdx(chId,mipId+1000*sizeClu); //set index of cluster return; } - - if(mipId==-1) { - pTrk->SetHMPIDcluIdx(chId,9999); //set index of cluster - pTrk->SetHMPIDsignal(kMipQdcCut); - return; - } //no clusters with QDC more the threshold at all - pTrk->SetHMPIDcluIdx(chId,mipId+1000*sizeClu); //set chamber, index of cluster + cluster size - if(dMin>fParam->DistCut()) {pTrk->SetHMPIDsignal(kMipDistCut); return;} //closest cluster with enough charge is still too far from intersection - + fMipPos.Set(mipX,mipY); - - + //PATTERN RECOGNITION STARTED: - Int_t iNrec=FlagPhot(HoughResponse()); //flag photons according to individual theta ckov with respect to most probable + 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(iNrec<1){ pTrk->SetHMPIDsignal(kNoPhotAccept); //no photon candidates are accepted return; } + 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(); }//CkovAngle() //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ @@ -359,7 +357,7 @@ Double_t AliHMPIDRecon::FindRingCkov(Int_t) 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() @@ -368,6 +366,7 @@ Int_t AliHMPIDRecon::FlagPhot(Double_t ckov) // 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 @@ -382,10 +381,25 @@ Int_t AliHMPIDRecon::FlagPhot(Double_t ckov) fPhotFlag[i] = 0; if(fPhotCkov[i] >= tmin && fPhotCkov[i] <= tmax) { fPhotFlag[i]=2; + PhotIndex[iInsideCnt]=fPhotClusIndex[i]; iInsideCnt++; } } + + for (Int_t iClu=0; iCluGetEntriesFast();iClu++){//clusters loop + AliHMPIDCluster *pClu=(AliHMPIDCluster*)pCluLst->UncheckedAt(iClu); //get pointer to current cluster + for(Int_t j=0; jAddCalibObject(pClus); + } + } + } + + delete [] PhotIndex; + return iInsideCnt; + }//FlagPhot() //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ TVector2 AliHMPIDRecon::TracePhot(Double_t ckovThe,Double_t ckovPhi)const @@ -480,3 +494,31 @@ Double_t AliHMPIDRecon::HoughResponse() 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;jIsInDead(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; +}