]>
Commit | Line | Data |
---|---|---|
fcefbac4 | 1 | /************************************************************************** |
41c52850 | 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 | **************************************************************************/ | |
fcefbac4 | 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> | |
ce8cd162 | 27 | #include <TObjArray.h> |
fcefbac4 | 28 | #include <TTree.h> |
29 | #include <TString.h> | |
30 | #include <Riostream.h> | |
fcefbac4 | 31 | #include <TRandom.h> |
32 | #include <TROOT.h> | |
41c52850 | 33 | #include <TMath.h> |
fcefbac4 | 34 | |
35 | // STEER includes | |
fcefbac4 | 36 | #include "AliESDEvent.h" |
fcefbac4 | 37 | #include "AliRecoParam.h" |
38 | #include "AliCDBManager.h" | |
41c52850 | 39 | #include "AliRunLoader.h" |
40 | #include "AliLoader.h" | |
a99c3449 | 41 | #include "AliHeader.h" |
fcefbac4 | 42 | |
43 | // MUON includes | |
41c52850 | 44 | #include "AliMpConstants.h" |
fcefbac4 | 45 | #include "AliMpSegmentation.h" |
46 | #include "AliMpVSegmentation.h" | |
47 | #include "AliMpPad.h" | |
a99c3449 | 48 | #include "AliMUONCDB.h" |
fcefbac4 | 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" | |
41c52850 | 58 | #include "AliMUONVClusterStore.h" |
fcefbac4 | 59 | #include "AliMUONTrack.h" |
fcefbac4 | 60 | #include "AliMUONTrackParam.h" |
fcefbac4 | 61 | #endif |
62 | ||
63 | const Int_t printLevel = 1; | |
64 | ||
fcefbac4 | 65 | TTree* GetESDTree(TFile *esdFile); |
41c52850 | 66 | UInt_t buildClusterMap(AliMUONTrack &track); |
fcefbac4 | 67 | |
68 | //----------------------------------------------------------------------- | |
a99c3449 | 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") | |
fcefbac4 | 71 | { |
41c52850 | 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 | } | |
fcefbac4 | 90 | |
91 | AliMUONClusterInfo* clusterInfo = new AliMUONClusterInfo(); | |
92 | AliMUONPadInfo padInfo; | |
93 | AliMUONCalibrationData* calibData = 0x0; | |
8468d8c6 | 94 | AliMUONESDInterface esdInterface; |
41c52850 | 95 | AliMUONVClusterStore* clusterStore = 0x0; |
96 | AliMUONVDigitStore* digitStore = 0x0; | |
fcefbac4 | 97 | |
fcefbac4 | 98 | // open the ESD file and tree and connect the ESD event |
41c52850 | 99 | TFile* esdFile = 0x0; |
100 | TTree* esdTree = 0x0; | |
101 | AliESDEvent* esd = 0x0; | |
a99c3449 | 102 | Int_t runNumber = -1; |
41c52850 | 103 | if (useESD) { |
104 | esdFile = TFile::Open(esdFileName); | |
105 | esdTree = GetESDTree(esdFile); | |
106 | esd = new AliESDEvent(); | |
107 | esd->ReadFromTree(esdTree); | |
a99c3449 | 108 | if (esdTree->GetEvent(0) <= 0) { |
109 | Error("MUONClusterInfo", "no ESD object found for event 0"); | |
110 | return; | |
111 | } | |
112 | runNumber = esd->GetRunNumber(); | |
41c52850 | 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"); | |
a99c3449 | 123 | rl->LoadHeader(); |
124 | if (runNumber < 0) runNumber = rl->GetHeader()->GetRun(); | |
41c52850 | 125 | } |
fcefbac4 | 126 | |
a99c3449 | 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 | ||
fcefbac4 | 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 | ||
41c52850 | 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 | } | |
fcefbac4 | 160 | for (Int_t iEvent = 0; iEvent < nevents; iEvent++) { |
161 | ||
162 | //----------------------------------------------// | |
163 | // -------------- process event --------------- // | |
164 | //----------------------------------------------// | |
165 | // get the ESD of current event | |
41c52850 | 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(); | |
fcefbac4 | 176 | } |
fcefbac4 | 177 | |
41c52850 | 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 | } | |
fcefbac4 | 200 | |
fcefbac4 | 201 | //----------------------------------------------// |
202 | // ------------- fill cluster info ------------ // | |
203 | //----------------------------------------------// | |
41c52850 | 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) { | |
fcefbac4 | 294 | |
41c52850 | 295 | TIter nextCluster(clusterStore->CreateIterator()); |
296 | AliMUONVCluster *cluster; | |
297 | while ( ( cluster = static_cast<AliMUONVCluster*>(nextCluster()) ) ) { | |
298 | ||
fcefbac4 | 299 | clusterInfo->Clear("C"); |
300 | ||
301 | // fill cluster info | |
41c52850 | 302 | clusterInfo->SetRunId(rl->GetRunNumber()); |
fcefbac4 | 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 | ||
41c52850 | 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 | ||
fcefbac4 | 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()-> | |
41c52850 | 329 | GetMpSegmentation(digit->DetElemId(),AliMp::GetCathodType(digit->Cathode())); |
168e9c4d | 330 | AliMpPad pad = seg->PadByIndices(digit->PadX(), digit->PadY()); |
fcefbac4 | 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(); | |
41c52850 | 336 | Int_t planeType = 0; |
337 | if ( digit->ManuId() & AliMpConstants::ManuMask(AliMp::kNonBendingPlane)) { | |
338 | planeType = 1; | |
339 | } | |
fcefbac4 | 340 | |
341 | // fill pad info | |
342 | padInfo.SetPadId(digit->GetUniqueID()); | |
41c52850 | 343 | padInfo.SetPadPlaneType(planeType); |
68828de2 | 344 | padInfo.SetPadXY(pad.GetPositionX(), pad.GetPositionY()); |
345 | padInfo.SetPadDimXY(pad.GetDimensionX(), pad.GetDimensionY()); | |
fcefbac4 | 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 | ||
fcefbac4 | 360 | } |
361 | ||
41c52850 | 362 | delete digitStore; |
363 | delete clusterStore; | |
364 | ||
fcefbac4 | 365 | } |
366 | ||
fcefbac4 | 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; | |
41c52850 | 385 | if (useRecPoints) { |
386 | MUONLoader->UnloadDigits(); | |
387 | MUONLoader->UnloadRecPoints(); | |
a99c3449 | 388 | rl->UnloadHeader(); |
41c52850 | 389 | delete rl; |
390 | } | |
391 | if (useESD) { | |
392 | esdFile->Close(); | |
393 | delete esd; | |
394 | } | |
fcefbac4 | 395 | cout<<endl<<"time to fill cluster/track info: R:"<<timer.RealTime()<<" C:"<<timer.CpuTime()<<endl<<endl; |
396 | } | |
397 | ||
398 | //----------------------------------------------------------------------- | |
fcefbac4 | 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 | ||
41c52850 | 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 |