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(const AliMpVSegmentation* segmentations[2],
 114                                        const AliMUONVDigitStore& digitStore)
 115 {
 116   /// Prepare for clustering
 117
 118   // FIXME: should we get the Mathieson from elsewhere ?
 119
 120   // Find out the DetElemId
 121   Int_t detElemId(-1);
 122
 123   TIter next(digitStore.CreateIterator());
 124   AliMUONVDigit* d = static_cast<AliMUONVDigit*>(next());
 125
 126   if (d)
 127   {
 128     detElemId = d->DetElemId();
 129   }
 130   else
 131   {
 132     AliWarning("Could not find DE. Probably no digits at all : here's the digitStore :");
 133     StdoutToAliWarning(digitStore.Print(););
 134     return kFALSE;
 135   }
 136
 137   AliMp::StationType stationType = AliMpDEManager::GetStationType(detElemId);
 138
 139   Float_t kx3 = AliMUONConstants::SqrtKx3();
 140   Float_t ky3 = AliMUONConstants::SqrtKy3();
 141   Float_t pitch = AliMUONConstants::Pitch();
 142
 143   if ( stationType == AliMp::kStation1 )
 144   {
 145     kx3 = AliMUONConstants::SqrtKx3St1();
 146     ky3 = AliMUONConstants::SqrtKy3St1();
 147     pitch = AliMUONConstants::PitchSt1();
 148   }
 149
 150   delete fMathieson;
 151   fMathieson = new AliMUONMathieson;
 152
 153   fMathieson->SetPitch(pitch);
 154   fMathieson->SetSqrtKx3AndDeriveKx2Kx4(kx3);
 155   fMathieson->SetSqrtKy3AndDeriveKy2Ky4(ky3);
 156
 157   return fClusterFinder->Prepare(segmentations,digitStore);
 158 }
 159
 160 //_____________________________________________________________________________
 161 AliMUONCluster*
 162 AliMUONClusterFinderSimpleFit::NextCluster()
 163 {
 164   /// Returns next cluster
 165
 166   if ( !fClusterFinder ) return 0x0;
 167   AliMUONCluster* cluster = fClusterFinder->NextCluster();
 168   if ( cluster )
 169   {
 170     ComputePosition(*cluster);
 171
 172     if ( cluster->Charge() < 7 )
 173     {
 174       // skip that one
 175       return NextCluster();
 176     }
 177   }
 178   return cluster;
 179 }
 180
 181 //_____________________________________________________________________________
 182 void
 183 AliMUONClusterFinderSimpleFit::ComputePosition(AliMUONCluster& cluster)
 184 {
 185   /// Compute the position of the given cluster, by fitting a Mathieson
 186   /// charge distribution to it
 187
 188   TVirtualFitter* fitter = TVirtualFitter::Fitter(0,2);
 189   fitter->SetFCN(FitFunction);
 190
 191   if ( cluster.Multiplicity() < 3 ) return;
 192
 193   // We try a Mathieson fit, starting
 194   // with the center-of-gravity estimate as a first guess
 195   // for the cluster center.
 196
 197   Double_t xCOG = cluster.Position().X();
 198   Double_t yCOG = cluster.Position().Y();
 199
 200   Float_t stepX = 0.01; // cm
 201   Float_t stepY = 0.01; // cm
 202
 203   Double_t arg(-1); // disable printout
 204
 205   fitter->ExecuteCommand("SET PRINT",&arg,1);
 206
 207   fitter->SetParameter(0,"cluster X position",xCOG,stepX,0,0);
 208   fitter->SetParameter(1,"cluster Y position",yCOG,stepY,0,0);
 209
 210   TObjArray userObjects;
 211
 212   userObjects.Add(&cluster);
 213   userObjects.Add(fMathieson);
 214
 215   fitter->SetObjectFit(&userObjects);
 216
 217   Int_t val = fitter->ExecuteCommand("MIGRAD",0,0);
 218   AliDebug(1,Form("ExecuteCommand returned value=%d",val));
 219   if ( val )
 220   {
 221     // fit failed. Using COG results, with big errors
 222     AliWarning("Fit failed. Using COG results for cluster=");
 223     StdoutToAliWarning(cluster.Print());
 224     cluster.SetPosition(TVector2(xCOG,yCOG),TVector2(TMath::Abs(xCOG),TMath::Abs(yCOG)));
 225     cluster.SetChi2(1E3);
 226   }
 227
 228   Double_t results[] = { fitter->GetParameter(0),
 229     fitter->GetParameter(1) };
 230
 231   Double_t errors[] = { fitter->GetParError(0),
 232     fitter->GetParError(1) };
 233
 234   cluster.SetPosition(TVector2(results[0],results[1]),
 235                       TVector2(errors[0],errors[1]));
 236
 237   Double_t amin, edm, errdef;
 238   Int_t nvpar, nparx;
 239
 240   fitter->GetStats(amin, edm, errdef, nvpar, nparx);
 241
 242   Double_t chi2 = amin;
 243
 244   AliDebug(1,Form("Cluster fitted to (x,y)=(%e,%e) (xerr,yerr)=(%e,%e) \n chi2=%e ndf=%d",
 245                   results[0],results[1],
 246                   errors[0],errors[1],chi2,fitter->GetNumberFreeParameters()));
 247   cluster.SetChi2(chi2);
 248 }
 249
 250
 251