]>
Commit | Line | Data |
---|---|---|
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 <TRandom.h> | |
32 | #include <TROOT.h> | |
33 | #include <TMath.h> | |
34 | ||
35 | // STEER includes | |
36 | #include "AliESDEvent.h" | |
37 | #include "AliRecoParam.h" | |
38 | #include "AliCDBManager.h" | |
39 | #include "AliRunLoader.h" | |
40 | #include "AliLoader.h" | |
41 | #include "AliHeader.h" | |
42 | ||
43 | // MUON includes | |
44 | #include "AliMpConstants.h" | |
45 | #include "AliMpSegmentation.h" | |
46 | #include "AliMpVSegmentation.h" | |
47 | #include "AliMpPad.h" | |
48 | #include "AliMUONCDB.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 "AliMUONVDigit.h" | |
56 | #include "AliMUONVDigitStore.h" | |
57 | #include "AliMUONVCluster.h" | |
58 | #include "AliMUONVClusterStore.h" | |
59 | #include "AliMUONTrack.h" | |
60 | #include "AliMUONTrackParam.h" | |
61 | #endif | |
62 | ||
63 | const Int_t printLevel = 1; | |
64 | ||
65 | TTree* GetESDTree(TFile *esdFile); | |
66 | UInt_t buildClusterMap(AliMUONTrack &track); | |
67 | ||
68 | //----------------------------------------------------------------------- | |
69 | void MUONClusterInfo(Int_t nevents = -1, const char* esdFileName = "AliESDs.root", const char* inFileName = "galice.root", | |
70 | const TString ocdbPath = "local://$ALICE_ROOT/OCDB", const char* outFileName = "clusterInfo.root") | |
71 | { | |
72 | /// 1) if (esdFileName != "") | |
73 | /// loop over ESD event and fill AliMUONClusterInfo object with cluster + corresponding track parameters; | |
74 | /// 2) if (inFileName != "") | |
75 | /// loop over RecPoints and fill AliMUONClusterInfo object with cluster not attached to a track; | |
76 | /// 3) write results in a new root file. | |
77 | /// | |
78 | /// ******************************************* WARNING ******************************************* /// | |
79 | /// track parameters at each cluster are recomputed by the interface using Kalman filter + Smoother /// | |
80 | /// (It can be changed by resetting the tracker in the interface with a new recoParam object) /// | |
81 | /// and the magnetic field set in the function prepare() /// | |
82 | /// ******************************************* WARNING ******************************************* /// | |
83 | ||
84 | Bool_t useESD = (strcmp(esdFileName,"")); | |
85 | Bool_t useRecPoints = (strcmp(inFileName,"")); | |
86 | if (!useESD && !useRecPoints) { | |
87 | Error("MUONClusterInfo","you must provide ESD and/or galice root file(s)"); | |
88 | return; | |
89 | } | |
90 | ||
91 | AliMUONClusterInfo* clusterInfo = new AliMUONClusterInfo(); | |
92 | AliMUONPadInfo padInfo; | |
93 | AliMUONCalibrationData* calibData = 0x0; | |
94 | AliMUONESDInterface esdInterface; | |
95 | AliMUONVClusterStore* clusterStore = 0x0; | |
96 | AliMUONVDigitStore* digitStore = 0x0; | |
97 | ||
98 | // open the ESD file and tree and connect the ESD event | |
99 | TFile* esdFile = 0x0; | |
100 | TTree* esdTree = 0x0; | |
101 | AliESDEvent* esd = 0x0; | |
102 | Int_t runNumber = -1; | |
103 | if (useESD) { | |
104 | esdFile = TFile::Open(esdFileName); | |
105 | esdTree = GetESDTree(esdFile); | |
106 | esd = new AliESDEvent(); | |
107 | esd->ReadFromTree(esdTree); | |
108 | if (esdTree->GetEvent(0) <= 0) { | |
109 | Error("MUONClusterInfo", "no ESD object found for event 0"); | |
110 | return; | |
111 | } | |
112 | runNumber = esd->GetRunNumber(); | |
113 | } | |
114 | ||
115 | // get the cluster from RecPoints | |
116 | AliRunLoader * rl = 0x0; | |
117 | AliLoader* MUONLoader = 0x0; | |
118 | if (useRecPoints) { | |
119 | rl = AliRunLoader::Open(inFileName,"MUONLoader"); | |
120 | MUONLoader = rl->GetDetectorLoader("MUON"); | |
121 | MUONLoader->LoadRecPoints("READ"); | |
122 | MUONLoader->LoadDigits("READ"); | |
123 | rl->LoadHeader(); | |
124 | if (runNumber < 0) runNumber = rl->GetHeader()->GetRun(); | |
125 | } | |
126 | ||
127 | // load necessary data from OCDB | |
128 | AliCDBManager::Instance()->SetDefaultStorage(ocdbPath); | |
129 | AliCDBManager::Instance()->SetRun(runNumber); | |
130 | if (!AliMUONCDB::LoadField()) return; | |
131 | if (!AliMUONCDB::LoadMapping(kTRUE)) return; | |
132 | AliMUONRecoParam* recoParam = AliMUONCDB::LoadRecoParam(); | |
133 | if (!recoParam) return; | |
134 | ||
135 | // reset tracker for track restoring initial track parameters at cluster | |
136 | AliMUONESDInterface::ResetTracker(recoParam); | |
137 | ||
138 | // prepare access to calibration data | |
139 | calibData = new AliMUONCalibrationData(runNumber); | |
140 | ||
141 | // prepare the output tree | |
142 | gROOT->cd(); | |
143 | TFile* clusterInfoFile = TFile::Open(outFileName, "RECREATE"); | |
144 | clusterInfoFile->SetCompressionLevel(1); | |
145 | ||
146 | TTree* clusterInfoTree = new TTree("clusterInfoTree","clusterInfoTree"); | |
147 | clusterInfoTree->Branch("clusterInfo", &clusterInfo, 32000, 99); | |
148 | ||
149 | // timer start... | |
150 | TStopwatch timer; | |
151 | ||
152 | // Loop over events | |
153 | if (useESD) { | |
154 | if (nevents > 0) nevents = TMath::Min(nevents,(Int_t)esdTree->GetEntries()); | |
155 | else nevents = (Int_t)esdTree->GetEntries(); | |
156 | } else { | |
157 | if (nevents > 0) nevents = TMath::Min(nevents,(Int_t)rl->GetNumberOfEvents()); | |
158 | else nevents = (Int_t)rl->GetNumberOfEvents(); | |
159 | } | |
160 | for (Int_t iEvent = 0; iEvent < nevents; iEvent++) { | |
161 | ||
162 | //----------------------------------------------// | |
163 | // -------------- process event --------------- // | |
164 | //----------------------------------------------// | |
165 | // get the ESD of current event | |
166 | if (useESD) { | |
167 | esdTree->GetEvent(iEvent); | |
168 | if (!esd) { | |
169 | Error("MUONClusterInfo", "no ESD object found for event %d", iEvent); | |
170 | return; | |
171 | } | |
172 | // load the current esd event | |
173 | esdInterface.LoadEvent(*esd); | |
174 | // get digit store | |
175 | if (!useRecPoints) digitStore = esdInterface.GetDigits(); | |
176 | } | |
177 | ||
178 | if (useRecPoints) { | |
179 | if (!(rl->GetEvent(iEvent) == 0)) { | |
180 | Error("MUONClusterInfo", "unable to load event %d", iEvent); | |
181 | return; | |
182 | } | |
183 | // get the clusters of current event | |
184 | TTree* treeR = MUONLoader->TreeR(); | |
185 | clusterStore = AliMUONVClusterStore::Create(*treeR); | |
186 | if ( clusterStore != 0x0 ) { | |
187 | clusterStore->Clear(); | |
188 | clusterStore->Connect(*treeR); | |
189 | treeR->GetEvent(0); | |
190 | } | |
191 | // get the digits of current event | |
192 | TTree* treeD = MUONLoader->TreeD(); | |
193 | digitStore = AliMUONVDigitStore::Create(*treeD); | |
194 | if ( digitStore != 0x0 ) { | |
195 | digitStore->Clear(); | |
196 | digitStore->Connect(*treeD); | |
197 | treeD->GetEvent(0); | |
198 | } | |
199 | } | |
200 | ||
201 | //----------------------------------------------// | |
202 | // ------------- fill cluster info ------------ // | |
203 | //----------------------------------------------// | |
204 | // --- loop over the refitted tracks --- | |
205 | if (useESD) { | |
206 | ||
207 | TIter nextTrack(esdInterface.CreateTrackIterator()); | |
208 | AliMUONTrack* track; | |
209 | while ((track = static_cast<AliMUONTrack*>(nextTrack()))) { | |
210 | ||
211 | UInt_t muonClusterMap = buildClusterMap(*track); | |
212 | ||
213 | // loop over clusters | |
214 | AliMUONTrackParam* trackParam = static_cast<AliMUONTrackParam*>(track->GetTrackParamAtCluster()->First()); | |
215 | while (trackParam) { | |
216 | clusterInfo->Clear("C"); | |
217 | ||
218 | // fill cluster info | |
219 | AliMUONVCluster* cluster = trackParam->GetClusterPtr(); | |
220 | clusterInfo->SetRunId(esd->GetRunNumber()); | |
221 | clusterInfo->SetEventId(iEvent); | |
222 | clusterInfo->SetZ(cluster->GetZ()); | |
223 | clusterInfo->SetClusterId(cluster->GetUniqueID()); | |
224 | clusterInfo->SetClusterXY(cluster->GetX(), cluster->GetY()); | |
225 | clusterInfo->SetClusterXYErr(cluster->GetErrX(), cluster->GetErrY()); | |
226 | clusterInfo->SetClusterChi2(cluster->GetChi2()); | |
227 | clusterInfo->SetClusterCharge(cluster->GetCharge()); | |
228 | ||
229 | // fill track info | |
230 | clusterInfo->SetTrackId(track->GetUniqueID()); | |
231 | clusterInfo->SetTrackXY(trackParam->GetNonBendingCoor(), trackParam->GetBendingCoor()); | |
232 | clusterInfo->SetTrackThetaXY(TMath::ATan(trackParam->GetNonBendingSlope()), TMath::ATan(trackParam->GetBendingSlope())); | |
233 | clusterInfo->SetTrackP(trackParam->P()); | |
234 | const TMatrixD paramCov = trackParam->GetCovariances(); | |
235 | clusterInfo->SetTrackXYErr(TMath::Sqrt(paramCov(0,0)), TMath::Sqrt(paramCov(2,2))); | |
236 | clusterInfo->SetTrackChi2(track->GetNormalizedChi2()); | |
237 | clusterInfo->SetTrackCharge((Short_t)trackParam->GetCharge()); | |
238 | clusterInfo->SetTrackNHits(track->GetNClusters()); | |
239 | clusterInfo->SetTrackChamberHitMap(muonClusterMap); | |
240 | ||
241 | // fill pad info if available | |
242 | for (Int_t i=0; i<cluster->GetNDigits(); i++) { | |
243 | AliMUONVDigit* digit = digitStore->FindObject(cluster->GetDigitId(i)); | |
244 | if (!digit) continue; | |
245 | ||
246 | // pad location | |
247 | const AliMpVSegmentation* seg = AliMpSegmentation::Instance()-> | |
248 | GetMpSegmentation(digit->DetElemId(),AliMp::GetCathodType(digit->Cathode())); | |
249 | AliMpPad pad = seg->PadByIndices(digit->PadX(), digit->PadY()); | |
250 | ||
251 | // calibration parameters | |
252 | AliMUONVCalibParam* ped = calibData->Pedestals(digit->DetElemId(), digit->ManuId()); | |
253 | AliMUONVCalibParam* gain = calibData->Gains(digit->DetElemId(), digit->ManuId()); | |
254 | Int_t manuChannel = digit->ManuChannel(); | |
255 | Int_t planeType = 0; | |
256 | if ( digit->ManuId() & AliMpConstants::ManuMask(AliMp::kNonBendingPlane)) { | |
257 | planeType = 1; | |
258 | } | |
259 | ||
260 | // fill pad info | |
261 | padInfo.SetPadId(digit->GetUniqueID()); | |
262 | padInfo.SetPadPlaneType(planeType); | |
263 | padInfo.SetPadXY(pad.GetPositionX(), pad.GetPositionY()); | |
264 | padInfo.SetPadDimXY(pad.GetDimensionX(), pad.GetDimensionY()); | |
265 | padInfo.SetPadCharge((Double_t)digit->Charge()); | |
266 | padInfo.SetPadADC(digit->ADC()); | |
267 | padInfo.SetSaturated(digit->IsSaturated()); | |
268 | padInfo.SetCalibrated(digit->IsCalibrated()); | |
269 | padInfo.SetPedestal(ped->ValueAsFloatFast(manuChannel,0), ped->ValueAsFloatFast(manuChannel,1)); | |
270 | padInfo.SetGain(gain->ValueAsFloatFast(manuChannel,0), gain->ValueAsFloatFast(manuChannel,1), | |
271 | gain->ValueAsFloatFast(manuChannel,2), gain->ValueAsFloatFast(manuChannel,3)); | |
272 | ||
273 | clusterInfo->AddPad(padInfo); | |
274 | } | |
275 | ||
276 | // remove clusters attached to a track | |
277 | if (useRecPoints) { | |
278 | AliMUONVCluster* cl = clusterStore->FindObject(cluster->GetUniqueID()); | |
279 | if (cl) clusterStore->Remove(*cl); | |
280 | } | |
281 | ||
282 | // fill cluster info tree | |
283 | clusterInfoTree->Fill(); | |
284 | ||
285 | trackParam = static_cast<AliMUONTrackParam*>(track->GetTrackParamAtCluster()->After(trackParam)); | |
286 | } | |
287 | ||
288 | } | |
289 | ||
290 | } | |
291 | ||
292 | // --- loop over clusters not attached to a track --- | |
293 | if (useRecPoints) { | |
294 | ||
295 | TIter nextCluster(clusterStore->CreateIterator()); | |
296 | AliMUONVCluster *cluster; | |
297 | while ( ( cluster = static_cast<AliMUONVCluster*>(nextCluster()) ) ) { | |
298 | ||
299 | clusterInfo->Clear("C"); | |
300 | ||
301 | // fill cluster info | |
302 | clusterInfo->SetRunId(rl->GetRunNumber()); | |
303 | clusterInfo->SetEventId(iEvent); | |
304 | clusterInfo->SetZ(cluster->GetZ()); | |
305 | clusterInfo->SetClusterId(cluster->GetUniqueID()); | |
306 | clusterInfo->SetClusterXY(cluster->GetX(), cluster->GetY()); | |
307 | clusterInfo->SetClusterXYErr(cluster->GetErrX(), cluster->GetErrY()); | |
308 | clusterInfo->SetClusterChi2(cluster->GetChi2()); | |
309 | clusterInfo->SetClusterCharge(cluster->GetCharge()); | |
310 | ||
311 | // fill dummy track info | |
312 | clusterInfo->SetTrackId(0); | |
313 | clusterInfo->SetTrackXY(0.,0.); | |
314 | clusterInfo->SetTrackThetaXY(0.,0.); | |
315 | clusterInfo->SetTrackP(0.); | |
316 | clusterInfo->SetTrackXYErr(0.,0.); | |
317 | clusterInfo->SetTrackChi2(0.); | |
318 | clusterInfo->SetTrackCharge(0); | |
319 | clusterInfo->SetTrackNHits(0); | |
320 | clusterInfo->SetTrackChamberHitMap(0); | |
321 | ||
322 | // fill pad info if available | |
323 | for (Int_t i=0; i<cluster->GetNDigits(); i++) { | |
324 | AliMUONVDigit* digit = digitStore->FindObject(cluster->GetDigitId(i)); | |
325 | if (!digit) continue; | |
326 | ||
327 | // pad location | |
328 | const AliMpVSegmentation* seg = AliMpSegmentation::Instance()-> | |
329 | GetMpSegmentation(digit->DetElemId(),AliMp::GetCathodType(digit->Cathode())); | |
330 | AliMpPad pad = seg->PadByIndices(digit->PadX(), digit->PadY()); | |
331 | ||
332 | // calibration parameters | |
333 | AliMUONVCalibParam* ped = calibData->Pedestals(digit->DetElemId(), digit->ManuId()); | |
334 | AliMUONVCalibParam* gain = calibData->Gains(digit->DetElemId(), digit->ManuId()); | |
335 | Int_t manuChannel = digit->ManuChannel(); | |
336 | Int_t planeType = 0; | |
337 | if ( digit->ManuId() & AliMpConstants::ManuMask(AliMp::kNonBendingPlane)) { | |
338 | planeType = 1; | |
339 | } | |
340 | ||
341 | // fill pad info | |
342 | padInfo.SetPadId(digit->GetUniqueID()); | |
343 | padInfo.SetPadPlaneType(planeType); | |
344 | padInfo.SetPadXY(pad.GetPositionX(), pad.GetPositionY()); | |
345 | padInfo.SetPadDimXY(pad.GetDimensionX(), pad.GetDimensionY()); | |
346 | padInfo.SetPadCharge((Double_t)digit->Charge()); | |
347 | padInfo.SetPadADC(digit->ADC()); | |
348 | padInfo.SetSaturated(digit->IsSaturated()); | |
349 | padInfo.SetCalibrated(digit->IsCalibrated()); | |
350 | padInfo.SetPedestal(ped->ValueAsFloatFast(manuChannel,0), ped->ValueAsFloatFast(manuChannel,1)); | |
351 | padInfo.SetGain(gain->ValueAsFloatFast(manuChannel,0), gain->ValueAsFloatFast(manuChannel,1), | |
352 | gain->ValueAsFloatFast(manuChannel,2), gain->ValueAsFloatFast(manuChannel,3)); | |
353 | ||
354 | clusterInfo->AddPad(padInfo); | |
355 | } | |
356 | ||
357 | // fill cluster info tree | |
358 | clusterInfoTree->Fill(); | |
359 | ||
360 | } | |
361 | ||
362 | delete digitStore; | |
363 | delete clusterStore; | |
364 | ||
365 | } | |
366 | ||
367 | } | |
368 | ||
369 | // ...timer stop | |
370 | timer.Stop(); | |
371 | printf("Writing Tree\n"); | |
372 | // write output tree | |
373 | clusterInfoFile->cd(); | |
374 | clusterInfoTree->Write(); | |
375 | printf("Deleting Tree\n"); | |
376 | delete clusterInfoTree; | |
377 | printf("Closing File\n"); | |
378 | clusterInfoFile->Close(); | |
379 | ||
380 | // free memory | |
381 | printf("Deleting calibData\n"); | |
382 | delete calibData; | |
383 | printf("Deleting clusterInfo\n"); | |
384 | delete clusterInfo; | |
385 | if (useRecPoints) { | |
386 | MUONLoader->UnloadDigits(); | |
387 | MUONLoader->UnloadRecPoints(); | |
388 | rl->UnloadHeader(); | |
389 | delete rl; | |
390 | } | |
391 | if (useESD) { | |
392 | esdFile->Close(); | |
393 | delete esd; | |
394 | } | |
395 | cout<<endl<<"time to fill cluster/track info: R:"<<timer.RealTime()<<" C:"<<timer.CpuTime()<<endl<<endl; | |
396 | } | |
397 | ||
398 | //----------------------------------------------------------------------- | |
399 | TTree* GetESDTree(TFile *esdFile) | |
400 | { | |
401 | /// Check that the file is properly open | |
402 | /// Return pointer to the ESD Tree | |
403 | ||
404 | if (!esdFile || !esdFile->IsOpen()) { | |
405 | Error("GetESDTree", "opening ESD file failed"); | |
406 | exit(-1); | |
407 | } | |
408 | ||
409 | TTree* tree = (TTree*) esdFile->Get("esdTree"); | |
410 | if (!tree) { | |
411 | Error("GetESDTree", "no ESD tree found"); | |
412 | exit(-1); | |
413 | } | |
414 | ||
415 | return tree; | |
416 | ||
417 | } | |
418 | ||
419 | ||
420 | //----------------------------------------------------------------------- | |
421 | UInt_t buildClusterMap(AliMUONTrack &track) | |
422 | { | |
423 | /// Build the map of clusters in tracking chambers | |
424 | ||
425 | UInt_t muonClusterMap = 0; | |
426 | ||
427 | AliMUONTrackParam* trackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->First()); | |
428 | while (trackParam) { | |
429 | ||
430 | muonClusterMap |= BIT(trackParam->GetClusterPtr()->GetChamberId()); | |
431 | ||
432 | trackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->After(trackParam)); | |
433 | } | |
434 | ||
435 | return muonClusterMap; | |
436 | ||
437 | } | |
438 |