Update of ACORDE-QA-Raw data histograms (now they go from -0.5 to 59.5)
[u/mrichter/AliRoot.git] / MUON / MUONClusterInfo.C
CommitLineData
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
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
63const Int_t printLevel = 1;
64
fcefbac4 65TTree* GetESDTree(TFile *esdFile);
41c52850 66UInt_t buildClusterMap(AliMUONTrack &track);
fcefbac4 67
68//-----------------------------------------------------------------------
a99c3449 69void 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 399TTree* 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//-----------------------------------------------------------------------
421UInt_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