1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
19 /// \file MUONClusterInfo.C
20 /// \brief Macro to fill AliMUONClusterInfo objects
22 /// \author Philippe Pillot, SUBATECH
24 #if !defined(__CINT__) || defined(__MAKECINT__)
25 #include <TStopwatch.h>
27 #include <TObjArray.h>
30 #include <Riostream.h>
36 #include "AliESDEvent.h"
37 #include "AliRecoParam.h"
38 #include "AliCDBManager.h"
39 #include "AliRunLoader.h"
40 #include "AliLoader.h"
41 #include "AliHeader.h"
44 #include "AliMpConstants.h"
45 #include "AliMpSegmentation.h"
46 #include "AliMpVSegmentation.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"
63 const Int_t printLevel = 1;
65 TTree* GetESDTree(TFile *esdFile);
66 UInt_t buildClusterMap(AliMUONTrack &track);
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")
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.
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 ******************************************* ///
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)");
91 AliMUONClusterInfo* clusterInfo = new AliMUONClusterInfo();
92 AliMUONPadInfo padInfo;
93 AliMUONCalibrationData* calibData = 0x0;
94 AliMUONESDInterface esdInterface;
95 AliMUONVClusterStore* clusterStore = 0x0;
96 AliMUONVDigitStore* digitStore = 0x0;
98 // open the ESD file and tree and connect the ESD event
100 TTree* esdTree = 0x0;
101 AliESDEvent* esd = 0x0;
102 Int_t runNumber = -1;
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");
112 runNumber = esd->GetRunNumber();
115 // get the cluster from RecPoints
116 AliRunLoader * rl = 0x0;
117 AliLoader* MUONLoader = 0x0;
119 rl = AliRunLoader::Open(inFileName,"MUONLoader");
120 MUONLoader = rl->GetDetectorLoader("MUON");
121 MUONLoader->LoadRecPoints("READ");
122 MUONLoader->LoadDigits("READ");
124 if (runNumber < 0) runNumber = rl->GetHeader()->GetRun();
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;
135 // reset tracker for track restoring initial track parameters at cluster
136 AliMUONESDInterface::ResetTracker(recoParam);
138 // prepare access to calibration data
139 calibData = new AliMUONCalibrationData(runNumber);
141 // prepare the output tree
143 TFile* clusterInfoFile = TFile::Open(outFileName, "RECREATE");
144 clusterInfoFile->SetCompressionLevel(1);
146 TTree* clusterInfoTree = new TTree("clusterInfoTree","clusterInfoTree");
147 clusterInfoTree->Branch("clusterInfo", &clusterInfo, 32000, 99);
154 if (nevents > 0) nevents = TMath::Min(nevents,(Int_t)esdTree->GetEntries());
155 else nevents = (Int_t)esdTree->GetEntries();
157 if (nevents > 0) nevents = TMath::Min(nevents,(Int_t)rl->GetNumberOfEvents());
158 else nevents = (Int_t)rl->GetNumberOfEvents();
160 for (Int_t iEvent = 0; iEvent < nevents; iEvent++) {
162 //----------------------------------------------//
163 // -------------- process event --------------- //
164 //----------------------------------------------//
165 // get the ESD of current event
167 esdTree->GetEvent(iEvent);
169 Error("MUONClusterInfo", "no ESD object found for event %d", iEvent);
172 // load the current esd event
173 esdInterface.LoadEvent(*esd);
175 if (!useRecPoints) digitStore = esdInterface.GetDigits();
179 if (!(rl->GetEvent(iEvent) == 0)) {
180 Error("MUONClusterInfo", "unable to load event %d", iEvent);
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);
191 // get the digits of current event
192 TTree* treeD = MUONLoader->TreeD();
194 digitStore = AliMUONVDigitStore::Create(*treeD);
195 if ( digitStore != 0x0 ) {
197 digitStore->Connect(*treeD);
203 //----------------------------------------------//
204 // ------------- fill cluster info ------------ //
205 //----------------------------------------------//
206 // --- loop over the refitted tracks ---
209 TIter nextTrack(esdInterface.CreateTrackIterator());
211 while ((track = static_cast<AliMUONTrack*>(nextTrack()))) {
213 UInt_t muonClusterMap = buildClusterMap(*track);
215 // loop over clusters
216 AliMUONTrackParam* trackParam = static_cast<AliMUONTrackParam*>(track->GetTrackParamAtCluster()->First());
218 clusterInfo->Clear("C");
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());
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);
243 // fill pad info if available
244 if (digitStore) for (Int_t i=0; i<cluster->GetNDigits(); i++) {
245 AliMUONVDigit* digit = digitStore->FindObject(cluster->GetDigitId(i));
246 if (!digit) continue;
249 const AliMpVSegmentation* seg = AliMpSegmentation::Instance()->
250 GetMpSegmentation(digit->DetElemId(),AliMp::GetCathodType(digit->Cathode()));
251 AliMpPad pad = seg->PadByIndices(digit->PadX(), digit->PadY());
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();
258 if ( digit->ManuId() & AliMpConstants::ManuMask(AliMp::kNonBendingPlane)) {
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));
275 clusterInfo->AddPad(padInfo);
278 // remove clusters attached to a track
280 AliMUONVCluster* cl = clusterStore->FindObject(cluster->GetUniqueID());
281 if (cl) clusterStore->Remove(*cl);
284 // fill cluster info tree
285 clusterInfoTree->Fill();
287 trackParam = static_cast<AliMUONTrackParam*>(track->GetTrackParamAtCluster()->After(trackParam));
294 // --- loop over clusters not attached to a track ---
297 TIter nextCluster(clusterStore->CreateIterator());
298 AliMUONVCluster *cluster;
299 while ( ( cluster = static_cast<AliMUONVCluster*>(nextCluster()) ) ) {
301 clusterInfo->Clear("C");
304 clusterInfo->SetRunId(rl->GetRunNumber());
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());
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);
324 // fill pad info if available
325 if (digitStore) for (Int_t i=0; i<cluster->GetNDigits(); i++) {
326 AliMUONVDigit* digit = digitStore->FindObject(cluster->GetDigitId(i));
327 if (!digit) continue;
330 const AliMpVSegmentation* seg = AliMpSegmentation::Instance()->
331 GetMpSegmentation(digit->DetElemId(),AliMp::GetCathodType(digit->Cathode()));
332 AliMpPad pad = seg->PadByIndices(digit->PadX(), digit->PadY());
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();
339 if ( digit->ManuId() & AliMpConstants::ManuMask(AliMp::kNonBendingPlane)) {
344 padInfo.SetPadId(digit->GetUniqueID());
345 padInfo.SetPadPlaneType(planeType);
346 padInfo.SetPadXY(pad.GetPositionX(), pad.GetPositionY());
347 padInfo.SetPadDimXY(pad.GetDimensionX(), pad.GetDimensionY());
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));
356 clusterInfo->AddPad(padInfo);
359 // fill cluster info tree
360 clusterInfoTree->Fill();
373 printf("Writing Tree\n");
375 clusterInfoFile->cd();
376 clusterInfoTree->Write();
377 printf("Deleting Tree\n");
378 delete clusterInfoTree;
379 printf("Closing File\n");
380 clusterInfoFile->Close();
383 printf("Deleting calibData\n");
385 printf("Deleting clusterInfo\n");
388 MUONLoader->UnloadDigits();
389 MUONLoader->UnloadRecPoints();
397 cout<<endl<<"time to fill cluster/track info: R:"<<timer.RealTime()<<" C:"<<timer.CpuTime()<<endl<<endl;
400 //-----------------------------------------------------------------------
401 TTree* GetESDTree(TFile *esdFile)
403 /// Check that the file is properly open
404 /// Return pointer to the ESD Tree
406 if (!esdFile || !esdFile->IsOpen()) {
407 Error("GetESDTree", "opening ESD file failed");
411 TTree* tree = (TTree*) esdFile->Get("esdTree");
413 Error("GetESDTree", "no ESD tree found");
422 //-----------------------------------------------------------------------
423 UInt_t buildClusterMap(AliMUONTrack &track)
425 /// Build the map of clusters in tracking chambers
427 UInt_t muonClusterMap = 0;
429 AliMUONTrackParam* trackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->First());
432 muonClusterMap |= BIT(trackParam->GetClusterPtr()->GetChamberId());
434 trackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->After(trackParam));
437 return muonClusterMap;