Write only detector coefficients from HLT (Raphaelle)
[u/mrichter/AliRoot.git] / MUON / AliMUONClusterFinderSimpleFit.cxx
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 "TClonesArray.h"
29 #include "TObjArray.h"
30 #include "TVector2.h"
31 #include "TVirtualFitter.h"
32 #include "TF1.h"
33 #include "AliMUONVDigitStore.h"
34 #include <Riostream.h>
35
36 //-----------------------------------------------------------------------------
37 /// \class AliMUONClusterFinderSimpleFit
38 ///
39 /// Basic cluster finder 
40 /// 
41 /// We simply use AliMUONPreClusterFinder to get basic cluster,
42 /// and then we try to fit the charge repartition using a Mathieson
43 /// distribution, varying the position.
44 ///
45 /// FIXME: this one is still at the developping stage...
46 ///
47 /// \author Laurent Aphecetche
48 //-----------------------------------------------------------------------------
49
50 /// \cond CLASSIMP
51 ClassImp(AliMUONClusterFinderSimpleFit)
52 /// \endcond
53
54 namespace
55 {
56   //___________________________________________________________________________
57   void 
58   FitFunction(Int_t& /*notused*/, Double_t* /*notused*/, 
59               Double_t& f, Double_t* par, 
60               Int_t /*notused*/)
61   {
62     /// Chi2 Function to minimize: Mathieson charge distribution in 2 dimensions
63     
64     TObjArray* userObjects = static_cast<TObjArray*>(TVirtualFitter::GetFitter()->GetObjectFit());
65     
66     AliMUONCluster* cluster = static_cast<AliMUONCluster*>(userObjects->At(0));
67     AliMUONMathieson* mathieson = static_cast<AliMUONMathieson*>(userObjects->At(1));
68     
69     f = 0.0;
70     Float_t qTot = cluster->Charge();
71 //    Float_t chargeCorrel[] = { cluster->Charge(0)/qTot, cluster->Charge(1)/qTot };
72 //    Float_t qRatio[] = { 1.0/par[2], par[2] };
73     
74     for ( Int_t i = 0 ; i < cluster->Multiplicity(); ++i )
75     {
76       AliMUONPad* pad = cluster->Pad(i);
77       // skip pads w/ saturation or other problem(s)
78       if ( pad->Status() ) continue; 
79       TVector2 lowerLeft = TVector2(par[0],par[1]) - pad->Position() - pad->Dimensions();
80       TVector2 upperRight(lowerLeft + pad->Dimensions()*2.0);
81       Float_t estimatedCharge = mathieson->IntXY(lowerLeft.X(),lowerLeft.Y(),
82                                                  upperRight.X(),upperRight.Y());
83 //      estimatedCharge *= 2/(1+qRatio[pad->Cathode()]);
84       Float_t actualCharge = pad->Charge()/qTot;
85       
86       Float_t delta = (estimatedCharge - actualCharge);
87       
88       f += delta*delta;    
89     }  
90   }
91 }
92
93 //_____________________________________________________________________________
94 AliMUONClusterFinderSimpleFit::AliMUONClusterFinderSimpleFit(AliMUONVClusterFinder* clusterFinder)
95 : AliMUONVClusterFinder(),
96 fClusterFinder(clusterFinder),
97 fMathieson(0x0),
98 fLowestClusterCharge(0)
99 {
100   /// ctor
101 }
102
103 //_____________________________________________________________________________
104 AliMUONClusterFinderSimpleFit::~AliMUONClusterFinderSimpleFit()
105 {
106   /// dtor
107   delete fClusterFinder;
108   delete fMathieson;
109 }
110
111 //_____________________________________________________________________________
112 Bool_t 
113 AliMUONClusterFinderSimpleFit::Prepare(Int_t detElemId,
114                                        TClonesArray* pads[2],
115                                        const AliMpArea& area)
116 {
117   /// Prepare for clustering
118
119   // FIXME: should we get the Mathieson from elsewhere ?
120   
121   // Find out the DetElemId
122   AliMq::Station12Type stationType = AliMpDEManager::GetStation12Type(detElemId);
123   
124   Float_t kx3 = AliMUONConstants::SqrtKx3();
125   Float_t ky3 = AliMUONConstants::SqrtKy3();
126   Float_t pitch = AliMUONConstants::Pitch();
127   
128   if ( stationType == AliMq::kStation1 )
129   {
130     kx3 = AliMUONConstants::SqrtKx3St1();
131     ky3 = AliMUONConstants::SqrtKy3St1();
132     pitch = AliMUONConstants::PitchSt1();
133   }
134   
135   delete fMathieson;
136   fMathieson = new AliMUONMathieson;
137   
138   fMathieson->SetPitch(pitch);
139   fMathieson->SetSqrtKx3AndDeriveKx2Kx4(kx3);
140   fMathieson->SetSqrtKy3AndDeriveKy2Ky4(ky3);
141
142   return fClusterFinder->Prepare(detElemId,pads,area);
143 }
144
145 //_____________________________________________________________________________
146 AliMUONCluster* 
147 AliMUONClusterFinderSimpleFit::NextCluster()
148 {
149   /// Returns next cluster
150   
151   if ( !fClusterFinder ) return 0x0;
152   AliMUONCluster* cluster = fClusterFinder->NextCluster();
153   if ( cluster )
154   {
155     ComputePosition(*cluster);
156
157     if ( cluster->Charge() < fLowestClusterCharge )
158     {
159       // skip that one
160       return NextCluster();
161     }    
162   }
163   return cluster;
164 }
165
166 //_____________________________________________________________________________
167 void 
168 AliMUONClusterFinderSimpleFit::ComputePosition(AliMUONCluster& cluster)
169 {
170   /// Compute the position of the given cluster, by fitting a Mathieson
171   /// charge distribution to it
172   
173   TVirtualFitter* fitter = TVirtualFitter::Fitter(0,2);
174   fitter->SetFCN(FitFunction);
175
176   if ( cluster.Multiplicity() < 3 ) return;
177   
178   // We try a Mathieson fit, starting
179   // with the center-of-gravity estimate as a first guess
180   // for the cluster center.
181   
182   Double_t xCOG = cluster.Position().X();
183   Double_t yCOG = cluster.Position().Y();
184   
185   Float_t stepX = 0.01; // cm
186   Float_t stepY = 0.01; // cm
187   
188   Double_t arg(-1); // disable printout
189   
190   fitter->ExecuteCommand("SET PRINT",&arg,1);
191   
192   fitter->SetParameter(0,"cluster X position",xCOG,stepX,0,0);
193   fitter->SetParameter(1,"cluster Y position",yCOG,stepY,0,0);
194   
195   TObjArray userObjects;
196   
197   userObjects.Add(&cluster);
198   userObjects.Add(fMathieson);
199   
200   fitter->SetObjectFit(&userObjects);
201   
202   Int_t val = fitter->ExecuteCommand("MIGRAD",0,0);
203   AliDebug(1,Form("ExecuteCommand returned value=%d",val));
204   if ( val ) 
205   {
206     // fit failed. Using COG results, with big errors
207     AliWarning("Fit failed. Using COG results for cluster=");
208     StdoutToAliWarning(cluster.Print());
209     cluster.SetPosition(TVector2(xCOG,yCOG),TVector2(TMath::Abs(xCOG),TMath::Abs(yCOG)));
210     cluster.SetChi2(1E3);
211   }
212   
213   Double_t results[] = { fitter->GetParameter(0),
214     fitter->GetParameter(1) };
215   
216   Double_t errors[] = { fitter->GetParError(0),
217     fitter->GetParError(1) };
218   
219   cluster.SetPosition(TVector2(results[0],results[1]),
220                       TVector2(errors[0],errors[1]));
221   
222   Double_t amin, edm, errdef;
223   Int_t nvpar, nparx;
224   
225   fitter->GetStats(amin, edm, errdef, nvpar, nparx);
226
227   Double_t chi2 = amin;
228   
229   AliDebug(1,Form("Cluster fitted to (x,y)=(%e,%e) (xerr,yerr)=(%e,%e) \n chi2=%e ndf=%d",
230                   results[0],results[1],
231                   errors[0],errors[1],chi2,fitter->GetNumberFreeParameters()));
232   cluster.SetChi2(chi2);
233 }
234
235
236