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