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