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