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