PMD info on da/amore updated
[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 "AliESDEvent.h"
36 #include "AliESDMuonTrack.h"
37 #include "AliCDBManager.h"
38 #include "AliGeomManager.h"
39
40 // MUON includes
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"
49 #endif
50
51 const Int_t printLevel = 1;
52 const Bool_t reconstructFromDigits = kTRUE; // kFALSE = reconstruct from clusters
53
54 TTree* GetESDTree(TFile *esdFile);
55
56 //-----------------------------------------------------------------------
57 void 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")
59 {
60   /// Example of muon refitting:
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   
66   // open the ESD file and tree
67   TFile* esdFile = TFile::Open(esdFileNameIn);
68   TTree* esdTree = GetESDTree(esdFile);
69   
70   // connect ESD event to the ESD tree
71   AliESDEvent* esd = new AliESDEvent();
72   esd->ReadFromTree(esdTree);
73   
74   // create the ESD output file and tree
75   TFile* newESDFile = TFile::Open(esdFileNameOut, "RECREATE");
76   newESDFile->SetCompressionLevel(2);
77   
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);
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   
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) {
131       Error("MUONRefit", "no ESD object found for event %d", iEvent);
132       return;
133     }
134     Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks();
135     if (nTracks < 1) continue;
136     
137     // load the current event
138     esdInterface.LoadEvent(*esd, kFALSE);
139     
140     // remove digits and clusters from ESD
141     esd->FindListObject("MuonClusters")->Clear("C");
142     esd->FindListObject("MuonPads")->Clear("C");
143     
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       
164     }
165     
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);
175       
176       // skip ghost tracks (leave them unchanged in the new ESD file)
177       if (!esdTrack->ContainTrackerData()) continue;
178       
179       // get the corresponding MUON track
180       AliMUONTrack* track = esdInterface.FindTrack(esdTrack->GetUniqueID());
181       
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
186       if (newTrack && (!recoParam->ImproveTracks() || newTrack->IsImproved())) {
187         
188         // fill the track info
189         Double_t vertex[3] = {esdTrack->GetNonBendingCoor(), esdTrack->GetBendingCoor(), esdTrack->GetZ()};
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);
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
231   newESDFile->cd();
232   newESDTree->Write();
233   delete newESDTree;
234   newESDFile->Close();
235   
236   // free memory
237   esdFile->Close();
238   delete esd;
239   delete recoParam;
240   
241   cout<<endl<<"time to refit: R:"<<timer.RealTime()<<" C:"<<timer.CpuTime()<<endl<<endl;
242 }
243
244 //-----------------------------------------------------------------------
245 TTree* 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