Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[u/mrichter/AliRoot.git] / MUON / AliMUONClusterFinderSimpleFit.cxx
CommitLineData
e036a7af 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"
e036a7af 22#include "AliMUONCluster.h"
23#include "AliMUONConstants.h"
3887d11e 24#include "AliMUONVDigit.h"
e036a7af 25#include "AliMUONMathieson.h"
26#include "AliMUONPad.h"
e036a7af 27#include "AliMpArea.h"
e036a7af 28#include "TObjArray.h"
29#include "TVector2.h"
30#include "TVirtualFitter.h"
31#include "TF1.h"
3887d11e 32#include "AliMUONVDigitStore.h"
33#include <Riostream.h>
e036a7af 34
3d1463c8 35//-----------------------------------------------------------------------------
e036a7af 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
3d1463c8 47//-----------------------------------------------------------------------------
e036a7af 48
78649106 49/// \cond CLASSIMP
e036a7af 50ClassImp(AliMUONClusterFinderSimpleFit)
51/// \endcond
52
53namespace
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();
32074b02 70// Float_t chargeCorrel[] = { cluster->Charge(0)/qTot, cluster->Charge(1)/qTot };
71// Float_t qRatio[] = { 1.0/par[2], par[2] };
e036a7af 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;
ff331318 78 TVector2 lowerLeft = TVector2(par[0],par[1]) - pad->Position() - pad->Dimensions();
e036a7af 79 TVector2 upperRight(lowerLeft + pad->Dimensions()*2.0);
32074b02 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;
e036a7af 84
85 Float_t delta = (estimatedCharge - actualCharge);
86
32074b02 87 f += delta*delta;
e036a7af 88 }
89 }
90}
91
92//_____________________________________________________________________________
b1a19e07 93AliMUONClusterFinderSimpleFit::AliMUONClusterFinderSimpleFit(AliMUONVClusterFinder* clusterFinder)
e036a7af 94: AliMUONVClusterFinder(),
b1a19e07 95fClusterFinder(clusterFinder),
110edb51 96fMathieson(0x0),
97fLowestClusterCharge(0)
e036a7af 98{
99 /// ctor
100}
101
102//_____________________________________________________________________________
103AliMUONClusterFinderSimpleFit::~AliMUONClusterFinderSimpleFit()
104{
105 /// dtor
106 delete fClusterFinder;
107 delete fMathieson;
108}
109
110//_____________________________________________________________________________
111Bool_t
24935e58 112AliMUONClusterFinderSimpleFit::Prepare(Int_t detElemId,
9e41a340 113 TObjArray* pads[2],
24935e58 114 const AliMpArea& area)
e036a7af 115{
116 /// Prepare for clustering
117
118 // FIXME: should we get the Mathieson from elsewhere ?
119
120 // Find out the DetElemId
4e51cfd2 121 AliMq::Station12Type stationType = AliMpDEManager::GetStation12Type(detElemId);
e036a7af 122
123 Float_t kx3 = AliMUONConstants::SqrtKx3();
124 Float_t ky3 = AliMUONConstants::SqrtKy3();
125 Float_t pitch = AliMUONConstants::Pitch();
126
4e51cfd2 127 if ( stationType == AliMq::kStation1 )
e036a7af 128 {
129 kx3 = AliMUONConstants::SqrtKx3St1();
130 ky3 = AliMUONConstants::SqrtKy3St1();
131 pitch = AliMUONConstants::PitchSt1();
132 }
133
134 delete fMathieson;
135 fMathieson = new AliMUONMathieson;
136
137 fMathieson->SetPitch(pitch);
138 fMathieson->SetSqrtKx3AndDeriveKx2Kx4(kx3);
139 fMathieson->SetSqrtKy3AndDeriveKy2Ky4(ky3);
140
24935e58 141 return fClusterFinder->Prepare(detElemId,pads,area);
e036a7af 142}
143
144//_____________________________________________________________________________
145AliMUONCluster*
146AliMUONClusterFinderSimpleFit::NextCluster()
147{
148 /// Returns next cluster
149
150 if ( !fClusterFinder ) return 0x0;
151 AliMUONCluster* cluster = fClusterFinder->NextCluster();
152 if ( cluster )
153 {
154 ComputePosition(*cluster);
155
110edb51 156 if ( cluster->Charge() < fLowestClusterCharge )
e036a7af 157 {
158 // skip that one
159 return NextCluster();
160 }
161 }
162 return cluster;
163}
164
165//_____________________________________________________________________________
166void
167AliMUONClusterFinderSimpleFit::ComputePosition(AliMUONCluster& cluster)
168{
169 /// Compute the position of the given cluster, by fitting a Mathieson
170 /// charge distribution to it
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
32074b02 187 Double_t arg(-1); // disable printout
e036a7af 188
32074b02 189 fitter->ExecuteCommand("SET PRINT",&arg,1);
e036a7af 190
191 fitter->SetParameter(0,"cluster X position",xCOG,stepX,0,0);
192 fitter->SetParameter(1,"cluster Y position",yCOG,stepY,0,0);
e036a7af 193
194 TObjArray userObjects;
195
196 userObjects.Add(&cluster);
197 userObjects.Add(fMathieson);
198
e036a7af 199 fitter->SetObjectFit(&userObjects);
200
32074b02 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 }
e036a7af 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
32074b02 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;
e036a7af 227
32074b02 228 AliDebug(1,Form("Cluster fitted to (x,y)=(%e,%e) (xerr,yerr)=(%e,%e) \n chi2=%e ndf=%d",
e036a7af 229 results[0],results[1],
32074b02 230 errors[0],errors[1],chi2,fitter->GetNumberFreeParameters()));
231 cluster.SetChi2(chi2);
e036a7af 232}
233
234
235