]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - MUON/AliMUONClusterSplitterMLEM.cxx
Added analysis of fakes (Francesco)
[u/mrichter/AliRoot.git] / MUON / AliMUONClusterSplitterMLEM.cxx
index dfc92a7e102575e2adab1edf2e320e60e9b7843c..a758b3474a91c37a4c1e8baa1840d16aef29970c 100644 (file)
 #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>
@@ -44,7 +45,6 @@
 #include <TMath.h>
 #include <TMatrixD.h>
 #include <TObjArray.h>
-#include <TROOT.h>
 #include <TRandom.h>
 #include <Riostream.h>
 
@@ -57,25 +57,31 @@ const Double_t AliMUONClusterSplitterMLEM::fgkCouplMin = 1.e-2; // threshold on
 
 //_____________________________________________________________________________
 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();
@@ -119,7 +125,7 @@ AliMUONClusterSplitterMLEM::AddBin(TH2 *mlem,
       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), 
@@ -164,7 +170,7 @@ AliMUONClusterSplitterMLEM::BinToPix(TH2 *mlem,
   // 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;
@@ -207,7 +213,8 @@ AliMUONClusterSplitterMLEM::Fcn1(const AliMUONCluster& cluster,
   {
     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(); 
@@ -231,7 +238,7 @@ AliMUONClusterSplitterMLEM::Fcn1(const AliMUONCluster& cluster,
 }
 
 //_____________________________________________________________________________
-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
   
@@ -246,15 +253,14 @@ Double_t AliMUONClusterSplitterMLEM::Param2Coef(Int_t icand, Double_t coef, Doub
 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);
@@ -270,7 +276,7 @@ AliMUONClusterSplitterMLEM::Fit(const AliMUONCluster& cluster,
     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;
@@ -334,18 +340,18 @@ AliMUONClusterSplitterMLEM::Fit(const AliMUONCluster& cluster,
   
   // 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;
@@ -354,7 +360,7 @@ AliMUONClusterSplitterMLEM::Fit(const AliMUONCluster& cluster,
   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;
   
@@ -363,6 +369,9 @@ AliMUONClusterSplitterMLEM::Fit(const AliMUONCluster& cluster,
   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) 
@@ -402,14 +411,18 @@ AliMUONClusterSplitterMLEM::Fit(const AliMUONCluster& cluster,
   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};
     
@@ -513,8 +526,8 @@ AliMUONClusterSplitterMLEM::Fit(const AliMUONCluster& cluster,
         {
          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) 
@@ -593,7 +606,7 @@ AliMUONClusterSplitterMLEM::Fit(const AliMUONCluster& cluster,
        }
       } // 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++;
         
@@ -633,8 +646,8 @@ AliMUONClusterSplitterMLEM::Fit(const AliMUONCluster& cluster,
       
     //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
@@ -693,7 +706,7 @@ AliMUONClusterSplitterMLEM::Fit(const AliMUONCluster& cluster,
   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;
       }
@@ -739,14 +752,20 @@ AliMUONClusterSplitterMLEM::Fit(const AliMUONCluster& cluster,
     //                                             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.
@@ -803,7 +822,7 @@ AliMUONClusterSplitterMLEM::Split(const AliMUONCluster& cluster,
       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));
@@ -945,9 +964,9 @@ AliMUONClusterSplitterMLEM::Split(const AliMUONCluster& cluster,
           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
@@ -956,7 +975,7 @@ AliMUONClusterSplitterMLEM::Split(const AliMUONCluster& cluster,
       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 ...
@@ -973,8 +992,8 @@ AliMUONClusterSplitterMLEM::Split(const AliMUONCluster& cluster,
         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)
@@ -1014,7 +1033,7 @@ AliMUONClusterSplitterMLEM::Split(const AliMUONCluster& cluster,
         {
           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];
@@ -1030,7 +1049,7 @@ AliMUONClusterSplitterMLEM::Split(const AliMUONCluster& cluster,
             }
           }
           //pad->SetStatus(-8);
-          pad->SetStatus(AliMUONClusterFinderMLEM::fgkOver);
+          pad->SetStatus(AliMUONClusterFinderMLEM::GetOverFlag());
         } // for (Int_t j=0; j<npad;
       } // if (nCoupled > 3)
     } // while (nCoupled > 0)
@@ -1051,7 +1070,7 @@ AliMUONClusterSplitterMLEM::Split(const AliMUONCluster& cluster,
 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)
 {
@@ -1103,7 +1122,7 @@ AliMUONClusterSplitterMLEM::Merge(const AliMUONCluster& cluster,
     {
       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;
     }
@@ -1113,8 +1132,8 @@ AliMUONClusterSplitterMLEM::Merge(const AliMUONCluster& cluster,
 
 //_____________________________________________________________________________
 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
   
@@ -1199,8 +1218,8 @@ AliMUONClusterSplitterMLEM::MinGroupCoupl(Int_t nCoupled, Int_t *clustNumb,
 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.
@@ -1231,9 +1250,9 @@ AliMUONClusterSplitterMLEM::SelectPad(const AliMUONCluster& cluster,
         ++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) 
@@ -1257,7 +1276,7 @@ AliMUONClusterSplitterMLEM::SelectPad(const AliMUONCluster& cluster,
     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; 
@@ -1278,7 +1297,7 @@ AliMUONClusterSplitterMLEM::UpdatePads(const AliMUONCluster& cluster,
   {
     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;
@@ -1294,9 +1313,9 @@ AliMUONClusterSplitterMLEM::UpdatePads(const AliMUONCluster& cluster,
     } // 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;
 }