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 "AliMagF.h"
  36 #include "AliTracker.h"
  37 #include "AliESDEvent.h"
  38 #include "AliESDMuonTrack.h"
  39 #include "AliRecoParam.h"
  40 #include "AliCDBManager.h"
  41 #include "AliCDBEntry.h"
  42 #include "AliCDBPath.h"
  43 #include "AliGeomManager.h"
  44
  45 // MUON includes
  46 #include "AliMpCDB.h"
  47 #include "AliMUONRecoParam.h"
  48 #include "AliMUONESDInterface.h"
  49 #include "AliMUONRefitter.h"
  50 #include "AliMUONVDigit.h"
  51 #include "AliMUONTrack.h"
  52 #include "AliMUONVTrackStore.h"
  53 #include "AliMUONTrackParam.h"
  54 #include "AliMUONTrackExtrap.h"
  55 #endif
  56
  57 const Int_t printLevel = 1;
  58
  59 void Prepare();
  60 TTree* GetESDTree(TFile *esdFile);
  61
  62 //-----------------------------------------------------------------------
  63 void MUONRefit(Int_t nevents = -1, const char* esdFileNameIn = "AliESDs.root", const char* esdFileNameOut = "AliESDs_New.root")
  64 {
  65   /// refit ESD tracks from ESD pads (i.e. re-clusterized the attached ESD clusters);
  66   /// reset the charge of the digit using their raw charge before refitting;
  67   /// compare results with original ESD tracks;
  68   /// write results in a new ESD file
  69
  70   // prepare the refitting
  71   gRandom->SetSeed(1);
  72   Prepare();
  73
  74   // reconstruction parameters for the refitting
  75   AliMUONRecoParam* recoParam = AliMUONRecoParam::GetLowFluxParam();
  76   Info("MUONRefit", "\n Reconstruction parameters for refitting:");
  77   recoParam->Print("FULL");
  78
  79   AliMUONESDInterface esdInterface;
  80   AliMUONRefitter refitter(recoParam);
  81   refitter.Connect(&esdInterface);
  82
  83   // open the ESD file and tree
  84   TFile* esdFile = TFile::Open(esdFileNameIn);
  85   TTree* esdTree = GetESDTree(esdFile);
  86
  87   // create the ESD output file and tree
  88   TFile* newESDFile = TFile::Open(esdFileNameOut, "RECREATE");
  89   newESDFile->SetCompressionLevel(2);
  90   TTree* newESDTree = esdTree->CloneTree(0);
  91
  92   // connect ESD event to the ESD tree
  93   AliESDEvent* esd = new AliESDEvent();
  94   esd->ReadFromTree(esdTree);
  95
  96   // timer start...
  97   TStopwatch timer;
  98
  99   // Loop over ESD events
 100   if (nevents > 0) nevents = TMath::Min(nevents,(Int_t)esdTree->GetEntries());
 101   else nevents = (Int_t)esdTree->GetEntries();
 102   for (Int_t iEvent = 0; iEvent < nevents; iEvent++) {
 103     if (printLevel>0) cout<<endl<<"            ****************event #"<<iEvent+1<<"****************"<<endl;
 104
 105     //----------------------------------------------//
 106     // -------------- process event --------------- //
 107     //----------------------------------------------//
 108     // get the ESD of current event
 109     esdTree->GetEvent(iEvent);
 110     if (!esd) {
 111       Error("MUONRefit", "no ESD object found for event %d", iEvent);
 112       return;
 113     }
 114     Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks();
 115     if (nTracks < 1) continue;
 116
 117     // load the current event
 118     esdInterface.LoadEvent(*esd);
 119
 120     // loop over digit to modify their charge
 121     AliMUONVDigit *digit;
 122     TIter next(esdInterface.CreateDigitIterator());
 123     while ((digit = static_cast<AliMUONVDigit*>(next()))) {
 124       digit->SetCharge(digit->ADC());
 125       digit->Calibrated(kFALSE);
 126     }
 127
 128     // refit the tracks from digits
 129     AliMUONVTrackStore* newTrackStore = refitter.ReconstructFromDigits();
 130
 131     //----------------------------------------------//
 132     // ------ fill new ESD and print results ------ //
 133     //----------------------------------------------//
 134     // loop over the list of ESD tracks
 135     TClonesArray *esdTracks = (TClonesArray*) esd->FindListObject("MuonTracks");
 136     for (Int_t iTrack = 0; iTrack <  nTracks; iTrack++) {
 137
 138       // get the ESD track
 139       AliESDMuonTrack* esdTrack = (AliESDMuonTrack*) esdTracks->UncheckedAt(iTrack);
 140
 141       // skip ghost tracks (leave them unchanged in the new ESD file)
 142       if (!esdTrack->ContainTrackerData()) continue;
 143
 144       // get the corresponding MUON track
 145       AliMUONTrack* track = esdInterface.FindTrack(esdTrack->GetUniqueID());
 146
 147       // Find the corresponding re-fitted MUON track
 148       AliMUONTrack* newTrack = (AliMUONTrack*) newTrackStore->FindObject(esdTrack->GetUniqueID());
 149
 150       // replace the content of the current ESD track or remove it if the refitting has failed
 151       if (newTrack) {
 152         Double_t vertex[3] = {esdTrack->GetNonBendingCoor(), esdTrack->GetBendingCoor(), esdTrack->GetZ()};
 153         AliMUONESDInterface::MUONToESD(*newTrack, *esdTrack, vertex, esdInterface.GetDigits());
 154       } else {
 155         esdTracks->Remove(esdTrack);
 156       }
 157
 158       // print initial and re-fitted track parameters at first cluster if any
 159       if (printLevel>0) {
 160         cout<<"            ----------------track #"<<iTrack+1<<"----------------"<<endl;
 161         cout<<"before refit:"<<endl;
 162         AliMUONTrackParam *param = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->First();
 163         param->Print("FULL");
 164         if (printLevel>1) param->GetCovariances().Print();
 165         if (!newTrack) continue;
 166         cout<<"after refit:"<<endl;
 167         param = (AliMUONTrackParam*) newTrack->GetTrackParamAtCluster()->First();
 168         param->Print("FULL");
 169         if (printLevel>1) param->GetCovariances().Print();
 170         cout<<"            ----------------------------------------"<<endl;
 171       }
 172
 173     }
 174
 175     // free memory
 176     delete newTrackStore;
 177
 178     // fill new ESD tree with new tracks
 179     esdTracks->Compress();
 180     newESDTree->Fill();
 181
 182     if (printLevel>0) cout<<"            ****************************************"<<endl;
 183   }
 184
 185   // ...timer stop
 186   timer.Stop();
 187
 188   // write output ESD tree
 189   newESDFile->cd();
 190   newESDTree->Write();
 191   delete newESDTree;
 192   newESDFile->Close();
 193
 194   // free memory
 195   esdFile->Close();
 196   delete esd;
 197
 198   cout<<endl<<"time to refit: R:"<<timer.RealTime()<<" C:"<<timer.CpuTime()<<endl<<endl;
 199 }
 200
 201 //-----------------------------------------------------------------------
 202 void Prepare()
 203 {
 204   /// Set the geometry, the magnetic field, the mapping and the reconstruction parameters
 205
 206   // Import TGeo geometry (needed by AliMUONTrackExtrap::ExtrapToVertex)
 207   if (!gGeoManager) {
 208     AliGeomManager::LoadGeometry("geometry.root");
 209     if (!gGeoManager) {
 210       Error("MUONRefit", "getting geometry from file %s failed", "generated/galice.root");
 211       return;
 212     }
 213   }
 214
 215   // set  mag field
 216   // waiting for mag field in CDB
 217   if (!TGeoGlobalMagField::Instance()->GetField()) {
 218     printf("Loading field map...\n");
 219     AliMagF* field = new AliMagF("Maps","Maps",2,1.,1., 10.,AliMagF::k5kG);
 220     TGeoGlobalMagField::Instance()->SetField(field);
 221   }
 222   // set the magnetic field for track extrapolations
 223   AliMUONTrackExtrap::SetField();
 224
 225   // Load mapping
 226   AliCDBManager* man = AliCDBManager::Instance();
 227   man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
 228   man->SetRun(0);
 229   if ( ! AliMpCDB::LoadDDLStore() ) {
 230     Error("MUONRefit","Could not access mapping from OCDB !");
 231     exit(-1);
 232   }
 233
 234   // Load initial reconstruction parameters from OCDB
 235   AliMUONRecoParam* recoParam = 0x0;
 236   AliCDBPath path("MUON","Calib","RecoParam");
 237   AliCDBEntry *entry=man->Get(path.GetPath());
 238   if(entry) {
 239     recoParam = dynamic_cast<AliMUONRecoParam*>(entry->GetObject());
 240     entry->SetOwner(0);
 241     AliCDBManager::Instance()->UnloadFromCache(path.GetPath());
 242   }
 243   if (!recoParam) {
 244     printf("Couldn't find RecoParam object in OCDB: create default one");
 245     recoParam = AliMUONRecoParam::GetLowFluxParam();
 246   }
 247   Info("MUONRefit", "\n initial recontruction parameters:");
 248   recoParam->Print("FULL");
 249   AliMUONESDInterface::ResetTracker(recoParam);
 250
 251 }
 252
 253 //-----------------------------------------------------------------------
 254 TTree* GetESDTree(TFile *esdFile)
 255 {
 256   /// Check that the file is properly open
 257   /// Return pointer to the ESD Tree
 258
 259   if (!esdFile || !esdFile->IsOpen()) {
 260     Error("GetESDTree", "opening ESD file failed");
 261     exit(-1);
 262   }
 263
 264   TTree* tree = (TTree*) esdFile->Get("esdTree");
 265   if (!tree) {
 266     Error("GetESDTree", "no ESD tree found");
 267     exit(-1);
 268   }
 269
 270   return tree;
 271
 272 }
 273