1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 #include "AliMUONClusterFinderSimpleFit.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"
33 #include "TVirtualFitter.h"
36 /// \class AliMUONClusterFinderSimpleFit
38 /// Basic cluster finder
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.
44 /// FIXME: this one is still at the developping stage...
46 /// \author Laurent Aphecetche
49 ClassImp(AliMUONClusterFinderSimpleFit)
54 //___________________________________________________________________________
56 FitFunction(Int_t& /*notused*/, Double_t* /*notused*/,
57 Double_t& f, Double_t* par,
60 /// Chi2 Function to minimize: Mathieson charge distribution in 2 dimensions
62 TObjArray* userObjects = static_cast<TObjArray*>(TVirtualFitter::GetFitter()->GetObjectFit());
64 AliMUONCluster* cluster = static_cast<AliMUONCluster*>(userObjects->At(0));
65 AliMUONMathieson* mathieson = static_cast<AliMUONMathieson*>(userObjects->At(1));
68 Float_t qTot = cluster->Charge();
69 Float_t chargeCorrel[] = { cluster->Charge(0)/qTot, cluster->Charge(1)/qTot };
71 for ( Int_t i = 0 ; i < cluster->Multiplicity(); ++i )
73 AliMUONPad* pad = cluster->Pad(i);
74 // skip pads w/ saturation or other problem(s)
75 if ( pad->Status() ) continue;
76 TVector2 lowerLeft = TVector2(par[0],par[1]) - pad->Position() - pad->Dimensions();
77 TVector2 upperRight(lowerLeft + pad->Dimensions()*2.0);
78 Float_t estimatedCharge =
79 qTot*mathieson->IntXY(lowerLeft.X(),lowerLeft.Y(),
80 upperRight.X(),upperRight.Y());
81 estimatedCharge *= chargeCorrel[pad->Cathode()];
82 Float_t actualCharge = pad->Charge();
84 Float_t delta = (estimatedCharge - actualCharge);
86 f += delta*delta/qTot;
91 //_____________________________________________________________________________
92 AliMUONClusterFinderSimpleFit::AliMUONClusterFinderSimpleFit()
93 : AliMUONVClusterFinder(),
100 //_____________________________________________________________________________
101 AliMUONClusterFinderSimpleFit::~AliMUONClusterFinderSimpleFit()
104 delete fClusterFinder;
108 //_____________________________________________________________________________
110 AliMUONClusterFinderSimpleFit::Prepare(const AliMpVSegmentation* segmentations[2],
111 TClonesArray* digits[2])
113 /// Prepare for clustering
115 // FIXME: should we get the Mathieson from elsewhere ?
117 // Find out the DetElemId
120 for ( Int_t i = 0; i < 2; ++i )
122 AliMUONDigit* d = static_cast<AliMUONDigit*>(digits[i]->First());
125 detElemId = d->DetElemId();
132 AliWarning("Could not find DE. Probably no digits at all ?");
136 AliMpStationType stationType = AliMpDEManager::GetStationType(detElemId);
138 Float_t kx3 = AliMUONConstants::SqrtKx3();
139 Float_t ky3 = AliMUONConstants::SqrtKy3();
140 Float_t pitch = AliMUONConstants::Pitch();
142 if ( stationType == kStation1 )
144 kx3 = AliMUONConstants::SqrtKx3St1();
145 ky3 = AliMUONConstants::SqrtKy3St1();
146 pitch = AliMUONConstants::PitchSt1();
150 fMathieson = new AliMUONMathieson;
152 fMathieson->SetPitch(pitch);
153 fMathieson->SetSqrtKx3AndDeriveKx2Kx4(kx3);
154 fMathieson->SetSqrtKy3AndDeriveKy2Ky4(ky3);
156 delete fClusterFinder;
157 fClusterFinder = new AliMUONClusterFinderCOG;
158 return fClusterFinder->Prepare(segmentations,digits);
161 //_____________________________________________________________________________
163 AliMUONClusterFinderSimpleFit::NextCluster()
165 /// Returns next cluster
167 if ( !fClusterFinder ) return 0x0;
168 AliMUONCluster* cluster = fClusterFinder->NextCluster();
171 ComputePosition(*cluster);
173 if ( cluster->Charge() < 7 )
176 return NextCluster();
182 //_____________________________________________________________________________
184 AliMUONClusterFinderSimpleFit::ComputePosition(AliMUONCluster& cluster)
186 /// Compute the position of the given cluster, by fitting a Mathieson
187 /// charge distribution to it
189 TVirtualFitter* fitter = TVirtualFitter::Fitter(0,2);
190 fitter->SetFCN(FitFunction);
192 if ( cluster.Multiplicity() < 3 ) return;
194 // We try a Mathieson fit, starting
195 // with the center-of-gravity estimate as a first guess
196 // for the cluster center.
198 Double_t xCOG = cluster.Position().X();
199 Double_t yCOG = cluster.Position().Y();
201 Float_t stepX = 0.01; // cm
202 Float_t stepY = 0.01; // cm
204 Double_t arg(1);//0);
206 fitter->ExecuteCommand("SET PRINT",&arg,1); // disable printout
208 fitter->SetParameter(0,"cluster X position",xCOG,stepX,0,0);
209 fitter->SetParameter(1,"cluster Y position",yCOG,stepY,0,0);
210 // fitter->SetParameter(2,"charge ratio",1.0,0.01,0.1,10);
212 TObjArray userObjects;
214 userObjects.Add(&cluster);
215 userObjects.Add(fMathieson);
217 // fitter->SetUserFunc(&userObjects);
218 fitter->SetObjectFit(&userObjects);
220 fitter->ExecuteCommand("MIGRAD",0,0);
222 Double_t results[] = { fitter->GetParameter(0),
223 fitter->GetParameter(1) };
225 Double_t errors[] = { fitter->GetParError(0),
226 fitter->GetParError(1) };
228 cluster.SetPosition(TVector2(results[0],results[1]),
229 TVector2(errors[0],errors[1]));
231 TF1* func = static_cast<TF1*>(fitter->GetUserFunc());
233 if ( func ) chi2 = func->GetChisquare();
235 AliDebug(1,Form("Cluster fitted to (x,y)=(%e,%e) (xerr,yerr)=(%e,%e) chi2=%e",
236 results[0],results[1],
237 errors[0],errors[1],chi2));
238 // cluster.SetChi2(chi2/fitter->GetNumberFreeParameters());