#include "AliMUONCluster.h"
#include "AliMUONPad.h"
#include "AliMUONPad.h"
-#include "AliMpStationType.h"
#include "AliMUONConstants.h"
#include "AliMpDEManager.h"
#include "AliMUONMathieson.h"
+#include "AliMpEncodePair.h"
+
#include "AliLog.h"
#include <TClonesArray.h>
#include <TMath.h>
#include <TMatrixD.h>
#include <TObjArray.h>
-#include <TROOT.h>
#include <TRandom.h>
#include <Riostream.h>
//_____________________________________________________________________________
AliMUONClusterSplitterMLEM::AliMUONClusterSplitterMLEM(Int_t detElemId,
- TObjArray* fPixArray)
+ TObjArray* pixArray,
+ Double_t lowestPixelCharge,
+ Double_t lowestPadCharge,
+ Double_t lowestClusterCharge)
: TObject(),
-fPixArray(fPixArray),
+fPixArray(pixArray),
fMathieson(0x0),
fDetElemId(detElemId),
fNpar(0),
fQtot(0),
fnCoupled(0),
-fDebug(0)
+fDebug(0),
+fLowestPixelCharge(lowestPixelCharge),
+fLowestPadCharge(lowestPadCharge),
+fLowestClusterCharge(lowestClusterCharge)
{
/// Constructor
- AliMp::StationType stationType = AliMpDEManager::GetStationType(fDetElemId);
+ AliMq::Station12Type stationType = AliMpDEManager::GetStation12Type(fDetElemId);
Float_t kx3 = AliMUONConstants::SqrtKx3();
Float_t ky3 = AliMUONConstants::SqrtKy3();
Float_t pitch = AliMUONConstants::Pitch();
- if ( stationType == AliMp::kStation1 )
+ if ( stationType == AliMq::kStation1 )
{
kx3 = AliMUONConstants::SqrtKx3St1();
ky3 = AliMUONConstants::SqrtKy3St1();
cont1 = mlem->GetCellContent(j,i);
if (mode && cont1 > cont) continue;
used[(i-1)*nx+j-1] = kTRUE;
- if (cont1 < 0.5) continue;
+ if (cont1 < fLowestPixelCharge) continue;
if (pix) pix->Add(BinToPix(mlem,j,i));
else {
pixPtr = new AliMUONPad (mlem->GetXaxis()->GetBinCenter(j),
// Compare pixel and bin positions
for (Int_t i = 0; i < nPix; ++i) {
pixPtr = (AliMUONPad*) fPixArray->UncheckedAt(i);
- if (pixPtr->Charge() < 0.5) continue;
+ if (pixPtr->Charge() < fLowestPixelCharge) continue;
if (TMath::Abs(pixPtr->Coord(0)-xc)<1.e-4 && TMath::Abs(pixPtr->Coord(1)-yc)<1.e-4)
{
//return (TObject*) pixPtr;
{
AliMUONPad* pad = cluster.Pad(j);
//if ( pad->Status() !=1 || pad->IsSaturated() ) continue;
- if ( pad->Status() != AliMUONClusterFinderMLEM::fgkUseForFit) continue;
+ if ( pad->Status() != AliMUONClusterFinderMLEM::GetUseForFitFlag() ||
+ pad->Charge() == 0 ) continue;
if (iflag == 0) {
if ( pad->IsReal() ) npads++; // exclude virtual pads
qTot += pad->Charge();
}
//_____________________________________________________________________________
-Double_t AliMUONClusterSplitterMLEM::Param2Coef(Int_t icand, Double_t coef, Double_t *par)
+Double_t AliMUONClusterSplitterMLEM::Param2Coef(Int_t icand, Double_t coef, Double_t *par) const
{
/// Extract hit contribution scale factor from fit parameters
Int_t
AliMUONClusterSplitterMLEM::Fit(const AliMUONCluster& cluster,
Int_t iSimple, Int_t nfit,
- Int_t *clustFit, TObjArray **clusters,
+ const Int_t *clustFit, TObjArray **clusters,
Double_t *parOk,
- TObjArray& clusterList)
+ TObjArray& clusterList, TH2 *mlem)
{
/// Steering function and fitting procedure for the fit of pad charge distribution
// AliDebug(2,Form("iSimple=%d nfit=%d",iSimple,nfit));
- TH2D *mlem = (TH2D*) gROOT->FindObject("mlem");
Double_t xmin = mlem->GetXaxis()->GetXmin() - mlem->GetXaxis()->GetBinWidth(1);
Double_t xmax = mlem->GetXaxis()->GetXmax() + mlem->GetXaxis()->GetBinWidth(1);
Double_t ymin = mlem->GetYaxis()->GetXmin() - mlem->GetYaxis()->GetBinWidth(1);
AliMUONPad* pad = cluster.Pad(i);
if ( !pad->IsReal() ) ++nVirtual;
//if ( pad->Status() !=1 || pad->IsSaturated() ) continue;
- if ( pad->Status() != AliMUONClusterFinderMLEM::fgkUseForFit ) continue;
+ if ( pad->Status() != AliMUONClusterFinderMLEM::GetUseForFitFlag() ) continue;
if ( pad->IsReal() )
{
++npads;
// Get number of pads in X and Y
//const Int_t kStatusToTest(1);
- const Int_t kStatusToTest(AliMUONClusterFinderMLEM::fgkUseForFit);
+ const Int_t kStatusToTest(AliMUONClusterFinderMLEM::GetUseForFitFlag());
- AliMpIntPair nofPads = cluster.NofPads(kStatusToTest);
- Int_t nInX = nofPads.GetFirst();
- Int_t nInY = nofPads.GetSecond();
+ Long_t nofPads = cluster.NofPads(kStatusToTest);
+ Int_t nInX = AliMp::PairFirst(nofPads);
+ Int_t nInY = AliMp::PairSecond(nofPads);
if (fDebug) {
Int_t npadOK = 0;
for (Int_t j = 0; j < cluster.Multiplicity(); ++j) {
AliMUONPad *pad = cluster.Pad(j);
//if (pad->Status() == 1 && !pad->IsSaturated()) npadOK++;
- if (pad->Status() == AliMUONClusterFinderMLEM::fgkUseForFit && !pad->IsSaturated()) npadOK++;
+ if (pad->Status() == AliMUONClusterFinderMLEM::GetUseForFitFlag() && !pad->IsSaturated()) npadOK++;
}
cout << " Number of pads to fit: " << npadOK << endl;
cout << " nInX and Y: " << nInX << " " << nInY << endl;
Int_t nfitMax = 3;
nfitMax = TMath::Min (nfitMax, (npads + 1) / 3);
if (nfitMax > 1) {
- if (nInX < 3 && nInY < 3 || nInX == 3 && nInY < 3 || nInX < 3 && nInY == 3) nfitMax = 1; // not enough pads in each direction
+ if (((nInX < 3) && (nInY < 3)) || ((nInX == 3) && (nInY < 3)) || ((nInX < 3) && (nInY == 3))) nfitMax = 1; // not enough pads in each direction
}
if (nfit > nfitMax) nfit = nfitMax;
AliMUONPad *pixPtr;
Int_t npxclu;
Double_t cont, cmax = 0, xseed = 0, yseed = 0, errOk[8], qq = 0;
+
+ for ( int i = 0; i < 8; ++i ) errOk[i]=0.0;
+
Double_t xyseed[3][2], qseed[3], xyCand[3][2] = {{0},{0}}, sigCand[3][2] = {{0},{0}};
for (Int_t ifit = 1; ifit <= nfit0; ++ifit)
if (fDebug) cout << xyCand[0][0] << " " << xyCand[0][1] << " " << sigCand[0][0] << " " << sigCand[0][1] << endl;
Int_t nDof, maxSeed[3];//, nMax = 0;
-
+
+ if ( nfit0 < 0 || nfit0 > 3 ) {
+ AliErrorStream() << "Wrong nfit0 value: " << nfit0 << endl;
+ return nfit;
+ }
TMath::Sort(nfit0, qseed, maxSeed, kTRUE); // in decreasing order
Double_t step[3]={0.01,0.002,0.02}, fmin, chi2o = 9999, chi2n;
- Double_t *gin = 0, func0, func1, param[8], step0[8];
+ Double_t *gin = 0, func0, func1, param[8]={0}, step0[8]={0};
Double_t param0[2][8]={{0},{0}}, deriv[2][8]={{0},{0}};
- Double_t shift[8], stepMax, derMax, parmin[8], parmax[8], func2[2], shift0;
- Double_t delta[8], scMax, dder[8], estim, shiftSave = 0;
+ Double_t shift[8]={0}, stepMax, derMax, parmin[8]={0}, parmax[8]={0}, func2[2]={0}, shift0;
+ Double_t delta[8]={0}, scMax, dder[8], estim, shiftSave = 0;
Int_t min, max, nCall = 0, nLoop, idMax = 0, iestMax = 0, nFail;
Double_t rad, dist[3] = {0};
param[j] = parOk[j];
}
param[fNpar] = 0.6;
- parmin[fNpar] = 0;
+ parmin[fNpar] = 1E-9;
parmax[fNpar++] = 1;
}
{
shift[j] = 0;
}
- else if (deriv[min][j]*deriv[!min][j] > 0 && TMath::Abs(deriv[min][j]) > TMath::Abs(deriv[!min][j])
- || TMath::Abs(deriv[0][j]-deriv[1][j]) < 1.e-3 || TMath::Abs(dder[j]) < 1.e-6)
+ else if (((deriv[min][j]*deriv[!min][j] > 0) && (TMath::Abs(deriv[min][j]) > TMath::Abs(deriv[!min][j])))
+ || (TMath::Abs(deriv[0][j]-deriv[1][j]) < 1.e-3) || (TMath::Abs(dder[j]) < 1.e-6))
{
shift[j] = -TMath::Sign (shift[j], (func2[0]-func2[1]) * (param0[0][j]-param0[1][j]));
if (min == max)
}
} // for (Int_t j=0; j<fNpar;
- if (estim < 1 && derMax < 2 || nLoop > 150) break; // minimum was found
+ if (((estim < 1) && (derMax < 2)) || nLoop > 150) break; // minimum was found
nLoop++;
//if (fNpar > 2) cout << param0[min][fNpar-3] << " " << chi2n * (1+TMath::Min(1-param0[min][fNpar-3],0.25)) << endl;
//if (chi2n*1.2+1.e-6 > chi2o )
- if (fNpar > 2 && (chi2n > chi2o || iseed == nfit-1
- && chi2n * (1+TMath::Min(1-param0[min][fNpar-3],0.25)) > chi2o))
+ if (fNpar > 2 && (chi2n > chi2o || ((iseed == nfit-1)
+ && (chi2n * (1+TMath::Min(1-param0[min][fNpar-3],0.25)) > chi2o))))
{ fNpar -= 3; break; }
// Save parameters and errors
if (fDebug) {
for (Int_t i=0; i<fNpar; ++i) {
if (i == 4 || i == 7) {
- if (i == 7 || i == 4 && fNpar < 7) cout << parOk[i] << endl;
+ if ((i == 7) || ((i == 4) && (fNpar < 7))) cout << parOk[i] << endl;
else cout << parOk[i] * (1-parOk[7]) << endl;
continue;
}
// Double_t /*sigy*/,
// Double_t /*dist*/)
- if ( coef*fQtot >= 14 )
+ if ( coef*fQtot >= fLowestClusterCharge )
{
//AZ AliMUONCluster* cluster1 = new AliMUONCluster();
AliMUONCluster* cluster1 = new AliMUONCluster(cluster);
cluster1->SetCharge(coef*fQtot,coef*fQtot);
cluster1->SetPosition(TVector2(parOk[indx],parOk[indx+1]),TVector2(sigCand[0][0],sigCand[0][1]));
- cluster1->SetChi2(dist[TMath::LocMin(nfit,dist)]);
+ //cluster1->SetChi2(dist[TMath::LocMin(nfit,dist)]);
+ Int_t idx = TMath::LocMin(nfit,dist);
+ if ( idx < 0 || idx > 2 ) {
+ AliErrorStream() << "Wrong index value: " << idx << endl;
+ return nfit;
+ }
+ cluster1->SetChi2(dist[idx]);
// FIXME: we miss some information in this cluster, as compared to
// the original AddRawCluster code.
indx = (i-1)*nx + j - 1;
if (used[indx]) continue;
cont = mlem->GetCellContent(j,i);
- if (cont < 0.5) continue;
+ if (cont < fLowestPixelCharge) continue;
pix = new TObjArray(20);
used[indx] = 1;
pix->Add(BinToPix(mlem,j,i));
// Flag clusters for fit
nForFit = 0;
- while (minGroup[nForFit] >= 0 && nForFit < 3)
+ while (nForFit < 3 && minGroup[nForFit] >= 0)
{
if (fDebug) cout << clustNumb[minGroup[nForFit]] << " ";
clustFit[nForFit] = clustNumb[minGroup[nForFit]];
AliMUONPad* pad = cluster.Pad(j);
//if ( pad->Status()==1 ) pad->SetStatus(0);
//if ( pad->Status()==-9) pad->SetStatus(-5);
- if ( pad->Status() == AliMUONClusterFinderMLEM::fgkUseForFit ||
- pad->Status() == AliMUONClusterFinderMLEM::fgkCoupled)
- pad->SetStatus(AliMUONClusterFinderMLEM::fgkZero);
+ if ( pad->Status() == AliMUONClusterFinderMLEM::GetUseForFitFlag() ||
+ pad->Status() == AliMUONClusterFinderMLEM::GetCoupledFlag())
+ pad->SetStatus(AliMUONClusterFinderMLEM::GetZeroFlag());
}
// Merge the failed cluster candidates (with too few pads to fit) with
// the one with the strongest coupling
else
{
// Do the fit
- nfit = Fit(cluster,0, nForFit, clustFit, clusters, parOk, clusterList);
+ nfit = Fit(cluster,0, nForFit, clustFit, clusters, parOk, clusterList, mlem);
if (nfit == 0) {
//cout << " (nfit == 0) " << fNpar << " " << cluster.Multiplicity() << endl;
fNpar = 0; // should be 0 by itself but just in case ...
AliMUONPad* pad = cluster.Pad(j);
//if ( pad->Status()==1 ) pad->SetStatus(-2);
//if ( pad->Status()==-9) pad->SetStatus(-5);
- if ( pad->Status() == AliMUONClusterFinderMLEM::fgkUseForFit )
- pad->SetStatus(AliMUONClusterFinderMLEM::fgkModified);
+ if ( pad->Status() == AliMUONClusterFinderMLEM::GetUseForFitFlag() )
+ pad->SetStatus(AliMUONClusterFinderMLEM::GetModifiedFlag());
}
// Sort the clusters (move to the right the used ones)
{
AliMUONPad* pad = cluster.Pad(j);
//if ( pad->Status() != -2) continue;
- if ( pad->Status() != AliMUONClusterFinderMLEM::fgkModified) continue;
+ if ( pad->Status() != AliMUONClusterFinderMLEM::GetModifiedFlag()) continue;
for (Int_t iclust=0; iclust<nCoupled; ++iclust)
{
indx = clustNumb[iclust];
}
}
//pad->SetStatus(-8);
- pad->SetStatus(AliMUONClusterFinderMLEM::fgkOver);
+ pad->SetStatus(AliMUONClusterFinderMLEM::GetOverFlag());
} // for (Int_t j=0; j<npad;
} // if (nCoupled > 3)
} // while (nCoupled > 0)
void
AliMUONClusterSplitterMLEM::Merge(const AliMUONCluster& cluster,
Int_t nForFit, Int_t nCoupled,
- Int_t *clustNumb, Int_t *clustFit,
+ const Int_t *clustNumb, const Int_t *clustFit,
TObjArray **clusters,
TMatrixD& aijcluclu, TMatrixD& aijclupad)
{
{
AliMUONPad* pad = cluster.Pad(j);
//if ( pad->Status() < 0 && pad->Status() != -5 ) continue;// exclude used pads
- if ( pad->Status() != AliMUONClusterFinderMLEM::fgkZero) continue;// exclude used pads
+ if ( pad->Status() != AliMUONClusterFinderMLEM::GetZeroFlag()) continue;// exclude used pads
aijclupad(imax,j) += aijclupad(indx,j);
aijclupad(indx,j) = 0;
}
//_____________________________________________________________________________
Double_t
-AliMUONClusterSplitterMLEM::MinGroupCoupl(Int_t nCoupled, Int_t *clustNumb,
- TMatrixD& aijcluclu, Int_t *minGroup)
+AliMUONClusterSplitterMLEM::MinGroupCoupl(Int_t nCoupled, const Int_t *clustNumb,
+ const TMatrixD& aijcluclu, Int_t *minGroup)
{
/// Find group of clusters with minimum coupling to all the others
Int_t
AliMUONClusterSplitterMLEM::SelectPad(const AliMUONCluster& cluster,
Int_t nCoupled, Int_t nForFit,
- Int_t *clustNumb, Int_t *clustFit,
- TMatrixD& aijclupad)
+ const Int_t *clustNumb, const Int_t *clustFit,
+ const TMatrixD& aijclupad)
{
/// Select pads for fit. If too many coupled clusters, find pads giving
/// the strongest coupling with the rest of clusters and exclude them from the fit.
++nOK; // pad to be used in fit
}
*/
- if ( pad->Status() != AliMUONClusterFinderMLEM::fgkZero
+ if ( pad->Status() != AliMUONClusterFinderMLEM::GetZeroFlag()
|| pad->IsSaturated() ) continue; // used pads and overflows
- pad->SetStatus(AliMUONClusterFinderMLEM::fgkUseForFit);
+ pad->SetStatus(AliMUONClusterFinderMLEM::GetUseForFitFlag());
++nOK; // pad to be used in fit
if (nCoupled > 3)
if (padpix[j] < fgkCouplMin) continue;
aaa += padpix[j];
//cluster.Pad(j)->SetStatus(-1); // exclude pads with strong coupling to the other clusters
- cluster.Pad(j)->SetStatus(AliMUONClusterFinderMLEM::fgkCoupled); // exclude pads with strong coupling to the other clusters
+ cluster.Pad(j)->SetStatus(AliMUONClusterFinderMLEM::GetCoupledFlag()); // exclude pads with strong coupling to the other clusters
nOK--;
}
delete [] padpix;
{
AliMUONPad* pad = cluster.Pad(j);
//if ( pad->Status() != -1 ) continue;
- if ( pad->Status() != AliMUONClusterFinderMLEM::fgkCoupled ) continue;
+ if ( pad->Status() != AliMUONClusterFinderMLEM::GetCoupledFlag() ) continue;
if (fNpar != 0)
{
charge = 0;
} // if (fNpar != 0)
//if (pad->Charge() > 6 /*fgkZeroSuppression*/) pad->SetStatus(0);
- if (pad->Charge() > 6 /*fgkZeroSuppression*/) pad->SetStatus(AliMUONClusterFinderMLEM::fgkZero);
+ if (pad->Charge() > fLowestPadCharge) pad->SetStatus(AliMUONClusterFinderMLEM::GetZeroFlag());
// return pad for further using // FIXME: remove usage of zerosuppression here
- else pad->SetStatus(AliMUONClusterFinderMLEM::fgkOver); // do not use anymore
+ else pad->SetStatus(AliMUONClusterFinderMLEM::GetOverFlag()); // do not use anymore
} // for (Int_t j=0;
}