From 9bbd7f601bc641bfee69069ded14d09a8020a215 Mon Sep 17 00:00:00 2001 From: ivana Date: Mon, 25 Jun 2007 12:35:53 +0000 Subject: [PATCH] Making MLEM cluster finder working properly (reproduces results from AZ cluster finder for single upsilons). (Sasha) --- MUON/AliMUONClusterFinderMLEM.cxx | 85 +++++++++++++++++-------- MUON/AliMUONClusterSplitterMLEM.cxx | 96 +++++++++++++++++------------ MUON/AliMUONClusterSplitterMLEM.h | 3 + 3 files changed, 116 insertions(+), 68 deletions(-) diff --git a/MUON/AliMUONClusterFinderMLEM.cxx b/MUON/AliMUONClusterFinderMLEM.cxx index 30808ea3f14..889e971aeff 100644 --- a/MUON/AliMUONClusterFinderMLEM.cxx +++ b/MUON/AliMUONClusterFinderMLEM.cxx @@ -71,7 +71,7 @@ fZpad(0.0), fReco(1), fCathBeg(0), fPixArray(new TObjArray(20)), -fDebug(1), +fDebug(0), fPlot(plot), fTimers(new TObjArray(kLast)), fSplitter(0x0), @@ -95,6 +95,7 @@ fNAddVirtualPads(0) t->Start(kTRUE); t->Stop(); } + if (fPlot) fDebug = 1; } //_____________________________________________________________________________ @@ -146,6 +147,7 @@ AliMUONClusterFinderMLEM::Prepare(const AliMpVSegmentation* segmentations[2], delete fSplitter; fSplitter = new AliMUONClusterSplitterMLEM(fDetElemId,fPixArray); + fSplitter->SetDebug(fDebug); // find out current event number, and reset the cluster number fEventNumber = AliRunLoader::GetRunLoader()->GetEventNumber(); @@ -163,15 +165,9 @@ AliMUONClusterFinderMLEM::NextCluster() { /// Return next cluster - ++fClusterNumber; - // 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(o); - } + TObject* o = fClusterList.At(++fClusterNumber); + if ( o != 0x0 ) return static_cast(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 @@ -189,6 +185,7 @@ AliMUONClusterFinderMLEM::NextCluster() } fClusterList.Delete(); // reset the list of clusters for this pre-cluster + fClusterNumber = -1; //AZ WorkOnPreCluster(); @@ -214,10 +211,28 @@ AliMUONClusterFinderMLEM::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; jMultiplicity(); ++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 ) @@ -309,8 +324,6 @@ AliMUONClusterFinderMLEM::CheckPrecluster(const AliMUONCluster& origCluster) AliMUONCluster* cluster = static_cast(origCluster.Clone()); - cluster->Sort(); - AliDebug(2,"Start of CheckPreCluster="); StdoutToAliDebug(2,cluster->Print("full")); @@ -382,6 +395,8 @@ AliMUONClusterFinderMLEM::CheckPreclusterTwoCathodes(AliMUONCluster* cluster) if (nFlags > 0) { // not all pads overlap. + if (fDebug) cout << " nFlags: " << nFlags << endl; + TObjArray toBeRemoved; for (Int_t i=0; iPad(i); @@ -391,16 +406,22 @@ AliMUONClusterFinderMLEM::CheckPreclusterTwoCathodes(AliMUONCluster* cluster) // 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(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(); @@ -843,6 +864,7 @@ AliMUONClusterFinderMLEM::Plot(const char* basename) { /// Make a plot and save it as png + return; //AZ if (!fPlot) return; TCanvas* c = new TCanvas("MLEM","MLEM",800,600); @@ -1245,6 +1267,7 @@ void AliMUONClusterFinderMLEM::Mlem(AliMUONCluster& cluster, // Do iterations for (Int_t ipix=0; ipixSetChargeBackup(0); // Correct each pixel if (probi[ipix] < 0.01) continue; // skip "invisible" pixel Double_t sum = 0; @@ -1281,9 +1304,14 @@ void AliMUONClusterFinderMLEM::Mlem(AliMUONCluster& cluster, 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; } @@ -1672,7 +1700,6 @@ void AliMUONClusterFinderMLEM::AddVirtualPad(AliMUONCluster& cluster) 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; @@ -1868,7 +1895,8 @@ void AliMUONClusterFinderMLEM::AddVirtualPad(AliMUONCluster& cluster) 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 @@ -1891,15 +1919,16 @@ void AliMUONClusterFinderMLEM::AddVirtualPad(AliMUONCluster& cluster) // 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; @@ -1914,7 +1943,8 @@ void AliMUONClusterFinderMLEM::AddVirtualPad(AliMUONCluster& cluster) // 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) { @@ -1934,12 +1964,13 @@ void AliMUONClusterFinderMLEM::AddVirtualPad(AliMUONCluster& cluster) // 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; } @@ -1991,7 +2022,7 @@ void AliMUONClusterFinderMLEM::Simple(AliMUONCluster& cluster) 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); } diff --git a/MUON/AliMUONClusterSplitterMLEM.cxx b/MUON/AliMUONClusterSplitterMLEM.cxx index 0f077dc3486..770feff792e 100644 --- a/MUON/AliMUONClusterSplitterMLEM.cxx +++ b/MUON/AliMUONClusterSplitterMLEM.cxx @@ -43,6 +43,7 @@ #include #include #include +#include /// \cond CLASSIMP ClassImp(AliMUONClusterSplitterMLEM) @@ -59,7 +60,8 @@ fMathieson(0x0), fDetElemId(detElemId), fNpar(0), fQtot(0), -fnCoupled(0) +fnCoupled(0), +fDebug(0) { /// Constructor @@ -196,8 +198,8 @@ AliMUONClusterSplitterMLEM::Fcn1(const AliMUONCluster& cluster, for (Int_t j=0; j< cluster.Multiplicity(); ++j) { AliMUONPad* pad = cluster.Pad(j); - if ( pad->Status() != 1 ) continue; - if ( pad->DX() > 0 ) npads++; // exclude virtual pads + if ( pad->Status() !=1 || pad->IsSaturated() ) continue; + if ( pad->IsReal() ) npads++; // exclude virtual pads qTot += pad->Charge(); // c.fXyq[2][j]; charge = 0; for (Int_t i=fNpar/3; i>=0; --i) @@ -217,7 +219,7 @@ AliMUONClusterSplitterMLEM::Fcn1(const AliMUONCluster& cluster, coef = i==1 ? coef*par[indx+2] : coef - par[7]; } coef = TMath::Max (coef, 0.); - charge += ChargeIntegration(par[indx],par[indx+1],*pad); + charge += ChargeIntegration(par[indx],par[indx+1],*pad) * coef; } charge *= fQtot; delta = charge - pad->Charge(); //c.fXyq[2][j]; @@ -251,12 +253,13 @@ AliMUONClusterSplitterMLEM::Fit(const AliMUONCluster& cluster, // Number of pads to use and number of virtual pads Int_t npads = 0, nVirtual = 0, nfit0 = nfit; + //cluster.Print("full"); for (Int_t i=0; iDX() < 0 ) ++nVirtual; - if ( pad->Status() !=1 ) continue; - if ( pad->DX() > 0 ) + if ( !pad->IsReal() ) ++nVirtual; + if ( pad->Status() !=1 || pad->IsSaturated() ) continue; + if ( pad->IsReal() ) { ++npads; if (yPad > 9999) @@ -325,7 +328,16 @@ AliMUONClusterSplitterMLEM::Fit(const AliMUONCluster& cluster, AliMpIntPair nofPads = cluster.NofPads(kStatusToTest); Int_t nInX = nofPads.GetFirst(); Int_t nInY = nofPads.GetSecond(); - //cout << " nInX and Y: " << nInX << " " << nInY << endl; + + 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++; + } + 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); @@ -375,7 +387,7 @@ AliMUONClusterSplitterMLEM::Fit(const AliMUONCluster& cluster, sigCand[0][0] = sigCand[0][0] > 0 ? TMath::Sqrt (sigCand[0][0]) : 0; sigCand[0][1] = sigCand[0][1]/qq - xyCand[0][1]*xyCand[0][1]; // - ^2 sigCand[0][1] = sigCand[0][1] > 0 ? TMath::Sqrt (sigCand[0][1]) : 0; -// if (fDebug) cout << xyCand[0][0] << " " << xyCand[0][1] << " " << sigCand[0][0] << " " << sigCand[0][1] << endl; + if (fDebug) cout << xyCand[0][0] << " " << xyCand[0][1] << " " << sigCand[0][0] << " " << sigCand[0][1] << endl; Int_t nDof, maxSeed[3];//, nMax = 0; Double_t fmin, chi2o = 9999, chi2n; @@ -386,7 +398,7 @@ AliMUONClusterSplitterMLEM::Fit(const AliMUONCluster& cluster, 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; - Int_t min, max, nCall = 0, memory[8] = {0}, nLoop, idMax = 0, iestMax = 0, nFail; + Int_t min, max, nCall = 0, nLoop, idMax = 0, iestMax = 0, nFail; Double_t rad, dist[3] = {0}; // Try to fit with one-track hypothesis, then 2-track. If chi2/dof is @@ -394,6 +406,7 @@ AliMUONClusterSplitterMLEM::Fit(const AliMUONCluster& cluster, for (Int_t iseed=0; iseed 1) - { - shift[j] *= 2; - } - memory[j]++; - } + if (min == max) + { + if (memory[j] > 1) + { + shift[j] *= 2; + } + memory[j]++; + } + } else { shift[j] = dder[j] != 0 ? -deriv[min][j] / dder[j] : 0; @@ -597,7 +610,7 @@ AliMUONClusterSplitterMLEM::Fit(const AliMUONCluster& cluster, nDof = npads - fNpar + nVirtual; if (!nDof) nDof++; chi2n = fmin / nDof; -// if (fDebug) cout << " Chi2 " << chi2n << " " << fNpar << endl; + if (fDebug) cout << " Chi2 " << chi2n << " " << fNpar << endl; if (chi2n*1.2+1.e-6 > chi2o ) { fNpar -= 3; break; } @@ -648,18 +661,18 @@ AliMUONClusterSplitterMLEM::Fit(const AliMUONCluster& cluster, chi2o = chi2n; if (fmin < 0.1) break; // !!!??? - } // for (Int_t iseed=0; - -// if (fDebug) { -// for (Int_t i=0; i= 14 ) { - AliMUONCluster* cluster = new AliMUONCluster(); + //AZ AliMUONCluster* cluster1 = new AliMUONCluster(); + AliMUONCluster* cluster1 = new AliMUONCluster(cluster); - cluster->SetCharge(coef*fQtot,coef*fQtot); - cluster->SetPosition(TVector2(parOk[indx],parOk[indx+1]),TVector2(sigCand[0][0],sigCand[0][1])); - cluster->SetChi2(dist[TMath::LocMin(nfit,dist)]); + 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)]); // FIXME: we miss some information in this cluster, as compared to // the original AddRawCluster code. AliDebug(2,Form("Adding RawCluster detElemId %4d mult %2d charge %5d (xl,yl)=(%9.6g,%9.6g)", - fDetElemId,cluster->Multiplicity(),(Int_t)cluster->Charge(), - cluster->Position().X(),cluster->Position().Y())); + fDetElemId,cluster1->Multiplicity(),(Int_t)cluster1->Charge(), + cluster1->Position().X(),cluster1->Position().Y())); - clusterList.Add(cluster); + clusterList.Add(cluster1); } // AddRawCluster (parOk[indx], // double x // parOk[indx+1], // double y @@ -1218,7 +1232,7 @@ AliMUONClusterSplitterMLEM::UpdatePads(const AliMUONCluster& cluster, for (Int_t j=0; jStatus() != 1 ) continue; + if ( pad->Status() != -1 ) continue; if (fNpar != 0) { charge = 0; @@ -1240,7 +1254,7 @@ AliMUONClusterSplitterMLEM::UpdatePads(const AliMUONCluster& cluster, coef = i==1 ? coef*par[indx+2] : coef - par[7]; } coef = TMath::Max (coef, 0.); - charge += ChargeIntegration(par[indx],par[indx+1],*pad); + charge += ChargeIntegration(par[indx],par[indx+1],*pad) * coef; } charge *= fQtot; pad->SetCharge(pad->Charge()-charge); diff --git a/MUON/AliMUONClusterSplitterMLEM.h b/MUON/AliMUONClusterSplitterMLEM.h index 1c7dc5bfdc9..96d9ac629fb 100644 --- a/MUON/AliMUONClusterSplitterMLEM.h +++ b/MUON/AliMUONClusterSplitterMLEM.h @@ -73,6 +73,8 @@ public: void UpdatePads(const AliMUONCluster& cluster, Int_t nfit, Double_t *par); + /// Set debug level + void SetDebug (Int_t debug) { fDebug = debug; } private: /// will not be implemented @@ -90,6 +92,7 @@ private: Int_t fNpar; //!< number of fit parameters Double_t fQtot; //!< total charge Int_t fnCoupled; //!< number of coupled pixels ? + Int_t fDebug; //!< debug level ClassDef(AliMUONClusterSplitterMLEM,1) // }; -- 2.43.0