Line data Source code
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 "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 : //-----------------------------------------------------------------------------
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 :
49 : /// \cond CLASSIMP
50 18 : 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 0 : TObjArray* userObjects = static_cast<TObjArray*>(TVirtualFitter::GetFitter()->GetObjectFit());
64 :
65 0 : AliMUONCluster* cluster = static_cast<AliMUONCluster*>(userObjects->At(0));
66 0 : AliMUONMathieson* mathieson = static_cast<AliMUONMathieson*>(userObjects->At(1));
67 :
68 0 : f = 0.0;
69 0 : 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 0 : for ( Int_t i = 0 ; i < cluster->Multiplicity(); ++i )
74 : {
75 0 : AliMUONPad* pad = cluster->Pad(i);
76 : // skip pads w/ saturation or other problem(s)
77 0 : if ( pad->Status() ) continue;
78 0 : TVector2 lowerLeft = TVector2(par[0],par[1]) - pad->Position() - pad->Dimensions();
79 0 : TVector2 upperRight(lowerLeft + pad->Dimensions()*2.0);
80 0 : Float_t estimatedCharge = mathieson->IntXY(lowerLeft.X(),lowerLeft.Y(),
81 0 : upperRight.X(),upperRight.Y());
82 : // estimatedCharge *= 2/(1+qRatio[pad->Cathode()]);
83 0 : Float_t actualCharge = pad->Charge()/qTot;
84 :
85 0 : Float_t delta = (estimatedCharge - actualCharge);
86 :
87 0 : f += delta*delta;
88 0 : }
89 0 : }
90 : }
91 :
92 : //_____________________________________________________________________________
93 : AliMUONClusterFinderSimpleFit::AliMUONClusterFinderSimpleFit(AliMUONVClusterFinder* clusterFinder)
94 0 : : AliMUONVClusterFinder(),
95 0 : fClusterFinder(clusterFinder),
96 0 : fMathieson(0x0),
97 0 : fLowestClusterCharge(0)
98 0 : {
99 : /// ctor
100 0 : }
101 :
102 : //_____________________________________________________________________________
103 : AliMUONClusterFinderSimpleFit::~AliMUONClusterFinderSimpleFit()
104 0 : {
105 : /// dtor
106 0 : delete fClusterFinder;
107 0 : delete fMathieson;
108 0 : }
109 :
110 : //_____________________________________________________________________________
111 : Bool_t
112 : AliMUONClusterFinderSimpleFit::Prepare(Int_t detElemId,
113 : TObjArray* pads[2],
114 : const AliMpArea& area)
115 : {
116 : /// Prepare for clustering
117 :
118 : // FIXME: should we get the Mathieson from elsewhere ?
119 :
120 : // Find out the DetElemId
121 0 : AliMq::Station12Type stationType = AliMpDEManager::GetStation12Type(detElemId);
122 :
123 0 : Float_t kx3 = AliMUONConstants::SqrtKx3();
124 0 : Float_t ky3 = AliMUONConstants::SqrtKy3();
125 0 : Float_t pitch = AliMUONConstants::Pitch();
126 :
127 0 : if ( stationType == AliMq::kStation1 )
128 : {
129 0 : kx3 = AliMUONConstants::SqrtKx3St1();
130 0 : ky3 = AliMUONConstants::SqrtKy3St1();
131 0 : pitch = AliMUONConstants::PitchSt1();
132 0 : }
133 :
134 0 : delete fMathieson;
135 0 : fMathieson = new AliMUONMathieson;
136 :
137 0 : fMathieson->SetPitch(pitch);
138 0 : fMathieson->SetSqrtKx3AndDeriveKx2Kx4(kx3);
139 0 : fMathieson->SetSqrtKy3AndDeriveKy2Ky4(ky3);
140 :
141 0 : return fClusterFinder->Prepare(detElemId,pads,area);
142 0 : }
143 :
144 : //_____________________________________________________________________________
145 : AliMUONCluster*
146 : AliMUONClusterFinderSimpleFit::NextCluster()
147 : {
148 : /// Returns next cluster
149 :
150 0 : if ( !fClusterFinder ) return 0x0;
151 0 : AliMUONCluster* cluster = fClusterFinder->NextCluster();
152 0 : if ( cluster )
153 : {
154 0 : ComputePosition(*cluster);
155 :
156 0 : if ( cluster->Charge() < fLowestClusterCharge )
157 : {
158 : // skip that one
159 0 : return NextCluster();
160 : }
161 : }
162 0 : return cluster;
163 0 : }
164 :
165 : //_____________________________________________________________________________
166 : void
167 : AliMUONClusterFinderSimpleFit::ComputePosition(AliMUONCluster& cluster)
168 : {
169 : /// Compute the position of the given cluster, by fitting a Mathieson
170 : /// charge distribution to it
171 :
172 0 : TVirtualFitter* fitter = TVirtualFitter::Fitter(0,2);
173 0 : fitter->SetFCN(FitFunction);
174 :
175 0 : 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 0 : Double_t xCOG = cluster.Position().X();
182 0 : Double_t yCOG = cluster.Position().Y();
183 :
184 : Float_t stepX = 0.01; // cm
185 : Float_t stepY = 0.01; // cm
186 :
187 0 : Double_t arg(-1); // disable printout
188 :
189 0 : fitter->ExecuteCommand("SET PRINT",&arg,1);
190 :
191 0 : fitter->SetParameter(0,"cluster X position",xCOG,stepX,0,0);
192 0 : fitter->SetParameter(1,"cluster Y position",yCOG,stepY,0,0);
193 :
194 0 : TObjArray userObjects;
195 :
196 0 : userObjects.Add(&cluster);
197 0 : userObjects.Add(fMathieson);
198 :
199 0 : fitter->SetObjectFit(&userObjects);
200 :
201 0 : Int_t val = fitter->ExecuteCommand("MIGRAD",0,0);
202 0 : AliDebug(1,Form("ExecuteCommand returned value=%d",val));
203 0 : if ( val )
204 : {
205 : // fit failed. Using COG results, with big errors
206 0 : AliWarning("Fit failed. Using COG results for cluster=");
207 0 : StdoutToAliWarning(cluster.Print());
208 0 : cluster.SetPosition(TVector2(xCOG,yCOG),TVector2(TMath::Abs(xCOG),TMath::Abs(yCOG)));
209 0 : cluster.SetChi2(1E3);
210 0 : }
211 :
212 0 : Double_t results[] = { fitter->GetParameter(0),
213 0 : fitter->GetParameter(1) };
214 :
215 0 : Double_t errors[] = { fitter->GetParError(0),
216 0 : fitter->GetParError(1) };
217 :
218 0 : cluster.SetPosition(TVector2(results[0],results[1]),
219 0 : TVector2(errors[0],errors[1]));
220 :
221 0 : Double_t amin, edm, errdef;
222 0 : Int_t nvpar, nparx;
223 :
224 0 : fitter->GetStats(amin, edm, errdef, nvpar, nparx);
225 :
226 0 : Double_t chi2 = amin;
227 :
228 0 : 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 0 : cluster.SetChi2(chi2);
232 0 : }
233 :
234 :
235 :
|
