#include "AliHMPIDReconHTA.h"//class header
#include "AliHMPIDCluster.h" //CkovHiddenTrk()
#include "AliHMPIDRecon.h" //FunMinPhot()
+#include <TFile.h> //Database()
#include <TMinuit.h> //FitFree()
#include <TClonesArray.h> //CkovHiddenTrk()
#include <AliESDtrack.h> //CkovHiddenTrk()
-#include <TH2I.h> //InitDatabase()
+#include <TH2F.h> //InitDatabase()
#include <TGraph.h> //ShapeModel()
#include <TSpline.h> //ShapeModel()
+#include <TCanvas.h> //ShapeModel()
#include "TStopwatch.h" //
-TH2I* AliHMPIDReconHTA::fgDatabase = 0x0;
-
+Int_t AliHMPIDReconHTA::fgDB[500][150]={{75000*0}};
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-AliHMPIDReconHTA::AliHMPIDReconHTA():TTask("RichRec","RichPat")
+AliHMPIDReconHTA::AliHMPIDReconHTA():
+ TTask("RichRec","RichPat"),
+ fMipX(-999),
+ fMipY(-999),
+ fMipQ(-999),
+ fRadX(-999),
+ fRadY(-999),
+ fIdxMip(0),
+ fNClu(0),
+ fXClu(0),
+ fYClu(0),
+ fPhiPhot(0),
+ fThetaPhot(0),
+ fClCk(0),
+ fThTrkIn(-999),
+ fPhTrkIn(-999),
+ fThTrkFit(-999),
+ fPhTrkFit(-999),
+ fCkovFit(-999),
+ fNCluFit(0),
+ fCkovSig2(0),
+ fFitStatus(0),
+ fParam(AliHMPIDParam::Instance())
{
//..
//hidden algorithm
//..
- fMipX=fMipY=fThTrkFit=fPhTrkFit=fCkovFit=fMipQ=fRadX=fRadY=-999;
- fIdxMip=fNClu=0;
- fCkovSig2=0;
- fXClu = 0x0;
- fYClu = 0x0;
- fClCk = 0x0;
-
- fParam=AliHMPIDParam::Instance();
-
fParam->SetRefIdx(fParam->MeanIdxRad()); // initialization of ref index to a default one
- if(!fgDatabase) InitDatabase();
+ InitDatabase();
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
AliHMPIDReconHTA::~AliHMPIDReconHTA()
//..
fXClu = new Double_t[n];
fYClu = new Double_t[n];
+ fPhiPhot = new Double_t[n];
+ fThetaPhot = new Double_t[n];
fClCk = new Bool_t[n];
for(Int_t i=0;i<n;i++) fClCk[i] = kTRUE;
//
//..
if(fXClu) delete fXClu;
if(fYClu) delete fYClu;
+ if(fPhiPhot) delete fPhiPhot;
+ if(fThetaPhot) delete fThetaPhot;
if(fClCk) delete fClCk;
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-Bool_t AliHMPIDReconHTA::CkovHiddenTrk(AliESDtrack *pTrk,TClonesArray *pCluLst,Double_t nmean, Double_t qthre)
+Bool_t AliHMPIDReconHTA::CkovHiddenTrk(AliESDtrack *pTrk,TClonesArray *pCluLst,Int_t index, Double_t nmean)
{
// Pattern recognition method without any infos from tracking:HTA (Hidden Track Algorithm)...
// The method finds in the chmuber the cluster with the highest charge
// compatibile with a MIP, then the strategy is applied
// Arguments: pTrk - pointer to ESD track
// pCluLs - list of clusters for a given chamber
-// nmean - mean freon ref. index
+// pNmean - pointer to ref. index
+// pQthre - pointer to qthre
// Returns: - 0=ok,1=not fitted
AliHMPIDParam *pParam = AliHMPIDParam::Instance();
- pParam->SetRefIdx(nmean);
- if(!CluPreFilter(pCluLst)) {return kFALSE;}
-
- Float_t mipX=-1,mipY=-1;Int_t mipId=-1,mipQ=-1;
- Double_t qRef = 0;
+ if(!CluPreFilter(pCluLst)) return kFALSE;
+
Int_t nCh=0;
- for (Int_t iClu=0;iClu<pCluLst->GetEntriesFast();iClu++){ //clusters loop
+ Int_t sizeClu=0;
+
+ fNClu = pCluLst->GetEntriesFast();
+
+ for (Int_t iClu=0;iClu<fNClu;iClu++){ //clusters loop
AliHMPIDCluster *pClu=(AliHMPIDCluster*)pCluLst->UncheckedAt(iClu); //get pointer to current cluster
- nCh = pClu->Ch();
fXClu[iClu] = pClu->X();fYClu[iClu] = pClu->Y(); //store x,y for fitting procedure
fClCk[iClu] = kTRUE; //all cluster are accepted at this stage to be reconstructed
- if(pClu->Q()>qthre) fClCk[iClu] = kFALSE; // not a good photon in any case (multiple MIPs)
- if(pClu->Q()>qRef){ //searching the highest charge to select a MIP
- qRef = pClu->Q();
- mipId=iClu; mipX=pClu->X();mipY=pClu->Y();mipQ=(Int_t)pClu->Q();
- }
-// Printf(" n. %d x %f y %f Q %f",iClu,pClu->X(),pClu->Y(),pClu->Q());
+
+ if(iClu == index) {
+
+ fMipX = pClu->X();
+ fMipY = pClu->Y();
+ fMipQ = pClu->Q();
+ sizeClu = pClu->Size();
+ nCh = pClu->Ch();
+ fClCk[index] = kFALSE;
+ fIdxMip = index;
+ AliDebug(1,Form(" MIP n. %i x %f y %f Q %f",iClu,pClu->X(),pClu->Y(),pClu->Q()));
+ }
}//clusters loop
-
- fNClu = pCluLst->GetEntriesFast();
- if(qRef>qthre){ //charge compartible with MIP clusters
- fIdxMip = mipId;
- fClCk[mipId] = kFALSE;
- fMipX = mipX; fMipY=mipY; fMipQ = qRef;
-// Printf(" mipId %d x %f y %f Q %f",fIdxMip,fMipX,fMipY,fMipQ);
- if(!DoRecHiddenTrk()) {
- pTrk->SetHMPIDsignal(pParam->kNoPhotAccept);
- return kFALSE;
- } //Do track and ring reconstruction,if problems returns 1
- pTrk->SetHMPIDtrk(fRadX,fRadY,fThTrkFit,fPhTrkFit); //store track intersection info
- pTrk->SetHMPIDmip(fMipX,fMipY,(Int_t)fMipQ,fNClu); //store mip info
- pTrk->SetHMPIDcluIdx(nCh,fIdxMip); //set cham number and index of cluster
- pTrk->SetHMPIDsignal(fCkovFit); //find best Theta ckov for ring i.e. track
- pTrk->SetHMPIDchi2(fCkovSig2); //errors squared
-// Printf(" n clusters tot %i accepted %i",pCluLst->GetEntriesFast(),fNClu);
-// Printf("CkovHiddenTrk: thetaC %f th %f ph %f",fCkovFit,fThTrkFit,fPhTrkFit);
- return kTRUE;
- }
- return kFALSE;
-}//CkovHiddenTrk()
-//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-Bool_t AliHMPIDReconHTA::DoRecHiddenTrk()
-{
-// Pattern recognition method without any infos from tracking...
-// First a preclustering filter to avoid part of the noise
-// Then only ellipsed-rings are fitted (no possibility,
-// for the moment, to reconstruct very inclined tracks)
-// Finally a fitting with (th,ph) free, starting by very close values
-// previously evaluated.
-// Arguments: none
-// Returns: none
- Double_t thTrkRec,phiTrkRec,thetaCRec;
+ pParam->SetRefIdx(nmean);
- if(!FindShape(thTrkRec,phiTrkRec,thetaCRec)) return kFALSE;
-
- if(!FitFree(thTrkRec,phiTrkRec)) {return kFALSE;}
+ //
+ Float_t xra,yra,th,ph; pTrk->GetHMPIDtrk(xra,yra,th,ph);
+ AliDebug(1,Form(" simulated phiTRK %6.2f thetaTRK %6.2f",ph*TMath::RadToDeg(),th*TMath::RadToDeg()));
+ //
+
+ if(!DoRecHiddenTrk()) {
+ pTrk->SetHMPIDsignal(pParam->kNoPhotAccept);
+ return kFALSE;
+ } //Do track and ring reconstruction,if problems returns 1
+ AliDebug(1,Form(" fitted phi %6.2f ",fPhTrkFit*TMath::RadToDeg()));
+
+ pTrk->SetHMPIDtrk(fRadX,fRadY,fThTrkFit,fPhTrkFit); //store track intersection info
+ pTrk->SetHMPIDmip(fMipX,fMipY,(Int_t)fMipQ,NCluFit()); //store mip info + n. phots of the ring
+ pTrk->SetHMPIDcluIdx(nCh,fIdxMip+1000*sizeClu); //set cham number, index of cluster + cluster size
+ pTrk->SetHMPIDsignal(fCkovFit); //find best Theta ckov for ring i.e. track
+ pTrk->SetHMPIDchi2(fCkovSig2); //errors squared
+ AliDebug(1,Form(" n clusters tot %i fitted to ring %i",fNClu,NCluFit()));
+ for(Int_t i=0;i<fNClu;i++) {
+ AliDebug(1,Form(" n.%3i ThetaCer %8.3f PhiCer %8.3f check %i",i,fThetaPhot[i],fPhiPhot[i],fClCk[i]));
+ }
+ AliDebug(1,Form("CkovHiddenTrk: thetaC %f th %f ph %f",fCkovFit,fThTrkFit*TMath::RadToDeg(),fPhTrkFit*TMath::RadToDeg()));
return kTRUE;
-}//DoRecHiddenTrk()
-/*
+
+}//CkovHiddenTrk()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-Bool_t AliHMPIDReconHTA::DoRecHiddenTrk(TClonesArray *pCluLst)
+Bool_t AliHMPIDReconHTA::DoRecHiddenTrk()
{
// Pattern recognition method without any infos from tracking...
// First a preclustering filter to avoid part of the noise
// Returns: none
Double_t thTrkRec,phiTrkRec,thetaCRec;
- if(!FindShape(thTrkRec,phiTrkRec,thetaCRec)) return kFALSE;
+ if(!FindShape(thTrkRec,phiTrkRec,thetaCRec)) {
+ AliDebug(1,Form(" FindShape failed...!"));
+ return kFALSE;
+ }
+ AliDebug(1,Form(" FindShape accepted...!"));
-// Printf("thTrkRec %f phiTrkRec %f ThetaCRec %f",thTrkRec*TMath::RadToDeg(),phiTrkRec*TMath::RadToDeg(),thetaCRec*TMath::RadToDeg());
- Int_t nClTmp1 = pCluLst->GetEntriesFast()-1; //minus MIP...
- Int_t nClTmp2 = 0;
-
- while(nClTmp1 != nClTmp2){
- SetNClu(pCluLst->GetEntriesFast());
- if(!FitFree(thTrkRec,phiTrkRec)) {return kFALSE;}
- nClTmp2 = NClu();
- if(nClTmp2!=nClTmp1) {nClTmp1=nClTmp2;nClTmp2=0;}
+ if(!FitRing(thTrkRec,phiTrkRec)) {
+ AliDebug(1,Form(" FitRing failed...!"));
+ return kFALSE;
}
+ AliDebug(1,Form(" FitRing accepted...!"));
- fNClu = nClTmp2;
+ if(!UniformDistrib()) {
+ AliDebug(1,Form(" UniformDistrib failed...!"));
+ return kFALSE;
+ }
+ AliDebug(1,Form(" UniformDistrib accepted...!"));
+
return kTRUE;
}//DoRecHiddenTrk()
-*/
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Bool_t AliHMPIDReconHTA::CluPreFilter(TClonesArray *pCluLst)
{
Bool_t status;
-// Sort in phi angle...
+// Sort in phi angle...
for(Int_t i=0;i<fNClu;i++) {
+ if(!fClCk[i]) {
+ AliDebug(1,Form(" n.%3i xMIP %8.3f yMIP %8.3f check %i",i,fMipX,fMipY,fClCk[i]));
+ phiphot[i] = 999.;
+ dist[i] = 999.;
+ continue;
+ }
phiphot[i] = (TMath::ATan2(fMipY-fYClu[i],fMipX-fXClu[i])+TMath::Pi())*TMath::RadToDeg();
dist[i]=TMath::Sqrt((fMipX-fXClu[i])*(fMipX-fXClu[i])+(fMipY-fYClu[i])*(fMipY-fYClu[i]));
+ AliDebug(1,Form(" n.%3i phiphot %8.3f dist %8.3f check %i",i,phiphot[i],dist[i],fClCk[i]));
}
TMath::Sort(fNClu,phiphot,indphi,kFALSE);
for(Int_t i=0;i<fNClu;i++) {
if(!fClCk[indphi[i]]) continue; // Check if a good photon candidate or not
- if(TMath::Abs(dMean-dist[indphi[i]]) > 2*rms) {
+ if(TMath::Abs(dMean-dist[indphi[i]]) > 1.5*rms) {
fClCk[indphi[i]] = kFALSE;
continue;
}
}
+ AliDebug(1,"Purification of photons...");
//
// purify vectors for good photon candidates
//
Double_t *phiphotP = new Double_t[fNClu+1];
Double_t *distP = new Double_t[fNClu+1];
for(Int_t i=0;i<fNClu;i++) {
+ AliDebug(1,Form(" n. %3i phiphot %8.3f dist %8.3f check %i",i,phiphot[indphi[i]],dist[indphi[i]],fClCk[indphi[i]]));
if(!fClCk[indphi[i]]) continue; // Check if a good photon candidate or not
phiphotP[npeff] = phiphot[indphi[i]];
distP[npeff] = dist[indphi[i]];
delete [] indphi;
if(npeff<3) {
+ AliDebug(1,Form("FindShape failed: no enough photons = %i...",npeff));
delete [] phiphotP;
delete [] distP;
return kFALSE;
}
-// for(Int_t i=0;i<npeff;i++) {Printf(" n. %d phiphot %f dist %f",i,phiphotP[i],distP[i]);}
-
Double_t xA,xB;
- if(ShapeModel(npeff,phiphotP,distP,xA,xB,phiTrkRec)) {
+ status = kFALSE;
+
+ if(!ShapeModel(npeff,phiphotP,distP,xA,xB,phiTrkRec)) {AliDebug(1,Form("ShapeModel failed ")); return kFALSE;}
+
+// if(xA > 50 || xB > 15) {AliDebug(1,Form("xA and xB failed out of range")); return kFALSE;}
+
+ Int_t binxDB,binyDB;
+ Int_t compactDB=-1;
+
+ if(xA > xB) { //geometrically not possible, try to recover on a mean values...
+
+ FindBinDB(xA,xA,binxDB,binyDB);
+ if(binxDB<0 || binyDB<0) {AliDebug(1,Form("bin < 0 ! failed ")); return kFALSE;}
+ Int_t compactDB1 = CompactDB(binxDB,binyDB);
+ FindBinDB(xB,xB,binxDB,binyDB);
+ if(binxDB<0 || binyDB<0) {AliDebug(1,Form("bin < 0 ! failed ")); return kFALSE;}
+ Int_t compactDB2 = CompactDB(binxDB,binyDB);
+ Double_t thetaCRec1 = (Double_t)(compactDB1%1000);
+ Double_t thetaCRec2 = (Double_t)(compactDB2%1000);
+ Double_t thTrkRec1 = (Double_t)(compactDB1/1000);
+ Double_t thTrkRec2 = (Double_t)(compactDB2/1000);
+ thetaCRec = 0.5*(thetaCRec1+thetaCRec2);
+ thTrkRec = 0.5*( thTrkRec1+ thTrkRec2);
-// Printf("FindShape: phi start %f",phiTrkRec*TMath::RadToDeg());
+ } else {
+
+ FindBinDB(xA,xB,binxDB,binyDB);
+ if(binxDB<0 || binyDB<0) {AliDebug(1,Form("bin < 0 ! failed ")); return kFALSE;}
- Int_t bin = fgDatabase->FindBin(xA,xB);
- Int_t compact = (Int_t)fgDatabase->GetBinContent(bin);
- thetaCRec = (Double_t)(compact%1000);
- thTrkRec = (Double_t)(compact/1000);
+ compactDB = CompactDB(binxDB,binyDB);
- thTrkRec *= TMath::DegToRad();
- thetaCRec *= TMath::DegToRad();
+ if(compactDB<0) {AliDebug(1,Form("compact< 0! failed ")); return kFALSE;}
+ //
+ //
+ thetaCRec = (Double_t)(compactDB%1000);
+ thTrkRec = (Double_t)(compactDB/1000);
-// Printf("FindShape: xA %f xB %f compact %d thTrk %f thC %f",xA,xB,compact,thTrkRec*TMath::RadToDeg(),thetaCRec*TMath::RadToDeg());
-
- status = kTRUE;
-
- } else {
-
- status = kFALSE;
-
}
+ AliDebug(1,Form(" CompactDB %i thTrkRec %8.3f thetaCRec %8.3f ",compactDB,thTrkRec,thetaCRec));
+
+ phiTrkRec *= TMath::DegToRad();
+ thTrkRec *= TMath::DegToRad();
+ thetaCRec *= TMath::DegToRad();
+
+ status = kTRUE;
+
delete [] phiphotP;
delete [] distP;
// phiStart- estimate of the track phi
TGraph *phigr = new TGraph(np,phiphot,dist);
- TSpline3 *sphi = new TSpline3("sphi",phigr);
- if(!sphi) {Printf("Spline not created!Bye.");return kFALSE;}
-
- Int_t locMin = TMath::LocMin(np,dist);
- Int_t locMax = TMath::LocMax(np,dist);
-
- Double_t minX = phiphot[locMin];
-// Double_t minY = dist[locMin];
- Double_t maxX = phiphot[locMax];
-// Double_t maxY = dist[locMax];
-
- Int_t ip[3] = {-1,0,1};
- if(locMin==0 ) {ip[0]= 0;ip[1]= 1;ip[2]= 2;}
- if(locMin==np-1) {ip[0]=-2;ip[1]=-1;ip[2]= 0;}
-
- Double_t minXf = VertParab(phiphot[locMin+ip[0]],dist[locMin+ip[0]],
- phiphot[locMin+ip[1]],dist[locMin+ip[1]],
- phiphot[locMin+ip[2]],dist[locMin+ip[2]]);
- if(minXf< phiphot[locMin+ip[0]] || minXf > phiphot[locMin+ip[2]]) minXf = minX;
+ phiStart = FindSimmPhi();
- ip[0]=-1;ip[1]=0;ip[2]=1;
- if(locMax==0 ) {ip[0]= 0;ip[1]= 1;ip[2]= 2;}
- if(locMax==np-1) {ip[0]=-2;ip[1]=-1;ip[2]= 0;}
-
- Double_t maxXf = VertParab(phiphot[locMax+ip[0]],dist[locMax+ip[0]],
- phiphot[locMax+ip[1]],dist[locMax+ip[1]],
- phiphot[locMax+ip[2]],dist[locMax+ip[2]]);
- if(maxXf< phiphot[locMax+ip[0]] || maxXf > phiphot[locMax+ip[2]]) maxXf = maxX;
-
-// Printf(" phi at mindist %f and found %f",minX,minXf);
-// Printf(" phi at maxdist %f and found %f",maxX,maxXf);
-//
- if(TMath::Abs(maxXf-minXf)>30) {
- xA = sphi->Eval(minXf);
- if(xA < 0) xA = dist[sphi->FindX(xA)];
- xB = sphi->Eval(minXf-90);
- if(xB < 0) xB = dist[sphi->FindX(xB)];
- phiStart = minXf-180; //open ring or acceptance effect...so believe to min phi angle!
- } else {
- phiStart = 0.5*(maxXf-180+minXf);
- xA = sphi->Eval(phiStart+180);
- if(xA < 0) xA = dist[sphi->FindX(xA)];
- xB = sphi->Eval(phiStart+90);
- if(xB < 0) xB = dist[sphi->FindX(xB)];
- }
- //
-// Printf("ShapeModel: phiStart %f xA %f xB %f",phiStart,xA,xB);
+ Double_t phiStart1 = phiStart;
+ if(phiStart1 > 360) phiStart1 -= 360;
+ Double_t phiStart2 = phiStart+90;
+ if(phiStart2 > 360) phiStart2 -= 360;
+ xA = phigr->Eval(phiStart1);
+ xB = phigr->Eval(phiStart2);
+ //---
+ AliDebug(1,Form("phiStart %f phiStart1 %f phiStart2 %f ",phiStart,phiStart1,phiStart2));
+ AliDebug(1,Form("xA %f xB %f",xA,xB));
+
+ phiStart += 180;
+ if(phiStart>360) phiStart-=360;
- phiStart*=TMath::DegToRad();
return kTRUE;
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
return -0.5*b/a;
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+Bool_t AliHMPIDReconHTA::FitRing(Double_t thTrkRec,Double_t phiTrkRec)
+{
+ Double_t th = thTrkRec;
+ Double_t ph = phiTrkRec;
+
+ FitFree(th,ph);
+ while(FitStatus()) {
+ th = fThTrkFit;
+ ph = fPhTrkFit;
+ FitFree(th,ph);
+ }
+ return kTRUE;
+}
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Bool_t AliHMPIDReconHTA::FitFree(Double_t thTrkRec,Double_t phiTrkRec)
{
// Fit performed by minimizing RMS/sqrt(n) of the
if(thTrkRec==0) thTrkRec = 3.*TMath::DegToRad(); // not to start from the edge...
+ AliDebug(1,Form(" Minimization STARTED with phiTRK %6.2f thetaTRK %6.2f",phiTrkRec,thTrkRec));
+
pMinuit->mnparm(0," thTrk ",thTrkRec ,parStep=0.01,parLow=0,parHigh=TMath::PiOver4(),iErrFlg);
pMinuit->mnparm(1," phiTrk ",phiTrkRec,parStep=0.01,parLow=0,parHigh=TMath::TwoPi(),iErrFlg);
par[0] = th;par[1] = ph;
pMinuit->Eval(2,grad,f,par,3);
-// Printf("FitFree: theta %f phi %f",th,ph);
-
SetTrkFit(th,ph);
return kTRUE;
}
Double_t thetaCer,phiCer;
Int_t nClAcc = 0;
Int_t nClTot=pRecHTA->NClu();
-
+
for(Int_t i=0;i<nClTot;i++) {
+
+ if(pRecHTA->IdxMip() == i) {
+ pRecHTA->SetPhotAngles(i,999.,999.);
+ continue;
+ }
+
if(!(pRecHTA->ClCk(i))) continue;
- pRec.FindPhotCkov(pRecHTA->XClu(i),pRecHTA->YClu(i),thetaCer,phiCer);
+
+ Bool_t status = pRec.FindPhotCkov(pRecHTA->XClu(i),pRecHTA->YClu(i),thetaCer,phiCer);
+ if(!status) {
+ pRecHTA->SetPhotAngles(i,999.,999.);
+ continue;
+ }
+ pRecHTA->SetPhotAngles(i,thetaCer,phiCer);
meanCkov += thetaCer;
meanCkov2 += thetaCer*thetaCer;
nClAcc++;
}
+
if(nClAcc==0) {f=999;return;}
meanCkov /=(Double_t)nClAcc;
meanCkov2 /=(Double_t)nClAcc;
f = rms/TMath::Sqrt((Double_t)nClAcc);
if(iflag==3) {
-/*
- Printf("FunMinPhot before: photons candidates %i used %i",nClTot,nClAcc);
+ Int_t nClAccStep1 = nClAcc;
nClAcc = 0;
Double_t meanCkov1=0;
- Double_t meanCkov2=0;
+ Double_t meanCkov3=0;
for(Int_t i=0;i<nClTot;i++) {
- if(!(pRec->ClCk(i))) continue;
- pRec->FindPhotCkov(pRec->XClu(i),pRec->YClu(i),thetaCer,phiCer);
+
+ if(!(pRecHTA->ClCk(i))) continue;
+ thetaCer = pRecHTA->PhotTheta(i);
if(TMath::Abs(thetaCer-meanCkov)<2*rms) {
meanCkov1 += thetaCer;
- meanCkov2 += thetaCer*thetaCer;
+ meanCkov3 += thetaCer*thetaCer;
nClAcc++;
- } else pRec->SetClCk(i,kFALSE);
+ } else pRecHTA->SetClCk(i,kFALSE);
+ }
+
+ if(nClAcc<3) {
+ pRecHTA->SetFitStatus(kFALSE);
+ pRecHTA->SetCkovFit(meanCkov);
+ pRecHTA->SetCkovSig2(rms*rms);
+ pRecHTA->SetNCluFit(nClAccStep1);
+ return;
}
+
meanCkov1/=nClAcc;
- Double_t rms2 = (meanCkov2 - meanCkov*meanCkov*nClAcc)/nClAcc;
- Printf("FunMinPhot after: photons candidates %i used %i thetaCer %f",nClTot,nClAcc,meanCkov1);
- pRec->SetCkovFit(meanCkov1);
- pRec->SetCkovSig2(rms2);
- pRec->SetNClu(nClAcc);
-*/
-// Printf("FunMinPhot: reconstructed theta Cerenkov %f with %d photons",meanCkov,nClAcc);
- pRecHTA->SetCkovFit(meanCkov);
- pRecHTA->SetCkovSig2(rms*rms);
- pRecHTA->SetNClu(nClAcc);
+ Double_t rms2 = (meanCkov3 - meanCkov*meanCkov*nClAcc)/nClAcc;
+
+ // get a logger instance
+ // what for??
+ //AliLog::GetRootLogger();
+
+ if(nClAcc!=nClAccStep1) pRecHTA->SetFitStatus(kTRUE); else pRecHTA->SetFitStatus(kFALSE);
+
+ pRecHTA->SetCkovFit(meanCkov1);
+ pRecHTA->SetCkovSig2(rms2);
+ pRecHTA->SetNCluFit(nClAcc);
}
}//FunMinPhot()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Construction a database of ring shapes on fly
// Arguments: none
// Returns : none
-// N.B. fgDatabase points to a TH2I with x-min dist from MIP
-// y-dist from the ring of the MIP perpendicular to major axis
+// N.B. fgDB is the distance with x-min from MIP
+// y-dist from the ring of the MIP perpendicular to major axis
// The content is the packed info of track theta and thetaC in degrees
// thetaC+1000*thTrk
//
+// TFile *pout = new TFile("./database.root","recreate");
+
+ static Bool_t isDone = kFALSE;
+
TStopwatch timer;
- timer.Start();
+ if(isDone) {
+ return;
+ }
+
+ if(!isDone) {
+ timer.Start();
+ }
+
AliInfo(Form("database HTA is being built.Please, wait..."));
//
- Double_t x[3],y[3];
+ Double_t x[3]={0,0,0},y[3];
AliHMPIDRecon rec;
Int_t nstepx = 1000;
Int_t nstepy = 1000;
- TH2I *deconv = new TH2I("deconv","database;d1;d2;thC+1000*thTrk",500,0,50,150,0,15);
//
Double_t xrad = 0;
Double_t yrad = 0;
//
//mip position
//
- Double_t sizeCh = fParam->RadThick()+fParam->WinThick()+fParam->GapThick();
+ Double_t sizeCh = 0.5*fParam->RadThick()+fParam->WinThick()+fParam->GapThick();
Double_t xmip = xrad + sizeCh*TMath::Tan(thTrk)*TMath::Cos(phTrk);
Double_t ymip = yrad + sizeCh*TMath::Tan(thTrk)*TMath::Sin(phTrk);
pos=rec.TracePhot(thetaC,0);
if(pos.X()==-999) {
- dist1 = 0; //open ring...anly the distance btw mip and point at 180 will be considered
+ dist1 = 0; //open ring...only the distance btw mip and point at 180 will be considered
} else {
x[0] = pos.X(); y[0] = pos.Y();
dist1 = TMath::Sqrt((x[0]-xmip)*(x[0]-xmip)+(y[0]-ymip)*(y[0]-ymip));
rec.SetTrack(xrad,yrad,thTrk,phTrk);
pos=rec.TracePhot(thetaC,TMath::Pi());
- if(pos.X()==-999) {Printf("it should not happens!Bye");return;}
+ if(pos.X()==-999) {AliDebug(1,Form("it should not happens!Bye"));return;}
x[1] = pos.X(); y[1] = pos.Y();
if((x[1]-xmip)*(x[0]-xmip)>0) continue; // to avoid circles out mips (for very low ThetaC)
dist2 = TMath::Sqrt((x[1]-xmip)*(x[1]-xmip)+(y[1]-ymip)*(y[1]-ymip));
Double_t distB = TMath::Sqrt((x[2]-xmip)*(x[2]-xmip)+(y[2]-ymip)*(y[2]-ymip));
// compact the infos...
Int_t compact = (Int_t)(thetaC*TMath::RadToDeg())+1000*(Int_t)(thTrk*TMath::RadToDeg());
- Int_t bin = deconv->FindBin(distA,distB);
- if(deconv->GetBinContent(bin)==0) deconv->Fill(distA,distB,compact);
+ Int_t binxDB,binyDB;
+ FindBinDB(distA,distB,binxDB,binyDB);
+ if(fgDB[binxDB][binyDB]==0) fgDB[binxDB][binyDB] = compact;
}
}
- FillZeroChan(deconv);
- fgDatabase = deconv;
+ FillZeroChan();
- timer.Stop();
- Double_t nSecs = timer.CpuTime();
- AliInfo(Form("database HTA successfully open in %3.1f sec.(CPU). Reconstruction is started.",nSecs));
-}
+ if(!isDone) {
+ timer.Stop();
+ Double_t nSecs = timer.CpuTime();
+ AliInfo(Form("database HTA successfully open in %3.1f sec.(CPU). Reconstruction is started.",nSecs));
+ isDone = kTRUE;
+ }
+
+// pout->Write();
+// pout->Close();
+
+}//InitDatabase()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-void AliHMPIDReconHTA::FillZeroChan(TH2I *deconv)const
+void AliHMPIDReconHTA::FillZeroChan()const
{
//If fills eventually channel without entries
//inthe histo "database" jyst interpolating the neighboring cells
// Arguments: histogram pointer of the database
// Returns: none
//
- Int_t nbinx = deconv->GetNbinsX();
- Int_t nbiny = deconv->GetNbinsY();
- for(Int_t i = 0;i<nbinx;i++) {
- for(Int_t j = 0;j<nbiny;j++) {
- if(deconv->GetBinContent(i,j) == 0) {
+ const Int_t nxDB = 500;
+ const Int_t nyDB = 150;
+
+ for(Int_t i = 0;i<nxDB;i++) {
+ for(Int_t j = 0;j<nyDB;j++) {
+ if(fgDB[i][j] == 0) {
+ fgDB[i][j] = -1;
Int_t nXmin = i-1; Int_t nXmax=i+1;
Int_t nYmin = j-1; Int_t nYmax=j+1;
Int_t nc = 0;
Double_t meanC =0;
Double_t meanTrk =0;
for(Int_t ix=nXmin;ix<=nXmax;ix++) {
- if(ix<0||ix>nbinx) continue;
+ if(ix<0||ix>=nxDB) continue;
for(Int_t iy=nYmin;iy<=nYmax;iy++) {
- if(iy<0||iy>nbiny) continue;
- meanC += (Int_t)deconv->GetBinContent(ix,iy)%1000;
- meanTrk+= (Int_t)deconv->GetBinContent(ix,iy)/1000;
+ if(iy<0||iy>=nyDB) continue;
+ meanC += (Int_t)(fgDB[ix][iy]%1000);
+ meanTrk+= (Int_t)(fgDB[ix][iy]/1000);
nc++;
}
meanC/=nc; meanTrk/=nc;
Int_t compact = (Int_t)meanC+1000*(Int_t)meanTrk;
- if(compact>0)deconv->SetCellContent(i,j,compact);
+ if(compact>0)fgDB[i][j] = compact;
}
}
}
}
+}//FillZeroChan()
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+Int_t AliHMPIDReconHTA::r2(Double_t *coef, Double_t &x1, Double_t &x2)
+{
+ //2nd deg. equation
+ //solution
+ // Arguments: coef 2 1 0: ax^2+bx+c=0
+ // Returns: n. of solutions
+ // x1= 1st sol
+ // x2= 2nd sol
+ Double_t a,b,c;
+ a = coef[2];
+ b = coef[1];
+ c = coef[0];
+ Double_t delta = b*b-4*a*c;
+ if(delta<0) {return 0;}
+ if(delta==0) {
+ x1=x2=-b/(2*a);
+ return 1;
+ }
+ // delta>0
+ x1 = (-b+TMath::Sqrt(delta))/(2*a);
+ x2 = (-b-TMath::Sqrt(delta))/(2*a);
+ return 2;
+}//r2()
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+Double_t AliHMPIDReconHTA::FindSimmPhi()
+{
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+// Reconstruction of phiTRK angle with two methods (in switching)
+//
+// - least square method (for closed rings)
+// - by minimum distance mip-photon (for open rings)
+
+ Float_t coeff1ord=0; Float_t coeff2ord=0; Float_t coeff0ord=0;
+ Float_t xrotsumm =0; Float_t yrotsumm =0; Float_t xx =0;
+ Float_t yy =0; Float_t xy =0;
+ Double_t xmin=0;
+ Double_t ymin=0;
+
+ Int_t np=0;
+
+ Double_t distMin = 999.;
+
+ for(Int_t i=0;i<fNClu;i++) {
+ if(!fClCk[i]) continue;
+ np++;
+ xrotsumm+=fXClu[i]; // summ xi
+ yrotsumm+=fYClu[i]; // summ yi
+ xx+=fXClu[i]*fXClu[i]; // summ xixi
+ yy+=fYClu[i]*fYClu[i]; // summ yiyi
+ xy+=fXClu[i]*fYClu[i]; // summ yixi
+ Double_t dist2= (fXClu[i]-fMipX)*(fXClu[i]-fMipX)+(fYClu[i]-fMipY)*(fYClu[i]-fMipY);
+ if(dist2<distMin) {
+ distMin = dist2;
+ xmin = fXClu[i];
+ ymin = fYClu[i];
+ }
+ }
+
+ Double_t AngM = TMath::ATan2(ymin-fMipY,(xmin-fMipX))*TMath::RadToDeg();
+ if (AngM<0) AngM+=360;
+
+ AliDebug(1,Form(" Simple angle prediction with MIN phi = %f",AngM));
+
+ //_____calc. met min quadr using effective distance _________________________________________________
+
+ coeff2ord = xy-xrotsumm*yrotsumm/np;
+ coeff1ord = yrotsumm*yrotsumm/np - xrotsumm*xrotsumm/np - yy + xx;
+ coeff0ord = -coeff2ord;
+
+ AliDebug(1,Form(" a = %f b = %f c = %f",coeff2ord,coeff1ord,coeff0ord));
+
+ Double_t m1=0, m2=0; Double_t n1=0, n2=0;
+ // c // b // a
+ Double_t coeff[3]={coeff0ord,coeff1ord,coeff2ord};
+
+ r2(coeff,m1,m2);
+
+ n1=(yrotsumm-m1*xrotsumm)/np;
+ n2=(yrotsumm-m2*xrotsumm)/np;
+ // 2 solutions.................
+
+ // negative angles solved...
+
+ Double_t d1 =(1/(m1*m1+1))*(yy+m1*m1*xx+np*n1*n1-2*m1*xy-2*n1*yrotsumm+2*m1*n1*xrotsumm);
+ Double_t d2 =(1/(m2*m2+1))*(yy+m2*m2*xx+np*n2*n2-2*m2*xy-2*n2*yrotsumm+2*m2*n2*xrotsumm);
+
+ AliDebug(1,Form(" predicted distance d1 %f for angle %6.2f d2 %f for angle %6.2f",d1,TMath::ATan(m1)*TMath::RadToDeg(),
+ d2,TMath::ATan(m2)*TMath::RadToDeg()));
+ Double_t mMin;
+ if(d1 > d2) mMin = m2; else mMin = m1;
+
+ Double_t phiTrk=0;
+ //
+ if(ymin < fMipY && xmin > fMipX) {phiTrk = TMath::ATan(mMin)*TMath::RadToDeg()+180;}
+ if(ymin > fMipY && xmin < fMipX) {phiTrk = TMath::ATan(mMin)*TMath::RadToDeg()+360;}
+ if(ymin > fMipY && xmin > fMipX) {phiTrk = TMath::ATan(mMin)*TMath::RadToDeg()+180;}
+ if(ymin < fMipY && xmin < fMipX) {phiTrk = TMath::ATan(mMin)*TMath::RadToDeg();}
+ if(ymin == fMipY && xmin > fMipX) {phiTrk = TMath::ATan(mMin)*TMath::RadToDeg()+180;}
+ if(ymin == fMipY && xmin < fMipX) {phiTrk = TMath::ATan(mMin)*TMath::RadToDeg();}
+ if(ymin < fMipY && xmin == fMipX) {phiTrk = 90;}
+ if(ymin > fMipY && xmin == fMipX) {phiTrk = 270;}
+
+ // ------------------------- choose the best-----------------------
+
+
+ if( AngM-40 <= phiTrk && AngM+40 >= phiTrk) return phiTrk; else return AngM;
+}
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+void AliHMPIDReconHTA::FindBinDB(Double_t x,Double_t y,Int_t &binX,Int_t &binY)
+{
+ const Int_t nxDB = 500;
+ const Int_t nyDB = 150;
+ const Double_t xlowDB = 0;
+ const Double_t xhigDB = 50;
+ const Double_t ylowDB = 0;
+ const Double_t yhigDB = 30;
+
+ binX = -1;
+ binY = -1;
+ if(x<xlowDB && x>xhigDB &&
+ y<ylowDB && y>yhigDB) return;
+ binX = Int_t((x-xlowDB)/(xhigDB-xlowDB)*nxDB);
+ binY = Int_t((y-ylowDB)/(yhigDB-ylowDB)*nyDB);
+ if(binX>=nxDB || binY>=nyDB) {
+ binX = -1;
+ binY = -1;
+ }
+
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+Bool_t AliHMPIDReconHTA::UniformDistrib()
+{
+ AliHMPIDParam *pParam=AliHMPIDParam::Instance();
+ AliHMPIDRecon pRec;
+
+ Double_t sizeCh = 0.5*pParam->RadThick()+pParam->WinThick()+pParam->GapThick();
+ Double_t xrad = MipX() - sizeCh*TMath::Tan(fThTrkFit)*TMath::Cos(fPhTrkFit);
+ Double_t yrad = MipY() - sizeCh*TMath::Tan(fThTrkFit)*TMath::Sin(fPhTrkFit);
+
+ pRec.SetTrack(xrad,yrad,fThTrkFit,fPhTrkFit);
+
+ Int_t npeff=0;
+ Int_t nPhotPerBin = 4;
+ for(Int_t i=0;i<fNClu;i++) {
+ if(!ClCk(i)) continue;
+ npeff++;
+ }
+
+ Int_t nPhiBins = npeff/nPhotPerBin+1;
+ if(nPhiBins<=1) nPhiBins = 2;
+
+ Double_t *iPhiBin;
+ iPhiBin = new Double_t[nPhiBins];
+
+ for(Int_t i=0;i<nPhiBins;i++) iPhiBin[i] =0;
+
+ for(Int_t i=0;i<fNClu;i++) {
+ if(!ClCk(i)) continue;
+ Double_t phiCer = PhotPhi(i);
+
+ Double_t PhiProva = phiCer*TMath::RadToDeg();
+ if(PhiProva<0) PhiProva+= 360;
+ Int_t index = (Int_t)((Float_t)nPhiBins*PhiProva/360.);
+ iPhiBin[index]++;
+ }
+
+ Double_t chi2 = 0;
+ for(Int_t i=0;i<nPhiBins;i++) {
+ Double_t theo = (Double_t)npeff/(Double_t)nPhiBins;
+ chi2+= (iPhiBin[i] - theo)*(iPhiBin[i] - theo)/theo;
+ }
+
+ delete iPhiBin;
+
+ Double_t prob = TMath::Prob(chi2, nPhiBins-1);
+ AliDebug(1,Form(" Probability for uniform distrib: %6f.3 %s",prob,(prob<0.05) ? "rejected" : "accepted"));
+ if(prob<0.05) return kFALSE;
+ return kTRUE;
+
+ }
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++