STEER/AliReconstruction.cxx

   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 ///////////////////////////////////////////////////////////////////////////////
  19 //                                                                           //
  20 // class for running the reconstruction                                      //
  21 //                                                                           //
  22 // Clusters and tracks are created for all detectors and all events by       //
  23 // typing:                                                                   //
  24 //                                                                           //
  25 //   AliReconstruction rec;                                                  //
  26 //   rec.Run();                                                              //
  27 //                                                                           //
  28 // The Run method returns kTRUE in case of successful execution.             //
  29 // The name of the galice file can be changed from the default               //
  30 // "galice.root" by                                                          //
  31 //                                                                           //
  32 //   rec.SetGAliceFile("...");                                               //
  33 //                                                                           //
  34 // The reconstruction can be switched on or off for individual detectors by  //
  35 //                                                                           //
  36 //   rec.SetRunReconstruction("...");                                        //
  37 //                                                                           //
  38 // The argument is a (case sensitive) string with the names of the           //
  39 // detectors separated by a space. The special string "ALL" selects all      //
  40 // available detectors. This is the default.                                 //
  41 //                                                                           //
  42 // The tracking in ITS, TPC and TRD and the creation of ESD tracks can be    //
  43 // switched off by                                                           //
  44 //                                                                           //
  45 //   rec.SetRunTracking(kFALSE);                                             //
  46 //                                                                           //
  47 // The filling of additional ESD information can be steered by               //
  48 //                                                                           //
  49 //   rec.SetFillESD("...");                                                  //
  50 //                                                                           //
  51 // Again, the string specifies the list of detectors. The default is "ALL".  //
  52 //                                                                           //
  53 // The reconstruction requires digits as input. For the creation of digits   //
  54 // have a look at the class AliSimulation.                                   //
  55 //                                                                           //
  56 ///////////////////////////////////////////////////////////////////////////////
  57
  58
  59 #include "AliReconstruction.h"
  60 #include "AliRunLoader.h"
  61 #include "AliRun.h"
  62 #include "AliModule.h"
  63 #include "AliDetector.h"
  64 #include "AliTracker.h"
  65 #include "AliESD.h"
  66 #include "AliHeader.h"
  67 #include "AliGenEventHeader.h"
  68 #include "AliESDpid.h"
  69 #include <TArrayF.h>
  70
  71
  72 ClassImp(AliReconstruction)
  73
  74
  75 //_____________________________________________________________________________
  76 AliReconstruction::AliReconstruction(const char* name, const char* title) :
  77   TNamed(name, title)
  78 {
  79 // create reconstruction object with default parameters
  80
  81   Init();
  82 }
  83
  84 //_____________________________________________________________________________
  85 AliReconstruction::AliReconstruction(const AliReconstruction& rec) :
  86   TNamed(rec)
  87 {
  88 // copy constructor
  89
  90   fRunReconstruction = rec.fRunReconstruction;
  91   fRunTracking = rec.fRunTracking;
  92   fStopOnError = rec.fStopOnError;
  93
  94   fGAliceFileName = rec.fGAliceFileName;
  95
  96   fRunLoader = NULL;
  97 }
  98
  99 //_____________________________________________________________________________
 100 AliReconstruction& AliReconstruction::operator = (const AliReconstruction& rec)
 101 {
 102 // assignment operator
 103
 104   this->~AliReconstruction();
 105   new(this) AliReconstruction(rec);
 106   return *this;
 107 }
 108
 109 //_____________________________________________________________________________
 110 AliReconstruction::~AliReconstruction()
 111 {
 112 // clean up
 113
 114 }
 115
 116 //_____________________________________________________________________________
 117 void AliReconstruction::Init()
 118 {
 119 // set default parameters
 120
 121   fRunReconstruction = "ALL";
 122   fRunTracking = kTRUE;
 123   fFillESD = "ALL";
 124   fStopOnError = kFALSE;
 125
 126   fGAliceFileName = "galice.root";
 127
 128   fRunLoader = NULL;
 129 }
 130
 131
 132 //_____________________________________________________________________________
 133 void AliReconstruction::SetGAliceFile(const char* fileName)
 134 {
 135 // set the name of the galice file
 136
 137   fGAliceFileName = fileName;
 138 }
 139
 140
 141 //_____________________________________________________________________________
 142 Bool_t AliReconstruction::Run()
 143 {
 144 // run the reconstruction
 145
 146   // open the run loader
 147   if (fRunLoader) delete fRunLoader;
 148   fRunLoader = AliRunLoader::Open(fGAliceFileName.Data());
 149   if (!fRunLoader) {
 150     Error("Run", "no run loader found in file %s",
 151           fGAliceFileName.Data());
 152     return kFALSE;
 153   }
 154   fRunLoader->LoadgAlice();
 155   gAlice = fRunLoader->GetAliRun();
 156   if (!gAlice) {
 157     Error("Run", "no gAlice object found in file %s",
 158           fGAliceFileName.Data());
 159     return kFALSE;
 160   }
 161
 162   // local reconstruction
 163   if (!fRunReconstruction.IsNull()) {
 164     if (!RunReconstruction(fRunReconstruction)) {
 165       if (fStopOnError) return kFALSE;
 166     }
 167   }
 168   if (!fRunTracking && fFillESD.IsNull()) return kTRUE;
 169
 170   // get loaders and trackers
 171   fITSLoader = fRunLoader->GetLoader("ITSLoader");
 172   if (!fITSLoader) {
 173     Error("Run", "no ITS loader found");
 174     if (fStopOnError) return kFALSE;
 175   }
 176   fITSTracker = NULL;
 177   if (gAlice->GetDetector("ITS")) {
 178     fITSTracker = gAlice->GetDetector("ITS")->CreateTracker();
 179   }
 180   if (!fITSTracker) {
 181     Error("Run", "couldn't create a tracker for ITS");
 182     if (fStopOnError) return kFALSE;
 183   }
 184
 185   fTPCLoader = fRunLoader->GetLoader("TPCLoader");
 186   if (!fTPCLoader) {
 187     Error("Run", "no TPC loader found");
 188     if (fStopOnError) return kFALSE;
 189   }
 190   fTPCTracker = NULL;
 191   if (gAlice->GetDetector("TPC")) {
 192     fTPCTracker = gAlice->GetDetector("TPC")->CreateTracker();
 193   }
 194   if (!fTPCTracker) {
 195     Error("Run", "couldn't create a tracker for TPC");
 196     if (fStopOnError) return kFALSE;
 197   }
 198
 199   fTRDLoader = fRunLoader->GetLoader("TRDLoader");
 200   if (!fTRDLoader) {
 201     Error("Run", "no TRD loader found");
 202     if (fStopOnError) return kFALSE;
 203   }
 204   fTRDTracker = NULL;
 205   if (gAlice->GetDetector("TRD")) {
 206     fTRDTracker = gAlice->GetDetector("TRD")->CreateTracker();
 207   }
 208   if (!fTRDTracker) {
 209     Error("Run", "couldn't create a tracker for TRD");
 210     if (fStopOnError) return kFALSE;
 211   }
 212
 213   fTOFLoader = fRunLoader->GetLoader("TOFLoader");
 214   if (!fTOFLoader) {
 215     Error("Run", "no TOF loader found");
 216     if (fStopOnError) return kFALSE;
 217   }
 218   fTOFTracker = NULL;
 219   if (gAlice->GetDetector("TOF")) {
 220     fTOFTracker = gAlice->GetDetector("TOF")->CreateTracker();
 221   }
 222   if (!fTOFTracker) {
 223     Error("Run", "couldn't create a tracker for TOF");
 224     if (fStopOnError) return kFALSE;
 225   }
 226
 227   // create the ESD output file
 228   TFile* file = TFile::Open("AliESDs.root", "RECREATE");
 229   if (!file->IsOpen()) {
 230     Error("Run", "opening AliESDs.root failed");
 231     if (fStopOnError) return kFALSE;
 232   }
 233
 234   // loop over events
 235   for (Int_t iEvent = 0; iEvent < fRunLoader->GetNumberOfEvents(); iEvent++) {
 236     Info("Run", "processing event %d", iEvent);
 237     AliESD* esd = new AliESD;
 238     fRunLoader->GetEvent(iEvent);
 239     esd->SetRunNumber(gAlice->GetRunNumber());
 240     esd->SetEventNumber(gAlice->GetEvNumber());
 241
 242     // barrel tracking
 243     if (fRunTracking) {
 244       if (!RunTracking(esd)) {
 245         if (fStopOnError) return kFALSE;
 246       }
 247     }
 248
 249     // fill ESD
 250     if (!fFillESD.IsNull()) {
 251       if (!FillESD(esd, fFillESD)) {
 252         if (fStopOnError) return kFALSE;
 253       }
 254     }
 255
 256     // combined PID
 257     AliESDpid::MakePID(esd);
 258
 259     // write ESD
 260     char name[100];
 261     sprintf(name, "ESD%d", iEvent);
 262     file->cd();
 263     if (!esd->Write(name)) {
 264       Error("Run", "writing ESD failed");
 265       if (fStopOnError) return kFALSE;
 266     }
 267   }
 268
 269   file->Close();
 270
 271   return kTRUE;
 272 }
 273
 274
 275 //_____________________________________________________________________________
 276 Bool_t AliReconstruction::RunReconstruction(const TString& detectors)
 277 {
 278 // run the reconstruction
 279
 280   TString detStr = detectors;
 281   TObjArray* detArray = gAlice->Detectors();
 282   for (Int_t iDet = 0; iDet < detArray->GetEntriesFast(); iDet++) {
 283     AliModule* det = (AliModule*) detArray->At(iDet);
 284     if (!det || !det->IsActive()) continue;
 285     if (IsSelected(det->GetName(), detStr)) {
 286       Info("RunReconstruction", "running reconstruction for %s",
 287            det->GetName());
 288       det->Reconstruct();
 289     }
 290   }
 291
 292   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
 293     Error("RunReconstruction", "the following detectors were not found: %s",
 294           detStr.Data());
 295     if (fStopOnError) return kFALSE;
 296   }
 297
 298   return kTRUE;
 299 }
 300
 301 //_____________________________________________________________________________
 302 Bool_t AliReconstruction::RunTracking(AliESD* esd)
 303 {
 304 // run the barrel tracking
 305
 306   // get the primary vertex (from MC for the moment)
 307   TArrayF vertex(3);
 308   fRunLoader->GetHeader()->GenEventHeader()->PrimaryVertex(vertex);
 309   Double_t vtxPos[3] = {vertex[0], vertex[1], vertex[2]};
 310   Double_t vtxCov[6] = {
 311     0.005,
 312     0.000, 0.005,
 313     0.000, 0.000, 0.010
 314   };
 315   Double_t vtxErr[3] = {vtxCov[0], vtxCov[2], vtxCov[5]}; // diag. elements
 316   esd->SetVertex(vtxPos, vtxCov);
 317   fITSTracker->SetVertex(vtxPos, vtxErr);
 318   fTPCTracker->SetVertex(vtxPos, vtxErr);
 319   fTRDTracker->SetVertex(vtxPos, vtxErr);
 320
 321   // TPC inward fit
 322   Info("RunTracking", "TPC inward fit");
 323   fTPCLoader->LoadRecPoints("read");
 324   TTree* tpcTree = fTPCLoader->TreeR();
 325   if (!tpcTree) {
 326     Error("RunTracking", "Can't get the TPC cluster tree");
 327     return kFALSE;
 328   }
 329   fTPCTracker->LoadClusters(tpcTree);
 330   if (fTPCTracker->Clusters2Tracks(esd) != 0) {
 331     Error("RunTracking", "TPC inward track fit failed");
 332     return kFALSE;
 333   }
 334
 335   // ITS inward fit
 336   Info("RunTracking", "ITS inward fit");
 337   fITSLoader->LoadRecPoints("read");
 338   TTree* itsTree = fITSLoader->TreeR();
 339   if (!itsTree) {
 340     Error("RunTracking", "Can't get the ITS cluster tree");
 341     return kFALSE;
 342   }
 343   fITSTracker->LoadClusters(itsTree);
 344   if (fITSTracker->Clusters2Tracks(esd) != 0) {
 345     Error("RunTracking", "ITS inward track fit failed");
 346     return kFALSE;
 347   }
 348
 349   // ITS back propagation
 350   Info("RunTracking", "ITS back propagation");
 351   if (fITSTracker->PropagateBack(esd) != 0) {
 352     Error("RunTracking", "ITS backward propagation failed");
 353     return kFALSE;
 354   }
 355
 356   // TPC back propagation
 357   Info("RunTracking", "TPC back propagation");
 358   if (fTPCTracker->PropagateBack(esd) != 0) {
 359     Error("RunTracking", "TPC backward propagation failed");
 360     return kFALSE;
 361   }
 362
 363   // TRD back propagation
 364   Info("RunTracking", "TRD back propagation");
 365   fTRDLoader->LoadRecPoints("read");
 366   TTree* trdTree = fTRDLoader->TreeR();
 367   if (!trdTree) {
 368     Error("RunTracking", "Can't get the TRD cluster tree");
 369     return kFALSE;
 370   }
 371   fTRDTracker->LoadClusters(trdTree);
 372   if (fTRDTracker->PropagateBack(esd) != 0) {
 373     Error("RunTracking", "TRD backward propagation failed");
 374     return kFALSE;
 375   }
 376
 377   // TOF back propagation
 378   Info("RunTracking", "TOF back propagation");
 379   fTOFLoader->LoadDigits("read");
 380   TTree* tofTree = fTOFLoader->TreeD();
 381   if (!tofTree) {
 382     Error("RunTracking", "Can't get the TOF digits tree");
 383     return kFALSE;
 384   }
 385   fTOFTracker->LoadClusters(tofTree);
 386   if (fTOFTracker->PropagateBack(esd) != 0) {
 387     Error("RunTracking", "TOF backward propagation failed");
 388     return kFALSE;
 389   }
 390   fTOFTracker->UnloadClusters();
 391   fTOFLoader->UnloadDigits();
 392
 393   // TRD inward refit
 394   Info("RunTracking", "TRD inward refit");
 395   if (fTRDTracker->RefitInward(esd) != 0) {
 396     Error("RunTracking", "TRD inward refit failed");
 397     return kFALSE;
 398   }
 399   fTRDTracker->UnloadClusters();
 400   fTRDLoader->UnloadRecPoints();
 401
 402   // TPC inward refit
 403   Info("RunTracking", "TPC inward refit");
 404   if (fTPCTracker->RefitInward(esd) != 0) {
 405     Error("RunTracking", "TPC inward refit failed");
 406     return kFALSE;
 407   }
 408   fTPCTracker->UnloadClusters();
 409   fTPCLoader->UnloadRecPoints();
 410
 411   // ITS inward refit
 412   Info("RunTracking", "ITS inward refit");
 413   if (fITSTracker->RefitInward(esd) != 0) {
 414     Error("RunTracking", "ITS inward refit failed");
 415     return kFALSE;
 416   }
 417   fITSTracker->UnloadClusters();
 418   fITSLoader->UnloadRecPoints();
 419
 420   return kTRUE;
 421 }
 422
 423 //_____________________________________________________________________________
 424 Bool_t AliReconstruction::FillESD(AliESD* esd, const TString& detectors)
 425 {
 426 // fill the event summary data
 427
 428   TString detStr = detectors;
 429   TObjArray* detArray = gAlice->Detectors();
 430   for (Int_t iDet = 0; iDet < detArray->GetEntriesFast(); iDet++) {
 431     AliModule* det = (AliModule*) detArray->At(iDet);
 432     if (!det || !det->IsActive()) continue;
 433     if (IsSelected(det->GetName(), detStr)) {
 434       Info("FillESD", "filling ESD for %s",
 435            det->GetName());
 436       det->FillESD(esd);
 437     }
 438   }
 439
 440   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
 441     Error("FillESD", "the following detectors were not found: %s",
 442           detStr.Data());
 443     if (fStopOnError) return kFALSE;
 444   }
 445
 446   return kTRUE;
 447 }
 448
 449
 450 //_____________________________________________________________________________
 451 Bool_t AliReconstruction::IsSelected(TString detName, TString& detectors) const
 452 {
 453 // check whether detName is contained in detectors
 454 // if yes, it is removed from detectors
 455
 456   // check if all detectors are selected
 457   if ((detectors.CompareTo("ALL") == 0) ||
 458       detectors.BeginsWith("ALL ") ||
 459       detectors.EndsWith(" ALL") ||
 460       detectors.Contains(" ALL ")) {
 461     detectors = "ALL";
 462     return kTRUE;
 463   }
 464
 465   // search for the given detector
 466   Bool_t result = kFALSE;
 467   if ((detectors.CompareTo(detName) == 0) ||
 468       detectors.BeginsWith(detName+" ") ||
 469       detectors.EndsWith(" "+detName) ||
 470       detectors.Contains(" "+detName+" ")) {
 471     detectors.ReplaceAll(detName, "");
 472     result = kTRUE;
 473   }
 474
 475   // clean up the detectors string
 476   while (detectors.Contains("  ")) detectors.ReplaceAll("  ", " ");
 477   while (detectors.BeginsWith(" ")) detectors.Remove(0, 1);
 478   while (detectors.EndsWith(" ")) detectors.Remove(detectors.Length()-1, 1);
 479
 480   return result;
 481 }