]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - MUON/AliMUONClusterFinderSimpleFit.cxx
code cleanup: renaming functions; adding prototype code for later development; no...
[u/mrichter/AliRoot.git] / MUON / AliMUONClusterFinderSimpleFit.cxx
index fa77e5e5e0cf2c224bfd8b4c78b5c13dbf703261..3886a5f18d954982b4ce44a3e798409d4d98cdbb 100644 (file)
 
 #include "AliLog.h"
 #include "AliMpDEManager.h"
-#include "AliMpStationType.h"
 #include "AliMUONCluster.h"
 #include "AliMUONConstants.h"
-#include "AliMUONDigit.h"
+#include "AliMUONVDigit.h"
 #include "AliMUONMathieson.h"
 #include "AliMUONPad.h"
-#include "AliMUONClusterFinderCOG.h"
 #include "AliMpArea.h"
-#include "TClonesArray.h"
 #include "TObjArray.h"
 #include "TVector2.h"
 #include "TVirtualFitter.h"
 #include "TF1.h"
+#include "AliMUONVDigitStore.h"
+#include <Riostream.h>
 
+//-----------------------------------------------------------------------------
 /// \class AliMUONClusterFinderSimpleFit
 ///
 /// Basic cluster finder 
@@ -44,8 +44,9 @@
 /// FIXME: this one is still at the developping stage...
 ///
 /// \author Laurent Aphecetche
+//-----------------------------------------------------------------------------
 
-/// \nocond CLASSIMP
+/// \cond CLASSIMP
 ClassImp(AliMUONClusterFinderSimpleFit)
 /// \endcond
 
@@ -66,33 +67,34 @@ namespace
     
     f = 0.0;
     Float_t qTot = cluster->Charge();
-    Float_t chargeCorrel[] = { cluster->Charge(0)/qTot, cluster->Charge(1)/qTot };
+//    Float_t chargeCorrel[] = { cluster->Charge(0)/qTot, cluster->Charge(1)/qTot };
+//    Float_t qRatio[] = { 1.0/par[2], par[2] };
     
     for ( Int_t i = 0 ; i < cluster->Multiplicity(); ++i )
     {
       AliMUONPad* pad = cluster->Pad(i);
       // skip pads w/ saturation or other problem(s)
       if ( pad->Status() ) continue; 
-      TVector2 lowerLeft(TVector2(par[0],par[1])-pad->Position()-pad->Dimensions());
+      TVector2 lowerLeft = TVector2(par[0],par[1]) - pad->Position() - pad->Dimensions();
       TVector2 upperRight(lowerLeft + pad->Dimensions()*2.0);
-      Float_t estimatedCharge = 
-        qTot*mathieson->IntXY(lowerLeft.X(),lowerLeft.Y(),
-                              upperRight.X(),upperRight.Y());
-      estimatedCharge *= chargeCorrel[pad->Cathode()];
-      Float_t actualCharge = pad->Charge();
+      Float_t estimatedCharge = mathieson->IntXY(lowerLeft.X(),lowerLeft.Y(),
+                                                 upperRight.X(),upperRight.Y());
+//      estimatedCharge *= 2/(1+qRatio[pad->Cathode()]);
+      Float_t actualCharge = pad->Charge()/qTot;
       
       Float_t delta = (estimatedCharge - actualCharge);
       
-      f += delta*delta/qTot;    
+      f += delta*delta;    
     }  
   }
 }
 
 //_____________________________________________________________________________
-AliMUONClusterFinderSimpleFit::AliMUONClusterFinderSimpleFit()
+AliMUONClusterFinderSimpleFit::AliMUONClusterFinderSimpleFit(AliMUONVClusterFinder* clusterFinder)
 : AliMUONVClusterFinder(),
-fClusterFinder(0x0),
-fMathieson(0x0)
+fClusterFinder(clusterFinder),
+fMathieson(0x0),
+fLowestClusterCharge(0)
 {
   /// ctor
 }
@@ -107,39 +109,22 @@ AliMUONClusterFinderSimpleFit::~AliMUONClusterFinderSimpleFit()
 
 //_____________________________________________________________________________
 Bool_t 
-AliMUONClusterFinderSimpleFit::Prepare(const AliMpVSegmentation* segmentations[2],
-                                       TClonesArray* digits[2])
+AliMUONClusterFinderSimpleFit::Prepare(Int_t detElemId,
+                                       TObjArray* pads[2],
+                                       const AliMpArea& area)
 {
   /// Prepare for clustering
 
   // FIXME: should we get the Mathieson from elsewhere ?
   
   // Find out the DetElemId
-  Int_t detElemId(-1);
-  
-  for ( Int_t i = 0; i < 2; ++i )
-  {
-    AliMUONDigit* d = static_cast<AliMUONDigit*>(digits[i]->First());
-    if (d)
-    {
-      detElemId = d->DetElemId();
-      break;
-    }
-  }
-  
-  if ( detElemId < 0 )
-  {
-    AliWarning("Could not find DE. Probably no digits at all ?");
-    return kFALSE;
-  }
-  
-  AliMpStationType stationType = AliMpDEManager::GetStationType(detElemId);
+  AliMq::Station12Type stationType = AliMpDEManager::GetStation12Type(detElemId);
   
   Float_t kx3 = AliMUONConstants::SqrtKx3();
   Float_t ky3 = AliMUONConstants::SqrtKy3();
   Float_t pitch = AliMUONConstants::Pitch();
   
-  if ( stationType == kStation1 )
+  if ( stationType == AliMq::kStation1 )
   {
     kx3 = AliMUONConstants::SqrtKx3St1();
     ky3 = AliMUONConstants::SqrtKy3St1();
@@ -153,9 +138,7 @@ AliMUONClusterFinderSimpleFit::Prepare(const AliMpVSegmentation* segmentations[2
   fMathieson->SetSqrtKx3AndDeriveKx2Kx4(kx3);
   fMathieson->SetSqrtKy3AndDeriveKy2Ky4(ky3);
 
-  delete fClusterFinder;
-  fClusterFinder = new AliMUONClusterFinderCOG;
-  return fClusterFinder->Prepare(segmentations,digits);
+  return fClusterFinder->Prepare(detElemId,pads,area);
 }
 
 //_____________________________________________________________________________
@@ -170,7 +153,7 @@ AliMUONClusterFinderSimpleFit::NextCluster()
   {
     ComputePosition(*cluster);
 
-    if ( cluster->Charge() < 7 )
+    if ( cluster->Charge() < fLowestClusterCharge )
     {
       // skip that one
       return NextCluster();
@@ -201,23 +184,30 @@ AliMUONClusterFinderSimpleFit::ComputePosition(AliMUONCluster& cluster)
   Float_t stepX = 0.01; // cm
   Float_t stepY = 0.01; // cm
   
-  Double_t arg(1);//0);
+  Double_t arg(-1); // disable printout
   
-  fitter->ExecuteCommand("SET PRINT",&arg,1); // disable printout
+  fitter->ExecuteCommand("SET PRINT",&arg,1);
   
   fitter->SetParameter(0,"cluster X position",xCOG,stepX,0,0);
   fitter->SetParameter(1,"cluster Y position",yCOG,stepY,0,0);
-//  fitter->SetParameter(2,"charge ratio",1.0,0.01,0.1,10);
   
   TObjArray userObjects;
   
   userObjects.Add(&cluster);
   userObjects.Add(fMathieson);
   
-//  fitter->SetUserFunc(&userObjects);
   fitter->SetObjectFit(&userObjects);
   
-  fitter->ExecuteCommand("MIGRAD",0,0);
+  Int_t val = fitter->ExecuteCommand("MIGRAD",0,0);
+  AliDebug(1,Form("ExecuteCommand returned value=%d",val));
+  if ( val ) 
+  {
+    // fit failed. Using COG results, with big errors
+    AliWarning("Fit failed. Using COG results for cluster=");
+    StdoutToAliWarning(cluster.Print());
+    cluster.SetPosition(TVector2(xCOG,yCOG),TVector2(TMath::Abs(xCOG),TMath::Abs(yCOG)));
+    cluster.SetChi2(1E3);
+  }
   
   Double_t results[] = { fitter->GetParameter(0),
     fitter->GetParameter(1) };
@@ -228,14 +218,17 @@ AliMUONClusterFinderSimpleFit::ComputePosition(AliMUONCluster& cluster)
   cluster.SetPosition(TVector2(results[0],results[1]),
                       TVector2(errors[0],errors[1]));
   
-  TF1* func = static_cast<TF1*>(fitter->GetUserFunc());
-  Double_t chi2 = 0;
-  if ( func ) chi2 = func->GetChisquare();
+  Double_t amin, edm, errdef;
+  Int_t nvpar, nparx;
+  
+  fitter->GetStats(amin, edm, errdef, nvpar, nparx);
+
+  Double_t chi2 = amin;
   
-  AliDebug(1,Form("Cluster fitted to (x,y)=(%e,%e) (xerr,yerr)=(%e,%e) chi2=%e",
+  AliDebug(1,Form("Cluster fitted to (x,y)=(%e,%e) (xerr,yerr)=(%e,%e) \n chi2=%e ndf=%d",
                   results[0],results[1],
-                  errors[0],errors[1],chi2));
-//  cluster.SetChi2(chi2/fitter->GetNumberFreeParameters());
+                  errors[0],errors[1],chi2,fitter->GetNumberFreeParameters()));
+  cluster.SetChi2(chi2);
 }