// * provided "as is" without express or implied warranty. *
// **************************************************************************
+#include "AliHMPIDCluster.h" //class header
#include <TVirtualFitter.h> //Solve()
#include <TMinuit.h> //Solve()
#include <TClonesArray.h> //Solve()
#include <TMarker.h> //Draw()
-#include "AliLog.h" //FitFunc()
-#include "AliHMPIDCluster.h" //class header
+#include "AliLog.h" //FindCusterSize()
Bool_t AliHMPIDCluster::fgDoCorrSin=kTRUE;
for(Int_t i=0;i<nPads;i++){ //loop on all pads of the cluster
for(Int_t j=0;j<iNshape;j++){ //Mathiesons loop as all of them may contribute to this pad
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)->Mathieson(par[3*j]-pClu->Dig(i)->LorsX()-0.5*AliHMPIDParam::SizePadX())-
- pClu->Dig(i)->Mathieson(par[3*j]-pClu->Dig(i)->LorsX()+0.5*AliHMPIDParam::SizePadX()))*
- pClu->Dig(i)->IntPartMathi(par[3*j+1],2);
- derivPart[3*j+1][i] += par[3*j+2]*(pClu->Dig(i)->Mathieson(par[3*j+1]-pClu->Dig(i)->LorsY()-0.5*AliHMPIDParam::SizePadY())-
- pClu->Dig(i)->Mathieson(par[3*j+1]-pClu->Dig(i)->LorsY()+0.5*AliHMPIDParam::SizePadY()))*
- pClu->Dig(i)->IntPartMathi(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<iNpars;i++) delete [] derivPart[i]; delete [] derivPart;
}
//---gradient calculations ended
+
+// fit ended. Final calculations
+
}//FitFunction()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
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;
}
if(fDigs) fDigs->Print();
}//Print()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-Int_t AliHMPIDCluster::Solve(TClonesArray *pCluLst,Bool_t isTryUnfold)
+Int_t AliHMPIDCluster::Solve(TClonesArray *pCluLst,Int_t *pSigmaCut, Bool_t isTryUnfold)
{
//This methode is invoked when the cluster is formed to solve it. Solve the cluster means to try to unfold the cluster
//into the local maxima number of clusters. This methode is invoked by AliHMPIDRconstructor::Dig2Clu() on cluster by cluster basis.
const Int_t kMaxLocMax=6; //max allowed number of loc max for fitting
//
CoG(); //First calculate CoG for the given cluster
+
Int_t iCluCnt=pCluLst->GetEntriesFast(); //get current number of clusters already stored in the list by previous operations
- if(isTryUnfold==kFALSE || Size()==1) { //if cluster contains single pad there is no way to improve the knowledge
- fSt = (isTryUnfold)? kSi1: kNot;
- if(fParam->GetInstType()) SetClusterParams(fXX,fYY,fCh);
+
+ Int_t rawSize = Size(); //get current raw cluster size
+
+ if(rawSize>100 || isTryUnfold==kFALSE || Size()==1) { //No deconv if: 1 - big cluster (also avoid no zero suppression!)
+ fSt = (isTryUnfold)? kSi1: kNot; // 2 - flag is set to FALSE
+ if(fParam->GetInstType()) SetClusterParams(fXX,fYY,fCh); // 3 - size = 1
new ((*pCluLst)[iCluCnt++]) AliHMPIDCluster(*this); //add this raw cluster
return 1;
}
//Phase 1. Find number of local maxima. Strategy is to check if the current pad has QDC more then all neigbours. Also find the box contaning the cluster
fNlocMax=0;
- for(Int_t iDig1=0;iDig1<Size();iDig1++) { //first digits loop
+ for(Int_t iDig1=0;iDig1<rawSize;iDig1++) { //first digits loop
AliHMPIDDigit *pDig1 = Dig(iDig1); //take next digit
Int_t iCnt = 0; //counts how many neighbouring pads has QDC more then current one
- for(Int_t iDig2=0;iDig2<Size();iDig2++) { //loop on all digits again
+ for(Int_t iDig2=0;iDig2<rawSize;iDig2++) { //loop on all digits again
if(iDig1==iDig2) continue; //the same digit, no need to compare
AliHMPIDDigit *pDig2 = Dig(iDig2); //take second digit to compare with the first one
//Phase 2. Fit loc max number of Mathiesons or add this current cluster to the list
// case 1 -> no loc max found
if ( fNlocMax == 0) { // case of no local maxima found: pads with same charge...
-
- ierflg = fitter->SetParameter(3*fNlocMax ,Form("x%i",fNlocMax),fXX,0.1,0,0); // Init values taken from CoG() -> fXX,fYY,fQRaw
- ierflg = fitter->SetParameter(3*fNlocMax+1,Form("y%i",fNlocMax),fYY,0.1,0,0); //
- ierflg = fitter->SetParameter(3*fNlocMax+2,Form("q%i",fNlocMax),fQRaw,0.1,0,10000); //
-
fNlocMax = 1;
fSt=kNoLoc;
+ if(fParam->GetInstType()) SetClusterParams(fXX,fYY,fCh); //need to fill the AliCluster3D part
+ new ((*pCluLst)[iCluCnt++]) AliHMPIDCluster(*this); //add new unfolded cluster
+
+ return fNlocMax;
}
// case 2 -> loc max found. Check # of loc maxima
Double_t dummy; char sName[80]; //vars to get results from Minuit
Double_t edm, errdef;
Int_t nvpar, nparx;
-
+
for(Int_t i=0;i<fNlocMax;i++){ //store the local maxima parameters
- fitter->GetParameter(3*i ,sName, fXX, fErrX , dummy, dummy); // X
- fitter->GetParameter(3*i+1 ,sName, fYY, fErrY , dummy, dummy); // Y
+ fitter->GetParameter(3*i ,sName, fXX, fErrX , dummy, dummy); // X
+ fitter->GetParameter(3*i+1 ,sName, fYY, fErrY , dummy, dummy); // Y
fitter->GetParameter(3*i+2 ,sName, fQ, fErrQ , dummy, dummy); // Q
fitter->GetStats(fChi2, edm, errdef, nvpar, nparx); //get fit infos
- if(fSt!=kAbn) {
- if(fNlocMax!=1)fSt=kUnf; // if unfolded
- if(fNlocMax==1&&fSt!=kNoLoc) fSt=kLo1; // if only 1 loc max
- if ( !IsInPc()) fSt = kEdg; // if Out of Pc
- if(fSt==kNoLoc) fNlocMax=0; // if with no loc max (pads with same charge..)
- }
- if(fParam->GetInstType()) SetClusterParams(fXX,fYY,fCh); //need to fill the AliCluster3D part
- new ((*pCluLst)[iCluCnt++]) AliHMPIDCluster(*this); //add new unfolded cluster
-
+
+ if(fNlocMax>1)FindClusterSize(i,pSigmaCut); //find clustersize for deconvoluted clusters
+ //after this call, fSi temporarly is the calculated size. Later is set again
+ //to its original value
+ if(fSt!=kAbn) {
+ if(fNlocMax!=1)fSt=kUnf; // if unfolded
+ if(fNlocMax==1&&fSt!=kNoLoc) fSt=kLo1; // if only 1 loc max
+ if ( !IsInPc()) fSt = kEdg; // if Out of Pc
+ if(fSt==kNoLoc) fNlocMax=0; // if with no loc max (pads with same charge..)
+ }
+ if(fParam->GetInstType()) SetClusterParams(fXX,fYY,fCh); //need to fill the AliCluster3D part
+ new ((*pCluLst)[iCluCnt++]) AliHMPIDCluster(*this); //add new unfolded cluster
+ if(fNlocMax>1)SetSize(rawSize); //Original raw size is set again to its proper value
}
}
}//Solve()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+void AliHMPIDCluster::FindClusterSize(Int_t i,Int_t *pSigmaCut)
+{
+
+//Estimate of the clustersize for a deconvoluted cluster
+ Int_t size = 0;
+ for(Int_t iDig=0;iDig<Size();iDig++) { //digits loop
+ AliHMPIDDigit *pDig = Dig(iDig); //take digit
+ Int_t iCh = pDig->Ch();
+ Double_t qPad = Q()*pDig->IntMathieson(X(),Y()); //pad charge
+ 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); //in case of size == 0, original raw clustersize used
+}