MUON/AliMUONClusterFinderSimpleFit.cxx

   1 /**************************************************************************
   2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15
  16 // $Id$
  17
  18 #include "AliMUONClusterFinderSimpleFit.h"
  19
  20 #include "AliLog.h"
  21 #include "AliMpDEManager.h"
  22 #include "AliMpStationType.h"
  23 #include "AliMUONCluster.h"
  24 #include "AliMUONConstants.h"
  25 #include "AliMUONVDigit.h"
  26 #include "AliMUONMathieson.h"
  27 #include "AliMUONPad.h"
  28 #include "AliMpArea.h"
  29 #include "TClonesArray.h"
  30 #include "TObjArray.h"
  31 #include "TVector2.h"
  32 #include "TVirtualFitter.h"
  33 #include "TF1.h"
  34 #include "AliMUONVDigitStore.h"
  35 #include <Riostream.h>
  36
  37 //-----------------------------------------------------------------------------
  38 /// \class AliMUONClusterFinderSimpleFit
  39 ///
  40 /// Basic cluster finder
  41 ///
  42 /// We simply use AliMUONPreClusterFinder to get basic cluster,
  43 /// and then we try to fit the charge repartition using a Mathieson
  44 /// distribution, varying the position.
  45 ///
  46 /// FIXME: this one is still at the developping stage...
  47 ///
  48 /// \author Laurent Aphecetche
  49 //-----------------------------------------------------------------------------
  50
  51 /// \cond CLASSIMP
  52 ClassImp(AliMUONClusterFinderSimpleFit)
  53 /// \endcond
  54
  55 namespace
  56 {
  57   //___________________________________________________________________________
  58   void
  59   FitFunction(Int_t& /*notused*/, Double_t* /*notused*/,
  60               Double_t& f, Double_t* par,
  61               Int_t /*notused*/)
  62   {
  63     /// Chi2 Function to minimize: Mathieson charge distribution in 2 dimensions
  64
  65     TObjArray* userObjects = static_cast<TObjArray*>(TVirtualFitter::GetFitter()->GetObjectFit());
  66
  67     AliMUONCluster* cluster = static_cast<AliMUONCluster*>(userObjects->At(0));
  68     AliMUONMathieson* mathieson = static_cast<AliMUONMathieson*>(userObjects->At(1));
  69
  70     f = 0.0;
  71     Float_t qTot = cluster->Charge();
  72 //    Float_t chargeCorrel[] = { cluster->Charge(0)/qTot, cluster->Charge(1)/qTot };
  73 //    Float_t qRatio[] = { 1.0/par[2], par[2] };
  74
  75     for ( Int_t i = 0 ; i < cluster->Multiplicity(); ++i )
  76     {
  77       AliMUONPad* pad = cluster->Pad(i);
  78       // skip pads w/ saturation or other problem(s)
  79       if ( pad->Status() ) continue;
  80       TVector2 lowerLeft = TVector2(par[0],par[1]) - pad->Position() - pad->Dimensions();
  81       TVector2 upperRight(lowerLeft + pad->Dimensions()*2.0);
  82       Float_t estimatedCharge = mathieson->IntXY(lowerLeft.X(),lowerLeft.Y(),
  83                                                  upperRight.X(),upperRight.Y());
  84 //      estimatedCharge *= 2/(1+qRatio[pad->Cathode()]);
  85       Float_t actualCharge = pad->Charge()/qTot;
  86
  87       Float_t delta = (estimatedCharge - actualCharge);
  88
  89       f += delta*delta;
  90     }
  91   }
  92 }
  93
  94 //_____________________________________________________________________________
  95 AliMUONClusterFinderSimpleFit::AliMUONClusterFinderSimpleFit(AliMUONVClusterFinder* clusterFinder)
  96 : AliMUONVClusterFinder(),
  97 fClusterFinder(clusterFinder),
  98 fMathieson(0x0)
  99 {
 100   /// ctor
 101 }
 102
 103 //_____________________________________________________________________________
 104 AliMUONClusterFinderSimpleFit::~AliMUONClusterFinderSimpleFit()
 105 {
 106   /// dtor
 107   delete fClusterFinder;
 108   delete fMathieson;
 109 }
 110
 111 //_____________________________________________________________________________
 112 Bool_t
 113 AliMUONClusterFinderSimpleFit::Prepare(Int_t detElemId,
 114                                        TClonesArray* pads[2],
 115                                        const AliMpArea& area)
 116 {
 117   /// Prepare for clustering
 118
 119   // FIXME: should we get the Mathieson from elsewhere ?
 120
 121   // Find out the DetElemId
 122   AliMp::StationType stationType = AliMpDEManager::GetStationType(detElemId);
 123
 124   Float_t kx3 = AliMUONConstants::SqrtKx3();
 125   Float_t ky3 = AliMUONConstants::SqrtKy3();
 126   Float_t pitch = AliMUONConstants::Pitch();
 127
 128   if ( stationType == AliMp::kStation1 )
 129   {
 130     kx3 = AliMUONConstants::SqrtKx3St1();
 131     ky3 = AliMUONConstants::SqrtKy3St1();
 132     pitch = AliMUONConstants::PitchSt1();
 133   }
 134
 135   delete fMathieson;
 136   fMathieson = new AliMUONMathieson;
 137
 138   fMathieson->SetPitch(pitch);
 139   fMathieson->SetSqrtKx3AndDeriveKx2Kx4(kx3);
 140   fMathieson->SetSqrtKy3AndDeriveKy2Ky4(ky3);
 141
 142   return fClusterFinder->Prepare(detElemId,pads,area);
 143 }
 144
 145 //_____________________________________________________________________________
 146 AliMUONCluster*
 147 AliMUONClusterFinderSimpleFit::NextCluster()
 148 {
 149   /// Returns next cluster
 150
 151   if ( !fClusterFinder ) return 0x0;
 152   AliMUONCluster* cluster = fClusterFinder->NextCluster();
 153   if ( cluster )
 154   {
 155     ComputePosition(*cluster);
 156
 157     if ( cluster->Charge() < 7 )
 158     {
 159       // skip that one
 160       return NextCluster();
 161     }
 162   }
 163   return cluster;
 164 }
 165
 166 //_____________________________________________________________________________
 167 void
 168 AliMUONClusterFinderSimpleFit::ComputePosition(AliMUONCluster& cluster)
 169 {
 170   /// Compute the position of the given cluster, by fitting a Mathieson
 171   /// charge distribution to it
 172
 173   TVirtualFitter* fitter = TVirtualFitter::Fitter(0,2);
 174   fitter->SetFCN(FitFunction);
 175
 176   if ( cluster.Multiplicity() < 3 ) return;
 177
 178   // We try a Mathieson fit, starting
 179   // with the center-of-gravity estimate as a first guess
 180   // for the cluster center.
 181
 182   Double_t xCOG = cluster.Position().X();
 183   Double_t yCOG = cluster.Position().Y();
 184
 185   Float_t stepX = 0.01; // cm
 186   Float_t stepY = 0.01; // cm
 187
 188   Double_t arg(-1); // disable printout
 189
 190   fitter->ExecuteCommand("SET PRINT",&arg,1);
 191
 192   fitter->SetParameter(0,"cluster X position",xCOG,stepX,0,0);
 193   fitter->SetParameter(1,"cluster Y position",yCOG,stepY,0,0);
 194
 195   TObjArray userObjects;
 196
 197   userObjects.Add(&cluster);
 198   userObjects.Add(fMathieson);
 199
 200   fitter->SetObjectFit(&userObjects);
 201
 202   Int_t val = fitter->ExecuteCommand("MIGRAD",0,0);
 203   AliDebug(1,Form("ExecuteCommand returned value=%d",val));
 204   if ( val )
 205   {
 206     // fit failed. Using COG results, with big errors
 207     AliWarning("Fit failed. Using COG results for cluster=");
 208     StdoutToAliWarning(cluster.Print());
 209     cluster.SetPosition(TVector2(xCOG,yCOG),TVector2(TMath::Abs(xCOG),TMath::Abs(yCOG)));
 210     cluster.SetChi2(1E3);
 211   }
 212
 213   Double_t results[] = { fitter->GetParameter(0),
 214     fitter->GetParameter(1) };
 215
 216   Double_t errors[] = { fitter->GetParError(0),
 217     fitter->GetParError(1) };
 218
 219   cluster.SetPosition(TVector2(results[0],results[1]),
 220                       TVector2(errors[0],errors[1]));
 221
 222   Double_t amin, edm, errdef;
 223   Int_t nvpar, nparx;
 224
 225   fitter->GetStats(amin, edm, errdef, nvpar, nparx);
 226
 227   Double_t chi2 = amin;
 228
 229   AliDebug(1,Form("Cluster fitted to (x,y)=(%e,%e) (xerr,yerr)=(%e,%e) \n chi2=%e ndf=%d",
 230                   results[0],results[1],
 231                   errors[0],errors[1],chi2,fitter->GetNumberFreeParameters()));
 232   cluster.SetChi2(chi2);
 233 }
 234
 235
 236