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