]>
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(); | |
bef83d1e | 193 | if (treeD) { |
194 | digitStore = AliMUONVDigitStore::Create(*treeD); | |
195 | if ( digitStore != 0x0 ) { | |
196 | digitStore->Clear(); | |
197 | digitStore->Connect(*treeD); | |
198 | treeD->GetEvent(0); | |
199 | } | |
41c52850 | 200 | } |
201 | } | |
fcefbac4 | 202 | |
fcefbac4 | 203 | //----------------------------------------------// |
204 | // ------------- fill cluster info ------------ // | |
205 | //----------------------------------------------// | |
41c52850 | 206 | // --- loop over the refitted tracks --- |
207 | if (useESD) { | |
208 | ||
209 | TIter nextTrack(esdInterface.CreateTrackIterator()); | |
210 | AliMUONTrack* track; | |
211 | while ((track = static_cast<AliMUONTrack*>(nextTrack()))) { | |
212 | ||
213 | UInt_t muonClusterMap = buildClusterMap(*track); | |
214 | ||
215 | // loop over clusters | |
216 | AliMUONTrackParam* trackParam = static_cast<AliMUONTrackParam*>(track->GetTrackParamAtCluster()->First()); | |
217 | while (trackParam) { | |
218 | clusterInfo->Clear("C"); | |
219 | ||
220 | // fill cluster info | |
221 | AliMUONVCluster* cluster = trackParam->GetClusterPtr(); | |
222 | clusterInfo->SetRunId(esd->GetRunNumber()); | |
223 | clusterInfo->SetEventId(iEvent); | |
224 | clusterInfo->SetZ(cluster->GetZ()); | |
225 | clusterInfo->SetClusterId(cluster->GetUniqueID()); | |
226 | clusterInfo->SetClusterXY(cluster->GetX(), cluster->GetY()); | |
227 | clusterInfo->SetClusterXYErr(cluster->GetErrX(), cluster->GetErrY()); | |
228 | clusterInfo->SetClusterChi2(cluster->GetChi2()); | |
229 | clusterInfo->SetClusterCharge(cluster->GetCharge()); | |
230 | ||
231 | // fill track info | |
232 | clusterInfo->SetTrackId(track->GetUniqueID()); | |
233 | clusterInfo->SetTrackXY(trackParam->GetNonBendingCoor(), trackParam->GetBendingCoor()); | |
234 | clusterInfo->SetTrackThetaXY(TMath::ATan(trackParam->GetNonBendingSlope()), TMath::ATan(trackParam->GetBendingSlope())); | |
235 | clusterInfo->SetTrackP(trackParam->P()); | |
236 | const TMatrixD paramCov = trackParam->GetCovariances(); | |
237 | clusterInfo->SetTrackXYErr(TMath::Sqrt(paramCov(0,0)), TMath::Sqrt(paramCov(2,2))); | |
238 | clusterInfo->SetTrackChi2(track->GetNormalizedChi2()); | |
239 | clusterInfo->SetTrackCharge((Short_t)trackParam->GetCharge()); | |
240 | clusterInfo->SetTrackNHits(track->GetNClusters()); | |
241 | clusterInfo->SetTrackChamberHitMap(muonClusterMap); | |
242 | ||
243 | // fill pad info if available | |
bef83d1e | 244 | if (digitStore) for (Int_t i=0; i<cluster->GetNDigits(); i++) { |
41c52850 | 245 | AliMUONVDigit* digit = digitStore->FindObject(cluster->GetDigitId(i)); |
246 | if (!digit) continue; | |
247 | ||
248 | // pad location | |
249 | const AliMpVSegmentation* seg = AliMpSegmentation::Instance()-> | |
250 | GetMpSegmentation(digit->DetElemId(),AliMp::GetCathodType(digit->Cathode())); | |
251 | AliMpPad pad = seg->PadByIndices(digit->PadX(), digit->PadY()); | |
252 | ||
253 | // calibration parameters | |
254 | AliMUONVCalibParam* ped = calibData->Pedestals(digit->DetElemId(), digit->ManuId()); | |
255 | AliMUONVCalibParam* gain = calibData->Gains(digit->DetElemId(), digit->ManuId()); | |
256 | Int_t manuChannel = digit->ManuChannel(); | |
257 | Int_t planeType = 0; | |
258 | if ( digit->ManuId() & AliMpConstants::ManuMask(AliMp::kNonBendingPlane)) { | |
259 | planeType = 1; | |
260 | } | |
261 | ||
262 | // fill pad info | |
263 | padInfo.SetPadId(digit->GetUniqueID()); | |
264 | padInfo.SetPadPlaneType(planeType); | |
265 | padInfo.SetPadXY(pad.GetPositionX(), pad.GetPositionY()); | |
266 | padInfo.SetPadDimXY(pad.GetDimensionX(), pad.GetDimensionY()); | |
267 | padInfo.SetPadCharge((Double_t)digit->Charge()); | |
268 | padInfo.SetPadADC(digit->ADC()); | |
269 | padInfo.SetSaturated(digit->IsSaturated()); | |
270 | padInfo.SetCalibrated(digit->IsCalibrated()); | |
271 | padInfo.SetPedestal(ped->ValueAsFloatFast(manuChannel,0), ped->ValueAsFloatFast(manuChannel,1)); | |
272 | padInfo.SetGain(gain->ValueAsFloatFast(manuChannel,0), gain->ValueAsFloatFast(manuChannel,1), | |
273 | gain->ValueAsFloatFast(manuChannel,2), gain->ValueAsFloatFast(manuChannel,3)); | |
274 | ||
275 | clusterInfo->AddPad(padInfo); | |
276 | } | |
277 | ||
278 | // remove clusters attached to a track | |
279 | if (useRecPoints) { | |
280 | AliMUONVCluster* cl = clusterStore->FindObject(cluster->GetUniqueID()); | |
281 | if (cl) clusterStore->Remove(*cl); | |
282 | } | |
283 | ||
284 | // fill cluster info tree | |
285 | clusterInfoTree->Fill(); | |
286 | ||
287 | trackParam = static_cast<AliMUONTrackParam*>(track->GetTrackParamAtCluster()->After(trackParam)); | |
288 | } | |
289 | ||
290 | } | |
291 | ||
292 | } | |
293 | ||
294 | // --- loop over clusters not attached to a track --- | |
295 | if (useRecPoints) { | |
fcefbac4 | 296 | |
41c52850 | 297 | TIter nextCluster(clusterStore->CreateIterator()); |
298 | AliMUONVCluster *cluster; | |
299 | while ( ( cluster = static_cast<AliMUONVCluster*>(nextCluster()) ) ) { | |
300 | ||
fcefbac4 | 301 | clusterInfo->Clear("C"); |
302 | ||
303 | // fill cluster info | |
41c52850 | 304 | clusterInfo->SetRunId(rl->GetRunNumber()); |
fcefbac4 | 305 | clusterInfo->SetEventId(iEvent); |
306 | clusterInfo->SetZ(cluster->GetZ()); | |
307 | clusterInfo->SetClusterId(cluster->GetUniqueID()); | |
308 | clusterInfo->SetClusterXY(cluster->GetX(), cluster->GetY()); | |
309 | clusterInfo->SetClusterXYErr(cluster->GetErrX(), cluster->GetErrY()); | |
310 | clusterInfo->SetClusterChi2(cluster->GetChi2()); | |
311 | clusterInfo->SetClusterCharge(cluster->GetCharge()); | |
312 | ||
41c52850 | 313 | // fill dummy track info |
314 | clusterInfo->SetTrackId(0); | |
315 | clusterInfo->SetTrackXY(0.,0.); | |
316 | clusterInfo->SetTrackThetaXY(0.,0.); | |
317 | clusterInfo->SetTrackP(0.); | |
318 | clusterInfo->SetTrackXYErr(0.,0.); | |
319 | clusterInfo->SetTrackChi2(0.); | |
320 | clusterInfo->SetTrackCharge(0); | |
321 | clusterInfo->SetTrackNHits(0); | |
322 | clusterInfo->SetTrackChamberHitMap(0); | |
323 | ||
fcefbac4 | 324 | // fill pad info if available |
bef83d1e | 325 | if (digitStore) for (Int_t i=0; i<cluster->GetNDigits(); i++) { |
fcefbac4 | 326 | AliMUONVDigit* digit = digitStore->FindObject(cluster->GetDigitId(i)); |
327 | if (!digit) continue; | |
328 | ||
329 | // pad location | |
330 | const AliMpVSegmentation* seg = AliMpSegmentation::Instance()-> | |
41c52850 | 331 | GetMpSegmentation(digit->DetElemId(),AliMp::GetCathodType(digit->Cathode())); |
168e9c4d | 332 | AliMpPad pad = seg->PadByIndices(digit->PadX(), digit->PadY()); |
fcefbac4 | 333 | |
334 | // calibration parameters | |
335 | AliMUONVCalibParam* ped = calibData->Pedestals(digit->DetElemId(), digit->ManuId()); | |
336 | AliMUONVCalibParam* gain = calibData->Gains(digit->DetElemId(), digit->ManuId()); | |
337 | Int_t manuChannel = digit->ManuChannel(); | |
41c52850 | 338 | Int_t planeType = 0; |
339 | if ( digit->ManuId() & AliMpConstants::ManuMask(AliMp::kNonBendingPlane)) { | |
340 | planeType = 1; | |
341 | } | |
fcefbac4 | 342 | |
343 | // fill pad info | |
344 | padInfo.SetPadId(digit->GetUniqueID()); | |
41c52850 | 345 | padInfo.SetPadPlaneType(planeType); |
68828de2 | 346 | padInfo.SetPadXY(pad.GetPositionX(), pad.GetPositionY()); |
347 | padInfo.SetPadDimXY(pad.GetDimensionX(), pad.GetDimensionY()); | |
fcefbac4 | 348 | padInfo.SetPadCharge((Double_t)digit->Charge()); |
349 | padInfo.SetPadADC(digit->ADC()); | |
350 | padInfo.SetSaturated(digit->IsSaturated()); | |
351 | padInfo.SetCalibrated(digit->IsCalibrated()); | |
352 | padInfo.SetPedestal(ped->ValueAsFloatFast(manuChannel,0), ped->ValueAsFloatFast(manuChannel,1)); | |
353 | padInfo.SetGain(gain->ValueAsFloatFast(manuChannel,0), gain->ValueAsFloatFast(manuChannel,1), | |
354 | gain->ValueAsFloatFast(manuChannel,2), gain->ValueAsFloatFast(manuChannel,3)); | |
355 | ||
356 | clusterInfo->AddPad(padInfo); | |
357 | } | |
358 | ||
359 | // fill cluster info tree | |
360 | clusterInfoTree->Fill(); | |
361 | ||
fcefbac4 | 362 | } |
363 | ||
41c52850 | 364 | delete digitStore; |
365 | delete clusterStore; | |
366 | ||
fcefbac4 | 367 | } |
368 | ||
fcefbac4 | 369 | } |
370 | ||
371 | // ...timer stop | |
372 | timer.Stop(); | |
373 | printf("Writing Tree\n"); | |
374 | // write output tree | |
375 | clusterInfoFile->cd(); | |
376 | clusterInfoTree->Write(); | |
377 | printf("Deleting Tree\n"); | |
378 | delete clusterInfoTree; | |
379 | printf("Closing File\n"); | |
380 | clusterInfoFile->Close(); | |
381 | ||
382 | // free memory | |
383 | printf("Deleting calibData\n"); | |
384 | delete calibData; | |
385 | printf("Deleting clusterInfo\n"); | |
386 | delete clusterInfo; | |
41c52850 | 387 | if (useRecPoints) { |
388 | MUONLoader->UnloadDigits(); | |
389 | MUONLoader->UnloadRecPoints(); | |
a99c3449 | 390 | rl->UnloadHeader(); |
41c52850 | 391 | delete rl; |
392 | } | |
393 | if (useESD) { | |
394 | esdFile->Close(); | |
395 | delete esd; | |
396 | } | |
fcefbac4 | 397 | cout<<endl<<"time to fill cluster/track info: R:"<<timer.RealTime()<<" C:"<<timer.CpuTime()<<endl<<endl; |
398 | } | |
399 | ||
400 | //----------------------------------------------------------------------- | |
fcefbac4 | 401 | TTree* GetESDTree(TFile *esdFile) |
402 | { | |
403 | /// Check that the file is properly open | |
404 | /// Return pointer to the ESD Tree | |
405 | ||
406 | if (!esdFile || !esdFile->IsOpen()) { | |
407 | Error("GetESDTree", "opening ESD file failed"); | |
408 | exit(-1); | |
409 | } | |
410 | ||
411 | TTree* tree = (TTree*) esdFile->Get("esdTree"); | |
412 | if (!tree) { | |
413 | Error("GetESDTree", "no ESD tree found"); | |
414 | exit(-1); | |
415 | } | |
416 | ||
417 | return tree; | |
418 | ||
419 | } | |
420 | ||
41c52850 | 421 | |
422 | //----------------------------------------------------------------------- | |
423 | UInt_t buildClusterMap(AliMUONTrack &track) | |
424 | { | |
425 | /// Build the map of clusters in tracking chambers | |
426 | ||
427 | UInt_t muonClusterMap = 0; | |
428 | ||
429 | AliMUONTrackParam* trackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->First()); | |
430 | while (trackParam) { | |
431 | ||
432 | muonClusterMap |= BIT(trackParam->GetClusterPtr()->GetChamberId()); | |
433 | ||
434 | trackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->After(trackParam)); | |
435 | } | |
436 | ||
437 | return muonClusterMap; | |
438 | ||
439 | } | |
440 |