- AliMUONRecoParam.cxx:
[u/mrichter/AliRoot.git] / MUON / MUONRefit.C
CommitLineData
7deb8eb0 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/// \ingroup macros
19/// \file MUONRefit.C
20/// \brief Macro for refitting ESD tracks from ESD pads
21///
22/// \author Philippe Pillot, SUBATECH
23
24#if !defined(__CINT__) || defined(__MAKECINT__)
25#include <TStopwatch.h>
26#include <TFile.h>
27#include <TTree.h>
28#include <TString.h>
29#include <Riostream.h>
30#include <TGeoManager.h>
31#include <TRandom.h>
32#include <TROOT.h>
33
34// STEER includes
35#include "AliMagFMaps.h"
36#include "AliTracker.h"
37#include "AliESDEvent.h"
38#include "AliESDMuonTrack.h"
39#include "AliRecoParam.h"
40#include "AliCDBManager.h"
41#include "AliGeomManager.h"
42
43// MUON includes
44#include "AliMpCDB.h"
45#include "AliMUONRecoParam.h"
46#include "AliMUONESDInterface.h"
47#include "AliMUONRefitter.h"
48#include "AliMUONVDigit.h"
49#include "AliMUONVCluster.h"
50#include "AliMUONVClusterStore.h"
51#include "AliMUONTrack.h"
52#include "AliMUONVTrackStore.h"
53#include "AliMUONTrackParam.h"
54#include "AliMUONTrackExtrap.h"
55#endif
56
57const Int_t printLevel = 1;
58
59void Prepare();
60TTree* GetESDTree(TFile *esdFile);
61
62//-----------------------------------------------------------------------
63void MUONRefit(Int_t nevents = -1, const char* esdFileNameIn = "AliESDs.root", const char* esdFileNameOut = "AliESDs_New.root")
64{
65 /// refit ESD tracks from ESD pads (i.e. re-clusterized the attached ESD clusters);
66 /// reset the charge of the digit using their raw charge before refitting;
67 /// compare results with original ESD tracks;
68 /// write results in a new ESD file
69
70 // prepare the refitting
71 gRandom->SetSeed(1);
72 Prepare();
ad3c6eda 73
74 // set reconstruction parameters used for refitting
75 AliMUONRecoParam *muonRecoParam = AliMUONRecoParam::GetLowFluxParam();
76 muonRecoParam->Print("FULL");
77
7deb8eb0 78 AliMUONESDInterface esdInterface;
ad3c6eda 79 AliMUONRefitter refitter(muonRecoParam);
7deb8eb0 80 refitter.Connect(&esdInterface);
81
82 // open the ESD file and tree
83 TFile* esdFile = TFile::Open(esdFileNameIn);
84 TTree* esdTree = GetESDTree(esdFile);
85
ad3c6eda 86 // create the ESD output file and tree
87 TFile* newESDFile = TFile::Open(esdFileNameOut, "RECREATE");
88 newESDFile->SetCompressionLevel(2);
7deb8eb0 89 TTree* newESDTree = esdTree->CloneTree(0);
90
91 // connect ESD event to the ESD tree
92 AliESDEvent* esd = new AliESDEvent();
93 esd->ReadFromTree(esdTree);
94
95 // timer start...
96 TStopwatch timer;
97
98 // Loop over ESD events
99 if (nevents > 0) nevents = TMath::Min(nevents,(Int_t)esdTree->GetEntries());
100 else nevents = (Int_t)esdTree->GetEntries();
101 for (Int_t iEvent = 0; iEvent < nevents; iEvent++) {
102 if (printLevel>0) cout<<endl<<" ****************event #"<<iEvent+1<<"****************"<<endl;
103
104 //----------------------------------------------//
105 // -------------- process event --------------- //
106 //----------------------------------------------//
107 // get the ESD of current event
108 esdTree->GetEvent(iEvent);
109 if (!esd) {
110 Error("CheckESD", "no ESD object found for event %d", iEvent);
111 return;
112 }
113 Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks();
114 if (nTracks < 1) continue;
115
116 // load the current event
117 esdInterface.LoadEvent(*esd);
118
119 // loop over digit to modify their charge
120 AliMUONVDigit *digit;
121 TIter next(esdInterface.CreateDigitIterator());
b1fea02e 122 while ((digit = static_cast<AliMUONVDigit*>(next()))) {
123 digit->SetCharge(digit->ADC());
124 digit->Calibrated(kFALSE);
125 }
7deb8eb0 126
127 // refit the tracks from digits
128 AliMUONVTrackStore* newTrackStore = refitter.ReconstructFromDigits();
129
130 //----------------------------------------------//
131 // ------ fill new ESD and print results ------ //
132 //----------------------------------------------//
133 // loop over the list of ESD tracks
134 TClonesArray *esdTracks = (TClonesArray*) esd->FindListObject("MuonTracks");
135 for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
136
137 // get the ESD track
138 AliESDMuonTrack* esdTrack = (AliESDMuonTrack*) esdTracks->UncheckedAt(iTrack);
b1fea02e 139
140 // skip ghost tracks (leave them unchanged in the new ESD file)
141 if (!esdTrack->ContainTrackerData()) continue;
142
7deb8eb0 143 // get the corresponding MUON track
630711ed 144 AliMUONTrack* track = esdInterface.FindTrack(esdTrack->GetUniqueID());
b1fea02e 145
7deb8eb0 146 // Find the corresponding re-fitted MUON track
147 AliMUONTrack* newTrack = (AliMUONTrack*) newTrackStore->FindObject(esdTrack->GetUniqueID());
148
149 // replace the content of the current ESD track or remove it if the refitting has failed
150 if (newTrack) {
151 Double_t vertex[3] = {esdTrack->GetNonBendingCoor(), esdTrack->GetBendingCoor(), esdTrack->GetZ()};
152 AliMUONESDInterface::MUONToESD(*newTrack, *esdTrack, vertex, esdInterface.GetDigits());
153 } else {
154 esdTracks->Remove(esdTrack);
155 }
156
157 // print initial and re-fitted track parameters at first cluster if any
158 if (printLevel>0) {
159 cout<<" ----------------track #"<<iTrack+1<<"----------------"<<endl;
160 cout<<"before refit:"<<endl;
161 AliMUONTrackParam *param = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->First();
162 param->Print("FULL");
163 if (printLevel>1) param->GetCovariances().Print();
164 if (!newTrack) continue;
165 cout<<"after refit:"<<endl;
166 param = (AliMUONTrackParam*) newTrack->GetTrackParamAtCluster()->First();
167 param->Print("FULL");
168 if (printLevel>1) param->GetCovariances().Print();
169 cout<<" ----------------------------------------"<<endl;
170 }
171
172 }
173
174 // free memory
175 delete newTrackStore;
176
177 // fill new ESD tree with new tracks
178 esdTracks->Compress();
179 newESDTree->Fill();
180
181 if (printLevel>0) cout<<" ****************************************"<<endl;
182 }
183
184 // ...timer stop
185 timer.Stop();
186
187 // write output ESD tree
ad3c6eda 188 newESDFile->cd();
7deb8eb0 189 newESDTree->Write();
ad3c6eda 190 delete newESDTree;
7deb8eb0 191 newESDFile->Close();
192
193 // free memory
194 esdFile->Close();
7deb8eb0 195 delete esd;
196
197 cout<<endl<<"time to refit: R:"<<timer.RealTime()<<" C:"<<timer.CpuTime()<<endl<<endl;
198}
199
200//-----------------------------------------------------------------------
201void Prepare()
202{
203 /// Set the geometry, the magnetic field, the mapping and the reconstruction parameters
204
205 // Import TGeo geometry (needed by AliMUONTrackExtrap::ExtrapToVertex)
206 if (!gGeoManager) {
207 AliGeomManager::LoadGeometry("geometry.root");
208 if (!gGeoManager) {
209 Error("MUONRefit", "getting geometry from file %s failed", "generated/galice.root");
210 return;
211 }
212 }
213
214 // set mag field
215 printf("Loading field map...\n");
216 AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 1, 1., 10., AliMagFMaps::k5kG);
217 AliTracker::SetFieldMap(field, kFALSE);
218 AliMUONTrackExtrap::SetField(AliTracker::GetFieldMap());
219
220 // Load mapping
221 AliCDBManager* man = AliCDBManager::Instance();
222 man->SetDefaultStorage("local://$ALICE_ROOT");
223 man->SetRun(0);
224 if ( ! AliMpCDB::LoadDDLStore() ) {
225 Error("MUONRefit","Could not access mapping from OCDB !");
226 exit(-1);
227 }
228
7deb8eb0 229}
230
231//-----------------------------------------------------------------------
232TTree* GetESDTree(TFile *esdFile)
233{
234 /// Check that the file is properly open
235 /// Return pointer to the ESD Tree
236
237 if (!esdFile || !esdFile->IsOpen()) {
238 Error("GetESDTree", "opening ESD file failed");
239 exit(-1);
240 }
241
242 TTree* tree = (TTree*) esdFile->Get("esdTree");
243 if (!tree) {
244 Error("GetESDTree", "no ESD tree found");
245 exit(-1);
246 }
247
248 return tree;
249
250}
251