a1b85eabe5fadad37b40380c0a8d7c56ed6c4f0f
[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 "AliMagF.h"
36 #include "AliTracker.h"
37 #include "AliESDEvent.h"
38 #include "AliESDMuonTrack.h"
39 #include "AliRecoParam.h"
40 #include "AliCDBManager.h"
41 #include "AliCDBEntry.h"
42 #include "AliCDBPath.h"
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"
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   
74   // reconstruction parameters for the refitting
75   AliMUONRecoParam* recoParam = AliMUONRecoParam::GetLowFluxParam();
76   Info("MUONRefit", "\n Reconstruction parameters for refitting:");
77   recoParam->Print("FULL");
78   
79   AliMUONESDInterface esdInterface;
80   AliMUONRefitter refitter(recoParam);
81   refitter.Connect(&esdInterface);
82   
83   // open the ESD file and tree
84   TFile* esdFile = TFile::Open(esdFileNameIn);
85   TTree* esdTree = GetESDTree(esdFile);
86   
87   // create the ESD output file and tree
88   TFile* newESDFile = TFile::Open(esdFileNameOut, "RECREATE");
89   newESDFile->SetCompressionLevel(2);
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) {
111       Error("MUONRefit", "no ESD object found for event %d", iEvent);
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());
123     while ((digit = static_cast<AliMUONVDigit*>(next()))) {
124       digit->SetCharge(digit->ADC());
125       digit->Calibrated(kFALSE);
126     }
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);
140       
141       // skip ghost tracks (leave them unchanged in the new ESD file)
142       if (!esdTrack->ContainTrackerData()) continue;
143       
144       // get the corresponding MUON track
145       AliMUONTrack* track = esdInterface.FindTrack(esdTrack->GetUniqueID());
146       
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
189   newESDFile->cd();
190   newESDTree->Write();
191   delete newESDTree;
192   newESDFile->Close();
193   
194   // free memory
195   esdFile->Close();
196   delete esd;
197   
198   cout<<endl<<"time to refit: R:"<<timer.RealTime()<<" C:"<<timer.CpuTime()<<endl<<endl;
199 }
200
201 //-----------------------------------------------------------------------
202 void 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   
215   // set  mag field 
216   // waiting for mag field in CDB 
217   if (!TGeoGlobalMagField::Instance()->GetField()) {
218     printf("Loading field map...\n");
219     AliMagF* field = new AliMagF("Maps","Maps",2,1.,1., 10.,AliMagF::k5kG);
220     TGeoGlobalMagField::Instance()->SetField(field);
221   }
222   // set the magnetic field for track extrapolations
223   AliMUONTrackExtrap::SetField();
224   
225   // Load mapping
226   AliCDBManager* man = AliCDBManager::Instance();
227   man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
228   man->SetRun(0);
229   if ( ! AliMpCDB::LoadDDLStore() ) {
230     Error("MUONRefit","Could not access mapping from OCDB !");
231     exit(-1);
232   }
233   
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   
251 }
252
253 //-----------------------------------------------------------------------
254 TTree* 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