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