#include "AliHMPIDCluster.h"
#include <TVirtualFitter.h> //Solve()
#include <TMinuit.h> //Solve()
#include <TClonesArray.h> //Solve()
#include <TMarker.h> //Draw()
#include "AliLog.h"
Bool_t AliHMPIDCluster::fgDoCorrSin=kTRUE;
ClassImp(AliHMPIDCluster)
void AliHMPIDCluster::SetClusterParams(Double_t xL,Double_t yL,Int_t iCh )
{
fParam = AliHMPIDParam::Instance();
if(!fParam->GetInstType())
{
new(this) AliCluster3D(); return;
}
UShort_t volId=AliGeomManager::LayerToVolUID(AliGeomManager::kHMPID,iCh);
const TGeoHMatrix *t2l= AliGeomManager::GetTracking2LocalMatrix(volId);
fParam = AliHMPIDParam::Instance();
xL -= 0.5*fParam->SizeAllX();
yL -= 0.5*fParam->SizeAllY();
Double_t posL[3]={xL, yL, 0.};
Double_t posT[3];
t2l->MasterToLocal(posL,posT);
Double_t covL[9] = {
0.8*0.8/12., 0., 0.0,
0., 0.84*0.84/12., 0.0,
0., 0., 0.1,
};
TGeoHMatrix m;
m.SetRotation(covL);
m.Multiply(t2l);
m.MultiplyLeft(&t2l->Inverse());
Double_t *covT = m.GetRotationMatrix();
new(this) AliCluster3D(volId,
posT[0],posT[1],posT[2],
covT[0],covT[1],covT[2],
covT[4],covT[5],
covT[8],
0x0);
}
AliHMPIDCluster::~AliHMPIDCluster()
{
if(fDigs) delete fDigs; fDigs=0;
}
void AliHMPIDCluster::CoG()
{
Int_t minPadX=999,minPadY=999,maxPadX=-1,maxPadY=-1;
if(fDigs==0) return;
fXX=fYY=fQRaw=0;
fCh = -1;
Int_t maxQpad=-1,maxQ=-1;
AliHMPIDDigit *pDig=0x0;
for(Int_t iDig=0;iDig<fDigs->GetEntriesFast();iDig++){
pDig=(AliHMPIDDigit*)fDigs->At(iDig);
if(!pDig) continue;
if(pDig->PadPcX() > maxPadX) maxPadX = pDig->PadPcX();
if(pDig->PadPcY() > maxPadY) maxPadY = pDig->PadPcY();
if(pDig->PadPcX() < minPadX) minPadX = pDig->PadPcX();
if(pDig->PadPcY() < minPadY) minPadY = pDig->PadPcY();
Float_t q=pDig->Q();
fXX += pDig->LorsX()*q;fYY +=pDig->LorsY()*q;
fQRaw+=q;
if(q>maxQ) {maxQpad = pDig->Pad();maxQ=(Int_t)q;}
fCh=pDig->Ch();
}
fBox=(maxPadX-minPadX+1)*100+maxPadY-minPadY+1;
if ( fQRaw != 0 ) {fXX/=fQRaw;fYY/=fQRaw;}
if(fDigs->GetEntriesFast()>1&&fgDoCorrSin)CorrSin();
fQ = fQRaw;
fMaxQpad = maxQpad; fMaxQ=maxQ;
fChi2=0;
fNlocMax=0;
fSt=kCoG;
SetClusterParams(fXX,fYY,fCh);
}
void AliHMPIDCluster::CorrSin()
{
Int_t pc,px,py;
fParam->Lors2Pad(fXX,fYY,pc,px,py);
Float_t x=fXX-fParam->LorsX(pc,px);
fXX+=3.31267e-2*TMath::Sin(2*TMath::Pi()/0.8*x)-2.66575e-3*TMath::Sin(4*TMath::Pi()/0.8*x)+2.80553e-3*TMath::Sin(6*TMath::Pi()/0.8*x)+0.0070;
}
void AliHMPIDCluster::Draw(Option_t*)
{
TMarker *pMark=new TMarker(X(),Y(),5); pMark->SetUniqueID(fSt);pMark->SetMarkerColor(kBlue); pMark->Draw();
}
void AliHMPIDCluster::FitFunc(Int_t &iNpars, Double_t* deriv, Double_t &chi2, Double_t *par, Int_t iflag)
{
AliHMPIDCluster *pClu=(AliHMPIDCluster*)TVirtualFitter::GetFitter()->GetObjectFit();
Int_t nPads = pClu->Size();
chi2 = 0;
Int_t iNshape = iNpars/3;
for(Int_t i=0;i<nPads;i++){
Double_t dQpadMath = 0;
for(Int_t j=0;j<iNshape;j++){
Double_t fracMathi = pClu->Dig(i)->IntMathieson(par[3*j],par[3*j+1]);
dQpadMath+=par[3*j+2]*fracMathi;
}
if(dQpadMath>0 && pClu->Dig(i)->Q()>0) {
chi2 +=TMath::Power((pClu->Dig(i)->Q()-dQpadMath),2)/pClu->Dig(i)->Q();
}
}
if(iflag==2) {
Double_t **derivPart;
derivPart = new Double_t*[iNpars];
for(Int_t j=0;j<iNpars;j++){
deriv[j] = 0;
derivPart[j] = new Double_t[nPads];
for(Int_t i=0;i<nPads;i++){
derivPart[j][i] = 0;
}
}
for(Int_t i=0;i<nPads;i++){
for(Int_t j=0;j<iNshape;j++){
Double_t fracMathi = pClu->Dig(i)->IntMathieson(par[3*j],par[3*j+1]);
derivPart[3*j ][i] += par[3*j+2]*(pClu->Dig(i)->MathiesonX(par[3*j]-pClu->Dig(i)->LorsX()-0.5*AliHMPIDParam::SizePadX())-
pClu->Dig(i)->MathiesonX(par[3*j]-pClu->Dig(i)->LorsX()+0.5*AliHMPIDParam::SizePadX()))*
pClu->Dig(i)->IntPartMathiY(par[3*j+1]);
derivPart[3*j+1][i] += par[3*j+2]*(pClu->Dig(i)->MathiesonY(par[3*j+1]-pClu->Dig(i)->LorsY()-0.5*AliHMPIDParam::SizePadY())-
pClu->Dig(i)->MathiesonY(par[3*j+1]-pClu->Dig(i)->LorsY()+0.5*AliHMPIDParam::SizePadY()))*
pClu->Dig(i)->IntPartMathiX(par[3*j]);
derivPart[3*j+2][i] += fracMathi;
}
}
for(Int_t i=0;i<nPads;i++){
Double_t dQpadMath = 0;
for(Int_t j=0;j<iNshape;j++){
Double_t fracMathi = pClu->Dig(i)->IntMathieson(par[3*j],par[3*j+1]);
dQpadMath+=par[3*j+2]*fracMathi;
if(dQpadMath>0 && pClu->Dig(i)->Q()>0) {
deriv[3*j] += 2/pClu->Dig(i)->Q()*(pClu->Dig(i)->Q()-dQpadMath)*derivPart[3*j ][i];
deriv[3*j+1] += 2/pClu->Dig(i)->Q()*(pClu->Dig(i)->Q()-dQpadMath)*derivPart[3*j+1][i];
deriv[3*j+2] += 2/pClu->Dig(i)->Q()*(pClu->Dig(i)->Q()-dQpadMath)*derivPart[3*j+2][i];
}
}
}
for(Int_t i=0;i<iNpars;i++) delete [] derivPart[i]; delete [] derivPart;
}
}
void AliHMPIDCluster::Print(Option_t* opt)const
{
const char *status=0;
switch(fSt){
case kFrm : status="formed " ;break;
case kUnf : status="unfolded (fit)" ;break;
case kCoG : status="coged " ;break;
case kLo1 : status="locmax 1 (fit)" ;break;
case kMax : status="exceeded (cog)" ;break;
case kNot : status="not done (cog)" ;break;
case kEmp : status="empty " ;break;
case kEdg : status="edge (fit)" ;break;
case kSi1 : status="size 1 (cog)" ;break;
case kNoLoc: status="no LocMax(fit)" ;break;
case kAbn : status="Abnormal fit " ;break;
case kBig : status="Big Clu(>100) " ;break;
default: status="??????" ;break;
}
Double_t ratio=0;
if(Q()>0&&QRaw()>0) ratio = Q()/QRaw()*100;
Printf("%sCLU: ch=%i (%7.3f,%7.3f) Q=%8.3f Qraw=%8.3f(%3.0f%%) Size=%2i DimBox=%i LocMax=%i Chi2=%7.3f %s",
opt,Ch(),X(),Y(),Q(),QRaw(),ratio,Size(),fBox,fNlocMax,fChi2,status);
if(fDigs) fDigs->Print();
}
Int_t AliHMPIDCluster::Solve(TClonesArray *pCluLst,Int_t *pSigmaCut, Bool_t isTryUnfold)
{
const Int_t kMaxLocMax=6;
CoG();
Int_t iCluCnt=pCluLst->GetEntriesFast();
Int_t rawSize = Size();
if(rawSize>100) {
fSt = kBig;
} else if(isTryUnfold==kFALSE) {
fSt = kNot;
} else if(rawSize==1) {
fSt = kSi1;
}
if(rawSize>100 || isTryUnfold==kFALSE || rawSize==1) {
SetClusterParams(fXX,fYY,fCh);
new ((*pCluLst)[iCluCnt++]) AliHMPIDCluster(*this);
return 1;
}
Double_t arglist[10];
Int_t ierflg = 0;
TVirtualFitter* fitter = TVirtualFitter::Fitter(this,3*6);
arglist[0] = -1;
ierflg = fitter->ExecuteCommand("SET PRI", arglist, 1);
ierflg = fitter->ExecuteCommand("SET NOW", arglist, 0);
arglist[0] = 1;
ierflg = fitter->ExecuteCommand("SET GRA", arglist, 1);
fitter->SetFCN(AliHMPIDCluster::FitFunc);
fNlocMax=0;
for(Int_t iDig1=0;iDig1<rawSize;iDig1++) {
AliHMPIDDigit *pDig1 = Dig(iDig1);
Int_t iCnt = 0;
for(Int_t iDig2=0;iDig2<rawSize;iDig2++) {
if(iDig1==iDig2) continue;
AliHMPIDDigit *pDig2 = Dig(iDig2);
Int_t dist = TMath::Sign(Int_t(pDig1->PadChX()-pDig2->PadChX()),1)+TMath::Sign(Int_t(pDig1->PadChY()-pDig2->PadChY()),1);
if(dist==1)
if(pDig2->Q()>=pDig1->Q()) iCnt++;
}
if(iCnt==0&&fNlocMax<kMaxLocMax){
Double_t xStart=pDig1->LorsX();Double_t yStart=pDig1->LorsY();
Double_t xMin=xStart-fParam->SizePadX();
Double_t xMax=xStart+fParam->SizePadX();
Double_t yMin=yStart-fParam->SizePadY();
Double_t yMax=yStart+fParam->SizePadY();
ierflg = fitter->SetParameter(3*fNlocMax ,Form("x%i",fNlocMax),xStart,0.1,xMin,xMax);
ierflg = fitter->SetParameter(3*fNlocMax+1,Form("y%i",fNlocMax),yStart,0.1,yMin,yMax);
ierflg = fitter->SetParameter(3*fNlocMax+2,Form("q%i",fNlocMax),pDig1->Q(),0.1,0,10000);
fNlocMax++;
}
}
if ( fNlocMax == 0) {
fNlocMax = 1;
fSt=kNoLoc;
SetClusterParams(fXX,fYY,fCh);
new ((*pCluLst)[iCluCnt++]) AliHMPIDCluster(*this);
return fNlocMax;
}
if ( fNlocMax >= kMaxLocMax) {
SetClusterParams(fXX,fYY,fCh);
fSt = kMax; new ((*pCluLst)[iCluCnt++]) AliHMPIDCluster(*this);
} else {
arglist[0] = 500;
arglist[1] = 1.;
ierflg = fitter->ExecuteCommand("SIMPLEX",arglist,2);
if (!ierflg)
fitter->ExecuteCommand("MIGRAD" ,arglist,2);
if(ierflg) {
Double_t strategy=2;
ierflg = fitter->ExecuteCommand("SET STR",&strategy,1);
if(!ierflg) {
ierflg = fitter->ExecuteCommand("SIMPLEX",arglist,2);
if (!ierflg)
fitter->ExecuteCommand("MIGRAD" ,arglist,2);
}
}
if(ierflg) fSt=kAbn;
Double_t dummy; char sName[80];
Double_t edm, errdef;
Int_t nvpar, nparx;
for(Int_t i=0;i<fNlocMax;i++){
fitter->GetParameter(3*i ,sName, fXX, fErrX , dummy, dummy);
fitter->GetParameter(3*i+1 ,sName, fYY, fErrY , dummy, dummy);
fitter->GetParameter(3*i+2 ,sName, fQ, fErrQ , dummy, dummy);
fitter->GetStats(fChi2, edm, errdef, nvpar, nparx);
if(fNlocMax>1)FindClusterSize(i,pSigmaCut);
if(fSt!=kAbn) {
if(fNlocMax!=1)fSt=kUnf;
if(fNlocMax==1&&fSt!=kNoLoc) fSt=kLo1;
if ( !IsInPc()) fSt = kEdg;
if(fSt==kNoLoc) fNlocMax=0;
}
SetClusterParams(fXX,fYY,fCh);
new ((*pCluLst)[iCluCnt++]) AliHMPIDCluster(*this);
if(fNlocMax>1)SetSize(rawSize);
}
}
return fNlocMax;
}
void AliHMPIDCluster::FindClusterSize(Int_t i,Int_t *pSigmaCut)
{
Int_t size = 0;
for(Int_t iDig=0;iDig<Size();iDig++) {
AliHMPIDDigit *pDig = Dig(iDig);
Int_t iCh = pDig->Ch();
Double_t qPad = Q()*pDig->IntMathieson(X(),Y());
AliDebug(1,Form("Chamber %i X %i Y %i SigmaCut %i pad %i qpadMath %8.2f qPadRaw %8.2f Qtotal %8.2f cluster n.%i",iCh,pDig->PadChX(),pDig->PadChY(),
pSigmaCut[iCh],iDig,qPad,pDig->Q(),QRaw(),i));
if(qPad>pSigmaCut[iCh]) size++;
}
AliDebug(1,Form(" Calculated size %i",size));
if(size>0) SetSize(size);
}