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