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