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 "AliMUONClusterFinderCOG.h"
  29 #include "AliMpArea.h"
  30 #include "TClonesArray.h"
  31 #include "TObjArray.h"
  32 #include "TVector2.h"
  33 #include "TVirtualFitter.h"
  34 #include "TF1.h"
  35 #include "AliMUONVDigitStore.h"
  36 #include <Riostream.h>
  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 /// \cond CLASSIMP
  51 ClassImp(AliMUONClusterFinderSimpleFit)
  52 /// \endcond
  53
  54 namespace
  55 {
  56   //___________________________________________________________________________
  57   void
  58   FitFunction(Int_t& /*notused*/, Double_t* /*notused*/,
  59               Double_t& f, Double_t* par,
  60               Int_t /*notused*/)
  61   {
  62     /// Chi2 Function to minimize: Mathieson charge distribution in 2 dimensions
  63
  64     TObjArray* userObjects = static_cast<TObjArray*>(TVirtualFitter::GetFitter()->GetObjectFit());
  65
  66     AliMUONCluster* cluster = static_cast<AliMUONCluster*>(userObjects->At(0));
  67     AliMUONMathieson* mathieson = static_cast<AliMUONMathieson*>(userObjects->At(1));
  68
  69     f = 0.0;
  70     Float_t qTot = cluster->Charge();
  71 //    Float_t chargeCorrel[] = { cluster->Charge(0)/qTot, cluster->Charge(1)/qTot };
  72 //    Float_t qRatio[] = { 1.0/par[2], par[2] };
  73
  74     for ( Int_t i = 0 ; i < cluster->Multiplicity(); ++i )
  75     {
  76       AliMUONPad* pad = cluster->Pad(i);
  77       // skip pads w/ saturation or other problem(s)
  78       if ( pad->Status() ) continue;
  79       TVector2 lowerLeft = TVector2(par[0],par[1]) - pad->Position() - pad->Dimensions();
  80       TVector2 upperRight(lowerLeft + pad->Dimensions()*2.0);
  81       Float_t estimatedCharge = mathieson->IntXY(lowerLeft.X(),lowerLeft.Y(),
  82                                                  upperRight.X(),upperRight.Y());
  83 //      estimatedCharge *= 2/(1+qRatio[pad->Cathode()]);
  84       Float_t actualCharge = pad->Charge()/qTot;
  85
  86       Float_t delta = (estimatedCharge - actualCharge);
  87
  88       f += delta*delta;
  89     }
  90   }
  91 }
  92
  93 //_____________________________________________________________________________
  94 AliMUONClusterFinderSimpleFit::AliMUONClusterFinderSimpleFit()
  95 : AliMUONVClusterFinder(),
  96 fClusterFinder(0x0),
  97 fMathieson(0x0)
  98 {
  99   /// ctor
 100 }
 101
 102 //_____________________________________________________________________________
 103 AliMUONClusterFinderSimpleFit::~AliMUONClusterFinderSimpleFit()
 104 {
 105   /// dtor
 106   delete fClusterFinder;
 107   delete fMathieson;
 108 }
 109
 110 //_____________________________________________________________________________
 111 Bool_t
 112 AliMUONClusterFinderSimpleFit::Prepare(const AliMpVSegmentation* segmentations[2],
 113                                        const AliMUONVDigitStore& digitStore)
 114 {
 115   /// Prepare for clustering
 116
 117   // FIXME: should we get the Mathieson from elsewhere ?
 118
 119   // Find out the DetElemId
 120   Int_t detElemId(-1);
 121
 122   TIter next(digitStore.CreateIterator());
 123   AliMUONVDigit* d = static_cast<AliMUONVDigit*>(next());
 124
 125   if (d)
 126   {
 127     detElemId = d->DetElemId();
 128   }
 129   else
 130   {
 131     AliWarning("Could not find DE. Probably no digits at all : here's the digitStore :");
 132     StdoutToAliWarning(digitStore.Print(););
 133     return kFALSE;
 134   }
 135
 136   AliMp::StationType stationType = AliMpDEManager::GetStationType(detElemId);
 137
 138   Float_t kx3 = AliMUONConstants::SqrtKx3();
 139   Float_t ky3 = AliMUONConstants::SqrtKy3();
 140   Float_t pitch = AliMUONConstants::Pitch();
 141
 142   if ( stationType == AliMp::kStation1 )
 143   {
 144     kx3 = AliMUONConstants::SqrtKx3St1();
 145     ky3 = AliMUONConstants::SqrtKy3St1();
 146     pitch = AliMUONConstants::PitchSt1();
 147   }
 148
 149   delete fMathieson;
 150   fMathieson = new AliMUONMathieson;
 151
 152   fMathieson->SetPitch(pitch);
 153   fMathieson->SetSqrtKx3AndDeriveKx2Kx4(kx3);
 154   fMathieson->SetSqrtKy3AndDeriveKy2Ky4(ky3);
 155
 156   delete fClusterFinder;
 157   fClusterFinder = new AliMUONClusterFinderCOG;
 158   return fClusterFinder->Prepare(segmentations,digitStore);
 159 }
 160
 161 //_____________________________________________________________________________
 162 AliMUONCluster*
 163 AliMUONClusterFinderSimpleFit::NextCluster()
 164 {
 165   /// Returns next cluster
 166
 167   if ( !fClusterFinder ) return 0x0;
 168   AliMUONCluster* cluster = fClusterFinder->NextCluster();
 169   if ( cluster )
 170   {
 171     ComputePosition(*cluster);
 172
 173     if ( cluster->Charge() < 7 )
 174     {
 175       // skip that one
 176       return NextCluster();
 177     }
 178   }
 179   return cluster;
 180 }
 181
 182 //_____________________________________________________________________________
 183 void
 184 AliMUONClusterFinderSimpleFit::ComputePosition(AliMUONCluster& cluster)
 185 {
 186   /// Compute the position of the given cluster, by fitting a Mathieson
 187   /// charge distribution to it
 188
 189   TVirtualFitter* fitter = TVirtualFitter::Fitter(0,2);
 190   fitter->SetFCN(FitFunction);
 191
 192   if ( cluster.Multiplicity() < 3 ) return;
 193
 194   // We try a Mathieson fit, starting
 195   // with the center-of-gravity estimate as a first guess
 196   // for the cluster center.
 197
 198   Double_t xCOG = cluster.Position().X();
 199   Double_t yCOG = cluster.Position().Y();
 200
 201   Float_t stepX = 0.01; // cm
 202   Float_t stepY = 0.01; // cm
 203
 204   Double_t arg(-1); // disable printout
 205
 206   fitter->ExecuteCommand("SET PRINT",&arg,1);
 207
 208   fitter->SetParameter(0,"cluster X position",xCOG,stepX,0,0);
 209   fitter->SetParameter(1,"cluster Y position",yCOG,stepY,0,0);
 210
 211   TObjArray userObjects;
 212
 213   userObjects.Add(&cluster);
 214   userObjects.Add(fMathieson);
 215
 216   fitter->SetObjectFit(&userObjects);
 217
 218   Int_t val = fitter->ExecuteCommand("MIGRAD",0,0);
 219   AliDebug(1,Form("ExecuteCommand returned value=%d",val));
 220   if ( val )
 221   {
 222     // fit failed. Using COG results, with big errors
 223     AliWarning("Fit failed. Using COG results for cluster=");
 224     StdoutToAliWarning(cluster.Print());
 225     cluster.SetPosition(TVector2(xCOG,yCOG),TVector2(TMath::Abs(xCOG),TMath::Abs(yCOG)));
 226     cluster.SetChi2(1E3);
 227   }
 228
 229   Double_t results[] = { fitter->GetParameter(0),
 230     fitter->GetParameter(1) };
 231
 232   Double_t errors[] = { fitter->GetParError(0),
 233     fitter->GetParError(1) };
 234
 235   cluster.SetPosition(TVector2(results[0],results[1]),
 236                       TVector2(errors[0],errors[1]));
 237
 238   Double_t amin, edm, errdef;
 239   Int_t nvpar, nparx;
 240
 241   fitter->GetStats(amin, edm, errdef, nvpar, nparx);
 242
 243   Double_t chi2 = amin;
 244
 245   AliDebug(1,Form("Cluster fitted to (x,y)=(%e,%e) (xerr,yerr)=(%e,%e) \n chi2=%e ndf=%d",
 246                   results[0],results[1],
 247                   errors[0],errors[1],chi2,fitter->GetNumberFreeParameters()));
 248   cluster.SetChi2(chi2);
 249 }
 250
 251
 252