- AliMUONRecoParam.cxx:
[u/mrichter/AliRoot.git] / MUON / AliMUONRefitter.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 <cstdlib>
19 #include "AliMUONRefitter.h"
20 #include "AliMUONGeometryTransformer.h"
21 #include "AliMUONClusterFinderCOG.h"
22 #include "AliMUONClusterFinderMLEM.h"
23 #include "AliMUONClusterFinderSimpleFit.h"
24 #include "AliMUONPreClusterFinder.h"
25 #include "AliMUONPreClusterFinderV2.h"
26 #include "AliMUONPreClusterFinderV3.h"
27 #include "AliMUONSimpleClusterServer.h"
28 #include "AliMUONTrackReconstructor.h"
29 #include "AliMUONTrackReconstructorK.h"
30 #include "AliMUONRecoParam.h"
31 #include "AliMUONESDInterface.h"
32 #include "AliMUONVClusterStore.h"
33 #include "AliMUONVTrackStore.h"
34 #include "AliMUONTrack.h"
35 #include "AliMUONTracker.h"
36 #include "AliLog.h"
37
38 //-----------------------------------------------------------------------------
39 /// \class AliMUONRefitter
40 ///
41 /// This class has been developped to simplify the re-reconstruction of the MUON tracks
42 /// stored into ESD with different recoParams and/or after having re-calibrated the digits.
43 /// It creates new MUON object from ESD objects given as input (through the ESDInterface) then:
44 ///
45 /// - re-clusterize the ESD clusters using the attached ESD pads
46 ///   (several new clusters can be reconstructed per ESD cluster)
47 /// - re-fit the ESD tracks using the attached ESD clusters
48 /// - reconstruct the ESD tracks from ESD pads (i.e. re-clusterized the attached clusters)
49 ///
50 /// note:
51 /// - connexion between an ESD cluster and corresponding MUON clusters from re-clustering
52 ///   can be made through the detection element ID
53 /// - connexion between an ESD track and the corresponding refitted MUON track
54 ///   can be made through their unique ID
55 ///
56 /// \author Philippe Pillot
57 //-----------------------------------------------------------------------------
58
59 /// \cond CLASSIMP
60 ClassImp(AliMUONRefitter)
61 /// \endcond
62
63 //_____________________________________________________________________________
64 AliMUONRefitter::AliMUONRefitter(const AliMUONRecoParam* recoParam)
65 : TObject(),
66   fRecoParam(recoParam),
67   fGeometryTransformer(0x0),
68   fClusterServer(0x0),
69   fTracker(0x0),
70   fESDInterface(0x0)
71 {
72   /// Default constructor
73   CreateGeometryTransformer();
74   CreateClusterServer(*fGeometryTransformer);
75   if (fClusterServer) fTracker = AliMUONTracker::CreateTrackReconstructor(recoParam,fClusterServer);
76   if (!fClusterServer || !fTracker) {
77     AliFatal("refitter initialization failed");
78     exit(-1);
79   }
80 }
81
82 //_____________________________________________________________________________
83 AliMUONRefitter::~AliMUONRefitter()
84 {
85   /// Destructor
86   delete fGeometryTransformer;
87   delete fClusterServer;
88   delete fTracker;
89 }
90
91 //_____________________________________________________________________________
92 AliMUONVTrackStore* AliMUONRefitter::ReconstructFromDigits()
93 {
94   /// re-reconstruct all tracks and attached clusters from the digits
95   /// it is the responsability of the user to delete the returned store
96   
97   if (!fESDInterface) {
98     AliError("the refitter must be connected to an ESDInterface containing the ESD event to reconstruct");
99     return 0x0;
100   }
101   
102   // prepare new track(s)
103   AliMUONVTrackStore* newTrackStore = AliMUONESDInterface::NewTrackStore();
104   if (!newTrackStore) return 0x0;
105   
106   // loop over tracks and refit them (create new tracks)
107   AliMUONTrack *track;
108   TIter next(fESDInterface->CreateTrackIterator());
109   while ((track = static_cast<AliMUONTrack*>(next()))) {
110     AliMUONTrack *newTrack = RetrackFromDigits(*track);
111     newTrackStore->Add(newTrack);
112     delete newTrack;
113   }
114   
115   return newTrackStore;
116 }
117
118 //_____________________________________________________________________________
119 AliMUONVTrackStore* AliMUONRefitter::ReconstructFromClusters()
120 {
121   /// refit all tracks from the attached clusters
122   /// it is the responsability of the user to delete the returned store
123   
124   if (!fESDInterface) {
125     AliError("the refitter must be connected to an ESDInterface containing the ESD event to reconstruct");
126     return 0x0;
127   }
128   
129   // prepare new track(s)
130   AliMUONVTrackStore* newTrackStore = AliMUONESDInterface::NewTrackStore();
131   if (!newTrackStore) return 0x0;
132   
133   // loop over tracks and refit them (create new tracks)
134   AliMUONTrack *track;
135   TIter next(fESDInterface->CreateTrackIterator());
136   while ((track = static_cast<AliMUONTrack*>(next()))) {
137     AliMUONTrack* newTrack = newTrackStore->Add(*track);
138     if (!fTracker->RefitTrack(*newTrack)) newTrackStore->Remove(*newTrack);
139   }
140   
141   return newTrackStore;
142 }
143
144 //_____________________________________________________________________________
145 AliMUONTrack* AliMUONRefitter::RetrackFromDigits(UInt_t trackId)
146 {
147   /// refit track "trackId" from the digits (i.e. re-clusterized the attached clusters)
148   /// it is the responsability of the user to delete the returned track
149   
150   if (!fESDInterface) {
151     AliError("the refitter must be connected to an ESDInterface containing the ESD event to reconstruct");
152     return 0x0;
153   }
154   
155   // get the track to refit
156   AliMUONTrack* track = fESDInterface->FindTrack(trackId);
157   
158   return track ? RetrackFromDigits(*track) : 0x0;
159 }
160
161 //_____________________________________________________________________________
162 AliMUONTrack* AliMUONRefitter::RetrackFromClusters(UInt_t trackId)
163 {
164   /// refit track "trackId" form the clusters (i.e. do not re-clusterize)
165   /// it is the responsability of the user to delete the returned track
166   
167   if (!fESDInterface) {
168     AliError("the refitter must be connected to an ESDInterface containing the ESD event to reconstruct");
169     return 0x0;
170   }
171   
172   // get the track to refit
173   AliMUONTrack* track = fESDInterface->FindTrack(trackId);
174   if (!track) return 0x0;
175   
176   // refit the track (create a new one)
177   AliMUONTrack* newTrack = new AliMUONTrack(*track);
178   if (!fTracker->RefitTrack(*newTrack)) {
179     delete newTrack;
180     return 0x0;
181   }
182   
183   return newTrack;
184 }
185
186 //_____________________________________________________________________________
187 AliMUONVClusterStore* AliMUONRefitter::ReClusterize(UInt_t trackId, UInt_t clusterId)
188 {
189   /// re-clusterize cluster numbered "clusterId" in track "trackId"
190   /// several new clusters may be reconstructed
191   /// it is the responsability of the user to delete the returned store
192   
193   if (!fESDInterface) {
194     AliError("the refitter must be connected to an ESDInterface containing the ESD event to reconstruct");
195     return 0x0;
196   }
197   
198   // get the cluster to re-clusterize
199   AliMUONVCluster* cluster = fESDInterface->FindCluster(trackId,clusterId);
200   if (!cluster) return 0x0;
201   
202   // check if digits exist
203   if (cluster->GetNDigits() == 0) {
204     AliError(Form("no digit attached to cluster #%d in track %d",clusterId,trackId));
205     return 0x0;
206   }
207   
208   // create the cluster store
209   AliMUONVClusterStore* clusterStore = AliMUONESDInterface::NewClusterStore();
210   if (!clusterStore) return 0x0;
211   
212   // re-clusterize
213   TIter next(fESDInterface->CreateDigitIterator(trackId, clusterId));
214   fClusterServer->UseDigits(next);
215   fClusterServer->Clusterize(cluster->GetChamberId(),*clusterStore,AliMpArea());
216   
217   return clusterStore;
218 }
219
220 //_____________________________________________________________________________
221 AliMUONVClusterStore* AliMUONRefitter::ReClusterize(UInt_t clusterId)
222 {
223   /// re-clusterize cluster "clusterId"
224   /// several new clusters may be reconstructed
225   /// it is the responsability of the user to delete the returned store
226   
227   if (!fESDInterface) {
228     AliError("the refitter must be connected to an ESDInterface containing the ESD event to reconstruct");
229     return 0x0;
230   }
231   
232   // get the cluster to re-clusterize
233   AliMUONVCluster* cluster = fESDInterface->FindCluster(clusterId);
234   if (!cluster) return 0x0;
235   
236   // check if digits exist
237   if (cluster->GetNDigits() == 0) {
238     AliError(Form("no digit attached to cluster %d",clusterId));
239     return 0x0;
240   }
241   
242   // create the cluster store
243   AliMUONVClusterStore* clusterStore = AliMUONESDInterface::NewClusterStore();
244   if (!clusterStore) return 0x0;
245   
246   // re-clusterize
247   TIter next(fESDInterface->CreateDigitIteratorInCluster(clusterId));
248   fClusterServer->UseDigits(next);
249   fClusterServer->Clusterize(cluster->GetChamberId(),*clusterStore,AliMpArea());
250   
251   return clusterStore;
252 }
253
254 //_____________________________________________________________________________
255 void AliMUONRefitter::CreateGeometryTransformer()
256 {
257   /// Create geometry transformer (local<->global)
258   /// and load geometry data
259   fGeometryTransformer = new AliMUONGeometryTransformer();
260   fGeometryTransformer->LoadGeometryData();
261 }
262
263 //_____________________________________________________________________________
264 void AliMUONRefitter::CreateClusterServer(AliMUONGeometryTransformer& transformer)
265 {
266   /// Create cluster server
267   AliMUONVClusterFinder* clusterFinder = AliMUONReconstructor::CreateClusterFinder(GetRecoParam()->GetClusteringMode());
268   fClusterServer = clusterFinder ? new AliMUONSimpleClusterServer(clusterFinder,transformer) : 0x0;
269 }
270
271 //_____________________________________________________________________________
272 AliMUONTrack* AliMUONRefitter::RetrackFromDigits(const AliMUONTrack& track)
273 {
274   /// refit the given track from the digits (i.e. re-clusterized the attached clusters):
275   /// several new clusters may be reconstructed per initial ESD cluster:
276   /// -> all the combinations of clusters are considered to build the new tracks
277   /// -> return the best track (largest number of clusters or best chi2 in case of equality)
278   
279   // check if digits exist
280   UInt_t trackId = track.GetUniqueID();
281   if (fESDInterface->GetNDigits(trackId) == 0) {
282     AliError(Form("no digit attached to track #%d",trackId));
283     return 0x0;
284   }
285   
286   // prepare new track(s)
287   AliMUONVTrackStore* newTrackStore = AliMUONESDInterface::NewTrackStore();
288   if (!newTrackStore) return 0x0;
289   newTrackStore->Add(track)->Clear("C");
290   
291   // prepare new cluster store
292   AliMUONVClusterStore* newClusterStore = AliMUONESDInterface::NewClusterStore();
293   if (!newClusterStore) {
294     delete newTrackStore;
295     return 0x0;
296   }
297   
298   // loop over clusters, re-clusterize and build new tracks
299   AliMUONVCluster* cluster;
300   TIter nextCluster(fESDInterface->CreateClusterIterator(trackId));
301   while ((cluster = static_cast<AliMUONVCluster*>(nextCluster()))) {
302     
303     // reset the new cluster store
304     newClusterStore->Clear();
305     
306     // re-clusterize current cluster
307     TIter nextDigit(fESDInterface->CreateDigitIterator(trackId, cluster->GetUniqueID()));
308     fClusterServer->UseDigits(nextDigit);
309     Int_t nNewClusters = fClusterServer->Clusterize(cluster->GetChamberId(),*newClusterStore,AliMpArea());
310     
311     // check that re-clusterizing succeeded
312     if (nNewClusters == 0) {
313       AliWarning(Form("refit gave no cluster (chamber %d)",cluster->GetChamberId()));
314       AliInfo("initial ESD cluster:");
315       cluster->Print("FULL");
316       continue;
317     }
318     
319     // add the new cluster(s) to the tracks
320     AddClusterToTracks(*newClusterStore, *newTrackStore);
321     
322   }
323   
324   // refit the tracks and pick up the best one
325   AliMUONTrack *currentTrack, *bestTrack = 0x0;
326   Double_t currentChi2, bestChi2 = 1.e10;
327   Int_t currentNCluster, bestNClusters = 0;
328   TIter next(newTrackStore->CreateIterator());
329   while ((currentTrack = static_cast<AliMUONTrack*>(next()))) {
330     
331     // set the track parameters at first cluster if any (used as seed in original tracking)
332     AliMUONTrackParam* param = (AliMUONTrackParam*) currentTrack->GetTrackParamAtCluster()->First();
333     if (param) *param = *((AliMUONTrackParam*) track.GetTrackParamAtCluster()->First());
334     
335     // refit the track
336     if (!fTracker->RefitTrack(*currentTrack)) break;
337     
338     // find best track (the one with the higher number of cluster or the best chi2 in case of equality)
339     currentNCluster = currentTrack->GetNClusters();
340     currentChi2 = currentTrack->GetGlobalChi2();
341     if (currentNCluster > bestNClusters || (currentNCluster == bestNClusters && currentChi2 < bestChi2)) {
342       bestTrack = currentTrack;
343       bestNClusters = currentNCluster;
344       bestChi2 = currentChi2;
345     }
346     
347   }
348   
349   // copy best track and free memory
350   AliMUONTrack* newTrack = bestTrack ? new AliMUONTrack(*bestTrack) : 0x0;
351   delete newClusterStore;
352   delete newTrackStore;
353   
354   return newTrack;
355 }
356
357 //_____________________________________________________________________________
358 void AliMUONRefitter::AddClusterToTracks(const AliMUONVClusterStore &clusterStore, AliMUONVTrackStore &trackStore)
359 {
360   /// add clusters to each of the given tracks
361   /// duplicate the tracks if there are several clusters and add one cluster per copy
362   
363   // create new track store if there are more than 1 cluster to add per track
364   Int_t nClusters = clusterStore.GetSize();
365   if (nClusters < 1) return;
366   
367   AliMUONTrackParam dummyParam;
368   AliMUONTrack *currentTrack, *track;
369   AliMUONVCluster *newCluster;
370   Int_t nTracks = trackStore.GetSize();
371   Int_t iTrack = 0;
372   Int_t iCluster = 0;
373   
374   // loop over existing tracks to add the cluster(s)
375   TIter nextTrack(trackStore.CreateIterator());
376   while ((currentTrack = static_cast<AliMUONTrack*>(nextTrack())) && (iTrack < nTracks)) {
377     
378     iTrack++;
379     
380     // add the new cluster(s) to the tracks
381     // duplicate the tracks if there are several clusters
382     // the loop after loading the last cluster which is added to the current track
383     iCluster = 0;
384     TIter nextCluster(clusterStore.CreateIterator());
385     while ((newCluster = static_cast<AliMUONVCluster*>(nextCluster())) && (iCluster < nClusters - 1)) {
386       
387       iCluster++;
388       
389       // add a copy of the current track to the store
390       track = trackStore.Add(AliMUONTrack(*currentTrack));
391       
392       // only set Z parameter to avoid error in AddTrackParamAtCluster()
393       // the rest will be recomputed during refit
394       dummyParam.SetZ(newCluster->GetZ());
395       
396       // add new cluster to the new track
397       track->AddTrackParamAtCluster(dummyParam, *newCluster, kTRUE);
398       
399     }
400     
401     // only set Z parameter to avoid error in AddTrackParamAtCluster()
402     // the rest will be recomputed during refit
403     dummyParam.SetZ(newCluster->GetZ());
404     
405     // add new cluster to the current track
406     currentTrack->AddTrackParamAtCluster(dummyParam, *newCluster, kTRUE);
407     
408   }
409   
410 }
411