AliRoot Core  edcc906 (edcc906)
AliMUONClusterFinderSimpleFit.cxx
Go to the documentation of this file.
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 
19 
20 #include "AliLog.h"
21 #include "AliMpDEManager.h"
22 #include "AliMUONCluster.h"
23 #include "AliMUONConstants.h"
24 #include "AliMUONVDigit.h"
25 #include "AliMUONMathieson.h"
26 #include "AliMUONPad.h"
27 #include "AliMpArea.h"
28 #include "TObjArray.h"
29 #include "TVector2.h"
30 #include "TVirtualFitter.h"
31 #include "TF1.h"
32 #include "AliMUONVDigitStore.h"
33 #include <Riostream.h>
34 
35 //-----------------------------------------------------------------------------
47 //-----------------------------------------------------------------------------
48 
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  {
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 //_____________________________________________________________________________
95 fClusterFinder(clusterFinder),
96 fMathieson(0x0),
97 fLowestClusterCharge(0)
98 {
100 }
101 
102 //_____________________________________________________________________________
104 {
106  delete fClusterFinder;
107  delete fMathieson;
108 }
109 
110 //_____________________________________________________________________________
111 Bool_t
113  TObjArray* pads[2],
114  const AliMpArea& area)
115 {
117 
118  // FIXME: should we get the Mathieson from elsewhere ?
119 
120  // Find out the DetElemId
122 
123  Float_t kx3 = AliMUONConstants::SqrtKx3();
124  Float_t ky3 = AliMUONConstants::SqrtKy3();
125  Float_t pitch = AliMUONConstants::Pitch();
126 
127  if ( stationType == AliMq::kStation1 )
128  {
131  pitch = AliMUONConstants::PitchSt1();
132  }
133 
134  delete fMathieson;
136 
137  fMathieson->SetPitch(pitch);
140 
141  return fClusterFinder->Prepare(detElemId,pads,area);
142 }
143 
144 //_____________________________________________________________________________
147 {
149 
150  if ( !fClusterFinder ) return 0x0;
152  if ( cluster )
153  {
154  ComputePosition(*cluster);
155 
156  if ( cluster->Charge() < fLowestClusterCharge )
157  {
158  // skip that one
159  return NextCluster();
160  }
161  }
162  return cluster;
163 }
164 
165 //_____________________________________________________________________________
166 void
168 {
171 
172  TVirtualFitter* fitter = TVirtualFitter::Fitter(0,2);
173  fitter->SetFCN(FitFunction);
174 
175  if ( cluster.Multiplicity() < 3 ) return;
176 
177  // We try a Mathieson fit, starting
178  // with the center-of-gravity estimate as a first guess
179  // for the cluster center.
180 
181  Double_t xCOG = cluster.Position().X();
182  Double_t yCOG = cluster.Position().Y();
183 
184  Float_t stepX = 0.01; // cm
185  Float_t stepY = 0.01; // cm
186 
187  Double_t arg(-1); // disable printout
188 
189  fitter->ExecuteCommand("SET PRINT",&arg,1);
190 
191  fitter->SetParameter(0,"cluster X position",xCOG,stepX,0,0);
192  fitter->SetParameter(1,"cluster Y position",yCOG,stepY,0,0);
193 
194  TObjArray userObjects;
195 
196  userObjects.Add(&cluster);
197  userObjects.Add(fMathieson);
198 
199  fitter->SetObjectFit(&userObjects);
200 
201  Int_t val = fitter->ExecuteCommand("MIGRAD",0,0);
202  AliDebug(1,Form("ExecuteCommand returned value=%d",val));
203  if ( val )
204  {
205  // fit failed. Using COG results, with big errors
206  AliWarning("Fit failed. Using COG results for cluster=");
207  StdoutToAliWarning(cluster.Print());
208  cluster.SetPosition(TVector2(xCOG,yCOG),TVector2(TMath::Abs(xCOG),TMath::Abs(yCOG)));
209  cluster.SetChi2(1E3);
210  }
211 
212  Double_t results[] = { fitter->GetParameter(0),
213  fitter->GetParameter(1) };
214 
215  Double_t errors[] = { fitter->GetParError(0),
216  fitter->GetParError(1) };
217 
218  cluster.SetPosition(TVector2(results[0],results[1]),
219  TVector2(errors[0],errors[1]));
220 
221  Double_t amin, edm, errdef;
222  Int_t nvpar, nparx;
223 
224  fitter->GetStats(amin, edm, errdef, nvpar, nparx);
225 
226  Double_t chi2 = amin;
227 
228  AliDebug(1,Form("Cluster fitted to (x,y)=(%e,%e) (xerr,yerr)=(%e,%e) \n chi2=%e ndf=%d",
229  results[0],results[1],
230  errors[0],errors[1],chi2,fitter->GetNumberFreeParameters()));
231  cluster.SetChi2(chi2);
232 }
233 
234 
235 
static AliMq::Station12Type GetStation12Type(Int_t detElemId)
Interface of a cluster finder.
void SetSqrtKx3AndDeriveKx2Kx4(Float_t SqrtKx3)
Mathieson sqrt{Kx3} and derived Kx2 and Kx4.
static Float_t Pitch()
Return wire pitch.
Float_t IntXY(Float_t xi1, Float_t yi1, Float_t xi2, Float_t yi2) const
Charge integration on region (x1,y1,x2,y2).
#define TObjArray
A group of adjacent pads.
static Float_t SqrtKy3St1()
Return SqrtKy3 for Station 1 & 2.
AliMUONClusterFinderSimpleFit(AliMUONVClusterFinder *clusterFinder)
static Float_t SqrtKx3()
Return SqrtKx3 for Slat.
Implementation of Mathieson response.
TVector2 Dimensions() const
Return half dimensions in x and y (cm)
Definition: AliMUONPad.h:56
virtual void Print(Option_t *opt="") const
AliMUONVClusterFinder * fClusterFinder
! the preclustering we use
A rectangle area positioned in plane..
Definition: AliMpArea.h:20
void ComputePosition(AliMUONCluster &cluster)
Int_t Multiplicity() const
virtual Bool_t Prepare(Int_t detElemId, TObjArray *pads[2], const AliMpArea &area)
void SetPosition(const TVector2 &pos, const TVector2 &errorOnPos)
Set (x,y) of that cluster and errors.
static Float_t PitchSt1()
Return wire pitch for Station 1 & 2.
#define AliWarning(message)
Definition: AliLog.h:541
static Float_t SqrtKy3()
Return SqrtKy3 for Slat.
static Float_t SqrtKx3St1()
Return SqrtKx3 for Station 1 & 2.
virtual Bool_t Prepare(Int_t detElemId, TObjArray *pads[2], const AliMpArea &area)
return clusterFinder
Double_t chi2
Definition: AnalyzeLaser.C:7
void SetPitch(Float_t p1)
AliMUONMathieson * fMathieson
! Mathieson to compute the charge repartition
#define StdoutToAliWarning(whatever)
Definition: AliLog.h:570
TVector2 Position() const
Return (x,y) of that cluster.
TF1 * f
Definition: interpolTest.C:21
void SetChi2(Float_t chi2)
Set chi2 of the RawCharge fit.
TVectorD errors
Definition: driftITSTPC.C:97
TVector2 Position() const
Return positions in x and y (cm)
Definition: AliMUONPad.h:85
Float_t Charge() const
#define AliDebug(logLevel, message)
Definition: AliLog.h:300
virtual AliMUONCluster * NextCluster()=0
Double_t fLowestClusterCharge
! minimum cluster charge we allow
Double_t Charge() const
Return pad charge.
Definition: AliMUONPad.h:48
Int_t Status() const
Return status word.
Definition: AliMUONPad.h:123
AliMUONPad * Pad(Int_t index) const
void SetSqrtKy3AndDeriveKy2Ky4(Float_t SqrtKy3)
Mathieson sqrt{Ky3} and derived Ky2 and Ky4.
Combination of digit and mppad informations.
Definition: AliMUONPad.h:25