02f8e796407d66f837f2414c5cc4bd7c518693c8
[u/mrichter/AliRoot.git] / MUON / MUONClusterInfo.C
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 <TGeoManager.h>
32 #include <TRandom.h>
33 #include <TROOT.h>
34
35 // STEER includes
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"
43
44 // MUON includes
45 #include "AliMpCDB.h"
46 #include "AliMpSegmentation.h"
47 #include "AliMpVSegmentation.h"
48 #include "AliMpPad.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"
60 #endif
61
62 const Int_t printLevel = 1;
63
64 void Prepare();
65 TTree* GetESDTree(TFile *esdFile);
66
67 //-----------------------------------------------------------------------
68 void MUONClusterInfo(Int_t nevents = -1, const char* esdFileName = "AliESDs.root", const char* outFileName = "clusterInfo.root")
69 {
70   /// load ESD event in the ESDInterface to recover MUON objects;
71   /// track parameters at each cluster are recomputed by the interface using Kalman filter + Smoother
72   /// (It can be changed by resetting the tracker in the interface with a new recoParam object);
73   /// fill AliMUONESDClusterInfo object with ESD cluster + corresponding track parameters;
74   /// write results in a new root file.
75   
76   AliMUONClusterInfo* clusterInfo = new AliMUONClusterInfo();
77   AliMUONPadInfo padInfo;
78   AliMUONCalibrationData* calibData = 0x0;
79   AliMUONESDInterface esdInterface;
80   
81   // prepare the refitting during ESD->MUON conversion
82   Prepare();
83   
84   // open the ESD file and tree and connect the ESD event
85   TFile* esdFile = TFile::Open(esdFileName);
86   TTree* esdTree = GetESDTree(esdFile);
87   AliESDEvent* esd = new AliESDEvent();
88   esd->ReadFromTree(esdTree);
89   
90   // prepare the output tree
91   gROOT->cd();
92   TFile* clusterInfoFile = TFile::Open(outFileName, "RECREATE");
93   clusterInfoFile->SetCompressionLevel(1);
94   
95   TTree* clusterInfoTree = new TTree("clusterInfoTree","clusterInfoTree");
96   clusterInfoTree->Branch("clusterInfo", &clusterInfo, 32000, 99);
97   
98   // timer start...
99   TStopwatch timer;
100   
101   // Loop over ESD events
102   if (nevents > 0) nevents = TMath::Min(nevents,(Int_t)esdTree->GetEntries());
103   else nevents = (Int_t)esdTree->GetEntries();
104   for (Int_t iEvent = 0; iEvent < nevents; iEvent++) {
105     
106     //----------------------------------------------//
107     // -------------- process event --------------- //
108     //----------------------------------------------//
109     // get the ESD of current event
110     esdTree->GetEvent(iEvent);
111     if (!esd) {
112       Error("CheckESD", "no ESD object found for event %d", iEvent);
113       return;
114     }
115     if (esd->GetNumberOfMuonTracks() < 1) continue;
116     
117     // prepare access to calibration data
118     if (!calibData) calibData = new AliMUONCalibrationData(esd->GetESDRun()->GetRunNumber());
119     
120     // load the current event
121     esdInterface.LoadEvent(*esd);
122     
123     // get digit store
124     AliMUONVDigitStore* digitStore = esdInterface.GetDigits();
125     
126     //----------------------------------------------//
127     // ------------- fill cluster info ------------ //
128     //----------------------------------------------//
129     // loop over the refitted tracks
130     TIter nextTrack(esdInterface.CreateTrackIterator());
131     AliMUONTrack* track;
132     while ((track = static_cast<AliMUONTrack*>(nextTrack()))) {
133       
134       // loop over clusters
135       AliMUONTrackParam* trackParam = static_cast<AliMUONTrackParam*>(track->GetTrackParamAtCluster()->First());
136       while (trackParam) {
137         clusterInfo->Clear("C");
138         
139         // fill cluster info
140         AliMUONVCluster* cluster = trackParam->GetClusterPtr();
141         clusterInfo->SetEventId(iEvent);
142         clusterInfo->SetZ(cluster->GetZ());
143         clusterInfo->SetClusterId(cluster->GetUniqueID());
144         clusterInfo->SetClusterXY(cluster->GetX(), cluster->GetY());
145         clusterInfo->SetClusterXYErr(cluster->GetErrX(), cluster->GetErrY());
146         clusterInfo->SetClusterChi2(cluster->GetChi2());
147         clusterInfo->SetClusterCharge(cluster->GetCharge());
148         
149         // fill track info
150         clusterInfo->SetTrackId(track->GetUniqueID());
151         clusterInfo->SetTrackXY(trackParam->GetNonBendingCoor(), trackParam->GetBendingCoor());
152         clusterInfo->SetTrackThetaXY(TMath::ATan(trackParam->GetBendingSlope()), TMath::ATan(trackParam->GetNonBendingSlope()));
153         clusterInfo->SetTrackP(trackParam->P());
154         const TMatrixD paramCov = trackParam->GetCovariances();
155         clusterInfo->SetTrackXYErr(TMath::Sqrt(paramCov(0,0)), TMath::Sqrt(paramCov(2,2)));
156         clusterInfo->SetTrackChi2(track->GetNormalizedChi2());
157         clusterInfo->SetTrackCharge((Short_t)trackParam->GetCharge());
158
159         // fill pad info if available     
160         for (Int_t i=0; i<cluster->GetNDigits(); i++) {
161           AliMUONVDigit* digit = digitStore->FindObject(cluster->GetDigitId(i));
162           if (!digit) continue;
163           
164           // pad location
165           const AliMpVSegmentation* seg = AliMpSegmentation::Instance()->
166             GetMpSegmentation(digit->DetElemId(),AliMp::GetCathodType(digit->Cathode()));
167           AliMpPad pad = seg->PadByIndices(AliMpIntPair(digit->PadX(), digit->PadY()));
168           
169           // calibration parameters
170           AliMUONVCalibParam* ped = calibData->Pedestals(digit->DetElemId(), digit->ManuId());
171           AliMUONVCalibParam* gain = calibData->Gains(digit->DetElemId(), digit->ManuId());
172           Int_t manuChannel = digit->ManuChannel();
173           
174           // fill pad info
175           padInfo.SetPadId(digit->GetUniqueID());
176           padInfo.SetPadXY(pad.Position().X(), pad.Position().Y());
177           padInfo.SetPadDimXY(pad.Dimensions().X(), pad.Dimensions().Y());
178           padInfo.SetPadCharge((Double_t)digit->Charge());
179           padInfo.SetPadADC(digit->ADC());
180           padInfo.SetSaturated(digit->IsSaturated());
181           padInfo.SetCalibrated(digit->IsCalibrated());
182           padInfo.SetPedestal(ped->ValueAsFloatFast(manuChannel,0), ped->ValueAsFloatFast(manuChannel,1));
183           padInfo.SetGain(gain->ValueAsFloatFast(manuChannel,0), gain->ValueAsFloatFast(manuChannel,1),
184                           gain->ValueAsFloatFast(manuChannel,2), gain->ValueAsFloatFast(manuChannel,3));
185           
186           clusterInfo->AddPad(padInfo);
187         }
188         
189         // fill cluster info tree
190         clusterInfoTree->Fill();
191         
192         trackParam = static_cast<AliMUONTrackParam*>(track->GetTrackParamAtCluster()->After(trackParam));
193       }
194       
195     }
196     
197   }
198   
199   // ...timer stop
200   timer.Stop();
201   printf("Writing Tree\n");
202   // write output tree
203   clusterInfoFile->cd();
204   clusterInfoTree->Write();
205   printf("Deleting Tree\n");
206   delete clusterInfoTree;
207   printf("Closing File\n");
208   clusterInfoFile->Close();
209   
210   // free memory
211   printf("Deleting calibData\n");
212   delete calibData;
213   printf("Deleting clusterInfo\n");
214   delete clusterInfo;
215   printf("Closing esdFile\n");
216   esdFile->Close();
217   printf("Deleting esd\n");
218   delete esd;
219   //  delete padInfo;
220   cout<<endl<<"time to fill cluster/track info: R:"<<timer.RealTime()<<" C:"<<timer.CpuTime()<<endl<<endl;
221 }
222
223 //-----------------------------------------------------------------------
224 void Prepare()
225 {
226   /// Set the geometry, the magnetic field, the mapping and the reconstruction parameters
227   
228   gRandom->SetSeed(0);
229   
230   // Import TGeo geometry (needed for track extrapolation)
231   if (!gGeoManager) {
232     AliGeomManager::LoadGeometry("geometry.root");
233     if (!gGeoManager) {
234       Error("MUONRefit", "getting geometry from file %s failed", "generated/galice.root");
235       return;
236     }
237   }
238   
239   // set mag field
240   printf("Loading field map...\n");
241   AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 1, 1., 10., AliMagFMaps::k5kG);
242   AliTracker::SetFieldMap(field, kFALSE);
243   
244   // Load mapping
245   AliCDBManager* man = AliCDBManager::Instance();
246   man->SetDefaultStorage("local://$ALICE_ROOT");
247   man->SetRun(0);
248   if ( ! AliMpCDB::LoadDDLStore() ) {
249     Error("MUONRefit","Could not access mapping from OCDB !");
250     exit(-1);
251   }
252   
253   // Reset the reconstruction parameters for track refitting if needed
254   // (by default will use Kalman filter + Smoother)
255 //  AliMUONRecoParam *muonRecoParam = AliMUONRecoParam::GetLowFluxParam();
256 //  AliMUONESDInterface::ResetTracker(muonRecoParam);
257   
258 }
259
260 //-----------------------------------------------------------------------
261 TTree* GetESDTree(TFile *esdFile)
262 {
263   /// Check that the file is properly open
264   /// Return pointer to the ESD Tree
265   
266   if (!esdFile || !esdFile->IsOpen()) {
267     Error("GetESDTree", "opening ESD file failed");
268     exit(-1);
269   }
270   
271   TTree* tree = (TTree*) esdFile->Get("esdTree");
272   if (!tree) {
273     Error("GetESDTree", "no ESD tree found");
274     exit(-1);
275   }
276   
277   return tree;
278   
279 }
280