]>
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> | |
27 | #include <TClonesArray.h> | |
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 | |
f7a1cc68 | 36 | #include "AliMagF.h" |
fcefbac4 | 37 | #include "AliTracker.h" |
38 | #include "AliESDEvent.h" | |
fcefbac4 | 39 | #include "AliRecoParam.h" |
40 | #include "AliCDBManager.h" | |
41c52850 | 41 | #include "AliRunLoader.h" |
42 | #include "AliLoader.h" | |
fcefbac4 | 43 | |
44 | // MUON includes | |
41c52850 | 45 | #include "AliMpConstants.h" |
fcefbac4 | 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" | |
fcefbac4 | 56 | #include "AliMUONVDigit.h" |
57 | #include "AliMUONVDigitStore.h" | |
58 | #include "AliMUONVCluster.h" | |
41c52850 | 59 | #include "AliMUONVClusterStore.h" |
fcefbac4 | 60 | #include "AliMUONTrack.h" |
fcefbac4 | 61 | #include "AliMUONTrackParam.h" |
fcefbac4 | 62 | #endif |
63 | ||
64 | const Int_t printLevel = 1; | |
65 | ||
66 | void Prepare(); | |
67 | TTree* GetESDTree(TFile *esdFile); | |
41c52850 | 68 | UInt_t buildClusterMap(AliMUONTrack &track); |
fcefbac4 | 69 | |
70 | //----------------------------------------------------------------------- | |
41c52850 | 71 | void MUONClusterInfo(Int_t nevents = -1, const char* esdFileName = "AliESDs.root", |
72 | const char* inFileName = "galice.root", const char* outFileName = "clusterInfo.root") | |
fcefbac4 | 73 | { |
41c52850 | 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 | } | |
fcefbac4 | 92 | |
93 | AliMUONClusterInfo* clusterInfo = new AliMUONClusterInfo(); | |
94 | AliMUONPadInfo padInfo; | |
95 | AliMUONCalibrationData* calibData = 0x0; | |
8468d8c6 | 96 | AliMUONESDInterface esdInterface; |
41c52850 | 97 | AliMUONVClusterStore* clusterStore = 0x0; |
98 | AliMUONVDigitStore* digitStore = 0x0; | |
fcefbac4 | 99 | |
8468d8c6 | 100 | // prepare the refitting during ESD->MUON conversion |
fcefbac4 | 101 | Prepare(); |
fcefbac4 | 102 | |
103 | // open the ESD file and tree and connect the ESD event | |
41c52850 | 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 | } | |
fcefbac4 | 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 | ||
41c52850 | 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 | } | |
fcefbac4 | 143 | for (Int_t iEvent = 0; iEvent < nevents; iEvent++) { |
144 | ||
145 | //----------------------------------------------// | |
146 | // -------------- process event --------------- // | |
147 | //----------------------------------------------// | |
148 | // get the ESD of current event | |
41c52850 | 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(); | |
fcefbac4 | 159 | } |
fcefbac4 | 160 | |
41c52850 | 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 | } | |
fcefbac4 | 183 | |
41c52850 | 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()); | |
fcefbac4 | 187 | |
fcefbac4 | 188 | //----------------------------------------------// |
189 | // ------------- fill cluster info ------------ // | |
190 | //----------------------------------------------// | |
41c52850 | 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) { | |
fcefbac4 | 281 | |
41c52850 | 282 | TIter nextCluster(clusterStore->CreateIterator()); |
283 | AliMUONVCluster *cluster; | |
284 | while ( ( cluster = static_cast<AliMUONVCluster*>(nextCluster()) ) ) { | |
285 | ||
fcefbac4 | 286 | clusterInfo->Clear("C"); |
287 | ||
288 | // fill cluster info | |
41c52850 | 289 | clusterInfo->SetRunId(rl->GetRunNumber()); |
fcefbac4 | 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 | ||
41c52850 | 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 | ||
fcefbac4 | 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()-> | |
41c52850 | 316 | GetMpSegmentation(digit->DetElemId(),AliMp::GetCathodType(digit->Cathode())); |
168e9c4d | 317 | AliMpPad pad = seg->PadByIndices(digit->PadX(), digit->PadY()); |
fcefbac4 | 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(); | |
41c52850 | 323 | Int_t planeType = 0; |
324 | if ( digit->ManuId() & AliMpConstants::ManuMask(AliMp::kNonBendingPlane)) { | |
325 | planeType = 1; | |
326 | } | |
fcefbac4 | 327 | |
328 | // fill pad info | |
329 | padInfo.SetPadId(digit->GetUniqueID()); | |
41c52850 | 330 | padInfo.SetPadPlaneType(planeType); |
68828de2 | 331 | padInfo.SetPadXY(pad.GetPositionX(), pad.GetPositionY()); |
332 | padInfo.SetPadDimXY(pad.GetDimensionX(), pad.GetDimensionY()); | |
fcefbac4 | 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 | ||
fcefbac4 | 347 | } |
348 | ||
41c52850 | 349 | delete digitStore; |
350 | delete clusterStore; | |
351 | ||
fcefbac4 | 352 | } |
353 | ||
fcefbac4 | 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; | |
41c52850 | 372 | if (useRecPoints) { |
373 | MUONLoader->UnloadDigits(); | |
374 | MUONLoader->UnloadRecPoints(); | |
375 | delete rl; | |
376 | } | |
377 | if (useESD) { | |
378 | esdFile->Close(); | |
379 | delete esd; | |
380 | } | |
fcefbac4 | 381 | cout<<endl<<"time to fill cluster/track info: R:"<<timer.RealTime()<<" C:"<<timer.CpuTime()<<endl<<endl; |
382 | } | |
383 | ||
384 | //----------------------------------------------------------------------- | |
385 | void Prepare() | |
386 | { | |
41c52850 | 387 | /// Set the magnetic field, the mapping and the reconstruction parameters |
fcefbac4 | 388 | |
8468d8c6 | 389 | gRandom->SetSeed(0); |
390 | ||
fcefbac4 | 391 | // set mag field |
f7a1cc68 | 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 | } | |
fcefbac4 | 397 | |
398 | // Load mapping | |
399 | AliCDBManager* man = AliCDBManager::Instance(); | |
162637e4 | 400 | man->SetDefaultStorage("local://$ALICE_ROOT/OCDB"); |
fcefbac4 | 401 | man->SetRun(0); |
402 | if ( ! AliMpCDB::LoadDDLStore() ) { | |
41c52850 | 403 | Error("Prepare","Could not access mapping from OCDB !"); |
fcefbac4 | 404 | exit(-1); |
405 | } | |
406 | ||
ad3c6eda | 407 | // Reset the reconstruction parameters for track refitting if needed |
408 | // (by default will use Kalman filter + Smoother) | |
41c52850 | 409 | // AliMUONRecoParam *muonRecoParam = AliMUONRecoParam::GetLowFluxParam(); |
410 | // AliMUONESDInterface::ResetTracker(muonRecoParam); | |
ad3c6eda | 411 | |
fcefbac4 | 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 | ||
41c52850 | 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 |