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