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