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 "AliMUONDigit.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
  36 /// \class AliMUONClusterFinderSimpleFit
  37 ///
  38 /// Basic cluster finder
  39 ///
  40 /// We simply use AliMUONPreClusterFinder to get basic cluster,
  41 /// and then we try to fit the charge repartition using a Mathieson
  42 /// distribution, varying the position.
  43 ///
  44 /// FIXME: this one is still at the developping stage...
  45 ///
  46 /// \author Laurent Aphecetche
  47
  48 /// \cond CLASSIMP
  49 ClassImp(AliMUONClusterFinderSimpleFit)
  50 /// \endcond
  51
  52 namespace
  53 {
  54   //___________________________________________________________________________
  55   void
  56   FitFunction(Int_t& /*notused*/, Double_t* /*notused*/,
  57               Double_t& f, Double_t* par,
  58               Int_t /*notused*/)
  59   {
  60     /// Chi2 Function to minimize: Mathieson charge distribution in 2 dimensions
  61
  62     TObjArray* userObjects = static_cast<TObjArray*>(TVirtualFitter::GetFitter()->GetObjectFit());
  63
  64     AliMUONCluster* cluster = static_cast<AliMUONCluster*>(userObjects->At(0));
  65     AliMUONMathieson* mathieson = static_cast<AliMUONMathieson*>(userObjects->At(1));
  66
  67     f = 0.0;
  68     Float_t qTot = cluster->Charge();
  69 //    Float_t chargeCorrel[] = { cluster->Charge(0)/qTot, cluster->Charge(1)/qTot };
  70 //    Float_t qRatio[] = { 1.0/par[2], par[2] };
  71
  72     for ( Int_t i = 0 ; i < cluster->Multiplicity(); ++i )
  73     {
  74       AliMUONPad* pad = cluster->Pad(i);
  75       // skip pads w/ saturation or other problem(s)
  76       if ( pad->Status() ) continue;
  77       TVector2 lowerLeft = TVector2(par[0],par[1]) - pad->Position() - pad->Dimensions();
  78       TVector2 upperRight(lowerLeft + pad->Dimensions()*2.0);
  79       Float_t estimatedCharge = mathieson->IntXY(lowerLeft.X(),lowerLeft.Y(),
  80                                                  upperRight.X(),upperRight.Y());
  81 //      estimatedCharge *= 2/(1+qRatio[pad->Cathode()]);
  82       Float_t actualCharge = pad->Charge()/qTot;
  83
  84       Float_t delta = (estimatedCharge - actualCharge);
  85
  86       f += delta*delta;
  87     }
  88   }
  89 }
  90
  91 //_____________________________________________________________________________
  92 AliMUONClusterFinderSimpleFit::AliMUONClusterFinderSimpleFit()
  93 : AliMUONVClusterFinder(),
  94 fClusterFinder(0x0),
  95 fMathieson(0x0)
  96 {
  97   /// ctor
  98 }
  99
 100 //_____________________________________________________________________________
 101 AliMUONClusterFinderSimpleFit::~AliMUONClusterFinderSimpleFit()
 102 {
 103   /// dtor
 104   delete fClusterFinder;
 105   delete fMathieson;
 106 }
 107
 108 //_____________________________________________________________________________
 109 Bool_t
 110 AliMUONClusterFinderSimpleFit::Prepare(const AliMpVSegmentation* segmentations[2],
 111                                        TClonesArray* digits[2])
 112 {
 113   /// Prepare for clustering
 114
 115   // FIXME: should we get the Mathieson from elsewhere ?
 116
 117   // Find out the DetElemId
 118   Int_t detElemId(-1);
 119
 120   for ( Int_t i = 0; i < 2; ++i )
 121   {
 122     AliMUONDigit* d = static_cast<AliMUONDigit*>(digits[i]->First());
 123     if (d)
 124     {
 125       detElemId = d->DetElemId();
 126       break;
 127     }
 128   }
 129
 130   if ( detElemId < 0 )
 131   {
 132     AliWarning("Could not find DE. Probably no digits at all ?");
 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,digits);
 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