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 "AliMUONVDigit.h"
26 #include "AliMUONMathieson.h"
27 #include "AliMUONPad.h"
28 #include "AliMpArea.h"
29 #include "TClonesArray.h"
30 #include "TObjArray.h"
32 #include "TVirtualFitter.h"
34 #include "AliMUONVDigitStore.h"
35 #include <Riostream.h>
37 //-----------------------------------------------------------------------------
38 /// \class AliMUONClusterFinderSimpleFit
40 /// Basic cluster finder
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.
46 /// FIXME: this one is still at the developping stage...
48 /// \author Laurent Aphecetche
49 //-----------------------------------------------------------------------------
52 ClassImp(AliMUONClusterFinderSimpleFit)
57 //___________________________________________________________________________
59 FitFunction(Int_t& /*notused*/, Double_t* /*notused*/,
60 Double_t& f, Double_t* par,
63 /// Chi2 Function to minimize: Mathieson charge distribution in 2 dimensions
65 TObjArray* userObjects = static_cast<TObjArray*>(TVirtualFitter::GetFitter()->GetObjectFit());
67 AliMUONCluster* cluster = static_cast<AliMUONCluster*>(userObjects->At(0));
68 AliMUONMathieson* mathieson = static_cast<AliMUONMathieson*>(userObjects->At(1));
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] };
75 for ( Int_t i = 0 ; i < cluster->Multiplicity(); ++i )
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;
87 Float_t delta = (estimatedCharge - actualCharge);
94 //_____________________________________________________________________________
95 AliMUONClusterFinderSimpleFit::AliMUONClusterFinderSimpleFit(AliMUONVClusterFinder* clusterFinder)
96 : AliMUONVClusterFinder(),
97 fClusterFinder(clusterFinder),
103 //_____________________________________________________________________________
104 AliMUONClusterFinderSimpleFit::~AliMUONClusterFinderSimpleFit()
107 delete fClusterFinder;
111 //_____________________________________________________________________________
113 AliMUONClusterFinderSimpleFit::Prepare(const AliMpVSegmentation* segmentations[2],
114 const AliMUONVDigitStore& digitStore)
116 /// Prepare for clustering
118 // FIXME: should we get the Mathieson from elsewhere ?
120 // Find out the DetElemId
123 TIter next(digitStore.CreateIterator());
124 AliMUONVDigit* d = static_cast<AliMUONVDigit*>(next());
128 detElemId = d->DetElemId();
132 AliWarning("Could not find DE. Probably no digits at all : here's the digitStore :");
133 StdoutToAliWarning(digitStore.Print(););
137 AliMp::StationType stationType = AliMpDEManager::GetStationType(detElemId);
139 Float_t kx3 = AliMUONConstants::SqrtKx3();
140 Float_t ky3 = AliMUONConstants::SqrtKy3();
141 Float_t pitch = AliMUONConstants::Pitch();
143 if ( stationType == AliMp::kStation1 )
145 kx3 = AliMUONConstants::SqrtKx3St1();
146 ky3 = AliMUONConstants::SqrtKy3St1();
147 pitch = AliMUONConstants::PitchSt1();
151 fMathieson = new AliMUONMathieson;
153 fMathieson->SetPitch(pitch);
154 fMathieson->SetSqrtKx3AndDeriveKx2Kx4(kx3);
155 fMathieson->SetSqrtKy3AndDeriveKy2Ky4(ky3);
157 return fClusterFinder->Prepare(segmentations,digitStore);
160 //_____________________________________________________________________________
162 AliMUONClusterFinderSimpleFit::NextCluster()
164 /// Returns next cluster
166 if ( !fClusterFinder ) return 0x0;
167 AliMUONCluster* cluster = fClusterFinder->NextCluster();
170 ComputePosition(*cluster);
172 if ( cluster->Charge() < 7 )
175 return NextCluster();
181 //_____________________________________________________________________________
183 AliMUONClusterFinderSimpleFit::ComputePosition(AliMUONCluster& cluster)
185 /// Compute the position of the given cluster, by fitting a Mathieson
186 /// charge distribution to it
188 TVirtualFitter* fitter = TVirtualFitter::Fitter(0,2);
189 fitter->SetFCN(FitFunction);
191 if ( cluster.Multiplicity() < 3 ) return;
193 // We try a Mathieson fit, starting
194 // with the center-of-gravity estimate as a first guess
195 // for the cluster center.
197 Double_t xCOG = cluster.Position().X();
198 Double_t yCOG = cluster.Position().Y();
200 Float_t stepX = 0.01; // cm
201 Float_t stepY = 0.01; // cm
203 Double_t arg(-1); // disable printout
205 fitter->ExecuteCommand("SET PRINT",&arg,1);
207 fitter->SetParameter(0,"cluster X position",xCOG,stepX,0,0);
208 fitter->SetParameter(1,"cluster Y position",yCOG,stepY,0,0);
210 TObjArray userObjects;
212 userObjects.Add(&cluster);
213 userObjects.Add(fMathieson);
215 fitter->SetObjectFit(&userObjects);
217 Int_t val = fitter->ExecuteCommand("MIGRAD",0,0);
218 AliDebug(1,Form("ExecuteCommand returned value=%d",val));
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);
228 Double_t results[] = { fitter->GetParameter(0),
229 fitter->GetParameter(1) };
231 Double_t errors[] = { fitter->GetParError(0),
232 fitter->GetParError(1) };
234 cluster.SetPosition(TVector2(results[0],results[1]),
235 TVector2(errors[0],errors[1]));
237 Double_t amin, edm, errdef;
240 fitter->GetStats(amin, edm, errdef, nvpar, nparx);
242 Double_t chi2 = amin;
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);