PMD info on da/amore updated
[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 <TObjArray.h>
28 #include <TTree.h>
29 #include <TString.h>
30 #include <Riostream.h>
31 #include <TRandom.h>
32 #include <TROOT.h>
33 #include <TMath.h>
34
35 // STEER includes
36 #include "AliESDEvent.h"
37 #include "AliRecoParam.h"
38 #include "AliCDBManager.h"
39 #include "AliRunLoader.h"
40 #include "AliLoader.h"
41 #include "AliHeader.h"
42
43 // MUON includes
44 #include "AliMpConstants.h"
45 #include "AliMpSegmentation.h"
46 #include "AliMpVSegmentation.h"
47 #include "AliMpPad.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"
61 #endif
62
63 const Int_t printLevel = 1;
64
65 TTree* GetESDTree(TFile *esdFile);
66 UInt_t buildClusterMap(AliMUONTrack &track);
67
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")
71 {
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   }
90   
91   AliMUONClusterInfo* clusterInfo = new AliMUONClusterInfo();
92   AliMUONPadInfo padInfo;
93   AliMUONCalibrationData* calibData = 0x0;
94   AliMUONESDInterface esdInterface;
95   AliMUONVClusterStore* clusterStore = 0x0;
96   AliMUONVDigitStore* digitStore = 0x0;
97   
98   // open the ESD file and tree and connect the ESD event
99   TFile* esdFile = 0x0;
100   TTree* esdTree = 0x0;
101   AliESDEvent* esd = 0x0;
102   Int_t runNumber = -1;
103   if (useESD) {
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");
110       return;
111     }
112     runNumber = esd->GetRunNumber();
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");
123     rl->LoadHeader();
124     if (runNumber < 0) runNumber = rl->GetHeader()->GetRun();
125   }
126   
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   
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   
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   }
160   for (Int_t iEvent = 0; iEvent < nevents; iEvent++) {
161     
162     //----------------------------------------------//
163     // -------------- process event --------------- //
164     //----------------------------------------------//
165     // get the ESD of current event
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();
176     }
177     
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       if (treeD) {
194         digitStore = AliMUONVDigitStore::Create(*treeD);
195         if ( digitStore != 0x0 ) {
196           digitStore->Clear();
197           digitStore->Connect(*treeD);
198           treeD->GetEvent(0);
199         }
200       }
201     }
202     
203     //----------------------------------------------//
204     // ------------- fill cluster info ------------ //
205     //----------------------------------------------//
206     // --- loop over the refitted tracks ---
207     if (useESD) {
208       
209       TIter nextTrack(esdInterface.CreateTrackIterator());
210       AliMUONTrack* track;
211       while ((track = static_cast<AliMUONTrack*>(nextTrack()))) {
212         
213         UInt_t muonClusterMap = buildClusterMap(*track);
214         
215         // loop over clusters
216         AliMUONTrackParam* trackParam = static_cast<AliMUONTrackParam*>(track->GetTrackParamAtCluster()->First());
217         while (trackParam) {
218           clusterInfo->Clear("C");
219           
220           // fill cluster info
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());
230           
231           // fill track info
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);
242           
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;
247             
248             // pad location
249             const AliMpVSegmentation* seg = AliMpSegmentation::Instance()->
250             GetMpSegmentation(digit->DetElemId(),AliMp::GetCathodType(digit->Cathode()));
251             AliMpPad pad = seg->PadByIndices(digit->PadX(), digit->PadY());
252             
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();
257             Int_t planeType = 0;
258             if ( digit->ManuId() & AliMpConstants::ManuMask(AliMp::kNonBendingPlane)) {
259               planeType = 1;
260             }
261             
262             // fill pad info
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));
274             
275             clusterInfo->AddPad(padInfo);
276           }
277           
278           // remove clusters attached to a track
279           if (useRecPoints) {
280             AliMUONVCluster* cl = clusterStore->FindObject(cluster->GetUniqueID());
281             if (cl) clusterStore->Remove(*cl);
282           }
283           
284           // fill cluster info tree
285           clusterInfoTree->Fill();
286           
287           trackParam = static_cast<AliMUONTrackParam*>(track->GetTrackParamAtCluster()->After(trackParam));
288         }
289         
290       }
291       
292     }
293     
294     // --- loop over clusters not attached to a track ---
295     if (useRecPoints) {
296       
297       TIter nextCluster(clusterStore->CreateIterator());
298       AliMUONVCluster *cluster;
299       while ( ( cluster = static_cast<AliMUONVCluster*>(nextCluster()) ) ) {
300         
301         clusterInfo->Clear("C");
302         
303         // fill cluster info
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());
312         
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);
323         
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;
328           
329           // pad location
330           const AliMpVSegmentation* seg = AliMpSegmentation::Instance()->
331           GetMpSegmentation(digit->DetElemId(),AliMp::GetCathodType(digit->Cathode()));
332           AliMpPad pad = seg->PadByIndices(digit->PadX(), digit->PadY());
333           
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();
338           Int_t planeType = 0;
339           if ( digit->ManuId() & AliMpConstants::ManuMask(AliMp::kNonBendingPlane)) {
340             planeType = 1;
341           }
342           
343           // fill pad info
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));
355           
356           clusterInfo->AddPad(padInfo);
357         }
358         
359         // fill cluster info tree
360         clusterInfoTree->Fill();
361         
362       }
363       
364       delete digitStore;
365       delete clusterStore;
366       
367     }
368     
369   }
370   
371   // ...timer stop
372   timer.Stop();
373   printf("Writing Tree\n");
374   // write output tree
375   clusterInfoFile->cd();
376   clusterInfoTree->Write();
377   printf("Deleting Tree\n");
378   delete clusterInfoTree;
379   printf("Closing File\n");
380   clusterInfoFile->Close();
381   
382   // free memory
383   printf("Deleting calibData\n");
384   delete calibData;
385   printf("Deleting clusterInfo\n");
386   delete clusterInfo;
387   if (useRecPoints) {
388     MUONLoader->UnloadDigits();
389     MUONLoader->UnloadRecPoints();
390     rl->UnloadHeader();
391     delete rl;
392   }
393   if (useESD) {
394     esdFile->Close();
395     delete esd;
396   }
397   cout<<endl<<"time to fill cluster/track info: R:"<<timer.RealTime()<<" C:"<<timer.CpuTime()<<endl<<endl;
398 }
399
400 //-----------------------------------------------------------------------
401 TTree* GetESDTree(TFile *esdFile)
402 {
403   /// Check that the file is properly open
404   /// Return pointer to the ESD Tree
405   
406   if (!esdFile || !esdFile->IsOpen()) {
407     Error("GetESDTree", "opening ESD file failed");
408     exit(-1);
409   }
410   
411   TTree* tree = (TTree*) esdFile->Get("esdTree");
412   if (!tree) {
413     Error("GetESDTree", "no ESD tree found");
414     exit(-1);
415   }
416   
417   return tree;
418   
419 }
420
421       
422 //-----------------------------------------------------------------------
423 UInt_t buildClusterMap(AliMUONTrack &track)
424 {
425   /// Build the map of clusters in tracking chambers
426   
427   UInt_t muonClusterMap = 0;
428   
429   AliMUONTrackParam* trackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->First());
430   while (trackParam) {
431     
432     muonClusterMap |= BIT(trackParam->GetClusterPtr()->GetChamberId());
433     
434     trackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->After(trackParam));
435   }
436   
437   return muonClusterMap;
438   
439 }
440