#include "AliHMPIDTracker.h"
#include "AliHMPIDtrack.h"
#include "AliHMPIDCluster.h"
#include "AliHMPIDParam.h"
#include "AliHMPIDPid.h"
#include "AliHMPIDRecon.h"
#include "AliHMPIDRecoParamV1.h"
#include "AliHMPIDReconstructor.h"//Recon()
#include "AliHMPIDReconHTA.h"
#include <AliLog.h>
#include <AliESDEvent.h>
#include <AliESDtrack.h>
#include <AliTracker.h>
#include <AliRun.h>
#include <AliTrackPointArray.h>
#include <AliAlignObj.h>
#include <AliCDBManager.h>
#include <AliCDBEntry.h>
#include "TTreeStream.h"
ClassImp(AliHMPIDTracker)
AliHMPIDTracker::AliHMPIDTracker():
AliTracker(),
fClu(new TObjArray(AliHMPIDParam::kMaxCh+1)),
fDebugStreamer(0)
{
fClu->SetOwner(kTRUE);
for(int i=AliHMPIDParam::kMinCh;i<=AliHMPIDParam::kMaxCh;i++) fClu->AddAt(new TClonesArray("AliHMPIDCluster"),i);
fDebugStreamer = new TTreeSRedirector("HMPIDdebug.root");
}
AliHMPIDTracker::~AliHMPIDTracker(){
delete fClu;
if (fDebugStreamer) delete fDebugStreamer;
}
Bool_t AliHMPIDTracker::GetTrackPoint(Int_t idx, AliTrackPoint& point) const
{
if(idx<0) return kFALSE;
Int_t iCham=idx/1000000; Int_t iClu=idx%1000000;
iClu = iClu%1000;
point.SetVolumeID(AliGeomManager::LayerToVolUID(AliGeomManager::kHMPID,iCham));
TClonesArray *pArr=(TClonesArray*)(*fClu)[iCham];
if(!pArr) return kFALSE;
AliHMPIDCluster *pClu=(AliHMPIDCluster*)pArr->UncheckedAt(iClu);
if(!pClu) return kFALSE;
Float_t xyz[3];
pClu->GetGlobalXYZ(xyz);
Float_t cov[6];
pClu->GetGlobalCov(cov);
point.SetXYZ(xyz,cov);
return kTRUE;
}
Int_t AliHMPIDTracker::IntTrkCha(AliESDtrack *pTrk,Float_t &xPc,Float_t &yPc,Float_t &xRa,Float_t &yRa,Float_t &theta,Float_t &phi)
{
AliHMPIDtrack *hmpTrk = new AliHMPIDtrack(*pTrk);
for(Int_t i=AliHMPIDParam::kMinCh;i<=AliHMPIDParam::kMaxCh;i++){
Int_t chInt = IntTrkCha(i,hmpTrk,xPc,yPc,xRa,yRa,theta,phi);
if(chInt>=0) {delete hmpTrk;hmpTrk=0x0;return chInt;}
}
delete hmpTrk; hmpTrk=0x0;
return -1;
}
Int_t AliHMPIDTracker::IntTrkCha(Int_t ch,AliHMPIDtrack *pTrk,Float_t &xPc,Float_t &yPc,Float_t &xRa,Float_t &yRa,Float_t &theta,Float_t &phi)
{
AliHMPIDParam *pParam=AliHMPIDParam::Instance();
Double_t p1[3],n1[3];
pParam->Norm(ch,n1);
pParam->Point(ch,p1,AliHMPIDParam::kRad);
Double_t p2[3],n2[3];
pParam->Norm(ch,n2);
pParam->Point(ch,p2,AliHMPIDParam::kPc);
if(pTrk->Intersect(p1,n1)==kFALSE) return -1;
if(pTrk->Intersect(p2,n2)==kFALSE) return -1;
pParam->Mars2LorsVec(ch,n1,theta,phi);
pParam->Mars2Lors (ch,p1,xRa,yRa);
pParam->Mars2Lors (ch,p2,xPc,yPc);
if(AliHMPIDParam::IsInside(xPc,yPc,pParam->DistCut())==kTRUE) return ch;
return -1;
}
Int_t AliHMPIDTracker::LoadClusters(TTree *pCluTree)
{
for(int i=AliHMPIDParam::kMinCh;i<=AliHMPIDParam::kMaxCh;i++) pCluTree->SetBranchAddress(Form("HMPID%d",i),&((*fClu)[i]));
pCluTree->GetEntry(0);
return 0;
}
Int_t AliHMPIDTracker::PropagateBack(AliESDEvent *pEsd)
{
AliCDBEntry *pNmeanEnt =AliCDBManager::Instance()->Get("HMPID/Calib/Nmean");
AliCDBEntry *pQthreEnt =AliCDBManager::Instance()->Get("HMPID/Calib/Qthre");
if(!pNmeanEnt) AliError("No Nmean C6F14 ");
if(!pQthreEnt) AliError("No Qthre");
return Recon(pEsd,fClu,(TObjArray*)pNmeanEnt->GetObject(),(TObjArray*)pQthreEnt->GetObject());
}
Int_t AliHMPIDTracker::Recon(AliESDEvent *pEsd,TObjArray *pClus,TObjArray *pNmean, TObjArray *pQthre)
{
const Double_t kMaxSnp=0.9;
const Double_t kdRadiator=10;
AliHMPIDRecon recon;
AliHMPIDParam *pParam = AliHMPIDParam::Instance();
Float_t xPc,yPc,xRa,yRa,theta,phi;
Double_t cluLORS[2]={0};
Int_t nMipClusTot=0;
Double_t qthre = 0; Double_t nmean=0; Int_t hvsec=0;
Bool_t tsRight = kTRUE;
UInt_t tsmin = (UInt_t)((TF1*)pQthre->At(0))->GetXmin();
UInt_t tsmax = (UInt_t)((TF1*)pQthre->At(0))->GetXmax();
UInt_t ts = pEsd->GetTimeStamp();
if(ts<tsmin || ts>tsmax) {
AliWarning(Form(" in HMPID time stamp out of range!!! Please check!!! ts = %i",ts));
tsRight = kFALSE;
}
for(Int_t iTrk=0;iTrk<pEsd->GetNumberOfTracks();iTrk++){
Double_t dmin=999999,distCut=1,distParams[5]={1};
Bool_t isOkDcut=kFALSE;
Bool_t isOkQcut=kFALSE;
Bool_t isMatched=kFALSE;
AliHMPIDCluster *bestHmpCluster=0x0;
AliESDtrack *pTrk = pEsd->GetTrack(iTrk);
if(!pTrk->IsOn(AliESDtrack::kTPCout)) continue;
if(pTrk->IsOn(AliESDtrack::kTPCrefit)) continue;
AliESDfriendTrack *ftrack= (AliESDfriendTrack *)pTrk->GetFriendTrack();
if (!ftrack) continue;
if (!ftrack->GetTPCOut()) continue;
AliHMPIDtrack *hmpTrk = new AliHMPIDtrack(*pTrk);
AliHMPIDtrack *hmpTrkConstrained = 0;
hmpTrk->Set(ftrack->GetTPCOut()->GetX(), ftrack->GetTPCOut()->GetAlpha(),ftrack->GetTPCOut()->GetParameter(), ftrack->GetTPCOut()->GetCovariance());
pTrk->SetHMPIDtrk(0,0,0,0);
pTrk->SetHMPIDmip(0,0,0,0);
pTrk->SetHMPIDcluIdx(99,99999);
pTrk->SetHMPIDsignal(AliHMPIDRecon::kNotPerformed);
Int_t ipCh=IntTrkCha(pTrk,xPc,yPc,xRa,yRa,theta,phi);
if(ipCh<0) {delete hmpTrk;hmpTrk=0x0;continue;}
pTrk->SetHMPIDtrk(xPc,yPc,theta,phi);
pTrk->SetHMPIDcluIdx(ipCh,9999);
TClonesArray *pMipCluLst=(TClonesArray *)pClus->At(ipCh);
nMipClusTot = pMipCluLst->GetEntries();
if(nMipClusTot==0) {delete hmpTrk;hmpTrk=0x0;continue;}
Int_t index=-1;
for (Int_t iClu=0; iClu<nMipClusTot;iClu++) {
AliHMPIDCluster *pClu=(AliHMPIDCluster*)pMipCluLst->UncheckedAt(iClu);
if(tsRight){
if(pQthre->GetEntriesFast()==pParam->kMaxCh+1) {
qthre=((TF1*)pQthre->At(pClu->Ch()))->Eval(ts);
} else {
hvsec = pParam->InHVSector(pClu->Y());
if(hvsec>=0) qthre=((TF1*)pQthre->At(6*ipCh+hvsec))->Eval(ts);
}
} else qthre = pParam->QCut();
if(pClu->Q()<qthre) continue;
isOkQcut = kTRUE;
cluLORS[0]=pClu->X(); cluLORS[1]=pClu->Y();
Double_t dist = TMath::Sqrt((xPc-cluLORS[0])*(xPc-cluLORS[0])+(yPc-cluLORS[1])*(yPc-cluLORS[1]));
if(dist<dmin) {
dmin = dist;
index=iClu;
bestHmpCluster=pClu;
}
}
if(!bestHmpCluster) {delete hmpTrk;hmpTrk=0x0; delete hmpTrkConstrained;hmpTrkConstrained=0x0; continue;}
TVector3 vG = pParam->Lors2Mars(ipCh,bestHmpCluster->X(),bestHmpCluster->Y());
Double_t gx = vG[0];
Double_t gy = vG[1];
Double_t gz = vG[2];
Double_t alpha=TMath::ATan2(gy,gx);
Double_t radiusH=TMath::Sqrt(gy*gy+gx*gx);
if (!(hmpTrk->Rotate(alpha,kTRUE))) continue;
if(!AliTrackerBase::PropagateTrackToBxByBz(hmpTrk,radiusH,pTrk->GetMass(),1,kFALSE,kMaxSnp,-1)) {delete hmpTrk;hmpTrk=0x0; delete hmpTrkConstrained;hmpTrkConstrained=0x0; continue;}
AliExternalTrackParam trackC(*hmpTrk);
Double_t posC[2]={0,gz};
Double_t covC[3]={0.1*0.1, 0, 0.1*0.1};
trackC.Update(posC,covC);
hmpTrkConstrained = new AliHMPIDtrack(*pTrk);
hmpTrkConstrained->Set(trackC.GetX(), trackC.GetAlpha(),trackC.GetParameter(), trackC.GetCovariance());
if(!AliTrackerBase::PropagateTrackToBxByBz(hmpTrkConstrained,radiusH-kdRadiator,pTrk->GetMass(),1,kFALSE,kMaxSnp,1)) {delete hmpTrk;hmpTrk=0x0; delete hmpTrkConstrained;hmpTrkConstrained=0x0;continue;}
Float_t xPc0=0, yPc0=0;
IntTrkCha(ipCh, hmpTrk, xPc0,yPc0,xRa,yRa,theta,phi);
IntTrkCha(ipCh, hmpTrkConstrained, xPc,yPc,xRa,yRa,theta,phi);
Int_t cluSiz = bestHmpCluster->Size();
pTrk->SetHMPIDmip(bestHmpCluster->X(),bestHmpCluster->Y(),(Int_t)bestHmpCluster->Q(),0);
pTrk->SetHMPIDcluIdx(ipCh,index+1000*cluSiz);
pTrk->SetHMPIDtrk(xPc0,yPc0,theta,phi);
if (AliHMPIDReconstructor::StreamLevel()>0) {
AliExternalTrackParam * trackTPC=new AliExternalTrackParam(*(ftrack->GetTPCOut()));
AliExternalTrackParam * trackCurrent=new AliExternalTrackParam(*pTrk);
if(!trackTPC->Rotate(alpha)) continue;
if(!trackCurrent->Rotate(alpha)) continue;
Bool_t statusTPC= AliTracker::PropagateTrackToBxByBz(trackTPC,radiusH,pTrk->GetMass(),1,kFALSE,kMaxSnp,-1);
Bool_t statusCurrent=AliTracker::PropagateTrackToBxByBz(trackCurrent,radiusH,pTrk->GetMass(),1,kFALSE,kMaxSnp,-1);
AliExternalTrackParam * trackTPCNB=new AliExternalTrackParam(*(ftrack->GetTPCOut()));
if(!trackTPCNB->Rotate(alpha)) continue;
Bool_t statusTPCNB=kTRUE;
Double_t bfield[3]={0,0,0};
for (Double_t radius=trackTPCNB->GetX(); radius<radiusH; radius+=1){
Double_t xyz[3];
trackTPCNB->GetXYZ(xyz);
GetBxByBz(xyz,bfield);
statusTPCNB&=trackTPCNB->PropagateToBxByBz(radius,bfield);
}
statusTPCNB&=trackTPCNB->PropagateToBxByBz(radiusH,bfield);
Double_t tanAlpha=TMath::Tan(TMath::ASin(trackTPC->GetSnp()));
Double_t deltaC= trackTPC->GetC(AliTrackerBase::GetBz())-ftrack->GetTPCOut()->GetC(AliTrackerBase::GetBz());
AliExternalTrackParam * trackTPCConstrained= new AliExternalTrackParam(*trackTPC);
Double_t pos[2]={0,gz};
Double_t cov[3]={0.1*0.1, 0, 0.1*0.1};
Double_t chi2C = trackTPCConstrained->GetPredictedChi2(pos,cov);
trackTPCConstrained->Update(pos,cov);
(*fDebugStreamer)<<"track"<<
"rH="<<radiusH<<
"angle="<<tanAlpha<<
"dC="<<deltaC<<
"trackTPC.="<<trackTPC<<
"trackTPCNB.="<<trackTPCNB<<
"chi2C="<<chi2C<<
"trackTPCC.="<<trackTPCConstrained<<
"trackCurrent.="<<trackCurrent<<
"sTPC="<<statusTPC<<
"sCurrent="<<statusCurrent<<
"cl.="<<bestHmpCluster<<
"t.="<<pTrk<<
"ft.="<<ftrack<<
"hmpTrk.="<<hmpTrk<<
"hmpTrkC.="<<hmpTrkConstrained<<
"gx="<<gx<<
"gy="<<gy<<
"gz="<<gz<<
"\n";
}
if(!isOkQcut) {
pTrk->SetHMPIDsignal(pParam->kMipQdcCut);
delete hmpTrk;hmpTrk=0x0;
delete hmpTrkConstrained;hmpTrkConstrained=0x0;
continue;
}
if(AliHMPIDReconstructor::GetRecoParam())
{
if(AliHMPIDReconstructor::GetRecoParam()->IsFixedDistCut()==kTRUE)
{
distCut=AliHMPIDReconstructor::GetRecoParam()->GetHmpTrackMatchingDist();
}
else
{
for(Int_t ipar=0;ipar<5;ipar++) distParams[ipar]=AliHMPIDReconstructor::GetRecoParam()->GetHmpTrackMatchingDistParam(ipar);
distCut=distParams[0]+distParams[1]*TMath::Power(distParams[2]*pTrk->GetP(),distParams[3]);
}
}
else {distCut=pParam->DistCut();}
dmin = TMath::Sqrt((xPc0-bestHmpCluster->X())*(xPc0-bestHmpCluster->X())+(yPc0-bestHmpCluster->Y())*(yPc0-bestHmpCluster->Y()));
if(dmin < distCut) {
isOkDcut = kTRUE;
}
if(!isOkDcut) {
pTrk->SetHMPIDsignal(pParam->kMipDistCut);
}
if(isOkQcut*isOkDcut) isMatched = kTRUE;
if(!isMatched) {delete hmpTrk;hmpTrk=0x0;delete hmpTrkConstrained;hmpTrkConstrained=0x0;continue;}
Bool_t isOk = kTRUE;
if(!isOk) {delete hmpTrk;hmpTrk=0x0; delete hmpTrkConstrained;hmpTrkConstrained=0x0; continue;}
pTrk->SetOuterHmpParam(hmpTrkConstrained,AliESDtrack::kHMPIDout);
FillResiduals(hmpTrk,bestHmpCluster,kFALSE);
Int_t iRad = pParam->Radiator(yRa);
if(tsRight){
if(pNmean->GetEntries()==21) {
nmean=((TF1*)pNmean->At(3*ipCh))->Eval(ts);
} else {
if(iRad < 0) {
nmean = -1;
} else {
Double_t tLow = ((TF1*)pNmean->At(6*ipCh+2*iRad ))->Eval(ts);
Double_t tHigh = ((TF1*)pNmean->At(6*ipCh+2*iRad+1))->Eval(ts);
Double_t tExp = pParam->FindTemp(tLow,tHigh,yRa);
nmean = pParam->NIdxRad(AliHMPIDParam::Instance()->GetEPhotMean(),tExp);
}
if(nmean < 0){
pTrk->SetHMPIDsignal(AliHMPIDRecon::kNoRad);
pTrk->SetHMPIDcluIdx(ipCh,index+1000*cluSiz);
delete hmpTrk;hmpTrk=0x0;
delete hmpTrkConstrained;hmpTrkConstrained=0x0;
continue;
}
}
} else nmean = pParam->MeanIdxRad();
recon.SetImpPC(xPc0,yPc0);
recon.CkovAngle(pTrk,(TClonesArray *)pClus->At(ipCh),index,nmean,xRa,yRa);
Double_t thetaCkov = pTrk->GetHMPIDsignal();
if (AliHMPIDReconstructor::StreamLevel()>0) {
AliExternalTrackParam * trackTPC=new AliExternalTrackParam(*(ftrack->GetTPCOut()));
AliExternalTrackParam * trackCurrent=new AliExternalTrackParam(*pTrk);
if(!trackTPC->Rotate(alpha)) continue;
if(!trackCurrent->Rotate(alpha)) continue;
Bool_t statusTPC= AliTracker::PropagateTrackToBxByBz(trackTPC,radiusH,pTrk->GetMass(),1,kFALSE,kMaxSnp,-1);
Bool_t statusCurrent=AliTracker::PropagateTrackToBxByBz(trackCurrent,radiusH,pTrk->GetMass(),1,kFALSE,kMaxSnp,-1);
Double_t tanAlpha=TMath::Tan(TMath::ASin(trackTPC->GetSnp()));
Double_t deltaC= trackTPC->GetC(AliTrackerBase::GetBz())-ftrack->GetTPCOut()->GetC(AliTrackerBase::GetBz());
AliExternalTrackParam * trackTPCNB=new AliExternalTrackParam(*(ftrack->GetTPCOut()));
if(!trackTPCNB->Rotate(alpha)) continue;
Bool_t statusTPCNB=kTRUE;
Double_t bfield[3]={0,0,0};
for (Double_t radius=trackTPCNB->GetX(); radius<radiusH; radius+=1){
Double_t xyz[3];
trackTPCNB->GetXYZ(xyz);
GetBxByBz(xyz,bfield);
statusTPCNB&=trackTPCNB->PropagateToBxByBz(radius,bfield);
}
statusTPCNB&=trackTPCNB->PropagateToBxByBz(radiusH,bfield);
AliExternalTrackParam * trackTPCConstrained= new AliExternalTrackParam(*trackTPC);
Double_t pos[2]={0,gz};
Double_t cov[3]={0.1*0.1, 0, 0.1*0.1};
Double_t chi2C = trackTPCConstrained->GetPredictedChi2(pos,cov);
trackTPCConstrained->Update(pos,cov);
(*fDebugStreamer)<<"track2"<<
"rH="<<radiusH<<
"angle="<<tanAlpha<<
"dC="<<deltaC<<
"trackTPC.="<<trackTPC<<
"trackTPCNB.="<<trackTPCNB<<
"chi2C="<<chi2C<<
"trackTPCC.="<<trackTPCConstrained<<
"trackCurrent.="<<trackCurrent<<
"sTPC="<<statusTPC<<
"sCurrent="<<statusCurrent<<
"cl.="<<bestHmpCluster<<
"t.="<<pTrk<<
"ft.="<<ftrack<<
"hmpTrk.="<<hmpTrk<<
"hmpTrkC.="<<hmpTrkConstrained<<
"gx="<<gx<<
"gy="<<gy<<
"gz="<<gz<<
"thetaCkov="<<thetaCkov<<
"\n";
}
if(pTrk->GetHMPIDsignal()<0) {
delete hmpTrk;hmpTrk=0x0;
delete hmpTrkConstrained; hmpTrkConstrained=0x0;
continue;}
AliHMPIDPid pID;
Double_t prob[5];
pID.FindPid(pTrk,nmean,5,prob);
pTrk->SetHMPIDpid(prob);
delete hmpTrk; hmpTrk=0x0;
delete hmpTrkConstrained; hmpTrkConstrained=0x0;
}
return 0;
}
Int_t AliHMPIDTracker::ReconFromKin(AliESDEvent *pEsd,TObjArray *pClus,TObjArray *pNmean, TObjArray *pQthre)
{
AliHMPIDRecon recon;
AliHMPIDParam *pParam = AliHMPIDParam::Instance();
Float_t xPc,yPc,xRa,yRa,theta,phi;
Double_t cluLORS[2]={0};
Int_t nMipClusTot=0;
Double_t qthre = 0; Double_t nmean=0; Int_t hvsec=0;
Bool_t tsRight = kTRUE;
UInt_t tsmin = (UInt_t)((TF1*)pQthre->At(0))->GetXmin();
UInt_t tsmax = (UInt_t)((TF1*)pQthre->At(0))->GetXmax();
UInt_t ts = pEsd->GetTimeStamp();
if(ts<tsmin || ts>tsmax) {
AliWarning(Form(" in HMPID time stamp out of range!!! Please check!!! ts = %i",ts));
tsRight = kFALSE;
}
for(Int_t iTrk=0;iTrk<pEsd->GetNumberOfTracks();iTrk++){
Double_t dmin=999999,distCut=1,distParams[5]={1};
Bool_t isOkDcut=kFALSE;
Bool_t isOkQcut=kFALSE;
Bool_t isMatched=kFALSE;
AliHMPIDCluster *bestHmpCluster=0x0;
AliESDtrack *pTrk = pEsd->GetTrack(iTrk);
AliHMPIDtrack *hmpTrk = new AliHMPIDtrack(*pTrk);
pTrk->SetHMPIDtrk(0,0,0,0);
pTrk->SetHMPIDmip(0,0,0,0);
pTrk->SetHMPIDcluIdx(99,99999);
pTrk->SetHMPIDsignal(AliHMPIDRecon::kNotPerformed);
Int_t ipCh=IntTrkCha(pTrk,xPc,yPc,xRa,yRa,theta,phi);
if(ipCh<0) {delete hmpTrk;hmpTrk=0x0;continue;}
pTrk->SetHMPIDtrk(xPc,yPc,theta,phi);
pTrk->SetHMPIDcluIdx(ipCh,9999);
TClonesArray *pMipCluLst=(TClonesArray *)pClus->At(ipCh);
nMipClusTot = pMipCluLst->GetEntries();
if(nMipClusTot==0) {delete hmpTrk;hmpTrk=0x0;continue;}
Int_t index=-1;
for (Int_t iClu=0; iClu<nMipClusTot;iClu++) {
AliHMPIDCluster *pClu=(AliHMPIDCluster*)pMipCluLst->UncheckedAt(iClu);
if(tsRight){
if(pQthre->GetEntriesFast()==pParam->kMaxCh+1) {
qthre=((TF1*)pQthre->At(pClu->Ch()))->Eval(ts);
} else {
hvsec = pParam->InHVSector(pClu->Y());
if(hvsec>=0) qthre=((TF1*)pQthre->At(6*ipCh+hvsec))->Eval(ts);
}
} else qthre = pParam->QCut();
if(pClu->Q()<qthre) continue;
isOkQcut = kTRUE;
cluLORS[0]=pClu->X(); cluLORS[1]=pClu->Y();
Double_t dist = TMath::Sqrt((xPc-cluLORS[0])*(xPc-cluLORS[0])+(yPc-cluLORS[1])*(yPc-cluLORS[1]));
if(dist<dmin) {
dmin = dist;
index=iClu;
bestHmpCluster=pClu;
}
}
if(!isOkQcut) {
pTrk->SetHMPIDsignal(pParam->kMipQdcCut);
delete hmpTrk;hmpTrk=0x0; continue;
}
Double_t radius = (pParam->Lors2Mars(ipCh,pParam->SizeAllX()/2,pParam->SizeAllY()/2)).Mag();
if(!AliTracker::PropagateTrackToBxByBz(hmpTrk,radius,pTrk->GetMass(),1,kFALSE)) {delete hmpTrk;hmpTrk=0x0;continue;}
if(!hmpTrk->PropagateTo(bestHmpCluster)) {delete hmpTrk;hmpTrk=0x0;continue;}
Int_t cluSiz = bestHmpCluster->Size();
pTrk->SetHMPIDmip(bestHmpCluster->X(),bestHmpCluster->Y(),(Int_t)bestHmpCluster->Q(),0);
pTrk->SetHMPIDcluIdx(ipCh,index+1000*cluSiz);
if(AliHMPIDReconstructor::GetRecoParam())
{
if(AliHMPIDReconstructor::GetRecoParam()->IsFixedDistCut()==kTRUE)
{
distCut=AliHMPIDReconstructor::GetRecoParam()->GetHmpTrackMatchingDist();
}
else
{
for(Int_t ipar=0;ipar<5;ipar++) distParams[ipar]=AliHMPIDReconstructor::GetRecoParam()->GetHmpTrackMatchingDistParam(ipar);
distCut=distParams[0]+distParams[1]*TMath::Power(distParams[2]*pTrk->GetP(),distParams[3]);
}
}
else {distCut=pParam->DistCut();}
if(dmin < distCut) {
isOkDcut = kTRUE;
}
if(!isOkDcut) {
pTrk->SetHMPIDsignal(pParam->kMipDistCut);
}
if(isOkQcut*isOkDcut) isMatched = kTRUE;
if(!isMatched) {delete hmpTrk;hmpTrk=0x0;continue;}
Bool_t isOk = hmpTrk->Update(bestHmpCluster,0.1,0);
if(!isOk) {delete hmpTrk;hmpTrk=0x0;continue;}
pTrk->SetOuterHmpParam(hmpTrk,AliESDtrack::kHMPIDout);
FillResiduals(hmpTrk,bestHmpCluster,kFALSE);
Int_t iRad = pParam->Radiator(yRa);
if(tsRight){
if(pNmean->GetEntries()==21) {
nmean=((TF1*)pNmean->At(3*ipCh))->Eval(ts);
} else {
if(iRad < 0) {
nmean = -1;
} else {
Double_t tLow = ((TF1*)pNmean->At(6*ipCh+2*iRad ))->Eval(ts);
Double_t tHigh = ((TF1*)pNmean->At(6*ipCh+2*iRad+1))->Eval(ts);
Double_t tExp = pParam->FindTemp(tLow,tHigh,yRa);
nmean = pParam->NIdxRad(AliHMPIDParam::Instance()->GetEPhotMean(),tExp);
}
if(nmean < 0){
pTrk->SetHMPIDsignal(AliHMPIDRecon::kNoRad);
pTrk->SetHMPIDcluIdx(ipCh,index+1000*cluSiz);
delete hmpTrk;hmpTrk=0x0;
continue;
}
}
} else nmean = pParam->MeanIdxRad();
recon.SetImpPC(xPc,yPc);
recon.CkovAngle(pTrk,(TClonesArray *)pClus->At(ipCh),index,nmean,xRa,yRa);
if(pTrk->GetHMPIDsignal()<0) {delete hmpTrk;hmpTrk=0x0;continue;}
AliHMPIDPid pID;
Double_t prob[5];
pID.FindPid(pTrk,nmean,5,prob);
pTrk->SetHMPIDpid(prob);
delete hmpTrk;hmpTrk=0x0;
}
return 0;
}
Int_t AliHMPIDTracker::ReconHiddenTrk(AliESDEvent *pEsd,TObjArray *pClus,TObjArray *pNmean, TObjArray *pQthre)
{
AliHMPIDReconHTA reconHTA;
AliHMPIDParam *pParam = AliHMPIDParam::Instance();
Bool_t tsRight = kTRUE;
UInt_t tsmin = (UInt_t)((TF1*)pQthre->At(0))->GetXmin();
UInt_t tsmax = (UInt_t)((TF1*)pQthre->At(0))->GetXmax();
UInt_t ts = pEsd->GetTimeStamp();
if(ts<tsmin || ts>tsmax) {
AliWarning(Form(" in HMPID time stamp out of range!!! Please check!!! ts = %i",ts));
tsRight = kFALSE;
}
for(Int_t iTrk=0;iTrk<pEsd->GetNumberOfTracks();iTrk++){
AliESDtrack *pTrk = pEsd->GetTrack(iTrk);
Int_t ipCh;
ipCh = pTrk->GetHMPIDcluIdx();ipCh/=1000000;
if(ipCh<0) continue;
TClonesArray *pMipCluLst=(TClonesArray *)pClus->At(ipCh);
Int_t nMipClusTot = pMipCluLst->GetEntries();
Double_t qMip=-1;
Int_t chMip=-1;
Double_t xMip = 0;
Double_t yMip = 0;
Int_t indMip = 0;
Int_t cluMipSiz = 0;
for (Int_t iClu=0; iClu<nMipClusTot;iClu++) {
AliHMPIDCluster *pClu=(AliHMPIDCluster*)pMipCluLst->UncheckedAt(iClu);
Double_t qClus = pClu->Q();
if(qClus>qMip) {
qMip = qClus;
chMip = pClu->Ch();
xMip = pClu->X();
yMip = pClu->Y();
indMip = iClu;
cluMipSiz = pClu->Size();
}
}
if(chMip<0) return 1;
Int_t hvsec;
Double_t qthre=0;
if(tsRight){
if(pQthre->GetEntriesFast()==pParam->kMaxCh+1) {
qthre=((TF1*)pQthre->At(chMip))->Eval(ts);
} else {
hvsec = pParam->InHVSector(yMip);
if(hvsec>=0) qthre=((TF1*)pQthre->At(6*chMip+hvsec))->Eval(ts);
}
} else qthre = pParam->QCut();
if(qMip<qthre) {
pTrk->SetHMPIDmip(xMip,yMip,(Int_t)qMip,0);
pTrk->SetHMPIDcluIdx(chMip,indMip+1000*cluMipSiz);
pTrk->SetHMPIDsignal(pParam->kMipQdcCut);
return 1;
}
pTrk->SetHMPIDmip(xMip,yMip,(Int_t)qMip,0);
pTrk->SetHMPIDcluIdx(chMip,indMip+1000*cluMipSiz);
Double_t yRa = yMip;
Double_t nmean;
Int_t iRad = pParam->Radiator(yRa);
if(tsRight){
if(pNmean->GetEntries()==21) {
nmean=((TF1*)pNmean->At(3*chMip))->Eval(ts);
} else {
if(iRad < 0) {
nmean = -1;
} else {
Double_t tLow = ((TF1*)pNmean->At(6*chMip+2*iRad ))->Eval(ts);
Double_t tHigh = ((TF1*)pNmean->At(6*chMip+2*iRad+1))->Eval(ts);
Double_t tExp = pParam->FindTemp(tLow,tHigh,yRa);
nmean = pParam->NIdxRad(AliHMPIDParam::Instance()->GetEPhotMean(),tExp);
}
if(nmean < 0){
pTrk->SetHMPIDsignal(AliHMPIDRecon::kNoRad);
return 1;
}
}
} else nmean = pParam->MeanIdxRad();
if(!reconHTA.CkovHiddenTrk(pTrk,(TClonesArray *)pClus->At(ipCh),indMip,nmean)) {
AliHMPIDPid pID;
Double_t prob[5];
pID.FindPid(pTrk,nmean,5,prob);
pTrk->SetHMPIDpid(prob);
}
}
return 0;
}
void AliHMPIDTracker::FillClusterArray(TObjArray* array) const {
for(Int_t iCh=AliHMPIDParam::kMinCh;iCh<=AliHMPIDParam::kMaxCh;iCh++){
TClonesArray *pCluArr=(TClonesArray*)(*fClu)[iCh];
for (Int_t iClu=0; iClu<pCluArr->GetEntriesFast();iClu++){
AliHMPIDCluster *pClu=(AliHMPIDCluster*)pCluArr->UncheckedAt(iClu);
array->AddLast(pClu);
}
pCluArr->Delete();
}
return;
}