/* $Id$ */
+//-----------------------------------------------------------------------------
/// \class AliMUONClusterFinderMLEM
///
/// Clusterizer class based on the Expectation-Maximization algorithm
///
/// \author Laurent Aphecetche (for the "new" C++ structure) and
/// Alexander Zinchenko, JINR Dubna, for the hardcore of it ;-)
+//-----------------------------------------------------------------------------
#include "AliMUONClusterFinderMLEM.h"
#include "AliLog.h"
#include "AliMUONCluster.h"
#include "AliMUONClusterSplitterMLEM.h"
-#include "AliMUONDigit.h"
+#include "AliMUONVDigit.h"
#include "AliMUONPad.h"
#include "AliMUONPreClusterFinder.h"
#include "AliMpPad.h"
#include "AliMpVPadIterator.h"
#include "AliMpVSegmentation.h"
#include "AliRunLoader.h"
-
+#include "AliMUONVDigitStore.h"
#include <Riostream.h>
#include <TH2.h>
#include <TMinuit.h>
#include <TStopwatch.h>
#include <TMath.h>
#include <TROOT.h>
+//#include "AliCodeTimer.h"
/// \cond CLASSIMP
ClassImp(AliMUONClusterFinderMLEM)
TMinuit* AliMUONClusterFinderMLEM::fgMinuit = 0x0;
//_____________________________________________________________________________
-AliMUONClusterFinderMLEM::AliMUONClusterFinderMLEM(Bool_t plot)
+AliMUONClusterFinderMLEM::AliMUONClusterFinderMLEM(Bool_t plot, AliMUONVClusterFinder* clusterFinder)
: AliMUONVClusterFinder(),
-fPreClusterFinder(new AliMUONPreClusterFinder),
+fPreClusterFinder(clusterFinder),
fPreCluster(0x0),
fClusterList(),
fEventNumber(0),
fReco(1),
fCathBeg(0),
fPixArray(new TObjArray(20)),
-fDebug(1),
+fDebug(0),
fPlot(plot),
fTimers(new TObjArray(kLast)),
fSplitter(0x0),
t->Start(kTRUE);
t->Stop();
}
+ if (fPlot) fDebug = 1;
}
//_____________________________________________________________________________
//_____________________________________________________________________________
Bool_t
AliMUONClusterFinderMLEM::Prepare(const AliMpVSegmentation* segmentations[2],
- TClonesArray* digits[2])
+ const AliMUONVDigitStore& digitStore)
{
/// Prepare for clustering
+// AliCodeTimerAuto("")
for ( Int_t i = 0; i < 2; ++i )
{
// Find out the DetElemId
fDetElemId = -1;
- for ( Int_t i = 0; i < 2; ++i )
+ TIter next(digitStore.CreateIterator());
+ AliMUONVDigit* d = static_cast<AliMUONVDigit*>(next());
+
+ if (d)
{
- AliMUONDigit* d = static_cast<AliMUONDigit*>(digits[i]->First());
- if (d)
- {
- fDetElemId = d->DetElemId();
- break;
- }
+ fDetElemId = d->DetElemId();
}
if ( fDetElemId < 0 )
delete fSplitter;
fSplitter = new AliMUONClusterSplitterMLEM(fDetElemId,fPixArray);
+ fSplitter->SetDebug(fDebug);
// find out current event number, and reset the cluster number
fEventNumber = AliRunLoader::GetRunLoader()->GetEventNumber();
fClusterNumber = -1;
fClusterList.Delete();
-// AliDebug(3,Form("EVT %d DE %d",fEventNumber,fDetElemId));
+ AliDebug(3,Form("EVT %d DE %d",fEventNumber,fDetElemId));
- return fPreClusterFinder->Prepare(segmentations,digits);
+ return fPreClusterFinder->Prepare(segmentations,digitStore);
}
//_____________________________________________________________________________
AliMUONClusterFinderMLEM::NextCluster()
{
/// Return next cluster
-
- ++fClusterNumber;
+// AliCodeTimerAuto("clustering + pre-clustering")
// if the list of clusters is not void, pick one from there
- if ( fClusterList.GetLast() >= 0 )
- {
- TObject* o = fClusterList.At(0);
- fClusterList.RemoveAt(0);
- return static_cast<AliMUONCluster*>(o);
- }
+ TObject* o = fClusterList.At(++fClusterNumber);
+ if ( o != 0x0 ) return static_cast<AliMUONCluster*>(o);
//FIXME : at this point, must check whether we've used all the digits
//from precluster : if not, let the preclustering know about those unused
}
fClusterList.Delete(); // reset the list of clusters for this pre-cluster
+ fClusterNumber = -1; //AZ
WorkOnPreCluster();
/// Starting from a precluster, builds a pixel array, and then
/// extract clusters from this array
+ // Set saturation flag - it is not set if working directly with MC digits (w/out
+ // creating raw data) !!!
+ for (Int_t j = 0; j < fPreCluster->Multiplicity(); ++j) {
+ AliMUONPad* pad = fPreCluster->Pad(j);
+ if (pad->IsSaturated()) break;
+ if (pad->Charge() > fgkSaturation-1) pad->SetSaturated(kTRUE); //FIXME : remove usage of fgkSaturation
+ }
+
+ if (fDebug) {
+ cout << " *** Event # " << AliRunLoader::GetRunLoader()->GetEventNumber()
+ << " det. elem.: " << fDetElemId << endl;
+ for (Int_t j=0; j<fPreCluster->Multiplicity(); ++j) {
+ AliMUONPad* pad = fPreCluster->Pad(j);
+ printf(" bbb %3d %1d %8.4f %8.4f %8.4f %8.4f %6.1f %3d %3d %2d %1d %1d \n",
+ j, pad->Cathode(), pad->Coord(0), pad->Coord(1), pad->DX()*2, pad->DY()*2,
+ pad->Charge(), pad->Ix(), pad->Iy(), pad->Status(), pad->IsReal(), pad->IsSaturated());
+ }
+ }
+
AliMUONCluster* cluster = CheckPrecluster(*fPreCluster);
-
if (!cluster) return kFALSE;
-
+
BuildPixArray(*cluster);
if ( fPixArray->GetLast() < 0 )
AliMUONCluster* cluster = static_cast<AliMUONCluster*>(origCluster.Clone());
- cluster->Sort();
-
AliDebug(2,"Start of CheckPreCluster=");
- StdoutToAliDebug(2,cluster->Print("full"));
+// StdoutToAliDebug(2,cluster->Print("full"));
// Check if one-cathode precluster
Int_t i1 = cluster->Multiplicity(0) ? 0 : 1;
{
/// Check single-cathode precluster
AliWarning("Reimplement me!");
- AliDebug(2,"End of CheckPreClusterOneCathode=");
- StdoutToAliDebug(2,cluster->Print("full"));
+ AliDebug(2,"End of CheckPreClusterOneCathode=");
+// StdoutToAliDebug(2,cluster->Print("full"));
return cluster;
}
if (nFlags > 0)
{
// not all pads overlap.
+ if (fDebug) cout << " nFlags: " << nFlags << endl;
+ TObjArray toBeRemoved;
for (Int_t i=0; i<npad; ++i)
{
AliMUONPad* pad = cluster->Pad(i);
// Check for edge effect (missing pads on the _other_ cathode)
AliMpPad mpPad = fSegmentation[cath1]->PadByPosition(pad->Position(),kFALSE);
if (!mpPad.IsValid()) continue;
+ if (nFlags == 1 && pad->Charge() < fgkZeroSuppression * 3) continue;
AliDebug(2,Form("Releasing the following pad : de,cath,ix,iy %d,%d,%d,%d charge %e",
fDetElemId,pad->Cathode(),pad->Ix(),pad->Iy(),pad->Charge()));
- cluster->RemovePad(pad);
+ toBeRemoved.AddLast(pad);
+ //AZ cluster->RemovePad(pad);
fPreCluster->Pad(i)->Release();
- --npad;
+ //AZ --npad;
+ }
+ for ( Int_t i = 0; i <= toBeRemoved.GetLast(); ++i )
+ {
+ cluster->RemovePad(static_cast<AliMUONPad*>(toBeRemoved.At(i)));
}
}
// Check correlations of cathode charges
- if ( !cluster->IsSaturated() && cluster->ChargeAsymmetry()*2 > 1 )
+ if ( !cluster->IsSaturated() && cluster->ChargeAsymmetry() > 1 )
{
// big difference
Int_t cathode = cluster->MaxRawChargeCathode();
delete[] flags;
AliDebug(2,"End of CheckPreClusterTwoCathodes=");
- StdoutToAliDebug(2,cluster->Print("full"));
+// StdoutToAliDebug(2,cluster->Print("full"));
return cluster;
}
Int_t nPix = fPixArray->GetLast()+1;
-// AliDebug(2,Form("nPix after BuildPixArray=%d",nPix));
+ AliDebug(2,Form("nPix after BuildPixArray=%d",nPix));
Double_t xPadMin(1E9);
Double_t yPadMin(1E9);
if (pixPtr->Charge() < 1)
{
AliDebug(2,Form("Removing pixel %d with charge<1 : ",i));
- StdoutToAliDebug(2,pixPtr->Print());
+// StdoutToAliDebug(2,pixPtr->Print());
RemovePixel(i);
}
}
fPixArray->Compress();
nPix = fPixArray->GetEntriesFast();
-// AliDebug(2,Form("nPix after AdjustPixel=%d",nPix));
+ AliDebug(2,Form("nPix after AdjustPixel=%d",nPix));
if ( nPix > cluster.Multiplicity() )
{
-// AliDebug(2,Form("Will trim number of pixels to number of pads"));
+ AliDebug(2,Form("Will trim number of pixels to number of pads"));
// Too many pixels - sort and remove pixels with the lowest signal
fPixArray->Sort();
// StdoutToAliDebug(2,cout << "End of BuildPixelArray:" << endl;
// fPixArray->Print(););
- CheckOverlaps();//FIXME : this is for debug only. Remove it.
+// CheckOverlaps();//FIXME : this is for debug only. Remove it.
}
//_____________________________________________________________________________
{
/// From a single-cathode cluster, build the pixel array
-// AliDebug(2,Form("cluster.Multiplicity=%d",cluster.Multiplicity()));
+ AliDebug(2,Form("cluster.Multiplicity=%d",cluster.Multiplicity()));
for ( Int_t j=0; j<cluster.Multiplicity(); ++j)
{
{
/// From a two-cathodes cluster, build the pixel array
-// AliDebug(2,Form("cluster.Multiplicity=%d",cluster.Multiplicity()));
+ AliDebug(2,Form("cluster.Multiplicity=%d",cluster.Multiplicity()));
Int_t i1 = cluster.Pad(0)->Cathode();
Int_t i2 = TMath::Even(i1);
{
/// Make a plot and save it as png
+ return; //AZ
if (!fPlot) return;
TCanvas* c = new TCanvas("MLEM","MLEM",800,600);
// coef is the charge (given by Mathieson integral) on pad, assuming
// the Mathieson is center at pixel.
coef[indx1] = fSplitter->ChargeIntegration(pixPtr->Coord(0), pixPtr->Coord(1), *pad);
-// AliDebug(2,Form("pad=(%d,%d,%e,%e,%e,%e) pix=(%e,%e,%e,%e) coef %e",
-// pad->Ix(),pad->Iy(),
-// pad->X(),pad->Y(),
-// pad->DX(),pad->DY(),
-// pixPtr->Coord(0),pixPtr->Coord(1),
-// pixPtr->Size(0),pixPtr->Size(1),
-// coef[indx1]));
+ AliDebug(2,Form("pad=(%d,%d,%e,%e,%e,%e) pix=(%e,%e,%e,%e) coef %e",
+ pad->Ix(),pad->Iy(),
+ pad->X(),pad->Y(),
+ pad->DX(),pad->DY(),
+ pixPtr->Coord(0),pixPtr->Coord(1),
+ pixPtr->Size(0),pixPtr->Size(1),
+ coef[indx1]));
probi[ipix] += coef[indx1];
}
Int_t nPix = fPixArray->GetLast()+1;
AliDebug(2,Form("nPix=%d iSimple=%d, precluster=",nPix,iSimple));
- StdoutToAliDebug(2,cluster.Print("full"););
+// StdoutToAliDebug(2,cluster.Print("full"););
if ( nPix < 0 )
{
Double_t* probi(0x0);
Int_t lc(0); // loop counter (for debug)
- Plot("mlem.start");
+// Plot("mlem.start");
while (1)
{
AliDebug(2,Form("lc %d nPix %d(%d) npadTot %d npadOK %d",lc,nPix,fPixArray->GetLast()+1,npadTot,npadOK));
AliDebug(2,Form("EVT%d PixArray=",fEventNumber));
- StdoutToAliDebug(2,fPixArray->Print("","full"));
+// StdoutToAliDebug(2,fPixArray->Print("","full"));
coef = new Double_t [npadTot*nPix];
probi = new Double_t [nPix];
{
AliMUONPad* pixel = Pixel(ipix);
AliDebug(2,Form("Setting the following pixel to invisible as its probi<0.01:"));
- StdoutToAliDebug(2,cout << Form(" -- ipix %3d --- "); pixel->Print(););
+// StdoutToAliDebug(2,cout << Form(" -- ipix %3d --- "); pixel->Print(););
pixel->SetCharge(0); // "invisible" pixel
}
}
nPix = fPixArray->GetEntriesFast();
AliDebug(2,Form("After shift:"));
- StdoutToAliDebug(2,fPixArray->Print("","full"););
- Plot(Form("mlem.lc%d",lc+1));
+// StdoutToAliDebug(2,fPixArray->Print("","full"););
+// Plot(Form("mlem.lc%d",lc+1));
AliDebug(2,Form(" xyCOG=%9.6f %9.6f xylim=%9.6f,%9.6f,%9.6f,%9.6f",
xyCOG[0],xyCOG[1],
j = TMath::Even (i/2);
p->SetCoord(j, xyCOG[j]);
AliDebug(2,Form("Adding pixel on the edge (i=%d) ",i));
- StdoutToAliDebug(2,cout << " ---- ";
- p->Print("corners"););
+// StdoutToAliDebug(2,cout << " ---- ";
+// p->Print("corners"););
fPixArray->Add(p);
++nPix;
}
} // while (1)
AliDebug(2,Form("At the end of while loop nPix=%d : ",fPixArray->GetLast()+1));
- StdoutToAliDebug(2,fPixArray->Print("","full"););
+// StdoutToAliDebug(2,fPixArray->Print("","full"););
// remove pixels with low signal or low visibility
// Cuts are empirical !!!
Mlem(cluster,coef,probi,2);
AliDebug(2,Form("Before splitting nPix=%d EVT %d DE %d",fPixArray->GetLast()+1,fEventNumber,fDetElemId));
- StdoutToAliDebug(2,fPixArray->Print("","full"););
- Plot("mlem.beforesplit");
+// StdoutToAliDebug(2,fPixArray->Print("","full"););
+// Plot("mlem.beforesplit");
// Update histogram
for (Int_t i=0; i<nPix; ++i)
// Do iterations
for (Int_t ipix=0; ipix<nPix; ++ipix)
{
+ Pixel(ipix)->SetChargeBackup(0);
// Correct each pixel
if (probi[ipix] < 0.01) continue; // skip "invisible" pixel
Double_t sum = 0;
continue;
}
- if (coef[indx] > 1.e-6)
+ if (sum1 > 1.e-6)
{
sum += pad->Charge()*coef[indx]/sum1;
}
AliMUONPad* pixPtr = Pixel(ipix);
if (probi1[ipix] > 1.e-6)
{
- pixPtr->SetCharge(pixPtr->Charge()*sum/probi1[ipix]);
+ //AZ pixPtr->SetCharge(pixPtr->Charge()*sum/probi1[ipix]);
+ pixPtr->SetChargeBackup(pixPtr->Charge()*sum/probi1[ipix]);
}
} // for (Int_t ipix=0;
+ for (Int_t i = 0; i < nPix; ++i) {
+ AliMUONPad* pixPtr = Pixel(i);
+ pixPtr->SetCharge(pixPtr->ChargeBackup());
+ }
} // for (Int_t iter=0;
delete [] probi1;
}
Int_t nInX = -1, nInY;
PadsInXandY(cluster,nInX, nInY);
++fNClusters;
- if (fDebug) cout << " Chamber: " << fDetElemId / 100 - 1 << " " << nInX << " " << nInY << endl;
// Add virtual pad only if number of pads per direction == 2
if (nInX != 2 && nInY != 2) return;
AliMpPad pad = fSegmentation[cath]->PadByIndices(AliMpIntPair(ix,iy),kTRUE);
// fXyq[0][npads] = pad.Position().X();
// fXyq[1][npads] = pad.Position().Y();
- AliMUONPad muonPad(fDetElemId, cath, ix, iy, pad.Position().X(), pad.Position().Y(), 0, 0, 0);
+ AliMUONPad muonPad(fDetElemId, cath, ix, iy, pad.Position().X(), pad.Position().Y(),
+ pad.Dimensions().X(), pad.Dimensions().Y(), 0);
// fSegmentation[cath]->GetPadC(ix, iy, fXyq[0][npads], fXyq[1][npads], zpad);
// if (fXyq[0][npads] > 1.e+5) continue; // temporary fix
if (muonPad.Coord(0) > 1.e+5) continue; // temporary fix
// fXyq[3][npads] = -2; // flag
// fPadIJ[2][npads] = ix;
// fPadIJ[3][npads] = iy;
- muonPad.SetSize(0,-2.); //flag
+ muonPad.SetReal(kFALSE);
// fnPads[1]++;
// iAddX = npads;
iAddX = 1;
- //AliDebug(1,Form("Add virtual pad in X %f %f %f %3d %3d \n",
- // fXyq[2][npads], fXyq[0][npads], fXyq[1][npads], ix, iy));
- //muonPad.Charge(), muonPad.Coord(0), muonPad.Coord(1), ix, iy));
- if (fDebug) printf(" ***** Add virtual pad in X ***** %f %f %f %3d %3d \n",
- muonPad.Charge(), muonPad.Coord(0), muonPad.Coord(1), ix, iy);
+ if (fDebug) printf(" ***** Add virtual pad in X ***** %f %f %f %3d %3d %f %f \n",
+ muonPad.Charge(), muonPad.Coord(0), muonPad.Coord(1), ix, iy,
+ muonPad.DX(), muonPad.DY());
cluster.AddPad(muonPad); // add pad to the cluster
ix1 = ix0;
continue;
// fXyq[0][npads] = pad.Position().X();
// fXyq[1][npads] = pad.Position().Y();
// fSegmentation[cath]->GetPadC(ix, iy, fXyq[0][npads], fXyq[1][npads], zpad);
- AliMUONPad muonPad(fDetElemId, cath, ix, iy, pad.Position().X(), pad.Position().Y(), 0, 0, 0);
+ AliMUONPad muonPad(fDetElemId, cath, ix, iy, pad.Position().X(), pad.Position().Y(),
+ pad.Dimensions().X(), pad.Dimensions().Y(), 0);
if (iy1 == iy0) continue;
//if (iPad && iy1 == iy0) continue;
if (maxpad[0][0] < 0 || mirror && maxpad[1][0] >= 0) {
// fXyq[3][npads] = -2; // flag
// fPadIJ[2][npads] = ix;
// fPadIJ[3][npads] = iy;
- muonPad.SetSize(0,-2.); //flag
+ muonPad.SetReal(kFALSE);
// fnPads[1]++;
// iAddY = npads;
iAddY = 1;
- if (fDebug) printf(" ***** Add virtual pad in Y ***** %f %f %f %3d %3d \n",
- muonPad.Charge(), muonPad.Coord(0), muonPad.Coord(1), ix, iy);
+ if (fDebug) printf(" ***** Add virtual pad in Y ***** %f %f %f %3d %3d %f %f \n",
+ muonPad.Charge(), muonPad.Coord(0), muonPad.Coord(1), ix, iy,
+ muonPad.DX(), muonPad.DY());
cluster.AddPad(muonPad); // add pad to the cluster
iy1 = iy0;
}
for (Int_t i = 0; i < cluster.Multiplicity(); ++i)
{
AliMUONPad* pad = cluster.Pad(i);
- if ( pad->Charge() > fgkSaturation-1) //FIXME : remove usage of fgkSaturation
+ if ( pad->IsSaturated())
{
pad->SetStatus(-9);
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff