changes in the MagF constructor
[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
f7a1cc68 35#include "AliMagF.h"
7deb8eb0 36#include "AliTracker.h"
37#include "AliESDEvent.h"
38#include "AliESDMuonTrack.h"
39#include "AliRecoParam.h"
40#include "AliCDBManager.h"
17e52d9f 41#include "AliCDBEntry.h"
42#include "AliCDBPath.h"
7deb8eb0 43#include "AliGeomManager.h"
44
45// MUON includes
46#include "AliMpCDB.h"
47#include "AliMUONRecoParam.h"
48#include "AliMUONESDInterface.h"
49#include "AliMUONRefitter.h"
50#include "AliMUONVDigit.h"
7deb8eb0 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
17e52d9f 74 // reconstruction parameters for the refitting
75 AliMUONRecoParam* recoParam = AliMUONRecoParam::GetLowFluxParam();
76 Info("MUONRefit", "\n Reconstruction parameters for refitting:");
77 recoParam->Print("FULL");
ad3c6eda 78
7deb8eb0 79 AliMUONESDInterface esdInterface;
17e52d9f 80 AliMUONRefitter refitter(recoParam);
7deb8eb0 81 refitter.Connect(&esdInterface);
82
83 // open the ESD file and tree
84 TFile* esdFile = TFile::Open(esdFileNameIn);
85 TTree* esdTree = GetESDTree(esdFile);
86
ad3c6eda 87 // create the ESD output file and tree
88 TFile* newESDFile = TFile::Open(esdFileNameOut, "RECREATE");
89 newESDFile->SetCompressionLevel(2);
7deb8eb0 90 TTree* newESDTree = esdTree->CloneTree(0);
91
92 // connect ESD event to the ESD tree
93 AliESDEvent* esd = new AliESDEvent();
94 esd->ReadFromTree(esdTree);
95
96 // timer start...
97 TStopwatch timer;
98
99 // Loop over ESD events
100 if (nevents > 0) nevents = TMath::Min(nevents,(Int_t)esdTree->GetEntries());
101 else nevents = (Int_t)esdTree->GetEntries();
102 for (Int_t iEvent = 0; iEvent < nevents; iEvent++) {
103 if (printLevel>0) cout<<endl<<" ****************event #"<<iEvent+1<<"****************"<<endl;
104
105 //----------------------------------------------//
106 // -------------- process event --------------- //
107 //----------------------------------------------//
108 // get the ESD of current event
109 esdTree->GetEvent(iEvent);
110 if (!esd) {
17e52d9f 111 Error("MUONRefit", "no ESD object found for event %d", iEvent);
7deb8eb0 112 return;
113 }
114 Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks();
115 if (nTracks < 1) continue;
116
117 // load the current event
118 esdInterface.LoadEvent(*esd);
119
120 // loop over digit to modify their charge
121 AliMUONVDigit *digit;
122 TIter next(esdInterface.CreateDigitIterator());
b1fea02e 123 while ((digit = static_cast<AliMUONVDigit*>(next()))) {
124 digit->SetCharge(digit->ADC());
125 digit->Calibrated(kFALSE);
126 }
7deb8eb0 127
128 // refit the tracks from digits
129 AliMUONVTrackStore* newTrackStore = refitter.ReconstructFromDigits();
130
131 //----------------------------------------------//
132 // ------ fill new ESD and print results ------ //
133 //----------------------------------------------//
134 // loop over the list of ESD tracks
135 TClonesArray *esdTracks = (TClonesArray*) esd->FindListObject("MuonTracks");
136 for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
137
138 // get the ESD track
139 AliESDMuonTrack* esdTrack = (AliESDMuonTrack*) esdTracks->UncheckedAt(iTrack);
b1fea02e 140
141 // skip ghost tracks (leave them unchanged in the new ESD file)
142 if (!esdTrack->ContainTrackerData()) continue;
143
7deb8eb0 144 // get the corresponding MUON track
630711ed 145 AliMUONTrack* track = esdInterface.FindTrack(esdTrack->GetUniqueID());
b1fea02e 146
7deb8eb0 147 // Find the corresponding re-fitted MUON track
148 AliMUONTrack* newTrack = (AliMUONTrack*) newTrackStore->FindObject(esdTrack->GetUniqueID());
149
150 // replace the content of the current ESD track or remove it if the refitting has failed
151 if (newTrack) {
152 Double_t vertex[3] = {esdTrack->GetNonBendingCoor(), esdTrack->GetBendingCoor(), esdTrack->GetZ()};
153 AliMUONESDInterface::MUONToESD(*newTrack, *esdTrack, vertex, esdInterface.GetDigits());
154 } else {
155 esdTracks->Remove(esdTrack);
156 }
157
158 // print initial and re-fitted track parameters at first cluster if any
159 if (printLevel>0) {
160 cout<<" ----------------track #"<<iTrack+1<<"----------------"<<endl;
161 cout<<"before refit:"<<endl;
162 AliMUONTrackParam *param = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->First();
163 param->Print("FULL");
164 if (printLevel>1) param->GetCovariances().Print();
165 if (!newTrack) continue;
166 cout<<"after refit:"<<endl;
167 param = (AliMUONTrackParam*) newTrack->GetTrackParamAtCluster()->First();
168 param->Print("FULL");
169 if (printLevel>1) param->GetCovariances().Print();
170 cout<<" ----------------------------------------"<<endl;
171 }
172
173 }
174
175 // free memory
176 delete newTrackStore;
177
178 // fill new ESD tree with new tracks
179 esdTracks->Compress();
180 newESDTree->Fill();
181
182 if (printLevel>0) cout<<" ****************************************"<<endl;
183 }
184
185 // ...timer stop
186 timer.Stop();
187
188 // write output ESD tree
ad3c6eda 189 newESDFile->cd();
7deb8eb0 190 newESDTree->Write();
ad3c6eda 191 delete newESDTree;
7deb8eb0 192 newESDFile->Close();
193
194 // free memory
195 esdFile->Close();
7deb8eb0 196 delete esd;
197
198 cout<<endl<<"time to refit: R:"<<timer.RealTime()<<" C:"<<timer.CpuTime()<<endl<<endl;
199}
200
201//-----------------------------------------------------------------------
202void Prepare()
203{
204 /// Set the geometry, the magnetic field, the mapping and the reconstruction parameters
205
206 // Import TGeo geometry (needed by AliMUONTrackExtrap::ExtrapToVertex)
207 if (!gGeoManager) {
208 AliGeomManager::LoadGeometry("geometry.root");
209 if (!gGeoManager) {
210 Error("MUONRefit", "getting geometry from file %s failed", "generated/galice.root");
211 return;
212 }
213 }
214
f7a1cc68 215 // set mag field
216 // waiting for mag field in CDB
217 if (!TGeoGlobalMagField::Instance()->GetField()) {
218 printf("Loading field map...\n");
4642ac4b 219 AliMagF* field = new AliMagF("Maps","Maps",1.,1.,AliMagF::k5kG);
f7a1cc68 220 TGeoGlobalMagField::Instance()->SetField(field);
221 }
222 // set the magnetic field for track extrapolations
223 AliMUONTrackExtrap::SetField();
7deb8eb0 224
225 // Load mapping
226 AliCDBManager* man = AliCDBManager::Instance();
162637e4 227 man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
7deb8eb0 228 man->SetRun(0);
229 if ( ! AliMpCDB::LoadDDLStore() ) {
230 Error("MUONRefit","Could not access mapping from OCDB !");
231 exit(-1);
232 }
233
17e52d9f 234 // Load initial reconstruction parameters from OCDB
235 AliMUONRecoParam* recoParam = 0x0;
236 AliCDBPath path("MUON","Calib","RecoParam");
237 AliCDBEntry *entry=man->Get(path.GetPath());
238 if(entry) {
239 recoParam = dynamic_cast<AliMUONRecoParam*>(entry->GetObject());
240 entry->SetOwner(0);
241 AliCDBManager::Instance()->UnloadFromCache(path.GetPath());
242 }
243 if (!recoParam) {
244 printf("Couldn't find RecoParam object in OCDB: create default one");
245 recoParam = AliMUONRecoParam::GetLowFluxParam();
246 }
247 Info("MUONRefit", "\n initial recontruction parameters:");
248 recoParam->Print("FULL");
249 AliMUONESDInterface::ResetTracker(recoParam);
250
7deb8eb0 251}
252
253//-----------------------------------------------------------------------
254TTree* GetESDTree(TFile *esdFile)
255{
256 /// Check that the file is properly open
257 /// Return pointer to the ESD Tree
258
259 if (!esdFile || !esdFile->IsOpen()) {
260 Error("GetESDTree", "opening ESD file failed");
261 exit(-1);
262 }
263
264 TTree* tree = (TTree*) esdFile->Get("esdTree");
265 if (!tree) {
266 Error("GetESDTree", "no ESD tree found");
267 exit(-1);
268 }
269
270 return tree;
271
272}
273