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 <TRandom.h>
  32 #include <TROOT.h>
  33 #include <TMath.h>
  34
  35 // STEER includes
  36 #include "AliMagF.h"
  37 #include "AliTracker.h"
  38 #include "AliESDEvent.h"
  39 #include "AliRecoParam.h"
  40 #include "AliCDBManager.h"
  41 #include "AliRunLoader.h"
  42 #include "AliLoader.h"
  43
  44 // MUON includes
  45 #include "AliMpConstants.h"
  46 #include "AliMpCDB.h"
  47 #include "AliMpSegmentation.h"
  48 #include "AliMpVSegmentation.h"
  49 #include "AliMpPad.h"
  50 #include "AliMUONCalibrationData.h"
  51 #include "AliMUONVCalibParam.h"
  52 #include "AliMUONPadInfo.h"
  53 #include "AliMUONClusterInfo.h"
  54 #include "AliMUONRecoParam.h"
  55 #include "AliMUONESDInterface.h"
  56 #include "AliMUONVDigit.h"
  57 #include "AliMUONVDigitStore.h"
  58 #include "AliMUONVCluster.h"
  59 #include "AliMUONVClusterStore.h"
  60 #include "AliMUONTrack.h"
  61 #include "AliMUONTrackParam.h"
  62 #endif
  63
  64 const Int_t printLevel = 1;
  65
  66 void Prepare();
  67 TTree* GetESDTree(TFile *esdFile);
  68 UInt_t buildClusterMap(AliMUONTrack &track);
  69
  70 //-----------------------------------------------------------------------
  71 void MUONClusterInfo(Int_t nevents = -1, const char* esdFileName = "AliESDs.root",
  72                      const char* inFileName = "galice.root", const char* outFileName = "clusterInfo.root")
  73 {
  74   /// 1) if (esdFileName != "")
  75   /// loop over ESD event and fill AliMUONClusterInfo object with cluster + corresponding track parameters;
  76   /// 2) if (inFileName != "")
  77   /// loop over RecPoints and fill AliMUONClusterInfo object with cluster not attached to a track;
  78   /// 3) write results in a new root file.
  79   ///
  80   /// ******************************************* WARNING ******************************************* ///
  81   /// track parameters at each cluster are recomputed by the interface using Kalman filter + Smoother ///
  82   /// (It can be changed by resetting the tracker in the interface with a new recoParam object)       ///
  83   /// and the magnetic field set in the function prepare()                                            ///
  84   /// ******************************************* WARNING ******************************************* ///
  85
  86   Bool_t useESD = (strcmp(esdFileName,""));
  87   Bool_t useRecPoints = (strcmp(inFileName,""));
  88   if (!useESD && !useRecPoints) {
  89     Error("MUONClusterInfo","you must provide ESD and/or galice root file(s)");
  90     return;
  91   }
  92
  93   AliMUONClusterInfo* clusterInfo = new AliMUONClusterInfo();
  94   AliMUONPadInfo padInfo;
  95   AliMUONCalibrationData* calibData = 0x0;
  96   AliMUONESDInterface esdInterface;
  97   AliMUONVClusterStore* clusterStore = 0x0;
  98   AliMUONVDigitStore* digitStore = 0x0;
  99
 100   // prepare the refitting during ESD->MUON conversion
 101   Prepare();
 102
 103   // open the ESD file and tree and connect the ESD event
 104   TFile* esdFile = 0x0;
 105   TTree* esdTree = 0x0;
 106   AliESDEvent* esd = 0x0;
 107   if (useESD) {
 108     esdFile = TFile::Open(esdFileName);
 109     esdTree = GetESDTree(esdFile);
 110     esd = new AliESDEvent();
 111     esd->ReadFromTree(esdTree);
 112   }
 113
 114   // get the cluster from RecPoints
 115   AliRunLoader * rl = 0x0;
 116   AliLoader* MUONLoader = 0x0;
 117   if (useRecPoints) {
 118     rl = AliRunLoader::Open(inFileName,"MUONLoader");
 119     MUONLoader = rl->GetDetectorLoader("MUON");
 120     MUONLoader->LoadRecPoints("READ");
 121     MUONLoader->LoadDigits("READ");
 122   }
 123
 124   // prepare the output tree
 125   gROOT->cd();
 126   TFile* clusterInfoFile = TFile::Open(outFileName, "RECREATE");
 127   clusterInfoFile->SetCompressionLevel(1);
 128
 129   TTree* clusterInfoTree = new TTree("clusterInfoTree","clusterInfoTree");
 130   clusterInfoTree->Branch("clusterInfo", &clusterInfo, 32000, 99);
 131
 132   // timer start...
 133   TStopwatch timer;
 134
 135   // Loop over events
 136   if (useESD) {
 137     if (nevents > 0) nevents = TMath::Min(nevents,(Int_t)esdTree->GetEntries());
 138     else nevents = (Int_t)esdTree->GetEntries();
 139   } else {
 140     if (nevents > 0) nevents = TMath::Min(nevents,(Int_t)rl->GetNumberOfEvents());
 141     else nevents = (Int_t)rl->GetNumberOfEvents();
 142   }
 143   for (Int_t iEvent = 0; iEvent < nevents; iEvent++) {
 144
 145     //----------------------------------------------//
 146     // -------------- process event --------------- //
 147     //----------------------------------------------//
 148     // get the ESD of current event
 149     if (useESD) {
 150       esdTree->GetEvent(iEvent);
 151       if (!esd) {
 152         Error("MUONClusterInfo", "no ESD object found for event %d", iEvent);
 153         return;
 154       }
 155       // load the current esd event
 156       esdInterface.LoadEvent(*esd);
 157       // get digit store
 158       if (!useRecPoints) digitStore = esdInterface.GetDigits();
 159     }
 160
 161     if (useRecPoints) {
 162       if (!(rl->GetEvent(iEvent) == 0)) {
 163         Error("MUONClusterInfo", "unable to load event %d", iEvent);
 164         return;
 165       }
 166       // get the clusters of current event
 167       TTree* treeR = MUONLoader->TreeR();
 168       clusterStore = AliMUONVClusterStore::Create(*treeR);
 169       if ( clusterStore != 0x0 ) {
 170         clusterStore->Clear();
 171         clusterStore->Connect(*treeR);
 172         treeR->GetEvent(0);
 173       }
 174       // get the digits of current event
 175       TTree* treeD = MUONLoader->TreeD();
 176       digitStore = AliMUONVDigitStore::Create(*treeD);
 177       if ( digitStore != 0x0 ) {
 178         digitStore->Clear();
 179         digitStore->Connect(*treeD);
 180         treeD->GetEvent(0);
 181       }
 182     }
 183
 184     // prepare access to calibration data
 185     if (useESD && !calibData) calibData = new AliMUONCalibrationData(esd->GetESDRun()->GetRunNumber());
 186     else if (!calibData) calibData = new AliMUONCalibrationData(rl->GetRunNumber());
 187
 188     //----------------------------------------------//
 189     // ------------- fill cluster info ------------ //
 190     //----------------------------------------------//
 191     // --- loop over the refitted tracks ---
 192     if (useESD) {
 193
 194       TIter nextTrack(esdInterface.CreateTrackIterator());
 195       AliMUONTrack* track;
 196       while ((track = static_cast<AliMUONTrack*>(nextTrack()))) {
 197
 198         UInt_t muonClusterMap = buildClusterMap(*track);
 199
 200         // loop over clusters
 201         AliMUONTrackParam* trackParam = static_cast<AliMUONTrackParam*>(track->GetTrackParamAtCluster()->First());
 202         while (trackParam) {
 203           clusterInfo->Clear("C");
 204
 205           // fill cluster info
 206           AliMUONVCluster* cluster = trackParam->GetClusterPtr();
 207           clusterInfo->SetRunId(esd->GetRunNumber());
 208           clusterInfo->SetEventId(iEvent);
 209           clusterInfo->SetZ(cluster->GetZ());
 210           clusterInfo->SetClusterId(cluster->GetUniqueID());
 211           clusterInfo->SetClusterXY(cluster->GetX(), cluster->GetY());
 212           clusterInfo->SetClusterXYErr(cluster->GetErrX(), cluster->GetErrY());
 213           clusterInfo->SetClusterChi2(cluster->GetChi2());
 214           clusterInfo->SetClusterCharge(cluster->GetCharge());
 215
 216           // fill track info
 217           clusterInfo->SetTrackId(track->GetUniqueID());
 218           clusterInfo->SetTrackXY(trackParam->GetNonBendingCoor(), trackParam->GetBendingCoor());
 219           clusterInfo->SetTrackThetaXY(TMath::ATan(trackParam->GetNonBendingSlope()), TMath::ATan(trackParam->GetBendingSlope()));
 220           clusterInfo->SetTrackP(trackParam->P());
 221           const TMatrixD paramCov = trackParam->GetCovariances();
 222           clusterInfo->SetTrackXYErr(TMath::Sqrt(paramCov(0,0)), TMath::Sqrt(paramCov(2,2)));
 223           clusterInfo->SetTrackChi2(track->GetNormalizedChi2());
 224           clusterInfo->SetTrackCharge((Short_t)trackParam->GetCharge());
 225           clusterInfo->SetTrackNHits(track->GetNClusters());
 226           clusterInfo->SetTrackChamberHitMap(muonClusterMap);
 227
 228           // fill pad info if available
 229           for (Int_t i=0; i<cluster->GetNDigits(); i++) {
 230             AliMUONVDigit* digit = digitStore->FindObject(cluster->GetDigitId(i));
 231             if (!digit) continue;
 232
 233             // pad location
 234             const AliMpVSegmentation* seg = AliMpSegmentation::Instance()->
 235             GetMpSegmentation(digit->DetElemId(),AliMp::GetCathodType(digit->Cathode()));
 236             AliMpPad pad = seg->PadByIndices(digit->PadX(), digit->PadY());
 237
 238             // calibration parameters
 239             AliMUONVCalibParam* ped = calibData->Pedestals(digit->DetElemId(), digit->ManuId());
 240             AliMUONVCalibParam* gain = calibData->Gains(digit->DetElemId(), digit->ManuId());
 241             Int_t manuChannel = digit->ManuChannel();
 242             Int_t planeType = 0;
 243             if ( digit->ManuId() & AliMpConstants::ManuMask(AliMp::kNonBendingPlane)) {
 244               planeType = 1;
 245             }
 246
 247             // fill pad info
 248             padInfo.SetPadId(digit->GetUniqueID());
 249             padInfo.SetPadPlaneType(planeType);
 250             padInfo.SetPadXY(pad.GetPositionX(), pad.GetPositionY());
 251             padInfo.SetPadDimXY(pad.GetDimensionX(), pad.GetDimensionY());
 252             padInfo.SetPadCharge((Double_t)digit->Charge());
 253             padInfo.SetPadADC(digit->ADC());
 254             padInfo.SetSaturated(digit->IsSaturated());
 255             padInfo.SetCalibrated(digit->IsCalibrated());
 256             padInfo.SetPedestal(ped->ValueAsFloatFast(manuChannel,0), ped->ValueAsFloatFast(manuChannel,1));
 257             padInfo.SetGain(gain->ValueAsFloatFast(manuChannel,0), gain->ValueAsFloatFast(manuChannel,1),
 258                             gain->ValueAsFloatFast(manuChannel,2), gain->ValueAsFloatFast(manuChannel,3));
 259
 260             clusterInfo->AddPad(padInfo);
 261           }
 262
 263           // remove clusters attached to a track
 264           if (useRecPoints) {
 265             AliMUONVCluster* cl = clusterStore->FindObject(cluster->GetUniqueID());
 266             if (cl) clusterStore->Remove(*cl);
 267           }
 268
 269           // fill cluster info tree
 270           clusterInfoTree->Fill();
 271
 272           trackParam = static_cast<AliMUONTrackParam*>(track->GetTrackParamAtCluster()->After(trackParam));
 273         }
 274
 275       }
 276
 277     }
 278
 279     // --- loop over clusters not attached to a track ---
 280     if (useRecPoints) {
 281
 282       TIter nextCluster(clusterStore->CreateIterator());
 283       AliMUONVCluster *cluster;
 284       while ( ( cluster = static_cast<AliMUONVCluster*>(nextCluster()) ) ) {
 285
 286         clusterInfo->Clear("C");
 287
 288         // fill cluster info
 289         clusterInfo->SetRunId(rl->GetRunNumber());
 290         clusterInfo->SetEventId(iEvent);
 291         clusterInfo->SetZ(cluster->GetZ());
 292         clusterInfo->SetClusterId(cluster->GetUniqueID());
 293         clusterInfo->SetClusterXY(cluster->GetX(), cluster->GetY());
 294         clusterInfo->SetClusterXYErr(cluster->GetErrX(), cluster->GetErrY());
 295         clusterInfo->SetClusterChi2(cluster->GetChi2());
 296         clusterInfo->SetClusterCharge(cluster->GetCharge());
 297
 298         // fill dummy track info
 299         clusterInfo->SetTrackId(0);
 300         clusterInfo->SetTrackXY(0.,0.);
 301         clusterInfo->SetTrackThetaXY(0.,0.);
 302         clusterInfo->SetTrackP(0.);
 303         clusterInfo->SetTrackXYErr(0.,0.);
 304         clusterInfo->SetTrackChi2(0.);
 305         clusterInfo->SetTrackCharge(0);
 306         clusterInfo->SetTrackNHits(0);
 307         clusterInfo->SetTrackChamberHitMap(0);
 308
 309         // fill pad info if available
 310         for (Int_t i=0; i<cluster->GetNDigits(); i++) {
 311           AliMUONVDigit* digit = digitStore->FindObject(cluster->GetDigitId(i));
 312           if (!digit) continue;
 313
 314           // pad location
 315           const AliMpVSegmentation* seg = AliMpSegmentation::Instance()->
 316           GetMpSegmentation(digit->DetElemId(),AliMp::GetCathodType(digit->Cathode()));
 317           AliMpPad pad = seg->PadByIndices(digit->PadX(), digit->PadY());
 318
 319           // calibration parameters
 320           AliMUONVCalibParam* ped = calibData->Pedestals(digit->DetElemId(), digit->ManuId());
 321           AliMUONVCalibParam* gain = calibData->Gains(digit->DetElemId(), digit->ManuId());
 322           Int_t manuChannel = digit->ManuChannel();
 323           Int_t planeType = 0;
 324           if ( digit->ManuId() & AliMpConstants::ManuMask(AliMp::kNonBendingPlane)) {
 325             planeType = 1;
 326           }
 327
 328           // fill pad info
 329           padInfo.SetPadId(digit->GetUniqueID());
 330           padInfo.SetPadPlaneType(planeType);
 331           padInfo.SetPadXY(pad.GetPositionX(), pad.GetPositionY());
 332           padInfo.SetPadDimXY(pad.GetDimensionX(), pad.GetDimensionY());
 333           padInfo.SetPadCharge((Double_t)digit->Charge());
 334           padInfo.SetPadADC(digit->ADC());
 335           padInfo.SetSaturated(digit->IsSaturated());
 336           padInfo.SetCalibrated(digit->IsCalibrated());
 337           padInfo.SetPedestal(ped->ValueAsFloatFast(manuChannel,0), ped->ValueAsFloatFast(manuChannel,1));
 338           padInfo.SetGain(gain->ValueAsFloatFast(manuChannel,0), gain->ValueAsFloatFast(manuChannel,1),
 339                           gain->ValueAsFloatFast(manuChannel,2), gain->ValueAsFloatFast(manuChannel,3));
 340
 341           clusterInfo->AddPad(padInfo);
 342         }
 343
 344         // fill cluster info tree
 345         clusterInfoTree->Fill();
 346
 347       }
 348
 349       delete digitStore;
 350       delete clusterStore;
 351
 352     }
 353
 354   }
 355
 356   // ...timer stop
 357   timer.Stop();
 358   printf("Writing Tree\n");
 359   // write output tree
 360   clusterInfoFile->cd();
 361   clusterInfoTree->Write();
 362   printf("Deleting Tree\n");
 363   delete clusterInfoTree;
 364   printf("Closing File\n");
 365   clusterInfoFile->Close();
 366
 367   // free memory
 368   printf("Deleting calibData\n");
 369   delete calibData;
 370   printf("Deleting clusterInfo\n");
 371   delete clusterInfo;
 372   if (useRecPoints) {
 373     MUONLoader->UnloadDigits();
 374     MUONLoader->UnloadRecPoints();
 375     delete rl;
 376   }
 377   if (useESD) {
 378     esdFile->Close();
 379     delete esd;
 380   }
 381   cout<<endl<<"time to fill cluster/track info: R:"<<timer.RealTime()<<" C:"<<timer.CpuTime()<<endl<<endl;
 382 }
 383
 384 //-----------------------------------------------------------------------
 385 void Prepare()
 386 {
 387   /// Set the magnetic field, the mapping and the reconstruction parameters
 388
 389   gRandom->SetSeed(0);
 390
 391   // set mag field
 392   if (!TGeoGlobalMagField::Instance()->GetField()) {
 393     printf("Loading field map...\n");
 394     AliMagF* field = new AliMagF("Maps","Maps",2,1.,1., 10.,AliMagF::k5kG);
 395     TGeoGlobalMagField::Instance()->SetField(field);
 396   }
 397
 398   // Load mapping
 399   AliCDBManager* man = AliCDBManager::Instance();
 400   man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
 401   man->SetRun(0);
 402   if ( ! AliMpCDB::LoadDDLStore() ) {
 403     Error("Prepare","Could not access mapping from OCDB !");
 404     exit(-1);
 405   }
 406
 407   // Reset the reconstruction parameters for track refitting if needed
 408   // (by default will use Kalman filter + Smoother)
 409   //  AliMUONRecoParam *muonRecoParam = AliMUONRecoParam::GetLowFluxParam();
 410   //  AliMUONESDInterface::ResetTracker(muonRecoParam);
 411
 412 }
 413
 414 //-----------------------------------------------------------------------
 415 TTree* GetESDTree(TFile *esdFile)
 416 {
 417   /// Check that the file is properly open
 418   /// Return pointer to the ESD Tree
 419
 420   if (!esdFile || !esdFile->IsOpen()) {
 421     Error("GetESDTree", "opening ESD file failed");
 422     exit(-1);
 423   }
 424
 425   TTree* tree = (TTree*) esdFile->Get("esdTree");
 426   if (!tree) {
 427     Error("GetESDTree", "no ESD tree found");
 428     exit(-1);
 429   }
 430
 431   return tree;
 432
 433 }
 434
 435
 436 //-----------------------------------------------------------------------
 437 UInt_t buildClusterMap(AliMUONTrack &track)
 438 {
 439   /// Build the map of clusters in tracking chambers
 440
 441   UInt_t muonClusterMap = 0;
 442
 443   AliMUONTrackParam* trackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->First());
 444   while (trackParam) {
 445
 446     muonClusterMap |= BIT(trackParam->GetClusterPtr()->GetChamberId());
 447
 448     trackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->After(trackParam));
 449   }
 450
 451   return muonClusterMap;
 452
 453 }
 454