1) New class "AliMUONRefitter" to:
[u/mrichter/AliRoot.git] / MUON / MUONRefit.C
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
57 const Int_t printLevel = 1;
58
59 void Prepare();
60 TTree* GetESDTree(TFile *esdFile);
61
62 //-----------------------------------------------------------------------
63 void 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();
73   AliMUONESDInterface esdInterface;
74   AliMUONRefitter refitter;
75   refitter.Connect(&esdInterface);
76   
77   // open the ESD file and tree
78   TFile* esdFile = TFile::Open(esdFileNameIn);
79   TTree* esdTree = GetESDTree(esdFile);
80   
81   // create the ESD output tree
82   gROOT->cd();
83   TTree* newESDTree = esdTree->CloneTree(0);
84   
85   // connect ESD event to the ESD tree
86   AliESDEvent* esd = new AliESDEvent();
87   esd->ReadFromTree(esdTree);
88
89   // timer start...
90   TStopwatch timer;
91   
92   // Loop over ESD events
93   if (nevents > 0) nevents = TMath::Min(nevents,(Int_t)esdTree->GetEntries());
94   else nevents = (Int_t)esdTree->GetEntries();
95   for (Int_t iEvent = 0; iEvent < nevents; iEvent++) {
96     if (printLevel>0) cout<<endl<<"            ****************event #"<<iEvent+1<<"****************"<<endl;
97     
98     //----------------------------------------------//
99     // -------------- process event --------------- //
100     //----------------------------------------------//
101     // get the ESD of current event
102     esdTree->GetEvent(iEvent);
103     if (!esd) {
104       Error("CheckESD", "no ESD object found for event %d", iEvent);
105       return;
106     }
107     Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks();
108     if (nTracks < 1) continue;
109     
110     // load the current event
111     esdInterface.LoadEvent(*esd);
112     
113     // loop over digit to modify their charge
114     AliMUONVDigit *digit;
115     TIter next(esdInterface.CreateDigitIterator());
116     while ((digit = static_cast<AliMUONVDigit*>(next()))) digit->SetCharge(digit->ADC());
117     
118     // refit the tracks from digits
119     AliMUONVTrackStore* newTrackStore = refitter.ReconstructFromDigits();
120     
121     //----------------------------------------------//
122     // ------ fill new ESD and print results ------ //
123     //----------------------------------------------//
124     // loop over the list of ESD tracks
125     TClonesArray *esdTracks = (TClonesArray*) esd->FindListObject("MuonTracks");
126     for (Int_t iTrack = 0; iTrack <  nTracks; iTrack++) {
127       
128       // get the ESD track
129       AliESDMuonTrack* esdTrack = (AliESDMuonTrack*) esdTracks->UncheckedAt(iTrack);
130       // get the corresponding MUON track
131       AliMUONTrack* track = esdInterface.GetTrackFast(iTrack);
132       // Find the corresponding re-fitted MUON track
133       AliMUONTrack* newTrack = (AliMUONTrack*) newTrackStore->FindObject(esdTrack->GetUniqueID());
134       
135       // replace the content of the current ESD track or remove it if the refitting has failed
136       if (newTrack) {
137         Double_t vertex[3] = {esdTrack->GetNonBendingCoor(), esdTrack->GetBendingCoor(), esdTrack->GetZ()};
138         AliMUONESDInterface::MUONToESD(*newTrack, *esdTrack, vertex, esdInterface.GetDigits());
139       } else {
140         esdTracks->Remove(esdTrack);
141       }
142       
143       // print initial and re-fitted track parameters at first cluster if any
144       if (printLevel>0) {
145         cout<<"            ----------------track #"<<iTrack+1<<"----------------"<<endl;
146         cout<<"before refit:"<<endl;
147         AliMUONTrackParam *param = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->First();
148         param->Print("FULL");
149         if (printLevel>1) param->GetCovariances().Print();
150         if (!newTrack) continue;
151         cout<<"after refit:"<<endl;
152         param = (AliMUONTrackParam*) newTrack->GetTrackParamAtCluster()->First();
153         param->Print("FULL");
154         if (printLevel>1) param->GetCovariances().Print();
155         cout<<"            ----------------------------------------"<<endl;
156       }
157       
158     }
159     
160     // free memory
161     delete newTrackStore;
162     
163     // fill new ESD tree with new tracks
164     esdTracks->Compress();
165     newESDTree->Fill();
166     
167     if (printLevel>0) cout<<"            ****************************************"<<endl;
168   }
169   
170   // ...timer stop
171   timer.Stop();
172   
173   // write output ESD tree
174   TFile* newESDFile = TFile::Open(esdFileNameOut, "RECREATE");
175   newESDFile->SetCompressionLevel(2);
176   newESDTree->Write();
177   newESDFile->Close();
178   
179   // free memory
180   esdFile->Close();
181   delete newESDTree;
182   delete esd;
183   
184   cout<<endl<<"time to refit: R:"<<timer.RealTime()<<" C:"<<timer.CpuTime()<<endl<<endl;
185 }
186
187 //-----------------------------------------------------------------------
188 void Prepare()
189 {
190   /// Set the geometry, the magnetic field, the mapping and the reconstruction parameters
191   
192   // Import TGeo geometry (needed by AliMUONTrackExtrap::ExtrapToVertex)
193   if (!gGeoManager) {
194     AliGeomManager::LoadGeometry("geometry.root");
195     if (!gGeoManager) {
196       Error("MUONRefit", "getting geometry from file %s failed", "generated/galice.root");
197       return;
198     }
199   }
200   
201   // set mag field
202   printf("Loading field map...\n");
203   AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 1, 1., 10., AliMagFMaps::k5kG);
204   AliTracker::SetFieldMap(field, kFALSE);
205   AliMUONTrackExtrap::SetField(AliTracker::GetFieldMap());
206   
207   // Load mapping
208   AliCDBManager* man = AliCDBManager::Instance();
209   man->SetDefaultStorage("local://$ALICE_ROOT");
210   man->SetRun(0);
211   if ( ! AliMpCDB::LoadDDLStore() ) {
212     Error("MUONRefit","Could not access mapping from OCDB !");
213     exit(-1);
214   }
215   
216   // set reconstruction parameters
217   AliMUONRecoParam *muonRecoParam = AliMUONRecoParam::GetLowFluxParam();
218   muonRecoParam->ImproveTracks(kFALSE);
219   muonRecoParam->Print("FULL");
220   AliRecoParam::Instance()->RegisterRecoParam(muonRecoParam);
221   
222 }
223
224 //-----------------------------------------------------------------------
225 TTree* GetESDTree(TFile *esdFile)
226 {
227   /// Check that the file is properly open
228   /// Return pointer to the ESD Tree
229   
230   if (!esdFile || !esdFile->IsOpen()) {
231     Error("GetESDTree", "opening ESD file failed");
232     exit(-1);
233   }
234   
235   TTree* tree = (TTree*) esdFile->Get("esdTree");
236   if (!tree) {
237     Error("GetESDTree", "no ESD tree found");
238     exit(-1);
239   }
240   
241   return tree;
242   
243 }
244