renaming function to avoid library conflict (discovered in test/generators/TUHKMgen)
[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
7deb8eb0 35#include "AliESDEvent.h"
36#include "AliESDMuonTrack.h"
7deb8eb0 37#include "AliCDBManager.h"
38#include "AliGeomManager.h"
39
40// MUON includes
a99c3449 41#include "AliMUONCDB.h"
7deb8eb0 42#include "AliMUONRecoParam.h"
43#include "AliMUONESDInterface.h"
44#include "AliMUONRefitter.h"
45#include "AliMUONVDigit.h"
7deb8eb0 46#include "AliMUONTrack.h"
47#include "AliMUONVTrackStore.h"
48#include "AliMUONTrackParam.h"
7deb8eb0 49#endif
50
51const Int_t printLevel = 1;
69877654 52const Bool_t reconstructFromDigits = kTRUE; // kFALSE = reconstruct from clusters
7deb8eb0 53
7deb8eb0 54TTree* GetESDTree(TFile *esdFile);
55
56//-----------------------------------------------------------------------
a99c3449 57void MUONRefit(Int_t nevents = -1, const char* esdFileNameIn = "AliESDs.root", const char* esdFileNameOut = "AliESDs_New.root",
58 const char* geoFilename = "geometry.root", const char* ocdbPath = "local://$ALICE_ROOT/OCDB")
7deb8eb0 59{
a99c3449 60 /// Example of muon refitting:
7deb8eb0 61 /// refit ESD tracks from ESD pads (i.e. re-clusterized the attached ESD clusters);
62 /// reset the charge of the digit using their raw charge before refitting;
63 /// compare results with original ESD tracks;
64 /// write results in a new ESD file
65
7deb8eb0 66 // open the ESD file and tree
67 TFile* esdFile = TFile::Open(esdFileNameIn);
68 TTree* esdTree = GetESDTree(esdFile);
69
fe0324de 70 // connect ESD event to the ESD tree
71 AliESDEvent* esd = new AliESDEvent();
72 esd->ReadFromTree(esdTree);
73
ad3c6eda 74 // create the ESD output file and tree
75 TFile* newESDFile = TFile::Open(esdFileNameOut, "RECREATE");
76 newESDFile->SetCompressionLevel(2);
7deb8eb0 77
fe0324de 78 // connect ESD event to the ESD tree (recreate track branch for backward compatibility)
79 esdTree->SetBranchStatus("*MuonTracks*",0);
80 TTree* newESDTree = esdTree->CloneTree(0);
81 esdTree->SetBranchStatus("*MuonTracks*",1);
82 esd->WriteToTree(newESDTree);
a99c3449 83
84 // get run number
85 if (esdTree->GetEvent(0) <= 0) {
86 Error("MUONRefit", "no ESD object found for event 0");
87 return;
88 }
89 Int_t runNumber = esd->GetRunNumber();
90
91 // Import TGeo geometry
92 if (!gGeoManager) {
93 AliGeomManager::LoadGeometry(geoFilename);
94 if (!gGeoManager) {
95 Error("MUONRefit", "getting geometry from file %s failed", "generated/galice.root");
96 exit(-1);
97 }
98 }
99
100 // load necessary data from OCDB
101 AliCDBManager::Instance()->SetDefaultStorage(ocdbPath);
102 AliCDBManager::Instance()->SetRun(runNumber);
103 if (!AliMUONCDB::LoadField()) return;
104 if (!AliMUONCDB::LoadMapping(kTRUE)) return;
105
106 // reconstruction parameters for the refitting
107 AliMUONRecoParam* recoParam = AliMUONRecoParam::GetLowFluxParam();
108 Info("MUONRefit", "\n Reconstruction parameters for refitting:");
109 recoParam->Print("FULL");
110
111 AliMUONESDInterface esdInterface;
112 AliMUONRefitter refitter(recoParam);
113 refitter.Connect(&esdInterface);
114 gRandom->SetSeed(1);
115
7deb8eb0 116 // timer start...
117 TStopwatch timer;
118
119 // Loop over ESD events
120 if (nevents > 0) nevents = TMath::Min(nevents,(Int_t)esdTree->GetEntries());
121 else nevents = (Int_t)esdTree->GetEntries();
122 for (Int_t iEvent = 0; iEvent < nevents; iEvent++) {
123 if (printLevel>0) cout<<endl<<" ****************event #"<<iEvent+1<<"****************"<<endl;
124
125 //----------------------------------------------//
126 // -------------- process event --------------- //
127 //----------------------------------------------//
128 // get the ESD of current event
129 esdTree->GetEvent(iEvent);
130 if (!esd) {
17e52d9f 131 Error("MUONRefit", "no ESD object found for event %d", iEvent);
7deb8eb0 132 return;
133 }
134 Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks();
135 if (nTracks < 1) continue;
136
137 // load the current event
a99c3449 138 esdInterface.LoadEvent(*esd, kFALSE);
7deb8eb0 139
fe0324de 140 // remove digits and clusters from ESD
141 esd->FindListObject("MuonClusters")->Clear("C");
142 esd->FindListObject("MuonPads")->Clear("C");
143
69877654 144 AliMUONVTrackStore* newTrackStore = 0x0;
145 if (reconstructFromDigits) {
146
147 // loop over digits to modify their charge
148 AliMUONVDigit *digit;
149 TIter next(esdInterface.CreateDigitIterator());
150 while ((digit = static_cast<AliMUONVDigit*>(next()))) {
151 digit->SetCharge(digit->ADC());
152 digit->Calibrated(kFALSE);
153 }
154
155 // refit the tracks from digits
156 refitter.SetFirstClusterIndex(0);
157 newTrackStore = refitter.ReconstructFromDigits();
158
159 } else {
160
161 // refit the tracks from clusters
162 newTrackStore = refitter.ReconstructFromClusters();
163
b1fea02e 164 }
7deb8eb0 165
7deb8eb0 166 //----------------------------------------------//
167 // ------ fill new ESD and print results ------ //
168 //----------------------------------------------//
169 // loop over the list of ESD tracks
170 TClonesArray *esdTracks = (TClonesArray*) esd->FindListObject("MuonTracks");
171 for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
172
173 // get the ESD track
174 AliESDMuonTrack* esdTrack = (AliESDMuonTrack*) esdTracks->UncheckedAt(iTrack);
b1fea02e 175
176 // skip ghost tracks (leave them unchanged in the new ESD file)
177 if (!esdTrack->ContainTrackerData()) continue;
178
7deb8eb0 179 // get the corresponding MUON track
630711ed 180 AliMUONTrack* track = esdInterface.FindTrack(esdTrack->GetUniqueID());
b1fea02e 181
7deb8eb0 182 // Find the corresponding re-fitted MUON track
183 AliMUONTrack* newTrack = (AliMUONTrack*) newTrackStore->FindObject(esdTrack->GetUniqueID());
184
185 // replace the content of the current ESD track or remove it if the refitting has failed
fe0324de 186 if (newTrack && (!recoParam->ImproveTracks() || newTrack->IsImproved())) {
187
188 // fill the track info
7deb8eb0 189 Double_t vertex[3] = {esdTrack->GetNonBendingCoor(), esdTrack->GetBendingCoor(), esdTrack->GetZ()};
fe0324de 190 AliMUONESDInterface::MUONToESD(*newTrack, *esdTrack, vertex);
191
192 // add the clusters (and the digits if previously there)
193 for (Int_t i = 0; i < newTrack->GetNClusters(); i++) {
194 AliMUONTrackParam *param = static_cast<AliMUONTrackParam*>(newTrack->GetTrackParamAtCluster()->UncheckedAt(i));
195 AliMUONESDInterface::MUONToESD(*(param->GetClusterPtr()), *esd, esdInterface.GetDigits());
196 }
197
198 } else esdTracks->Remove(esdTrack);
7deb8eb0 199
200 // print initial and re-fitted track parameters at first cluster if any
201 if (printLevel>0) {
202 cout<<" ----------------track #"<<iTrack+1<<"----------------"<<endl;
203 cout<<"before refit:"<<endl;
204 AliMUONTrackParam *param = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->First();
205 param->Print("FULL");
206 if (printLevel>1) param->GetCovariances().Print();
207 if (!newTrack) continue;
208 cout<<"after refit:"<<endl;
209 param = (AliMUONTrackParam*) newTrack->GetTrackParamAtCluster()->First();
210 param->Print("FULL");
211 if (printLevel>1) param->GetCovariances().Print();
212 cout<<" ----------------------------------------"<<endl;
213 }
214
215 }
216
217 // free memory
218 delete newTrackStore;
219
220 // fill new ESD tree with new tracks
221 esdTracks->Compress();
222 newESDTree->Fill();
223
224 if (printLevel>0) cout<<" ****************************************"<<endl;
225 }
226
227 // ...timer stop
228 timer.Stop();
229
230 // write output ESD tree
ad3c6eda 231 newESDFile->cd();
7deb8eb0 232 newESDTree->Write();
ad3c6eda 233 delete newESDTree;
7deb8eb0 234 newESDFile->Close();
235
236 // free memory
237 esdFile->Close();
7deb8eb0 238 delete esd;
a99c3449 239 delete recoParam;
7deb8eb0 240
241 cout<<endl<<"time to refit: R:"<<timer.RealTime()<<" C:"<<timer.CpuTime()<<endl<<endl;
242}
243
244//-----------------------------------------------------------------------
7deb8eb0 245TTree* GetESDTree(TFile *esdFile)
246{
247 /// Check that the file is properly open
248 /// Return pointer to the ESD Tree
249
250 if (!esdFile || !esdFile->IsOpen()) {
251 Error("GetESDTree", "opening ESD file failed");
252 exit(-1);
253 }
254
255 TTree* tree = (TTree*) esdFile->Get("esdTree");
256 if (!tree) {
257 Error("GetESDTree", "no ESD tree found");
258 exit(-1);
259 }
260
261 return tree;
262
263}
264