1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
20 /// \brief Macro for refitting ESD tracks from ESD pads
22 /// \author Philippe Pillot, SUBATECH
24 #if !defined(__CINT__) || defined(__MAKECINT__)
25 #include <TStopwatch.h>
29 #include <Riostream.h>
30 #include <TGeoManager.h>
35 #include "AliESDEvent.h"
36 #include "AliESDMuonTrack.h"
37 #include "AliCDBManager.h"
38 #include "AliGeomManager.h"
41 #include "AliMUONCDB.h"
42 #include "AliMUONRecoParam.h"
43 #include "AliMUONESDInterface.h"
44 #include "AliMUONRefitter.h"
45 #include "AliMUONVDigit.h"
46 #include "AliMUONTrack.h"
47 #include "AliMUONVTrackStore.h"
48 #include "AliMUONTrackParam.h"
51 const Int_t printLevel = 1;
53 TTree* GetESDTree(TFile *esdFile);
55 //-----------------------------------------------------------------------
56 void MUONRefit(Int_t nevents = -1, const char* esdFileNameIn = "AliESDs.root", const char* esdFileNameOut = "AliESDs_New.root",
57 const char* geoFilename = "geometry.root", const char* ocdbPath = "local://$ALICE_ROOT/OCDB")
59 /// Example of muon refitting:
60 /// refit ESD tracks from ESD pads (i.e. re-clusterized the attached ESD clusters);
61 /// reset the charge of the digit using their raw charge before refitting;
62 /// compare results with original ESD tracks;
63 /// write results in a new ESD file
65 // open the ESD file and tree
66 TFile* esdFile = TFile::Open(esdFileNameIn);
67 TTree* esdTree = GetESDTree(esdFile);
69 // connect ESD event to the ESD tree
70 AliESDEvent* esd = new AliESDEvent();
71 esd->ReadFromTree(esdTree);
73 // create the ESD output file and tree
74 TFile* newESDFile = TFile::Open(esdFileNameOut, "RECREATE");
75 newESDFile->SetCompressionLevel(2);
77 // connect ESD event to the ESD tree (recreate track branch for backward compatibility)
78 esdTree->SetBranchStatus("*MuonTracks*",0);
79 TTree* newESDTree = esdTree->CloneTree(0);
80 esdTree->SetBranchStatus("*MuonTracks*",1);
81 esd->WriteToTree(newESDTree);
84 if (esdTree->GetEvent(0) <= 0) {
85 Error("MUONRefit", "no ESD object found for event 0");
88 Int_t runNumber = esd->GetRunNumber();
90 // Import TGeo geometry
92 AliGeomManager::LoadGeometry(geoFilename);
94 Error("MUONRefit", "getting geometry from file %s failed", "generated/galice.root");
99 // load necessary data from OCDB
100 AliCDBManager::Instance()->SetDefaultStorage(ocdbPath);
101 AliCDBManager::Instance()->SetRun(runNumber);
102 if (!AliMUONCDB::LoadField()) return;
103 if (!AliMUONCDB::LoadMapping(kTRUE)) return;
105 // reconstruction parameters for the refitting
106 AliMUONRecoParam* recoParam = AliMUONRecoParam::GetLowFluxParam();
107 Info("MUONRefit", "\n Reconstruction parameters for refitting:");
108 recoParam->Print("FULL");
110 AliMUONESDInterface esdInterface;
111 AliMUONRefitter refitter(recoParam);
112 refitter.Connect(&esdInterface);
118 // Loop over ESD events
119 if (nevents > 0) nevents = TMath::Min(nevents,(Int_t)esdTree->GetEntries());
120 else nevents = (Int_t)esdTree->GetEntries();
121 for (Int_t iEvent = 0; iEvent < nevents; iEvent++) {
122 if (printLevel>0) cout<<endl<<" ****************event #"<<iEvent+1<<"****************"<<endl;
124 //----------------------------------------------//
125 // -------------- process event --------------- //
126 //----------------------------------------------//
127 // get the ESD of current event
128 esdTree->GetEvent(iEvent);
130 Error("MUONRefit", "no ESD object found for event %d", iEvent);
133 Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks();
134 if (nTracks < 1) continue;
136 // load the current event
137 esdInterface.LoadEvent(*esd, kFALSE);
139 // remove digits and clusters from ESD
140 esd->FindListObject("MuonClusters")->Clear("C");
141 esd->FindListObject("MuonPads")->Clear("C");
143 // loop over digits to modify their charge
144 AliMUONVDigit *digit;
145 TIter next(esdInterface.CreateDigitIterator());
146 while ((digit = static_cast<AliMUONVDigit*>(next()))) {
147 digit->SetCharge(digit->ADC());
148 digit->Calibrated(kFALSE);
151 // refit the tracks from digits
152 refitter.SetFirstClusterIndex(0);
153 AliMUONVTrackStore* newTrackStore = refitter.ReconstructFromDigits();
155 //----------------------------------------------//
156 // ------ fill new ESD and print results ------ //
157 //----------------------------------------------//
158 // loop over the list of ESD tracks
159 TClonesArray *esdTracks = (TClonesArray*) esd->FindListObject("MuonTracks");
160 for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
163 AliESDMuonTrack* esdTrack = (AliESDMuonTrack*) esdTracks->UncheckedAt(iTrack);
165 // skip ghost tracks (leave them unchanged in the new ESD file)
166 if (!esdTrack->ContainTrackerData()) continue;
168 // get the corresponding MUON track
169 AliMUONTrack* track = esdInterface.FindTrack(esdTrack->GetUniqueID());
171 // Find the corresponding re-fitted MUON track
172 AliMUONTrack* newTrack = (AliMUONTrack*) newTrackStore->FindObject(esdTrack->GetUniqueID());
174 // replace the content of the current ESD track or remove it if the refitting has failed
175 if (newTrack && (!recoParam->ImproveTracks() || newTrack->IsImproved())) {
177 // fill the track info
178 Double_t vertex[3] = {esdTrack->GetNonBendingCoor(), esdTrack->GetBendingCoor(), esdTrack->GetZ()};
179 AliMUONESDInterface::MUONToESD(*newTrack, *esdTrack, vertex);
181 // add the clusters (and the digits if previously there)
182 for (Int_t i = 0; i < newTrack->GetNClusters(); i++) {
183 AliMUONTrackParam *param = static_cast<AliMUONTrackParam*>(newTrack->GetTrackParamAtCluster()->UncheckedAt(i));
184 AliMUONESDInterface::MUONToESD(*(param->GetClusterPtr()), *esd, esdInterface.GetDigits());
187 } else esdTracks->Remove(esdTrack);
189 // print initial and re-fitted track parameters at first cluster if any
191 cout<<" ----------------track #"<<iTrack+1<<"----------------"<<endl;
192 cout<<"before refit:"<<endl;
193 AliMUONTrackParam *param = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->First();
194 param->Print("FULL");
195 if (printLevel>1) param->GetCovariances().Print();
196 if (!newTrack) continue;
197 cout<<"after refit:"<<endl;
198 param = (AliMUONTrackParam*) newTrack->GetTrackParamAtCluster()->First();
199 param->Print("FULL");
200 if (printLevel>1) param->GetCovariances().Print();
201 cout<<" ----------------------------------------"<<endl;
207 delete newTrackStore;
209 // fill new ESD tree with new tracks
210 esdTracks->Compress();
213 if (printLevel>0) cout<<" ****************************************"<<endl;
219 // write output ESD tree
230 cout<<endl<<"time to refit: R:"<<timer.RealTime()<<" C:"<<timer.CpuTime()<<endl<<endl;
233 //-----------------------------------------------------------------------
234 TTree* GetESDTree(TFile *esdFile)
236 /// Check that the file is properly open
237 /// Return pointer to the ESD Tree
239 if (!esdFile || !esdFile->IsOpen()) {
240 Error("GetESDTree", "opening ESD file failed");
244 TTree* tree = (TTree*) esdFile->Get("esdTree");
246 Error("GetESDTree", "no ESD tree found");