]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - MUON/AliMUONClusterSplitterMLEM.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff
? search:
summary | shortlog | log | commit | commitdiff | tree
raw | inline | side by side
Add a protection to avoid array boundary error in IntSpecGeant
[u/mrichter/AliRoot.git] / MUON / AliMUONClusterSplitterMLEM.cxx
diff --git a/MUON/AliMUONClusterSplitterMLEM.cxx b/MUON/AliMUONClusterSplitterMLEM.cxx
index 44883c1abd865d8a4adfe397fc71e0f2d958f283..acc31137566d2b9bd089460e984c7332e7771fe0 100644 (file)
--- a/MUON/AliMUONClusterSplitterMLEM.cxx
+++ b/MUON/AliMUONClusterSplitterMLEM.cxx
@@ -15,6 +15,7 @@
 
 // $Id$
 
+//-----------------------------------------------------------------------------
 /// \class AliMUONClusterSplitterMLEM
 /// 
 /// Splitter class for the MLEM algorithm...
@@ -23,6 +24,7 @@
 ///
 /// \author Laurent Aphecetche (for the "new" C++ structure) and 
 /// Alexander Zinchenko, JINR Dubna, for the hardcore of it ;-)
+//-----------------------------------------------------------------------------
 
 #include "AliMUONClusterSplitterMLEM.h"
 
@@ -43,6 +45,7 @@
 #include <TObjArray.h>
 #include <TROOT.h>
 #include <TRandom.h>
+#include <Riostream.h>
 
 /// \cond CLASSIMP
 ClassImp(AliMUONClusterSplitterMLEM)
@@ -59,7 +62,8 @@ fMathieson(0x0),
 fDetElemId(detElemId),
 fNpar(0),
 fQtot(0),
-fnCoupled(0)
+fnCoupled(0),
+fDebug(0)
 {
   /// Constructor
   
@@ -87,6 +91,8 @@ fnCoupled(0)
 //_____________________________________________________________________________
 AliMUONClusterSplitterMLEM::~AliMUONClusterSplitterMLEM()
 {
+  /// Destructor
+  
   delete fMathieson;
 }
 
@@ -194,8 +200,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)
@@ -215,7 +221,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];
@@ -249,12 +255,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; i<cluster.Multiplicity(); ++i ) 
   {
     AliMUONPad* pad = cluster.Pad(i);
-    if ( pad->DX() < 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) 
@@ -323,7 +330,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);
@@ -373,7 +389,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]; // <y^2> - <y>^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;
@@ -384,7 +400,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 
@@ -392,6 +408,7 @@ AliMUONClusterSplitterMLEM::Fit(const AliMUONCluster& cluster,
     for (Int_t iseed=0; iseed<nfit; iseed++) 
     {
       
+      Int_t memory[8] = {0};
       if (iseed) 
       { 
         for (Int_t j=0; j<fNpar; j++) 
@@ -484,15 +501,15 @@ AliMUONClusterSplitterMLEM::Fit(const AliMUONCluster& cluster,
                    || 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) 
-          { 
-            if (memory[j] > 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;
@@ -595,7 +612,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; }
       
@@ -646,18 +663,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<fNpar; i++) {
-//        if (i == 4 || i == 7) {
-//          if (i == 7 || i == 4 && fNpar < 7) cout << parOk[i] << endl;
-//          else cout << parOk[i] * (1-parOk[7]) << endl;
-//          continue;
-//        }
-//        cout << parOk[i] << " " << errOk[i] << endl;
-//      }
-//    }
+    } // for (Int_t iseed=0; 
+   
+    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;
+          else cout << parOk[i] * (1-parOk[7]) << endl;
+          continue;
+        }
+        cout << parOk[i] << " " << errOk[i] << endl;
+      }
+    }
     nfit = (fNpar + 1) / 3;
     dist[0] = dist[1] = dist[2] = 0;
     
@@ -706,20 +723,21 @@ AliMUONClusterSplitterMLEM::Fit(const AliMUONCluster& cluster,
       
       if ( coef*fQtot >= 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
@@ -1216,7 +1234,7 @@ AliMUONClusterSplitterMLEM::UpdatePads(const AliMUONCluster& cluster,
   for (Int_t j=0; j<cluster.Multiplicity(); ++j) 
   {
     AliMUONPad* pad = cluster.Pad(j);
-    if ( pad->Status() != 1 ) continue;
+    if ( pad->Status() != -1 ) continue;
     if (fNpar != 0) 
     {
       charge = 0;
@@ -1238,7 +1256,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);