Macro to fill new AliMUONClusterInfo objects from ESDs (Philippe P.)
[u/mrichter/AliRoot.git] / MUON / MUONClusterInfo.C
CommitLineData
fcefbac4 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 MUONClusterInfo.C
20/// \brief Macro to fill AliMUONClusterInfo objects
21///
22/// \author Philippe Pillot, SUBATECH
23
24#if !defined(__CINT__) || defined(__MAKECINT__)
25#include <TStopwatch.h>
26#include <TFile.h>
27#include <TClonesArray.h>
28#include <TTree.h>
29#include <TString.h>
30#include <Riostream.h>
31#include <TGeoManager.h>
32#include <TRandom.h>
33#include <TROOT.h>
34
35// STEER includes
36#include "AliMagFMaps.h"
37#include "AliTracker.h"
38#include "AliESDEvent.h"
39#include "AliESDMuonTrack.h"
40#include "AliRecoParam.h"
41#include "AliCDBManager.h"
42#include "AliGeomManager.h"
43
44// MUON includes
45#include "AliMpCDB.h"
46#include "AliMpSegmentation.h"
47#include "AliMpVSegmentation.h"
48#include "AliMpPad.h"
49#include "AliMUONCalibrationData.h"
50#include "AliMUONVCalibParam.h"
51#include "AliMUONPadInfo.h"
52#include "AliMUONClusterInfo.h"
53#include "AliMUONRecoParam.h"
54#include "AliMUONESDInterface.h"
55#include "AliMUONRefitter.h"
56#include "AliMUONVDigit.h"
57#include "AliMUONVDigitStore.h"
58#include "AliMUONVCluster.h"
59#include "AliMUONVClusterStore.h"
60#include "AliMUONTrack.h"
61#include "AliMUONVTrackStore.h"
62#include "AliMUONTrackParam.h"
63#include "AliMUONTrackExtrap.h"
64#endif
65
66const Int_t printLevel = 1;
67
68void Prepare();
69TTree* GetESDTree(TFile *esdFile);
70
71//-----------------------------------------------------------------------
72void MUONClusterInfo(Int_t nevents = -1, const char* esdFileName = "AliESDs.root", const char* outFileName = "clusterInfo.root")
73{
74 /// refit ESD tracks from ESD clusters to get track parameters at each attached cluster;
75 /// fill AliMUONESDClusterInfo object with ESD cluster + corresponding track parameters;
76 /// write results in a new root file.
77
78 AliMUONClusterInfo* clusterInfo = new AliMUONClusterInfo();
79 AliMUONPadInfo padInfo;
80 AliMUONCalibrationData* calibData = 0x0;
81
82 // prepare the refitting
83 gRandom->SetSeed(0);
84 Prepare();
85 AliMUONESDInterface esdInterface;
86 AliMUONRefitter refitter;
87 refitter.Connect(&esdInterface);
88
89 // open the ESD file and tree and connect the ESD event
90 TFile* esdFile = TFile::Open(esdFileName);
91 TTree* esdTree = GetESDTree(esdFile);
92 AliESDEvent* esd = new AliESDEvent();
93 esd->ReadFromTree(esdTree);
94
95 // prepare the output tree
96 gROOT->cd();
97 TFile* clusterInfoFile = TFile::Open(outFileName, "RECREATE");
98 clusterInfoFile->SetCompressionLevel(1);
99
100 TTree* clusterInfoTree = new TTree("clusterInfoTree","clusterInfoTree");
101 clusterInfoTree->Branch("clusterInfo", &clusterInfo, 32000, 99);
102
103 // timer start...
104 TStopwatch timer;
105
106 // Loop over ESD events
107 if (nevents > 0) nevents = TMath::Min(nevents,(Int_t)esdTree->GetEntries());
108 else nevents = (Int_t)esdTree->GetEntries();
109 for (Int_t iEvent = 0; iEvent < nevents; iEvent++) {
110
111 //----------------------------------------------//
112 // -------------- process event --------------- //
113 //----------------------------------------------//
114 // get the ESD of current event
115 esdTree->GetEvent(iEvent);
116 if (!esd) {
117 Error("CheckESD", "no ESD object found for event %d", iEvent);
118 return;
119 }
120 Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks();
121 if (nTracks < 1) continue;
122
123 // prepare access to calibration data
124 if (!calibData) calibData = new AliMUONCalibrationData(esd->GetESDRun()->GetRunNumber());
125
126 // load the current event
127 esdInterface.LoadEvent(*esd);
128
129 // get digit store
130 AliMUONVDigitStore* digitStore = esdInterface.GetDigits();
131
132 // refit the tracks from clusters
133 AliMUONVTrackStore* newTrackStore = refitter.ReconstructFromClusters();
134
135 //----------------------------------------------//
136 // ------------- fill cluster info ------------ //
137 //----------------------------------------------//
138 // loop over the refitted tracks
139 TIter nextTrack(newTrackStore->CreateIterator());
140 AliMUONTrack* newTrack;
141 while ((newTrack = static_cast<AliMUONTrack*>(nextTrack()))) {
142
143 // loop over clusters
144 AliMUONTrackParam* trackParam = static_cast<AliMUONTrackParam*>(newTrack->GetTrackParamAtCluster()->First());
145 while (trackParam) {
146 clusterInfo->Clear("C");
147
148 // fill cluster info
149 AliMUONVCluster* cluster = trackParam->GetClusterPtr();
150 clusterInfo->SetEventId(iEvent);
151 clusterInfo->SetZ(cluster->GetZ());
152 clusterInfo->SetClusterId(cluster->GetUniqueID());
153 clusterInfo->SetClusterXY(cluster->GetX(), cluster->GetY());
154 clusterInfo->SetClusterXYErr(cluster->GetErrX(), cluster->GetErrY());
155 clusterInfo->SetClusterChi2(cluster->GetChi2());
156 clusterInfo->SetClusterCharge(cluster->GetCharge());
157
158 // fill track info
159 clusterInfo->SetTrackId(newTrack->GetUniqueID());
160 clusterInfo->SetTrackXY(trackParam->GetNonBendingCoor(), trackParam->GetBendingCoor());
161 clusterInfo->SetTrackThetaXY(TMath::ATan(trackParam->GetBendingSlope()), TMath::ATan(trackParam->GetNonBendingSlope()));
162 clusterInfo->SetTrackP(trackParam->P());
163 const TMatrixD paramCov = trackParam->GetCovariances();
164 clusterInfo->SetTrackXYErr(TMath::Sqrt(paramCov(0,0)), TMath::Sqrt(paramCov(2,2)));
165 clusterInfo->SetTrackChi2(newTrack->GetNormalizedChi2());
166 clusterInfo->SetTrackCharge((Short_t)trackParam->GetCharge());
167
168 // fill pad info if available
169 for (Int_t i=0; i<cluster->GetNDigits(); i++) {
170 AliMUONVDigit* digit = digitStore->FindObject(cluster->GetDigitId(i));
171 if (!digit) continue;
172
173 // pad location
174 const AliMpVSegmentation* seg = AliMpSegmentation::Instance()->
175 GetMpSegmentation(digit->DetElemId(),AliMp::GetCathodType(digit->Cathode()));
176 AliMpPad pad = seg->PadByIndices(AliMpIntPair(digit->PadX(), digit->PadY()));
177
178 // calibration parameters
179 AliMUONVCalibParam* ped = calibData->Pedestals(digit->DetElemId(), digit->ManuId());
180 AliMUONVCalibParam* gain = calibData->Gains(digit->DetElemId(), digit->ManuId());
181 Int_t manuChannel = digit->ManuChannel();
182
183 // fill pad info
184 padInfo.SetPadId(digit->GetUniqueID());
185 padInfo.SetPadXY(pad.Position().X(), pad.Position().Y());
186 padInfo.SetPadDimXY(pad.Dimensions().X(), pad.Dimensions().Y());
187 padInfo.SetPadCharge((Double_t)digit->Charge());
188 padInfo.SetPadADC(digit->ADC());
189 padInfo.SetSaturated(digit->IsSaturated());
190 padInfo.SetCalibrated(digit->IsCalibrated());
191 padInfo.SetPedestal(ped->ValueAsFloatFast(manuChannel,0), ped->ValueAsFloatFast(manuChannel,1));
192 padInfo.SetGain(gain->ValueAsFloatFast(manuChannel,0), gain->ValueAsFloatFast(manuChannel,1),
193 gain->ValueAsFloatFast(manuChannel,2), gain->ValueAsFloatFast(manuChannel,3));
194
195 clusterInfo->AddPad(padInfo);
196 }
197
198 // fill cluster info tree
199 clusterInfoTree->Fill();
200
201 trackParam = static_cast<AliMUONTrackParam*>(newTrack->GetTrackParamAtCluster()->After(trackParam));
202 }
203
204 }
205
206 // free memory
207 delete newTrackStore;
208 }
209
210 // ...timer stop
211 timer.Stop();
212 printf("Writing Tree\n");
213 // write output tree
214 clusterInfoFile->cd();
215 clusterInfoTree->Write();
216 printf("Deleting Tree\n");
217 delete clusterInfoTree;
218 printf("Closing File\n");
219 clusterInfoFile->Close();
220
221 // free memory
222 printf("Deleting calibData\n");
223 delete calibData;
224 printf("Deleting clusterInfo\n");
225 delete clusterInfo;
226 printf("Closing esdFile\n");
227 esdFile->Close();
228 printf("Deleting esd\n");
229 delete esd;
230 // delete padInfo;
231 cout<<endl<<"time to fill cluster/track info: R:"<<timer.RealTime()<<" C:"<<timer.CpuTime()<<endl<<endl;
232}
233
234//-----------------------------------------------------------------------
235void Prepare()
236{
237 /// Set the geometry, the magnetic field, the mapping and the reconstruction parameters
238
239 // Import TGeo geometry (needed by AliMUONTrackExtrap::ExtrapToVertex)
240 if (!gGeoManager) {
241 AliGeomManager::LoadGeometry("geometry.root");
242 if (!gGeoManager) {
243 Error("MUONRefit", "getting geometry from file %s failed", "generated/galice.root");
244 return;
245 }
246 }
247
248 // set mag field
249 printf("Loading field map...\n");
250 AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 1, 1., 10., AliMagFMaps::k5kG);
251 AliTracker::SetFieldMap(field, kFALSE);
252 AliMUONTrackExtrap::SetField(AliTracker::GetFieldMap());
253
254 // Load mapping
255 AliCDBManager* man = AliCDBManager::Instance();
256 man->SetDefaultStorage("local://$ALICE_ROOT");
257 man->SetRun(0);
258 if ( ! AliMpCDB::LoadDDLStore() ) {
259 Error("MUONRefit","Could not access mapping from OCDB !");
260 exit(-1);
261 }
262
263 // set reconstruction parameters
264 AliMUONRecoParam *muonRecoParam = AliMUONRecoParam::GetLowFluxParam();
265 muonRecoParam->Print("FULL");
266 AliRecoParam::Instance()->RegisterRecoParam(muonRecoParam);
267
268}
269
270//-----------------------------------------------------------------------
271TTree* GetESDTree(TFile *esdFile)
272{
273 /// Check that the file is properly open
274 /// Return pointer to the ESD Tree
275
276 if (!esdFile || !esdFile->IsOpen()) {
277 Error("GetESDTree", "opening ESD file failed");
278 exit(-1);
279 }
280
281 TTree* tree = (TTree*) esdFile->Get("esdTree");
282 if (!tree) {
283 Error("GetESDTree", "no ESD tree found");
284 exit(-1);
285 }
286
287 return tree;
288
289}
290