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 <TClonesArray.h>
30 #include <Riostream.h>
31 #include <TGeoManager.h>
36 #include "AliMagFMaps.h"
37 #include "AliTracker.h"
38 #include "AliESDEvent.h"
39 #include "AliESDMuonTrack.h"
40 #include "AliRecoParam.h"
41 #include "AliCDBManager.h"
42 #include "AliGeomManager.h"
46 #include "AliMpSegmentation.h"
47 #include "AliMpVSegmentation.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 "AliMUONTrack.h"
59 #include "AliMUONTrackParam.h"
62 const Int_t printLevel = 1;
65 TTree* GetESDTree(TFile *esdFile);
67 //-----------------------------------------------------------------------
68 void MUONClusterInfo(Int_t nevents = -1, const char* esdFileName = "AliESDs.root", const char* outFileName = "clusterInfo.root")
70 /// load ESD event in the ESDInterface to recover MUON objects;
71 /// fill AliMUONESDClusterInfo object with ESD cluster + corresponding track parameters;
72 /// write results in a new root file.
74 AliMUONClusterInfo* clusterInfo = new AliMUONClusterInfo();
75 AliMUONPadInfo padInfo;
76 AliMUONCalibrationData* calibData = 0x0;
77 AliMUONESDInterface esdInterface;
79 // prepare the refitting during ESD->MUON conversion
82 // open the ESD file and tree and connect the ESD event
83 TFile* esdFile = TFile::Open(esdFileName);
84 TTree* esdTree = GetESDTree(esdFile);
85 AliESDEvent* esd = new AliESDEvent();
86 esd->ReadFromTree(esdTree);
88 // prepare the output tree
90 TFile* clusterInfoFile = TFile::Open(outFileName, "RECREATE");
91 clusterInfoFile->SetCompressionLevel(1);
93 TTree* clusterInfoTree = new TTree("clusterInfoTree","clusterInfoTree");
94 clusterInfoTree->Branch("clusterInfo", &clusterInfo, 32000, 99);
99 // Loop over ESD events
100 if (nevents > 0) nevents = TMath::Min(nevents,(Int_t)esdTree->GetEntries());
101 else nevents = (Int_t)esdTree->GetEntries();
102 for (Int_t iEvent = 0; iEvent < nevents; iEvent++) {
104 //----------------------------------------------//
105 // -------------- process event --------------- //
106 //----------------------------------------------//
107 // get the ESD of current event
108 esdTree->GetEvent(iEvent);
110 Error("CheckESD", "no ESD object found for event %d", iEvent);
113 if (esd->GetNumberOfMuonTracks() < 1) continue;
115 // prepare access to calibration data
116 if (!calibData) calibData = new AliMUONCalibrationData(esd->GetESDRun()->GetRunNumber());
118 // load the current event
119 esdInterface.LoadEvent(*esd);
122 AliMUONVDigitStore* digitStore = esdInterface.GetDigits();
124 //----------------------------------------------//
125 // ------------- fill cluster info ------------ //
126 //----------------------------------------------//
127 // loop over the refitted tracks
128 TIter nextTrack(esdInterface.CreateTrackIterator());
130 while ((track = static_cast<AliMUONTrack*>(nextTrack()))) {
132 // loop over clusters
133 AliMUONTrackParam* trackParam = static_cast<AliMUONTrackParam*>(track->GetTrackParamAtCluster()->First());
135 clusterInfo->Clear("C");
138 AliMUONVCluster* cluster = trackParam->GetClusterPtr();
139 clusterInfo->SetEventId(iEvent);
140 clusterInfo->SetZ(cluster->GetZ());
141 clusterInfo->SetClusterId(cluster->GetUniqueID());
142 clusterInfo->SetClusterXY(cluster->GetX(), cluster->GetY());
143 clusterInfo->SetClusterXYErr(cluster->GetErrX(), cluster->GetErrY());
144 clusterInfo->SetClusterChi2(cluster->GetChi2());
145 clusterInfo->SetClusterCharge(cluster->GetCharge());
148 clusterInfo->SetTrackId(track->GetUniqueID());
149 clusterInfo->SetTrackXY(trackParam->GetNonBendingCoor(), trackParam->GetBendingCoor());
150 clusterInfo->SetTrackThetaXY(TMath::ATan(trackParam->GetBendingSlope()), TMath::ATan(trackParam->GetNonBendingSlope()));
151 clusterInfo->SetTrackP(trackParam->P());
152 const TMatrixD paramCov = trackParam->GetCovariances();
153 clusterInfo->SetTrackXYErr(TMath::Sqrt(paramCov(0,0)), TMath::Sqrt(paramCov(2,2)));
154 clusterInfo->SetTrackChi2(track->GetNormalizedChi2());
155 clusterInfo->SetTrackCharge((Short_t)trackParam->GetCharge());
157 // fill pad info if available
158 for (Int_t i=0; i<cluster->GetNDigits(); i++) {
159 AliMUONVDigit* digit = digitStore->FindObject(cluster->GetDigitId(i));
160 if (!digit) continue;
163 const AliMpVSegmentation* seg = AliMpSegmentation::Instance()->
164 GetMpSegmentation(digit->DetElemId(),AliMp::GetCathodType(digit->Cathode()));
165 AliMpPad pad = seg->PadByIndices(AliMpIntPair(digit->PadX(), digit->PadY()));
167 // calibration parameters
168 AliMUONVCalibParam* ped = calibData->Pedestals(digit->DetElemId(), digit->ManuId());
169 AliMUONVCalibParam* gain = calibData->Gains(digit->DetElemId(), digit->ManuId());
170 Int_t manuChannel = digit->ManuChannel();
173 padInfo.SetPadId(digit->GetUniqueID());
174 padInfo.SetPadXY(pad.Position().X(), pad.Position().Y());
175 padInfo.SetPadDimXY(pad.Dimensions().X(), pad.Dimensions().Y());
176 padInfo.SetPadCharge((Double_t)digit->Charge());
177 padInfo.SetPadADC(digit->ADC());
178 padInfo.SetSaturated(digit->IsSaturated());
179 padInfo.SetCalibrated(digit->IsCalibrated());
180 padInfo.SetPedestal(ped->ValueAsFloatFast(manuChannel,0), ped->ValueAsFloatFast(manuChannel,1));
181 padInfo.SetGain(gain->ValueAsFloatFast(manuChannel,0), gain->ValueAsFloatFast(manuChannel,1),
182 gain->ValueAsFloatFast(manuChannel,2), gain->ValueAsFloatFast(manuChannel,3));
184 clusterInfo->AddPad(padInfo);
187 // fill cluster info tree
188 clusterInfoTree->Fill();
190 trackParam = static_cast<AliMUONTrackParam*>(track->GetTrackParamAtCluster()->After(trackParam));
199 printf("Writing Tree\n");
201 clusterInfoFile->cd();
202 clusterInfoTree->Write();
203 printf("Deleting Tree\n");
204 delete clusterInfoTree;
205 printf("Closing File\n");
206 clusterInfoFile->Close();
209 printf("Deleting calibData\n");
211 printf("Deleting clusterInfo\n");
213 printf("Closing esdFile\n");
215 printf("Deleting esd\n");
218 cout<<endl<<"time to fill cluster/track info: R:"<<timer.RealTime()<<" C:"<<timer.CpuTime()<<endl<<endl;
221 //-----------------------------------------------------------------------
224 /// Set the geometry, the magnetic field, the mapping and the reconstruction parameters
228 // Import TGeo geometry (needed for track extrapolation)
230 AliGeomManager::LoadGeometry("geometry.root");
232 Error("MUONRefit", "getting geometry from file %s failed", "generated/galice.root");
238 printf("Loading field map...\n");
239 AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 1, 1., 10., AliMagFMaps::k5kG);
240 AliTracker::SetFieldMap(field, kFALSE);
243 AliCDBManager* man = AliCDBManager::Instance();
244 man->SetDefaultStorage("local://$ALICE_ROOT");
246 if ( ! AliMpCDB::LoadDDLStore() ) {
247 Error("MUONRefit","Could not access mapping from OCDB !");
251 // eventually set reconstruction parameters for refit (otherwise read from OCDB)
252 /* AliMUONRecoParam *muonRecoParam = AliMUONRecoParam::GetLowFluxParam();
253 muonRecoParam->Print("FULL");
254 AliRecoParam::Instance()->RegisterRecoParam(muonRecoParam);
258 //-----------------------------------------------------------------------
259 TTree* GetESDTree(TFile *esdFile)
261 /// Check that the file is properly open
262 /// Return pointer to the ESD Tree
264 if (!esdFile || !esdFile->IsOpen()) {
265 Error("GetESDTree", "opening ESD file failed");
269 TTree* tree = (TTree*) esdFile->Get("esdTree");
271 Error("GetESDTree", "no ESD tree found");