MUON/MUONRefit.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 MUONRefit.C
  20 /// \brief Macro for refitting ESD tracks from ESD pads
  21 ///
  22 /// \author Philippe Pillot, SUBATECH
  23
  24 #if !defined(__CINT__) || defined(__MAKECINT__)
  25 #include <TStopwatch.h>
  26 #include <TFile.h>
  27 #include <TTree.h>
  28 #include <TString.h>
  29 #include <Riostream.h>
  30 #include <TGeoManager.h>
  31 #include <TRandom.h>
  32 #include <TROOT.h>
  33
  34 // STEER includes
  35 #include "AliESDEvent.h"
  36 #include "AliESDMuonTrack.h"
  37 #include "AliCDBManager.h"
  38 #include "AliGeomManager.h"
  39
  40 // MUON includes
  41 #include "AliMUONCDB.h"
  42 #include "AliMUONRecoParam.h"
  43 #include "AliMUONESDInterface.h"
  44 #include "AliMUONRefitter.h"
  45 #include "AliMUONVDigit.h"
  46 #include "AliMUONTrack.h"
  47 #include "AliMUONVTrackStore.h"
  48 #include "AliMUONTrackParam.h"
  49 #endif
  50
  51 const Int_t printLevel = 1;
  52 const Bool_t reconstructFromDigits = kTRUE; // kFALSE = reconstruct from clusters
  53
  54 TTree* GetESDTree(TFile *esdFile);
  55
  56 //-----------------------------------------------------------------------
  57 void MUONRefit(Int_t nevents = -1, const char* esdFileNameIn = "AliESDs.root", const char* esdFileNameOut = "AliESDs_New.root",
  58                const char* geoFilename = "geometry.root", const char* ocdbPath = "local://$ALICE_ROOT/OCDB")
  59 {
  60   /// Example of muon refitting:
  61   /// refit ESD tracks from ESD pads (i.e. re-clusterized the attached ESD clusters);
  62   /// reset the charge of the digit using their raw charge before refitting;
  63   /// compare results with original ESD tracks;
  64   /// write results in a new ESD file
  65
  66   // open the ESD file and tree
  67   TFile* esdFile = TFile::Open(esdFileNameIn);
  68   TTree* esdTree = GetESDTree(esdFile);
  69
  70   // connect ESD event to the ESD tree
  71   AliESDEvent* esd = new AliESDEvent();
  72   esd->ReadFromTree(esdTree);
  73
  74   // create the ESD output file and tree
  75   TFile* newESDFile = TFile::Open(esdFileNameOut, "RECREATE");
  76   newESDFile->SetCompressionLevel(2);
  77
  78   // connect ESD event to the ESD tree (recreate track branch for backward compatibility)
  79   esdTree->SetBranchStatus("*MuonTracks*",0);
  80   TTree* newESDTree = esdTree->CloneTree(0);
  81   esdTree->SetBranchStatus("*MuonTracks*",1);
  82   esd->WriteToTree(newESDTree);
  83
  84   // get run number
  85   if (esdTree->GetEvent(0) <= 0) {
  86     Error("MUONRefit", "no ESD object found for event 0");
  87     return;
  88   }
  89   Int_t runNumber = esd->GetRunNumber();
  90
  91   // Import TGeo geometry
  92   if (!gGeoManager) {
  93     AliGeomManager::LoadGeometry(geoFilename);
  94     if (!gGeoManager) {
  95       Error("MUONRefit", "getting geometry from file %s failed", "generated/galice.root");
  96       exit(-1);
  97     }
  98   }
  99
 100   // load necessary data from OCDB
 101   AliCDBManager::Instance()->SetDefaultStorage(ocdbPath);
 102   AliCDBManager::Instance()->SetRun(runNumber);
 103   if (!AliMUONCDB::LoadField()) return;
 104   if (!AliMUONCDB::LoadMapping(kTRUE)) return;
 105
 106   // reconstruction parameters for the refitting
 107   AliMUONRecoParam* recoParam = AliMUONRecoParam::GetLowFluxParam();
 108   Info("MUONRefit", "\n Reconstruction parameters for refitting:");
 109   recoParam->Print("FULL");
 110
 111   AliMUONESDInterface esdInterface;
 112   AliMUONRefitter refitter(recoParam);
 113   refitter.Connect(&esdInterface);
 114   gRandom->SetSeed(1);
 115
 116   // timer start...
 117   TStopwatch timer;
 118
 119   // Loop over ESD events
 120   if (nevents > 0) nevents = TMath::Min(nevents,(Int_t)esdTree->GetEntries());
 121   else nevents = (Int_t)esdTree->GetEntries();
 122   for (Int_t iEvent = 0; iEvent < nevents; iEvent++) {
 123     if (printLevel>0) cout<<endl<<"            ****************event #"<<iEvent+1<<"****************"<<endl;
 124
 125     //----------------------------------------------//
 126     // -------------- process event --------------- //
 127     //----------------------------------------------//
 128     // get the ESD of current event
 129     esdTree->GetEvent(iEvent);
 130     if (!esd) {
 131       Error("MUONRefit", "no ESD object found for event %d", iEvent);
 132       return;
 133     }
 134     Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks();
 135     if (nTracks < 1) continue;
 136
 137     // load the current event
 138     esdInterface.LoadEvent(*esd, kFALSE);
 139
 140     // remove digits and clusters from ESD
 141     esd->FindListObject("MuonClusters")->Clear("C");
 142     esd->FindListObject("MuonPads")->Clear("C");
 143
 144     AliMUONVTrackStore* newTrackStore = 0x0;
 145     if (reconstructFromDigits) {
 146
 147       // loop over digits to modify their charge
 148       AliMUONVDigit *digit;
 149       TIter next(esdInterface.CreateDigitIterator());
 150       while ((digit = static_cast<AliMUONVDigit*>(next()))) {
 151         digit->SetCharge(digit->ADC());
 152         digit->Calibrated(kFALSE);
 153       }
 154
 155       // refit the tracks from digits
 156       refitter.SetFirstClusterIndex(0);
 157       newTrackStore = refitter.ReconstructFromDigits();
 158
 159     } else {
 160
 161       // refit the tracks from clusters
 162       newTrackStore = refitter.ReconstructFromClusters();
 163
 164     }
 165
 166     //----------------------------------------------//
 167     // ------ fill new ESD and print results ------ //
 168     //----------------------------------------------//
 169     // loop over the list of ESD tracks
 170     TClonesArray *esdTracks = (TClonesArray*) esd->FindListObject("MuonTracks");
 171     for (Int_t iTrack = 0; iTrack <  nTracks; iTrack++) {
 172
 173       // get the ESD track
 174       AliESDMuonTrack* esdTrack = (AliESDMuonTrack*) esdTracks->UncheckedAt(iTrack);
 175
 176       // skip ghost tracks (leave them unchanged in the new ESD file)
 177       if (!esdTrack->ContainTrackerData()) continue;
 178
 179       // get the corresponding MUON track
 180       AliMUONTrack* track = esdInterface.FindTrack(esdTrack->GetUniqueID());
 181
 182       // Find the corresponding re-fitted MUON track
 183       AliMUONTrack* newTrack = (AliMUONTrack*) newTrackStore->FindObject(esdTrack->GetUniqueID());
 184
 185       // replace the content of the current ESD track or remove it if the refitting has failed
 186       if (newTrack && (!recoParam->ImproveTracks() || newTrack->IsImproved())) {
 187
 188         // fill the track info
 189         Double_t vertex[3] = {esdTrack->GetNonBendingCoor(), esdTrack->GetBendingCoor(), esdTrack->GetZ()};
 190         AliMUONESDInterface::MUONToESD(*newTrack, *esdTrack, vertex);
 191
 192         // add the clusters (and the digits if previously there)
 193         for (Int_t i = 0; i < newTrack->GetNClusters(); i++) {
 194           AliMUONTrackParam *param = static_cast<AliMUONTrackParam*>(newTrack->GetTrackParamAtCluster()->UncheckedAt(i));
 195           AliMUONESDInterface::MUONToESD(*(param->GetClusterPtr()), *esd, esdInterface.GetDigits());
 196         }
 197
 198       } else esdTracks->Remove(esdTrack);
 199
 200       // print initial and re-fitted track parameters at first cluster if any
 201       if (printLevel>0) {
 202         cout<<"            ----------------track #"<<iTrack+1<<"----------------"<<endl;
 203         cout<<"before refit:"<<endl;
 204         AliMUONTrackParam *param = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->First();
 205         param->Print("FULL");
 206         if (printLevel>1) param->GetCovariances().Print();
 207         if (!newTrack) continue;
 208         cout<<"after refit:"<<endl;
 209         param = (AliMUONTrackParam*) newTrack->GetTrackParamAtCluster()->First();
 210         param->Print("FULL");
 211         if (printLevel>1) param->GetCovariances().Print();
 212         cout<<"            ----------------------------------------"<<endl;
 213       }
 214
 215     }
 216
 217     // free memory
 218     delete newTrackStore;
 219
 220     // fill new ESD tree with new tracks
 221     esdTracks->Compress();
 222     newESDTree->Fill();
 223
 224     if (printLevel>0) cout<<"            ****************************************"<<endl;
 225   }
 226
 227   // ...timer stop
 228   timer.Stop();
 229
 230   // write output ESD tree
 231   newESDFile->cd();
 232   newESDTree->Write();
 233   delete newESDTree;
 234   newESDFile->Close();
 235
 236   // free memory
 237   esdFile->Close();
 238   delete esd;
 239   delete recoParam;
 240
 241   cout<<endl<<"time to refit: R:"<<timer.RealTime()<<" C:"<<timer.CpuTime()<<endl<<endl;
 242 }
 243
 244 //-----------------------------------------------------------------------
 245 TTree* GetESDTree(TFile *esdFile)
 246 {
 247   /// Check that the file is properly open
 248   /// Return pointer to the ESD Tree
 249
 250   if (!esdFile || !esdFile->IsOpen()) {
 251     Error("GetESDTree", "opening ESD file failed");
 252     exit(-1);
 253   }
 254
 255   TTree* tree = (TTree*) esdFile->Get("esdTree");
 256   if (!tree) {
 257     Error("GetESDTree", "no ESD tree found");
 258     exit(-1);
 259   }
 260
 261   return tree;
 262
 263 }
 264