MUON/AliMUONReconstructor.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 /* $Id$ */
  16
  17 ///////////////////////////////////////////////////////////////////////////////
  18 //                                                                           //
  19 // class for MUON reconstruction                                             //
  20 //                                                                           //
  21 ///////////////////////////////////////////////////////////////////////////////
  22
  23 #include "AliMUONReconstructor.h"
  24
  25 #include "AliESD.h"
  26 #include "AliESDMuonTrack.h"
  27 #include "AliLog.h"
  28 #include "AliMUON.h"
  29 #include "AliMUONCalibrationData.h"
  30 #include "AliMUONClusterFinderAZ.h"
  31 #include "AliMUONClusterFinderVS.h"
  32 #include "AliMUONClusterReconstructor.h"
  33 #include "AliMUONData.h"
  34 #include "AliMUONDigitCalibrator.h"
  35 #include "AliMUONEventRecoCombi.h"
  36 #include "AliMUONRawReader.h"
  37 #include "AliMUONTrack.h"
  38 #include "AliMUONTrackParam.h"
  39 #include "AliMUONTrackReconstructor.h"
  40 #include "AliMUONTriggerTrack.h"
  41 #include "AliRawReader.h"
  42 #include "AliRun.h"
  43 #include "AliRunLoader.h"
  44 #include "TTask.h"
  45 #include "TStopwatch.h"
  46
  47 ClassImp(AliMUONReconstructor)
  48
  49 //_____________________________________________________________________________
  50 AliMUONReconstructor::AliMUONReconstructor()
  51   : AliReconstructor(), fCalibrationData(0x0)
  52 {
  53 /// Default constructor
  54
  55     AliDebug(1,"");
  56 }
  57
  58 //______________________________________________________________________________
  59 AliMUONReconstructor::AliMUONReconstructor(const AliMUONReconstructor& right)
  60   : AliReconstructor(right)
  61 {
  62 /// Protected copy constructor (not implemented)
  63
  64   AliFatal("Copy constructor not provided.");
  65 }
  66
  67 //_____________________________________________________________________________
  68 AliMUONReconstructor::~AliMUONReconstructor()
  69 {
  70 /// Destructor
  71
  72   AliDebug(1,"");
  73   delete fCalibrationData;
  74 }
  75
  76 //______________________________________________________________________________
  77 AliMUONReconstructor&
  78 AliMUONReconstructor::operator=(const AliMUONReconstructor& right)
  79 {
  80 /// Protected assignement operator (not implemented)
  81
  82   // check assignement to self
  83   if (this == &right) return *this;
  84
  85   AliFatal("Assignement operator not provided.");
  86
  87   return *this;
  88 }
  89
  90 //_____________________________________________________________________________
  91 TTask*
  92 AliMUONReconstructor::GetCalibrationTask(AliMUONData* data) const
  93 {
  94 /// Create the calibration task(s).
  95
  96   const AliRun* run = fRunLoader->GetAliRun();
  97
  98   // Not really clean, but for the moment we must check whether the
  99   // simulation has decalibrated the data or not...
 100   const AliMUON* muon = static_cast<AliMUON*>(run->GetModule("MUON"));
 101   if ( muon->DigitizerType().Contains("NewDigitizer") )
 102   {
 103     AliInfo("Calibration will occur.");
 104     Int_t runNumber = run->GetRunNumber();
 105     fCalibrationData = new AliMUONCalibrationData(runNumber);
 106     if ( !fCalibrationData->IsValid() )
 107     {
 108       AliError("Could not retrieve calibrations !");
 109       delete fCalibrationData;
 110       fCalibrationData = 0x0;
 111       return 0x0;
 112     }
 113     TTask* calibration = new TTask("MUONCalibrator","MUON Digit calibrator");
 114     calibration->Add(new AliMUONDigitCalibrator(data,fCalibrationData));
 115     //FIXME: calibration->Add(something about dead channels should go here).
 116     return calibration;
 117   }
 118   else
 119   {
 120     AliInfo("Detected the usage of old digitizer (w/o decalibration). "
 121             "Will not calibrate then.");
 122     return 0x0;
 123   }
 124 }
 125
 126 //_____________________________________________________________________________
 127 void
 128 AliMUONReconstructor::Init(AliRunLoader* runLoader)
 129 {
 130 /// Initialize
 131
 132   fRunLoader = runLoader;
 133 }
 134
 135 //_____________________________________________________________________________
 136 void AliMUONReconstructor::Reconstruct(AliRunLoader* runLoader) const
 137 {
 138 /// Reconstruct
 139 /// \todo add more
 140
 141   AliLoader* loader = runLoader->GetLoader("MUONLoader");
 142   Int_t nEvents = runLoader->GetNumberOfEvents();
 143
 144   AliMUONData* data = new AliMUONData(loader,"MUON","MUON");
 145
 146 // passing loader as argument.
 147   AliMUONTrackReconstructor* recoEvent = new AliMUONTrackReconstructor(loader, data);
 148
 149   if (strstr(GetOption(),"Original"))
 150     recoEvent->SetTrackMethod(1); // Original tracking
 151   else if (strstr(GetOption(),"Combi"))
 152     recoEvent->SetTrackMethod(3); // Combined cluster / track
 153   else
 154     recoEvent->SetTrackMethod(2); // Kalman
 155
 156   AliMUONClusterReconstructor* recoCluster = new AliMUONClusterReconstructor(data);
 157
 158   AliMUONClusterFinderVS *recModel = recoCluster->GetRecoModel();
 159
 160   if (!strstr(GetOption(),"VS")) {
 161     recModel = (AliMUONClusterFinderVS*) new AliMUONClusterFinderAZ();
 162     recoCluster->SetRecoModel(recModel);
 163   }
 164   recModel->SetGhostChi2Cut(10);
 165
 166   loader->LoadDigits("READ");
 167   loader->LoadRecPoints("RECREATE");
 168   loader->LoadTracks("RECREATE");
 169
 170   TTask* calibration = GetCalibrationTask(data);
 171
 172   Int_t chBeg = recoEvent->GetTrackMethod() == 3 ? 6 : 0;
 173   //   Loop over events
 174   for(Int_t ievent = 0; ievent < nEvents; ievent++) {
 175
 176     AliDebug(1,Form("Event %d",ievent));
 177
 178     runLoader->GetEvent(ievent);
 179
 180     //----------------------- digit2cluster & Trigger2Trigger -------------------
 181     if (!loader->TreeR()) loader->MakeRecPointsContainer();
 182
 183     // tracking branch
 184     if (recoEvent->GetTrackMethod() != 3) {
 185       data->MakeBranch("RC");
 186       data->SetTreeAddress("D,RC");
 187     } else {
 188       data->SetTreeAddress("D");
 189       data->SetTreeAddress("RCC");
 190     }
 191     // Important for avoiding a memory leak when reading digits ( to be investigated more in detail)
 192     // In any case the reading of GLT is needed for the Trigger2Tigger method below
 193     data->SetTreeAddress("GLT");
 194
 195     data->GetDigits();
 196
 197     if ( calibration )
 198     {
 199       calibration->ExecuteTask();
 200     }
 201
 202     recoCluster->Digits2Clusters(chBeg);
 203
 204     if (recoEvent->GetTrackMethod() == 3) {
 205       // Combined cluster / track finder
 206       AliMUONEventRecoCombi::Instance()->FillEvent(data, (AliMUONClusterFinderAZ*)recModel);
 207       ((AliMUONClusterFinderAZ*) recModel)->SetReco(2);
 208     }
 209     else data->Fill("RC");
 210
 211     // trigger branch
 212     data->MakeBranch("TC");
 213     data->SetTreeAddress("TC");
 214     recoCluster->Trigger2Trigger();
 215     data->Fill("TC");
 216
 217     //AZ loader->WriteRecPoints("OVERWRITE");
 218
 219     //---------------------------- Track & TriggerTrack ---------------------
 220     if (!loader->TreeT()) loader->MakeTracksContainer();
 221
 222     // trigger branch
 223     data->MakeBranch("RL"); //trigger track
 224     data->SetTreeAddress("RL");
 225     recoEvent->EventReconstructTrigger();
 226     data->Fill("RL");
 227
 228     // tracking branch
 229     data->MakeBranch("RT"); //track
 230     data->SetTreeAddress("RT");
 231     recoEvent->EventReconstruct();
 232     data->Fill("RT");
 233
 234     loader->WriteTracks("OVERWRITE");
 235
 236     if (recoEvent->GetTrackMethod() == 3) {
 237       // Combined cluster / track
 238       ((AliMUONClusterFinderAZ*) recModel)->SetReco(1);
 239       data->MakeBranch("RC");
 240       data->SetTreeAddress("RC");
 241       AliMUONEventRecoCombi::Instance()->FillRecP(data, recoEvent);
 242       data->Fill("RC");
 243     }
 244     loader->WriteRecPoints("OVERWRITE");
 245
 246     //--------------------------- Resetting branches -----------------------
 247     data->ResetDigits();
 248     data->ResetRawClusters();
 249     data->ResetTrigger();
 250
 251     data->ResetRawClusters();
 252     data->ResetTrigger();
 253     data->ResetRecTracks();
 254     data->ResetRecTriggerTracks();
 255
 256   }
 257   loader->UnloadDigits();
 258   loader->UnloadRecPoints();
 259   loader->UnloadTracks();
 260
 261   delete recoCluster;
 262   delete recoEvent;
 263   delete data;
 264   delete calibration;
 265 }
 266
 267 //_____________________________________________________________________________
 268 void AliMUONReconstructor::Reconstruct(AliRunLoader* runLoader, AliRawReader* rawReader) const
 269 {
 270 /// Recontruct
 271 /// \todo add more
 272
 273   //  AliLoader
 274   AliLoader* loader = runLoader->GetLoader("MUONLoader");
 275   AliMUONData data(loader,"MUON","MUON");
 276
 277   // passing loader as argument.
 278   AliMUONTrackReconstructor recoEvent(loader, &data);
 279
 280   AliMUONRawReader rawData(&data);
 281
 282   AliMUONClusterReconstructor recoCluster(&data);
 283
 284   if (strstr(GetOption(),"Original"))
 285   {
 286     recoEvent.SetTrackMethod(1); // Original tracking
 287   }
 288   else
 289   {
 290     recoEvent.SetTrackMethod(2); // Kalman
 291   }
 292
 293   AliMUONClusterFinderVS *recModel = recoCluster.GetRecoModel();
 294   if (!strstr(GetOption(),"VS"))
 295   {
 296     recModel = (AliMUONClusterFinderVS*) new AliMUONClusterFinderAZ();
 297     recoCluster.SetRecoModel(recModel);
 298   }
 299   recModel->SetGhostChi2Cut(10);
 300
 301   TTask* calibration = GetCalibrationTask(&data);
 302
 303   loader->LoadRecPoints("RECREATE");
 304   loader->LoadTracks("RECREATE");
 305   loader->LoadDigits("READ");
 306
 307   //   Loop over events
 308   Int_t iEvent = 0;
 309
 310   TStopwatch totalTimer;
 311   TStopwatch rawTimer;
 312   TStopwatch calibTimer;
 313   TStopwatch clusterTimer;
 314   TStopwatch trackingTimer;
 315
 316   rawTimer.Start(kTRUE); rawTimer.Stop();
 317   calibTimer.Start(kTRUE); calibTimer.Stop();
 318   clusterTimer.Start(kTRUE); clusterTimer.Stop();
 319   trackingTimer.Start(kTRUE); trackingTimer.Stop();
 320
 321   totalTimer.Start(kTRUE);
 322
 323   while (rawReader->NextEvent())
 324   {
 325     AliDebug(1,Form("Event %d",iEvent));
 326
 327     runLoader->GetEvent(iEvent++);
 328
 329     //----------------------- raw2digits & raw2trigger-------------------
 330     if (!loader->TreeD())
 331     {
 332       AliDebug(1,Form("Making Digit Container for event %d",iEvent));
 333       loader->MakeDigitsContainer();
 334     }
 335
 336     data.SetTreeAddress("D,GLT");
 337     rawTimer.Start(kFALSE);
 338     rawData.Raw2Digits(rawReader);
 339     rawTimer.Stop();
 340
 341     if ( calibration )
 342     {
 343       calibTimer.Start(kFALSE);
 344       calibration->ExecuteTask();
 345       calibTimer.Stop();
 346     }
 347
 348     //----------------------- digit2cluster & Trigger2Trigger -------------------
 349     clusterTimer.Start(kFALSE);
 350
 351     if (!loader->TreeR()) loader->MakeRecPointsContainer();
 352
 353     // tracking branch
 354     data.MakeBranch("RC");
 355     data.SetTreeAddress("RC");
 356     recoCluster.Digits2Clusters();
 357     data.Fill("RC");
 358
 359     // trigger branch
 360     data.MakeBranch("TC");
 361     data.SetTreeAddress("TC");
 362 //    recoCluster.Trigger2Trigger();
 363     data.Fill("TC");
 364
 365     loader->WriteRecPoints("OVERWRITE");
 366
 367     clusterTimer.Stop();
 368
 369     //---------------------------- Track & TriggerTrack ---------------------
 370     trackingTimer.Start(kFALSE);
 371     if (!loader->TreeT()) loader->MakeTracksContainer();
 372
 373     // trigger branch
 374     data.MakeBranch("RL"); //trigger track
 375     data.SetTreeAddress("RL");
 376     recoEvent.EventReconstructTrigger();
 377     data.Fill("RL");
 378
 379     // tracking branch
 380     data.MakeBranch("RT"); //track
 381     data.SetTreeAddress("RT");
 382     recoEvent.EventReconstruct();
 383     data.Fill("RT");
 384
 385     loader->WriteTracks("OVERWRITE");
 386     trackingTimer.Stop();
 387
 388     //--------------------------- Resetting branches -----------------------
 389     data.ResetDigits();
 390     data.ResetRawClusters();
 391     data.ResetTrigger();
 392
 393     data.ResetRawClusters();
 394     data.ResetTrigger();
 395     data.ResetRecTracks();
 396     data.ResetRecTriggerTracks();
 397
 398   }
 399
 400   totalTimer.Stop();
 401
 402   loader->UnloadRecPoints();
 403   loader->UnloadTracks();
 404   loader->UnloadDigits();
 405
 406   AliInfo(Form("Execution time for converting RAW data to digits in MUON : R:%.2fs C:%.2fs",
 407                rawTimer.RealTime(),rawTimer.CpuTime()));
 408   AliInfo(Form("Execution time for calibrating MUON : R:%.2fs C:%.2fs",
 409                calibTimer.RealTime(),calibTimer.CpuTime()));
 410   AliInfo(Form("Execution time for clusterizing MUON : R:%.2fs C:%.2fs",
 411                clusterTimer.RealTime(),clusterTimer.CpuTime()));
 412   AliInfo(Form("Execution time for tracking MUON : R:%.2fs C:%.2fs",
 413                trackingTimer.RealTime(),trackingTimer.CpuTime()));
 414   AliInfo(Form("Total Execution time for Reconstruct(from raw) MUON : R:%.2fs C:%.2fs",
 415                totalTimer.RealTime(),totalTimer.CpuTime()));
 416 }
 417
 418 //_____________________________________________________________________________
 419 void AliMUONReconstructor::FillESD(AliRunLoader* runLoader, AliESD* esd) const
 420 {
 421 /// Fill ESD
 422 /// \todo add more
 423
 424   TClonesArray* recTracksArray = 0;
 425   TClonesArray* recTrigTracksArray = 0;
 426
 427   AliLoader* loader = runLoader->GetLoader("MUONLoader");
 428   loader->LoadTracks("READ");
 429   AliMUONData* muonData = new AliMUONData(loader,"MUON","MUON");
 430
 431    // declaration
 432   Int_t iEvent;// nPart;
 433   Int_t nTrackHits;// nPrimary;
 434   Double_t fitFmin;
 435
 436   Double_t bendingSlope, nonBendingSlope, inverseBendingMomentum;
 437   Double_t xRec, yRec, zRec, chi2MatchTrigger;
 438   Bool_t matchTrigger;
 439
 440   // setting pointer for tracks, triggertracks & trackparam at vertex
 441   AliMUONTrack* recTrack = 0;
 442   AliMUONTrackParam* trackParam = 0;
 443   AliMUONTriggerTrack* recTriggerTrack = 0;
 444 //   TParticle* particle = new TParticle();
 445 //   AliGenEventHeader* header = 0;
 446   iEvent = runLoader->GetEventNumber();
 447   runLoader->GetEvent(iEvent);
 448
 449   // Get vertex
 450   Double_t vertex[3] = {0};
 451   const AliESDVertex *esdVert = esd->GetVertex();
 452   if (esdVert) esdVert->GetXYZ(vertex);
 453
 454   //  nPrimary = 0;
 455 //   if ( (header = runLoader->GetHeader()->GenEventHeader()) ) {
 456 //     header->PrimaryVertex(vertex);
 457 //   } else {
 458 //     runLoader->LoadKinematics("READ");
 459 //     runLoader->TreeK()->GetBranch("Particles")->SetAddress(&particle);
 460 //     nPart = (Int_t)runLoader->TreeK()->GetEntries();
 461 //     for(Int_t iPart = 0; iPart < nPart; iPart++) {
 462 //       runLoader->TreeK()->GetEvent(iPart);
 463 //       if (particle->GetFirstMother() == -1) {
 464 //      vertex[0] += particle->Vx();
 465 //      vertex[1] += particle->Vy();
 466 //      vertex[2] += particle->Vz();
 467 //      nPrimary++;
 468 //       }
 469 //       if (nPrimary) {
 470 //      vertex[0] /= (double)nPrimary;
 471 //      vertex[1] /= (double)nPrimary;
 472 //      vertex[2] /= (double)nPrimary;
 473 //       }
 474 //     }
 475 //   }
 476   // setting ESD MUON class
 477   AliESDMuonTrack* theESDTrack = new  AliESDMuonTrack() ;
 478
 479   //-------------------- trigger tracks-------------
 480   Long_t trigPat = 0;
 481   muonData->SetTreeAddress("RL");
 482   muonData->GetRecTriggerTracks();
 483   recTrigTracksArray = muonData->RecTriggerTracks();
 484
 485   // ready global trigger pattern from first track
 486   if (recTrigTracksArray)
 487     recTriggerTrack = (AliMUONTriggerTrack*) recTrigTracksArray->First();
 488   if (recTriggerTrack) trigPat = recTriggerTrack->GetGTPattern();
 489
 490   //printf(">>> Event %d Number of Recconstructed tracks %d \n",iEvent, nrectracks);
 491
 492   // -------------------- tracks-------------
 493   muonData->SetTreeAddress("RT");
 494   muonData->GetRecTracks();
 495   recTracksArray = muonData->RecTracks();
 496
 497   Int_t nRecTracks = 0;
 498   if (recTracksArray)
 499     nRecTracks = (Int_t) recTracksArray->GetEntriesFast(); //
 500
 501   // loop over tracks
 502   for (Int_t iRecTracks = 0; iRecTracks <  nRecTracks;  iRecTracks++) {
 503
 504     // reading info from tracks
 505     recTrack = (AliMUONTrack*) recTracksArray->At(iRecTracks);
 506
 507     trackParam = (AliMUONTrackParam*) (recTrack->GetTrackParamAtHit())->First();
 508     trackParam->ExtrapToVertex(vertex[0],vertex[1],vertex[2]);
 509
 510     bendingSlope            = trackParam->GetBendingSlope();
 511     nonBendingSlope         = trackParam->GetNonBendingSlope();
 512     inverseBendingMomentum = trackParam->GetInverseBendingMomentum();
 513     xRec  = trackParam->GetNonBendingCoor();
 514     yRec  = trackParam->GetBendingCoor();
 515     zRec  = trackParam->GetZ();
 516
 517     nTrackHits       = recTrack->GetNTrackHits();
 518     fitFmin          = recTrack->GetFitFMin();
 519     matchTrigger     = recTrack->GetMatchTrigger();
 520     chi2MatchTrigger = recTrack->GetChi2MatchTrigger();
 521
 522     // setting data member of ESD MUON
 523     theESDTrack->SetInverseBendingMomentum(inverseBendingMomentum);
 524     theESDTrack->SetThetaX(TMath::ATan(nonBendingSlope));
 525     theESDTrack->SetThetaY(TMath::ATan(bendingSlope));
 526     theESDTrack->SetZ(zRec);
 527     theESDTrack->SetBendingCoor(yRec); // calculate vertex at ESD or Tracking level ?
 528     theESDTrack->SetNonBendingCoor(xRec);
 529     theESDTrack->SetChi2(fitFmin);
 530     theESDTrack->SetNHit(nTrackHits);
 531     theESDTrack->SetMatchTrigger(matchTrigger);
 532     theESDTrack->SetChi2MatchTrigger(chi2MatchTrigger);
 533
 534     // storing ESD MUON Track into ESD Event
 535     if (nRecTracks != 0)
 536       esd->AddMuonTrack(theESDTrack);
 537   } // end loop tracks
 538
 539   // reset muondata
 540   muonData->ResetRecTracks();
 541   muonData->ResetRecTriggerTracks();
 542
 543   //} // end loop on event
 544   loader->UnloadTracks();
 545  //  if (!header)
 546 //     runLoader->UnloadKinematics();
 547   delete theESDTrack;
 548   delete muonData;
 549   // delete particle;
 550 }//_____________________________________________________________________________
 551 void AliMUONReconstructor::FillESD(AliRunLoader* runLoader, AliRawReader* /*rawReader*/, AliESD* esd) const
 552 {
 553 /// Fill ESD
 554 /// \todo add more
 555
 556   // don't need rawReader ???
 557   FillESD(runLoader, esd);
 558 }